0 引言

飞船返回舱再入到130 km左右进行推返分离，需要在飞行高度约125 km到105 km期间根据预装的配平迎角进行调姿配平；当飞船返回舱跨流区飞行到中间大气层区域85~60 km附近，为了控制着陆地面落点精度，将根据预装的升阻比等参数实施升力控制。如果地面预测的高超声速再入段的配平迎角不准确，将会引起反作用控制系统(Reaction Control System，RCS)脉冲发动机的多次点火。这将会浪费过多的燃料，从而影响RCS发动机在跨声速流域的稳定控制^[1]。因此准确预测配平迎角随再入高度的变化对控制系统以及返回舱落点精度都是非常重要的。其它深空探测返回器与飞船返回舱有相似的钝体外形，以近第二宇宙速度、半弹道跳跃式再入大气层^[2-3]，由于是两次再入，高空稀薄段飞行时间显著增加，对返回器在稀薄流区域配平特性的准确预测对于落点控制至关重要。另外这类飞行器的再入速度超过7.5 km/s，再入飞行具有极高马赫数、低雷诺数特征，表面摩阻系数大幅度升高引发强黏性效应，而高超声速飞行绕流产生强烈的脱体激波，激波后的温度可达上万度，高温环境将导致气体分子发生结构上的变化，并会发生剧烈的化学反应。再入飞行器周围流场出现强烈的热、化学非平衡流动，又会对飞行器的力、热及配平特性产生影响^[4]。传统空气动力学研究手段预测得到的飞船再入配平迎角与实际飞行遥测结果偏差较大，如何建立可靠的模拟手段并正确分析评估返回舱从外层空间再入飞行过程特别是高空稀薄气体流域的配平特性是本文研究重点。

文献[5]从质心位置、壁面反射模型、马赫数和高温真实气体效应等多个方面对返回舱在稀薄流域的配平特性就行了初步研究，文献[6]则考察了烧蚀外形的变化对配平特性的影响。文献[7-10]分别对连续流域和稀薄流域高温真实气体效应对流场和气动特性的影响进行了数值模拟，缺乏高温和稀薄气体耦合效应下的跨流域流动机理分析。本文基于直接模拟蒙特卡罗(DSMC)方法，通过发展流场直角与表面三角形非结构混合网格及网格自适应技术，构造适于高稀薄流到近连续滑移流多流区共存的变时间步长模拟策略及DSMC区域分解并行算法，计算分析了返回舱再入跨流域(稀薄过渡流、近连续滑移流)激波与边界层的演变机理以及气动力系数变化规律，重点研究了高温真实气体效应对滑移流和稀薄过渡流区域气动力特性以及配平特性的影响，考察了壁面反射模型、质心纵横向位置变化对配平特性的影响规律。

1 数值模拟方法 1.1 计算网格及网格自适应

在DSMC方法中，流场中的网格是用来选取可能的碰撞分子对以及对宏观流动参数取样。流场采用均匀的直角坐标网格，追踪分子的效率非常高，计算区域内的模拟分子可以直接根据分子的位置坐标来确定分子所属的网格，而不必跟踪分子从一个网格运动到另一个网格。其缺点是无法精确地描述物面边界。本文采用直角与表面非结构混合网格的DSMC数值方法，在描述物面几何形状的非结构网格建立以后，直接将其嵌入到直角网格的流场中，使DSMC计算对流场网格的依赖程度大大降低。同时通过判断分子运动轨迹方程和物面三角形面元上任一点的位置方程，唯一确定出分子与物面的碰撞点坐标^[11]，解决了这种混合网格流场分子运动与物面碰撞的难题。对分子在物体三角形面元上碰撞、反射前后的流场参数进行统计取样就可以获得飞行器的整体气动力特性以及表面力、热载荷分布。

在流动梯度变化较大的区域，碰撞对应该通过有效的最相邻搜索方法选取。本文在高度非平衡区域采用基于密度梯度的动态自适应网格，避免在网格内随机选取的碰撞分子之间的距离远大于当地的分子平均自由程。即在背景网格的基础上，根据流场中密度梯度的变化分别对碰撞网格和取样网格进行细化，碰撞分子则是在自适应后最小的亚网格内选取，保证了计算的空间精度。

1.2 变时间步长格式

当计算区域中流场密度有较大的变化时(如强烈的激波压缩流动或激波干扰流动)，按照DSMC方法的模拟要求，时间步长要小于当地的分子平均碰撞时间，如果整个计算区域统一采用较小的时间步长，那么整个流场要达到稳定状态需要花费非常多的计算时间，计算效率很低。如果流场中每个模拟粒子的权重相同，则在高密度区域模拟的粒子数非常少，而在低密度区域粒子数又大大超过求解的需要^[12]。过多的时间浪费在低密度区域的计算上，而在高密度区域取样严重不足。由于每个网格内的粒子数与气体密度的平方呈反比，为了使计算区域每个网格的模拟分子数分布更加均匀，并且不会被分子“复制”产生有害的影响，必须将当地的时间步长和分子的权重相匹配。为了保证模拟粒子穿越网格交界面时通量守恒(包括质量、动量和动力学能量守恒)，W_i/Δt_i在整个计算区域都保持相同。模拟分子穿过网格交界面时的剩余时间，需要按照分子运动的初始网格和目的网格的时间步长的比值重新调整。经过时间步长的自适应过程后，网格i的时间步长可用下式计算：

$ \Delta {t_i} = \Delta {t_\infty }\left[ {\frac{{{n_\infty }}}{{{n_i}}}} \right] $

(1)

式中, Δt_∞和n_∞是时间步长和数密度的参考值，[·]表示截断。如果某个网格内的数密度大于参考数密度值，则其时间步长就等于参考时间步长。网格内的模拟粒子数仅根据其自身的权重自适应，而粒子权重又与时间步长成正比，因此在高密度区域粒子数会增加，密度低的区域粒子数则减少。这就有可能使自适应后的碰撞网格尺度接近当地的分子平均自由程。同时由于流场中总的模拟粒子数减少，也会大幅提高计算效率。

作者在文献[13]中针对返回舱外形对比分析了网格尺度、网格自适应和采用变时间步长对计算精度和计算效率的影响，其中的算例3结果表明采用本文给出的变时间步长技术可使流场中总的模拟粒子数减少至基准算例的1/4左右，计算所花时间也仅为基准算例的1/4，计算效率提高明显。

1.3 DSMC并行算法

DSMC并行算法采用区域分解的策略，根据计算的处理器数(或CPU核心数)的多少将计算区域划分为等量的子区域，每个处理器在其分配的子区域内部独立地计算模拟分子的碰撞和迁移，离开子区域的模拟分子把携带的信息传递给对应子区域的处理器。并行计算总的时间包括每个处理器计算碰撞、迁移的时间、处理器之间的通讯时间以及各处理器之间为同步而等待的时间。提高并行计算的效率主要是通过减少通讯和同步等待的时间来实现。

为了减少不同处理器计算时间的差别，通常采用负载平衡技术^[14-16]，本文采用静态随机负载平衡方法用于解决不同处理器之间的计算时间同步问题^[17]，该方法基于概率近似原理，将计算区域的全部网格平均分配给指定的所有处理器。由于采用相同的概率随机选取流场网格，按照均分后的数量分配给每个处理器，因此当计算区域的网格数量较多时，每个处理器包含近似相等的物体边界、高密度流动区域和稀薄气体区域的网格数，这样每个处理器的计算负载也非常接近。

2 计算结果分析 2.1 再入过程的跨流域流动机理分析

飞船返回舱再入过程跨越了自由分子流、稀薄过渡流、滑移流和连续流等多个流动区域，在稀薄过渡流和滑移流区域覆盖的飞行高度范围内，高超声速再入飞行环境下的黏性干扰、稀薄气体效应、高温真实气体效应相互耦合相互影响，采用DSMC数值方法研究了热化学非平衡流场结构和气动力系数随高度的演变规律。为了分析方便，本文计算中全部采用正迎角，这时的轴向力和法向力系数为正，抬头的俯仰力矩系数为正，升阻比为负。图 1是迎角20°、再入速度7.5 km/s条件下，将高度70~100 km的流场马赫数等值线分布画在一起，可以非常直观地看到随着高度的增加，激波层的厚度从高度70 km接近间断到高度100 km，激波演变成一个渐进的压缩过程。

从图 2、图 3不同高度的压力和密度沿流场中心轴线的变化可以看出，从激波开始压缩到激波压缩结束的两个平台之间，激波层由薄到厚的演变过程。图 3的密度变化同时反映出随着高度的增加，激波层和边界层都在不断增大，到100 km以上高空，激波层已经与边界层融合在一起，流动演变为渐进的压缩过程。

从激波层和边界层随高度(或稀薄度)的演变规律可以知道，随着高度增加，边界层增厚，黏性的影响逐渐加大，稀薄气体效应的影响逐渐明显。对于飞船返回舱这类大钝体飞行器，稀薄气体效应对法向力系数、俯仰力矩系数影响显著。图 4展示了迎角20°时返回舱轴向力、法向力和对质心俯仰力矩系数随高度的变化曲线。图中分别给出了压力和摩擦力的贡献，对于70~80 km高度而言，流动处于滑移流区域，稀薄气体效应影响微弱，以黏性干扰为主，90~130 km流动处于稀薄过渡流区，稀薄气体效应影响显著。而130 km以上高度则是自由分子流，气动力系数基本不再变化。对于轴向力系数来说，压力的贡献占主要部分，摩擦力的贡献非常微弱。而压力对法向力系数和俯仰力矩系数的贡献随高度几乎不变，摩擦力的贡献在90 km高度以上近似呈线性增长，这也是在过渡流区稀薄气体效应对返回器法向力系数、俯仰力矩系数影响显著的主要原因。

2.2 高温真实气体效应对气动特性影响分析

飞船返回舱或探月返回器都以极高的速度再入，飞行器周围的超高速气流通过激波的压缩和加热，会发生剧烈的化学反应。反应后的气体分子通过能量的再分配以及产生新的化学组元等改变了当地的流场结构，进而影响气动力系数以及配平特性。以返回舱为例，分别选取80 km和90 km两个高度来考察高温气体效应对流场和气动特性的影响。考虑到7.5 km/s的再入速度下80 km以上高度的气体密度相对稀薄，气体的电离较弱，因此化学反应仅考虑了五组元空气的离解和置换反应。在DSMC模拟中，TCE化学反应模型^[18]获得了最广泛的应用，非常可靠，本文亦在计算中采用该模型。作为对比，分别计算了考虑分子转动、振动激发和五组元空气化学反应(真实气体模型)^[9]以及将振动激发和化学反应冻结(完全气体模型)的返回舱流场及气动力系数随迎角的变化。图 5给出了迎角20°下80 km和90 km高度的流场压力等值线分布云图，可以看出在80 km高度头部激波的压缩非常强烈，激波层很薄。随着高度增加，激波压缩逐渐减弱，激波厚度也逐渐增加。图 6则是两个高度下的化学组元摩尔分数沿中心轴线的分布，图中显示出较强的化学非平衡流动特征。在激波层区域，分子的高密度产生高碰撞率，使离解活化能较低的O₂分子离解速率远大于N₂的离解速率，在到达物面前，O₂全部离解完毕。在五个化学组元中NO的含量是最低的，并且是先上升再下降接近于零，说明发生交换反应生成NO的几率并不高，同时也表明生成的NO在靠近物面前也基本离解完毕。随着高度的增加，化学反应的强度不断减弱，来流速度7.5 km/s时，在100 km以上几乎不再有化学反应的发生。

图 7分别给出了两个高度下两种气体模型计算的气动力系数随迎角的变化。对于90 km高度，高温真实气体效应对轴向力系数、法向力系数、俯仰力矩系数的影响相对较小，仅使配平迎角减小1.6°。而80 km高度处于滑移流，激波后高温气体的化学反应更加剧烈，两种气体模型的轴向力系数有明显差别，在0°迎角时的相对偏差最大达到5.2%，并且随着迎角的增大两者的偏差逐渐减小，但法向力系数两者相差较小，对配平迎角的影响较大，使配平迎角减小2.6°。同样高温真实气体效应使压心位置后移并随着高度增加影响减弱。

计算中还发现对于90 km以下高度在小迎角范围(≤15°)高温真实气体影响对轴向力系数的影响起主导作用，70~90 km高度小迎角范围的轴向力系数随高度的变化趋势与90 km以上高空稀薄气体效应占主导的轴向力系数变化规律正好相反，但在20°迎角下70~90 km的轴向力系数已非常接近。

2.3 壁面反射模型对配平特性影响分析

由于飞船返回舱在再入的稀薄气体流域，早期的地面试验和数值计算预测的配平迎角与实际飞行遥测结果有较大的差异，因此本文开展了不同壁面反射模型对120~90 km高度配平迎角和配平升阻比的影响规律研究。最经典的壁面反射模型是镜面反射和完全漫反射，这是两个极限状态。镜面反射时，气体分子碰到壁面后，只有法向速度改变方向，切向速度不变，因此在壁面没有黏性，计算的配平迎角最小。完全漫反射则是在壁面达到壁面温度下的Maxwell平衡分布，是壁面黏性影响最大的一种情况，而黏性增加会使返回舱的低头力矩增大，而要达到姿态配平则需要增大迎角以增加相应的抬头力矩，这样计算的配平迎角最大。在这两种壁面反射模型的基础上，又选取了Maxwell反射模型(镜面反射和漫反射的混合)以及CLL(Cercignani-Lampis-Lord)反射模型，通过改变调节系数，考察不同调节系数下的这两种壁面反射模型对返回舱配平特性的影响。对于Maxwell反射模型，调节系数分别取0.2和0.5，即20%的漫反射和80%的镜面反射以及漫反射和镜面反射各占50%两种情况，为了和Maxwell模型对应，CLL模型的法向动量调节系数和切向动量调节系数也都同时分别取0.2和0.5。为了避免过多的因素交叉影响，所有的反射模型计算中都不考虑气体分子的能量调节。

图 8给出了不同的壁面反射模型下计算的返回舱配平迎角随高度的变化。在镜面反射模型下，壁面没有黏性，返回舱配平迎角随高度变化不大，配平迎角从90 km和100 km高度的不到20°，在110 km和120 km配平迎角上升至接近22°和23°。完全漫反射模型下的配平迎角从25°不断上升至44°多，反映出漫反射模型的黏性影响使配平迎角增大较快。Maxwell反射模型和CLL反射模型在所取的调节系数下对90 km配平迎角的影响不到1°，在0.2的调节系数下，随高度增加，配平迎角增加的比较缓慢，从90 km到120 km仅增加了3°，与镜面反射时的增量近似。0.5的调节系数对90 km和100 km高度配平迎角的影响都比较小，与完全漫反射模型下的配平迎角非常接近，而在110 km和120 km配平迎角相比完全漫反射增加的相对较缓一些。因此90 km以下高空，壁面反射模型对返回舱配平迎角的影响非常微弱，到100 km及以上高空，调节系数的变化对配平迎角的影响明显加大。

图 9是不同壁面反射模型对应配平迎角下的配平升阻比随高度的变化曲线，尽管镜面反射模型、0.2调节系数下的Maxwell模型和CLL的配平迎角相差近5°，但三者的配平升阻比却比较接近，在100 km以上高度配平升阻比还有所增加。完全漫反射模型的配平升阻比随高度增大一直是下降的，0.5调节系数下的Maxwell模型和CLL模型在100 km以上高度变化比较小。

上述计算结果表明壁面反射模型对返回舱这类钝体外形90 km高度以上的配平迎角和配平升阻比影响较大，对于90 km以下高度影响迅速减弱。作为验证，本文选取了探月试验返回器第二次再入时100~80 km范围内的配平迎角进行计算并与飞行试验数据对比。图 10是计算和飞行试验数据的对比情况。计算中采用的壁面反射模型是Maxwell反射模型，首先将所有壁面全部采用0.8的调节系数(图中三角符号)，壁面温度为1500K等温壁。除了100 km高度由于飞行试验对配平迎角进行控制偏差较大外，只有95 km高度计算和飞行数据的偏差超过1°, 其它高度两者的偏差都在1°以内。通过分析不同部位的气动力对配平迎角的贡献发现，头部迎风锥面对配平迎角的影响较大，因此又计算了只是头部迎风锥面采用0.8的调节系数，其它部位全部采用完全漫反射(图中的菱形符号)，壁面反射模型的改变对90 km及以下高度配平迎角的影响微弱，与95 km高度飞行试验数据的偏差不到1°，并且使100 km高度的配平迎角也大幅降度。因此对于返回舱类钝体外形在90 km以上高空不宜采用完全漫反射模型计算气动力特性。

2.4 质心位置变化对配平迎角的影响分析

在返回舱的设计中，一般采用横偏质心位置的方法，来提供返回舱的再入配平迎角和实现飞行轨迹机动控制所需的配平升阻比。所谓飞行器的气动配平状态，就是指绕飞行器质心的俯仰力矩等于零。由于俯仰力矩系数本身是一个小量，又受外界因素的影响较大，对准确预测配平迎角造成一定的困难。飞船返回舱尤其是深空探测返回器由于调姿、控制中的燃料消耗以及再入过程中的烧蚀，总会引起质心位置的持续变化，因此需要关注质心位置变化对配平特性的影响。在“阿波罗”飞船研制过程中，美国开展了大量的地面试验和飞行试验来研究飞船质心位置的变化对配平特性的影响^[18-19]。本文对80~100 km高度范围飞船返回舱横向和纵向质心位置的变化对配平迎角的影响进行了计算分析。图 11给出了不同高度横向和纵向质心系数分别正向和负向变化1%时引起的配平迎角的变化量，可以看出无论增加或减少1%的横向质心系数，都会引起6°左右的配平迎角变化，并且随着飞行高度的增加，横向质心系数的改变对配平迎角的影响量有所减小。而对于纵向质心系数±1%的变化仅引起配平迎角约1°的改变，并且随着飞行高度的增加其变化量有所增大。因此横向质心位置的变化对配平迎角的影响较大，而纵向质心位置的变化对配平迎角的影响相对较小。

3 结论

本文建立了基于当地流动参数自适应的可考虑热化学非平衡流动的DSMC数值模拟方法，通过模拟返回器再入高度130~70 km范围跨越稀薄过渡流和近连续流的热化学非平衡流动，研究分析了跨流域激波层和边界层的演变规律以及稀薄气体效应、高温真实气体效应对气动力特性和配平特性的影响规律；对比分析了返回舱120~90 km高度范围不同壁面反射模型对配平特性的影响规律，计算的探月试验返回器的配平迎角与飞行试验数据一致；计算分析了高度100~80 km返回舱质心位置在横向和纵向的偏移对配平迎角的影响，有以下几点初步结论：

(1) 在滑移过渡流(高度70~80 km)，返回舱再入时的激波层和边界层都很薄，随着高度增加，激波脱体距离增大，激波层和边界层厚度相应增加，激波的间断效应减弱，到100 km以上高空激波层和边界层完全融合，激波演变为一个渐进的强扰动压缩过程。

(2) 高温真实气体效应使返回舱的轴向力系数增大、法向力系数和配平迎角减小、压心位置后移，并且随着高度增加影响迅速减弱。

(3) 压力对轴向力系数的贡献占主要部分，而俯仰力矩系数则以摩擦力的贡献为主；压力对法向力系数和俯仰力矩系数的贡献随高度变化微弱，摩擦力的贡献在90 km以上高度近似呈线性增长，这也是稀薄气体效应对法向力系数、俯仰力矩系数影响显著的主要原因。

(4) Maxwell反射模型和CLL反射模型在相同调节系数情况下，100 km以下高度对返回舱配平迎角的影响非常接近，改变调节系数Maxwell反射模型和CLL反射模型对90 km以下高度配平迎角的影响不超过1°。随着高度增加，反射模型对配平迎角的影响越来越大。

(5) 返回舱增加或减少1%的横向质心系数，都会引起约6°的配平迎角变化。随着飞行高度增加，横向质心系数的改变对配平迎角的影响量有所减小；增加或减少1%的纵向质心系数，大约会引起1°左右配平迎角的变化，随着飞行高度增加，纵向质心系数的改变对配平迎角的影响量略有增加。