呂 達(dá)，張維桐，張魯民，趙俊波，張石玉，蘇浩秦

(1.中國航天空氣動(dòng)力技術(shù)研究院，北京 100074；2.北京航空航天大學(xué) 宇航學(xué)院，北京 100191；3.國家計(jì)算流體力學(xué)實(shí)驗(yàn)室，北京 100191)

0 引言

隨著打擊時(shí)敏目標(biāo)和機(jī)動(dòng)突防的需求日益增長，具備高升阻比、高機(jī)動(dòng)性能的面對稱升力式構(gòu)型已成為高超聲速飛行器氣動(dòng)布局主流發(fā)展方向。2010 年美國HTV-2、X-51A、X-37B 的試飛[1-4]，標(biāo)志著繼航天飛機(jī)以后的新型面對稱高超聲速飛行技術(shù)的工程化成為可能。同時(shí)，以進(jìn)一步提升氣動(dòng)與飛行性能為設(shè)計(jì)目標(biāo)的無尾布局是高超聲速飛行器未來的重要發(fā)展方向。無尾布局外形簡潔，具備全向?qū)掝l隱身特性好、升阻特性優(yōu)、氣動(dòng)熱載荷低等性能優(yōu)勢，但在通道耦合和考慮耦合的操穩(wěn)設(shè)計(jì)方面存在難題，急需技術(shù)攻關(guān)。主要體現(xiàn)在：對于無尾布局高超聲速飛行器，由于不存在垂尾，俯仰、滾轉(zhuǎn)和偏航三個(gè)通道的穩(wěn)定特性只能依賴機(jī)體的設(shè)計(jì)，而機(jī)體的形態(tài)和幾何參數(shù)又受容積率、升阻比等指標(biāo)的約束，穩(wěn)定性調(diào)節(jié)范圍較小，且其橫航向穩(wěn)定性較差，因此綜合多重指標(biāo)的優(yōu)化設(shè)計(jì)難度極大。同時(shí)，無尾布局飛行器通常采用欠驅(qū)動(dòng)控制，其通道間耦合效應(yīng)較強(qiáng)，耦合操穩(wěn)特性的設(shè)計(jì)難度較大。因此，發(fā)展以三通道判據(jù)為核心的精細(xì)化耦合穩(wěn)定性分析思想對于先進(jìn)的無尾布局高超聲速飛行器的設(shè)計(jì)具有重要意義。

20 世紀(jì)50 年代中期以來，科技工作者對飛行器耦合失穩(wěn)開展了大量研究。文獻(xiàn)[5]系統(tǒng)總結(jié)了各類判據(jù)，如：Moul、Weissman、Kalviste、Johnston、Bihrle、Pelikan 等的研究成果。歷經(jīng)多年發(fā)展，如今上述各類判據(jù)已經(jīng)成為航空航天飛行器穩(wěn)定性分析的重要工具。文獻(xiàn)[6]詳細(xì)總結(jié)了飛行器耦合失穩(wěn)各類判據(jù)在現(xiàn)階段國內(nèi)外飛行器設(shè)計(jì)工作中的應(yīng)用。

近年來，隨著無尾布局高超聲速飛行器技術(shù)的發(fā)展，國內(nèi)外對于耦合穩(wěn)定性的研究進(jìn)一步深入。例如，Park 等[7]使用一系列氣動(dòng)分析方法評估了無尾飛行器的穩(wěn)定性和控制特性，利用縱向和橫向運(yùn)動(dòng)的系統(tǒng)矩陣推導(dǎo)了無尾飛行器的穩(wěn)定性準(zhǔn)則，對幾種飛行條件進(jìn)行了穩(wěn)定性分析。Fu 等[8-10]提出了一種飛行試驗(yàn)，用以研究具備荷蘭橫滾模式穩(wěn)定性但不具備偏航穩(wěn)定性的飛行器的耦合穩(wěn)定特性。Shen[11-12]為了準(zhǔn)確預(yù)測偏離行為，在控制系統(tǒng)的設(shè)計(jì)過程中考慮偏航-滾轉(zhuǎn)耦合的影響，并進(jìn)行了不同耦合比下的偏航-滾轉(zhuǎn)耦合風(fēng)洞試驗(yàn)。此外，還有許多文獻(xiàn)在飛行器氣動(dòng)設(shè)計(jì)[13-16]和飛行器控制系統(tǒng)設(shè)計(jì)[17-23]的過程中考慮了耦合穩(wěn)定性偏離特性的影響。

本文通過理論推導(dǎo)以及數(shù)值仿真試驗(yàn)，引入了開環(huán)及閉環(huán)情況下的3 個(gè)靜穩(wěn)定性耦合判據(jù)和3 個(gè)動(dòng)穩(wěn)定性耦合判據(jù)。與傳統(tǒng)的判據(jù)相比，更為全面地考慮了飛行器三通道耦合對其穩(wěn)定性的影響，可用于精細(xì)化地分析近代先進(jìn)無尾布局高超聲速飛行器三通道的耦合穩(wěn)定特性，且對其設(shè)計(jì)過程具有指導(dǎo)意義。

1 飛行器的姿態(tài)動(dòng)力學(xué)方程組

利用剛體假設(shè)建立飛行器的六自由度動(dòng)力學(xué)模型，忽略地球自轉(zhuǎn)的影響，根據(jù)動(dòng)量、動(dòng)量矩定理，可以導(dǎo)出飛行器機(jī)體坐標(biāo)系下的六自由度動(dòng)力學(xué)方程組[24]。

飛行器運(yùn)動(dòng)可分解為平動(dòng)和轉(zhuǎn)動(dòng)，而確保轉(zhuǎn)動(dòng)運(yùn)動(dòng)穩(wěn)定性是保證飛行器飛行穩(wěn)定性的最根本條件。根據(jù)飛行器的六自由度動(dòng)力學(xué)方程組，忽略平動(dòng)速度影響，可得：

同時(shí)，對于角加速度方程組，由小擾動(dòng)線性化方法，可得到角加速度方程組：

飛機(jī)相對于速度矢量的方向可由3 個(gè)角度定義，其中兩個(gè)分別為攻角和側(cè)滑角 (α,β)，現(xiàn)增加側(cè)傾角μ。與體軸系滾轉(zhuǎn)角繞體軸ox轉(zhuǎn)動(dòng)的定義不同，側(cè)傾角 μ定義的是飛行器的升力矢量相對于速度矢量的旋轉(zhuǎn)角度，具體見文獻(xiàn)[24]，側(cè)傾角速率與體坐標(biāo)系中角速率的關(guān)系為：

略去位移運(yùn)動(dòng)的影響，略去慣性積，由小擾動(dòng)線性化方法，得到簡化的飛行器姿態(tài)動(dòng)力學(xué)方程組為：

式中的算子定義如下：

2 開環(huán)三通道靜動(dòng)穩(wěn)定性耦合判據(jù)

2.1 俯仰靜動(dòng)穩(wěn)定性耦合判據(jù)

利用小擾動(dòng)線性化方法，選擇式（4）中第一個(gè)方程和式（5）組成俯仰通道系統(tǒng)的線化動(dòng)力學(xué)模型：

對于該類飛行器，其橫航向阻尼較小，僅考慮由俯仰通道 Δα和 Δq產(chǎn)生的耦合力矩，忽略式（7）中的滾轉(zhuǎn)和偏航阻尼項(xiàng)，可得簡化的系統(tǒng)矩陣如下：

則其特征方程為：

由于忽略了橫航向的阻尼，故特征方程失去了一個(gè)根，則分別由勞斯判據(jù)的c＞0 和bc-ad＞0，有俯仰靜穩(wěn)定性開環(huán)耦合判據(jù)和俯仰動(dòng)穩(wěn)定性開環(huán)耦合判據(jù)：

進(jìn)一步得其無量綱的形式為：

2.2 偏航靜動(dòng)穩(wěn)定性耦合判據(jù)

選擇式（4）中第二個(gè)方程以及式（5）組成偏航通道系統(tǒng)的線化動(dòng)力學(xué)模型：

僅考慮由偏航通道 Δβ 和 Δr產(chǎn)生的耦合力矩，故忽略式（12）中的滾轉(zhuǎn)阻尼項(xiàng)，保留俯仰阻尼系數(shù)，可得偏航通道簡化的系統(tǒng)矩陣如下：

其特征方程為：

分別由勞斯判據(jù)的c＞0 和bc-ad＞0，從而有偏航靜穩(wěn)定性開環(huán)耦合判據(jù)和偏航動(dòng)穩(wěn)定性開環(huán)耦合判據(jù)：

無量綱化得到下式，式中Q為動(dòng)壓，S為參考面積，l為參考長度。

2.3 滾轉(zhuǎn)靜動(dòng)穩(wěn)定性耦合判據(jù)

選擇式（4）中第三個(gè)方程以及式（5）組成滾轉(zhuǎn)通道系統(tǒng)的線化動(dòng)力學(xué)模型。此時(shí)由于后三個(gè)方程不顯式地含有μ，但μ變化時(shí)會(huì)同時(shí)引起β的變化，μ對力矩的影響可以通過力矩對于β的導(dǎo)數(shù)來體現(xiàn)，故應(yīng)尋找在滾轉(zhuǎn)通道運(yùn)動(dòng)下μ和β的近似關(guān)系。

由于面對稱的無尾布局高超聲速飛行器的側(cè)滑角通常希望被控制在0°，圖1 為當(dāng)β接近于0°時(shí)由于飛行器滾轉(zhuǎn)運(yùn)動(dòng)帶來的β及μ變化的近似關(guān)系。圖中OAC平面為飛行器受擾前的對稱面，OAB為受擾后的對稱面，OXs為受擾后的穩(wěn)定軸系的X軸，A為體軸系上一點(diǎn)，AB為過點(diǎn)A向OXs軸做的垂線，AC為過點(diǎn)A向風(fēng)軸系X軸做的垂線，由于OXs為穩(wěn)定軸系的X軸，故OAB平面與OBC平面垂直，且AB垂直于OXs，故AB垂直于OBC平面，從而AB垂直于OC，又由于AC垂直于OC，且AB和AC又分別在飛行器受擾后與受擾前的對稱面中，故AB和AC可以認(rèn)為是飛行器的升力方向，則根據(jù)側(cè)傾角的定義，AB和AC之間的夾角即為側(cè)傾角μ，故由圖中幾何關(guān)系有：

圖1 當(dāng)β 接近于0 時(shí)由于滾轉(zhuǎn)通道運(yùn)動(dòng)帶來的β 及μ 變化的近似關(guān)系Fig.1 Approximate relationship between variations of β and μ due to roll channel motion when β is close to 0

即當(dāng)研究飛行器的滾轉(zhuǎn)通道運(yùn)動(dòng)時(shí)，有如下近似關(guān)系：

式（19）即μ的力矩系數(shù)的近似表達(dá)，注意，μ的力矩系數(shù)在實(shí)際中并不存在，本文只是利用了飛行器滾轉(zhuǎn)運(yùn)動(dòng)時(shí)μ和β的近似關(guān)系導(dǎo)出了由于μ的變化而造成了β的伴隨變化而導(dǎo)致的力矩改變。

則滾轉(zhuǎn)通道系統(tǒng)的線化動(dòng)力學(xué)模型近似為：

與其他兩通道類似，僅考慮由滾轉(zhuǎn)通道 Δμ和Δp產(chǎn)生的耦合力矩，忽略式（20）中的偏航阻尼項(xiàng)，保留俯仰阻尼系數(shù)，可得簡化的系統(tǒng)矩陣如下：

分別由勞斯判據(jù)的c＞0 和bc-ad＞0，有滾轉(zhuǎn)靜穩(wěn)定性開環(huán)耦合判據(jù)和滾轉(zhuǎn)動(dòng)穩(wěn)定性開環(huán)耦合判據(jù)：

其無量綱簡化的形式如下：

2.4 本文判據(jù)與現(xiàn)有判據(jù)的對比闡釋

與傳統(tǒng)的氣動(dòng)穩(wěn)定性判據(jù)相比，本文判據(jù)在攻角為正，且攻角和側(cè)滑角接近于零時(shí)，可退化為前者，即ClβDYN=Clβ，CmαDYN=Cmα，CnβDYN=Cnβ。此外，若忽略本文判據(jù)中氣動(dòng)力矩系數(shù)的小量，則本文判據(jù)同樣可退化為傳統(tǒng)判據(jù)。與傳統(tǒng)判據(jù)的臨界值皆為0 不同，由于本文考慮了俯仰通道動(dòng)導(dǎo)數(shù)對橫航向的影響，故橫航向判據(jù)的臨界值不為0。因此，本文所述判據(jù)與傳統(tǒng)的氣動(dòng)穩(wěn)定性判據(jù)仍然有著不可分割的聯(lián)系，本文所述判據(jù)是現(xiàn)有判據(jù)的進(jìn)一步完善，對于穩(wěn)定性可調(diào)節(jié)范圍較小的無尾布局高速飛行器的氣動(dòng)穩(wěn)定性設(shè)計(jì)具有指導(dǎo)意義。

3 閉環(huán)三通道靜動(dòng)穩(wěn)定性耦合判據(jù)

3.1 俯仰閉環(huán)靜動(dòng)穩(wěn)定性耦合判據(jù)

采用升降舵作為俯仰通道控制力來源。考慮攻角及俯仰角速度鎮(zhèn)定，利用攻角及俯仰角速度反饋至升降舵的方式，反饋控制律采用以下形式：

式中 δe為升降舵的舵偏角，則有閉環(huán)的俯仰系統(tǒng)矩陣為：

由勞斯判據(jù)c＞0，有俯仰通道閉環(huán)靜穩(wěn)定性耦合判據(jù)為：

3.2 偏航閉環(huán)靜動(dòng)穩(wěn)定性耦合判據(jù)

采用副翼作為偏航方向控制力來源?？紤]鎮(zhèn)定側(cè)滑角及偏航角速度，即實(shí)施側(cè)滑角及偏航角速度反饋至副翼，反饋控制律采用以下形式：

式中，δa為副翼的差動(dòng)偏角，忽略副翼引起的側(cè)向力以及俯仰力矩的變化，則有閉環(huán)的偏航系統(tǒng)矩陣為：

則其閉環(huán)特征方程為：

由勞斯判據(jù)c＞0，有偏航通道閉環(huán)靜穩(wěn)定性耦合判據(jù)為：

3.3 滾轉(zhuǎn)閉環(huán)靜動(dòng)穩(wěn)定性耦合判據(jù)

同樣采用副翼作為滾轉(zhuǎn)通道控制力來源?？紤]鎮(zhèn)定側(cè)傾角μ及滾轉(zhuǎn)角速度p，即實(shí)施側(cè)傾角及滾轉(zhuǎn)角速度反饋至副翼 δa，反饋控制律采用以下形式：

式中，δa為副翼的差動(dòng)偏角，忽略副翼引起的俯仰力矩的變化，則有閉環(huán)的滾轉(zhuǎn)系統(tǒng)矩陣為：

由勞斯判據(jù)c＞0，有滾轉(zhuǎn)靜穩(wěn)定性閉環(huán)耦合判據(jù)：

其無量綱簡化的形式為：

4 三通道穩(wěn)定性耦合判據(jù)數(shù)值仿真驗(yàn)證

本節(jié)通過一個(gè)數(shù)值仿真案例來驗(yàn)證所提出的三通道穩(wěn)定性耦合判據(jù)的有效性。采用某型面對稱無尾布局高超聲速飛行器作為研究對象[15]，其基礎(chǔ)外形1 如圖2 所示，將基礎(chǔ)外形的前緣后掠角增大，且將飛行器底部截面高度增加，側(cè)緣切削量增大，得到外形 2；將基礎(chǔ)外形的前緣后掠角減小，并將飛行器底部截面高度增加，且將其質(zhì)心適當(dāng)向后調(diào)整，得到外形3；最后將基礎(chǔ)外形的前緣后掠角減小并適當(dāng)飛行器底部截面高度，得到外形 4。

圖2 基礎(chǔ)外形的三視圖[15]Fig.2 Three views of basic configuration[15]

首先選取基礎(chǔ)外形1，對其氣動(dòng)穩(wěn)定性進(jìn)行評估。圖3 為計(jì)算域以及對稱面的網(wǎng)格劃分，采用結(jié)構(gòu)網(wǎng)格，壁面第一層網(wǎng)格高度設(shè)置為0.001 mm，以保證y+值小于1。采用基于三維積分形式的雷諾平均N-S方程求解，湍流模型采用S-A 模型，無黏通量采用差分分裂的Roe 格式，黏性項(xiàng)采用二階中心差分格式進(jìn)行離散，時(shí)間推進(jìn)格式采用隱式近似因子分解法。選取類似外形的少量試驗(yàn)數(shù)據(jù)對該計(jì)算方法進(jìn)行驗(yàn)證，風(fēng)洞試驗(yàn)?zāi)Ｐ腿鐖D4[25]所示，采用該方法得到的結(jié)果與試驗(yàn)得到的結(jié)果對比如圖5，可見兩者結(jié)果在量值及趨勢上吻合較好。其他三種外形類似，不再贅述。

圖4 風(fēng)洞試驗(yàn)?zāi)Ｐ蚚25]Fig.4 Wind tunnel test model[25]

圖5 數(shù)值計(jì)算及風(fēng)洞試驗(yàn)結(jié)果對比Fig.5 Comparison of CFD and experimental results

利用CFD 計(jì)算得到的氣動(dòng)數(shù)據(jù)進(jìn)行六自由度動(dòng)力學(xué)仿真。對于基礎(chǔ)外形1，其在攻角2°時(shí)由CFD 仿真得到的核心氣動(dòng)穩(wěn)定性數(shù)據(jù)以及其他三種外形的核心氣動(dòng)穩(wěn)定性數(shù)據(jù)如表1 所示。為了方便判斷其偏離特性，仿真的判據(jù)臨界值近似取0，可見外形1 不同時(shí)滿足三通道穩(wěn)定性耦合判據(jù)，由于沒有垂尾，導(dǎo)致CnβDYN＜0，偏航通道不穩(wěn)定，且ClβDYN的絕對值偏小，即滾轉(zhuǎn)通道穩(wěn)定性偏小。取初始配平狀態(tài)：α0=2°，β0=0°，μ0=0°，p0=0(°)/s，q0=0(°)/s，r0=0(°)/s，假設(shè)飛行器的滾轉(zhuǎn)角速度在仿真初始時(shí)刻受到擾動(dòng)變?yōu)閜=p0+Δp=1(°)/s，則可得側(cè)傾角、攻角和側(cè)滑角以及姿態(tài)角速度的變化曲線如圖6 和圖7 所示。由圖可得，飛行器的側(cè)滑角以及側(cè)傾角迅速發(fā)散，攻角在一段時(shí)間之后也發(fā)散，飛行器處于不穩(wěn)定狀態(tài)。

表1 本文仿真涉及到的4 種外形在攻角2°時(shí)的核心氣動(dòng)穩(wěn)定性數(shù)據(jù)Table 1 Core aerodynamic data at 2° angle of attack of the four configurations involved in the simulation

圖6 飛行器的姿態(tài)角速度的變化曲線Fig.6 Variation of aircraft angular velocity

圖7 飛行器的側(cè)傾角、攻角和側(cè)滑角的變化曲線Fig.7 Variations of aircraft roll angle,angle of attack,and sideslip angle

對于外形2，使其在同樣的初始攻角情況下滿足CnβDYN＞0，CmαDYN＜0，但ClβDYN＞0，如表1 所示，即滾轉(zhuǎn)通道不穩(wěn)定，偏航通道弱穩(wěn)定，給定初始狀態(tài)α0=2°，β0=0°，μ0=0°，p0=0(°)/s，q0=0(°)/s，r0=0(°)/s，在仿真初始時(shí)刻給相同的滾轉(zhuǎn)角速度擾動(dòng)Δp=1(°)/s，仿真結(jié)果如圖8 和圖9所示。由圖可得，飛行器側(cè)傾角迅速發(fā)散，且由于耦合效應(yīng)，攻角以及側(cè)滑角隨后也發(fā)散，飛行器處于不穩(wěn)定狀態(tài)。

圖8 飛行器的姿態(tài)角速度的變化曲線Fig.8 Variation of aircraft angular velocity

圖9 飛行器的側(cè)傾角、攻角和側(cè)滑角的變化曲線Fig.9 Variations of aircraft roll angle,angle of attack,and sideslip angle

對于外形3，使其在同樣的初始攻角情況下滿足CnβDYN＞0、ClβDYN＜0，但CmαDYN＞0，如表1 所示，即俯仰通道不穩(wěn)定，滾轉(zhuǎn)和偏航通道弱穩(wěn)定，給定初始狀態(tài)α0=2°，β0=0°，μ0=0°，p0=0(°)/s，q0=0(°)/s，r0=0(°)/s，仿真初始時(shí)刻給定相同的擾動(dòng) Δp=1(°)/s，仿真結(jié)果如圖10 和圖11 所示。由圖可得，飛行器攻角迅速發(fā)散，且由于耦合效應(yīng)，側(cè)傾角以及側(cè)滑角隨后發(fā)散，飛行器處于不穩(wěn)定狀態(tài)。

圖10 飛行器的姿態(tài)角速度的變化曲線Fig.10 Variation of aircraft angular velocity

圖11 飛行器的側(cè)傾角、攻角和側(cè)滑角的變化曲線Fig.11 Variations of aircraft roll angle,angle of attack,and sideslip angle

對于外形4，使其在同樣的初始攻角情況下滿足CmαDYN＜0、ClβDYN＜0，CnβDYN＞0，如表1 所示，即滿足三通道的穩(wěn)定性耦合判據(jù)。給定初始狀態(tài)α0=2°，β0=0°，μ0=0°，p0=0(°)/s，q0=0(°)/s，r0=0(°)/s，同樣在仿真初始時(shí)刻施加相同的擾動(dòng) Δp=1(°)/s，仿真結(jié)果如圖12 和圖13 所示。由圖可得，飛行器側(cè)傾角、攻角以及側(cè)滑角均收斂到穩(wěn)定的值，飛行器處于穩(wěn)定狀態(tài)。

圖12 飛行器的姿態(tài)角速度的變化曲線Fig.12 Variation of aircraft angular velocity

圖13 飛行器的側(cè)傾角、攻角和側(cè)滑角的變化曲線Fig.13 Variations of aircraft roll angle,angle of attack,and sideslip angle

綜上所述，當(dāng)無尾布局高超聲速飛行器同時(shí)滿足本文所提出的三通道穩(wěn)定性耦合判據(jù)時(shí)，飛行器處于穩(wěn)定狀態(tài)。反之，則會(huì)出現(xiàn)各種耦合失穩(wěn)問題。初步說明本文提出的三通道穩(wěn)定性耦合判據(jù)是有效的。

5 結(jié)論

針對無尾布局高超聲速飛行器飛行中的通道耦合偏離失穩(wěn)問題，本文推導(dǎo)了面對稱無尾布局高超聲速飛行器三通道的氣動(dòng)穩(wěn)定性耦合判據(jù)，并通過幾組數(shù)值仿真算例初步驗(yàn)證了本文判據(jù)的有效性，完善了判斷此類飛行器靜動(dòng)穩(wěn)定性的判別方法。與傳統(tǒng)判據(jù)的臨界值皆為0 不同，由于本文考慮了俯仰通道動(dòng)導(dǎo)數(shù)對橫航向的影響，故橫航向判據(jù)的臨界值不為0，本文判據(jù)相比傳統(tǒng)判據(jù)考慮的因素更加全面。

本文判據(jù)在推導(dǎo)過程中應(yīng)用了小擾動(dòng)線化理論且省略了許多因素，因此，該判據(jù)仍具有一定局限性，未來需要通過數(shù)值仿真、地面試驗(yàn)及飛行試驗(yàn)，進(jìn)一步驗(yàn)證判據(jù)的有效性。同理，對本文建立的閉環(huán)穩(wěn)定性耦合判據(jù)也需要通過進(jìn)一步地研究以驗(yàn)證其可靠性。