助推-滑翔飞行器弹道分段优化研究

李 柯,聂万胜,冯必鸣

(装备学院,北京 101416)

助推-滑翔飞行器具有远程快速精确打击、大范围区域覆盖、机动性好等优点,已成为国内外研究的热点。其整个飞行过程可描述为:固体火箭从地面垂直发射,经过垂直上升和转弯过程后发动机关机,飞行器与火箭和整流罩分离,这段过程称为主动段;飞行器与箭体分离后,依靠气动控制实现远距离无动力滑翔,完成对目标的打击,这段过程称为滑翔段。

整个飞行过程中,飞行器可以达到的极限距离为最大射程,是评价助推-滑翔飞行器性能的一项重要指标。文献[1]利用解析法近似分析了最大射程和升阻比及关机点速度的关系。文献[2]利用Gauss伪谱法对 CAV再入滑翔段的轨迹进行了优化,突出了伪谱法在计算效率和精度上的优势。文献[3]基于Gauss伪谱法对固体火箭上升段进行了优化研究。文献[4-5]设定少量的约束条件,通过序列二次规划法对整体过程进行了优化。总的来讲,目前对全程弹道优化研究的仍然很少,尤其对于主动段关机点指标的选取仍不明确。本文在上述文献的基础上,利用Radau伪谱法分别对主动段和滑翔段进行优化,然后将不同优化方案求得的最大射程进行了对比分析,明确了主动段的性能指标,给出了最佳方案。

1 数学模型

1.1 动力学模型

在方案论证阶段可忽略次要影响因素,因此不考虑地球自转、扁率和弹体自身滚转,其简化的主动段运动方程为

其中,V为飞行速度, T为发动机推力,D为阻力, L为升力,α为攻角,θ为弹道倾角,s为航程,h为高度, m为火箭总质量,R0为地球平均半径,g为重力加速度,I为发动机比冲。飞行器处于滑翔段时,自身不携带发动机,质量为定值,所以将T＝0,d m / d t ＝ 0 代入式(1)即得到滑翔段的运动方程。

式中,密度ρ、重力加速度g,阻力D和升力L可由式(2)计算得出:

其中,0ρ为海平面大气密度,H为常量,μ为引力常数,CD为阻力系数,CL为升力系数,S为气动参考面积。

1.2 气动模型

固体火箭采用轴对称结构,其气动参数可按照文献[6]选取:

飞行器采用升力体构型,气动外形如图1所示。以飞行高度30km为例,压强1170Pa,温度216.5K,大气密度 1.9×10-5kg/m2,计算得到不同马赫数状态下,升力系数、阻力系数和升阻比随攻角的变化,如图 2至图4所示。

图1 飞行器的气动外形

图2 升力系数随攻角变化情况

图3 阻力系数随攻角变化情况

图4 升阻比随攻角变化情况

由图2、图3可知,升力系数随攻角的变化曲线近似为一次函数,阻力系数随攻角的变化曲线近似为二次函数。在马赫数高于9时,气动特性曲线非常接近,这与高马赫数无关原理相吻合。因此,气动系数可拟合成攻角的函数,即

1.3 优化模型

不失一般性,考虑 Bolza形式的最优控制问题,设状态变量 x （ τ）∈ Rn,控制变量 u （τ）∈ Rm,初始时间t0,终端时间tf,指标函数为[7]

受材料和自身结构限制,必须使驻点热流、动压和过载始终处于飞行器允许的范围之内,即路径约束:

状态变量的初始值和终端值要根据初始条件和末端指标来确定,即边界条件:

其中包括:初始参数h0、0θ、s0、V0、m0,末端参数hf、fθ、sf、Vf、mf。

攻角作为控制量,应处于飞行器所允许的最大可控范围之内,即控制量约束:

2 Radau优化算法

Radau伪谱法将状态变量和控制变量在一系列LGR(Legendre-Gauss-Radau)点上离散,并以离散点为节点构造Largrange插值多项式来拟合状态变量和控制变量。通过对全局插值多项式求导来近似状态变量对时间的导数,将动力学微分方程约束转化成一组代数约束,将轨迹优化问题最终转化为非线性规划问题[8]。

由于LGR点所定义的时间区间 τ ＝ [-1 ,1),若要把动力学模型在LGR点上进行离散,需要把优化问题的时间区间 t ＝[t0, tf)通过式(11)转换:

转换后,τ取代t成为独立变量,则指标函数,微分方程,约束条件和边界条件分别为

状态变量和控制变量用Lagrange插值多项式近似表示:

经过上述离散化,指标函数、微分方程、路径约束和边界条件分别为

其中, wk,（ k ＝ 1 ,… … ,N -1）为LGR权重。

最终将连续问题转化为非线性规划问题进行求解。

3 仿真及分析

在主动段,初始条件设为h0=0km,V0=0m/s,θ0= 90°,由于整流罩的存在,可不考虑热流密度约束,仅考虑动压和法向过载约束,qmax＝ 1 00kPa,nymax＝ 1 .5。洲际固体导弹的攻角一般不超过7°,所以控制量约束设为 - 7°≤ α ≤ 7°。假设发动机关机时速度为5000m/s,以此设计出两级固体火箭的性能参数如表1所示,具体设计方法参考文献[9]。

由于发动机的关机点参数(即滑翔段的初始参数)直接决定了滑翔弹道的形状和射程,因此,根据不同的性能指标,主动段采取不同的方案进行优化,然后从中挑选出最佳方案。

方案1:主动段飞行距离最远 J ＝ m axsf;

方案3:关机速度最大 J ＝ m axVf。

不同方案优化出的关机点参数如表2所示。

表1 固体火箭性能参数

表2 关机点参数

在滑翔段,设飞行器气动参考面积 S ＝ 0 .48387m2,质量m ＝ 9 07.18kg,将主动段的关机点参数作为初始值,仿真终止条件为 hf=0km 。需要考虑路径约束包括动压约束、过载约束和驻点热流密度约束,即控制量约 束 为0≤ α ≤ 3 0°。 性 能 指 标 为 滑 翔 距 离 最 远

如图5和图6所示,3种方案的弹道形状有很大差异。方案 1中,滑翔段的初始倾角较小,弹道平滑,跳跃幅度小,飞行高度都保持在 30～50km,该高度为大部分防空导弹的拦截盲区,利于突防,攻角变化较为平滑,易于控制,但射程最小;方案2中,滑翔段的初始倾角较大,弹道的跳跃幅度较大,自由飞行弹道(近似于弹道导弹的中段飞行)较长,不利于突防,而且在200～543s这段时间内,攻角出现剧烈波动,不利于控制,射程比方案1略大;方案3中,弹道跳跃幅度适中,攻角最大值为 16.8º,没有出现剧烈波动的情况,且射程最远。由此可见,主动段的终端倾角决定滑翔段的跳跃幅度,主动段的终端速度决定最大射程。

图5 不同方案的弹道曲线

图6 不同方案的控制量变化曲线

图7 主动段路径约束变化情况

图8 滑翔段路径约束变化情况

表3 路径约束最大值

由表3可知,3种方案的路径约束最大值都满足本文所设的约束条件,验证了方案的可行性。对比方案2和方案3,在热流密度最大值相同的情况下,方案3的射程最远,此外,方案2在滑翔段的路径约束最大值更接近极限值,这对飞行器的自身承受能力要求较高,因此方案3优于方案2。对比方案1和方案3,方案1在滑翔段的路径约束最大值皆小于方案3,但在主动段的法向最大过载接近极限值,这对固体火箭的结构强度要求较高,同时射程也小于方案3,可见方案3优于方案1。综合对比发现,方案3射程最远,实现难度适中,可选为最佳方案。

文献[1]利用解析法,得到的最大射程计算公式为

式中,L/ D为飞行器的升阻比,Vf为主动段的关机点速度,Vs为第一宇宙速度,可见最大射程与飞行器的升阻比成正比。令 L / D ＝ 3 .681,即最大升阻比,求得方案3的最大射程解析值如表4所示。

由表 4可知,通过本文的设计方案,优化后最大射程比解析值增加了9.8%,优化方案可行。

4 结束语

本文利用Radau伪谱法,将轨迹优化问题转化为非线性规划问题,分别优化了主动段和滑翔段的弹道,仿真计算了不同方案时的飞行器最大射程,得出了以下结论:

1)性能指标选取为主动段关机点速度最大和滑翔段射程最远的组合方式,可以使总射程达到最大,实现难度适中;

2)主动段的终端倾角决定滑翔段的跳跃幅度,主动段的终端速度决定最大射程;

3)采用分段优化的方法处理助推滑翔问题,便于分析各段的弹道特征,为技战术指标的制定提供参考。

[1]Eggers A J, Allen J H. A comparative analysis of the performance of long-range hypervelocity vehicles[R].NACA Technical Report, March l954.

[2]雍恩米.高超声速滑翔式再入飞行器轨迹优化与制导方法研究[D].长沙:国防科技大学, 2008.

[3]杨希祥.基于Gauss伪谱法的固体运载火箭上升段轨迹快速优化研究[J].宇航学报,2011,32(1):15-21.

[4]李瑜.助推-滑翔导弹弹道优化研究[J].宇航学报,2008,29(1):66-71.

[5]李瑜.助推-滑翔导弹最大射程优化[J].弹道学报,2008,20(4):53-56.

[6]何麟书.弹道导弹和运载火箭设计[M].北京:北京航空航天大学出版社,2002.

[7]胡寿松.自动控制原理[M].北京:科学出版社,2004.

[8]Divya Garg, Michael A Patterson. Direct trajectory optimization and costate estimation of general optimal problems using a radau pseudospectral method[C].AIAA Guidance, Navigation, and Control Conference.AIAA, August 2009.

[9][苏]A.M.西纽科夫.固体弹道式导弹[M].北京:国防工业出版社,1984.