水下高速运动航行体表面附近的水域会因压力低于水的饱和蒸汽压而发生空化现象，进而在航行体表面形成空泡，当速度较高时会形成完全包裹航行体表面的超空泡^[1].超空泡的形成能够明显降低航行体水下运动时的阻力，可大幅提高航行体的水下运动速度^[2]，但是由于空泡流动本身的复杂性使得超空泡技术在实际应用中面临诸多挑战，以超空泡航行体的稳定性和控制问题尤为突出^[3-4].

高速超空泡射弹在水下运动时，任何微小的扰动(如射弹在发射时受到的扰动或者水的不稳定流动等)都可能使射弹尾部在运动过程中与空泡壁面发生连续往复撞击，这一现象称为尾拍运动^[5].近年来，超空泡射弹的尾拍运动引起了国内外相关学者的关注和研究，孟庆昌等^[5]对超空泡射弹的尾拍运动进行了数学建模及分析；Ruzzene等^[6-7]通过初步研究建立了高速运动的超空泡射弹在空泡内部往复振动的动力学模型，研究了射弹在尾拍撞击作用力下的结构响应问题；Savchenko等^[8-12]通过在开放式水洞中进行的模型实验研究了圆柱体模型在超空泡内部滑行时的滑行力；金大桥等^[13]设计并实施了通气超空泡水下射弹试验，对射弹水下产生的超空泡形态进行了试验研究；曹伟等^[14]通过射弹试验对超空泡的发展规律和几何特征进行了深入研究；何乾坤等^[15]对超空泡射弹尾拍运动结构动力学响应和空泡摆动对射弹的尾拍撞击影响进行了研究；冯光等^[16]应用细长体理论计算得到水下超空泡航行体的流体动力，进而得到超空泡状态下的水下航行体弹道；而初始扰动对超空泡射弹尾拍运动和弹道特性影响的相关研究较少.

随着计算机技术的飞速发展，基于N-S方程的计算流体动力学(CFD)方法已成为揭示一些复杂流动现象的有效途径.本文通过对流体动力计算软件CFX进行二次开发，运用程序自带的CEL语言将射弹在垂直平面内运动的刚体运动方程离散化并添加到程序的计算求解过程中，结合动网格技术，实现了多相流动非定常雷诺平均N-S方程(URANS)与刚体运动学方程的耦合求解，建立了水下超空泡射弹纵向平面的尾拍运动动力学计算模型并分析了初始扰动对空泡形态及射弹尾拍运动和弹道特性的影响.本文建立的计算模型考虑了射弹与空泡表面复杂的作用力、射弹的自由运动、空泡的时变特征及湍流对射弹的流体动力影响.

1 数值计算模型 1.1 基本控制方程

本文研究水下运动的超空泡射弹，涉及的是自然空化流动问题，流体域的两相介质为水蒸汽和水.连续性方程为

$\frac{\partial {{\rho }_{m}}}{\partial t}+({{\rho }_{m}}{{v}_{j}}){{x}_{j}}=0.$

动量方程为

$\begin{align} & \frac{\partial \left( {{\rho }_{m}}{{v}_{i}} \right)}{\partial t}+\frac{\partial \left( {{\rho }_{m}}{{v}_{i}}{{v}_{j}} \right)}{\partial {{x}_{j}}}=-\frac{\partial p}{\partial {{x}_{i}}}+ \\ & \frac{\partial }{\partial {{x}_{j}}}\left[ ({{\mu }_{m}}+{{\mu }_{t}})\left( \frac{\partial {{v}_{i}}}{\partial {{x}_{j}}}+\frac{\partial {{v}_{j}}}{\partial {{x}_{i}}} \right) \right]. \\ \end{align}$

式中：${{\rho }_{m}}=\sum\limits_{k=1}^{n}{{{\alpha }_{k}}{{\rho }_{k}}}$、${{\mu }_{m}}=\sum\limits_{k=1}^{n}{{{\alpha }_{k}}{{\mu }_{k}}}$分别为按体积分数确定的混合物的密度和动力黏性系数；v_j为速度.

1.2 湍流模型

湍流是一种强非线性流体状态，是空化流动的主要特征之一.空泡发展及闭合等均伴随着湍流现象.本文采用的湍流模型是基于Baseline k-ω模型的SST(shear stress transport)湍流模型.该湍流模型能够较好地解决湍流剪切应力的传输问题，对于逆压梯度下的流动分离问题有更精确的预测，基本方程为^[17]

$\begin{align} & \frac{\partial \left( \rho k \right)}{\partial t}+\frac{\partial \left( \rho {{v}_{j}}k \right)}{\partial {{x}_{j}}}= \\ & \frac{\partial }{\partial {{x}_{j}}}\left[ \frac{\partial k}{\partial {{x}_{j}}}\left( \mu +\frac{{{\mu }_{t}}}{{{\sigma }_{k3}}} \right) \right]+{{P}_{k}}-\beta \prime \rho k\omega +{{P}_{kb}}, \\ & \frac{\partial \left( \rho \omega \right)}{\partial t}+\frac{\partial (\rho {{v}_{j}}\omega )}{\partial {{x}_{j}}}= \\ & \frac{\partial }{\partial {{x}_{j}}}\left[ \frac{\partial \omega }{\partial {{x}_{j}}}\left( \mu +\frac{{{\mu }_{t}}}{{{\sigma }_{\omega 3}}} \right) \right]+{{\alpha }_{3}}\frac{\omega }{k}{{P}_{k}}-{{\beta }_{3}}\rho \omega +{{P}_{\omega b}}+ \\ & 2(1-{{F}_{1}})\rho \frac{1}{{{\sigma }_{\omega 2}}\omega }\frac{\partial k}{\partial {{x}_{j}}}\frac{\partial \omega }{\partial {{x}_{j}}}. \\ \end{align}$

SST湍流模型的涡黏度限制方程为

${{\nu }_{t}}=\frac{{{a}_{1}}k}{max\left( {{a}_{1}}\omega ,S{{F}_{2}} \right)}.$

式中：μ为流体黏性系数；μ_t为湍流黏度；ν_t=μ_t/ρ；P_k为湍动能中由黏性力产生的部分；P_kb为湍动能中由浮力产生的部分；F₁、F₂分别为混合函数；常数β′=0.09、α₁=5/9；S为应变率的不变测度.

1.3 空化模型

用来描述流场空化程度强弱的特征参数为空化数，为

$\sigma =\frac{{{p}_{\infty }}-{{p}_{c}}}{12\rho {{U}^{2}}_{\infty }}.$

式中：p_∞为流场无穷远处压力；p_c为当地温度下水的饱和蒸汽压；ρU²_∞/2为流场的参考动压力.由空化数的定义可以看出，空化数越小，流场的空化程度越强.

不计热传输和非平衡相变效应，空化流动中水蒸汽的体积分数输运方程为

$\frac{\partial ~({{\rho }_{v}}{{\alpha }_{v}})}{\partial t}+\frac{\partial ({{\rho }_{v}}{{\alpha }_{v}}{{v}_{i}})}{\partial {{x}_{i}}}={{{\dot{m}}}^{+}}-{{{\dot{m}}}^{-}}.$

式中：${{{\dot{m}}}^{+}}$ 、${{{\dot{m}}}^{-}}$分别为蒸发项和凝结项.本文采用Rayleigh-Plesset气泡动力学模型来描述液体中气泡的生成和溃灭，蒸发项和凝结项分别为

$\begin{align} & {{{\dot{m}}}^{+}}={{F}_{e}}\frac{3{{r}_{nuc}}\left( 1-{{\alpha }_{v}} \right){{\rho }_{v}}}{{{R}_{B}}}\sqrt{\frac{2}{3}\frac{|{{p}_{v}}-p|}{{{\rho }_{l}}}}, \\ & {{{\dot{m}}}^{-}}={{F}_{c}}\frac{3{{\alpha }_{v}}{{\rho }_{v}}}{{{R}_{B}}}\sqrt{\frac{2}{3}\frac{|{{p}_{v}}-p|}{{{\rho }_{l}}}}. \\ \end{align}$

式中：p_v、p分别为气泡内压力和环境压力；α_v为蒸汽的体积分数；r_nuc=5×10^-4为成核点的体积分数；R_B为气泡半径；F_e=50，F_c=0.01均为经验常数.

2 模型尺寸及网格划分

本文计算采用的射弹模型为具有大长细比的锥柱结合体，如图 1所示.

射弹结构包括3段，分别为头部柱段、肩部锥段和后体柱段.头部柱段直径4 mm，长度3 mm；肩部锥段的半锥角为4°；后体柱段直径D=10 mm；弹体总长L= 170 mm；质量为0.09 kg；质心距离弹体头部中心点95.8 mm；射弹绕通过质心且垂直于弹体轴线的坐标轴的转动惯量为1.72×10^-4 kg·m².

计算流域为圆柱体，入口距弹体头部长度取2L=340 mm，出口距弹体尾部取20L=3 400 mm，直径取60D=600 mm，x轴位于弹体轴线上，计算中考虑重力影响.整个计算域均采用六面体结构化网格离散，单元数量约为163万，弹体附近网格如图 2所示.计算中采用CFX中CEL语言将刚体运动学方程离散后并嵌入程序包中与U-RNAS方程耦合求解，瞬态计算开始前，首先通过调整流场参数得到覆盖整个弹体的超空泡流动稳态数值解作为后续计算的初始值，以利于计算收敛.射弹初始位置位于水面下5 m深度，初始速度v₀ = 150 m/s为水平方向，初始自然空化数σ_v=0.013 1.

针对图 1所示的射弹模型，建立了粗糙(101万)、中等精细(163万)、精细(210万)3种不同精度的网格，3种网格数量下射弹头部附近压力系数对比如图 3所示.

从图 3中可以看出，3种网格精度下射弹头部附近压力系数变化差异较小，其中精细网格(210万)和中等精细网格(163万)的计算结果基本一致，而粗糙网格(101万)在射弹头部与精细网格和中等精细网格的计算结果有一定差异，综合考虑，为了节省计算资源，并且达到所需精度，本文选用中等精细网格(163万)进行计算.

为验证本文所建立数值计算模型的有效性，进行了v=150 m/s，水深5 m时的仿真计算，空泡形态与Lognovinch独立膨胀原理所得的空泡形态对比如图 4所示.

由图 4可以看出，本文数值模型所得空泡形态与独立膨胀原理所得空泡形态基本一致，在空泡尾部由于二者的处理方法不同而略有差异，这说明了本文所建立数值计算模型的有效性.

3 结果与分析 3.1 空泡形态分析

初始空化数σ_v= 0.013 1，初始扰动角速度为θ·₀=6 rad/s的射弹尾拍运动在一个周期的空泡形态变化如图 5所示.

由图 5可以看出，射弹的尾拍运动引起空泡形态的改变，破坏空泡原有的光滑对称壁面.弹体尾部与空泡壁面碰撞的位置，空泡形态发生显著改变，碰撞点处发生动量交换，空泡壁面凸出，造成空泡形态不对称.

以射弹尾部截面为0点，图 5(c) 中发生撞击时弹体尾部沿轴向的压力分布如图 6所示.

由图 6可以看出，当射弹尾部撞击空泡壁面时，射弹尾部表面的压力由原来空泡内较低且均匀的空化压力改变为图 6所示的在弹体尾部0<L/D<2.5范围内产生一个压力峰值，此压力峰值的产生主要由于弹体尾部与空泡壁面的撞击使得弹体尾部的压力与空泡外撞击流域的压力相近，压力峰值位于弹体尾部与空泡壁面的撞击点附近.

3.2 尾拍运动特性分析

初始水平速度v₀ = 150 m/s，初始空化数σ_v= 0.013 1，初始扰动角速度分别为为θ·₀= 2、4、6、8 rad/s的射弹尾拍运动的俯仰角变化历程如图 7所示.

从图 7中可以看出，不同初始扰动下射弹俯仰角的变化不同，具体表现为俯仰角的幅值和变化周期均不相同.初始扰动角速度越大，俯仰角变化的最大幅值越大，变化频率较快.在t=0.055 s后，俯仰角不再呈现周期性变化，因为此时空泡开始溃灭于弹体尾部并逐渐减小至消失，如图 8所示.

不同初始扰动角速度下射弹尾拍运动的俯仰角速度和俯仰角加速度变化历程分别如图 9,10所示.

尾拍振动的周期越来越短，频率越来越高，随着后期射弹速度的衰减，空化数变化，超空泡长度变短并逐渐溃灭于弹体上，射弹尾拍振动时刺透空泡壁面浸入流体的深度增加，造成尾拍反力的增加，进而影响尾拍运动俯仰角加速度的增加；当空泡开始溃灭于弹体上并逐渐消失时，射弹由尾拍运动变为双空泡运动(如图 8所示)，尾拍反力急剧减小并消失，故射弹的俯仰角加速度的变化也呈现先增加后急剧减小的趋势；且初始扰动角速度越大，射弹尾拍运动的俯仰角速度和角加速度的幅值也越大，变化频率越高.

3.3 弹道特性分析

图 11给出了射弹在不同初始扰动下以相同初速度自由运动的水平速度变化情况.初始扰动越大，射弹的速度衰减越快，但整体而言射弹水平速度变化的差别并不明显.

不同初始扰动下射弹的阻力系数变化如图 12所示.阻力系数的变化主要包括两个比较明显的阶段，当t<0.06 s，阻力系数主要呈现为周期振动特性，这是由于射弹尾拍运动时弹体尾部刺透空泡壁面并浸入水中增加了流体阻力；当t>0.06 s，阻力系数开始急剧增加，主要是由于空泡长度变短接触到弹体尾部并开始逐渐溃灭于弹体上，因此阻力系数急剧增加.不同初始扰动下尾拍运动时弹体尾部浸入水中的深度不同，所以引起的阻力变化幅值也不同，初始扰动越大，阻力系数变化的振幅越大，频率越高.

图 13,14分别给出了不同初始扰动下射弹自由衰减运动的水平位移和垂直位移的变化历程.从图中可以看出，不同初始扰动对射弹的水平位移影响较小，而对射弹的垂直位移影响较大.

初始扰动分别为2 rad/s和8 rad/s的两种情况下，在相同运动时间内射弹的垂直位移相差约1倍；随着初始扰动的增加，射弹的垂直位移增大；射弹垂直位移的增加是重力和扰动的共同作用.

4 结 论

1) 不同初始扰动下射弹尾拍运动的俯仰角、俯仰角速度、俯仰角加速度、阻力系数都表现出周期性变化的特点，且初始扰动越大，变化的幅值越大，频率越高，尾拍振动的频率随着时间发展越来越高.

2) 不同初始扰动对射弹速度衰减和水平位移变化的影响较小，但对射弹垂直位移的影响较大.

3) 超空泡射弹尾拍运动时弹体尾部与空泡壁面发生碰撞并浸入流体中，破坏空泡形态的对称性，射弹尾部碰撞点附近产生较大的压力峰值.