0 引　言

近年来，流固耦合（Fluid-Structure Interaction, FSI）问题在海洋工程领域备受关注，特别是三维浮式结构物与系泊装置耦合系统的动态响应研究。随着计算力学的发展，学者们通过耦合有限元法（FEM）^[1]、边界元法（BEM）^{[2 - 3]}与计算流体动力学（CFD）技术^[4]，在波浪荷载预测、系泊缆非线性动力特性分析等方面取得显著进展。

然而，传统网格方法在处理强非线性自由液面、结构物大变形及复杂接触边界时仍面临网格畸变与重构的挑战。在此背景下，光滑粒子流体动力学（Smoothed Particle Hydrodynamics, SPH）方法因其无网格、拉格朗日特性的优势，在三维流固耦合研究中崭露头角。SPH方法通过离散粒子描述连续介质，天然适应自由表面流动、流体破碎等强非线性现象，并能有效模拟浮体六自由度运动与系泊缆索的柔性变形耦合过程，已成功应用于海上浮式平台^[5]、波浪能装置等复杂系统的动力学分析以及浮体与波浪相互作用^{[5 - 7]}、漂浮结构水动力响应^[8]以及系泊系统耦合建模^{[9 - 10]}等方面。He等^[11]建立SPH数值模型研究波与多浮子结构的相互作用，对波况、间隙间距、模块数等不同因素的影响进行了全面讨论。Salis等^[12]基于SPH方法开发3D数值模型，研究聚焦波对系泊浮动结构的影响，波高误差约为5%，可以对该物理过程进行全面和准确的预测。Liu等^[13]建立DualSPHysics与MoorDyn之间的耦合模型，研究几何参数对系泊浮式防波堤波浪衰减性能的影响，发现将每个单独测试组中的优化参数组合起来可以有效地优化浮式防波堤。Wan等^[14]应用SPH方法研究了混合筏式系统地水动力性能，揭示了波与筏子之间地相互作用和干扰效应。此外，Domínguez等^[15]通过SPH模拟系泊平台，验证其对系泊缆-浮体耦合运动的捕捉能力优于传统频域方法。Jiao等^[16]基于DualSPHysics模拟了并排LNG船舶与液舱晃荡耦合作用，发现液舱装载率和船间间隙对船体运动具有显著影响，验证了该方法在并靠工况下的工程适用性。

尽管SPH方法在波浪砰击以及自由液面破碎变形的模拟中展现出独特优势，但其应用场景的单一性严重影响了算法的发展。为此，本文提出了一种耦合SPH数值模型，采用SPH算法模拟流场中的波浪载荷，基于集中质量法求解系泊缆张力，并通过顺序交错算法，将二者进行耦合求解，同时利用GPU加速技术，显著提高了SPH算法在计算三维问题中的计算效率，并通过模型试验，验证了该耦合模型的数值精度和收敛性，最终成功的模拟了非线性波浪载荷与单点系泊船舶的耦合作用过程，为海洋浮式结构物与系泊系统耦合动力学的精细化研究提供了高效数值工具。通过系统分析短周期波浪下的砰击现象与结构响应变化，揭示其对单点系泊船舶稳定性与极限工况的影响机制，为复杂海况下的系泊系统设计与评估提供理论支持。

1 数值模型 1.1 SPH基本理论

SPH方法用粒子对问题域进行离散化。通过SPH中的核近似和粒子近似方法，具有密度扩散项和动量方程的拉格朗日连续性方程可以离散为：

$\left\{ {\begin{aligned} & {\frac{{{\rm{d}}{\rho _i}}}{{{\rm{d}}t}} = \sum\limits_k {{m_k}} {\nu _{ik}} \cdot {\nabla _i}W\left( {{x_i} - {x_k},h} \right) + } \\ & {0.1h{c_0}\sum\limits_k {{\psi _{ik}}} \cdot {\nabla _i}W\left( {{x_i} - {x_k},h} \right)\frac{{{m_k}}}{{{\rho _k}}}}，\\ & {{\psi _{ik}} = 2\left( {\rho _k^D - \rho _i^D} \right)\frac{{{x_i} - {x_k}}}{{{{\left\| {{x_i} - {x_k}} \right\|}^2}}}} 。\end{aligned}} \right. $

(1)

$ \frac{{{\rm{d}}{\nu _i}}}{{{\rm{d}}t}} = - \sum\limits_k {{m_k}} \left[ {\left( {\frac{{{p_i}}}{{\rho _i^2}} + \frac{{{p_k}}}{{\rho _i^2}}} \right) + {\prod _{ik}}} \right]{\nabla _i}W({x_i} - {x_k},h) + g 。$

(2)

式中：${\rho _i}$和${\rho _k}$分别为流体粒子$i$和$k$的密度；${m_k}$为流体粒子$k$的质量；${v_{ik}}$为粒子$i$和粒子$k$之间的速度矢量差，即$ {\nu _{ik}} = {\nu _i} - {v_k} $；$W$为核函数；$x$为粒子的位置；$h$为$W$的光滑长度，本文取$h = 1.25\Delta x$，其中$\Delta x$为粒子间距；${c_0}$为人工声速；上标$D$表示动密度或动压力；${p_i}$和${p_k}$分别为粒子$i$和粒子$k$处的压强；${\Pi _{ik}}$为人工粘度项；$g$为重力加速度。研究使用的光滑函数为：

$ W(r,h) = \frac{{21{\text{π}} {h^3}}}{{16}}{\left( {1 - \frac{q}{2}} \right)^4}(2q + 1),0 \leqslant q \leqslant 2。$

(3)

式中：$q = r/h$，$r$为支持域内目标粒子与其相邻粒子的距离。

人工粘度${\Pi _{ik}}$用于抑制激波引起的数值非物理振荡，表示为：

$ \Pi _{ik} = \displaystyle\left\{\begin{aligned} &{\displaystyle\frac{{ - \alpha {{\bar c}_{ik}}{\mu _{ik}}}}{{{{\bar \rho }_{ik}}}},{\nu _{ik}} \cdot {x_{ik}} \lt 0}，\\ & {0,\;\;\; {\nu _{ik}} \cdot {x_{ik}} \geqslant 0}。\end{aligned} \right. $

(4)

$\quad{\mu _{ik}} = \frac{{h{\nu _{ik}} \cdot {x_{ik}}}}{{{{\left| {{x_{ik}}} \right|}^2} + {\eta ^2}}} 。$

(5)

式中：$\alpha $为系数；${\bar c_{ik}} = ({c_i} + {c_k})/2$为平均声速；$ {x_{ik}} = {x_i} - {x_k} $；${\eta ^2} = 0.01{h^2}$。

为满足弱可压假设条件，在SPH法计算过程中，其密度的变化率应控制在1%以为，同时本文采用状态方程描述流体压力与密度的关系，表示为：

$ P = b\left[ {{{\left( {\frac{\rho }{{{\rho _0}}}} \right)}^\gamma } - 1} \right];\;\;\; b = c_0^2{\rho _0}/\gamma。$

(6)

式中：ρ₀=1 000 kg/m³为水的参考密度；$\gamma $为流体的多向系数，对于水流通常取7；${c_0} = \sqrt {(\partial P/\partial \rho )} {|_{{\rho _0}}}$为流体的人工声速。

1.2 集中质量法基本理论

对于集中质量法，是将一根连续的系泊缆离散成$ N $段尺寸均匀的线段，为方便叙述，用$ {P_n} $表示张力；$ {F_n} $表示阻尼力；用$ {(n + 1)}{ \mathord{\left/ {\vphantom {1 2}} \right. } 2} $表示某一段系泊缆的编号，并将施加在某个节点$ i $处的系泊缆重量被定义为$ {M_n} $。

对于一段编号为$ {(n + 1)}{ \mathord{\left/ {\vphantom {1 2}} \right. } 2} $表示的系泊缆，其受到的浮力表示为：

$ {W_{n + \displaystyle\frac{1}{2}}} = \frac{\text π}{4}{d^2}l\left( {{\rho _w} - {\rho _{{\rm{moorings}}}}} \right)g。$

(7)

式中：$ d $为水深；$ {\rho _w} $为水密度（本文中，ρ_w=1 000 kg/m³）。由于2个连接的节点在承受着相同大小的浮力作用，所以，每个节点受到的浮力矢量可表示为：

$ {W_n} = \frac{1}{2}\left( {{W_{n + \displaystyle\frac{1}{2}}} + {W_{n - \displaystyle\frac{1}{2}}}} \right)\mathop {{k_z}}\limits^ \wedge。$

(8)

式中：$ \mathop {{k_z}}\limits^ \wedge $为正Z方向的单位向量。

对于作用在$ n + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2} $段系泊缆的张力可表示为：

$ {P_{n + \displaystyle\frac{1}{2}}} = E\frac{\text π}{4}{R^2}\left( {\frac{1}{l} - \frac{1}{{\left\| {{r_{n + 1}} - {r_n}} \right\|}}} \right)\left( {{r_{n + 1}} - {r_n}} \right)。$

(9)

式中：$ l $为这段系泊缆的长度；$ E $为杨氏模量；$ R $为系泊缆的直径；$ {r_n} $为第$ n $处节点位移矢量；$ {r_{n + 1}} $为第$ n + 1 $处节点位移矢量。

通过为每个节点分配2个相邻系泊缆段质量的一半，可以实现以点质量来代替系泊缆的质量。因此，节点$ i $的质量矩阵可以表示为：

$ {m_n} = \frac{\text{π} }{4}{R^2}l{\rho _{\rm moorings}}{\boldsymbol{E}}。$

(10)

式中：$ {\boldsymbol{E}} $为单位矩阵。

节点$ i $的附加质量矩阵可以表示为：

$ {{\boldsymbol{k}}_n} = {{\boldsymbol{k}}_{pn}} + {{\boldsymbol{k}}_{qn}} = \frac{\text{π} }{4}{R^2}{\rho _w}l\left[ {{{\boldsymbol{C}}_{ap}}\left( {{\boldsymbol{E}} - \mathop {{h_n}}\limits^ \wedge {{\mathop {{h_n}}\limits^ \wedge }^{\mathrm{T}}}} \right) + {{\boldsymbol{C}}_{aq}}\left( {\mathop {{h_n}}\limits^ \wedge {{\mathop {{h_n}}\limits^ \wedge }^{\mathrm{T}}}} \right)} \right]。$

(11)

式中：$ {{\boldsymbol{k}}_{pn}} $为横向附加质量矩阵；$ {{\boldsymbol{k}}_{qn}} $为切向附加质量矩阵；$ {{\boldsymbol{C}}_{ap}} $为横向附加质量系数；$ {{\boldsymbol{C}}_{aq}} $为切向附加质量系数。

最终，对于某一节点$ n $，其完整的运动方程可表示为：

$\begin{split} \underbrace {\left[ {{m_n} + {k_n}} \right]\mathop {{r_n}}\limits^{ \cdot \cdot } }_{\rm{mass \; and \; added \; mass }} = &\underbrace {{P_{n + \displaystyle \frac{1}{2}}} - {P_{n - \displaystyle \frac{1}{2}}} + {F_{n +\displaystyle \frac{1}{2}}} - {F_{n -\displaystyle \frac{1}{2}}}}_{\rm {internal \; stiffness \; and \; damping \; force}} + \\ &\underbrace {{W_n} + {B_n}}_{\rm {weigth \; and \; contact}} + \underbrace {{D_{pn}} + {D_{qn}}}_{\rm drag} 。\end{split} $

(12)

式中：$ {B_n} = Rl\left[ {\left( {{z_{bot}} - {z_n}} \right){k_b} - \mathop {{z_n}}\limits^ \cdot {c_b}} \right]\mathop {{k_z}}\limits^ \wedge $；$ {D_{pn}} $为横向拖曳力；$ {D_{qn}} $为切向拖曳力；$ {k_b} $为刚度系数；$ {c_b} $为阻尼系数。

最后，采用常时间步二阶龙格-库塔（RK2）积分算法求解整个系泊系统。

1.3 耦合SPH数值模型

本文采用集中质量法实现对系泊缆张力的计算，其主要思想为它将系泊缆离散化为由线性弹簧阻尼器段连接的点质量^[17]。系泊缆模型结合了内部轴向刚度、阻尼力与重量、浮力等，这些力构成了系泊缆每个节点的运动方程。同时利用顺序交错流固耦合算法将系泊缆数值模型与SPH算法模型分区求解，其中在SPH算法中计算流场载荷，并将载荷作为外力传递至刚体六自由度运动方程，随后将刚体数据传递至系泊缆张力求解器中，更新求解导缆孔处的物理信息，并在下一个流场时间步开始之前，重新将结构信息传递回流场。

2 数值验证 2.1 收敛性分析

考虑到SPH算法在计算效率上的短板，为了兼顾计算成本和计算精度，本文首先选取了5种不同的粒子分辨率进行收敛性分析，其中测试工况为规则波，而计算分辨率$ L/dp $分别取200、250、300、500和1000。入射波高设置为0.15 m，入射波周期设置为1.6 s，图1所示为不同分辨率下测点处波浪高程的时域曲线对比情况。可以看出，当$ L/dp = 300 $时，粒子精度已经满足收敛性要求，随着分辨率继续增大，变化较小，但高分辨率易带来较高的计算成本，例如物理计算时间为20 s，粒子总数1767846，采用GPU加速时，计算设备采用 NVIDIA GeForce RTX 3090 GPU，计算时长为4 h，采用CPU进行并行加速时，计算设备为Intel Core i9-14900KF CPU，所需计算时长约为128 h，因此在后续的模拟中本文采用的计算分辨率为$ L/dp = 300 $，并采用GPU加速完成所有工况计算。

2.2 数值验证

前面章节完成了数值模型的收敛性分析，并确定了最优计算分辨率，接下来本节将针对耦合SPH模型开展数值精度验证。

在本节中，选择了Tan等^[18]波浪作用下的浮式方箱水动力试验。试验布置如图2所示。系泊缆采用304不锈钢制成，弹性模量为193 GPa，密度为7.93 g/m³，线重为0.0305 kg/m，等效直径为0.002213 m。系泊缆长度为0.6 m。系泊点和锚固点布置如图所示。需要注意的是，在本案例计算使用的的粒子分辨率收敛性分析一节中保持一致。

本节选择波高0.1 m，周期1.0 s的规则波作为载荷输入。图3所示为本文建立的三维系泊浮体运动响应与实验结果之间的对比。可以看出，数值模拟的计算结果与实验数据吻合良好。图4所示为典型时刻系泊浮体的运动姿态以及系泊缆运动姿态，同时与实验进行了对比。可以看出，基于SPH方法建立的三维系泊浮体模型可以准确地捕捉到浮体在波浪作用下的运动响应和受力特性，系泊张力求解模型也可以准确地接受浮体运动参数的输入，并提供准确的拖曳力，实现自由浮体的动态锚泊。因此，本文将应用该模型对波浪与系泊船舶之间的耦合效应进行数值求解，并分析波浪引起的砰击效应和船舶的运动特性具有较高的可行性。

3 结果分析 3.1 不同波浪周期

对波浪作用下的三维单点系泊船舶的水动力学特性问题，为提高计算效率，采用耦合SPH-GPU计算模型对其进行数值求解，具体而言，整个计算流程分为3个步骤：1）在程序迭代开始之前，需要设置基本参数，并通过cudaMemcpyToSymbol函数将粒子信息拷贝至GPU内存；2）计算时在GPU端给每一个粒子单独分配一个线程，并计算与所有相邻粒子的相互作用；3）完成计算后将粒子信息从CPU端拷贝至GPU端。

船舶的具体参数如表1所示，系泊缆的参数如表2所示。与模型缩尺类似，系泊缆也需要根据缩尺比1∶10进行等效缩尺。具体数值布置如图5所示，其中水深为2 m，选择的粒子分辨率为$ L/dp = 300 $。

为研究不同波浪条件对系泊船舶运动响应的影响，本文将选择规则波浪作为波浪条件，首先考虑的是不同波浪周期对船舶运动响应的影响。为此，本文选择了5种波浪周期，如表3所示。

图6所示为不同入射波周期下船舶纵荡运动响应的变化。在初始阶段，由于系泊缆预张力的影响，初始位置的船舶将发生横向位移并到达平衡位置。然而，当入射波与系泊船舶发生耦合作用时，船舶将在初始阶段被推离其原始平衡位置。因此，在初始阶段，船舶的纵荡运动逐渐加剧。由于在此过程中入射波高几乎等于船舶吃水深度，波浪和船舶将发生剧烈耦合作用。船舶出现最大纵荡值的时刻并不是系泊缆被拉直的时刻，这是由于船舶特殊的外形构造，因此与常规的规则物体相比，其重心位置并不在结构中心位置，因此导致了系泊缆被拉直时船舶的纵荡幅值并未达到峰值。

从图7可以看出，随着入射波周期的增加，船舶的垂荡运动增加，但纵摇运动减小。以周期1.58 s为例，垂荡运动由初始的0.107 m增加至0.175 m，而纵摇运动由19.1 °下降至15.8 °。

图8所示为船舶所受波浪砰击力的变化曲线。可以看出，周期越短，波浪垂向砰击力对船舶的作用越大，砰击力峰值相差可达200 N，波浪周期越小，这意味着波长越短。因此，当波浪传播到船舶附近时，更难完全通过，在这种情况下，波浪必然会对船舶产生强烈的砰击作用。

图9所示为2个典型周期下船舶运动的流场云图。从云图上可以清楚地看到，当短周期波的波峰位于船中时，长周期波的波峰刚刚到达船艏位置。然而，在船舶经历第1个波峰后，由于波浪周期较短，将出现第2个波峰。对于长周期波浪，第1个波峰刚刚越过船尾，第2个波峰尚未出现。这2种情况导致波浪和船舶之间的运动相位不同，导致船舶在短周期波浪条件下的运动响应更强烈，而长周期波浪则相反。砰击力曲线的变化进一步证实了这样一个结论，即周期越小，船舶的纵摇运动就越强烈。

3.2 不同浪向

在实际工程中，系泊船舶与入射波之间的角度并不总是0 °。因此，本文将对不同波浪方向下波浪和系泊船舶的冲击响应特性进行研究。因此，所选的波浪条件为波高0.3 m、周期1.58 s、波向角分别为0 °、30 °、45 °、60 °和90 °。

需要注意的是，90 °浪向意味着波浪为横浪。此外，如图10所示，为了简化计算，本文将波浪浪向角转换为船舶与入射波之间的初始偏转角。通过这种方法，模拟了波浪方向的变化，系泊缆的相应位置也需要相应地改变。

图11所示为不同浪向下船舶六自由度运动曲线的比较。可以看出，随着入射波的波向角逐渐增大，船舶纵荡运动的幅值逐渐增大，最大值在90°浪向（横浪）条件下达到0.298 m，0 °与90 °的纵荡运动幅值相差2 m。此时，由于系泊缆位置的变化，船舶在设定的计算期内不会被系泊缆拉回其初始位置。此外，可以看出，浪向角度的变化对船舶的垂荡运动影响相对较小，并在整个计算周期内保持有规律的变化状态。然而，随着波浪方向角的增加，垂荡振幅也逐渐增加，这表明随着船舶与波浪接触面积的增加，船体将逐渐承受更多的波浪载荷。从图11(b)所示的船舶横荡运动曲线可以看出，在0 °条件下，船舶没有横荡运动，随着浪向角度的增加，船舶的横荡运动逐渐增大。这是因为随着浪向角度的变化，船舶逐渐开始沿Y方向受到波浪砰击力。

从图12所示的船舶所受波浪砰击力的比较曲线可以看出，当浪向角为0时，船舶几乎不受Y方向的波浪砰击力。然而，当浪向角不为0时，Y方向的波浪力开始出现，导致船舶横荡运动的变化。

从图13所示的横浪作用前期船舶运动姿态云图中还可以看出，船舶仍能保持一定的稳性，保证横摇角在较小范围内变化。然而，从图11(e)所示的船舶纵摇曲线可以看出，在横浪作用下，船舶的纵摇并不像其他浪向那样呈规律性变化，而是在计算期间逐渐增大。

图14所示为船舶在横浪后期的运动姿态。可以观察到，随着横浪的持续作用，船舶首尾出现了明显的甲板上浪现象，船舶的纵摇角度不断增大，随时都有倾覆的风险。因此，在进行单点系泊时船舶应避免横浪作用。

4 结　语

针对SPH方法在求解自由表面流与三维系泊浮式结构耦合作用时计算效率较低的问题，建立了求解三维波浪与系泊结构物耦合作用问题的SPH-GPU数值计算模型，弥补SPH方法在计算效率上的缺陷，同时将该耦合模型应用于船舶在波浪环境中的单点系泊问题。结果表明，短周期波浪导致更大的垂向砰击力与纵摇运动，而垂荡运动随周期增加而增强，纵摇则减弱；初始阶段系泊缆预张力引发船舶横向位移，但波浪耦合作用加剧纵荡响应，最大纵荡值因船舶重心偏移未出现在系泊缆拉直时，同时，短周期波浪因相位差导致更频繁波峰冲击，加剧运动响应；随着浪向角增大，纵荡运动幅值显著增加，横荡运动仅在非0°浪向下出现且随角度增大而加剧，垂荡振幅因接触面积增加逐渐上升。横浪条件下，船舶初期稳性较好，但后期纵摇角度持续增大，甲板上浪现象明显，存在倾覆风险。