叢戎飛，葉友達，趙忠良

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

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

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

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

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

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

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

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

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

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

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

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

本文研究的吸氣式高超聲速飛行器模型是參照美國SR-72[27-28]概念飛行器設(shè)計的，采用了翼身組合體氣動布局，與HTV-2和X-51A等高超聲速飛行器普遍采用的“乘波體”構(gòu)型有較大的不同。圖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導(dǎo)彈標模，利用數(shù)值強迫振蕩的方法，計算導(dǎo)彈的俯仰、滾轉(zhuǎn)動導(dǎo)數(shù)，并與試驗結(jié)果對比，驗證程序?qū)τ跉鈩?運動耦合狀態(tài)下非定常氣動力的計算能力。

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

圖5 Basic Finner滾轉(zhuǎn)動導(dǎo)數(shù)

3 滾轉(zhuǎn)穩(wěn)定性

3.1 靜穩(wěn)定性

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

圖7 滾轉(zhuǎn)靜導(dǎo)數(shù)

3.2 動穩(wěn)定性

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

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

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

圖9 滾轉(zhuǎn)動導(dǎo)數(shù)

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

圖10 自由滾轉(zhuǎn)運動時滾轉(zhuǎn)角時間歷程曲線

圖11 自由滾轉(zhuǎn)運動相圖

圖12 滾轉(zhuǎn)運動頻譜圖

4 強迫俯仰/自由滾轉(zhuǎn)兩自由度耦合運動

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

4.1 轉(zhuǎn)動慣量影響

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

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

表1 模型轉(zhuǎn)動慣量

(1)

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

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

圖14 慣性耦合力矩與氣動滾轉(zhuǎn)力矩時間歷程曲線

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

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

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

(2)

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

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

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

4.2 俯仰運動頻率影響

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

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

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

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

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

圖19 滾轉(zhuǎn)角速度-滾轉(zhuǎn)角相圖

4.3 耦合運動機理分析

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

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

(3)

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

(4)

式(4)可進一步寫為

(5)

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

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

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

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

5 結(jié) 論

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

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

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

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