乘波體荷蘭滾模態(tài)特性研究

劉 文, 張陳安, 王發(fā)民

(1.西北工業(yè)大學(xué) 翼型葉柵國防科技重點實驗室, 陜西 西安 710072;2.中國科學(xué)院力學(xué)研究所 高溫氣體動力學(xué)國家重點實驗室, 北京 100190)

?

乘波體荷蘭滾模態(tài)特性研究

劉 文1,2, 張陳安2,*, 王發(fā)民2

(1.西北工業(yè)大學(xué) 翼型葉柵國防科技重點實驗室, 陜西 西安 710072;2.中國科學(xué)院力學(xué)研究所 高溫氣體動力學(xué)國家重點實驗室, 北京 100190)

乘波體非軸對稱、扁平、大長細(xì)比的幾何外形特點決定了其存在嚴(yán)重的橫航向耦合動穩(wěn)定性問題。目前對乘波體橫航向穩(wěn)定性的研究還相對較少，且一般只針對單一外形，對飛行器設(shè)計未能獲得有指導(dǎo)意義的定性、定量結(jié)論。以冪次乘波體為研究對象，首先引入設(shè)計參數(shù)kw和φ描述其外形特點，然后結(jié)合CFD數(shù)值模擬和Kriging代理模型，獲得了整個設(shè)計參數(shù)空間內(nèi)乘波體的靜/動導(dǎo)數(shù)，進(jìn)而通過求解飛行動力學(xué)耦合方程特征根獲得乘波體的橫航向耦合動穩(wěn)定荷蘭滾模態(tài)特性。定義了荷蘭滾動穩(wěn)定性導(dǎo)數(shù)的概念，推導(dǎo)了荷蘭滾阻尼近似表達(dá)式，解釋了不同攻角下荷蘭滾模態(tài)發(fā)散/收斂的成因，獲得了荷蘭滾阻尼隨設(shè)計參數(shù)和攻角的分布規(guī)律。推導(dǎo)了荷蘭滾頻率近似表達(dá)式，獲得了荷蘭滾頻率隨設(shè)計參數(shù)和攻角的分布規(guī)律。

乘波體；高超聲速流動；動穩(wěn)定；荷蘭滾

0 引 言

高超聲速飛行方式正由傳統(tǒng)的彈道式飛行和太空返回向近空間大氣層內(nèi)做長時間高速機動飛行的方向發(fā)展。以美國的“HTV-2”為代表的高超聲速滑翔飛行器是其典型代表之一，該類飛行器以單級或多級火箭為動力，助推到一定高度和馬赫數(shù)后箭體分離，之后滑翔器在近空間進(jìn)行高超聲速滑翔機動飛行。美國國防預(yù)先研究局對“HTV-2”進(jìn)行過兩次試飛，均發(fā)生了飛行器因橫航向動態(tài)失穩(wěn)而墜毀[1-2]。

乘波體的概念由Nonweiler[3]在1959年提出，近30年來得到了廣泛發(fā)展。其原理是在乘波體的前緣產(chǎn)生附體激波，阻止了下表面的高壓氣體上溢，使波后高壓氣體都限制在下表面，從而獲得在相同攻角下比普通外形更高的升力以及升阻比。乘波體的高升阻比優(yōu)勢，使其成為高超聲速滑翔飛行器的一種非常理想的潛在構(gòu)型。

與HTV-2類似，乘波體非軸對稱、扁平、大長細(xì)比的幾何外形特點也決定了乘波體飛行器將面臨類似的橫航向穩(wěn)定性問題。首先，側(cè)滑對橫向滾轉(zhuǎn)力矩的影響不可忽略，同時大長細(xì)比外形使其橫向轉(zhuǎn)動慣量相比縱向和航向要小的多，因而在出現(xiàn)側(cè)滑時會引起明顯的橫滾效應(yīng)；第二，具有尖前緣特征的乘波體飛行器在高超聲速飛行時同樣可能出現(xiàn)由于氣動燒蝕引起防熱層脫落，進(jìn)而誘導(dǎo)出激波引起附加的橫航向氣動力/力矩，乘波體下表面可能出現(xiàn)的非對稱轉(zhuǎn)捩會引發(fā)下游流動的不對稱，形成橫向的環(huán)流和橫航向的非對稱氣動力[4]，這些影響因素難以通過計算或?qū)嶒炦M(jìn)行評估，因而對飛行器本身的橫航向穩(wěn)定性提出了更苛刻的要求；第三，扁平的幾何外形使得其由側(cè)滑產(chǎn)生的側(cè)向力較小，難以通過直接側(cè)滑拐彎方式直接獲得較大的側(cè)向加速度，所以通常采用傾斜拐彎方式，也會面臨橫航向運動耦合問題[5]；第四，飛行器在高超聲速飛行時氣動阻尼很小，橫航向穩(wěn)定性較低，容易發(fā)生失穩(wěn)現(xiàn)象，由于飛行速度極快，即使平靜、緩慢發(fā)展的耦合不穩(wěn)定也可能演化成迅速發(fā)展的耦合不穩(wěn)定，使飛行器產(chǎn)生強烈的不可控運動[6]。綜上，橫航向穩(wěn)定性問題是高超聲速乘波體飛行器設(shè)計過程中面臨的一個嚴(yán)峻考驗。

研究高超聲速飛行器橫航向耦合運動機理的最常用方法是對橫航向小擾動線化模型進(jìn)行簡化，獲得各個模態(tài)的近似表達(dá)式，進(jìn)而分析影響各個模態(tài)的主要氣動參數(shù)。已有文獻(xiàn)[7-12]對飛行器在亞、跨、超聲速的橫航向模態(tài)進(jìn)行了簡化分析，但是高超聲速飛行器所處的高空高馬赫數(shù)飛行條件、大后掠氣動布局特點以及慣量沿縱軸集中的質(zhì)量分布特點，與傳統(tǒng)飛行器有很大不同，因而無法通過傳統(tǒng)結(jié)論來直接分析高超聲速飛行器的橫航向運動模態(tài)。例如，傳統(tǒng)上認(rèn)為荷蘭滾模態(tài)僅包含偏航和側(cè)滑運動，忽略滾轉(zhuǎn)運動的物理量，即忽略與滾轉(zhuǎn)角速度p和滾轉(zhuǎn)角φ相關(guān)項，只保留與偏航角速度r和側(cè)滑角β相關(guān)項。Phillips[13]通過推導(dǎo)指出，這種傳統(tǒng)上忽略滾轉(zhuǎn)運動的荷蘭滾近似是滾轉(zhuǎn)阻尼無窮大時的漸近解，實際上大多數(shù)情況滾轉(zhuǎn)阻尼沒有大到可以忽略荷蘭滾模態(tài)中的滾轉(zhuǎn)運動。在高超聲速飛行條件下，滾轉(zhuǎn)阻尼較小，滾轉(zhuǎn)運動在荷蘭滾模態(tài)中甚至起主導(dǎo)作用，傳統(tǒng)的荷蘭滾模態(tài)分析將得到明顯錯誤的結(jié)論。針對上述問題，近年來已有少量文獻(xiàn)對高超聲速飛行時的橫航向耦合模態(tài)機理進(jìn)行了研究。

Heller和Breitsamter等[14-15]推導(dǎo)了高超聲速飛行器HYTEX R-A3的模態(tài)近似表達(dá)式，分析了橫航向模態(tài)隨馬赫數(shù)和質(zhì)心的變化規(guī)律，并推導(dǎo)了滾轉(zhuǎn)和螺旋模態(tài)耦合為橫向長周期模態(tài)的條件，最后指出滾轉(zhuǎn)阻尼小、荷蘭滾不穩(wěn)定和滾轉(zhuǎn)-偏航耦合度高是導(dǎo)致高超聲速飛行器面臨飛行品質(zhì)缺陷的主要原因。

Yin等[16]推導(dǎo)了荷蘭滾模態(tài)中滾轉(zhuǎn)角速度與偏航角速度振幅比值|p/r|的表達(dá)式，結(jié)果表明由于高超聲速飛行器較大的飛行攻角和細(xì)長體的外形特征，該比值一般較大，使得傳統(tǒng)荷蘭滾阻尼比近似失效。同時由于|p/r|較大，因而他們忽略了r的相關(guān)項得到了荷蘭滾阻尼比的近似表達(dá)式，結(jié)果表明其大小主要受滾轉(zhuǎn)阻尼導(dǎo)數(shù)和橫向靜導(dǎo)數(shù)的影響。

從目前對高超聲速飛行器橫航向模態(tài)的少量研究來看，還存在幾個方面的不足：1) 對荷蘭滾模態(tài)頻率的分析理論較為完善，但缺乏對荷蘭滾阻尼有效而簡潔的近似表達(dá)式，因而難以分析荷蘭滾阻尼特性的主要影響因素； 2) 已有研究大多數(shù)針對某一個而不是某一類飛行器來研究，結(jié)論缺乏普適性；3) 缺乏對高超聲速乘波體類飛行器橫航向耦合運動機理性研究。

針對上述問題，本文嘗試對高超聲速乘波體飛行器橫航向耦合動穩(wěn)定模態(tài)中最為復(fù)雜的荷蘭滾模態(tài)特性進(jìn)行研究。建立幾何特征可參數(shù)化描述的乘波體模型，通過CFD定常/非定常氣動力計算和代理模型獲得整個樣本設(shè)計空間內(nèi)所有乘波體的靜/動導(dǎo)數(shù)，結(jié)合橫航向線化小擾動模型及模態(tài)簡化，得到荷蘭滾模態(tài)頻率和阻尼特性隨外形參數(shù)的變化規(guī)律及主要影響因素。

1 研究對象和計算方法

1.1 冪次乘波體

乘波體是從已知流場中通過反設(shè)計方法獲得的一種高升阻比飛行器。由于乘波體前緣線附著在激波面上，將波后高壓氣體限制在下表面，從而在小攻角下即可獲得較高升力。通過冪次軸對稱流場獲得的冪次乘波體是一種容易設(shè)計成縱向靜穩(wěn)定的構(gòu)型[17]，因而本文以冪次乘波體為研究對象,其外形由冪次曲線y=Cxn中的參數(shù)C和n決定，此處選用的外形參數(shù)為：C=1，n=0.7；設(shè)計馬赫數(shù)選為15，繞冪次軸對稱體的基準(zhǔn)流場通過CFD求解Euler方程獲得。 乘波體通常可通過前緣在底面的投影進(jìn)行定義。在底面上定義一條前緣線的基準(zhǔn)曲線，將該曲線沿自由來流方向在激波面上投影所獲得的交線即為乘波體前緣線，對前緣線上的點進(jìn)行流線追蹤即可獲得乘波體，生成過程如圖1所示。

圖1 乘波體生成示意圖Fig.1 Generation of waverider

采用三次多項式定義基準(zhǔn)曲線：

(1)

令基準(zhǔn)曲線與激波圓的交點處斜率為0，設(shè)激波圓半徑為Rs，基準(zhǔn)曲線縱截距絕對值為R0，方位角為φ，令參數(shù)kw=R0/Rs，化簡可得：

(2)

選取距離冪次體頂點200 m處的平面為基準(zhǔn)面，則根據(jù)冪次曲線方程可確定Rs，然后給定設(shè)計參數(shù)kw和φ的值即可確定基準(zhǔn)曲線方程，從而生成冪次乘波體。

1.2 數(shù)值模擬方法

采用格心格式有限體積法求解RANS(Reynolds-Averaged Navier-Stokes)方程，空間離散采用AUSM+格式，時間離散采用隱式LU-SGS方法。在計算非定常氣動力時，采用雙時間法求解控制方程，對子迭代采用四階龍格-庫塔法進(jìn)行推進(jìn)，使用基于RBF內(nèi)插的網(wǎng)格變形技術(shù)[18]。選用層流和完全氣體模型，忽略真實氣體效應(yīng)的影響。

文獻(xiàn)[19]采用上述計算方法獲得了高空高馬赫數(shù)條件下航天飛機“OV-102”的升阻比，并與試飛結(jié)果進(jìn)行對比，表明了計算程序在強粘性干擾飛行條件下的可靠性；文獻(xiàn)[20]對比了跨音速條件下NACA0012翼型的非定常氣動力計算與實驗結(jié)果，驗證了求解器對非定常氣動力計算的可靠性，此處不再贅述。

1.3 橫航向耦合飛行力學(xué)模型

本文通過求解線化橫航向小擾動方程特征矩陣的特征根，來獲得乘波體的開環(huán)橫航向耦合動穩(wěn)定模態(tài)特性[21]。在體軸系下的方程為：

(3)其中Yβ=CyβqS/(mV)，Lβ=ClβqSb/Ix，Nβ=CnβqSb/Iz，Li=CliqSb2/(2IxV)，Ni=CniqSb2/(2IzV)，i∈{p,r}。式中α、θ、β、p、r、φ分別表示攻角、俯仰角、側(cè)滑角、滾轉(zhuǎn)角速度、偏航角速度和滾轉(zhuǎn)角；Y、L、N分別表示側(cè)向力、滾轉(zhuǎn)力矩和偏航力矩，Yβ、Lβ、Nβ為相應(yīng)氣動力對側(cè)滑角的靜導(dǎo)數(shù)；Lp、Np、Lr、Nr為相應(yīng)氣動力對相應(yīng)角速度的動導(dǎo)數(shù)，此處不考慮航跡角，因而有α=θ。

靜導(dǎo)數(shù)通過求解靜態(tài)氣動力即可獲得，動導(dǎo)數(shù)采用簡諧振蕩法進(jìn)行計算[22]。所有外形采用相同的質(zhì)量(m=500 kg)和慣量(Ix=75 kg·m2，Iz=750 kg·m2，忽略Ixz)，乘波體長度L=4 m，質(zhì)心位置為(0.59L，-0.05L，0)，其中坐標(biāo)原點位于乘波體頭部頂點處。研究的飛行工況為：Ma=15，H=50 km。

1.4 代理模型

由于需要研究乘波體橫航向動穩(wěn)定模態(tài)特性在一定設(shè)計參數(shù)空間內(nèi)的分布，而對于每個點都通過CFD求得其相應(yīng)的靜/動導(dǎo)數(shù)是不現(xiàn)實的。為了解決該問題，采用Kriging代理模型，在求得采樣點的靜/動導(dǎo)數(shù)后，其他點通過代理模型插值獲得相應(yīng)的靜/動導(dǎo)數(shù)。

表征乘波體外形的設(shè)計參數(shù)kw和φ的取值范圍分別選為[0.3,0.5]，[35°,55°]。為了提高插值精度，對于每個研究工況，采用全析因試驗設(shè)計方法對kw和φ各取9個水平。通過CFD獲得81個設(shè)計樣本點的橫航向靜/動導(dǎo)數(shù)，設(shè)計空間內(nèi)其它點通過代理模型插值獲得。

為了驗證所構(gòu)建代理模型計算靜/動導(dǎo)數(shù)的精度，在設(shè)計空間隨機選取五個樣本點，通過CFD求解氣動力獲得靜/動導(dǎo)數(shù)作為基準(zhǔn)值，與通過代理模型獲得的靜/動導(dǎo)數(shù)進(jìn)行對比。以零攻角的Lβ和Lp為例，相對誤差如表2所列，最大誤差僅為2.54%?？梢?，由于訓(xùn)練樣本點取得較密，構(gòu)建的代理模型有著較高的精度。

表1 測試樣本點Lβ和Lp相對誤差Table 1 Relative error for Lβ and Lp of test samples

2 橫航向靜穩(wěn)定性

在獲得橫航向耦合模態(tài)特性分布之前，根據(jù)CFD計算和代理模型的結(jié)果，我們首先可以得到橫航向靜穩(wěn)定性隨參數(shù)的分布規(guī)律。

對無尾飛行器來說，影響橫向靜穩(wěn)定性的主要幾何特征為上反角和后掠角。由于乘波體外形的特殊性，無法按照常規(guī)飛機那樣精確定義其上反角Γ和后掠角Λ，為此我們將其按圖2所示定義為：

(4)

乘波體三維外形隨kw和φ的變化特征如圖3和圖4所示，其上反角和后掠角分布如圖5和圖6所示。結(jié)合圖3～圖6可知，上反角隨著kw和φ的增大而增大，后掠角隨著kw和φ的增大而減小。

圖2 上反角和后掠角示意圖Fig.2 Sketch of dihedral and sweep angle

圖3 外形隨φ變化示意圖，kw=0.4Fig.3 Variation of shape with φ，kw=0.4

圖4 外形隨kw變化示意圖，φ=45°Fig.4 Variation of shape with kw，φ=45°

圖5 上反角分布Fig.5 Distribution of dihedral angle

圖6 后掠角分布Fig.6 Distribution of sweep angle

圖7和圖8分別給出了三個攻角的橫航向靜導(dǎo)數(shù)分布。可以看到在靠近參數(shù)空間區(qū)域的右上角，橫向和航向靜穩(wěn)定性都較強，靠近左下角則較弱。

對于橫向靜穩(wěn)定性來說，隨著攻角的增加，穩(wěn)定性都逐漸增強，且分布規(guī)律與上反角一致，其可能原因是乘波體本身后掠角已經(jīng)很大，其隨參數(shù)變化的相對變化幅度較小，因而影響橫向靜穩(wěn)定性的主要幾何特征為上反角；0°攻角時靠近左下角虛線以下區(qū)域Lβ>0，即橫向靜不穩(wěn)定，其主要原因是上反角較?。辉黾庸ソ强梢云鸬缴戏醋饔?，因而隨攻角增加橫向靜穩(wěn)定性顯著增強。

對于航向靜穩(wěn)定性來說，三個攻角所有區(qū)域都是靜穩(wěn)定的。在定常側(cè)滑直線飛行時，偏航力矩主要取決于機身側(cè)面積的大小以及航向壓心與質(zhì)心的距離。乘波體的下表面由流線追蹤而成，上表面與自由來流平行，因此乘波體靠近尾部的側(cè)面積比靠近頭部的側(cè)面積明顯要大，其接近尾部的側(cè)向力就大，只要質(zhì)心位置不是過于靠后，是容易實現(xiàn)航向靜穩(wěn)定性的。

3 橫航向耦合動穩(wěn)定特性

橫航向耦合動穩(wěn)定性一般可以由三個模態(tài)表征：螺旋模態(tài)、滾轉(zhuǎn)模態(tài)和荷蘭滾模態(tài)。其中荷蘭滾模態(tài)是運動形式最復(fù)雜、控制難度最大的模態(tài)。目前對高超聲速飛行器的橫航向模態(tài)研究很少，且已有的研究一般在0°攻角下進(jìn)行，不能考慮攻角的影響；同時對荷蘭滾阻尼的近似較為復(fù)雜，不易分析影響阻尼特性的主要因素。因而，本文結(jié)合高超聲速乘波體類飛行器靜/動導(dǎo)數(shù)的特點，通過模態(tài)簡化來獲得荷蘭滾模態(tài)特性隨外形參數(shù)的分布規(guī)律及主要影響因素。

圖9所示為三個攻角荷蘭滾模態(tài)的分布，紅色區(qū)域表示模態(tài)發(fā)散，綠色區(qū)域表示模態(tài)收斂?？梢姡煌ソ窍潞商m滾模態(tài)呈現(xiàn)出明顯不同的發(fā)散/收斂特性。

3.1 荷蘭滾阻尼

3.1.1 模態(tài)簡化

我們根據(jù)各個模態(tài)特征根的特點，對特征方程進(jìn)行分解以獲得各個模態(tài)特征根的近似表達(dá)式。將式(3)中特征矩陣(用A表示)的特征方程展開可得：

det(sI-A)=s4+a3s3+a2s2+a1s+a0

(5)

(a) α=0°

(b) α=5°

(c) α=10°

圖7 橫向靜導(dǎo)數(shù)Lβ分布

Fig.7 Distribution ofLβ

(a) α=0°

(b) α=5°

(c) α=10°

圖8 航向靜導(dǎo)數(shù)Nβ分布

Fig.8 Distribution ofNβ

(a) α=0°

(b) α=5°

(c) α=10°

圖9 荷蘭滾模態(tài)分布

Fig.9 Distribution of Dutch roll

其中:

a0= (LβNr-LrNβ)g/V·cosθ+

(LpNβ-LβNp)g/V·sinθ

a1= (LβNp-LpNβ-Lβg/V)cosα+

(LβNr-LrNβ-Nβg/V)sinα

a2=Nβcosα-Lβsinα+LpNr-

(6)

式中ξ和ωn分別表示荷蘭滾阻尼和無阻尼自振頻率，λR和λS分別為滾轉(zhuǎn)和螺旋模態(tài)特征根。

對于高超聲速乘波體飛行器來說，荷蘭滾模態(tài)特征根的虛部遠(yuǎn)大于其實部以及滾轉(zhuǎn)與螺旋模態(tài)的實根大小，因而可以只保留特征方程的二階及更低階數(shù)項來獲得滾轉(zhuǎn)和螺旋模態(tài)的近似表達(dá)式，然后將式(5)展開，根據(jù)對應(yīng)項相等即可得

(7)

進(jìn)而可得荷蘭滾實部η和虛部ω的表達(dá)式

(8)

(9)

荷蘭滾實部表征荷蘭滾的阻尼大小，通過式(8)所求η幾乎與精確值完全相等。但是，由于a1、a2和a3所包含的物理量較多，通過式(8)仍然難以直觀地獲得定性的結(jié)論。因此，需要對a1、a2和a3進(jìn)一步簡化，以損失一定程度的精度為代價來獲得更簡潔的表達(dá)式，進(jìn)而分析荷蘭滾阻尼的主要影響因素。

以10°攻角時點(0.4,45°)所對應(yīng)的冪次乘波體為例，其靜/動導(dǎo)數(shù)如表3所列。

對于a1，

兩項值的大小相差很大，忽略第二項，即：

(10)

對于a2，Lβ和Nβ遠(yuǎn)大于其他導(dǎo)數(shù)大小，因而可簡化為：

(11)

對于a3，Yβ+Nr與Lp相差接近一個數(shù)量級，可做如下近似：

(12)

根據(jù)上述簡化，可得：

(13)

其中

(14)

Nβ DYN即為有量綱形式的荷蘭滾靜穩(wěn)定性導(dǎo)數(shù)[23]。參考荷蘭滾靜穩(wěn)定性導(dǎo)數(shù)的定義，此處定義一個荷蘭滾動穩(wěn)定性導(dǎo)數(shù)的概念——NpDYN，其表達(dá)式為：

(15)

則有

(16)

根據(jù)式(16)得到的三個攻角下近似荷蘭滾阻尼如圖10所示，根據(jù)特征根得到的精確荷蘭滾阻尼如圖11所示。從圖中可以看到，三個攻角下的荷蘭滾阻尼近似解與精確解的分布趨勢幾乎完全一致，僅僅在定量上存在較小的誤差，因而通過簡化后的式(16)可以較準(zhǔn)確的對荷蘭滾阻尼進(jìn)行定性分析。

3.1.2 荷蘭滾模態(tài)邊界判斷

通過式(16)可以直接判斷荷蘭滾模態(tài)是否穩(wěn)定，而對于具有橫航向靜穩(wěn)定性的乘波體，荷蘭滾模態(tài)收斂的近似條件可進(jìn)一步簡化為：

(a) α=0°

(b) α=5°

(c) α=10°

圖10 近似荷蘭滾阻尼分布

Fig.10 Distribution of approximate Dutch roll damping

(a) α=0°

(b) α=5°

(c) α=10°

圖11 荷蘭滾阻尼分布

Fig.11 Distribution of the Dutch roll damping

(17)

荷蘭滾模態(tài)邊界如圖12所示。圖中“DR” 表示通過求解特征矩陣獲得的與圖9對應(yīng)的荷蘭滾精確邊界，“DR_Approx.”表示通過式(16)得到的荷蘭滾近似邊界，“A”表示攻角。可見，每個攻角下的荷蘭滾模態(tài)近似邊界與精確邊界都比較接近，表明通過式(16)來近似判斷荷蘭滾模態(tài)的發(fā)散/收斂特性是可靠的。

圖12 荷蘭滾模態(tài)發(fā)散/收斂邊界Fig.12 Boundaries of the Dutch roll mode

接下來根據(jù)式(16)來分析三個攻角荷蘭滾模態(tài)發(fā)散/收斂的原因。荷蘭滾動穩(wěn)定性導(dǎo)數(shù)NpDYN在三個攻角下設(shè)計參數(shù)空間內(nèi)的分布如圖13所示。對于Ma=15、H=50 km的飛行工況，g/V=0.002。

當(dāng)α=0°時，整個區(qū)域都有NpDYN0(因為Nβ>0)，此時荷蘭滾模態(tài)邊界由Lβ的正負(fù)決定：當(dāng)Lβ>0時，荷蘭滾模態(tài)收斂；反之，荷蘭滾模態(tài)發(fā)散。由圖7(a)可知，圖中虛線左下角區(qū)域橫向靜不穩(wěn)定，即Lβ?0，因而該部分區(qū)域荷蘭滾模態(tài)收斂，其余部分區(qū)域荷蘭滾模態(tài)發(fā)散。 當(dāng)α=5°時，整個區(qū)域橫航向都是靜穩(wěn)定的，即Lβ/NβDYN<0。在靠近右上角的部分區(qū)域，NpDYN>gcosα/V，因而此區(qū)域荷蘭滾模態(tài)收斂，其余區(qū)域荷蘭滾模態(tài)發(fā)散。 當(dāng)α=10°時，整個區(qū)域同樣為橫航向靜穩(wěn)定的，靠近左下角的小部分區(qū)域NpDYN較小，有NpDYN

3.1.3 荷蘭滾阻尼分布規(guī)律

對荷蘭滾模態(tài)來說，當(dāng)特征根實部為負(fù)且絕對值越大時，飛行器受擾后越能夠快速地回復(fù)到平衡狀態(tài)，荷蘭滾阻尼特性越好，反之則越差。

在0°攻角時，由圖13和圖14(Lβ/NβDYN在參數(shù)設(shè)計空間內(nèi)的分布)可知，Lβ/NβDYN與NpDYN兩項在整個參數(shù)空間內(nèi)的變化范圍都比較大，二者共同作用導(dǎo)致荷蘭滾阻尼特性越靠近左下角越好，越靠近右上

(a)α=0° (b)α=5° (c)α=10°

圖13NpDYN分布角越差，如圖11所示。

Fig.13 Distribution ofNpDYN

(a)α=0°

(b)α=5°

(c)α=10°

圖14Lβ/Nβ DYN分布

Fig.14 Distribution ofLβ/Nβ DYN

在5°和10°攻角時，參數(shù)kw和φ越大，荷蘭滾阻尼特性越好，而且每個點的阻尼也隨攻角增大而增大，其原因可以通過式(16)來解釋。在兩個攻角下，整個參數(shù)空間的Lβ/Nβ DYN的相對變化范圍較小，阻尼特性的變化主要由NpDYN決定：參數(shù)kw和φ越大，NpDYN越大，荷蘭滾阻尼特性越好；攻角越大，NpDYN越大，因而荷蘭滾阻尼特性越好。

根據(jù)上述分析，對于高超聲速乘波體飛行器來說，荷蘭滾動穩(wěn)定性導(dǎo)數(shù)NpDYN對荷蘭滾模態(tài)的阻尼特性起著關(guān)鍵作用。在飛行器初步設(shè)計與控制律設(shè)計時，應(yīng)著重考慮NpDYN的影響。

3.1.4 橫航向耦合運動的CFD/RBD數(shù)值模擬

采用計算流體動力學(xué)方程和剛體動力學(xué)方程(CFD/RBD)的耦合求解來模擬受擾后的橫航向方程可參考文獻(xiàn)[21]，由于此處只針對橫航向耦合運動進(jìn)行數(shù)值模擬，因而不考慮與飛行器縱向運動相關(guān)的自由度，只考慮與橫航向運動相關(guān)的自由度，即只保留繞x軸滾轉(zhuǎn)、繞z軸偏航和沿y軸平動三個自由度上的位移和氣動力。非定常計算的物理時間步長取0.0004 s。

以10°攻角為例，選取兩個典型外形，對應(yīng)的(kw,φ)分別為：外形A(0.3,35°)，外形B(0.3,50°)。根據(jù)圖9和圖12(c)可知，兩個外形對應(yīng)的荷蘭滾模態(tài)分別為發(fā)散和收斂。兩個外形受擾后的橫航向耦合運動的CFD/RBD數(shù)值模擬結(jié)果(以滾轉(zhuǎn)角速度p為例)如圖15所示。由于荷蘭滾模態(tài)特征根實部很小，因而荷蘭滾模態(tài)的發(fā)散或者收斂速度較為緩慢，但仍然可以看出外形A呈現(xiàn)出荷蘭滾模態(tài)發(fā)散特性，而外形B呈現(xiàn)出荷蘭滾模態(tài)收斂特性，與小擾動理論和判據(jù)式(17)的結(jié)果相吻合。

3.2 荷蘭滾頻率

荷蘭滾虛部ω表征荷蘭滾頻率，由于其值明顯大于η、λR和λS，因而可忽略特征矩陣中的重力項和阻尼項，根據(jù)式(5)可得：

(18)

此時兩個零根對應(yīng)滾轉(zhuǎn)模態(tài)和螺旋模態(tài)[16]，則根據(jù)對應(yīng)項相等可得：

(19)

圖16所示為該式計算所得的荷蘭滾頻率相對誤差，可見誤差幾乎為0。由此可知，荷蘭滾頻率特性由荷蘭滾靜穩(wěn)定性導(dǎo)數(shù)Nβ DYN決定。荷蘭滾頻率和Nβ DYN分布分別如圖17和圖18所示。從圖中可見，同一攻角下，kw和φ越大，荷蘭滾頻率越大；隨攻角增大，每個參數(shù)點的荷蘭滾頻率也越大。因而，在飛行攻角較大時，乘波體會面臨荷蘭滾振蕩速度較快的不利現(xiàn)象，需要在飛行器設(shè)計初期合理選擇外形參數(shù)和適當(dāng)?shù)目刂坡稍O(shè)計來解決這一飛行品質(zhì)問題。

(a) 外形A

(b) 外形B圖15 CFD/RBD數(shù)值模擬結(jié)果，α=10°Fig.15 Simulation results of CFD/RBD, α=10°

(a)α=0°

(b)α=5°

(c)α=10°

圖16 荷蘭滾頻率誤差(%)

Fig.16 Error of the Dutch roll frequency

(a)α=0°

(b)α=5°

(c)α=10°

圖17 荷蘭滾頻率分布

Fig.17 Distribution of the Dutch roll frequency

(a)α=0°

(b)α=5°

(c)α=10°

圖18Nβ DYN分布

Fig.18 Distribution ofNβ DYN

4 結(jié) 論

橫航向靜/動穩(wěn)定性分析是飛行器設(shè)計中的關(guān)鍵一步，本文引入?yún)?shù)kw和φ對冪次乘波體外形進(jìn)行參數(shù)化處理，結(jié)合CFD計算和代理模型得到了三個攻角下橫航向靜穩(wěn)定性和荷蘭滾模態(tài)隨參數(shù)的分布規(guī)律，主要內(nèi)容和結(jié)論如下：

1) 乘波體橫向靜穩(wěn)定性主要受上反角影響，隨著設(shè)計參數(shù)kw和φ的增大，上反角變大，橫向靜穩(wěn)定性增強；乘波體的外形特點使其容易設(shè)計成航向靜穩(wěn)定的。

2) 得到了三個攻角下的荷蘭滾模態(tài)分布，不同攻角下呈現(xiàn)出明顯不同的發(fā)散/收斂特性。

3) 推導(dǎo)了荷蘭滾阻尼的近似表達(dá)式，定義了一個新的概念——荷蘭滾動穩(wěn)定性導(dǎo)數(shù)NpDYN。分析表明，NpDYN對荷蘭滾阻尼隨參數(shù)和攻角的變化規(guī)律起著主導(dǎo)作用；根據(jù)阻尼近似表達(dá)式，可以清晰地分析出荷蘭滾模態(tài)發(fā)散/收斂的主要成因。

4) 荷蘭滾靜穩(wěn)定性導(dǎo)數(shù)Nβ DYN對荷蘭滾頻率起決定性作用。設(shè)計參數(shù)kw和φ越大，荷蘭滾頻率越大；攻角越大，荷蘭滾頻率也越大，因而在較大攻角時會面臨荷蘭滾振蕩較快的不利現(xiàn)象。

5) 在實際工程應(yīng)用中，需要對原始乘波體進(jìn)行局部修型以滿足工程需求來獲得最終的高超聲速滑翔飛行器方案，但是經(jīng)過合理得修型，對原有的乘波氣動布局特征不會產(chǎn)生較大的影響。因而，本文關(guān)于荷蘭滾阻尼和頻率的研究結(jié)論對乘波布局高超聲速滑翔飛行器的開環(huán)穩(wěn)定性設(shè)計、評估及閉環(huán)控制律設(shè)計有一定參考價值。

[1]DARPA concludes review of falcon HTV-2 flight anomaly[OL].http://www.darpa.mil/WorkAera/DownloadAsset.aspx?id=2147484134[2]Engineering review board concludes review of HTV-2 second test flight[OL].http://www.darpa.mil/NewsEvents/Releases/2012/04/20.aspx[3]Nonweiler T R F.Aerodynamic problems of manned space vehicles[J].Journal of Royal Aeronautical Society, 1959, 63(585): 521-528.

[4]Zhu L G, Wang Y F, Zhuang F G, et al.The derivation , development of weissman chart and applications on configuration design of reentry vehicle[J].Journal of Astronautics, 2009, 30(1): 13-17.(in Chinese)祝立國, 王永豐, 莊逢甘, 等.Weissman圖的產(chǎn)生、發(fā)展及其在再入航天飛行器氣動布局設(shè)計中的應(yīng)用[J].宇航學(xué)報, 2009, 30(1): 13-17.

[5]Gao Q, Li Q, Chen N, et al.The influence of asymmetric transition on stability of hypersonic aircrafts[J].Tactical Missile Technology, 2012, (6): 12-15.(in Chinese)高清, 李潛, 陳農(nóng), 等.高超聲速飛行器非對稱轉(zhuǎn)捩對穩(wěn)定性的影響[J].戰(zhàn)術(shù)導(dǎo)彈技術(shù), 2012, (6): 12-15.

[6]Gao Q, Li Q.The lateral stability of hypersonic aircraft[J].Aerodynamic Missiles Journal, 2012, (12): 14-18.高清, 李潛.美國高超聲速飛行器橫側(cè)向穩(wěn)定性研究[J].飛航導(dǎo)彈, 2012, (12): 14-18.

[7]Nelson R C.Flight stability and automatic control[M].New York: McGraw-Hill, 1990.

[8]Livneh R.Improved literal approximation for lateral-directional

dynamics of rigid aircraft[J].Journal of Guidance, Control, and Dynamics, 1995, 18(4): 925-927.

[9]Lutze F H, Durham W C, Mason W H.Unified development of lateral-directional departure criteria[J].Journal of Guidance, Control, and Dynamics, 1996, 19(2): 489-493.

[10]Schmidt L V.Introduction to aircraft flight dynamics[M].New York: AIAA Education Series, AIAA, 1998.

[11]Chakravarty M K.Prediction, modeling, and mechanism of aircraft wing rock[D].Indian Inst.Of Technology.Bombay, India, 1999.

[12]Ananthkrishnan N, Unnikrishnan S.Literal approximations to aircraft dynamic modes[J].Journal of Guidance, Control, and Dynamics, 2001, 24(6): 1196-1203.

[13]Phillips W F.Improved closed form approximation for dutch roll[J].Journal of Aircraft, 2000, 37(3): 484-490.

[14]Heller M, Holzapfel F, Sachs G.Robust lateral control of hypersonic vehicles[R].AIAA 2000-4248, 2000.

[15]Breitsamter C, Cvrlje T, Laschka B, et al.Lateral-directional coupling and unsteady aerodynamic effects of hypersonic vehicles[J].Journal of Spacecraft and Rockets, 2001, 38(2): 159-167.

[16]Yin L L, Huang Y, Sun C Z, et al.Improved dutch roll approximation for hypersonic vehicle[J].Sensors & Transducers, 2014, 173(6): 28-33.

[17]賈子安, 張陳安, 王柯穆, 等.乘波布局高超聲速飛行器縱向靜穩(wěn)定特性分析[J].中國科學(xué): 技術(shù)科學(xué), 2014, 44(10): 1114-1122.

[18]Boer A, Schoot M S, Faculty H B.Mesh deformation based on radial basis function interpolation[J].Computers & Structures, 2007, 85: 784-795.

[19]Liu W, Zhang C A, Han H Q, et al.Local piston theory with viscous correction and its application[J].AIAA Journal, 2017, 55(3): 942-954.

[20]Kou J Q, Zhang W W.An approach to enhance the generalization capability of nonlinear aerodynamic reduced-order models[J].Aerospace Science and Technology, 2015, 49: 197-208.

[21]方振平, 陳萬春, 張曙光.航空飛行器飛行動力學(xué)[M].北京: 北京航空航天大學(xué)出版社, 2005.

[22]Liu W.Nonlinear dynamics analysis for mechanism of slender wing rock and study of numerical simulation method[D].National University of Defense Technology, 2004.劉偉.細(xì)長機翼搖滾機理的非線性動力學(xué)分析及數(shù)值模擬方法研究[D].國防科學(xué)技術(shù)大學(xué), 2004.

[23]Moul M T, Paulson J.Dynamic lateral behavior of high performance aircraft[R].NACA RML58E16, 1958.

Study on characteristics of Dutch roll mode for hypersonic waverider

Liu Wen1,2, Zhang Chen′an2,*, Wang Famin2

(1.NationalKeyLaboratoryofAerodynamicDesignandResearch,NorthwesternPolytechnicalUniversity,Xi’an710072,China;2.InstituteofMechanics,ChineseAcademyofSciences,Beijing100190,China)

A waverider aircraft faces the problem regarding lateral-directional coupled dynamic stability due to its non-axisymmetric, flat and slender geometry.The study of the lateral-directional stability of a waverider is relatively limited in the current literature, and instructive conclusions can hardly be found.To investigate the lateral-directional stability, a power-law waverider was taken as the objective, whose Dutch roll mode characteristic was studied in detail.First, the design parameterskwand φ were introduced to describe the geometry of the waverider.The static and dynamic derivatives of the waverider for the whole design parameter space were computed by using CFD and Kriging surrogate model.Then, the characteristics of the Dutch roll mode were obtained according to the solution of linearized small-disturbance equations for lateral-directional motions.An approximate expression was derived for the damping of the Dutch roll mode.A new concept of dynamic stability derivative of the Dutch roll mode was defined.The reason for the divergence or convergence of the Dutch roll mode was analyzed with the help of the approximate expression.The distribution of the damping was obtained at different angles of attack.For the frequency of the Dutch roll mode, an approximate expression was derived, and the variation of the frequency was also analyzed followed by the changing design parameters and angles of attack.

waverider; hypersonic flow; dynamic stability; Dutch roll

0258-1825(2017)03-0444-10

2017-01-16;

2016-03-14

國家自然科學(xué)基金(91216205&11272319)

劉文(1988-),男，山東煙臺人，博士，研究方向：高超聲速空氣動力學(xué)和氣動布局設(shè)計.E-mail: 576825638@qq.com

張陳安*(1982-)，男，廣西南寧人，高級工程師，研究方向：高超聲速空氣動力學(xué)和氣動布局設(shè)計.E-mail: zhch_a@imech.ac.cn

劉文, 張陳安, 王發(fā)民.乘波體荷蘭滾模態(tài)特性研究[J].空氣動力學(xué)學(xué)報, 2017, 35(3): 444-453.

10.7638/kqdlxxb-2017.0024 Liu W, Zhang C A, Wang F M.Study on characteristics of Dutch roll mode for hypersonic waverider[J].Acta Aerodynamica Sinica, 2017, 35(3): 444-453.

V212.12

A doi: 10.7638/kqdlxxb-2017.0024