修正型缓坡方程的有限元模型

倪云林，滕斌，丛龙飞

(1.大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116024； 2. 浙江海洋大学 港航与交通运输工程学院，浙江 舟山 316022)



修正型缓坡方程的有限元模型

倪云林1,2，滕斌1*，丛龙飞1

(1.大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116024； 2. 浙江海洋大学 港航与交通运输工程学院，浙江 舟山 316022)

与缓坡方程相比，修正型缓坡方程增加了地形曲率项和坡度平方项，从而提高了数值求解的复杂性。本文将计算域划分为内域和外域，内域为水深变化区域，使用修正型缓坡方程，其中的地形曲率项和坡度平方项可用有限单元各节点的水深信息和单元插值函数表示，外域为水深恒定区，速度势满足Helmholtz方程，通过内外域的边界匹配建立有限元方程，并用高斯消去法求解。进而分别模拟了波浪传过Homma岛和圆形浅滩的变形，其结果与相关的解析解和实验数据吻合良好，证明了本文有限元模型的正确性。同时，通过与实验数据的对比也明显看出，在地形坡度较陡的情况下，修正型缓坡方程较缓坡方程具有更高的计算精度。

修正型缓坡方程；有限元；Homma岛；圆形浅滩

1 引言

波浪在由深海向浅水传播的过程中，受海底地形的影响，会发生浅化、折射、绕射等一系列的变形。为了研究这一过程中波浪要素的变化规律，许多学者建立和改进了各种数学模型。缓坡方程(mild slope equation, MSE)就是其中的一种，通过沿水深积分的方法将三维问题转化为二维问题，具有方程形式简单且应用范围较广的优点。

缓坡方程最初由Berkhoff(1972年)提出，又称为折射-绕射联合模型，是在缓坡假设的基础上，用线性波浪理论研究缓变、不可渗透地形上波浪的传播变形问题[1]。关于该方程的使用条件，Booij(1983年)证明缓坡方程在底坡小于1∶3的情况下具有足够的精度[2]，但在底坡大于1∶3或地形波动的情况下，缓坡方程则不能准确描述波浪的传播变形。为了改善“缓坡假设”的条件限制，Kirby(1986年)最先对缓坡方程进行了改进，他将剧变地形分解成一个缓变地形和一个小振幅波动地形的叠加，导出了扩展型缓坡方程(extended mild slope equation，EMSE)，提高了地形波动情况下缓坡方程的计算精度[3]。之后，Chamberlain和Porter[4](1995年)，Chandrasekera和Cheung[5](1997年)通过引入地形曲率项和坡度平方项对缓坡方程进行了修正，得到了修正型缓坡方程(modified mild slope equation, MMSE)，提高了陡变地形和剧变地形情况下缓坡方程的计算精度，从而大大拓宽了缓坡方程的应用范围。

目前，有不少文献对修正型缓坡方程的求解方法进行了论述。例如，Suh等将修正型缓坡方程转化为一阶双曲型方程组，采用有限差分法进行求解[6]。Silva等采用中心差分格式离散控制方程和边界条件来求解修正型缓坡方程[7]。有限差分法具有数学概念直观、表达简单、发展成熟等优点，但差分网格通常为矩形，在边界不规则或形状复杂时精度降低。因此，本文将采用有限元法直接对修正型缓坡方程进行求解，解决地形曲率项和坡度平方项的空间离散问题。在求解过程中，将整个计算域划分为内域和外域，内域是水深变化的有限区域，可用有限元离散，外域是水深恒定的无限区域，速度势满足Helmholtz方程，其解为级数形式，通过内外域公共边界上的匹配条件建立并求解有限元方程[8—11]。

2 控制方程

针对水深变化海床上的波浪运动问题，Chamberlain和Porter[4]采用变分原理推导了修正型缓坡方程：

▽·(ccg▽φ)+[k2ccg+f1(kh)g▽2h

+f2(kh)gk(▽h)2]φ=0.

(1)

式中，▽为水平梯度算子；φ为速度势的空间分量；c和cg分别为波速和群速度；k为波数；h为水深；g为重力加速度；f1(kh)和f2(kh)分别为地形曲率项系数和坡度平方项系数，其表达式如式(2)、(3)所示[5]。

f1(kh)=[-4khcosh(kh)+sinh(3kh)+sinh(kh)

+8(kh)2sinh(kh)]/{8cosh3(kh)[2kh+sinh(2kh)]}

(2)

16(kh)3sinh(2kh)-9sinh2(2kh)cosh(2kh)+

12(kh)[1+2sinh4(kh)][kh+sinh(2kh)]}.

(3)

上述两个系数均为kh的函数。从图1可以看出，f1(kh)和f2(kh)在有限水深区数值较大，而随着水深的增大，两者均趋于0。

图1 f1和f2随相对水深kh的变化Fig.1 Variation of coefficients f1 and f2 with kh

当忽略地形曲率项和坡度平方项时，修正型缓坡方程就简化为Berkhoff的缓坡方程：

▽·(ccg▽φ)+k2ccgφ=0.

(4)

当水深不变时，则转换为Helmholtz方程：

▽2φ+k2φ=0.

(5)

3 有限元方程的建立

3.1 计算域的划分

如图2所示，将计算域划分为内域和外域。其中，内域Ω1是地形变化区域，划分六节点三角形单元网格，用单元插值函数作为权函数以建立有限元方程；外域Ω2是水深恒定区域，不划分网格，用Hankel函数作为权函数，并应用格林公式将外域的面积分转化为边界Ba上的线积分以建立方程式；Ba为内外域的公共边界，通常情况下取为圆周，通过Ba上的匹配条件可进行有限元方程的求解；Bb为物面边界；B∞为无穷远边界。

图2 计算域划分Fig.2 Division of computational domain

3.2 边界条件

设内域Ω1的速度势为φ，则在物面边界Bb上满足

(6)

(7)

(8)

式中，r为Ω2域内x-y平面上某点到原点O的距离；φ为r与x轴正方向的夹角；Hn(kr)为n阶第1类Hankel函数；An(n=0,1,2，…)和Bn(n=1,2，…)为待定系数。在n=N处截断，并令G1=H0(kr)，G2i=Hi(kr)cosiφ，G2i+1=Hi(kr)siniφ，α1=A0，α2i=Ai，α2i+1=Bi(i=1,2,…，N)，式(8)可简写成：

(9)

式中，M=2N+1。则无穷远处辐射边界条件可以表示为:

(10)

在公共边界Ba上满足：

(11)

3.3 内域有限元方程

设φj为单元各节点速度势，Nj为单元插值函数，则单元内的速度势：

(12)

利用单元各节点的水深信息和单元插值函数，坡度平方项可表示为:

(13)

(14)

(15)

类似的，地形曲率项

(16)

取单元插值函数Ni为权函数，利用Galerkin方法在单元内建立积分方程

[▽·(ccg▽φ)+(k2ccg+f1(kh)g▽2h

+f2(kh)gk(▽h)2)φ]NidΩ=0.

(17)

引入三节点线单元插值函数Li，对边界Ba上的单元而言，Li为该单元外侧边界上的单元插值函数。将边界Ba上的速度势分解为入射势φI和绕射势φS，应用格林公式，式(17)可以重新写成:

(18)

结合式(12)、(15)和(16)，可以写成单元矩阵形式

(19)

其中,

集合内域各单元矩阵可得:

(20)

式中，q为总节点数；l为边界Ba上的节点数。

3.4 外域方程式[12—13]

(21)

应用第二格林公式，式(21)可以写成:

(22)

(23)

(24)

3.5 总矩阵方程

根据式(20)和式(24)，可以形成总矩阵方程

[K](q+M)×(q+M){φ}(q+M)×1={B}(q+M)×1.

(25)

式中，

；

由于Ba为圆周，所以[K22]可用解析表达式表示

[K22]=2πkRccgdiag{H0′(kR)H0(kR),

H1′(kR)H1(kR),…,HM′(kR)HM(kR)},

(26)

式中,R为圆周Ba的半径。

最后，本文采用列主元高斯消去法求解总矩阵方程式(25)。

4 结果验证和讨论

为了验证本文修正型缓坡方程的有限元模型，首先计算了长波经过Homma岛的绕射问题，计算结果与Zhai等[14]的解析解进行了比较；接着计算了圆形浅滩上波浪的绕射问题，计算结果与Suh等[6]的实验数据进行了比较，并进一步探讨了较陡地形对不同周期波浪的影响。

4.1 长波过Homma岛的绕射问题

Homma岛由旋转抛物体及上部圆柱体组合而成。如图3所示，上部圆柱体的半径a=10 km，旋转抛物体的底面半径b=30 km。距岛屿中心r(单位：km)处的水深(单位：km)为:

(27)

图3 Homma岛地形Fig.3 Sketch of Homma Island

图4 岸线圆周上修正型缓坡方程计算的相对波高和解析解的对比Fig.4 Comparison of calculated wave height by MMSE with analytical solutions along the shoreline

4.2 波浪过圆形浅滩的绕射问题

为了进一步验证较陡地形上修正型缓坡方程有限元模型的正确性和应用范围，本文将修正型缓坡方程有限元解、缓坡方程的有限元解与Sub等[6]的实验结果进行比较。实验设计的圆形浅滩地形如图5所示，距浅滩中心r(0≤r≤R，单位：m)处的水深(单位：m)为:

(28)

图5 圆形浅滩地形Fig.5 Sketch of circular shoal

图6是修正型缓坡方程和缓坡方程计算的y=0断面上相对波高与实验结果的对比。从图中可以看出，相比缓坡方程，修正型缓坡方程的计算结果与实验数据更加吻合，这说明在地形较陡的情况下，地形曲率项和坡度平方项对波浪的传播变形会起到十分重要的作用。通过比较图6a～c，并结合表1可以发现，在相对水深较大时，波高最大值发生在滩顶后方，而随着相对水深的减小，波高的最大值有所增加，发生最大值的位置向前方移动，缓坡方程的计算结果在滩顶后方也表现出更大的误差，这说明长波较短波更易受到地形曲率和坡度平方的影响。

5 结论

(1)本文建立了修正型缓坡方程的有限元模型，并应用该模型计算了波浪过Homma岛和圆形浅滩的绕射问题，计算结果与翟西媛等的解析解和Sub等的实验结果吻合良好，证明了本文有限元模型的正确性。

图6 y=0断面修正型缓坡方程和缓坡方程计算结果与实验结果的对比Fig.6 Comparison of MMSE and MSE results with experimental data along y=0 section

波浪参数最大相对波高Hmax/H0最大相对波高发生位置x/R实验结果MMSE计算结果MSE计算结果实验结果MMSE计算结果MSE计算结果T=1.259s,kh0=0.9651.2621.2561.239-0.220-0.178-0.022T=0.791s,kh0=2.0031.2521.2351.3130.6670.7110.800T=0.636s,kh0=3.0031.1971.1821.2311.3251.0000.889

(2)通过本文修正型缓坡方程计算结果、缓坡方程计算结果和实验数据的三者对比，可以看出，对于较陡地形上的波浪变形问题，修正型缓坡方程较缓坡方程具有更高的计算精度，这从侧面证明了本文的有限元模型具有更宽的适用范围。

(3)针对长波入射的情况，Homma岛上部圆柱周围的最大波高不是发生在迎浪点，而是对称分布于迎浪点左右两侧，这与直立圆柱的情形不同。

(4)在地形坡度较陡的情况下，地形曲率项和坡度平方项会对波浪的传播变形起到不可忽略的影响，并且随着波浪周期的增大，波浪受地形曲率和坡度平方的影响也随之增大。

[1] Berkhoff J C W. Computation of Combined Refraction-diffraction[M]. Delft Hydraulics Laboratory, 1974.

[2] Booij N. A note on the accuracy of the mild-slope equation[J]. Coastal Engineering, 1983, 7(3): 191-203.

[3] Kirby J T. A general wave equation for waves over rippled beds[J]. Journal of Fluid Mechanics, 1986, 162: 171-186.

[4] Chamberlain P G, Porter D. The modified mild-slope equation[J]. Journal of Fluid Mechanics, 1995, 291: 393-407.

[5] Chandrasekera C N, Cheung K F. Extended linear refraction-diffraction model[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 1997, 123(5): 280-286.

[6] Suh K D, Lee C, Park Y H, et al. Experimental verification of horizontal two-dimensional modified mild-slope equation model[J]. Coastal Engineering, 2001, 44(1): 1-12.

[7] Silva R, Borthwick A G L, Taylor R E. Numerical implementation of the harmonic modified mild-slope equation[J]. Coastal engineering, 2005, 52(5): 391-407.

[8] 李孟国, 蒋德才. 关于波浪缓坡方程的研究[J]. 海洋通报, 1999, 18(4): 70-92.

Li Mengguo, Jiang Decai. A review on the study of mild-slope equation[J]. Marine Science Bulletin, 1999,18(4): 70-92.

[9] Chen H S, Mei C C. Oscillations and wave forces in a man-made harbor in the open sea[C]//Symposium on Naval Hydrodynamics, 10th, Proceedings, Pap and Discuss, Cambridge, Mass, June 24-28, 1974. 1976.

[10] Tsay T K, Liu P L F. A finite element model for wave refraction and diffraction[J]. Applied Ocean Research, 1983, 5(1): 30-37.

[11] Houston J R. Combined refraction and diffraction of short waves using the finite element method[J]. Applied Ocean Research, 1981, 3(4): 163-170.

[12] 赵明. 波浪作用下建筑物周围的泥沙冲刷及海床演变[D]. 大连: 大连理工大学, 2002.

Zhao Ming. The local scour and topographical change around offshore structures under wave action[D]. Dalian: Dalian University of Technology, 2002.

[13] 赵明, 滕斌. 缓坡方程的有限元解[J]. 大连理工大学学报, 2000, 40(1): 117-119.

Zhao Ming, Teng Bin. Finite element solutions for mild slope equation[J]. Journal of Dalian University of Technology, 2000, 40(1): 117-119.

[14] Zhai X Y, Liu H W, Xie J J. Analytic study to wave scattering by a general Homma island using the explicit modified mild-slope equation[J]. Applied Ocean Research, 2013, 43: 175-183.

[15] Chau F P, Taylor R E. Second-order wave diffraction by a vertical cylinder[J]. Journal of Fluid Mechanics, 1992, 240(1): 571-599.

FEM model of the modified mild slope equation

Ni Yunlin1,2，Teng Bin1，Cong Longfei1

(1.StateKeyLaboratoryofCoastalandOffshoreEngineering,DalianUniversityofTechnology,Dalian116024,China；2.SchoolofPortandTransportationEngineering,ZhejiangOceanUniversity,Zhoushan316022,China)

Compared with the mild slope equation, the bottom curvature and slope-squared terms are contained in the modified mild slope equation, which increases the complexity of solving the equation numerically. In this paper, the physical domain is divided into a finite inner region and an infinite outer region. The inner region of a variable depth is studied with the modified mild slope equation, and the outer region has a constant water depth, the velocity potential satisfying the Helmholtz equation. The bottom curvature and slope-squared terms in the equation are evaluated from the input bathymetry at each node of every element using the interpolation functions. With the boundary matching method and Gaussian elimination technique, the finite element equations are established and solved. Then, wave transformation over a Homma Island and a circular shoal is simulated respectively, and the results agree well with analytical solutions and experimental data, verifying the validity of the FEM model in this paper. Meanwhile, the comparison between the numerical and experimental results indicates the modified mild slope equation can provide more accurate predictions than the mild slope equation for wave propagation over relatively steep bathymetry.

the modified mild slope equation; finite element method; Homma Island; circular shoal

10.3969/j.issn.0253-4193.2017.01.011

2016-03-07；

2016-05-17。

国家自然科学基金(51379032，51490672)。

倪云林(1986—)，男，浙江省舟山市人，讲师，博士生，主要从事波浪对海上建筑物作用的研究工作。E-mail:nylzjou@126.com

*通信作者：滕斌，教授，主要从事波浪对海上建筑物的作用研究。E-mail: bteng@dlut.edu.cn

P753

A

0253-4193(2017)01-0104-07

倪云林, 滕斌, 丛龙飞. 修正型缓坡方程的有限元模型[J]. 海洋学报, 2017, 39(1):104-110,

Ni Yunlin，Teng Bin，Cong Longfei. FEM model of the modified mild slope equation[J]. Haiyang Xuebao, 2017, 39(1): 104-110, doi: 10.3969/j.issn.0253-4193.2017.01.011