雙層Boussinesq模型非線性波浪模擬研究

饒永紅，劉忠波，梁書秀，李 欣

(1.中國人民解放軍91053部隊，青島 266100； 2.大連海事大學 交通運輸工程學院，大連 116026； 3.大連理工大學 海岸和近海工程國家重點實驗室，大連 116023)

波浪是近岸水域主要的水動力，掌握波浪資料對建筑物結構設計、泥沙運動、海床演變十分重要。作為一種水工構筑物，潛堤常用于促淤保灘，潛堤的存在使局部水深減少，淺水更有利于產生高次諧波，過淺時波浪將發生破碎，損耗大量波能，減小了堤后波高。波浪與潛堤相互作用問題的研究是海岸工程界和學界關注的課題，工程界多關注采用什么樣的結構型式更有利于減少堤后波浪，學界多關注如何精準地描述波浪在潛堤上的演變過程。

潛堤上的波浪傳播變形是復雜的水動力過程，其涉及波浪反射、變淺效應、反變淺效應甚至波浪破碎等現象。為了研究潛堤上的波浪傳播變形特性，眾多國內外學者利用理論解析、物理模型試驗或數值模型模擬等多種方式開展研究。Beji和Battjes[1]、Luth等[2]、Ohyama等[3]、鄒志利等[4]開展了不同的物理模型試驗來研究波浪傳播變形，試驗數據被用來驗證各類數值模型。時瑩等[5]基于前坡1:15、后坡1:5潛堤的物理模型試驗，考察了SWAN 模型對潛堤上波浪演化的模擬能力，發現譜平均周期計算值明顯偏小，頻譜模擬出現高頻高估、低頻低估現象，三波非線性項的計算精確度影響了模型在近岸淺水區譜平均周期和頻譜的模擬能力，需更多改進。史寶凱等[6]基于δ-SPH方法，通過在控制方程中引入介質阻力項，建立波浪與多孔潛堤相互作用的數值波浪水槽，數值模擬與物理模型試驗結果較好。基于典型的勢流理論，Boussinesq方程是一種常用的波浪數值模擬工具，然而不同的Boussinesq方程模型，存在多種表達形式，理論、性能存在明顯差異，所以，模擬波浪水動力的精度也不同(Madsen和Fuhrman[7]；Brocchini[8]；Kirby[9]；張堯等[10]；琚海軍等[11]；孫家文等[12])。波浪演化研究的重點為波面位移和波浪速度，前者用于計算波高、平均水位(增減水)和非線性波浪參數，后者是計算水工建筑物波浪力荷載的依據。精確模擬波面演變的全過程，需要數值模型具有優良的色散性能(相位)、非線性(描述和差頻和高次諧波波幅)和變淺性能(變淺波幅)，而精確模擬波浪速度場則要求模型具有良好的速度垂向分布特征。Liu等[13]給出的多層Boussinesq方程，在色散性、非線性以及速度分布特征等方面理論適用水深得到較大的拓展，該模型已被應用到滑坡興波、淺水非線性長波、聚焦波、深水雙色波、流函數強非線性波浪和陡坡潛堤上的波浪傳播變形等[14-17]。

然而，多數Boussinesq方程模型以波面位移為主要研究內容，對波浪水平速度和垂向速度未進行深入研究，在研究潛堤上波浪傳播問題時，也僅以相關試驗來考察模型的適用性。Liu等的多層Boussinesq方程[16]，理論上具有良好的線性和非線性性能，但研究僅限于用一些理想情況或以試驗為依據討論方程的適用性。本文以Liu等的雙層Boussinesq方程模型[13]為研究對象，考察模擬1:4坡度潛堤上非破碎波浪演變的能力；為驗證模型精度，開展了相應的物理模型試驗。研究內容包括：(1)該模型對潛堤上強非線性波浪波面演化的模擬能力；(2)方程能否精確模擬波浪速度場。

1 數值模型

1.1 控制方程

(1)在自由水面處，滿足以下運動學和動力學邊界條件

(1)

(2)

(2)自由面處與靜止水位處的速度關系表達式為

(3)

(4)

(3)靜止水位處的速度與上層計算速度之間滿足關系式

(5)

(6)

(4)在連接上下兩層水體處，滿足速度相等條件，即

(7)

(8)

(5)在水底處，海底為不可穿透的運動學邊界條件，即

(9)

(10)

式中：下標x表示變量對空間x求導；g為重力加速度；h是靜止水平面下的水深；色散系數α1和α2分別取值為0.172和0.328。

1.2 任一水深處的水平速度和垂向速度

(11)

(12)

(13)

(14)

(15)

(16)

式中：zα1=-α1h；zα2=-(2α1+α2)h

表達式(11)、(12)至(15)、(16)給出了水平和垂向速度的表達式。當水深z位于區間[-h, -2α1h]時，采用表達式(11)、(12)進行求解；當z位于區間[-2α1h,0]，采用(13)、(14)求解；當波面高于靜止水位z=0時，采用(15)、(16)計算速度。

1.3 數值模型的求解

本文模型最高空間導數是2，比開源程序FUNWAVE對應方程的空間導數少1[18]。為此，本文模型在時間格式上采用與FUNWAVE源程序一樣的混合4階Adams-Bashforth-Moulton(ABM)格式。在預報階段，先利用顯示3階Adams-Bashforth格式計算方程(1)與(2)，可獲得水面位移和自由水面處水平速度的預報值；進一步依次求解方程(3)、(5)和(7)可得3個水平速度的預報值，進而求解方程(9)、(8)、(6)和(4)得到4個垂向速度的預報值。在校正階段，采用4階Adams-Moulton時間差分格式求解方程(1)和(2)，可求得水面位移和自由表面處水平速度的校正值，而其他速度校正值的求解過程類似預報階段。當所有變量的校正值與預報值的誤差控制在0.000 1內，當前計算結束，否則重新賦值給各個變量后繼續進行迭代求解。

圖1 方程的無因次相速度Fig.1 Non-dimensional phase celerity of the model

為了避免地形反射波浪引起造波板的二次反射，在數值模型中，引入了Hsiao等[19]的內波造波技術，且在計算域兩端邊界設置海綿層來吸收波浪[18]。

1.4 模型的理論性能

為了展示方程相速度與速度場的模擬精度，圖1給出方程的無因次相速度，圖2給出kh=0.5、1、3和6時方程解析速度場與Stokes線性波解析解的比較，對應的沿水深積分速度誤差見圖3。其他非線性特征的計算精度可參閱Liu等的文獻[16]。

圖2 方程的速度分布Fig.2 Velocity profile of the model

圖3 沿水深積分的速度誤差Fig.3 Error of the integrated velocity along vertical water depth

表1 波浪要素Tab.1 Wave parameters

2 物理模型簡介

物理模型試驗在長50 m、寬1.2 m、深1.2 m的波浪水槽中進行，圖4為模型布置示意圖。試驗采用規則波，研究迎浪面和背浪面坡度均為1:4潛堤上的波浪傳播變形。水槽水深h為0.5 m，潛堤上方水深h1為0.15 m，潛堤前趾距離造波板15 m。浪高儀設置在距離造波板x=12 m、15.7 m、16.7 m、18.1 m、19.1 m、20.1 m、21.1 m、22.1 m、23.1 m、24.1 m和 25.1 m處。試驗中采用超聲多普勒流速儀(ADV)測量瞬時流速，垂向方向上布設在靜止水位下0.06 m處，水平位置至造波板的距離依次為x=17.4 m、19.1 m和20.232 m。波面位移采樣間隔0.017～0.024 s，波浪速度的采樣間隔為0.025 s，設置約采集10個波的樣本長度。試驗中，通過多次重復不同波高(以0.005 m為間隔)入射條件，觀察潛堤上波浪是否會發生破碎，從而確定出非破碎強非線性波浪工況，具體見表1，表中k=2π/L，L為波長，A為波幅。

圖4 試驗布置示意圖(單位：m)Fig.4 Sketch of the experimental layout

3 結果比較和討論

數值模擬中，整個計算長度為60 m，采用0.01 s的時間步長和0.03 m的空間步長。內部造波的中心線設在x=10 m處，對應于試驗造波機的0點位置，模擬中前后邊界處均設置9 m長海綿邊界層吸收波浪。

3.1 波面位移和諧波幅值

以x=12 m處測得的波面位移作為表1入射條件的校核標準，數值模擬中的波浪要素與表1給定的目標值一致。波面位移計算值與試驗結果的對比見圖5、圖6和圖7，兩者符合程度高。觀察其中x=15.7 m和x=16.7 m處的波面位移過程可知，規則波在潛堤前坡上傳播時，波浪以鎖相波形態存在；而在平坡段(潛堤上方)常水深傳播時，淺水容易釋放出更多高次諧波，它們在潛堤后坡(x=19.1 m)由于水深增加變為自由波，并以各自相速度傳播，因此，坡后出現更為復雜的波浪形態。此外，伴隨波浪周期的增大，潛堤上方無因次水深減小，更容易激發出更多的高次諧波，圖7中x=18.1 m處出現明顯的三次峰就是佐證。

圖5 波面時間歷程對比 (T=1.4 s, H=0.055 m)Fig.5 Comparisons of the computed surface elevation and the experimental data (T=1.4 s, H=0.055 m)

圖6 波面時間歷程對比 (T=1.8 s, H=0.05 m)Fig.6 Comparisons of the computed surface elevation and the experimental data (T=1.8 s, H=0.05 m)

圖7 波面時間歷程對比 (T=2.2 s, H=0.045 m)Fig.7 Comparisons of the computed surface elevation and the experimental data (T=2.2 s, H=0.045 m)

對三組工況的波面位移進行傅里葉分析，得到各次諧波幅值見圖8、圖9和圖10。由圖可見，整體上和試驗結果吻合良好，且隨著波浪周期逐漸增大，高次諧波的幅值快速增加。在潛堤后方(x≥19.1 m)，工況1幅值大小順序是基頻、2倍頻、3倍頻和4倍頻；工況2幅值大小順序是2倍頻、基頻、3倍頻和4倍頻；而工況3幅值大小順序是3倍頻、基頻、2倍頻和4倍頻。進一步統計潛堤后方(x≥19.1 m)7個測點處的各次諧波幅值，并取各點平均值，結果顯示，周期1.4 s工況中2倍、3倍和4倍頻諧波的幅值平均值依次為基頻的0.74、0.11和0.08倍，周期1.8 s工況中該倍數依次為1.50、0.49和0.13，周期2.2 s工況中該倍數依次為0.71、1.31和0.38。其中，工況3的波浪周期和波長最大。堤上相對水深最淺，說明淺水更有利于波-波能量由低頻向高頻傳遞，導致堤后3次諧波幅值最大。

圖13 計算流速與試驗結果的比較 (T=2.2 s, H=0.045 m)Fig.13 Comparisons of computed velocity and the experimental data(T=2.2 s, H=0.045 m)

3.2 波浪速度

物理模型試驗測量了波浪水平速度和垂向速度過程，利用三點光滑(中間值權重0.5，相鄰點權重為0.25)進行濾波處理。將計算波浪流速與濾波后試驗數據進行比較，詳見圖11～圖13。由圖可見，無論是水平速度還是垂向速度，兩者流速過程的相位和峰谷值均吻合良好。

3.3 波浪非線性特征分析

對表1的三組工況，采用數值模型計算了最大波峰面高度(ηmax)、最大波高(Hmax)、波長(L1，L2，L3)和最大波峰處的水平速度(uη)。這里L1為淺水線性波長，用Stokes線性波色散關系求解得到；L2為淺水非線性波長，用數值模擬結果分析得到；L3為與流函數波浪對應的波長，分別將1.4 s、1.8 s和2.2 s周期對應的Hmax輸入流函數波浪波長數值逼近程序計算得到。計算結果見表2。表中，Hmax/h1為淺水非線性參數，U=HmaxL2/h13為厄塞爾數。由表可見，在潛堤上的平坡淺水水域，波高與水深比在0.539～0.599范圍內，屬于強非線性波浪，厄塞爾數在68.46～229.79，當每組工況的特征參數超過表2所列臨界值之后，波浪將發生破碎。

表2 波浪特征和非線性參數臨界值Tab.2 Computed results of wave features and wave nonlinear parameters

4 結論

本文針對雙層高階Boussinesq方程，建立了基于預報-校正-迭代的有限差分數值模型。利用數值模型，模擬了迎浪面和背浪面坡度均為1:4潛堤上強非線性波浪的傳播變形，給出了不同位置點的波面位移和速度(傳播方向上的水平速度和垂向速度)時間歷程。為探究模型的模擬能力，開展了相應的波浪物理模型試驗，并確定了強非線性波浪工況。通過以上研究，主要結論如下：

(1)在潛堤上淺水區域，波高與水深之比在0.539～0.599，屬于強非線性波浪。本文試驗工況的厄塞爾數在68.46～229.79。

(2)針對不同周期的強非線性波浪工況，雙層Boussinesq方程模型計算了堤前、堤頂和堤后的波面位移、水平速度、垂向速度，與試驗結果吻合程度較高，充分展現了數值模型良好的非線性波浪模擬能力。

(3)隨著波浪周期逐漸增大，高次諧波的幅值增加較快，在堤后超過了基頻波的波幅，尤其是周期2.2 s工況，3倍頻諧波幅值最大，超過了基頻和2倍頻，這充分說明了淺水更有利于高次諧波的產生。