一種極坐標系升力式再入六自由度仿真模型

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

(中國空氣動力研究與發展中心 計算空氣動力學研究所, 四川 綿陽 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