邹文峰，王 平, ，刘忠波，房克照，孙家文，张宁川
(1. 大连理工大学 海岸及近海工程国家重点实验室，辽宁 大连 116024; 2. 国家海洋环境监测中心 国家环境保护海洋生态环境整治修复重点实验室，辽宁 大连 116023; 3. 大连海事大学 交通运输工程学院，辽宁 大连 116026)
近几十年来，Boussinesq水波理论得到了长足发展，Madsen和Sørensen[2]、Nwogu[3]、Wei等[4]、邹志利[5]、Madsen等[6]、Lynett和Liu[7]、Liu和Fang[8]都对其进行了研究。基于上述方程的模型也被成功用于模拟近岸区域的波浪运动，最新的Boussinesq水波方程的研究进展可详见Kirby[9]、张尧等[10]、孙家文等[11]的综述。
传统型式的Boussinesq水波方程将三维波浪问题简化成二维问题，较大程度上提高了数值模型的计算效率。当前，基于Wei等[4]方程的FUNWAVE模型、基于Madsen和Sørensen[2]方程的MIKE21-BW模块以及基于Lynett和Liu[7]方程的COULWAVE已得到较大程度的应用。与这些国外流行的模型相比，国内学者在一些基础数值模型方面也取得了一些进展，如柳淑学等[12]基于Beji和Nadaoka[13]的Boussinesq方程建立了有限元模型；
刘忠波[14]、房克照[15]针对邹志利[5]的Boussinesq水波方程建立了一系列有限差分数值计算模型。不同于传统型式的Boussinesq水波方程，Fuhrman[16]、王本龙[17]基于Madsen等[6]的三维Boussinesq方程开发了近岸波浪计算模型。
Liu等[1]给出了最高空间导数为2、3和5的多层Boussinesq水波方程，此多层方程具有最大的稳定空间(Madsen和Fuhrman[18])。最近Liu和Fang[19]、Liu等[20]利用有限差分法开发了最高空间导数为2和3的垂直二维模型，并将数值模型应用到波浪与潜堤相互作用、双色波群演化、非线性聚焦波演化以及海床运动兴波问题(Fang等[21])。
近期，Sun等[22]综合分析了最高导数为2的双层Boussinesq水波方程(Liu等[1])，发现方程中参数取值在0.13～0.25任意选择时，水深在0 Liu等[1]将水体分为两层，分别考虑自由面处的变量η,uη,wη(z=η)、静止水面处的变量u0,w0(z=0)、第一层中间位置处的变量uα,wα(z=zα)、第二层中间位置处的变量uβ,wβ(z=zβ)，水体分层及各变量如图1所示。 图1 双层Boussinesq方程的坐标示意Fig. 1 Schematic diagram of the two-layer Boussinesq model 在自由面上，满足连续性条件及动力学边界条件，对应的方程为： (1) (2) 式中：η为波面升高，uη和wη为自由水平面上的水平速度和垂向速度,uη=(uη，vη)，Ut为U在时间上的偏导，g为重力加速度，∇为水平梯度算子。 自由面上的速度可利用静水位处的速度来表达(Madsen等[6]； (3) (4) 式中：u0和w0表示定义在静水位处的速度，可以通过将z=0代入第一层的速度表达式来确定。 结合两层连接处的速度相等条件、水底运动学边界条件、每层水体速度的泰勒展开式，并引入伪速度来表达，得到： u0=uα*-σ1∇(∇·uα*)+σ2∇wα*-σ3∇h(∇·uα*)-σ4∇h∇2wα* (5) w0=wα*-σ1∇2wα*-σ2∇·uα*-σ3∇h·∇wα*+σ4∇h·∇(∇·uα*) (6) (7) (8) (9) 方程(1)～(9)组成了双层Boussinesq模型的控制方程，具体方程推导过程见文献[1]。1%色散误差下，方程适应无因次水深到kh=19.7，见图2。当α1=0.5时水体变为一层，方程(7)和(8)取消，上述双层模型可退化为单层Boussinesq模型。 图2 方程无因次相速度变化Fig. 2 Phase velocity of the two-layer Boussinesq model 上述控制方程由η,uη,wη,u0,w0,uα*,wα,uβ*,wβ这9个变量和9个方程组成，方程的空间离散采用规则网格下的有限差分方法，在时间步内采用预报—校正—迭代的形式进行求解，直到计算误差落入设定的误差范围内，一次迭代过程分6步求解，数值计算流程详见图3。参数η,uη由方程(1)和(2)求解，在预报时采用三阶Adams-Bashforth的时间格式，校正时采用四阶Adams-Moulton格式，具体构造为： 图3 三维波浪模型数值计算流程Fig. 3 Numerical simulation process of the two-layer Boussinesq model (10) (11) 式中：上标n为第n个时间步的标记，即对应n×Δt时刻的参数值； 针对u0,uα*,uβ*,wβ,wα则分别由式(3)、(5)、(7) 、(9)、(8)求解，将u0,uα*,uβ*,wβ,wα作为未知变量移至方程的左侧，其余项移至方程右侧作为已知量。在空间离散上，对于方程左端未知量的一次、二次导数采用式(12)、(14)的二阶精度格式； fx=(fi+1-fi-1)/(2Δx)+O(Δx)2 (12) fx=(-fi+2+8fi+1-8fi-1+fi-2)/(12Δx)+O(Δx)4 (13) fxx=(fi+1-2fi+fi-1)/(Δx)2+O(Δx)2 (14) (15) (16) (17) 对于w0,wη由方程(6)和(4)直接求解，完成1～6步求解后即得到更新后的变量值，分析校正值和预报值误差，若不满足要求，则对预报值进行更新，并重复1～6步计算调整后的校正值，直至其满足误差要求。 边界造波法和内部造波法是各类Boussinesq数值模型采用的两种主要造波技术，边界造波通常可以选择Stokes线性波理论、二阶波理论或其他波浪理论等，给出边界上的波面位移及速度值。 内部造波法可避免在初始边界采用固定造波板造波时，地形反射给造波板引起的二次反射问题，模型采用了Hsiao等[23]的内部造波法，其源项表达式为： (18) 式中：Δt为数值模型中的时间步长，T为入射波浪周期，Cg为波群速度，H为波高，C为波速，β=20/L2，L为波长，x为网格点坐标，x0为造波源点坐标，ω为波浪角频率，k为波数，t为数值模型当前运行时间。 对于出口边界在本文所有的数值计算中，为提高数值模型稳定性，数值水槽两端均采用海绵层吸收条件，而两侧边界采用反对称的可滑移边界条件； 为了验证本文建立的模型，开展了在深水规则波传播演化、波浪在椭圆形浅滩地形的演化以及波浪在凹形浅滩上的演化三种情况的模拟研究，并与相关解析解和试验结果进行了比较，验证模型的适应性。 深水波浪传播试验中波浪周期为 10 s，波高为 1.0 m，水深分别为 156 m和 468 m，其对应的无因次水深kh分别为2π和6π。在数值模拟中采用边界造波方式，空间网格大小为5 m，时间步长为0.2 s。 图4给出了在空间上模型和波面解析解的对比，图5给出了x=1 000 m处波面历时过程线对比，对比结果显示，无论是kh=2π还是kh=6π，模型给出的波面演变结果均能与解析解较好地吻合，这反映出模型能够精确再现规则波演化过程。 图4 沿波浪传播方向上的数值模拟与解析结果对比Fig. 4 Comparisons of the surface elevations between the numerical model and the analytical solution 图5 x=1 000 m处波面过程数值模拟与解析结果对比Fig. 5 Comparisons of the surface elevations between the numerical model and the analytical solution at x=1 000 m 为考察三维模型对水平速度u及垂向速度w的模拟精度，针对kh=2π工况，图6给出了x=1 000 m处波面位置的u和w，速度的模拟值与解析解吻合较好； 图6 x=1 000 m处波面速度u和w数值模拟与解析结果对比(kh=2π)Fig. 6 Comparisons of the horizontal and vertical velocities at the wave surface between the numerical model and the analytical solution (kh=2π) 图7 x=1 000 m处沿水深速度u和w数值模拟与解析结果对比(kh=2π)Fig. 7 Comparisons of the horizontal and vertical velocities along the water depth between the numerical model and the analytical solution (kh=2π) 图8 不同水深条件下的速度累积误差情况Fig. 8 Accuracy of the horizontal and vertical velocity profiles 为研究波浪传播过程中的折射、绕射、浅化以及非线性作用，Berkhoff等[25]设计了波浪在斜坡上叠加椭圆形浅滩的物理模型试验，相关试验数据已广泛用于验证各种波浪数值模型，如基于N-S方程、各类缓坡方程和不同形式的Boussinesq方程的模型。 在数值模拟中，计算区域为29 m×20 m，斜坡坡度为1∶50，波浪周期为1.0 s，波高为0.046 4 m，波浪正向入射，斜坡梯度方向与波浪入射方向的夹角为20°，内部造波源位于x=4 m处的位置。斜坡旋转后的坐标与计算坐标的关系为： (19) 椭圆形浅滩中心坐标为(0, 0)，浅滩边界定义为： (x1/3.0)2+(y1/4.0)2=1.0 (20) 平底区域及斜坡上的水深hs(单位：mm)为: (21) 浅滩上水深h(单位：mm)为: h=hs-0.5[1-(x1/3.75)2-(y1/5.0)2]0.5+0.3 (22) 数值模拟空间步长为0.1 m×0.1 m，时间步长为0.01 s，为避免过小水深引起数值异常，设置最小水深为0.07 m。分别采用了双层和单层Boussinesq模型计算，得到各个断面的计算结果与试验值对比见图9。 从图9可以看出，单层和双层Boussinesq(BT)数值模型在8个断面的计算结果与试验结果都基本吻合，在D3～D5断面上模拟结果略小于试验数据，在Wei等[4]的数值研究中也出现类似现象，可能是由于强非线性作用，降低了波浪汇聚区的波幅值(Kirby和Dalrymple[26])。多层Boussinesq模型的非线性适用范围更广，其在D4～D5断面上计算的波高值略小于单层模型。 图9 数值结果与实测数据的比较Fig. 9 Comparisons of the wave heights between the numerical model and the measurement results Whalian[27]设计了波浪在凹形浅滩上的非线性折射和绕射物理试验，在该地形上波浪会发生变浅汇聚作用，且浅水非线性作用产生了高次谐波过程，因此该试验被广泛地用来验证各类波浪模型。试验水槽长25.6 m，宽为6.096 m，试验的地形由下式给出： (23) 其中，G=[y(6.096-y)]0.5，0≤y≤6.096。 数值模型中考虑不同周期及波高条件下的波浪演化过程，工况1～3的波浪周期为1 s、2 s、3 s，对应的波高为0.019 4 m、0.015 0 m、0.013 6 m； 在数值模拟中，数值水槽设置长35～45 m，采用内部造波时，周期为1 s和2 s、3 s工况的源位置分别位于x=5 m和x=10 m处，时间步长为0.02 s，空间步长为0.1 m×0.101 6 m。图10给出了双层Boussinesq模型下工况1～3计算结果和试验结果的比较，图11给出了工况4～6计算结果和试验结果对比，图12给出了t=40 s时工况1～3的波面分布。 图10 工况1～3计算波幅与实测值的比较Fig. 10 Comparisons of the computed and measurement results of wave amplitudes under case 1～3 图11 工况4～6计算波幅与实测值的比较Fig. 11 Comparison of the computed and the measurement results of wave amplitudes under case 4～6 图12 不同工况下的波面分布Fig. 12 Wave surface distribution under different cases 从图10和11可见，对于波浪周期为1 s的工况，不同波高下的一次谐波和二次谐波的模拟结果与试验值吻合均较好。对于波浪周期为2 s工况，在非线性较弱的工况2中，一次谐波的模拟值试验值基本一致，而二次谐波和三次谐波的模拟值均略小于试验值； 对Liu等[1]给出的最高导数为2的双层Boussinesq水波方程在矩形网格上进行了空间离散，时间积分上采用混合4阶Adams-Bashforth-Moulton的预报—校正格式，得到了双层Boussinesq方程下的三维波浪数值模型，模型具有较好的色散性和非线性。 针对水深在kh<10深水波浪，模型可以很好地给出波面、波面处的速度以及速度的垂向分布，具有较好的计算精度。针对椭圆形浅滩的物理试验，除波浪汇聚区的波幅值略小于实测值外，模型很好地刻画了复杂地形上的波浪折射、绕射、浅化以及非线性作用。针对凹形浅滩地形上波浪浅水非线性作用产生的高次谐波过程，在周期1 s和2 s的工况中，模型均能很好地得到各次谐波值，且随着非线性的增大，高次谐波更为显著。通过不同典型算例的验证，本文建立的双层Boussinesq模型对强非线性波浪的演化具有较好的模拟精度。1.1 双层Boussinesq方程
Liu和Fang[8])，以获取更准确的非线性特征，二者的关系式为：
1.2 数值求解
n+1对应(n+1)×Δt时刻的参数值；
其他依此类推。
对于方程右端已知量的一次导数采用四阶精度式(13)求解，二次导数采用二阶精度式(14)求解，具体格式为：
1.3 入射边界与出口边界
为使程序运行稳定，模型采用了Kirby等[24]给出的光滑滤波技术。2.1 深水中规则波传播过程
图7给出u和w的模拟值与解析解的对比，除第一层和第二层连接处附近略有差异外，其余深度的速度均能较好再现。针对kh=6π工况，尽管可以较好再现波面及波面处速度的演化过程，但由于其水深已超过本文双层Boussinesq模型在速度分布特性上的适用范围，其垂向速度累积误差已超过6%，见图8，这里不再给出对比结果。
2.2 椭圆形浅滩上的波浪传播变形
2.3 凹形浅滩上的波浪传播变形
工况4～6的波浪周期为1 s、2 s、3 s，对应的波高为0.039 0 m、0.029 8 m、0.029 2 m。随着入射波高的增大，波浪的非线性逐渐增强。
当非线性增强时(工况5)，一次谐波和二次谐波的模拟值均略大于试验值，三次谐波的模拟值与试验值吻合较好。对于波浪周期为3 s工况，模型结果和试验结果差异较大，目前所见的模型结果多存在这一现象(Chen和Liu[28])。总体而言，模型非线性会对数值结果产生一定影响，本文建立的双层Boussinesq模型对强非线性波浪的演化具有较好的模拟精度。