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

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

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

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

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

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

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

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

1 數(shù)學模型

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

1.1 控制方程

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

(1)

(2)

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

(3)

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

1.2 VOF方法

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

ρ=α1ρ1+α2ρ2

(4)

μ=α1μ1+α2μ2

(5)

(6)

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

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

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

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

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

(7)

(8)

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

(9)

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

1.3 松弛區(qū)方法

松弛區(qū)方法的原理為：引入權(quán)重因子w，對松弛區(qū)內(nèi)流體速度和體積分數(shù)進行加權(quán)，實現(xiàn)目標解析解和CFD計算解的過渡，從而實現(xiàn)速度入口造波邊界和CFD計算域間的過渡，以及CFD計算域與壓力出口邊界的過渡，對應(yīng)的表達式如下：

(10)

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

(11)

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

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

w=1-w*Co/Comax

(12)

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

(13)

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

(14)

w*=1-γp

(15)

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

2 結(jié)果與討論

2.1 計算模型

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

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

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

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

2.2 收斂性分析

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

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

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

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

2.2.2庫朗數(shù)收斂性分析 庫朗數(shù)定義如下：

(16)

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

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

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

2.3 造波分析

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

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

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

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

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

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

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

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

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

2.4 消波分析

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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