基于缓坡方程在岛礁地形上波浪破碎的模拟研究

方亚冰，柳淑学，李金宣，刘思

（大连理工大学海岸和近海工程国家重点实验室，大连116024）

基于缓坡方程在岛礁地形上波浪破碎的模拟研究

方亚冰，柳淑学，李金宣，刘思

（大连理工大学海岸和近海工程国家重点实验室，大连116024）

波浪从深海传至近岸岛礁时，一般需要经过落差较大的礁坪边缘，水深急剧变化，导致波浪在传播过程中发生破碎，因此准确模拟波浪的破碎过程以及破碎后的波高大小，对于岛礁海岸工程建设具有重要的意义。缓坡方程是描述近岸波浪传播变形较好的数值模型之一，文章在采用自适应有限元求解缓坡方程所建立的数值计算模型的基础上，引入描述波浪破碎的模型，建立可以描述波浪破碎影响的近岸波浪数值模型。基于二维岛礁地形上的波浪实验，比较分析了4种不同的波浪破碎能量损失因子，给出了适合于岛礁地形条件下波浪传播破碎模拟的模型。

二维岛礁地形；缓坡方程；波浪破碎

中国海域幅员辽阔、珊瑚岛礁众多，波浪从深海传至近岸岛礁时，一般需经过落差较大的礁坪边缘。而由于水深的急剧变化，导致与在大多数缓变地形上波浪传播相比，其规律有所不同。因此建立起能够模拟波浪在岛礁地形上传播的数值计算模型，对波浪变形进行较为准确的预报是十分有必要的。

考虑到在波浪传播的过程中，会受到折射、绕射、反射等因素的影响，很多学者在二维椭圆型缓坡方程的基础上，建立起了多个数值模型（例如：Mei［1］；Tsay and Liu［2］；Chen and Houston［3］；Mattioli［4］），并运用到了实际工程当中。但是，大部分模型都没有考虑到如岛礁地形这种水深差别较大情况下波浪破碎的影响。本文在Berkhoff［5］发展的缓坡方程的基础上，基于自适应有限元方法来求解缓坡方程［6］，着重考虑波浪破碎因素的影响，针对不同学者提出的4种不同的波浪破碎能量损失因子，对二维岛礁地形条件下波浪的传播进行模拟，通过与实验结果进行对比、分析，建立起适用于岛礁地形下波浪传播的数值计算模型。

1 数值计算模型

1.1控制方程

考虑到波浪在传播过程中波浪破碎的效应，有些学者认为可以把波浪破碎作为一项耗散项加入到Berk⁃hoff的原始缓坡方程当中［7-9］。

式中

其中C（x,y）=ω/k为相速度，Cg(x,y)=为群速度，γ为波浪破碎能量损失因子,可有多种计算模型，将在下文给出。k为波数，与波浪频率满足如下色散关系

式中：h为水深。

1.2边界条件

在近岸波浪传播模拟过程中，一般需考虑如下边界条件。

（1）全反射边界。

（2）部分反射边界。

式中：α=α1+iα2为复常数，与复反射系数有关，复反射系数可表示为ρ=Reiε；ε为入射波和反射波之间相位差；R为振幅衰减；α与ρ的关系为

（3）入射边界条件。

对于一般波浪传播模拟，可认为波浪是从一直线上开始向计算区域内入射传播的，而在这条直线上有入射势（ΦI=a）和散射势ΦS，且不考虑入射边界对于散射势的反射，因此有

1.3方程求解

假设将要求解的区域划分为有限个单元，取其中一个具有代表性的单元，假定权函数为N，采用加权余量法、分部积分和格林公式等，可以建立求解单元上的有限元方程，最后将单元有限元方程在全区域上进行叠加，即可得到如下的有限元方程组

其中[] K具有对称和稀疏的特点，计算中采用索引存储法，只储存上三角的非零元素。方程的求解则采用双共轭梯度法，模型的建立及具体求解过程参见Liu S X等［6］。

2 波浪破碎能量损失因子的选取及其模拟结果

本文选取了Battjes和Janssen［10］，Dally et al.［11］，Massel［12］，Chawla et al.［13］等建立的4种不同的波浪破碎能量损失因子进行数值计算模拟。这些破碎因子参数都依赖于波高而变化。随着地形的变化，波高将产生明显的变化，当波高超过了一定的限制条件，波浪发生破碎，能量发生损失，直到产生适于当地水深条件的波浪。

本文将基于岛礁地形波浪传播的实验结果，针对上述常用的4种破碎模式进行研究，确定可用于岛礁地形条件下波浪破碎模拟的模型。

2.1物理实验模型

为了模拟岛礁地形条件下波浪的传播特性，刘宁［14］进行了相应的物理模型试验研究，实验布置如图1，水槽长69 m，宽2 m，深1.8 m。实验将岛礁地形简化为1：5的斜坡，斜坡段水平距离2.5 m，礁坪部分简化为水平地形，表面为光滑混凝土面，高度0.5 m，水平距离30 m，水槽后端布置斜坡式消能器用于吸收波浪。

模型中布置了18根浪高仪，具体位置及间距详见图2。

图2试验浪高仪布置图（单位：cm）Fig.2Sketch of wave height gauge

实验中坡前水深分别为h1=0.625 m、0.715 m和0.835 m，为研究入射波浪大小对于波浪传播破碎的影响，在特定周期T条件下，逐渐增大入射波浪波高H0进行实验，具体实验参数和过程可参考刘宁［14］的论文。

2.2波浪破碎能量损失因子的数值对比分析

2.2.1Battjes和Jassen［10］参数（BJ）

Battjes和Jassen［10］建议波浪破碎能量损失因子可由下式计算

式中：α为可调常数，一般可以取α=1，

在式（13）中，H为当地波高；Hm为最大允许波高，定义如下

式中：γ0为波浪破碎参数。对于浅水，式（14）则可以简化为Hm=γ0h（γ0取0.8，h表示当地水深）。

事实上，当式（13）中的b=0.3的时候，Qb=1.5×10-5，即Qb→0，γ→0。因此式（13）可以给出破碎波高条件

即当波高H≤Hb时，γ的值等于0（未破碎），否则γ将按照式（11）进行计算。

将式（11）代入式（1）针对上述实验条件下的波浪传播进行模拟计算，图3给出了3种坡前水深条件下、相同周期不同入射大小波浪传播的模拟结果与物理模拟实验结果的比较，图4分别给出了同一水深条件下相同入射波高H0、不同周期波浪传播的模拟结果与实验结果的比较，从图中可以看出：数值模拟给出的波浪破碎点位置与实验结果基本一致，发生在图2所示的10～14#浪高仪之间，即距造波边界距离为30.75～32.5 m；相同水深、相同周期不同入射波高情况下，入射波高越大，波浪破碎后稳定的波高稍小，与实验结果基本一致；但是相同水深、同一周期不同入射波高情况下，入射波高越大，波浪破碎后的稳定点位置相对于实验结果前移，而且稳定点之间的距离较大，与实验结果差别较大，其原因在于式（11）中的Qb是通过迭代求解的，迭代过程甚至达到上千次，迭代过程中会造成破碎后收敛的不一致性，对于落差较大的岛礁地形适应性较差。因此，总体来讲，虽然该模型可以给出与实验数据比较一致的结果，但是该方法不是模拟岛礁地形条件下破碎波浪的最佳选择。

2.2.2Dally et al.［11］参数（DDD）

以Horikawa和Kuo［15］的实验数据为基础，Dallyet al.［11］建议波浪破碎能量损失因子由下式计算

图3三种水深下相同周期不同入射波高H0情况下基于BJ破碎因子模型模拟波浪结果与实验结果的比较Fig.3Comparison of the simulated results based on BJ wave breaking factor and experimental ones for the same periods and different incident wave heightH0under three different depths

式中：H和h分别为当地波高和水深，参数Γ和χ分别表示稳定波浪因子和波浪延迟因子。刘宁［14］亦根据实验波浪结果给出了不同情况下的Γ和χ值（表1）,其所得结果与Dally et al.［11］给出相同实验条件下的结果是一致的，因此采用刘宁［14］所给出的Γ和χ值代入式（16）进行计算。

同样，式（16）中也包含了一个最低破碎限制条件，即

因此，在实际的计算过程中，如果在之前的迭代过程中所得到的波高H≤Hb，那么γ的值将取0，除此之外γ将按照式（16）进行计算。

与前述类似，图5给出了采用式（16）所定义的破碎因子，针对3种坡前水深条件下、相同周期不同入射波高波浪传播的模拟结果与物理模拟实验结果的比较，而图6分别给出了同一水深条件下相同入射波高、不同周期波浪传播的模拟结果与实验结果的比较。从图中可以看出：与前述Battjes和Jassen［10］模型模拟结果类似，数值模拟给出的波浪破碎点位置与实验结果基本一致；相同水深、相同周期不同入射波高情况下，入射波高越大，同样模拟所得波浪破碎后稳定的波高略微减小，同时，入射波高越大，波浪破碎后的稳定点位置前移，与实验结果一致。而且，与实验结果类似，相同水深、相同入射波高但不同周期情况下，波浪破碎后稳定的波高变化不大。因此，综合不同工况情况下的计算结果，式（16）所定义的波浪破碎因子可以较好地描述岛礁地形条件下破碎波浪及波浪传播过程。

图4h1=0.715 m水深下入射波高H0=0.14 m不同周期情况下基于BJ破碎因子模型模拟波浪结果与实验结果的比较Fig.4Comparison of the simulated results based on BJ wave breaking factor and experimental ones for the incident wave heightH0=0.14 mand different periods under the depthh1=0.715 m

表1不同水深比下规则波试验所得Γ和χ参数取值Tab.1TheΓandχvalues for regular waves with different ε

图5三种水深下相同周期不同入射波高H0情况下基于DDD破碎因子模型模拟波浪结果与实验结果的比较Fig.5Comparison of the simulated results based on DDD wave breaking factor and experimental ones for the same periods and different incident wave heightH0under three different depths

2.2.3Massel［12］的参数（SD）

Massel［12］给出波浪破碎能量损失因子的计算模式为

但是在实际计算过程中，发现当H＞2.85h的时候，γ＜0，因此在计算γ时，当地波高H应包含一个最高限制条件，即

另一方面，在数值计算过程中，同样发现当波浪发生破碎后，波高会随着波浪的传播，逐渐递减至0，这不符合实际波浪的传播情况，其原因是，不同于前述2种破碎波浪因子，式（18）未给出最低限制条件，导致该式不能给出有效的结果，因此为了更好地模拟实际波浪传播，仍取式（17）作为采用式（18）进行波浪模拟时的限制条件。

与上述结果类似，图7和图8分别给出了采用式（18）所定义的波浪破碎因子数值模拟的结果与实验结果的对比，从图中可以看出，采用该模型数值模拟所得波浪破碎点所在位置普遍比实验结果前移，即波浪较实验波浪提前破碎，但破碎波高大致相同；而在相同水深、相同周期和不同入射波高情况下，入射波高越大，波浪破碎后模拟所得稳定的波高比实验结果偏大，同时入射波高增大时，波浪破碎后的稳定点位置没有明显的前移，甚至还有退后的趋势，与实验结果差别较大。另外，相同水深、相同入射波高不同周期情况下，随着周期的减小，波浪破碎后稳定的波高增大，且增大幅度比较大，与实验结果差别也较大。这说明式（18）描述的破碎因子不能较好地模拟岛礁地形条件下波浪破碎的情况。

图6h1=0.715 m水深下入射波高H0=0.14 m不同周期情况下基于DDD破碎因子模型模拟波浪结果与实验结果的比较Fig.6Comparison of the simulated results based on DDD wave breaking factor and experimental ones for the incident wave heightH0=0.14 mand different periods under the depthh1=0.715 m

图7三种水深下同周期不同入射波高H0情况下基于SD破碎因子模型模拟波浪结果与实验结果的比较Fig.7Comparison of the simulated results based on SD wave breaking factor and experimental ones for the same periods and different incident wave heightH0under three different depths

2.2.4Chawla et al.［13］的参数（COK）

Chawla et al.［13］依据Thornton和Guza［16］所统计的波浪破碎数据，结合基于方程（1）式的椭圆形近似模型，归纳推导出了下面的破碎因子计算公式：

式中:λ和B建议分别取1.0和0.6。该式与式（18）类似，缺少一个限制条件（即γ=0时Hb的取值），导致波高出现不合理的结果，所以仍然采取式（17）作为相应的限制条件。

同样，图9和图10分别给出了数值计算结果与实验结果的比较，从图中可以看出：受所引入的破碎因子影响，波浪破碎之前的波高波动比上述模型（BJ、DDD、SD）给出的结果大，而且波浪破碎点的位置与实验结果亦有很大的差别，而且随入射波高增大，数值计算所得波浪破碎后的稳定波高稍小于实验结果，而且波浪破碎后的稳定点位置前移，幅度较大，与实验结果差别较大。因此，式（20）也不能较好地描述岛礁地形条件下波浪的传播情况。

图8h1=0.715 m水深下入射波高H0=0.14 m不同周期情况下基于SD破碎因子模型模拟波浪结果与实验结果的比较Fig.8Comparison of the simulated results based on SD wave breaking factor and experimental ones for the incident wave heightH0=0.14 m and different periods under the depthh1=0.715 m

2.2.5波浪破碎能量损失因子的对比

为进一步清楚的比较上述4种波浪破碎能量损失因子对岛礁地形条件下波浪传播的适应性，图11给出了这4种波浪破碎能量损失因子在3种坡前水深条件下、相同周期、相同入射大小波浪传播的模拟结果的比较。由图可以看出，如前所述，BJ破碎因子模拟波浪破碎后稳定点的位置偏后，且稳定波高明显偏高；COK破碎因子虽然模拟所得稳定波高与试验结果差别不大，但是亦存在波浪破碎后稳定点位置偏后的问题；SD破碎因子有时波浪破碎后模拟所得稳定波高偏大；综合比较，DDD破碎因子各方面模拟结果均与实验结果一致，可以较好地描述岛礁地形条件下的破碎波浪。

图9三种水深下同周期不同入射波高H0情况下基于COK破碎因子模型模拟波浪结果与实验结果的比较Fig.9Comparison of the simulated results based on COK wave breaking factor and experimental ones for the same periods and different incident wave heightH0under three different depths

图10h1=0.715 m水深下入射波高H0=0.14 m入射波高不同周期情况下基于COK破碎因子模型模拟波浪结果与实验结果的比较Fig.10Comparison of the simulated results based on COK wave breaking factor and experimental ones for the incident wave heightH0=0.14 mand different periods under the depthh1=0.715 m

图11三种水深下同周期同入射波高H0情况下四种破碎因子的比较Fig.11Comparison of the four wave breaking factors for the same periods and incident wave heightH0under three different depths

3 结语

波浪在向岛礁传播时，水深急剧变化，与其他缓变地形上的波浪传播运动不同。准确地模拟岛礁地形条件下波浪的破碎，对于岛礁地形上波浪传播的模拟具有重要的意义，本文在采用自适应有限元方法来求解缓坡方程的基础上，着重考虑波浪破碎因子对于波浪传播模拟结果的影响，建立起基于缓坡方程适用于求解岛礁地形下波浪传播的数值计算模型。

本文选用了Battjes和Janssen［10］，Dally et al.［11］，Massel［12］，Chawla et al.［13］等人提出的不同的波浪破碎能量损失因子分别建立数值计算模型，针对二维岛礁地形条件下波浪的传播进行模拟。为了考虑入射波高大小对于波浪破碎的影响，在模拟过程中，逐渐增大入射波高，对于波浪的传播进行计算。通过计算结果与实验结果进行对比分析，综合考虑波浪破碎点的位置和波高、波浪破碎后稳定点的位置以及稳定后的波高大小等因素，对比分析结果表明，由Dally et al.［11］建议的式（16）作为描述波浪破碎能量损失因子，结合通用的缓坡方程所建立的数值计算模型可以较好地描述典型岛礁地形条件下波浪的传播计算。

参考文献：

［1］Mei C C.The Applied Dynamics of Ocean Surface Waves［M］.New York:Wiley,1983.

［2］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.

［3］Chen H S,Houston J R.Calculation of Water Oscillations in Coastal Harbors;HARBS and HARBD User′s Manual［R］.Vicks⁃burg MS:Coastal Engineering Research Center,1987.

［4］Mattioli F.Dynamic response of the Lido channel to wavemotion in the presence of movable barriers［J］.IlNuovoCimento C,1996, 19（1）:177-194.

［5］Berkhoff J C W.Mathematical models for simple harmonic liner water waves⁃wave diffraction and refraction［R］.Delft,the Nether⁃lands:Delft Hydraulic Lab.,1976.

［6］LIU S X,SUN B,SUN Z B,et al.Self⁃adaptive FEM numerical modelling of the mild⁃slope equation［J］.Applied Mathematical Modelling,2008,32:2 775-2 791.

［7］Liuzhi Zhao,Vijay Panchang,Chen W,et al.Simulation of wave breaking effects in two⁃dimensional elliptic harbor wave models［J］.Coastal Engineering,2001,42:359-373.

［8］唐军,沈永明,郑永红，等.结合椭圆型缓坡方程模拟近岸波流场［J］.海洋学报,2006,28（1）:146-151. TANG J,SHEN Y M,ZHENG Y H,et al.The numerical simulation of nearshore current combined with the elliptic mild⁃slope equation［J］.Acta Oceanologica Sinica,2006,28（1）:146-151.

［9］祁峰,曹宏生.应用改进的椭圆型缓坡方程数值模拟波浪传播变形［C］//中国海洋工程学会.第十三届中国海洋（岸）工程学术讨论会论文集.北京：海洋出版社，2007:316-318.

［10］Battjes J A,Janssen J.Energy loss and set⁃up due to breaking of random waves［J］.Coastal Engineering Proceedings,1978,1（16）:569-587.

［11］Dally W R,Dean R G,Dalrymple R A.Wave height variation across beaches of arbitrary profile［J］.Journal of Geophysical Re⁃search:Oceans（1978-2012）,1985,90（C6）:11 917-11 927.

［12］Massel S R.Inclusion of wave⁃breaking mechanism in a modified mild⁃slope model［M］//Banner M L,Grimshaw R H J.Breaking Waves.Springer⁃Verlag Berlin and Heidelberg GmbH&Co.K,1992:319-324.

［13］Chawla A,Özkan⁃Haller H T,Kirby J T.Spectral model for wave transformation and breaking over irregular bathymetry［J］.Jour⁃nal of waterway,port,coastal,and ocean engineering,1998,124（4）:189-198.

［14］刘宁.波浪在岛礁地形上传播特性的试验研究［D］.大连:大连理工大学,2014.

［15］Horikawa K,Kuo C T.A study on wave transformation inside surf zone［J］.Coastal Engineering Proceedings,1966,1（10）:217-233.

［16］Thornton E B,Guza R T.Transformation of wave height distribution［J］.J.Geophys.Res,1983,88（10）:5 925-5 938.

Simulation of wave breaking on the reefs terrain using mild⁃slope equation

FANG Ya⁃bing,LIU Shu⁃xue,LI Jin⁃xuan,LIU Si
（State Key Laboratory of Coastal and Offshore Engineering,Dalian University of Technology,Dalian 116024,China）

When the waves propagate from the deep water to the inshore reefs,a large gap reef flat must be overcome,leading to dramatic changes in water depth.This will cause the wave breaking.It is important for the engi⁃neering design on the terrain to accurately simulate the wave breaking and calculate the wave heights after wave breaking.Mild⁃slope equation is one of the best numerical models to describe near shore wave transformation. Based on a self⁃adaptive finite element numerical model to solve the mild⁃slope equation,a near shore numerical model to simulate wave breaking was developed by introducing the wave breaking model in this paper.Four kinds of wave breaking model were used and compared based on the experimental results of the wave propagation on the two⁃dimensional reefs terrain.A numerical model for solving the wave transformation on the two⁃dimensional reefs ter⁃rain was proposed.

two⁃dimensional reefs terrain;mild⁃slope equation;wave breaking

TV142；TV131.6

A

1005-8443（2015）04-0290-07

2014-12-05；

2015-01-04

国家重点基础研究发展（973）计划资助项目（2013CB036101，2011CB013703）；国家自然科学基金创新研究群体基金（51221961）

方亚冰（1989-），男，湖南省郴州市桂东县人，硕士研究生，主要从事近岸波浪传播特性的模拟研究。

Biography：FANG Ya⁃bing(1989-),male,master student.