劉必勁，張振偉，劉忠波

（1.浙江大學海南研究院，海南三亞 572025；2.廈門理工學院土木工程與建筑學院，福建廈門 361024；3.大連海事大學交通運輸工程學院，遼寧大連 116024）

0 引 言

畸形波是非線性極強的瞬態(tài)大浪，受實際海洋中各種不確定的眾多環(huán)境因素影響，導致人類對這種畸形波的產(chǎn)生機制并沒有完全掌握。當這種瞬態(tài)大浪作用在過往的船只或海上建筑物時，將會嚴重影響它們的安全。為了探求畸形波的生成機制，國內外學者嘗試從相位聚焦、波-流非線性相互作用、地形聚焦作用、表面風的作用、波浪調制不穩(wěn)定性等方面開展畸形波的研究。

對畸形波的研究主要方法有物理模型實驗室模擬[1-5]、數(shù)值模擬[6]和理論分析[7-8]等。近年來，伴隨計算機性能的大幅度提升和各類數(shù)值模型的不斷完善，采用數(shù)值模型開展畸形波的研究成為比較常見且主要的手段，其中時域計算模型多采用勢流理論、歐拉方程或粘性流理論等。在勢流理論方面，有高階邊界元求解Laplace 方程的數(shù)值模型[6]、高階譜數(shù)值模型[9-11]和Boussinesq 型水波數(shù)值模型[12-13]。歐拉方程或粘性Navier-Stokes方程基礎上建立起來的非靜壓模型[14]，將內部造波法引入CIP求解對流項的Navier-Stokes 方程的數(shù)值模型[15]。相關研究表明，采用相位聚焦方式可以形成畸形波，由于畸形波的強非線性特征，線性理論和二階Stokes波理論均無法精準預報聚焦波面和波峰面下的速度場[1]。

聚焦波演化形成畸形波的過程中，多種頻率組成成分的波浪非線性相互作用引起的多波共振或調制不穩(wěn)定等是生成畸形波的最有可能的原因之一。精確預報畸形波的生成與發(fā)展過程勢必要求數(shù)值模型具有足夠精準的色散性和非線性特征。一方面，聚焦波中波浪頻率組成成分越多，在同一水深下，這些波浪具有不同的無因次水深，要求模型在更大水深范圍內具備優(yōu)良的線性色散特性，才能保證每一種頻率波浪的計算結果與線性解析解具有良好的相位匹配關系[16]；另一方面，不同頻率的波浪之間存在較為明顯的非線性相互作用，這將產(chǎn)生低頻、高頻等其他頻率的差頻與和頻波浪，精準的非線性特征（二階和差頻、三階非線性波幅、波幅色散）則可以保證數(shù)值模型在模擬非線性波浪傳播時具有較好的適用性[16-19]。當方程在描述速度分布的精度不準確時，數(shù)值計算得到的速度場與相關實驗結果會出現(xiàn)一定的差異[13]。劉必勁等[13]采用Liu 和Fang[20]的單層Boussinesq 型水波數(shù)值模型，在開展聚焦波的數(shù)值模擬時發(fā)現(xiàn)，聚焦波面與試驗結果的吻合程度良好，但聚焦波峰值下的水平速度場計算值呈S 形，與Baldock 等[1]試驗的速度分布對比時僅基本吻合，這表明受限于單層Boussinesq 水波型方程在速度分布方面的精度，該數(shù)值模型無法勝任深水聚焦波速度場的模擬。為此，本文選用了Liu 和Fang[20]的雙層Boussinesq 型方程，建立相應的數(shù)值模型，研究該數(shù)值模型在應用于模擬深水聚焦波問題時，波面位移特別是聚焦波峰下水平速度的準確性。

1 雙層Boussinesq數(shù)值模型

1.1 控制方程

在無旋、無粘的假設下，Liu 和Fang[20]保留了垂向速度，經(jīng)推導給出了適合緩變地形的雙層Boussinesq 方程，立面二維模型由9 個方程組成，含有9 個獨立待求變量：波面位移η、自由表面處水平和垂向速度uη和wη、靜止水位處的速度u0和w0，以及計算速度（u*i，w*i）（i=1，2）。方程的具體表達形式如下：

在自由面處，滿足精確的運動學邊界條件和動力學邊界條件：

自由面處速度與靜止水位處速度之間的表達式為

靜止水位處速度與第一層計算速度的關系式為

在連接水深位置處，滿足速度連續(xù)，則有

在水底處，滿足運動學邊界條件：

式中，下標x表示對空間x求導，g是重力加速度。

方程（1）～（9）構成了本文數(shù)值模型對應的雙層Boussinesq水波方程，該方程具有較好的線性與非線性性能。對波浪沿水深分布的速度進行優(yōu)化，Liu等[18]給出了方程中的系數(shù)（α1，α2，β12，β13，β22，β23）=（0.129，0.371，0.941，0.896，0.948，0.608）。

1.2 數(shù)值模型的求解

雙層Boussinesq 水波方程的最高空間導數(shù)是3，其與目前常用于海岸波浪模擬的開源程序FUNWAVE[21]對應的Boussinesq 水波[22]方程最高空間導數(shù)的階數(shù)上是一致的。為此，該方程在時空差分格式的離散上均采用與該開源代碼一樣的數(shù)值格式。在時間差分格式上，采用混合4 階Adams-Bashforth-Moulton 格式。其中，預報時利用3 階Adams-Bashforth 格式求解方程（1）和（2），可以得到波面位移以及波面處水平速度的預報值；分別求解方程（3）、（5）和（7）可得到3個水平速度的預報值；進而求解方程（9）、（8）、（6）和（4），可依次得到4 個垂向速度的預報值。校正時則利用4 階Adams-Moulton 時間步進格式，得到波面位移以及波面處水平速度的校正值，其他速度變量校正值的求解過程類似于預報階段，不同點在于所有方程右端項均是當前預報層的值。當9 個變量的校正值與預報值之間的誤差控制在0.000 1 內，當前迭代結束，否則更新變量的預報值，其采用校正值和預報值的加權值，繼續(xù)進行當前步的迭代過程。校正值的加權系數(shù)R為0.05～0.1，而預報值加權系數(shù)為1-R。

2 數(shù)值模型的驗證

本章將重點研究聚焦波的一些特性，并與相關理論解析解和實驗結果進行比較，驗證本文數(shù)值模型的適用性。首先開展線性聚焦波的模擬，考察數(shù)值模型在具有多個不同周期組成下線性色散的適用性問題，其最重要的表現(xiàn)是疊加后波面的相位問題，當多波組成的波面與解析波面一致時，則表明數(shù)值模型在更大無因次水深范圍內有較好的色散精度；其次，開展非線性深水聚焦波的模擬，考察數(shù)值模型波-波非線性相互作用和非線性色散的精確性問題。

2.1 線性聚焦波的模擬

雙層Boussinesq 方程具有確切的色散關系表達式，其與Stokes 波色散關系表達式的對比在Liu 和Fang[20]以及Liu 等[18]中均明確給出，本文的雙層Boussinesq 數(shù)值模型中采用了Liu 等[18]給出的參數(shù)值，原圖可詳見Liu 等[18]中的圖3，這里為了展現(xiàn)其色散性能對圖進行了放大，重新繪制得到的無因次相速度見圖1。由圖1 可見，方程的相速度在無因次水深0＜kh＜20范圍內與解析解的誤差僅為0.003 68%，這表明雙層Boussinesq 水波方程具有極為精確的色散特征。盡管方程具有精確的色散關系式，但數(shù)值模型因離散誤差以及數(shù)值積累誤差等，一般來說，其計算出來的色散方面的性能低于理論分析得到的色散性能。

當忽略非線性項后，雙層Boussinesq 模型退化為線性模型，利用該模型對線性聚焦波進行模擬，目的是考察建立的雙層Boussinesq 數(shù)值模型的精準度，這里選擇Baldock 等[1]給出的波浪頻率進行數(shù)值計算。Baldock等[1]的聚焦波深水實驗含有不同波浪頻率段范圍，在不同周期段內等分周期成29份，每個頻率波幅取值相等，其中的B 組（T=0.6～1.4 s，寬譜）和D 組（T=0.8～1.2 s，窄譜）對應的線性無因次水深kh分別為1.568～7.825 和2.026～4.403。在數(shù)值模擬中，水深為0.7 m，計算域設置為24 m，末端設置4 m長的海綿邊界層吸收波浪，時間步長為0.01 s，空間步長為0.04 m。設置線性聚焦位置在8 m，線性聚焦時間為60 s，聚焦波幅為0.055 m，選擇這一聚焦時間點的原因是避開數(shù)值模型在模擬波群演化中波前（wave front）的影響。在聚焦時刻t=60 s，B 組和D 組的計算結果與一階Stokes 波解析解的比較見圖2 和圖3。由圖可見，數(shù)值模擬結果與解析解吻合程度高，這一方面表明數(shù)值模型可以精準捕捉不同頻率波浪的相位，另一方面也從數(shù)值角度佐證了如圖1展示的雙層Boussinesq水波方程在較大頻率范圍內精準的相速度。此外，在聚焦位置處，兩組工況的計算波峰面均是0.054 7 m，與解析波峰面0.055 m的誤差是0.5%。

圖1 無因次相速度Fig.1 Non-dimensional phase velocity

圖2 B55計算結果與Stokes波解析解的比較Fig.2 Comparison of the simulated surface elevation with the analytical solution for Case B55

圖3 D55計算結果與Stokes波解析解的比較Fig.3 Comparison of the simulated surface elevation with the analytical solution for Case D55

2.2 非線性聚焦波的模擬

在2.1 節(jié)線性聚焦波的模擬中，不同成分的線性波浪滿足疊加關系，以上僅是表征方程在線性波群中相位捕捉的能力。實際中，非線性波浪并不滿足線性疊加原理。一方面，聚焦波組成成分之間存在明確的非線性相互作用，這會導致波群演化過程中出現(xiàn)多個頻率之間產(chǎn)生和頻項與差頻項，因此，要求方程必須具有足夠精確的和頻與差頻性能；另一方面，根據(jù)三階Stokes波理論可知，波幅越大，波長將越長（同等條件下），在相位聚焦生成大波的過程中，波峰面越大，將會導致其與設置的線性聚焦位置發(fā)生更多偏離[1，14]。

在進行數(shù)值模擬之前，首先確認方程的非線性和速度分布特征方面的適用范疇。在1%誤差下，雙層Boussinesq 水波方程在二階非線性、三階非線性和波幅色散三個特性方面的最大適用水深kh是35.7、33.1和36.7，而線性水平速度和垂向速度的最大適用水深kh均是17.6，關于這些特性對應的詳圖可見參考文獻[18]。這表明雙層Boussinesq 數(shù)值模型可以勝任深水聚焦波演化的模擬。其次，對所采用的聚焦波物理試驗進行說明。第一，水槽中開展的波浪試驗，波高在沿程方向存在一定程度的衰減現(xiàn)象；第二，造波機驅動信號產(chǎn)生的波浪與理論解存在差異；第三，實際水槽無法克服邊界的反射問題，這種反射波會“污染”試驗結果。

針對Baldock 等[1]試驗中B 組和D 組的系列工況開展數(shù)值模擬研究，在計算中，數(shù)值水槽長24 m，末端設置4 m 的海綿邊界層消波?？臻g步長和時間步長分別是0.04 m 和0.01 s，計算時間為40 s。對B 組和D 組分別含有三組線性聚焦波峰面Af（分別是22 mm、38 mm 和55 mm）工況進行數(shù)值模擬。在入射邊界，造波采用線性波浪條件，設置的線性聚焦位置是8 m，通過反算可以確定入射波浪的情況，聚焦時間設置為30 s，其他條件均與Baldock 等[1]實驗描述的一致。實際計算得到的最大聚焦位置處，波面位移與實驗結果[1]的對比見圖4。相應地，最大聚焦波峰面下的水平速度與實驗結果的對比見圖5。從圖4中可見：非線性較強的工況B55和D55的計算結果比實驗結果大；作為非線性相對弱一些的D38 工況，計算波面比實驗結果小，其他非線性工況計算結果與實驗結果相差不大。從圖5 對比的水平速度場來分析，除D22 外，其他模擬結果與實驗結果的吻合程度較好，同時雙層模型計算的速度剖面比較光滑，不會出現(xiàn)單層模型給出的S 形[13]，這表明雙層Boussinesq 水波方程具有良好的速度分布特性。

圖4 計算波面位移與Baldock等[1]（1996）實驗結果的對比Fig.4 Comparison of the calculated surface elevation with the experimental results of Baldock et al[1](1996)

圖5 水平速度與Baldock等[1]（1996）實驗結果的對比Fig.5 Comparison of the calculated horizontal velocity with the experimental results of Baldock et al[1](1996)

此外，實驗中流速儀無法捕捉到波峰面處的速度，而數(shù)值模型可以提供波峰面速度，因此，速度場的最大值高于試驗值是可以接受的。所有工況算例表明，雙層Boussinesq數(shù)值模型可以較為準確地模擬深水聚焦波演化和垂向分布的速度場。

根據(jù)前面提及水槽實驗中的波高衰減等問題，對兩組非線性最強的兩個算例（B55和D55）減少入射波高進行模擬，其中，B55工況對應的線性聚焦波峰面Af設置為0.053 m,而D55工況對應的Af設置為0.053 5 m，具體計算得到的波面位移和水平速度結果見圖6。由該圖可見，計算的波面位移和水平速度均與實驗結果吻合。

圖6 計算的波面位移和水平速度與Baldock等（1996）實驗結果的對比Fig.6 Comparison of the calculated surface elevation and horizontal velocity with the experimental results of Baldock et al(1996)

3 數(shù)值模擬結果的進一步討論與分析

從已有的實驗研究[1]和相關數(shù)值模擬研究[6，13]可確知，保持中心頻率不變，在同樣線性聚焦峰面情況下，周期范圍越寬，則非線性越弱。因此本文以D組為例，首先與Baldock等[1]的實驗條件類似，保持線性聚焦位置與聚焦時間不變化，討論本文建立的雙層Boussinesq數(shù)值模型計算得到的最大波面位移ηf和最大波面位移處的水平速度uf隨線性聚焦波峰幅值Af的變化情況；其次，以Baldock 等[1]的D55 為例，討論線性聚焦位置（xf=6 m,7 m,8 m,9 m,10 m）對ηf和uf的影響，同時為避免短波沒有傳播到聚焦位置等對數(shù)值結果的影響，數(shù)值模型中將線性聚焦時間tf設置為35 s。最后，仍以Baldock 等[1]的D55為例，討論了計算域內每個空間點的最大波峰面位移η和對應η處的水平速度uη。

在數(shù)值模擬中，保持水深為0.7 m，計算域仍然設置為24 m，末端設置4 m 長的海綿吸收邊界，時間步長為0.01 s，空間步長為0.04 m，計算時間設置為40 s。

3.1 線性聚焦波峰幅值Af對ηf和uf的影響

圖7和圖8分別給出了D組系列工況中的最大波面位移ηf和最大波面位移處的水平速度uf隨線性聚焦波峰幅值Af（0.01～0.055 m，每個工況間隔0.005 m）的變化情況。由圖7 可見，數(shù)值模型計算的最大波峰面在不同工況時與實驗結果[1]趨勢基本一致且比較接近；進一步分析可知，伴隨線性聚焦波峰面（入射波高）的增大，最大波峰面呈現(xiàn)非線性變化。類似地，圖8給出的最大波峰面處的水平速度也具有同樣變化趨勢。

圖7 最大波面位移隨線性聚焦波峰幅值的變化Fig.7 Variation of the maximum surface elevation with the amplitude of the linearly focused wave crest

圖8 最大波峰面處水平速度隨線性聚焦波峰幅值的變化Fig.8 Variation of horizontal velocity at the maximum surface elevation with the amplitude of the linearly focused wave crest

3.2 線性聚焦位置xf對ηf和uf的影響

表1 給出了Af=55 mm 工況（tf=35 s）最大波面位移ηf和最大波面位移處的水平速度uf隨線性聚焦位置（xf=6 m,7 m,8 m,9 m,10 m）的變化情況。在表中也給出了相鄰兩次大波峰面之間的傳播速度cf，其計算公式可寫為

表1 ηf和uf隨xf的變化情況Tab.1 Variation of ηf and uf with xf

式中：xf1和xf2分別為相鄰兩次的最大波峰面位置，對應著時間tf1和tf2；xf1為設置線性聚焦位置后的第一個聚焦波峰面位置，xf2為相鄰的最大波峰面位置。需要強調的是，從數(shù)值計算結果來看，xf2可以在xf1之前，也可以在xf1之后，根據(jù)具體的數(shù)值模擬情況確定。

由表1 可知：一方面，伴隨線性聚焦位置的增大，線性聚焦位置后的第一個聚焦波峰面位置xf1處的波面和對應的水平速度值均在增大，且變化范圍較大；另一方面，較短的距離范圍內，不同頻率之間的非線性相互作用還未充分發(fā)展。因此，在線性聚焦位置為xf（5m,6m 和7m）的第1次聚焦之后，又出現(xiàn)了更大的聚焦波峰面，在本文計算中并未發(fā)現(xiàn)其后有第3 次聚焦大峰。而在其他算例（xf=8m，9m，10 m 和12 m）中，均在最大聚焦位置前，也出現(xiàn)了較大的聚焦波峰面（可以近似認為是局部最大波峰）。綜合來看，所有計算工況全域內的最大波峰面在74.42～81.27 mm 之間，這說明線性聚焦位置的設置對計算的最大波峰面產(chǎn)生了較大影響；類似地，所有計算工況中最大波峰面處的水平速度變化范圍為0.596～0.702 m/s之間，說明線性聚焦位置的設置對uf也產(chǎn)生了較大影響。此外，有趣的是式（10）給出的相鄰兩次最大波峰面間的傳播速度cf卻變化很小，其變化范圍在0.79～0.81 m/s之間。所有工況中，波峰面的最大速度為0.702 m/s，其小于最大波峰間的傳播速度0.81 m/s，這間接表明本文計算工況中波浪均未發(fā)生破碎。

3.3 波峰面最大值及水平速度的沿程變化

針對Baldock等[1]的D55工況，計算得到的波峰面及水平速度的沿程變化見圖9和圖10。從圖9中可知，以70 mm為界，存在較多的空間點波峰面超過該值，且存在一個較大的共振區(qū)域。類似地，圖10也表明，波峰面處的水平速度也存在一個較大區(qū)域內的共振現(xiàn)象。

圖9 波峰面最大值的沿程變化Fig.9 Maximum crest elevations along the flume

圖10 波峰面最大值處水平速度的沿程變化Fig.10 Horizontal velocity of maximum crest elevations along the flume

4 結 論

本文基于雙層Boussinesq 水波方程建立了聚焦波立面二維數(shù)值模型，將其應用于深水聚焦波的數(shù)值模擬中，得出的主要結論如下：

（1）數(shù)值模型計算線性聚焦波面處的波面位移與一階Stokes波解析解吻合程度很高，表明本文建立的雙層Boussinesq數(shù)值模型能夠精確捕捉線性相位問題。

（2）非線性聚焦演化中，波浪非線性起重要作用，數(shù)值模擬波面與實驗結果吻合較好，表明雙層Boussinesq數(shù)值模型能勝任聚焦波波面演化的模擬，也印證了數(shù)值模型對應的方程具有良好的非線性特性。此外，對波峰面下的水平速度場，本文建立的雙層Boussinesq數(shù)值模型也能較好地模擬。

（3）在強非線性波浪聚焦過程中，將會出現(xiàn)一個較大范圍的波浪聚焦區(qū)。以Baldock 等（1996）的D55 為例，在設置不同線性聚焦位置的情況下，最大波峰面的傳播速度保持在0.79～0.81 m/s 之間，這表明不同線性聚焦位置的設置對最大波峰面的傳播速度影響不大，但其對波峰面下的速度影響較大。