0 引　言

高马赫数弹体入水，是一个气液固三相耦合问题，涉及到不同介质之间的相互转换以及由于高速弹体运动而产生的激波现象，其中多相流场之间的相互转换、流体可压缩性等影响因素，一直是入水问题研究的难点^[1]。高速运动的弹体入水后，其表面出现大量水蒸气泡，将弹体紧紧围绕，形成超空泡现象^{[2 − 3]}。弹体入水过程具有瞬态非定常及大载荷特点，对弹体的运动轨迹和受力特性产生严重影响。

国外对弹体入水问题的研究较早，从实验和数值模拟2个方面进行了大量研究。

Grumstrup等^[4]结合数值模拟和实验研究2种方法对弹体入水过程中超空泡壁的振动规律及其内部的特征进行了研究。Sudo等^[5]采用实验方法研究了涂有磁流体的永磁弹体入水时的空泡动力学特性，用高速摄像机观察到了永磁体弹体入水过程中空泡的产生过程，获得了磁流体产生的交变磁场对空泡动力学特性的影响规律。Bodily等^[6]利用高速摄像机对细长轴对称弹体的弹头几何形状、润湿方法、入水倾斜角如何影响弹体的受力和运动轨迹开展了实验研究。实验结果表明，在垂直入水和相同入水深度的条件下，弹头对称润湿的弹体侧向位移比弹头半疏水-半亲水的弹体侧向位移大，锥形弹头的弹体比椭圆形弹头和扁平形弹头的弹体侧向位移大；垂直入水时，3种弹头几何构型中，扁平形弹头的弹体受到的冲力最大；斜入水时，扁平形弹头的弹体相对于垂直入水时，弹体受到冲击力会减小。Erfanian等^[7]采用流-固耦合数值模拟方法计算了弹体以不同攻角低速入水的过程，获得了弹体入水过程中的超空泡形状和弹体的运动轨迹等结果，并与实验进行了对比。Derakhshanian等^[8]利用钝性体垂直入水和斜入水的实验数据，验证了3种数值模拟方法的可靠性，研究结果表明采用ALE的流固耦合方法是一种解决钝性体入水问题的可靠方法。Mirzaei等^[9]建立了一个圆柱弹体入水时预测超空泡形状、弹体姿态和弹道的瞬态动力学数学模型，并以几种不同入水角度的弹体为对象，通过数值仿真和物理实验对瞬态动力学模型进行了验证。Akbari等^[10]利用CFD方法建立了在入水初速度为280 m/s时的3种长径比的圆柱弹体，在3种不同入水角度条件下的数值计算模型，获得了3种入水角度下防止圆柱弹体翻转的临界长径比，并通过前人的实验结果进行验证。Yoo等^[11]针对弹体的入水问题，基于多相流弱可压缩SPH方法建立了弹体入水多相流固耦合数值计算模型，研究了球体润湿特性对空泡的形成和形状的影响规律，并通过物理实验对所获得的结论进行验证。

在国内，自21世纪初期才开始对超空泡高速弹体的相关技术进行深入研究，尤其是最近十几年来，越来越多的专家学者在该领域跟踪国内外前沿技术，开展了相关实验研究和数值模拟。路中磊等^[12]通过实验开展了圆柱空心壳体低速垂直入水过程的研究，发现入水速度和空腔壳体的几何结构对入水空泡形状有较大影响。周杰等^[13]采用二级轻气炮和流场显示系统通过实验方法研究了高速弹丸斜侵彻入水过程中气/水界面变形破碎规律、入水空泡和水中冲击波传播的特性，获得了弹丸斜侵彻入水流场的演化图像，弹丸头部的构型对气泡的体积大小和水下弹道的稳定性有较大影响的结论。孙钊等^[14]采用数值仿真模拟的方法，得出了半疏水-半亲水球体垂直入水后的运动轨迹会偏离竖直的运动轨道。孙铁志等^[15]通过数值模拟和实验对圆柱体倾斜入水空泡的演化过程和机理进行了研究，获得了围绕圆柱体存在多尺度流体漩涡的结论；该团队还开展了多相流的可压缩性对空泡演化过程的影响研究，得出了多相流的可压缩性对空泡的演化和空泡内的压力有重要影响 ^[16]。

目前，由于实验条件的限制，针对弹体入水的研究，很难实现弹体以高马赫数的初始速度条件入水。但弹体高速入水时，空泡的形成、演化和溃灭的发展过程以及弹体所受到的载荷与弹体低速入水时的情况完全不同^[17]。本文考虑了不同初始速度、流体的可压缩性、相变影响、空泡演化等问题，在建立了多相流固耦合计算模型的基础上，仿真计算了多种初始速度条件下的锥形弹体垂直入水过程中弹体的动力学特性，包括弹体的运动轨迹、空泡的演化规律及流体速度和压力的变化规律等。

1 高马赫数弹体入水过程数学模型 1.1 流体控制方程

高速弹体入水过程是典型的多相流固耦合过程，不仅需建立高速弹体的控制方程，而且也需建立流体的控制方程。

能量守恒方程：

$ \frac{\partial }{{\partial t}}\int_\Omega {\rho {e_t}{\mathrm{d}}V} + \int_\Gamma {\rho {e_t}{{v}}{\mathrm{d}}S} = \int_\Omega {\rho {{bv}}{\mathrm{d}}V} + \int_\Gamma {\sigma v{\mathrm{d}}S} 。$

(1)

式中：$ {e_t} $为单位质量的总能量。

质量守恒方程：

$ \frac{\partial }{{\partial t}}\int_\Omega {\rho {\mathrm{d}}V} + \int_\Gamma {\rho v{\mathrm{d}}S} = 0 。$

(2)

式中：t为时间；ρ为流体密度；v为速度。

动量守恒方程：

$ \frac{\partial }{{\partial t}}\int_\Omega {\rho v{\mathrm{d}}V} + \int_\Gamma {\rho v{{(v}}{\mathrm{d}}S{\text{)}}} = \int_\Omega {\rho {{b}}{\mathrm{d}}V} + \int_\Gamma {{t^{{{(n)}}}}{\mathrm{d}}S} 。$

(3)

式中：$ {t^{{{(n)}}}} $为边界上单位面积上的力；b为体积力。

1.2 弹体控制方程

计算中需考虑弹体与流体间的相互作用，因此需建立弹体的控制方程。弹体入水过程中承受的外载荷是随时间变化的，为考虑其动力学特性，平衡方程需考虑惯性力和阻尼力，具体形式如下：

$ {\boldsymbol{M{{a}} }}+ {\boldsymbol{C{{v}} }}+{\boldsymbol{K{{d}} = }}{\boldsymbol{P}}。$

(4)

式中：M为弹体的质量矩阵；a为弹体的加速度矢量矩阵；C为阻尼矩阵；v为弹体速度矢量矩阵；K为刚度矩阵；d为弹体位移矢量矩阵；P为随时间变化的载荷矢量矩阵，可通过直接积分法对方程进行求解。

1.3 流体可压缩性方程

弹体高速入水过程中，高速弹体会带动弹体表面附近的流体高速流动，因此，高马赫数弹体入水过程中的流体必考虑可压缩性。Tait状态方程能较好描述流体高速流动和高压下的流体行为，其方程的简化形式为：

$ {\left( {\frac{\rho }{{{\rho _0}}}} \right)^n} = \frac{{\boldsymbol{K}}}{{{{\boldsymbol{K}}_0}}}，$

(5)

$ {\boldsymbol{K}} = {{\boldsymbol{K}}_0} + n\Delta {\boldsymbol{p}} ，$

(6)

$ \Delta {\boldsymbol{p}} = {\boldsymbol{p}} - {{\boldsymbol{p}}_0}。$

(7)

式中：$ {p_0} $为大气压力，作为参考压力值；$ p $为当前压力值，Pa；$ {\rho _0} $为水的密度，取值为1000 kg/m³；$ {K_0} $为水的体积模量，取值为2.2 GPa；n为密度指数常数，取值为7.15。

声速$ c $为压力扰动传播速度：

$ {\boldsymbol{c}} = \sqrt {\frac{{\boldsymbol{K}}}{\rho }}。$

(8)

式中：$ \rho $为液体密度，kg/m³；$ K $为液体体积模量，Pa。

1.4 流-固耦合求解方法

显式时间积分求解法即可用来求解拉格朗日方程，也可用来求解欧拉方程。

设当前时间步为第n步，显式求解法会将式(9)改写为式(10)：

$ {\boldsymbol{M{{{a}}}}_{{n}}}{{ + }}{\boldsymbol{C{{{v}}}}_{{n}}}{{ + }}{\boldsymbol{K{{{d}}}}_{{n}}}{{ = }}F，$

(9)

$ {\boldsymbol{M{{{a}}}}_{{n}}}{\text{ = }} F - T 。$

(10)

式中：F为外载荷；T为内力，$ T{{ = }}{\boldsymbol{C{{{v}}}}_{{n}}}{{ + }}{\boldsymbol{K{{{d}}}}_{{n}}} $。

弹体加速度可通过质量矩阵求逆再乘上剩余力矢量$ {F^0}{\text{ = }}F - T $求得：

$ {{{{\boldsymbol{a}}}}_{{n}}}{{ = }}{{\boldsymbol{M}}^{ - 1}}{F_{{n}}}^0。$

(11)

式中：M为对角阵，由于质量分布在节点上，式(11)可写成不同自由度下独立的一元一次方程，因而加速度可表示为：

$ {{{{\boldsymbol{a}}}}_{{{ni}}}}{{ = }}{F_{{{ni}}}}^0{{/}}{{\boldsymbol{M}}_{{i}}}。$

(12)

2 高速弹体流-固耦合数值模型验证

为了验证所采用的流固耦合计算方法的可靠性，以文献[1]所给出的弹体为例，对本文采用的数值计算方法进行验证。建立直径为10 mm的球头回转体二维固体几何模型和周围流体几何模型，划分流体网格和弹体网格，流固耦合边界施加在弹体与流体的交界面上，其余流体边界条件与文献[1]相同，弹体约束所有位移自由度，沿x正向给定20 m/s水流速度，观察超空泡形态，计算结果如图1所示。

通过对比可知，图1(b)与图1(a)中的空泡外形几乎一致。图1(c)为通过仿真得到的空泡轮廓和实验得出的空泡轮廓对比，这里以弹体顶部圆心为原点，弹体轴向为x方向，径向为y方向，绘制弹体表面超空泡的轮廓。如图1(c)中曲线所示，仿真结果与文献[1]中的弹体表面空泡轮廓几乎重合，平均误差仅为2.7%。在此情况下认为该计算模型有效，能有效模拟弹体入水形成的超空泡现象。

3 高马赫数弹体入水过程流-固耦合数值计算 3.1 数值计算模型

弹体几何模型如图2所示。

流体域模型为长方体，水上方的气体区域尺寸：0.8 m×0.8 m×0.4 m；水区域尺寸：0.8 m×0.8 m×7.6 m。

3.2 计算边界条件与材料属性

弹体初始速度为500 m/s、750 m/s、1000 m/s、1250 m/s、1500 m/s；弹体为钨合金，密度为17250 kg/m³；气体介质为空气；液体介质为清水。

3.3 流体域网格划分

气-液-固耦合下的接触问题是建立弹体流-固耦合数值仿真的关键，弹体入水涉及到弹体大位移运动，采用流体网格重构的方法不仅会耗费大量的计算资源，而且可能由于网格重构的失效而导致无法顺利完成数值计算。因而，本文采用基于规整的笛卡尔网格的浸入边界法，其核心思想是将物体形状和位置显式表示为一个函数，该函数定义了等值曲面，曲面上的点为物体的表面，可得到物体形状随着时间的变化，该方法在模拟流体与固体相互作用问题上非常适用，对弹体入水过程及空泡形态演化过程可很好地呈现，同时流场的计算是基于直角坐标网格，极大地提高仿真计算效率，无需重划网格，只需更新耦合界面。在本模型中流体域的网格为5000万六面体网格。

3.4 流-固耦合分析结果

图3为弹体在空气中高速运动即将入水瞬间的液体表面状态。高速弹体运动在行进方向压缩空气，进而对液体表面产生作用，形成了液面局部凹陷。

弹体入水时，因弹体头部首先跟水接触，弹头表面附近水受到弹体的挤压，形成一个压缩区，压缩区中的压力以弹性波的形式向周围扩散，压缩区的大小及形状与弹头端部的形状以及弹体的速度密切相关。由于弹头受到水的阻力作用，压缩波不仅会在水中传播，同时也会在弹体结构中传播。在液体形成压缩区的瞬间过程中，弹体会承受很大的液体冲击力，可能会造成弹体局部承受很高的应力，从而可能造成弹体局部结构的损坏。

由于在弹体入水的瞬间，气体与初始形成的微小空泡并没有完全断开，即存在气体通道，水表面的空气会不断被卷入到水中，引起空泡的不断膨胀和形状变化，此阶段称为开空泡阶段。由图4可知，空泡分布在弹体周围的薄层内，随着弹体一起高速运动，由于弹体入水速度快，弹体表面空泡来不及闭合。在弹体冲击和压力作用下，初始速度为1500 m/s的弹体入水后，在其顶端产生强大的压力波对液体进行压缩，弹头局部水域的密度增大至1254 kg/m³，该部水域被高度压缩。

图5为弹体入水前和入水后压力等值线图，图5(a)为入水前流场的压力分布图，最高压力集中在弹体头部，压力值为0.814 MPa；图5(b)为弹体入水后在0.00096 s时，弹体在水中的压力等值线分布图。

由图5(b)可知，高速运动的弹体入水后，弹体对水具有强烈的压迫，弹体头部附近液体的压强在波阵面上发生突跃变化，压强等值线十分密集，变化梯度很大。可知在弹体顶部的运动方向上，流体压强最大为2443 MPa，从顶部向弹身方向周边流体的压力逐渐降低。

图6为弹体入水后，不同时刻的水面变化过程。在弹体初始入水瞬间，大气与空泡相连，即开空泡阶段，随着时间变化，弹体入水深度不断增加，如图6(b)和图6(c)所示，高速度运动的弹体带动附着在弹体表面的流体一起高速运动，从而造成弹体周围薄层内的流体压力急剧减小，空泡不断膨胀，体积不断变大，直至空泡发生溃灭。当空泡溃灭后，弹体周围薄层的流体压力很小，在弹体周围远区水压力的作用下，水向薄层内急剧汇聚，造成水液面的喷溅，如图6(c)和图6(d)所示，

图7为弹体随时间变化的位移曲线，随着时间的增加，弹体在阻力作用下，速度虽有降低，但与初始速度相比，变化值相对较小，因此在入水初始阶段还可将位移-时间曲线近似看成直线。图8为弹体运动速度随时间变化曲线。可以看出，弹体入水前，可视为匀速运动，入水后受到阻力作用速度逐渐降低。弹体初始速度越高，入水后速度下降越快。通过观察可看出，弹体入水后的速度几乎为线性变化，说明弹体受力稳定，速度逐渐降低。

图9为弹体运动加速度与时间之间的关系曲线，当弹体初始速度为500 m/s 时，弹体的加速度峰值为−12120 m/s²；当弹体初始速度为1500 m/s时，加速度峰值为−35980 m/s²，入水载荷最大。弹体入水后，加速度迅速增大，当增大到某一值后加速度数值的变化不再剧烈，入水载荷相对平稳。弹体初始速度越大，弹体入水后的加速度和加速度振幅越大，初速度越低，弹体的加速度和加速度振幅越小。

图10为弹体表面压强与弹体的初始速度关系曲线。可知，弹体表面的压强峰值会随着弹体初始速度的增加而线性增加。当弹体初始速度为500 m/s时，弹体压强峰值为720 MPa，弹体初始速度为1500 m/s时，弹体表面压强峰值为2440 MPa。因此，超高的初始速度对弹体的强度也提出了更高要求。

4 结　语

1） 建立了高马赫数弹体垂直入水流固耦合的数值仿真方法，该方法可得出弹体入水过程中运行轨迹、承受压力的变化情况，同时流体的速度、压力、密度变化以及超空泡的演化过程；

2） 高速弹体入水加速度峰值和弹体表面压强峰值会随着弹体初始速度的增加而增加，而弹体的运动速度反而逐渐下降；

3）弹体高速入水瞬间，液体表面压力快速增加，流体密度随之增大；由于入水速度快，弹体表面空泡不易闭合；

4）弹体以初速度1500 m/s入水时，在其顶端产生强大的压力波对液体进行压缩，形成空泡，水的密度增大至1254 kg/m³，弹体压强峰值为2440 MPa。