一种极坐标系升力式再入六自由度仿真模型

雍恩米, 汪清, 钱炜祺, 何开锋

(中国空气动力研究与发展中心 计算空气动力学研究所, 四川 绵阳 621000)

一种极坐标系升力式再入六自由度仿真模型

雍恩米, 汪清, 钱炜祺, 何开锋

(中国空气动力研究与发展中心 计算空气动力学研究所, 四川 绵阳 621000)

建立了升力式再入飞行器在极坐标系下的六自由度模型。首先,定义了建模所需坐标系及转换关系;然后,建立了极坐标系下的升力式再入质心运动方程以及面对称飞行器在体坐标系下的转动方程,并利用坐标转换关系建立欧拉角转换辅助方程;最后,对欧拉角关系、开环六自由度进行了仿真验证。仿真结果表明,该模型可为同时关注位置变量和速度变量的升力式再入飞行器的飞行力学问题研究,如六自由度的轨迹优化、制导控制一体化设计等提供模型参考。

升力式再入; 六自由度; 飞行仿真; 建模; 极坐标

0 引言

自航天飞机成功研制以来,各种升力式再入飞行器成为各国的研究热点,如美国验证飞行器系列中的X-33,X-37以及高升阻比滑翔飞行器HTV-2,均为从轨道或亚轨道高度升力式再入飞行器。

从飞行力学的角度看,关于升力式再入飞行器的研究活跃在两个方向,一是以质点运动模型为基础的再入轨迹优化与制导技术研究,二是以包括转动运动的动力学模型为基础的飞行控制方法研究。前者一般建立极坐标系下的三自由度模型[1-2],因为在极坐标系下建立质心运动模型的优势是描述位置变量的方程相对简单,特别是考虑地球模型为旋转圆球或椭球模型时;后者一般在笛卡尔坐标系(如速度系)描述质心运动,在体坐标系描述绕质心转动[3-4],且研究中主要关注的变量是与控制系统设计相关的速度和角速度变量,然而在这类笛卡尔坐标系建立的运动模型,当考虑地球球形与自转模型时,其位置变量方程会变得很复杂[5],有潜在的奇异性问题且难以无量纲化,因此一般在设计控制律时采用平面大地假设建立位置变量方程[6],然而平面大地假设对航程较远的升力式再入飞行器不再适用。

本文的研究目的是建立一种极坐标系下的六自由度模型,并综合上述两种模型的优势,为同时关注位置、速度及角速度的研究,如六自由度的轨迹优化、制导控制一体化设计等提供模型基础。

1 坐标系定义及转换关系

考虑到本文所用部分坐标系与常规定义有所不同,因此对建模中所用坐标系及转换关系进行说明。

1.1 坐标系定义

(1)地心球面固连坐标系OEXYZ。该坐标系原点位于地心OE,OEZ轴垂直于赤道平面指向北极,OEX和OEY在赤道平面内,其中OEX沿赤道平面与Greenwich子午面的相交线,OEY则由右手法则确定。该坐标系用于确定极坐标系描述的再入飞行器运动参数。

(2)弹体系Obxbybzb。弹体系原点Ob位于飞行器质心,Obxb沿纵轴指向头部,Obyb在飞行器对称面内垂直于Obxb,Obzb则与Obxb和Obyb构成右手坐标系。该坐标系与常规体坐标系定义相同。

(3)地心系OExEyEzE。定义地心固连坐标系原点OE位于地心,OExE垂直于赤道平面指向北极,OEyE沿赤道平面与经度为λ0的子午面的交线,OEzE与OExE和OEyE构成右手坐标系。该坐标系与常规地心系定义略有不同,即OEyE指向再入点经度而非零经度,该定义是为了简化建模中的坐标转换。

(4)再入系Oexeyeze。再入系原点为再入点星下点Oe,Oexe在水平面内指向目标点,Oeye在水平面内与Oexe垂直,Oeze与Oexe和Oeye构成右手坐标系。该坐标系即再入动力学中通常采用再入坐标系。

(5)当地北天东坐标系Odxdydzd。该坐标原点位于飞行器质心或其星下点Od,Odxd指向正北方向,Odyd与Odxd垂直并指向天,Odzd与Odxd和Odxd构成右手坐标系。该坐标系定义与一般北天东坐标系定义略有区别,该定义是为了减少建模中的坐标转换次数。

(6)速度坐标系O1xvyvzv。速度系原点位于再入飞行器质心O1,O1xv沿飞行器速度方向,O1xv在飞行器主对称面内垂直于O1xv,O1zv与O1xv和O1xv构成右手坐标系。

(7)半速度坐标系O1xhyhzh。该坐标系原点位于飞行器质心O1,O1xh与O1xv方向重合,O1yh在再入坐标系的xeOye平面内垂直于O1xh轴,O1zh与O1xh和O1yh构成右手坐标系。

1.2 坐标转换关系

根据上述坐标系定义,可以得出以下坐标系转换关系,其中所用再入点参数分别为再入点星下点经度λ0、纬度μ0和再入方位角A0。

(1)当地北天东坐标系Odxdydzd与地心系OExEyEzE的转换矩阵:

(1)

式中:Δλ=λ0-λ;λ,μ分别为再入飞行器当前位置星下点经度和纬度。

(2)地心系OExEyEzE与再入系Oexeyeze的转换矩阵:

(2)

(3)再入系Oexeyeze与弹体系Obxbybzb的转换矩阵。考虑升力式再入飞行器有较大升阻比,且多采用侧倾转弯的控制方式,姿态欧拉角按“2-3-1”的顺序定义,即按偏航ψ、滚转φ和俯仰θ顺序进行坐标转换,因此从再入系到弹体系的转换矩阵为:

φ)My(ψ)

(3)

以下公式中绕各坐标轴的基元变换矩阵同上,进行简化书写。

(4)再入系到半速度系的转换矩阵:

(4)

式中:Θ为速度倾角,是速度方向在Oxeye平面上的投影与Oxe间的夹角;υ为航迹偏航角,是速度方向与Oexeye平面的夹角。

(5)速度系到弹体系的转换矩阵:

(5)

式中：α,β分别为迎角和侧滑角。

(6)半速度系到速度系的转换矩阵:

(6)

式中:σ为侧倾角。

上述坐标系转换关系如图1所示。

图1 坐标转换关系Fig.1 Relationship of coordinate transformation

2 六自由度仿真模型及辅助方程

(1)极坐标系下建立再入质心运动方程

在地心球面固连坐标系下,基于极坐标建立升力式再入飞行器质心运动方程,其中位置参数包括地心距r、经度λ和纬度μ,速度参数包括速度大小V、航迹角θv和航向角ψv,如图2所示。图中：θv为速度向量与当地水平面的夹角,向上为正;ψv为速度向量在当地水平面投影与正北方向的夹角,顺时针旋转为正。

图2 升力式再入运动参数示意图Fig.2 Kinetic parameters of lifting reentry

考虑地球为旋转圆球的再入运动方程[7]为:

(7)

(8)

式中:ω,g,β,L,D,C分别为地球自转角速度、引力加速度、侧滑角、升力、阻力和侧向力。

气动力计算表达式为:

(9)

式中:q=ρV2/2为动压;CL,CD,CC,Sref分别为升力系数、阻力系数、侧向力系数和气动参考面积。

(2)体坐标系下建立转动动力学方程

设体坐标系相对再入系的角速度矢量为ω=[ωx,ωy,ωz]T,面对称再入飞行器的转动惯量为:

(10)

根据动量矩定理建立转动方程:

Idω/dt+ω×(I·ω)=M

(11)

根据式(12)有:

(12)

分解得到标量形式的转动方程为:

(13)

根据坐标转换关系,还可建立姿态运动学方程[6]:

(14)

气动力矩计算表达式为:

(15)

(3)辅助方程

气动升力L、阻力D以及气动力矩(Mx1,My1,Mz1)是飞行状态以及迎角α、侧滑角β的函数。另外,采用倾侧转弯控制方式的升力式再入飞行器,制导方程一般设计侧倾角σ,因此需要建立动力学方程所涉及的欧拉角的转换关系,作为六自由度仿真模型的辅助方程。已包含于质心运动方程及绕质心转动方程状态变量中的欧拉角包括:航迹角θv、航向角ψv、俯仰角θ、偏航角ψ和滚转角φ,待确定的欧拉角包括迎角α、侧滑角β和侧倾角σ。

根据定义,迎角和侧滑角可由体坐标系的速度分量求解,即:

(16)

为获得计算迎角和侧滑角所需的体坐标系的速度分量,首先,根据描述质心运动方程的变量,即速度大小V、航迹角θv和航向角ψv,计算出本文定义的当地北天东坐标系中的速度分量,即:

(17)

根据图1所示的坐标转换关系,可由当地北东天坐标系速度得到体坐标系下的速度,即:

(18)

结合式(16)和式(18),则可利用极坐标描述的运动参数求解迎角α和侧滑角β。

为求解侧倾角,还需用到两个欧拉角,即速度倾角Θ和航迹偏航角υ。根据定义,速度倾角和航迹偏航角可由再入系的速度分量求解,即:

(19)

而再入系的速度也可由当地北天东坐标系速度求得,即:

(20)

从而,利用半速度系到速度系的转换矩阵可求解侧倾角,即:

(21)

(22)

根据图1所示的坐标转换关系,从半速度系到速度系转换矩阵还可以通过再入系、弹体系再到速度系进行转换,即:

(23)

因此,结合式(16)、式(19)、式(20)、式(22)和式(23)可根据极坐标系下的运动参数求出侧倾角σ。

3 仿真验证及结果分析

3.1 欧拉角转换关系验证

采用文献[1]关于欧拉角的转换关系表达式对上述模型进行局部验证,利用迎角、侧滑角以及姿态角求解侧倾角的表达式为:

sinσ=[ cosαsinβsinθ-

sinαsinβcosφcosθ+

cosβsinφcosθ]/cosΘ

(24)

表1为三组算例计算状态,表2为计算所得欧拉角,其中α,β,σ为本文的辅助方程模型计算的欧拉角,σ′为根据式(24)计算的侧倾角。可以看出,两组几何方程计算出的侧倾角一致。

表1 计算状态Table 1 States of the computation

表2 欧拉角计算结果Table 2 Results of Euler angles

3.2 开环六自由度仿真验证

采用X-33的基本数据进行开环六自由度仿真计算。总体参数为:质量为35 828 kg;参考面积为149.39 m2;参考长度为10.97 m;转动惯量Ixx=588 790 kg·m2,Iyy=1 534 200 kg·m2,Izz=1.303 200 kg·m2,Ixy=24 242 kg·m2。气动数据采用拟合公式[7]拟合为迎角、侧滑角和舵偏角的函数。考虑X-33从亚轨道高度再入的情况,取初始高度为50 km,初始速度为3 680 m/s,初始航迹角为0°,舵偏角均为零,仿真结果如图3～图5所示。该算例仅为开环六自由度仿真模型计算,即再入处于无控状态,因此姿态角呈发散状态,可以在此基础上增加制导控制模型，形成闭环六自由度仿真模型。

图4 马赫数-航程变化曲线Fig.4 Curve of Ma vs R

图5 姿态角变化曲线Fig.5 Variation curves of altitude angle

4 结束语

本文建立了一种极坐标下的升力式再入六自由度仿真模型。极坐标系下建立的考虑地球为旋转圆球的质心运动方程形式相对简单,便于升力式再入飞行器的轨迹优化与制导方法研究。极坐标系下的

运动方程和体系坐标系下的转动方程的耦合关系源于隐含于气动力和气动力矩系数的迎角和侧滑角。通过坐标系转换,可以获得描述欧拉角转换关系的辅助方程。下一步工作是开展基于该六自由度模型的制导控制方法研究,建立升力式再入飞行器闭环六自由度仿真模型。

[1] Harpold J C,Jr Graves C A.Shuttle entry guidance[J].The Journal of the Astronautical Sciences,1979,27(3):239-268.

[2] Shen Zuojun,Lu Ping.Onboard generation of three-dimensional constrained entry trajectories[J].Journal of Guidance,Control,and Dynamics,2003,26(1):110-121.

[3] 李菁菁,任章,宋剑爽.高超声速再入滑翔飞行器的模糊变结构控制[J].上海交通大学学报,2011,45(2):295-300.

[4] 叶振信,李传峰,傅维贤,等.大气层内无动力滑翔炸弹倾斜转弯控制设计[J].导弹与航天运载技术,2010(6):6-9.

[5] 贾沛然,陈克俊,何力.远程火箭弹道学[M].长沙:国防科学技术大学出版社,1993:57-67.

[6] 钱杏芳,林瑞雄,赵亚男.导弹飞行力学[M].北京:北京理工大学出版社,2008:36-48.

[7] Bollino K P.High-fidelity real-time trajectory optimization for reusable launch vehicles[D].California:Naval Postgraduate School,2006.

(编辑：李怡)

Six DOF simulation model in polar coordinate for lifting reentry vehicles

YONG En-mi, WANG Qing, QIAN Wei-qi, HE Kai-feng

(Computational Aerodynamics Institute, CARDC, Mianyang 621000, China)

A six degree-of-freedom (DOF) simulation model in polar coordinate for lifting reentry vehicles is presented. Firstly, some reference frames used in this model were defined. Secondly, the kinetic equations in polar coordinate and the equation in body axis coordinate were both constructed. Then the auxiliary equations for Euler angles were built by coordinate transformation. Finally, the transformation relationship between Euler angles was verified and an open loop six DOF reentry flight simulation was implemented. Simulation results show that the reentry dynamic model is potential to be used in the reentry flight dynamic research which is focused on both of the position and velocity variables, such as six DOF trajectory optimization and integrated guidance and control problems.

lifting reentry; six DOF; fight simulation; modeling; polar coordinate

2014-06-16;

2014-09-18;

时间:2014-11-04 08:28

雍恩米(1979-)，女，四川绵阳人，博士，研究方向为飞行器轨迹优化与制导控制、飞行建模与仿真。

V412

A

1002-0853(2015)01-0052-05