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

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

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

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

為了改善質量源造波的效果，上述文獻均對質量源形狀進行了限制.質量源形狀影響造波的主要原因是推導源強方程時沒有完整地考慮質量輸出在各個方向對造波所做的貢獻，從而導致在一些情況下模擬值較理論解偏小.尤其在長寬比不滿足條件時，僅考慮單方向質量輸出不能滿足精度要求，因此才有了使源項強度線性增大來彌補波能損失這一措施，但對于不同的波況，源強增大倍數(shù)也各有不同，每次都需要試算來確定顯然費時費力.

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

1 數(shù)值模型

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

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

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

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

1.2 改進的質量源造波理論推導

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

式中：u和v分別是質量源區(qū)內質點的水平和垂向速度；x和y為位置，s（x，y，t）為質量源源強.不同的目標波浪對應不同的源強表達式，其推導過程如下：

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

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

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

在質量源區(qū)，假設每個網(wǎng)格應該滿足以下關系：

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

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

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

最終，源強的表達式為

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

1.3 消波方法

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

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

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

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

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

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

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

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

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

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

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

表2 時間步長的驗證Tab.2 Verification of time step effects

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

數(shù)值造波對非線性參數(shù)ε要求相對比較苛刻，對于非線性較強的波浪，往往模擬效果不理想[19].為了檢驗改進后質量源造波方法的可行性和適用范圍，本文探究了相同水深下不同波高時Stokes二階波的模擬效果.設置時間步長為0.02s，網(wǎng)格0.004 m，d＝0.2m，T＝1s，H＝0.01，0.04，0.05m，相應ε＝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時，本文方法模擬效果比較理想.雖然波谷處與理論波面存在微小誤差，但仍在精度范圍內，波峰與理論波面吻合很好，整個波形的模擬效果非常理想.當H／d＝0.25時，波谷處數(shù)值解與理論波形的偏差增大，且整個波形有上漂趨勢.因此，本文方法對于ε＝0.05～0.2的Stokes二階波模擬效果較好.

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

ε為衡量數(shù)值水槽造波能力的重要參數(shù)，為了驗證本文方法對造波效果的改善作用，以ε大小為標準與前人模擬結果進行了比較.文獻[20]基于質量源造波方法對Stokes二階波進行了數(shù)值模擬，波浪參數(shù)中H＝0.1m，d＝1.8m，文獻[6]利用質量源造波方法模擬Stokes二階波，在源強推導中考慮了水平方向質量輸出，選取H／d為0.2和0.3.針對3種不同ε值，本文的數(shù)值模擬結果與文獻進行對比，結果如圖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可見：本文改進后的方法在ε＝0.05～0.20時，對波幅誤差有較好的改善效果；尤其在ε＝0.20時，波幅誤差及相位誤差相對于傳統(tǒng)質量源造波方法分別降低21.7%和46.5%；當ε＝0.30時，波幅誤差與文獻中模擬方法相比結果接近，但相位誤差減小了63.6%.

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

同前人工作類似，文獻[7]通過設置多組波況擬合得到源區(qū)寬度的經(jīng)驗公式，不同的相對水深（kd）代入經(jīng)驗公式可得到不同的源區(qū)寬度，而本文在介紹改進的質量源源項推導過程時，發(fā)現(xiàn)改進方法的造波效果不依賴于源區(qū)形狀.為了探究這一結論的正確性，本文設置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可知，改進的質量源造波方法在相對水深為0.30～1.08均適用，無需嚴格按照經(jīng)驗公式[7]計算源區(qū)寬度，這證實本文方法不僅降低了造波時對源區(qū)形狀的限制，并且能取得較好的模擬效果.

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

為了驗證本文方法的實用性，利用改進后的質量源方法對孤立波的傳播進行模擬.孤立波的自由表面方程如下：

將其代入式（12）即可得到孤立波對應的源強表達式.

3.1 單個孤立波的傳播模擬

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

左右兩列孤立波W1和W2為鏡像對稱，因此選取任意1列作為研究對象即可.圖10展示了不同時刻下W2的波面圖，H＝0.005m，將模擬結果無因次化，可見孤立波在傳播過程中波形穩(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)了擾動，導致一部分能量損失，并且隨著ε值的增大，擾動出現(xiàn)了增大趨勢.分析發(fā)現(xiàn)，ε達到0.2時，模擬精度在可接受范圍內.

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

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

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

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

圖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

文獻[21]通過理論推導，得出孤立波疊加的理論最大波高計算公式為

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

4 結語

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

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

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

（3）本文方法提高了對強非線性波浪的模擬精度，對于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ù)值解與理論解吻合很好，并且同樣適用于孤立波碰撞疊加過程的模擬.

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