基于OpenFOAM幾何流體體積方法的波浪數值模擬

田 康，張 堯,李金龍，張新曙,3，尤云祥,3

(1.上海交通大學 海洋工程國家重點實驗室，上海 200240；2.自然資源部海洋減災中心， 北京 100194；3.上海交通大學 高新船舶與深海開發裝備協同創新中心，上海 200240)

在對波浪和結構物相互作用問題進行模擬時，如研究海洋平臺立柱的波浪高階力，船舶在波浪中運動的操縱性、耐波性等問題，為了保證數值結果的可靠性，需對波浪進行準確模擬.此外，在對波浪進行機理性研究時，如波浪非線性、波浪的傳播與演化、波浪破碎等，也對波浪模擬提出了很高要求.

波浪模擬需要解決3個關鍵問題：造波、自由液面捕捉和消波.造波方法主要有3種：速度入口造波[1]、動邊界造波[2]和源造波[3].速度入口法又稱為邊界輸入法，即根據波浪理論，在每一時間步指定造波邊界上速度場分布；動邊界法又名仿物理造波法，即模仿模型試驗中的推板造波法，通過邊界運動推動計算域內的流體產生波浪；源造波法是一種在計算域內部產生波浪的方法，根據源項添加方式不同，可以分為動量源法和質量源法.

自由液面捕捉方法可以分為拉格朗日法和歐拉法.拉格朗日法從自由面處粒子運動入手，通過描述粒子運動軌跡來刻畫自由液面變化，Harlow等[4]在1965年提出的MAC (Mark and Cell)方法就是其中一種代表方法.歐拉法關注自由液面附近某一固定空間區域流體流動變化，代表方法是Hirt 等[5]在1981年提出的流體體積(VOF)方法.根據求解過程的不同，VOF方法可細分為代數VOF方法和幾何VOF方法，兩者的區別在于：幾何VOF方法會對自由液面進行幾何重構，然后利用重構的自由液面更新體積分數率分布；代數VOF方法直接對體積分數率輸運方程進行求解，不涉及任何幾何對象的操作.傳統幾何VOF方法只適用于結構化網格，在推廣到三維非結構化網格時遇到了困難，為了能夠將幾何VOF方法擴展到三維任意多面體網格，Roenby等[6]于2016年提出了一種新型幾何VOF方法——isoAdvector方法.

消波主要有兩種思路：出流邊界消波和數值消波.出流邊界是指采用人工截斷邊界，讓流體自然地流出計算域邊界而對計算內的求解沒有任何影響，在眾多方法中，Sommerfield輻射邊界條件的效果最好[7].但是由于該條件在時間空間上是局部的，僅限于頻率已知的規則入射波和長波情況[8].目前更為流行的方法是數值消波方法，數值消波方法分為3種：阻尼消波、主動式吸收消波和松弛區方法消波.阻尼消波通過在消波區域對自由邊界條件或動量方程添加阻尼耗散項，從而達到消波效果；主動式吸收消波的原理是在水池末端設置推板，通過推板運動制造一個與反射波波幅相同但相位相反的波去抵消反射波，與此同時對入射波的傳播沒有任何影響；松弛區方法[9]消波是近幾年提出的一種新型消波方法，其原理是在計算值和理論值之間合理分配權重，通過解析解與計算流體動力學(CFD)解之間的過渡實現消波.Chen等[10]在用OpenFOAM研究波浪-結構物相互作用問題時，使用了松弛區方法進行消波；此外，Jensen等[11]在對近海結構物多孔介質方程及阻力系數問題進行研究時也用到了松弛區消波方法.

本文研究基于開源流體力學計算平臺OpenFOAM對波浪進行數值模擬.利用OpenFOAM已有幾何VOF方法求解器InterIsoFoam(內置isoAdvector方法[6])，通過二次開發將InterIsoFoam與松弛區方法(Relaxation Zone Scheme)[12]消波結合，生成新的求解器WaveIsoFoam，對波浪進行模擬.首先采用幾何VOF模擬波浪，并進行網格和庫朗數收斂性分析，然后將計算結果分別與代數VOF方法及理論值進行了比較和分析，最后分析了松弛區方法中消波區長度和消波權重對波浪模擬的影響.

1 數學模型

在本文波浪模擬問題中，將氣液兩相視為單一不可壓縮連續介質，并考慮流體黏性.控制方程、VOF方法和松弛區消波方法的介紹如下.

1.1 控制方程

在計算域中流體流動同時滿足質量守恒定律和動量守恒定律，對應控制方程中連續性方程和Navier-Stokes方程的積分形式如下：

(1)

(2)

式中：V為計算域中任意網格；?V代表組成網格的面的集合；ni,nj為網格表面對應單位法向量分量，其方向為指向網格外部；ui,uj為流體速度分量；ρ為流體密度；xi,xj,xk為笛卡爾坐標系下3個分量；μ為動力黏度；S為網格單元表面面積；pd為剔除靜壓之后的壓力，總壓力p的表達式為p=pd+ρgkxk；重力加速度g=[gigjgk]T=[0 0 -9.81]Tm/s2；源項fσ,i為表面張力，根據連續表面張力理論[13-14]，fσ,i可以用下式進行表達：

(3)

式中：表面張力系數一般取σ=0.078 2 kg/s2[14]；κ為界面曲率；α1為液體體積分數.

1.2 VOF方法

VOF方法通過引入流體體積分數對自由液面進行捕捉.流體體積分數定義為目標流體的體積與網格體積比值，基于這個函數在每個網格體心上的值，可以實現對運動的自由液面捕捉.在本文中，將液體和氣體分別用下標1和2表示.兩相流體域中單一流體密度ρ，動力黏度μ和運動黏度ν可以用下列公式表達：

ρ=α1ρ1+α2ρ2

(4)

μ=α1μ1+α2μ2

(5)

(6)

式中：ρ1,ρ2分別為液體和氣體的密度；μ1,μ2分別為液體和氣體的動力黏度；α1,α2分別為液體和氣體的體積分數.

1.2.1isoAdvector方法 isoAdvector方法是Roenby等[6]在2016年提出的一種新型幾何VOF方法.在每個時間步長內，isoAdvector方法的求解過程可以分為以下兩步.

(1) 界面幾何重構.當某一個網格單元所包含的液體體積分數在0和1之間時，該網格單元就稱為界面網格單元，即自由液面經過這個網格單元.在界面幾何重構中，關鍵是找出每個界面網格單元中，最適合劃分液體和氣體的平面(Isoface).Isoface需滿足的條件為：從該網格單元(母網格單元)劃分出的兩個子網格單元的幾何體積比值應與該網格單元體心所儲存的液體體積分數匹配.搜尋最佳劃分Isoface的具體方法可以參考文獻[6].如圖1所示，在相鄰時間步長中，母網格內的Isoface會隨著流體的輸運而向某一方向移動.這些由Isoface所重構的界面將用于下一步體積分數更新.由于在尋找切分Isoface時是單獨對每一個網格單元進行計算處理，所以重構界面可能并不連續[6].

圖1 Isoface在任意多面體網格單元中移動示意圖Fig.1 Sketch of movement of Isoface within a polyhedral mesh

(2) 體積分數率更新.下一時刻體積分數場可由下式[15]得出：

(7)

(8)

1.2.2代數VOF方法 代數VOF方法需要對體積分數輸運方程進行求解，為獲取更薄的界面厚度(定義液體體積分數為0.01～0.99的厚度)，在代數VOF方法中引入了人工壓縮項.在OpenFOAM中，代數VOF方法采用MULES (Multidimensional Universal Limiter for Explicit Solution) 方法進行求解.MULES方法引入了FCT算法(Flux-corrected Transport Algorithm) 求解體積分數輸運方程：

(9)

式中：等號左端最后一項為人工壓縮項；uc,j的處理詳見文獻[15].

1.3 松弛區方法

松弛區方法的原理為：引入權重因子w，對松弛區內流體速度和體積分數進行加權，實現目標解析解和CFD計算解的過渡，從而實現速度入口造波邊界和CFD計算域間的過渡，以及CFD計算域與壓力出口邊界的過渡，對應的表達式如下：

(10)

α=(1-w)αA+wαC

(11)

在入口處的松弛區，w從0以指定形式增長到1，實現解析解到CFD解的過渡；出口處的松弛區，w從1以指定的形式減小到0，實現CFD解到解析解的過渡，達到消波的效果，如圖2所示.權重w的基本表達式為

圖2 松弛區方法示意圖Fig.2 Sketch of relaxation zone scheme

w=1-w*Co/Comax

(12)

式中：Co,Comax分別為局部庫朗數和最大庫朗數；w*有3種不同的形式，分別對應3種不同類型的權重分布，即指數權重、三階多項式權重和自由多項式權重，

(13)

w*=-2(1-γ)3+3(1-γ)2

(14)

w*=1-γp

(15)

γ為松弛區域內的局部坐標，坐標位置取決于松弛區域的形狀.需要說明的是，在所有情況中，松弛權重w*都是局部坐標γ下的函數，滿足w*(γ=1)=0，w*(γ=0)=1；式(13)中指數p一般取3.5；式(15)中指數p一般取為整數.

2 結果與討論

2.1 計算模型

本文基于OpenFOAM對波浪進行數值模擬.波浪模擬計算域如圖3所示，全局笛卡爾坐標系(xi,xj,xk)原點O位于靜水面、入口處垂直壁面和中縱剖面三面的相交處，xi軸方向與波浪傳播方向一致，xk軸方向垂直向上，xj軸方向滿足右手定則.計算域入口處及頂部邊界條件為速度入口，出口處邊界條件為壓力出口；計算域底部為固壁邊界條件，前后兩側為對稱邊界條件.在計算域入口處和出口處均設置了松弛區，其中入口處松弛區長度為兩個波長；出口處的松弛區長度設置對消波效果有影響，具體長度選取見后文關于消波問題的分析.在計算域自由液面處沿xi軸均勻分布320個波高監測點，用于獲取波浪在數值模擬過程中的空間數據和時歷數據.計算域長度為20 m,寬度為1 m，高度為3 m，計算域水深為2.5 m,實際計算水深為5 m.

圖3 三維計算域示意圖(m)Fig.3 3D sketch of computational domain (m)

本研究數值模擬應用開源軟件ESI OpenFOAM-v1812，基于層流理論[16-17]，使用WaveIsoFoam對兩相流進行求解.WaveIsoFoam是在interIsoFoam求解器基礎上引入松弛區功能修改得到.InterIsoFoam內置了isoAdvector方法對自由液面進行捕捉，因此代表了幾何VOF方法求解兩相流.

時間離散格式采用二階Crank-Nicolson方法，對流項采用高分辨率Gauss vanLeer求解，Laplace項采用高斯線性修正方法進行求解；控制方程解耦采用SIMPLE方法，動量方程用Gauss-Seidel方法求解.在每個時間步長內，有兩次壓力修正，分別應用GAMG求解器和PCG求解器對這兩次壓力修正時的壓力泊松方程進行求解.

2.2 收斂性分析

在實際對波浪相關問題進行模擬研究時，為了盡可能減小由于數值計算產生的誤差，同時盡可能提高計算效率，需先進行網格和庫朗數收斂性分析.為了探究幾何VOF方法在收斂性方面的特性，本文選用了無限水深Stokes五階波浪進行分析，波浪模擬參數見表1.

表1 Stokes五階波波浪模擬參數Tab.1 Simulation parameters of Stokes fifth-order waves

2.2.1網格收斂性分析 為了能夠較好地模擬波浪傳播，需要在自由表面處進行網格加密.一般來說,在橫向、縱向和垂向3個方向中，垂向網格密度對波浪模擬的影響尤為明顯.因此，在自由液面處，固定垂向網格尺寸與縱向(波浪傳播方向)網格尺寸比為1∶2，在垂向選取4組網格H/Δz=6,8,12,16(H為波高，Δz為垂向網格尺寸).圖4所示為模擬時間為30 s時不同網格密度下波高的空間分布，圖中：xi為沿計算域縱向方向的距離；η為波高.由圖可知，當網格密度H/Δz>8時，幾何VOF方法能較好地對波浪進行模擬；而當H/Δz<8時，幾何VOF方法模擬波浪在波高和相位上隨著波浪傳播逐漸與理論解產生偏離.因此，在使用幾何VOF方法對波浪進行模擬時，單個波高垂向分布網格層數至少應為8層.

圖4 不同網格密度下波高的空間分布Fig.4 Spatial distributions of wave elevations at different mesh densities

2.2.2庫朗數收斂性分析 庫朗數定義如下：

(16)

式中：Δt為時間步長；|U|為某個網格單元內速度矢量的模，Δx為速度方向網格長度.在計算域網格已劃定的情況下，網格長度Δx確定；在同一工況下，固定位置的速度矢量的模|U|也確定.因此，可以用庫朗數Co對CFD計算的時間步長Δt進行反映.用于收斂性分析所選取的庫朗數Co=0.1,0.2,0.4,1.0，圖5為模擬時間為30 s時不同庫朗數下波高空間分布.

圖5 不同庫朗數下波高的空間分布Fig.5 Spatial distributions of wave elevations at different Courant numbers

計算結果顯示，不同庫朗數下波浪模擬結果均能較好地與理論結果相吻合.即在使用幾何VOF方法進行波浪模擬時，選取較大的庫朗數或時間步長，也能得到較為滿意的結果.

2.3 造波分析

圖6為模擬時間為30 s時不同波陡下幾何和代數VOF方法波高模擬空間分布，其中波高為H、波長為λ、波陡為δ.圖7為模擬時間為30 s時圖6對應的波高模擬時歷曲線.從圖6和圖7中可以看出：

圖6 不同波陡下幾何和代數VOF方法波高模擬空間分布Fig.6 Spatial dstributions of wave elevations by using geometrical and algebraic VOF methods at varied wave steepnesses

圖7 不同波陡下幾何和代數VOF方法波高模擬時歷曲線Fig.7 Time histories of wave elevations by using geometrical and algebraic VOF method at varied wave steepnesses

(1) 對于波高的空間分布和時歷模擬，在不同波陡情況下，相比于代數VOF方法，幾何VOF方法都能得到更準確的結果；在相同設置下，代數VOF方法模擬得到的波高則稍大于理論值得到的波高；

(2) 波浪在空間域傳播時，無論是幾何VOF方法還是代數VOF方法，產生的波浪均會與理論結果產生微小的空間相位差，且隨著模擬波浪的波高增大，空間相位差也逐漸累積變大；波浪在時間域中，隨著模擬時間增加，模擬結果與理論結果無相位差產生，同時不受模擬波浪波高變化的影響.

圖8為模擬時間為30 s時不同波頻下幾何和代數VOF方法波高模擬空間分布，其中波頻用ω表示.比較結果顯示，對于不同的波浪頻率，幾何VOF方法均能給出較好的波浪模擬結果.相比之下代數VOF方法在不同波浪頻率情況下模擬波浪的空間相位與理論值吻合較好，但幅值略大于理論值.

圖8 不同波頻下幾何和代數VOF方法波高模擬空間分布Fig.8 Spatial dstributions of wave elevations by using geometrical and algebraic VOF method at varied wave frequencies

幾何VOF方法相較于代數VOF能更加精確地模擬波浪，其中最主要的原因在于幾何VOF方法在求解體積分數之前對網格進行了重構，使其能夠更加精確地自由液面進行捕捉.圖9為幾何VOF方法和代數VOF方法分別對應的自由液面處液體體積分數分布情況，藍色部分代表α1=0 (氣體)，紅色部分代表α1=1 (液體)；黑線α1=0.01，黃線α1=0.99.用液體體積分數為0.01和0.99兩個分界之間的厚度以表示自由液面捕捉精度，可以看到幾何VOF方法所得的自由液面處液體分布更為集中，即界面厚度比代數VOF方法更小.這體現了幾何VOF方法較代數VOF方法有更高的自由液面捕捉精度，即能更好地對波浪進行模擬.

圖9 自由液面處液體體積分數α1分布情況Fig.9 Distributions of the liquid volume fraction α1 around free surface

2.4 消波分析

為了能夠模擬出質量較高的波浪，除了需要選擇合適的網格密度、庫朗數和自由液面捕捉方法以外，消波方法的選擇也非常重要.若消波效果不好，可能會導致波浪反射從而對入射波造成干擾，使得波浪模擬結果與理論結果發生偏離.

在前面介紹中已經提到，本次研究所采用的消波方法為松弛區方法.在松弛區方法中，消波區長度和消波加權方法會對消波效果有較為明顯的影響.因此，在本次研究中，分別對這兩個因素進行研究和分析.

2.4.1消波區長度分析 在消波區長度對消波效果影響的研究分析中，選取1倍波長(2 m)、2倍波長(4 m)和3倍波長(6 m)進行模擬計算，消波權重采用指數方法，即式(13).圖10為模擬時間為30 s時不同消波區長度(用Lwaz表示)下波高空間分布，兩條紅色垂向虛線間的距離為波高削減至零所需最短空間長度.

圖10 不同消波區長度下波高空間分布 Fig.10 Spatial distributions of wave elevations at different lengths of relaxation zone

從圖10中可以看出，隨著消波區長度增加，波幅削減至0所需最短空間長度也隨之增加，這是由計算解至解析解過渡區域延長所致.但是最短空間長度占總消波區長度比值較小，說明松弛區方法可以快速地對波高進行抑制.

而在波浪反射方面，從圖10(a)中看出，當消波區長度為一個波長時，靠近消波區附近的波高模擬結果明顯高于理論值，說明波浪反射未得到完全遏制.而2倍和3倍波長的消波區長度，不僅能夠實現快速消波，并且能很好地遏制波浪反射.

2.4.2松弛區消波權重分析 松弛區消波方法有3種權重控制流體速度和體積分數率從邊界到CFD計算域的過渡，分別是：指數方法(式(13))、三階多項式方法(式(14))以及自由多項式方法(式(15)).圖11為選取消波區長度為2倍波長、模擬時間為30 s時不同消波權重下波高空間分布.比較發現，相比于指數方法進行權重分配，三階多項式及自由多項式方法，均能在非常短的區域內，將波幅削減至0；但波幅模擬結果都大于理論值，說明兩種權重分配方法都不是最佳過渡方式，不能很好地遏制波浪反射.雖然指數方法將波幅削減至0所需最短空間長度略大于后兩種方法，但能更充分遏制波浪反射，得到最為精確的波浪模擬結果.綜合考慮，在使用松弛區方法進行波浪消波處理時，應當優先選擇指數法方法作為權重分配方法.

圖11 不同消波權重下波高空間分布Fig.11 Spatial distributions of wave elevations at different weights of wave absorption

3 結論

本次研究基于OpenFOAM，對Stokes五階波進行了模擬分析，著重研究了幾何VOF方法及松弛區方法在波浪模擬中的應用.

在網格和庫朗數收斂性分析中發現，若要模擬出較為滿意的波浪，在一個波高范圍內至少垂向設置8層網格，庫朗數可取Co=1.在不同波陡的波浪模擬中，引入了代數VOF方法，將之與幾何VOF方法的模擬結果和理論值進行比較發現：兩種VOF方法模擬的波浪均會產生較小的空間相位差，而幾何VOF能夠更好地對波高進行模擬.在不同波頻模擬中，幾何VOF方法能較好地對不同波頻下波高及空間相位進行模擬；代數VOF能較好地對相位進行模擬，而波高結果則略大于理論值.通過對比幾何和代數VOF方法下自由液面處液體體積分數分布情況，可知幾何VOF之所以較代數VOF方法有更高的波浪模擬精度，是因為幾何VOF方法通過網格重構有效減少了偽擴散現象.

在松弛區方法消波的研究中，對消波區長度及權重方法進行了分析.對于消波區長度，較短消波區長度雖足夠以波幅削減至0，但不能很好地遏制波浪反射.因此，需要選取適當的消波區長度(2倍波長及以上)，以確保能遏制波浪反射，從而避免對入射波的模擬產生影響.在權重分布的選擇上，相比于三階多項式權重及自由多項式權重，采用指數方法進行權重分配的消波效果最好.