顧杰，張曙光,*，楊帆，王保印

1.北京航空航天大學(xué) 交通科學(xué)與工程學(xué)院，北京 100083 2.成都飛機(jī)設(shè)計(jì)研究所，成都 610091

再入飛行器沉浮特性近似解析及應(yīng)用

顧杰1，張曙光1,*，楊帆2，王保印1

1.北京航空航天大學(xué) 交通科學(xué)與工程學(xué)院，北京 100083 2.成都飛機(jī)設(shè)計(jì)研究所，成都 610091

采用基于平衡滑翔的數(shù)值或解析預(yù)測(cè)-校正再入制導(dǎo)方法的再入飛行器，從初始下降段到平衡滑翔段過渡或出現(xiàn)較大預(yù)測(cè)偏差時(shí)易產(chǎn)生沉浮振蕩，且隨著近年來所研究飛行器升阻比的增加，沉浮振蕩更加明顯，從而引起了研究者對(duì)高超聲速沉浮特性的重新審視。首先,通過三階縱向動(dòng)態(tài)方程及平衡滑翔條件推導(dǎo)出了形式簡(jiǎn)潔、能直觀表達(dá)主要影響因素的再入飛行器高超聲速沉浮特性近似解。在此基礎(chǔ)上，分析發(fā)現(xiàn)高超聲速沉浮阻尼特性隨高度的變化規(guī)律主要由軌道速度比和沉浮修正參數(shù)主導(dǎo)，澄清了以往對(duì)大氣密度梯度參數(shù)影響的猜測(cè)。最后,推導(dǎo)出再入軌跡振蕩抑制器設(shè)計(jì)的近似解析關(guān)系，進(jìn)一步完善了基于平衡滑翔的數(shù)值或解析預(yù)測(cè)-校正再入制導(dǎo)方法，仿真驗(yàn)證表明該方法能夠有效抑制再入軌跡的沉浮振蕩。

再入；高超聲速飛行器；沉??；阻尼；解析解

近十多年來，基于平衡滑翔的數(shù)值或解析預(yù)測(cè)-校正再入制導(dǎo)方法，以其良好的自適應(yīng)性得到了越來越多的關(guān)注[1-11]。該方法設(shè)計(jì)過程簡(jiǎn)潔，不需要像傳統(tǒng)的航天飛機(jī)再入制導(dǎo)方法[12]一樣預(yù)先精心設(shè)計(jì)依賴于任務(wù)的參考剖面，因而能夠節(jié)省大量人力物力。但是預(yù)測(cè)-校正制導(dǎo)方法在出現(xiàn)較大的模型偏差或遇到黑障區(qū)較大的導(dǎo)航誤差[13-14]時(shí)，將產(chǎn)生較大的校正偏差，從而使得再入軌跡偏離平衡滑翔狀態(tài)，進(jìn)而引起沉浮振蕩，特別是從初始下降段到平衡滑翔段過渡時(shí)最為明顯[3-4]。這些不確定性較大的軌跡振蕩可能會(huì)導(dǎo)致飛行器超出氣動(dòng)熱、動(dòng)壓和過載等約束限制，從而使得熱防護(hù)結(jié)構(gòu)和乘員面臨危險(xiǎn)[2,15]。

為消除以上再入飛行過程中因沉浮振蕩產(chǎn)生的安全風(fēng)險(xiǎn)，Lu[1-4]、Shen[10-11]和Xue[7]等在尋找有效消除再入軌跡沉浮振蕩的自適應(yīng)方法上進(jìn)行了大量工作，并得到了有益的結(jié)果[4]。Shen通過預(yù)測(cè)搜索得到了再入初始下降段與平衡滑翔段平滑過渡的轉(zhuǎn)換點(diǎn)信息[10-11]。Xue則在Shen的基礎(chǔ)上引入與地球自轉(zhuǎn)相關(guān)的修正因子進(jìn)一步提高了再入軌跡的平滑性[7]。但以上方法均是基于模型的預(yù)測(cè)控制，其控制效果依賴于模型計(jì)算的準(zhǔn)確性。而Lu[1-4]則在多年工作的基礎(chǔ)上采用簡(jiǎn)潔的反饋控制抑制了再入軌跡振蕩，但是其反饋增益是依據(jù)經(jīng)驗(yàn)設(shè)定的，缺乏較為直接的理論依據(jù)。因而為了從理論上完善Lu的沉浮振蕩抑制方法，且隨著近年來對(duì)易于產(chǎn)生沉浮振蕩的中高升阻比高超聲速飛行器研究的增多[4,16-18]和相關(guān)問題的顯現(xiàn)[3-4]，以及近期在火星滑翔再入研究中升力飛行器易于跳躍問題的出現(xiàn)[19]，使得高超聲速飛行器沉浮振蕩特性的研究重新受到了關(guān)注和審視[20-21]。

航空發(fā)展早期，“沉浮”(Phugoid)一詞最早是英國學(xué)者Lanchester[22]于1908年在研究低速飛機(jī)飛行穩(wěn)定性時(shí)提出的，其將低速飛機(jī)水平飛行時(shí)相對(duì)較慢的軌跡俯沖上仰運(yùn)動(dòng)命名為“沉浮”運(yùn)動(dòng)。Lanchester在迎角近似為常值和切向力平衡的假設(shè)下，依據(jù)機(jī)械能守恒推導(dǎo)出了低速飛機(jī)沉浮模態(tài)周期的簡(jiǎn)潔解析表達(dá)式，表明低速飛機(jī)沉浮模態(tài)的周期僅與飛行速度有關(guān)。但由于Lanchester方法未考慮能量耗散，故無法得到任何阻尼信息。此后，文獻(xiàn)[23-26]分別推導(dǎo)出了較為準(zhǔn)確的低速飛機(jī)沉浮模態(tài)特征根或特征方程系數(shù)的解析表達(dá)式；而Etkin[23]則得到了最有直觀指導(dǎo)意義的結(jié)果，其所得解析解表明低速飛機(jī)的沉浮阻尼比完全由升阻比大小和動(dòng)力系統(tǒng)類型確定。

研究低空低速飛行時(shí)，可認(rèn)為大氣密度近似不變以及大氣具有不可壓縮性。但對(duì)于中高空高速飛行則需要考慮大氣密度和壓縮性的影響。Scheubel[27]于1942年在Lanchester研究基礎(chǔ)上考慮了大氣密度隨高度變化的影響；Neumark[28]還考慮了壓縮性影響，推導(dǎo)出復(fù)雜但較為準(zhǔn)確的超聲速飛行沉浮模態(tài)特征根近似解析表達(dá)式。

隨著高超聲速飛行器的發(fā)展，Etkin[29]于1961年首次發(fā)表采用直接數(shù)值方法的高超聲速縱向飛行動(dòng)力學(xué)研究結(jié)果，指出該研究必須考慮大氣密度隨高度變化和由地球曲率產(chǎn)生的“離心力”的影響，且在接近軌道速度時(shí)重力隨高度變化的影響會(huì)顯著起來。由Etkin的數(shù)值結(jié)果可以看出高超聲速飛行的沉浮周期隨高度增加而增加，并趨近于對(duì)應(yīng)高度的軌道周期，但不受推力規(guī)律的影響；而沉浮阻尼特性則對(duì)推力規(guī)律較為敏感。Etkin通過對(duì)沉浮阻尼特性隨高度變化曲線和大氣密度梯度參數(shù)隨高度變化曲線的對(duì)比，認(rèn)為高超聲速飛行器在90多千米高度處出現(xiàn)“重阻尼”是由于大氣密度梯度參數(shù)的絕對(duì)值在該高度出現(xiàn)峰值導(dǎo)致的，后面將對(duì)這一看法作較深入討論。

為了顯式表達(dá)各種參數(shù)對(duì)高超聲速飛行器沉浮特性的影響，Laitone和Chou[30]推導(dǎo)出了常值推力下高超聲速飛行器沉浮周期和阻尼項(xiàng)的近似解析表達(dá)式；Vinh和Dobrzelecki[31]則通過引入一個(gè)由拉格朗日展開法求解出的無窮級(jí)數(shù)，對(duì)Laitone和Chou的結(jié)果進(jìn)行了修正，表達(dá)式更加準(zhǔn)確，但是也變得更加復(fù)雜。

20世紀(jì)80年代末，隨著以超燃沖壓發(fā)動(dòng)機(jī)為關(guān)鍵技術(shù)的美國國家空天飛機(jī)計(jì)劃(National Aero-Space Plane，NASP)的啟動(dòng)，部分學(xué)者開始了關(guān)于推力規(guī)律對(duì)高超聲速飛行器縱向動(dòng)力學(xué)特性影響的深入研究。Berry[32]認(rèn)為空天飛機(jī)二階變化的推力對(duì)沉浮模態(tài)和高度模態(tài)具有顯著影響；Markopoulos等[33]則進(jìn)行了隨速度、高度變化而變化的推力對(duì)高超聲速沉浮特性的影響研究，結(jié)果表明當(dāng)推力規(guī)律函數(shù)為相關(guān)參數(shù)一階及一階以下變化律時(shí)，高超聲速飛行器沉浮周期不受推力規(guī)律影響，這與Etkin和Berry的結(jié)論一致。

對(duì)于高超聲速滑翔再入研究，F(xiàn)erreira[34]以及Chen等[20]分別采用Liouville變換和通用多尺度理論求解出了滑翔再入沉浮運(yùn)動(dòng)軌跡的解析表達(dá)式(由平衡滑翔項(xiàng)和阻尼振蕩項(xiàng)組成)和整個(gè)再入過程軌跡振蕩次數(shù)的表達(dá)式，兩者結(jié)果基本相同；胡錦川和陳萬春[21]則在進(jìn)行平穩(wěn)滑翔彈道研究過程中通過高度偏差的二階微分方程，得到了滑翔再入沉浮振蕩的自然頻率和阻尼比的近似解析表達(dá)式，但是滑翔再入縱向沉浮運(yùn)動(dòng)為三階動(dòng)態(tài)系統(tǒng)，故其近似解析解的精度還有提高的空間。

綜上所述，迄今對(duì)飛行器沉浮運(yùn)動(dòng)特性的研究已經(jīng)相對(duì)全面。但是已有的高超聲速沉浮特性解析式均較復(fù)雜，還沒有像Etkin等所推導(dǎo)出的低速飛機(jī)沉浮特性解析式一樣具有直觀指導(dǎo)意義的簡(jiǎn)潔結(jié)果，因而不易理解再入飛行器沉浮特性隨再入軌跡的變化規(guī)律及影響因素。為此，本文在已有研究基礎(chǔ)上，利用平衡滑翔條件，推導(dǎo)出再入飛行器沉浮運(yùn)動(dòng)特性的簡(jiǎn)潔近似解析式，并在沉浮特性規(guī)律解釋和軌跡振蕩抑制兩方面進(jìn)行了應(yīng)用，其對(duì)于再入運(yùn)動(dòng)分析與控制設(shè)計(jì)均具有指導(dǎo)意義。

1 再入運(yùn)動(dòng)方程

1.1 非線性再入運(yùn)動(dòng)方程

為進(jìn)行較為精確的飛行器再入運(yùn)動(dòng)分析和仿真計(jì)算，考慮地球自轉(zhuǎn)引起的向心加速度和科氏加速度、地球曲率引起的向心加速度等相關(guān)項(xiàng)，且認(rèn)為地球大氣是靜止的，即地球大氣隨地球以相同角速度自轉(zhuǎn)，則均質(zhì)圓地球條件下的三自由度再入質(zhì)點(diǎn)運(yùn)動(dòng)方程為[7]

(1)

(2)

(3)

(4)

(5)

(6)

式中：m為再入飛行器的質(zhì)量；r為地球球心到飛行器質(zhì)心的距離；g=μ/r2為重力加速度，μ=3.986 031 954×1014m3/s2為地球引力常數(shù)[35-36]；Ω為地球自轉(zhuǎn)角速度；θ和φ分別為地球經(jīng)度和緯度；V為飛行器相對(duì)地球表面的速度；γ為航跡角(飛行器速度矢量高于當(dāng)?shù)厮矫鏋檎?；ψ為航向角(飛行器速度矢量從北向開始順時(shí)針為正)；X、Y和Z為除重力外作用于飛行器上的合力在3個(gè)方向上的分量，分別為切向力(沿速度方向?yàn)檎?、側(cè)向力(垂直于速度矢量所在的鉛垂面向右為正)和法向力(垂直于速度矢量向下為正)。對(duì)于無動(dòng)力再入飛行器來說，X、Y和Z可以表示為

(7)

ρ=ρreβ(r-Re)

(8)

式中：ρr為參考大氣密度；β=(?ρ/?r)/ρ為大氣密度高度變化率與當(dāng)?shù)卮髿饷芏鹊谋戎?，?jiǎn)稱大氣密度梯度參數(shù)；Re=6 371.2×103m為地球平均半徑。

圖1直觀給出了0～180 km高度下幾種大氣密度梯度參數(shù)β的值，其中實(shí)線依據(jù)1966年美國標(biāo)準(zhǔn)大氣擬合數(shù)據(jù)得到[31,34]，其他3種線型是為了便于對(duì)比分析β的影響而給出的,β通常取為常值，本文中若無特別說明，取以下常值[29,35-36](圖1點(diǎn)劃線，最接近真實(shí)情況)：

(9)

式中：hr為大氣密度梯度參數(shù)的參考高度。

本文選用航天飛機(jī)軌道器作為算例飛行器，其再入質(zhì)量和氣動(dòng)參考面積分別為80 000 kg和249.91 m2[36]?？紤]到高超聲速飛行氣動(dòng)系數(shù)的Oswatitsch馬赫數(shù)無關(guān)性原理[29],其高超聲速升力系數(shù)和阻力系數(shù)可近似擬合為[36]

(10)

式中：α為迎角。

圖1 地球大氣密度梯度參數(shù) Fig.1 Density-gradient parameters of Earth’s atmosphere

1.2 線性再入沉浮動(dòng)態(tài)方程

飛行器再入時(shí)地球自轉(zhuǎn)引起的向心加速度特征量?jī)H占重力加速度的0.35%，隨速度增加而增加的科氏加速度特征量在軌道速度時(shí)也不超過重力加速度的10%，因而在進(jìn)行初步再入運(yùn)動(dòng)分析時(shí)忽略地球自轉(zhuǎn)不會(huì)對(duì)分析結(jié)果產(chǎn)生較大影響[29]，且數(shù)值分析也表明忽略地球自轉(zhuǎn)對(duì)沉浮動(dòng)態(tài)特性影響很小，故可令Ω=0,代入式(1)、式(4)和式(5)得到再入沉浮動(dòng)態(tài)方程為

(11)

式中：x=[Vγr]T和u=[kσ]=[cosσ]分別為狀態(tài)向量和控制向量。

Ferreira[34]和Chen等[20]的研究已經(jīng)表明，滑翔再入沉浮運(yùn)動(dòng)由平衡滑翔項(xiàng)和阻尼振蕩項(xiàng)組成，因而可取平衡滑翔狀態(tài)為基準(zhǔn)狀態(tài)，并可由Ferreira[34]和Chen等[20]推導(dǎo)的平衡滑翔近似解得到。在基準(zhǔn)狀態(tài)對(duì)再入沉浮動(dòng)態(tài)方程進(jìn)行李雅普諾夫線性化[37-38]，即可得到線性再入沉浮動(dòng)態(tài)方程為

(12)

(13)

B=[0 -Zkσ/(mV) 0]T

(14)

式中：前綴Δ表示擾動(dòng)量；Δx=[ΔVΔγΔr]T為狀態(tài)向量；Δu=[Δkσ]為控制向量；A為狀態(tài)矩陣；B為控制矩陣；A13、A21和A23分別為

(15)

飛行器高超聲速再入飛行時(shí)，由于氣動(dòng)系數(shù)具有馬赫數(shù)無關(guān)性，因而氣動(dòng)力主要由迎角決定。而迎角為短周期參數(shù)，其相對(duì)于緩慢變化的沉浮運(yùn)動(dòng)來說可以認(rèn)為是常值。故高超聲速飛行時(shí)的縱向氣動(dòng)導(dǎo)數(shù)可表示為

(16)

2 再入沉浮特性參數(shù)一次近似

2.1 特征方程

鑒于中高升阻比飛行器平衡滑翔再入過程中航跡角及航跡角變化率都很小[34]，故在對(duì)以平衡滑翔狀態(tài)為基準(zhǔn)狀態(tài)的線性再入運(yùn)動(dòng)方程推導(dǎo)分析時(shí)均將其近似取為零，于是可得到特征方程為

|λI-A|=λ3+C1λ2+C2λ+C3=0

(17)

式中：λ為特征值;

(18)

根據(jù)Etkin[29]的數(shù)值研究，高超聲速飛行具有振蕩沉浮模態(tài)和單調(diào)縱向螺旋模態(tài)(高度模態(tài))。其特征方程具有一對(duì)共軛復(fù)根和一個(gè)離原點(diǎn)很近的實(shí)根，共軛復(fù)根可表示為

(19)

式中：η為阻尼項(xiàng)；ζ為阻尼比；ω為阻尼振蕩頻率；ωn為無阻尼振蕩頻率，即自然頻率。

特征方程式(17)還可表示為

(λ2+p1λ+p2)(λ-λ3)=0

(20)

式中：

(21)

2.2 阻尼振蕩頻率

特征方程式(17)中常數(shù)項(xiàng)C3的主導(dǎo)項(xiàng)表征阻力的變化，也即切向力變化率項(xiàng)。Etkin[29]關(guān)于火箭動(dòng)力(推力為常值)和吸氣式動(dòng)力(推力與密度成正比)高超聲速飛行器沉浮周期的數(shù)值研究表明，切向力變化對(duì)沉浮周期影響微弱[29-30]，因而可以忽略C3項(xiàng)求解近似沉浮周期[30]，則有λ3=0、p1=C1和p2=C2,于是特征根虛部(即阻尼振蕩頻率)可表示為

(22)

考慮到平衡滑翔條件Z/m+g=V2/r,經(jīng)適當(dāng)變形整理，得沉浮模態(tài)阻尼振蕩頻率解析式為

(23)

式中：

(24)

2.3 阻尼項(xiàng)

阻力是再入飛行器沉浮運(yùn)動(dòng)中消耗能量的唯一來源，因而也是沉浮運(yùn)動(dòng)阻尼的根本來源，故不能通過忽略含阻力相關(guān)項(xiàng)的特征方程得到準(zhǔn)確的阻尼項(xiàng)解析式。下面對(duì)原三階特征方程進(jìn)行推導(dǎo)，將式(20)展開后與式(17)對(duì)比可得

C1=p1-λ3=-2η-λ3=2ζωn-λ3

(25)

(26)

(27)

由式(25)、式(27)可得到沉浮模態(tài)特征根實(shí)部(即阻尼項(xiàng))的解析式為

(28)

(29)

式中：無量綱沉浮修正參數(shù)f的解析式為

(30)

2.4 自然頻率、阻尼比和半衰期內(nèi)振蕩次數(shù)

由式(23)和式(29)可依次得到沉浮模態(tài)自然頻率、阻尼比、半衰期內(nèi)振蕩次數(shù)的解析式為

(31)

(32)

(33)

式中：

(34)

至此，基于平衡滑翔條件和阻力變化量對(duì)周期影響微弱的特點(diǎn)，建立了表征再入沉浮運(yùn)動(dòng)特性的自然頻率、阻尼比兩個(gè)基本參數(shù)的一次近似解析式(31)和式(32)，以及阻尼振蕩頻率、阻尼項(xiàng)和半衰期內(nèi)振蕩次數(shù)的一次近似解析式(23)、式(29)和式(33)。

2.5 高超聲速沉浮近似解析式的低空低速適用性

研究低空低速飛行時(shí)，來流可認(rèn)為是不可壓縮流，地球曲率引起的“離心力”可忽略，大氣密度梯度參數(shù)可取為β=0,低速情況下軌道速度比也可近似取為U≈0。則將高超聲速再入沉浮阻尼振蕩頻率解析式(23)和自然頻率解析式(31)展開，并將β=0和U≈0代入，可推導(dǎo)出低空低速飛行沉浮運(yùn)動(dòng)阻尼振蕩頻率和自然頻率的近似簡(jiǎn)化解析式為

(35)

(36)

由式(30)可知，當(dāng)U≈0時(shí)，有f≈0。于是由式(29)可得到低空低速飛行沉浮運(yùn)動(dòng)阻尼項(xiàng)(特征根實(shí)部)的近似簡(jiǎn)化解析式為

η≈-g/(KVV)

(37)

展開式(32)和式(33)，并將β=0、U≈0和f≈0代入，可得低空低速飛行沉浮運(yùn)動(dòng)阻尼比和半衰期內(nèi)振蕩次數(shù)的近似簡(jiǎn)化解析式為

(38)

(39)

以上低空低速飛行沉浮特性參數(shù)近似簡(jiǎn)化解析式的推導(dǎo)結(jié)果均與Lanchester[22]、Etkin[23]以及方振平[38]等的經(jīng)典結(jié)果一致，這是表明所推導(dǎo)的高超聲速再入沉浮特性近似解析式正確性的必要條件，但并不充分。后續(xù)將在對(duì)解析式二次近似的推導(dǎo)分析過程中，與由狀態(tài)矩陣式(13)得到的數(shù)值解進(jìn)行統(tǒng)一數(shù)值對(duì)比，以充分驗(yàn)證所推導(dǎo)的高超聲速再入沉浮特性近似解析式的正確性。

3 再入沉浮特性參數(shù)二次近似

第2節(jié)推導(dǎo)的一次近似解析式，包含的影響因素關(guān)系仍顯繁瑣，本節(jié)在數(shù)值比較的基礎(chǔ)上，通過忽略高超聲速段的小量進(jìn)行二次近似簡(jiǎn)化。

3.1 簡(jiǎn)化的數(shù)值分析基礎(chǔ)

圖2直觀表示出不同縱向升阻比下W2、W3、W4與W1的相對(duì)大小關(guān)系，結(jié)合式(23)、式(31)和式(32)，可以看出在大部分高超聲速情況下，W1為阻尼振蕩頻率、自然頻率和阻尼比的主導(dǎo)項(xiàng)之一，因而可以忽略W2、W3和W4。

圖3直觀表示出不同縱向升阻比下W6與W5的相對(duì)大小關(guān)系，KV=1.5和KV=100的兩條曲線表明當(dāng)KV>1.5時(shí),W6對(duì)縱向升阻比不敏感，并且結(jié)合式(33)可以看出在大部分高超聲速情況下，W5為表征阻尼特性的半衰期內(nèi)振蕩次數(shù)的主導(dǎo)項(xiàng)之一，因而可忽略W6。

圖4直觀表示出了作為軌道速度比U、大氣密度梯度參數(shù)β和高度r的函數(shù)的沉浮修正參數(shù)f在不同β取值下隨速度的變化情況，可以看出當(dāng)軌道速度比U<0.98,即飛行速度V<0.98Vc時(shí)，有0

圖2 Wi/W1(i=1,2,3,4)隨軌道速度比的變化Fig.2 Wi/W1(i=1,2,3,4) vs the ratio of flight speed to orbital speed

圖3 Wi/W5(i=5,6)隨軌道速度比的變化Fig.3 Wi/W5(i=5,6) vs the ratio of flight speed to orbital speed

圖4 沉浮修正參數(shù)隨軌道速度比的變化Fig.4 Phugoid correction factor vs the ratio of flight speed to orbital speed

3.2 阻尼項(xiàng)二次近似解析式

由3.1節(jié)數(shù)值分析結(jié)果，忽略式(29)中的f,可得到高超聲速再入沉浮阻尼項(xiàng)(特征根實(shí)部)的二次近似簡(jiǎn)化解析式為

η≈-g/(KVV)

(40)

式(40)與低空低速式(37)完全相同，表明與低空低速飛行器一樣，高超聲速飛行器縱向升阻比越大其沉浮阻尼越小，并且高超聲速沉浮阻尼項(xiàng)與速度的反比關(guān)系從形式上決定了弱沉浮阻尼是中高升阻比高超聲速飛行器的固有特性。

不過從本質(zhì)上看，這一固有特性源于高空很小的大氣密度。因?yàn)樽枘嶂饕獊碓从跉鈩?dòng)阻力隨速度的變化，而由其表達(dá)式

dD/dV=?D/?V+(?D/?ρ)(dρ/dr)(dr/dV)≈

?D/?V=ρVSCD

(41)

可知，由于大氣密度隨高度呈指數(shù)減小，使ρV值沿著平衡滑翔軌跡隨著高度和速度增加反而快速減小，從而導(dǎo)致高超聲速飛行器具有了弱阻尼的固有特性。

3.3 自然頻率二次近似解析式

根據(jù)數(shù)值分析結(jié)果，忽略式(23)、式(31)中的W2、W3和W4,則高超聲速再入沉浮阻尼振蕩頻率與自然頻率近似相等(間接表明沉浮阻尼比很小)，再考慮到平衡滑翔條件可得出高超聲速再入沉浮自然頻率的二次近似簡(jiǎn)化解析式為

(42)

圖5中對(duì)自然頻率一次近似解式(31)、二次近似解式(42)、低空低速近似解式(36)和由狀態(tài)矩陣式(13)得到的數(shù)值解進(jìn)行了對(duì)比，驗(yàn)證了所推導(dǎo)解析式的正確性，并表明二次近似解在高超聲速飛行段(U>0.2)的精度很好。由式(42)可知高超聲速再入沉浮運(yùn)動(dòng)的自然頻率主要由軌道速度比U和與飛行器特性無關(guān)的大氣密度梯度參數(shù)β、重力加速度g決定，因而不同飛行器高超聲速沉浮自然頻率隨軌道速度比的變化情況相差不大。

由式(42)可得到沉浮運(yùn)動(dòng)等效彈性系數(shù)為

(43)

式(43)表明高超聲速再入沉浮振蕩的恢復(fù)力主要來源于隨高度變化的縱向平面升力變化量。

圖5 沉浮自然頻率隨軌道速度比的變化Fig.5 Phugoid natural frequency vs the ratio of flight speed to orbital speed

3.4 阻尼比和半衰期振蕩次數(shù)二次近似解析式

依據(jù)簡(jiǎn)化的數(shù)值分析結(jié)果，忽略式(32)中f、W2和W4,以及式(33)中的f和W6,可分別得到高超聲速再入沉浮阻尼比和半衰期內(nèi)振蕩次數(shù)的二次近似簡(jiǎn)化解析式為

(44)

(45)

圖6中將阻尼比一次近似解式(32)、二次近似解式(44)、低空低速近似解式(38)、由狀態(tài)矩陣式(13)得到的數(shù)值解和由胡錦川和陳萬春[21]通過高度偏差的二階微分方程推導(dǎo)出的近似解進(jìn)行了對(duì)比，驗(yàn)證了所推導(dǎo)解析式的正確性，且可以看出本文的阻尼比二次近似解式(44)在除軌道速度外的高超聲速飛行段(0.2

二次近似解式(44)中，由于-βr≈900,當(dāng)0.2≤U≤0.98時(shí),有0.066 7/KV≤ζ≤0.17/KV,可見縱向升阻比KV和軌道速度比U是沉浮阻尼特性的主導(dǎo)影響參數(shù)，同時(shí)表明低阻尼比是中高升阻比高超聲速飛行器再入沉浮運(yùn)動(dòng)的固有特性，且縱向升阻比越大阻尼比越小。

綜合式(40)和式(44)的分析可知，再入飛行器在高超聲速滑翔時(shí)具有弱阻尼和低阻尼比的固有特性，且作為飛行器本體參數(shù)的升阻比對(duì)其也具有重要的主導(dǎo)作用(與頻率不同)，同時(shí)其反比關(guān)系也解釋了隨升阻比增大高超聲速飛行器越易產(chǎn)生沉浮振蕩的原因。

圖6 沉浮阻尼比隨軌道速度比的變化Fig.6 Phugoid damping ratio vs the ratio of flight speed to orbital speed

4 再入沉浮特性近似解的應(yīng)用

本節(jié)對(duì)以上所推導(dǎo)的沉浮特性近似解及近似關(guān)系分別在沉浮特性規(guī)律的解釋和軌跡振蕩的抑制方面進(jìn)行了應(yīng)用。

4.1 高超聲速沉浮阻尼特性規(guī)律解釋

表征高超聲速沉浮阻尼特性的半衰期內(nèi)振蕩次數(shù)N1/2隨高度變化的曲線(圖7圓圈)和大氣密度梯度參數(shù)β隨高度變化的曲線(圖1實(shí)線)具有相似的變化規(guī)律，很容易讓人產(chǎn)生兩者具有強(qiáng)相關(guān)性的認(rèn)識(shí)。最早進(jìn)行高超聲速飛行動(dòng)力學(xué)數(shù)值研究的Etkin認(rèn)為，高超聲速飛行器在90多千米高度處出現(xiàn)“重阻尼”(N1/2出現(xiàn)極小值)是由于大氣密度梯度參數(shù)β的絕對(duì)值在該高度出現(xiàn)峰值導(dǎo)致的[29]。但通過解析分析和數(shù)值驗(yàn)證，均表明此觀點(diǎn)不準(zhǔn)確。以下進(jìn)行詳細(xì)說明，并應(yīng)用所推導(dǎo)的近似解分析其變化規(guī)律的主導(dǎo)因素。

首先，由所推導(dǎo)的高超聲速沉浮運(yùn)動(dòng)近似式(45)可以看出，β絕對(duì)值越大，N1/2越大(即沉浮阻尼越弱)，而不是Etkin推測(cè)的阻尼越重。

其次，為便于同Etkin的數(shù)值結(jié)果進(jìn)行對(duì)比，圖7選用了與Etkin相同的飛行器模型。圖中圓圈為通過狀態(tài)矩陣式(13)得到的數(shù)值解，其他各曲線則依據(jù)一次近似解式(33)得到?？梢钥闯?，當(dāng)β取常值-1/7 250時(shí)，N1/2同樣會(huì)在90多千米高度出現(xiàn)極小值，這足以說明N1/2在90多千米高度出現(xiàn)極小值與β的絕對(duì)值在該高度有極大值不相關(guān)或相關(guān)性很低。

為了進(jìn)一步確定高超聲速沉浮阻尼特性變化規(guī)律的主要影響因素，圖8給出了依據(jù)一次近似解式(33)得到的與圖7中曲線相對(duì)應(yīng)的沉浮運(yùn)動(dòng)半衰期振蕩次數(shù)N1/2隨軌道速度比U的變化曲線。從圖中可看出每條曲線均在U≈0.707附近出現(xiàn)峰值，在接近U=1附近出現(xiàn)谷值，這與圖7中的峰值和谷值相對(duì)應(yīng)。

N1/2在U≈0.707附近出現(xiàn)峰值的原因可以由一次近似解式(33)中的主導(dǎo)項(xiàng)W5中的軌道速度比項(xiàng)解釋，其有如下關(guān)系：

(1-U2)U2≤[(1-U2+U2)/2]2=0.25

(46)

式中：當(dāng)1-U2=U2,即U2=0.5或U≈0.707時(shí)，等號(hào)成立，可取得最大值。

N1/2在接近U=1附近出現(xiàn)谷值的原因可由一次近似解式(33)中的1/(1-f)項(xiàng)來解釋。由圖4可看出，當(dāng)U接近于1時(shí)，沉浮修正參數(shù)f會(huì)突然由0快速趨近于1，進(jìn)而使得1/(1-f)項(xiàng)急劇增大，從而改變之前由(1-U2)U2項(xiàng)主導(dǎo)的N1/2隨U增加而遞減的趨勢(shì)。于是在接近U=1時(shí)N1/2出現(xiàn)了谷值，并且出現(xiàn)谷值的高度必然在f由0快速趨近于1的高度區(qū)間內(nèi)。圖9給出了依據(jù)式(30)得到的沉浮修正參數(shù)f隨高度的變化曲線，可以看出圖7中不同大氣密度梯度參數(shù)下N1/2出現(xiàn)谷值的高度均在圖9中f由0快速趨近于1的高度區(qū)間內(nèi)。

圖7 沉浮運(yùn)動(dòng)半衰期內(nèi)振蕩次數(shù)隨高度的變化Fig.7 Cycles to half amplitude of phugoid oscillations vs altitude

圖8 沉浮運(yùn)動(dòng)半衰期內(nèi)振蕩次數(shù)隨軌道速度比的變化Fig.8 Cycles to half amplitude of phugoid oscillations vs the ratio of flight speed to orbital speed

圖9 沉浮修正參數(shù)隨高度的變化Fig.9 Phugoid correction factor vs altitude

綜上分析，高超聲速沉浮阻尼特性的變化規(guī)律主要由軌道速度比U和沉浮修正參數(shù)f主導(dǎo)，而與大氣密度梯度參數(shù)β的相關(guān)性較低。其中軌道速度比U的影響關(guān)系式(1-U2)U2表征縱向升力與“離心力”的乘積，且當(dāng)縱向升力與“離心力”相等，即軌道速度比U為0.707時(shí)，沉浮阻尼比取極小值；而沉浮修正參數(shù)f則對(duì)接近軌道速度時(shí)沉浮阻尼特性的急劇變化具有主導(dǎo)作用。

4.2 再入沉浮振蕩軌跡的抑制

由前述可知高超聲速滑翔再入飛行具有弱阻尼和低阻尼比的固有特性，因而當(dāng)初始再入條件偏離平衡滑翔條件過多或各種不確定性過大導(dǎo)致預(yù)測(cè)控制誤差較大時(shí)，飛行器容易產(chǎn)生縱向沉浮振蕩，故需要設(shè)計(jì)相應(yīng)的控制器進(jìn)行軌跡振蕩抑制。本節(jié)對(duì)抑制振蕩的控制關(guān)系進(jìn)行了推導(dǎo)。

采用以下反饋控制結(jié)構(gòu)調(diào)節(jié)縱向沉浮運(yùn)動(dòng)阻尼比，進(jìn)而對(duì)再入縱向軌跡振蕩進(jìn)行抑制。

Δu=Δv-FΔx

(47)

式中：Δv=[Δkσ]為閉環(huán)反饋系統(tǒng)的控制變量；F=[kVkγkr]為反饋增益矩陣。將式(47)代入再入飛行器縱向動(dòng)態(tài)方程式(11)，可得閉環(huán)反饋系統(tǒng)狀態(tài)矩陣為

(48)

對(duì)比在不同速度下單獨(dú)反饋速度、航跡角、高度的根軌跡，如圖10所示，可知反饋航跡角γ時(shí)具有最有效的阻尼比調(diào)節(jié)效果，且反饋增益需隨速度增加而增加。進(jìn)一步由圖11中不同傾側(cè)角σ下反饋航跡角γ的根軌跡圖可以看出，若要系統(tǒng)保持在最佳阻尼比附近，反饋增益需隨傾側(cè)角增加而減小。

下面推導(dǎo)航跡角反饋增益的近似解析表達(dá)式。令kV=kr=0,則航跡角反饋系統(tǒng)的閉環(huán)特征方程可表示為

(49)

式中：

(50)

(51)

考慮到高超聲速飛行縱向螺旋模態(tài)的根離S平面原點(diǎn)很近，則由式(25)和式(26)可得以下近似關(guān)系：

(52)

(53)

分別將式(50)和式(51)代入式(52)和式(53)，并根據(jù)式(23)的推導(dǎo)和圖2中所示出的W1為主導(dǎo)項(xiàng)的特性，可得到航跡角反饋增益近似表達(dá)式為

(54)

圖10 反饋不同參數(shù)下的根軌跡圖Fig.10 Root locus for different parameters feedback

圖11 反饋航跡角的根軌跡圖Fig.11 Root locus for flight path angle feedback

圖12 再入狀態(tài)量與控制量曲線圖Fig.12 Reentry state parameters and control parameters curves

圖13 速度-高度空間中的再入走廊和再入軌跡Fig.13 Reentry corridor and trajectories in velocity-altitude space

為了對(duì)比軌跡振蕩抑制器的效果，設(shè)置初始再入條件使得基準(zhǔn)制導(dǎo)律再入軌跡產(chǎn)生振蕩，如圖12和圖13中虛線所示；而嵌入軌跡振蕩抑制器后的仿真結(jié)果如圖中實(shí)線所示，可以看出所設(shè)計(jì)的軌跡振蕩抑制器在未過多增加控制負(fù)擔(dān)的情況下有效抑制了再入軌跡的振蕩。

5 結(jié) 論

1) 所推導(dǎo)出的高超聲速再入沉浮特性近似解析解簡(jiǎn)潔直觀地表明：自然頻率隨軌道速度比U的變化與飛行器本體參數(shù)關(guān)系不大；高超聲速再入沉浮運(yùn)動(dòng)具有弱阻尼和低阻尼比的固有特性，且與縱向升阻比KV成反比關(guān)系，故升阻比越大的飛行器越易產(chǎn)生再入沉浮振蕩。

2) 驗(yàn)證了高超聲速沉浮阻尼特性隨高度的變化規(guī)律的主要影響因素是軌道速度比U和所推導(dǎo)的沉浮修正參數(shù)f,而不是早期研究者推測(cè)的大氣密度梯度參數(shù)β,其中軌道速度比的影響關(guān)系式表征縱向升力與“離心力”的乘積，且當(dāng)縱向升力與“離心力”相等，即軌道速度比U為0.707時(shí)，沉浮阻尼比取極小值。

3) 建立了再入軌跡振蕩抑制的閉環(huán)反饋系統(tǒng)近似關(guān)系，彌補(bǔ)了依靠經(jīng)驗(yàn)確定參數(shù)的不足，從而進(jìn)一步完善了平衡滑翔預(yù)測(cè)-校正再入制導(dǎo)方法。數(shù)值仿真表明，所采用的近似控制關(guān)系能夠有效抑制由初始偏差等引起的再入沉浮振蕩。

[1] JOHNSON W R, LU P, STACHOWIAK S J. Automated re-entry system using FNPEG: AIAA-2017-1899[R]. Reston， VA: AIAA, 2017.

[2] LU P, BRUNNER C W, STACHOWIAK S J, et al. Verification of a fully numerical entry guidance algorithm: AIAA-2016-0377[R]. Reston， VA: AIAA, 2016.

[3] LU P. Entry guidance: A unified method[J]. Journal of Guidance, Control, and Dynamics, 2014, 37(3): 713-728.

[4] LU P, FORBES S, BALDWIN M. Gliding guidance of high L/D hypersonic vehicles: AIAA-2013-4648[R]. Reston， VA: AIAA, 2013.

[5] WEBB K, LU P. Entry guidance by onboard trajectory planning and tracking: AIAA-2016-0279[R]. Reston, VA: AIAA, 2016.

[6] PUTNAM Z R, GRANT M J, KELLY J R, et al. Feasibility of guided entry for a crewed lifting body without Angle-of-Attack control[J]. Journal of Guidance, Control, and Dynamics, 2014, 37(3): 729-740.

[7] XUE S B, LU P. Constrained predictor-corrector entry guidance[J]. Journal of Guidance, Control, and Dynamics, 2010, 33(4): 1273-1281.

[8] JOSHI A, SIVAN K, AMMA S S. Predictor-corrector reentry guidance algorithm with path constraints for atmospheric entry vehicles[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(5): 1307-1318.

[9] LU P. Asymptotic analysis of quasi-equilibrium glide in lifting entry flight[J]. Journal of Guidance, Control, and Dynamics, 2006, 29(3): 662-670.

[10] SHEN Z J, LU P. Onboard generation of Three-Dimensional constrained entry trajectories[J]. Journal of Guidance, Control, and Dynamics, 2003, 26(1): 111-121.

[11] SHEN Z J. On-board three-dimensional constrained entry flight trajectory generation[D]. Ames, Iowa: Iowa State University, 2002: 19-64.

[12] HARPOLD J C, GRAVES C A. Shuttle entry guidance[J]. Journal of Astronautical Sciences, 1979, 37(3): 239-268.

[13] 胡軍, 張釗. 載人登月飛行器高速返回再入制導(dǎo)技術(shù)研究[J]. 控制理論與應(yīng)用, 2014, 31(12): 1678-1685.

HU J, ZHANG Z. A study on the reentry guidance for a manned lunar return vehicle[J]. Control Theory & Applications, 2014, 31(12): 1678-1685 (in Chinese).

[14] 崔乃剛, 黃榮, 傅瑜, 等. 基于匹配漸進(jìn)展開的跳躍式再入解析預(yù)測(cè)-校正制導(dǎo)律設(shè)計(jì)[J]. 航空學(xué)報(bào), 2015, 36(8): 2764-2772.

CUI N G, HUANG R, FU Y, et al. Design of analytical prediction-correction skip entry guidance law based on matched asymptotic expansions[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(8): 2764-2772 (in Chinese).

[15] 楊勇, 張輝, 鄭宏濤. 有翼高超聲速再入飛行器氣動(dòng)設(shè)計(jì)難點(diǎn)問題[J]. 航空學(xué)報(bào), 2015, 36(1): 49-57.

YANG Y, ZHANG H, ZHENG H T. Difficulties in aerodynamic design problems of the winged hypersonic reentry vehicle[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 49-57 (in Chinese).

[16] JORRIS T R. Common aero vehicle autonomous reentry trajectory optimization satisfying waypoint and no-fly zone constraints[D]. Wright-Patterson Air Force Base: Air University, 2007: 34-41, 89.

[17] GUO J, WU X Z, TANG S J. Autonomous gliding entry guidance with geographic constraints[J]. Chinese Journal of Aeronautics, 2015, 28(5): 1343-1354.

[18] ZHAO J, ZHOU R. Particle swarm optimization applied to hypersonic reentry trajectories[J]. Chinese Journal of Aeronautics, 2015, 28(3): 822-831.

[19] MALL K, GRANT M J. High mass Mars exploration using slender entry vehicles: AIAA-2016-0019[R]. Reston, VA: AIAA, 2016.

[20] CHEN X Q, HOU Z X, LIU J X, et al. Phugoid dynamic characteristic of hypersonic gliding vehicles[J]. Science China (Information Sciences), 2011, 54(3): 542-550.

[21] 胡錦川, 陳萬春. 高超聲速飛行器平穩(wěn)滑翔彈道設(shè)計(jì)方法[J]. 北京航空航天大學(xué)學(xué)報(bào), 2015, 41(8): 1464-1475.

HU J C, CHEN W C. Steady glide trajectory planning method for hypersonic reentry vehicle[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(8): 1464-1475 (in Chinese).

[22] LANCHESTER F W. Aerodonetics[M]. London: Archibald Constable & Co.Ltd., 1908: 18-82.

[23] ETKIN B. Dynamics of atmospheric flight[M]. New York: John Wiley & Sons, Inc., 1972: 323, 329-332, 337-339.

[24] BAIRSTOW L. Applied aerodynamics[M]. London: Longmans, Green and Co., 1939: 447-551.

[25] ASHKENAS I L, MCRUER D T. Approximate airframe transfer functions and application to single sensor control systems: WADC-TR-58-82[R]. Wright-Patterson Air Force Base: Wright Air Development Center, 1958.

[26] 徐瑞娟. 對(duì)飛機(jī)縱向沉浮模態(tài)解析法的探討[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 1995, 13(4): 457-462.

XU R J. An exploration on the analytical-solution for the longitudinal phugoid mode of aircraft[J]. Acta Aerodynamica Sinica, 1995, 13(4): 457-462 (in Chinese).

[27] SCHEUBEL F N. Der einfluss des dichtegradienten der atmosphere auf die langsbewegung des flugzeugs[J]. Luftfahrtforschung, 1942, 19(4): 132-136.

SCHEUBEL F N. The influence of the density gradient of the atmosphere on the slow motion of the airplane[J]. Aviation Research, 1942, 19(4): 132-136 (in German).

[28] NEUMARK S. Longitudinal stability, speed and height: an examination of dynamic longitudinal stability in level flight, including the effects of compressibility and changes in atmospheric phenomena with height[J]. Aircraft Engineering and Aerospace Technology, 1950, 22(11): 323-334.

[29] ETKIN B. Longitudinal dynamics of a lifting vehicle in orbital flight[J]. Journal of the Aerospace Sciences, 1961, 28(28): 779-788, 832.

[30] LAITONE E V, CHOU Y S. Phugoid oscillations at hypersonic speeds[J]. AIAA Journal, 1965, 3(4): 732-735.

[31] VINH N X, DOBRZELECKI A J. Nonlinear longitudinal dynamics of an orbital lifting vehicle[J]. Celestial Mechanics, 1971, 3(4): 427-460.

[32] BERRY D T. Longitudinal long-period dynamics of aerospace craft: AIAA-1988-4358[R]. Reston, VA: AIAA, 1988.

[33] MARKOPOULOS N, MEASE K D, VINH N X. Thrust law effects on the long-period modes of aerospace craft: AIAA-1989-3379[R]. Reston, VA: AIAA, 1989.

[34] FERREIRA L D O. Nonlinear dynamics and stability of hypersonic reentry vehicles[D]. Michigan: University of Michigan, 1995: 9-82, 138-139.

[35] REGAN F J, ANANDAKRISHNAN S M. Dynamics of atmospheric reentry[M]. Reston, VA: AIAA, 1993: 37-39, 56.

[36] BETTS J T. Practical methods for optimal control using nonlinear programming[M]. Philadelphia: Society for Industrial and Applied Mathematics, 2001: 133-134.

[37] SLOTINE J E, LI W P. 應(yīng)用非線性控制[M]. 程代展, 譯. 北京: 機(jī)械工業(yè)出版社, 2006: 28-38.

SLOTINE J E, LI W P. Applied nonlinear control[M]. CHENG D Z, translated. Beijing: China Machine Press, 2006: 28-38 (in Chinese).

[38] 方振平, 陳萬春, 張曙光. 航空飛行器飛行動(dòng)力學(xué)[M]. 北京: 北京航空航天大學(xué)出版社, 2005: 186-204, 299.

FANG Z P, CHEN W C, ZHANG S G. Aircraft flight dynamics[M]. Beijing: Beihang University Press, 2005: 186-204, 299 (in Chinese).

Approximateanalyticalanalysisforphugoidcharacteristicofreentryvehiclesanditsapplications

GUJie1,ZHANGShuguang1,*,YANGFan2,WANGBaoyin1

1.SchoolofTransportationScienceandEngineering,BeihangUniversity,Beijing100083,China2.ChengduAircraftDesignandResearchInstitute,Chengdu610091,China

Forthehypersonicreentryvehicleadoptingtheequilibriumglide-basednumericaloranalyticalpredictor-correctorguidance,thephugoidoscillationofthereentrytrajectoryispronetobecausedatthetransitionpointfromtheinitialdescenttotheequilibriumglideorbylargepredictiondeviations.Withtheincreaseofthelift-to-dragratioofreentryvehiclesunderresearchinrecentyears,thephugoidoscillationbecomesmorenoticeabletocauseresearcherstore-examinethehypersonicphugoidcharacteristic.Aconciseapproximateanalyticalsolutionforintuitivelyexpressingthedominatefactorsforhypersonicphugoidcharacteristicofreentryvehiclesisderived,accordingtothreeorderlongitudinalflightdynamicsequationsandtheequilibriumglidecondition.Itisfoundthatthehypersonicphugoiddampingcharacteristicmainlydependsontheratioofflightspeedtoorbitalspeedandthephugoidcorrectionfactor,ratherthanontheatmosphericdensity-gradientparameterconsideredbypreviousresearchers.Acontrollerisdesignedtoeliminatethetrajectoryoscillationbyusingtheapproximaterelationderived,furtherimprovingtheequilibriumglidebasednumericaloranalyticalpredictor-correctorguidance.Effectivenessofthemethodisverifiedbythesimulationresults.

reentry;hypersonicvehicle;phugoid;damping;analyticalsolution

2017-02-15；Revised2017-04-06；Accepted2017-05-09;Publishedonline2017-05-121059

URL：http://hkxb.buaa.edu.cn/CN/html/20171004.html

.E-mailgnahz@buaa.edu.cn

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2017.121174

V212.1

A

1000-6893(2017)10-121174-12

2017-02-15；退修日期2017-04-06；錄用日期2017-05-09；< class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2017-05-121059

http://hkxb.buaa.edu.cn/CN/html/20171004.html

.E-mailgnahz@buaa.edu.cn

顧杰，張曙光，楊帆，等．再入飛行器沉浮特性近似解析及應(yīng)用J．航空學(xué)報(bào),2017,38(10):121174.GUJ,ZHANGSG,YANGF,etal.ApproximateanalyticalanalysisforphugoidcharacteristicofreentryvehiclesanditsapplicationsJ.ActaAeronauticaetAstronauticaSinica,2017,38(10):121174.

(責(zé)任編輯：李明敏)