0 引　言

在海洋工程领域中，船舶航行时的局部结构入水以及波浪对结构物相互作用时产生的上浪砰击是属于较为常见的流体动力学现象。在这些过程中，均牵扯到了自由表面流与结构的强非线性耦合作用，尤其是产生的水体高速飞溅，将会对船体结构以及内部人员产生威胁。

目前国内外已有许多学者对结构物出入水问题进行研究。提出了许多数值方法，包括流体体积法（VOF）^[1]、边界元法（BEM）^[2]、光滑粒子流体动力学法（SPH）^[3]等，Feng等^[4]基于BEM方法提出了一种弹性楔型体在浅水域入水问题的耦合方案。Cui等^[5]基于SPH方法建立了用于碰撞的HBP本构数值模型，研究了楔型体对非牛顿流体表面的撞击。Ping等^[6]对带有底部附体和不带有底部附体的2种船体模型进行了入水实验。Wang等^[7]设计了2种不同主船体形状的模型，一种能够延长主船体的三角帆船截面模型，另一种是一个主船体与原始三体船横截面相似的模型，都揭示了砰击载荷和流体力学的行为差异。上述方法可以简单描述二维结构物的入水砰击问题，但三维效应对于进水问题的研究也很重要，Wang等^[8]开展了全尺寸3D纤维增强复合材料在船首部位的系列砰击实验。Xie等^[9]通过附带俯仰角的方式进行入水试验，阐明了3D船尾模型的斜向砰击特性，揭示了三维效应和不同下落高度对峰值压力的影响。Sui等^[10]对具有俯仰角的高速弹丸进行了入水实验，以探究弹丸头部形状对空腔演变和冲击载荷的影响。基于此，对进水过程中的砰击载荷的研究几乎都基于简单的结构物，对复杂结构（如船体模型）入水问题的研究很少。Li等^[11]在具有倾斜角的三维船头上进行了跌落实验，并提供了对3D船头实际砰击的初步效果。Chen等^[12]通过将虚拟粒子边界和法向通量法相结合，提出了一种改进的固体边界处理技术，应用边界处理方法模拟救生艇的进水情况。Lyu等^[13]基于SPH方法，建立多分辨率3D-SPH模型，系统研究和讨论了不同入水角度下救生艇的水动力学行为。

在采用这些数值手段对这些问题进行研究时，由于涉及到了自由液面的翻卷大变形，若采用网格法进行模拟，则必须结合多种复杂的网格处理技术和自由液面捕捉技术，而SPH方法在对于自由表面流与三维结构的耦合求解上计算效率较低，难以适用。为此，针对上述问题，本文将基于优化的GPU加速技术构建三维SPH数值求解模型，对船舶入水砰击问题进行数值模拟研究。

1 SPH理论 1.1 控制方程

SPH方法中用于计算相应密度变化的动量方程和连续性方程为：

$ \frac{{{\rm{d}}v}}{{{\rm{d}}t}} = - \frac{1}{\rho }\nabla P + g + \Gamma，$

(1)

$ \frac{{\rm{d}} \rho_a}{{\rm{d}}t}=-\sum_b m_b v_{ab}\nabla _aW_{ab}。$

(2)

式中：${p_a}$和${p_a}$分别为粒子$a$和$b$的压强；$W$为光滑核函数。$\prod $为人工黏度，作用是取代流体的真实黏度，并提供剪切力，保证数值积分的稳定性。${\rho _a}$和${\rho _b}$分别为粒子$a$和$b$的密度，$g$为重力加速度。$W$和${\prod _{ab}}$的表达式可表示为：

$ W(r,h) = {\alpha _D}{(1 - \frac{q}{2})^4}(2q + 1)，0 \leqslant q \leqslant 2，$

(3)

$ \prod_{ab}= \left\{\begin{aligned} &\frac{-\alpha\mu_{ab}\overline{c_{ab}}}{\overline{\rho_{ab}}}，{\boldsymbol{v}}_{ab} {\boldsymbol{r}}_{ab}< 0 ，\\ &0，{\boldsymbol{v}}_{ab}{\boldsymbol{r}}_{ab}> 0。\end{aligned}\right. $

(4)

式中：无因次距离$q = r/h$为任意2个给定粒子之间的距离，记为$a$和$b$，其中$r$为实际距离，$h$为光滑长度。二维${\alpha _D}$值为$7/(4\text{π} {h^2})$，三维${\alpha _D}$值为$21/(16\text{π} {h^3})$。如式(4)所示，支持域的大小为$2\;h$。另外，令${r_{ab}} = {r_a} - {r_b}$，${v_{ab}} = {v_a} - {v_b}$。其中：${v_a}$为粒子$a$的速度，${r_a}$为粒子$a$的位移。平均声速为${c_{ab}} = 0.5 ({c_a} + {c_b} ) $，${\mu _{ab}} = h{v_{ab}}.{r_{ab}}/(r_{ab}^2 + {\eta ^2})$，${\eta ^2} = 0.01{h^2}$。粘性系数$\alpha = 0.01$。

由于SPH方法考虑流体为弱可压缩（WC），因此状态方程允许建立粒子密度与流体压力之间的关系，使用Cole状态方程可以得到：

$ P = {B_f}\left[ {{{\left( {\frac{\rho }{{{\rho _0}}}} \right)}^\gamma } - 1} \right]。$

(5)

流体的参考密度为${\rho _0} = 1\;000\;{\rm{kg}}/{{\rm{m}}^3}$，$\gamma $为几何常数，取值范围为1～7。参数${B_f}$与流体的可压缩性有关，可表示为${B_f} = {c^2}{\rho _0}/\gamma $。其中：$c$为参考密度下的声速，本文将计算模型的CFL值设为0.2。

1.2 三维问题GPU并行加速数值技术 1.2.1 改进的链表搜索算法

本文采用CLL方法来保存邻近粒子链表。图1为该方法的简化示意图。在该方法中，使用了4个单独的数组：Id、Cell、IdSort和CellBegin，上标*表示已排序的数组。先将粒子的信息储存在相应的单元网格中；根据Id数组中粒子索引建立Cell数组，并排序生成Cell*数组；使用Nvidia的radixsort算法，按照6个背景网格的顺序对Cell数组重新排序，并使用Cell*数组根据粒子所属的背景网格对IdSort进行重新排序，形成排序后的IdSort*数组；生成IdSort*后，考虑到Id_new[i]＝Id[IdSort*[i]]，所有描述粒子属性（Id、位置、速度、密度…）的数组都需要排序，从而产生新数组（Id_new、Position _new、Velocity_new、Density_new…）。例如，Id_new[4]=Id[IdSort*[4]]=Id[6]=2，如图1所示，蓝色圆圈标记粒子2，红色圆圈标记索引位置6。最后，根据每个背景网格中第一个粒子的索引位置创建数组CellBegin。

1.2.2 GPU并行加速优化技术

使用GPU并行加速时容易出现共享内存大小不足、计算线程工作负载不平衡以及程序执行出现分歧等问题。因此，需重点考虑GPU内部核心的占用率。对GPU并行计算的具体优化策略包括：GPU使用效率的最大化。如图2(a)所示，通过将背景网格的划分尺寸进一步减小，通过提高GPU线程束（Warp）的使用率来实现对更多计算单元的调用；减少全局内存的访问。舍去能从其他变量计算得到的物理量，同时将部分数组进行合并；简化邻近粒子搜寻的方法。如图2(b)所示，根据粒子编号在各方向的连续性，将邻近背景网格分成9组。以三维问题为例，对于一个给定的粒子，其相关的邻近单元总数为27个。然而，这个数量可以进行优化。当与目标粒子相互作用的粒子范围已知时，由于粒子按照单元进行重排，单元按照x、y、z轴的顺序排列，所以X轴上3个连续单元(cellx,y,z;cellx+1,y,z;cellx+2,y,z)的粒子范围等于cellx,y,z第一个粒子到cellx+2,y,z最后一个粒子。因此，可以将这27个单元定义为9个连续单元，如图2(b)所示。在这一过程中，减少了对内存访问的需求。

2 数值验证 2.1 收敛性分析

针对SPH-GPU模型的数值收敛性问题，本文选择了Jalalisendi等^[14]进行的三维船舶入水试验作为数值案例，船体质量为0.60 1kg，初始速度为2.98 m/s，仅受重力作用，做自由落体运动，选取3种不同的粒子分辨率（$ L/dp = 100,300,350 $，其中$ L $为船长参数，$ dp $为粒子间距），代码采用Nvidia Geforce3090显卡进行加速计算。如图3所示，基于GPU加速的三维SPH数值模型能够准确捕捉船舶入水问题中的动力响应特性，数值收敛性较好同时随着分辨率的提高，数值模拟的结果与试验的误差逐渐减小，综合考虑计算成本和求解精度，最终确定粒子分辨率$ L/dp = 300 $。图4为船舶入水砰击在典型时刻的速度云图，由于SPH方法具有独特的优势，无需采用复杂的技术便可直接对其进行精确捕捉，这也进一步说明了SPH方法在处理自由液面翻卷大变形领域的优势。

与此同时，本文采用CPU并行计算对该工况进行了同步计算，CPU设备为：Intel(R)Core(TM) i9-12900K,采用GPU加速计算时长比只采用CPU进行并行加速的时长缩短了约54倍，涉及到的其他计算时间对比如表1所示。

2.2 数值验证

为了进一步说明所建立的SPH模型计算在大尺度入水砰击问题上的精确性，本文选取具有一定初速度的船模结构作为测试案例，船模长度为1 m，宽度为0.25 m，初始纵倾角为30°，质量为10 kg，由于没有相关的模型试验，所以本文将采用国际主流CFD商业软件STAR-CCM+对模型进行同步计算，具体布置如图5所示。

在开始对比计算之前，本文需要首先验证STAR-CCM+数值计算模型的收敛性，为此选择了3种不同的网格尺寸进行收敛性分析，网格参数如表2所示。图6为不同网格分辨率下的船模速度和位移曲线对比图，可以看出STAR−CCM+数值计算模型具有较好的收敛性和可靠性。接下来，本文将使用三维SPH-GPU数值计算模型求解船模入水过程。

图7为SPH数值计算模型与STAR-CCM+数值计算软件计算结果的位移和速度曲线对比结果。可以看出，所建立的三维SPH数值计算模型得到的船模运动响应与STAR−CCM+结果吻合较好，说明SPH模型具有较好的可靠性。仅在0.84 s时，纵摇值最大相对误差达到12.46%。图8为基于STAR−CCM+和SPH数值模型在典型时刻下船舶模型姿态以及入水过程中产生的自由液面飞溅现象对比。可知，通过SPH数值计算模型和STAR−CCM+数值计算软件求解器得到的船模姿态几乎相同，而且基于拉格朗日特性的SPH数值计算模型能够准确捕捉自由表面流与结构耦合作用产生的飞溅射流和自由面破碎变形现象。

因此，本文所建立的三维SPH-GPU数值计算模型具有较高的计算精度和更好的适用性。

3 计算结果分析 3.1 数值设置

本文选取的船型参数如表3所示，计算域尺寸如图9所示，本文选取的缩尺比为1∶10，计算的粒子分辨率为$ L/dp = 300 $。

因为选择的船模具有复杂的曲面，而初始化后的船模粒子分布对计算结果影响较大，所以本文对于船模的离散采用了Zhu等^[15]开发的前处理数值技术，如图10所示。

由于砰击属于瞬时过程，本节设置的计算时间为1.5 s，由于粒子数量较多，所有程序都在A100高性能计算显卡上进行，为了提高计算精度，采用双精度计算。本节主要考虑船舶质量以及出水时的速度对船舶在跃出水面后，对水体的砰击作用，因此设置了3种不同的质量，分别为200、300、400 kg，对应船舶空载、压载和满载状态，初始速度的水平分量${v_x}$分别设置为2、4、6 m/s，对应船舶合同航速、设计航速和结构航速，本文将船舶初始纵倾角设置为10°。

3.2 不同初始速度

图11为不同初始航速下不同质量船舶的纵荡运动响应曲线，可知，对于相同质量的船舶，当初始纵倾角给定时，船舶的纵荡运动随着速度的增加而增大，当船舶形状和初始纵倾角保持不变时，可以认定附连水质量的影响以及初始水体阻力的影响均一样，速度对于船舶出入水过程中纵荡运动的影响远大于质量的影响，可很清晰看出，对于相同速度的船舶，随着质量的增加，其纵荡运动的幅值会有所不同，但差异不大。然而，对于相同质量的船舶，当速度变化时，其纵荡幅值之间的差异较为明显。

图12为不同初始速度下流体对船舶的横向砰击力，明显可知，速度越大，横向砰击力越大，特别是在初始阶段（t<0.15 s），由于初始速度越高带来的初始动能越大，产生的瞬间砰击力也越大。

图13为不同初始速度下船舶受到的垂向砰击力，可知，在后期，速度越高，流体对船体的垂向砰击力越大，在这种情况下，船舶出水再重新入水时，会被水抬高，其垂荡幅值自然会小一些。在前期，低速船受到的砰击力要大于高速船，因为低速船初始速度较小，其初始动能可以在短时间内转化为系统势能，因此可以更早地与水面发生全面接触，而此时高速船仍保持斜向上的运动姿态，只有尾部与水面相互作用，承受着水体砰击力，因此，低速时，流体对船舶的砰击力较大。

图11描述了船舶纵倾角的时域变化曲线，可知，在砰击的初始阶段（t=0.12 s），速度越大，纵倾角越大，但随着时间的推移，速度越大，纵倾角越小，而速度越小时，其船舶纵倾角越大。以质量200 kg为例，图14为船舶在不同初始速度下入水过程中几个典型时刻流场变化云图。可知，初始速度越大，纵倾角越大，当t=0.25 s时，可清楚看到不同速度下船舶纵倾角度的变化，由于初速度的垂直分量较大，因此在初始阶段船体的纵倾角度较大。随着时间的推移，初始速度较低的船舶由于其动能较小，会在很短的时间内重新砸回水面，且由于附着在船体表面的自由流体较多，总质量大于其他2种情况，因此船首入水深度更大。随着船舶航速的增加，一方面船舶与水面的接触面积同时变小，更有利于调整到水平位置入水，另一方面由于初始速度较大，附着在船表面的水的质量较小，使得总质量较低。因此，在发生入水过程时，船首入水深度相对较小。

图15为不同质量船舶在入水过程中流场变化的三维云图，可知，船舶质量越大，入水后期船首入水后产生的甲板上浪现象越明显。

3.3 不同初始倾角

研究不同初始倾角下船舶出入水过程。设置船舶质量为400 kg，初始纵倾角分别为10°、20°和30°，覆盖了船舶遭遇的常见海况和极端海况条件。初始速度的水平分量设为6 m/s。

图16为2个典型时刻的流场变化云图，可知，在前期由于纵倾角较大，重心位置较高，而入水后期，由于速度较快，入水过程中对水体的挤压和飞溅现象更加明显，不同纵倾角度下入水深度的差异也可清晰捕捉。与传统救生艇入水过程不同的是，在整个入水过程中，由于其倾角较小，速度慢，因此船舶自始至终都没有完全脱离水体，因此，由此产生的沃辛顿射流对船体结构安全的威胁较小。

如图17所示，当船舶以给定的初始速度运动时，纵倾角越大，垂直速度分量越大，垂向位移也就越大。在前期，能够避免过早地接触流体，当倾角为30°时，与其他2种工况相比，随着时间的推移，船舶受到的主要砰击力集中分布开始逐渐从船底向船艏过渡。同时，船舶在另外2种工况下已经实现了整个底部的入水过程，砰击力的分布更加均匀，从而减小受到的砰击力。但随着时间的推移，当船体落回水面时，倾角越大，积累的重力势能越大，砰击现象就会越剧烈，船体受到的砰击力也就越大。

4 结　语

本文针对SPH方法在求解自由表面流与三维结构耦合作用时计算效率低的问题，建立了基于GPU加速框架的三维SPH数值计算模型，对船舶出入水问题进行分析研究。本文的主要研究结论有：

1） 本文提出了一种基于GPU加速技术的高性能SPH数值模型，通过与已发表数据对比，验证了SPH-GPU模型在三维砰击入水问题中的准确性，并与传统CPU并行进行对比，表明GPU加速技术的引入可以大大提高SPH方法计算效率。

2） 在船舶出入水问题中，当船舶的质量给定时，速度越大，船舶的纵荡幅值越大。而在保持船速不变，仅改变船舶质量的情况下，当船舶出水后再入水时，船舶初始质量越小，其入水的深度越小，而初始质量越大时，在船首入水以后产生的甲板上浪现象越明显。

3） 在保证船舶质量和速度不变的情况下，仅改变初始纵倾角时，可以发现其对船舶的纵荡运动影响较小；此外，初始纵倾角越大，船舶在前期所受的砰击力越小，但随着时间的推移，船体落回水面时，由于前期较大的纵倾角导致产生了较大垂向位移，积累了更多重力势能，导致产生的砰击力更大。