刘忠波， 刘泽， 王铖龙， 王彦， 房克照

(1.大连海事大学交通运输工程学院， 大连 116026； 2.大连理工大学海岸和近海国家重点实验室， 大连 116024)

Boussinesq型水波方程是一类时域水波模型，此类模型最显著的基本特征是控制方程中变量不对垂向坐标z求导，等同于将三维波浪问题降至水平二维问题，同时Boussinesq模型也无需追踪自由面问题，从而大大降低了计算时间，这促进了其在海岸波浪水动力模拟研究的广泛应用。

考虑到模型在实际工程应用的可能性，Boussinesq型水波方程的变量对空间坐标求导一般不超过3，关于这一点，很多理论和数值应用文献均有陈述[1-4]。吴昕炜等[5]利用Boussinesq波浪(Boussinesq wave，BW)模型对波浪在实际地形中的传播变形规律进行了数值模拟研究，并利用相关试验结果进行了验证，该模型最高空间导数为3。最高空间导数为2的弱非线性Boussinesq型水波方程由Zhao等[6]给出，该方程具有Padé(2，2)的色散性能，一些学者研究泥沙动力学问题前选择了这个模型开展相关的波浪水动力研究，但弱非线性是该模型的主要限制。Liu等[2]、Chazel等[7]给出最高空间导数为2的多层三维Boussinesq水波方程；Liu等[8]进一步分析了Liu等[2]多层方程的理论性能，多层模型显著提高了方程的线性和非线性性能。作为代价，层数越多，计算时间也会相应增加。Liu等[8]模拟水槽聚焦波算例的结果表明，3层(11个方程)比2层(9个方程)计算时间增约为23%。因此，选择单层模型也具有研究价值，尽管在1%色散误差下原方程的适用最大因次水深kh=3.8，但精确的色散关系仅在kh<1.23左右(0.1%误差)。与色散性能更为优良的多层模型不同，单层模型在实际工程中被选择应用的概率更大，因此有限水深时的性能更应予以重视。

除色散性外，浅化和非线性等性能是推导Boussinesq型水波方程必须检测的重要指标[2,8-10]。以变浅性能为例，刘忠波等[11]对水平二维Boussinesq方程[9,12-13]中的速度用变换速度取代，并通过分析变浅性能获取了相关参数值。事实上，变换速度中引入水深变化项使得控制方程含有的独立色散参数和独立变浅参数的个数相等，这是方程具有更优良变浅性能的关键因素。与水平二维方程不同，三维Boussinesq型方程含垂向速度，确保了方程具有更优良的非线性性能[2,7-8,14]。

为了平衡效率和模型精度，现选择Liu等[2]给出最高空间导数为2的单层三维Boussinesq水波方程为研究对象，对其理论性能的改进方面进行初探，并建立相应的数值模型，通过相关试验来验证理论模型的改进。

1 三维Boussinesq型方程的改进

Liu等[2]在无旋、无黏假设下，首先沿垂向方向将水体分为N层，并规定每层中点水深处的速度为特征速度，在每一层中，对该特征速度关于z做泰勒展开，并利用连续方程和无旋假定，则可消除掉变量对垂向坐标z的导数，可以得到一组以N层特征速度表达的Boussinesq型水波方程。其次，对N个特征速度变量引入L变换算子得到N个变换速度(又称为伪速度或计算速度)，引入变换速度的核心目的是提高方程的色散性，最终推导出具有更高性能的三组N层Boussinesq型水波方程。选取最高导数为2的单层方程作为研究对象，方程中的主要变量为波面位移η、自由表面处水平和垂向速度uη和wη、静止水位处的速度u10和w10以及变换速度uE和wE。

在自由表面处，精确的运动学和动力学边界条件为

(1)

(2)

U=uη+wη∇η

(3)

式中：U为定义在自由表面处速度势的水平梯度；g为重力加速度；∇为水平梯度算子。

在海底，边界满足条件为

ub·∇h+wb=0

(4)

式(4)中：ub和wb分别为海底水平速度和垂向速度；h为静止水位下的水深。

自由表面处速度与静止水位处速度满足的关系表达式为

(5)

从水底到静止水平面的速度场可表达为

(z-zα)∇zα∇·uE+(z-zα)∇wE+

(6)

(z-zα)∇zα·∇wE-(z-zα)∇·uE-

(7)

将z=0和z=-h分别代入式(6)和式(7)可得到静止水位处和海底处的速度。在式(6)和式(7)中，引入的变换速度(uE,wE)与静止的中间水深处速度(uα,wα)之间的L变换表达式为

(8)

(9)

与式(8)和式(9)给出的L变换表达式不同，引入的变换如下。

uα=uE+α1(αh)2∇(∇·uE)+β1α2h∇h∇·uE

(10)

wα=wE+α2(αh)2∇·(∇wE)+β2α2h∇h·(∇wE)

(11)

式中：α1、α2为色散参数；β1、β2为变浅参数。

由此，速度场表达式为

(z-zα)∇zα∇·uE+(z-zα)∇wE+

[-2α2(z-zα)αh+(z-zα)2]∇zα∇·

(∇wE)+β1α2h∇h∇·uE+

β2α2h(z-zα)∇h∇2wE

(12)

(z-zα)∇zα·∇wE-(z-zα)∇·uE-

[-2α1(z-zα)αh+(z-zα)2]∇zα·

∇(∇·uE)+β2α2h∇h·(∇wE)-

(z-zα)β1α2h∇h·∇2uE

(13)

类似于文献[10]，在水底方程中引入变浅系数β3，最终水底运动学方程可写为

(14)

式(1)、式(2)、式(5)、式(14)以及静止水位处速度与变换速度的表达式[z=0代入式(12)和式(13)]，构成了新三维Boussinesq型水波方程，当α1=α2=1/6，β1=β2=0时，则可以转化为原来的方程。

2 方程的基本性能

在立面二维和常水深情况下，通过傅里叶分析，则可分析出方程的相速度、和差频、速度分布特征等，分析的详细过程参见文献[10]。

2.1 相速度

方程的色散关系表达式可写为

(15)

式(15)中：ω为波浪圆频率；k为波数。

当表达式(15)确定的相速度与解析相速度在0

图1 无因次相速度Fig.1 Non-dimensional wave phase celerity

2.2 速度分布

为了表征速度分布的精度，采用Madsen等[14]给出的公式考察本文方程与解析解之间的误差，公式为

(16)

式(16)中：us(0) 和ws(0)分别为z=0时的Stokes波线性解。图2给出了改进值与原参数下的速度分布误差的比较，整体来看，改进参数后的水平速度和垂向速度误差更小。

图2 速度分布误差Fig.2 The error of the velocity profiles

2.3 不规则波中的和差频

一般来说，对于不含垂向速度的Boussinesq水波方程，差频性能是核心短板，关于这一点，文献[9,12,15]中均有表述。当选取色散参数α1=0.170，α2=0.158(改进值)和α1=α2=1/6(原参数)情况下的，Boussinesq方程的和差频解析结果如图3所示，其中，k1和k2是两个不同频率波浪的波数。无因次和频和差频均以对角线为对称线，对角线上面部分为和频，下面部分是差频，图3仅给出一半的图形。改进后，1%误差下和频、差频适用范围均超过了kh=4。1%误差下原参数方程的和频、差频适用范围仅约等kh=1.7和kh=2.2。这表明，本文方程具有更优良的和差频特性。此外，无论是原方程还是本文方程，在非线性性能方面特别是差频性能方面均优于传统型的Boussinesq水波方程，以Gobbi等[16]的方程为例进行说明，该方程精确到4阶保留了所有非线性项，色散适用水深约为kh=6，其1%误差下差频的适用范围仅为kh=1左右。

图3 无因次和差频Fig.3 Non-dimensional amplitudes of sub-harmonics and super-harmonics

2.4 变浅性能

变浅性能代表着线性波浪由深水到浅水演化过程中的浅化特征，在0

图4 方程的变浅梯度Fig.4 The shoaling gradient of the present model

3 数值验证

对Boussinesq型水波方程可构建相应的垂向二维模型，数值模型中对时间差分采用3阶Adams-Bashforth格式和4阶Adams-Moulton格式，分别用于预报和校正阶段。当预报和校正的所有变量的误差控制在一定范围内，则当前计算步结束。波浪的产生采用Hsiao等[17]的域内造波法，并在两边界上采用2倍波长的海绵边界层进行消波。

为对比改进前后模型的差别，选择了对Luth等[18]潜堤上规则波传播变形的试验数值模拟。先简要分析波浪在潜堤上的演化过程，波浪的传播变形是一个复杂的物理过程，其演化过程中主要包括波-波相互作用、变浅与反变浅、相位锁定与解锁甚至波浪破碎等。为了精准模拟演化过程，数值模型应具有足够精度的色散性、变浅性和非线性性能。

图5 计算波面位移与文献[14]试验结果的比较Fig.5 The comparisons of the computed surface elevations with reference[14]

数值模拟中，选取了周期T=2.02 s，波高为0.02 m的工况，时间步长采用0.01 s，空间步长采用0.03 m。图5所示为改进前与改进后的计算波面位移与实验结果的比较。改进后的本文方程和原方程的模拟结果与实验结果吻合程度均较高，但在浪高仪位置x=19.0 m和x=21.0 m处，本文方程比原方程在3次峰波形的准确性上更高。具有Padé (2, 2)色散关系的各种Boussinesq数值模型无法准确捕捉到最后3个位置点的波面位移。具有Padé (4, 4)色散关系的方程，若变浅性和非线性性能足够，则可以较为精确捕捉到这3个位置点的波面，倘若变浅性能不佳，则也无法精准模拟这一过程[19]。

4 结论

通过引入带参数的变换速度，导出了一组最高空间导数为2的三维Boussinseq型方程。并进行了理论分析和潜堤上波浪传播变形模拟研究，主要得出以下结论。

(1)本文方程比原文献在线性(相速度和速度分布)、和差频等方面均有不同程度的改善，特别是相速度与解析解在0

(2)以实验结果为标准，潜堤后方高次谐波释放，原方程对三次峰的捕捉不够精确，这反映出原方程在色散性精度不足，但值得注意的是，本文方程模拟的主峰与实验结果吻合程度在最后两个浪高仪位置与实验结果的吻合程度略次于原方程，关于这一点，更广泛的数值模拟验证将会在下一步研究中进行。