基于改進(jìn)質(zhì)量源造波方法的非線性波數(shù)值模擬

馬 哲， 周 婷， 孫家文， 房克照， 翟鋼軍

（1.大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連116024；2.大連理工大學(xué) 深海工程研究中心，遼寧 大連116024；3.國家海洋環(huán)境監(jiān)測中心，遼寧 大連116023）

基于計算流體力學(xué)（CFD）的數(shù)值造波方法是海洋工程領(lǐng)域的研究熱點(diǎn)，常用的有推板造波法、速度入口法和源造波法.文獻(xiàn)[1]基于速度入口法實(shí)現(xiàn)了二階Stokes的數(shù)值模擬.文獻(xiàn)[2－3]針對推板造波方法開展了相關(guān)研究，成功模擬了線性波和瞬時極限波.源造波法根據(jù)源項(xiàng)添加位置的不同可分為質(zhì)量源和動量源兩種[4－5]，作為一種新型造波方法，因能充分利用數(shù)值模擬的優(yōu)勢而受到越來越多的關(guān)注.以往研究[5]表明，質(zhì)量源造波較動量源法適用的水深范圍更廣，且可以實(shí)現(xiàn)線性波和非線性波等不同目標(biāo)波浪的模擬，但源區(qū)形狀會對目標(biāo)波浪的模擬產(chǎn)生明顯影響.文獻(xiàn)[6]通過試算得出質(zhì)量源造波源區(qū)的長度和寬度只有滿足特定條件時才能達(dá)到比較理想的模擬效果；文獻(xiàn)[7]在內(nèi)域造波中通過多組波況的模擬分析擬合得到質(zhì)量源長度和寬度的經(jīng)驗(yàn)公式；文獻(xiàn)[8]采用水平形式的質(zhì)量源，考慮垂向質(zhì)量輸出，模擬了孤立波的行進(jìn)過程；同樣對孤立波傳播進(jìn)行數(shù)值模擬，文獻(xiàn)[9]選擇了垂向造波源，并指出以一個網(wǎng)格的水平長度Δx作為源區(qū)寬度可保證波浪模擬的精確性.然而，采用質(zhì)量源進(jìn)行數(shù)值造波得到的波高幅值往往小于理論值，僅僅依靠調(diào)節(jié)質(zhì)量源的長度和寬度對模擬效果的改善存在局限性.文獻(xiàn)[10]認(rèn)為質(zhì)量源在造波過程中的部分能量未能發(fā)揮作用，可以通過在源項(xiàng)表達(dá)式中乘一個放大系數(shù)以彌補(bǔ)損失的能量；文獻(xiàn)[11]針對淺水情況，在源強(qiáng)方程中考慮了垂向源區(qū)上方的質(zhì)量輸出，避免了造波過程中源區(qū)上方水體高度值多次試算帶來的不便.

為了改善質(zhì)量源造波的效果，上述文獻(xiàn)均對質(zhì)量源形狀進(jìn)行了限制.質(zhì)量源形狀影響造波的主要原因是推導(dǎo)源強(qiáng)方程時沒有完整地考慮質(zhì)量輸出在各個方向?qū)υ觳ㄋ龅呢暙I(xiàn)，從而導(dǎo)致在一些情況下模擬值較理論解偏小.尤其在長寬比不滿足條件時，僅考慮單方向質(zhì)量輸出不能滿足精度要求，因此才有了使源項(xiàng)強(qiáng)度線性增大來彌補(bǔ)波能損失這一措施，但對于不同的波況，源強(qiáng)增大倍數(shù)也各有不同，每次都需要試算來確定顯然費(fèi)時費(fèi)力.

基于上述問題，本文提出了一種無需限制源區(qū)形狀的改進(jìn)方法.首先從理論分析的角度，揭示了傳統(tǒng)質(zhì)量源造波法非線性波模擬出現(xiàn)波幅偏小的原因，在此基礎(chǔ)上修改源項(xiàng)方程，針對矩形質(zhì)量源，充分考慮源區(qū)質(zhì)量輸出在各個方向?qū)υ觳ㄋ龅呢暙I(xiàn)，得到了新的質(zhì)量輸出的源項(xiàng)函數(shù).然后，利用改進(jìn)方法模擬了不同相對水深、不同非線性條件下的Stokes二階波和孤立波.在本文的模擬過程中，改進(jìn)方法不依賴于源區(qū)形狀的長寬比，因此對于不同的波況無需多次進(jìn)行源區(qū)形狀的調(diào)試工作，大大提高了工作效率.

1 數(shù)值模型

1.1 控制方程與數(shù)值求解方法

數(shù)值模型基于開源CFD類庫OpenFOAM中不可壓縮氣液二相流求解器interFoam[12]，采用Volume of Fluid（VOF）方法[13－14]進(jìn)行自由面的捕捉，控制方程為連續(xù)性方程、引入相體積分?jǐn)?shù)α的動量方程以及VOF方程：

式中：u為平均速度；t為時間；ρτ為雷諾應(yīng)力項(xiàng)，其中τ表征由于Δ脈動值引起的雷諾應(yīng)力張量，為二階對稱張量；Cκα為表面張力項(xiàng)，C為表面張力系數(shù)，一般取0.07kg／s2，κ為自由面曲率；g為重力加速度；X為位置矢量；ur代表相對速度，含有ur的項(xiàng)為人工壓縮項(xiàng)[15]，只在二相交界面處起作用.ρ和μ分別為二相流的平均密度和平均動力黏度：

對控制方程的求解需要對時間和空間進(jìn)行離散，OpenFOAM采用有限體積法離散框架，時間積分項(xiàng)采用隱式歐拉格式，體積分中的對流項(xiàng)采用Guass limitedLinearV 1格式離散，黏性擴(kuò)散項(xiàng)采用linear corrected，其余用線性插值.速度壓力場采用分離式（segregated）算法PISO.離散后的線性方程組采用預(yù)處理共軛梯度法求解.

1.2 改進(jìn)的質(zhì)量源造波理論推導(dǎo)

為了實(shí)現(xiàn)質(zhì)量源造波，需要在造波源區(qū)對式（1）添加源項(xiàng)，如下所示.

式中：u和v分別是質(zhì)量源區(qū)內(nèi)質(zhì)點(diǎn)的水平和垂向速度；x和y為位置，s（x，y，t）為質(zhì)量源源強(qiáng).不同的目標(biāo)波浪對應(yīng)不同的源強(qiáng)表達(dá)式，其推導(dǎo)過程如下：

式中：c為目標(biāo)波浪的相速度；η（t）為自由表面方程；Ω為質(zhì)量源所在的區(qū)域.式（7）的意義為在某一時間段內(nèi)，質(zhì)量源區(qū)質(zhì)量的輸出用于產(chǎn)生目標(biāo)波浪[16].右端系數(shù)2表示質(zhì)量源區(qū)的質(zhì)量輸入向左右兩邊都產(chǎn)生了波浪.簡化式（7）可得

式中：A為質(zhì)量源區(qū)的面積.

上述推導(dǎo)過程考慮了質(zhì)量源在相速度方向的質(zhì)量輸出，忽略了垂向，而實(shí)際上源區(qū)的質(zhì)量輸出在4個方向同時存在.基于此，本文進(jìn)行了如下改進(jìn).

在質(zhì)量源區(qū)，假設(shè)每個網(wǎng)格應(yīng)該滿足以下關(guān)系：

式中：Ls和Hs分別是質(zhì)量源的長和寬.對式（9）進(jìn)行化簡可得

根據(jù)文獻(xiàn)[17]，質(zhì)量源造波在原理上類似于水下爆炸，爆炸瞬間質(zhì)點(diǎn)向各個方向的運(yùn)動速度相同，質(zhì)量源向外輸出流體沿各方向的速度也相等，即u＝v.最終源強(qiáng)簡化為

為了保證數(shù)值的穩(wěn)定性，源強(qiáng)函數(shù)的初始值應(yīng)為0.在源強(qiáng)中添加緩啟動因子（SSF），定義在很短的初始時間段t0內(nèi)，源強(qiáng)函數(shù)的值從0開始增加，具體表達(dá)如下：

最終，源強(qiáng)的表達(dá)式為

由式（11）可見：當(dāng)Ls和Hs不在同等量級時，以橫向質(zhì)量源（Ls?Hs）為例，此時可以忽略式（10）右端第一項(xiàng)，僅需要考慮質(zhì)量源區(qū)垂向方向質(zhì)量輸出就能滿足造波的精度，所以利用式（8）推導(dǎo)出的質(zhì)量源源項(xiàng)在橫向質(zhì)量源中可以取得較好模擬結(jié)果；同理，豎向質(zhì)量源（Ls?Hs）只需考慮橫向質(zhì)量輸出也能達(dá)到數(shù)值造波的精度要求.但當(dāng)質(zhì)量源的長寬比介于0.1～10時，只考慮單方向質(zhì)量輸出時得到的數(shù)值結(jié)果普遍較差，因?yàn)檫@種近似本身存在理論上的缺陷，需要對質(zhì)量源形狀限制來彌補(bǔ).式（11）考慮了質(zhì)量源區(qū)4個方向的質(zhì)量輸出，其造波精度不依賴于矩形源區(qū)的形狀變化，提高了源區(qū)形狀變化的靈活性.

1.3 消波方法

為了防止波浪反射對模擬效果的影響，本文在數(shù)值水槽的左右兩端分別添加阻尼層.常用的阻尼消波方式有指數(shù)、線性、根號、平方等形式，合理選擇消波函數(shù)會對反射率產(chǎn)生影響[18].本文算例均采用線性阻尼消波形式，衰減系數(shù)取8，可取得較好的消波效果.

2 Stokes二階波的數(shù)值模擬

以二維情形為例進(jìn)行模擬，數(shù)值波浪水槽總長8.46m，高0.26m，水深0.2m，水槽左右兩端和底端為固壁邊界類型，上端為壓力出口邊界.質(zhì)量源區(qū)位于水槽中間，根據(jù)文獻(xiàn)[6]源區(qū)位置設(shè)置經(jīng)驗(yàn)法則，源區(qū)中心距靜水面的距離設(shè)為1／3水深，阻尼消波段首尾各為1.21m，如圖1所示.

圖1 二維數(shù)值水槽示意圖Fig.1 Diagram of 2Dnumerical wave tank

Stokes二階波的自由表面函數(shù)為

式中：k為波數(shù)；d為水深；ω為圓頻率.將式（14）代入式（12），得到Stokes二階波的質(zhì)量源源強(qiáng)表達(dá)式.

2.1 網(wǎng)格、時間步長對造波效果的影響

既往研究[2]表明，網(wǎng)格和時間步長尺度選取帶來的數(shù)值黏性對數(shù)值水槽的模擬效果具有重要影響.以Stokes二階波為例，水深d＝0.2m，波高H＝0.04m，并定義無量綱參數(shù)ε表示波浪的非線性（ε＝H／d）.分別選取3組不同網(wǎng)格大小和4組不同時間步長的工況，對數(shù)值模擬結(jié)果的影響進(jìn)行比較分析.

3組不同網(wǎng)格的設(shè)置和模擬結(jié)果如表1所示（其中T為波浪周期），橫向網(wǎng)格始終保持為0.012 m，可見加密垂向網(wǎng)格有利于數(shù)值模擬精度的提高.從算例M1到算例M2，數(shù)值模擬波幅與理論值的誤差明顯減小.從算例M2到算例M3，網(wǎng)格數(shù)量翻倍，但誤差結(jié)果相差不大.模擬結(jié)果證明1個波高內(nèi)有10個網(wǎng)格時可實(shí)現(xiàn)理想的造波精度，因此采用0.004m作為垂向網(wǎng)格的長度.

表1 網(wǎng)格的驗(yàn)證Tab.1 Verification of mesh

4組不同時間步長的設(shè)置和模擬結(jié)果如表2所示，時間步長對數(shù)值模擬的影響程度不大，當(dāng)時間步長取0.02s時，模擬結(jié)果與理論值的誤差為3.15%，滿足精度要求.算例T1和T2時間步長分別取0.002和0.001s，同算例T3相比，誤差改善效果不明顯，但是計算時間大大增加.綜合考慮后，選取時間步長為0.02s，即為0.02T時，計算代價最小.

表2 時間步長的驗(yàn)證Tab.2 Verification of time step effects

2.2 不同波高時二階Stokes波的模擬

數(shù)值造波對非線性參數(shù)ε要求相對比較苛刻，對于非線性較強(qiáng)的波浪，往往模擬效果不理想[19].為了檢驗(yàn)改進(jìn)后質(zhì)量源造波方法的可行性和適用范圍，本文探究了相同水深下不同波高時Stokes二階波的模擬效果.設(shè)置時間步長為0.02s，網(wǎng)格0.004 m，d＝0.2m，T＝1s，H＝0.01，0.04，0.05m，相應(yīng)ε＝0.05，0.2，0.25.距離源區(qū)中心1倍波長處的自由液面高度隨時間的變化過程以及與理論值對比如圖2～4所示.

圖2 ε＝0.05時數(shù)值解和理論解的對比 （H＝0.01m）Fig.2 Comparison of numerical results and theoretical analysis（H＝0.01m，ε＝0.05）

圖3 ε＝0.2時數(shù)值解和理論解的對比 （H＝0.04m）Fig.3 Comparison of numerical results and theoretical analysis（H＝0.04m，ε＝0.2）

圖4 ε＝0.25時數(shù)值解和理論解的對比 （H＝0.05m）Fig.4 Comparison of numerical results and theoretical analysis（H＝0.05m，ε＝0.25）

由圖可見，ε＝0.05，0.2時，本文方法模擬效果比較理想.雖然波谷處與理論波面存在微小誤差，但仍在精度范圍內(nèi)，波峰與理論波面吻合很好，整個波形的模擬效果非常理想.當(dāng)H／d＝0.25時，波谷處數(shù)值解與理論波形的偏差增大，且整個波形有上漂趨勢.因此，本文方法對于ε＝0.05～0.2的Stokes二階波模擬效果較好.

2.3 數(shù)值模擬精度的改善

ε為衡量數(shù)值水槽造波能力的重要參數(shù)，為了驗(yàn)證本文方法對造波效果的改善作用，以ε大小為標(biāo)準(zhǔn)與前人模擬結(jié)果進(jìn)行了比較.文獻(xiàn)[20]基于質(zhì)量源造波方法對Stokes二階波進(jìn)行了數(shù)值模擬，波浪參數(shù)中H＝0.1m，d＝1.8m，文獻(xiàn)[6]利用質(zhì)量源造波方法模擬Stokes二階波，在源強(qiáng)推導(dǎo)中考慮了水平方向質(zhì)量輸出，選取H／d為0.2和0.3.針對3種不同ε值，本文的數(shù)值模擬結(jié)果與文獻(xiàn)進(jìn)行對比，結(jié)果如圖5～7和表3所示.

圖5 二階Stokes波時歷曲線 （ε＝0.05）Fig.5 Time history of second－order Stokes（ε＝0.05）

圖6 二階Stokes波時歷曲線 （ε＝0.2）Fig.6 Time history of second－order Stokes（ε＝0.2）

圖7 二階Stokes波時歷曲線 （ε＝0.3）Fig.7 Time history of second－order Stokes（ε＝0.3）

表3 數(shù)值誤差對比Table 3 Comparison of numerical errors

由表3可見：本文改進(jìn)后的方法在ε＝0.05～0.20時，對波幅誤差有較好的改善效果；尤其在ε＝0.20時，波幅誤差及相位誤差相對于傳統(tǒng)質(zhì)量源造波方法分別降低21.7%和46.5%；當(dāng)ε＝0.30時，波幅誤差與文獻(xiàn)中模擬方法相比結(jié)果接近，但相位誤差減小了63.6%.

2.4 不同相對水深時二階Stokes波的模擬

同前人工作類似，文獻(xiàn)[7]通過設(shè)置多組波況擬合得到源區(qū)寬度的經(jīng)驗(yàn)公式，不同的相對水深（kd）代入經(jīng)驗(yàn)公式可得到不同的源區(qū)寬度，而本文在介紹改進(jìn)的質(zhì)量源源項(xiàng)推導(dǎo)過程時，發(fā)現(xiàn)改進(jìn)方法的造波效果不依賴于源區(qū)形狀.為了探究這一結(jié)論的正確性，本文設(shè)置4組波況，選取H＝0.01m，d＝0.2m （ε＝0.05），T＝1.0，1.5，2.0，3.0s，則距離源區(qū)中心1倍波長處的自由液面高度隨時間的變化過程與理論值對比如圖8所示.

圖8 不同周期時數(shù)值解和理論解的對比Fig.8 Comparison of numerical results and theoretical analysis with different periods

由圖8可知，改進(jìn)的質(zhì)量源造波方法在相對水深為0.30～1.08均適用，無需嚴(yán)格按照經(jīng)驗(yàn)公式[7]計算源區(qū)寬度，這證實(shí)本文方法不僅降低了造波時對源區(qū)形狀的限制，并且能取得較好的模擬效果.

3 孤立波傳播的數(shù)值模擬

為了驗(yàn)證本文方法的實(shí)用性，利用改進(jìn)后的質(zhì)量源方法對孤立波的傳播進(jìn)行模擬.孤立波的自由表面方程如下：

將其代入式（12）即可得到孤立波對應(yīng)的源強(qiáng)表達(dá)式.

3.1 單個孤立波的傳播模擬

孤立波為一種淺水波，有其獨(dú)特性.數(shù)值水槽尺寸如圖9所示.時間步長為0.001s，水平網(wǎng)格邊長為0.02m，垂向網(wǎng)格邊長為0.01m.

左右兩列孤立波W1和W2為鏡像對稱，因此選取任意1列作為研究對象即可.圖10展示了不同時刻下W2的波面圖，H＝0.005m，將模擬結(jié)果無因次化，可見孤立波在傳播過程中波形穩(wěn)定且振幅衰減小.

圖9 單源數(shù)值水槽示意圖Fig.9 Diagram of a single－source numerical wave tank

圖10 不同時刻下孤立波的波面Fig.10 Wave surface elevation of solitary waves

圖11分別展示了在d＝0.2m，t＝5s時H＝0.005，0.01，0.04，0.05m的孤立波波面曲線.由圖11可見，3種不同波高情況下基于本文方法得到的波面與理論解吻合很好，隨著ε的增大，波峰處誤差逐漸增大，主要是由于兩列孤立波中間出現(xiàn)了擾動，導(dǎo)致一部分能量損失，并且隨著ε值的增大，擾動出現(xiàn)了增大趨勢.分析發(fā)現(xiàn)，ε達(dá)到0.2時，模擬精度在可接受范圍內(nèi).

圖11 不同波高時孤立波的數(shù)值解與理論解的對比 （t＝5s）Fig.11 Comparison between numerical results and theoretical analysis with different wave heights（t＝5s）

3.2 孤立波“碰撞”過程的數(shù)值模擬

為了檢驗(yàn)本文方法的有效性，選取兩列孤立波相向運(yùn)動及“碰撞”后演化過程作為數(shù)值算例.為了產(chǎn)生兩列相向運(yùn)動的孤立波，在水槽中設(shè)置2個質(zhì)量源區(qū)，數(shù)值水槽尺寸如圖12所示.

左右設(shè)置的2個質(zhì)量源區(qū)S1和S2都產(chǎn)生了沿左右傳播的孤立波，左側(cè)S1產(chǎn)生的波標(biāo)記為L1，L2，右側(cè)S2產(chǎn)生的波標(biāo)記為R1，R2.圖13給出了孤立波產(chǎn)生和傳播的過程.當(dāng)t＝5s時，S1和S2分別產(chǎn)生兩列孤立波，在t＝6s時，L2和R1剛好重疊，此時能量達(dá)到最大.

圖12 水槽布置示意圖Fig.12 Layout of the numerical wave tank

圖13 孤立波碰撞過程中不同時刻下的波面圖Fig.13 Wave surface elevation at different time during the process of collisions in solitary waves

文獻(xiàn)[21]通過理論推導(dǎo)，得出孤立波疊加的理論最大波高計算公式為

式中：右側(cè)第3項(xiàng)和第4項(xiàng)分別為2階和3階小量，在此忽略其量綱的貢獻(xiàn).本文取H1和H2都為0.01m，計算得到理論最大波高為0.020 05m，數(shù)值解和理論解的誤差小于1%.因此本文方法同樣適用于孤立波“碰撞”過程模擬.

4 結(jié)語

本文針對傳統(tǒng)質(zhì)量源造波法對強(qiáng)非線性波模擬效果不佳的問題，考慮雙向質(zhì)量輸出，提出改進(jìn)方法.通過理論分析與推導(dǎo)得出新的質(zhì)量源源項(xiàng)函數(shù)，利用該源項(xiàng)公式對Stokes二階波、孤立波傳播過程進(jìn)行數(shù)值模擬，計算結(jié)果表明：

（1）質(zhì)量源對于垂向網(wǎng)格尺寸比較敏感，1個波高內(nèi)含10個網(wǎng)格時模擬效果較好；時間步長越短，數(shù)值黏性越小，取時間步長為0.02T時，在保證造波精度的前提下減少了計算代價.

（2）本文方法避免了對質(zhì)量源區(qū)形狀的過多限制（傳統(tǒng)線形造波源長寬比小于0.1或大于10），省去了源區(qū)形狀調(diào)試工作，提高了工作效率.

（3）本文方法提高了對強(qiáng)非線性波浪的模擬精度，對于H／d＝0.05～0.3的Stokes二階波取得很好的模擬效果；尤其在H／d＝0.2時，波幅誤差及相位誤差相對于傳統(tǒng)方法分別降低21.7%和46.5%.

（4）對于H／d＝0.025～0.25的孤立波，本文方法的數(shù)值解與理論解吻合很好，并且同樣適用于孤立波碰撞疊加過程的模擬.

圍繞質(zhì)量源造波的相關(guān)問題還需要進(jìn)一步開展相關(guān)研究，如改善質(zhì)量源造波在強(qiáng)非線性范圍的造波精度；參照前人的工作，將質(zhì)量源造波方法與動量源造波方法進(jìn)行比較和分析；在現(xiàn)有二維算例的基礎(chǔ)上，拓展三維質(zhì)量源造波方法，并開展斜向波的數(shù)值模擬.