0 引言

导弹等水下航行体从水下一定深度的发射筒内，利用燃气向水面垂直发射。航行体出筒后发射筒内高温高压的燃气会涌出发射筒，在海水中形成脉动压力场，并与筒盖结构以及航行体尾部交互耦合作用。这种脉动的筒口压力场将会对发射筒盖等结构强度形成不利影响^{[1, 2]}。

针对筒口气泡脉动形态和压力预示，目前常采用理论计算和数值CFD模拟2种方法开展研究。在理论计算方面，王诚等^[3]采用一维非定常气流场模型和轴对称理想水流场模型，对水下发射导弹的尾部非定常燃气泡内外流场进行耦合数值求解，计算了燃气泡发展对导弹所受水动力的影响。李杰等^[4]采用Rayleigh-Plesset方法，建立压球状尾空泡模型，研究了空泡对弹道的影响。程少华等^[5]结合空泡独立膨胀原理和Rayleigh-Plesset方程，建立了导弹尾空泡非定长发展数学模型，对弹尾空泡和压力进行预示。

在数值CFD模拟方面，王汉平等^[6-7]基于VOF模型和动网格技术对潜射导弹出筒筒口压力场特性进行了二维和三维仿真分析。邱海强、王亚东等^[8-9]利用Mixture多项流模型和动网格技术，对导弹弹射出筒筒口气泡的发展情况进行研究，并重点关注了筒口气泡对发射艇表面压力脉动的影响。

本文采取刚体Lagrange结构网格和Euler流场网格（L/E）耦合的方法^[10-11]开展筒口压力场三维数值仿真计算。计算模型考虑航行体出筒过程对筒口流场的扰动和横向来流对筒口压力场的影响。

1 控制方程及数值方法 1.1 控制方程

流体运动采用Euler方法模拟，建立流体欧拉域。控制方程包括质量、动量和能量守恒。

质量守恒：

$ \frac{\partial }{{\partial t}}\iiint\limits_{vol} {\rho {\rm d}V} = - \iint\limits_{surf} {\rho (u \cdot n){\rm d}A} \text{；} $

动量守恒：

$ \begin{array}{l} \displaystyle\int \!\!\!{\int\limits_{vol}\!\!\! {\int {\rho {u_i}{\rm d}V} } } + \int \!\!\! {\int\limits_{surf} {\rho {u_i}(u \cdot n){\rm d}A = - } } \\ \quad \quad \quad \quad \quad \displaystyle\int\!\!\! {\int\limits_{surf} {p{n_i}{\rm d}A + \int \!\!\! {\int\limits_{surf} {{s_{ij}}{n_i}{\rm d}A} } } } \text{；} \end{array} $

能量守恒：

$ \frac{\partial }{{\partial t}}\iiint\limits_{vol} {\rho e{\rm d}V} + \iint\limits_{surf} {\rho e(u \cdot n)} {\rm d}A = - \iint\limits_{surf} {{u_i}p{n_i}} {\rm d}A\text{。} $

式中：ρ为材料密度；u_i为速度分量；p为压力；s_ij为粘性偏应力张量；e为比内能。

1.2 状态方程

状态方程定义各种不同状态下的压力和密度以及比内能之间的函数关系。对于海水，采用多项式状态方程。

压缩状态（μ > 0）：

$ p \!=\! {a_1}\mu \!+\! {a_2}{\mu ^2} \!+\! {a_3}{\mu ^3} \!+\! ({b_0} + {b_1}\mu \!+\! {b_2}{\mu ^2} \!+\! {b_3}{\mu ^3}){\rho _0}e\text{，} $

(1)

拉伸状态（μ < 0）：

$ p = {a_1}\mu + ({b_0} + {b_1}\mu ){\rho _0}e\text{。} $

(2)

其中$\mu = \frac{\rho }{{{\rho _0}}} - 1$。取海水的初始密度ρ₀=1 024 kg/m³，体积模量a₁=2.2 × 10⁹ Pa，对于高阶项大多数情况取b₀=b₁=b₂=b₃=0，比内能e取83 950 J/kg^[11]。

对于燃气采用Gamma状态方程：

$ p = (\gamma - 1)\rho e\text{。} $

(3)

式中：e为燃气比内能；ρ为燃气密度；比热比γ=1.2。

1.3 耦合算法

本计算中包括固体域和流体域2个计算域。固体域主要模拟航行体，采用Lagrange方法模拟，Lagrange网格设置为刚体模拟航行体运动，但不发生变形。流体域包括海水和燃气所填充的Euler网格区域，采用Euler方法模拟。2个计算域相互重叠。为了使固体域和流体域间产生相互作用，在Lagrange网格结构外部定义一封闭耦合面，用于传递2种解算器之间的力。此耦合面对Euler网格来说是一种流动边界，材料在Euler网格中流动产生的压力通过耦合面传到Lagrange结构上，使结构发生运动或变形；Lagrange结构运动或变形也反过来影响Euler材料的流动和压力值。

2 计算模型 2.1 航行体模型

航行体沿轴线建立一半对称模型，头部采用双曲线头型，如图 1所示。航行体设置为刚体，除了垂直向上的运动外，其余5个自由度固定。航行体筒内和出筒运动的轨迹已有较为精确的计算数据，直接加载到航行体上。航行体圆柱段直径为D_m，总长为7D_m。

2.2 流场模型

由于流体域相对于过轴线的海水横向流流动方向对称，为减小计算量，建立对称流体域计算模型。取发射筒口中心为坐标原点，海水流场区域为半圆柱，半径5.6D_m，深22D_m。海水区域下为发射筒流场，按照航行体尺寸建立。

计算区域及边界条件如图 2所示。流场区域上下为压力边界，左右为海水速度入口和出口边界，最底部为燃气入口边界，其他为固体边界。

计算流场区域网格如图 3所示。对航行体运动区域和筒口压力场作用区域网格进行加密。初始时刻模拟航行体的Lagrange计算网格与发射筒流场区域的Euler网格重叠（见图 4）。Lagrange固体网格表面为流体-固体耦合面，计算中耦合面内重叠区域Euler网格材料自动为空。

3 计算结果分析 3.1 含横向流筒口压力计算

航行体离开发射筒前，发射筒内燃气封闭于发射筒与航行体组成的空腔内，压力高于周围海水压力。航行体尾部离开发射筒后，由于受压力及惯性等作用，筒内燃气向外扩张（见图 5），形成第1个压力峰。随着燃气在海水中的膨胀，燃气泡内的压力不断下降。由于惯性作用，在到达压力与外界压力平衡时，燃气泡将继续膨胀，气泡内的压力将小于气泡外界压力（见图 6），直到燃气泡压力达到最小极值停止膨胀。随后，由于燃气泡内外压力不平衡，燃气泡开始收缩。

由于惯性作用，燃气泡压力将超越与海水环境压力等值点继续升高，直到达到第2个压力峰值，在此过程中部分燃气随航行体尾部与主脉动气泡分离（见图 7）。此后燃气泡持续脉动，一部分能量在海水环境中耗散，一部分能量随燃气泡上升过程带走。此过程与水下爆炸压力场脉动特征相近。由计算结果可以看出，由于计算中考虑了横向海流，筒口气泡在脉动过程中出现了向X正向的飘移。

3.2 无横向流筒口压力场计算

将海水速度入口和海水速度出口边界条件改为压力边界条件，并将海水初始速度赋值为0，开展无横向流筒口压力场计算。在此工况下，筒口气泡在脉动过程中没有出现向X正向的飘移（见图 8）。

3.3 试验结果对比

图 9给出了试验与含横向流计算的筒口压力曲线。图 10为试验与无横向流计算的筒口压力曲线。曲线X时间轴从航行体尾部出筒时刻开始。Y轴为测点压力与航行体出筒时刻筒内燃气压力的比值。由2组对比曲线可知，计算筒口压力曲线和试验压力曲线均呈现剧烈的压力脉动，随着时间推移压力脉动快速衰减。试验与计算的2个主要压力波峰曲线形状基本一致。

取2个主要压力波峰极值作为特征点，对比计算与试验结果偏差（见表 1）。由表 1可知，计算结果相对于试验值均略为偏高。其中含横向流工况2个特征点的计算偏差为3.81%和14.3%，无横向流工况2个特征点的计算偏差为4.59%和30.1%。含横向流计算压力波峰值的偏差明显小于无横向流计算，特别是第2个波峰计算偏差减小至0.5倍以下。说明计算中考虑横向流作用可以有效减少筒口压力场计算偏差。

4 结语

本文采用Lagrange结构网格和Euler流场网格耦合的方法，对潜射航行器出筒后筒口压力场进行三维数值仿真分析。由于计算模型考虑了海水的横向来流，筒口气泡在脉动中出现了随来流偏转的特征。试验与计算的2个主要压力波峰曲线形状基本一致。2个主要压力峰值特征点的计算偏差为3.81%和14.3%，可以满足工程研制需要。含横向流工况和无横向流工况计算偏差对比表明，计算中考虑横向流作用可以有效减少筒口压力场计算偏差。本文的建模方法和仿真结果对潜射航行器筒盖系统的受载分析以及筒盖系统结构设计，具有工程参考价值。