劉 緒，趙云飛，王東方，劉 偉

（國防科學(xué)技術(shù)大學(xué) 航天與材料工程學(xué)院，長沙410073）

高超聲速內(nèi)外流一體化飛行器的動態(tài)穩(wěn)定性參數(shù)（工程上常稱為“動導(dǎo)數(shù)”）研究是控制系統(tǒng)設(shè)計和動態(tài)品質(zhì)分析的重要參數(shù)，其工程需求主要體現(xiàn)在以下3個方面［1］：①飛行器軌道設(shè)計的重要參數(shù)；②飛行器姿態(tài)控制系統(tǒng)設(shè)計的重要參數(shù)；③飛行器縱、橫向動態(tài)穩(wěn)定性分析的重要依據(jù)。

隨著美國 HyperX［2］、FALCON（獵鷹）［3］和“黑雨燕”等項目以及在歐洲、日本等項目的相繼開展，吸氣式高超聲速飛行器的研究工作逐漸進(jìn)入到工程預(yù)發(fā)展階段［4］。內(nèi)外流一體化飛行器姿態(tài)控制精度要求很高以確保進(jìn)氣道啟動，其操穩(wěn)特性分析評估也更重要，需要準(zhǔn)確預(yù)測動導(dǎo)數(shù)。X－51“乘波者”高超聲速飛行器第2次飛行試驗，由于超燃沖壓發(fā)動機(jī)的進(jìn)氣道未能啟動而失敗。HTV－2首飛試驗失敗，在BTT（Bank to Turn）方式實現(xiàn)偏航過程中飛行器沿縱軸偏轉(zhuǎn)，在橫滾速度達(dá)到極限后超出可控范圍。美國工程審查委員會（ERB）分析認(rèn)為：對飛行過程中的若干空氣動力學(xué)關(guān)鍵參數(shù)認(rèn)識有限，在缺少主動控制的條件下，HTV－2將進(jìn)入彈道螺旋飛行（典型的橫側(cè)不穩(wěn)定問題）。

目前國內(nèi)外對動導(dǎo)數(shù)辨識方面的工作主要還是針對傳統(tǒng)的以火箭發(fā)動機(jī)為動力的彈箭類飛行器，對包含內(nèi)流的吸氣式飛行器動態(tài)特性研究較少。因此高超聲速內(nèi)外流一體化飛行器的研制亟需開展動態(tài)特性模擬技術(shù)研究。其核心是三方向（俯仰／偏航／滾轉(zhuǎn)）的直接阻尼導(dǎo)數(shù)模擬方法研究；同時還迫切需要開展加速度導(dǎo)數(shù)、交叉導(dǎo)數(shù)、交叉耦合導(dǎo)數(shù)的數(shù)值算法研究，為縱橫向耦合運(yùn)動特性研究提供技術(shù)支撐。動導(dǎo)數(shù)預(yù)測及穩(wěn)定性預(yù)示將為飛行器控制系統(tǒng)設(shè)計、高超聲速飛行器動不穩(wěn)定發(fā)生的邊界分析及相應(yīng)的動態(tài)穩(wěn)定性判據(jù)研究提供關(guān)鍵氣動參數(shù)。

動導(dǎo)數(shù)計算方法分為小擾動線化理論、牛頓理論與氣動力工程模型相結(jié)合的動導(dǎo)數(shù)工程近似方法、模擬動態(tài)實驗的數(shù)值自由振蕩法、數(shù)值強(qiáng)迫振蕩法。工程近似方法的最大特點(diǎn)在于快捷，但依賴于經(jīng)驗性，考慮氣流分離、再附和尾流等效應(yīng)較為困難，對于包含內(nèi)流的一體化復(fù)雜外形流動難以適應(yīng)。因此考慮了流動非線性的模擬動態(tài)實驗的全數(shù)值計算是動導(dǎo)數(shù)研究發(fā)展的重要方向。

本文在Etkin非定常氣動力模型基礎(chǔ)上，給出了小振幅強(qiáng)迫振動下的俯仰、偏航和滾轉(zhuǎn)三方向直接阻尼導(dǎo)數(shù)的強(qiáng)迫簡諧分析方法，并發(fā)展了包括交叉導(dǎo)數(shù)、交叉耦合導(dǎo)數(shù)在內(nèi)的多種動導(dǎo)數(shù)辨識方法。在采用Finner標(biāo)模驗證的基礎(chǔ)上，開展了高超聲速內(nèi)外流一體化飛行器的動態(tài)特性分析。

1 數(shù)值方法

在貼體坐標(biāo)系ξ－η－ζ下，對完全氣體、忽略質(zhì)量

力下的三維無量綱Navier－Stokes方程形式如下：

式中：U為守恒變量；E，F(xiàn)，G為無粘通量；Ev，F(xiàn)v，Gv為粘性通量。

本文在空間上采用二階精度的Roe格式。Roe格式是基于Godunov方法的基本思路發(fā)展起來的，是通量差分分裂（FDS）格式的一種，由于其優(yōu)秀的間斷分辨率和較小的數(shù)值耗散性，目前得到廣泛的應(yīng)用。采用引入“雙時間步”（dual－time－step）方法的LU－SGS隱式格式離散流體運(yùn)動方程時間項。動網(wǎng)格生成采用剛性網(wǎng)格生成方法。遠(yuǎn)場入流邊界采用適用于動態(tài)邊界條件下的一維Riemann不變量的無反射邊界條件。對于定常／非定常超音速出流，邊界點(diǎn)上的值由內(nèi)流場計算結(jié)果外推得到。對于壁面邊界，采用速度無滑移、絕熱壁條件。奇性軸采用外插后周向平均處理。

2 動導(dǎo)數(shù)計算方法

在飛行器姿態(tài)控制系統(tǒng)設(shè)計及軌道（彈道）設(shè)計中，所需要的動態(tài)穩(wěn)定性參數(shù)有數(shù)十個之多。表1列出了一些重要的動態(tài)阻尼導(dǎo)數(shù)［5］。表中：Cl，Cm，Cn分別為滾轉(zhuǎn)、俯仰、偏航力矩系數(shù)；α，β分別為攻角和側(cè)滑角；p，q，r分別為滾動軸、俯仰軸、偏航軸的角速度分量。目前國內(nèi)外公開文獻(xiàn)主要是研究俯仰、偏航或滾轉(zhuǎn)三方向的直接阻尼導(dǎo)數(shù)，而對交叉導(dǎo)數(shù)、交叉耦合導(dǎo)數(shù)的數(shù)值計算較少涉及。本文采用小振幅強(qiáng)迫簡諧分析法給出了多種導(dǎo)數(shù)的數(shù)值辨識方法。

表1 動態(tài)阻尼導(dǎo)數(shù)

飛行器做單自由度的俯仰運(yùn)動，如果其質(zhì)心速度不變，則確定運(yùn)動的獨(dú)立狀態(tài)變量只有攻角α和俯仰角速度q。對俯仰力矩系數(shù)Cm，根據(jù)Etkin假定［6］，俯仰力矩系數(shù)可寫成：

給定強(qiáng)迫振動：α（t）＝θ＝α0＋αmsinkt，其中，α為瞬時攻角，θ為俯仰角，α0為基準(zhǔn)狀態(tài)的攻角，αm為攻角振幅，k為減縮頻率。在k不很大且忽略高階導(dǎo)數(shù)的影響時，當(dāng)非定常振動過程達(dá)到諧振解，通過數(shù)值積分可以求出俯仰阻尼導(dǎo)數(shù)為

式中：ts為積分起始時間，Tc為振動周期。

類似地，通過給定不同的擾動方式可以推導(dǎo)出其他類型的動態(tài)穩(wěn)定性參數(shù)的計算公式。其中強(qiáng)迫滾轉(zhuǎn)運(yùn)動的形式為φ＝φmsinkt，φ為滾轉(zhuǎn)角，φm為滾轉(zhuǎn)角振幅；強(qiáng)迫偏航運(yùn)動的形式為ψ＝ψmsinkt，ψ為偏航角，ψm為偏航角振幅。

3 驗證算例

本文采用Finner標(biāo)模作為驗證算例，F(xiàn)inner標(biāo)模是外形為“十”字翼的導(dǎo)彈，其動態(tài)特性有標(biāo)準(zhǔn)實驗數(shù)據(jù)和基于準(zhǔn)定常歐拉方程的求解結(jié)果。圖1給出了“十”字翼導(dǎo)彈的外形尺寸［7］。導(dǎo)彈全長L為10倍的彈體直徑。本文采用塊間完全對接的多塊結(jié)構(gòu)化網(wǎng)格，對稱面網(wǎng)格與分區(qū)情況見圖2。網(wǎng)格總量為200萬。近壁第一層網(wǎng)格到壁面的距離為1×10－4L。

圖1 “十”字翼導(dǎo)彈外形尺寸

圖2 對稱面網(wǎng)格與分區(qū)圖

3.1 俯仰阻尼導(dǎo)數(shù)辨識

給定強(qiáng)迫俯仰振動形式α（t）＝θ＝α0＋αmsinkt，其中初始攻角α0＝1.5°，振幅αm＝1.5°，減縮頻率k＝0.05，質(zhì)心位置Xcg／L＝0.5。

圖3給出了俯仰阻尼導(dǎo)數(shù)的計算結(jié)果與實驗數(shù)據(jù)［8］和基于準(zhǔn)定常歐拉方程結(jié)果［9］的比較。從中可以看出，本文預(yù)測的俯仰阻尼導(dǎo)數(shù)比文獻(xiàn)中采用準(zhǔn)定常歐拉方程的求解結(jié)果更加接近風(fēng)洞實驗的數(shù)據(jù)。

圖3 俯仰阻尼導(dǎo)數(shù)計算結(jié)果

3.2 滾轉(zhuǎn)阻尼導(dǎo)數(shù)辨識

給定強(qiáng)迫滾轉(zhuǎn)振動形式φ＝φmsinkt，取φm＝2.5°，k＝0.05，Xcg／L＝0.5。

圖4將強(qiáng)迫滾轉(zhuǎn)辨識出的動導(dǎo)數(shù)與文獻(xiàn)［10－11］中的實驗數(shù)據(jù)進(jìn)行了比較。滾轉(zhuǎn)阻尼導(dǎo)數(shù)的實驗精度一般不高，主要由于實驗中存在支架、洞壁干擾等不確定因素，這使對同一模型在不同風(fēng)洞中的實驗結(jié)果也存在一定的差異。本文計算得到的滾轉(zhuǎn)阻尼導(dǎo)數(shù)介于2個實驗結(jié)果之間，驗證了方法和程序的可靠性。

圖4 滾轉(zhuǎn)阻尼導(dǎo)數(shù)計算結(jié)果

4 動態(tài)特性分析

高超聲速巡航導(dǎo)彈X－51是美國空軍研究實驗室（AFRL）聯(lián)合波音公司（負(fù)責(zé)機(jī)身）和普惠公司（負(fù)責(zé)發(fā)動機(jī)）研制的一架超燃沖壓發(fā)動機(jī)驗證機(jī)－乘波器（SED－WR）飛行試驗平臺。本文以內(nèi)外流一體化飛行器X－51為背景，利用公開的圖像數(shù)據(jù)資料進(jìn)行反向建模，生成的類X－51三維實體外形見圖5，其進(jìn)氣道內(nèi)部結(jié)構(gòu)見圖6。通過開展小振幅強(qiáng)迫簡諧運(yùn)動的非定常流場數(shù)值模擬，對進(jìn)氣道冷流狀態(tài)下的類X－51動態(tài)穩(wěn)定性參數(shù)進(jìn)行數(shù)值辨識。

圖5 類X－51三維實體外形

圖6 進(jìn)氣道內(nèi)部結(jié)構(gòu)

4.1 直接阻尼導(dǎo)數(shù)辨識

計算采用塊間完全對接的多塊結(jié)構(gòu)化網(wǎng)格，網(wǎng)格總量為460萬，共劃分為62塊。不同位置的結(jié)構(gòu)化網(wǎng)格與拓?fù)浣Y(jié)構(gòu)見圖7和圖8。計算條件為：馬赫數(shù)6.5，高度27km。俯仰／偏航／滾轉(zhuǎn)3個方向的強(qiáng)迫簡諧運(yùn)動的振幅均為1°，k＝0.1，無量綱時間步長Δt＝0.005。圖9給出了攻角為0°時的強(qiáng)迫俯仰振動一個周期內(nèi)不同時刻的壓力等值線，其中ti＝ts＋（i－1）Tc／4，i＝1，2，3，4。頭部斜激波、膨脹波以及尾翼激波在圖中清楚地反映出來，內(nèi)部流動已經(jīng)建立。同時可以看出強(qiáng)迫振動的振幅很小導(dǎo)致各時刻流場狀態(tài)差別不大。在攻角0°～6°范圍之間本文給出了3個方向強(qiáng)迫振動的力矩系數(shù)曲線，如圖10～圖12所示，圖中，Cmm，Cml，Cmn分別為俯仰、滾轉(zhuǎn)、偏航力矩系數(shù)。

圖7 對稱面網(wǎng)格與分區(qū)圖

圖8 舵面附近的結(jié)構(gòu)化網(wǎng)格

圖9 一個周期內(nèi)不同時刻的壓力等值線

圖10 （a）、圖11（a）和圖12（a）分別給出了俯仰、滾轉(zhuǎn)和偏航運(yùn)動對應(yīng)力矩系數(shù)的遲滯環(huán)曲線。這是非定常運(yùn)動時物體上漩渦運(yùn)動和物面運(yùn)動之間存在的時間延遲現(xiàn)象在氣動力系數(shù)上的反映。遲滯環(huán)形狀飽滿，非定常滯后效應(yīng)比較明顯。在計算啟動后一個振蕩周期內(nèi)即進(jìn)入遲滯環(huán)，遲滯環(huán)重復(fù)性較好。不同攻角下的強(qiáng)迫振蕩遲滯環(huán)的轉(zhuǎn)動方式相互一致。從與3個遲滯環(huán)相對應(yīng)的時間歷程曲線圖10（b）、圖11（c）和圖12（d）可以看到，非定常氣動力矩收斂情況很好，從第2個周期開始已經(jīng)完全達(dá)到諧振。3個振動方向上的俯仰力矩系數(shù)、滾轉(zhuǎn)力矩系數(shù)、偏航力矩系數(shù)隨時間按正弦規(guī)律變化。

圖10 強(qiáng)迫俯仰振動力矩系數(shù)曲線

圖11 強(qiáng)迫滾轉(zhuǎn)振動力矩系數(shù)曲線

通過積分強(qiáng)迫簡諧振動的遲滯環(huán)曲線，表2給出了俯仰／偏航／滾轉(zhuǎn)3個方向的直接阻尼導(dǎo)數(shù)辨識結(jié)果。各方向下的直接阻尼導(dǎo)數(shù)均為負(fù)值，說明飛行器在俯仰／偏航／滾轉(zhuǎn)方向受到擾動后處于動穩(wěn)定狀態(tài)。俯仰阻尼導(dǎo)數(shù)和偏航阻尼導(dǎo)數(shù)量級一致，但滾轉(zhuǎn)阻尼導(dǎo)數(shù)要小1～2個量級。

圖12 強(qiáng)迫偏航振動力矩系數(shù)曲線

表2 直接阻尼導(dǎo)數(shù)辨識結(jié)果

4.2 交叉阻尼導(dǎo)數(shù)辨識

圖11（d）和圖12（c）的偏航和滾轉(zhuǎn)力矩系數(shù)隨時間變化能夠達(dá)到諧振，周期性變化明顯，說明滾轉(zhuǎn)方向和偏航方向的交叉性較強(qiáng)。但圖10（c）和圖11（b）中的力矩系數(shù)基本為一個常數(shù)，受周期性簡諧運(yùn)動的影響不大，說明俯仰方向和滾轉(zhuǎn)方向的交叉耦合性很弱。類似地，圖10（d）和圖12（b）顯示出俯仰方向和偏航方向的交叉耦合性同樣很弱。綜合上述規(guī)律，在本狀態(tài)下飛行器偏航與滾轉(zhuǎn)橫向之間的影響明顯，但縱橫向交叉耦合性不強(qiáng)，故本文對交叉導(dǎo)數(shù)開展了數(shù)值辨識工作，其結(jié)果列于表3。

表3 交叉阻尼導(dǎo)數(shù)辨識結(jié)果

從辨識結(jié)果來看，強(qiáng)迫滾轉(zhuǎn)運(yùn)動時滾轉(zhuǎn)阻尼導(dǎo)數(shù)在各計算狀態(tài)下均為負(fù)值，而滾轉(zhuǎn)－偏航力矩交叉導(dǎo)數(shù)均為正值，二者量級相差不大，說明飛行器在滾轉(zhuǎn)時是動穩(wěn)定的，但可能引起偏航方向的動不穩(wěn)定問題，需要控制系統(tǒng)對其姿態(tài)進(jìn)行控制。強(qiáng)迫偏航運(yùn)動與此類似，偏航阻尼導(dǎo)數(shù)在各計算狀態(tài)下均為負(fù)，而偏航－滾轉(zhuǎn)力矩交叉導(dǎo)數(shù)均為正，但其量級比偏航阻尼導(dǎo)數(shù)小2個量級。計算時涉及表面摩阻等復(fù)雜的計算問題，增加了計算難度，因此交叉導(dǎo)數(shù)的計算還需做深入研究。

4.3 強(qiáng)迫振蕩計算參數(shù)選取比較

強(qiáng)迫振動法辨識動態(tài)穩(wěn)定性參數(shù)的精度依賴于3個重要的計算參數(shù)：時間步長、子迭代步數(shù)、減縮頻率。時間步長和子迭代步數(shù)這2個參數(shù)的組合實際上決定了雙時間步方法計算非定常流動的收斂程度。頻率相似問題是影響強(qiáng)迫振動法辨識動導(dǎo)數(shù)精度的重要影響因素。本文選取了幾組不同的參數(shù)值，分別考察時間步長、子迭代步數(shù)、減縮頻率對包含內(nèi)流的高超聲速飛行器動導(dǎo)數(shù)辨識的影響。

減小時間步長、加大子迭代步數(shù)都可以提高非定常流動的收斂程度，從而提高動態(tài)穩(wěn)定性參數(shù)的辨識精度，但計算量也會大幅增加。因此有必要先進(jìn)行強(qiáng)迫振蕩的試算，確定保證準(zhǔn)確度的適宜計算時間步長和子迭代步數(shù)。本文無量綱時間步長Δt取0.005，子迭代步數(shù)n為5步，在保證計算精度的前提下有較好的計算效率。表4給出了更小的時間步長和更大的子迭代步數(shù)時的動導(dǎo)數(shù)辨識結(jié)果，可以看到計算結(jié)果基本一致。

表4 不同時間步長和子迭代步數(shù)下的動導(dǎo)數(shù)辨識結(jié)果

減縮頻率實際上是強(qiáng)迫振動頻率快慢的反映，其大小會影響動導(dǎo)數(shù)的量值乃至符號。減縮頻率的選取必須像馬赫數(shù)、攻角一樣考慮到與實驗的相似性。本文以強(qiáng)迫滾轉(zhuǎn)運(yùn)動為例，考察了不同減縮頻率對動導(dǎo)數(shù)辨識的影響。圖13顯示出遲滯環(huán)隨減縮頻率增加由內(nèi)而外形成嵌套。減縮頻率越小，遲滯環(huán)曲線越狹長，非定常滯后效應(yīng)越弱，動導(dǎo)數(shù)辨識難度高，非定常流場的計算量增大。但考慮到與實際情況的相似性，k不應(yīng)取得過大，否則影響計算精度。表5給出了不同k值的滾轉(zhuǎn)阻尼導(dǎo)數(shù)和滾轉(zhuǎn)－偏航交叉導(dǎo)數(shù)，可以看出減縮頻率對動導(dǎo)數(shù)辨識結(jié)果的影響，隨著k的增大，橫向之間的影響加強(qiáng)。

圖13 強(qiáng)迫滾轉(zhuǎn)振動不同減縮頻率下的遲滯環(huán)曲線

表5 不同減縮頻率下的動導(dǎo)數(shù)辨識結(jié)果

5 結(jié)論

本文基于三維非定常N－S方程，采用小振幅強(qiáng)迫簡諧分析法開展了高超聲速內(nèi)外流一體化飛行器動態(tài)氣動參數(shù)的計算研究。本文給出的動導(dǎo)數(shù)辨識技術(shù)不僅適用于傳統(tǒng)的以火箭發(fā)動機(jī)為動力的彈箭類飛行器，對包含內(nèi)流的吸氣式高超聲速飛行器適用性良好，各攻角下的非定常滯后效應(yīng)明顯，遲滯環(huán)重復(fù)性較強(qiáng)。辨識得到的高超聲速飛行器3個方向的直接阻尼導(dǎo)數(shù)均為負(fù)值。對橫側(cè)向動穩(wěn)定性、縱橫向耦合動穩(wěn)定性問題的分析顯示，其縱橫向耦合性不強(qiáng)，但橫向之間的交叉影響明顯，各交叉導(dǎo)數(shù)均為正值，但可能引起偏航／滾轉(zhuǎn)方向的動不穩(wěn)定問題，需要控制系統(tǒng)對其姿態(tài)進(jìn)行控制。

［1］劉偉.細(xì)長機(jī)翼搖滾機(jī)理的非線性動力學(xué)分析及數(shù)值模擬方法研究［D］.長沙：國防科學(xué)技術(shù)大學(xué)，2004.LIU Wei.Nonlinear dynamics analysis for mechanism of slender wing rock and study of numerical simulation method［D］.Changsha：National University of Defense Technology，2004.（in Chinese）

［2］PEEBLES C.The X－43Aflight research program：lessons learned on the road to Mach 10［M］.Washington，D.C.：AIAA Inc.，2007.

［3］WALKER S，RODGERS F.Falcon hypersonic technology overview，AIAA－2005－3253［R］.2005.

［4］楊超，許赟，謝長川.高超聲速飛行器氣動彈性力學(xué)研究綜述［J］.航空學(xué)報，2010，31（1）：1－11.YANG Chao，XU Yun，XIE Chang－chuan.Review of studies on aeroelasticity of hypersonic vehicles［J］.Acta Aeronautica et Astronautica Sinica，2010，31（1）：1－11.（in Chinese）

［5］童秉綱，陳強(qiáng).關(guān)于非定?？諝鈩恿W(xué)［J］.力學(xué)進(jìn)展，1983，13（4）：377－394.TONG Bing－gang，CHEN Qiang.Some remarks on unsteady aerodynamics［J］.Advances in Mechanics，1983，13（4）：377－394.（in Chinese）

［6］劉偉，楊小亮，趙云飛.高超聲速飛行器加速度導(dǎo)數(shù)數(shù)值模擬［J］.空氣動力學(xué)學(xué)報，2010，28（4）：426－429.LIU Wei，YANG Xiao－liang，ZHAO Yun－fei.Numerical simulation of acceleration derivative of hypersonic aircraft［J］.Acta Aerodynamica Sinica，2010，28（4）：426－429.（in Chinese）

［7］MIKHAIL A G.Roll damping for projectiles including wraparound，offset，and arbitrary number of fins［J］.Journal of Spacecraft and Rockets，1995，32（6）：929－937.

［8］SHANTZ I，GROVES R T.Dynamic and static stability measurements of the basic finner at supersonic speeds，NAVORD Report 4516［R］.1960.

［9］OKTAY E.CFD predictions of dynamic derivatives for missiles，AIAA 2002－0276［R］.2002.

［10］REGAN F J.Roll damping moment measurements for the basic finner at subsonic and supersonic speeds，NAVORD Report 6652［R］.1964.

［11］WHYTE R H.Spinner－a computer program for predicting the aerodynamic coefficients of spin stabilized projectiles，Genera1 Electric Class 2Reports［R］.1969.