基于Boussinesq方程近岸波浪演变的数值研究

王培涛，于福江,2

（1.国家海洋环境预报中心，北京 100081；2.国家海洋局海洋灾害预报技术研究重点实验室，北京 100081）

基于Boussinesq方程近岸波浪演变的数值研究

王培涛1，于福江1,2

（1.国家海洋环境预报中心，北京 100081；2.国家海洋局海洋灾害预报技术研究重点实验室，北京 100081）

首先对目前描述近岸波浪传播变形的数学模型进行了回顾与总结；对不同数学模型的特点、适用范围和发展情况进行了阐述与对比。应用基于Boussinesq方程的Coulwave模式针对几个经典实验地形进行了数值实验，数值结果和实验实测数据吻合较好。此外，分别采用不同的近岸波浪模型模拟了某渔港附近波浪的传播变形，结果表明：当考虑波浪的折射、绕射、反射联合作用时，Coulwave模式计算结果明显较缓坡方程及SWAN模型计算结果更加合理。

Boussinesq型方程；近岸波浪模型；缓坡方程SWAN模型；Coulwave模型

1 引言

近岸水域是与人类生活发生联系最多的区域，大气和海洋每年总交换热量的近一半分布于近岸水域，如港口航道工程、近海取排水设施、海水养殖设施等，都集中分布在近岸海域内。在海岸动力场中，海浪一直是施加在近岸建筑物上的最重要的环境荷载之一。当波浪由外海深水区传播到海岸附近的浅水区时，受到水深、地形、底摩擦、障碍物（如建筑物、岛屿、岬角、沙洲、岸壁）等的影响，波浪的传播速度、波长、波高、波面状况、及其产生的流场和压力场都会发生较大的变化。对于复杂地形和复杂边界上波浪传播变形的研究方法主要包括现场观测、物理模型、数值模拟。然而伴随着计算机软件、硬件技术以及计算方法的快速发展，利用数值模型来模拟近岸波浪的传播变形变得简单有效，已逐渐成为近岸波浪研究的主要手段。

目前研究波浪在近岸区域传播变形的模型主要有：基于光学snell定律的射线理论（ray theory）、基于微幅波理论的缓坡方程（mild-slope equation）、基于波能和波作用量守恒方程的谱模型、基于holmheltz方程的波浪绕射模型、基于直接求解Navier-stokes方程模型以及考虑波浪的非线性和色散效应的Boussinesq方程模型[1-7]。然而，上述几种数学模型都依赖于特定的物理假设和数学近似，所以各自的适用范围就受到了限制。

射线理论本身为线性理论，适用于缓变地形，且假定无能量跨过波向线，在海底坡度较大，地形较复杂时应用射线理论会出现焦散现象和盲区现象（波向线不能进入的区域）[8-9],但射线理论具有简单省时省力的特点，计算精度为一阶精度；抛物型缓坡方程仅考虑沿波峰线方向的绕射作用，而忽略波向线方向的绕射作用，且没有考虑地形和障碍物的反射作用，在反射作用比较强的区域不适用。此外，也不适用于对于入射波向偏离主方向太大的情况；椭圆型缓坡方程变量为复变量，方程本身具有不可分离的性质，对它的离散求解不能应用迭代法，对于大的计算域求解比较困难；双曲型缓坡方程没有忽略波浪的反射特性，且计算精度与椭圆形缓坡方程相当，但其为时变型方程，为使模型收敛，需要用小的时间步长；谱模型用源函数的方式分别考虑了风应力作用、波浪破碎、底摩擦耗散、波-波相互作用等各种物理过程，可用于较大范围较长时间的计算，但由于其在波浪折射、绕射等方面的不足，对不规则地形下近岸波浪计算，特别是港池、航道等复杂地形下的波浪传播变形并不适用；基于holmheltz方程的波浪绕射模型是基于线性简谐波及等水深水域的波浪绕射方程，可用于任意形状港湾和防波堤掩护水域的波场计算[10-11]；基于不可压缩完整的三维Navier-stokes方程和连续方程的波浪数学模型可以描述全水深、任意坡度的线性和非线性的波浪传播变形，并可计算速度的垂向分布，但将完整的三维Navier-stokes方程应用于实际还有很多工作要做。

Boussinesq方程用流速和水位描述水体变化[12]，假定垂向速度沿水深为线性分布，由此得出水平速度和压力沿水深为抛物型分布，考虑了波浪的非线性及色散性特征，运动保持质量和动量守恒，因此能较好地模拟波浪传播过程中多种因素的相互作用，其适用性广泛。随着计算方法的逐渐成熟，Boussinesq方程已得到越来越广泛的应用与发展，特别是在拓展和改进模型的非线性及频散特征方面取得了较大的进步。本文将重点介绍基于Boussinesq方程的Coulwave模型对近岸波浪传播变形的模拟。

2 Boussinesq方程模型研究进展

Boussinesq（1872）[13]假定水平速度上下均匀，垂向速度由底到自由表面是线性分布，得出了一维非线性控制方程，称为Boussinesq方程。Peregrine（1967）[14]推导出了变水深条件下的二维Boussinesq方程，称为经典Boussinesq方程：

式中：x、y为与静止水面重合的直角坐标系坐标，u(x,y,t)、v（x,y,t）分别为水深平均的质点速度矢量沿x、y方向的分量，ζ为波面到静止水面的距离。

由于经典的Boussinesq方程包含非线性色散性，能够反映和描述近岸区域波浪的变形现象，被认为是短波数值模拟领域的一个重大突破。但Boussinesq方程具有弱色散性、弱非线性以及方程本身没有考虑底摩擦、波浪破碎和环境水流的影响，其适用性也受到了限制，适用水深范围h/L≤0.12。近20年来，Boussinesq方程的理论和应用都得到了快速发展，特别是在色散性及变浅作用、非线性、波浪破碎作用、环境流速等方面改进显著，从而拓展了其水深适用范围。

Madsen等采用在动量方程中加入含有1个待定系数常数Boussinesq三阶导数项的方法得到了1个常水深情况下的用单宽流量（或称水深积分速度）和波面表达的二维改进型Boussinesq方程，在进一步的工作中，Madsen等将其推广到适用缓变坡度变水深情况[15]。当B＝1/15时，改进型Boussinesq方程的线性色散性能达到Airy波精确解的Pade[2,2]阶近似。Madsen等同时提出了线性变浅梯度（linear shoaling gradient）概念作为衡量改进型Boussinesq方程变浅作用性能的指标。Madsen等的改进型Boussinesq方程的变浅作用性能的精度达到O(μ4)，能够适用于h/L≤0.5。在这一工作中，Madsen还详细讨论了方程色散性与方程速度变量选择的关系，指出了方程用不同的速度变量表达时，所得的色散关系是不同的，其中，用水深平均速度给出结果最好，而用自由表面速度结果最差。这一思想在Nwogu（1993）的研究中得到了很好的发展，他取任意水深处速度表达Boussinesq方程[16]，并通过与线性色散方程的Pade’，展开拟合，获得最优速度，所得方程色散关系与Witting的方程和Madsen的方程一致。Liu（1994）[17]和Wei等（1995）[18]扩展Nwogu方法于完全非线性波，导致模型不仅可以适用于中等水深，也适用于模拟强非线性相互作用。Gobbi（2000）[19]使用四阶多项式推导的模型，线性色散性精度到kh=6。Madsen等（2002）基于Agnon（1999）的方法提出了一个Boussinesq类模型，准确性达到深水（kh=40）。本方法采用Laplace方程的最优展开，通过使用水柱不同层的多级展开，只需要五阶空间导数就可以达到深水。Lynett[20-23]和Philip Liu提出了基于Boussinesq方程的波浪分层数学模型，假定在每一层中的垂向流速按线性分布，分布在各层对垂向流速进行垂向积分，求其平均流速，用这一平均流速来代替每一层的垂向流速。这样的处理既提高了方程的精度，又改进了方程的色散性，使方程能应用于更大的水深范围。图1是Boussinesq方程各阶段主要研究阶段色散关系精度对比情况。

图1 几种改进后的Boussinesq方程线性色散关系与Airy波色散关系精确解比较

3COULWAVE模型的数学物理基础及验证

Coulwave模型在垂直方向上将水柱分解成多层，采用深度积分模型来代替对垂直流场的高阶多项式近似，这种方法不仅提高了Boussinesq方程的计算精度，而且对方程的线性色散精度也进一步改善，更重要的是方程的物理阶数却没有提高，避免了数值求解高阶导数项带来的数值振荡问题。经拓展后的Coulwave模型在源函数造波、底摩擦效应、波浪破碎、海绵层滤波方面得到了完善，模型可以用来模拟水波浅化、折射、绕射和反射等综合现象，特别是在刻画结构物附近的波浪场，有其显著的优势。

3.1 COULWAVE模型的数学物理基础

本节将以两层模型为例进行介绍及应用。二层的Boussinesq方程模型在一层模型基础上，对水柱面再进行分层，在每一层中的垂向流速都按线性分布，分别对其进行沿水深方向积分求其平均流速，用此平均流速代替各层的流速（见图2）。

图2 二层模型示意图

对于两层模型来说，定义其水平流速向量为：

模型的连续方程及动量方程为：

其中u2可以用u1的函数关系式表达如：

底层流速可以通过上层流速计算求解。式（6）—（8）即为二层系统的联合控制方程。

3.2 COULWAVE模型验证

为了验证Coulwave模型的可用性，本文针对2个经典试验地形下波浪演变进行了数值研究，并将计算结果与实验数据进行了对比。

3.2.1波浪跨越潜堤的传播和变形

Dingemans（1994）的实验模型（见图3），左端x=0(h0=0.86 m)的位置设置造波机，取空间步长dx=0.10 m，时间步长Δt=0.0189 s，模拟时间120 s。在模拟中左右两端分别设置了1.25 L的海绵层。在图示的四个靶位放置测波仪测量波面时间序列。x1=26.04 m，x2=33.64 m，x3=37.04 m，x4=41.04 m

图3 Dingemans实验模型示意图

在实验中，入射波(波幅H=0.02 m，波周期T= 2.86，ε=a0/h0=0.023kh=0.7，L=7.7 m，h=0.86 m)传播经过＃1位置后，这时由于水深变浅，波浪发生浅水效应，使波浪变陡，同时由于非线性作用的影响，波浪的能量不断从主频传向高频，波浪中的高频谐波将不断增加。由于各谐波的相速度不同，以至于在堤后，高次谐波被释放出来，波浪的形状及相速度发生改变。我们分别用一层模型和两层模型对这次传播过程进行了数值模拟。从模拟结果图4来看：在坡顶（＃1位置处）出现了明显的非线性特征，表现为波峰变尖，波谷变平。当波浪到达坡顶附近时，波幅达到最大值，出现明显的次峰，次峰随着地形的变化而变化。从图4可以看到对于＃1位置，两种方法模拟结果与实验结果吻合良好。两层模型与一层模型模拟的效果相当，并没有显示出两层模型的优势所在。这是因为此时的kh＜2.0，一层、两层模式的色散性都能保证合理的精度。当波列行进到坡后，特别是＃3，＃4位，实验结果与模拟结果有了一些差异，这是由于当高次谐波释放后形成较短的波，μ=k0h0值较大，在坡后μ=kh≈4，超出了一层模型的预报精度，同时波浪在堤后的变形更加剧烈，这时对模型非线性和色散性精度要求更高，两层模型模拟的结果与实验数据更加吻合，说明采用两层处速度近似可以大大提高方程的色散性能（见图4）。

图4 波面高度历时变化与实验数据对比

3.2.2 波浪在Berkhoff实验地形上的传播演变

1982年Berkhoff等设计和进行了著名的椭圆缓坡地形波浪传播、变形的物理模型实验[24]。本实验可反映出波浪在浅水变形、折射、反射和非线性弥散关系等诸多因素影响下的变形和传播情况，已经成为数值模型计算验证的经典实验之一。为了验证本文模型在波浪浅水效应、反射、折射和非线性弥散等诸多因素共同作用下的适用性，本小节将分别用本文模型和抛物型缓坡方程对实验地形进行了数值模拟，并与实验数据进行了比较。

具体水下地形和实验测量断面分布见图5，计算区域为25 m×20 m，整个水域底坡为1∶50，斜坡梯度方向与波浪入射边界的法向夹角为20°。以x0=10m，y0=10m为中心，有一长半轴为4 m，短半轴为3 m的椭圆浅滩。斜坡旋转后的坐标为：

图5 Berkhoff地形实验设置

图6 计算所得T=10 s、30 s时刻三维波面

(x0，y0)位于浅滩中心，水深为0.1332 m。浅滩的边界为：

浅滩外水深如下：

浅滩上水深由下式给出

入射波波高为0.0464 m，波浪周期T＝1.0 s，入射方向从x正方向入射。计算空间步长为0.05 m，时间步长dt=0.01s，总时间步数为3000，即30个波周期。模拟启动的初始条件为“冷启动”，在x＝1.5 m处设置波源，左、右两端设置l＝1.5 m的海绵滤波层，上下两个侧边界为全反射边界（见图5）。

从模拟的结果（见图6）可以看出，当入射波前行遇到斜坡和椭圆形浅滩后开始变形，在斜坡作用下波浪发生折射，波峰线开始变得不再平行，一侧波长变短波高增加，一侧波长变长波高减小。当入射波传播至浅滩后方，由于两侧波浪绕过浅滩后的叠加，波能集中，波高显著增加。从断面比波高对比来看（见图7），缓坡方程及改进的Boussinesq方程模型对于波浪在缓坡地形上的传播变形都能得到较满意的结果；但对于断面6及断面8后半段而言，本文的模型跟实验数据吻合更好，究其原因可能为斜坡和浅滩引起波浪复杂的波-波相互作用所致。

4COULWAVE模型的工程应用对比

上一节所述的几个例子都是在边界相对简单，地形不太复杂的情况下进行的，实际的海底地形和岸边界情况往往比这复杂的多，计算区域也大得多。本节将利用本文模式和缓坡模式、SWAN模式对某渔港（见图8）附近区域波浪传播进行模拟对比。同时也将本文模型模拟波浪在渔港内的演变（见图9），以及利用上述三种近岸模式模拟外海波浪在特定工程点处的统计波高情况进行了研究,选取经常侵袭港口水域E（X=0）方向的入射波进行模拟。计算区域1500 m×1500 m，波长（L）150 m，波高（H）2.5m，周期（T）12.542s,计算网格步长Δx=Δy=5 m，时间步长Δt=0.2s。侧边界采用海绵层滤波。

图7 各个断面比波高对比

图8 渔港三维地形图及控制点分布（单位：m）

图9 T=30 s、90 s、120 s、200 s波浪演进三维立体图

表1三种模型工程测点波高计算对比（单位/m）

不同时刻波面变化（见图9）中可以看到：波浪传播到口门防波堤前趾时发生绕射，部分波能向港池内扩散。防波堤前方很明显看到，入射波和反射波叠加，波幅明显增大。随着时间推移，反射波的影响范围越来越大。同时由于水深和地形的影响波向发生折射，波面变得越来越不规则，在堤前区域形成短峰波。由于防波堤的存在，在港池内波幅明显减小。在波浪折射、绕射、反射联合作用下，可以看到在槽道内出现了明显的波幅。港池内，波浪传播遇到直立堤，同样发生反射与进行波叠加，使得波幅增大。

对比表1可以知道，在开阔水域（1、2、3、4）三类模型计算的波高相差不大，在反射作用较强的区域（5、6、7、8），只有Coulwave模式模拟的结果比较真实，缓坡方程次之，SWAN模式基本没有表现出任何反射效应。当波浪绕射作用较强时（9、10、11），Coulwave模式的优势就更加明显，其他两类模型未能刻画出波浪在堤后的绕射效应。而对于能够反映波浪折射、绕射及反射联合作用的控制点（12—18），缓坡模型及swan模型计算结果在量级上与Boussinesq方程模型偏差较大。

Coulwave模式全面考虑了波浪的折射、绕射、反射效应，而抛物型缓坡方程没有考虑波浪的反射效应，测点处得到的波高小于Coulwave的结果就在意料之中了。而Swan模型对波浪的绕射、反射效应模拟能力相对比较差，得到的结果小于其他两类模式的结果，并且Swan模型在模拟规则波的反射、折射、绕射方面均不如前两类模式模拟的效果好，特别是波高在一两个波长水平范围内变化比较大时，模型的计算结果将出现更大偏差。因此，Swan模式在模拟由近岸海底地形和建筑物引起的波浪绕射和反射效应时会产生较大的误差，这时一般不采用此类模式。

5 结语

Coulwave模型采用分层Boussinesq方程方法来改进经典的Boussinesq方程的色散性，使得改进后的模型的水深范围得到了极大拓展。从本文的模拟结果可以看出：基于Boussinesq方程的Coulwave模型对复杂边界问题处理能力与缓坡方程模型、Swan模型相比具有明显优势，同时Coulwave模型很好的模拟了外海波浪传播到近岸过程中出现的折射、反射、绕射、浅水效应和破碎等波浪在近岸发生的变形现象，模拟结果明显优于另外两种波浪模型，特别是对波浪反射过程的模拟，模型更加适用于复杂地形水域(有人工岛、防波堤等建筑物存在的区域)附近波浪传播变形的模拟研究。

由于海洋中实际波浪一般呈现不规则方向谱形式，波浪不仅具有不规则性，而且具有多向性。本研究仅限于对规则入射波，且与地形的联合作用采用二维两层Boussinesq模型进行了讨论，对不规则方向谱的模拟未能涉及。研究表明分层Boussinesq波浪模型中层数分得越多，其色散性和非线性的效果体现越好，计算模型的精度越高，适用水深范围越广。因此对不规则波入射传播及更高层数的Boussinesq方程模型的探讨将是今后研究的重点。

[1]邹华志.近岸波浪的Boussinesq型方程模型及其应用[J].广东水利水电,2009（4）:1-5.

[2]李孟国,王正林,蒋德才.近岸波浪传播变形数学模型的研究与进展[J].海洋工程,2002,20(4):43-57.

[3]文圣常,余宙文.海浪理论与计算原理[M].北京:科学出版社, 1984.

[4]李孟国,蒋德才.近岸波浪传播折射变形的数学模型概述[J].海岸工程,1999,18（4）:100-109.

[5]李春颖.基于Boussinesq方程的海岸地区波浪数学模型研究[D].天津:天津大学,2003.

[6]王培涛.基于Boussinesq方程的近岸波浪模型的研究及应用[D].北京:国家海洋环境预报中心,2007.

[7]王海珑.基于Boussinesq方程波浪分层数值模拟方法[D].天津:天津大学,2004.

[8]Arthur R S.Refraction of water waves by islands and shoals with circular bottom contours[J].Trans.Amer.Geophys.Un.,1946,27 (2):168-177.

[9]Munk W H,Arthur R S.Wave intensity along a refacted ray.In: Gravity waves[J].Natl.Bur.Standards,1952,Circular521:98-108.

[10]龚崇准,戴功虎.双突堤波浪绕射的数学模型[J].海洋学报, 1984,6(1):99-116.

[11]龚崇准,张长宽.任意反射率边界条件的港湾波浪绕射数学模型[J].河海大学学报,1987,15(6):72-79.

[12]邹志利.水波理论及其应用[M].北京:科学出版社，2005

[13]Boussinesq J.Theory of wave and swells propagated in long horizontal rectangular canal and imparting to the liquid contained in this canal[J].Journal de Mathmatiques Pures et Appliquees,1872,17(2):55-108.

[14]Peregrine D H.Long waves on a beach[J].J Fluid Mech.,1967, 27:815-827.

[15]Madsen P A,Sorensen O R.A new form of Boussinesq equations with improved linear dispersion characteristics:PartⅡ.A slowly varying bathymetry[J].Coast Eng,1992,18:183-204.

[16]Nwogu O.Alternative from of Boussinesq equations for nearshore wave propagation[J].Journal of Waterway,Port,Coastal and Ocean Engineering,1993,119(6):618-639.

[17]Liu P L F,Model equations for wave propagation from deep to shallow water[J].Advances in Coastal Engineering,1994,1:125-157.

[18]Wei G,Kirby J T,Grilli S T,et al.A fully nonlinear Boussinesq model for surface waves:Part I.Highly nonlinear unsteady waves[J].J.Fluid Mech.1995,294:71-92.

[19]Gobbi M F,Kirby J T.Wave evolution over submerged sills: tests of a high-order Boussinesq model[J].Coast Eng,1999,37: 57-96.

[20]Lynett P J,Wu T R,Liu P L F,Modelling wave runup with depth-integrated equations[J],CoastalEngineering,2002,46: 89-107.

[21]Lynett P,Liu P L F,Hwung H H,et al.Multi-layer modeling of wave groups from deep to shallow water,Proceedings of OMAE 2003[R].Ocean Engineering Symposium,2003,13:8-13.

[22]Lynett P,Liu P L F.A two-dimensional,depth-integrated model for internal wave propagation[J].Wave Motion in press.2002,36: 221-240.

[23]Lynett P.A multi-layer approach to modeling generation propagation and interaction of water waves[D].USA:CornellUniversity,2002.

[24]Berkhoff J C W,Booy N,Radder A C.Verification of numerical wave propagation models for simple harmonic linear water waves[J].Coastal Engineering,1982,6:255-279.

Numerical study on the nearshore wave evolution based on Boussinesq-type equations

WANG Pei-tao1，YU Fu-jiang1,2
（1.National Marine Environmental Forecasting Center,Beijing 100081,China；2.Key Laboratory of Research on Marine Hazards Forecasting, State Oceanic Administration,Beijing 100081 China)

The paper reviewed several advanced mathematical models that can be used for the description of wave deformation in port and coastal engineering.The different features of mathematical models,the applicable scale and the evaluation were discussed.The Coulwave(Cornell University Long and Intermediate Wave Modeling)model,based on Boussinesq-type equations,is used to simulate the wave propagation and deformation for several classic topography experiments with specific incident waves.The numerical results agree well with the experimental data.In addition,we compared the results of Coulwave with the SWAN and models with mild-slope equation.The results show that:when considering the interaction of the wave refraction,diffraction, reflection,Coulwave model shows high performance than the model based on mild-slope equations and SWAN model.

Boussinesq-type equations;nearshore wave models;mild-slope equations;SWAN model；Coulwave model

P731

A

1003-0239（2012）04-0007-17

2011-12-29

海洋公益性行业科研专项（200905013）

王培涛（1981-），男，助理研究员，主要从事海啸、风暴潮理论及预警报技术研究。E-mail：wpt@nmefc.gov.cn