船舶在恶劣海况下航行易发生砰击现象，即船体与波浪发生剧烈冲击，该现象主要发生于船艏。尤其是在迎浪高速航行时^[1]，受高等级风浪的影响，船舶发生剧烈的纵摇运动，导致船艏露出水面，然后以较大的相对速度落入水中，其产生的冲击压力常可使船艏结构严重变形或破坏。随着海况的恶劣，船舯和船艉也难以幸免。砰击使船体受到很大的冲击压力，其大小与船体局部的相对速度及剖面形状有关。入水砰击通常可分为2类：波浪对船底结构产生的底部砰击以及对舷侧外飘产生的外张砰击。吃水较浅且底部平坦的船舶(如登陆舰等)以底部砰击为主；对于艏部有较大外张的船舶(如航空母舰，集装箱船等)以外张砰击为主。

关于入水砰击的研究最早可以追溯到20世纪30年代，如Von Karman^[2]和Wagner^[3]的相关研究。然后，关于砰击的理论研究和试验有了一定的进展，随着计算机技术的发展，近些年数值模拟也成为了预报砰击载荷的有效手段之一。Zhao等^[4]对二维楔形体的入水砰击问题进行了研究。陈震等^[5]通过Dytran软件对平底结构的砰击压力进行了模拟。陈小平等^[6]研究了流固耦合效应对弹性结构入水响应的影响。Qin等^[7]基于改进的Logvinovich法对带横摇运动非对称楔形体的入水砰击力进行了计算。骆寒冰等^[8]基于显式有限元方法对二维楔形体入水问题进行了研究。张岳青等^[9]发现楔形体密闭空腔会造成加速度最大值增大近一倍。Hu等^[10]用CIP方法模拟了双楔形体的入水问题。Bihnam等^[11]基于OpenFOAM研究了流体可压缩性对二维楔形体入水过程的影响。更多关于入水砰击问题的研究可参见文献[12-14]等。

以往文献的主要研究对象为楔形体或简化后的球鼻艏，对形状复杂的实船型剖面的研究较少，本研究将对实船的3个典型剖面展开研究。研究基于显式有限元方法软件LS-DYNA建立了二维楔形体模型，并与已有试验数据进行对比分析，验证该数值模拟在处理入水砰击问题上的可靠性和精确性。然后，根据S-175集装箱船的横剖线图建立具有代表性的第19站、第10站和第1站的实尺度模型，分别作为二维球鼻艏、平行舯体和船艉剖面模型。最后，根据何广华等^[15]对S-175船型在波浪中强非线性运动的模拟结果，计算出上述各剖面处的下沉速度，以此为参考来设置本数值模拟的入水速度，分析各剖面在入水过程中的速度和压力变化等情况。

1 任意拉格朗日-欧拉算法

基于显式有限元方法，采用能够有效跟踪物质结构边界运动，同时可以避免网格出现严重畸变的任意拉格朗日-欧拉算法(arbitrary lagrangian-eulerian，ALE)算法。ALE算法的控制方程可由质量、动量以及能量守恒方程给定：

$ \frac{\partial \rho}{\partial t}=-\rho \frac{\partial v_{i}}{\partial x_{i}}-w_{i} \frac{\partial \rho}{\partial x_{i}} $

(1)

$ v \frac{\partial v_{i}}{\partial t}=\sigma_{i j, j}+\rho b_{i}-\rho w_{i} \frac{\partial v_{i}}{\partial x_{j}} $

(2)

$ \rho \frac{\partial E}{\partial t}=\sigma_{i j} v_{i, j}+\rho b_{i} v_{i}-\rho w_{j} \frac{\partial E}{\partial x_{j}} $

(3)

式中：t为时间；ρ为密度；相对速度w为物质速度v与网格速度u之差；w_{i, j}、v_{i, j}分别为速度分量；x_{i, j}为欧拉坐标；σ_ij为应力张量。

LS-DYNA数据处理器采用中心差分法按时间递增进行求解，采用时间显式法，提供二阶时间精度。每个节点的速度和位移进行更新：

$ u^{n+1 / 2}=u^{n-1 / 2}+\Delta t \cdot \boldsymbol{M}^{-1} \cdot\left(\boldsymbol{F}_{\mathrm{eut}}^{n}+\boldsymbol{F}_{\mathrm{int}}^{n}\right) $

(4)

$ x^{n+1}=x^{n-1}+\Delta t u^{n+1 / 2} $

(5)

式中：F_int为内力矢量；F_ext为外力矢量；M为质量对角矩阵。

模型采用罚函数约束的方法实现流固耦合，罚函数耦合系数追踪拉格朗日节点和欧拉流体节点间的相对位移d，检查节点是否出现贯穿现象。如果发生贯穿现象，界面力F会分布到欧拉流体的节点上，其大小与贯穿数量成正比：

$ F=k_{i} \cdot d $

(6)

式中k_i表示刚度系数。

经多次模拟和验证，最终水和空气均采用空白材料，选用线性多项式方程作为其状态方程，其压力可表示为：

$ P=C_{0}+C_{1} \mu+C_{2} \mu^{2}+C_{3} \mu^{3}+\left(C_{4}+C_{5} \mu+C_{6} \mu^{2}\right) E $

(7)

$ \mu=\frac{\rho}{\rho_{0}}-1 $

(8)

式中：C_i为多项式方程系数；ρ为当前密度；ρ₀为初始密度；E为内能；C₀=C₂=C₃=C₆=0。参数设置如表 1所示。

2 数值模型及其验证 2.1 可靠性验证

为了验证数值模型的可靠性，首先基于Zhao等^[4]的模型试验，对二维楔形体的入水问题展开了研究，并将计算结果与试验数据进行了验证。模型入水速度为6.15 m/s，试验楔形体模型参数如下：模型宽0.5 m，总长1 m，斜升角30°，总重241 kg。压强测量点示意图如图 1所示。

图 2给出了不同网格尺寸的数值计算结果与试验数据的对比验证，图中实线是本数值模拟结果，三角点为Zhao等^[4]试验数据；P₁~P₅为楔形体表面各压强测量点处的压强。从图 2中可以看出，最小网格尺寸0.001 25 m时的模拟结果整体砰击力波动变化更小，与试验数据吻合较好。

2.2 船体模型

模型选取S-175集装箱船的横剖线图中具有代表性的第19站、第10站和第1站型线建立实尺度模型，作为二维球鼻艏、平行舯体和船艉的剖面模型。船体横剖线图和3种剖面模型如图 3和图 4所示，单位为m。

S-175船型船长175 m，宽25.4 m，型深15.4 m，吃水9.5 m，排水体积24 742 m³，重心纵向位置据艉垂线84.97 m，垂向位置据龙骨线8.5 m。各剖面的入水速度如下^[15]：球鼻艏剖面入水速度7.3 m/s，平行舯体剖面2.7 m/s，船艉剖面3.8 m/s。

考虑到船体的对称性，文中均采用半剖面模型进行模拟。考虑研究对象为刚体，且研究重点为底部和外张砰击，因此在保证面质量不变的情况下对球鼻艏的内部进行了掏空处理(图 4(b))，从而保证了模型的网格质量。

2.3 网格划分

由于采用半剖面模型，因此计算域也做出了相应的设置。计算域左边界面(x=0)为对称边界，四周边界均采用无反射边界。经过系统的数值试验，计算域最终确定为水深30 m，水面上40 m，半宽100 m；流体区域采用非均匀网格，近入水区域网格均匀布置，局部最小网格尺寸为0.05 m，根据库朗数设置时间步长为0.01 s。在半剖面模型表面的x为1.5、3、6 m处设置压强监测点(P₁, P₂, P₃)。球鼻艏模型的网格布置示意如图 5所示。

3 数值模拟与结果分析 3.1 球鼻艏剖面

船艏是遭遇砰击的主要部位，因此对形状复杂的球鼻艏进行入水砰击研究至关重要。根据S-175集装箱船横剖线图中的第19站型线和已有研究中的运动响应，建立了实尺度的数值模型。

图 6和图 7分别为球鼻艏的外张砰击过程以及在入水过程中的速度变化曲线。从图 7中可以看到速度变化存在3个转折点：第1个速度转折点在0.75 s左右发生，此时球鼻艏的底部开始接触水面；第2个转折点在1.83 s附近，球鼻艏底部溅起的波浪到达舷侧从而产生了砰击现象；第3个转折点在约2.01 s时，由外张砰击产生的气泡消失，液面变化趋于稳定。

考虑球鼻艏以外张砰击为主，因此各记录点(x为1.5、3、6 m)均位于球鼻之上的舷侧外表面。图 8展示了球鼻艏剖面的入水砰击力和2个记录点处的压力。

与图 7中3个速度转折点对应，球鼻艏入水过程中，砰击力也存在3个变化阶段，如图 8(a)所示：第1阶段球鼻入水，船艏剖面受到一个较小的砰击力；至1.83 s附近时，舷侧开始受到波浪的砰击，受力逐渐增大；在船波之间的气泡消失殆尽，舷侧波浪逐渐稳定时，会形成一次峰值很大的冲击。从图 8(b)可看到，P₁处的压强存在2个同等量级的峰值，第1次出现在波浪开始拍击舷侧，第2次出现在船波之间气体消失时。记录点P₂处的峰值时间位于P₁ 2次峰值之间，波浪拍击舷侧后逐渐漫延到P₂处，其压强峰值大于P₁，随后逐渐减小并在P₁第2次峰值的相同时间产生一个略小的波动。

3.2 平行舯体剖面

平行舯体剖面形状接近矩形，主要以底部砰击为主。根据S-175集装箱船的第10站型线和已有研究中的运动响应，建立了实尺度的数值模型，为保证网格质量并减少计算量，在保证平行舯体面质量不变的情况下，对模型进行了薄壳化处理。

与球鼻艏不同，平行舯体的剖面形状和入水过程相对简单，其速度变化整体更为平滑。从图 9可以看到，平行舯体在很短的时间内就达到了底部与水的全面接触，水花沿船底向外侧溅射。如图 10所示，平行舯体在0.3 s左右开始接触水面，在整个底面的入水过程中，下沉速度会稍有缓和，待底面充分接触自由面后，速度开始平稳变化，呈现出先增大后减小的趋势。

从图 11可以看到，平行舯体的砰击载荷主要在接触水面时产生，因为底升高较小，在剖面入水时砰击压力较大且不稳定，压力峰值由船底向舷侧递减，该趋势与陈震等^[5]对平底结构砰击压力研究的结论一致。当底面全部入水后，剖面受力随着下沉深度逐渐增加。平行舯体入水速度峰值虽远小于球鼻艏，但其砰击压力情况因剖面形状而更加恶劣，且其压力峰值分布沿中线面向舷侧逐渐递减，发生时刻也依次滞后，入水过程稳定后，各记录点处的压力基本一致。

3.3 船艉剖面

船艉剖面与球鼻艏剖面类似，主要以外张砰击为主。根据S-175集装箱船横剖线图中的第1站型线和已有研究中的运动响应，建立了实尺度的数值模型，为保证网格质量并减少计算量，在保证船艉面质量不变的情况下，对模型进行了薄壳化处理。

在船波相对运动中，船艉底部较船艏不易发生砰击现象，从图 12和图 13可以看出，船艉在底部入水时速度并未发生明显改变，在整个底面的入水过程中，速度趋势几乎不变；至外张部分开始接触水面后，速度开始有了明显的转变。

图 14展示了船艉剖面的入水砰击力及记录点1、2处的压力。从图 14可以看出，在船艉底部开始入水时，砰击力逐渐缓慢增长，至外张部分接触水面时，艉部砰击力迅速增加达到峰值，随后回落至稳定状态。

3.4 3个典型剖面

通过对3个最具代表性剖面的研究，有利于对整个船体砰击响应的分析，剖面的形状对入水速度的变化和压力曲线有着重要影响。

图 15展示了3个典型剖面的入水速度、砰击力和压力曲线的对比。在同一海况下，船艏可达到的下沉速度最大；通过图 15(a)可以看出，船艉底部入水对下沉速度的影响可忽略不计，但球鼻的入水可以明显减缓下沉速度的增大。图 15(b)、(c)揭示了：速度大小不能直接反映压力情况，速度最小的平行舯体剖面受到的砰击力最大，且无论是整体砰击力还是记录点处的压力峰值都远大于球鼻艏和船艉剖面。另外，船艉剖面的砰击载荷较球鼻艏和平行舯体剖面的小。

4 结论

1) 球鼻艏剖面入水过程中的速度和压力存在3个变化阶段，在外张砰击中，气泡消失时存在二次砰击现象；

2) 平行舯体虽然在同一海况下的入水速度最小，但一旦发生砰击，其剖面形状使得其砰击力较大，底面的压力分布由船底向舷侧递减；

3) 船艉剖面在3个剖面入水中砰击力最小，且压力变化最简单；速度变化和砰击载荷主要发生在外张部分的入水过程中，底部对速度和砰击力影响很小；

4) 剖面形状对下沉速度和砰击载荷具有较大影响，且速度不能直接反映压力情况，船艉部分的砰击载荷较小。

该研究重点在于探索S-175集装箱船3个典型剖面在入水过程中的速度、砰击力和压力变化，尤其是形状复杂的球鼻艏剖面。本研究结果对S-175集装箱船的底部砰击载荷和外张砰击载荷的预报以及船型优化具有一定的参考价值。