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

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

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

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

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

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

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

1 數(shù)學(xué)模型

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

1.1 控制方程

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

(1)

(2)

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

(3)

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

1.2 VOF方法

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

ρ=α1ρ1+α2ρ2

(4)

μ=α1μ1+α2μ2

(5)

(6)

式中：ρ1,ρ2分別為液體和氣體的密度；μ1,μ2分別為液體和氣體的動(dòng)力黏度；α1,α2分別為液體和氣體的體積分?jǐn)?shù).

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

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

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

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

(7)

(8)

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

(9)

式中：等號(hào)左端最后一項(xiàng)為人工壓縮項(xiàng)；uc,j的處理詳見文獻(xiàn)[15].

1.3 松弛區(qū)方法

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

(10)

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

(11)

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

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

w=1-w*Co/Comax

(12)

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

(13)

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

(14)

w*=1-γp

(15)

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

2 結(jié)果與討論

2.1 計(jì)算模型

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

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

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

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

2.2 收斂性分析

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

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

2.2.1網(wǎng)格收斂性分析 為了能夠較好地模擬波浪傳播，需要在自由表面處進(jìn)行網(wǎng)格加密.一般來說,在橫向、縱向和垂向3個(gè)方向中，垂向網(wǎng)格密度對(duì)波浪模擬的影響尤為明顯.因此，在自由液面處，固定垂向網(wǎng)格尺寸與縱向(波浪傳播方向)網(wǎng)格尺寸比為1∶2，在垂向選取4組網(wǎng)格H/Δz=6,8,12,16(H為波高，Δz為垂向網(wǎng)格尺寸).圖4所示為模擬時(shí)間為30 s時(shí)不同網(wǎng)格密度下波高的空間分布，圖中：xi為沿計(jì)算域縱向方向的距離；η為波高.由圖可知，當(dāng)網(wǎng)格密度H/Δz>8時(shí)，幾何VOF方法能較好地對(duì)波浪進(jìn)行模擬；而當(dāng)H/Δz<8時(shí)，幾何VOF方法模擬波浪在波高和相位上隨著波浪傳播逐漸與理論解產(chǎn)生偏離.因此，在使用幾何VOF方法對(duì)波浪進(jìn)行模擬時(shí)，單個(gè)波高垂向分布網(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為時(shí)間步長；|U|為某個(gè)網(wǎng)格單元內(nèi)速度矢量的模，Δx為速度方向網(wǎng)格長度.在計(jì)算域網(wǎng)格已劃定的情況下，網(wǎng)格長度Δx確定；在同一工況下，固定位置的速度矢量的模|U|也確定.因此，可以用庫朗數(shù)Co對(duì)CFD計(jì)算的時(shí)間步長Δt進(jìn)行反映.用于收斂性分析所選取的庫朗數(shù)Co=0.1,0.2,0.4,1.0，圖5為模擬時(shí)間為30 s時(shí)不同庫朗數(shù)下波高空間分布.

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

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

2.3 造波分析

圖6為模擬時(shí)間為30 s時(shí)不同波陡下幾何和代數(shù)VOF方法波高模擬空間分布，其中波高為H、波長為λ、波陡為δ.圖7為模擬時(shí)間為30 s時(shí)圖6對(duì)應(yīng)的波高模擬時(shí)歷曲線.從圖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方法波高模擬時(shí)歷曲線Fig.7 Time histories of wave elevations by using geometrical and algebraic VOF method at varied wave steepnesses

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

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

圖8為模擬時(shí)間為30 s時(shí)不同波頻下幾何和代數(shù)VOF方法波高模擬空間分布，其中波頻用ω表示.比較結(jié)果顯示，對(duì)于不同的波浪頻率，幾何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方法在求解體積分?jǐn)?shù)之前對(duì)網(wǎng)格進(jìn)行了重構(gòu)，使其能夠更加精確地自由液面進(jìn)行捕捉.圖9為幾何VOF方法和代數(shù)VOF方法分別對(duì)應(yīng)的自由液面處液體體積分?jǐn)?shù)分布情況，藍(lán)色部分代表α1=0 (氣體)，紅色部分代表α1=1 (液體)；黑線α1=0.01，黃線α1=0.99.用液體體積分?jǐn)?shù)為0.01和0.99兩個(gè)分界之間的厚度以表示自由液面捕捉精度，可以看到幾何VOF方法所得的自由液面處液體分布更為集中，即界面厚度比代數(shù)VOF方法更小.這體現(xiàn)了幾何VOF方法較代數(shù)VOF方法有更高的自由液面捕捉精度，即能更好地對(duì)波浪進(jìn)行模擬.

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

2.4 消波分析

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

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

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

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

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

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

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

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

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