航行体水下作业过程涉及问题综合性高，技术难度大^[1]，试验装置往往设计复杂，为了提高试验过程的可重复性，经济性，并尽可能形成试验装置的有效保护措施，提高试验装置的可回收性，对航行体下落入水过程进行研究具有极大意义。试验航行体的回收重点是实现航行体可打捞性的提高以及尽可能减小对水下发射平台的破坏威胁。目前航行体下落入水回收主要采用的方式是在下落过程中弹出助浮装置^[2-3]，携带助浮装置的航行体在下落入水过程中不仅仅可以依靠助浮装置进行冲击缓冲，而且助浮装置提供的浮力可以使得航行体下沉深度减小，进一步保护水下发射平台，助浮装置的存在使得航行体最终上浮至水面，因此带助浮装置航行体的设计也减小了试验航行体的搜寻与打捞难度。

航行体入水过程涉及结构冲击问题、气液两相跨介质问题以及流固动力学耦合问题，其冲击过程复杂，非线性显著，研究难度大，研究综合性高，近些年航行体入水问题也一直是学者们研究的重点问题^[4-7]。祁晓斌等^[8]基于VOF多相流模型和动网格技术对带尾翼航行体在高速小角度入水时的空泡发展和弹体姿态变化进行了研究。宋武超等^[9]利用数值模拟方法对不同头型回转体低速倾斜入水过程进行了分析，得到了不同头型条件下回转体入水空泡形态发展规律、回转体运动特性及流体动力特性变化规律。Dong等^[10]对无人潜航器入水问题进行了研究，其结果表明初始倾角对航行体入水过程具有一定影响，无人潜航器初始倾角对水平速度的影响随桨距角增大而减小，对水平力的影响在入水后期则更为显著。Yuan等^[11]对半球形弹头竖直入水阻力特性进行了数值研究，研究表明航行体头部平面扩张半角越大，弹头头部冲击速度越小，最大阻力系数值越大。Huang等^[12]研究了救生艇以不同倾斜角度与高度释放的下落入水过程，结果表明，较大的落水角可以减小入水压力冲击，但入水角也决定了救生艇入水时的运动模式，入水角度过低或过高都会导致救生艇入水后出现不良状态。

除了静水环境下的航行体入水问题，研究者们还对航行体在波浪环境下的入水进行了研究。史崇镔^[13]通过建立数值预报模型对不同波浪环境下飞潜器入水过程进行了模拟，结果表明波浪会对飞潜器入水运动特性产生影响，波浪工况下飞潜器具有更大的水平速度和更小的竖直速度，升、阻力系数也发生了变化。杨晓光等^[14]通过实验对波浪环境下的高速入水物体进行了研究，分别对静水、不同波面位置、不同波高等条件下的模型入水过程中的加速度、姿态角等物理量进行了测量，结果表明波浪对入水过程的影响主要体现在入水时入水角度相对于水面的变化。李治涛等^[15]分析了高速旋转射弹在不同波浪相位下的入水过程，发现波面的存在使得弹体实际有效入水角发生了变化，进而影响了入水运动特征。Zhang等^[16]在波浪环境下对楔形断面体的入水过程进行了数值模拟，重点分析了不同断面曲率楔形体在不同波浪相位处的入水流场发展情况以及入水物体的力学响应特性。Zhao等^[17]采用了试验与数值相结合的方法探究了周期性发展波浪作用下的球体入水空泡演化过程以及流动结构特性，具体分析了不同波浪相位入水位置和不同弗劳德数对入水小球所受的水动力特性影响。Cheng等^[18]基于高阶边界元方法研究了夹带气半圆航行器斜入五阶Stokes波的入水问题，通过数值模拟的方法研究了波幅、波长、初始相位、初始气压以及航行器速度等影响参数对航行器入水空泡的影响，发现波浪的存在延缓了入水气穴的形成，随着波浪非线性的增加，气穴的形状不对称性增加，气穴形成后流体压力急剧增加。

由于助浮装置在航行体下落入水过程中可以起到较好的缓冲与提供浮力作用，对带助浮装置的航行体下落入水冲击过程进行研究可以进一步形成并优化助浮装置结构设计。赵蛟龙等^[19]在考虑可压缩性的基础上推导出了带空腔复杂结构的落水数学模型，并研究了不同落水角度下结构表面水动力载荷的变化规律。陈洋等^[20]基于LS-DYNA对带环形密闭气囊弹体入水冲击过程进行了数值分析，计算了不同初始气囊内压条件下弹体的入水深度、减速时间及连接绳的拉力峰值。陈开颜等^[21]通过试验与数值结合的方式辨识出了带囊回转体落水过程的影响因素，并基于大尺度带囊回转体落水试验研究了不同囊-体连接结构形式的落水特征。包健等^[22]基于耦合欧拉-拉格朗日方法对不同姿态角和不同初始囊压条件下的带气囊航行体低速落水后的运动过程、囊压变化以及连接带的受力情况进行了研究。

综合上述，带助浮装置航行体落水相关研究为航行体落水回收的方案设计与结构参数设计提供了一定的参考依据，但是相关研究均为考虑静水环境下航行体的下落过程，而现实航行体下落环境往往伴随着波浪的发展，因此本文基于CFD数值模拟技术构建波浪环境下带助浮装置航行体下落冲击的数值模型，基于该模型，重点针对航行体动力学响应以及载荷变化情况，对波浪环境下带助浮装置航行体的下落冲击与下沉-上浮过程进行研究。

1.   数值理论基础

1.1   控制方程

本文中的流体为不可压缩的黏性流体，且过程中不涉及热力学能量的相互交换，所以采用的数值模拟控制方程不引入能量方程，仅包括质量守恒方程（连续性方程）和动量守恒方程：

式中：$\rho$为流体混合相密度；${u}_{i}$为流体速度在$i$方向上的分量；$t$为时间，${x}_{i}$为流体在$i$方向上的位置坐标；$p$为压力；$g$为重力加速度；${\tau }_{ij}$为黏性剪应力，具体定义如下：

式中：$\mu$为混合动力黏度系数项。

1.2   VOF多相流模型

本文通过VOF多相流模型来捕捉不同相流体之间的交界面，该模型以各相不同的的体积分数来描述不同的流体，并假设流体共享速度场和压力场，${\alpha }_{\mathrm{l}}$、${\alpha }_{\mathrm{v}}$、${\alpha }_{\mathrm{g}}$满足下式关系：

计算过程中每个流体单元的动态混合密度项和动态混合动力粘度项分别由以下两式确定：

式中：${\rho }_{\mathrm{l}}$、${\rho }_{\mathrm{v}}$、${\rho }_{\mathrm{g}}$分别为液相物质、水蒸气、其他气体的密度，${\mu }_{\mathrm{l}}$、${\mu }_{\mathrm{v}}$、${\mu }_{\mathrm{g}}$分别为对应的动力黏度。

1.3   湍流模型

本文采用的k-ω SST湍流模型考虑到湍流剪切应力的输运，综合了近壁面k-ω模型的稳定性及边界层外部k-$\varepsilon$模型独立性的优点，不但能够对各种来流进行准确的预测，还能在各种压力梯度下精确模拟分离现象。其输运方程为^[23]：

式中：${P}_{k}$为方程生成项，$k$为湍动能，$\omega$为湍动能比耗散率，${\mu }_{\mathrm{t}}$为湍流黏度，${\sigma }_{k}$、${\sigma }_{\omega }$、$\gamma$、${a}_{1}$、$\beta$、${\beta }^{*}$均为常数，$S$为应变速率，${F}_{1}$与${F}_{2}$均为控制常数的加权函数。

1.4   空化模型

由于航行体下落入水过程中往往伴随着空化现象，所以本文采用通过简化Rayleigh-Plesset空化动力学方程^[24]得到的Schnerr-Sauer空化模型来模拟空化现象，其水蒸气输运方程为：

式中：${\dot{m}}^{+}$与${\dot{m}}^{-}$分别为蒸发源项与冷凝源项，${R}_{\mathrm{b}}$为气泡半径，$p$为环境压力，${p}_{\mathrm{v}}$为饱和蒸气压。

1.5   数值波浪水池模型

为了得到与真实海浪接近的波浪环境，本文采用Stokes五阶非线性波理论构建数值波浪水池，具体采用Fenton^[25]通将流函数展开为傅里叶级数形式并将波陡作为摄动参数进行展开所得到的五阶波的解，其波面方程与速度势表达式分别如下：

式中：$\eta \left(x\right)$为波面高程，$\varphi \left(x,y,t\right)$为速度势，$h$、$c$、$\kappa$、$\varepsilon$分别为、水深、波速、波数、无量纲波高，$\bar{u}$为任意高程处一个波长范围内的平均流速，${C}_{i}$、${A}_{ij}$、${B}_{ij}$为无量纲系数，$O\left({\varepsilon }^{6}\right)$为${\varepsilon }^{6}$的同阶无穷小量。

基于上述波浪理论，本文通过wave forcing^[26]方法实现波浪水池的数值造波与消波，即通过在控制方程中加入一源项，使得从距离航行体一定距离处开始，由控制方程得到的流动变量逐渐趋近于波浪模型的理论解，直到在流域边界处，所得的流动变量完全等于波浪模型的理论解。所添加的源项为

式中：$\rho$为流体密度，$S$为输运方程的当前解，${S}^{*}$为wave forcing求解所接近的值，$\alpha$为力系数，定义如下：

式中：$d$为距边界距离，${d}_{1}$为应用wave forcing的区域宽度，${\alpha }_{0}$为边界处的力系数。

2.   数值方法验证与模型建立

2.1   下落入水冲击数值方法验证

本文通过CFD方法构建了带助浮装置航行体下落入水冲击的数值模型：利用VOF多相流模型捕捉不同流体之间的交界面，利用k-ω SST湍流模型捕捉湍流流动细节，利用Schnerr-Sauer空化模型实现空化模拟，利用重叠网格与滑移网格相结合的方法实现航行体下落冲击的六自由度运动与助浮装置的转动展开过程。为验证以上数值方法的准确性，本文通过数值方法对Wei等^[27]的下落圆柱试验进行复现，并将数值模拟得到的下落圆柱的入水空泡形态与下落位移时历曲线结果与试验结果进行对比。

图1给出了试验与数值得到的圆柱入水过程中空泡形态的发展情况，可以看出本文数值方法能有效捕捉到物体入水过程中空泡的生成、发展与拉断现象，对比可知试验结果和数值模拟结果基本吻合。

图  1  圆柱入水空泡形态试验与数值结果对比

Figure  1.  Comparison of experimental and numerical results of cylindrical water-entry cavitation

下载: 全尺寸图片 幻灯片

图2给出了试验与数值得到的圆柱入水冲击过程中质心位置时历变化情况。经过计算本文数值算法结果与试验结果之间在不同时刻的误差均保持在5%以下，因此本文构建的下落入水冲击数值方法可以较好预示物体下落冲击过程。

图  2  圆柱下落冲击过程质心位置时历曲线对比

Figure  2.  Comparison of the time history curves of centroid position during cylinder falling impact

下载: 全尺寸图片 幻灯片

2.2   数值波浪水池方法验证

本文基于Stokes五阶非线性波理论，通过wave forcing方法进行了波浪水池的构建，为了验证数值水池的有效性与稳定性，本文进行了数值波浪水池的方法有效性验证。在有效性验证中构建了如图3所示的长60 m，宽7 m，高15 m的立方体数值造波水池，水池入口，出口以及底面均设为速度入口边界，顶部设置为压力出口边界，前后设置为对称平面。

图  3  数值造波水池几何模型

Figure  3.  Numerical wave pool geometry model

下载: 全尺寸图片 幻灯片

基于构建完成的数值波浪水池对波高H=1 m、周期T=3 s、水深h=10 m的海况进行数值模拟，模拟中保证波高方向单位网格尺度为H/32，波长方向与垂直波长方向单位网格尺度为λ/120，λ为波长。设定距离入口或出口不足15 m的区域为消波区，中间30 m的区域为有效计算区，每隔5 m设置一个浪高仪对波高时历进行监测，具体布置如图4所示。

图  4  浪高仪布置情况

Figure  4.  The layout of waveprobes

下载: 全尺寸图片 幻灯片

图5给出了t=35 s时刻计算域中波浪波高分布的数值与理论对比情况，图6给出了x=20 m处与x=40 m处的波高时历曲线数值与理论对比情况，对比结果表明本文数值波浪水池方法在造波稳定性与准确性方面都满足后期开展波浪环境下带助浮装置航行体冲击入水研究的条件。

图  5  t=35 s时刻计算域中波浪波高分布情况

Figure  5.  Wave height distribution in the calculation domain at time t=35 s

下载: 全尺寸图片 幻灯片

图  6  x=20 m处与x=40 m处波高时历曲线

Figure  6.  Time history curve of wave height at x=20 m and x=40 m

下载: 全尺寸图片 幻灯片

2.3   带助浮装置航行体下落入水冲击数值模型构建

本文通过将航行体与助浮装置作为同一个六自由度体来模拟助浮装置与航行体的虚拟连接，实现了助浮装置与航行体运动绑定与受力绑定，助浮装置与航行体绑定运动的同时利用滑移网格的方式实现了助浮装置在冲击入水过程中的转动展开过程。计算中对航行体与助浮装置之间的连接杆进行了忽略，模拟实现了助浮装置与航行体之间的钢性连接，固定了航行体与助浮装置之间的相对距离，由于连接杆相较于助浮装置与航行体体积来说其直径较小，其对流场的扰动可以忽略不计，所以计算过程中忽略了连接杆的实际体积。计算采用的航行体模型与助浮装置模型如图7所示，航行体总长为L，第一级直径为D。助浮装置采用椭球体结构，椭球长轴半径为R，短轴为长轴的一半。助浮装置分布在航行体两侧，助浮装置椭球中心距离航行体头部为2.34D，距离航行体中轴线为2.56D。

图  7  计算模型几何尺寸

Figure  7.  The geometric dimensions of computational model

下载: 全尺寸图片 幻灯片

图8给出了波浪环境下带助浮装置航行体的下落冲击过程计算域设置与网格划分，计算域为长60D、宽24D、高24D、深14D的方形波浪水池，计算域边界设置与数值造波水池边界设置保持一致，航行体与助浮装置壁面设置为无滑移壁面边界条件。本文计算设置有大地坐标系xyz与航行体坐标系x₀y₀z₀，其中大地坐标系z轴方向为波高方向，y轴为波长方向，x轴为波浪展向，航行体坐标系原点设定在航行体质心处，x₀轴与y₀轴为航行体径向方向，z₀轴为航行体轴向方向。距离入口20D区域与距离出口20D区域设定为消波区，航行体质心前后10D区域为有效计算区域。计算域整体网格分为背景区域与重叠网格区域，背景区域中设置波面加密区与航行体下落路径加密区，重叠区与背景区通过重叠网格交界面连接，重叠区内航行体与助浮装置壁面设置棱柱层网格，助浮装置转动展开过程通过滑移网格方法实现，助浮装置区域网格通过滑移网格交界面与外部网格进行求解过程中的信息交互。

图  8  计算域设置与网格划分方案

Figure  8.  Computing domain setup and meshing scheme

下载: 全尺寸图片 幻灯片

助浮装置在航行体下落过程中进行转动展开，转动围绕助浮装置旋转中心进行，助浮装置初始触水时开始转动，转动角度按照正弦曲线规律变化，当旋转角度达到60°时停止转动，图9给出了助浮装置转动方向示意与转动角度的时历曲线变化，整个转动展开过程从0.405 s时开始，0.804 s时结束，展开完成后助浮装置始终保持与航行体中轴线成60°夹角。

图  9  助浮装置转动方向示意与转动角度变化时历曲线

Figure  9.  Boost floatation aids rotation direction and time history curve of rotation angle

下载: 全尺寸图片 幻灯片

本文在波浪环境下的带助浮装置航行体下落冲击模拟中设定航行体初始时刻中轴线与z轴轴线成15°夹角，中轴线与航行体尾部交点距离波峰点0.2 m距离，下落初始速度设定为沿z轴逆向7.5 m/s，如图10所示。航行体水下作业试验中，由于航行体尾部助推动力耗尽，航行体在空中进入自由下落过程，该过程中航行体头部在上，尾部在下，其尾部往往最先接触水面，因此本文研究航行体尾部向下入水过程对工程指导更具有针对性，并且航行体尾部构型相比于头部流线型来说更加具有复杂性，其流场演化更加剧烈，入水冲击过程更加复杂，这使得选取尾部入水冲击过程开展研究可以对入水冲击流场进行更加深入的挖掘分析。随着初始下落航行体中轴线与z轴所成夹角的减小，航行体最大入水深度逐渐增加，并在竖直落水时达到峰值^[22]，因此为了综合考虑航行体姿态对入水过程的影响以及捕捉较极端危险工况，本文计算设定航行体初始时刻轴线与z轴轴线呈15°夹角。

图  10  航行体初始下落位置

Figure  10.  Initial drop position of vehicle

下载: 全尺寸图片 幻灯片

为了便于后续结果的分析比对，对相关参数进行无量化：

式中：${v}_{0}$为航行体初始下落速度，${v}_{0}$=7.5 m/s；$L$=3.25 m；$y$为实际横向位移；z为实际垂向位移；$D$为航行体第一级直径；${v}_{{\textit{z}}}$为实际垂向速度；$g$为重力加速度，本文取9.81 m/s²；${v}_{y}$为实际横向速度；${a}_{{\textit{z}}}$为实际垂向加速度；${a}_{y}$为实际横向加速度；$p$为实际压力；${p}_{0}$为一个大气压压强值，${p}_{0}$=101325 Pa。

通过改变重叠区网格与背景加密区网格大小，形成总体网格数量分别为7.5×10⁶、6.3×10⁶、5.3×10⁶的细、中、粗三套网格划分方案，第一层网格大小设置为2.5×10⁻⁴ m，计算得到实际y+值保持在0～150之间，计算采用全y+壁面处理模型。基于三套网格划分方案对波高H=1 m，周期T=3 s的特定海况带助浮装置航行体下落入水过程进行计算，计算结果如图11所示，对比分析可知三套网格计算结果最大相对误差均保持在5%以下。综合考虑模拟计算精度，稳定性以及计算效率，本文模拟研究选择网格数量为6.3×10⁶的中网格划分方案。

图  11  网格无关性计算结果

Figure  11.  Grid independence calculation results

下载: 全尺寸图片 幻灯片

3.   计算结果与分析

3.1   波浪环境下航行体下落运动过程分析

对波高H=1 m，周期T=3 s的特定海况下的航行体下落入水过程进行分析。如图12所示，航行体带助浮装置下落冲击过程具体经历了下落冲击段、下沉段、上浮段以及下落-上浮过程衰减段，其中衰减段中航行体竖直位移与竖直速度逐渐趋近于稳定值，不再发生较大变化，这也表征着带助浮装置航行体的下落冲击过程结束。第一次的下落-下沉-上浮过程由于航行体具有一定的初始速度，在重力的作用下进入短暂加速阶段，加速结束后其下落速度达到最大值${\bar{v}}_{{\textit{z}}}$=−1.465，所以波浪对其运动状态没有明显的影响。当$\bar{t}$=2.52时，航行体下落速度减为零，达到下落最大位移$\bar{{\textit{z}}}$=−1.351，此后航行体进入首次上浮阶段。在上浮阶段中向上的浮力使得航行体具有向上的加速度，在浮力加速度的作用下航行体再次越出水面并于$\bar{t}$=5.47时达到最大上浮位置，紧接着航行体在重力的作用下再次下落，即进入衰减运动段。进入运动衰减过程后航行体在竖直方向速度较小，因此衰减阶段波浪对其的运动影响显著。

图  12  航行体下落-下沉-上浮过程

Figure  12.  Vehicle fall-sink-rise process

下载: 全尺寸图片 幻灯片

图13与图14分别给出了航行体的位移与速度曲线，图15给出了航行体冲击下落过程的速度方向矢量发展情况。可以发现当航行体进入运动衰减段后其运动状态与波浪具有强相关，特别是航行体横向运动状态与波浪传播相位直接相关，波浪的作用使得航行体的运动状态偏移了物体静浮运动的震荡衰减模式，其横向运动呈现随波性的摇荡模式。在航行体下落初期的短暂加速过程中，竖直速度曲线初期出现了两个速度突变时刻，分别为$\bar{t}$=0.08与$\bar{t}$=0.19：$\bar{t}$=0.08为航行体底部触水时刻，底部接触水面形成了短暂的反向砰击力，这使得航行体下落加速阶段出现了短暂的减速；$\bar{t}$=0.08为航行体助浮装置触水时刻，触水形成的反向砰击力导致了航行体第二次减速，并且由于触水面积大，其触水砰击力相较于航行体尾部触水时更大，持续时间更长，这也使得航行体下落加速阶段结束，进入了下落减速阶段。航行体第二次下落进入运动衰减段并于$\bar{t}$=6.23时航行体下落速度达到最大，此时下落速度仅为${\bar{v}}_{{\textit{z}}}$=−0.189 m/s，在助浮装置浮力的作用下，航行体开始减速下沉，当$\bar{t}$=7.59时航行体下落速度再次减为零并开始上浮，此时由于上浮位置处于波谷与波峰之间的波后相位，航行体与助浮装置受到了具有竖直向下速度的水质点的作用，上浮速度受到极大抑制，上浮仅持续到$\bar{t}$=8.55后航行体再次进入下落阶段，此后航行体的浮动主要受波浪影响，即进入随浪漂移浮沉状态。

图  13  航行体位移时历曲线图

Figure  13.  Time history curves of vehicle’s displacement

下载: 全尺寸图片 幻灯片

图  14  航行体速度时历曲线

Figure  14.  Time history curves of vehicle’s velocity

下载: 全尺寸图片 幻灯片

图  15  航行体冲击下落过程速度方向矢量

Figure  15.  Velocity vectors of vehicle during impact fall process

下载: 全尺寸图片 幻灯片

由于航行体下落初始不具有水平速度，所以航行体在水平方向上的运动受到初始释放角度与波浪水质点速度的影响较为显著。航行体下落时尾部砰击水面，入水存在一定的倾斜角度，所以沿航行体轴线的砰击力不仅仅产生了竖直方向上的冲击力，而且还产生了水平+y方向的推动力，加之航行体初始下落时会受到波峰处具有水平向右速度的水质点的作用，航行体初始下落时将产生+y方向位移，而随着航行体下落深度的增加，航行体周围的波浪相位发展为了波谷相位，波谷处具有−y方向速度的水质点使得航行体水平位移方向变为了−y方向。当航行体第二次下落进入运动衰减段时，由于下落位置正好处于波浪波峰相位，航行体尾部受到了与波浪传播方向相同即具有+y方向运动速度的水质点的冲击作用，其航行体−y方向速度开始减小并产生了+y方向速度，此后航行体水平速度随着所处波浪相位的不同而改变，而水平位移则呈现整体向+y方向运动，局部短时向−y方向运动的发展趋势。

图16给出了航行体冲击下落过程中的绕x₀轴转角变化情况，规定绕x₀轴逆时针转动为正向转角。初始下落过程中航行体尾部砰击力与航行体重力形成了沿x₀轴顺时针方向的偏转力矩；随着波浪的发展，航行体与助浮装置完全入水后，波谷处具有水平逆向速度的水质点作用于航行体与助浮装置，由于助浮装置布置位置相较于航行体质心偏上，所以波浪水质点的冲击使得航行体发生沿x₀轴逆时针方向的偏转；航行体在助浮装置浮力作用下第一次上浮，由于惯性，助浮装置再次出水；在下落过程中，航行体尾部段与波浪波峰遭遇，这使得航行体沿x₀轴逆时针偏转加剧，当$\bar{t}$=7.28时航行体偏转角度达到了沿x₀轴逆时针58.048°；此后航行体尾部与头部先后经历了波浪的波谷与波峰相位，使得航行体发生了沿x₀轴顺时针方向的偏转；最终航行体速度衰减达到一个较小水平，运动主要受波浪影响，转角随着波浪传播呈现周期性变化。

图  16  航行体转角时历曲线

Figure  16.  Time history curve of vehicle’s angle

下载: 全尺寸图片 幻灯片

3.2   波浪环境下航行体下落冲击过程分析

当$\bar{t}$=2.52时航行体下沉至最深位置，此后开始上浮，本节重点对$\bar{t}$=2.52之前的带助浮装置航行体下落冲击过程进行分析。图17为航行体下落冲击过程中竖直方向与水平方向上受到的冲击加速度。

图  17  带助浮装置航行体加速度时程曲线

Figure  17.  Acceleration time history curve of vehicle with boost floatation aids

下载: 全尺寸图片 幻灯片

图18给出了下落冲击过程中入水空泡的发展情况，由于冲击过程持续时间仅为波浪周期的1/3，所以冲击阶段波浪因素对航行体影响较小。

图  18  入水空泡发展过程

Figure  18.  Development process of water-entry cavity

下载: 全尺寸图片 幻灯片

航行体倾斜下落过程中，尾部喷口先发生入水砰击，喷口尾部端面砰击力分解到竖直与水平方向，体现为航行体下落过程中的下落阻力与水平正向推力，该砰击过程从尾部接触水面开始持续到了$\bar{t}$=0.10。尾部砰击结束后，当$\bar{t}$=0.46之前航行体两侧助浮装置与水面发生砰击，助浮装置冲击由于与水面的接触面积大，所以其竖直冲击力相较于尾部竖直冲击力要大得多。图19为下落过程中助浮装置的受力情况，其中助浮装置1与助浮装置2分别位于航行体两侧的−x方向与+x方向。结合下落加速度与助浮装置受力情况进行分析，带助浮装置航行体的竖直下落加速度变化情况与助浮装置z轴方向上的受力均出现了两次峰值，这说明航行体下落过程中竖直方向上的冲击主要是由于助浮装置的冲击与转动造成的，第一次峰值是因为助浮装置与水面发生砰击，产生了冲击加速度，第二次峰值是因为助浮装置展开过程中呈现先加速后减速的运动变化，展开速度在$\bar{t}$=0.69达到最大，这也是第二次加速度峰值出现的时间段，在加速展开过程中助浮装置加速拍水，水介质的反作用力直接导致了带助浮装置航行体在竖直与水平方向上的加速度高峰，$\bar{t}$=0.69后助浮装置展开速度减小，展开的角加速度方向逆转，航行体的加速度也相应减小。与竖直方向不同的是，在水平方向上带助浮装置航行体的加速度变化情况与助浮装置受力情况明显不同，这说明在水平方向上带助浮装置航行体的加速度变化情况不仅仅受到助浮装置冲击过程影响，而且与入水空泡的发展过程直接相关。入水过程中，带助浮装置航行体的水平加速度出现了3次正向峰值：当$\bar{t}$=0.57时，助浮装置砰击水面形成了第一次峰值；当$\bar{t}$=0.71时，助浮装置加速拍水形成了第二次峰值，也是水平正向最大的加速度峰值；当$\bar{t}$=0.93左右，带助浮装置航行体入水空泡脱离航行体背流侧进入大面积溃灭阶段，溃灭射流冲击航行体背流侧，这直接造成了带助浮装置航行体水平方向上的第三个正向加速度峰值。助浮装置1与助浮装置2在竖直与水平方向上的受力变化情况大致相同，在x轴方向上，两侧助浮装置均受到了背离航行体向外的作用力，并且该作用力峰值相较于其他方向上的峰值要大得多，这使得航行体与助浮装置之间的连接机构在设计强度上必须考虑该极限载荷情况。

图  19  助浮装置受力时历曲线

Figure  19.  Time history curves of boost floatation aids’ force

下载: 全尺寸图片 幻灯片

为了进一步对航行体冲击入水空泡发展情况进行研究，设定如图20所示的压力监测点，测点分为迎流端，背流端和侧流端三组，其中迎流端为航行体倾斜入水时的先触水侧，侧流端设定在助浮装置与航行体内侧，每组各包含5个压力监测点（G1、G2、···、G5），每个监测点之间距离0.5 m，距离航行体壁面0.00125 m，第1监测点G1与航行体质心位置的轴向距离为1.25 m。图21给出了各监测点的压力时历变化情况。

图  20  压力监测点布置位置

Figure  20.  Location of pressure monitor points

下载: 全尺寸图片 幻灯片

图  21  监测点压力时程曲线

Figure  21.  Time history curves of the pressure at monitor points

下载: 全尺寸图片 幻灯片

上述结果表明迎流端压力测点最先出现压力峰值，出现峰值的时刻对应了带助浮装置航行体水平加速度曲线的下落谷底出现时刻，这进一步说明了迎流侧入水空泡的溃灭是航行体出现水平加度逆向峰值的主要原因。从图21(a)中不难发现每个压力测点均出现了两次压力局部峰值波动，在$\bar{t}$=0.5到$\bar{t}$=0.8之间的第一次压力峰值是迎流侧入水空泡的溃灭所引起的；在$\bar{t}$=0.95的第二次压力峰值则是因为入水空泡逐渐溃灭至背流端，流体充满原来迎侧流端的空泡位置所引起的压力冲击。压力监测点越靠近航行体尾部其压力峰值越大且出现时刻越早，说明迎流端入水空泡溃灭发展是自下而上进行的，并且溃灭强度逐渐减弱。背流端与侧流端压力峰值均出现在$\bar{t}$=0.95左右，但是侧流端压力曲线较背流端出现扰动的时间要更早，这说明入水空泡沿迎流端-侧流端-背流端方向发展溃灭，并且在背流端溃灭阶段其溃灭最剧烈，溃灭压力最大，溃灭压力随深度增加而增加。

3.3   不同海况环境对航行体下落冲击过程的影响分析

为了研究不同海况环境下带助浮装置航行体冲击下落过程，本文设定了如表1所示的计算工况，其中工况1到工况4表示不同的海况环境，工况5为静水工况。

表  1  工况设置

Table  1.  Working conditions setting

工况	波高/m	波周期/s
1	0.50	3.0
2	0.75	3.0
3	1.00	3.0
4	1.25	3.0
5	0.00	/

下载: 导出CSV 

| 显示表格

各工况下带助浮装置航行体以相同的水面高度与初速度释放，图22给出了不同工况下航行体的下落-下沉-上浮过程，图23与图24分别给出了不同工况下带助浮装置航行体的位移与速度情况。从不同海况环境下的航行体运动状态来看，波浪环境对航行体下落冲击过程的影响主要体现在了运动衰减阶段，而在首次下落冲击阶段相对较小。从竖直位移来看，海况浪级越高，航行体在下沉过程中达到的深度越大，进而在上浮过程中达到的高度也越小，这是因为在下沉过程中航行体遭遇到了波峰与波谷之间的相位，该相位下水质点主要具有竖直向下的速度，并且海况等级越高，水质点速度越大，所以航行体下沉深度也越大，由于海况浪级高的环境下航行体下沉深度越大，所以导致了航行体在上浮过程中消耗的能量也越多，上浮高度也越低，在运动衰减段，海况浪级越低其航行体浮沉运动越接近静浮衰减规律。从竖直速度来看，航行体下落冲击过程中海况浪级越高其竖直速度也越大，而进入运动衰减段后，海况浪级越低，航行体速度反而越大，这是因为在带助浮装置航行体第二次下落过程中遭遇到了波谷到波峰之间的波浪相位，水质点竖直向上的速度导致了高海况浪级下的航行体速度大幅衰减。相较于竖直方向的运动响应，波浪环境对航行体的水平运动响应影响更明显，海况浪级越高的环境，航行体水平位移的幅值和逆向漂浮位移越明显，不同时刻的水平速度也越大。

图  22  不同工况下航行体下落-下沉-上浮过程

Figure  22.  Vehicle fall-sink-rise process under different working conditions

下载: 全尺寸图片 幻灯片

图  23  不同工况下航行体位移时历曲线

Figure  23.  Displacement history of the vehicle under different working conditions

下载: 全尺寸图片 幻灯片

图  24  不同工况下航行体速度时历曲线

Figure  24.  Velocity history of vehicle under different working conditions

下载: 全尺寸图片 幻灯片

图25给出了带助浮装置航行体下落冲击过程中的冲击加速度情况，不同海况环境下的竖直加速度总体变化浮动较小，尾部砰击结束到助浮装置砰击开始阶段，海况浪级与竖直加速度保持正相关关系，该区间波浪水质点的竖直分量加快了航行体的下落过程。不同工况航行体尾部砰击加速度变化不大，助浮装置砰击与助浮装置加速展开形成的砰击力随着海况等级的提高而略微减小，以静浮工况加速度为基准，工况1助浮装置砰击加速度占比静浮工况加速度为93.76%，随着工况波高每增加0.25 m，助浮装置砰击加速度按平均4.28%的占比速度减小；工况1助浮装置加速展开形成的砰击加速度占比静浮工况加速度为97.42%，随着工况波高每增加0.25 m，助浮装置加速展开砰击加速度按平均0.94%的占比速度减小。水平方向上的砰击加速度受到波浪影响较大，随着海况浪级的增加，助浮装置触水砰击后会形正向水平加速度的短时震荡，当助浮装置入水开始加速展开时，水平加速度的峰值随着海况浪级的增加而增加，波高每增加0.25 m，加速度峰值平均增加15.31%，而由于背流端入水空泡溃灭所造成的第二个正向水平加速度峰值则会随着海况浪级的增加而减小。

图  25  不同工况下带助浮装置航行体加速度时历曲线

Figure  25.  Acceleration time history curves of vehicle with boost floatation aids under different working conditions

下载: 全尺寸图片 幻灯片

不同海况环境对带助浮装置航行体水平方向的加速度冲击主要体现在了对入水空泡的影响上。图26为$\bar{t}$=0.92 s时不同海况环境下的计算域中纵剖面水的体积分数云图以及速度方向矢量图，入水空泡从迎流端向背流端进行溃灭，当$\bar{t}$=0.92 s时航行体入水位置处于波浪波峰位置，可以发现海况浪级越高其航行体入水空泡溃灭闭合时间越晚，这是因为波浪水质点的水平速度与入水空泡的溃灭方向相反，这极大减缓了入水空泡从迎流端到背流端的溃灭进程。

图  26  不同工况下入水空泡流场中截面示意图

Figure  26.  Schematic diagram of mid-profile in the flow field of water-entry cavity under different working conditions

下载: 全尺寸图片 幻灯片

4.   结　论

本文通过CFD数值模拟方法对波浪环境下的带助浮装置航行体下落冲击过程进行了探究，主要研究了带助浮装置航行体下落冲击过程中的运动学以及动力学参数变化特性以及空化流场演化过程，以此为基础具体分析了航行体以及助浮装置在下落入水冲击过程中的冲击特性以及运动演化，具体研究了带助浮装置航行体下落冲击过程在不同海况等级下的具体表现，主要得到了以下结论。

(1)按照运动状态的不同，带助浮装置航行体下落冲击过程具体分为了下落冲击段、下沉段、上浮段以及下落-上浮过程衰减段，其中当航行体进入运动衰减段后其运动状态体现出了与波浪的强相关性。

(2)带助浮装置航行体下落冲击过程中主要经历了航行体尾部入水砰击、助浮装置入水砰击以及助浮装置旋转展开过程。竖直方向的冲击中助浮装置加速展平阶段产生的加速度峰值最大，脉宽最大；水平方向的冲击中，助浮装置加速展平阶段产生的水平加速度峰值最大，并且入水空泡背流端溃灭冲击也将产生较大的加速度峰值。在助浮装置入水冲击过程中，助浮装置x轴方向会产生较大冲击力，威胁助浮装置与航行体之间的连接结构。

(3)不同海况对带助浮装置航行体下落冲击过程的影响主要体现在运动衰减阶段，而在首次下落冲击阶段相对较小。相对于竖直方向的冲击，水平方向的冲击受到不同海况下波浪环境的影响较大，不同海况对航行体的水平冲击造成的影响主要是通过影响入水空泡的形成与溃灭过程实现的。

(4)不同海况对航行体尾部砰击加速度影响远远小于助浮装置砰击影响，其中在竖直方向上，随着海况波高每增加0.25 m，助浮装置砰击加速度按平均4.28%占比的速度减小，助浮装置加速展开砰击加速度按平均0.94%占比的速度减小；在水平方向上，随着海况波高每增加0.25 m，助浮装置加速展开砰击加速度峰值平均增加15.31%。