基于改进缓坡方程的波浪传播数值模拟

王红川，周正萍

(南京水利科学研究院 水文水资源与水利工程国家重点实验室，江苏 南京 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］给出: