火星大氣高溫光譜建模與非平衡輻射熱流預測

呂俊明，李飛，李齊，程曉麗

1.中國航天空氣動力技術研究院，北京 100074 2.中國航天科技集團有限公司 航天飛行器氣動熱防護實驗室，北京 100074 3.中國科學院 力學研究所 高溫氣體動力學國家重點實驗室，北京 100190 4.北京空間飛行器總體設計部，北京 100094

火星是距離地球最近的行星，也是太陽系中被探測最頻繁的行星。以往的火星登陸任務中，失敗多發(fā)生在進入-下降-著陸(Entry-Descent-Landing，EDL)過程，所以這個過程常被稱為“死亡7分鐘”，是關系任務成敗的關鍵因素。近年來，隨著火星登陸成功率提升，對火星大氣環(huán)境的了解不斷增加，對EDL過程、尤其是進入階段的氣體動力學和氣動熱力學預測能力隨之不斷提高，各種預測模型和方法也得到了驗證。Wright等對火星進入氣動熱力學模型的發(fā)展和應用進行了綜述，指出與對流加熱和表面催化模型相比，高速進入引起的激波層高溫氣體輻射加熱對于未來的火星探索任務可能非常重要。Johnston等的研究也表明，輻射加熱是火星進入熱防護設計中不確定性的主要來源之一。隨后，火星科學實驗室(Mars Science Laboratory，MSL)成功登陸火星，其飛行重建數(shù)據(jù)證明了這一點，同時證明火星進入輻射加熱預測存在較大不確定性。

隨著物理模型、計算技術和硬件能力的進步，數(shù)值模擬已成為氣動輻射預測的重要和有效手段之一。氣動輻射模擬一般由4個模塊組成：熱化學非平衡流動、氣體粒子激發(fā)態(tài)布居、氣體光譜特性和輻射傳輸過程。目前，火星進入氣動輻射預測的不確定性主要來自于非平衡流動與光譜特性的物理與計算模型。

高溫非平衡流動模擬為輻射光譜計算提供氣體密度、多模態(tài)溫度和組分濃度等參數(shù)信息，其計算準確性直接決定光譜計算精度。Hollis和Prabhu以對流熱通量預測精度為對象，綜述了相關試驗研究以及同數(shù)值模擬的對比情況，評估了Park化學反應動力學模型，發(fā)現(xiàn)對流加熱受壁面催化的影響更大，對化學反應和兩溫度模型不敏感。Johnston等對氣動輻射進行了不確定性分析，主要關注化學反應模型和平動-振動松弛模型，發(fā)現(xiàn)在進入速度為6.3～7.7 km/s的范圍內，駐點輻射加熱關于上述模型參數(shù)的不確定性為50%～200%，主要由CO離解速率和CO重粒子激發(fā)速率造成?；诖朔治龊碗娀〖げü?Electric Arc Shock Tube，EAST)的試驗結果，Johnston和Brandis提出了新的火星大氣化學反應速率模型，CO的離解反應速率被大大提高，結果表明新模型與試驗值吻合較好，傳統(tǒng)的Park模型明顯高估輻射。

光譜特性計算為輻射傳輸提供發(fā)射系數(shù)和吸收系數(shù)等輻射特性。董士奎等以HITEMP為基礎計算了CO窄譜帶模型參數(shù)。Palmer和Cruden使用NEQAIR(Nonequilibrium Radiative Transport and Spectra Program)計算了中低溫條件的CO紅外輻射，采用簡化CDSD-4000光譜數(shù)據(jù)庫，結果驗證了該溫度范圍內的CO光譜輻射模型。Lemal等結合數(shù)值模擬和地面試驗預測了CO的非平衡輻射加熱，分析了平衡和非平衡條件下的光譜輻射，特別在平衡狀態(tài)下對比了HITEMP-2000和CDSD-4000等不同數(shù)據(jù)庫，發(fā)現(xiàn)使用CDSD的計算結果與試驗值符合更好。針對高速非平衡流動的CO反應氣體輻射特性尚有大量基于流動場景的模型與試驗對比研究需要開展。

精細化光譜輻射強度測量試驗是氣動輻射研究的重要組成部分，試驗獲得的基礎數(shù)據(jù)是支撐計算模型驗證與確認的關鍵。Cruden在EAST上開展了速度為6～8 km/s的發(fā)射光譜定量測試，和相對低速條件下CO振動紅外輻射測量，發(fā)現(xiàn)在低密度條件下，絕大多數(shù)試驗均未達到平衡態(tài)。CO真空紫外輻射貢獻很大，還有CN輻射，CO振動躍遷引起的中紅外輻射則較小。Takayanagi等在HVST(Hyper-Velocity Shock Tube)的模擬火星大氣環(huán)境中，通過兩臺光譜儀測量了真空紫外到近紅外的發(fā)射光譜。還開展了速度2.5～8 km/s范圍的CO振動中紅外輻射測量，發(fā)現(xiàn)激波速度小于6 km/s時，紅外輻射隨著速度降低而增加。

未來的火星探測進入器具有更大的尺寸和重量，需要采用新的低冗余熱防護系統(tǒng)設計，以實現(xiàn)降低發(fā)射成本、增加載荷和增強系統(tǒng)穩(wěn)定性與安全性的目標，那么，更為準確和可靠的氣動輻射預測就尤為重要。目前的火星進入氣動輻射計算模型尚存在較大的不確定性，用于支撐模型驗證與改進的地面試驗數(shù)據(jù)在多樣性和廣泛性方面也顯不足，因此，有必要利用不同地面設備開展更多的精細光譜測量試驗，加強試驗和數(shù)值研究的相互協(xié)作，完善和驗證計算模型，以此為基礎構建準確的火星進入氣動輻射預測手段。

1 控制方程與計算方法

1.1 非平衡流動

1.1.1 守恒方程

火星進入流動求解多組分Navier-Stokes方程，考慮熱化學非平衡模型。質量守恒方程為

(1)

動量守恒方程為

(2)

式中：為密度;為壓力；為應力張量分量；為Kronecker符號。

考慮高速流動下內能與平動能處于非平衡狀態(tài)，能量方程包括內能守恒和總能守恒。內能守恒方程為

(3)

式中：為質量平均內能；int,為內能梯度引起的方向能量通量；源項代表混合氣體中由粒子碰撞過程引起的內能松弛速率。平均內能由各組分內能得到：

(4)

總能量守恒方程為

(5)

式中：為能量梯度引起的方向總能量通量；為組分的焓;為總能量，其表達式為

(6)

其中：為組分的總能量，包括平動能、轉動能、振動能和電子能，本文采用雙溫度模型，即平動能與轉動能平衡，振動能與電子能平衡：

(7)

(8)

的表達式為

(9)

式中:為組分的壓力。

對于質量、動量和能量的輸運過程，有

(10)

式中：為組分的摩爾分數(shù)；組分的等效擴散系數(shù)、黏性系數(shù)、平動轉動熱傳導系數(shù)和振動電子熱傳導系數(shù)的計算公式與取值參見文獻[22]。

1.1.2 源 項

化學反應源項由有限速率化學反應模型計算。將化學反應寫為

(11)

式中：為化學反應數(shù)；,和,分別為反應中組分的化學計量系數(shù)。那么單位體積內組分的質量生成率源項為

(12)

式中:為組分的摩爾質量;f,和b,為反應的正向和逆向反應速率：

(13)

其中:f,和b,分別為反應的正向和逆向反應速率系數(shù)，由Arrhenius反應速率模型得到：

(14)

式中：f,、f,和b,、b,分別為反應的正向和逆向反應速率常數(shù)；f,、b,為正向和逆向反應的活化能;為玻爾茲曼常數(shù)。

對于火星大氣環(huán)境，采用Park的8組分(CO、CO、O、O、C、N、N、NO)、9反應模型，具體反應見表1，反應速率常數(shù)參考文獻[6-7]。

內能松弛源項對于電離氣體而言形式非常復雜，對于典型火星進入場景，離子化反應程度不強，可以不考慮離子和電子等造成的內能轉換，僅考慮由碰撞過程引起的振動能松弛，該源項一般采用Landau-Teller模型：

表1 火星大氣的Park反應模型Table 1 Park’s reaction model of Mars atmosphere

(15)

1.1.3 數(shù)值方法

采用基于多塊結構網格的有限體積法離散計算域，對流項使用van Leer格式，界面值由MUSCL (Monotonic Upwind Scheme for Conservation Laws)方法使用van Albda限制器計算；黏性項采用二階中心格式；時間推進采用LU-SGS (Lower Upper Symmetric Gauss Siedel scheme)方法；并行使用MPICH實現(xiàn)。

1.2 光譜輻射模型

氣體輻射包括原子束縛-束縛躍遷(原子線譜)、原子束縛-自由躍遷(光致電離)、原子自由-自由躍遷(軔致輻射)和分子束縛-束縛躍遷(分子帶譜)，并可能包括化學反應發(fā)光。火星進入氣體輻射的主要來源包括：分子帶譜，如CO振動躍遷，CO(4+，真空紫外)，CN(B-X，紫色)、C(Swan，藍紫色)和CN(A-X，紅色)，以及原子輻射，如C和O(紫外-可見-紅外)。

對于電子躍遷，采用窄帶法計算輻射特性：

(16)

對于振動-轉動躍遷，以CDSD數(shù)據(jù)庫為基礎進行計算。氣體吸收譜線之間會發(fā)生部分重疊，在波數(shù)處的光譜吸收系數(shù)等于相互重疊譜線在波數(shù)處的吸收系數(shù),之和：

(17)

式中：為標準化為單個分子的譜線積分強度，計算方法和數(shù)據(jù)庫參見文獻[11,13]；(-0,)為譜線線型函數(shù),0,為第條譜線的中心波數(shù)；為分子數(shù)密度。非平衡條件下應用流場模擬獲得的平動轉動溫度和振動電子溫度計算配分函數(shù)等參數(shù)，具體參考文獻[7,26]。

1.3 輻射傳輸方程

存在吸收、發(fā)射和散射的非灰氣體輻射傳輸方程為

(18)

式中：為輻射強度；為空間位置；為立體角；為譜帶發(fā)射系數(shù)；為吸收系數(shù)；s為散射系數(shù)；(,′)為散射相函數(shù)，其中和′分別為入射角和散射角;下標表示光譜離散的譜帶。散射項對計算資源需求非常大，所幸進入氣動環(huán)境中高溫氣體的散射非常弱，可忽略，這樣極大簡化了方程求解。劉林華等通過離散坐標法進行了三維空腔的輻射計算，牛青林等采用視在光線法對滑翔飛行器紅外輻射特性進行了模擬。本文基于進入器輻射熱流預測需求，選擇光線追跡法和有限體積法解算輻射傳輸方程，以獲得特定方向輻射強度和飛行器表面輻射熱流分布。

光線追跡法針對一條光線，能夠得到流場內沿該傳播方向的詳細光譜結構和數(shù)據(jù)，適用于光譜模型的驗證。在建立沿一定方向的光線和離散點后，如圖1(a)所示，通過流場插值獲得節(jié)點處的密度、溫度和濃度等信息，進而得到發(fā)射和吸收系數(shù)等，由節(jié)點開始計算輻射傳輸，至節(jié)點結束，即得到節(jié)點的光譜輻射強度。具體的，根據(jù)Beer定律：

(19)

式中：為頻率；為光學厚度：

(20)

有限體積法則能夠直接獲取進入器表面所有位置與空間的輻射通量，并直接利用流場網格與信息，避免插值引起的誤差，借助CFD求解框架實現(xiàn)并行計算，適用于進入器熱環(huán)境評估的輻射強度和輻射熱流計算。圖1(b)顯示了有限體積法在位置空間和角度空間離散輻射傳輸方程得到的控制體，為內部格心節(jié)點，、、、、和分別為表面積分點，角度離散將4π空間分為個立體角，向量為立體角中心矢量。通過簡化和數(shù)值變換，可得不計散射的輻射傳輸離散方程：

(21)

圖1 光線追跡法和有限體積法離散Fig.1 Ray tracing and finite volume discretization

2 CO2輻射計算與驗證

2.1 光譜輻射強度

采用激波管試驗的測量結果驗證中高溫條件CO光譜輻射計算。試驗分別選自NASA Ames研究中心于2014年和JAXA于2018年開展的地面試驗。數(shù)值模擬使用一個長2 m的激波管為對象，其中充滿純CO氣體，利用平衡反應產物作為激波管虛擬驅動氣體，進行一維非定常流動模擬，劃分1 000個網格點，時間推進采用3階Runge-Kutta方法。

圖2是被驅段壓力為133 Pa、激波速度為3.0 km/s 激波管試驗狀態(tài)的模擬結果，分別為壓力、溫度和各組分質量分數(shù)。波后非平衡區(qū)比較窄，平動溫度和振動溫度迅速平衡，平衡溫度約3 500 K。此條件下，約10%的CO離解，產物主要為CO和少量的O與O。

圖3為激波到達1.55 m位置時，波后10 cm處光譜發(fā)射強度的文獻測量值和數(shù)值預測結果的對比。結果表明，激波速度為3 km/s時，計算得到的光譜輻射結構與文獻相符，輻射主要集中于CO的4.7 μm波段，輻射峰值和譜帶寬度均與測量結果符合較好，驗證了流動及輻射特性計算模型。

圖2 被驅壓力為133 Pa、激波速度為3 km/s的激波管流動參數(shù)Fig.2 Shock tube parameters at 133 Pa driven pressure and 3 km/s shock velocity

圖3 CO2光譜發(fā)射強度的計算值與參考對比Fig.3 Comparison of CO2 spectral emission intensity between calculation and reference

2.2 相對輻射強度測量與驗證

以CO為試驗介質，在典型進入速度范圍內，基于燃燒驅動激波管開展了一系列氣體發(fā)射強度測量試驗。相關試驗設備和測量方法詳見文獻[32-33]。被驅動段充入純CO氣體，壓力設為300 Pa，通過調節(jié)驅動段氣體壓力以產生所需的激波速度，測得的激波速度范圍為3～6.6 km/s。采用濾光片和光電探測器測量氣體發(fā)射強度，測量位置均位于波后平衡區(qū)內。

圖4為以4 km/s激波速度的發(fā)射強度為參考值得到的歸一化發(fā)射強度隨激波速度的變化，包括計算值、試驗測量結果和Cruden等的試驗結果。3套數(shù)據(jù)之間符合較好，進一步驗證了火星進入非平衡流動和輻射特性計算模型在不同狀態(tài)下的準確性。同時，結果也表明，當進入速度小于6 km/s時，發(fā)射強度隨速度降低而增加，與地球再入的輻射變化規(guī)律截然不同，此結論與文獻[19]一致。

圖4 不同激波速度下計算、試驗和文獻的發(fā)射強度E對比Fig.4 Comparison of emission intensity E among computations, tests and reference data under different velocities of shock

3 進入器輻射熱流

3.1 模型與計算條件

選擇火星探路者號(Mars Pathfinder，MPF)為進入器模型進行數(shù)值分析。MPF的前體為直徑2.65 m的70°鈍球錐，后體為單錐，外形和網格如圖5所示。0°側滑下流場關于=0 m平面對稱，因此網格僅考慮半模。周向采用C型拓撲，流向采用C-H網格，表面駐點區(qū)為非奇性H型網格。根據(jù)初算結果，對網格在邊界層和激波附近進行加密，以確保網格正交性、波系平行和分辨率。計算條件見表2，表中為高度，為來流速度，為來流密度，為來流溫度，為攻角。

圖5 MPF外形和網格Fig.5 Configuration and mesh demonstration of MPF

表2 MPF計算條件Table 2 Computational conditions for MPF

3.2 網格收斂性

針對MPF前體外形，生成了3套網格，具有不同的壁面法向最小網格長度和法向網格數(shù)目，法向網格最小距離()分別為10m、10m和10m，法向網格數(shù)目分別為121、121和201。圖6顯示了數(shù)值計算得到的表面對流熱流()沿徑向的變化，圖7是輻射強度沿駐點線的變化?？梢?，3套網格的計算結果相差不大，最小距離10m和10m網格的計算結果非常接近，綜合考慮計算精度和速度，以下分析均使用10m網格進行計算。

圖6 不同網格的表面對流熱流Fig.6 Surface convective heating flow for different grids

圖7 不同網格的駐點線輻射強度Fig.7 Radiation intensity along stagnation line for different grids

3.3 流動與輻射加熱

圖8和圖9分別是兩種進入速度的對稱面流動參量分布，包括溫度和CO質量分數(shù)。不同的激波強度導致激波層內氣體溫度和組分濃度差異明顯。高速條件下，激波層氣體溫度超過7 000 K，在這種高溫條件下化學反應劇烈，CO完全離解，進入器周圍主要為CO、O和O。尾流中氣體由于快速膨脹，溫度迅速降低至約2 000 K，此處氣體主要為由上游流至，由于流動速度很快，復合反應沒有足夠時間展開，因此，CO仍維持完全離解狀態(tài)。在尾流中心處，旋渦結構及其上方的剪切結構導致氣體溫度升高至4 000 K左右，此處氣體較少來自上游，而當?shù)貧怏w由于缺少激波壓縮，密度較低，化學反應程度較弱，因此CO離解程度相對稍低。低速條件下，激波層厚度明顯增加，尾流呈現(xiàn)更為明顯的剪切和反射激波，激波層和尾流核心區(qū)內的氣體溫度都在1 500 K左右，雖然低空的氣流密度更高，但溫度較低導致CO基本沒有離解，進入器周圍主要為CO。

圖8 MPF對稱面溫度分布Fig.8 Temperature distribution in symmetry plane of MPF

圖9 MPF對稱面CO2質量分數(shù)分布Fig.9 CO2 mass fraction distribution in symmetry plane of MPF

考慮沿駐點線傳播的光線，建立光線追跡網格，獲取相關流場參數(shù)，如圖10(a)所示，分別為2個狀態(tài)對應的溫度和組分質量分數(shù)沿駐點線的變化。對該光線進行輻射特性和輻射傳輸計算，得到如圖10(b)所示的輻射強度變化，可見隨著光線進入激波層并向壁面?zhèn)鞑?，高溫氣體不斷發(fā)射輻射能量，輻射強度逐漸增大，最終傳輸至邊界層，考慮到邊界層內氣體溫度較低，對輻射傳輸?shù)挠绊戄^小。

圖10 流動參量和輻射強度沿駐點線變化Fig.10 Variations of flow parameters and radiation intensity along stagnation line

圖11對比了2個速度狀態(tài)的駐點總光譜輻射強度，結果顯示速度為2.1 km/s時，雖然氣體溫度遠低于高速條件，但紅外段的輻射強度遠大于高速條件。速度為6.6 km/s時，短波光譜有明顯升高，對輻射產生更大的影響。

圖11 不同速度的駐點光譜輻射強度Fig.11 Spectral radiation intensity at stagnation point for different velocities

圖12為速度6.6 km/s狀態(tài)的駐點光譜輻射強度隨波長的變化，由于對CO(4+)的真空紫外輻射尚未完成驗證，因此譜帶并未包括真空紫外段。輻射能量主要集中在波長2.5～3.5 μm和4～6 μm的紅外譜段，由CO和CO疊加形成，O、N等原子和O、N等分子的電子躍遷譜線集中于波長2 μm以下的短波區(qū)域。

圖12 速度6.6 km/s狀態(tài)駐點光譜輻射強度Fig.12 Spectral radiation intensity at stagnation point under 6.6 km/s

圖13為速度2.1 km/s時CO和CO產生的駐點光譜輻射強度。可見低速時，由于流場內組分基本是中等溫度的CO，故輻射基本完全由CO產生，其強度遠高于其他組分。對比圖12，高速時高溫及離解產生的CO成為重要的輻射來源，因此，對于火星進入，高速和低速條件產生的非平衡流動及其輻射機制是完全不同的。

圖13 速度2.1 km/s狀態(tài)駐點光譜輻射強度Fig.13 Spectral radiation intensity at stagnation point under 2.1 km/s

基于流場計算結果，在所有網格點進行發(fā)射系數(shù)和吸收系數(shù)的計算，根據(jù)Planck平均公式得到0.2～10 μm光譜范圍內的平均吸收系數(shù)，對稱面分布如圖14所示。吸收系數(shù)分布和流場結構高度相似，發(fā)射較強的區(qū)域一般為高密度且高溫度區(qū)域，例如，高速條件下，吸收系數(shù)較大的區(qū)域為激波層和尾流剪切層，核心區(qū)密度較低，因此吸收系數(shù)也較低。低速條件下激波層內的氣體吸收系數(shù)明顯高于其他區(qū)域。

圖15分別為進入速度6.6 km/s和2.1 km/s的前體和后體表面輻射熱流()分布。高速時，前體激波層較薄，雖然氣體溫度很高，但CO完全離解，輻射熱流較低，峰值熱流出現(xiàn)在后體鄰近肩部處，主要來自尾流中大范圍的中高溫氣體。同時，由于球錐外形的駐點區(qū)激波層較錐體更薄，導致駐點輻射熱流低于錐體。低速時，CO未離解，氣體輻射主要來自CO的貢獻。激波層和尾流的溫度接近，但尾流密度明顯小于激波層，導致前體輻射熱流高于后體。同樣的，雖然駐點區(qū)與錐體的激波脫體距離差異不如高速時明顯，但駐點區(qū)輻射熱流仍小于錐體。另外，駐點區(qū)輻射熱流在低速條件下高于高速條件，這與激波管試驗的結論一致。

圖14 對稱面平均吸收系數(shù)分布Fig.14 Average absorption coefficient distribution in symmetric plane

圖15 前體和后體表面輻射熱流Fig.15 Radiative heating rate on fore and aft surface

4 結 論

對典型火星進入條件的氣體光譜輻射結構和輻射加熱開展了數(shù)值和試驗研究，得到以下結論：

1) 建立了火星大氣氣體光譜特性模型，計算結果與文獻數(shù)據(jù)在光譜結構和輻射強度上均符合較好，驗證了數(shù)值模型，CDSD光譜數(shù)據(jù)庫適用于火星進入光譜計算。

2) 通過激波管地面試驗，獲得了不同激波速度的輻射強度測量結果，與計算值和文獻數(shù)據(jù)符合良好，結果表明火星大氣進入速度低于6 km/s時，氣體輻射隨進入速度減小而增大，和地球再入規(guī)律不同。

3) 針對MPF典型進入狀態(tài)，開展了高速和低速進入的非平衡流動、輻射特性和輻射傳輸模擬。發(fā)現(xiàn)駐點區(qū)輻射熱流小于錐體，低速進入的輻射主要由CO形成，高速進入的輻射由CO和CO共同主導，低速進入的駐點輻射熱流高于高速條件的，高速時輻射峰值熱流出現(xiàn)在后體的肩部附近，低速時則出現(xiàn)在前體的錐體，以上規(guī)律可輔助火星進入器的熱防護系統(tǒng)設計。