Boussinesq水波方程新型数值解法

房克照，孙家文，刘忠波，尹 晶，张 哲

(1.大连理工大学海岸和近海工程国家重点实验室，辽宁大连 116024;2.国家海洋环境监测中心，辽宁大连 116023)

由于在模型稳定性、波浪破碎处理和海岸动边界问题方面较传统有限差分方法具有优势，近十年来，将有限体积方法应用于Boussinesq类水波方程的数值求解成为研究热点，形成了一类有限体积和有限差分方法混合求解的数值格式［1-9］。

建立该类混合数值格式的主要思路如下，将Boussinesq类方程写为守恒格式，对通量项采用有限体积方法计算，剩余项采用有限差分方法求解。当波浪破碎在水-陆交界附近，色散性和高阶非线性对流体运动贡献较小，不参与计算，则Boussinesq方程退化为完全非线性浅水方程。完全非线性浅水方程属于典型双曲守恒律方程，有许多基于有限体积方法的高分辨率数值格式［10］(多以黎曼间断解问题为理论基础构筑)，因此其足以胜任破碎波浪(视为水跃)和水-陆交界(视为间断解)的数值计算。

由于Boussinesq方程中大多数项仍通过有限差分方法求解，有限体积方法的数值实现重点体现在网格界面处单调性数值通量的计算。通常，采用中心格式和迎风格式进行［10-12］。中心格式不需要详细了解特征波传播状态，构造简单，但精度低。迎风格式(如HLL格式和Roe格式)精度高，但需要精确或者近似求解黎曼问题，计算量大，构造过程相对繁琐。现有混合格式大都采用HLL格式进行数值通量计算［1-4，6-9］。也有部分学者采用Roe格式计算［5］，但其每个时间步都涉及到特征矩阵的求解和分裂计算，十分繁琐。HLL和Roe格式均建立于详细的黎曼解状态基础上，并且均存在一些不足，譬如HLL格式在水平二维问题时需要扩展至HLLC格式以弥补不能考虑接触波的不足，Roe格式某些情况下不满足熵增条件，对于实际应用而言，较准确地求解黎曼解不是一件易事［10］。Toro等近年来发展了MUSTA格式［11-12］，与传统中心格式和迎风格式不同，其通过控制方程直接构造数值通量，是一种既具有中心格式的简单性同时具备迎风格式精确性的新型数值格式，其精确性和有效性在求解非线性浅水方程时得到了验证［11-13］。

将MUSTA格式应用于Boussinesq类水波方程的求解，建立有限体积/有限差分混合求解数值格式。除对所建立模型进行验证外，重点比较MUSTA和HLL两种格式在计算精度、计算效率以及程序编制等方面的表现。

1 数学模型

1.1 适用于快速变化地形的Boussinesq水波方程

以波面η和通量q表达的二阶Boussinesq方程［14］为:

式中:h为静水深;d=h+η为当地水深;g为重力加速度;q=du为断面流量(u为水深平均速度)，下标x和t分别表示变量对空间和时间的偏导数;B=1/15为自由参数，通过优化方程的色散性得到［15］。上述方程为Madsen和Sorensen方程［15］的扩展，更适用于在地形快速变化水域中应用。忽略式(2)中色散项ψx，方程退化为完全非线性浅水方程。

考虑海底摩擦，并将上述方程写成守恒形式如下:

式中:f为底摩擦系数。

1.2 方程空间离散

将计算域在空间、时间上做如下离散:xi=i(i=1，…，N)，tn=nΔt，其中Δx、Δt分别为空间、时间网格步长。在有限体积［xi-1/2，xi+1/2］×［tn，tn+1］内对控制方程(4)进行积分并应用格林定理，可得

1.3 数值通量计算

高分辨率数值通量Fi+1/2的计算均通过求解黎曼解问题进行，即求解初始值问题:

第一步:通量计算

若k=K，流程结束。

第二步:通过控制方程求近似黎曼解

第三步:返回第一步。

图1 MUSTA格式构造示意图(上标表示第k步值)Fig.1 The sketch of MUSTA scheme(subscript k denotes the value after kth stage)

文中，k=1，2，3时对应的格式分别简称为MUSTA-1，MUSTA-2和MUSTA-3。MUSTA格式具有简洁、直观的表达形式，不需要求解局部黎曼解的特征波状态。而传统的HLL格式［1-4，6-9］需要判断黎曼解特征波状态，需要计算特征波波速，并以此进行通量计算，计算繁琐。因此，MUSTA格式在计算效率和程序编制方面均要比HLL格式具有优势。

上述MUSTA格式仅为一阶精度，为提高精度，采用四阶精度的状态插值方法［17］对界面左右变量进行重构，利用重构后的变量进行上述计算。

1.4 波浪破碎和水-陆动边界处理

在混合格式中，波浪破碎发生时，Boussinesq方程退化为完全非线性浅水方程，处理为水跃。根据建议［1］，当波面升高同水深比达到0.8时，认为波浪发生破碎，控制方程(1)～(3)中的Boussinesq高阶非线性项和色散性不参与计算，方程退化为浅水非线性方程。对于水-陆动边界，采用简单有效的薄层水体法处理［13］，即定义一极小水深，本文取0.001 m，若网格水深大于该值视为水域;若网格水深小于该值，视为干网格，水深强制赋值0.001 m，流速为零。

1.5 时间积分以及边界条件

时间积分通过具有TVD性质的三阶龙格-库塔方法进行［1］。计算域两端设置为固壁边界，模型采用在质量方程中增加源项的方法产生波浪，同时视需要在计算域末端设置海绵吸收层吸收波浪［1］。

2 模型验证

本节将通过几个典型算例的数值模拟，对所建立模型进行验证。除采用MUSTA格式进行计算外，在保持其他计算条件一致情况下，也采用HLL格式进行了计算，以便两种格式的比较。

2.1 孤立波在等水深水槽中的传播

孤立波是波浪色散性和非线性平衡制约的典型代表，其在常水深水槽中长距离传播的数值模拟是检验混合类格式的有力工具。计算时，计算域长度500 m，水深1.0 m。在计算域内给出孤立波作为初始条件，波幅a=0.6 m，幅值中心位于x=60 m处，网格尺寸Δx=0.10 m，计算域两端设置5 m长海绵层。采用MUSTA-1格式的计算结果在图2中给出。可见，经过长距离的传播，波形和幅值保持不变，表明所建立的数值格式没有引入数值耗散以及伪色散。

针对这一问题，文献［8］给出了解析解。t=100 s时，数值解同解析解的对比在图3给出，可见数值解同解析解高度吻合。此外，数值计算的传播速度为4.038 m/s，同解析解4.037 3 m/s吻合。

图2 孤立波传播过程中不同时刻波面升高Fig.2 The surface elevations at different moments during the solitary wave propagation

图3 t=100 s时计算波面同解析解的比较Fig.3 Comparison of computed and analytical surface profiles at t=100 s

为定量考证MUSTA格式同HLL格式的区别，引入以下误差判断标准［13］:

式中:Y代表要比较的变量，下标i为网格标记，上标sim和exact分别为数值解和解析解。计算结果在表1中给出，同时给出各种格式的CPU耗时。对比可见，对MUSTA格式而言，随计算步数增加，模型精度提高，但计算时间也显著增加。

表1 不同格式误差汇总Tab.1 A summary of error norms for various MUSTA and HLL schemes

当K=3时，计算精度的提高几乎可以忽略。MUSTA格式计算效率均高于HLL格式。其中，同HLL格式相比，MUSTA-1精度略有降低，但计算耗时显著减少。综合考虑计算精度、计算效率以及程序编制难易程度，本文以下计算均采用MUSTA-1格式进行。这同其他学者［12-13］推荐在求解非线性浅水方程时使用MUSTA-1格式一致。

2.2 规则波在潜堤地形上的传播

规则波跨越潜堤传播是非常复杂的波浪演化过程，涉及波浪的变浅、高次谐波产生释放等，是检验模型色散性和非线性综合性能的有力工具。这里针对Beji和Battjes实验［17］中工况(a)进行模拟，入射波波高H=0.02 m，周期T=2.02 s，实验布置见图4。计算时，空间步长Δx=0.02 m，计算域两端设置1.5倍波长海绵吸收层。图5给出6个浪高仪位置上计算结果同实验数据的对比。MUSTA-1和HLL格式给出的数值结果几乎相同。在x=13.5 m以前，数值结果同实验数据吻合。在此之后，由于控制方程色散性的不足，模拟波面同实验数据差别较大。经统计发现，MUSTA-1较HLL格式CPU耗时减少约5%。

图4 规则波跨越潜堤传播实验布置示意(Beji和Battjes［17］)Fig.4 The sketch of experimental setup of regular waves propagation over submerged bar(Beji and Battjes［17］)

2.3 潜礁地形上孤立波传播

为验证模型，针对孤立波在潜礁上传播过程进行数值模拟，其中涉及到波浪破碎、水-陆动边界、水跃形成传播等复杂过程，是检验模型较为苛刻的算例。这里，采用Roeber实验［2］中第三组进行模拟，该组实验地形设置类似自然界中带潟湖潜礁，非常具有代表性。潜礁峰露出水面0.06 m，而潜礁平台淹没于水下0.14 m［2］。

图5 波面升高数值结果和实验数据对比Fig.5 The comparison of surface elevations between numerical model and experimental data

图6 波面升高数值结果和实验数据对比Fig.6 The comparison of surface elevations between numerical model and experimental data

计算时，底摩擦系数取0.007 5，左侧设置海绵层5 m，右侧为完全反射边界，Δx=0.10 m。入射孤立波波高0.75 m，造波板处水深2.5 m。数值模拟结果在图6中给出。可见，孤立波在潜礁前坡传播过程中，波前变陡，破碎发生形成水跃。t’=59时刻起(t’=t/(g/h0)1/2为无因次时间，h0为造波板处水深)，波浪开始跨越潜礁峰传播至潟湖区域。在潜礁峰后激起水跃，并保持尖锐形状向前传播(至t’=80.11)。随后(t’=91.26，100.61和117.11)，被完全反射后反向传播，至深水区域由于色散性作用加强也出现诸多短波。数值结果同实验数据的良好吻合，表明模型也适用于这种特殊地形条件下孤立波传播和破碎的数值模拟。同本文其他算例一致，MUSTA-1和HLL格式给出的数值结果接近，但前者计算速度提高约6%。

3 结语

将MUSTA格式应用建立求解Boussinesq水波方程的混合数值格式。针对一维守恒格式的控制方程，采用有限体积方法求解数值通量，剩余项利用有限差分方法求解。时间积分、边界处理等同现有的混合格式基本类似，但本文采用MUSTA格式计算数值通量。除进行模型验证外，重点对MUSTA和HLL格式进行了对比研究。综合考虑计算精度、计算效率、程序编制和实际应用，MUSTA-1格式较普遍应用的HLL格式具有优势，值得在Boussinesq类水波方程的求解方面进行推广和应用。

［1］ Shi F，Kirby J T，Harris J C，et al.A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation ［J］.Ocean Modelling，2012，43-44:36-51.

［2］ Roeber V，Cheung K F.Boussinesq-type model for energetic breaking waves in fringing reef environments［J］.Coastal Engineering，2012，70:1-20.

［3］ Roeber V，Cheung K F，Kobayashi M H.Shock-capturing Boussinesq-type model for nearshore wave processes［J］.Coastal Engineering，2010，57(4):407-423.

［4］ Bradford S F，Sanders B F.Finite volume schemes for the Boussinesq equations［J］.Ocean Wave Measurement and Analysis，2001:953-962.

［5］ Erduran K S，Llic S，Kutijia V.Hybrid-finite volume finite-difference scheme for the solution of Boussinesq equations［J］.International Journal for Numerical Methods in Fluids，2005，49:1213-1232.

［6］ Cienfuegos R，Barthelemy E，Bonneton P.A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations.Part II:boundary conditions and validation ［J］.International Journal for Numerical Methods in Fluids，2007，53:1423-1455.

［7］ Shiach J B，Mingham C G.A temporally second-order accurate Godunov-type scheme for solving the extended Boussinesq equations［J］.Coastal Engineering，2009，59(1):32-45.

［8］ Orszaghova J，Borthwick A G L，Taylor P H.From the paddle to the beach-A Boussinesq shallow water numerical wave tank based on Madsen and Sorensen’s equations［J］.Journal of Computational Physics，2012，231(2):328-344.

［9］ Dutykh D，Katsaounis T，Mistotakis D.Finite volume scheme for dispersive wave propagation and runup［J］.Journal of Computational Physics，2011，230(8):3035-3061.

［10］ Toro E F.Riemann solvers and numerical methods for fluid dynamics［M］.Springer Verlag，2009.

［11］ Toro E F.Musta fluxes for systems of conservation laws［J］.Journal of Computational Physics，2006，216(2):403-429.

［12］ Toro E F.Musta:A multi-stage numerical flux［J］.Applied Numerical Mathematics，2006，56(10-11):1464-1479.

［13］ Guo W D，Lai J S，Lin G F.Finite-volume multi-stage schemes for shallow-water flow simulation［J］.International Journal for Numerical Methods in Fluids，2008，57(2):177-204.

［14］ Kim G，Lee C，Suh K.Extended Boussinesq equations for rapidly varying topography［J］.Ocean Engineering，2009，36:842-851.

［15］ Madsen P A，Sorensen O R.A new form of the Boussinesq equations with improved linear dispersion characteristics.Part 2.A slowly-varying bathymetry［J］.Coastal Engineering，1992，18:183-204.

［16］ Yamamoto S，Kano S，Daiguji H.An efficient CFD approach for simulating unsteady hypersonic shock-shock interference flows［J］.Computers＆ Fluids，1998，27(5-6):571-580.

［17］ Beji S，Battjes J A.Experimental investigation of wave propagation over a bar［J］.Coastal Engineering，1993，19(1-2):151-162.