改進的全水深非線性波傳播模型及其模擬

張 俞,洪廣文,馮衛(wèi)兵

(1.河海大學海岸與海洋工程研究所,江蘇南京 210098;2.河海大學港口海岸與近海工程學院,江蘇南京 210098)

目前波浪傳播模型主要沿著完全頻散性的線性緩坡方程(MSE)[1]及其擴展(EMSE)[2-5]與弱頻散性的非線性Boussinesq型方程[6]及其擴展[7-8]2種典型模型不斷發(fā)展.完全頻散性的線性緩坡方程及其擴展模型主要應用于大范圍水域波浪傳播計算,其優(yōu)點是計算簡便,可計算大范圍水域,缺點是對淺水域、非線性作用強的區(qū)域特別是破碎區(qū)域難以有效模擬.Boussinesq型方程及其擴展模型主要用于淺水水域波浪的變形計算,能夠很好地反映其非線性特性,但是計算難度大,較難應用于深水域的波浪計算.為此,能夠同時兼顧頻散特性和淺水非線性特性,較全面、精確模擬波浪傳播的全水深方程波浪傳播模型成為新型的波浪模擬途徑.

全水深方程波浪傳播模型最早由Nadaoka等[9-10]提出,但其只適用于無流情況且未考慮能耗與地形因子影響.洪廣文等[11-12]基于格林公式和變分法提出了含有能量系數的水流中聯合折射-繞射線性模型,隨后擴展為緩變流場、水位和水深的非線性模型,并證明在淺水域可以轉化為Boussinesq型方程、KdV方程、Airy淺水波方程,在深水域可化為非線性緩坡方程.吳中等[13]利用該非線性模型進行了較系統的計算驗證;張俞等[14]也利用其計算了考慮波浪破碎及摩阻損耗的單坡地形上波浪傳播,考慮了反映地形坡度與曲率因子J的影響.

本文在長波上非線性重力表面波傳播數學模型基礎上建立了保留了三階項的數學模型及其求解模式,并且在方程中(kη)2項近似取為2.采用Berkhoff試驗地形和用于驗證Boussinesq型方程模型的單坡地形進行模擬計算,并將結果與試驗值及線性、二階模型結果相比較,驗證三階全水深非線性數學模型的計算精度及適用性.

1 全水深非線性波傳播數值模型

1.1 三階非線性數學模型

基于長波上非線性重力表面波傳播數學模型[12],保留三階項可得無流情況下三階非線性數學模型,控制方程如下:

其中

式中:x,y——平面坐標,置于靜水面處;z——垂直坐標,豎直向上為正;h——當地水深;Φ——二維復值波動勢函數;η——復值波面函數;W*——能量系數,包括能量輸入與能量損耗(包括摩阻損耗和破波損耗);J——反映水底地形坡度與曲率有關的參數[12,14];k——波數,～c,～cg,～σ——波速、波群速和考慮能耗的相對圓頻率;ω——絕對圓頻率;g——重力加速度.

1.2 二階非線性模型

1.3 線性傳播數值模型

略去高階項,由方程(1),(2)可得如下線性模式控制方程:

1.4 邊界條件

為了使得入射條件與非線性模型控制方程組相匹配,本文的入射邊界采用Stokes二階形式:

式中:y=y0——入射邊界;θ0——入射波的波向與 x軸的夾角;a——入射波的振幅;k0——波數.出流邊界采用消波層消波與統一出流邊界相結合的方式[14].

1.5 初值條件

初始情況(t=0)除入射邊界外,其余水域 Φ,η均為0.

1.6 控制方程的數值解法

采用可變網格的預測-校正-迭代的Crank-Nicolson算法[14].為了增加計算的穩(wěn)定性,本文程序設計上采用將線性項與非線性項分離,先由線性模型計算至穩(wěn)定,再將線性結果作為非線性初值進行迭代計算.

2 模型驗證與分析

2.1 橢圓淺灘地形

1982年,Berkhoff[1]在水力實驗室針對均勻斜坡上布置一橢圓形淺灘的地形進行了波浪物理模型試驗,以驗證線性緩坡方程的計算結果.試驗采用的波周期為T=1.0s,入射波高為H0=0.0464m,并取得了8個斷面的實測資料.

Berkhoff試驗地形如圖1所示,本文的計算區(qū)域為 0 m＜x＜21.5m,0m＜y＜20m,波浪垂直y斷面正向入射.將該計算區(qū)域坐標轉換到橢圓中心,可得

在新坐標系中斜坡上水深:

斜坡加上橢圓淺灘后的地形為

圖1 Berkhoff橢圓淺灘地形試驗水深等值線及斷面分布Fig.1 Sketch of Berkhoff's test bathymetry

模型中空間網格步長為Δx=Δy=0.1,時間步長取Δt=T/60,計算總時長為32T.非線性模式入射邊界 Φ采用Stokes二階波的形式給定,線性模式 Φ采用Stokes一階波的形式 Φ1給定,即 Φ=Φ1;左右邊界及出流邊界均采用消波層與統一邊界條件進行處理,消波層長度為波長的2倍;計算達到穩(wěn)定后對各計算點波面時間過程線采用上跨零點法求取波高.各斷面的相對波高結果見圖2(圖中H為波高,H0為各組況下的入射波高).圖3和圖4為整個水域的波高等值線圖及波面影像圖.

圖2 1～8號斷面相對波高計算值與試驗值比較Fig.2 Comparison between measured and computed relative wave heights for section No.1—No.8

圖3 數值模擬的相對波高等值線(單位:m)Fig.3 Contours of simulated relative waveheights(unit:m)

圖4 t=32T瞬時波面影像Fig.4 Instantaneous wavesurface elevation for t=32T

由圖2可見,本文模型的結果優(yōu)于二階模型及線性模型結果,特別是在水深較淺的5號以及6～8號的淺水處.由圖3和圖4明顯可以看出,隨著水深變淺,波向朝著垂直于等深線的方向傳播,且在橢圓地形附近發(fā)生輻聚,符合波浪的變形趨勢.

2.2 單坡地形

采用文獻[7]中為驗證Boussinesq型方程擴展模型進行的單坡水槽波浪傳播試驗資料作為驗證依據,其試驗模型布置如圖5所示.該水槽有效長度30m、水深0.36m,在水槽末端設1∶34.26的斜坡.入射波采用表1中給出的5種破波形態(tài)的波要素.

圖5 Hansen與 Svendesen試驗地形Fig.5 Layout of Hansen and Svendesen's test bathymetry

表1 Hansen與Svendesn試驗入射波參數Table1 Parameters of incident waves for Hansen and Svendesen's tests

本文模型中入射處采用非線性入射條件,出流邊界處采用消波層與統一邊界條件相結合的方式處理.各模型中的空間步長隨該水深處對應的波長呈連續(xù)性變化,Li/Δxi≥7,Δt=T/50～T/100.

波浪破碎采用試算的方法先確定最初的破碎點.破碎點確定后,破碎點后采用公式

式中 ω為絕對圓頻率;系數S參照文獻[15]和Bettjes等[16]破碎公式調試確定.對于表1中所取5組波要素,S依次取0.6496,0.6496,0.4974,0.2210,0.1240,破碎點前=0,重新進行整個水域 Φ,η的計算,得到波高值.破碎后的波參數與試驗結果吻合良好.

由圖6結果可以看出,斜坡上水深較深段,線性模型結果尚能滿足要求,但隨著水深變淺,線性模型結果明顯偏小,不能有效體現波浪的傳播變形規(guī)律;本文模型與二階模型計算結果相當,與試驗值吻合良好,說明該三階修正是合理的.值得一提的是,Kirby對線性緩坡方程進行非線性擴展[2],其模型對Berkhoff試驗地形驗證后得到了很好的結果,但是該修正用于計算單坡地形,結果并不理想.

圖6 計算波高與試驗值比較Fig.6 Comparison between measured and computed wave heights

3 結 論

a.本文模型能同時兼顧非線性與頻散性特性,模型中W*能考慮能量的輸入與耗散.

b.本文模型對“長波上非線性重力表面波傳播數學模型”保留至3階精度,并用2替代 η2,采用含有松弛因子的Crank-Nicklson算法,入射邊界采用Stokes非線性波輸入改善與模式的匹配性,出流條件采用消波層與統一邊界條件相結合的方法,模型計算后獲得穩(wěn)定的波高結果.

c.Berkhoff試驗地形驗證結果表明,本文模型與試驗值吻合度有所提高,其效果類似于Kirby修正;單坡地形驗證結果表明,本文模型能模擬波浪由深水至淺水直至破碎的過程,所作的修正是合理的.

[1]BERKHOFF J C W.Computation of combined refraction-diffraction[C]//Proceeding 13thConference on Coastal Engineering.Vancouver:ASCE,1972:471-490.

[2]KIRBY JT,DALRYMPLE RA.Verification of a parabolic equationfor propagation of weakly nonlinear waves[J].Coastal Eng,1984,8:219-232.

[3]HSUT W,WEN CC.A parabolic equation extended to account for rapidly varying topography[J].Ocean Engineering,2001,28(11):1479-1498.

[4]PAN Jun-ning,HONGGuang-wen,ZUOQi-hua,et al.Verification of a numerical harbor wavemodel[J].China Ocean Engineering,2007,21(4):647-658.

[5]張揚,李瑞杰,張素香.波浪Boussinesq方程與緩坡方程的比較[J].河海大學學報:自然科學版,2005,33(5):588-591.(ZHANG Yang,LI Rui-jie,ZHANG Su-xiang.Comparasion of boussinesq equation with mild slope equation[J].Journal of Hohai University:Natural Science,2005,33(5):588-591.(in Chinese))

[6]PEREGRINEDH.Long waves on a beach slope[J].JFluid Mech,1967,27(4):815-827.

[7]KENNEY A B,CHEN Q,KIRBY JM,et al.Boussinesq modeling of wave transformation,breaking and runup.I:1D[J].Journal of Waterway,Port,Coastal and Ocean Engineering,2000,126(1):39-47.

[8]LI B.Wave equations for regular and irregular water wave propagation[J].Journal of Waterway,Port,Coastal and Ocean Engineering,2008,134(2):121-142.

[9]NADAOKA K,BEJI S,NAKAGAWA Y.A fully-dispersive nonlinear wave model and its numerical solutions[C]//Proceeding of the 24thInternational Conference on Coastal Engineeing.Kobe:ASCE,1994:427-441.

[10]ISOBE M.Time-dependent mild slope equation for random waves[C]//Proceeding of the 24thInternational Conference on Coastal Engineering.Kobe:ASCE,1994:285-299.

[11]HONGGuang-wen.Mathematical models for combined refraction-diffraction of wave on non-uniform current and depth[J].China Ocean Engineering,1996,10(4):433-454.

[12]洪廣文,吳中,張俞.長波上非線性重力表面波傳播數學模型[C]//第14屆全國海岸工程學術討論會論文集:上.呼和浩特:海洋出版社,2009:21-34.

[13]吳中,洪廣文,張洪生,等.完全頻散性非線性波數值模型及應用[C]//第14屆全國海岸工程學術討論會論文集:上.呼和浩特:海洋出版社,2009:448-457.

[14]張俞,洪廣文,馮衛(wèi)兵,等.考慮能量損耗的全水非線性波傳播[C]//第14屆全國海岸工程學術討論會論文集:上.呼和浩特:海洋出版社,2009:465-472.

[15]DALLY W R,DEAN RG,DALRYMPLE R A.Wave height variation across beaches of arbitrary profile[J].JGeophys Res,1985,90(C6):11917-11927.

[16]BATTJESJA,JANSSEN J PF.Energy loss and set-up due to breaking of random waves[C]//Proceeding 16thConference on Coastal Engineering.Hamburg:ASCE,1978:569-587.