土卫六探测器高超声速化学非平衡特性数值计算

徐一航, 陈少松, 周 航

(南京理工大学能源与动力工程学院, 江苏南京 210094)

引 言

土卫六是土星最大的卫星, 也是太阳系中唯一有大气的卫星.土卫六的半径为2 575 km, 比半径为2 430 km的水星大, 略小于半径为3 400 km的火星.土卫六的大气质量更大, 表面压力也比地球大.土卫六大气中的有机化学可能类似于地球上生命开始之前的原始大气, 因此发射探测器对研究土卫六大气有着重要的科学意义.

1943年Kuiper探测到甲烷时, 第一次在土卫六上探测到大气[1].其大气主要由分子氮、甲烷和氢组成, 而且可能含有甲烷云[2].1997年10月, “Cassini”号星际探测器被发射到飞往土星的轨道, 经过3.5×109km的漫长太空旅行之后, “Cassini”号成功进入土星轨道, 并向土卫六发射“Huygens”号探测器[3-6].Olejniczak等[7]对以12°攻角飞行的土卫六探测器进行了数值模拟, 得到了不同速度下探测器的阻力系数、升力系数和表面热流密度的变化规律.Bouilly[8]对“Huygens”号探测器在最后测试阶段出现的问题进行了总结.发现与设计之初相比, 飞行器表面热通量增加显著, 对热防护系统提出了更高的要求.Brandis等[9]通过实验证实了土卫六探测器激波层内形成的氰化物是一种高辐射物质, 忽略辐射耦合引起的非绝热效应会导致对表面温度的过度预测, 其峰值约为不考虑辐射效应的两倍.West等[10]对土卫六探测器表面激波层辐射加热预测的不确定性进行了研究, 结果表明辐射加热不确定性的最大来源是激波层内的化学非平衡流场.

图1是土卫六探测器进入大气过程的示意图.由于土卫六的大气环境与地球相似, 因此常采用在地球上进行再入飞行实验来检验探测器的综合性能[8].由于地球的大气成分和土卫六的大气成分存在差异, 可能会造成探测器表面激波层内的化学非平衡流场、热流密度的不同.因此, 探测器在地球大气和土卫六大气中飞行时气动特性的区别须进一步明确.采用数值模拟的方法分别对地球大气成分和土卫六大气成分下探测器前方流场内气动热环境进行了研究, 分析其表面热流密度存在的差异, 对比不同攻角下土卫六探测器前方热流密度的不同.

图1 土卫六探测器进入过程示意图

1 数值计算方法

1.1 流动控制方程

对于高温非平衡气体, 通常采用两温度模型进行数值模拟, 控制方程为三维热化学非平衡Navier-Stokes方程

其中, 对于两温度模型

Q=(ρi,ρEv,ρ,ρu,ρv,ρw,ρE)T

W=(wi,wv,0,0,0,0,0)T

ρi是组分i的密度,u，v，w是直角坐标系下3个方向的分速度,E为总能,Ev为分子组分的总振动能,F，G，H，Fv，Gv，Hv为对流项和黏性项,W为非平衡源项.

数值计算采用有限体积法对流动控制方程进行离散.采用 Roe′s FDS迎风差分格式对无黏通量进行离散, 限制器选用2阶精度的min mod限制器, 时间推进采用Euler隐式格式, 采用非催化壁面.

1.2 热力学模型

对于两温度热力学模型, 假设分子的平转动能和振动动能并不平衡, 其表述形式用平转动温度T和振动温度Tv来描述.混合气体总能为

气体的内能e为

e=etr+ero+eve+eo

气体的平动能和转动能分别表示为

其中，eo为混合气体零点能，ero和etr分别为混合气体平动能和转动能，R为摩尔气体常数，Ttr和Tro分别为平动温度和转动温度，Ms为组元s的摩尔质量.

1.3 化学反应模型

正向化学反应速率采用Arrhenius公式[11]

其中,Af，Bf和Cf分别为正向化学反应频率因子、温度指数和活化温度.

逆向化学反应速率通过化学反应平衡常数确定

Keq=exp(A0+A1z+A2z2+A3z3+A4z4)

其中,z=104/T,A0，A1，A2，A3，A4为常数.

对于探测器在地球大气中飞行, 采用5组元化学反应对探测器进行热化学非平衡数值计算, 其反应方程如表1所示.

对于探测器在土卫六大气中飞行, 其化学反应模型要比地球大气中更为复杂, 目前常见的是由Gokcen等[13]提出的N2-CH4-Ar模型, 该模型包含35个化学反应方程式, 这并不利于工程计算, 因此有必要对该反应模型进行简化.对于空间内参数数量较少的、以及不必要的中间反应过程进行简化.但是对于已知的强辐射物质(如: CN, C2), 无论它们是否是显著量, 都应当考虑.其次, 有些物质并不大量存在, 但是它为生成强辐射物质提供了重要反应途径(如: HCN), 因此也应当考虑[13].表2为简化后的土卫六大气化学反应模型.

1.4 数值模型验证

为了验证本文数值计算的合理性, 采用5组元化学反应模型对RAM-C球锥飞行器进行数值计算, 并将计算结果与实验值[12]进行了对比, 如图2所示.其计算结果与实验结果吻合, 因此采用本文的数值计算方法计算高超声速飞行器化学非平衡特性具有一定的可信度.

图2 壁面热流密度分布

2 数值计算分析

2.1 不同大气成分表面气动热环境分析

为了分析不同大气化学成分对探测器表面激波层内气动热环境的影响, 对典型外型的土卫六探测器进行了数值模拟, 模型与网格如图3, 4所示, 由于主要关注防热大底前激波层内的化学非平衡反应和表面热流密度, 因此采用图4的网格划分形式.计算条件如表3所示, 攻角为0°～20°; 采用非催化壁面, 壁温为300 K.地球大气成分的质量比例为: N2为74.7%, O2为25.3%.土卫六大气成分的质量比例为: N2为99.2%, CH4为0.8%[14].

表3 数值计算条件[15]

图3 模型示意图

图4 网格示意图

图5是不同速度下探测器在地球大气成分下和土卫六大气成分下表面热流密度的对比, 探测器在土卫六大气成分下表面的热流密度要明显高于其在地球大气成分下, 峰值可达1倍以上, 这对热防护提出了更高的要求.探测器在进入土卫六大气之后, 摩擦导致速度下降，但是由于大气逐渐稠密、压强增大, 探测器表面热流密度并没有随着速度的减小而减小，而是先随着高度的减小而增加, 随后当速度小于一定程度时热流密度开始减小.

图5 不同速度下探测器表面热流密度对比

图6是不同速度下探测器前方流场温度云图对比, 云图上方为地球大气成分, 下方为土卫六大气成分.当速度较高时, 土卫六大气成分下探测器的激波层厚度要比地球成分下的薄.并且土卫六大气成分下激波层内温度虽然没有地球成分下的温度高, 但是其激波层内的温度梯度小、整体温度偏高.图7是不同速度下头部轴线温度的对比, 地球大气成分下激波之后温度下降的速度较快; 而土卫六大气成分下激波之后温度几乎不变, 只有在接近壁面的时候突然降低, 这就导致了壁面的热流密度较大.图8是v=5 973 m/s时CH4组分云图, 激波之后温度较高, CH4解离程度达到80%.置换反应生成的强辐射物CN在脱体激波之前形成一道小的激波薄层(如图9所示), 这在一定程度上会降低壁面的热流密度[16].

图6 不同速度下探测器流场温度云图对比

图7 不同速度下头部轴线温度对比

图8 v=5 973 m/s时CH4组分云图

图10是CH4等质量分数沿头部轴线的分布曲线, 土卫六大气中的CH4解离得非常快, 导致激波后面碳、氢以及含有碳和氢原子的物质迅速增加.在地球大气中, 氧气的含量比较高, O-O的键能为498 kJ/mol, N≡N的键能为946 kJ/mol.在土卫六大气中虽然CH4的含量相对较低, 但是其有4个键能为414 kJ/mol的C-H.这就使得土卫六大气中气体解离所需要的能量比地球大气多, 在近壁面释放的能量也较多.地球大气中激波层内自由基反应生成物的化学键键能都相对较低(如N-O为607 kJ/mol), 释放能量相对较少.但土卫六大气中自由基反应生成物的化学键键能都相对较大(如C≡N为891 kJ/mol，C≡C为837 kJ/mol), 释放能量相对较多, 这导致了土卫六大气条件下探测器前方激波层的不同.

图10 v=5 973 m/s时CH4等质量分数沿头部轴线分布

2.2 不同攻角表面气动热环境分析

图11是不同攻角、土卫六大气成分下探测器前方流场温度云图, 随着攻角的增大, 脱体激波逐渐向上移动, 上方激波层内温度有所降低.图12是v=5 231 m/s时不同攻角下探测器表面热流密度的对比, 随着攻角的增大下表面热流密度增加, 上表面热流密度减小.因此, 当探测器以一定攻角进入大气时, 下表面的热防护十分重要.

图11 不同攻角下探测器流场温度云图对比

图12 v=5 231 m/s时不同攻角下探测器表面热流密度

3 结论

采用数值模拟方法分别对地球大气成分和土卫六大气成分下探测器前方流场内气动热环境进行了数值模拟, 得到以下结论:

(1)与地球大气条件相比, 土卫六大气条件下探测器表面的热流密度较高, 其峰值可达到地球大气条件的1倍以上.探测器表面的热流密度先随着高度的降低而增加, 而后当速度小于一定程度时热流密度开始减小.

(2)CH4在激波层内发生较大程度的解离反应, 置换反应生成的强辐射物CN在脱体激波之前形成一道小的激波薄层.

(3)随着攻角的增大下表面热流密度增加, 上表面热流密度减小.因此, 当探测器以一定攻角进入大气时, 下表面的热防护十分重要.