0 引　言

结构物入水现象在工程领域广泛存在，如船体冲击、海上救生艇投放、水上飞机降落以及空投鱼雷入水等。入水冲击是典型的强非线性问题，涉及流体飞溅、破碎及空气效应等物理过程，作用机制极为复杂，已成为航空、航海等多领域亟待解决的核心问题。

当结构物入水发生砰击时，压力的时空分布受入水速度、几何形状、结构弹性、流体压缩性等多种因素综合影响，呈现出高度复杂性和不确定性。目前尚无完全适用的解析方法来分析其压力特性，相关研究主要集中于实验和数值分析。部分学者虽针对特定入水问题提出了局部解或特定数值方法，如Zhao等^[1]针对折角处的流动分离问题建立了局部解；Sun等^[2]采用非线性边界元方法求解了机翼剖面的入水波浪问题，但这些方法在处理复杂结构和流动分离现象时存在局限性，难以满足实际需求。

近年来，为准确模拟入水过程中涉及的复杂气穴效应，计算流体动力学（CFD）方法得到了广泛应用。沈雁鸣等^[3]采用光滑粒子水动力（SPH）方法分析了二维楔形体入水冲击的压力分布特征。Wen 等^[4]进一步将楔形体入水问题扩展到三维情况，分析了不同低升角下的压力响应特征。Yang 等^[5] 将自适应粒子细化（APR）和粒子移动技术（PST）相结合，进一步扩展到了结构的出水问题中。针对复杂的曲面三维结构，Shi 等^[6]基于笛卡尔网格多相流模型，探讨了三维外飘结构入水的载荷特性，并发现了射流砰击效应。Xie 等^[7]研究三体船入水的气泡流动特性等。然而，现有研究多集中于简化结构，对带附体的复杂结构，如带球鼻艏的外飘结构、头部带导航装备的鱼雷等的入水冲击研究较少。Liu等^[8]在尾部入水冲击实验中发现了尾轴对砰击载荷的重要影响。实际上，附体的存在会使流体在爬升过程中出现流动分离、破碎及明显的气穴效应。这些现象模拟的合理性直接影响砰击载荷预报的准确性，因此对附体结构入水的研究显得十分必要。

为进一步阐明附体对结构物入水时流体动力特性和冲击载荷的影响，本研究基于纳维 - 斯托克斯（N - S）方程建立了附体结构物入水砰击问题的气 - 液两相流模型，充分考虑气液两相的相互作用及流体的复杂流动特性。以带球鼻艏附体的复杂外飘结构为研究对象，采用流体体积函数（VOF）方法进行数值模拟。重点研究球鼻艏引起的自由面演变和压力特征，并讨论不同入水速度对砰击载荷的影响，为实际工程结构的入水设计和安全评估提供参考。

1 计算模型及研究对象 1.1 控制方程

基于N-S方程建立入水冲击问题的两相流模型。假设流体不可压，并忽略流体的粘性和重力效应，则连续性方程和动量方程简化为

$ \nabla \cdot {\boldsymbol{u}} = 0，$

(1)

$ \frac{{\partial {\boldsymbol{u}}}}{{\partial t}} + ({\boldsymbol{u}} \cdot \nabla ){\boldsymbol{u}} = - \frac{1}{\rho }\nabla p。$

(2)

式中：p为压力；u为速度矢量；ρ为流体的密度。应用VOF方法求解上述方程，假设水为基本相，引入相方程：

$ \frac{{\partial \alpha }}{{\partial t}} + \nabla \cdot(\alpha {\boldsymbol{u}}) = 0。$

(3)

式中：α为基本相的体积分数。则流体的密度可进一步表示为

$ \rho = \alpha {\rho _1} + (1 - \alpha ){\rho _2}。$

(4)

式中：${\rho _1}$、${\rho _2}$分别为水和空气的密度，取为1.0×10³ kg/m³、1.25 kg/m³。式(1)～式(3)构成了应用VOF法求解两相流模型的连续性方程、动量方程和相方程。

流动方程采用Fluent中标准的k - ε湍流模型进行求解。控制方程的离散化通过有限体积法实现。针对时间项，采用显式欧拉格式进行处理。对流项则选用中心差分格式进行近似。速度场与压力场的离散均采用二阶精度的迎风格式。入水过程中耦合六自由度运动求解器进行运动速度的指定。采用重叠网格进行运动的处理。对于自由曲面的捕捉，采用分段线性界面重构（PLIC）格式。对于利用VOF方法解决曲面截面入水问题的详细数值描述，可参考Xie等^[9]的研究。

1.2 研究对象

本研究以船体砰击中常见的带球鼻艏附体的外飘结构作为研究对象。选取某船典型外飘结构按1∶35的缩尺比进行建模，几何尺寸为1.52 m×0.93 m，入水问题示意图如图1所示。模型以垂向速度V匀速入水。气-液两相的交界面上侧为空气，下侧为水。结构本身划分为3个区域：底部的球鼻艏区域、球鼻艏与外飘相连的外飘区域1和远离球鼻艏的外飘区域2。

整个流体域尺寸为10 m×6 m，水深为7 m。基于ICEM CFD软件对流域网格进行划分。网格类型以四边形网格为主。模型位于计算域的中心。整体网格尺寸控制在0.1 m。为了准确捕捉自由表面流动的动态特征，网格在模型上下3倍于模型高度的区域内进行细化。该区域的细化网格尺寸统一设置为0.02 m。同样，为了保证结构运动方向上的流动细节，网格也在模型的主要运动区域进行细化。此网格尺寸为0.02 m。在内部的重叠区域使用不同的网格配置对模型周围的中间区域进行细化。图2为流域网格细节。

1.3 求解过程

步骤1　对相方程(3)进行求解，更新α和流体密度$\rho $；

步骤2　基于有限体积法对动量守恒方程进行空间离散化，结合二阶迎风格式完成对流项与扩散项的数值离散，通过线性方程组迭代求解获得速度场的预测值；

步骤3　采用PISO算法^[10]对压力-速度耦合修正策略实现压力场和速度场更新。

步骤4　利用动网格更新技术实现物体的入水过程^[11]，回到步骤1重新进行下次迭代直到计算终止。

2 数值验证

为了验证本文数值方法的可靠性，选取Aarsnes^[12]中的外飘剖面进行建模，模型尺寸为0.32 m×0.24 m，总质量为261 kg，落体高度为0.318 m，数值模拟中入水速度来自于实验值。

首先进行网格的收敛性研究。表1中列出了3种不同的网格尺寸：Grid1、Grid2和Grid3。图3给出了不同网格参数下压力和总垂向砰击力的对比结果。可知，Grid1的预测结果有所偏小，而Grid2和Grid3的结果基本相同。这表明当网格尺寸小于Grid2时，数值结果收敛。从计算时间来看，Grid3的计算时间是Grid2的2.1倍，因此可采用Grid2的网格尺寸来平衡计算时间和计算精度。

下面分析数值结果的误差，定义相对误差$ \varPhi $来评估数值结果的差异性，即

$ \varPhi = |{Q_{{\text{num}}}} - {Q_{\exp }}|/{Q_{\exp }} 。$

(5)

其中：一般项Q_num和Q_exp分别为本文数值解（Grid2结果）和实验中相应物理量（力和压力）的最大值。对于总的砰击力，数值和实验的相对误差为4.8%。需要注意的是，实验结果中存在一定振荡，其控制频率约为100 Hz，来自于实验过程中落体架的振动^[13]。对于局部的砰击压力，相对误差为3.1%。事实上，总的砰击力为计算方法在整个流场上的模拟精度，而局部压力为局部流动的模拟精度。从对比结果来看，本文方法在整体和局部流动的模拟中均小于5%，满足工程精度。

3 结果分析

应用数值方法对带有球鼻艏附体的外飘结构入水过程进行研究，模型存在对称性，选取1/2模型进行数值模拟。流域尺寸取为10 m×15 m，水深设为10 m。流域的上部为压力出口边界，底部和两侧为无滑移边界条件，物体表面为无滑移的动边界。物体下落过程中假定为匀速入水，入水速度V分别为2、4、6、8 m/s。

3.1 自由面特征分析

图4给出了不同时刻下自由面的演变过程。不同于常规的楔形体入水，带附体外飘结构入水时，球鼻艏的存在会产生一些复杂的流动现象。在初始阶段，球鼻艏首先没入水中，由于球鼻艏和外飘区域交界处的凹面，爬升的流体从物体处分离，如图4(a)所示。这一阶段出现的流动形式类似于圆柱入水时的流动分离^[14]，分离点并非固定，临界条件为分离点处压力为0，受流速和物体几何形状影响。随着物体入水深度的增加，分离的流体重新撞击到物体表面，产生射流砰击现象，如图4(b)所示。在这个阶段，空气将被捕获，形成气穴，并依附于物体的表面，如图4 (c)所示。流体上升过程中，气穴受流场压力影响产生了变形，如图4(d)所示。最后流体在折角处再一次发生分离，分离后的流体逐渐上升，如图4 (e)和图4(f)所示。不同于图4(a)中发生的分离现象，这种流动分离发生在物体固定的折角点处，仅取决于物体的几何形状，此时折角处的压力并非为0。

3.2 压力特性分析

图5给出了图4对应时刻的无因次化压力系数Cp沿纵轴z的分布，其中$ Cp{\text{ = }}2p/{\rho _w}{V^2} $。可知在首次流动分离时（见图4(b)），此时流场速度剧烈变化，导致球鼻艏与外飘交界的内凹处出现了负压。随后该处压力逐渐增加，由负压转变为正压。当流体在折角处发生分离后（见图4(d)），整个物面上的压力分布逐渐趋于稳定，这表明流体从折角分离后，整个流场区域稳定。图6给出了对应时刻的压力云图。在t=0.04 s，发生流动分离的区域存在负压。球鼻艏的凸起导致内部的压力不足以维持其贴体，从而发生分离现象^[15]。在射流砰击发生时（t=0.06 s），射流位置出现局部的高压。随后气泡封闭后，受流体的挤压高压区域扩散。图7给出了入水冲击速度8 m/s时的压力分布特征。与6 m/s的入水情况相比，两者的压力分布演变特性基本相同。这说明压力的规律在不同速度下具有一定的相似性。

3.3 二次砰击特性分析

为了进一步分析分离的流体撞击物面时冲击载荷的时历特征，在外飘区域定义4个典型压力监测点，如表2所示。其中，P₁、P₂点位于外飘区域1；P₃、P₄点位于外飘区域2。图8给出了初始入水速度6 m/s时不同测点处无因次化压力系数Cp的时历特征。可以看到P₁、P₂处在t=0.064 s处出现了历时极短（$\Delta t \approx 0.005\;{\text{s}}$）而极大（$Cp \gt 3.1$）的砰击压力，之后在t=0.08 s时又出现了一个压力峰值。这表明P₁和P₂处发生了2次砰击作用，这种现象称之为二次砰击现象^[16]。在图4(b)中给出了P₂达到第一次压力峰值时自由面的形状。可知，此时分离的流体刚好撞击到P₁和P₂所在区域，由于此时物体和流体的相对速度较大且有效砰击角也小于此处的斜率，导致砰击几乎在一瞬间达到峰值。P₁和P₂处的第2个峰值附近的压力时历特征与P₃和P₄处类似，是入水过程中结构排挤水产生的压力，类似于常规的楔形体砰击压力特征^[17]。

3.4 入水速度的影响

在3.3节的讨论中，入水速度对冲击载荷尤其是首次砰击压力峰值具有重要影响。为进一步明确这种影响，图9给出了不同速度入水下无因次化压力系数Cp的对比结果。可知随着入水速度的增加，P₁和P₂处的Cp随入水速度增加而增加，而P₃和P₄处随速度的变化并不明显。说明球鼻艏附体在一定程度上会增大与球鼻艏连接区域（外飘区域1）附近的压力（P₁和P₂），而且这种影响会随着入水速度的增加而增大。图10给出了不同入水速度下的峰值时刻压力云图。可知在球鼻艏上方的外飘区域1出现了远大于外飘区域2的压力。这与该处射流冲击封闭气泡逃逸的路径后引起的挤压有关。对于远离球鼻艏的区域（外飘区域2），入水速度对砰击压力（P₃和P₄）的影响并不明显。因此，通过应用无因次压力系数进行结构设计时，存在附体的结构应当格外关注由此引发的射流冲击压力（图8中第1个峰值）。在这种情况下，虽然仍可使用无量纲的Cp，但需要注意Cp并非独立于入水速度。

4 结　语

基于气液两相流模型对带附体的外飘结构入水过程进行了数值模拟。通过分析入水时流体的演变，压力的特性及不同速度对其影响，得出以下结论：

1）由于球鼻艏的影响，入水中伴随着流动分离和气穴效应，模拟中需要合理对待；外飘剖面入水中出现了明显的二次砰击现象，一定程度上增加了与附体连接处（外飘区域1）的砰击压力。

2）靠近附体区域处压力会出现双峰值，第1个峰值来源于射流冲击，明显高于第2个，且作用时间极短，这种瞬时的高脉冲力在结构设计中值得注意。

3）入水速度对带附体结构的砰击压力具有重要影响。入水速度越大，靠近附体区域的无因次化压力系数越大。工程应用中需要尤为注意入水速度对该系数的影响。