王凡瑜，權(quán)曉波，魏海鵬，孔德才

（北京宇航系統(tǒng)工程研究所，北京 100076）

波浪是水利工程、船舶與海洋工程等領(lǐng)域關(guān)心的重要環(huán)境因素之一。研究波浪及其影響主要有實(shí)驗(yàn)與數(shù)值模擬兩種手段。數(shù)值模擬手段具有成本低、可研究的空間尺度大等優(yōu)點(diǎn)，并能夠較簡(jiǎn)便地模擬更接近真實(shí)水文環(huán)境的各種不規(guī)則波，在工程實(shí)踐中發(fā)揮著重要作用。

數(shù)值消波方法用于消除計(jì)算域出口附近產(chǎn)生的反射波，避免反射波影響上游流場(chǎng)，是開展波浪數(shù)值模擬的關(guān)鍵點(diǎn)之一。其中人工阻尼消波法通過修改控制方程模擬自然界中沙灘或?qū)嶒?yàn)水槽中消波端的阻尼作用實(shí)現(xiàn)消波，具有實(shí)現(xiàn)簡(jiǎn)便、適用于非線性波數(shù)值模擬的優(yōu)點(diǎn)，在波浪數(shù)值模擬中得到了廣泛應(yīng)用。ISRAELI M等[1]首先提出人工阻尼消波法，其后LARSEN J等[2]以及BAKER G R等[3]分別在Boussinesq方程和二維邊界元模型中實(shí)現(xiàn)了人工阻尼消波。國(guó)內(nèi)方面，高學(xué)平等[4]結(jié)合人工阻尼消波法和輻射邊界條件消波法模擬了二階Stokes波、駐波與不規(guī)則波，董志等[5]基于人工阻尼消波法比較了不同造波方法的二階Stokes波模擬效果，査晶晶[6]在開源計(jì)算流體力學(xué)軟件OpenFOAM中植入了人工阻尼消波法，王巧莎等[7]研究了適用于不同重力加速度條件與不同波陡條件的人工阻尼消波法，穆澤楠[8]基于人工阻尼消波法比較了不同造波方法的不規(guī)則波模擬效果。但是一些研究[6,9-10]指出，波浪數(shù)值模擬中存在如圖1所示的平均自由面緩慢抬升現(xiàn)象，數(shù)值解波形逐漸向上偏離解析解，使波浪模擬精度略有損失，給開展波浪數(shù)值模擬帶來一定困擾。

圖1 平均自由面抬升現(xiàn)象

本文提出了一種基于單方向阻尼源項(xiàng)的數(shù)值消波方法，采用VOF模型與邊界條件造波法模擬了大、小兩種波陡條件的線性波，通過比較獲得的數(shù)值波面并定量分析平均自由面抬升程度，驗(yàn)證了單方向阻尼源項(xiàng)抑制平均自由面抬升的效果。在此基礎(chǔ)上，進(jìn)一步研究了阻尼系數(shù)與密度依賴性對(duì)波浪數(shù)值模擬精度的影響。

1 數(shù)值方法

1.1 控制方程

對(duì)常黏性系數(shù)牛頓流體的二維不可壓縮流動(dòng)，有連續(xù)性方程：

動(dòng)量方程：

式中：u、v分別為水平速度分量與豎直速度分量；ρ為密度；p為壓強(qiáng)；μ為動(dòng)力黏性系數(shù)；g為重力加速度。

氣水自由面采用VOF模型[11]模擬，相分?jǐn)?shù)隱式求解；瞬態(tài)項(xiàng)采用一階隱式格式離散，壓力項(xiàng)采用PRESTO!格式離散，對(duì)流項(xiàng)采用二階迎風(fēng)格式離散；控制方程組采用SIMPLE算法[12]求解。

1.2 造波方法

本文采用邊界條件法造波。邊界條件造波法根據(jù)波浪解析解或數(shù)值解設(shè)置入口邊界的速度條件與相分?jǐn)?shù)條件，具有計(jì)算量小、易生成各種波浪的優(yōu)點(diǎn)。

以線性波為例，波浪理論指出，自由面高度為：

液相速度分量為：

式中：H為波高；k為波數(shù)；ω為角頻率；T為周期；d為靜水水深。入口邊界的相分?jǐn)?shù)條件可由式(3)確定，速度條件可由式(4)確定。

1.3 消波方法

1.3.1 人工阻尼消波法 人工阻尼消波法以多孔介質(zhì)模擬自然界中的沙灘或?qū)嶒?yàn)水槽中的消波端。假設(shè)消波區(qū)填充有多孔介質(zhì)，為模擬對(duì)流體運(yùn)動(dòng)的阻礙作用，向式(2)右端添加阻尼源項(xiàng)。

式中：μu/α和μv/α反映多孔介質(zhì)導(dǎo)致的黏性損失，α為多孔介質(zhì)滲透率；和反映多孔介質(zhì)導(dǎo)致的慣性損失，C為慣性阻力系數(shù)。實(shí)際使用中一般只取μu/α和μv/α，并改寫為[6,8,13]：

式中：μ0為阻尼系數(shù)；f(x)為阻尼分布函數(shù)。通常采用線性阻尼分布函數(shù)：

使阻尼源項(xiàng)自消波區(qū)起始端從零逐漸增大，以避免在消波區(qū)起始端產(chǎn)生反射波。式中xs、xe分別為消波區(qū)起止端水平坐標(biāo)。

1.3.2 單方向阻尼源項(xiàng) 多孔介質(zhì)對(duì)流體運(yùn)動(dòng)的阻尼作用使消波區(qū)內(nèi)流體趨于靜止，在避免出口處產(chǎn)生反射波的同時(shí)，也阻礙流體通過出口邊界流出。采用邊界條件造波法時(shí)，液相流體在入口邊界存在凈流入量，計(jì)算域內(nèi)液相流體總量不斷增加，使得平均自由面高度不斷抬升[6]。

平均自由面抬升的根本原因是計(jì)算域內(nèi)液相流體總量不斷增加，故考慮取消水平方向阻尼，即令

以減少流體水平速度衰減，使液相流體能夠通過出口邊界流出，降低計(jì)算域內(nèi)液相流體凈增量。為盡可能避免出口處產(chǎn)生反射波影響上游流場(chǎng)，仍保留豎直方向阻尼，使自由面波動(dòng)在消波區(qū)內(nèi)逐漸衰減。因僅有豎直方向阻尼源項(xiàng)，所以稱作基于單方向阻尼源項(xiàng)的數(shù)值消波方法，相應(yīng)地將式(6)所示的阻尼源項(xiàng)稱作兩方向阻尼源項(xiàng)。1.3.3 密度依賴性 一些研究[5,7]采用常系數(shù)代替阻尼源項(xiàng)中的流體密度，本文稱之為密度無(wú)關(guān)型阻尼源項(xiàng)，該系數(shù)數(shù)值上等于液相流體密度。相應(yīng)地將式(6)所示的阻尼源項(xiàng)稱作密度依賴型阻尼源項(xiàng)。

2 算例設(shè)置

計(jì)算域?yàn)殚L(zhǎng)400 m、高35 m的矩形區(qū)域，如圖2所示，其中靜水深25 m。左側(cè)邊界設(shè)置為速度入口，右側(cè)與頂部邊界設(shè)置為壓力出口，底部邊界設(shè)置為無(wú)滑移壁面。時(shí)間步長(zhǎng)取0.005 s。

圖2 計(jì)算域示意圖

分別在大、小波陡條件下采用各種形式的阻尼源項(xiàng)開展波浪數(shù)值模擬。小波陡條件下，波高H=1 m，波長(zhǎng)L=50 m，消波區(qū)長(zhǎng)度Ld=50 m；大波陡條件下，波高H=2 m，波長(zhǎng)L=40 m，消波區(qū)長(zhǎng)度Ld=80 m。計(jì)算條件如表1所示。

表1 計(jì)算條件

3 結(jié)果與分析

3.1 單方向阻尼源項(xiàng)改進(jìn)效果

3.1.1 平均自由面抬升 比較采用不同方向性阻尼源項(xiàng)獲得的數(shù)值波面，如圖3所示，黑色虛線對(duì)應(yīng)兩方向阻尼源項(xiàng)，紅色實(shí)線對(duì)應(yīng)單方向阻尼源項(xiàng)。與兩方向阻尼源項(xiàng)相比，采用單方向阻尼源項(xiàng)獲得的數(shù)值波面在小波陡條件下略有降低，在大波陡條件下降低明顯。

圖3 不同方向性算例的數(shù)值波面對(duì)比

為定量分析，取觀察區(qū)間內(nèi)平均波峰高度與平均波谷高度的算術(shù)平均值表征平均自由面高度：

式中：m、n分別為觀察區(qū)間內(nèi)波峰或波谷的個(gè)數(shù)；yi,cj、yi,tj分別表示算例i中第j個(gè)波峰或波谷的高度。觀察區(qū)間定義為2～5倍波長(zhǎng)范圍，即小波陡條件下取為100～250 m，大波陡條件下取為80～200 m。

比較采用不同方向性阻尼源項(xiàng)獲得的平均自由面高度，如表2所示。相較于兩方向阻尼源項(xiàng)，采用單方向阻尼源項(xiàng)可使平均自由面高度下降60%～80%，抑制平均自由面抬升效果顯著。

表2 不同方向性算例的平均自由面高度對(duì)比

3.1.2 波浪模擬精度 為評(píng)價(jià)波浪模擬精度，定義自由面高度數(shù)值解與解析解的偏差平方和為：

式中：η表示自由面高度；n表示數(shù)值解；a表示解析解；i=1,2…表示觀察區(qū)間內(nèi)沿水平方向分布的各離散點(diǎn)。偏差平方和越小，表示模擬精度越高。

不同方向性阻尼源項(xiàng)的偏差平方和對(duì)比如圖4所示，黑色實(shí)心條柱對(duì)應(yīng)兩方向阻尼源項(xiàng)，紅色斜紋條柱對(duì)應(yīng)單方向阻尼源項(xiàng)。保持其他條件不變，采用單方向阻尼源項(xiàng)的偏差平方和普遍較小，表明采用單方向阻尼源項(xiàng)抑制平均自由面抬升，可提高波浪模擬精度，獲得更接近波浪理論解析解的數(shù)值波面。

圖4 不同方向性算例的偏差平方和對(duì)比

3.2 阻尼系數(shù)影響

比較采用不同系數(shù)阻尼源項(xiàng)的偏差平方和，如圖5所示。圖中黑色實(shí)心條柱對(duì)應(yīng)μ0=1算例，紅色斜紋條柱對(duì)應(yīng)μ0=10算例。小波陡條件下，取μ0=1的偏差平方和較??；大波陡條件下，取μ0=10的偏差平方和較小。阻尼系數(shù)對(duì)波浪模擬精度的影響在大、小波陡條件下相反，原因在于：大波陡條件下，波浪能量較強(qiáng)，消波所需阻尼較大，故取μ0=10模擬精度較高；小波陡條件下相反，波浪能量較弱，取μ0=1產(chǎn)生的阻尼較合適，模擬精度較高。

圖5 不同阻尼系數(shù)算例的偏差平方和對(duì)比

3.3 密度依賴性影響

比較采用不同密度依賴性阻尼源項(xiàng)的偏差平方和，如圖6所示。圖中黑色實(shí)心條柱對(duì)應(yīng)密度依賴型阻尼源項(xiàng)，紅色斜紋條柱對(duì)應(yīng)密度無(wú)關(guān)型阻尼源項(xiàng)。與采用不同阻尼系數(shù)相比，采用不同密度依賴性的偏差平方和相差不大，說明密度依賴性對(duì)波浪模擬精度的影響小于阻尼系數(shù)。

圖6 不同密度依賴性算例的偏差平方和對(duì)比

盡管兩者差異較小，但仍可以看出：小波陡條件下，密度依賴型阻尼源項(xiàng)的偏差平方和較??；大波陡條件下，密度無(wú)關(guān)型阻尼源項(xiàng)的偏差平方和較小。阻尼源項(xiàng)依賴密度與否主要影響消波區(qū)內(nèi)氣相阻尼大小，其中密度無(wú)關(guān)型阻尼源項(xiàng)產(chǎn)生較大的氣相阻尼。如前所述，大波陡波浪能量較強(qiáng)，消波所需阻尼較大，故大波陡條件下密度無(wú)關(guān)型阻尼源項(xiàng)波浪模擬精度較高。相反，小波陡條件下波浪能量較弱，消波所需阻尼較小，密度依賴型阻尼源項(xiàng)模擬精度較高。

4 結(jié) 論

基于VOF模型與邊界條件造波法，采用不同方向性阻尼源項(xiàng)在大、小波陡條件下開展了波浪數(shù)值模擬，定量比較了平均自由面抬升程度，并進(jìn)一步研究了阻尼系數(shù)與密度依賴性對(duì)波浪模擬精度的影響，結(jié)果分析表明：（1）采用基于單方向阻尼源項(xiàng)的數(shù)值消波方法，可有效抑制平均自由面抬升，提高數(shù)值模擬精度，獲得更接近理論解的數(shù)值波面；（2）阻尼系數(shù)的優(yōu)取值與波陡有關(guān)，波陡較大時(shí)，消波所需阻尼較強(qiáng)，阻尼系數(shù)取值偏大，波浪模擬精度較高；相反，波陡較小時(shí)，阻尼系數(shù)取值也應(yīng)較小；（3）密度依賴性對(duì)波浪模擬精度的影響弱于阻尼系數(shù)，小波陡條件下，密度依賴型阻尼源項(xiàng)模擬精度較高；大波陡條件下，密度無(wú)關(guān)型阻尼源項(xiàng)模擬精度較高。