徐敏, 白斌, 祝小平

(1.西北工業(yè)大學(xué) 航天學(xué)院, 陜西 西安 710072;2.西北工業(yè)大學(xué) 無人機(jī)研究所, 陜西 西安 710072)

0 引言

目前大部分高超聲速飛行器都具有典型的乘波體機(jī)身和矩形進(jìn)氣道超燃沖壓發(fā)動(dòng)機(jī)一體化設(shè)計(jì)外形[1]。由于吸氣式高超聲速飛行器建模涉及到彈性、氣動(dòng)和推進(jìn)的耦合作用,屬于多學(xué)科交叉問題,建模難度較大,傳統(tǒng)的飛行動(dòng)力學(xué)模型將不再適用。文獻(xiàn)[2-4]分別建立了簡(jiǎn)化的二維高超聲速飛行器模型,并進(jìn)行了相應(yīng)研究。

在30多年的探索中,作為非線性動(dòng)力學(xué)特性分析的重要工具,分支分析的理論發(fā)展經(jīng)歷了兩個(gè)重要的階段。1977 年,Mehra等[5]首次提出分支分析與突變理論方法(簡(jiǎn)稱常規(guī)分支分析方法),并對(duì)飛機(jī)大迎角快速滾轉(zhuǎn)的操縱性和穩(wěn)定性進(jìn)行了非線性動(dòng)力學(xué)分析。2001年,Ananthkrishnan和 Sinha提出擴(kuò)展分支分析方法[6],利用此方法可以有效地分析約束飛行狀態(tài)下的非線性動(dòng)力學(xué)特性,從而彌補(bǔ)常規(guī)分支分析方法只能連續(xù)變化單一控制參數(shù)達(dá)不到約束飛行條件這一缺陷。

本文對(duì)Bolender和Doman提出的二維高超聲速飛行器模型[3]進(jìn)行縱向分析。首先利用空氣動(dòng)力學(xué)相關(guān)理論得到氣動(dòng)力和推力數(shù)據(jù),進(jìn)而完成彈性模型計(jì)算,最后得到了巡航狀態(tài)下飛行器的分支圖,并針對(duì)平衡曲線對(duì)系統(tǒng)穩(wěn)定性和機(jī)身彈性的影響進(jìn)行了簡(jiǎn)要分析。

1 幾何外形描述

如文獻(xiàn)[3],縱向剖面視圖及機(jī)體坐標(biāo)系Oxyz如圖1所示。

圖1 吸氣式高超聲速飛行器幾何外形Fig.1 Geometry of hypersonic air-breathing vehicle

飛行器在 30 km的高度以Ma0=10的速度進(jìn)行巡航運(yùn)動(dòng)。聲速計(jì)算方法如下:

(1)

式中:γ為比熱比,這里為固定值1.4;Ra為空氣氣體常數(shù)287.0041 J/(kg·K)。利用式(1)可以方便地計(jì)算出飛行速度V0。

2 氣動(dòng)計(jì)算

自由來流流過機(jī)體前部下表面ad時(shí),與下表面的夾角為δ=τ1,l+α。當(dāng)δ>0時(shí),采用Oblique-Shock理論來計(jì)算激波后的馬赫數(shù)Ma1、溫度T1和壓力P1;當(dāng)δ<0時(shí),則采用Prandtl-Meyer理論來計(jì)算膨脹波后的氣體參數(shù)。

當(dāng)δ>0時(shí),氣流分布如圖2所示。為了實(shí)現(xiàn)空氣流量的最大化,發(fā)動(dòng)機(jī)的進(jìn)氣道唇口需設(shè)計(jì)為可變結(jié)構(gòu),使得在不同迎角和馬赫數(shù)下,機(jī)體前端的激波經(jīng)過唇口前端后產(chǎn)生的激波剛好打在進(jìn)氣口上端前部,以保證進(jìn)入發(fā)動(dòng)機(jī)的氣流為一維流。

圖2 機(jī)體下表面氣流分布圖Fig.2 Air flow distribution at the lower surface of the vehicle

同理,利用Oblique-Shock理論或Prandtl-Meyer理論可以分別計(jì)算出機(jī)體上表面ab的壓力Pf、發(fā)動(dòng)機(jī)下表面ef的壓力Pn、機(jī)體尾部下表面bc的壓力Pa、升降舵δe的壓力Pδe和前端鴨舵δc的壓力Pδc,進(jìn)而計(jì)算相應(yīng)的氣動(dòng)力,詳見文獻(xiàn)[3]。各個(gè)初始參數(shù)取值范圍如表1所示。

表1 各參數(shù)取值范圍Table 1 Range of the parameters

表1中的ψ為燃油當(dāng)量比。當(dāng)Ma0=10,ψ=0.4,δc=0°時(shí),得到升力、阻力和俯仰力矩隨迎角、升降舵偏角變化的三維圖如圖3所示。

圖3 升力、阻力和俯仰力矩隨迎角、升降舵偏角變化曲線Fig.3 Changes of lift, drag and pitch moment with AOA and elevator deflection

對(duì)比文獻(xiàn)[7]相應(yīng)的結(jié)果,可以發(fā)現(xiàn)兩者在走勢(shì)上是一致的。出現(xiàn)差異的主要原因是本文在機(jī)體前部加入了鴨舵δc,同時(shí)飛行器的幾何參數(shù)也有一定的不同。

3 推力計(jì)算

發(fā)動(dòng)機(jī)模型如圖4所示,分為以下3個(gè)階段:A1～A2為氣體壓縮減速階段;A2～A3為燃燒室;A3～Ae為氣體膨脹加速階段。

圖4 發(fā)動(dòng)機(jī)模型Fig.4 Engine model

在A1～A2段,利用連續(xù)方程(質(zhì)量守恒)可以計(jì)算出燃燒室入口A2截面處的氣體參數(shù)Ma3,T3,P3。

A2～A3燃燒室內(nèi),為避免碳?xì)淙剂先紵耶a(chǎn)生復(fù)雜的物理化學(xué)變化給推力的計(jì)算帶來困難,采用氫氣為燃料,根據(jù)一維Rayleigh流理論可知:

(2)

(3)

(4)

(5)

式中:qin為燃燒室內(nèi)增加的熱量;cp為氫氣定壓比熱容;ηc為燃燒效率;LHV為氫氣低發(fā)熱值;fst為氫氣理論化學(xué)當(dāng)量比。

為簡(jiǎn)化起見,本文引用文獻(xiàn)[8]給出的燃燒效率與燃油當(dāng)量比的簡(jiǎn)單函數(shù)關(guān)系進(jìn)行計(jì)算。T04,T03為總溫,P04,P03為總壓,存在以下關(guān)系:

(6)

(7)

在A3～Ae段,同樣利用連續(xù)方程可以得到發(fā)動(dòng)機(jī)出口處的氣體參數(shù)Mae,Te,Pe,則發(fā)動(dòng)機(jī)的推力為:

(8)

當(dāng)Ma0=10時(shí),得到推力隨α和ψ的變化曲線如圖5所示。采用擬合方法即可得到推力關(guān)于Ma0,α,ψ的解析表達(dá)式。

圖5 推力隨α和ψ的變化曲線Fig.5 Change of thrust with α and ψ

4 彈性模型

如圖6所示,可以將飛行器看成兩個(gè)固定于重心處的懸臂梁模型,忽略扭轉(zhuǎn)變形和沿x軸方向的伸縮變形,只考慮彎曲變形,由于變形位移較小,故滿足胡克定律。

圖6 HSV彈性模型Fig.6 HSV aeroelastic model

彎曲振動(dòng)微分方程為:

(9)

φf,k(x)=

(10)

φa,k(x)=

(11)

(12)

(13)

式(13)有無限多組解βk,前幾組解如下:

βkx=1.875 1, 4.694 1, 7.854 8, 10.995 5，…

整個(gè)梁的一階模態(tài)及其導(dǎo)數(shù)變化如圖7所示。

圖7 一階模態(tài)曲線及其導(dǎo)數(shù)Fig.7 First order mode curve and its derivatives

當(dāng)梁受到分布力和集中力作用時(shí),式(9)變?yōu)?

P(x,t)+Fj(t)δ(x-xj)

(14)

上式的解可以表示為:

(15)

式中:ηk(t)為廣義坐標(biāo),滿足以下方程:

(16)

式中:Nk(t)為k階模態(tài)廣義力。整個(gè)飛行器受到的集中力有升降舵產(chǎn)生的Fe和前部鴨舵產(chǎn)生的Fc。為了簡(jiǎn)化,直接將以上用工程算法得到的壓力分布作為彈性機(jī)身的壓力分布:

(17)

對(duì)于前部梁,廣義力可以表示為:

(18)

對(duì)于后部梁,廣義力可以表示為;

(19)

當(dāng)Ma0=10,δe=δc=0,ψ=0.4 rad時(shí),得到機(jī)體前部廣義力Nf,1和機(jī)體后部廣義力Na,1隨迎角的變化如圖8所示。由于飛行器幾何參數(shù)有差異,本文增加了鴨舵,并且文獻(xiàn)[7]給出的是簡(jiǎn)化解析擬合式,因此曲線有一定的差異,但是在走勢(shì)上是一致的,結(jié)果具有一定的可信度。

圖8 廣義力隨迎角變化曲線Fig.8 Change of generalized force with AOA

5 非線性動(dòng)態(tài)特性分析

對(duì)于非線性系統(tǒng)可用以下方程來描述:

(20)

在標(biāo)準(zhǔn)分支分析方法中,x為n維狀態(tài)變量,μ為可變控制參數(shù),λ為m維固定控制參數(shù),通過連續(xù)的改變?chǔ)碳纯汕蟮孟到y(tǒng)的平衡曲線圖。而實(shí)際情況中往往需要多個(gè)控制參數(shù)的協(xié)調(diào)配合才能夠達(dá)到約束飛行的條件,此時(shí)就需要運(yùn)用擴(kuò)展分支分析方法才能達(dá)到目的。擴(kuò)展分支分析方法通常分為兩步。首先用標(biāo)準(zhǔn)分支分析方法求解:

(21)

式中:g(x)為k維約束方程(k

(22)

然后采用標(biāo)準(zhǔn)分支分析方法再次計(jì)算式(22),即可得到約束條件下的非線性動(dòng)力學(xué)特性。

在高超聲速飛行器縱向剛體-彈性動(dòng)力學(xué)模型[3]中,令:

(23)

則飛行器縱向剛體-彈性動(dòng)力學(xué)模型轉(zhuǎn)化為:

(24a)

(24b)

(25)

圖9 控制舵對(duì)應(yīng)關(guān)系Fig.9 Control surface schedule

從圖9可以看到,存在兩條平衡曲線,δc(δe)不滿足一一對(duì)應(yīng)關(guān)系,此時(shí)需要將這兩條曲線分離,分別建立δc(δe)關(guān)系,然后代入式(22)再次分別利用常規(guī)分支分析方法進(jìn)行平衡狀態(tài)的求解。

計(jì)算過程中保持飛行高度不變,重心xcg=0.55l。將所有計(jì)算結(jié)果匯總得到V,α,ηf和ηa隨升降舵δe變化的平衡曲線,如圖10所示。因?yàn)楸疚膶?duì)Ma0,α,δe,δc的取值有范圍限制,范圍之外的平衡解在圖10中沒有顯示。

圖10 V, α, ηf和ηa隨δe變化的平衡曲線Fig.10 Changes of V, α, ηf and ηa with δe

圖11 偏導(dǎo)數(shù)隨迎角變化曲線Fig.11 Change of partial derivative with AOA

由于彈性的存在,當(dāng)機(jī)身受到氣動(dòng)力的作用會(huì)發(fā)生彎曲,引起的迎角變化[9]可近似表示為:

(26)

令ω(x,t)=φ(x)η(t),則式(26)可以轉(zhuǎn)化為:

(27)

由圖10可知,在x=0和x=l彈性引起的迎角變化最大,利用上面求得的所有平衡狀態(tài)可以得到飛行器處于巡航平衡狀態(tài)時(shí)機(jī)體前端由于彈性引起的迎角變化最大值Δαf,max=0.109 2°,機(jī)體尾部由于彈性引起的迎角變化最大值Δαa,max=0.120 8°。

6 結(jié)論

本文對(duì)具有典型乘波體結(jié)構(gòu)的吸氣式高超聲速飛行器縱向進(jìn)行了相關(guān)分析,得到以下結(jié)論:

(1)針對(duì)高空巡航狀態(tài),利用擴(kuò)展分支分析方法,得到了系統(tǒng)的分支圖。分析發(fā)現(xiàn)該飛行器所有的平衡狀態(tài)均是不穩(wěn)定的,這主要是由于飛行器重心位于壓心之前引起的。為了實(shí)現(xiàn)穩(wěn)定,需要進(jìn)一步利用相關(guān)控制方法來實(shí)現(xiàn)。

(2)根據(jù)本文建立的彈性模型,當(dāng)該飛行器處于高空巡航平衡狀態(tài)時(shí),由于二維機(jī)身彈性引起的迎角變化較小,在初步進(jìn)行分析時(shí)可以不考慮彈性對(duì)系統(tǒng)的影響。

參考文獻(xiàn)：

[1] 唐碩,祝強(qiáng)軍.吸氣式高超聲速飛行器動(dòng)力學(xué)建模研究進(jìn)展[J].力學(xué)進(jìn)展,2011,41(2):187-200.

[2] Chavez F,Schmidt D.Analytical aeropropulsive/aeroelastic hypersonic-vehicle model with dynamic analysis[J].Journal of Guidance,Control,and Dynamics,1994,17(6):1308-1319.

[3] Bolender M A,Doman D B.A non-linear model for the longitudinal dynamics of a hypersonic air-breathing vehicle[C]//Guidance,Navigation,and Control Conference and Exhibit.San Francisco,California,2005.

[4] Clark A,Chivey Wu,Mirmirani M,et al.Development of an airframe-propulsion integrated generic hypersonic vehicle model[C]//AIAA Aerospace Sciences Meeting and Exhibit.Reno,Nevada,2006.

[5] Mehra R K,Kessel W C,Carroll J V.Global stability and control analysis of aircraft at high angle of attack[R].AD A051-850,1977.

[6] Ananthkrishnan N,Sinha N K.Level flight trim and stability analysis using an extended bifurcation and continuation procedure [J].Journal of Guidance,Control,and Dynamics,2001,24(6):1225-1228.

[7] Jason T P,Andrea S,Stephen Y,et al.Control-oriented modeling of an air-breathing hypersonic vehicle[J].Journal of Guidance,Control,and Dynamics, 2007,30(3):856-869.

[8] 姚照輝.考慮飛/推耦合特性的超然沖壓發(fā)動(dòng)機(jī)控制方法研究[D].哈爾濱:哈爾濱工業(yè)大學(xué),2010.

[9] Clark A D,Mirmirani M D,Chivey W,et al.An aero-propulsion integrated elastic model of a generic airbreathing hypersonic vehicle[C]//AIAA Guidance,Navigation,and Control Conference and Exhibit.Keystone,Colorado,2006.