PHAROS求解火星進入熱化學非平衡流場的測試及應用

楊星鏈，王京盈，郝佳傲，孫柯

1.山東大學，山東 濟南 250061

2.香港理工大學，香港 999077

火星登陸對于人類探測月球以遠空間、開發(fā)太空資源甚至未來外星移民具有重要意義，已成為當前世界主要航天強國積極發(fā)展的熱點。2021年2月18日，美國“毅力號”探測器在火星杰澤羅隕石坑安全著陸，美國國家航空航天局（NASA）再次成功執(zhí)行了火星進入、下降和著陸（EDL）任務[1]。2021 年5 月15 日，我國“天問一號”探測器在火星烏托邦平原南部預選著陸區(qū)著陸，標志著我國首次火星探測任務著陸火星取得圓滿成功，中國成為繼美國之后第二個成功登陸火星的國家[2]。

探測器在EDL 過程中以極高速度進入火星大氣層（約5km/s），經(jīng)受著劇烈氣動加熱的嚴酷考驗[3]，此時有效的熱防護系統(tǒng)成為決定任務成敗的關鍵。然而，EDL進入段流場復雜的高溫熱化學非平衡效應顯著影響氣動熱環(huán)境的準確預測，給火星探測器的熱防護設計帶來很大的困難。因此，開展火星進入高溫熱化學非平衡流動模型和算法研究、計算流體力學（CFD）求解器的開發(fā)以及氣動熱環(huán)境精準預測分析，具有重要的科學意義和工程價值[4-5]。

在NASA 近60 年系列火星探測計劃的推動下，美國的火星進入段流場CFD研究一直在穩(wěn)步推進。G.Candler等[6]較早實現(xiàn)了相對完整的火星大氣高超聲速熱化學非平衡流場數(shù)值求解。C.Park 等[7]給出了進入條件下高溫火星大氣的16 組分35 反應化學動力學模型，為之后火星進入段流場熱化學非平衡效應的模擬奠定了重要的模型基礎[8]。國內(nèi)呂俊明等[9]采用三維Navier-Stokes方程求解程序結(jié)合5 組分8 反應模型研究了MSL 流場中化學非平衡效應的影響，并利用5 組分5 反應模型分析了火星大氣參數(shù)對MSL 氣動特性的影響[10]。楊肖峰等[11]基于自主研發(fā)的FL-CAPTER 軟件平臺，采用有效比熱比方法模擬分析了MSL 升力-彈道式進入火星大氣層的氣動特性；而后，楊肖峰等又以5 組分6 反應模型研究了“探路者號”激波層流場的非平衡特性[12]，以及壁面催化條件對氣動加熱預測的影響[13]。劉慶宗等[14]利用自研的CFD 軟件平臺AEROPH_Flow 結(jié)合兩溫度和8 組分12 反應模型對球錐大鈍頭外形的火星進入段流場進行了求解，分析了表面材料催化特性對氣動加熱規(guī)律的影響。綜上所述，國內(nèi)針對火星進入熱化學非平衡流場的CFD 模型、方法及軟件研究起步較晚，既往諸多研究未考慮熱力學非平衡效應，亦或未采用更新的化學反應動力學數(shù)據(jù)。

本文基于自研的高超聲速氣動熱力學及熱輻射優(yōu)化并行求解器PHAROS[15-17]，采用Park 兩溫度模型考慮熱力學非平衡[7]，使用Jaffe 等修正的反應速率系數(shù)[18]更新了原始的Park 火星大氣組分化學反應動力學數(shù)據(jù)，能更加準確地考慮N2離解與CO 離解和置換反應，同時配合流場高精度算法，實現(xiàn)對火星進入熱化學非平衡流場的求解。本文首先通過MSL 風洞試驗算例對PHAROS 模擬能力進行可靠性測試，而后利用PHAROS對MSL真實飛行工況開展三維應用模擬研究。

1 控制方程與數(shù)值方法

1.1 控制方程

CFD求解器PHAROS仍然基于連續(xù)流假設進行研發(fā)。本文暫不考慮湍流和熱輻射效應，采用兩溫度模型描述火星進入流場的熱力學非平衡狀態(tài)：分子轉(zhuǎn)動模態(tài)完全激發(fā)且與重粒子平動模態(tài)平衡，對應于一個平動-轉(zhuǎn)動溫度Ttr，分子振動能以及電子平動能對應于一個振動-電子溫度Tve。

流場中各組分密度、動量、能量以及振動-電子能量方程[19]在笛卡兒坐標系{x,y,z}下的守恒型形式如下

守恒變量U和源項Ω分別為

式中：ρi為組分i的密度；e和eve分別為混合物的總比能和比振動-電子能；ωi為化學反應源項；ωve為振動-電子能量方程源項。x方向?qū)α魍亢宛ば酝烤唧w表達式為

式中：指標mol.表示分子組分；p為混合物壓強；Jx,i為組分i質(zhì)量擴散通量的x方向分量；hs和eve,s分別表示組分s的比焓和比振動-電子能；qtr,x和qve,x分別為對應于平動-轉(zhuǎn)動模態(tài)與振動-電子模態(tài)的熱傳導通量x方向分量。G和Gv以及H和Hv形式與F和Fv類似，分別為對應y和z方向的對流通量和黏性通量。

1.2 熱化學非平衡模型

重粒子平動能和轉(zhuǎn)動能以及電子平動能由能均分定理直接給出。諧振子假設下組分s比振動能有

式中：Rs表示組分s的氣體常數(shù)；gr,s和θv,r,s為組分s振動模態(tài)r的簡并度和振動特征溫度?；鹦谴髿饨M分中，CO2為線性三原子分子，具有彎曲、對稱及反對稱三個振動模態(tài)[20]，其余分子組分均為單一振動模態(tài)，本文研究中振動能級數(shù)據(jù)取自參考文獻[6]。

本文選用Micheltree 和Gnoffo 的8 組分（C、O、N、O2、N2、NO、CO、CO2）14 反應模型[20]計算控制方程的化學反應源項。該模型同樣源自Park 等的16 組分35 反應化學動力學模型[7]，同時引入了Jaffe 等對N2和CO 離解反應速率系數(shù)的修正[18]，并已被D.Bose等[21]證實可實現(xiàn)對火星進入段氣動熱環(huán)境的準確預測。同時，PHAROS采用平動-轉(zhuǎn)動溫度Ttr與振動-電子溫度Tve的加權平均作為化學反應的控制溫度，以及化學反應與熱力學非平衡過程的耦合效應。

控制方程中的振動-電子能方程源項ωve可進一步分解為

其中，平動-振動能量輸運項ωt-v采用Landau?Teller形式[22]計算

式（4）中平動-振動松弛時間τv,s由Millikan?White關系式[23]計算。M.Camac[24]指出CO2組分的三個模態(tài)間松弛時間較短，可采用統(tǒng)一的平動-振動松弛時間計算式

源項(3)中的ωchem，v表示化學反應對振動能的影響，表達式如下

1.3 數(shù)值方法

采用有限體積法求解1.1 節(jié)中控制方程組。對流通量采用修正的Steger?Warming 格式[25]計算，利用MUSCL 插值[26]提升至二階精度。黏性通量離散格式采用二階中心差分。時間推進使用線松弛迭代[27]。

2 PHAROS測試及應用

2.1 HET風洞MSL模型試驗

選取伊利諾伊大學的HET 風洞零迎角MSL 模型試驗[28]作為PHAROS的驗證測試算例。模型幾何示意如圖1所示，Rn=12.7mm，Rs=1.27mm，Rb=25.4mm，前體傾角α=70°，后體傾角αf=40°。風洞試驗來流條件u∞= 3058m/s（Ma=5.5），ρ∞=1.44×10?2kg/m3，T∞=1172K，來流為純CO2氣體，壁面溫度取Tw=300K，分別考慮完全非催化和超催化壁面條件。計算網(wǎng)格量100×120，為獲得準確的氣動熱數(shù)據(jù)保證壁面附近第一層網(wǎng)格雷諾數(shù)Recell小于5[29]，壁面附近第一層網(wǎng)格法向間距統(tǒng)一設置為5×10?7m（Recell=0.59），網(wǎng)格劃分如圖2所示。

圖1 HET風洞試驗MSL模型Fig.1 MSL testing model geometry in the HET wind tunnel

圖2 MSL模型計算網(wǎng)格Fig.2 Computational grids for MSL testing model

圖3 首先對比了PHAROS 給出的MSL 模型激波脫體距離和風洞試驗數(shù)據(jù)，二者十分一致。圖4 進一步對比了PHAROS 計算與風洞測量的模型表面熱流，同時也比較了M.Sharma 等[28]的CFD 數(shù)據(jù)，本文還分別考慮了完全非催化和超催化兩種壁面條件。注意圖4 中在同一個位置有兩個風洞試驗數(shù)據(jù)，這是由于M.Sharma 等在試驗中進行了重復測量[28]。由圖4 可見，PHAROS 熱流預測值與風洞數(shù)據(jù)仍然一致，與Sharma等的數(shù)值結(jié)果相近。超催化較完全非催化壁面條件的熱流值更高，這是超催化條件下離解組分在壁面處再次復合放熱所致。本節(jié)算例測試表明PHAROS預測能力準確可信。

圖3 MSL模型壓強分布Fig.3 Pressure distribution around MSL model

圖4 MSL模型表面熱流Fig.4 Wall heat flux of MSL model

2.2 HYPULSE風洞試驗模型

利用PHAROS 對B.R.Hollis 和J.N.Perkins[30]在NASA HYPULSE膨脹風洞開展的70°球錐火星進入器模型開展軸對稱流場模擬研究。試驗為純CO2來流，u∞=4772m/s（Ma=8.9），ρ∞=5.79×10?3kg/m3，T∞=1088K，壁面溫度取Tw=300K，采用完全非催化壁面條件。模型幾何尺寸見參考文獻[30]，計算總網(wǎng)格量約4.8 萬，為保證熱流預測精度壁面第一層網(wǎng)格法向間距1×10?7m（Recell=0.08），如圖5所示。

圖5 HYPULSE模型計算網(wǎng)格Fig.5 Computational grids for HYPULSE model

圖6展示了該火星進入器模型流場的壓強分布和流線結(jié)果，而圖7給出了馬赫數(shù)分布，可見模型頭部附近激波層很薄，后體臺階處則出現(xiàn)了大尺度分離，且分離旋渦誘導出尾跡區(qū)內(nèi)的二次壓縮激波。圖8 和圖9 分別展示了流場的平動-轉(zhuǎn)動溫度Ttr和振動-電子溫度Tve分布，Ttr總體高于Tve，激波后Ttr高達13000K，高溫導致了顯著的熱化學非平衡效應。特別地，PHAROS 結(jié)果表明模型肩部誘導的膨脹區(qū)內(nèi)Ttr和Tve也存在明顯差異，反映了此區(qū)域經(jīng)歷著相當水平的熱力學非平衡過程。

圖6 HYPULSE模型壓強分布（單位：Pa）Fig.6 Pressure distribution around HYPULSE model

圖7 HYPULSE模型馬赫數(shù)分布Fig.7 Mach number distribution around HYPULSE model

圖10和圖11分別給出了CO2和CO質(zhì)量分數(shù)分布。高溫激波層內(nèi)CO2大量離解為CO，而低濃度CO2分布一直延續(xù)至下游尾跡區(qū)。值得注意的是，分離區(qū)內(nèi)部CO 濃度有所降低，這是由于此區(qū)域溫度相對較低，部分CO 重新復合為CO2所致。圖12 進一步給出了駐點線溫度分布，清晰展示了激波層內(nèi)存在較大尺度的熱力學非平衡區(qū)域，此時能量在平動-轉(zhuǎn)動模態(tài)和振動-電子模態(tài)間發(fā)生交換。圖13展示了模型表面熱流，PHAROS 數(shù)值預測結(jié)果與風洞試驗數(shù)據(jù)具有很好的一致性，其中模型駐點熱流qw,stag量級接近10MW/m2。

圖10 HYPULSE模型CO2質(zhì)量分數(shù)分布Fig.10 CO2 mass fraction distribution around HYPULSE model

圖11 HYPULSE模型CO質(zhì)量分數(shù)分布Fig.11 CO mass fraction distribution around HYPULSE model

圖12 HYPULSE模型沿駐點線溫度變化Fig.12 Temperature variation along the stagnation line of HYPULSE model

圖13 HYPULSE模型表面熱流Fig.13 Wall heat flux of HYPULSE model

2.3 MSL典型飛行工況三維流場模擬

本節(jié)使用PHAROS 對MSL 探測器真實飛行工況點的熱化學非平衡流場進行三維模擬，飛行工況數(shù)據(jù)取自參考文獻[31]，僅考慮MSL 前體部分，幾何構型仍如圖1 所示，但實際尺寸為Rb=2.25m。飛行參數(shù)為u∞=5660m/s（Ma=27.54），ρ∞=2.7×10?4kg/m3，T∞=157K，迎角α=15.7°?；鹦谴髿饨M分取質(zhì)量分數(shù)為97%的CO2和3%的N2。壁面初始溫度設置為Tw=300K，采用完全催化壁面條件，設置為輻射平衡壁，壁面發(fā)射率取ε=0.8。總網(wǎng)格量約95萬，為保證熱流預測精度壁面第一層網(wǎng)格法向間距2×10?5m（Recell=3.8），網(wǎng)格劃分及分塊如圖14所示。

圖14 MSL計算網(wǎng)格Fig.14 Computational grids for MSL

圖15和圖16分別給出了MSL流場的壓強和馬赫數(shù)分布，可見激波層十分薄，在有迎角飛行狀態(tài)下，駐點位于MSL 前體下半表面。與純球頭型前體的表面壓強圍繞駐點形成單環(huán)式分布不同[17]，MSL 球錐形前體在有迎角飛行時上半表面亦會出現(xiàn)一個低壓環(huán)狀分布，與下半表面圍繞駐點的高壓環(huán)狀分布對應，形成雙環(huán)式壓強分布。

圖15 MSL壓強分布Fig.15 Pressure distribution of MSL

圖16 MSL馬赫數(shù)分布Fig.16 Mach number distribution of MSL

圖17和圖18分別展示了MSL流場Ttr和Tve分布。激波后Ttr可突破8500K，Ttr與Tve分布差異明顯，意味著顯著的熱力學非平衡效應；輻射平衡壁面導致的MSL前體表面高溫分別出現(xiàn)在球頭和下肩部，約1300K，同時表面的高低溫差近300K。圖19 給出了流場的CO2質(zhì)量分數(shù)分布，可見CO2在高溫激波層內(nèi)大量分解。圖20展示了MSL前體表面對流熱流分布，峰值位于球錐鈍頭中心附近，為0.53MW/m2；在有迎角飛行狀態(tài)下，MSL下肩部同樣具有較高的氣動加熱水平。

圖17 MSL平動-轉(zhuǎn)動溫度分布Fig.17 Translational-rotational temperature distribution of MSL

圖18 MSL振動-電子溫度分布Fig.18 Vibrational-electron-electronic temperature distribution of MSL

圖19 MSL CO2質(zhì)量分數(shù)分布Fig.19 CO2 mass fraction distribution of MSL

圖20 MSL前體表面熱流Fig.20 Wall heat flux of MSL forebody

3 結(jié)論

本文基于自研的有限體積CFD求解器PHAROS，引入兩溫度和8 組分火星大氣化學反應動力學模型，成功實現(xiàn)了對火星進入高溫熱化學非平衡流場的數(shù)值模擬，并利用風洞試驗和MSL 典型飛行工況算例對PHAROS 進行了可靠性測試和應用研究，研究表明，PHAROS有望為火星進入熱化學非平衡流場求解和氣動熱環(huán)境預測研究提供一定的工具支持。本文研究結(jié)論如下：

（1）PHAROS 預測的激波脫體距離和表面熱流與HET風洞MSL 模型試驗數(shù)據(jù)一致，證實了PHAROS 求解器的可靠性。同時，超催化壁面比完全非催化壁面的熱流值更高。

（2）HYPULSE 風洞試驗模型高溫激波層和肩部下游膨脹區(qū)內(nèi)存在顯著的熱力學非平衡，HYPULSE風洞模型前體表面氣動加熱水平均在5MW/m2以上，駐點熱流高達10MW/m2。

（3）三維MSL 有迎角飛行工況計算表明，MSL 輻射平衡壁面高溫分布在球頭和下肩部，約1300K，表面高低溫差近300K；MSL 前體表面具有雙環(huán)式壓強分布，壁面熱流峰值出現(xiàn)在球頭附近，為0.53MW/m2，下肩部同樣具有較高的氣動加熱水平。