大型航天器再入解体气动力热特性模拟的直接模拟蒙特卡洛方法研究

梁 杰，李志辉，2，李绪国，杜波强

(1. 中国空气动力研究与发展中心超高速空气动力研究所， 绵阳621000；2. 北京航空航天大学北京前沿创新中心国家计算流力学实验室，北京100191)

1 引言

服役期满大型航天器如类天宫飞行器离轨再入大气层跨越自由分子流、稀薄过渡流和连续流等多个流动区域，一般在连续流区域是通过求解Euler 或N-S 方程的计算流体力学(CFD)准确模拟，在稀薄流区域则采用直接模拟蒙特卡洛(Direct Simulation Monte Carlo，DSMC)方法进行计算。 飞行器作为一个整体在连续流域飞行时，在高超声速流场中也可能存在局部的稀薄流动区域或高度非平衡流区，如激波、边界层或尾迹区。 在这些区域，传统的CFD 方法失效，而DSMC 方法在计算高密度区域的粒子碰撞时花费巨大，因此发展了CFD-DSMC 混合方法来求解这类多尺度复杂流动问题[1-2]。 这种混合方法的计算精度和计算效率都比较高，但程序代码复杂，2 种不同类型方法之间的信息交换困难，DSMC 方法的统计散布或统计波动极大地影响CFD 计算的稳定性。同时也发展了一些基于粒子模拟方法的混合技术(如平衡取样技术[3-4]、基于DSMC 的低扩散粒子方法[5]、eDSMC[6])来模拟近连续流，这类混合方法大多建立在DSMC 碰撞限制器的概念基础之上[7-9]。 在一个限制碰撞的DSMC 模拟中，数值碰撞频率被限制在一个较低的值，假设使用碰撞限制器的流场区域无粘，扩散输运的影响可以忽略。 如果将DSMC 碰撞限制器应用在流场中的强激波和边界层区域，明显的数值扩散和耗散会导致模拟结果有较大的误差[10]。 当采用DSMC 方法模拟高超声速近连续流时，可发现计算区域中有许多网格，其中的流动处于或接近于热力学平衡状态。 网格内的分子经过有限的碰撞次数就可以使分子速度达到Maxwell 平衡分布，更多的碰撞并不能改变当地的分布函数。 如果在这些区域采用限制网格内每个分子碰撞次数(碰撞限制器)的方法，并与大时间步长和大网格尺度相结合，可以节省大量的计算时间。

因此，针对大型航天器再入解体过程气动力热问题模拟，本文提出一种基于碰撞限制器技术的DSMC 混合方法，根据高超声速流动中连续流假设失效的判断准则，将流场区域分解为遵循气体分子速度函数平衡态分布的连续流区和强热非平衡的稀薄流区。 在连续流区域，采用限制碰撞的DSMC 方法，通过采用网格自适应和变时间步长技术来减少数值扩散的误差[11]。 同时采用大规模并行算法[12]，实现大型复杂飞行器二次解体过程高超声速近连续过渡流气动力热特性的模拟，以解决类天宫飞行器首次、二次解体过程中的高超声速气动热问题。 通过计算类天宫飞行器两舱结构体在低密度风洞试验状态的气动力/热系数，并与低密度风洞试验数据进行对比，以验证该算法的可靠实用性。 最后计算分析带太阳电池帆板的类天宫飞行器激波/激波、激波/边界层干扰及流动分离的复杂流场结构，为数值预报天宫飞行器首次、二次解体的气动力热特性可靠求解提供平台基础。

2 数值模拟方法

2.1 碰撞限制器技术的DSMC 方法实现

对于大型航天器再入解体过程中的绕流，因其大尺度非规则，前后不同区域的局部高超声速流场中，混合有气体分子速度分布函数呈现平衡与非平衡的流动，根据基于DSMC 碰撞限制器的混合方法，利用连续流失效参数来指定采用限制碰撞的DSMC 流场区域和当地自适应DSMC 计算区域。 连续流失效参数根据Boyd 等的研究[13]，采用式(1)所示基于当地流动梯度长度(GLL)的克努森数(Kn)来判断:

式中，梯度取流动参数的最大梯度方向， Q是某些量如密度、压力、温度或速度。 本文采用基于密度的KnGLL来进行区域分解，KnGLL＞0.02 的区域定义为非平衡区域，其它区域则为平衡区域。采用碰撞限制器的DSMC 方法计算。 由于整个流场都统一采用基于DSMC 的方法进行计算，因此在穿越连续流失效边界时并不存在额外的信息交换。

对于精确模拟，标准的DSMC 方法要求足够小的计算网格、足够小的时间步长和每个网格内足够多的模拟粒子数。 其网格尺寸应该小于当地分子平均自由程，时间步长间隔应小于当地分子平均碰撞时间，才能得到有效的结果。 但在高超声速流场的平衡区域，由于穿过每个网格的流动都为常量，所以时间步长和网格尺寸可以大许多倍。 网格中的每个粒子经过一定次数的碰撞后，当地的速度分布函数就能接近Maxwell 平衡态分布。 碰撞限制器就是把网格中每个粒子在一个时间步长之内的碰撞次数限制为一个常数值χ，即在一个时间步长内把一个网格内粒子总的碰撞次数限制为χ·N 次，其中N 是这个网格内模拟的粒子数。 不少文献对碰撞限制器DSMC 中χ 值的影响进行了检验[14-16]，本文算例中选取χ ＝2。与在连续流平衡区域采用真实的高碰撞频率相比，每个时间步长内每个网格总的粒子碰撞次数限制为2N 次，可以减少粒子碰撞处理过程的计算时间。

2.2 计算网格及网格自适应

本文采用直角与表面非结构混合网格的DSMC 数值方法。 在描述物面几何形状的非结构网格建立以后，直接将其嵌入直角网格的流场中，使DSMC 计算对流场网格的依赖程度大大降低。同时通过判断分子运动轨迹方程和物面三角形面元上任一点的位置方程，唯一确定出分子与物面的碰撞点坐标[17]，解决了这种混合网格流场分子运动与物面碰撞的难题。 对分子在物体三角形面元上碰撞、反射前后的流场参数进行统计取样，就可以获得飞行器的整体气动力特性以及表面力、热载荷分布。

在流动梯度较大的区域，碰撞分子应该通过有效的最相邻搜索方法选取。 本文在高度非平衡区采用基于密度梯度的动态自适应碰撞网格和取样网格，避免网格内随机选取的碰撞分子之间的距离远大于当地的分子平均自由程。 即在背景网格的基础上，根据流场中密度梯度的变化分别对碰撞网格和取样网格进行细化，碰撞分子则是在自适应后最小的亚网格内选取，以保证计算的空间精度。 但是为保证高密度区域精细的空间分辨率而增加总的模拟分子数，会导致低密度区域的粒子过度增多，严重降低计算效率。 为了使计算区域每个网格更均匀地分布模拟粒子，解决方法之一就是采用变时间步长格式。

2.3 变时间步长格式

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

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

3 计算结果分析

3.1 类天宫飞行器风洞试验状态气动力特性模拟验证

类天宫飞行器两舱结构模型的气动力试验是在中国空气动力研究与发展中心的高超声速低密度风洞完成的，试验来流气体是M12 型面喷管喷出的氮气，采用2%的模型缩比，模型总长为0.202 89 m，模型参考面积为3.525 66×10-3m2，俯仰力矩的参考点位于头部顶点。 表1 给出了试验状态参数，试验的迎角范围为-2 ° ～25 °。DSMC 计算中采用VHS 分子模型、L-B 内能交换模型，壁面采用固定壁温300 K 的完全漫反射模型。

表1 类天宫飞行器气动力试验状态参数Table 1 Test conditions for Tiangong-type vehicle

虽然试验状态的克努森数Kn 非常低，按传统流区划分，其模拟高度处于连续流区域，但是采用本文构建的碰撞限制器技术和当地参数自适应的混合DSMC 方法进行计算，既能保持计算的高精度又能大幅度提高计算效率。 图1 分别给出了计算的轴向力、法向力和对头部顶点的俯仰力矩系数与试验结果的对比情况，可看出计算与试验得到的气动力系数随迎角的变化规律完全一致。本文DSMC 模拟得到的轴向力系数与试验结果在0°迎角的最大偏差为5%，法向力系数在20°迎角的最大偏差7%左右，俯仰力矩系数在20°迎角的最大偏差在5%附近。 对于相对细长的类天宫飞行器模型，计算的气动力系数与低密度风洞试验数据的最大偏差均在7%以内，表明本文算法可以初步实现大型复杂航天器跨流域气动力特性的模拟能力。

图1 类天宫飞行器气动力系数DSMC 计算与试验比较Fig.1 Comparison of Tiangong-type aerodynamic coefficients between DSMC results and experimental data

图2～3 分别给出了0°和20°迎角下试验状态流场的密度、压力和马赫数等值线分布云图，可以看出复杂绕流现象及头部区域的激波/边界层流动干扰特征，表明基于当地流动参数自适应的DSMC 方法对于捕捉复杂的流动干扰细节是有效的。

图2 类天宫飞行器0°迎角流场等值线分布Fig.2 Flowfield contours distribution of Tiangongtype vehicle at 0 degree angle of attack

3.2 带太阳电池帆板的类天宫飞行器气动力/热复杂流动模拟分析

图3 类天宫飞行器20°迎角流场等值线分布Fig.3 Flowfield contours distribution of Tiangongtype vehicle at 20 degree angle of attack

为了研究寿命末期的天宫1 号目标飞行器无控再入过程中的熔融烧蚀和解体情况，需要准确模拟跨流域不同高度的气动力、气动热详细分布及流场干扰特征。 图4 绘出了带太阳电池帆板的类天宫飞行器在克努森数Kn ＝6.7×10-4(特征长度10.4 m)、0°迎角下的压力等值线分布，计算的来流条件是:速度7500 m/s、高度85 km。 图中可以清晰地看到复杂激波干扰下的流动特征:头部简化的对接机构产生的脱体激波，在锥段压缩激波以及与实验舱连接处膨胀波的干扰下向外延伸，并与太阳电池帆板前方的脱体激波相交，产生的反射激波与帆板边界层相互干扰，使帆板局部表面压力急剧上升。 另有多道反射激波从帆板脱体激波内部与之相交，使帆板脱体激波形成波浪形状。 在舱体表面附近区域，由于气体的膨胀，压力迅速下降，因此在头部脱体激波、舱体和太阳帆板表面之间形成较大的逆压梯度，引起较强的流动分离，产生一个较大的分离主涡和两个较小的分离附着涡，如图5 所示。 还可看出，从反射激波作用在帆板表面边界层的位置开始，沿帆板左右两侧各形成一股射流，朝向帆板外侧的射流在帆板表面附近仅引起微弱的流动分离(参见图5)；而朝向舱体的射流在较大的压力梯度下加速至超音速并冲向资源舱的表面，形成局部的脱体激波并在资源舱表面引起局部高达120 kW/m2左右的热流分布，如图6 所示。 太阳电池帆板与舱体间构成的这种特殊非规则绕流冲击形成的强气动热，导致的结构响应变形最终撕裂了太阳电池帆板，使之发生首次解体。

图4 带帆板的类天宫飞行器0°迎角压力等值线分布Fig.4 Pressure contours of Tiangong-type vehicle with solar panels at 0 degree angle of attack

图5 带帆板的类天宫飞行器0°迎角流线分布Fig.5 Streamline distribution of Tiangong-type vehicle with solar panels at 0 degree angle of attack

图6 激波干扰引起的局部表面高热流Fig.6 High heat flux on local surface induced by shock wave interactions

4 结论

1)本文针对大尺度复杂结构航天器再入绕流不同位置呈现多流区局部流场特点，构建了基于碰撞限制器技术的全粒子混合方法，并进一步发展了基于密度梯度的动态自适应混合网格处理技术与变时间步长方案。 根据连续流失效参数，将流场区域分解为平衡的连续流区域和非平衡的稀薄流区域，在平衡区域采用大尺度网格和大时间步长的碰撞限制器DSMC 方法，而在非平衡区域则采用网格自适应和变时间步长的DSMC 方法，在提高计算精度的同时，节省了大量计算时间；

2)使用本文算法对类天宫飞行器两舱体风洞试验状态气动力特性进行不同迎角的气动特性计算，解得的气动力/力矩系数与低密度风洞试验数据对比一致性较好，偏差7%以内，验证了本文提出的碰撞限制器和当地参数自适应的全粒子混合算法对大型复杂航天器再入气动特性高精度模拟能力；

3)模拟了带太阳电池帆板的类天宫飞行器85 km 高度再入时的高超声速复杂气动热及流场激波干扰结构，揭示了头部脱体激波、舱体和太阳电池帆板非规则绕流激波引起流动分离、产生强气动热和结构响应变形导致太阳电池帆板毁坏发生首次解体的机制。