雍恩米, 汪清, 錢煒祺, 何開鋒

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力學(xué)研究所, 四川 綿陽(yáng) 621000)

一種極坐標(biāo)系升力式再入六自由度仿真模型

雍恩米, 汪清, 錢煒祺, 何開鋒

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力學(xué)研究所, 四川 綿陽(yáng) 621000)

建立了升力式再入飛行器在極坐標(biāo)系下的六自由度模型。首先,定義了建模所需坐標(biāo)系及轉(zhuǎn)換關(guān)系;然后,建立了極坐標(biāo)系下的升力式再入質(zhì)心運(yùn)動(dòng)方程以及面對(duì)稱飛行器在體坐標(biāo)系下的轉(zhuǎn)動(dòng)方程,并利用坐標(biāo)轉(zhuǎn)換關(guān)系建立歐拉角轉(zhuǎn)換輔助方程;最后,對(duì)歐拉角關(guān)系、開環(huán)六自由度進(jìn)行了仿真驗(yàn)證。仿真結(jié)果表明,該模型可為同時(shí)關(guān)注位置變量和速度變量的升力式再入飛行器的飛行力學(xué)問題研究,如六自由度的軌跡優(yōu)化、制導(dǎo)控制一體化設(shè)計(jì)等提供模型參考。

升力式再入; 六自由度; 飛行仿真; 建模; 極坐標(biāo)

0 引言

自航天飛機(jī)成功研制以來(lái),各種升力式再入飛行器成為各國(guó)的研究熱點(diǎn),如美國(guó)驗(yàn)證飛行器系列中的X-33,X-37以及高升阻比滑翔飛行器HTV-2,均為從軌道或亞軌道高度升力式再入飛行器。

從飛行力學(xué)的角度看,關(guān)于升力式再入飛行器的研究活躍在兩個(gè)方向,一是以質(zhì)點(diǎn)運(yùn)動(dòng)模型為基礎(chǔ)的再入軌跡優(yōu)化與制導(dǎo)技術(shù)研究,二是以包括轉(zhuǎn)動(dòng)運(yùn)動(dòng)的動(dòng)力學(xué)模型為基礎(chǔ)的飛行控制方法研究。前者一般建立極坐標(biāo)系下的三自由度模型[1-2],因?yàn)樵跇O坐標(biāo)系下建立質(zhì)心運(yùn)動(dòng)模型的優(yōu)勢(shì)是描述位置變量的方程相對(duì)簡(jiǎn)單,特別是考慮地球模型為旋轉(zhuǎn)圓球或橢球模型時(shí);后者一般在笛卡爾坐標(biāo)系(如速度系)描述質(zhì)心運(yùn)動(dòng),在體坐標(biāo)系描述繞質(zhì)心轉(zhuǎn)動(dòng)[3-4],且研究中主要關(guān)注的變量是與控制系統(tǒng)設(shè)計(jì)相關(guān)的速度和角速度變量,然而在這類笛卡爾坐標(biāo)系建立的運(yùn)動(dòng)模型,當(dāng)考慮地球球形與自轉(zhuǎn)模型時(shí),其位置變量方程會(huì)變得很復(fù)雜[5],有潛在的奇異性問題且難以無(wú)量綱化,因此一般在設(shè)計(jì)控制律時(shí)采用平面大地假設(shè)建立位置變量方程[6],然而平面大地假設(shè)對(duì)航程較遠(yuǎn)的升力式再入飛行器不再適用。

本文的研究目的是建立一種極坐標(biāo)系下的六自由度模型,并綜合上述兩種模型的優(yōu)勢(shì),為同時(shí)關(guān)注位置、速度及角速度的研究,如六自由度的軌跡優(yōu)化、制導(dǎo)控制一體化設(shè)計(jì)等提供模型基礎(chǔ)。

1 坐標(biāo)系定義及轉(zhuǎn)換關(guān)系

考慮到本文所用部分坐標(biāo)系與常規(guī)定義有所不同,因此對(duì)建模中所用坐標(biāo)系及轉(zhuǎn)換關(guān)系進(jìn)行說(shuō)明。

1.1 坐標(biāo)系定義

(1)地心球面固連坐標(biāo)系OEXYZ。該坐標(biāo)系原點(diǎn)位于地心OE,OEZ軸垂直于赤道平面指向北極,OEX和OEY在赤道平面內(nèi),其中OEX沿赤道平面與Greenwich子午面的相交線,OEY則由右手法則確定。該坐標(biāo)系用于確定極坐標(biāo)系描述的再入飛行器運(yùn)動(dòng)參數(shù)。

(2)彈體系Obxbybzb。彈體系原點(diǎn)Ob位于飛行器質(zhì)心,Obxb沿縱軸指向頭部,Obyb在飛行器對(duì)稱面內(nèi)垂直于Obxb,Obzb則與Obxb和Obyb構(gòu)成右手坐標(biāo)系。該坐標(biāo)系與常規(guī)體坐標(biāo)系定義相同。

(3)地心系OExEyEzE。定義地心固連坐標(biāo)系原點(diǎn)OE位于地心,OExE垂直于赤道平面指向北極,OEyE沿赤道平面與經(jīng)度為λ0的子午面的交線,OEzE與OExE和OEyE構(gòu)成右手坐標(biāo)系。該坐標(biāo)系與常規(guī)地心系定義略有不同,即OEyE指向再入點(diǎn)經(jīng)度而非零經(jīng)度,該定義是為了簡(jiǎn)化建模中的坐標(biāo)轉(zhuǎn)換。

(4)再入系Oexeyeze。再入系原點(diǎn)為再入點(diǎn)星下點(diǎn)Oe,Oexe在水平面內(nèi)指向目標(biāo)點(diǎn),Oeye在水平面內(nèi)與Oexe垂直,Oeze與Oexe和Oeye構(gòu)成右手坐標(biāo)系。該坐標(biāo)系即再入動(dòng)力學(xué)中通常采用再入坐標(biāo)系。

(5)當(dāng)?shù)乇碧鞏|坐標(biāo)系Odxdydzd。該坐標(biāo)原點(diǎn)位于飛行器質(zhì)心或其星下點(diǎn)Od,Odxd指向正北方向,Odyd與Odxd垂直并指向天,Odzd與Odxd和Odxd構(gòu)成右手坐標(biāo)系。該坐標(biāo)系定義與一般北天東坐標(biāo)系定義略有區(qū)別,該定義是為了減少建模中的坐標(biāo)轉(zhuǎn)換次數(shù)。

(6)速度坐標(biāo)系O1xvyvzv。速度系原點(diǎn)位于再入飛行器質(zhì)心O1,O1xv沿飛行器速度方向,O1xv在飛行器主對(duì)稱面內(nèi)垂直于O1xv,O1zv與O1xv和O1xv構(gòu)成右手坐標(biāo)系。

(7)半速度坐標(biāo)系O1xhyhzh。該坐標(biāo)系原點(diǎn)位于飛行器質(zhì)心O1,O1xh與O1xv方向重合,O1yh在再入坐標(biāo)系的xeOye平面內(nèi)垂直于O1xh軸,O1zh與O1xh和O1yh構(gòu)成右手坐標(biāo)系。

1.2 坐標(biāo)轉(zhuǎn)換關(guān)系

根據(jù)上述坐標(biāo)系定義,可以得出以下坐標(biāo)系轉(zhuǎn)換關(guān)系,其中所用再入點(diǎn)參數(shù)分別為再入點(diǎn)星下點(diǎn)經(jīng)度λ0、緯度μ0和再入方位角A0。

(1)當(dāng)?shù)乇碧鞏|坐標(biāo)系Odxdydzd與地心系OExEyEzE的轉(zhuǎn)換矩陣:

(1)

式中:Δλ=λ0-λ;λ,μ分別為再入飛行器當(dāng)前位置星下點(diǎn)經(jīng)度和緯度。

(2)地心系OExEyEzE與再入系Oexeyeze的轉(zhuǎn)換矩陣:

(2)

(3)再入系Oexeyeze與彈體系Obxbybzb的轉(zhuǎn)換矩陣?？紤]升力式再入飛行器有較大升阻比,且多采用側(cè)傾轉(zhuǎn)彎的控制方式,姿態(tài)歐拉角按“2-3-1”的順序定義,即按偏航ψ、滾轉(zhuǎn)φ和俯仰θ順序進(jìn)行坐標(biāo)轉(zhuǎn)換,因此從再入系到彈體系的轉(zhuǎn)換矩陣為:

φ)My(ψ)

(3)

以下公式中繞各坐標(biāo)軸的基元變換矩陣同上,進(jìn)行簡(jiǎn)化書寫。

(4)再入系到半速度系的轉(zhuǎn)換矩陣:

(4)

式中:Θ為速度傾角,是速度方向在Oxeye平面上的投影與Oxe間的夾角;υ為航跡偏航角,是速度方向與Oexeye平面的夾角。

(5)速度系到彈體系的轉(zhuǎn)換矩陣:

(5)

式中：α,β分別為迎角和側(cè)滑角。

(6)半速度系到速度系的轉(zhuǎn)換矩陣:

(6)

式中:σ為側(cè)傾角。

上述坐標(biāo)系轉(zhuǎn)換關(guān)系如圖1所示。

圖1 坐標(biāo)轉(zhuǎn)換關(guān)系Fig.1 Relationship of coordinate transformation

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

(1)極坐標(biāo)系下建立再入質(zhì)心運(yùn)動(dòng)方程

在地心球面固連坐標(biāo)系下,基于極坐標(biāo)建立升力式再入飛行器質(zhì)心運(yùn)動(dòng)方程,其中位置參數(shù)包括地心距r、經(jīng)度λ和緯度μ,速度參數(shù)包括速度大小V、航跡角θv和航向角ψv,如圖2所示。圖中：θv為速度向量與當(dāng)?shù)厮矫娴膴A角,向上為正;ψv為速度向量在當(dāng)?shù)厮矫嫱队芭c正北方向的夾角,順時(shí)針旋轉(zhuǎn)為正。

圖2 升力式再入運(yùn)動(dòng)參數(shù)示意圖Fig.2 Kinetic parameters of lifting reentry

考慮地球?yàn)樾D(zhuǎn)圓球的再入運(yùn)動(dòng)方程[7]為:

(7)

(8)

式中:ω,g,β,L,D,C分別為地球自轉(zhuǎn)角速度、引力加速度、側(cè)滑角、升力、阻力和側(cè)向力。

氣動(dòng)力計(jì)算表達(dá)式為:

(9)

式中:q=ρV2/2為動(dòng)壓;CL,CD,CC,Sref分別為升力系數(shù)、阻力系數(shù)、側(cè)向力系數(shù)和氣動(dòng)參考面積。

(2)體坐標(biāo)系下建立轉(zhuǎn)動(dòng)動(dòng)力學(xué)方程

設(shè)體坐標(biāo)系相對(duì)再入系的角速度矢量為ω=[ωx,ωy,ωz]T,面對(duì)稱再入飛行器的轉(zhuǎn)動(dòng)慣量為:

(10)

根據(jù)動(dòng)量矩定理建立轉(zhuǎn)動(dòng)方程:

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

(11)

根據(jù)式(12)有:

(12)

分解得到標(biāo)量形式的轉(zhuǎn)動(dòng)方程為:

(13)

根據(jù)坐標(biāo)轉(zhuǎn)換關(guān)系,還可建立姿態(tài)運(yùn)動(dòng)學(xué)方程[6]:

(14)

氣動(dòng)力矩計(jì)算表達(dá)式為:

(15)

(3)輔助方程

氣動(dòng)升力L、阻力D以及氣動(dòng)力矩(Mx1,My1,Mz1)是飛行狀態(tài)以及迎角α、側(cè)滑角β的函數(shù)。另外,采用傾側(cè)轉(zhuǎn)彎控制方式的升力式再入飛行器,制導(dǎo)方程一般設(shè)計(jì)側(cè)傾角σ,因此需要建立動(dòng)力學(xué)方程所涉及的歐拉角的轉(zhuǎn)換關(guān)系,作為六自由度仿真模型的輔助方程。已包含于質(zhì)心運(yùn)動(dòng)方程及繞質(zhì)心轉(zhuǎn)動(dòng)方程狀態(tài)變量中的歐拉角包括:航跡角θv、航向角ψv、俯仰角θ、偏航角ψ和滾轉(zhuǎn)角φ,待確定的歐拉角包括迎角α、側(cè)滑角β和側(cè)傾角σ。

根據(jù)定義,迎角和側(cè)滑角可由體坐標(biāo)系的速度分量求解,即:

(16)

為獲得計(jì)算迎角和側(cè)滑角所需的體坐標(biāo)系的速度分量,首先,根據(jù)描述質(zhì)心運(yùn)動(dòng)方程的變量,即速度大小V、航跡角θv和航向角ψv,計(jì)算出本文定義的當(dāng)?shù)乇碧鞏|坐標(biāo)系中的速度分量,即:

(17)

根據(jù)圖1所示的坐標(biāo)轉(zhuǎn)換關(guān)系,可由當(dāng)?shù)乇睎|天坐標(biāo)系速度得到體坐標(biāo)系下的速度,即:

(18)

結(jié)合式(16)和式(18),則可利用極坐標(biāo)描述的運(yùn)動(dòng)參數(shù)求解迎角α和側(cè)滑角β。

為求解側(cè)傾角,還需用到兩個(gè)歐拉角,即速度傾角Θ和航跡偏航角υ。根據(jù)定義,速度傾角和航跡偏航角可由再入系的速度分量求解,即:

(19)

而再入系的速度也可由當(dāng)?shù)乇碧鞏|坐標(biāo)系速度求得,即:

(20)

從而,利用半速度系到速度系的轉(zhuǎn)換矩陣可求解側(cè)傾角,即:

(21)

(22)

根據(jù)圖1所示的坐標(biāo)轉(zhuǎn)換關(guān)系,從半速度系到速度系轉(zhuǎn)換矩陣還可以通過再入系、彈體系再到速度系進(jìn)行轉(zhuǎn)換,即:

(23)

因此,結(jié)合式(16)、式(19)、式(20)、式(22)和式(23)可根據(jù)極坐標(biāo)系下的運(yùn)動(dòng)參數(shù)求出側(cè)傾角σ。

3 仿真驗(yàn)證及結(jié)果分析

3.1 歐拉角轉(zhuǎn)換關(guān)系驗(yàn)證

采用文獻(xiàn)[1]關(guān)于歐拉角的轉(zhuǎn)換關(guān)系表達(dá)式對(duì)上述模型進(jìn)行局部驗(yàn)證,利用迎角、側(cè)滑角以及姿態(tài)角求解側(cè)傾角的表達(dá)式為:

sinσ=[ cosαsinβsinθ-

sinαsinβcosφcosθ+

cosβsinφcosθ]/cosΘ

(24)

表1為三組算例計(jì)算狀態(tài),表2為計(jì)算所得歐拉角,其中α,β,σ為本文的輔助方程模型計(jì)算的歐拉角,σ′為根據(jù)式(24)計(jì)算的側(cè)傾角?？梢钥闯?兩組幾何方程計(jì)算出的側(cè)傾角一致。

表1 計(jì)算狀態(tài)Table 1 States of the computation

表2 歐拉角計(jì)算結(jié)果Table 2 Results of Euler angles

3.2 開環(huán)六自由度仿真驗(yàn)證

采用X-33的基本數(shù)據(jù)進(jìn)行開環(huán)六自由度仿真計(jì)算?？傮w參數(shù)為:質(zhì)量為35 828 kg;參考面積為149.39 m2;參考長(zhǎng)度為10.97 m;轉(zhuǎn)動(dòng)慣量Ixx=588 790 kg·m2,Iyy=1 534 200 kg·m2,Izz=1.303 200 kg·m2,Ixy=24 242 kg·m2。氣動(dòng)數(shù)據(jù)采用擬合公式[7]擬合為迎角、側(cè)滑角和舵偏角的函數(shù)?？紤]X-33從亞軌道高度再入的情況,取初始高度為50 km,初始速度為3 680 m/s,初始航跡角為0°,舵偏角均為零,仿真結(jié)果如圖3～圖5所示。該算例僅為開環(huán)六自由度仿真模型計(jì)算,即再入處于無(wú)控狀態(tài),因此姿態(tài)角呈發(fā)散狀態(tài),可以在此基礎(chǔ)上增加制導(dǎo)控制模型，形成閉環(huán)六自由度仿真模型。

圖4 馬赫數(shù)-航程變化曲線Fig.4 Curve of Ma vs R

圖5 姿態(tài)角變化曲線Fig.5 Variation curves of altitude angle

4 結(jié)束語(yǔ)

本文建立了一種極坐標(biāo)下的升力式再入六自由度仿真模型。極坐標(biāo)系下建立的考慮地球?yàn)樾D(zhuǎn)圓球的質(zhì)心運(yùn)動(dòng)方程形式相對(duì)簡(jiǎn)單,便于升力式再入飛行器的軌跡優(yōu)化與制導(dǎo)方法研究。極坐標(biāo)系下的

運(yùn)動(dòng)方程和體系坐標(biāo)系下的轉(zhuǎn)動(dòng)方程的耦合關(guān)系源于隱含于氣動(dòng)力和氣動(dòng)力矩系數(shù)的迎角和側(cè)滑角。通過坐標(biāo)系轉(zhuǎn)換,可以獲得描述歐拉角轉(zhuǎn)換關(guān)系的輔助方程。下一步工作是開展基于該六自由度模型的制導(dǎo)控制方法研究,建立升力式再入飛行器閉環(huán)六自由度仿真模型。

[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] 李菁菁,任章,宋劍爽.高超聲速再入滑翔飛行器的模糊變結(jié)構(gòu)控制[J].上海交通大學(xué)學(xué)報(bào),2011,45(2):295-300.

[4] 葉振信,李傳峰,傅維賢,等.大氣層內(nèi)無(wú)動(dòng)力滑翔炸彈傾斜轉(zhuǎn)彎控制設(shè)計(jì)[J].導(dǎo)彈與航天運(yùn)載技術(shù),2010(6):6-9.

[5] 賈沛然,陳克俊,何力.遠(yuǎn)程火箭彈道學(xué)[M].長(zhǎng)沙:國(guó)防科學(xué)技術(shù)大學(xué)出版社,1993:57-67.

[6] 錢杏芳,林瑞雄,趙亞男.導(dǎo)彈飛行力學(xué)[M].北京:北京理工大學(xué)出版社,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;

時(shí)間:2014-11-04 08:28

雍恩米(1979-)，女，四川綿陽(yáng)人，博士，研究方向?yàn)轱w行器軌跡優(yōu)化與制導(dǎo)控制、飛行建模與仿真。

V412

A

1002-0853(2015)01-0052-05