張振偉，劉必勁，鄒志利，傅丹娟

(1.廈門理工學院土木工程與建筑學院，福建 廈門 361024 ； 2.海南浙江大學研究院，海南 三亞 572025；3.大連理工大學海岸和近海工程國家重點實驗室，遼寧 大連 116024)

波浪在近岸區(qū)域的傳播運動是海岸動力學研究的重要內容之一，也是研究近岸環(huán)流、海岸演變等的基礎。波浪傳播運動過程中因受地形影響會產生折射、繞射、反射等物理現(xiàn)象。Boussinesq[1]在研究淺水波運動時建立了一維平底非線性Boussinesq方程；Peregrine[2]推導建立了適用于非平底的經典二維Boussinesq方程；Abbott等[3]首次應用Peregrine方程建立了二維Boussinesq方程計算模型，研究了港口內波浪傳播運動；Witting[4]針對經典Boussinesq方程存在的缺點，提出了Boussinesq型方程的適用范圍。Madsen等[5-6]建立了二階色散精度的Boussinesq方程，因該形式的方程具有一定精度且計算量相對較小，在實際中得到了一定的應用[7-8]。此后，許多學者不斷對Boussinesq型方程進行改進[9-16]，進一步提高了方程的適用范圍。國內有關學者也開展了Boussinesq型方程發(fā)展及應用的研究。張堯等[17]綜述了Boussinesq型方程的理論發(fā)展、實際應用和理論瓶頸，提出了Boussinesq型方程可能的科學突破方向。孫家文等[18]對1967—2018年Boussinesq型方程的理論研究進行了綜述，將Boussinesq型方程進一步分為水平二維和三維兩種情況。劉必勁等[19]應用簡化的雙層Boussinesq模型研究了聚焦波波面位移、波面處的水平速度和垂向速度，結果顯示模型計算結果與試驗測量結果符合很好。

Boussinesq型方程色散參數是通過分析方程色散關系與線性波理論色散關系的Padé展開給出的，而后通過線性變淺分析得到變淺參數。劉忠波等[20-22]針對Madsen等[10]提出的一組四階Boussinesq型方程，注意到色散參數取值不同可以影響方程的色散精度，于是使方程相速度與線性波相速度在無因次水深kh=15內誤差平方累計之和最小以優(yōu)化色散參數，計算表明這一方法可以顯著提高色散適用范圍；同時對四階Boussinesq型方程變淺系數重新進行了優(yōu)化，優(yōu)化后的參數顯著提高了方程的色散性；另外通過引入含參數的變換速度變量取代原方程中的速度變量，對3組著名的Boussinesq型方程[9，10，13]進行了改進，顯著提高了方程的變淺性能。Simarro等[23]對一些經典的Boussinesq型方程[10，13，24]進行了參數優(yōu)化，參數優(yōu)化后的方程控制了線性波的波幅誤差，變淺性能得到了改善。

文獻[11]的Boussinesq型方程是具有四階色散精度的水平二維方程，方程最高空間導數為三階，對數值計算要求相對較低。目前，工程上應用最廣的是二階Boussinesq型方程[7-8]，隨著計算能力的提高，四階Boussinesq型方程必將在實際工程中得到應用。另一方面，沒有針對文獻[11]中Boussinesq型方程參數優(yōu)化的相關研究，這里不注重提高文獻[11]中方程的適用水深范圍，而是保持原方程的適用水深范圍內，通過優(yōu)化參數提高方程模擬近岸水波運動的精度。本文在無因次水深kh=6.28條件下，將文獻[11]中Boussinesq型方程的相速度與線性波相速度之間誤差平方累計之和最小以優(yōu)化方程色散參數，而后根據變淺分析確定變淺參數，并分析了不同參數對方程色散和變淺性能的影響。最后通過數值模擬潛堤上波浪運動變形問題來驗證該方法的合理、有效性。

1 Boussinesq水波方程及其數值模型

1.1 Boussinesq水波方程

(1)

其中μ=kh

最終方程如下：

(2)

(3)

其中d=h+ηa=a1+a2b=b1+b2c=c1+c2

式中：k為波數；h為水深；a、b、c為待定參數；h為靜水水深；η為波面升高；d為總水深。

為使方程具有較好的色散和變淺作用性能，文獻[11]通過將方程色散關系與線性Stokes色散關系的Padé四階逼近一致得到色散參數，而后應用變淺性能的優(yōu)化得到變淺參數，最終參數(a1，a2，b1，b2，c1，c2)=(-0.006 7，-0.022 0，-0.001 7，0.000 5，0.039 3，0.061 3)。劉忠波[16]開發(fā)了六參數Boussinesq型方程的計算程序，應用六參數方程開展了潛堤上波浪傳播運動的模擬研究，在其計算程序中六參數取值為(a1，a2，b1，b2，c1，c2)=(-0.006 645 873，-0.022，-0.001 691，0.000 547，0.039 291，0.061 3)，與文獻[11]相比，文獻[16]中的部分參數保留了更多的有效數字。經重新計算和確認，這與文獻[16]中Luth等[26]算例采用的參數是一致的。

六參數方程的色散關系表達式、Stokes線性色散關系的Padé四階逼近表達式分別為

(4)

ω2=gks2h[1+1/9(ksh)2+1/945(ksh)4][1+4/9(ksh)2+1/63(ksh)4]-1

(5)

式中ks為Stokes線性波波數。

采用文獻[20]中的方法對六參數方程的色散性參數進行改進。對于從深水傳播到近岸淺水的波浪，當波浪傳播到近岸之后kh值不斷減小，取kh=6.28可以使方程考慮從深水到淺水的波浪傳播過程，即在kh≤6.28時，保證方程相速度與線性波相速度之間誤差平方累計之和最小。由此可以得到(a，b，c)= (-0.024 494 62，-0.000 764 2，0.095 389)。

圖1為不同參數的無因次相速度(cs為Stokes線性波相速度，c為Boussinesq型方程的相速度)對比。由圖1可見，kh=6.28時，最新優(yōu)化參數、文獻[11]參數及文獻[16]參數的色散誤差分別為0.069%、5.21%和3.09%。由此可知，優(yōu)化參數可以較大地提高波浪相速度的模擬精度。

圖1 無因次相速度比較Fig.1 Non-dimensional wave phase celerity comparisons

1.2 變淺性能分析

忽略方程的非線性項、(hx)2和hxx項，一維方程可寫成：

(6)

(7)

其中γ1=-b/aλ1=-(a+b/a)λ2=-(3a+5b/a+2a2+2b2/a)

η=A(x)exp[i(ωt-kx)]u=U(x)(1+iσhx)exp[i(ωt-kx)]

(8)

式中：A(x)為波幅；U(x)為速度振幅；ω為頻率；σ為考慮地形變化的小變量參數。將式(8)代入式(6)和式(7)并分離實部和虛部，可得關系式：

h(Ux+σkhxU)+γ1h2(2ωkAx+ωkxA)+λ1h3(-3k2Ux-3kkxU-σk3hxU)+

hx(U-λ2k2h2U+γ2ωkhA)=0

(9)

ωA-khU-γ1h2ωk2A+λ1k3h3U=0

(10)

-ωσhxU+λ3ωh2(2kUx+kxU+σk2hxU)+gAx+3cgh2(k2Ax+kkxA)+

hhx(4c2gk2A+λ4ωkU)=0

(11)

ωU-λ3ωk2h2U-gkA-cgk3h2A=0

(12)

式(10)和式(12)存在非零解的條件是A和U的系數矩陣行列式為0，可以得到六參數方程的色散關系。而后應用式(9)和式(11)推導得到變淺梯度的表達式。將Stokes線性波變淺系數與六參數方程的均方差從0至6.28積分，通過積分誤差最小確定變淺系數，由此可以得到六參數方程的6個參數(a1，a2，b1，b2，c1，c2)=(0.003 605 38，-0.028 1，-0.000 764 3，0.0，0.029 389，0.066 0)，線性變淺的詳細過程可以參看文獻[21，27]。kh=6.28時，最新優(yōu)化參數變淺與理論值偏差為0.008 4，較文獻[11]有了較大提高。

圖2為不同參數計算得到的變淺梯度。由圖2可見，優(yōu)化參數給出的變淺梯度在kh為0～6.28的范圍內與理論變淺符合較好。為進一步分析參數對方程變淺性能的影響，應用不同參數模擬了線性波在斜坡上的傳播運動：地形坡度為19∶600，斜坡深水水深10 m，另一端水深0.5 m，波浪周期為3.581 s，入射波高為0.1 m，kh=3.14、0.407。模擬時，去掉方程所有的非線性項，空間步長為0.05 m，時間步長dt=0.05 s，圖3為優(yōu)化參數、文獻[16]參數和文獻[11]參數的計算結果，由圖3可知，優(yōu)化參數給出的結果最好，其次是文獻[16]參數，再次是文獻[11]參數，這從線性變淺的數值方面證明了優(yōu)化參數可以更好地模擬波浪傳播運動。

圖2 不同參數的變淺梯度系數Fig.2 Linear shoaling gradient for different paramerers

圖3 波浪在斜坡上的變淺Fig.3 Wave shoaling along the slope

為了更詳細地研究參數的影響，分析六參數方程的二階諧波波幅，并將其與二階Stokes波二階波幅進行比較。利用Madsen等[10]給出的方法進一步考察參數的影響。

η=A1cosθ+εA2cos2θu=U1cosθ+εU2cos2θ

(13)

其中θ=ωt-kx

式中：A1、A2分別為一階和二階振幅；U1、U2分別為一階和二階速度振幅。

將式(13)代入式(1)和式(2)，分析ε的一階問題,則可以得到二階諧波的波幅，二階諧波方程為

(14)

其中

由式(14)可求得A2：

(15)

相應的Stokes二階非線性解二階諧波振幅為

(16)

圖4為六參數方程中不同參數對應的二階諧波波幅對比。從圖4可見，優(yōu)化參數對二階諧波基本沒有改善。Kennedy等[28]指出，方程參數優(yōu)化使方程線性性能在某一相對水深范圍內是可以接受的，但非線性誤差隨著相對水深的增加而快速變大。本文僅優(yōu)化了方程的色散和變淺參數，提高了方程的線性性能，沒有改進方程的非線性性能，所以二階非線性基本沒有得到改善。

圖4 方程的二階諧波波幅與二階Stokes波幅的比值Fig.4 Ratio of second order hammonic wave amplitude to second order stokes wave amplitude of equation

1.3 數值方法

六參數方程求解時時間差分為四階Adams-Bashforth-Moulton混合格式，空間差分為中心差分格式對式和式進行求解。預報時采用3階Adams-Bashforth預報格式計算式和式中波面和速度的預報值，然后利用4階Adams-Moulton時間步進格式，計算方程波面和速度的校正值。當變量的校正值與預報值之間的誤差小于等于0.000 1時，迭代終止。若不滿足迭代終止條件則用校正值和預報值的加權作為變量的新預報值，預報值和校正值的加權系數分別取0.05～0.1和0.95～0.90，可以不斷重復校正，直至滿足條件為止[16]。在數值模擬時，采用內部造波法產生波浪，波浪產生位置與模型試驗中造波機的位置一致，另外在計算域的左端和右端設置2倍波長的海綿層吸收波浪，阻止反射波的產生。

2 算 例

采用六參數方程數值模型，針對潛堤存在時波浪的傳播運動進行數值模擬，并與試驗的測量結果進行了比較，給出了采用不同參數的計算結果。

Luth等[26]和Ohyama等[29]分別進行了梯形潛堤上波浪傳播變形試驗，試驗布置見圖5。表1為Luth等[26]試驗和Ohyama等[29]試驗的入射波浪要素，數值計算中空間步長取0.025 m，時間步長取0.005 s，數值波浪水槽長55 m。

圖5 Luth等[26]和Ohyama等[29]試驗布置(單位：m)Fig.5 Experiments setup of Luth et al. and Ohyama et al. (Unit：m)

表1 Luth等[26]試驗和Ohyama等[29]試驗波浪要素

圖6和圖7分別為Luth等[26]試驗數值計算的波面升高與實測數據的對比。由圖6的計算結果可以看出，從x=2.0 m到x=14.5 m，采用3組參數計算得到的波面升高差別不大，但從x=15.7 m之后，采用3組參數計算得到的波面升高的差別顯現(xiàn)出來。采用本文優(yōu)化參數的計算結果略好于采用文獻[16]參數的計算結果。采用文獻[11]參數計算的結果最差，這一點可以從x=17.3 m到x=21.0 m波面升高的計算結果看出來，波浪經過潛堤時，不同頻次的鎖相波被釋放成為自由波，波面時間過程線會出現(xiàn)多種形態(tài)的次峰。若方程的色散性能較好，則可以較為準確地捕捉到各個鎖相波的位置，而方程變淺性能好則可以較準確地模擬波浪在斜坡上運動的波高演化。在x=19.0 m處，采用本文優(yōu)化參數計算的結果與波面升高符合最好，其次是采用文獻[16]計算的結果，采用文獻[11]參數計算得到的波面升高最差。

圖6 不同浪高儀處數值模擬結果與試驗數據的對比(Luth等[26]試驗算例(a))Fig.6 Comparisons between the computed results and experimental data at different locations(Luth et al.[26] experiment example(a))

圖7 不同浪高儀處數值模擬結果與試驗數據的對比(Luth等[26]試驗算例(c))Fig.7 Comparisons between the computed results and experimental data at different locations(Luth et al.[26] experiment example(c))

圖7為Luth等[26]試驗算例(c)的計算結果。由圖7可知，自x=2.0 m到x=19.0 m，采用3組參數計算的結果存在明顯差別，這是由于該算例對方程的色散性能要求較高。因文獻[11]參數的色散性最差，文獻[16]參數的色散性稍好，采用本文優(yōu)化參數的Boussinesq型方程的色散性最好。特別是越過潛堤之后(x=14.5 m至x=19.0 m)，采用文獻[11]參數計算結果與試驗結果差別最大，采用本文優(yōu)化的參數與文獻[16]參數計算的結果偏差相對較小，總體上采用本文參數波面升高的計算結果與試驗數據符合最好。

圖8為Ohyama等[29]試驗波面升高計算結果與試驗結果的對比。由Ohyama等[29]試驗的入射波要素說明3組試驗的非線性參數一致，色散性參數不一致，其中算例(2)的色散性參數最大。圖8給出了3號浪高儀和5號浪高儀位置處的波面升高計算結果與試驗結果的對比。由圖8可以看出，對于算例(2)，本文優(yōu)化參數所得的計算結果明顯優(yōu)于采用文獻[11]和文獻[16]參數的計算結果。在5號浪高儀位置處，3組參數給出的計算結果差別最大，采用文獻[11]參數的計算結果與試驗數據各個諧波的相位不相符合。對于3號浪高儀，隨著試驗工況波浪色散性參數的減小，采用3組參數計算的結果差別也會變小。算例(6)3號浪高儀處的計算結果差別最小，這一點在5號浪高儀處也有體現(xiàn)，即隨著試驗波況色散性參數的降低，采用3組參數給出的計算結果差別變小。

圖8 不同浪高儀處數值模擬結果與試驗數據的對比(Ohyama等[29]試驗)Fig.8 Comparisons between the computed results and experimental data at different locations(Ohyama et al.[29] experiment )

3 結 論

a.通過優(yōu)化方程的參數，可以提高方程的色散性。在kh=6.28時，本文優(yōu)化參數、文獻[11]參數及文獻[16]參數給出的色散性誤差分別為0.069%、5.21%和3.09%，可見本文優(yōu)化參數顯著提高了六參數方程的色散性。參數對Boussinesq型方程計算結果影響很大，因此在給出參數時要保留足夠的有效數字。變淺參數優(yōu)化后，當kh=6.28時，最新優(yōu)化參數變淺與理論值偏差為0.008 4，精度也有了提高。線性變淺的計算表明，采用本文優(yōu)化參數的Boussinesq型方程變淺性能得到了提高。二階諧波分析表明，優(yōu)化參數對二階諧波基本沒有改善。二階非線性性能的提高需要對方程的非線性進行改進。

b.采用3組參數針對潛堤上波浪的傳播變形進行了計算研究，結果表明采用本文優(yōu)化參數給出的計算結果與試驗數據更為符合，這也進一步說明參數對方程色散和變淺的重要性。