超声速流场中子母弹分离与子弹姿态变化规律的数值模拟

陈时通,于 勇

(北京理工大学 宇航学院,北京 100081)

子母弹是现代武器中很重要的一种大规模杀伤性武器,它在飞行的过程中从母弹中弹射出多发子弹,造成大规模、大面积的毁伤[1]。子母弹的作战效能由多个子弹的散布效果决定。子母弹研制过程中面临的重要问题之一是使多个子弹正常分离,达到预期的散布效果。而要解决这一问题,首先要保证子弹分离过程中姿态稳定。

国外关于子母弹的研究资料较少。Edge等[2]自编CFD程序,采用嵌套网格,研究了子母弹抛撒过程中的激波相互作用。Panneerselvam等[3]通过风洞实验研究了子母弹在弹舱内和弹舱外的法向力系数、力矩系数和压心位置的变化。Corder等[4]使用自编程序研究了带弹翼和不带弹翼的子弹在不同攻角下的法向力系数和俯仰力矩系数的变化。Dietz[5]使用软件Overflow模拟了24枚子弹抛撒前的马赫数、压力分布以及抛撒后各子弹的位移和姿态变化。 Deep等[6]通过自编程序研究了有、无弹舱的情况下,子弹在不同马赫数及不同攻角情况下的法向力系数及俯仰力矩系数并与实验对比,吻合良好。

国内对子母弹的研究起步较晚。邹德坤等[7]使用CFX模拟了子弹在不同攻角、不同马赫数下的升力系数和阻力系数。王金龙等[8]将气囊实验的结果运用于FLUENT数值仿真,研究了在燃气干扰下不同来流马赫数对子母弹的流场结构和子弹的运动特性的影响。也有学者研究母弹的壳片在分离过程的运动特性。蒋增辉等[9]使用风洞投放模型试验技术研究了母弹处于大迎角(25°)状态下壳片的分离运动规律。王巍[10]利用弹簧近似和网格重构相结合的非结构动网格技术,耦合求解N-S方程及弹道方程,模拟了不同的壳片分离速度和角速度的抛壳过程。陶如意等利用AUSM+格式求解采用k-ωSST湍流模型的雷诺平均N-S方程,研究了子母弹时序抛撒过程中的子弹与子弹和子弹与母弹间的激波相互作用[11];通过风洞实验研究了子弹在母弹流场中不同位置、不同攻角下的阻力系数、升力系数和力矩系数的变化[12]。

实验研究可以给出子母弹抛撒过程中的气动参数及流场状况,但实验成本大,周期长。与实验相比,数值模拟成本低,周期短。已有大量的学者通过数值模拟研究子母弹抛撒过程,研究多采用动网格技术,使用非结构网格进行网格重构。但是,所带来的问题是网格数量多,网格重构计算量大,而且重构的网格质量会逐渐降低,严重时甚至使计算中断。因此,本文在计算方法上试图使用嵌套动网格技术,采用结构化网格。在网格边界运动时,不需要重新生成网格,从而提高网格质量,减少计算时间。同时耦合动网格技术和6DOF方程计算分离过程中的子弹轨迹与姿态,利用这种非定常的模拟方法研究了母弹上不同弹仓位置对子弹抛撒后的4枚子弹姿态的影响。

1 数值计算方法

1.1 流场求解

雷诺平均N-S方程在笛卡尔坐标系下的描述为

(1)

式中:Ω为控制体体积;S为控制体表面积;n为表面S的外法向单位矢量;dV为体积分微元;dS为面积分微元;Q为守恒变量,F(Q)为对流项,H(Q)为黏性通量。

1.2 子弹运动方程

子母弹运动6DOF方程为

(2)

(3)

(4)

(5)

本文使用的是Fluent18.0,采取k-ω湍流模型,耦合6DOF方程(弱耦合),使用的求解方法是隐式的Roe-FDS方法,离散方程是二阶迎风格式,选取密度求解器。

1.3 嵌套网格技术

随着工业技术的发展,工程项目模型更加复杂,网格质量的好坏成为了工程模拟能否进行的关键。工程项目中,网格的生成与修改往往要占据80%以上的时间,项目模型越复杂,网格生成越困难,需要的时间也越多。而且,在模拟的过程中,经常因为网格质量差导致模拟失败,最终工程项目只能暂停。因此,计算网格生成技术的发展成为了解决问题的关键[16]。

嵌套网格技术(Fluent中称为Overset)能够很好地解决模型复杂、网格质量不高的问题,同时还具有部件网格生成容易、网格数量少,网格质量高且高可移植性的特点。嵌套网格技术将模型拆解成一个背景网格和多个简单的组件网格。每套网格需单独划分结构化网格,并在Fluent中进行插值,网格重叠的边界处采用三线性差值的方式传递数据。组件网格运动的过程中,每运动一次都需要重新生成新的嵌套网格,保证生成网格具有高质量。同时,背景网格和组件网格单独划分,可根据工程需要随时更改组件形状,提高网格生成效率[13]。

1.4 嵌套网格技术验证

机弹分离问题是多体分离中的经典问题,本节将对机弹分离问题进行模拟,并与风洞试验进行对比,以此验证计算方法的正确性。算例来流静压p=20 657 Pa,静温T=216.65 K,Ma=1.2,攻角α=0°,仿真时间步长Δt=0.001 s。表1为机弹分离问题参数表。

表1 机弹分离问题参数

表1中,ms,Ls,Ds分别为弹体的质量、长度和直径;k为弹翼数;ds为质心位置;Ix,Iy,Iz为分别x,y,z方向转动惯量;FFE,FAE为弹射力。详细模型尺寸和实验结果可参考文献[14]。

图1为仿真模型图。

图1 机弹分离模型

图2、图3分别为导弹质心位移与姿态角随时间的变化规律,图中,D为质心位移,λ为姿态角(姿态角λ代表了滚转角γ、偏航角ψ和俯仰角θ)。图中导弹滚转角与实验值偏差较大,因为导弹轴向转动惯量较小,模拟较困难,需增加导弹周向的网格节点,但会极大地降低计算效率。总体来说,仿真数据与实验数据符合较好。机弹分离的经典算例模拟结果证明了计算方法的有效性。

图2 导弹质心位移随时间的变化

图3 导弹的姿态角随时间的变化

2 子母弹分离过程的数值模拟

本文计算模型为1枚母弹和8枚子弹装配,由于子弹对称分布,因此模型只考虑右边部分。本文采用嵌套网格对子母弹分离过程进行数值模拟,子母弹工作高度为8 km,母弹攻角为0°,且无旋转,飞行马赫数为2.43,静压为35 652 Pa,静温为236.2 K,子弹弹射速度为50 m/s(弹射沿径向方向均匀作用在子弹质心上),仿真时间为0.03 s。初始状态下弹舱已打开,子弹暴露在气流中。非定常计算时,给定子弹初始弹射速度,子弹做无控运动。

2.1 子母弹模型

子母弹结构设计分为2种方案:弹舱靠前设计(工况A),弹舱靠后设计(工况B)。

弹舱靠前设计及母弹尺寸如图4所示,子弹布局如图5所示,子弹尺寸如图6所示。弹舱靠后设计及母弹尺寸如图7所示。

图4 弹舱靠前设计母弹尺寸(单位:mm)

图5 子弹布局(单位:mm)

图6 子弹尺寸(单位:mm)

图7 弹舱靠后设计及母弹尺寸(单位:mm)

建立2个坐标系:母弹坐标系为定坐标系,以母弹质心为原点,以母弹弹轴为X1轴,弹尾指向弹头方向为X1轴正方向。Y1轴在母弹模型纵向对称面且与X1轴垂直,指向上与重力方向相反为Y1轴正方向。Z1轴由右手法则确定。子弹坐标系为动坐标系,以子弹质心为原点,子弹弹轴为X2轴,弹尾指向弹头方向为X2轴正方向。Y2轴始终垂直于X2轴且指向上为正方向,初始时刻Y2轴方向与重力方向相反。Z2轴由右手法则确定。4枚子弹分别位于母弹的不同周向位置,图8为子母弹装配图。

图8 子母弹装配图

2.2 边界条件

数值计算中,母弹速度为超声速,外部空气边界条件为压力远场,母弹设置为壁面,外部空气边界条件设置为Overset,子弹设置为壁面。由于母弹的弹头形状是小圆弧,弹尾是圆柱形,因此,计算域形状与母弹形状类似,计算域头部为半圆形,尾部为圆柱形,计算域半径为2 m,长度为11 m。计算域形状如图9所示。

图9 计算域形状

2.3 计算网格

图10是定常状态和非定常状态下某一时刻的Overset网格。

图10 子母弹Overset网格

由于计算采用的是嵌套网格,因此母弹和子弹需分别划分结构化网格。子弹和母弹靠近壁面的第1层网格高度为1 mm,子弹部件网格数量为21万,母弹网格数量为170万。母弹周围的网格为背景网格,子弹周围的网格为组件网格。在Fluent中选定背景网格和组件网格进行插值,重叠区域插值方式为最小二乘插值。计算网格壁面y+约为200,时间步长为10-4s,残差收敛达到10-3。在非定常状态下,每一个时间步组件网格要更新移动,嵌套网格也会进行重构,保证每一次迭代的网格质量。

3 结果分析

本文分别对子弹弹舱在前部(工况A)或后部(工况B)进行数值模拟。

3.1 子母弹压力场

子弹弹射出舱之前,需要计算定常状态下子弹在母弹舱内的绕流流场,并以此为初始流场计算子母弹的分离过程。图11(a)为定常状态下子母弹的压力场,图中,p为静压。从图中可以看出,在定常状态下,母弹弹头处压力较高,头部形成马赫波,从母弹头部到尾部,压力降低,在母弹弹舱后部位置形成马赫波。子弹尾部位于弹舱马赫波处,压力较高。

图11 子母弹压力场

图11(b)为工况A非定常状态下子弹流场压力分布,仿真时间t=0.03 s。图中截面为子弹2、子弹3的纵向截面,下同。从图中可以看出,子弹弹头处压力较高,形成激波且子弹形成的激波与母弹形成的激波相交。图11(c)为工况B非定常状态下子弹流场压力分布,仿真时间t=0.03 s。从图中可以看出,子弹弹头处也形成激波,但由于弹舱位置远离母弹激波,没有出现与母弹激波相交的状况。图11(d)为工况B非定常状态下子弹流场压力分布,仿真时间t=0.01 s。从图中可以看出,由于仿真时间较短,子弹位移较小,因此发生了子弹间的激波相交。子母弹分离的压力场表明,若弹舱位置在后部,子母弹分离初期子弹间发生激波相交;子母弹分离后期,由于子弹间距较大,没有发生激波相交。

3.2 弹舱位置对子弹姿态的影响

3.2.1 子弹1、子弹4的姿态变化

图12为子弹1的姿态角随时间的变化曲线。从图12可知,工况A时子弹1的姿态角变化较工况B更大。因为弹舱位置靠前,接近母弹激波,母弹激波附近压力较大,对子弹1的姿态影响较大。子弹1的近激波面压力较大,使子弹1的姿态发生偏转。在母弹激波作用下,子弹1的偏航角正向增大到0.7°,俯仰角负向增大到-1.7°。工况B时,由于子弹1位置远离母弹激波,在子弹2激波和空气动力的作用下,子弹1的偏航角和俯仰角负向增大到-0.5°左右,小于工况A。

图12 子弹1的姿态变化

图13为子弹4的姿态角随时间变化曲线。从图13中能够得出与子弹1类似的结论。

3.2.2 子弹2、子弹3的姿态变化

图14为子弹2的姿态角随时间变化曲线。子弹2的位移变化规律与子弹1相同。从图14可知,工况A时母弹激波对子弹2的姿态影响较大。因此,在母弹激波的作用下,子弹2的偏航角正向增大到1.7°,俯仰角负向增大到-0.7°。工况B时,子弹2受到子弹1、子弹3的激波和空气动力的作用,俯仰角正向增大到1.8°。

图13 子弹4的姿态变化

图14 子弹2的姿态变化

图15为子弹3的姿态角随时间变化曲线。子弹3的位移变化规律与子弹1相同。从图15中可以看出,在工况A时子弹3变化规律与子弹2相同。在工况B时子弹3偏航角与俯仰角变化都较大,与子弹2不同。工况A时,母弹激波对子弹3的姿态影响较大,使子弹3的偏航角正向增大到1.7°,俯仰角正向增大到0.7°。工况B时,子弹3受到子弹2和子弹4的激波的相交作用和空气动力作用,导致子弹3的俯仰角正向增大到1.8°,偏航角负向增大到-1.8°。

图15 子弹3的姿态变化

图16(a)、16(b)为工况A子母弹在t=0.01 s和t=0.02 s时的压力场,截面经过子弹质心,方向为X1轴方向。从图中可以看出,子母弹压力场对称,因此,子弹1、子弹4和子弹2、子弹3的姿态运动规律对称。

图16(c)、16(d)为工况B子母弹在t=0.01 s和t=0.02 s时的压力场。从图中可以看出,子弹1、子弹4的压力场轴对称,子弹2、子弹3的压力场不是轴对称。子弹2、子弹3的背激波面的低压区方向不同,造成了子弹2、子弹3的姿态变化不同。

图16 工况A、工况B不同时刻子母弹压力场

4 结论

采用嵌套网格的方法对子母弹的分离进行仿真与分析,由本文分析可知,采用嵌套网格技术和k-ω湍流模型,同时耦合6DOF方程,能够很好地模拟子母弹分离的姿态变化过程及其干扰流场。

子母弹分离过程中,若弹舱位置在前部,子弹分离过程始终出现子弹激波与母弹激波的相交,子弹受到母弹激波影响姿态变化比弹舱位置在后部更大。若弹舱位置在后部,子母弹分离初期出现子弹激波的相交,中间位置对称分布的子弹在抛撒过程中受到子弹激波和不对称压力场的作用,姿态会发生不对称的变化。无论弹舱位于前部还是后部,上下位置对称分布的子弹的压力场是对称的,运动姿态也对称。