孤立波与带有窄缝结构作用的数值模拟研究

苏晓杰，宁德志，滕斌

(大连理工大学海岸和近海工程国家重点实验室，辽宁大连116024)



孤立波与带有窄缝结构作用的数值模拟研究

苏晓杰，宁德志，滕斌

(大连理工大学海岸和近海工程国家重点实验室，辽宁大连116024)

摘要:针对孤立波式极值波浪与带有窄缝两结构相互作用问题，采用活塞推板造波法以及时域高阶边界元方法建立了自由水面满足完全非线性边界条件的二维数值波浪水槽模型。求解中采用混合欧拉－拉格朗日方法追踪流体瞬时水面，运用四阶龙格库塔方法更新下一时间步的波面和速度势，利用加速度势的方法来求得作用结构上的瞬时波浪荷载。通过与已发表的数值与实验结果对比，验证了所建模型的准确性。同时通过大量的数值计算研究了窄缝及结构尺寸对迎浪侧波浪爬高，窄缝内波面升高，背浪侧透射浪高以及箱体所受波浪荷载的影响规律。

关键词:窄缝;数值波浪水槽;孤立波;波浪力;高阶边界元

随着海洋工程的不断发展，大型海上结构物被大量使用，这些结构物多由多模块组成，模块之间存在尺度很小的窄缝［1］。多体船以及多船进行联合作业时，船体之间也会有相比船体尺度很小的窄缝出现。在某些频率波浪的作用下，共振现象随之发生，窄缝内发生很大的波面升高，对结构物产生很大的波浪荷载。很多学者研究了波浪与带窄缝浮体的相互作用，但大多研究局限于一般的规则正弦波。Miao等［1］采用渐近匹配法研究了带狭缝二维双箱的共振现象，指出了共振频率与方箱的吃水深度和狭缝宽度的关系。滕斌等［2］采用比例边界有限元方法研究了规则波浪与两箱结构作用下箱体的受力情况，发现随着箱体宽度的减小，共振频率向高频移动。Saitoh等［3］对不同入射波浪作用下两个方箱之间窄缝的波高变化进行了试验研究，发现窄缝内最大共振波高可以达到入射波高的五倍。Zhao等［4］采用势流理论研究了FLNG船与LNG船并排联合作业，发现窄缝内水体出现共振现象时，对船体的运动响应及荷载有很大的影响。Zhang等［5］采用频域方法，研究了两层流中双箱窄缝共振问题，发现共振现象受结构尺寸、流体相对密度的影响。

海洋环境极其复杂，经常出现破坏性更大的极值波浪。沿海地区经常出现的海啸就是极值波浪的一种。极值波浪一般含有惊人的能量，破坏力极大，对海洋结构物的安全产生巨大的危害。为了方便有效地研究极值波浪，可以将其近似为只有单一波峰，周期无限大的孤立波。

过去的几十年里，国内外学者做了很多关于孤立波的研究。Su等［6］给出了孤立波在直墙前爬高的三阶理论公式。Katell等［7］根据Rayleigh的孤立波理论解，推导出一种新的造孤立波的方法。齐鹏等［8］应用湍流数学模型和流体体积法模拟了孤立波翻越防波堤的流动和自由表面变化。刘长根等［9］用基于雷诺方程的二维数学模型研究孤立波与距离水面不同位置处的水平圆柱体的相互作用过程。He等［10］应用高阶边界元法模拟了孤立波与竖直弹性板相互作用，发现弹性板变形与板的刚度以及顶点处的约束情况有很大关系。房克照等［11］建立了基于高阶Boussinesq水波方程的波浪传播数学模型，模拟孤立波在潜礁地形上的传播。Chen等［12］通过模型实验发现当无量纲波高H/d＞0.5时，入射浪的非线性作用增强，并且出现Rayleigh－Taylor不稳定现象。

关于孤立波与带窄缝结构的相互作用问题的研究成果少见发表，因此本文通过采用高阶边界元方法建立自由水面满足完全非线性边界条件的时域数值波浪水槽，采用推板造波方法造波并在水槽另一侧布置阻尼层吸收出流波浪，进而求解孤立波与具有窄缝两结构的相互作用问题。进一步通过大量数值计算研究窄缝及结构尺寸对反射波高、透射波高、窄缝内波高和结构所受波浪荷载等的影响规律。

1　数学模型

考虑孤立波与具有窄缝的两方箱作用布置如图1所示，建立笛卡尔坐标系oxz，坐标原点位于静水面上，且z轴向上为正，波浪沿x轴正向传播。

图1　水槽示意图Fig.1 Definition sketch of the wave flume

图中:h为水槽静水深，W为箱体宽度，D为箱体吃水深度，Wg为两箱体间窄缝的宽度。在箱体1迎浪侧点A和箱体2背浪侧点B，窄缝中间位置水面点G分别布置测点记录波面波动历程。假定流体无粘并且流动无旋，这样整个流域可用速度势来描述。上述问题的控制方程为由速度势满足的拉普拉斯方程［13］，即:

在自由水面上，满足完全非线性动力学和运动学边界条件，在水槽一侧添加人工阻尼层来吸收向右传播的波浪，自由水面边界条件添加阻尼项后可以写成以下形式:式中粘性系数为:

式中:g是重力加速度，x0为阻尼层起始位置坐标，L为阻尼层长度。T为造波板运动时间，表达式如下［10］:

根据Rayleigh－Boussinesq理论，孤立波的波面升高公式如下［10］:

式中，K为outskirts衰减系数，c为孤立波波速，表达式如下:

在水槽的入射边界，造波板的运动方程以及在各瞬时位置处的运动速度可写为［10］:

在水底以及物面上，满足固壁不可渗透条件，即

假定初始时自由水面是静止的，即

在整个流域内对速度势应用格林第二定理，可得到如下边界积分方程［13］:

式中:p=(x0，z0)为源点，q=(x，z)为场点，C为固角系数，G是简单格林函数，考虑到水底镜像，可以表示为如下形式:

式中:r1为p和q两点距离，r2为p和q关于水底镜像之间距离。

求解作用在结构上的波浪力F={ fx，fz}可通过在瞬时物体湿表面上做压强积分得到:

式中:速度势时间导数通过加速度势方法求解可以得到稳定且准确的解。也即，速度势时间导数φt满足Laplace方程:

在自由水面上，φt由Bernoulli方程给出:

在固定边界上满足:

进而通过求解积分方程:

可以求得φt，其中系数矩阵与式(10)中相同，不用重新建立。最后通过式(12)求得作用在物体上的波浪力。

本文用三节点高阶边界元离散计算域成一些曲线单元，单元内任一点的几何坐标和速度势等物理量可以用二次形状函数插值得到。积分方程经高阶边界元离散后，可通过求解线性方程组得到未知量。计算中认为当前时刻物面上的速度势法向导数和自由水面上的速度势是已知的，根据积分方程计算当前时刻物面上的速度势和自由水面上的速度势法向导数，然后应用四阶Runga－Kutta法，根据自由水面条件式计算下一时刻的水质点位置和自由水面上的速度势，再用二次形状函数在旧单元上插值求得新节点上的物理量来对自由水面网格重新划分，重新应用积分方程计算下一时刻物面上的速度势和自由水面上的速度势法向导数。这样计算周而复始，直到计算结束［14］。

2　数值计算及讨论

2.1模型准确性和稳定性

首先验证本数值模型生成孤立波的功能。在数值模型中，水槽水深设定为h=0.5 m，孤立波波高H=0.012 m。计算域长度取70 m，在水槽的右端布置直墙，造波板位于x=0处。通过开展数值收敛性实验，自由水面上共布置140个单元，直墙上布置20个单元，造波板上布置20个单元。时间步长，共模拟200 s。

图2给出了几个不同时刻，孤立波的传播历程图。从图中可以看出在长时间传播过程中，孤立波波形保持不变，没有衰减，说明本模型造出的孤立波可以长时间传播且有很好的稳定性。

为了验证本文模型模拟孤立波与结构物的相互作用问题，图3给出了孤立波在直墙前的无量纲爬高Hmax/h与无量纲入射波高H/h的关系。并且给出了两种孤立波与直墙相互作用的模型实验的测量结果(Chen等［12］，Maxworthy等［15］)，一种数值模拟结果(Cooker等［16］)，以及一种孤立波在直墙前爬高的三阶理论公式(Su等［6］)，公式表示为:

式中，Hmax表示孤立波在直墙前的最大爬高，即波峰距离静水面的高度。从图中可以看出本文结果与以上四种结果稳合的很好，说明本模型可以准确模拟孤立波爬高问题。

图2　孤立波传播历程图Fig.2 Time series of solitary wave

图3　波浪爬高与无量高入射波高的关系Fig.3 Normalized maximum runup Hmax/h against H/h

图4　作用在直墙上波浪荷载随无量纲波高的分布Fig.4 Distribution of dimensionless wave forces on the wall against H/h

图4给出了孤立波作用于直墙的无量纲波浪力Fx/ρgh2(其中F为静水压力与波浪力之和)与无量纲入射波高H/h的关系，同时给出了Maiti等［17］的数值模型模拟的结果以及Fenton等［18］给出的孤立波作用于直墙水平方向波浪力理论公式对比。理论公式表示为:

从图中可以看出，当孤立波高较小时，本文模拟结果与理论公式［18］以及Maiti的数值结果［17］都吻合很好;当波高增大，非线性作用增强时，本文结果相比于Maiti的数值结果［17］更加接近Fenton的理论公式结果［18］。

利用势流理论，国内外学者采用时域［19］和频域［5］的方法模拟了规则波与带窄缝箱体相互作用问题，与物理模型实验［3］结果拟合的很好。通过这些对比研究，说明本文建立的模型可以准确模拟孤立波与直墙以及带窄缝结构物相互作用问题。

2.2数值结果

本文利用所建数学模型模拟研究孤立波与带有窄缝的两固定箱体相互作用问题。

首先选取水深h=0.5 m，入射波高分别为H=0.012、0.05、0.1 m，箱体宽度W=0.5 m，箱体吃水D=0.252 m，以及一系列的窄缝宽度(0.04～0.5 m)研究窄缝宽度不同时，水动力的特征以及变化趋势。计算域长度取100 m，在水槽的右端布置长为20 m的阻尼层，造波板位于x=0，箱体1侧面边界位于距离造波板70 m的位置，然后依次按Wg调整箱体2的位置。自由水面上共布置200个单元，窄缝宽度Wg=0.04 m时布置4个单元，随着窄缝宽度增大，单元数相应增加。箱体侧面边界均布置6个单元，底面边界布置12个单元，每个算例模拟80 s。

图5　窄缝内波浪爬高与窄缝宽度的关系Fig.5 Normalized maximum runup Hmax/h against Wg

图5给出了窄缝内无因次波高与窄缝宽度的关系。从图中可以看出，孤立波与带窄缝两箱体作用时，窄缝内并不会发生很明显的共振现象，与在某些频率规则波作用下会发生共振，产生很大的波浪升高的现象不同。这是因为孤立波的周期与窄缝内水体的共振频率相差较大，不能发生共振现象;另外一个原因是孤立波只有一个峰值，能量难以在窄缝内聚集。随着入射波高的增大，窄缝内的波面升高也增大。当Wg=0.05 m时，3种入射波高对应的窄缝波面升高分别为:H/H0=1.02，1.14和1.25。这是由于波高增大，非线性作用增强导致的结果。

图6给出了上述工况下，H=0.012 m，Wg=0.05 m时，不同位置测点处波面的时间历程。从图中可以看出，23 s时，孤立波峰到达结构物前40 m的位置，40 s到达结构物前A点处，通过结构物反射的波浪以一个类似单一周期的正弦波向远处传播，并在58 s时到达结构物前40 m的位置。孤立波在窄缝处有小幅度的爬高，随后水面有一个小的谷值，结构迎浪侧和背浪侧也有类似的现象，由于孤立波相当于长波，结构物对波形影响较小，因此A、G、B 3点处波形相差不大。

图6　不同位置测点处波面时间历程Fig.6 Time series of wave elevation

图7　波浪爬高与箱体宽度的关系Fig.7 Normalized maximum runup Hmax/h against Wg

图7给出了水深h=0.5 m，入射波高H=0.012 m，窄缝宽度Wg=0.05 m，箱体吃水D=0.252 m时，不同位置无因次波高与箱体宽度的关系。从图中可以发现，随着箱体宽度的增大，窄缝内的波面升高基本不变;而孤立波在结构物迎浪侧的爬高增大，透射浪逐渐减小。这是由于随着箱体宽度增大孤立波在箱体上的爬高增大，而透过箱体的波浪减小。随着入射波高的增大，迎浪测与窄缝内的波面升高增大，背浪侧的透射浪减小，这符合能量守恒的原理。

图8　箱体受力与箱体宽度的关系Fig.8 Dimensionless wave force against W

从图8给出的箱体受力与箱体宽度的关系图中可以看出，随着箱体宽度的增大，水平力和垂向力都有增大的趋势，但垂向力增大的程度更明显。这是由于随着箱体宽度增大，箱体1对孤立波的反射作用增强，透射浪减小，因此箱体2上的垂向力小于箱体1上的垂向力。随着入射波高的增大，两箱体所受水平力增大;箱体1所受垂向力基本不变;箱体2所受垂向力减小。

图9给出了在上述工况箱体宽度W=0.5 m时作用在两箱体上的无量纲波浪力时间历程。从图9(a)给出的水平力时间历程可以看出，两箱体水平力基本同时达到最大值，且相差不大。图9(b)给出的垂向力时间历程相比于水平力时间历程，差距较大，作用于箱体2上的垂向力达到最大值滞后箱体1上的垂向力，但由于箱体宽度相比于孤立波波长尺度较小，因此差距较小。由于透射到箱体2背浪侧的波浪小于箱体1迎浪侧的波浪，导致箱体1上垂向力大于箱体2上垂向力。

图9　箱体上波浪力时间历程Fig.9 The history of wave forces on two boxes

3　结论

本文基于时域高阶边界元方法建立孤立波与具有窄缝的两箱体结构相作用的完全非线性数值水槽模型，对窄缝内流体共振条件下反射波高、透射波高、窄缝内波高、作用在箱体上的波浪荷载等进行了模拟研究。通过与已发表实验数据和数值结果进行对比验证，表明本文所建立数学模型可以准确模拟孤立波与具有窄缝结构相作用过程。通过大量模拟发现:

1)相比于规则波，孤立波与带窄缝箱体相互作用，在窄缝内不会发生明显的共振现象。

2)随着箱体宽度的增大，箱体对孤立波的反射作用增强，透射浪逐渐减小;箱体上的水平力增长幅度较小，垂向力增长程度较大。

3)由于部分入射波被迎浪测箱体1反射，并且透过箱体1的波浪被背浪侧箱体2再次反射，导致作用在箱体二上的垂向力小于箱体1。

4)随着入射波高的增大，孤立波在结构物前的爬高增大，导致透过两箱体的透射浪波高减小，窄缝内的波面升高随入射波高增大而增大。

参考文献:

［1］MIAO Guoping，ISHIDA H，SAITOH T.Influence of gaps between multiple floating bodies on wave forces［J］.China Ocean Engineering，2000，14(4):407－422.

［2］滕斌，何广华，李博宁，等.应用比例边界有限元法求解狭缝对双箱水动力的影响［J］.海洋工程，2006，24(2):29－37.TENG Bin，HE Guanghua，LI Boning，et al.Research on the hydrodynamic influence from the gap between twin caissons by a scaled boundary finite element method［J］.The O－cean Engineering，2006，24(2):29－37.

［3］SAITOH T，MIAO G P，ISHIDA H.Theoretical analysis on appearance condition of fluid resonance in a narrow gap between two modules of very large floating structure［C］//Proceedings of the Third Asia－Pacific Workshop on Marine Hydrodynamics.Shanghai，China，2006:170－175.

［4］ZHAO Wenhua，YANG Jianmin，HU Zhiqiang.Hydrodynamic interaction between FLNG vessel and LNG carrier in side by side configuration［J］.Journal of Hydrodymics，Ser.B，2012，24(5):648－657.

［5］ZHANG Xinshu，BANDYK P.On two－dimensional moonpool resonance for twin bodies in a two－layer fluid［J］.Applied Ocean Research，2013，40:1－13.

［6］SU C H，MIRIE R M.On head－on collision between two solitary－waves［J］.Journal of Fluid Mechanics，1980，98:509－525.

［7］KATELL G，ERIC B.Accuracy of solitary wave generation by a piston wave maker［J］.Journal of Hydraulic Research，2002，40(3):321－331.

［8］齐鹏，王永学，侯一筠.孤立波翻越防波堤流动的湍流数值模拟［J］.水动力学研究与进展，2004，19(S):884－889.QI Peng，WANG Yongxue，HOU Yijun.Numerical simulation of solitary waves overtopping a breakwater［J］.Journal of Hydrodynamics，2004，19(S):884－889.

［9］刘长根，陶建华.孤立波与不同淹没深度水平圆柱体的相互作用［J］.水动力学研究与进展，2006，21(2):155－160.LIU Changgen，TAO Jianhua.Solitary waves interaction with horizontal cylinder submerged in different depth［J］.Journal of Hydrodynamics，2006，21(2):155－160.

［10］HE Guanghua，KASHIWAGI M.Numerical analysis of the hydroelastic behavior of a vertical plate due to solitary waves［J］.Journal of Marine Science and Technology，2012，17(2):154－167.

［11］房克照，刘忠波，唐军，等.潜礁上孤立波传播的数值模拟［J］.哈尔滨工程大学学报，2014，35(3):295－300.FANG Kezhao，LIU Zhongbo，TANG Jun，et al.Simulation of solitary wave transformation over reef profile［J］.Journal of Harbin Engineering University，2014，35(3):295－300.

［12］CHEN Y Y，KHARIF C，YANG J H，et al.An experimental study of steep solitary wave reflection at a vertical wall［J］.European Journal of Mechanics－B/Fluids，2015，49(Part A):20－28.

［13］周斌珍，宁德志，滕斌.造波板运动造波的实时模拟［J］.水动力学研究与进展，2009，24(4):1－12.ZHOU Binzhen，NING Dezhi，TENG Bin.Real－time simulation of waves generated by a wave maker［J］.Journal of Hydrodynamics，2009，24(4):1－12.

［14］宁德志，卓晓玲，滕斌，等.水流对聚焦波特性影响的数值模拟［J］.哈尔滨工程大学学报，2013，34(5):555－561.NING Dezhi，ZHUO Xiaoling，TENG Bin，et al.Numerical investigation of the current influence on the characteristics of focused waves［J］.Journal of Harbin Engineering U－niversity，2013，34(5):555－561.

［15］MAXWORTHY T.Experiments on collisions between solitary waves［J］.Journal of Fluid Mechanics，1976，76(1):177－185.

［16］COOKER M J，WEIDMAN P D，BALE D S.Reflection of a high－amplitude solitary wave at a vertical wall［J］.Journal of Fluid Mechanics，1997，342:141－158.

［17］MAITI S，SEN D.Computation of solitary waves during propagation and runup on a slope［J］.Ocean Engineering，1999，26(11):1063－1083.

［18］FENTON J D，RIENECKER M M.A Fourier method for solving nonlinear water－wave problems:application to solitary－wave interactions［J］.Journal of Fluid Mechanics，1982，118:411－443.

［19］NING D Z，SU X J，TENG B.Numerical study of wave interaction with two bodies in close proximity［C］//Proceedings of the 11th International Conference on Hydrodynamics，Nanyang Technological University，Singapore，2014:1－8.

Numerical simulation of the interaction of solitary waves and structures with narrow slits

SU Xiaojie，NING Dezhi，TENG Bin
(State Key Laboratory of Coastal and Offshore Engineering，Dalian University of Technology，Dalian 116024，China)

Abstract:This study investigated the interaction of extreme solitary waves and structures with narrow slits using a push plate wave－making method and a time－domain higher－order boundary element method.We established a twodimensional numerical wave flume model whose free water surface satisfied fully nonlinear boundary conditions.A mixed Eulerian－Lagrangian method was adopted to track the transient water surface，and the fourth－order Runge－Kutta method was used to refresh the wave surface and velocity potential at each time step.The acceleration potential method was used to calculate the transient wave loads on the structure.The superior accuracy of the proposed model in comparison with published experimental and numerical results was confirmed.Numerical simulations were performed to study the effect of the narrow slit and structure dimensions on wave runup at the head wave side，wave runup inside the narrow slit，wave height at the back wave side，and wave loads on the box.

Keywords:narrow slit;numerical wave flume;solitary wave;wave force;higher－order boundary element

通信作者:宁德志，E－mail:dzning@ dlut.edu.cn.

作者简介:苏晓杰(1989－)，男，硕士研究生;宁德志(1975－)，男，教授，博士生导师.

基金项目:国家自然科学基金资助项目(51222902，51490672和51209107);教育部新世纪优秀人才支持计划(NCET－13－0076).

收稿日期:2014－11－06.网络出版时间:2015－12－21.

中图分类号:O353.2

文献标志码:A

文章编号:1006－7043(2016)01－0086－06

doi:10.11990/jheu.201411027

网络出版地址:http://www.cnki.net/kcms/detail/23.1390.u.20151221.1603.032.html