火星進入器壁面脈動壓力環境數值模擬

石小潘，趙 瑞，榮吉利，袁 武

(1. 北京理工大學宇航學院，北京 100081；2. 中國科學院計算機網絡信息中心超級計算中心，北京 100190)

0 引 言

隨著航天技術的不斷進步和發展，人類已經展開了200多次的深空探測任務[1]。而火星作為距離地球最近的地外行星，和其他行星相比其自然環境更接近地球，因此，火星成為人類探索宇宙重要的一部分。火星和地球一樣存在大氣層，火星大氣主要由95.7%的CO2,2.7%的N2和1.6%的Ar組成，密度只有地球的1%，與地球大氣組分顯著不同。文獻[2]指出大氣中的主要氣體成分的物理化學性質對氣動特性有顯著影響，因此，進入器進入火星大氣與再入地球大氣時受到的氣動環境會有所差異。目前，各國開展火星無人進入或軟著陸探測任務時，發射的進入器均采用鈍體大底-倒錐布局的外形[3]，優點是在高超聲速進入時其布局外形的高阻力特性有利于進入器的氣動減速。而進入器在進入火星大氣層時，由于鈍體繞流效應，流體繞過肩部后出現分離，在艙體尾段及配平翼迎風面底端形成了分離區，分離區內渦流脈動劇烈，將會在進入器壁面產生振幅較高的脈動壓力，其頻率范圍包含了普通金屬面板的共振頻率。一方面，脈動壓力可使進入器壁面出現較大的局部載荷，容易激起飛行器結構的抖振響應，縮短材料的疲勞壽命；另一方面，脈動壓力作為一種隨機激勵以噪聲的形式通過透射及結構共振轉變為艙內噪聲，直接影響到內部儀器的可靠性，甚至導致任務失敗[4-5]。因此，研究火星進入器壁面的脈動壓力環境，對火星進入器的結構設計及艙內聲振環境的預測都具有重要的意義。

盡管若干火星著陸任務已經完成，但是與火星進入器相關的飛行試驗數據依然很匱乏，并且模擬火星大氣環境進行風洞試驗難度大且成本高，因此對火星進入器的脈動壓力環境的研究手段主要靠數值計算[6]。國外對飛行器脈動壓力環境作了諸多研究，Plotkin等[7]通過對各種實驗數據的分析，總結歸納了一套適用于跨聲速范圍內非定常流動脈動壓力環境的經驗公式，給出了均方根脈動壓力系數、功率譜及空間相關性的計算方法。Tsutsumi等[8]采用改進的延遲脫體渦方法對典型整流罩進行非定常數值模擬，通過與實驗值進行對比證明該方法的可行性。文獻[9-11]采用脫體渦方法對亞聲速、跨聲速及超聲速的鈍體繞流流場作了數值模擬，模擬結果整體和實驗數據相近，且在跨超聲速精度最好。Ross等[12]針對在亞跨聲速范圍內鈍體繞流非定常分離難以精確計算的問題，通過風洞實驗分析了雷諾數、來流馬赫數、攻角等對標準艙體外形脈動壓力的影響。在國內，徐立功等[13]對機動再入飛行器壁面所受的空氣動力脈動環境進行了分析，用統一的參變量給出了高超聲速范圍內飛行器壁面脈動壓力統計特性的工程計算公式。張志成等[14]使用工程估算公式分析了來流攻角和來流馬赫數對球頭雙錐體表面脈動壓力分布的影響。蔣華兵等[15]通過對典型球頭單錐飛行器定常時均流場的求解，獲得湍流附面層厚度及外緣流動參數，并代入經驗公式，對飛行器噪聲環境進行計算并指出了飛行器表面噪聲環境惡劣的區域。趙瑞等[16]分析了火箭整流罩外氣動噪聲環境，用隱式大渦模擬的方法對跨聲速條件下火箭整流罩外部氣動噪聲進行數值分析，并改進了經驗公式在分離區的預測精度；此外，趙瑞等[5]對跨聲速旋成體火箭脈動壓力產生的流場結構進行重新分區歸納，針對各個分區脈動流場提出改進的經驗公式，提高預測精度。

綜上所述，各國學者對飛行器脈動壓力環境作了不同的研究，并取得了顯著的成果，但對火星大氣下進入器的脈動壓力環境研究較少。進入器著陸火星表面要經歷四個階段：進入準備階段、高超聲速氣動減速下降階段、打開降落傘減速下降階段和終點下降著陸階段。本文將對進入器在開傘之前跨超聲速氣動減速下降階段的脈動壓力環境展開研究分析。基于脫體渦模擬(Detached eddy simulaion, DES)方法，通過對不同工況下進入器繞流進行非定常數值模擬，分析來流馬赫數、來流攻角及配平翼展開角的變化對進入器壁面脈動壓力環境的影響，綜合評價進入器的脈動壓力環境。

1 計算方法和模型

1.1 控制方程

流動控制方程為三維非定常可壓縮N-S方程組，其在直角坐標系下的守恒積分形式可表示為：

(1)

式中：U為控制體， ?U為控制面，Q為守恒變矢量，F為矢通量，n為邊界的外法向量。

采用有限體積法對控制方程進行數值離散，在數值離散過程中，空間對流項離散使用5階的WENO格式，黏性項使用4階的中心差分方法[17]。時間推進采用二階隱式雙時間法進行迭代[18]，內迭代選擇常用的LU-SGS(lower-upper symmetric-Gauss-Seidel )方法[19]。定常計算采用一方程S-A湍流模型并將其計算結果作為非定常DES計算的初始流場。

目前對于火星大氣介質的模擬，有兩種方法：一種方法是采用真實氣體模型[2]，考慮氣體介質中的各種化學組分；另一種方法是等效比熱比模型[20]，根據統計熱力學的相關理論[21]，標準狀態的各種氣體介質的熱力學參數都可用峰值能量分配原理來描述，也就能用有效比熱比、氣體常數和溫度等變量來表征。離解、電離等真實氣體效應亦可看作能量分配方式的增加，進而帶來更低的有效比熱比，因此，非空氣介質、真實氣體效應等影響可以用有效比熱比來等效和表征。本文使用的即是基于等效比熱比方法進行數值計算，其等效比熱比取為1.17[20]。另外，對于火星大氣環境，以CO2為主，其黏性系數隨溫度而改變，使用Sutherland公式來表征

(2)

式中：參考黏性系數μ0=1.48×10-5kg/m·s，參考溫度T0=293.15 K。和黏性力一樣，氣體微團之間的熱傳導也是由分子運動造成，導熱系數和黏性系數之間存在相似關系，在處理高速黏性流動問題，導熱系數由Pr=μcp/k確定，Pr=0.71。

1.2 基于S-A模型的DES方法

(3)

方程(3)右端的生成項變量定義為：

其中，Ω是渦量，d為到物面的最近距離,函數fW如下所示：

在S-A模型中定義湍流黏性系數為：

(4)

DES方法通過對Spalart-Allmaras湍流模型進行長度尺度修正，使其在近壁面采用雷諾平均方法(RANS)數值模擬附面層的流動，即用湍流模型來模擬近壁面的小尺度脈動；遠離物面附面層的區域DES轉換至Smagorinsky大渦模擬，成為亞格子雷諾應力模型。

(5)

(6)

1.3 計算模型和網格

本文主要研究進入器進入火星大氣層之后，在減速降落傘打開之前的跨超聲階段來流馬赫數(Ma)、來流攻角(α)及配平翼展開角(φ)對火星進入器壁面脈動壓力環境的影響。計算工況如表1所示。

表1 進入器計算條件Table 1 Computation condition for the capsule

計算模型為類似文獻[23]中所示的具有配平翼的進入器外形，基本尺寸如圖1所示。計算網格為六面體結構化網格，如圖2所示，網格總量約600萬，并在分離區附近(配平翼及艙體后方)進行加密，用于捕捉脫體渦結構。壁面第一層網格保證y+<1，用于捕捉黏性邊界層流動。

對進入器z=0截面進行測點采樣，如圖3所示，采樣點共計29個，采樣頻率為5000 Hz，從瞬時流場中提取壁面特征點壓力信號隨時間變化歷程曲線，獲得均方根脈動壓力系數，并采用快速傅里葉變換(FFT)，研究頻域空間中脈動壓力的功率譜密度分布。

1.4 數值方法驗證

為驗證本文使用的計算方法，以具有風洞試驗數據的神舟返回艙模型作為驗證算例。神舟返回艙風洞外形及壁面測點的位置見圖4，計算網格約為1000萬，計算參數選自風洞實驗參數見表2。

T0/KP0/PaMaU∞/(m·s-1)p∞/Paq∞/Pa(Re/L)/m-12871012090.6200.379297200471.213×107

2 結果與分析

2.1 來流馬赫數對脈動壓力環境的影響

均方根脈動壓力系數是衡量壁面脈動壓力總強度的關鍵統計參數。圖7為在來流攻角為0°,配平翼展開180°時不同來流馬赫數下進入器壁面各測點的均方根脈動壓力系數對比曲線。

從圖7可以看出，進入器整體的脈動壓力環境隨來流馬赫數的增加呈降低的趨勢，在馬赫數為1.2時，均方根脈動壓力系數第一峰值出現在艙體背風面下側，第二峰值出現在配平翼背風面；另外，配平翼迎風面翼根區及艙體背風面上側的脈動壓環境也十分嚴峻。圖8為不同來流馬赫數下進入器對稱面的馬赫數云圖。從圖8可以看出，由于鈍體效應，流動繞過肩部后出現分離，進入器艙體段淹沒在分離區中。由于分離區內渦流的非定常脈動，在艙體壁面誘導產生了脈動壓力。同時可以看到，在配平翼的迎風面翼根處也存在小的分離區。隨著馬赫數的增加，流動可壓縮效應增強，使得艙體分離區逐漸減小。當馬赫數較低時(見圖8(a))分離區較大，渦流脈動強烈，使得艙體背風面及配平翼背風面的均方根脈動壓力系數較大，脈動壓力環境惡劣。相反，隨著來流馬赫數的增加，來流氣體的可壓縮性逐漸增強，分離區也越來越小，分離區較短，導致渦流脈動空間不足、發展不充分(見圖8(b)和圖8(c))。因此，在配平翼背風區及艙體背風區的均方根脈動壓力系數較小，脈動壓力環境并不嚴峻。

在進入器艙體下側靠近肩部位置安裝有配平翼，下側氣流繞過肩部后經膨脹激波加速后沖擊在配平翼迎風面上，迎風面氣流在配平翼駐點處分為上下兩部分，上邊的氣流在艙翼連接的夾角范圍內形成了一個低壓分離區。在來流馬赫數為1.2時，鈍體前邊形成的脫體激波強度較小且不穩定，在肩部膨脹激波和翼尖高壓區的作用下，分離區極不穩定、運動劇烈，誘導出強烈的脈動壓力環境，致使配平翼迎風面翼根處均方根脈動壓力系數較大。而在來流馬赫數為2和3時，鈍體前方的脫體激波強度較大且穩定，并且從圖8可以看出，脫體激波隨著馬赫數的增加逐漸向大底靠近。配平翼迎風面的分離區受脫體激波抑制作用逐漸增強，分離區的運動逐漸保持穩定。因此，分離區在配平翼迎風面誘導的脈動壓力環境并不明顯。下邊的氣流繞過配平翼之后，經配平翼誘導在其背風面產生一個分離區。在來流馬赫數為1.2時，渦流運動劇烈，致使配平翼背風面的脈動壓力環境劇烈。該渦流流場向下游發展，和艙體脫體渦區域相互干擾，使艙體背風面下側靠近配平翼的區域脈動壓力達到最大，即均方根脈動壓力系數第一峰值所在的區域。相比艙體背風面下側，艙體背風面上側僅受艙體脫體渦的影響，其脈動壓力環境相對較小。

如上所述，進入器脈動壓力環境在低馬赫數下最為嚴峻，以下對攻角和配平翼展開角的研究主要基于Ma=1.2來流條件。

2.2 來流攻角對脈動壓力環境的影響

圖9為在來流馬赫數為1.2，配平翼展開180°時不同來流攻角下進入器壁面各測點的脈動壓力系數變化曲線。從圖9可以看出， 來流攻角從0°改變到15°，進入器整體的脈動壓力環境變化較大。

圖10為進入器在不同來流攻角下的馬赫數云圖。從圖10可以看出，隨著來流攻角的改變，進入器前方的脫體激波的位置向配平翼迎風面壓進，使得來流在配平翼迎風面上的再附點位置向翼根靠近，如圖11所示，再附點的位置改變使得分離區作用的區域向翼根縮進，從而使得在配平翼迎風面上誘導的脈動壓力峰值向翼根區轉移。另外，配平翼誘導的剪切層隨著攻角的改變向艙體下側面靠近，使得配平翼背風面及艙體背風面下側的分離區減小，渦流脈動空間減小、脈動減弱，使得均方脈動壓力系數相對減小；對于艙體背風面上側，分離區增大，渦流脈動空間增大、脈動增強，使得艙體上側面的均方根脈動壓力系數相對增大。

2.3 配平翼展開角度對脈動壓力環境的影響

圖12為在來流馬赫數為1.2,來流攻角0°時不同配平翼展開角度下進入器壁面各測點的均方根脈動壓力系數變化曲線。從圖12可以看出，配平翼角度變化對配平翼迎風面影響較大，對進入器其他區域影響較小。

圖13為不同配平翼展開角度下艙翼連接段的局部放大圖。從圖13可以看出，流動繞過下側肩部后出現分離，在配平翼展開角120°范圍內，配平翼和艙體段淹沒在分離區中，隨著配平翼的展開，配平翼迎風面的分離區逐漸減小，渦流脈動空間不足，脈動強度逐漸減弱，因此在配平翼迎風面上誘導產生的脈動壓力逐漸減小。相反，背風區的脈動壓力稍有增加。配平翼展開150°時，部分來流沖擊在配平翼迎風面上，來流受配平翼的影響，在艙翼連接段形成了一個低壓分離區，分離區再附點位置(見圖13(c))極不穩定，前后移動劇烈，將在配平翼迎風面上誘導產生脈動壓力峰值，均方根脈動壓力系數最大達到0.21。另外，再附點的前后移動加劇了分離區運動，將會在配平翼迎風面靠近翼根區域誘導嚴峻的脈動壓力環境。配平翼展開180°時，配平翼迎風面分離區再附點位置基本固定，相比配平翼展開150°，分離區的形態較穩定。因此，分離區在配平翼上誘導的脈動壓力得到減緩。

2.4 脈動壓力頻譜特性

功率譜密度表征脈動能量隨頻率的分布特性。通過對火星進入器脈動壓力的頻譜特性分析，可以確定脈動能量集中的頻率區間，進而對進入器進行針對的結構設計。針對進入器脈動壓力環境最惡劣的工況(來流馬赫數1.2，配平翼展開150°)，選取三個不同區域特征點進行功率譜密度分析，其中測點11位于配平翼迎風面分離區靠近翼根處，測點13位于配平翼再附點位置附近，測點25位于艙體段分離區。圖14為上述3個典型位置測點的功率譜密度頻譜變化曲線。從圖14可以看出，在測點11位置(見圖14(a))，該處位于艙翼連接段這個狹小的空間內，脈動頻率主要受分離區脈動頻率影響，在200 Hz左右達到峰值，脈動能量達10 Pa2·Hz-1。測點13與測點11類似，脈動能量集中頻率受分離區脈動頻率影響(見圖14(b))，但該位置處于再附點位置前后移動范圍，能量幅值在200 Hz左右達100 Pa2·Hz-1，相比其他區域高出一個量級，在結構設計中需要著重考慮。在測點25位置(見圖14(c))，艙體分離區內脈動壓力主要由艙體誘導的大分離渦流脈動所致，能量主要集中在低頻區域(10 Hz左右)。

3 結 論

本文使用DES方法研究來流馬赫數、來流攻角及配平翼展開角等因素對火星進入器壁面脈動壓力環境的影響規律，并得出以下結論：

1)采用DES方法對進入器壁面脈動壓力環境進行預測，能夠較為準確地獲得壁面脈動壓力的均方根系數以及頻譜特性。

2)在跨超聲速來流條件下，隨著來流馬赫數增加，艙體背風面分離區的可壓縮性逐漸增強，分離區逐漸減小，渦流脈動空間不足，致使艙體背風面的脈動壓力環境逐漸減弱。因此，在跨聲速時火星進入器壁面脈動壓力環境最為惡劣。

3)配平翼迎風面分離區受脫體激波的影響顯著：來流馬赫數較小時，脫體激波強度較小且不穩定，對翼根分離區抑制作用較弱，分離區運動劇烈，將會誘導強烈的脈動壓力環境；來流馬赫數較大時，脫體激波強度較大且穩定，對翼根分離區抑制作用較強，分離區運動平緩，誘導的脈動壓力環境較弱。

4)來流攻角改變使來流在配平翼迎風面上分離區再附點的位置向翼根方向轉移，使得當地的脈動壓力趨于惡劣。翼根區位于艙翼連接段，其結構穩定性在設計中需要著重注意。

5)配平翼展開150°時，配平翼迎風面分離區再附點的位置極不穩定、前后移動劇烈，將在配平翼迎風面誘導產生脈動壓力峰值，脈動能量主要集中在中頻區域；配平翼展開180°時，配平翼迎風面分離區再附點的位置基本固定，在配平翼迎風面上誘導的脈動壓力環境得到減緩。

參 考 文 獻

[1] 葉培建, 楊孟飛, 彭兢, 等. 中國深空探測進入/再入返回技術的發展現狀和展望[J]. 中國科學, 2015, 45(3):229-238. [Ye Pei-jian, Yang Meng-fei, Peng Jing, et al. Review and prospect of atmospheric entry and earth reentry technology of China deep space exploration [J]. Science China Press, 2015, 45(3): 229-238.]

[2] 呂俊明, 苗文博, 程小麗, 等. 火星大氣模型參數對MSL氣動特性的影響[J]. 空間科學學報, 2014, 34(4):377-383. [Lv Jun-ming, Miao Wen-bo, Cheng Xiao-li, et al. Impact of martian atmosphere model parameters on aerodynamic characteristics of Mars science laboratory [J]. Chinese Journal of Space Science, 2014, 34(4): 377-383.]

[3] 陳冰雁, 詹惠玲, 周偉江. 防熱大底外形對火星探測器氣動特性的影響分析[J]. 宇航學報, 2016, 37(4):388-396. [Chen Bing-yan, Zhan Hui-ling, Zhou Wei-jiang. Investigation on influence of heatshield shape changes on aerodynamic design for Mars entry capsule [J]. Journal of Astronautics, 2016, 37(4): 388-396.]

[4] 龍萬花, 陳偉芳, 宋松和. 旋成體跨聲速脈動壓力環境分析與預測[J]. 國防科技大學學報, 2004, 26(1):17-20. [Long Wan-hua, Chen Wei-fang, Song Song-he. Analysis and prediction of the fluctuating pressure induced by rotated aircraft at transonic Mach number [J]. Journal of National University Technology, 2004, 26(1): 17-20.]

[5] 趙瑞, 榮吉利, 任方, 等. 一種改進的跨聲速旋成體壁面脈動壓力經驗預測公式[J]. 宇航學報, 2016, 37(10):1179-1184. [Zhao Rui, Rong Ji-li, Ren Fang, et al. Improvement of the fluctuating pressure empirical functions for a body of revolution at transonic Mach numbers [J]. Journal of Astronautics, 2016, 37(10): 1179-1184.]

[6] Wright M J, Tang C Y, Edquist K T, et al. A review of aerothermal modeling for Mars entry mission[C]. The 48th AIAA Aerospace Science Meeting Including the New Horizons Forum and Aerospace Exposition, Orlando, USA, Jan 4-7, 2010.

[7] Plotkin K J, Robertson J E. Prediction of space shuttle fluctuating pressure environment, including rocket plume effects [R]. NASA-CR-124347, June 1973.

[8] Tsutsumi S, Takama Y, et al. Hybrid LES/RANS simulations of transonic flowfield around a rocket fairing[C]. The 30th AIAA Applied Aerodynamics Conference, New Orleans，Louisiana, USA, Jun 25-28, 2012.

[9] Alan M, Graham V. Validation of DES for capsule aerodynamics using 05-CA wind tunnel test data[C]. The 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Grapevine, Texas, USA,Jan 7-10, 2013.

[10] Alan M, Graham V. Detached-eddy simulation of capsule wake flows and comparison to wind-tunnel test data [J]. Journal of Spacecraft and Rockets, 2015, 52 (2):439-449.

[11] Josepf M, Pramod K, Graham V. Detached-eddy simulations of hypersonic capsule wake flow [J]. AIAA Journal, 2015, 53 (1):70-80.

[12] Ross J C, Burnside N J. Unsteady pressures on a generic capsule shape[C]. The 53rd AIAA Aerospace Sciences Meeting, Kissimmee, Florida, USA, Jan 5-9, 2015.

[13] 徐立功, 劉振寰． 再人飛行器脈動壓力環境的分析與預測[J]. 空氣動力學學報, 1991, 9(4):457-464. [Xu Li-gong, Liu Zhen-huan. Prediction of maneuvering reentry vehicles fluctuating pressure environments [J]. Acta Aerodynamic Sinica, 1991, 9(4): 457-464.]

[14] 張志成, 陳偉芳, 石于中, 等. 球頭雙錐再入體表面脈動壓力環境的分析與預測[J]. 宇航學報, 2002, 23(2):19-22. [Zhang Zhi-cheng, Chen Wei-fang, Shi Yu-zhong, et al. The analysis and prediction of fluctuating pressure on the spherebiconic reentry vehicle [J]. Journal of Astronautics, 2002, 23(2): 19-22.]

[15] 蔣華兵, 李春麗, 陳強洪. 再入飛行器脈動壓力環境特性分析[J]. 航天器環境工程, 2010, 27(3):378-382. [Jiang Hua-bing, Li Chun-li, Chen Qiang-hong. The characteristics of the fluctuating pressure environment for a re-entry vehicle [J]. Spacecraft Environment Technology, 2010, 27(3): 378-382.]

[16] 趙瑞, 榮吉利, 任方, 等. 火箭整流罩外氣動噪聲環境的大渦模擬研究[J]. 宇航學報, 2015, 36(9):988-994. [Zhao Rui, Rong Ji-li, Ren Fang, et al. Large eddy simulation of the aeroacoustic environment of a rocket fairing [J]. Journal of Astronautics, 2015, 36(9): 988-994.]

[17] Shen Y Q,Zha G G, Chen X Y. High order conservative differencing for viscous terms and the application to vortex-induced vibration flows [J]. Journal Computational Physics, 2008, 228(22):8283-8300.

[18] Dubuc L, Cantariti F, Woodgate M, et al. Solution of unsteady Euler equations using an implicit dual time step method [J]. AIAA Journal, 1998, 36(8): 1417-1424.

[19] Jameson A, Yoon S. Lower-upper implicit schemes with multiple grids for the Euler equations[J]. AIAA Journal, 1987, 25(7):929-935.

[20] 楊肖峰, 唐偉, 桂業偉. MSL火星探測器高超聲速流場預測及氣動性能分析[J]. 宇航學報, 2015, 36(4):384-389. [Yang Xiao-feng, Tang Wei, Gui Ye-wei. Hypersonic flow field prediction and aerodynamics analysis for MSL entry capsule [J]. Journal of Astronautics, 2015, 36(4): 384-389.]

[21] Anderson J D. Hypersonic and high-temperature gas dynamics[M]. Reston, Virginia: American Institute of Aeronautics and Astronautics, Inc. , 2006.

[22] Spalart P R, Allmaras S R. A one-equation turbulence model for aerodynamic flow [C]. The 30th Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, Jan 6-9, 1992.

[23] Andersen B M, Whitmore S A. Aerodynamic control on a lunar return capsule using trim-flaps[C].The 45th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada,USA, Jan 8-11, 2007.