超/高超聲速飛行器動態(tài)穩(wěn)定性導數(shù)極快速預測方法

李正洲，高昌，肖天航，馬自成，肖濟良，朱建輝

1. 中國空氣動力研究與發(fā)展中心 高超聲速沖壓發(fā)動機技術重點實驗室，綿陽 621000

2. 南京航空航天大學 航空學院，南京 210016

目前獲取動導數(shù)的主要方法有飛行試驗、風洞試驗[4-5]、數(shù)值計算[6-10]及工程近似方法等。飛行試驗和風洞試驗是飛行器動態(tài)特性判定的主要依據(jù)，但存在難度大、周期長、費用高的特點，需要與數(shù)值計算相互補充、互相驗證；采用非定常CFD數(shù)值模擬方法，可以考慮到流場的非線性特性，便于開展復雜外形的氣動力計算，但主要局限性在于計算量大，難以快速獲得定性的結(jié)論和定量判據(jù)；以牛頓撞擊理論[11]為代表的工程近似法考慮了線化空氣動力學理論和經(jīng)驗關系，具有快捷高效的優(yōu)勢，但對復雜外形的大迎角非線性流動可能在數(shù)值上存在量級甚至符號的差別。

動導數(shù)的高效預測關鍵在于飛行器非定常氣動力的快速、準確計算。在非定常氣動力快速計算方面，近年來的熱門方向是降階模型(Reduced Order Modeling, ROM)[12]，這一方法在計算精度和計算效率上取得了很好的平衡。但由于降階模型仍然部分依賴耦合非定常CFD數(shù)值計算，對非定常數(shù)值求解的魯棒性有較高要求，因此一些學者轉(zhuǎn)而采用定常CFD和工程方法相結(jié)合的思路構(gòu)建非定常氣動力模型，其中最典型的是基于CFD技術的當?shù)亓骰钊碚?Local Piston Theory, LPT)[13]。陳勁松和曹軍[14]、張偉偉[15-16]等基于定常歐拉(Euler)方程計算流場，利用獲得的物面當?shù)亓鲃訁?shù)結(jié)合活塞理論來計算非定常氣動力，這一方法解決了活塞理論只能計算一定馬赫數(shù)范圍內(nèi)小迎角、尖頭薄翼的缺點；劉溢浪[17]、秦之軒[18]等成功地將該方法應用于氣動導數(shù)快速預測。

結(jié)合定常CFD與當?shù)亓骰钊碚摰姆椒▉眍A測動導數(shù)，從結(jié)果來看，能夠較好地揭示超聲速/高超聲速流動下飛行器動導數(shù)的變化規(guī)律，相比于完全非定常CFD時域推進方法可以節(jié)約大量計算時間。然而，該方法仍然依賴定常CFD流場數(shù)據(jù)，無法在飛行器設計早期階段快速地給出定性的結(jié)果。本文針對上述問題，發(fā)展了一種結(jié)合當?shù)乇砻嫘倍确ê彤數(shù)亓骰钊碚摰母咚亠w行器非定常氣動力快速計算方法，再對非定常氣動力進行提取、辨識，可得到動導數(shù)。該方法不依賴CFD流場結(jié)果，在預測動導數(shù)方面具有極高的效率，同時具有較高的精度。

1 動導數(shù)高效預測策略及流程

本文動導數(shù)高效預測策略及流程如圖1所示：① 所需輸入為飛行器物面網(wǎng)格、來流參數(shù)及物面隨時間的變形、振動歷程；② 飛行器的非定常氣動力根據(jù)當?shù)亓骰钊碚摲譃闉槭茏杂蓙砹饕鸬臒o附加擾動項以及飛行器物面變形或運動引起的附加擾動項；③ 采用基于牛頓撞擊理論的當?shù)乇砻嫘倍确ㄇ蟪鰵鈩恿o附加擾動項，再根據(jù)激波后的等熵關系求解當?shù)亓鲃訁?shù)，結(jié)合物面變形或振動可求出當?shù)匚锩嫦孪此俣?，繼而求出氣動力附加擾動項隨時間變化；④ 對氣動力無附加擾動項和擾動項沿物面積分，即可求出高速飛行器受到的非定常氣動力。圖中：Ma、H、α、β分別馬赫數(shù)、高度、迎角、側(cè)滑角；Vl、Vb為飛行器物面的當?shù)厮俣仁噶?；?、αm、k、t分別為強迫旋轉(zhuǎn)運動的初始迎角、迎角振幅、縮減頻率、時刻；h、hm分別為強迫沉浮運動的高度、高度振幅。

圖1 高速飛行器動導數(shù)高效預測流程

本文方法的優(yōu)勢在于：采用當?shù)乇砻嫘倍确凹げê蟮牡褥仃P系式等求解飛行器表面流場參數(shù)，從而克服了傳統(tǒng)方法對CFD流場參數(shù)的依賴和耦合(圖1中虛線所示)。由于該方法不需要CFD方法求解定常流場，因此具有極高的效率；后文的算例也表明該方法具有較高的精度。因此本文方法可作為高速飛行器總體設計階段布局選型的工具。

值得說明的是，雖然本文在預測飛行器氣動導數(shù)過程中采用了基于牛頓撞擊理論的當?shù)乇砻嫘倍确?，但與“采用牛頓撞擊理論直接計算動導數(shù)”的工程近似法，兩者的原理不同：后者是將壓力系數(shù)對角速度(或迎角變化量)求導，再積分求出整個飛行器的動導數(shù)；本文方法是對飛行器施加強迫簡諧運動求出周期非定常氣動力，再對飛行器的非定常氣動力進行提取、辨識氣動導數(shù)。本文對飛行器施加強迫簡諧運動為風洞試驗、數(shù)值模擬計算動導數(shù)通用方法[19]。

2 非定常氣動力快速預測

2.1 當?shù)亓骰钊碚摰耐茖?/h3>

本文基于當?shù)亓骰钊碚搶Τ?高超聲速飛行器非定常氣動力進行快速預測，首先對當?shù)亓骰钊碚撨M行推導。

采用動量定理推導當?shù)亓骰钊碚揫20]，如圖2所示：考慮Ma?1的氣流經(jīng)過飛行器表面，此時擾動可近似為沿物面法向傳播，如同氣缸中活塞隊氣流的擾動一樣。圖中：V、P、a分別為速度、壓強、聲速；下標∞、l、n分別表示自由來流、當?shù)?、法向的流動參?shù)。

假設氣流受到物面變形或振動引起的擾動后速度在dt時間內(nèi)變化dW，擾動傳播的距離為a·dt，則受到擾動氣體質(zhì)量為ρaSdt(ρ為氣體密度，S為活塞面積)，受到擾動氣體的總動量變化為ρaSdt·dW。另外，活塞對氣體的沖量為dP·S·dt。根據(jù)動量定理可得

dP·S·dt=ρaSdt·dW

(1)

式(1)在當?shù)亓鲃訁?shù)下可簡化為

dP=ρlaldW

(2)

對式(2)積分，可得當?shù)亓骰钊碚摰膲毫τ嬎愎綖?/p>

P=ρlalW+C

(3)

式中：常數(shù)C為附加擾動為0時的壓力；ρlalW為物面擾動時的壓力。

將附加擾動為0時的物面壓力用Psteady表示，又物面擾動可分為變形或振動，則當?shù)亓骰钊碚摰膲毫τ嬎愎娇蓪憺?/p>

(4)

式中：n0為物面變形前的外法線單位矢量；n為物面變形后的外法線單位矢量；Vl和Vb分別為當?shù)亓鲃铀俣群臀锩嬲駝铀俣?；物面擾動速度W由物面變形速度Vl·δn和振動速度Vb·n組成。

圖2 當?shù)亓骰钊碚撌疽鈭D

與經(jīng)典活塞理論相比，基于動量定理推導出的當?shù)亓骰钊碚摚孩?推導過程中不需要級數(shù)展開，也就不存在忽略高階項的近似，因此精度比忽略高階項的活塞理論高；② 沒有活塞理論所必須的擾動速度小于來流聲速(W/a∞<1)的要求，因此當?shù)亓骰钊碚撝恍铦M足“擾動可近似為沿物面法向傳播”的假設即可，拓寬了馬赫數(shù)適用范圍的上限；③ 不需要等熵假設，這表明當?shù)亓骰钊碚摬粌H可以用于大迎角、大厚度翼面問題，也可用于三維翼面、機身以及鈍前緣或鈍頭體存在弓形激波的情況。上述特點是當?shù)亓骰钊碚摽赏茝V用于機身和其他復雜外形的基礎[21]。

式(4)表明超/高超聲速飛行器非定常氣動力可分為無擾動項和物面變形/振動后的擾動項。下面分別介紹這兩項氣動力的計算方法。

2.2 氣動力無附加擾動項計算

對于無物面變形或振動擾動的高速飛行器，可看做定常流動，再用基于牛頓撞擊理論的當?shù)乇砻嫘倍确ㄓ嬎銡鈩恿o附加擾動項Psteady。具體方法為[22]：將飛行器表面離散、劃分為三角形面元網(wǎng)格，對每個面元區(qū)分為迎風面或背風面，再分別采取不同的公式進行計算。迎風面、背風面的劃分是根據(jù)自由來流條件及面元法矢共同確定的；此外，在面元判斷時加入了對“遮擋面元”的判斷，即采用光線投射算法(Ray-casting Algorithm)[23]判斷當前面元是否被遮擋。若為遮擋面元，則將該面元作為背風面處理。

對于迎風面，采用活塞理論與修正牛頓理論結(jié)合的方法進行計算[24]：

(5)

式中：B=(2/π)arctan[(Ma2+6)/10]；Cp,max為駐點壓力系數(shù)。

對于背風面，則采用普朗特-邁耶膨脹波方法：

(6)

式中：θ和γ分別為物面偏折角、氣體比熱比。

2.3 氣動力附加擾動項計算

求解氣動力附加擾動項時，需要用到當?shù)亓鲃訁?shù)。為求解方便，本文采用求解邊界層外緣的擾動壓力系數(shù)而不直接求解物面擾動壓力系數(shù)。又根據(jù)邊界層內(nèi)的壓力梯度較小，可近似認為邊界層外緣的擾動壓力等于飛行器表面的壓力，因此所求出有效外形的非定常氣動力即為飛行器的非定常氣動力。

當?shù)亓鲃訁?shù)求解原理如下：對于完全氣體或平衡氣體，所有狀態(tài)參數(shù)(壓力、密度、溫度、焓、聲速、熵、流速、黏性系數(shù))中只有兩個是相互獨立的，其他狀態(tài)參數(shù)均可根據(jù)相應的熱力學關系由這兩個參數(shù)求出。

2.3.1 正激波完全氣體邊界層外緣參數(shù)

對于大鈍頭體飛行器形成的正激波完全氣體，可認為進入邊界層外緣的流線基本是通過弓形激波的正激波部分。由于正激波后的熵在來流條件給定后是唯一確定的，因此，利用所得出的物面定常壓強和正激波后的熵這兩個獨立的狀態(tài)參數(shù)，可以確定其他邊界層外緣參數(shù)。

激波后的氣體壓力P2和密度ρ2可利用正激波前后關系式求得，即

(7)

(8)

式中：P∞為來流靜壓；ρ∞為來流密度。

利用等熵關系，求邊界層外緣密度ρe：

(9)

再用邊界層外緣的壓力和密度求出邊界層外緣的聲速：

(10)

求出邊界層外緣的當?shù)孛芏群吐曀僖院?，結(jié)合物面的變形或振動的下洗速度，即可根據(jù)式(4)計算出氣動力的附加擾動項。

2.3.2 斜激波完全氣體邊界層外緣參數(shù)

對于尖銳前緣、細長體飛行器外形，通過其頭部的激波是一道或數(shù)道斜激波，斜激波完全氣體后的壓力和密度可分別為

(11)

(12)

式中：β為激波角，可由θ-β-Ma圖表查出[25]，也可通過代數(shù)變換后的激波角顯示表達式求出[26]。

求出激波后的氣體壓力P2和密度ρ2后，再采用正激波完全氣體相同的方法求出邊界外緣的密度ρe和聲速ae。

對高溫平衡氣體而言，邊界層參數(shù)的確定要比完全氣體復雜。波后壓力P2、焓H2和密度ρ2相互依賴，不能用波前參數(shù)簡單求出，需要用迭代方法求解[27]。

3 氣動導數(shù)提取及辨識

以求解高速飛行器俯仰方向組合動導數(shù)為例，求解組合動導數(shù)，通常對飛行器施加小幅強迫簡諧定軸振動的運動形式。強迫簡諧運動方程為

α=α0+αmsin(ωt)

(13)

根據(jù)線化氣動力理論，俯仰力矩系數(shù)的泰勒展開表達式為

(14)

(15)

(16)

式中：Lref為參考長度。

由于式(15)中包含間接量Cm0、Cmα，因此無法直接求出組合動導數(shù)。文獻[28]給出了一種求解方法，將俯仰力矩系數(shù)表達式改寫為

Cm=A+Bsin(ωt)+Ccos(ωt)

(17)

式中：A、B、C為待定系數(shù)。在俯仰振蕩非定常周期計算后，可得到俯仰力矩系數(shù)隨時間的變化曲線，在曲線上任取3組不同時刻的數(shù)據(jù)即可求出A、B、C這3個系數(shù)。結(jié)合式(15)與式(17)，從而，俯仰組合動導數(shù)可表示為

(18)

式(18)即為本文結(jié)合當?shù)亓骰钊碚撆c當?shù)乇砻嫘倍确ǎ捎脧娖日駝拥倪\動方式求解俯仰組合動導數(shù)的公式。

為比較本文方法與“采用牛頓撞擊理論直接計算動導數(shù)”的工程近似法之間的差異，下面給出工程近似法求動導數(shù)的方法[24,29]。

飛行器在體軸坐標系下的來流速度為

V∞=V∞[cosαcosβ,sinβ,cosβsinα]

(19)

物面的無量綱外法向速度投影可通過向量積求出，即

(V⊥/V∞)=nxcosαcosβ+nysinβ+

nzcosβsinα

(20)

式中：nx、ny、nz為物面法向矢量的三分量；下標⊥表示速度在垂直物面方向的分量。

在俯仰方向，由角速度引起物面與氣流相對運動的導數(shù)為

(V⊥/V∞)q=(xnz-znx)/l

(21)

式中：l表示物面與飛行器參考點之間的距離。

則根據(jù)牛頓撞擊理論，物面壓力系數(shù)的角速度導數(shù)為

Cpq=Cp,max[1/Ma∞+1.2(V⊥/V∞)]·

(V⊥/V∞)q

(22)

面元ds上力矩系數(shù)的角速度導數(shù)為

dMyq=dFaq·z-dFnq·x

(23)

式中：

dFaq=Cpq·nxds,dFnq=Cpq·nzds

(24)

(25)

式(25)即為“采用牛頓撞擊理論直接計算動導數(shù)”的工程近似法。

4 典型算例驗證

分別在超聲速和高超聲速來流條件下選取典型算例對本文方法進行驗證，以考核本文方法在不同速域范圍內(nèi)的適應性。

4.1 超聲速工況算例

在超聲速工況下，對有翼導彈(Basic Finner Missile，BFM)的俯仰組合動導數(shù)進行分析，目的是驗證本文方法超聲速工況預測動導數(shù)的精度以及效率。BFM是國際上驗證動導數(shù)計算的標準模型[30-31]，有比較成熟的實驗結(jié)果和工程估算結(jié)果。如圖3所示，BFM為尖錐形頭部、圓柱形彈身并帶有4個矩形尾翼的外形，尾翼為十字布局，圖中d為頭部直徑。

圖4為計算所用的有翼導彈面元網(wǎng)格，總網(wǎng)格量約為3.5×104。該算例的計算條件為：Ma=1.96，以頭部直徑為參考長度的雷諾數(shù)Red=0.187×106，質(zhì)心位置取外形縱向總長的61%處。

圖5為不同方法對有翼導彈的動導數(shù)預測結(jié)果對比，其中風洞試驗數(shù)據(jù)來源于文獻[30]；CFD數(shù)據(jù)為文獻[32]采用非定常CFD差分法的結(jié)果；準定常方法為對飛行器施加準定常定軸旋轉(zhuǎn)運動的工程近似結(jié)果。從圖中可以看出：在超聲速工況下本文方法能夠較好地反映有翼導彈動導數(shù)隨迎角變化規(guī)律，精度也比工程近似方法有較大的提高。

圖3 有翼導彈幾何外形示意圖

圖4 有翼導彈面元網(wǎng)格

圖5 有翼導彈動導數(shù)預測結(jié)果對比

表1為不同方法動導數(shù)計算效率對比。其中：非定常周期計算、非定常差分法預測動導數(shù)的耗時來源于文獻[32]；本文方法在Intel Core i7-8700 3.20 GHz單核上運行，僅需約5 s就能獲得所有工況的結(jié)果，表明本文方法具有極高的效率。特別的是，本文方法僅需要提供面元網(wǎng)格、不需要生成體網(wǎng)格，這節(jié)約了復雜外形體網(wǎng)格生成工作量與時間。

表1 不同方法動導數(shù)計算效率對比

Table 1 Efficiency comparison of different methods for dynamic derivatives computation

方法耗時非定常CFD周期計算[32]61 h非定常CFD差分法[32]12 h本文方法～5 s

4.2 高超聲速工況算例

在高超聲速工況下，選擇彈道外形(HyperBallistic Shape, HBS)進行俯仰組合動導數(shù)預測，目的是驗證本文方法高超聲速工況下預測動導數(shù)的精度。

HBS標模是AGARD用來測定高超聲速飛行器穩(wěn)定性的典型外形。如圖6所示，HBS由三段長均為1.5d′(d′為球體直徑)的球柱、5°半錐角的圓錐段和15°半錐角的圓錐段組成，其俯仰動導數(shù)具有較為精確的試驗結(jié)果[33]，通常被用來衡量高超聲速工況下動導數(shù)計算結(jié)果的準確程度。

圖7為彈道外形面元網(wǎng)格，總網(wǎng)格量約為2×104。 該算例的計算條件為：Ma=6.85，以頭部直徑為參考長度的雷諾數(shù)為Red′=0.72×106，質(zhì)心位置取外形縱向總長的72%處。

圖6 彈道外形示意圖

圖7 彈道外形面元網(wǎng)格

圖8為不同方法對彈道外形的動導數(shù)預測結(jié)果對比，其中風洞試驗數(shù)據(jù)來源于文獻[33]；CFD數(shù)據(jù)為文獻[34]采用非定常CFD方法的結(jié)果；準定常方法為對飛行器施加準定常定軸旋轉(zhuǎn)運動的工程近似結(jié)果。從圖中可以看出：在高超聲速工況下本文方法能夠較好地反映動導數(shù)隨迎角變化規(guī)律，而工程近似方法只能反映動導數(shù)的量級。

為了考察網(wǎng)格密度對本文方法的影響，將本算例模型分別劃分為粗糙、中等、加密3種不同密度的網(wǎng)格，以比較不同網(wǎng)格量下的動導數(shù)預測結(jié)果。不同網(wǎng)格在0°迎角下的動導數(shù)對比如表2所示。

圖8 彈道外形動導數(shù)預測結(jié)果對比

表2 不同網(wǎng)格密度對計算結(jié)果的影響

從表2可以看出，在面元網(wǎng)格能夠描述基本飛行器氣動外形的前提下，網(wǎng)格量變化對本文方法影響十分微小(不超過0.8%)。

在計算效率方面：使用非定常的雙時間方法計算本算例單個工況的俯仰動導數(shù)耗時為16 h，結(jié)合CFD和當?shù)亓骰钊碚摰姆椒ê臅r約為40 min[18]，本文方法計算單個工況的動導數(shù)平均不超過1 s，表明了本文方法具有極高的效率。

5 復雜外形飛行器氣動導數(shù)預測

為了驗證本文方法對復雜外形飛行器動導數(shù)預測的適用程度，以類X-37B飛行器為研究對象，分別采用本文方法和基于定常CFD的當?shù)亓骰钊碚摲椒A測類X-37B飛行器的動導數(shù)，并與施加強迫振動的非定常CFD方法作對比。

X-37B是美國為了驗證可重復使用空間技術和在軌空間飛行任務而啟動的項目[35]，并計劃在X-37B飛行器技術基礎上繼續(xù)進行能夠投送6名宇航員進入太空的X-37C計劃[36]。因此，以類X-37B飛行器為對象研究高速飛行器氣動導數(shù)具有典型的意義。圖9為類X-37B飛行器示意圖。

圖10為類X-37B飛行器動導數(shù)預測結(jié)果對比，其中CFD方法為文獻[37]對飛行器施加小幅強迫振蕩的非定常CFD結(jié)果；“無黏CFD+當?shù)亓骰钊骼碚摗睘榛诙ǔFD流場數(shù)據(jù)，采用當?shù)亓骰钊碚撚嬎惴嵌ǔ鈩恿?，再對動導?shù)進行提取、辨識方法得到的結(jié)果。通過圖10可以看出：本文方法很好地預測出了類X-37B飛行器的動導數(shù)隨迎角變化規(guī)律，精度與“無黏CFD+當?shù)亓骰钊骼碚摗狈椒ㄏ喈?，但本文方法效率遠高于后者。

圖9 類X-37B飛行器示意圖

圖10 類X-37B飛行器動導數(shù)預測結(jié)果對比

由于“無黏CFD+當?shù)亓骰钊骼碚摗碧崛?、辨識動導數(shù)方法與本文方法相同，因此兩者結(jié)果的誤差來源可能是無黏CFD流場與本文“當?shù)乇砻嫘倍确?激波后等熵關系”計算出的流場不一致造成。如式(4)所示：在使用當?shù)亓骰钊碚撚嬎惴嵌ǔ鈩恿r，用到了兩個重要的流場參數(shù)：邊界層外緣密度及邊界層外緣聲速。因此，本文對邊界層外緣密度及邊界層外緣聲速進行對比分析。

圖11、圖12分別為上述兩種方法得出的邊界層外緣密度與邊界層外緣聲速對比(α=30°)。從圖中可以看出，兩種方法計算出的流場大面積分布較為一致，但在駐點、翼前緣附近有較大差別。這應是兩種方法預測動導數(shù)誤差的來源。

圖11 邊界層外緣密度對比

圖12 邊界層外緣聲速對比

6 結(jié) 論

1) 本文動導數(shù)預測方法是基于“當?shù)亓骰钊碚摗奔啊芭ｎD撞擊理論”等相關理論推導而出，因此“當?shù)亓骰钊碚摗奔啊芭ｎD撞擊理論”的成立條件決定了本文方法適用的速域范圍，即擾動可近似為沿物面法向傳播的超聲速/高超聲速流動。

2) 當?shù)亓骰钊碚撏茖н^程中不需要等熵假設，拓寬了經(jīng)典活塞理論對飛行器外形和迎角的適用范圍，因此本文方法可用于復雜外形飛行器的動導數(shù)預測。

3) 本文方法避免了傳統(tǒng)動導數(shù)預測方法對CFD流場的依賴、耦合，從而大幅提高了計算效率；同時在面對復雜外形時也能較好地得出動導數(shù)隨迎角的變化規(guī)律，精度能夠滿足飛行器總體設計階段的要求，可作為飛行器布局選型階段的工具。