基于單方向阻尼源項的數值消波方法

王凡瑜，權曉波，魏海鵬，孔德才

（北京宇航系統工程研究所，北京 100076）

波浪是水利工程、船舶與海洋工程等領域關心的重要環境因素之一。研究波浪及其影響主要有實驗與數值模擬兩種手段。數值模擬手段具有成本低、可研究的空間尺度大等優點，并能夠較簡便地模擬更接近真實水文環境的各種不規則波，在工程實踐中發揮著重要作用。

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

圖1 平均自由面抬升現象

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

1 數值方法

1.1 控制方程

對常黏性系數牛頓流體的二維不可壓縮流動，有連續性方程：

動量方程：

式中：u、v分別為水平速度分量與豎直速度分量；ρ為密度；p為壓強；μ為動力黏性系數；g為重力加速度。

氣水自由面采用VOF模型[11]模擬，相分數隱式求解；瞬態項采用一階隱式格式離散，壓力項采用PRESTO!格式離散，對流項采用二階迎風格式離散；控制方程組采用SIMPLE算法[12]求解。

1.2 造波方法

本文采用邊界條件法造波。邊界條件造波法根據波浪解析解或數值解設置入口邊界的速度條件與相分數條件，具有計算量小、易生成各種波浪的優點。

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

液相速度分量為：

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

1.3 消波方法

1.3.1 人工阻尼消波法 人工阻尼消波法以多孔介質模擬自然界中的沙灘或實驗水槽中的消波端。假設消波區填充有多孔介質，為模擬對流體運動的阻礙作用，向式(2)右端添加阻尼源項。

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

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

使阻尼源項自消波區起始端從零逐漸增大，以避免在消波區起始端產生反射波。式中xs、xe分別為消波區起止端水平坐標。

1.3.2 單方向阻尼源項 多孔介質對流體運動的阻尼作用使消波區內流體趨于靜止，在避免出口處產生反射波的同時，也阻礙流體通過出口邊界流出。采用邊界條件造波法時，液相流體在入口邊界存在凈流入量，計算域內液相流體總量不斷增加，使得平均自由面高度不斷抬升[6]。

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

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

2 算例設置

計算域為長400 m、高35 m的矩形區域，如圖2所示，其中靜水深25 m。左側邊界設置為速度入口，右側與頂部邊界設置為壓力出口，底部邊界設置為無滑移壁面。時間步長取0.005 s。

圖2 計算域示意圖

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

表1 計算條件

3 結果與分析

3.1 單方向阻尼源項改進效果

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

圖3 不同方向性算例的數值波面對比

為定量分析，取觀察區間內平均波峰高度與平均波谷高度的算術平均值表征平均自由面高度：

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

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

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

3.1.2 波浪模擬精度 為評價波浪模擬精度，定義自由面高度數值解與解析解的偏差平方和為：

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

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

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

3.2 阻尼系數影響

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

圖5 不同阻尼系數算例的偏差平方和對比

3.3 密度依賴性影響

比較采用不同密度依賴性阻尼源項的偏差平方和，如圖6所示。圖中黑色實心條柱對應密度依賴型阻尼源項，紅色斜紋條柱對應密度無關型阻尼源項。與采用不同阻尼系數相比，采用不同密度依賴性的偏差平方和相差不大，說明密度依賴性對波浪模擬精度的影響小于阻尼系數。

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

盡管兩者差異較小，但仍可以看出：小波陡條件下，密度依賴型阻尼源項的偏差平方和較小；大波陡條件下，密度無關型阻尼源項的偏差平方和較小。阻尼源項依賴密度與否主要影響消波區內氣相阻尼大小，其中密度無關型阻尼源項產生較大的氣相阻尼。如前所述，大波陡波浪能量較強，消波所需阻尼較大，故大波陡條件下密度無關型阻尼源項波浪模擬精度較高。相反，小波陡條件下波浪能量較弱，消波所需阻尼較小，密度依賴型阻尼源項模擬精度較高。

4 結 論

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