吸氣式內(nèi)外流一體化飛行器動(dòng)導(dǎo)數(shù)數(shù)值模擬

劉 緒，劉 偉，周云龍，柴振霞

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

吸氣式內(nèi)外流一體化飛行器動(dòng)導(dǎo)數(shù)數(shù)值模擬

劉 緒＊，劉 偉，周云龍，柴振霞

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

動(dòng)導(dǎo)數(shù)是飛行器軌道及姿態(tài)控制系統(tǒng)設(shè)計(jì)時(shí)的重要參數(shù)，對飛行器開環(huán)系統(tǒng)受到擾動(dòng)時(shí)振蕩的斂散特性起重要作用。在標(biāo)模驗(yàn)證的基礎(chǔ)上開展吸氣式內(nèi)外流一體化飛行器WR-A動(dòng)態(tài)特性分析，采用N-S方程模擬強(qiáng)迫簡諧振動(dòng)的非定常流場，獲得了飛行器的直接阻尼導(dǎo)數(shù)、加速度導(dǎo)數(shù)和旋轉(zhuǎn)導(dǎo)數(shù)，并對WR-A外形大振幅振動(dòng)時(shí)的進(jìn)氣道性能參數(shù)進(jìn)行了分析。研究表明，在某些情況下反映流動(dòng)時(shí)滯效應(yīng)的加速度導(dǎo)數(shù)占直接阻尼導(dǎo)數(shù)的比例最高可以達(dá)到40%以上，而在某些情況與直接阻尼導(dǎo)數(shù)符號相反。此外，在動(dòng)導(dǎo)數(shù)研究基礎(chǔ)上，采用參數(shù)化的動(dòng)導(dǎo)數(shù)概念初步開展了非定常氣動(dòng)參數(shù)建模研究。針對WR-A大振幅強(qiáng)迫振動(dòng)的研究顯示，采用所建立的模型預(yù)測的氣動(dòng)參數(shù)與完全非定常計(jì)算吻合較好。

吸氣式高超聲速飛行器；動(dòng)導(dǎo)數(shù)；強(qiáng)迫振動(dòng)；數(shù)值模擬；N-S方程

0 引 言

吸氣式內(nèi)外流一體化飛行器因具有強(qiáng)突防、遠(yuǎn)航程及省時(shí)等優(yōu)勢而成為當(dāng)前各國重點(diǎn)發(fā)展項(xiàng)目。相比其他傳統(tǒng)飛行器而言，它能夠提高推進(jìn)效率、增強(qiáng)載荷能力、擴(kuò)展發(fā)射窗口以及節(jié)省使用成本［1］。吸氣式內(nèi)外流一體化飛行器在助推器分離、進(jìn)氣道打開/關(guān)閉、燃燒室點(diǎn)火/熄火等過程面臨擾動(dòng)力矩帶來的飛行姿態(tài)振蕩問題［2］。進(jìn)氣道作為高超聲速飛行器的關(guān)鍵部件，飛行姿態(tài)的非定常變化必然使進(jìn)氣道內(nèi)部流動(dòng)偏離設(shè)計(jì)工況，此時(shí)進(jìn)氣道是否還能正常啟動(dòng)，流場特征及性能如何變化，給下游燃燒室提供氣流的品質(zhì)如何，均是飛行器設(shè)計(jì)中需要探索的重要問題［3］。目前對高超聲速內(nèi)外流一體化飛行器的研究大多基于穩(wěn)態(tài)，對非定常的進(jìn)氣道性能研究較為缺乏。為了實(shí)現(xiàn)動(dòng)態(tài)運(yùn)動(dòng)的飛行器氣動(dòng)性能與發(fā)動(dòng)機(jī)推進(jìn)性能之間的耦合研究，復(fù)雜內(nèi)流流動(dòng)的非定常氣動(dòng)特性是飛行器設(shè)計(jì)中必須考慮的重要因素。

動(dòng)導(dǎo)數(shù)廣泛應(yīng)用于飛行器在擾動(dòng)（包括控制面作用）下的動(dòng)態(tài)穩(wěn)定性研究之中。通過飛行器繞定軸的振動(dòng)可以獲得俯仰阻尼導(dǎo)數(shù)等直接阻尼導(dǎo)數(shù)。直接阻尼導(dǎo)數(shù)是旋轉(zhuǎn)導(dǎo)數(shù)和加速度導(dǎo)數(shù)的組合。目前動(dòng)導(dǎo)數(shù)的計(jì)算和實(shí)驗(yàn)對直接阻尼導(dǎo)數(shù)關(guān)注較多。在工程應(yīng)用中旋轉(zhuǎn)導(dǎo)數(shù)和加速度導(dǎo)數(shù)一般采用設(shè)定比例從直接阻尼導(dǎo)數(shù)中得到，但現(xiàn)代飛行器外形設(shè)計(jì)和運(yùn)動(dòng)方式比傳統(tǒng)飛行器復(fù)雜，流動(dòng)的非線性及非定常性更強(qiáng)，上述方法是否適用，值得進(jìn)一步討論。研究表明，對某些飛行器來說反映流動(dòng)時(shí)滯效應(yīng)的加速度導(dǎo)數(shù)占直接阻尼導(dǎo)數(shù)的比例可以達(dá)到40%～50%，而在某些情況加速度導(dǎo)數(shù)所占比例很小，甚至與直接阻尼導(dǎo)數(shù)符號相反［4-5］。因此，作為飛行穩(wěn)定性研究的重要參數(shù)，旋轉(zhuǎn)導(dǎo)數(shù)和加速度導(dǎo)數(shù)也是動(dòng)導(dǎo)數(shù)預(yù)測研究的重要內(nèi)容。

本文在標(biāo)模驗(yàn)證的基礎(chǔ)上開展吸氣式內(nèi)外流一體化飛行器動(dòng)態(tài)特性分析。采用三維雷諾平均的RANS方程模擬飛行器強(qiáng)迫簡諧振動(dòng)的非定常流場，開展非定常運(yùn)動(dòng)對動(dòng)態(tài)氣動(dòng)特性的影響研究。針對飛行器直接阻尼導(dǎo)數(shù)，特別是旋轉(zhuǎn)導(dǎo)數(shù)和加速度導(dǎo)數(shù)，開展多種模擬方法研究。在此基礎(chǔ)上還對基于動(dòng)導(dǎo)數(shù)的非定常參數(shù)擬合方法進(jìn)行了分析對比。

1 控制方程及離散方法

在貼體坐標(biāo)系下，完全氣體、忽略質(zhì)量力下的三維無量綱Navier-Stokes方程形式如下：

式中：Q為守恒變量，E、F、G為對流通量，Ev、Fv、Gv為粘性通量。

本文采用二階迎風(fēng)型NND格式有限體積離散流動(dòng)控制方程組的空間導(dǎo)數(shù)項(xiàng)，采用含雙時(shí)間步的LU-SGS方法來提高時(shí)間求解的效率和精度。湍流模型采用基于工程經(jīng)驗(yàn)和量綱分析的SA模型。遠(yuǎn)場邊界采用適用于動(dòng)態(tài)邊界條件下的一維Riemann不變量的無反射邊界條件。壁面邊界、速度采用無滑移條件，溫度采用絕熱壁條件，壓力條件計(jì)入離心力的影響。

2 動(dòng)導(dǎo)數(shù)數(shù)值預(yù)測方法

通過繞定軸的非定常強(qiáng)迫簡諧振動(dòng)可以獲得直接阻尼導(dǎo)數(shù)。加速度導(dǎo)數(shù)代表飛行器做沉浮運(yùn)動(dòng)或洗流（上洗或下洗）的延遲效應(yīng)及非定常渦的時(shí)間遲滯特性［4，6］，通過飛行器的沉浮或側(cè)滑運(yùn)動(dòng)可以獲得加速度導(dǎo)數(shù)。旋轉(zhuǎn)導(dǎo)數(shù)表示姿態(tài)角速度的變化引起飛行器表面與來流夾角的局部變化而產(chǎn)生的附加氣動(dòng)特性，具有準(zhǔn)定常特性［7］，既可以通過非定常強(qiáng)迫拍動(dòng)運(yùn)動(dòng)獲得，也可以通過準(zhǔn)定常方法求解。

2.1 強(qiáng)迫簡諧振動(dòng)法

圖1是飛行器三種非定常強(qiáng)迫簡諧振動(dòng)（俯仰/沉浮/拍動(dòng)）的飛行姿態(tài)和運(yùn)動(dòng)軌跡示意圖。

式中，α0為初始攻角，θm為俯仰角振幅，k為無量綱減縮頻率（，其中為圓頻率，Lref為參考長度，V∞為來流速度）。

圖1 非定常運(yùn)動(dòng)示意圖Fig.1 Sketch of unsteady motion

采用文獻(xiàn)［4］發(fā)展的強(qiáng)迫沉浮簡諧振動(dòng)（圖1（b））方法計(jì)算俯仰力矩加速度導(dǎo)數(shù)。通過給定垂直方向沉浮運(yùn)動(dòng)速度的形式，實(shí)現(xiàn)強(qiáng)迫沉浮振動(dòng)：

強(qiáng)迫拍動(dòng)簡諧振動(dòng)（圖1（c））用于計(jì)算俯仰力矩旋轉(zhuǎn)導(dǎo)數(shù)Cmq。通過繞定軸的強(qiáng)迫振動(dòng)疊加沉浮振動(dòng)可以得到強(qiáng)迫拍動(dòng)振動(dòng)方程為：

基于Etkin模型［6］，通過上式對計(jì)算產(chǎn)生的時(shí)域數(shù)據(jù)進(jìn)行后處理來辨識得到動(dòng)導(dǎo)數(shù)，具體的推導(dǎo)過程見文獻(xiàn)［8］。

2.2 準(zhǔn)定常方法

旋轉(zhuǎn)導(dǎo)數(shù)除了采用強(qiáng)迫簡諧振動(dòng)法模擬外，還可以采用準(zhǔn)定常方法得到。旋轉(zhuǎn)導(dǎo)數(shù)表示姿態(tài)角速度的變化引起飛行器表面與來流夾角的局部變化而產(chǎn)生的附加氣動(dòng)特性［9］。飛行器在攻角不變時(shí)繞z軸以恒定俯仰角速度轉(zhuǎn)動(dòng)時(shí)的平衡方程為：

式中，Cm為飛行器以勻角速率轉(zhuǎn)動(dòng)時(shí)攻角α處的俯仰力矩系數(shù)，Cm0為俯仰角速度為0時(shí)攻角α處的俯仰力矩系數(shù)，Iz為繞z軸的轉(zhuǎn)動(dòng)慣量，為角加速度（勻角速度旋轉(zhuǎn)時(shí)為0）。

Cm0可以通過定常計(jì)算獲得。在準(zhǔn)定常假設(shè)下，Cm的計(jì)算需要考慮俯仰角速度引起的牽連速度的影響以反映局部迎角的變化。在空間格式及壁面邊界條件中加入牽連速度的貢獻(xiàn)，經(jīng)時(shí)間推進(jìn)得到收斂的俯仰力矩系數(shù)Cm。由式（6）可得

準(zhǔn)定常方法能夠快速方便地求解旋轉(zhuǎn)導(dǎo)數(shù)，其計(jì)算效率與定常方法大致相等。本文針對動(dòng)導(dǎo)數(shù)標(biāo)模外形和吸氣式內(nèi)外流一體化外形，采用準(zhǔn)定常方法求解旋轉(zhuǎn)導(dǎo)數(shù)，并與完全非定常計(jì)算結(jié)果進(jìn)行了對比驗(yàn)證。

3 驗(yàn)證算例

目前，國內(nèi)外高超聲速動(dòng)導(dǎo)數(shù)標(biāo)模外形主要包括HBS（Hyper Ballistic Shape）導(dǎo)彈標(biāo)模、Finner“十”字翼標(biāo)模以及ANSR（Army-Navy spinner rocket）旋轉(zhuǎn)穩(wěn)定式火箭彈標(biāo)模。HBS和Finner標(biāo)模有組合導(dǎo)數(shù)的風(fēng)洞實(shí)驗(yàn)結(jié)果。文獻(xiàn)［8］中給出了上一節(jié)強(qiáng)迫俯仰簡諧振動(dòng)計(jì)算的組合導(dǎo)數(shù)，實(shí)驗(yàn)與計(jì)算結(jié)果取得了較好的一致。國外針對ANSR火箭彈外形開展了組合導(dǎo)數(shù)、加速度導(dǎo)數(shù)及旋轉(zhuǎn)導(dǎo)數(shù)等一系列縱向動(dòng)態(tài)穩(wěn)定性導(dǎo)數(shù)的實(shí)驗(yàn)與計(jì)算［10-13］。本節(jié)將采用ANSR標(biāo)模針對俯仰力矩導(dǎo)數(shù)及法向力導(dǎo)數(shù)開展計(jì)算驗(yàn)證。

5口徑與9口徑ANSR標(biāo)模尺寸如圖2所示，1口徑長度代表20mm。分別對兩個(gè)標(biāo)模選擇了三個(gè)不同的質(zhì)心位置（見表1），考察了質(zhì)心位置對各種阻尼導(dǎo)數(shù)的影響。

圖2 ANSR標(biāo)模尺寸（單位：口徑，1口徑＝20mm）Fig.2 ANSR dimensions（unit：caliber，1cal.＝20mm）

表1 ANSR質(zhì)心位置（單位：口徑，1口徑＝20mm）Table 1 Position of ANSR center of mass（unit：caliber，1cal.＝20mm）

三維計(jì)算網(wǎng)格采用對稱面網(wǎng)格旋轉(zhuǎn)生成。彈身附近的對稱面網(wǎng)格及拓?fù)浣Y(jié)構(gòu)如圖3所示。5口徑與9口徑ANSR標(biāo)模網(wǎng)格量分別為27.5萬和29.6萬。壁面處第一層法向網(wǎng)格尺度為10-5D，D為彈身直徑，最小網(wǎng)格雷諾數(shù)在10左右。頭部附近的網(wǎng)格拉伸比控制在1.1以下。所有的算例的馬赫數(shù)為2.5，初始攻角和側(cè)滑角均為0°，自由來流條件取標(biāo)準(zhǔn)大氣下的海平面條件（101.325kPa，288K）。強(qiáng)迫振動(dòng)的振幅為1°，無量綱減縮頻率k取0.1。

圖3 ANSR對稱面網(wǎng)格拓?fù)銯ig.3 ANSR symmetric plane mesh topology

表2給出了ANSR標(biāo)模俯仰力矩系數(shù)的旋轉(zhuǎn)導(dǎo)數(shù)。表中Cmq（1）為采用式（7）準(zhǔn)定常計(jì)算結(jié)果，Cmq（2）為采用式（5）完全非定常的強(qiáng)迫簡諧振動(dòng)法的計(jì)算結(jié)果。Cmq（3）為采用完全非定常的強(qiáng)迫簡諧振動(dòng)法計(jì)算得到的與加速度導(dǎo)數(shù)相減得到的結(jié)果。表中采用diff.表示三者中的最大值和最小值之差與三者平均值比值的絕對值，可以看出三者之間的最大差別不超過4%。說明不同方式計(jì)算得到的旋轉(zhuǎn)導(dǎo)數(shù)差別不大，分別計(jì)算的加速度導(dǎo)數(shù)和旋轉(zhuǎn)導(dǎo)數(shù)之和與直接計(jì)算得到的組合導(dǎo)數(shù)相等。由于準(zhǔn)定常方法計(jì)算效率與定常計(jì)算相同，計(jì)算開銷明顯小于非定常方法，因此下文中均采用準(zhǔn)定常方法計(jì)算得到旋轉(zhuǎn)導(dǎo)數(shù)。

表2 ANSR旋轉(zhuǎn)導(dǎo)數(shù)對比Table 2 Comparison of ANSR rotary derivatives

圖4給出了針對ANSR標(biāo)?？v向動(dòng)導(dǎo)數(shù)本文計(jì)算結(jié)果（CFD）、文獻(xiàn)［13］計(jì)算結(jié)果（PNS、SBT）以及文獻(xiàn)［11］實(shí)驗(yàn)結(jié)果（EXP）的比較。其中PNS為英國的Qin及美國軍隊(duì)研究實(shí)驗(yàn)室（ARL）的Paul Weinacht等人采用拋物化的NS方程在非慣性系中求解錐運(yùn)動(dòng)和螺旋運(yùn)動(dòng)得到的縱向動(dòng)態(tài)穩(wěn)定性導(dǎo)數(shù)［13］。這種方法好處是效率較高，可直接得到組合導(dǎo)數(shù)的分量形式（加速度導(dǎo)數(shù)和旋轉(zhuǎn)導(dǎo)數(shù)），但該方法有其局限性，只適用于軸對稱或面對稱外形的小攻角線性狀態(tài)，且僅能夠計(jì)算俯仰和偏航通道。圖中SBT（Slender Body Theory）為采用工程近似方法中細(xì)長體理論得到的結(jié)果［13］，除了圖4（a）中質(zhì)心位置取2.5口徑時(shí)各項(xiàng)俯仰力矩導(dǎo)數(shù)與CFD計(jì)算一致外，采用SBT方法的其余狀態(tài)均有量級甚至符號的差別。圖中EXP為美國陸軍彈道研究實(shí)驗(yàn)室俯仰力矩組合導(dǎo)數(shù)的實(shí)驗(yàn)結(jié)果［11］，同條件下的實(shí)驗(yàn)的數(shù)據(jù)之間具有一定的散布度。工程方法SBT在動(dòng)導(dǎo)數(shù)的預(yù)測趨勢上與實(shí)驗(yàn)符合，但數(shù)值上差別較大。PNS方法和本文CFD計(jì)算與實(shí)驗(yàn)數(shù)據(jù)基本一致，說明了本文動(dòng)導(dǎo)數(shù)計(jì)算方法的合理性。

對ANSR標(biāo)模俯仰力矩導(dǎo)數(shù)來說，反映流動(dòng)時(shí)滯效應(yīng)的加速度導(dǎo)數(shù)占組合導(dǎo)數(shù)的比例變化范圍很大，圖4（c）中時(shí)的比例能夠達(dá)到42%。隨著質(zhì)心后移加速度導(dǎo)數(shù)所占比例逐漸減小，圖4（c）中時(shí)的比例僅為3%；對ANSR標(biāo)模法向力導(dǎo)數(shù)來說，旋轉(zhuǎn)導(dǎo)數(shù)與質(zhì)心位置成近似線性關(guān)系，加速度導(dǎo)數(shù)則不受質(zhì)心位置的影響，加速度導(dǎo)數(shù)占組合導(dǎo)數(shù)的比例大于旋轉(zhuǎn)導(dǎo)數(shù)。

圖4 ANSR標(biāo)?？v向動(dòng)導(dǎo)數(shù)隨質(zhì)心變化規(guī)律Fig.4 ANSR longitudinal dynamic derivatives variation as a function of center of mass

4 WR-A外形靜、動(dòng)態(tài)氣動(dòng)特性研究

4.1 WR-A外形

吸氣式高超聲速巡航導(dǎo)彈WR-A（Wave Rider）的氣動(dòng)外形是本文仿照美國氣動(dòng)/推進(jìn)一體化設(shè)計(jì)方案X-51A乘波飛行器［14］設(shè)計(jì)而成，其外形如圖5（a）所示。表3給出了WR-A和X-51A兩種外形幾何尺寸的對比。X-51A彈道設(shè)計(jì)包括助推段、爬升段、巡航段、俯沖段，本文重點(diǎn)針對X-51A轉(zhuǎn)級點(diǎn)（助推器分離：Ma＝4.8，H＝18 288km）和設(shè)計(jì)點(diǎn)（巡航狀態(tài)：Ma＝6.5，H＝24 384km）兩個(gè)典型彈道位置開展動(dòng)態(tài)導(dǎo)數(shù)研究。

WR-A乘波飛行器采用乘波體反向設(shè)計(jì)方法［15-16］設(shè)計(jì)出的楔錐形乘波構(gòu)型［17-18］。外壓縮段采用多道斜激波壓縮加等熵波壓縮的波系配置形式，內(nèi)壓縮段設(shè)計(jì)采用混壓式矩形截面進(jìn)氣道。隔離段設(shè)計(jì)采用略有擴(kuò)張的矩形管道。燃燒室采用開式火焰穩(wěn)定凹腔，長度0.1193m，深度0.0178m，長深比6.7，后緣角30°。WR-A內(nèi)流道設(shè)計(jì)如圖5（b）所示。

圖5 WR-A外形示意圖Fig.5 WR-A layout

表3 乘波外形幾何尺寸對比Table 3 Dimension comparison of waveriders

4.2 WR-A靜態(tài)氣動(dòng)特性

針對乘波體WR-A的轉(zhuǎn)級點(diǎn)和設(shè)計(jì)點(diǎn)開展了冷流狀態(tài)下靜態(tài)氣動(dòng)特性計(jì)算，并與X-51A氣動(dòng)性能進(jìn)行比較。此外，靜態(tài)流場計(jì)算結(jié)果作為動(dòng)態(tài)非定常計(jì)算的初始流場，可以加速非定常流場的收斂，提高計(jì)算效率與精度。

4.2.1 網(wǎng)格劃分

根據(jù)吸氣式內(nèi)外流一體化外形特點(diǎn)，本文采用“O-H”型拓?fù)渖扇S結(jié)構(gòu)化網(wǎng)格（見圖6），網(wǎng)格總量為4 271萬，共劃分為208塊。內(nèi)外流一體化外形表面存在很強(qiáng)的壓縮、剪切現(xiàn)象及邊界層干擾等因素引起的較大摩擦阻力，乘波體WR-A的數(shù)值模擬顯示，該部分摩阻在總阻力中所占比重能夠達(dá)到50%以上。乘波體WR-A壁面處第一層法向網(wǎng)格尺度量級為10-7L，網(wǎng)格拉伸比控制在1.1以下。邊界層網(wǎng)格的細(xì)致刻化是壁面附近速度梯度即粘性計(jì)算準(zhǔn)確度的重要保證。

4.2.2 內(nèi)外流耦合流場特性

對于乘波體構(gòu)型機(jī)身而言，前體幾何的精確設(shè)計(jì)和激波系的結(jié)構(gòu)布置是非常重要的。本文WR-A乘波體前體壓縮面按馬赫數(shù)6.5設(shè)計(jì)點(diǎn)的斜激波交于唇口處設(shè)計(jì)，如圖7所示。來流在進(jìn)氣道唇口處減速至馬赫數(shù)4左右，下表面產(chǎn)生一道斜激波和上表面的邊界層發(fā)生干擾，產(chǎn)生邊界層分離。分離與反射激波打到進(jìn)氣道下表面后不斷向下游折射發(fā)展。圖8給出了燃燒室和噴管段流動(dòng)的發(fā)展過程。燃燒室凹腔處產(chǎn)生流動(dòng)分離，幾道斜激波之間相互交叉干擾。噴管處流動(dòng)開始膨脹加速，馬赫數(shù)最高達(dá)到6.8左右。

圖6 WR-A對稱面網(wǎng)格及拓?fù)浣Y(jié)構(gòu)Fig.6 Symmetric plane mesh and topological structure of WR-A

圖7 頭部對稱面壓力云圖Fig.7 Pressure cloud of head symmetric plane

圖8 燃燒室及尾噴管對稱面壓力云圖Fig.8 Pressure cloud of combustion chamber and nozzle symmetric plane

4.2.3 靜態(tài)氣動(dòng)力特性

轉(zhuǎn)級點(diǎn)與設(shè)計(jì)點(diǎn)的升阻比隨攻角的變化曲線如圖9所示。圖中還給出了美國空軍研究室會(huì)議報(bào)告公布的X-51A氣動(dòng)特性的數(shù)值結(jié)果［19］，做為本文計(jì)算結(jié)果的參考。WR-A轉(zhuǎn)級點(diǎn)與設(shè)計(jì)點(diǎn)的升阻比在攻角范圍-4°至4°內(nèi)升阻比近似線性，攻角超過此范圍后升阻比存在極值。通過與美國X-51A氣動(dòng)特性的比較，本文建立的WR-A模型氣動(dòng)性能符合典型高升阻比乘波體/升力體氣動(dòng)特征。

圖9 WR-A乘波體升阻比與X-51A對比Fig.9 Comparison of lift-drag ratio between WR-A waverider and X-51A

4.3 WR-A動(dòng)態(tài)氣動(dòng)特性

針對乘波體WR-A的設(shè)計(jì)點(diǎn)在0°至8°攻角范圍內(nèi)開展俯仰方向直接導(dǎo)數(shù)、加速度導(dǎo)數(shù)及旋轉(zhuǎn)導(dǎo)數(shù)數(shù)值模擬。在此基礎(chǔ)上開展基于參數(shù)化動(dòng)導(dǎo)數(shù)概念下的非定常氣動(dòng)參數(shù)建模研究。

4.3.1 設(shè)計(jì)點(diǎn)縱向動(dòng)態(tài)攻角特性分析

采用繞定軸的強(qiáng)迫俯仰簡諧振動(dòng)和沉浮運(yùn)動(dòng)分別獲得直接阻尼導(dǎo)數(shù)和加速度導(dǎo)數(shù)。強(qiáng)迫振動(dòng)的振幅αm為1°，振動(dòng)頻率f為5Hz。當(dāng)在初始瞬值衰變之后，流場參數(shù)及氣動(dòng)荷載達(dá)到諧振狀態(tài)。圖10給出了0°攻角的俯仰力矩系數(shù)Cm及進(jìn)氣道性能參數(shù)的時(shí)間歷程曲線，其中壓升比Rp、總壓恢復(fù)系數(shù)σ和流量系數(shù)φ均為隔離段入口截面上的質(zhì)量加權(quán)平均值。由于非定常計(jì)算采用靜態(tài)流場作為初場，收斂速度很快，在不到半個(gè)周期內(nèi)初始瞬值衰變完成，流場達(dá)到諧振狀態(tài)。

飛行器的非定常運(yùn)動(dòng)使得進(jìn)氣道內(nèi)部流場處于動(dòng)態(tài)變化的狀態(tài)。此時(shí)針對內(nèi)流流場特征及性能變化的預(yù)測非常重要。例如壓升比Rp反映了進(jìn)氣道抗反壓的能力，如果動(dòng)態(tài)非定常運(yùn)動(dòng)過程中Rp下降，則有可能導(dǎo)致燃燒室反壓對進(jìn)氣道產(chǎn)生干擾而進(jìn)入不啟動(dòng)狀態(tài)。2011年X-51高超聲速飛行器第二次飛行實(shí)驗(yàn)中就是因?yàn)槌紱_壓發(fā)動(dòng)機(jī)的進(jìn)氣道未能啟動(dòng)而導(dǎo)致失?。?0］。目前對進(jìn)氣道性能的計(jì)算主要還是針對穩(wěn)態(tài)，本文則根據(jù)Etkin理論，采用與力矩導(dǎo)數(shù)相同的方式對進(jìn)氣道性能參數(shù)導(dǎo)數(shù)進(jìn)行求解，可以用于動(dòng)態(tài)非定常時(shí)進(jìn)氣道性能的預(yù)測，為飛行器進(jìn)排氣系統(tǒng)與控制系統(tǒng)耦合問題的解決提供相應(yīng)的氣動(dòng)參數(shù)。

圖10 俯仰力矩系數(shù)及進(jìn)氣道性能參數(shù)的時(shí)間歷程曲線Fig.10 Time history of pitch moment coefficient and air inlet performance parameters

圖11給出了俯仰力矩系數(shù)及進(jìn)氣道性能參數(shù)的動(dòng)導(dǎo)數(shù)預(yù)測結(jié)果。通過繞定軸的強(qiáng)迫俯仰振動(dòng)，可以得到不同氣動(dòng)參數(shù)λ的俯仰導(dǎo)數(shù)。通過沉浮運(yùn)動(dòng)，得到加速度導(dǎo)數(shù)。通過準(zhǔn)定常計(jì)算，得到旋轉(zhuǎn)導(dǎo)數(shù)λq。將準(zhǔn)定常計(jì)算出的旋轉(zhuǎn)導(dǎo)數(shù)與加速度導(dǎo)數(shù)相加得到。其中λ代表俯仰力矩系數(shù)Cm、壓升比Rp、總壓恢復(fù)系數(shù)σ和流量系數(shù)φ。圖11中不同攻角下基本相等。對圖11（a）中俯仰力矩系數(shù)來說，在攻角4°和6°時(shí)加速度導(dǎo)數(shù)與組合導(dǎo)數(shù)符號相反，起不穩(wěn)定的作用。圖11（a）中旋轉(zhuǎn)導(dǎo)數(shù)曲線與組合導(dǎo)數(shù)曲線十分靠近，加速度導(dǎo)數(shù)占組合導(dǎo)數(shù)的比例較低，最大不超過20%，反映出WR-A飛行器做沉浮運(yùn)動(dòng)的延遲效應(yīng)及非定常渦的時(shí)間遲滯特性較弱。進(jìn)氣道性能參數(shù)的加速度導(dǎo)數(shù)則與之相反，如圖11（b）～圖11（d）所示，加速度導(dǎo)數(shù)曲線更加接近組合導(dǎo)數(shù)曲線，占組合導(dǎo)數(shù)的比例普遍高于旋轉(zhuǎn)導(dǎo)數(shù)，說明非定常效應(yīng)對進(jìn)氣道性能參數(shù)變化的影響主要來自于非定常流場的時(shí)滯效應(yīng)。

不同攻角下繞定軸的強(qiáng)迫俯仰振動(dòng)和沉浮運(yùn)動(dòng)俯仰力矩系數(shù)的遲滯環(huán)曲線如圖12所示。遲滯環(huán)的傾斜形狀、面積及旋轉(zhuǎn)方向決定了動(dòng)導(dǎo)數(shù)的大小與符號。由于WR-A乘波體定常狀態(tài)的俯仰力矩系數(shù)具有較強(qiáng)的非線性特征，因此圖12中動(dòng)態(tài)條件下的遲滯環(huán)并不全是傾斜光滑的橢圓形狀。遲滯環(huán)面積表示乘波體振動(dòng)過程中環(huán)境對其做功的大小，WR-A乘波體沉浮運(yùn)動(dòng)的遲滯環(huán)面積明顯小于俯仰振動(dòng)，這與圖11（b）中阻尼大小的規(guī)律相一致。圖12中還用箭頭標(biāo)出了沉浮運(yùn)動(dòng)遲滯環(huán)的旋轉(zhuǎn)方向。遲滯環(huán)的旋轉(zhuǎn)方向表示做功的方向，逆時(shí)針旋轉(zhuǎn)表示乘波體對環(huán)境做功，消耗乘波體動(dòng)能，對動(dòng)導(dǎo)數(shù)的貢獻(xiàn)為負(fù)值，起動(dòng)穩(wěn)定作用。順時(shí)針表示環(huán)境對乘波體做功，乘波體動(dòng)能增加，對動(dòng)導(dǎo)數(shù)的貢獻(xiàn)為正值，起不穩(wěn)定作用。WR-A飛行器沉浮運(yùn)動(dòng)的遲滯環(huán)在攻角2°和4°時(shí)變換方向，在6°和8°時(shí)呈8字環(huán)運(yùn)動(dòng)，說明飛行器的加速度導(dǎo)數(shù)經(jīng)歷了三次變號。從遲滯環(huán)圖中可以看出WR-A飛行器的加速度導(dǎo)數(shù)在攻角3°至6.2°和大于7.8°攻角的范圍起不穩(wěn)定作用。

圖11 俯仰力矩系數(shù)及進(jìn)氣道性能參數(shù)的動(dòng)導(dǎo)數(shù)預(yù)測結(jié)果Fig.11 Dynamic derivative predictions of pitch moment coefficient and air inlet performance parameters

圖12 強(qiáng)迫俯仰振動(dòng)和沉浮運(yùn)動(dòng)俯仰力矩系數(shù)遲滯環(huán)Fig.12 Pitch moment coefficient hysteresis loop of forced pitch vibration and plunge motion

圖13給出了總壓恢復(fù)系數(shù)在不同攻角下繞定軸的強(qiáng)迫俯仰振動(dòng)和沉浮運(yùn)動(dòng)的遲滯環(huán)曲線。隨著攻角的增加，進(jìn)氣道前緣壓縮激波的強(qiáng)度增加，外壓段壓縮能力增強(qiáng)，外壓波系后的馬赫數(shù)降低。當(dāng)攻角從負(fù)攻角上升至1.5°時(shí)，內(nèi)流道激波邊界層干擾效應(yīng)降低，進(jìn)氣道出口總壓恢復(fù)性能增加。當(dāng)攻角從1.5°繼續(xù)上升時(shí)，多波系壓縮的外壓段產(chǎn)生激波干擾作用，進(jìn)氣道邊界層增厚，從而在一定程度上降低進(jìn)氣道流量捕獲能力并使出口總壓恢復(fù)性能降低。

4.3.2 基于參數(shù)化動(dòng)導(dǎo)數(shù)概念的非定常氣動(dòng)參數(shù)建模研究

圖13 強(qiáng)迫俯仰振動(dòng)和沉浮運(yùn)動(dòng)總壓恢復(fù)系數(shù)遲滯環(huán)Fig.13 Total pressure recovery coefficient hysteresis loop of force pitch vibration and plunge motion

非定常氣動(dòng)力的建模方法一直是空氣動(dòng)力學(xué)的重要研究內(nèi)容。當(dāng)飛行器在小攻角條件下飛行時(shí)，流動(dòng)以附著流型為主，氣動(dòng)力與運(yùn)動(dòng)狀態(tài)參數(shù)常常成線性關(guān)系，因此動(dòng)導(dǎo)數(shù)可看成是常數(shù)，可簡單得到非定常氣動(dòng)力的擬合結(jié)果?，F(xiàn)代飛行器外形設(shè)計(jì)和運(yùn)動(dòng)方式比傳統(tǒng)飛行器復(fù)雜，大攻角和非零側(cè)滑角的飛行狀態(tài)通常會(huì)引起非常復(fù)雜的氣動(dòng)現(xiàn)象，激波誘導(dǎo)分離、旋渦運(yùn)動(dòng)與破裂以及它們的相互作用使得流體運(yùn)動(dòng)呈現(xiàn)強(qiáng)烈的非定常、非線性特征。動(dòng)導(dǎo)數(shù)不再是常數(shù)，而是和攻角、側(cè)滑角、減縮頻率等參數(shù)相關(guān)，具有參數(shù)化的性質(zhì)。為了更好地描述動(dòng)態(tài)特性在大攻角時(shí)出現(xiàn)的非線性特征，本節(jié)采用參數(shù)化的動(dòng)導(dǎo)數(shù)概念建立非定常氣動(dòng)參數(shù)模型。通過對比氣動(dòng)參數(shù)模型與采用完全非定常計(jì)算得到的大振幅強(qiáng)迫振動(dòng)的非線性氣動(dòng)力結(jié)果，開展WR-A飛行器非定常氣動(dòng)參數(shù)擬合方法研究。

WR-A乘波體在設(shè)計(jì)點(diǎn)位置大振幅強(qiáng)迫運(yùn)動(dòng)的基準(zhǔn)攻角為0°，振幅為1°至8°，頻率為5Hz。對于作俯仰運(yùn)動(dòng)的飛行器，決定繞俯仰軸轉(zhuǎn)動(dòng)運(yùn)動(dòng)的狀態(tài)變量為。動(dòng)導(dǎo)數(shù)是狀態(tài)變量α的函數(shù)。本文首先在前述工作基礎(chǔ)上計(jì)算出i（i＝0，±1，±2，…，±8）攻角下的靜態(tài)與動(dòng)態(tài)氣動(dòng)參數(shù)，其中λ代表氣動(dòng)力系數(shù)或進(jìn)氣道性能參數(shù)。當(dāng)攻角α處于［i-0.5，i＋0.5］區(qū)間內(nèi)時(shí)，λ的表達(dá)式為：

式（8）考慮了動(dòng)導(dǎo)數(shù)的參數(shù)化特征，是一種局部擬合方式。如果考慮更加復(fù)雜的運(yùn)動(dòng)，還可以保留動(dòng)導(dǎo)數(shù)的高階項(xiàng)、交叉導(dǎo)數(shù)和交叉耦合導(dǎo)數(shù)項(xiàng)等。

圖14中的曲線為采用完全非定常計(jì)算得到的俯仰力矩系數(shù)和進(jìn)氣道性能參數(shù)的遲滯環(huán)。圖14中的數(shù)據(jù)點(diǎn)表示采用參數(shù)化的動(dòng)導(dǎo)數(shù)概念建立的氣動(dòng)模型計(jì)算得到的俯仰力矩系數(shù)Cms、壓升比Rps、總壓恢復(fù)系數(shù)σs和流量系數(shù)φs。當(dāng)攻角小于0°時(shí)，采用式（8）計(jì)算得到的Cms全部準(zhǔn)確地落在不同振幅下的俯仰力矩遲滯曲線上。當(dāng)攻角大于0°時(shí)，俯仰力矩系數(shù)曲線開始出現(xiàn)較為明顯的非線性特征，圖14（a）中4°振幅以內(nèi)的Cms與遲滯曲線符合較好。振幅6°和8°的Cms與遲滯曲線有一定差別，但趨勢上與遲滯曲線相同。說明對俯仰力矩系數(shù)而言，當(dāng)運(yùn)動(dòng)處于非線性特性較弱的區(qū)域時(shí)（攻角-8°至-2°），參數(shù)化的動(dòng)導(dǎo)數(shù)概念建立的氣動(dòng)模型預(yù)測大攻角的氣動(dòng)數(shù)據(jù)與非定常計(jì)算符合較好。當(dāng)運(yùn)動(dòng)處于非線性特性較強(qiáng)的區(qū)域時(shí)（0°攻角以上），大攻角的氣動(dòng)數(shù)據(jù)與非定常計(jì)算有一定偏差，總體吻合較好。

圖14 大攻角振動(dòng)遲滯環(huán)曲線Fig.14 Hysteresis loop curve of large-amplitude vibration

圖14（b）～圖14（d）給出了動(dòng)導(dǎo)數(shù)模型預(yù)測進(jìn)氣道性能參數(shù)時(shí)的特性曲線。與俯仰力矩系數(shù)不同的是，無論運(yùn)動(dòng)處于正攻角或負(fù)攻角范圍，采用參數(shù)化的動(dòng)導(dǎo)數(shù)概念建立的氣動(dòng)模型預(yù)測的進(jìn)氣道性能參數(shù)與非定常計(jì)算一致。圖14（b）和圖14（d）中的流量系數(shù)和壓升比隨著攻角的上升呈線性增加的趨勢。圖14（c）中總壓恢復(fù)系數(shù)則在1.5°攻角處增大到極值位置，然后逐漸減小。在負(fù)攻角范圍內(nèi)，總壓恢復(fù)系數(shù)的滯后特性微弱，遲滯環(huán)的上行程和下行程幾乎重合。當(dāng)飛行器攻角大于0°時(shí)，總壓恢復(fù)系數(shù)的滯后效應(yīng)開始變得明顯，遲滯環(huán)曲線在攻角越大的位置越顯得“豐滿”，總壓恢復(fù)系數(shù)導(dǎo)數(shù)的絕對值也逐漸增大，與圖11（c）所反映的規(guī)律一致。

5 結(jié) 論

本文在標(biāo)模驗(yàn)證的基礎(chǔ)上開展吸氣式內(nèi)外流一體化飛行器動(dòng)態(tài)特性分析。采用三維雷諾平均的RANS方程模擬飛行器強(qiáng)迫簡諧振動(dòng)的非定常流場，通過建立不同的強(qiáng)迫簡諧振動(dòng)模式，獲得了直接阻尼導(dǎo)數(shù)、加速度導(dǎo)數(shù)和旋轉(zhuǎn)導(dǎo)數(shù)的辨識方法。并在此基礎(chǔ)上得到的主要結(jié)論如下：

對俯仰力矩系數(shù)的加速度導(dǎo)數(shù)來說，WR-A乘波體的加速度導(dǎo)數(shù)占組合導(dǎo)數(shù)的比例較低，最大不超過20%，最小時(shí)與組合導(dǎo)數(shù)符號相反。而ANSR火箭彈的加速度導(dǎo)數(shù)所占比例變化范圍很大，最高能夠達(dá)到40%以上，最低僅為3%。

基于Etkin理論對進(jìn)氣道性能參數(shù)動(dòng)態(tài)導(dǎo)數(shù)進(jìn)行求解，用于分析非定常內(nèi)流流場特性。WR-A乘波體進(jìn)氣道性能參數(shù)的加速度導(dǎo)數(shù)占組合導(dǎo)數(shù)的比例大于旋轉(zhuǎn)導(dǎo)數(shù)，說明非定常效應(yīng)對進(jìn)氣道性能參數(shù)變化的影響主要來自于非定常流場的時(shí)滯效應(yīng)。

初步開展了參數(shù)化動(dòng)導(dǎo)數(shù)概念下的非定常氣動(dòng)參數(shù)建模研究。采用所建立的模型預(yù)測WR-A乘波體在大攻角時(shí)的俯仰力矩系數(shù)，在非線性特性較弱的區(qū)域與非定常計(jì)算符合較好，在非線性特性較強(qiáng)的區(qū)域則產(chǎn)生一定偏差，非定常氣動(dòng)參數(shù)模型預(yù)測的大攻角條件下的進(jìn)氣道性能參數(shù)與非定常計(jì)算吻合較好。

［1］ Bilaedo V J，Curran F M，Hunt J L，et al.The benefits of hypersonic airbreathing launch systems for access to space［R］.AIAA 2003-5265.

［2］ Mirmirani M D，Wu C，Clark A，et al.Modeling for control of ageneric airbreathing hypersonic vehicle［C］//Proceedings of AIAA Conference on Guidance，Navigation，and Control Conference，2005.

［3］ Bahm C，Baumann E，Martin J，et al.The X-43AHyper-X Mach 7 flight 2guidance，navigation，and control overview and flight test results［R］.AIAA 2005-3275.

［4］ Liu W，Yang X L，Zhao Y F.Numerical simulation of acceleration derivative of hypersonic aircraft［J］.Acta Aerodynamica Sinica，2010，（04）：426-429.（in Chinese）劉偉，楊小亮，趙云飛.高超聲速飛行器加速度導(dǎo)數(shù)數(shù)值模擬［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2010，（04）：426-429.

［5］ Sun H S.A measurement technique for derivatives of aircraft due to acceleration in heave and sideslip at high angle of attack［J］.Experiments and Measurements in Fluid Mechanics，2001，15（4）：15-19.（in Chinese）孫海生.飛行器大攻角升沉平移加速度導(dǎo)數(shù)測量技術(shù)［J］.流體力學(xué)實(shí)驗(yàn)與測量，2001，15（4）：15-19.

［6］ Etkin B.Dynamics of atmospheric flight［M］.Courier Dover Publications，2012.

［7］ Bergmann A，Hübner A.Integrated experimental and numerical research on the aerodynamics of unsteady moving aircraft：proceedings of the 3rd international symposium on integrating CFD and experiments in aerodynamics［R］.U.S.Air Force Academy，CO，USA，2007.

［8］ Liu X.Investigation of dynamic characteristics of hypersonic airframe/propulsion integrative vehicle［D］.Changsha：National University of Defense Technology，2011.劉緒.高超聲速內(nèi)外流一體化飛行器動(dòng)態(tài)特性研究［D］.長沙：國防科學(xué)技術(shù)大學(xué)，2011.（in Chinese）

［9］ Xi K，Yan C，Huang Y，et al.Numerical simulation of individual components of pitch-damping coefficient sum［J］.Journal of Beijing University of Aeronautics and Astronautics，2015，41（2）：222-227.（in Chinese）席柯，閻超，黃宇，等.俯仰阻尼導(dǎo)數(shù)分量的CFD數(shù)值模擬［J］.北京航空航天大學(xué)學(xué)報(bào)，2015，41（2）：222-227.

［10］DeSpirito J，Silton S I，Weinacht P.Navier-Stokes predictions of dynamic stability derivatives：evaluation of steady-state methods［J］.Journal of Spacecraft and Rockets，2009，46（6）：1142-1154.

［11］Murphy C H，Schmidt L E.The effect of length on the aerodynamics characteristics of bodies of revolution in supersonic flight［R］.U.S.Army Ballistic Research Lab.，Rept.876，1953.

［12］Schmidt L，Murphy C.The aerodynamic properties of the 7-caliber army-navy spinner rocket in transonic flight［R］.U.S.Army Ballistic Research Lab.，BRL-MR-775，1954.

［13］Weinacht P，Sturek W B，Schiff L B.Navier-Stokes predictions of pitch damping for axisymmetric projectiles［J］.Journal of Spacecraft and Rockets，1997，34（6）：753-761.

［14］Hank J M，Murphy J S，Mutzman R C.The X-51Ascramjet engine flight demonstration program［R］.AIAA 2008-2540.

［15］Kuchemann D.The aerodynamic design of aircraft［M］.AIAA Education Series，Amer Inst of Aeronautics and Astronautics，2012.

［16］Bauer S.Analysis of two viscous optimized waveriders［C］//Proceeding of First International Hypersonic Waverider Symposium，1990.

［17］Takashima N，Lewis M J.Waverider configurations based on non-axisymmetric flow fields for engine-airframe integration［C］//Proceedings of 32nd Aerospace Sciences Meeting and Exhibit.Reno，NV，1994.

［18］Wang Zhuo.High speed flow design using the theory of osculating cones and axisymmetric flows［J］.Chinese Journal of Aeronautics，1999，（01）：3-10.

［19］Mutzman R，Murphy S.X-51development：a chief engineer’s perspective［C］//Proceedings of the 17th AIAA International Space Planes and Hypersonic Systems and Technologies Conference.2011.

［20］Bai Y L，Bai Y.Failure analysis of X-51Aaircraft flight test［J］.Cruise Missile，2012，（3）：27-31.（in Chinese）白延隆，白云.X-51A飛行器飛行試驗(yàn)的故障分析［J］.飛航導(dǎo)彈，2012，（3）：27-31.

Numerical simulation of dynamic derivatives for air-breathing hypersonic vehicle

Liu Xu，Liu Wei，Zhou Yunlong，Chai Zhenxia
（CollegeofAerospaceScienceandEngineering，NationalUniversityofDefenseTechnology，Changsha410073，China）

Dynamic derivatives are important parameters for designing aircraft trajectory and attitude control system，and for deciding the divergence behavior of vibration under disturbance.After calibration model validation，the dynamic behavior of an air-breathing hypersonic vehicle，namely WR-A（Wave Rider）is characterized.The unsteady flow field around the under aircraft forced simple harmonic vibration（SHV）condition is simulated using Navier Stokes equation.The direct damping derivatives，acceleration derivatives and rotary derivatives of this air-breathing hypersonic vehicle are obtained.The air inlet performance parameter derivatives are solved using Etkin predictive aerodynamic model.The air inlet performance parameters under large-amplitude vibration are successfully predicted using the dynamic derivative model.This offers a guideline for characterizing the dynamic unsteady internal flow field and predicting air inlet performance variation.The proportion of acceleration derivative，which represents the flow time lag effect，in the direct damping derivative can be as high as forty percent but is opposite to the damping derivative direction in some cases，contributing to dynamic instability adversely.It is reasonable to using this dynamic derivative model to express the aerodynamic behavior of air-breathing hypersonic vehicle at large angles of attack according to this WR-A large-amplitude vibration simulation.

air-breathing hypersonic vehicle；dynamic derivative；forced vibration；numerical simulation；N-S equation

V211.3

：Adoi：10.7638/kqdlxxb-2015.0020

0258-1825（2015）02-0147-09

2015-01-02；

：2015-03-13

國家自然科學(xué)基金（90716015，11172325）

劉緒＊（1986-），男，博士研究生，研究方向：計(jì)算流體力學(xué)與應(yīng)用.E-mail：liuxuqd＠126.com

劉緒，劉偉，周云龍，等.吸氣式內(nèi)外流一體化飛行器動(dòng)導(dǎo)數(shù)數(shù)值模擬［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2015，33（2）：147-155.

10.7638/kqdlxxb-2015.0020 Liu X，Liu W，Zhou Y L，et al.Numerical simulation of dynamic derivatives for air-breathing hypersonic vehicle［J］.Acta Aerodynamica Sinica，2015，33（2）：147-155.