叢戎飛，葉友達，趙忠良

1. 中國空氣動力研究與發(fā)展中心 高速空氣動力研究所，綿陽 621000

2. 中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000

3. 國家計算流體力學實驗室，北京 100083

2004年X-43A實驗飛行器的成功試飛標志著吸氣式高超聲速飛行器已經(jīng)從實驗室階段走向工程研制階段[1-4]。2013年，美國洛克希德·馬丁公司宣布正在研制代號為“黑燕”的SR-72高超聲速無人偵察機作為SR-71的繼承機型[5]。目前以美俄為代表的世界軍事強國正在積極推動吸氣式高超聲速飛行器在軍事領域的應用。

由于高超聲速飛行器的飛行環(huán)境以及機體構型與傳統(tǒng)飛行器不同，在穩(wěn)定性研究中面臨著很多新的問題。高超聲速飛行器進行機動運動時，其氣動特性通常是非線性、非定常變化的，并存在氣動/運動耦合以及慣性耦合作用。由于高超聲速飛行器通常采用細長體構型，其繞體軸X方向的轉動慣量Ix通常較小，且高空高速飛行中的氣動阻尼較小，使得飛行器在受到擾動時容易發(fā)生非指令性滾轉[6-7]。

X-43A在2001年首次試飛時，在機彈分離后出現(xiàn)了滾轉振蕩發(fā)散現(xiàn)象，過大的載荷導致氣動舵失效，從而導致試飛失敗。研究人員分析認為是動力學耦合造成飛行器失穩(wěn)，并將問題指向了氣動模型的不可靠性[8]。在2010年4月和2011年8月，美國兩次進行HTV-2的飛行試驗，均以失敗告終。其中第1次試飛失敗是由于在飛行器爬升階段發(fā)生滾轉與側滑交叉耦合，從而引起滾轉發(fā)散，超過飛行器的滾轉控制極限，最終導致飛行器失控[9]。由此可見，高超聲速飛行器的動穩(wěn)定性問題十分突出，給飛行器的運動穩(wěn)定性和飛行安全帶來了嚴峻的挑戰(zhàn)[10-11]，因此高超聲速飛行器橫側向穩(wěn)定性及其動力學耦合是必須研究的課題。

國內(nèi)外針對傳統(tǒng)飛行器的動穩(wěn)定性開展的研究較早。Phillips[12]在研究XS-1模型時發(fā)現(xiàn)了由飛行器滾轉引起的慣性耦合現(xiàn)象，從而提出了用于研究急滾動力學的Phillips判據(jù)。20世紀70年代，Mehra和Carroll[13]將分叉分析方法應用于飛行器穩(wěn)定性分析，并指出機翼搖滾是由Hopf分叉引起的。Kandil[14]和劉偉[15]等分別開展了機翼單自由度滾動的非線性動力學分析研究，楊小亮[16]開展了后掠三角翼強迫俯仰、自由滾轉耦合運動特性研究，張涵信等[17]針對球冠倒錐外形返回艙在跨、超聲速配平攻角區(qū)出現(xiàn)極限環(huán)振蕩的現(xiàn)象，提出了動穩(wěn)定性判據(jù)及實現(xiàn)分叉的條件和分叉臨界馬赫數(shù)。袁先旭等采用耦合求解非定常Navier-Stokes方程和飛行器運動方程，數(shù)值模擬了飛船返回艙俯仰振蕩隨來流馬赫數(shù)變化的Hopf分叉過程，驗證了分析結論[18]。

針對高超聲速飛行器的動穩(wěn)定性研究則起步較晚，作者團隊[10-11]在張涵信已發(fā)展的正滾判則的基礎上，研究了類HTV-2高超飛行器的單自由度滾轉以及強迫俯仰/滾轉運動，完善了飛行器偏滾運動的穩(wěn)定性判則，并指出高超聲速飛行器高空飛行時，由于大氣的密度低，飛行器自身的慣性力矩不能被忽略，慣性力矩會迅速加快滾動速度，誘導偏離振蕩及慣性耦合發(fā)散。李乾等[19-20]對類HTV-2高超飛行器在兩自由度耦合運動條件下的氣動/運動耦合現(xiàn)象進行了研究，發(fā)現(xiàn)飛行器在俯仰機動過程中，滾轉方向會在氣動力矩以及慣性力矩共同作用下發(fā)生失穩(wěn)運動。

目前針對吸氣式高超聲速飛行器動態(tài)特性開展的研究較少，以單自由度動穩(wěn)定性研究為主，較少涉及多自由度耦合運動研究。Adamczak和Bolender[21]對HIFiRE-6一體化飛行器的縱向以及橫航向運動模態(tài)進行了分析，指出飛行器在小迎角下橫側向通道動不穩(wěn)定。劉緒等[22-23]等通過數(shù)值模擬針對類X-51機體/推進一體化飛行器在給出了進氣道堵塞和通流兩種外形下小振幅強迫振動的俯仰、偏航、滾轉阻尼導數(shù)和交叉導數(shù)，并對兩種外形下的計算結果進行了分析。趙忠良等[24-25]在0.5 m高超聲速風洞針對類X-51A外形的帶進氣道通流模型進行了動導數(shù)試驗，研究了模型進氣道喉道高度對吸氣式高超聲速飛行器動導數(shù)的影響。

對吸氣式飛行器而言，由于機體外形以及氣動特性與滑翔式飛行器有較大區(qū)別，其動穩(wěn)定性值得進行進一步的研究。本文針對一種類似SR-72構型的吸氣式高超聲速飛機開展了動穩(wěn)定性研究。通過數(shù)值模擬獲得了滾轉單自由度動穩(wěn)定性以及強迫俯仰/自由滾轉運動下的兩自由度耦合動穩(wěn)定性。

1 數(shù)值方法和計算模型

本文的數(shù)值模擬研究在課題組研發(fā)的軟件VFNS[26]上完成,該軟件采用并行結構網(wǎng)格求解器，求解非定常Navier-Stokes方程，采用高階Adams預估校正法在雙時間步內(nèi)實現(xiàn)空氣動力/飛行力學的緊耦合計算,控制方程空間離散采用的是基于多塊拼接網(wǎng)格的有限體積法，黏性項采用中心差分，無黏項采用Van Leer格式，非定常時間推進采用雙時間步法，湍流模型為Spalart-Allmaras模型。

本文研究的吸氣式高超聲速飛行器模型是參照美國SR-72[27-28]概念飛行器設計的，采用了翼身組合體氣動布局，與HTV-2和X-51A等高超聲速飛行器普遍采用的“乘波體”構型有較大的不同。圖1和圖2分別給出了模型幾何外形以及內(nèi)流道形狀。計算網(wǎng)格拓撲如圖3所示。

為了模擬吸氣式高超聲速飛行器典型的巡航狀態(tài)，計算工況選取為高度為27 km的標準大氣參數(shù)，馬赫數(shù)Ma=6。計算參考長度為模型全長，參考面積為機翼面積，力矩參考點為1/2全長的位置。以模型全長為參考長度的雷諾數(shù)Re=7.1×107。模型內(nèi)流道為通流狀態(tài)，忽略發(fā)動機燃燒等情況。

圖1 模型幾何外形

圖2 內(nèi)流道形狀

圖3 計算網(wǎng)格拓撲

2 驗證算例

本節(jié)選取Basic Finner導彈標模，利用數(shù)值強迫振蕩的方法，計算導彈的俯仰、滾轉動導數(shù)，并與試驗結果對比，驗證程序對于氣動/運動耦合狀態(tài)下非定常氣動力的計算能力。

圖4 Basic Finner俯仰動導數(shù)

圖5 Basic Finner滾轉動導數(shù)

3 滾轉穩(wěn)定性

3.1 靜穩(wěn)定性

圖6 不同俯仰角下滾轉力矩系數(shù)隨滾轉角變化

圖7 滾轉靜導數(shù)

3.2 動穩(wěn)定性

由上文可知模型在0°～15°俯仰角下，在一定滾轉角范圍內(nèi)具備滾轉靜穩(wěn)定性，即模型在平衡狀態(tài)下受到滾轉方向的瞬時擾動時具有恢復到平衡狀態(tài)的趨勢。而模型在受到擾動后的恢復過程以及在平衡位置的斂散特性則是由動穩(wěn)定性描述的。下面使用數(shù)值強迫滾轉的方式求得不同俯仰角下模型的滾轉動導數(shù)，從而評估模型的動穩(wěn)定性。

強迫滾轉振蕩運動的滾轉力矩系數(shù)遲滯曲線如圖8所示，遲滯曲線的傾斜程度與圖6中靜態(tài)滾轉力矩系數(shù)曲線相對應。遲滯曲線的方向如圖中箭頭所示，為逆時針方向，說明氣動力起阻尼作用，模型具備滾轉動穩(wěn)定性。滾轉動導數(shù)為圖9所示，在0°～15°俯仰角下，滾轉動導數(shù)均小于零，同樣說明模型具備滾轉動穩(wěn)定性。其中在俯仰角為5°時滾轉動導數(shù)的絕對值最小，說明此時滾轉運動受到的阻尼最小，滾轉動穩(wěn)定性最差。

圖8 滾轉力矩系數(shù)遲滯曲線

圖9 滾轉動導數(shù)

為了與強迫滾轉運動的計算結果進行對比驗證，數(shù)值模擬了模型在不同俯仰角下的自由滾轉運動，模型分別在5°、10°、15°俯仰角下，從5°初始滾轉角釋放，在氣動滾轉力矩的作用下自由滾轉，逐漸收斂至平衡位置。模型的滾轉方向轉動慣量Ix為0.1×106kg·m2。數(shù)值模擬得到的自由滾轉運動曲線如圖10所示，運動相圖如圖11所示。ωx為滾轉角速度。隨著俯仰角增大，滾轉收斂速度變快，滾轉頻率增大。滾轉運動頻譜圖如圖12所示。由圖可知θ=5°時，滾轉振蕩頻率f約為0.134 Hz；θ=10°時，滾轉振蕩頻率f約為0.213 Hz；θ=15°時，滾轉振蕩頻率f約為0.372 Hz。

圖10 自由滾轉運動時滾轉角時間歷程曲線

圖11 自由滾轉運動相圖

圖12 滾轉運動頻譜圖

4 強迫俯仰/自由滾轉兩自由度耦合運動

在第3節(jié)單自由度滾轉穩(wěn)定性研究的基礎上，通過數(shù)值模擬強迫俯仰/自由滾轉運動，進一步研究吸氣式高超聲速飛行器的氣動/運動耦合作用機理，考察在飛行器俯仰運動的過程中滾轉通道的穩(wěn)定性。

4.1 轉動慣量影響

為了研究飛行器的轉動慣量對強迫俯仰/自由滾轉運動帶來的影響，參照相似尺寸飛行器的轉動慣量分布，針對本文的模型設計了如表1所示的轉動慣量工況。對表中工況進行了存在初始滾轉角時的強迫俯仰/自由滾轉運動數(shù)值模擬。模型初始滾轉角為5°，強迫俯仰運動為正弦振蕩，平均迎角為5°，振幅為5°，頻率為0.025 Hz。

首先保持Ix不變，研究Iy、Iz變化帶來的影響。由于高超聲速飛行器通常采用細長體構型，其滾轉方向轉動慣量Ix通常遠小于偏航和俯仰方向轉動慣量Iy、Iz。在飛行器初始滾轉角不為0°時，俯仰方向的運動引起的慣性力矩在滾轉方向上的分量，也就是慣性耦合力矩可能會對滾轉方向產(chǎn)生擾動。由飛行器運動學方程和動力學方程得到慣性耦合力矩M′x的表達式為

表1 模型轉動慣量

(1)

針對Case1、Case3，分析不同工況下飛行器滾轉角隨時間的變化，如圖13所示。圖中黑色虛線為俯仰角，彩色曲線為不同工況下的滾轉角變化曲線。圖中3種工況的曲線全部重疊在一起。這是由于俯仰角速度較小，慣性耦合力矩很小，慣性耦合作用對滾轉運動的影響可以忽略不計。以Case1工況為例，慣性耦合力矩M′x(黑色曲線)和氣動滾轉力矩Mx(紅色曲線)隨時間的變化曲線如圖14所示，可以看到Mx比M′x大了3個數(shù)量級。這說明滾轉方向的運動是由氣動力主導的。

圖13 滾轉角和俯仰角時間歷程曲線(Cases1，2，3)

圖14 慣性耦合力矩與氣動滾轉力矩時間歷程曲線

在圖13中可以看到，隨著俯仰角運動區(qū)間的不同，滾轉角的振動模態(tài)呈現(xiàn)擬周期性的變化。如果把俯仰角大于平均振幅的周期稱為上半周期，如t=40～60 s(圖中藍背景區(qū)域)，小于平均振幅的稱為下半周期，如t=60～80 s(圖中黃背景區(qū)域)，可以發(fā)現(xiàn)，除去前20 s,在俯仰角上半周期存在滾轉角約為15°的中等振幅振蕩，頻率較高，而在下半周期，滾轉角以超過60°振幅進行大振幅振蕩，且振動頻率與俯仰角振蕩基本頻率一致。

盡管上文對滾轉穩(wěn)定性的研究表明模型在滾轉方向具有靜穩(wěn)定性和動穩(wěn)定性，但在強迫俯仰/自由滾轉耦合運動的過程中，模型在滾轉方向出現(xiàn)了中等振幅的振蕩，并在θ=0°的位置附近出現(xiàn)了滾轉角約70°的振蕩。這說明單自由度運動的動導數(shù)不能對兩自由度耦合運動進行準確預測。

葉友達、田浩等[10-11]在研究類HTV-2高超聲速飛行器兩自由度耦合運動時提出了兩自由度穩(wěn)定性判據(jù):

(2)

圖15 判據(jù)P和滾轉角時間歷程曲線

接下來保持Iy、Iz不變，研究Ix變化帶來的影響。針對Case2、4、5，分析不同工況下飛行器滾轉角隨時間的變化如圖16所示?？梢钥闯鯥x變化對滾轉角的運動規(guī)律影響較為明顯。在俯仰角上半周期，滾轉角振動頻率隨Ix的增大而減小，滾轉角的振幅隨著Ix的增大而 明顯增大，Case5(Ix=0.05×106kg·m2)最大振幅只有13°，Case4(Ix=0.2×106kg·m2)最大振幅達到40°。而在俯仰角下半周期，滾轉角的振幅也隨著Ix的增大而增大，但變化相對較小。

圖16 滾轉角和俯仰角時間歷程曲線(Cases2，4，5)

4.2 俯仰運動頻率影響

下面針對俯仰運動頻率的影響展開進一步的研究。保持模型初始滾轉角5°，強迫俯仰運動為正弦振蕩，平均迎角5°，振幅5°，模型慣量分布與Case2保持一致，分別對f=0.012 5 Hz和f=0.037 5 Hz的工況進行數(shù)值模擬，并與f=0.025 Hz的工況結果進行對比。

f=0.012 5 Hz和f=0.037 5 Hz工況下滾轉角隨時間的變化如圖17和圖18所示。f=0.012 5 Hz工況下，在俯仰角的上半周期滾轉振蕩的振幅小于10°，在下半周期超過50°。f=0.037 5 Hz 工況下，在俯仰角的上半周期滾轉振蕩的振幅約為20°，在下半周期超過60°。結合前文Case2工況俯仰角振蕩頻率為f=0.025 Hz的滾轉角變化曲線可以看出，在俯仰運動的上半周期，滾轉角振蕩頻率隨俯仰角振蕩頻率的變化較小，而在下半周期的大幅振蕩頻率則與俯仰角振蕩頻率基本一致。

f=0.012 5，0.025，0.037 5 Hz 3種工況的滾轉相圖如圖19所示，縱軸為滾轉角速度?？梢钥吹皆?種工況下，相圖均為類似極限環(huán)的形態(tài)，一個半徑較小的圓形極限環(huán)與一個扁平的橢圓形極限環(huán)交替出現(xiàn)，兩個極限環(huán)分別對應上半周期的中小幅度振蕩以及下半周期的大振幅振蕩。在圖中可以看出圓環(huán)半徑隨俯仰運動頻率的增加而增大，表明上半周期振幅受俯仰運動頻率影響較大；而橢圓形右端位置變化較小，表明下半周期振幅受俯仰運動頻率影響較小。

圖17 f=0.012 5 Hz工況滾轉角和俯仰角時間歷程曲線

圖18 f=0.037 5 Hz工況滾轉角和俯仰角時間歷程曲線

圖19 滾轉角速度-滾轉角相圖

4.3 耦合運動機理分析

上文研究了飛行器轉動慣量以及俯仰頻率對飛行器耦合運動的影響，可以發(fā)現(xiàn)在各個工況下滾轉方向的運動都呈現(xiàn)為高頻中小幅振蕩與低頻大幅振蕩交替出現(xiàn)的趨勢。下面對耦合運動進行穩(wěn)定性分析[24-25]。

飛行器的氣動滾轉力矩Mx在平衡態(tài)附近可以展開為

(3)

由于飛行器的慣性耦合作用可以忽略，滾轉方向的運動學方程和動力學方程可表示為

(4)

式(4)可進一步寫為

(5)

選取Case1作為典型工況進行具體分析，該工況下滾轉力矩系數(shù)以及俯仰角、側滑角、滾轉角隨時間的變化曲線如圖20所示?？梢钥吹剑跐L轉角高頻振蕩的區(qū)間，例如40～60 s范圍內(nèi)，飛行器的俯仰角超過5°，此時飛行器的側滑角隨著飛行器的滾轉角而劇烈振蕩，使得飛行器在滾轉的過程中獲得了較大的滾轉回復力矩，維持飛行器的高頻滾轉振蕩。而在60～80 s范圍內(nèi)，飛行器的俯仰角逐漸減小至0°后回到5°，在此期間，雖然飛行器出現(xiàn)了大幅滾轉，但側滑角始終小于2°，滾轉力矩系數(shù)非常小且變化緩慢，靜穩(wěn)定性差，無法在短時間內(nèi)改變飛行器的滾轉方向。飛行器在俯仰角接近0°的位置出現(xiàn)滾轉角大幅震蕩，這種現(xiàn)象與Adamczak和Bolender[21]對同樣采用翼身組合體構型的HIFiRE-6一體化飛行器的動態(tài)特性研究時指出的飛行器受到的氣動阻尼較小，在小迎角下橫側向通道動不穩(wěn)定的結論有一定相似性。說明本文的研究結論對采用類似構型的飛行器有一定的參考價值。

飛行器在60～80 s這個滾轉角大幅振蕩的區(qū)間內(nèi)滾轉力矩系數(shù)隨滾轉角的變化曲線如圖21所示。曲線為雙“8”字型，曲線左右兩端呈逆時針旋轉，表明飛行器的滾轉運動受到氣流的阻尼作用，向氣流釋放能量，而曲線的中間部分呈順時針旋轉，表明氣流向飛行器的滾轉運動輸入能量。這種能量輸入輸出與俯仰方向的運動有關。圖20中紅色虛線位置為俯仰角為0°的時刻，可以看到由于滾轉角達到峰值的時刻滯后于俯仰角達到0°的時刻，在滾轉角增大的過程中對應的俯仰角小于滾轉角減小的過程中對應的俯仰角，這種差異導致了滾轉過程中能量的輸入和耗散。

圖20 滾轉力矩系數(shù)、俯仰角、側滑角和滾轉角時間歷程曲線

圖21 滾轉力矩系數(shù)隨滾轉角的變化

5 結 論

1) 雖然模型具有滾轉靜穩(wěn)定性和動穩(wěn)定性，但是在強迫俯仰/自由滾轉運動過程中，滾轉通道卻出現(xiàn)了中小幅度振蕩與大振幅振蕩交替出現(xiàn)的情況。中小幅度振蕩出現(xiàn)在俯仰運動的上半周期，大振幅振蕩出現(xiàn)在俯仰運動的下半周期，俯仰角接近0°的位置，最大滾轉角超過70°，這說明單自由度運動的動導數(shù)不能對兩自由度耦合運動進行準確預測。

2) 由于俯仰角速度較小，慣性耦合力矩很小，因而Iy、Iz變化對耦合運動帶來的影響可以忽略。而Ix對俯仰運動上半周期出現(xiàn)的滾轉角中小幅振蕩影響較為明顯，滾轉角振動頻率隨Ix的增大而減小，振幅隨著Ix的增大而明顯增大。

3) 俯仰運動頻率對上半周期出現(xiàn)的滾轉振蕩振幅影響較為明顯，對滾轉振蕩頻率影響較小。振幅隨俯仰運動頻率的增大而增大。

4) 滾轉角增大的過程中對應的俯仰角小于滾轉角減小的過程中對應的俯仰角，這種差異導致了大振幅滾轉振蕩過程中能量的輸入和耗散。