0 引言
对于大部分运载火箭而言,助推器在工作结束后需要进行分离以降低结构重量,其中绝大多数飞行器的助推器不具备姿态控制能力,将在无控状态下自由再入。如今随着航天发射任务的高密度化,助推器再入后带来的安全性问题愈发受到关注,助推器落点范围的精确预示及落点范围缩减方法成为重点研究的方向。
传统工程上进行助推器落点范围预示的方法主要有极限偏差组合法和六自由度打靶法2种,极限偏差组合法不需要进行大量的气动数据计算,可利用少量典型状态弹道快速进行落区范围评估,但是该方法假设助推器飞行中全程按照阻力系数的极大或极小2种姿态飞行,工况过于保守,落点预示范围过大。六自由度打靶法是目前实际工程中最常用的方法,需要三通道大姿态角范围的助推器气动数据,往往需要大量的数值计算和风洞试验数据,并考虑各项偏差组合工况进行蒙特卡洛打靶仿真,解决了极限偏差组合法预示范围过于保守的问题,但是条件保障要求高、仿真周期长、难以辨识运动规律、难以指导落点范围缩减的正向设计。
助推器无控再入过程中,仅在重力和气动力的作用下进行自由运动,其动力学系统是一个自治非线性系统,理论上存在动态平衡状态即固有运动模态。例如美国Ares I火箭助推器安控区设计过程中,研究人员通过运动仿真分析辨识出助推器再入姿态运动存在的不同运动模态并对应不同落点散布规律。在自旋稳定弹丸的研究中,工程人员也通过分析弹丸的运动模态(动力平衡角)实现了对弹丸弹道散布程度的预估。可见自由状态的飞行器在一定条件下存在固有的运动模态,不同运动模态对应特定的质心受力条件,对应再入过程可持续保持的极限运动状态,因此可通过运动模态的预示实现落点范围的预示。
为了准确、高效的实现助推器无控再入范围的预示,并通过辨识其再入过程运动规律探索经济的落点范围缩减设计方案,本文中采用基于粘性扰流和细长体理论的工程方法快速获取助推器大范围气动特性,基于非线性动平衡点理论辨识助推器动姿态运动动平衡状态即飞行模态,进一步通过不同模态的质心等效受力模型快速实现落点范围的预示,并对助推器关键设计参数对落点范围的影响进行了敏感度分析,为助推器落点范围快速精细化设计提供了有益的参考。
1 助推器无控再入运动学模型
本文中建立以经典欧拉角及其角速度描述的姿态动力学模型,开展相关分析工作。助推器体坐标系(Oxbybzb)原点为O位于质心,Oxb轴沿助推器纵轴向前为正,Oyb在对称平面内指向上为正,Ozb轴与其余两轴满足右手定则。准速度坐标系(Ox5y5z5)坐标系原点为O位于质心,Ox5轴沿速度方向向前为正,Oy5位于包含Ox5的竖直平面内向上为正,Oz5轴与其余两轴满足右手定则。经典欧拉角的定义如图1所示。ζ定义为进动角,δ定义为章动角,λ定义为自转角。
图1 经典欧拉角定义
Fig.1 Classic Euler angle definition
考虑到助推器姿态运动相对质心运动为快时变运动,在分析姿态运动时可以认为助推器速度的大小及方向保持不变,准速度系可视为惯性坐标系,则再入助推器的角速度可以用章动角δ、进动角ζ及自转角λ的角速度表示:
(1)
对上式进行求导并进行转换,可以得到以经典欧拉角表示的助推器非线性动力学模型:
(2)
式(2)中,
为角速度
的3个分量,可根据欧拉方程求得,矢量形式为
(3)
取状态变量为
则可得五维非线性动力学模型:
(4)
2 基于非线性稳定平衡点辨识运动模态
对于助推器再入飞行无控姿态动力学系统,为了得到系统在固定参数下的平衡点,需要将模型中的所以状态变量关于时间的导数设为0。以上标“P”表示平衡点处的各状态变量,则得到平衡点处的运动方程:
Mp-ωp×(Ip·ωp)=0
(5)
上式即为平衡状态的运动方程,可展开得到平衡方程的标量形式为
(6)
式中,
由于
可通过式(1)由经典欧拉角表示,且由式(4)可知平衡条件下
及
均为0,则平衡方程组(6)为包括3个未知数
的3个互相独立的方程,因此可以利用演化搜索等数值计算方法求得系统平衡点,结果见表1。
表1 利用数值方法求解得到特征Ma下对应非线性模型平衡点
Table 1 Numerical methods are used to solve the equilibrium points of the corresponding nonlinear model under the Ma
Ma序号ζ·p/radδP/radλp/rad0.21000202.344 280302.080 690.960 98941.729 251.682 890.099 5410.81000202.369 870302.110 260.948 94347.006 211.598 230.063 418
续表(表1)
Ma序号ζ·p/radδP/radλp/rad31000202.745 210302.588 560.825 71648.385 541.589 80.008 29871000202.841 740302.720 760.808 224425.928 11.594 410.016 204
由表1可见,助推器无控再入动力学模型存在4组平衡点。针对四组平衡点,可通过李雅普诺夫第一法对局部稳定性进行分析,以章动角δP为例得到不同组平衡点随马赫数变化的稳定性如图2所示。
图2 δ(rad)解曲线的稳定性
Fig.2 Stability of the curve to δ(rad)
由图2可以看出,第1组平衡点全速域下不稳定,第2、4组平衡点在Ma∈(0.76,0.84)域内可保持稳定,第3组平衡点在全速域下稳定。稳定的平衡点意味着在合适的初始条件下系统会被锁定在这一平衡状态下,形成一种固定的运动模态,4个平衡点即对应4种运动模态(见图3—图6),其中模态一、二、三对应进动角速度为零,模态四特定马赫数下进动角速度为一常值。
图3 运动模态一(不稳定)
Fig.3 Movement Mode 1
图4 运动模态二(静态)
Fig.4 Motion Mode 2 (Static)
图5 运动模态三(静态)
Fig.5 Motion Mode 3 (Static)
图6 运动模态四(动态)
Fig.6 Motion Mode 4 (Dynamic)
3 基于运动模态受力模型的落点预示
助推器再入过程中,进动角速度为0的模态可称为静模态,进动角速度不为0的模态可称为动模态,在解姿态动力学平衡方程时,未考虑进动角,因此对于静模态稳定状态而言,进动角可以为值域内任意值,并且保持不变,而动模态稳定状态的以角速度呈周期变化。
通过受力分析可知2类模态的质心运动受力模型存在差异,详见表2。
表2 不同模态受力特点分析
Table 2 Analysis of the stress characteristics
模态Maδ/radζ/radζ·/rad等效阻力等效升力静模态2～3任意固定≈0动模态≈1.57周期变化5～40×
其中,静模态对应常值章动角和常值进动角,根据图1可知在速度系下即表现为固定方向常值总攻角。而对于动模态而言,章动角短周期内不变,进动角随时间快速周期变化,在速度系下体现为方向周期变化的常值总攻角。
由于上述特点,静模态质心运动受到固定方向的升力和阻力。动模态稳定状态助推器保持一定进动速度做进动运动,升力在周期内可近似抵消,质心运动在长周期上只体现出阻力特性。已知助推器再入质心运动方程组如下所示:
(7)
式中:X为阻力;Y为升力;Z为侧力。
由姿态转化关系可知,每一模态确定一组δp、λp进而确定一组αp=tan-1(tanδpcosλp)、 βp=sin-1(sinδpsinλp)。同样总升力L以及阻力X仅与总攻角相关,总攻角即对应姿态平衡模态章动角δp。至此建立了质心运动模型和姿态运动模态之间的关系。
综上,基于模态预测的助推器质心动力学模型逻辑如图7所示。首先通过非线性动平衡理论对助推器进行稳定性分析,辨识稳定运动模态并提取各模态不同马赫数下对应经典欧拉角,求取不同欧拉角对应总升力、阻力,将稳定平衡状态总升力、阻力代入等效质心3DOF运动模型,结合分离点参数初值,计算得到不同模态对应运动包络轨迹,选取6条轨迹落点作为控制点绘制落点预示包络。
图7 基于模态预示的助推器落点计算流程
Fig.7 Booster landing point calculation process based on modal prediction
4 助推器无控再入落点范围预示结果仿真
本文中以某飞行器固体火箭助推器为例,代入上述落点计算模型和方程,计算落点范围,并与极限阻力系数法和6DOF打靶仿真方法进行对比分析。助推器本身特性及分离点参数详见表3。
表3 某助推器特征参数
Table 3 A characteristic parameter of a booster
名称符号数值结构参数转动惯量/(kg·m2)Ix、Iy、Iz300、9 000、9 000参考长度/mL8参考面积/m2S0.9分离点参数分离点海拔高度/kmH30初始俯仰角/(°)φ010初始偏航、滚转角/(°)ψ0、γ00姿态角速度/((°)·s-1φ0ψ0、λ·00速度倾角偏差范围/(°)Δη±10
通过本文中的方法,计算得到助推器的6个边界控制点坐标见表4。空间及投影曲线见图8、图9。
表4 某助推器再入边界控制点
Table 4 A booster re-entry boundary control point
序号X/kmZ/m1113.998586.6052114.406-4642.933114.3825 818.094198.792842.9575197.53910 232.96197.621-8 556.26
图8 基于模态预示的再入轨迹
Fig.8 Re-entry trajectory based on modal prediction
图9 基于模态预示的落点区域
Fig.9 Drop point area based on modal prediction
将本文中的方法与目前工程上常用的极限阻力系数法和6DOF打靶仿真方法进行对比结果如图10所示。
图10 3种助推器落区预示方法对比分析
Fig.10 Comparative analysis of three booster drop prediction methods
经分析,采用本文中提出的方法助推器落点预示范围仅为极限阻力系数法的16.9%,相较6DOF打靶法模态法设计安全控制区纵向上包络了98.51%的打靶工况,侧向上包络了99.26%的打靶工况,但计算量仅为打靶法的0.12%(打靶法进行5 000次弹道计算,模态预示法仅计算6条弹道),实现了高效率与准确性的结合。
同时利用本文中方法高效、准确、落区预示范围清晰的特点,在放开速度倾角约束的条件下,进一步针对某型助推器分析了其质心系数对于落区范围的影响,结果如图11所示。
图11 质心系数变化对落点预示范围的影响分析
Fig.11 Analysis of the influence of the change of centroid coefficient on the predictive range of the landing point
由图11的仿真结果可见,随着质心系数的减小,利用模态预示法估计得到的子级落点散布范围会逐渐减小,且估计区域将逐渐向发射点靠近。分析原因为,质心系数越小,平衡状态对应的平衡章动角越小,且由于平衡章动角均大于90°,小的章动角带来大的阻力系数,使得弹道曲率变大,落点更接近于发射点。同时由于大的阻力系数使得再入时间缩短,因此升力在各向上的作用时间也越短,使得弹道散布更小。因此在助推器设计中,可以设计较小的质心系数,使得弹道散布范围更小。
5 结论
本文中为实现助推器安全控制区的高效精细化设计,提出基于助推器无控再入稳定性分析的落点预示方法,相比于打靶法5 000条弹道的计算量,本文中方法仅需进行6条弹道计算即可实现纵向98.51%、侧向99.26%的落点预示覆盖,在确保一定预示精度的同时大幅提升了预示效率。同时,本文中方法能够实现对影响落区范围大小的助推器设计参数进行灵敏度分析,辨识出助推器质心系数对落点预示范围的重要影响,以某助推器为例,质心系数提高5%,落点预示范围可缩减约63%,为助推器安控区的正向设计提供了参考。
[1]陈彬,邓舞燕,高家一.运载火箭助推器再入姿态稳定性研究[J].导弹与航天运载技术,2015(3):13-15.CHEN Bin,DENG Wuyan,GAO Jiayi.Stability study of reentry attitude for carrier rocket boosters[J].Journal of Missile and Space Launch Technology,2015(3):13-15.
[2]李艳,康卓,刘溥.等.算法及其应用[J].武汉汽车工业大学学报,2000(3):101-104.LI Yan,KANG Zhuo,LIU Pu.et al.Algorithm and its applications[J].Journal of Wuhan University of Automotive Technology,2000(3):101-104.
[3]向占宏,向永红.基于变搜索区间的郭涛算法求解非线性方程组[J].科学技术与工程,2007(13):3117-3120.XIANG Zhanhong,XIANG Yonghong.Algorithm for solving nonlinear equation systems based on variable search interval[J].Science,Technology and Engineering,2007(13):3117-3120.
[4]陈传森,谢资清.非线性微分方程多解计算的搜索延拓法[M].北京:科学出版社,2005.CHEN Chuansen,XIE Ziqing.Search extension method for calculating multiple solutions of nonlinear differential equations[M].Beijing:Science Press,2005.
[5]王亚飞.再入飞行器动力学特性分析与控制方法研究[D].北京:北京理工大学,2015.WANG Yafei.Research on dynamic characteristics analysis and control methods of reentry vehicles[D].Beijing:Beijing Institute of Technology,2015.
[6]PLATUS D H.Ballistic reentry vehicle flight dynamics[J].Journal of Guidance,Control,and Dynamics,1982,5(1):4-16.
[7]穆山,蒲婷.航天发射场选址条件与场址勘选方法[J].飞行器测控学报,2011,30(1):11-15.MU Shan,PU Ting.Site selection conditions and site survey methods for aerospace launch sites[J].Journal of Aerospace Instrumentation,2011,30(1):11-15.
[8]肖裴,沈怀荣,李强.航天器发射场选址决策准则分析与建模[J].装备指挥技术学院学报,2003(2):56-60.XIAO Pei,SHEN Huairong,LI Qiang.Analysis and modeling of decision criteria for aerospace launch site selection[J].Journal of Equipment Command,Technology Institute,2003(2):56-60.
[9]王景国,卞韩城,陈学林,等.CZ-2F火箭整流罩残骸落点预报方法研究[J].载人航天,2014,20(5):457-460.WANG Jingguo,BIAN Hancheng,CHEN Xuelin,et al.Research on debris landing point prediction method of CZ-2F rocket fairing[J].Manned Spaceflight,2014,20(5):457-460.
[10]吕斐凯,贺卫亮.运载火箭助推器分离后的姿态和轨迹分析[J].导弹与航天运载技术,2015(1):13-16.LYU Feikai,HE Weiliang.Attitude and trajectory analysis of carrier rocket boosters after separation[J].Journal of Missile and Space Launch Technology,2015(1):13-16.
[11]洪金森,洪诗权.再入飞行器极限环运动分析[J].空气动力学学报,2005(2):204-209.HONG Jinsen,HONG Shiquan.Analysis of limit cycle motion of reentry vehicle[J].Acta Aerodynamica Sinica,2005(2):204-209.
[12]TARTABINI P,STARR B,GUMBERT C,et al.Ares IX separation and reentry trajectory analyses[C]//AIAA Atmospheric Flight Mechanics Conference.2011:6462.
[13]STARR B R,GUMBERT C R,TARTABINI P V.Ares IX test flight reference trajectory development[C]//58th JANNAF Propulsion Meeting.2011 (NF1676L-11608).
[14]PINIER J T.Ares I and Ares IX stage separation aerodynamic testing[J].Journal of Spacecraft and Rockets,2012,49(5):842-852.
[15]KING R,HENGEL J,WOLF D.Ares I first stage booster deceleration system-an overview[C]//20th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar.2009:2984.
[16]STARR B R,GOWAN JR J W,THOMPSON B G,et al.Ares IX range safety analyses overview[C]//58th JANNAF Propulsion Meeting.2011 (NF1676L-11849).
[17]KARLGAARD C D,BECK R E,STARR B R,et al.Ares IX best estimated trajectory analysis and results[C]//44th Combustion Meeting.2011 (NF1676L-11837).
[18]VISHNYAK A.Parachute simulations for ariane-5 booster recovery system design[C]//13th Aerodynamic Decelerator Systems Technology Conference.1995:1593.