火星环绕器羽流效应仿真研究

陈 杰，贺碧蛟，蔡国飙

（北京航空航天大学宇航学院，北京 100191）

火星环绕器羽流效应仿真研究

陈 杰，贺碧蛟∗，蔡国飙

（北京航空航天大学宇航学院，北京 100191）

火星探测器发动机羽流效应影响飞行任务的安全，采用差分求解NS方程与直接模拟蒙特卡洛耦合的方法对火星环绕器姿轨控发动机羽流力、热以及污染效应进行了仿真研究。结果表明：发动机羽流对火星环绕器气动力、热及污染效应影响主要在发动机储箱上，最大热流值可以达到1789 W／m2，需要注意；发动机羽流对展开太阳能电池翼的气动力和气动力矩相对于发动机推力和推力力矩来说，小3～4个量级，影响较小；已分离着陆器在距发动机出口300 m以内时，羽流对已分离着陆器的影响较大，气动热流在100 W／m2以上；在100 m以内时，气动热流在1000 W／m2以上，需要对敏感设备采取一定的防护措施。不同热适应系数气动热流仿真对比计算结果表明，随着热适应系数的增加，热流密度值呈线性增加。数值模拟结果可以为火星环绕器工程设计提供参考。

火星探测；火星环绕器；羽流效应；DSMC

1 引言

2020年我国将进行首次火星探测，火星环绕器作为火星探测系统的一部分，将要完成地火转移飞行和轨道修正、完成火星捕获、进入火星停泊轨道、与巡视器可靠分离、携带载荷完成环火遥感等科学探测任务。由于发动机羽流会对火星环绕器及其上敏感设备、太阳能帆板等产生气动力和力矩、气动热以及污染效应，对火星环绕器的正常工作造成影响。为了使火星环绕器能够顺利完成任务，需要对火星环绕器发动机羽流效应进行仿真研究。

直接模拟蒙特卡洛（DSMC）方法［1］是目前模拟发动机羽流效应较为成熟的方法。基于DSMC方法，北京航空航天大学的蔡国飙等人开发了PWS（Plume Work Station）软件［2-3］，验证算例表明软件具有较高的置信度［4］。PWS软件已成功应用在神舟飞船、嫦娥3号着陆器等航天器的羽流仿真研究中［5-6］。

本文采用差分求解NS（Navier-Stokes）方程与DSMC相耦合的方法［7］，使用PWS对不同工况下发动机羽流对环绕器表面及敏感元件的气动力、热和污染效应进行仿真分析，对比不同热适应系数对环绕器气动热流仿真结果的影响。对于火星环绕器与火星着陆巡视器分离的过程，由于已释放的着陆巡视器与环绕器的距离不断改变，着陆巡视器与羽流作用是一个复杂的非稳态过程。本文采用将羽流流场与已释放着陆巡视器解耦的办法对羽流对已释放着陆巡视器的影响进行估算。

2 仿真模型

PWS软件采用二级构造法来构建复杂仿真模型［2］。首先在不影响羽流效应仿真分析结论的前提下，对实际模型进行适当简化。然后对PWS软件定义的几种常用的平面（直线）和曲面（曲线）等基本元素，使用“与”和“或”运算最终组合成拥有复杂边界的仿真模型。图1为使用PWS软件建立的火星环绕器羽流仿真简化模型。

3 数值方法

考虑到发动机内流场以及羽流核心区流场具有较高的压强和密度，使用DSMC方法求解，计算量过大。本文采用差分求解NS方程与DSMC相耦合的方法，对火星环绕器发动机羽流效应进行仿真研究。

对于火星环绕器与火星着陆巡视器分离的过程，本文采用解耦的办法对羽流对已释放着陆巡视器的影响进行估算。即首先计算发动机羽流场，得到羽流场压强、密度、温度、速度分布，然后根据着陆巡视器与环绕器距离和与发动机推力夹角，确定火星探测器分离轨迹，并得到轨迹上羽流压强、密度、温度和速度参数，最后根据局地平衡假设及分子动力学理论，计算分离轨迹上的羽流质量流率、压强和热流密度，估算羽流效应。

1）差分求解NS方程

目前计算连续流域流动的最常用方法是数值求解NS方程的方法。矢量形式的NS方程如式（1）所示［8］，其中，Q为守恒变量，F、G为无粘通量，Fv、Gv为粘性通量，S为源项。

一般真实的流场都是以湍流状态存在的，因此，控制方程中的粘性项和热传导项中的系数需由层流和湍流共同确定［9］，如式（2）、式（3）：

式中，μ、μl和μt分别为粘性系数、层流和湍流粘性系数，Pr、Prl和Prt分别为普朗特数、层流和湍流普朗特数。层流粘性系数一般随温度而变化，可由Sutherland公式较为精确地给出［9］。在求解雷诺平均控制方程时，为了使其封闭，必须引入计算湍流粘性系数的湍流模型［8］，计算中采用SST k-ω两方程模型。

2）DSMC方法

由于DSMC计算域的流场密度和气体温度较低，计算时做以下假设［10］：（1）假设流场中分子间的碰撞为二体碰撞；（2）仅考虑分子的转动内能；（3）气体流动为定常流动；（4）假设分子为可变径硬球模型（VHS）。

对于分子与壁面的作用模型，目前应用较为广泛的是Maxwell反射模型。Maxwell反射模型采用热适应系数α来表示分子在与壁面作用的过程中，镜面反射与漫反射各占的比例［10］，计算公式如式（4）：

其中，qi和qr是来流和反射流的能流量，qw是反射温度为壁温时完全Maxwell漫反射下的能流量。在完全Maxwell漫反射时，α＝1，完全镜面反射时，α＝0。热适应系数表征了反射分子的温度在多大程度上已经“适应”了壁面的温度。

4 羽流效应计算

火星环绕器在不同阶段将执行不同的飞行任务，需要不同数量、不同推力的发动机组合工作。对每个发动机工况，首先使用差分求解NS方程的方法得到发动机内流场及部分外流场，再使用得到的发动机出口截面参数作为发动机二维自由羽流场DSMC计算的入口条件，然后从得到的二维流场计算结果中截取流场参数作为三维羽流场入口条件，对三维羽流场进行仿真，得到不同工况下发动机羽流对火星环绕器及敏感元件的气动力、热及污染效应。采用不同热适应系数对发动机储箱的气动热流密度进行了对比计算。

1）入口条件

对火星环绕器进行全流场计算时，计算域较大，发动机推力较大以及发动机数量较多。在不影响计算结果的情况下，为减小计算量，将发动机喷管向外延伸一个母线倾角为15°的圆锥，以此圆锥面上的流场参数作为三维DSMC计算入口条件。

图2为分别以发动机出口和发动机喷管向外延伸一个母线倾角为15°的圆锥为DSMC入口条件得到的羽流场压力对比云图。可以看出，两种入口条件得到的羽流场基本吻合，只有在离发动机较远处、压力较小的流场区域有些许偏差，且偏差在工程允许范围内。结果说明，将发动机喷管向外延伸一个母线倾角为15°的圆锥作为DSMC入口条件，对羽流流场计算结果影响不大，此做法较为合理。对120 N和25 N发动机羽流场分析可以得到类似的结论。

为了得到此圆锥面上的流场参数，首先差分求解NS方程，得到发动机内流场及部分外流场结果，然后以发动机出口截面参数作为入口条件用DSMC方法对二维发动机羽流场进行计算，最后在得到的二维羽流场中截取圆锥面母线上的流场参数，作为三维羽流场DSMC计算入口条件。

图3为从3000 N发动机出口开始的5 m×5 m的外流场压力分布云图，图中红线表示作为三维DSMC计算入口的圆锥面位置。圆锥面所在位置的粒子压强和密度与发动机出口处相比，已大大降低，达到了减小计算量的目的。

2）边界条件

计算域边界：当粒子经计算域边界运动出计算域时，按逃逸处理，即将粒子删除。

模型壁面：当粒子撞击到模型边界时，将粒子运动按Maxwell反射处理。模型壁面的温度设置为常温300 K（模型壁面温度的变化，对羽流仿真结果影响不大）［2］，发动机壁面温度设置为1000 K，热适应系数设定为1。在对比算例中，分别计算了热适应系数为0、0.25、0.5、0.75和1的对比算例。

入口条件边界：将入口条件圆锥面设置为吸收边界条件，即在入口条件边界内部不存在粒子，粒子运动穿过吸收边界时即被删除。

5 结果分析

5.1 差分求解NS方程仿真结果

首先对发动机内流场及部分外流场使用差分求解NS方程的方法进行仿真计算。

图4为3000 N发动机内流场及部分外流场流场温度分布云图。可以看出，温度在燃烧室内最高，喷管扩张段为特形面，气体充分扩张。高温高速喷流喷入真空环境中，羽流边界形成了抛物线形的膨胀波，喷管型面作用使喷流形成两个连在一起的锥形激波。3000 N发动机流场压强、密度、Ma数流场分布云图也呈现出类似的规律。

5.2 羽流气动力、热以及污染效应的影响

在差分求解NS结果的基础上，使用DSMC方法对发动机二维羽流场和三维羽流场进行仿真计算。对所有工况仿真结果分析总结，发现在3000 N+4×120 N+2×25 N发动机组合工作时，火星环绕器及其上设备所受到的羽流效应影响最大。

图5为该工况下的三维流场切面温度分布图。可以看出，在发动机储箱上方流场温度最高，这是由于3000 N与120 N发动机喷流相互阻滞的缘故，这与发动机储箱上的热流密度最大的仿真结果相一致。流场密度、压强分布图等也呈现出相似的规律。

图6为该工况下环绕器底板及其上设备羽流影响分布云图。可以看出，在该工况下，发动机羽流对发动机储箱的羽流效应影响最大，储箱上的最大热流值可以达到1789 W／m2，CO2质量到达率可以达到9.12×10-5kg／m2／s。对其它结果云图分析发现，羽流气动力压强、H2O质量到达率等分布呈现出类似的规律，各最大值均在发动机储箱上；而对于离3000 N主发动机距离相对较远的火星环绕器上的其它敏感元器件来说，其附近的羽流流场密度压强、温度、密度等变得相对较小，对敏感元器件的羽流气动力、热以及污染效应也因此变得较小。因此，在该火星环绕器设计方案下，主要考虑发动机羽流对发动机储箱所造成的气动力、热及污染效应。

5.3 太阳能电池翼上羽流气动力和力矩的影响

由于太阳能电池翼展开半径较大，需要考虑各工况下展开太阳能电池翼上的羽流气动力以及其对探测器中心点坐标（0，0，0.705）的力矩值。图7为0°、45°和90°的太阳翼所受到的发动机羽流气动力压强分布云图。

对不同角度的太阳能电池翼上的力和力矩的计算结果总结分析，得到不同角度（°）下的太阳能电池翼所受到的气动力在坐标轴上的分力值如表1所示。

表1 不同角度太阳能电池翼上羽流气动力Table 1 Aerodynamic force on solar wings with differrent angles

表2为工作的两台发动机A1和A2的各自推力在坐标轴上的分力值，以及两台发动机推力合力在坐标轴上的推力值。表3为不同角度下的太阳能电池翼所受到的羽流气动力相对于参考点坐标的力矩值。表4为A1、A2发动机推力相对于参考点坐标的力矩值以及它们的和推力相对于参考坐标的力矩值。

表2 发动机推力Table 2 Thrust of motors

表3 太阳能电池翼上羽流气动力相对于参考点的力矩Table 3 Torque of aerodynamic force on solar wings relative to the reference point

表4 发动机推力相对于参考点的力矩Table 4 Torque of motor thrust relative to the reference point

由表1和表2可以看出，发动机羽流对太阳能电池翼的羽流气动力在0.001 N量级，而发动机的推力则在10 N量级，即与发动机推力相比，太阳能电池翼上的羽流气动力影响可以忽略。

由表3和表4可以看出，发动机羽流对太阳能电池翼的羽流气动力相对于参考点（0，0，0.705）的力矩在0.01 N·m量级，而发动机推力相对于参考点的力矩值在10 N·m量级，与发动机推力矩相比，羽流气动力所造成的力矩影响也较小。

图8为太阳能电池翼上的气动力随着角度变化的变化趋势图。图9为太阳能电池翼上的气动力相对于参考点（0，0，0.705）的气动力矩随着角度变化的变化趋势图。

由图8可以看出，45°太阳能电池翼所受到的x轴向的气动力最大，即45°太阳翼受到的绕火星环绕器中心轴的扭转力最大。随着太阳能电池翼偏转角度的增加，太阳能电池翼上所受到的沿z轴的气动力逐渐增大，即90°时最大，0°时最小；沿y轴的气动力先减小后增大，在45°时最小，90°时最大。由图9可以看出，沿x轴方向的羽流气动力矩随着太阳能电池翼转角的增大而增大，这是因为随着转角的增大，太阳能电池翼受到的沿x轴方向的力越来越大，而力臂基本相同；沿y轴方向的羽流气动力矩在各角度下都基本为0；沿z轴方向的力矩随着太阳能电池翼转角的增大先减小后增大，在转角为45°时绝对值最大，方向与z轴正向相反；合力矩随着太阳能电池翼转角的增大而增大。

5.4 发动机羽流对已分离着陆器的影响

采用解耦的方法估算不同推力发动机羽流对距离发动机出口不同距离的已分离火星着陆器的影响。图10为距发动机出口不同距离下，3000 N发动机羽流场气动热流分布对比图。可以看出，在与发动机轴线夹角相同时，热流密度随着与发动机出口距离的增加而减小；而在距离发动机出口相同距离时，随着与发动机轴线夹角度数的增加，热流密度也呈下降趋势。仿真结果符合羽流场结构特点，即随着流场位置与发动机轴线夹角的不断增大，流场逐渐由核心流域过渡到稀薄流域，密度、压强、温度等流场参数值不断降低；而在与发动机轴线同一夹角下的位置上，喷流逐渐向真空膨胀，密度、压强、温度等流场参数值随着喷流的膨胀而降低。对其它结果云图分析发现，羽流场质量通量分布以及正应力分布也呈现出类似的规律。表5为3000 N发动机工作时，在已分离着陆器距发动机出口不同距离下的羽流对着陆器表面的最大质量通量、最大正应力和最大热流密度值。可以看出，当已分离着陆器距发动机出口300 m时，质量通量达到10-5kg／m2量级，正应力达到0.01 Pa量级，热流达到100 W／m2量级；当已分离着陆器距发动机出口100 m时，质量通量达到10-4kg／m2量级，正应力达到0.1 Pa量级，热流达到1000 W／m2量级，即对于分离着陆器来说，发动机羽流对已分离着陆器的影响在300 m内应该考虑，在100 m内应该对相关敏感元器件采取一定的防护措施。

表5 3000 N发动机各计算结果的最大值Table 5 Maximum value of plume effects of 3000 N motor

5.5 热适应系数对羽流气动热流仿真结果的影响

图11为不同热适应系数下发动机储箱上最大热流密度值的变化趋势图。可以看出，随着热适应系数的增加，发动机储箱上最大热流密度值基本呈线性增加。

当热适应系数为0时，发动机储箱上的热流密度为0，即此时粒子与发动机储箱表面没有发生能量交换，仿真结果符合理论实际。当热适应系数由0.25增大到1时，发动机储箱上热流密度分布基本相同，只是热流密度值随着热适应系数的增大而增大。发动机储箱上的热流密度最大值由热适应系数α＝0.25时的464 W／m2，增大到热适应系数α＝1时的1789 W／m2。

6 结论

对火星环绕器在入轨及器器分离过程中的羽流效应进行了仿真分析；对同一发动机工况下不同角度太阳能电池翼羽流气动热流进行了对比计算；对热适应系数对羽流气动热流仿真结果的影响进行了研究，结果表明：

1）受羽流气动力、热及污染影响较大的为发动机储箱，气动热流值可以达到1789 W／m2；羽流对环绕器表面及其他敏感元件的气动力、热及污染影响相对较小，气动热流值基本都在100 W／m2量级和10 W／m2量级。

2）发动机羽流对太阳能电池翼的气动力以及气动力相对于参考点的力矩数值较小，气动力在0.001 N量级；相对于参考点的气动力矩值在0.01 N·m量级，比发动机推力和发动机推力相对于参考点的力矩小3～4个量级；太阳能电池翼上羽流气动力和力矩相对于发动机合力和合力矩可以忽略。

3）对于分离着陆器来说，发动机羽流对距发动机出口100 m的已分离着陆器的质量通量达到10-4kg／m2量级，正应力达到0.1 Pa量级，热通量达到1000 W／m2量级；对距发动机出口300 m的已分离着陆器质量通量达到10-5kg／m2量级，正应力达到0.01 Pa量级，热通量达到100 W／m2量级，即对于分离着陆器来说，发动机羽流对已分离着陆器的影响在300 m内应该考虑，在100 m内应该对相关敏感元器件采取一定的防护措施。

4）当太阳能电池翼处于不同角度时，45°太阳能电池翼所受到的x轴向的气动力最大。随着太阳能电池翼偏转角度的增加，太阳能电池翼上所受到的沿z轴的气动力逐渐增大；沿y轴的气动力先减小后增大，在45°时最小。对于相对于参考点的气动力矩来说，沿x轴方向的羽流气动力矩随着太阳能电池翼转角的增大而增大；沿y轴方向的羽流气动力矩在各角度下都基本为0；沿z轴方向的力矩随着太阳能电池翼转角的增大先减小后增大，在转角为45°时绝对值最大，方向与z轴正向相反；合力矩随着太阳能电池翼转角的增大而增大。

5）热适应系数对羽流气动热流值的仿真结果的影响是线性的，即随着热适应系数的增大，羽流气动热流值仿真结果呈线性增加。

（References）

［1］ Bird G.Molecular Gas Dynamics［M］.Oxford：Clarendon Press，1976.

［2］ 蔡国飙，贺碧蛟.PWS软件应用于探月着陆器羽流效应数值模拟研究［J］.航天器环境工程，2010，27（01）：18-23.Cai G B，He B J.PWS software applied to numerical simulation on plume effect of lunar exploration lander［J］.Spacecraft Environment Engineering，2010，27（01）：18-23.（in Chinese）

［3］ 侯凤龙，蔡国飙.应用于DSMC方法的直角网格技术［J］.北京航空航天大学学报，2010，36（06）：694-699.Hou F L，Cai G B.Cartesian grid techmology applied in DSMC method［J］.Journal of Beijing University of Aeronautics and Astronautics，2010，36（06）：694-699.（in Chinese）

［4］ He B，He X，Zhang M，et al.Plume aerodynamic effects of cushion engine in lunar landing［J］.Chinese Journal of Aeronautics，2013，26（2）：269-278.

［5］ 蔡国飙，贺碧蛟.神舟飞船平移发动机对太阳能帆板羽流效应数值模拟研究［J］.载人航天，2011，17（04）：48-53.Cai G B，He B J.DSMC simulation for thee plume effect of the shenzhou spacecraft attitude control thrusters on the solar board［J］.Manned Spaceflight，2011，17（04）：48-53.（in Chinese）

［6］ 张熇，蔡国飙，许映乔，等.嫦娥三号着陆器软着陆过程中羽流仿真分析及试验研究［J］.中国科学：技术科学，2014，44（04）：344-352.Zhang H，Cai G B，Xu Y Q，et al.Simulation and experimental study of the plume during the Chang’E-3 lunar landing［J］.Scientia Sinica Technologica，2014，44（04）：344-352.（in Chinese）

［7］ 唐振宇，贺碧蛟，蔡国飙.解耦N-S／DSMC方法计算推力器真空羽流的边界条件研究［J］.推进技术，2014，35（07）：897-904.Tang Z Y，He B J，Cai G B.Investigation on boundary conditions in decoupled N-S／DSMC method for vacuum plume simulation of thrusters［J］.Journal of Propulsion Technology，2014，35（07）：897-904.（in Chinese）

［8］ 约翰D.安德森.计算流体力学基础及应用［M］.吴颂平，刘赵淼，译.北京：机械工业出版社，2015.JohnD.Anderson.Foundation and Application of Computational Fluid Dynamics［M］.Wu S P，Liu Z M，translate.Beijing：China Machine Press，2015.（in Chinese）

［9］ 李斌.超燃冲压发动机进气道的数值模拟研究［D］.湘潭：湘潭大学，2013.Li B.Numerical Simulation of Scramjet Inlet／Isolator［D］.Xiangtan：Xiangtan University，2013.（in Chinese）

［10］ 沈青.稀薄气体动力学［M］.北京：国防工业出版社，2003.Shen Q.Rarefied Gas Dynamics［M］.Beijing：National Defense Industry Press，2003.（in Chinese）

Simulation Study of Plume Effect on Mars Probe

CHEN Jie，HE Bijiao∗，CAI Guobiao
（School of Astronautics，Beihang University，Beijing 100191，China）

To mitigate the threat of the motor plume on the flight mission of the Mars probe，the hybrid NS equations／DSMC method was used to simulate the plume effect on the Mars probe.The simulation results showed that the plume of the motor mainly affected the storage box and the maximum heat flow reached 1789 W／m2.The effects of the aerodynamic force and the torque relative to the reference point on the expended solar wing were 3～4 orders’smaller than that of the motor.The plume effect on the separated Mars lander within 300 m was large，and the value of the aerodynamic heat flux was more than 100 W／m2.Within 100 m，the value of the aerodynamic heat flux was more than 1000 W／m2.So some protective measures for the sensitive equipment was necessary.The simulation with different thermal accommodation coefficients showed that with the increase of the thermal accommodation coefficient，the heat flow value increased linearly.The numerical simulation results can provide references for the engineering design of the Mars probs.

Mars exploration；Mars orbiters；plume effect；DSMC

V476.4

A

1674-5825（2017）06-0743-08

2017-04-17；

2017-10-02

陈杰，男，硕士研究生，研究方向为姿轨控发动机羽流效应仿真。E-mail：fyaycj＠163.com

∗通讯作者：贺碧蛟，男，博士，讲师，研究方向为航天器真空羽流效应。E-mail：hbj＠buaa.edu.cn

（责任编辑：庞迎春）