彈道靶火星探測器試驗模型動態特性研究

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

(1.空氣動力學國家重點實驗室,綿陽 621000;2.中國空氣動力研究與發展中心 計算空氣動力研究所,綿陽 621000)

0 引 言

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

國內外火星環境的動態特性研究方法與地球環境的動態特性研究方法類似,主要有半經驗理論方法、數值模擬方法和火星風洞實驗方法。

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

火星探測的飛行數據稀缺,地面試驗難度大,工程計算方法進展緩慢,目前,數值模擬可作為火星探測氣動問題研究的重要工具。其中,近似比熱比數值模擬方法很大程度地簡化了計算模型[10-11],但國內外在近似比熱比的取值上有所不同,且研究尚顯不足。

因此,著眼于我國火星探測項目需求、亟需解決探測器進入火星大氣環境中所遇到的氣動問題,開展火星探測器氣動穩定性的研究。本文采用近似比熱比數值模擬方法,基于軟件in-house programm FLY3DMars,針對火星探測器彈道靶自由飛試驗外形,進行了火星大氣模擬環境下俯仰氣動特性的數值模擬研究。

1 數值模擬方法

1.1 控制方程

直角坐標系中,不計體積力和外加熱量的無量綱化、守恒形式Navier-Stokes方程組為:

其中,Q為守恒變量,E、F和G為對流項、無黏通矢量,Ev、Fv和Gv為黏性通矢量,具體形式分別為:

通過坐標變換,可將方程組(1)轉化為一般曲線坐標系中的形式:

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

1.2 動網格技術

動態非定常數值模擬與定態數值模擬的區別之一是需要建立相對于地面系的動坐標系,以描述飛行器的姿態和運動。本文使用的地面系與飛行力學習慣建立的地面系之間的差別在于:y軸方向不變,x、z軸方向相反。同樣,本文使用的體軸系的差別在于:yb軸方向不變,xb、zb軸方向相反。

一般俯仰運動計算時,可使用半流場網格;偏航、滾轉運動計算,使用全流場網格。對于網格運動要求滿足幾何守恒律,即:

由于計算非定常動態流場需要每一時間步都更新網格,因此動網格生成技術的計算量小就顯得尤為重要。目前,剛性運動網格生成和超限插值生成動網格是兩種常用的方法。剛性運動網格生成技術就是計算網格隨物體一起作剛體運動,其計算量小,并可保持初始網格的品質,但不適于生成較復雜的網格,也使運動開邊界的處理不易。超限插值生成動網格方法就是外邊界保持靜止,物面邊界由物體運動規律或運動方程得到,內場網格由超限插值的方法代數生成,其計算量也較小,能夠生成復雜網格,但不易保證網格品質,尤其是不能保證物面網格的正交性?？捎眉訖嗉夹g來綜合這兩種方法的優點[3]。

首先,由初始網格Xn生成剛性運動網格Xref1,再由n時刻的開邊界和n＋1時刻的物面邊界代數插值得到n＋1時刻的網格Xref2,最后由兩者加權得到新時刻的動網格:

只要適當選擇加權因子ω,就既可保證物面附近網格質量,又可使外邊界保持靜止而易于處理?？蓸嬙烊缦录訖嘁蜃?計算表明效果較好,令j為法向網格指標,則有:

1.3 靜、動態氣動穩定性參數計算辨識方法

在預測靜、動態氣動穩定性參數方面,數值計算明顯優于地面模擬試驗方法:數值計算費用低,可提供較長的非定常動態氣動力周期,振心位置可與實際重心位置一致[6]。靜、動導數參數辨識的思路在火星環境中仍然適用,主要是將動態氣動力進行Taylor級數展開,保留線性部分,以俯仰運動為例,對于單自由度俯仰運動,根據Etkin經典動態氣動力模型,有

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

本文采用強迫振動法求取動導數。以單自由度俯仰強迫振動為例,給定簡諧振蕩(無量綱化)形式如下:

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

其中,T＝π/k為強迫振蕩的周期。強迫振蕩的振幅對辨識結果的影響不大,一般取為1°。由于飛行器在空中的真實頻率并不知道,需要考察多個頻率對動穩定性方面的影響,才能得到有效結果。單自由度強迫偏航/滾轉及其靜、動導數的辨識類似。

2 算例驗證

第2節針對非軸對稱M3外形,開展俯仰強迫振動數值模擬,通過將數值模擬結果與彈道靶試驗結果對比,考核in-house programm FLY3D-MARS軟件在火星環境動態計算中的適用性。

2.1 MSL彈道靶試驗動態結果

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

圖1 M3模型動態試驗數據[9]Fig.1 Dynamic test data[9]of M3 model

2.2 M3模型計算條件及網格

數值計算模擬彈道靶試驗中的馬赫數與雷諾數,表1給出了M3模型的彈道靶試驗狀態,采用三種不同的氣體介質,分別為純CO2、CO2混合少量空氣以及空氣介質,研究試驗氣體成分對氣動特性特別是俯仰振蕩的影響。表中,帶?號的為空氣介質,其余狀態的試驗氣體介質采用純CO2。

表1 M3模型彈道靶試驗狀態Table 1 Test condition of ballistic range of M3 model

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

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

由于動態試驗條件的缺失,本文根據彈道靶試驗狀態,考察了多種減縮頻率以及時間步長對動態結果的影響。具體計算條件見表3。

表3 M3模型動態數值計算條件Table 3 Dynamic numerical conditions of M3 model

M3外形采用含奇性軸的網格進行數值計算,M3網格如圖2所示,半場網格規模為90萬,壁面最小網格間距0.03 mm,網格雷諾數范圍為4~10。

2.3 M3外形數值驗證

圖3展示了本文數值模擬結果與國外數值模擬結果以及經驗公式所得升力系數的對比驗證,結果吻合。

圖2 計算網格Fig.2 Computing mesh

圖3 升力系數對比圖Fig.3 Comparison of lift coefficients

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

針對表3中的計算狀態,開展M3模型俯仰強迫振動數值模擬計算,獲得了俯仰角-俯仰力矩遲滯曲線,如圖6(a,b)所示。

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

圖4 不同馬赫數下的流場云圖Fig.4 Flow field contour under different Mach number

圖5 氣動力系數與彈道靶試驗結果以及文獻數值模擬結果對比圖Fig.5 Comparison of ballistic range test results and dynamic force coefficient

圖6 M3模型俯仰強迫振動數值模擬計算結果Fig.6 Numerical simulation results of M3 model forced pitching motion

3 模型俯仰強迫振蕩及動態流場分析

本節通過數值模擬M1模型在火星環境中的俯仰強迫振蕩,研究探測器前、后體以及壓縮性效應等在俯仰動態穩定性中所起的作用和流動機理。為了保證時間精度,采用雙時間步隱式方法求解非定常繞流問題。振蕩幅值為1°,減縮頻率為0.13。遠場采用適用于動態邊界的無反射邊界條件,壁面無滑移邊界,物面溫度采用等溫壁假設,Twall＝1500 K。

3.1 壓縮性效應對高超階段靜、動導數影響規律

圖7為火星環境中,不同初始迎角下的俯仰角和俯仰力矩時間歷程曲線,圖8為馬赫數17.59,近似比熱比取值1.15時,不同迎角下動態流場分布。

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

上述結果表明,俯仰力矩曲線出現相位滯后。探測器在振蕩過程中,每一次位置的調整都是對周圍流場的一種擾動,流場達到穩定需要一定的延遲時間。前體大部分為超聲速或高超聲速流動,所需的調整時間很短;后體流場中,頸部點附近的流場受到剪切層、二次壓縮波、回流區內流場以及尾流流場之間的相互作用,流場十分復雜,受到擾動后,需要一定的調整時間,所以后體流場存在較大的滯后。隨著迎角的增加,滯后相位越大。

圖8 動態流場分布(Ma＝17.59)Fig.8 Dynamic flow field contour(Ma＝17.59)

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

圖10是在火星環境中,不同初始迎角下的迎角-俯仰力矩遲滯曲線,遲滯曲線的形狀存在明顯的差異。遲滯曲線的面積和旋轉方向決定了動導數的數值和符號,進而影響動導數的結果。圖11為初始迎角為15°時,不同比熱比條件下迎角-俯仰力矩遲滯曲線,結果表明,相同計算狀態下,不同的比熱比條件的遲滯圈的面積不同,靜、動導數存在差異。

根據遲滯曲線,通過參數計算辨識方法獲得了俯仰靜、動導數。圖12為相同比熱比條件下,靜、動導數隨迎角的變化曲線。結果表明:隨著迎角、比熱比等變化,靜、動導數只存在數值大小上的改變,符號并未發生變化;在所考察的范圍內,探測器均為靜穩定、動穩定狀態。

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

圖10 火星環境不同初始迎角下的遲滯圈(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 俯仰靜、動導數隨迎角的變化曲線(Ma＝17.59)Fig.12 Changing curves of pitching static/dynamic derivative according to attack angle(Ma＝17.59)

圖13為相同迎角條件下,靜、動導數隨比熱比的變化曲線。由靜導數結果(圖13a)可知,比熱比對靜導數結果存在影響,0°、2°迎角下,靜導數隨比熱比的增加而減小,靜穩定性增強;5°迎角下,靜導數隨比熱比先增大后減小,靜穩定性先減弱后增強;8°~20°迎角下,靜導數隨比熱比的增加而增加,靜穩定性減弱。由動導數隨比熱比變化曲線(圖13b)可知,5°~15°迎角除外,其余狀態相同的計算條件下,動導數隨著比熱比的增加,呈單調遞增變化,比熱比影響著動導數的結果;隨著比熱比的增加,動導數增加,動穩定性減弱。分析認為,5°~15°迎角時,由于不同比熱比狀態下的參考時間不一致,導致動導數的單調性不明顯。

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

圖14為馬赫數17.59、迎角20°、不同比熱比條件下的動態流場結果,圖中顯示,不同比熱比條件下脫體激波位置、通過激波后的壓強存在差異:比熱比越小,脫體激波越靠近物面,波后壓強越大。圖15為不同比熱比條件下,前、后體子午線動態壓力分布曲線,結果表明,前體壁面駐點以下部分,比熱比越小,壓力越大,體現了壓縮性效應的影響。

圖16為馬赫數17.59、迎角20°,兩種不同比熱比條件下的俯仰強迫振蕩瞬態壁面極限流線分布。由圖中可知,迎風面流線穿過脫體激波后改變流動方向,繞過物面向四周發散,流經物面的流線在第二段倒錐尾段物面匯合,后體下部形成半結點。

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

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

圖16 動態流線分布(Ma＝17.59,α＝20°)Fig.16 Dynamic streamline distribution(Ma＝17.59,α＝20°)

3.2 壓縮性效應對超聲速階段靜、動導數影響規律

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

圖20為不同比熱比條件的動態流場分布,結果表明,超聲速狀態中,比熱比越小,脫體激波越靠近物面,波后壓強越大,與高超聲速狀態下所得結論一致。

圖17 火星環境不同初始迎角下的遲滯圈(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 俯仰靜、動導數隨迎角的變化曲線(Ma＝2.5)Fig.19 Changing curves of pitching static/dynamic derivative according to attack angle(Ma＝2.5)

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

4 總 結

本文考核了FLY3D-MARS軟件在火星環境動態特性計算中的適用性,進一步確保流動模型的準確性;通過開展M1模型俯仰強迫振蕩數值模擬,系統分析壓縮性效應對動態特性的影響規律,結果表明,在所考察范圍內,比熱比越小,壓縮性增強,動導數減小,動穩定性增強。導致動穩定性隨γ變化的原因是不同γ流場壓縮性不同,與彈道靶試驗結論一致。