PHAROS求解火星進(jìn)入熱化學(xué)非平衡流場(chǎng)的測(cè)試及應(yīng)用

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

1.山東大學(xué)，山東 濟(jì)南 250061

2.香港理工大學(xué)，香港 999077

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

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

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

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

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

1.1 控制方程

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

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

守恒變量U和源項(xiàng)Ω分別為

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

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

1.2 熱化學(xué)非平衡模型

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

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

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

控制方程中的振動(dòng)-電子能方程源項(xiàng)ωve可進(jìn)一步分解為

其中，平動(dòng)-振動(dòng)能量輸運(yùn)項(xiàng)ωt-v采用Landau?Teller形式[22]計(jì)算

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

源項(xiàng)(3)中的ωchem，v表示化學(xué)反應(yīng)對(duì)振動(dòng)能的影響，表達(dá)式如下

1.3 數(shù)值方法

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

2 PHAROS測(cè)試及應(yīng)用

2.1 HET風(fēng)洞MSL模型試驗(yàn)

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

圖1 HET風(fēng)洞試驗(yàn)MSL模型Fig.1 MSL testing model geometry in the HET wind tunnel

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

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

圖3 MSL模型壓強(qiáng)分布Fig.3 Pressure distribution around MSL model

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

2.2 HYPULSE風(fēng)洞試驗(yàn)?zāi)Ｐ?/h3>

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

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

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

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

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

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

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

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

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

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

2.3 MSL典型飛行工況三維流場(chǎng)模擬

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

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

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

圖15 MSL壓強(qiáng)分布Fig.15 Pressure distribution of MSL

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

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

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

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

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

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

3 結(jié)論

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

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

（2）HYPULSE 風(fēng)洞試驗(yàn)?zāi)Ｐ透邷丶げ▽雍图绮肯掠闻蛎泤^(qū)內(nèi)存在顯著的熱力學(xué)非平衡，HYPULSE風(fēng)洞模型前體表面氣動(dòng)加熱水平均在5MW/m2以上，駐點(diǎn)熱流高達(dá)10MW/m2。

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