水下航行体垂直发射环境流场与弹道耦合数值模拟研究

2019-10-15

节能技术 2019年4期

(北京宇航系统工程研究所，北京 100076)

0 引言

潜射航行体的发射是一个涉及多因素干扰、多相介质参与的复杂物理过程，目前为应对各种复杂环境需求，航行体朝着变深度发射、简化结构、提高发射可靠性等趋势发展[1]。其中发射可靠性一直是其重点研究方向之一，由于航行体需经历穿越水层、空中点火等对其结构和强度有极高要求的过程，且以海洋为发射环境使得诸多不可预知因素成为抑制航行体弹道稳定性的技术瓶颈，又考虑水下发射点火升空时间极短，需要完成多个复杂的控制动作及过程，因此如何能在可预知和不可预知的多因素作用下开展航行体发射可靠性、提高其弹道稳定性的控制技术研究显得极为必要。

水下平台垂直冷发射[2]的航行体水下运动过程可分为三个阶段：出筒段、水中航行段和出水段[3]。发射时平台须保持一定的航行速度来维持其操纵性，平台航速的存在破坏了流场的对称性，使得航行体所受流体动力呈现显著的三维特征，进而出现较大的俯仰姿态[4-5]。变深度发射有利于航行体发射的灵活性，将是水下发射技术的发展趋势之一，因而探讨有平台航速条件下航行体的流体动力特征，研究平台航速和发射水深对航行体流体动力和运动速度的影响具有重要的工程价值。

随着多相流模型、动网格技术、湍流模型等CFD技术的不断发展，目前数值模拟已成为研究水下发射技术的常用手段。多相流模型根据界面处理的不同大致可分为两类，一类是自由界面追踪模型，如VOF[6]、MAC[7]模型，可以处理自由界面变化较剧烈但界面基本清晰的流动问题；另一类是混合模型，如Mixture模型等，可模拟各相之间强烈掺混、界面基本不清晰的流动问题。Dyment[8]，刘志勇等[9]分别采用VOF方法和改进的MAC方法模拟了尾空泡发展过程；李杰等[10]采用VOF方法对细长回转体出水过程进行了数值模拟；张红军等[11]采用Mixture模型对航行体出筒过程进行了三维非定常数值模拟研究；杨晓光[12]等采用VOF模型和动网格技术对航行体水下运动及出水过程进行了三维流场仿真。Ma 等采用MIX模型对出口航行体排气气膜孔参数进行了三维仿真[13]，后续又进行了考虑平台速度和偏航角不确定性的水下航行器弹道和姿态的数值研究[14]。

本文研究将采用Mixture模型、标准k-ε湍流模型以及动网格技术求解水下航行体非定常发射过程，通过对流场计算获得航行体所受水动力，而后将水动力作为力学参数输入以求解航行体的运动过程，从而实现航行体水下发射多相流场和弹道耦合求解。通过分析航行体表面压力以及流体动力的三维特性，以探讨发射平台航速和发射水深对航行体所受流体动力和运动速度的影响。

1 坐标系及无量纲化定义

图1所示为发射系和航行体体系坐标系示意图，本文将在发射系中讨论航行体的运动速度特征，在航行体体系中讨论流体动力特征。坐标系的定义如下:

发射系：以航行体初始时刻质心位置为原点O，以发射平台航速方向为X轴(横向)，以垂直出筒方向为Y轴(纵向)，Z轴方向由右手法则确定。航行体体系：以航行体实时的质心位置为原点O1，法向方向为X1轴，体轴方向为Y1轴，Z1轴方向由右手法则确定。

本文中压力p和时间t分别以发射水深静压P和出筒时刻T进行无量纲化；各速度分量均以水深H、平台航速U下计算获得的最大值为参考进行无量纲化；对力和力矩的无量纲化定义为

Cd=F/(0.5ρV2S)

(1)

Cm=M/(0.5ρV2SL)

(2)

式中ρ——水的密度;

V——选取水深H、平台航速U下的航行体最大纵向速度;

S——航行体特征横截面积;

L——航行体轴向长度。

2 数值计算模型

本文采用动网格技术将多相流流场和航行体弹道进行耦合计算。首先利用多相流模型计算当前时刻瞬态流场，通过对航行体表面压力积分，获得当前时刻航行体所受流体力和力矩，在此基础上计算航行体的加速度、速度和位移等；然后利用动网格技术获得下一时刻的航行体边界和网格，进而求解下一时刻的瞬态流场，接着计算下一时刻航行体的加速度、速度和位移等。上述过程不断循环，直至航行体出水为止。

2.1 流体控制方程

本文采用商用数值计算软件Ansys Fluent,采用Mixture模型模拟气液两相流动，忽略相间的滑移速度，计算中考虑水体为不可压缩相，气体为可压缩相。采用有限体积法离散控制方程，结合Simple算法来求解下列方程。

连续性方程

(3)

动量方程

(4)

能量方程

(5)

次相体积分数方程

(6)

状态方程

pg=ρgRT

(7)

ρw=const

(8)

αk——第k相的体积分数；

keff——有效传热系数；

SE——包含了所有的体积热源。

上述方程采用标准k-ε湍流模型封闭。

2.2 航行体动力学方程

本文研究中将航行体当作刚体，尾出筒后开始计算其俯仰运动(忽略出筒过程的俯仰运动)。通过下式计算航行体质心的合力和绕质心的合力矩

(9)

(10)

2.3 计算域、边界条件、网格划分及更新方法

本文建立如图2所示的三维对称模型，计算域包括发射筒、水域和大气域。发射筒底燃气入口为压力入口边界，航行体尾出筒后改为壁面边界；水域和空气域外围均为压力出口边界。为计算尾出筒后航行体的俯仰姿态，在航行体运动区域内划分非结构网格，采用网格当地重构更新动网格，为保证航行体壁面附近网格质量，对航行体表面附近进行了网格加密，该区域内的网格质量在航行体运动过程不发生变化。数值计算网格采用Gambit生成，初始计算网格数量为600万。

3 航行体表面压力和流体动力三维特征分析

为获得航行体表面压力的三维特征，计算时在航行体迎背水面布置了多个压力监测点，测点1为航行体驻点，测点2～4分布于头锥段，测点5～8分布在柱段，测点9位于尾段。各测点距离航行体实际顶点的长度如表1所示(以航行体轴向长度L为参考)。

图3对比了测点1和迎水面测点3、5、7、9处的计算与试验压力值。从曲线分布可知，计算值与试验值吻合较好，验证了计算模型的合理性。航行体出筒及水中航行过程中，驻点压力随着航行体速度增加而升高、下降而递减；头锥段压力因为绕流作用和环境压力的下降而呈现出下降趋势；柱段脱离发射筒后，压力随着环境水压的降低而降低，柱段靠近尾段部分(如测点7)在尾出筒时受燃气回射作用形成了明显的压力脉冲，且尾出筒后受尾空泡生成演化的影响，呈现出明显的压力振荡特性。

表1压力监控点布置

测点序号123456789测点位置00.034 L0.099 L0.198 L0.297 L0.481 L0.741 L0.916 L0.962 L

如图4所示，发射平台航速导致航行体迎背水面流场不对称，航行体以一定的俯仰姿态出水。图5给出了测点2～9迎背水面压差时变曲线，从中可以看出，迎水面压力高于背水面压力，压差主要集中在航行体头锥段。同时头锥段压差在出水前随着航行体速度增大而升高、下降而递降；柱段压差主要出现在水中航行段，两级上明显小于头锥段最大压差，且受燃气后效影响，靠近尾段部分(如测点8)压差在出筒时刻附近形成了脉冲变化；尾段压差则主要受尾空泡演化影响，除在尾空泡拉断后一段时间形成了一定量值的压差外，大部分时间几乎为零。

迎背水面压差的存在导致航行体承受较大的法向力和俯仰力矩。从图7给出了航行体所受流体法向力和俯仰力矩的分布可知，整个发射过程航行体所受法向力、俯仰力矩主要由头锥段和柱段组成，尾段贡献几乎为零：头锥段法向力量值在出筒过程中不断增大，并成为法向合力的主要组成部分，在尾出筒后逐渐减小；随着柱段进入水中，柱段法向力量值逐渐增大，水中航行段继续保持增大趋势，柱段出水过程逐渐减小至零；出水前，头锥段力矩在出筒前后呈现先增大后减小的变化趋势，并且决定了此阶段俯仰合力矩的变化趋势；柱段力矩在出水前量值较小，随着柱段出水逐渐减小至一个较大的负值后逐渐增大至零，本文中称之为俯仰力矩的“反号”变化，并决定了此阶段俯仰合力矩的变化趋势，这是由于法向合力作用点移动到质心之后所致。

4 发射平台航速影响分析

本文在相同发射水深H条件下，研究了发射平台航速(0.4U、0.7U和1.0U)对航行体流体动力和运动速度的影响。图7所示为，不同航速下航行体所受流体动力的时变曲线。由图可知，流体轴向力与尾段压力的变化规律基本一致。如图3中所示测点9压力在尾出筒后呈现周期性振荡变化；随着发射平台航速增加，流体法向力和俯仰力矩量值增大，但轴向力几乎不变。

图8比较了不同平台航速下航行体运动速度的时变曲线。航行体尾出筒后，横向速度逐渐减小，且出现了反号变化(即沿相反的方向运动)，尾出水后不再变化；纵向速度逐渐减小；俯仰角速度逐渐增大，受俯仰力矩“反号”变化的影响，头出水后逐渐减小，尾出水后不再变化；随着发射平台航速增大，横向速度下降加快，俯仰角速度增大，但纵向速度几乎不受影响。图9所示为尾出水后横向速度和俯仰角速度与初始运动时刻的绝对变化量，该变化量与平台航速基本呈线性关系。

5 发射水深影响分析

下文在相同发射平台航速U条件下，探讨了不同发射水深(0.8H、1.0H和1.2H)对航行体流体动力和运动速度的影响。图10所示为不同发射水深下航行体所受流体动力的时变曲线。由图可知不同发射水深条件下，航行体所受外力变化趋势基本一致；随发射水深增加，航行体水中段运动时间增长，轴向力经历周期性振荡次数增多，导致出水引起的法向力和俯仰力矩变化发生时间延迟。

由图11所示不同发射水深下航行体的运动速度时变曲线可知，不同发射水深条件下，航行体运动速度的变化历程几乎一致；随着发射水深增加，出水后横向速度和俯仰角速度量值增加，而纵向速度变化较小。