畢 林,袁先旭,,吳波佼,陳 浩,唐志共

(1.空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,綿陽(yáng) 621000;2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所,綿陽(yáng) 621000)

0 引 言

近年來(lái),國(guó)際火星探測(cè)迎來(lái)了新的高潮。由于進(jìn)入火星大氣時(shí)減速的需要,火星探測(cè)器為球冠倒錐構(gòu)型的鈍體飛行器[1-3]。對(duì)于此類飛行器,其后體流動(dòng)較為復(fù)雜,存在大范圍的分離與再附、激波/旋渦干擾等,相應(yīng)地,其動(dòng)態(tài)特性也較為復(fù)雜,是設(shè)計(jì)關(guān)鍵問(wèn)題之一。對(duì)于再入地球大氣的球冠倒錐鈍體飛船返回艙,其動(dòng)態(tài)特性問(wèn)題已有較多研究。但由于火星大氣與地球大氣差異較大,平均大氣密度約為地球的1.6%,表面大氣壓力僅為700 Pa,相當(dāng)于地球30~40 km高度時(shí)大氣壓力[4-5]。這么稀薄的大氣環(huán)境下,火星探測(cè)器進(jìn)入火星大氣環(huán)境的動(dòng)態(tài)特性問(wèn)題,已成為火星探測(cè)項(xiàng)目發(fā)展面臨的挑戰(zhàn)之一。

國(guó)內(nèi)外火星環(huán)境的動(dòng)態(tài)特性研究方法與地球環(huán)境的動(dòng)態(tài)特性研究方法類似,主要有半經(jīng)驗(yàn)理論方法、數(shù)值模擬方法和火星風(fēng)洞實(shí)驗(yàn)方法。

由于探測(cè)器在火星環(huán)境中的流場(chǎng)十分復(fù)雜,理論分析方法往往難以奏效。目前對(duì)返回艙的動(dòng)穩(wěn)定性研究主要依靠地面模擬試驗(yàn),包括強(qiáng)迫振蕩法、自由振動(dòng)法、自由翻滾法和模型自由飛試驗(yàn)法等[6-7]。然而,由于地面風(fēng)洞設(shè)備試驗(yàn)氣體以模擬空氣為主,主要針對(duì)地球環(huán)境的返回艙動(dòng)穩(wěn)定性問(wèn)題。目前關(guān)于火星環(huán)境探測(cè)器動(dòng)態(tài)特性預(yù)測(cè)的公開數(shù)據(jù)相對(duì)匱乏,其中,又以彈道靶中進(jìn)行的模型自由飛試驗(yàn)數(shù)據(jù)為主。如美國(guó)空軍ARF在常溫、常壓的空氣環(huán)境下進(jìn)行了MER火星探測(cè)器模型的彈道靶試驗(yàn)[8]。試驗(yàn)表明,在迎角較小的狀態(tài)下,這種短鈍頭體外形是不穩(wěn)定的,并且隨著速度的降低,穩(wěn)定性進(jìn)一步減弱;隨著質(zhì)心的前移,模型穩(wěn)定性增加。然而,常規(guī)空氣介質(zhì)風(fēng)洞所得的數(shù)據(jù)需要進(jìn)行地火相關(guān)性修正,才能模擬火星大氣相關(guān)氣動(dòng)問(wèn)題。美國(guó)阿姆斯研究中心的HFFAF(The Hypervelocity Free Flight Aerodynamics Facility)彈道靶中,進(jìn)行了二氧化碳環(huán)境下MSL(The Mars Smart Lander)探測(cè)器模型自由飛試驗(yàn)[9],獲得了靜、動(dòng)導(dǎo)數(shù)等動(dòng)態(tài)數(shù)據(jù),但公開的試驗(yàn)數(shù)據(jù)很少。

火星探測(cè)的飛行數(shù)據(jù)稀缺,地面試驗(yàn)難度大,工程計(jì)算方法進(jìn)展緩慢,目前,數(shù)值模擬可作為火星探測(cè)氣動(dòng)問(wèn)題研究的重要工具。其中,近似比熱比數(shù)值模擬方法很大程度地簡(jiǎn)化了計(jì)算模型[10-11],但國(guó)內(nèi)外在近似比熱比的取值上有所不同,且研究尚顯不足。

因此,著眼于我國(guó)火星探測(cè)項(xiàng)目需求、亟需解決探測(cè)器進(jìn)入火星大氣環(huán)境中所遇到的氣動(dòng)問(wèn)題,開展火星探測(cè)器氣動(dòng)穩(wěn)定性的研究。本文采用近似比熱比數(shù)值模擬方法,基于軟件in-house programm FLY3DMars,針對(duì)火星探測(cè)器彈道靶自由飛試驗(yàn)外形,進(jìn)行了火星大氣模擬環(huán)境下俯仰氣動(dòng)特性的數(shù)值模擬研究。

1 數(shù)值模擬方法

1.1 控制方程

直角坐標(biāo)系中,不計(jì)體積力和外加熱量的無(wú)量綱化、守恒形式Navier-Stokes方程組為:

其中,Q為守恒變量,E、F和G為對(duì)流項(xiàng)、無(wú)黏通矢量,Ev、Fv和Gv為黏性通矢量,具體形式分別為:

通過(guò)坐標(biāo)變換,可將方程組(1)轉(zhuǎn)化為一般曲線坐標(biāo)系中的形式:

將一般曲線坐標(biāo)系中的Navier-Stokes方程組(2)寫作:

1.2 動(dòng)網(wǎng)格技術(shù)

動(dòng)態(tài)非定常數(shù)值模擬與定態(tài)數(shù)值模擬的區(qū)別之一是需要建立相對(duì)于地面系的動(dòng)坐標(biāo)系,以描述飛行器的姿態(tài)和運(yùn)動(dòng)。本文使用的地面系與飛行力學(xué)習(xí)慣建立的地面系之間的差別在于:y軸方向不變,x、z軸方向相反。同樣,本文使用的體軸系的差別在于:yb軸方向不變,xb、zb軸方向相反。

一般俯仰運(yùn)動(dòng)計(jì)算時(shí),可使用半流場(chǎng)網(wǎng)格;偏航、滾轉(zhuǎn)運(yùn)動(dòng)計(jì)算,使用全流場(chǎng)網(wǎng)格。對(duì)于網(wǎng)格運(yùn)動(dòng)要求滿足幾何守恒律,即:

由于計(jì)算非定常動(dòng)態(tài)流場(chǎng)需要每一時(shí)間步都更新網(wǎng)格,因此動(dòng)網(wǎng)格生成技術(shù)的計(jì)算量小就顯得尤為重要。目前,剛性運(yùn)動(dòng)網(wǎng)格生成和超限插值生成動(dòng)網(wǎng)格是兩種常用的方法。剛性運(yùn)動(dòng)網(wǎng)格生成技術(shù)就是計(jì)算網(wǎng)格隨物體一起作剛體運(yùn)動(dòng),其計(jì)算量小,并可保持初始網(wǎng)格的品質(zhì),但不適于生成較復(fù)雜的網(wǎng)格,也使運(yùn)動(dòng)開邊界的處理不易。超限插值生成動(dòng)網(wǎng)格方法就是外邊界保持靜止,物面邊界由物體運(yùn)動(dòng)規(guī)律或運(yùn)動(dòng)方程得到,內(nèi)場(chǎng)網(wǎng)格由超限插值的方法代數(shù)生成,其計(jì)算量也較小,能夠生成復(fù)雜網(wǎng)格,但不易保證網(wǎng)格品質(zhì),尤其是不能保證物面網(wǎng)格的正交性?？捎眉訖?quán)技術(shù)來(lái)綜合這兩種方法的優(yōu)點(diǎn)[3]。

首先,由初始網(wǎng)格Xn生成剛性運(yùn)動(dòng)網(wǎng)格Xref1,再由n時(shí)刻的開邊界和n＋1時(shí)刻的物面邊界代數(shù)插值得到n＋1時(shí)刻的網(wǎng)格Xref2,最后由兩者加權(quán)得到新時(shí)刻的動(dòng)網(wǎng)格:

只要適當(dāng)選擇加權(quán)因子ω,就既可保證物面附近網(wǎng)格質(zhì)量,又可使外邊界保持靜止而易于處理?？蓸?gòu)造如下加權(quán)因子,計(jì)算表明效果較好,令j為法向網(wǎng)格指標(biāo),則有:

1.3 靜、動(dòng)態(tài)氣動(dòng)穩(wěn)定性參數(shù)計(jì)算辨識(shí)方法

在預(yù)測(cè)靜、動(dòng)態(tài)氣動(dòng)穩(wěn)定性參數(shù)方面,數(shù)值計(jì)算明顯優(yōu)于地面模擬試驗(yàn)方法:數(shù)值計(jì)算費(fèi)用低,可提供較長(zhǎng)的非定常動(dòng)態(tài)氣動(dòng)力周期,振心位置可與實(shí)際重心位置一致[6]。靜、動(dòng)導(dǎo)數(shù)參數(shù)辨識(shí)的思路在火星環(huán)境中仍然適用,主要是將動(dòng)態(tài)氣動(dòng)力進(jìn)行Taylor級(jí)數(shù)展開,保留線性部分,以俯仰運(yùn)動(dòng)為例,對(duì)于單自由度俯仰運(yùn)動(dòng),根據(jù)Etkin經(jīng)典動(dòng)態(tài)氣動(dòng)力模型,有

其中,θ＝α－α0,q下標(biāo)“0”表示在迎角α0處的取值,α0一般為平衡迎角αt,此時(shí)線性假設(shè)成立。工程上稱Cmα為靜導(dǎo)數(shù)＋Cmq為俯仰阻尼導(dǎo)數(shù)[3]。

本文采用強(qiáng)迫振動(dòng)法求取動(dòng)導(dǎo)數(shù)。以單自由度俯仰強(qiáng)迫振動(dòng)為例,給定簡(jiǎn)諧振蕩(無(wú)量綱化)形式如下:

其中,α0為振蕩的中心迎角,Am為振蕩的幅值,ω為振蕩的圓頻率,k＝ωLref/2V∞為減縮頻率(無(wú)量綱量),θ為俯仰振蕩角。對(duì)動(dòng)態(tài)俯仰力矩系數(shù)在振蕩中心迎角處展開如下:

其中,T＝π/k為強(qiáng)迫振蕩的周期。強(qiáng)迫振蕩的振幅對(duì)辨識(shí)結(jié)果的影響不大,一般取為1°。由于飛行器在空中的真實(shí)頻率并不知道,需要考察多個(gè)頻率對(duì)動(dòng)穩(wěn)定性方面的影響,才能得到有效結(jié)果。單自由度強(qiáng)迫偏航/滾轉(zhuǎn)及其靜、動(dòng)導(dǎo)數(shù)的辨識(shí)類似。

2 算例驗(yàn)證

第2節(jié)針對(duì)非軸對(duì)稱M3外形,開展俯仰強(qiáng)迫振動(dòng)數(shù)值模擬,通過(guò)將數(shù)值模擬結(jié)果與彈道靶試驗(yàn)結(jié)果對(duì)比,考核in-house programm FLY3D-MARS軟件在火星環(huán)境動(dòng)態(tài)計(jì)算中的適用性。

2.1 MSL彈道靶試驗(yàn)動(dòng)態(tài)結(jié)果

在美國(guó)NASA阿姆斯研究中心,針對(duì)MSL彈道靶模型進(jìn)行的彈道靶模型自由飛試驗(yàn)中,采用44.5 mm口徑的噴氣槍[9],在純CO2、CO2混合少量空氣以及空氣介質(zhì)中進(jìn)行對(duì)比試驗(yàn),得到了探測(cè)器M3模型的動(dòng)態(tài)數(shù)據(jù)。彈道靶試驗(yàn)動(dòng)態(tài)結(jié)果并未全部公開,并且缺乏動(dòng)態(tài)試驗(yàn)中的減縮頻率等條件,相關(guān)動(dòng)態(tài)結(jié)果十分匱乏,文獻(xiàn)僅僅給出了M3外形在不同介質(zhì)中的少量動(dòng)導(dǎo)數(shù)結(jié)果,如圖1所示,圖中為彈道靶試驗(yàn)中,三種不同的氣體介質(zhì)試驗(yàn)中,M3模型動(dòng)導(dǎo)數(shù)結(jié)果。由圖可知,在考察的彈道靶試驗(yàn)狀態(tài)中,空氣介質(zhì)以及空氣與CO2混合介質(zhì)中進(jìn)行試驗(yàn),所得動(dòng)導(dǎo)數(shù)為正,M3模型不穩(wěn)定;在CO2介質(zhì)中進(jìn)行試驗(yàn),動(dòng)導(dǎo)數(shù)為負(fù),模型穩(wěn)定。三種介質(zhì)中,CO2分子質(zhì)量最大,空氣介質(zhì)分子質(zhì)量最小。分析認(rèn)為,隨著分子質(zhì)量的增大,氣體的可壓縮性增強(qiáng),通過(guò)試驗(yàn)得到的動(dòng)導(dǎo)數(shù)越小,動(dòng)穩(wěn)定性越強(qiáng)。

圖1 M3模型動(dòng)態(tài)試驗(yàn)數(shù)據(jù)[9]Fig.1 Dynamic test data[9]of M3 model

2.2 M3模型計(jì)算條件及網(wǎng)格

數(shù)值計(jì)算模擬彈道靶試驗(yàn)中的馬赫數(shù)與雷諾數(shù),表1給出了M3模型的彈道靶試驗(yàn)狀態(tài),采用三種不同的氣體介質(zhì),分別為純CO2、CO2混合少量空氣以及空氣介質(zhì),研究試驗(yàn)氣體成分對(duì)氣動(dòng)特性特別是俯仰振蕩的影響。表中,帶?號(hào)的為空氣介質(zhì),其余狀態(tài)的試驗(yàn)氣體介質(zhì)采用純CO2。

表1 M3模型彈道靶試驗(yàn)狀態(tài)Table 1 Test condition of ballistic range of M3 model

根據(jù)質(zhì)心位置的不同,將彈道靶試驗(yàn)分為M3A、M3B、M3C、M3Ax以及M3A(high-M)五類狀態(tài),其中,前四類狀態(tài)的質(zhì)心位置有所不同,M3Ax狀態(tài)中,質(zhì)心在M3A的基礎(chǔ)上再偏離10%,質(zhì)心的具體位置如表2所示。M3A(high-M)與M3A狀態(tài)中質(zhì)心位置相同,馬赫數(shù)存在一定的差異,high-M狀態(tài)中馬赫數(shù)高于M3A狀態(tài)中馬赫數(shù)。

表2 彈道靶模型的慣性特征Table2 Inertia characteristics of ballistic range model

由于動(dòng)態(tài)試驗(yàn)條件的缺失,本文根據(jù)彈道靶試驗(yàn)狀態(tài),考察了多種減縮頻率以及時(shí)間步長(zhǎng)對(duì)動(dòng)態(tài)結(jié)果的影響。具體計(jì)算條件見表3。

表3 M3模型動(dòng)態(tài)數(shù)值計(jì)算條件Table 3 Dynamic numerical conditions of M3 model

M3外形采用含奇性軸的網(wǎng)格進(jìn)行數(shù)值計(jì)算,M3網(wǎng)格如圖2所示,半場(chǎng)網(wǎng)格規(guī)模為90萬(wàn),壁面最小網(wǎng)格間距0.03 mm,網(wǎng)格雷諾數(shù)范圍為4~10。

2.3 M3外形數(shù)值驗(yàn)證

圖3展示了本文數(shù)值模擬結(jié)果與國(guó)外數(shù)值模擬結(jié)果以及經(jīng)驗(yàn)公式所得升力系數(shù)的對(duì)比驗(yàn)證,結(jié)果吻合。

圖2 計(jì)算網(wǎng)格Fig.2 Computing mesh

圖3 升力系數(shù)對(duì)比圖Fig.3 Comparison of lift coefficients

圖4為不同馬赫數(shù)狀態(tài)下的流場(chǎng)云圖,其中圖4(a)為M3A類對(duì)應(yīng)的一個(gè)計(jì)算狀態(tài),馬赫數(shù)為2.69,迎角為15.8°,圖4(b)為M3A(high-M)類對(duì)應(yīng)的一個(gè)計(jì)算狀態(tài),馬赫數(shù)為3.56,迎角為15.9°。由圖可知,流線穿過(guò)脫體激波,壁面極限流線由駐點(diǎn)向四周散開,前體上半部流線繞過(guò)肩部耳片部件,出現(xiàn)內(nèi)折角,產(chǎn)生膨脹波,向后流動(dòng)加速,與繞過(guò)底部的尾流匯合,形成與來(lái)流方向平行的二次壓縮波,耳片后側(cè)形成大范圍回流區(qū)。馬赫數(shù)對(duì)M3模型的脫體激波位置存在影響,馬赫數(shù)越大,脫體激波越靠近物面,尾部回流區(qū)域變窄。圖5展示了本文數(shù)值模擬結(jié)果與彈道靶試驗(yàn)結(jié)果以及文獻(xiàn)數(shù)值模擬結(jié)果進(jìn)行對(duì)比,其中,test為不確定度為3倍標(biāo)準(zhǔn)差的彈道靶試驗(yàn)結(jié)果,ref-CFD為文獻(xiàn)中采用GASP程序進(jìn)行靜態(tài)數(shù)值預(yù)測(cè)的結(jié)果。結(jié)果表明,本文的數(shù)值結(jié)果與彈道靶試驗(yàn)結(jié)果基本吻合。

針對(duì)表3中的計(jì)算狀態(tài),開展M3模型俯仰強(qiáng)迫振動(dòng)數(shù)值模擬計(jì)算,獲得了俯仰角-俯仰力矩遲滯曲線,如圖6(a,b)所示。

根據(jù)遲滯曲線,通過(guò)參數(shù)計(jì)算辨識(shí)方法獲得了俯仰動(dòng)導(dǎo)數(shù)結(jié)果,將結(jié)果與彈道靶試驗(yàn)結(jié)果進(jìn)行對(duì)比,試驗(yàn)不確定度為3倍標(biāo)準(zhǔn)差(3σ),結(jié)果如圖6(c,d)所示。結(jié)果表明,在所考察的彈道靶試驗(yàn)狀態(tài)下,M3模型的動(dòng)導(dǎo)數(shù)均小于0,飛行穩(wěn)定;減縮頻率對(duì)參數(shù)辨識(shí)的結(jié)果影響較小,時(shí)間步長(zhǎng)對(duì)動(dòng)導(dǎo)數(shù)計(jì)算辨識(shí)結(jié)果存在一定的影響,時(shí)間步長(zhǎng)取無(wú)量綱時(shí)間0.02時(shí),動(dòng)導(dǎo)數(shù)的結(jié)果更接近試驗(yàn)值。數(shù)值結(jié)果與彈道靶試驗(yàn)的結(jié)果基本吻合,本文發(fā)展的FLY3D-MARS軟件適用于火星大氣環(huán)境的動(dòng)態(tài)數(shù)值模擬計(jì)算。

圖4 不同馬赫數(shù)下的流場(chǎng)云圖Fig.4 Flow field contour under different Mach number

圖5 氣動(dòng)力系數(shù)與彈道靶試驗(yàn)結(jié)果以及文獻(xiàn)數(shù)值模擬結(jié)果對(duì)比圖Fig.5 Comparison of ballistic range test results and dynamic force coefficient

圖6 M3模型俯仰強(qiáng)迫振動(dòng)數(shù)值模擬計(jì)算結(jié)果Fig.6 Numerical simulation results of M3 model forced pitching motion

3 模型俯仰強(qiáng)迫振蕩及動(dòng)態(tài)流場(chǎng)分析

本節(jié)通過(guò)數(shù)值模擬M1模型在火星環(huán)境中的俯仰強(qiáng)迫振蕩,研究探測(cè)器前、后體以及壓縮性效應(yīng)等在俯仰動(dòng)態(tài)穩(wěn)定性中所起的作用和流動(dòng)機(jī)理。為了保證時(shí)間精度,采用雙時(shí)間步隱式方法求解非定常繞流問(wèn)題。振蕩幅值為1°,減縮頻率為0.13。遠(yuǎn)場(chǎng)采用適用于動(dòng)態(tài)邊界的無(wú)反射邊界條件,壁面無(wú)滑移邊界,物面溫度采用等溫壁假設(shè),Twall＝1500 K。

3.1 壓縮性效應(yīng)對(duì)高超階段靜、動(dòng)導(dǎo)數(shù)影響規(guī)律

圖7為火星環(huán)境中,不同初始迎角下的俯仰角和俯仰力矩時(shí)間歷程曲線,圖8為馬赫數(shù)17.59,近似比熱比取值1.15時(shí),不同迎角下動(dòng)態(tài)流場(chǎng)分布。

圖7 俯仰角、俯仰力矩時(shí)間歷程曲線(Ma＝17.59)Fig.7 Time history of pitching angle and moment(Ma＝17.59)

上述結(jié)果表明,俯仰力矩曲線出現(xiàn)相位滯后。探測(cè)器在振蕩過(guò)程中,每一次位置的調(diào)整都是對(duì)周圍流場(chǎng)的一種擾動(dòng),流場(chǎng)達(dá)到穩(wěn)定需要一定的延遲時(shí)間。前體大部分為超聲速或高超聲速流動(dòng),所需的調(diào)整時(shí)間很短;后體流場(chǎng)中,頸部點(diǎn)附近的流場(chǎng)受到剪切層、二次壓縮波、回流區(qū)內(nèi)流場(chǎng)以及尾流流場(chǎng)之間的相互作用,流場(chǎng)十分復(fù)雜,受到擾動(dòng)后,需要一定的調(diào)整時(shí)間,所以后體流場(chǎng)存在較大的滯后。隨著迎角的增加,滯后相位越大。

圖8 動(dòng)態(tài)流場(chǎng)分布(Ma＝17.59)Fig.8 Dynamic flow field contour(Ma＝17.59)

圖9為初始迎角為15°時(shí),不同比熱比條件下俯仰角和俯仰力矩的時(shí)間歷程曲線。相同迎角下,隨著比熱比的增加,俯仰角和俯仰力矩系數(shù)的滯后相位越大。

圖10是在火星環(huán)境中,不同初始迎角下的迎角-俯仰力矩遲滯曲線,遲滯曲線的形狀存在明顯的差異。遲滯曲線的面積和旋轉(zhuǎn)方向決定了動(dòng)導(dǎo)數(shù)的數(shù)值和符號(hào),進(jìn)而影響動(dòng)導(dǎo)數(shù)的結(jié)果。圖11為初始迎角為15°時(shí),不同比熱比條件下迎角-俯仰力矩遲滯曲線,結(jié)果表明,相同計(jì)算狀態(tài)下,不同的比熱比條件的遲滯圈的面積不同,靜、動(dòng)導(dǎo)數(shù)存在差異。

根據(jù)遲滯曲線,通過(guò)參數(shù)計(jì)算辨識(shí)方法獲得了俯仰靜、動(dòng)導(dǎo)數(shù)。圖12為相同比熱比條件下,靜、動(dòng)導(dǎo)數(shù)隨迎角的變化曲線。結(jié)果表明:隨著迎角、比熱比等變化,靜、動(dòng)導(dǎo)數(shù)只存在數(shù)值大小上的改變,符號(hào)并未發(fā)生變化;在所考察的范圍內(nèi),探測(cè)器均為靜穩(wěn)定、動(dòng)穩(wěn)定狀態(tài)。

圖9 俯仰角、俯仰力矩時(shí)間歷程曲線(Ma＝17.59,α＝15°)Fig.9 Time history of pitching angle and moment(Ma＝17.59,α＝15°)

圖10 火星環(huán)境不同初始迎角下的遲滯圈(Ma＝17.59,γ＝1.15)Fig.10 Hysteresis loops of different initial angles under Mars environment(Ma＝17.59,γ＝1.15)

圖11 不同比熱比條件的遲滯圈(Ma＝17.59,α＝15°)Fig.11 Hysteresis loops of different specific heat ratio(Ma＝17.59,α＝15°)

圖12 俯仰靜、動(dòng)導(dǎo)數(shù)隨迎角的變化曲線(Ma＝17.59)Fig.12 Changing curves of pitching static/dynamic derivative according to attack angle(Ma＝17.59)

圖13為相同迎角條件下,靜、動(dòng)導(dǎo)數(shù)隨比熱比的變化曲線。由靜導(dǎo)數(shù)結(jié)果(圖13a)可知,比熱比對(duì)靜導(dǎo)數(shù)結(jié)果存在影響,0°、2°迎角下,靜導(dǎo)數(shù)隨比熱比的增加而減小,靜穩(wěn)定性增強(qiáng);5°迎角下,靜導(dǎo)數(shù)隨比熱比先增大后減小,靜穩(wěn)定性先減弱后增強(qiáng);8°~20°迎角下,靜導(dǎo)數(shù)隨比熱比的增加而增加,靜穩(wěn)定性減弱。由動(dòng)導(dǎo)數(shù)隨比熱比變化曲線(圖13b)可知,5°~15°迎角除外,其余狀態(tài)相同的計(jì)算條件下,動(dòng)導(dǎo)數(shù)隨著比熱比的增加,呈單調(diào)遞增變化,比熱比影響著動(dòng)導(dǎo)數(shù)的結(jié)果;隨著比熱比的增加,動(dòng)導(dǎo)數(shù)增加,動(dòng)穩(wěn)定性減弱。分析認(rèn)為,5°~15°迎角時(shí),由于不同比熱比狀態(tài)下的參考時(shí)間不一致,導(dǎo)致動(dòng)導(dǎo)數(shù)的單調(diào)性不明顯。

圖13 俯仰靜、動(dòng)導(dǎo)數(shù)隨比熱比的變化曲線(Ma＝17.59)Fig.13 Changing curves of pitching static/dynamic derivative according to attack angle(Ma＝17.59)

圖14為馬赫數(shù)17.59、迎角20°、不同比熱比條件下的動(dòng)態(tài)流場(chǎng)結(jié)果,圖中顯示,不同比熱比條件下脫體激波位置、通過(guò)激波后的壓強(qiáng)存在差異:比熱比越小,脫體激波越靠近物面,波后壓強(qiáng)越大。圖15為不同比熱比條件下,前、后體子午線動(dòng)態(tài)壓力分布曲線,結(jié)果表明,前體壁面駐點(diǎn)以下部分,比熱比越小,壓力越大,體現(xiàn)了壓縮性效應(yīng)的影響。

圖16為馬赫數(shù)17.59、迎角20°,兩種不同比熱比條件下的俯仰強(qiáng)迫振蕩瞬態(tài)壁面極限流線分布。由圖中可知,迎風(fēng)面流線穿過(guò)脫體激波后改變流動(dòng)方向,繞過(guò)物面向四周發(fā)散,流經(jīng)物面的流線在第二段倒錐尾段物面匯合,后體下部形成半結(jié)點(diǎn)。

圖14 不同比熱比流場(chǎng)結(jié)果(Ma＝17.59,α＝20°)Fig.14 Flow field contour under different specific heat ratio(Ma＝17.59,α＝20°)

圖15 不同比熱比條件下的動(dòng)態(tài)壓力分布(Ma＝17.59,α＝20°)Fig.15 Pressure distribution under different specific heat ratio(Ma＝17.59,α＝20°)

圖16 動(dòng)態(tài)流線分布(Ma＝17.59,α＝20°)Fig.16 Dynamic streamline distribution(Ma＝17.59,α＝20°)

3.2 壓縮性效應(yīng)對(duì)超聲速階段靜、動(dòng)導(dǎo)數(shù)影響規(guī)律

由于超聲速環(huán)境下的近似比熱比取值與高超聲速環(huán)境中有所不同,本節(jié)針對(duì)超聲速計(jì)算狀態(tài)(馬赫數(shù)2.5),根據(jù)不同迎角下的比熱比(1.15、1.3、1.3755)狀態(tài)展開了數(shù)值模擬,得到了俯仰角-俯仰力矩遲滯曲線隨迎角、比熱比的變化,如圖17和18所示。通過(guò)參數(shù)計(jì)算辨識(shí)方法所得的俯仰靜、動(dòng)導(dǎo)數(shù)結(jié)果如圖19所示,結(jié)果表明,超聲速狀態(tài),相同比熱比條件下,隨著迎角的增加,靜導(dǎo)數(shù)增大,靜穩(wěn)定性減弱。在所考察的幾種計(jì)算狀態(tài)下,動(dòng)導(dǎo)數(shù)的符號(hào)均為負(fù),彈道靶模型為動(dòng)穩(wěn)定狀態(tài)。火星環(huán)境中的近似比熱比取值1.3,與此同時(shí),γ取1.3755最接近于地球大氣比熱比1.4,兩種比熱比條件下,通過(guò)參數(shù)計(jì)算辨識(shí)方法所得結(jié)果最接近,由此可知,相同迎角下,地球環(huán)境所得動(dòng)導(dǎo)數(shù)的數(shù)值小于火星環(huán)境中動(dòng)導(dǎo)數(shù)的數(shù)值,即,在超聲速階段相同計(jì)算狀態(tài)下,火星環(huán)境中探測(cè)器具有較強(qiáng)的動(dòng)穩(wěn)定性。由此可以得出結(jié)論:在所考察范圍內(nèi),隨著比熱比的減小,流體壓縮性增強(qiáng),動(dòng)穩(wěn)定性增強(qiáng),與彈道靶試驗(yàn)所得結(jié)論一致(圖1)。

圖20為不同比熱比條件的動(dòng)態(tài)流場(chǎng)分布,結(jié)果表明,超聲速狀態(tài)中,比熱比越小,脫體激波越靠近物面,波后壓強(qiáng)越大,與高超聲速狀態(tài)下所得結(jié)論一致。

圖17 火星環(huán)境不同初始迎角下的遲滯圈(Ma＝2.5,γ＝1.3)Fig.17 Hysteresis loops of different initial angles under Mars environment(Ma＝2.5,γ＝1.3)

圖18 不同比熱比條件的遲滯圈(Ma＝2.5,α＝15°)Fig.18 Hysteresis loops of different specific heat ratio(Ma＝2.5,α＝15°)

圖19 俯仰靜、動(dòng)導(dǎo)數(shù)隨迎角的變化曲線(Ma＝2.5)Fig.19 Changing curves of pitching static/dynamic derivative according to attack angle(Ma＝2.5)

圖20 不同比熱比的動(dòng)態(tài)流場(chǎng)分布(Ma＝2.5,α＝15°)Fig.20 Dynamic flow field distribution under different specific heat ratio(Ma＝2.5,α＝15°)

4 總 結(jié)

本文考核了FLY3D-MARS軟件在火星環(huán)境動(dòng)態(tài)特性計(jì)算中的適用性,進(jìn)一步確保流動(dòng)模型的準(zhǔn)確性;通過(guò)開展M1模型俯仰強(qiáng)迫振蕩數(shù)值模擬,系統(tǒng)分析壓縮性效應(yīng)對(duì)動(dòng)態(tài)特性的影響規(guī)律,結(jié)果表明,在所考察范圍內(nèi),比熱比越小,壓縮性增強(qiáng),動(dòng)導(dǎo)數(shù)減小,動(dòng)穩(wěn)定性增強(qiáng)。導(dǎo)致動(dòng)穩(wěn)定性隨γ變化的原因是不同γ流場(chǎng)壓縮性不同,與彈道靶試驗(yàn)結(jié)論一致。