SPH 方法的孤立波与部分淹没结构物相互作用数值计算研究

林金波，毛鸿飞*，田正林，纪然

(1.广东海洋大学 海洋工程学院，广东 湛江 524088；2.广东海洋大学 海运学院，广东 湛江 524088)

1 引言

波浪与结构物相互作用涉及海洋与海岸结构物的安全生产和灾害防护问题，在海洋工程领域受到广泛关注[1-3]。在风暴潮、台风浪或海啸等极端海况下，结构物周围流场和所受波浪力特征的研究尤其得到重视。风暴潮或海啸引发的波浪波长极长，有时可达上百千米，且当类似极端波浪从深海传播到近岸后，波高急剧增大，会形成波高几十米的拍岸巨浪，同时波浪的传播速度也极快，具备极大的能量，容易导致海洋平台和近岸结构物的破坏失效。在垂向尺度上能够对极端波况进行缩比的孤立波具备常速传播且波形不变的性质，此性质可为波浪与结构物相互作用的机理性研究提供很大便利。孤立波能够在波高较大时表现出极端波况的强非线性，因而在开展波浪与海上结构相互作用的研究中，孤立波常被用于替代极端波况[4-6]。因此，利用孤立波开展极端波况条件下的波浪与结构物相互作用的研究对海洋结构物和近海工程设计与安全运行具有重要的工程指导意义。

目前，国内外关于孤立波与结构物相互作用的物理模型实验已取得一定成果。刚傲等[7]利用物理模型试验，采用粒子图像测速（Particle Image Velocimetry，PIV）技术及小波变换涡结构识别方法，获得了下潜平板迎浪侧和背浪侧的涡结构，从而对孤立波与下潜水平板相互作用的水平板周围流场演化特征进行了研究分析。Arabi 等[8]结合物理模型试验与FLOW3D 计算结果，考虑非破碎波浪，对孤立波与矩形截面不透水结构相互作用进行研究分析，总结了孤立波与矩形障碍物相互作用过程中湍流拟序结构的生成和演化规律。王千等[9]在物理模型试验中引入多目视觉立体重构技术实现三维波面测量，并利用三分力测力传感器测量孤立波对淹没平板的作用力和力矩，得到了俯仰力矩和垂向力极值点的三维波面形状。邵奇等[10]利用物理模型试验，进行孤立波作用下淹没圆盘波浪荷载实验，针对非破碎波面工况，分析孤立波经过淹没圆盘时的垂向力、水平力和倾覆力矩的时间变化，对圆盘上下表面的压力分布特征及其对波浪荷载的影响规律进行了阐述。

随着计算机技术的快速发展，研究人员利用数值模拟方法开展了大量孤立波与结构物的相互作用研究。宋帅[11]基于粘性流理论，建立二维和三维数值波浪水槽，对孤立波与多孔介质及外壁透空的圆柱浮体相互作用进行了数值模拟研究，分析了孤立波通过多孔介质后的波面变化情况以及外壁透空圆柱结构内圆柱体上的水动力特性。Mo 和Liu[12]利用有限体积方法离散Navier-Stokes（N-S）方程，并结合流体体积（Volume of Fluid，VOF）自由表面追踪方法，对孤立波与细长圆柱的相互作用进行了数值模拟研究。Aghili 等[13]基于弱可压光滑粒子流体动力学（Smoothed Particle Hydrodynamics，SPH）模型，对孤立波与不同垂向位置的水平淹没板的相互作用进行数值模拟，研究分析了淹没深度对自由表面、压力和波浪力垂直分量的影响。高俊亮[14]釆用FUNWAVE2.0 数值模型对孤立波或波群诱发的港湾共振现象进行了讨论分析。Hu 等[15]在OpenFOAM 中引入新的波浪边界条件，对孤立波与固定和浮式截断圆柱以及简化的浮式采油平台的相互作用进行了数值模拟研究。Xu 等[16]利用人工神经网络方法建立一种新的预测海岸桥面板上孤立波作用力的模型，对孤立波作用下的桥面板波浪荷载进行了计算研究。耿添和刘桦[17]采用边界元方法建立二维数值模型，对相对波高、相对板长和相对浸没深度对浸没平板水动力荷载的影响规律进行了研究分析。邹丽等[18]根据势流理论，应用边界元法对非线性孤立波的传播和直墙反射进行模拟研究，分析了孤立波作用于直墙的滞留时间、最大爬高和直墙瞬时压力。殷铭简等[19]基于开源代码Open-FOAM，开展了孤立波与密排桩防波堤相互作用研究，对密排桩防波堤的消浪机理和波浪爬升进行了重点分析。许彦章和万德成[20]体运动，采用FEM 模型模拟垂直板的弹性变形，对孤立波与水下垂直刚性板和弹性板的相互作用进行数值模拟研究，讨论了弹性模量对垂直板与孤立波相互作用的影响。虽然，国内外相关研究较多，但多是关于孤立波与板式或圆柱式结构物的相互作用研究，且由于波浪与浮式结构物相互作用过程中剧烈自由面变化（例如波浪的翻卷和破碎等）的存在，使得相关研究多考虑的是淹没式结构的非破碎波况。

由于SPH 方法的无网格特性，使得SPH 模型在处理大变形和自由表面破碎上具有独特优势，因而近年来被越来越多的研究人员应用于波浪与结构相互作用的研究。Pan 等[21]应用SPH 模型对孤立波与浮式海洋结构物的相互作用进行模拟研究，分析了张力腿平台在孤立波冲击下的漂移运动。温鸿杰等[22]利用SPH 模型对波浪在珊瑚礁上的传播演变进行了数值模拟研究。He 等[23]利用SPH 模型建立数值波流水槽，对波流相互作用进行了数值模拟研究。贺铭等[24]基于冲击式造波理论，对关联楔体的浸没深度、楔角及生成波高的隐式方程组进行推导，并建立了楔体运动和波浪生成的SPH 孤立波数值模型，利用试验数据对方程组和模型的可靠性进行了验证。虽然研究人员基于SPH 模型开展的波浪与结构相互作用研究较多，但基于SPH 模型的孤立波与结构相互作用研究尚不完善。

本文利用开源代码DualSPHysics，基于无网格SPH 方法，结合Rayleigh 孤立波生成理论，建立能够处理大变形自由表面以及强非线性现象的孤立波与海洋结构相互作用无网格数值模型，并通过将模型计算自由表面与解析解进行对比验证模型精度。利用验证后的数值模型对不同相对波高的孤立波与部分淹没矩形结构物相互作用过程进行数值模拟研究，分析相对波高对结构物周围波面、流速、涡量及结构荷载的影响，总结相对波高对孤立波与部分淹没矩形结构物相互作用流场的影响规律。

2 数值模型

2.1 控制方程

黏性流体控制方程为连续性方程和动量方程：

式中，ρ代表流体密度；u代表速度；P代表压力；g=(0,0,-9.81) m/s2，代表重力加速度；Г表示黏性项。

应用SPH 方法离散以上两式，则拉格朗日形式N-S 方程[25-26]表示为：

式中，m代表粒子质量；uij=ui-uj，下标i代表插值点而下标j代表插值点周围的邻近粒子；Wij=W(rij,hs)代表核函数；rij=ri-rj表示粒子间距。

在连续性方程中引入一个修正项[27]以减少方程中的密度波动，可得到delta-SPH 方程：

模型光滑核函数采用能够提供高阶插值特性，并维持中等计算量的五次核函数[28-29]。

式中，r、r'分别表示计算粒子及其邻近粒子位置矢量；R=r-r'表示计算粒子与其邻近粒子间距离；hs代表光滑长度，后文计算中光滑长度取值为2；二维模型中 αd取值为。

2.2 黏性处理方法

动量方程中的黏性项有几种不同的计算方式，其中人工黏性由于形式简单且能够防止粒子相互接近时的非物理穿透，因此在SPH 方法中被广泛应用。

人工黏性项形式如下[30-31]：

式中，vij=vi-vj，rij=ri-rj，vi和ri分别代表粒子的速度和位置向量;α代表人工黏性系数，取值为0.01；=0.5(ci+cj)代表平均数值声速。

2.3 状态方程

为了保持SPH 方法显式特征，提高计算效率，引入状态方程代替求解压力泊松方程，从而基于粒子密度值来计算压力值。Tait 状态方程[32]可表示为

2.4 孤立波造波

模型中，孤立波的生成是基于Rayleigh 方法[34]，利用造波板来制造孤立波。孤立波造波方法的基本假设是波峰下水粒子的平均水平速度与造波板速度相同。

式中，xs为造波板位移；u(xs,t)为水粒子水深平均速度，可表达为

式中，cw代表波速；d为水深；η代表自由表面高程。

将式（10）与式（11）合并，积分可得造波板位移公式：

式中，k为边缘系数，其描述了自由表面高程在无穷远处趋于平均水面的方式；H为波高；c为波速。孤立波分布可以表达为

式（13）为隐式方程，其存在几种不同的求解方式，采用Rayleigh 方法进行求解，孤立波传播时振幅损失较小[35]，上式表达的理论自由表面高程可以重写为

模型中固壁边界采用动力边界方法[36]处理，动力边界法无需进行显式的边界处理，而是将边界处理耦合在了控制方程的求解过程中，因此代码实现简单的同时，计算量也相对较小，适用于具有复杂边界的问题。方程求解采用辛（Symplectic）数值积分方法[37]，该方法在没有摩擦或黏性影响的情况下是时间可逆的，且具有显式二阶精度O(Δt2)，时间步长采用可变时间步长方法。

3 模型验证

为了验证模型计算精度，利用数值模型对初始水深d0=0.143 m，相对波高H0/d0为0.1 和0.4 的孤立波在水槽中的生成与传播进行数值模拟，将数值计算自由表面结果与理论解进行对比。其中，孤立波波面的理论解通过Green-Naghdi（G-N）方程[4]进行计算得到。数值水槽长度为12 m，高0.2 m。粒子间距Δx/d0=0.014，输出时间间隔0.02 s，计算时间共4 s。

观察模型计算孤立波波面形态，孤立波沿水平方向传播稳定，孤立波波峰通过后无明显尾波，波面及自由表面均较为清晰，并且孤立波剖面全部在静水面以上，波长近似为无，计算结果符合孤立波的实际物理特性。为了进一步分析模型计算精度，图1a 和图1b 分别为相对波高0.1 和0.4 时，x=2 m 位置数值计算结果与理论解对比，图中纵坐标d/d0为相对水深。

图1 模型计算结果与理论解对比Fig.1 Comparison between calculated results and analytical solution

相对波高为0.1 时，模型计算结果的相位、波峰与理论解基本一致，计算结果与理论解吻合良好。相对波高为0.4 时，模型计算波高与理论解基本一致，但相位稍有差异，且孤立波通过后测点处水位稍有波动，但总体来说误差不大。为了定量分析模型计算精度，表1 给出了模型计算结果与理论解之间的L2误差，相对波高为0.1 和0.4 时，L2误差分别为0.005 和0.014，虽然随相对波高增大，误差增大，但整体误差较小。

表1 模型计算结果L2 误差Table 1 L2 error of the calculated results

式中，N为数据总量；表示位置i处的数值计算结果；表示位置i处的理论解。

从上述结果可知，模型计算结果与理论解具有较好的一致性，数值模型能够比较理想地模拟孤立波的生成与传播过程。

4 孤立波与部分淹没矩形结构物相互作用研究

4.1 模型布置

为了分析孤立波波高对部分淹没固定矩形结构物相互作用的自由表面高程、流速、涡量及水动力荷载的影响，利用长度为100 m 的矩形数值水槽[38]，对孤立波与结构物相互作用过程进行数值模拟计算。矩形水槽高度为2 m，静水深d0=1.0 m。一个尺寸为5.0 m× 0.6 m 的矩形结构固定在静水面以上，结构物中心坐标为（52.5 m，0.9 m），初始淹没高度0.4 m，模型布置见图2。相对粒子间距Δx/d0取0.01，相对波高H0/d0分别取0.025、0.05、0.1、0.15、0.2、0.3 和0.5。

图2 模型布置示意图Fig.2 Layout of the model

图3 展示了相对波高为0.1 的模型计算测点水位及结构物受力计算结果。图3a 至图3c 测点水位和结构物x方向受力计算结果较为光滑，而图3d 结构物z方向受力计算结果初始时刻存在较大震荡，12 s后z方向受力计算结果趋于稳定。因此，孤立波与部分淹没固定矩形结构物相互作用计算中设置初始稳定时间为15 s，15 s 后造波板开始运动产生孤立波，实际计算时间40 s，总计算时间为55 s。

图3 相对波高为0.1 模型计算测点水位及结构物受力Fig.3 Water level at the measurement point and force on the obstacle with relative wave height of 0.1

4.2 水位及水动力荷载计算结果分析

图4 给出了不同相对波高的孤立波与结构物相互作用过程中，结构物上、下游测点的无量纲自由表面高程变化情况。

图4 不同相对波高测点相对自由表面高程变化Fig.4 History of the relative free surface elevation with different relative wave height

图4 中，随着相对波高增大，波速逐渐增大，孤立波到达测点时间逐渐向前移动，特别是当相对波高为0.5 时，37 s 后均有反射波到达两测点位置，而其他相对波高条件下反射波未到达测点位置，水位无抬升。相对波高较小时水面变化较为光滑，而随着相对波高增大，当相对波高大于0.2 后，孤立波波峰部分爬升至结构物顶部，一部分水体继续传播进入下游与下游水体碰撞，而部分水体从结构物上游掉落入水体，都会引起水面发生震荡，并随着相对波高增大振幅逐渐增大。图4a 中能够明显观察到两个衰减的波峰，第一个波峰为右行孤立波到达结构物上游测点引起水位抬升从而形成孤立波，第二个波峰伴随一个波谷形成完整周期的波浪是由于入射孤立波与结构物作用后形成左行反射波到达测点位置形成的完整波浪；而图4b 中右行孤立波与结构物作用波峰降低后继续向下游传播，直到测点位置形成孤立波并继续向下游传播。

图5 为结构物中间位置x=52.5 m 处测点的无量纲自由表面高程变化情况。由于该测点位于结构物中央，结构物部分淹没在水中，所以初始时刻水面高程为结构物底高程，相对自由表面高程为0.6；相对波高小于0.2 时，由于波高较小，并未发生越浪，结构物顶部无过流，故测点相对自由表面高程均为0.6；当相对波高大于0.2 后，结构物顶部发生越浪，t＞15 s 后，测点相对自由表面高程大于1.2，超过结构物顶高程。因此，当结构物未淹没比（高出水面的高度与总高的比值）为0.33 时，相对波高大于0.2，结构物上会发生越浪。

图5 不同相对波高x=52.5 m 测点相对自由表面高程变化Fig.5 History of the relative free surface elevation with different relative wave height at x=52.5 m

图6 中给出了矩形结构物水平和垂直方向所受的水动力荷载系数Cx和Cz随时间变化情况。Cx和Cz的计算公式为

图6 不同相对波高结构物水动力荷载系数变化Fig.6 History of hydrodynamic load coefficient of obstacle with different relative wave height

式中，A=0.6 m× 5 m 为矩形结构物面积；Fx、Fz分别为结构物水平和垂直方向所受合力。Fx、Fz的值是首先计算构成结构物的单个粒子的受力，再对构成结构物的全部粒子进行求和，计算公式为

式中，Fa为构成结构物的任意粒子a所受的力；F（Fx、Fz）为结构物所受合力。

与图4 相似，图6 中波速随相对波高增大而增大，结构物开始受力时刻逐渐提前，且结构物水动力荷载随相对波高增大而增大。相对波高较小时结构物水动力荷载变化较为平缓，相对波高大于0.2 后，孤立波波峰部分从顶部越过结构物并与水体发生碰撞破碎，造成水面震荡，引起水动力荷载发生波动，相对波高越大震荡越剧烈，结构物水平和垂向负压也越大。水平和垂向受力均为单波峰，说明结构物上下游均无反射波到达结构物处。

为了分析相对波高与结构物荷载的定量关系，图7给出了结构物水平和垂直方向所受的水动力荷载幅度随相对波高变化情况。图中，结构物所受荷载幅度随着相对波高的增大而逐渐上升，且结构物所受垂向荷载幅度较水平荷载幅度更大，相对波高大于0.3 后，荷载幅度上升速度下降，荷载增加速度变慢。

图7 结构物受力荷载幅度与相对波高关系曲线Fig.7 Relation between structure load and relative wave height

4.3 自由表面形态计算结果分析

图8 为不同相对波高的孤立波在7 个不同时刻数值水槽内的自由表面形态。计算开始后，计算域左侧造波板开始运动，生成向右传播的孤立波。t=14.75 s 时，相对波高为0.5 的孤立波到达结构物处，并爬升至结构物顶部。相对波高为0.5 的孤立波在t=16.25 s 从上部和底部越过结构物，到达结构物下游并发生碰撞破碎，且部分孤立波受到结构物阻碍形成左行反射波，同时相对波高为0.3 的孤立波到达结构物处并爬升至结构物顶部。t=17.50 s 时，相对波高为0.2 的孤立波到达结构物处，同样爬升至结构物顶部；相对波高为0.3 的孤立波从底部传播至结构物下游引起下游水面抬升，而结构物顶部水体尚未越过结构物，上游同样形成了左行反射波；相对波高为0.5 的孤立波结构物上下游波浪均向前传播一定距离，水面震荡较大。t=18.50 s，相对波高为0.15 的孤立波到达结构物处，但并未爬升至结构物顶部，仅从底部通过结构物向下游继续传播；相对波高为0.2 和0.3 的孤立波均从上部和底部越过结构物，到达结构物下游并与下游水体碰撞；相对波高为0.5 孤立波的上下游波浪继续向前传播，自由表面震荡减弱。直到t=21.00 s，相对波高为0.025、0.05 和0.1 的孤立波均已到达结构物处，与结构物作用后继续向下游传播，并形成左行反射波。随后，左行反射波和右行孤立波继续向前传播，t=26.25 s 时，上下游波浪均已远离结构物。t=40.00 s时，上下游波浪与边界作用后反射波重新接近结构物，其中，仅相对波高为0.5 的孤立波反射波浪再次经过x=35 m 和x=60 m 水位测点，引起测点水位抬升。总之，随着相对波高的增大，波速逐渐增大，孤立波传播速度随之增加；相对波高较小时自由表面较为光滑，随着相对波高增加，相对波高大于0.2 后，孤立波波峰分为两个部分，分别从上部和底部越过结构物，结构物上游水体传播进入下游，与下游水体碰撞破碎，而部分水体从结构物前缘掉落进入上游水体，引起结构物附近上下游水面发生震荡，且振幅随着相对波高增大而增大。

图8 不同相对波高瞬时自由表面形态Fig.8 Free surface at instants with different relative wave height

4.4 流速及涡量分布计算结果分析

不同相对波高的孤立波与矩形结构物相互作用流速场及涡量场在图9 中给出，随着相对波高的增大，结构物周围流速逐渐增加，相对波高为0.5 时最大流速达2 m/s 以上；相对波高大于0.2 以后，部分水体从顶部越过结构物，且流速较大，越过结构物顶部的水体流入下游，与下游水体碰撞破碎，使得流场更加复杂。相对波高为0.025 时，结构物周围无涡生成；随着相对波高增大，结构物周围涡逐渐增大且逐渐远离结构物向周围发展；相对波高较小时，由于边界效应，仅结构物壁面处存在微小涡量；相对波高大于0.15 后，结构物前缘底部形成一较小顺时针涡及一较大的逆时针涡，并随着相对波高增大而逐渐增大，且向深度方向发展；结构物下游涡分布较为复杂，相对波高大于0.15 后，壁面附近逆时针涡随着相对波高增大逐渐增大，但并未向下游发展，而结构物顶部水体流入下游与下游水体碰撞形成顺时针及逆时针涡，并随着相对波高增大涡场更加复杂，涡场逐渐向下游和底部发展。

图9 不同相对波高孤立波与结构物相互作用流速与涡量分布Fig.9 Velocity and vorticity of the interaction between solitary wave with different relative wave height and obstacle

5 结论

（1）基于Rayleigh 孤立波生成理论及弱可压SPH方法，建立了孤立波与海洋结构相互作用无网格数值模型。将模型不同相对波高孤立波计算自由表面与解析解进行对比，基于人工黏性和Rayleigh 孤立波生成理论的弱可压delta-SPH 模型计算结果与理论解具有较好的一致性，数值模型能够准确模拟孤立波的生成与传播过程。

（2）利用验证后的数值模型对不同相对波高的孤立波与部分淹没矩形结构物相互作用过程进行数值模拟研究，由于受到边界处理方法的影响，计算初始时刻结构物竖直方向受力存在非物理震荡，需给出初始稳定时间，本文设置初始稳定时间为15 s。

（3）当结构物未淹没比为0.33 时，相对波高大于0.2，结构物上会发生越浪，爬升至结构物顶部的水体一部分向下游传播并与下游水体碰撞，一部分从结构物上游掉落入水体，引起水面发生震荡，水动力荷载同样产生波动，且随着相对波高增大，水面及水动力荷载震荡逐渐增大，结构物水平和垂向负压也越大。结构物所受荷载幅度随着相对波高的增大而逐渐上升，且结构物所受垂向荷载幅度较水平荷载幅度更大，相对波高大于0.3 后，荷载幅度上升速度下降，荷载增加速度变慢。结构物周围流速和涡量均随相对波高增大而增大，特别是当孤立波爬升至结构物顶部，并越过结构物与下游水体发生碰撞后，流场及涡量场变得更加复杂，并且结构物周围涡量随相对波高增大逐渐向深度方向和下游方向发展。