渗透介质中波浪传播的二阶Boussinesq水波数值模拟

琚海军，刘忠波，刘 勇，房克照

(1.中国海洋大学 工程学院，青岛 266100；2.大连海事大学 交通运输工程学院，大连 116026；3.大连理工大学 海岸和近海工程国家重点实验室，大连 116023)

为了保护海岸，国内外修建了不少带有多孔介质的抛石防波堤等海岸工程。这些建筑物本身带有很大孔隙，它们一方面可以将远海传来的波浪反射回去，另一方面也可通过建筑物自身内部的摩擦阻力消耗掉一部分波浪，对波浪运动产生明显的衰减影响，从而起到了保护海岸的效果。考虑到近岸海域的海床上覆盖着大片透水程度不同的颗粒介质，这些实际孔隙介质的渗透性也直接影响到岸滩上的波浪水流等水动力特性，所以渗透性对波浪运动产生的影响也是不能忽略的。

关于孔隙海床的波浪水动力特性研究，多集中于孔隙海床上波高衰减问题上，并且常借助于实验或理论分析等手段进行展开。过去一些研究者对渗透海床上的波浪水动力问题进行了初步探究和分析，但是这些模型不能直接应用于研究较大范围内可渗海床上的波浪信息，而模拟近岸海域附近的数学模型主要有基于缓坡假定的缓坡方程和基于小振幅假定的Boussinesq型水波方程。缓坡方程的模型则多以线性波浪为研究对象且仅能适应于考虑弱非线性的情况，无法用于非线性较强的波浪与波浪、波浪与水流等相互作用问题的研究。而高阶Boussinesq型水波方程在模拟强非线性波浪、波浪与波浪相互作用等相关问题上则显示出较好的适用性。

但是国内外研究者在推导这些方程的时候，为了简化问题，通常假设海底不透水[1-7]，仅有少数研究者进行可渗海床上这方面的研究[8-13]。而后一类模型也只能考虑上部水体为自由水体，下部为可渗多孔介质的波浪运动情况，这使得它们很难直接应用于模拟波浪在多孔介质中传播变形。Hsiao等[12]通过欧拉方程给出了适应于多孔介质中波浪传播变形的高阶Boussinesq型水波方程，模型是由可渗透多孔介质某一深度处速度和波面升高表达的。刘忠波等推导了两层渗透介质上的Boussinesq型水波方程[13]，我们通过对其一组方程进行了简化，我们导出了一组高阶Boussinesq型水波方程，对其进行了理论研究，并在此基础上进行了数值研究工作。

1 高阶Boussinesq型水波方程的理论推导

波浪在多孔介质中的传播见图1，η(x,y,t)表示波浪在多孔介质中的波面升高(以静水面为起算面)，多孔介质的厚度为h(x,y)，ψ(x,y,z,t)表示速度势。

图1 波浪在多孔介质中传播的示意图Fig.1 Sketch of wave propagation in porous media

在无因次化中常用的三个特征量为：a、h0和l0，为了确认方程中不同项的相对重要性，引入以下无因次

(1)

在Boussinesq类水波方程推导过程中，常引入非线性特征参数和色散性特征参数

(2)

为了简便起见，忽略了上标，则对应的控制方程和边界条件可写为

μ2▽2ψ+ψzz=0，-hb

(3)

(4)

μ2(ηt+ε▽ψ·▽η)=ψz，z=εη

(5)

μ2▽ψ·▽h+ψz=0，z=-h

(6)

这里α=α1+α2|U|，其中α1和α2分别代表层流和紊流阻力系数；cr=1+(1-n)cm/n为惯性系数，cm为附加质量系数，详细可参见文献[8-9]。

将速度势展开如下

(7)

分别对速度势求1次和2次导数，并将结果带入控制方程(3)中可得

(8)

所以速度势表达式可以表示为

(9)

将(9)代入(6)中，经整理可得

(10)

式中：w0和u0分别表示在静止水位初的波浪垂向和水平速度。

对方程(4)求导，并将(10)代入，整理可得

(11)

类似地，整理(4)得到的连续方程如下

ηt+▽·Q=0

(12)

引入水深积分平均速度

(13)

对(13)进行整理，可得

(14)

连续方程可写为

(15)

将(14)代入(11)中，整理可得

(16)

方程(15)和方程(16)组成了一组Boussinesq方程。

沿水深积分平均速度表达的方程色散关系式与经典Boussinesq方程的一致，因此该方程的适用水深有限。为了拓展方程的适用水深，在动量方程(16)中引入刘忠波等[7]建议的表达式

(17)

式(16)中的两个参数可以改善方程的色散性，当增加(17)后，方程(16)可以写成

(18)

由(15)和(18)组成的方程，当不考虑α和cr时，方程转化为刘忠波等给出的适合非渗透Boussinesq水波方程[7]。同时，当文献[13]不考虑第二层渗透时，也可以转化为本文的研究方程。

2 方程的色散性与衰减率特征分析

在不考虑海底地形变化情况下，忽略非线性和忽略水深变量对x的导数，对由方程(15)和(18)组成的模型进行色散性分析，可分别得到如下色散关系式

(19)

式中:ω代表频率；Kb为复波数，Kb=Krb+iKib，虚数i满足i2=-1,β=β1+β2。

为了对比不同方程的相速度和ki，采用了Hsiao等[12]文献中给出的解析色散关系表达式

ω2(cr+α1i/ω)=gKtanh(Kh)

(20)

这里K=krs+ikis为复波数。

对比(19)和(20)，可以测定不同情况下方程的相速度与衰减率，反推参数的取值。在这里, 从0.1 s/m变化到4 s/m，krh从0.1变化到3，通过比较，以下面两种判断式为标准我们给出第一组解和第二组解。

[1-(krs/krb)]2+[1-(kis/kib)]2=min (21)

由于方程中有两个参数，而根据以上对比，只能确定一个参数，因此，另外一个参数则只能通过非渗透性情况下的变浅性能获取，具体可参照文献[7]，这里直接给出结果，见表1。对应的无因次相速度、衰减率以及复波数见图2和图3。从衰减率和相速度来看，第一组参数的适用范围看略大一些。从复波数变化来看，第二组的适用范围更大。整体来看，二者差异性并不十分明显，在后文的数值模拟中采用了第二组参数。

2-a 无因次相速度2-b 无因次衰减率2-c 无因次复波数图2 方程的无因次相速度(a)、无因次衰减率(b)和无因次复波数(c)(第一组)Fig. 2 The phase velocity (a), the dimensionless attenuation rate (b), and the dimensionless complex wave number (c) of the equation. (Group 1)

3-a 无因次相速度3-b 无因次衰减率3-c 无因次复波数图3 方程的无因次相速度(a)、无因次衰减率(b)和无因次复波数(c)(第二组)Fig. 3 The phase velocity (a), the dimensionless attenuation rate (b), and the dimensionless complex wave number (c) of the equation. (Group 2)

3 算例

3.1 孤立波穿越多孔直立防波堤

针对上面给出的模型在非交错网格下进行了离散，时间步进采用了三阶Adams-Bashforth格式进行预报，进而采用四阶的Adams-Moulton格式进行校正[7]。当预报值和校正值满足一定误差时，当前步结束，否则重新更新变量，重新利用校正步进行校正。

图4 实验布置Fig.4 Experimental arrangement

用本文提出的理论模拟Vidal等进行的孤立波与多孔直立式防波堤之间的水槽实验[14]，测量反射系数和透射系数。在实验中，防波堤的的宽度为20 cm或40 cm，实验水深从25 cm到32 cm不等，防波堤材料的孔隙率为n=0.44，粒径采用两种d50=1.43 cm和2.43 cm，实验中采用的两个浪高仪分别安置在防波堤前3.9 m和后2 m的地方，如图4。

在数值模拟中，模型采用了空间步长△x= 0.01 m，时间步长为△t=0.005 s，图5给出了本文模型模拟结果与试验结果的比较。在本文计算中，针对d50=1.43 cm，参数α1和α2取值为18 s-1和234.45 m；针对d50=2.43 cm，参数α1和α2取值为6.23 s-1和137.97 m，其中α1的值参照Gent给出的系数确定[15]，α2根据Lin和Karunarathna给出的系数确定[16]。由图5可见，其反射系数和透射系数的数值结果与试验结果相对吻合。

5-a a=20 cm,h=30 cm,d50=1.43 cm5-b a=20 cm,h=30.2 cm,d50=2.43 cm5-c a=40 cm,h=30 cm,d50=1.43 cm5-d a=40 cm,h=31.7 cm,d50=2.43 cm注:符号为Vidal等(1988)试验数据;实线:数值反射系数;虚线:数值透射系数图5 孤立波与多孔直立式防波堤相互作用的反射系数和透射系数比较Fig.5 Comparison of reflection and transmission coefficients of solitary wave interaction with porous breakwaters

图6 实验布置Fig.6 Experimental arrangement

3.2 渗透地形的数值实验

与文献[17]物理模型实验类似，仅在8～10 m(堤顶)上面放置出水渗透介质。在这里只考虑0.4 m水深，波高为0.02 m，周期T=2.02 s。水槽长23 m、宽0.8 m、高0.8 m，水槽中设置一梯形障碍物(设比尺为l:50时，对应实际地形中的沙坝)。梯形顶端长2.0 m，顶端的中心位置距离水槽左端为9.0 m。水深为0.4 m，其实验装置如图6。

…… (改进的Boussinesq模型) ——— (对照模型)

图7 改进Boussinesq方程数值波面时间历程图对比Fig.7 Comparison of the numerical wave surface of the improved Boussinesq equation

假设本文方程参数α2=0，以α1取2.0的计算(以α1=0为对照模型)，T= 2.02 s，附加水体质量为0.34，孔隙率为n=0.5，得出对比结果。图7对应参数α1取2.0时各个位置的波形图。通过与对照模型的比较，其结果验证了改进的Boussinesq型水波方程的有效性，波浪衰减非常明显，反映出渗透对波浪的影响。

图8和图9是α1分别等于0.2，0.5，1.0，2.0的均方波高(统计最后一个周期)和t=40 s时的空间波形比较的结果。我们也可以明显得到参数α1会影响波浪的衰减，并且在同等条件下，随着α1的系数值增加，波高衰减增大。

图8 均方波高(对照组:α1=0)Fig.8 Mean square wave height (control group:α1=0)图9 t=40 s时的空间波形(对照组:α1=0)Fig.9 Spatial waveform (t= 40 s,control group:α1=0)

4 结论

在考虑多孔拖曳阻力和惯性阻力的基础上，推导了一组高阶Boussinesq型水波方程。在常水深情况下，分析了模型对应的色散关系式，并与解析解进行了比较。建立了Boussinesq型水波方程数值模型，并模拟了孤立波穿越渗透多孔直立式防波堤和波浪在渗透地形下传播。通过以上研究，得到结果如下：

(1)Boussinesq型水波方程中的参数标准不同时，将给出不同的参数值；其中以复波数比为标准下，得到的参数与文献[7]一致。

(2)本文推导的Boussinesq型水波方程可以看成是对文献[7]的拓展，相当于在原方程中多考虑多孔拖曳阻力和惯性阻力。

(3)数值结果与实验结果吻合良好，说明本文建立的Boussinesq型水波方程是有效的。