雍恩米，劉深深，程艷青，錢(qián)煒祺,2

1. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000 2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

升力式再入飛行器以高超聲速再入地球大氣層時(shí)，面臨著復(fù)雜的臨近空間力熱環(huán)境，各類(lèi)大氣環(huán)境參數(shù)以及氣動(dòng)參數(shù)的不確定性對(duì)其再入運(yùn)動(dòng)具有較大影響，較小的擾動(dòng)可使其偏離預(yù)定軌跡較大甚至逸出大氣層，導(dǎo)致航天任務(wù)的失敗。飛行器自由擾動(dòng)運(yùn)動(dòng)的穩(wěn)定性分析，可為飛行器總體設(shè)計(jì)、控制律研究等方面提供重要參考。已開(kāi)展的穩(wěn)定性研究工作主要集中在航空器的動(dòng)力學(xué)穩(wěn)定性分析以及航天器中彈道式再入飛行器如導(dǎo)彈、返回飛船的姿態(tài)穩(wěn)定性分析[1-3]，而對(duì)于升力式高超聲速再入飛行器在臨近空間的動(dòng)力學(xué)系統(tǒng)穩(wěn)定性的研究相當(dāng)有限。

從公開(kāi)發(fā)表的文獻(xiàn)看，關(guān)于升力體飛行器的高超聲速動(dòng)力學(xué)穩(wěn)定性研究最早出現(xiàn)在20世紀(jì)中葉。Etkin研究了再入航天器的動(dòng)力學(xué)穩(wěn)定性問(wèn)題,將經(jīng)典的飛行器縱向動(dòng)力學(xué)理論擴(kuò)展到高速軌道飛行，研究了有動(dòng)力升力體飛行器在大氣層邊緣的圓軌道上的攝動(dòng)運(yùn)動(dòng),該項(xiàng)研究中保留了運(yùn)動(dòng)方程中的行星曲率、引力梯度和以高度為自變量的大氣密度項(xiàng)，這些項(xiàng)在低速低高度飛行器的動(dòng)力學(xué)穩(wěn)定性研究中通常被忽略[4]。Vinh和Laitone研究了航天飛機(jī)滑翔再入的縱向動(dòng)力學(xué)穩(wěn)定性問(wèn)題，引入無(wú)量綱時(shí)間變量，推導(dǎo)出了迎角的二階線性微分方程，針對(duì)兩種特殊情況，一是彈道式再入，二是小彈道傾角滑翔再入，求得了近似解，分析得出俯仰力矩導(dǎo)數(shù)是影響穩(wěn)定性的主要因素[5]。文獻(xiàn)[6]研究了帶動(dòng)力的吸氣式高超聲速飛行器長(zhǎng)周期動(dòng)力學(xué)穩(wěn)定性，發(fā)現(xiàn)推力-馬赫數(shù)曲線的斜率對(duì)長(zhǎng)周期模態(tài)的影響較大，另外還提出了基本反饋控制策略，即采用速度和高度反饋回路來(lái)改善不穩(wěn)定特性。Stich等研究了高超聲速巡航飛行器的長(zhǎng)周期特性和相關(guān)的飛行品質(zhì)問(wèn)題，在NASA的飛行仿真器上分析了馬赫數(shù)10條件下的長(zhǎng)周期模態(tài)的不穩(wěn)定特性和發(fā)動(dòng)機(jī)時(shí)間延遲對(duì)操縱品質(zhì)的影響。針對(duì)不穩(wěn)定的長(zhǎng)周期動(dòng)力學(xué)特性，推導(dǎo)了新的操縱品質(zhì)判斷準(zhǔn)則[7-9]。Robert研究了超聲速巡航縱向動(dòng)力學(xué)穩(wěn)定性，分析得出推力和氣動(dòng)外形對(duì)超聲速巡航飛行穩(wěn)定性影響較大，采用高度和速度反饋有利于巡航高度保持[10]。Lael和Darryll針對(duì)上升、穩(wěn)定巡航和周期性巡航3種任務(wù)，研究了乘波體高超聲速巡航飛行器縱向穩(wěn)定性[11]。

國(guó)內(nèi)關(guān)于升力式高超聲速飛行器模態(tài)穩(wěn)定性分析的研究較少。賈子安等研究了乘波體飛行器縱向靜穩(wěn)定特性，給出了有利于縱向靜穩(wěn)定的乘波體形狀特征[12]。陳琛等針對(duì)升力體構(gòu)型的高超聲速飛行器穩(wěn)定性問(wèn)題，在高超聲速飛行器的建模、靜穩(wěn)定性分析、自由擾動(dòng)模態(tài)穩(wěn)定性3個(gè)方面進(jìn)行較深入的研究，其中模態(tài)穩(wěn)定性分析針對(duì)特征點(diǎn)，且縱向運(yùn)動(dòng)主要關(guān)注短周期模態(tài)特性[13]。李銳研究了高超聲速飛行器軌跡優(yōu)化方法以及橫航向穩(wěn)定性判據(jù)，分析了給定外形高超聲速飛行器在整個(gè)速度范圍內(nèi)的橫航向穩(wěn)定性[14]。蘇二龍和羅建軍研究了滑翔式高超聲速飛行器大迎角橫側(cè)向失穩(wěn)問(wèn)題，采用延拓算法和分岔理論求解，對(duì)平衡分支的穩(wěn)定性和突變點(diǎn)進(jìn)行了分析[15]。高清和李潛對(duì)美國(guó)典型高超聲速飛行器橫側(cè)向穩(wěn)定性進(jìn)行了分析[16]。葉友達(dá)等研究了類(lèi)HTV-2高超聲速飛行器滾轉(zhuǎn)運(yùn)動(dòng)的穩(wěn)定性以及滾轉(zhuǎn)運(yùn)動(dòng)的失穩(wěn)判據(jù)[17]。

本文研究了面向彈道優(yōu)化問(wèn)題的高超聲速飛行器縱向模態(tài)穩(wěn)定性，即針對(duì)無(wú)動(dòng)力再入的高超聲速飛行器，考慮大氣密度隨高度的變化和引力梯度，分析高超聲速飛行器自由擾動(dòng)運(yùn)動(dòng)的縱向穩(wěn)定性，即縱向模態(tài)穩(wěn)定性，并基于不同最優(yōu)彈道的穩(wěn)定性分析結(jié)果，提出對(duì)高超聲速飛行器彈道優(yōu)化設(shè)計(jì)的建議。

1 高超聲速飛行器再入縱向運(yùn)動(dòng)模型

針對(duì)再入高超聲速飛行器運(yùn)動(dòng)特性和彈道穩(wěn)定性分析的需求，考慮地球?yàn)榉切D(zhuǎn)圓球時(shí)，飛行器無(wú)動(dòng)力再入縱向動(dòng)力學(xué)方程為[5]

(1)

(2)

(3)

(4)

(5)

θ=γ+α

(6)

式中：V、γ、q、r、θ、α分別為飛行器速度、航跡角、俯仰角速率、地心距、俯仰角和迎角；m為飛行器質(zhì)量；Ix、Iy、Iz為主慣量矩；S為氣動(dòng)參考面積；l為參考長(zhǎng)度。主要運(yùn)動(dòng)參數(shù)如圖1所示，其中L、D、W分別為飛行器所受升力、阻力和重力，x軸為飛行器體軸縱向?qū)ΨQ(chēng)軸。大氣密度ρ、引力加速度g是隨高度變化的函數(shù)。氣動(dòng)力系數(shù)包括升力系數(shù)CL、阻力系數(shù)CD和俯仰力矩系數(shù)Cm。上述氣動(dòng)系數(shù)是高度h、馬赫數(shù)Ma和迎角α的函數(shù),即

(7)

圖1 縱向運(yùn)動(dòng)主要參數(shù)示意圖Fig.1 Sketch map of main parameters of longitudinal motion

2 縱向動(dòng)力學(xué)方程小擾動(dòng)線性化

為分析高超聲速飛行器的模態(tài)穩(wěn)定性，即自由擾動(dòng)運(yùn)動(dòng)的穩(wěn)定性，引入基準(zhǔn)運(yùn)動(dòng)和小擾動(dòng)運(yùn)動(dòng)的概念，建立小擾動(dòng)運(yùn)動(dòng)方程，使運(yùn)動(dòng)方程線性化，然后獲得狀態(tài)方程和轉(zhuǎn)移矩陣，在此基礎(chǔ)上進(jìn)行沿彈道的縱向模態(tài)分析。

首先根據(jù)式(2)、式(5)、式(6)得到關(guān)于迎角的微分方程，即

(8)

對(duì)式(1)～式(4)和式(8)進(jìn)行小擾動(dòng)線性化，得到

(9)

(10)

(11)

(12)

(13)

縱向運(yùn)動(dòng)方程小擾動(dòng)線性化過(guò)程中采用的大氣密度模型為指數(shù)模型，即

ρ=ρ0e[-β(r-R0)]

(14)

式中：β為指數(shù)大氣密度常數(shù)；R0為地球平均半徑。因而有大氣密度對(duì)地心距的導(dǎo)數(shù)為

(15)

引力模型采用平方反比模型，即

(16)

從而得到引力對(duì)地心距的導(dǎo)數(shù)為

(17)

(18)

(19)

那么縱向擾動(dòng)運(yùn)動(dòng)的穩(wěn)定性可以由求解特征方程D(s)=|sI-A|的特征根進(jìn)行分析。

3 高超聲速飛行器氣動(dòng)特性分析

為開(kāi)展高超聲速飛行器模態(tài)穩(wěn)定性分析，還需要獲得飛行器氣動(dòng)力數(shù)據(jù)。本文以美國(guó)高超聲速飛行器HX(Hypersonic X-plane)為背景，進(jìn)行類(lèi)HX氣動(dòng)布局設(shè)計(jì)，利用二次曲線及基于類(lèi)型函數(shù)和形狀函數(shù)的CST(Class function and Shape function Transformation technique)方法[18-19]，提出類(lèi)HX升力體飛行器氣動(dòng)外形布局，其計(jì)算網(wǎng)格如圖2所示。最后采用工程估算方法獲得飛行器的氣動(dòng)力特性數(shù)據(jù)。

圖2 飛行器計(jì)算網(wǎng)格Fig.2 Computational grid for vehicle

由于高超聲速飛行器飛行空域范圍較大，因此對(duì)于不同流區(qū)，采用了不同氣動(dòng)力工程計(jì)算方法。其中連續(xù)流區(qū)采用了修正牛頓理論。背風(fēng)面修正利用Prandtl-Meyer公式，該類(lèi)方法已經(jīng)在高超聲速飛行器計(jì)算中得到了廣泛應(yīng)用[20-21]，其精度基本滿(mǎn)足彈道設(shè)計(jì)與穩(wěn)定性分析的精度要求。自由分子流區(qū)的氣動(dòng)特性則采用碰撞模型進(jìn)行計(jì)算，而稀薄流區(qū)氣動(dòng)特性則采用當(dāng)?shù)鼗椒?，連續(xù)流和自由分子流氣動(dòng)力系數(shù)通過(guò)橋函數(shù)得到。

基于上述方法，得到了高超聲速飛行器在不同高度、馬赫數(shù)、迎角和側(cè)滑角下的氣動(dòng)特性。本文研究高超聲速飛行器的彈道穩(wěn)定性，因此主要考慮縱向氣動(dòng)力特性，圖3給出了典型狀態(tài)下的升力系數(shù)、阻力系數(shù)和俯仰力矩系數(shù)隨迎角和舵偏角的變化曲面。

基于工程估算方法獲得飛行器的氣動(dòng)力特性數(shù)據(jù)，進(jìn)行縱向氣動(dòng)數(shù)據(jù)配平，可以進(jìn)一步分析該高超聲速飛行器的升阻特性。圖4為給定迎角為10°，不同高度條件下升阻比隨馬赫數(shù)的變化曲線。圖5為給定馬赫數(shù)條件下，不同高度升阻比隨迎角的變化曲線。分析圖4、圖5的結(jié)果可知：① 該飛行器在典型狀態(tài)下，最大升阻比在2.4～2.5左右，屬于高升阻比高超聲速飛行器；② 該飛行器的最大升阻比迎角約為10°,在第4節(jié)模態(tài)分析設(shè)計(jì)典型彈道的標(biāo)準(zhǔn)迎角剖面時(shí)，可以此為依據(jù)，即給定標(biāo)準(zhǔn)迎角剖面應(yīng)盡量發(fā)揮升阻比的優(yōu)勢(shì)；③ 在 給定高度條件下，該布局的高超聲速飛行器升阻特性隨迎角變化較大，而隨馬赫數(shù)變化較小，特別是飛行速度在馬赫數(shù)10～20這個(gè)范圍內(nèi)。在設(shè)計(jì)高超聲速飛行器再入典型彈道時(shí)，可以利用該特性，對(duì)所采用的氣動(dòng)力模型進(jìn)行簡(jiǎn)化和建模，即主要考慮迎角變量對(duì)升阻特性的影響。

圖3 典型狀態(tài)縱向氣動(dòng)力特性Fig.3 Longitude aerodynamic characteristic on typical states

圖4 升阻比隨馬赫數(shù)的變化Fig.4 Variation of lift-to-drag ratio with Mach number

圖5 升阻比隨迎角變化Fig.5 Variation of lift-to-drag ratio with angle of attack

4 典型彈道模態(tài)穩(wěn)定性分析

4.1 優(yōu)化彈道模態(tài)

優(yōu)化計(jì)算得到最大射程彈道高度和馬赫數(shù)變化曲線如圖6(a)所示，其飛行時(shí)間為3 423 s，射程為18 343 km。最小射程彈道如圖6(b)所示，其飛行時(shí)間為1 500 s，射程為7 606 km。

為從穩(wěn)定性角度分析不同極限彈道的特性，這里分別計(jì)算了最大、最小射程彈道的模態(tài)特性，即沿彈道求解特征方程的特征根。早期關(guān)于高超聲速飛行器縱向長(zhǎng)周期運(yùn)動(dòng)的研究表明，當(dāng)考慮大氣密度隨高度的變化時(shí)，高超聲速飛行器相較于一般速域的飛行器，增加了一個(gè)長(zhǎng)周期模態(tài)，即高度模態(tài)，同時(shí)也存在兩個(gè)常規(guī)模態(tài)，即沉浮模態(tài)和短周期模態(tài)。本文通過(guò)模態(tài)分析，同樣也得到3種模態(tài)和對(duì)應(yīng)的5個(gè)特征根。

首先，分析最大射程彈道模態(tài)穩(wěn)定性。由于5個(gè)特征根中的一個(gè)實(shí)根與其他4個(gè)根的數(shù)量級(jí)相差較大，因此模態(tài)圖分別在兩幅圖中畫(huà)出，分別對(duì)應(yīng)一個(gè)實(shí)根與其他4個(gè)根。最大射程彈道對(duì)應(yīng)的根軌跡如圖6(c)和圖6(e)所示。圖6(e)為一個(gè)大實(shí)根對(duì)應(yīng)的模態(tài)，即高度模態(tài)，根軌跡基本位于負(fù)實(shí)軸上，這表明高度模態(tài)是穩(wěn)定的且呈現(xiàn)非周期、快速衰減特性。圖6(c)中一對(duì)位于復(fù)平面左半平面的共軛復(fù)根對(duì)應(yīng)短周期模態(tài)，即表征迎角和俯仰角速率的自由擾動(dòng)穩(wěn)定性變化。該短周期模態(tài)頻率沿彈道呈現(xiàn)由低到高變化的趨勢(shì)，且在再入初始階段，由于動(dòng)壓較低，短周期模態(tài)的頻率與長(zhǎng)周期模態(tài)頻率相當(dāng)。圖6(c)中位于復(fù)平面實(shí)軸附近的兩個(gè)根為沉浮模態(tài)，該模態(tài)對(duì)應(yīng)的根軌跡在實(shí)軸附近。在再入初段，沉浮模態(tài)為一對(duì)穩(wěn)定的共軛復(fù)根，進(jìn)入中段和末段，沉浮模態(tài)變?yōu)閮蓚€(gè)實(shí)根，且趨于不穩(wěn)定。圖6(d)和圖6(f) 為最小射程彈道模態(tài)分析圖。類(lèi)似的，圖6(f)為最小射程彈道高度模態(tài)圖，與最大射程高度模態(tài)相比最大射程具有更大的負(fù)實(shí)根，說(shuō)明對(duì)應(yīng)的狀態(tài)衰減更快。圖6(d)為最小射程彈道短周期模態(tài)和沉浮模態(tài)圖：其中短周期模態(tài)也是穩(wěn)定的，其頻率由小變大再變小；最小射程彈道的沉浮模態(tài)與最大射程彈道有相似特性，不同之處在于其位于左半平面的實(shí)根數(shù)目明顯增多，即沿最小射程彈道，穩(wěn)定的沉浮模態(tài)特征根增多。

表1 彈道條件和過(guò)程約束Table 1 Conditions and process constraints of trajectory

圖6 優(yōu)化彈道模態(tài)分析Fig.6 Modal analysis of optimal trajectories

4.2 無(wú)約束跳躍彈道模態(tài)

前面分析了考慮過(guò)程約束的最大、最小射程彈道的模態(tài)特性，為了對(duì)比分析更多形式的高超聲速再入彈道模態(tài)特性，本節(jié)計(jì)算了不考慮路徑約束的最大射程彈道和無(wú)約束跳躍彈道，并進(jìn)行模態(tài)穩(wěn)定性分析。

上述兩種彈道的計(jì)算方法為：不考慮過(guò)程約束的最大射程彈道仍采用文獻(xiàn)[22]的優(yōu)化方法，只是無(wú)動(dòng)壓等過(guò)程約束，初始條件和終端約束同表1，計(jì)算得到該彈道飛行時(shí)間為3 926 s，射程為22 144 km，其高度、馬赫數(shù)曲線如圖7(a)所示。而無(wú)約束跳躍彈道直接生成，初始高度速度同表1，初始彈道傾角由0°改為5°，以實(shí)現(xiàn)彈道在高度空間的大幅度跳躍，其側(cè)傾角為取0°，即升力朝上，迎角仍采用標(biāo)準(zhǔn)彈道迎角剖面，結(jié)束條件為高度下降到30 km。跳躍彈道高度、馬赫數(shù)曲線如圖7(b)所示，飛行時(shí)間為4 859 s，射程為27 732 km。

兩種彈道的模態(tài)特性如圖7(c)～圖7(f)所示。兩種彈道的高度模態(tài)均為負(fù)的大實(shí)根，且跳躍彈道的高度模態(tài)衰減更快。無(wú)過(guò)程約束最大射程彈道與考慮約束的最大射程彈道相比較，射程增大了約18%，但沉浮模態(tài)的不穩(wěn)定狀態(tài)增多。跳躍彈道與考慮約束的最大射程彈道相比，射程增加了約33%，但短周期模態(tài)出現(xiàn)了不穩(wěn)定點(diǎn)，位于正實(shí)軸上，對(duì)應(yīng)彈道為初始上升段。跳躍彈道的沉浮模態(tài)沒(méi)有出現(xiàn)共軛復(fù)根，均位于實(shí)軸上，且大部分都位于復(fù)平面右平面不穩(wěn)定區(qū)域。

圖7 兩種無(wú)約束彈道模態(tài)分析Fig.7 Modal analysis of two types of unconstrained trajectories

綜合分析4條典型彈道(圖6、圖7)的模態(tài)特性可以得出以下結(jié)論：① 當(dāng)考慮大氣密度隨高度變化時(shí)，高超聲速再入飛行器具有兩個(gè)長(zhǎng)周期模態(tài)，即高度模態(tài)和沉浮模態(tài)，且高度模態(tài)為穩(wěn)定的非周期模態(tài);② 考慮過(guò)程約束的最大射程彈道和最小射程彈道都具有穩(wěn)定的短周模態(tài)，且最小射程彈道穩(wěn)定的沉浮模態(tài)特征根更多。

從無(wú)約束最大射程彈道和跳躍彈道的模態(tài)分析結(jié)果來(lái)看，高超聲速再入飛行器，在高度范圍內(nèi)跳躍變化雖然能夠在很大程度上增大射程，但同時(shí)會(huì)使短周期模態(tài)和沉浮模態(tài)產(chǎn)生更多不穩(wěn)定特征根。因此，從穩(wěn)定性的角度，可以對(duì)彈道設(shè)計(jì)提出建議：應(yīng)避免所設(shè)計(jì)的彈道產(chǎn)生太大的跳躍，即使是犧牲一些射程上的性能。

4.3 最大射程彈道典型狀態(tài)

前面分析了典型彈道的模態(tài)特性，現(xiàn)選取其中最大射程彈道上的再入初段、滑翔段和再入末段的3個(gè)狀態(tài)點(diǎn)作為典型狀態(tài)(見(jiàn)表2)，給出具體的特征根并進(jìn)行穩(wěn)定性分析。典型狀態(tài)點(diǎn)對(duì)應(yīng)的模態(tài)特征根和無(wú)阻尼自振頻率見(jiàn)表3。分析表3 可知：再入初段短周期和長(zhǎng)周期模態(tài)的頻率都比較低，沒(méi)有明顯的短周期和長(zhǎng)周期區(qū)別；滑翔段短周期模態(tài)頻率明顯增大，沉浮模態(tài)由穩(wěn)定的復(fù)根變?yōu)閮蓚€(gè)實(shí)根，且有一個(gè)為不穩(wěn)定實(shí)根，高度模態(tài)變?yōu)榇蟮呢?fù)實(shí)根；再入末段，短周期模態(tài)和高度模態(tài)無(wú)阻尼自振頻率繼續(xù)增大，且都是穩(wěn)定的。

表2 典型飛行時(shí)刻狀態(tài)參數(shù)Table 2 State parameters of typical flight time

表3 典型飛行時(shí)刻特征根及頻率Table 3 Eigen values and frequencies of typical flight time

5 結(jié) 論

本文分析了高超聲速再入飛行器的縱向模態(tài)穩(wěn)定性。引入基準(zhǔn)運(yùn)動(dòng)和小擾動(dòng)運(yùn)動(dòng)概念，建立無(wú)動(dòng)力再入縱向動(dòng)力學(xué)小擾動(dòng)運(yùn)動(dòng)方程，然后獲得狀態(tài)方程和轉(zhuǎn)移矩陣，在此基礎(chǔ)上進(jìn)行沿彈道的縱向模態(tài)分析。針對(duì)最大射程彈道、最小射程彈道和跳躍彈道以及最大射程彈道典型狀態(tài)的模態(tài)分析，得出以下結(jié)論：

1) 再入初段短周期模態(tài)和長(zhǎng)周期模態(tài)的頻率都比較低，沒(méi)有明顯的短周期和長(zhǎng)周期區(qū)別。

2) 高超聲速再入飛行器具有兩個(gè)長(zhǎng)周期模態(tài)，即高度模態(tài)和沉浮模態(tài)，且高度模態(tài)為穩(wěn)定的非周期模態(tài)。

3) 考慮過(guò)程約束的最大射程彈道和最小射程彈道都具有穩(wěn)定的短周模態(tài)，且最小射程彈道穩(wěn)定的沉浮模態(tài)特征根更多。

4) 高超聲速再入飛行器，在高度范圍內(nèi)跳躍變化雖然能夠在很大程度上增大射程但同時(shí)會(huì)使短周期模態(tài)和沉浮模態(tài)產(chǎn)生更多不穩(wěn)定特征根。

高超聲速飛行器縱向模態(tài)分析，即自由擾動(dòng)運(yùn)動(dòng)的穩(wěn)定性分析可為飛行器總體設(shè)計(jì)、制導(dǎo)控制律研究等方面提供重要參考。