基于改進緩坡方程的波浪傳播數值模擬

王紅川，周正萍

(南京水利科學研究院 水文水資源與水利工程國家重點實驗室，江蘇 南京 210029)

波浪在由深海向近岸的傳播過程中，由于受到復雜地形、障礙物和水流等因素的影響，將發生淺化、折射、繞射、反射、底摩阻能量耗散以及破碎等一系列復雜現象，即波浪傳播變形，準確確定波浪傳播變形的波浪場是十分必要的。

20 世紀70年代以前波浪傳播問題都是把波浪折射和繞射分開研究的。Berkhoff［1］導出一類考慮波浪折射和繞射問題的緩坡方程(Mild-Slope Equation 簡寫MSE)方法，當kh ＜＜1 時MSE 簡化為淺水波方程，當kh ＞＞1 或kh=cons 時可以得到Helmholtz 方程。MSE 是在假定海底坡度為緩坡(μ=O(▽h/kh)=O(ε)＜＜1)基礎上將勢波理論的三維問題化為二維問題，使問題得到簡化，其精度為O(ε2)。Booij［2］詳細討論了緩坡方程在不同水底坡度條件下應用的精度，通過與三維模型計算結果比較，認為緩坡方程在海底坡度小于1/3 時可以保證使用的準確性。

如果海底地形變陡，這種現象在沙質海岸沿岸沙壩附近常常出現，有時地形變化是顯著的，海底地形曲率(≈▽2h)和底坡度的二階項((▽h)2)等將會對波浪傳播產生影響，而在Berkhoff 方程中這些項被忽略了。Kirby［3-4］將水深定義為緩變水深和復雜波動地形，緩變水深滿足緩坡方程的假設，將海底波動地形直接放入底邊界條件，在此基礎上發展了一類依時間的擴展緩坡方程，得到波浪在較陡地形上的傳播方程。Massel［5］使用Galerkin 方法構造問題的控制方程，進一步將勢函數描述為N 個與水深有關的正交函數Zn(x，y，z)上的分量，通過復雜推導得到描述波浪傳播變形的擴展的緩坡方程:

式中第三項括號內包含了海底坡度平方項和二階導數項。Chamberlain［6］采用變分原理導出包含海底坡度平方項和二階導數修正的緩坡方程，修正的緩坡方程表述為

式中:K=2kh。Porter［7］在Chamberlain 修正緩坡方程基礎上對海底地形不連續時通過考慮質量流的連續性，對此作進一步修正，提高了計算精度。

洪廣文等［8］根據無粘性、無旋流體動力學方程導出了一類基于重力表面波與長波非線性相互作用波浪傳播數學模型，適用于波浪在深水到極淺水、有長波流場與水位變化的水域傳播，控制方程中含有反映能量攝入、摩阻與破波損耗的能量系數因子和反映水底局部地形的因子。模型中海底地形變化的地形因子表達式用式(7)表示。

更多的學者在改進和擴展緩坡方程時考慮了較陡地形(▽h)2和▽2h 的影響，對復雜地形下的緩坡方程作了研究分析和數值計算［9-14］。研究結果表明，在緩坡方程中考慮水深一階導數平方和二階導數項后，對提高復雜地形和較陡坡地形(如海底波紋地形、Booij 實驗室斜坡地形等)的波浪場計算精度有明顯效果。

為更好地研究這些模型在復雜地形下的模擬精度和數值計算有效性，Kyung［14］對一個特殊淺灘地形進行了模型試驗，該地形淺灘上水深二階導數▽2h 的值為3.55，底坡度平方(▽h)2最大值達0.64，試驗中測量了淺灘附近多個不同斷面的波高，為數值模擬研究提供了較為詳細的參考數據。

從水動力學基本方程利用變分原理，推導出一種可考慮復雜地形下的改進的緩坡方程。考慮波動勢函數的沿水深分布是一個與未知波面有關的函數，對含自由面高度(η)的雙曲函數采用泰勒展式，略去O(η3)以上高階項，利用變分原理得到描述含有底坡度高階項和底坡曲率項的考慮波浪折射繞射傳播的緩坡方程，模型中描述陡變地形(▽h)2和▽2h 的系數僅與水深有關，因而計算時簡單方便，由于推導和簡化的方法不同，與一些前人得出的表達式也不同。導出的方程如略去底坡度高階項后與Berkhoff 線性假設導出的MSE是一致的，本模型是MSE 的改進模型。

1 方程的推導

1967年Luke 把變分法拓廣到具有自由水面的流體運動，后來很多學者也采用變分原理方法研究波浪問題。對于二維波動問題，考慮如下波能密度函數:

式中:φ 是波動勢函數，η 是自由面高度函數，h 是水深。假設:

這里f(η)=1，k 是特征值。當式(11)中忽略η，勢函數f(z)分布與線性情況下的分布函數一致，特征值k 代表波數。

將式(10)、(11)代入式(9)，考慮到以下關系式:

經過復雜計算，得到:

上式各項的系數中含有自由面參數η，為推導線性緩坡方程，將上式中雙曲函數泰勒展開，令

式中:σ=th(kh)，根據如下變分方程

式(19)、(20)中消去η，得

令FH= －F/(k2CCg)，~F=1 +FH，式(22)為

式(25)為導出的考慮復雜地形的改進緩坡方程(Modified Mild-Slope Equation，簡寫MMSE)，略去地形變化高階項FH，方程與MSE 方程有相同的形式。上述推導過程中對含自由面高度(η)的雙曲函數采用泰勒展式，導出的考慮地形高階項的參數與文獻［8］相比有所改進，兩者在形式上的不同參見附錄。

2 MMSE 數值計算

為數值求解方程(25)，假設

略去ε2項，將寫成t，上式變成如下依時間的拋物型方程，可將時間變量看作迭代參數，最后得:

3 水力模型試驗驗證

為驗證陡坡地形上波浪傳播變形，Kyung［14］進行了水力模型試驗，試驗在海岸工程實驗室進行，試驗港池尺寸為23 m×11 m×1 m(長×寬×高)。試驗地形為平底地形上放置一個圓形淺灘，如圖1。波浪自左側邊界(X= －6 m 處)形成，右側下游邊界(X=10.75 m 處)采用吸收邊界，圓形淺灘中心距離造波邊界6 m處，淺灘半徑R=0.45 m，離淺灘中心的灘面水深:

式中:h0是平底地形水深，h0= 0.3 m，淺灘中心處水深b=0.18 m。波高測量在5 個側向斷面(見圖1 中斷面X0～X4)和一個中心線斷面(見圖1 中斷面Y0)進行。

試驗采用規則波進行，入射波高為3 cm，波周期分別采用1.259 s、0.791 s 和0.636 s 三個波周期，對應k0h0=1、2、3，k0、h0分別是深水處波數和水深。

淺灘水深函數中有平方項，灘上水深的二階偏導數▽2h 為常數4b/R2=3.55，灘上一階導數平方(▽h)2的值為在r=0 處為0，在r=R 處達0.64。將式(25)寫成下列形式:

為分析MMSE 中▽2h、(▽h)2項的影響程度，可考察M1= －F1/k2CCg和M2= －F2/k2CCg隨水深變化的關系。圖2 顯示了當T=1.259 s 時M1和M2相對kh/π 的變化。從圖中可以看出在相對水深較大(kh/π＞1.5)時，M1和M2趨于0，此時▽2h、(▽h)2的影響很小。在相對水深較小(kh/π ＜0.1)時，M2趨近于0，(▽h)2項的影響較小;此時M1趨近于－0.167，▽2h 項不可忽略。在中等水深情況下，M1和M2均表現為不可忽略，M1的值在－0.167 ～0.047 之間，M2的值為－0.043 ～0.000 之間。

利用MMSE 模型對Kyung 地形上的波浪傳播進行數值計算。計算中空間步長取0.025 m，時間步長取0.002 s。計算波周期分別為1.259 s、0.791 s、0.636 s 三種，波高比等值線列于圖3。圖中粗虛線為地形中淺灘邊界，淺灘中心坐標為(0.0 m，0.0 m)，半徑為0.45 m。

圖1 Kyung 模型試驗布置及斷面位置示意Fig.1 Model layout and section location

圖2 M1、M2 隨kh 的變化曲線Fig.2 M1，M2 －kh curve

從波高比圖中可以看出，淺灘上或淺灘后波能集中、淺灘兩側的折射繞射、淺灘地形上的波浪反射等現象在數值模擬中都可以清晰表現。圖中還反映波能集中點隨波周期的減小向下游移動，其原因可能是由于波周期減小，淺灘上的波浪折射減弱。

圖3 波高比等值線Fig.3 Wave height contour

圖4 列出了淺灘中心線斷面(如圖1 中Y0)MMSE 模式和傳統MSE 模式的波高計算結果與實測資料的比較。分別對3 種不同的波周期進行了試驗，圖中MMSE 模式和MSE 模式計算結果分別用實線和虛線表示，實圓點代表試驗結果。從圖中可以看出，MMSE 模式計算結果與試驗結果有相當好的一致性，特別是淺灘附近出現的最大波高及其位置與實驗結果吻合得相當好，而MSE 模式則嚴重偏離試驗結果。從圖4(a)可以看出，當T=1.259 s 時MSE 模式計算的最大波高位置向下游漂移。

圖4 中心線斷面(Y0)數值模擬與試驗結果比較Fig.4 Comparison between simulated and tested results at center section (Y0)

圖5 側向斷面數模計算結果與試驗結果的比較(T=1.259 s)Fig.5 Comparison between simulated and tested results at side section (T=1.259 s)

還對Kyung 地形中幾個側向斷面的波高計算結果作了比較，圖5 ～圖7 分別是波周期為1.259 s、0.791 s、0.636 s 五個斷面(如圖1 中X0、X1、X2、X3、X4)的計算值與試驗結果的比較。從圖中可以看出，MMSE 模式計算的各斷面的波高值與試驗結果吻合較好，而MSE 模式計算結果與試驗結果偏差明顯，特別是長周期波底坡效應變得十分明顯，MSE 模式在波能集中的位置計算結果偏大，在發散區計算波高偏小。波浪在淺灘開始處(X= －R)波高側向變化較小，隨著波浪向淺灘傳播，波高的側向變化劇烈，其影響直接延伸到淺灘下游水域。

圖6 側向斷面數模計算結果與試驗結果的比較(T=0.791 s)Fig.6 Comparison between simulated and tested results at side section (T=0.791 s)

圖7 側向斷面數模計算結果與試驗結果的比較(T=0.636 s)Fig.7 Comparison between simulated and tested results at side section (T=0.636 s)

4 結 語

通過理論分析從水動力學基本方程出發，利用變分原理推導出復雜地形上的波浪傳播變形的改進的緩坡方程數學模型，主要結論如下:

1)考慮波動勢函數的沿水深分布是一個與未知波面有關的函數，對含自由面高度(η)的雙曲函數采用泰勒展開，略去O(η3)以上高階項，利用變分原理得到描述含有底坡高階項和底坡曲率項的波浪傳播方程。

2)與實驗室復雜地形實測資料驗證對比可以看出，MMSE 模型可以清晰地表現淺灘上或淺灘后波能集中、淺灘兩側的折射繞射、淺灘地形上的波浪反射等現象，計算的各斷面的波高值與試驗結果吻合較好，能夠準確地反映淺灘附近出現的最大波高及其位置。

3)由于波浪在淺灘開始處波高側向變化較小，隨著波浪向淺灘傳播，波高的側向變化劇烈，其影響直接延伸到淺灘下游水域，MMSE 模式的模擬精度可明顯提高。改進后的模型可以提高復雜地形下的波浪場計算精度。

［1］Berkhoff J C W.Computation of combined refraction-diffraction［C］//Proceedings of the 13th International Conference on Coastal Engineering，Vancouver.1972:745-747.

［2］Booij N.A note on the accuracy of the mild-slope equation［J］.Coastal Engineering，1983，7:191-203.

［3］James T Kirby.A general wave equation for wave over rippled beds［J］.Journal of Fluid Mechanics，1986，162:171-186.

［4］James T Kirby.On the gradual refraction of weakly nonlinear stokes waves in regions with varying topography［J］.Journal of Fluid Mechanics.1986，162:187-209.

［5］Massel S R.Extended refraction-diffraction equation for surface waves［J］.Coastal Engineering，1993，19:97-126.

［6］Chamberlain P G，Porter D.The modified mild-slope equation［J］.Journal of Fluid Mechanics，1995，291:393-407.

［7］Porter D，Staziker D J.Extensions of the mild-slope equation［J］.Journal of Fluid Mechanics，1995，300:367-382.

［8］洪廣文.長波上非線性重力表面波傳播數學模型［C］//第十四屆中國海洋(岸)工程學術討論會論文集.北京:海洋出版社，2009.

［9］Suh K D，Lee C，Park W S.Time-dependent equations for wave propagation on rapidly varying topography［J］.Coastal Engineering，1997，32:91-117.

［10］Ting-Kuei Tsay，Philip L F Liu，Nan-Jing Wu.A nonlinear model for wave propagation［C］//Proceedings of the International Conference on Coastal Engineering.1996.

［11］Zhang L，Edge B L.A uniform mild-slope model for waves over varying bottom［C］//Proceedings of the 25th International Conference on Coastal Engineering，ASCE.1996:941-954.

［12］Changhonn Lee，Gunwoo Kim，Kyung-Duck Suh.Extended mild-slope equation for random waves［J］.Coastal Engineering，2003，48:277-287.

［13］Chandrasekera C N，Cheung K F.Extended linear refraction-diffraction model［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，1997，123:280-286.

［14］Kyung D S，Changhoon Lee，Young-Hyun park，et al.Experimental verification of horizontal two-dimensional modified mildslope equation model［J］.Coastal Engineering，2001，44:1-12.

［15］李孟國，蔣德才.關于波浪緩坡方程的研究［J］.海洋通報，1999，18(4):70-91.

［16］李瑞杰.考慮非線性彌散影響的波浪變形數學模型［J］.海洋學報，2001，23(1):102-107.

［17］王紅川，左其華，潘軍寧.波浪傳播數值模型波向角計算［J］.水動力學研究與進展:A 輯，2006，21(1):139-144.

［18］林 鋼，邱大洪.波浪在雙灘地形上的傳播［J］.水利學報，1999(8):57-60.

［19］潘軍寧，洪廣文，左其華.一種推廣的緩坡方程［J］.海洋工程，2001，19(1):24-31.

［20］潘軍寧，左其華，王登婷.港域波浪數學模型的改進與驗證［J］.海洋工程，2008，26(2):34-42.

附錄:

與文獻［8］推導的地形二階項▽h 系數比較

這里所導出的地形參數表達式用F=F1▽2h+F2(▽h)2

由彌散關系ω2=gkth(kh)，可得

F1的表達式在形式上與文獻［8］推導的▽2h 系數R1相同，文獻［8］給出: