基于气动耦合的导弹适配器初始分离弹道数值分析

牛钰森,姜　毅,李　静,陈　苗,史少岩,于邵祯

(1.北京理工大学 宇航学院,北京 100081;2.海军装备研究院,北京 100161)



基于气动耦合的导弹适配器初始分离弹道数值分析

牛钰森1,姜毅1,李静1,陈苗1,史少岩1,于邵祯2

(1.北京理工大学 宇航学院,北京 100081;2.海军装备研究院,北京 100161)

摘要:为了准确预测与导弹分离后导弹适配器的初始分离弹道,使用基于局部重构方法的动网格技术,将适配器的运动与流场的变化耦合。通过CFD方法实时求解适配器飞行过程中的气动载荷,并对适配器的六自由度运动微分方程进行求解。使用基于尺寸函数的局部重构方法对适配器周围网格进行更新。通过数值计算得到了适配器在初始分离弹道阶段的气动载荷系数曲线,以及适配器的分离速度、角速度曲线和分离飞行轨迹,并与试验结果进行了对比,两者高度一致。结果表明,在分离初速的作用下，适配器与导弹的距离逐渐增大，在气动力作用下适配器发生偏转，直到初始分离弹道结束适配器未与导弹发生碰撞。

关键词:适配器;CFD;六自由度;动网格;气动耦合

安装在导弹与发射箱之间的适配器能够为导弹提供支撑,并且在发射过程中起到减震与导向的作用。对于安装在导弹外表面的实心适配器,在导弹发射过程中随着导弹一起运动,出箱后通过爆炸螺栓或者弹簧机构与导弹解锁分离。在初始分离弹道阶段,适配器的飞行轨迹和姿态受气动载荷的作用时刻发生变化,有可能与导弹发生碰撞造成发射事故。因此,在发射之前有必要对适配器的初始分离弹道进行计算分析。

传统的非耦合方法计算数个攻角下适配器的气动系数,得到一系列离散的气动载荷数据,拟合出气动载荷随姿态的变化曲线,然后将曲线作为外力输入,通过求解刚体的六自由度运动微分方程计算适配器的飞行姿态与轨迹。以上方法虽具有可操作性,但是为了得到精确的计算结果,需要计算大量的工况以得到更多的数据点,耗费时间和精力。并且,在真实的环境中适配器周围的流场会随着适配器与导弹间的距离发生变化,这种变化过程是传统的非耦合计算方法不能实现的。

本文使用基于重构方法的六自由度动网格技术,将适配器的运动状态与流场的变化相耦合。每一时刻适配器受到的气动载荷由CFD程序直接求解,然后由刚体的六自由度运动微分方程计算出适配器的位移与角度变化。适配器的运动状态更新后成为新的初始条件影响下一时刻气动载荷的计算。如此循环迭代可以真正实时求解适配器在初始分离弹道阶段的飞行姿态与轨迹。使用耦合方法计算刚体在气动力驱动作用下的运动过程,已在国内外关于飞机外挂物的分离弹道研究中得到了广泛应用[2-4],其正确性得到了验证。

1适配器模型

适配器安装在导弹中部靠近重心的位置处,总共4个,并以导弹中心线为轴线按90°阵列。发射筒与地面垂直,导弹被弹射出筒,整个发射系统模型如图1所示。

图1　发射系统示意图

适配器整体呈长方形,靠近弹头一端有向内的楔形切面。内外表面均呈圆弧形,以便与导弹和发射筒壁面贴合。适配器中心位置有放置弹簧结构的盲孔,其外形如图2所示。

图2　导弹适配器几何外形

发射前适配器被发射筒和弹体约束,弹簧处于压缩状态。导弹弹射出筒的过程中适配器与导弹一起运动,当适配器离开发射筒,失去了筒壁的约束作用,弹簧开始伸展使适配器产生沿导弹径向的速度。适配器的弹出距离x与弹簧工作时间t的关系为

(1)

式中:l为弹簧的压缩量,k为弹簧的刚度,m为适配器的质量。由此计算出弹簧完全伸展开需要的时间为10 ms,在此时间内适配器的速度很小,因此气动载荷对适配器产生的冲量很小,可以忽略不计。因为适配器弹出的时间极短,所以在此时间内导弹与适配器沿轴线方向的速度也变化很小,可以认为两者的速度保持相同。以弹簧完全展开时为计算的初始时刻,此时导弹与适配器沿轴线方向的速度为20 m/s,适配器沿径向的速度为5.8 m/s。将适配器下落到导弹尾部高度时指定为初始分离弹道的结束时刻。

计算适配器的平移运动选择全局坐标系较为方便,计算适配器的旋转运动选择随体坐标系较为方便。全局坐标系以大地为基准,垂直于地面向上为X轴正方向,垂直于X轴且在流场模型的对称平面内为Y轴,Z轴由右手螺旋定则得到。适配器的随体坐标系以适配器中心孔为基准;经过重心垂直于中心孔旋转轴且在适配器对称平面内为X′轴,指向适配器楔形面为正方向;垂直于X′轴经过重心与中心孔旋转轴平行为Y′轴,指向中心孔开口侧为正方向;Z′轴由右手螺旋定则确定。

设随体坐标系由全局坐标系先绕Z轴旋转φ角,再绕Y轴旋转θ角,然后绕X轴旋转ψ角得到。则从全局坐标系到随体坐标系的转换矩阵R为

(2)

r11=cosθcosψ, r12=cosθsinψ, r13=-sinθ

r21=sinφsinθcosψ-cosφsinψ

r22=sinφsinθcosψ+cosφcosψ, r23=sinφcosθ

r31=cosφsinθcosψ+sinφsinψ

r32=cosψsinθsinψ-sinφcosψ, r33=cosφcosθ

使用Gambit建模工具以四面体非结构网格对整个流场计算域进行网格划分,网格单元数量为200万。为保证适配器气动载荷计算的精度,对适配器周围区域进行局部加密。并且在计算过程加密区域随适配器一起运动,以确保适配器壁面附近的网格密度如图3所示。

图3　适配器与导弹周围网格

2理论模型

2.1运动微分方程

适配器在空中飞行的过程中,除受到的气动力和重力外没有其他约束条件。因此,具有6个方向上的自由度。将适配器的运动分解为两部分,即质心的质点运动以及刚体的转动。

在全局坐标系下建立质心的质点运动微分方程:

(3)

式中:FX、FY、FZ分别为X、Y、Z方向上适配器受到的作用力,vX、vY、vZ分别为这3个方向上的速度分量。在适配器的随体坐标系下建立刚体转动的运动微分方程:

L·ω=∑M-ω×(L·ω)

(4)

式中:L为适配器的惯性张量,ω为角速度向量,M为作用在适配器上的力矩。Admas-Moulton多步积分法非常适合求解运动微分方程,因此使用Admas-Moulton隐式三步四阶格式对以上微分方程进行数值求解。计算过程中式(3)与式(4)中的力与力矩通过积分适配器表面的压强与切应力得到,而压强与切应力由CFD程序解算流场控制方程得到。

2.2网格重构方法

目前常用的动网格方法有:光滑方法、分层方法和重构方法。光滑方法适用于网格变形较小的情况,分层方法适用于网格运动方向恒定的情况。而适配器在空中飞行的过程中,速度方向时刻发生变化,且有滚转、俯仰和偏航运动,周围网格的变形量较大,因此选用重构动网格方法最为合适。

基于恒定尺寸阈值的局部网格重构方法灵活性较差,有可能造成刚体运动方向前方的网格尺寸过小,而运动方向后方的网格尺寸过大,降低计算精度甚至引起迭代发散。并且随着刚体的运动可能导致计算域内网格单元不断增多,使计算速度逐渐降低。而采用基于尺寸函数的网格重构方法可以有效避免以上问题。其基本原理是,根据网格单元中心点与边界的距离以及边界网格单元的尺寸,计算此单元网格重构时的尺寸基准。重构的网格尺寸为边界网格单元尺寸的加权平均,其权值是重构网格单元与边界网格单元距离的反比,即:

sp=sb×γ

(5)

式中:

(6)

(7)

式中:FI与边界上的单元格尺寸相关;LI为网格单元中心点距离节点的距离;db是一个无量纲量,表征网格单元与边界单元之间的距离;α为尺寸函数变化率,控制网格单元尺寸相对于边界网格单元的大小;β为尺寸函数比例,控制网格单元尺寸随db的变化程度。当运动边界靠近时,db减小,sp也减小,此时重构单元网格加密,当运动边界远离时,db增大,sp也增大,此时重构单元网格稀疏。

2.3流场控制方程

在基于欧拉描述的流场中,任意空间位置处的标量φ在当地的时间变化率是由:流入流出此控制体的对流流量,温度梯度、浓度梯度等带来的扩散流量以及其他现象产生的源项引起的。因此标量φ的输运方程可表示为

(8)

p=ρRT

(9)

h=∫cp(T)dT

(10)

式中:ρ为当地流场的密度,u为流场的速度矢量,为拉普拉斯算子,Γ为扩散系数,Sφ为标量φ的产生源项。当φ=1时,可以得到流场的质量守恒方程;当φ分别取X、Y、Z方向上的速度标量u、v、w时,可以得到动量守恒方程;当φ取比焓h时,可以得到能量守恒方程。以上5个方程再加上描述气体状态的理想气体状态的式(9)以及描述比焓与温度之间关系的式(10),就可以构成封闭的控制方程组,对基本标量ρ、u、v、w、p、T与h进行求解。

2.4Spalart-Allmaras湍流模型

在适配器的气动流场中,适配器的楔形面与立面的交接过渡处存在曲率较大的流线以及较强的压强梯度,当适配器的攻角较大时还会存在贴体流动分离现象。对于这样的气动湍流问题,Spalart-Allmaras模型经过一定的优化,是专门为气动流场设计的湍流模型。且在本文研究的适配器初始分离弹道问题中导弹发动机没有点火,因此不存在燃气射流这样的自由剪切流,选用此模型作为湍流粘性模型对于本文研究的问题十分合适。Spalart-Allmaras模型是一方程湍流粘性模型,通过输运方程求解湍流粘性νT,其表达式为

(11)

式中:σν为常数,Sν为湍流粘性源项。

3计算结果

初始时刻流场内为标准大气环境,压强为101.325kPa,温度为293.15K。流场中有沿Y轴正方向速度为3m/s的横风。按照初始时刻全局坐标系下适配器质心位置给适配器命名,质心在Y轴正方向的为“SPQ-1”,质心在Y轴负方向的为“SPQ-2”,质心在Z轴正方向的为“SPQ-3”,质心在Z轴负方向的为“SPQ-4”。

3.1适配器气动载荷变化

初始时刻适配器“SPQ-1”、“SPQ-2”迎风面的压力分布云图如图4所示。

图4　初始时刻适配器迎风面压力分布

由图4可以看出，“SPQ-2”迎风面平均压强要高于“SPQ-1”。选取适配器楔形面在随体坐标系O′Y′Z′平面上的投影面积S为参考面积,选取适配器最长边的一半为参考长度l,由公式(12)可得适配器无量纲参数气动载荷系数。

(12)

式中:CF为作用在适配器上的气动力系数,CM为作用在适配器上的气动力矩系数,q∞为无穷远处自由流场的动压。适配器气动力系数CFX、CFY、CFZ在初始分离弹道阶段随时间变化曲线如图5所示。

图5　气动力系数曲线

适配器有沿X轴正方向的速度,因此受到沿X轴负方向的气动力作用。可以看出4个适配器X方向的气动力变化趋势一致。从初始时刻到0.15s这段时间内适配器气动力系数CFX幅值不断增大,因为受图5(b)、图5(c)所示的气动力矩作用,在此时间内适配器绕O′Z′转动,随着旋转角度的增大,适配器在OYZ平面内的投影面积不断增大,使X方向上的气动力作用面积增大,使FX幅值上升。0.15s左右4个适配器绕O′Z′旋转角度达到90°,适配器处于水平状态如图11(c)所示,之后适配器继续旋转,OYZ平面内的投影面积不断减小,同时受FX作用vX不断减小,2个因素共同作用下FX幅值下降如图5(a)所示。在0.3s左右时适配器的旋转角度达到180°,如图11(d)所示,此时OYZ平面内的投影面积达到最小,且vX小于初始时刻,因此FX幅值最小。在Y、Z方向上气动力的变化过程与此类似。

不同时刻适配器表面的压力分布如图6所示。初始时刻,适配器处于竖直状态,此时X方向气动力主要作用在适配器前端的楔形面上,如图6(a)所示。楔形面外法线与O′X′轴有夹角,因此楔形面受到的气动力有沿O′Y′轴正方向分量,在适配器质心处产生的力矩沿O′Z′轴正向,这是致使适配器旋转的主要因素。随体坐标系下适配器受到的气动力矩系数CMZ′随时间变化曲线如图7所示。

图6　适配器表面压力云图

图7　气动载荷系数CMZ′曲线

可以看出，在最初的一段时间内适配器受到的气动力矩有所上升,这是因为适配器绕O′Z′轴旋转后增大了X方向的气动力作用面积,适配器迎风面位于质心之前的区域提供了更多的沿O′Z′轴正向的气动力矩,但是位于质心之后的区域产生的是沿O′Z′轴负向的气动力矩,两者叠加气动力矩幅值有较小的上升。当适配器转过一定的角度时,质心后部区域产生的负向气动力矩超过前部正向气动力矩。并且随着适配器转角的增大，楔形面外法线与来流之间的夹角不断增大,由此产生的正向气动力矩逐渐减小,在2个因素的共同作用下致使适配器沿O′Z′轴正向的气动力矩随转角的增大而减小。在0.15s作用时适配器旋转至水平位置,此时CMZ′达到最小值,之后适配器继续旋转适配器背面和底部成为迎风面,从而产生了正向气动力矩使CMZ′有所上升,但是此时适配器的运动速度已经减小很多,因此上升幅度不大。

3.2适配器分离弹道轨迹

由图5、图7可以看出，适配器“SPQ-3”与“SPQ-4”的气动载荷系数体现出较好的对称性与一致性,而适配器“SPQ-1”与“SPQ-2”的气动载荷系数有较大的差别。这是因为适配器“SPQ-1”的分离速度沿Y轴正方向与风的方向相同,而适配器“SPQ-2”的分离速度沿Y轴负方向与风的方向相反,因此两者相对于空气的运动速度不同,造成两者的气动载荷也不相同。而气动载荷的非对称性又加剧了两者Y方向上速度的差异。适配器质心运动速度随时间变化曲线如图8所示。

可以看出,4个适配器在X方向上速度变化趋势保持一致。适配器“SPQ-2”直接受到从Y轴负方向吹来的风影响,vY的减小速度快于适配器“SPQ-1”。同时受风力作用,本来在Y方向没有速度的适配器“SPQ-3”与“SPQ-4”也逐渐产生了速度。在Z方向上没有风力作用,因此适配器“SPQ-3”与“SPQ-4”的速度vZ变化趋势完全一致,且适配器“SPQ-1”与“SPQ-2”在Z方向上的速度始终为0。适配器质心运动轨迹呈抛物线形如图9所示。

图8　适配器质心速度曲线

图9　适配器质心轨迹

适配器与导弹分离0.575s后导弹底部在X方向上超过所有适配器,初始分离弹道结束。设此时适配器质心与导弹轴线的垂直距离为s,结束时所有适配器的速度与位移见表1。

表1　初始分离弹道终点适配器运动参数

分离过程中绕随体坐标系O′Z′轴正方向的旋转是主要的旋转运动,致使适配器前端楔形面向弹体方向翻转,旋转角度ψ与角速度ωZ随时间的变化曲线如图10所示。

初始分离弹道结束时所有适配器的旋转角速度与角度见表1。数值计算得到的适配器运动姿态与实验中高速摄影记录的不同时刻适配器姿态变化对比如图11所示。

由数值计算结果与试验结果均可以看出,初始分离弹道结束时适配器未与导弹发生碰撞,适配器的分离是安全的。此外,数值计算预测的适配器飞行姿态和轨迹与高速摄影记录的实验结果保持一致。

图10　适配器旋转运动变化曲线

图11　适配器飞行姿态

5结束语

本文建立了适配器与导弹的初始分离弹道模型,划分了非结构网格,使用基于六自由度动网格技术将适配器的分离运动与气动流场相耦合,得到了初始分离弹道阶段适配器的运动轨迹与姿态,并与实验结果进行了对比,可以得出以下结论:

①计算结果与实验结果高度一致,表明此种研究方法具有较高的可信度和准确性;

②在本文研究的工况中, 适配器在初始分离阶

段没有与导弹发生碰撞,分离过程是安全的;

③受气动载荷的作用适配器产生旋转运动,且适配器前端楔形面向弹体方向翻转;

④沿着风速方向对称分布的2个适配器,受到的气动载荷不同,两者的运动状态也不同。

本文使用的研究方法对于研究适配器与导弹的分离过程,预防适配器与导弹的碰撞,以及地面设备和人员的防护有指导作用。

参考文献

赵华,王敏杰,杨为,等.箱式发射导弹适配器.战术导弹技术,2007(4):42-50.

ZHAO Hua,WANG Min-jie,YANG Wei,et al.Adapters for canister-launched missile.Tactical Missile Technology,2007(4):42-50.(in Chinese)

MAHMOOD T,AIZUD N,ZAHIR S.CFD simulations of the store separation and its effect on the attitude of the parent vehicle//International Bhurban Conference on Applied Sciences.Piscataway,NJ,USA:IEEE,2012:198-202.

PANAGIOTOPOULOS E E,KYPARISSIS S D.CFD transonic store separation trajectory predictions with comparison to wind tunnel investigations.International Journal of Engineering,2010(3):538-553.

张培红,王明,邓有奇,等.武器分离及舱门开启过程数值模拟研究.空气动力学学报,2013,31(3):277-293.

ZHANG Pei-hong,WANG Ming,DENG You-qi,et al.Numerical simulation of store separation and door operation.Acta Aerodynamic Sincia,2013,31(3):277-293.(in Chinese)

BURDEN R L,FAIRES J D.Numerical analysis.Youngstown:Youngstown State University,2010.

POPE S B.Turbulent flows.UK:Cambridge University Press,2000.

ANDERSON J D.Fundamentals of aerodynamics.5th ed.New York:McGraw-Hill,2009.

Numerical Analysis of Initial Separation Trajectory of

Missile Adapter Based on Aerodynamic Coupling

NIU Yu-sen1,JIANG Yi1,LI Jing1,CHEN Miao1,SHI Shao-yan1,YU Shao-zhen2

(1.School of Aerospace Engineering,BIT,Beijing 100081,China;

2.Naval Acadme of Armament,Beijing 100161,China)

Abstract:In order to accurately predict the initial separation trajectory of the adapter after releasing from the missile,the motion of the adapter was coupled with the variation of the flow field by dynamic mesh technology based on local reconstruction.The aerodynamic loads of the adapter were computed in real time by CFD routine.The 6DOF motion differential equations of the adapter were solved.The models of missile and four adapters were built,and the computational domain was meshed by using tetrahedron cells.The grid cells around the adapter were updated by using local remeshing method based on size function.The initial separation trajectory of the adapter was calculated.The results accord with the experimental result.The curves of motion velocities and aerodynamic loads coefficients were drawn.Results show that the distance between adapter and missile increases because of the initial separation velocity,and the adapters rotate because of the aerodynamic force.The adapter doesn’t collide with missile untile the end of initial separation trajectory.

Key words:adapter;CFD;six degree of freedom;dynamic mesh;aerodynamic coupling.

中图分类号:TJ303.4

文献标识码:A

文章编号:1004-499X(2015)04-0052-07

作者简介:牛钰森(1987- ),男,博士研究生,研究方向为兵器发射理论与技术。E-mail:shenzhou1987@aliyun.com。

收稿日期:2015-06-11