PHAROS求解火星进入热化学非平衡流场的测试及应用

杨星链，王京盈，郝佳傲，孙柯

1.山东大学，山东 济南 250061

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

火星登陆对于人类探测月球以远空间、开发太空资源甚至未来外星移民具有重要意义，已成为当前世界主要航天强国积极发展的热点。2021年2月18日，美国“毅力号”探测器在火星杰泽罗陨石坑安全着陆，美国国家航空航天局（NASA）再次成功执行了火星进入、下降和着陆（EDL）任务[1]。2021 年5 月15 日，我国“天问一号”探测器在火星乌托邦平原南部预选着陆区着陆，标志着我国首次火星探测任务着陆火星取得圆满成功，中国成为继美国之后第二个成功登陆火星的国家[2]。

探测器在EDL 过程中以极高速度进入火星大气层（约5km/s），经受着剧烈气动加热的严酷考验[3]，此时有效的热防护系统成为决定任务成败的关键。然而，EDL进入段流场复杂的高温热化学非平衡效应显著影响气动热环境的准确预测，给火星探测器的热防护设计带来很大的困难。因此，开展火星进入高温热化学非平衡流动模型和算法研究、计算流体力学（CFD）求解器的开发以及气动热环境精准预测分析，具有重要的科学意义和工程价值[4-5]。

在NASA 近60 年系列火星探测计划的推动下，美国的火星进入段流场CFD研究一直在稳步推进。G.Candler等[6]较早实现了相对完整的火星大气高超声速热化学非平衡流场数值求解。C.Park 等[7]给出了进入条件下高温火星大气的16 组分35 反应化学动力学模型，为之后火星进入段流场热化学非平衡效应的模拟奠定了重要的模型基础[8]。国内吕俊明等[9]采用三维Navier-Stokes方程求解程序结合5 组分8 反应模型研究了MSL 流场中化学非平衡效应的影响，并利用5 组分5 反应模型分析了火星大气参数对MSL 气动特性的影响[10]。杨肖峰等[11]基于自主研发的FL-CAPTER 软件平台，采用有效比热比方法模拟分析了MSL 升力-弹道式进入火星大气层的气动特性；而后，杨肖峰等又以5 组分6 反应模型研究了“探路者号”激波层流场的非平衡特性[12]，以及壁面催化条件对气动加热预测的影响[13]。刘庆宗等[14]利用自研的CFD 软件平台AEROPH_Flow 结合两温度和8 组分12 反应模型对球锥大钝头外形的火星进入段流场进行了求解，分析了表面材料催化特性对气动加热规律的影响。综上所述，国内针对火星进入热化学非平衡流场的CFD 模型、方法及软件研究起步较晚，既往诸多研究未考虑热力学非平衡效应，亦或未采用更新的化学反应动力学数据。

本文基于自研的高超声速气动热力学及热辐射优化并行求解器PHAROS[15-17]，采用Park 两温度模型考虑热力学非平衡[7]，使用Jaffe 等修正的反应速率系数[18]更新了原始的Park 火星大气组分化学反应动力学数据，能更加准确地考虑N2离解与CO 离解和置换反应，同时配合流场高精度算法，实现对火星进入热化学非平衡流场的求解。本文首先通过MSL 风洞试验算例对PHAROS 模拟能力进行可靠性测试，而后利用PHAROS对MSL真实飞行工况开展三维应用模拟研究。

1 控制方程与数值方法

1.1 控制方程

CFD求解器PHAROS仍然基于连续流假设进行研发。本文暂不考虑湍流和热辐射效应，采用两温度模型描述火星进入流场的热力学非平衡状态：分子转动模态完全激发且与重粒子平动模态平衡，对应于一个平动-转动温度Ttr，分子振动能以及电子平动能对应于一个振动-电子温度Tve。

流场中各组分密度、动量、能量以及振动-电子能量方程[19]在笛卡儿坐标系{x,y,z}下的守恒型形式如下

守恒变量U和源项Ω分别为

式中：ρi为组分i的密度；e和eve分别为混合物的总比能和比振动-电子能；ωi为化学反应源项；ωve为振动-电子能量方程源项。x方向对流通量和黏性通量具体表达式为

式中：指标mol.表示分子组分；p为混合物压强；Jx,i为组分i质量扩散通量的x方向分量；hs和eve,s分别表示组分s的比焓和比振动-电子能；qtr,x和qve,x分别为对应于平动-转动模态与振动-电子模态的热传导通量x方向分量。G和Gv以及H和Hv形式与F和Fv类似，分别为对应y和z方向的对流通量和黏性通量。

1.2 热化学非平衡模型

重粒子平动能和转动能以及电子平动能由能均分定理直接给出。谐振子假设下组分s比振动能有

式中：Rs表示组分s的气体常数；gr,s和θv,r,s为组分s振动模态r的简并度和振动特征温度。火星大气组分中，CO2为线性三原子分子，具有弯曲、对称及反对称三个振动模态[20]，其余分子组分均为单一振动模态，本文研究中振动能级数据取自参考文献[6]。

本文选用Micheltree 和Gnoffo 的8 组分（C、O、N、O2、N2、NO、CO、CO2）14 反应模型[20]计算控制方程的化学反应源项。该模型同样源自Park 等的16 组分35 反应化学动力学模型[7]，同时引入了Jaffe 等对N2和CO 离解反应速率系数的修正[18]，并已被D.Bose等[21]证实可实现对火星进入段气动热环境的准确预测。同时，PHAROS采用平动-转动温度Ttr与振动-电子温度Tve的加权平均作为化学反应的控制温度，以及化学反应与热力学非平衡过程的耦合效应。

控制方程中的振动-电子能方程源项ωve可进一步分解为

其中，平动-振动能量输运项ωt-v采用Landau−Teller形式[22]计算

式（4）中平动-振动松弛时间τv,s由Millikan−White关系式[23]计算。M.Camac[24]指出CO2组分的三个模态间松弛时间较短，可采用统一的平动-振动松弛时间计算式

源项(3)中的ωchem，v表示化学反应对振动能的影响，表达式如下

1.3 数值方法

采用有限体积法求解1.1 节中控制方程组。对流通量采用修正的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，分别考虑完全非催化和超催化壁面条件。计算网格量100×120，为获得准确的气动热数据保证壁面附近第一层网格雷诺数Recell小于5[29]，壁面附近第一层网格法向间距统一设置为5×10−7m（Recell=0.59），网格划分如图2所示。

图1 HET风洞试验MSL模型Fig.1 MSL testing model geometry in the HET wind tunnel

图2 MSL模型计算网格Fig.2 Computational grids for MSL testing model

图3 首先对比了PHAROS 给出的MSL 模型激波脱体距离和风洞试验数据，二者十分一致。图4 进一步对比了PHAROS 计算与风洞测量的模型表面热流，同时也比较了M.Sharma 等[28]的CFD 数据，本文还分别考虑了完全非催化和超催化两种壁面条件。注意图4 中在同一个位置有两个风洞试验数据，这是由于M.Sharma 等在试验中进行了重复测量[28]。由图4 可见，PHAROS 热流预测值与风洞数据仍然一致，与Sharma等的数值结果相近。超催化较完全非催化壁面条件的热流值更高，这是超催化条件下离解组分在壁面处再次复合放热所致。本节算例测试表明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]，计算总网格量约4.8 万，为保证热流预测精度壁面第一层网格法向间距1×10−7m（Recell=0.08），如图5所示。

图5 HYPULSE模型计算网格Fig.5 Computational grids for HYPULSE model

图6展示了该火星进入器模型流场的压强分布和流线结果，而图7给出了马赫数分布，可见模型头部附近激波层很薄，后体台阶处则出现了大尺度分离，且分离旋涡诱导出尾迹区内的二次压缩激波。图8 和图9 分别展示了流场的平动-转动温度Ttr和振动-电子温度Tve分布，Ttr总体高于Tve，激波后Ttr高达13000K，高温导致了显著的热化学非平衡效应。特别地，PHAROS 结果表明模型肩部诱导的膨胀区内Ttr和Tve也存在明显差异，反映了此区域经历着相当水平的热力学非平衡过程。

图6 HYPULSE模型压强分布（单位：Pa）Fig.6 Pressure distribution around HYPULSE model

图7 HYPULSE模型马赫数分布Fig.7 Mach number distribution around HYPULSE model

图10和图11分别给出了CO2和CO质量分数分布。高温激波层内CO2大量离解为CO，而低浓度CO2分布一直延续至下游尾迹区。值得注意的是，分离区内部CO 浓度有所降低，这是由于此区域温度相对较低，部分CO 重新复合为CO2所致。图12 进一步给出了驻点线温度分布，清晰展示了激波层内存在较大尺度的热力学非平衡区域，此时能量在平动-转动模态和振动-电子模态间发生交换。图13展示了模型表面热流，PHAROS 数值预测结果与风洞试验数据具有很好的一致性，其中模型驻点热流qw,stag量级接近10MW/m2。

图10 HYPULSE模型CO2质量分数分布Fig.10 CO2 mass fraction distribution around HYPULSE model

图11 HYPULSE模型CO质量分数分布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典型飞行工况三维流场模拟

本节使用PHAROS 对MSL 探测器真实飞行工况点的热化学非平衡流场进行三维模拟，飞行工况数据取自参考文献[31]，仅考虑MSL 前体部分，几何构型仍如图1 所示，但实际尺寸为Rb=2.25m。飞行参数为u∞=5660m/s（Ma=27.54），ρ∞=2.7×10−4kg/m3，T∞=157K，迎角α=15.7°。火星大气组分取质量分数为97%的CO2和3%的N2。壁面初始温度设置为Tw=300K，采用完全催化壁面条件，设置为辐射平衡壁，壁面发射率取ε=0.8。总网格量约95万，为保证热流预测精度壁面第一层网格法向间距2×10−5m（Recell=3.8），网格划分及分块如图14所示。

图14 MSL计算网格Fig.14 Computational grids for MSL

图15和图16分别给出了MSL流场的压强和马赫数分布，可见激波层十分薄，在有迎角飞行状态下，驻点位于MSL 前体下半表面。与纯球头型前体的表面压强围绕驻点形成单环式分布不同[17]，MSL 球锥形前体在有迎角飞行时上半表面亦会出现一个低压环状分布，与下半表面围绕驻点的高压环状分布对应，形成双环式压强分布。

图15 MSL压强分布Fig.15 Pressure distribution of MSL

图16 MSL马赫数分布Fig.16 Mach number distribution of MSL

图17和图18分别展示了MSL流场Ttr和Tve分布。激波后Ttr可突破8500K，Ttr与Tve分布差异明显，意味着显著的热力学非平衡效应；辐射平衡壁面导致的MSL前体表面高温分别出现在球头和下肩部，约1300K，同时表面的高低温差近300K。图19 给出了流场的CO2质量分数分布，可见CO2在高温激波层内大量分解。图20展示了MSL前体表面对流热流分布，峰值位于球锥钝头中心附近，为0.53MW/m2；在有迎角飞行状态下，MSL下肩部同样具有较高的气动加热水平。

图17 MSL平动-转动温度分布Fig.17 Translational-rotational temperature distribution of MSL

图18 MSL振动-电子温度分布Fig.18 Vibrational-electron-electronic temperature distribution of MSL

图19 MSL CO2质量分数分布Fig.19 CO2 mass fraction distribution of MSL

图20 MSL前体表面热流Fig.20 Wall heat flux of MSL forebody

3 结论

本文基于自研的有限体积CFD求解器PHAROS，引入两温度和8 组分火星大气化学反应动力学模型，成功实现了对火星进入高温热化学非平衡流场的数值模拟，并利用风洞试验和MSL 典型飞行工况算例对PHAROS 进行了可靠性测试和应用研究，研究表明，PHAROS有望为火星进入热化学非平衡流场求解和气动热环境预测研究提供一定的工具支持。本文研究结论如下：

（1）PHAROS 预测的激波脱体距离和表面热流与HET风洞MSL 模型试验数据一致，证实了PHAROS 求解器的可靠性。同时，超催化壁面比完全非催化壁面的热流值更高。

（2）HYPULSE 风洞试验模型高温激波层和肩部下游膨胀区内存在显著的热力学非平衡，HYPULSE风洞模型前体表面气动加热水平均在5MW/m2以上，驻点热流高达10MW/m2。

（3）三维MSL 有迎角飞行工况计算表明，MSL 辐射平衡壁面高温分布在球头和下肩部，约1300K，表面高低温差近300K；MSL 前体表面具有双环式压强分布，壁面热流峰值出现在球头附近，为0.53MW/m2，下肩部同样具有较高的气动加热水平。