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

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

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

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

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

1 數(shù)值模型

1.1 控制方程

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

(1)

(2)

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

(3)

(4)

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

(5)

(6)

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

(7)

(8)

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

(9)

(10)

式中：下標x表示變量對空間x求導；g為重力加速度；h是靜止水平面下的水深；色散系數(shù)α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位于區(qū)間[-h, -2α1h]時，采用表達式(11)、(12)進行求解；當z位于區(qū)間[-2α1h,0]，采用(13)、(14)求解；當波面高于靜止水位z=0時，采用(15)、(16)計算速度。

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

本文模型最高空間導數(shù)是2，比開源程序FUNWAVE對應方程的空間導數(shù)少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內，當前計算結束，否則重新賦值給各個變量后繼續(xù)進行迭代求解。

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

為了避免地形反射波浪引起造波板的二次反射，在數(shù)值模型中，引入了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為模型布置示意圖。試驗采用規(guī)則波，研究迎浪面和背浪面坡度均為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為間隔)入射條件，觀察潛堤上波浪是否會發(fā)生破碎，從而確定出非破碎強非線性波浪工況，具體見表1，表中k=2π/L，L為波長，A為波幅。

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

3 結果比較和討論

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

3.1 波面位移和諧波幅值

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

圖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倍頻。進一步統(tǒng)計潛堤后方(x≥19.1 m)7個測點處的各次諧波幅值，并取各點平均值，結果顯示，周期1.4 s工況中2倍、3倍和4倍頻諧波的幅值平均值依次為基頻的0.74、0.11和0.08倍，周期1.8 s工況中該倍數(shù)依次為1.50、0.49和0.13，周期2.2 s工況中該倍數(shù)依次為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)進行濾波處理。將計算波浪流速與濾波后試驗數(shù)據(jù)進行比較，詳見圖11～圖13。由圖可見，無論是水平速度還是垂向速度，兩者流速過程的相位和峰谷值均吻合良好。

3.3 波浪非線性特征分析

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

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

4 結論

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

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

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

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