徐敏, 白斌, 祝小平

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

0 引言

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

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

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

1 幾何外形描述

如文獻[3],縱向剖面視圖及機體坐標系Oxyz如圖1所示。

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

飛行器在 30 km的高度以Ma0=10的速度進行巡航運動。聲速計算方法如下:

(1)

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

2 氣動計算

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

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

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

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

表1 各參數取值范圍Table 1 Range of the parameters

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

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

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

3 推力計算

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

圖4 發(fā)動機模型Fig.4 Engine model

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

A2～A3燃燒室內,為避免碳氫燃料燃燒室產生復雜的物理化學變化給推力的計算帶來困難,采用氫氣為燃料,根據一維Rayleigh流理論可知:

(2)

(3)

(4)

(5)

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

為簡化起見,本文引用文獻[8]給出的燃燒效率與燃油當量比的簡單函數關系進行計算。T04,T03為總溫,P04,P03為總壓,存在以下關系:

(6)

(7)

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

(8)

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

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

4 彈性模型

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

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

彎曲振動微分方程為:

(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，…

整個梁的一階模態(tài)及其導數變化如圖7所示。

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

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

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

(14)

上式的解可以表示為:

(15)

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

(16)

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

(17)

對于前部梁,廣義力可以表示為:

(18)

對于后部梁,廣義力可以表示為;

(19)

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

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

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

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

(20)

在標準分支分析方法中,x為n維狀態(tài)變量,μ為可變控制參數,λ為m維固定控制參數,通過連續(xù)的改變μ即可求得系統(tǒng)的平衡曲線圖。而實際情況中往往需要多個控制參數的協(xié)調配合才能夠達到約束飛行的條件,此時就需要運用擴展分支分析方法才能達到目的。擴展分支分析方法通常分為兩步。首先用標準分支分析方法求解:

(21)

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

(22)

然后采用標準分支分析方法再次計算式(22),即可得到約束條件下的非線性動力學特性。

在高超聲速飛行器縱向剛體-彈性動力學模型[3]中,令:

(23)

則飛行器縱向剛體-彈性動力學模型轉化為:

(24a)

(24b)

(25)

圖9 控制舵對應關系Fig.9 Control surface schedule

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

計算過程中保持飛行高度不變,重心xcg=0.55l。將所有計算結果匯總得到V,α,ηf和ηa隨升降舵δe變化的平衡曲線,如圖10所示。因為本文對Ma0,α,δe,δc的取值有范圍限制,范圍之外的平衡解在圖10中沒有顯示。

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

圖11 偏導數隨迎角變化曲線Fig.11 Change of partial derivative with AOA

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

(26)

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

(27)

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

6 結論

本文對具有典型乘波體結構的吸氣式高超聲速飛行器縱向進行了相關分析,得到以下結論:

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

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

參考文獻：

[1] 唐碩,祝強軍.吸氣式高超聲速飛行器動力學建模研究進展[J].力學進展,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].哈爾濱:哈爾濱工業(yè)大學,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.