再入飞行器烧蚀热防护一体化计算方法

2021-08-03周印佳张志贤付新卫阿嵘

航空学报 2021年7期

周印佳，张志贤，付新卫，阿嵘

中国空间技术研究院钱学森空间技术实验室，北京 100094

烧蚀热防护是再入飞行器最成功的防热方法之一，但烧蚀热防护系统响应的计算是也非常复杂的，涉及到平衡或者非平衡化学反应、烧蚀机理、烧蚀外形变化、烧蚀率和烧蚀量的预测以及结构热响应等诸多因素[1-2]。国外使用最广泛的烧蚀预测数值方法涉及到计算流体动力学(Computational Fluid Dynamics，CFD)的边界层法(BLIMP程序)和材料热响应的一维内部热传导、CMA程序(Aerotherm’s Charring Material Thermal Response and Ablation Program)。随着计算机硬件的发展，近年来基于Navier-Stokes(N-S)方程方法耦合材料热响应的预测逐渐增多。文献[3]利用N-S方程求解轴对称再入流场与材料热响应的耦合，基于表面能量平衡方法计算出表面温度和烧蚀率，流场计算和材料响应计算以准静态或瞬态方式在每个时间步上进行耦合。计算中表面会发生后退和网格动态变化，并且流场网格在固定的弹道点上进行重画以适应新的表面形状。文献[4-5]将轴对称N-S程序与CMA程序耦合，采用松耦合的方法迭代得到CFD计算的热流、CMA计算的表面温度以及烧蚀率引起的表面温度和流率变化低于1%。计算中允许烧蚀产物与非平衡气体反应，但是表面不产生后退。文献[6]采用类似的计算方法对轴对称流场进行计算，并且计算中允许表面后退。文献[7-8]将轴对称热化学非平衡流场计算和二维固体热传导计算耦合。流体和固体的耦合通过表面质量和能量平衡实现，并假设了一个包含表面化学非平衡的热化学烧蚀模型。文献[9]采用松耦合算法耦合了热化学非平衡N-S程序和考虑表面烧蚀和形变的多维内部热传导。将N-S程序得到的气动参数结果视为常数，材料响应程序按时间推进计算烧蚀、形变表面温度和内部温度。但是该方法中烧蚀过程的高度非线性导致热传导和烧蚀率计算的严重不稳定性，为了解决该问题，文献[10]开发了一种迭代技术以消除前面所述的不稳定性，并显著提高烧蚀过程的精度。

国外的气动热预测和材料响应计算通常采用在SACCARA[9]、GIANTS[11]、COYOTE II、FIAT等内部代码的基础上进行修改和应用的方式。而近年来得益于商用CFD和有限元软件的快速发展，采用商业软件耦合预测气动热和材料热响应的研究也逐渐增多。文献[12]利用Fluent与LS-DYNA商用软件，以松耦合的方式求解了零攻角轴对称飞行器的非平衡化学反应流动，并采用网格运动算法实现再入飞行器烧蚀导致的表面后退模拟。

以上数值方法通过求解N-S方程和高精度的CMA程序或有限元程序，提高了计算气动参数和烧蚀响应的精度，但是该方法计算负担较重，对计算能力要求较高，并且CMA高精度程序很复杂，尤其是考虑到程序输入所需要的大量数据更是如此。另外，生成表面化学输入和热解气体焓需要运行额外的平衡或者非平衡化学程序。化学程序本身需要大量的输入，并需要材料组分以及热解气体组分的相关知识。为了完成这种分析，需要开发热响应模型并与电弧风洞试验数据相关联，还需要一组可用的包含大量壁面条件的表面化学表，还要开发热解气体焓模型。如果这些模型都没有，热解气体和热防护材料的确切成分都未知，建立和运行像CMA这样的高精度模型无论是从准确度还是从结果收敛性上都是有问题的。

当高保真的热烧蚀响应模型无法实现时，尤其是在需要进行大量快速计算的初样设计阶段，需要多次迭代以求得防热层厚度，相对不那么复杂的计算方法显然更实用。这些近似方法通常采用稳态烧蚀假设并采用烧蚀热Q*来预测再入过程中的烧蚀后退量。这些方法常采用近似解析方程来确定材料内部的瞬态一维热传导。采用近似热传导方程求解在指定的热防护系统结构界面处温度条件下的所需绝热层厚度。这些方法能够给设计者提供大概的答案，但是其缺点在于这些方法的稳态烧蚀假设在飞行中通常是不成立的，一般会高估后退量。另外一个不足是采用的近似瞬态热传导方程一般只在半无限固体或平板假设成立时才有效。

本文所提出的计算方法旨在提供一种具有较高精度，并能够快速应用的烧蚀热防护一体化计算手段，适用于缺少高精度材料热响应模型的初步设计阶段。该方法利用定义好的飞行器几何数据，采用工程算法估算驻点的对流热流和辐射热流，整合具有较高精度的烧蚀模型，通过Landau变换简化烧蚀后退带来的节点删除过程并保证空间离散精度，最后求解瞬态有限差分热传导来获得不同尺寸、形状和材料的热防护内部响应。

1 物理模型与计算方法

1.1 气动热环境

采用工程算法预测驻点的对流热流和辐射热流。Sutton-Graves冷壁对流热流计算公式[13]为

(1)

式中：K为大气化学的函数，针对地球大气和火星大气分别取为K=1.741 5×10-4和K=1.902 7×10-4;Rn为球头半径；ρ∞为自由来流密度；V为飞行速度。式(1)已经过众多地面和飞行验证，对于钝头体的热流预测偏差通常在5%～10%。

通常气体温度在8 000 K以下时辐射热流可以忽略不计，但是随着再入速度的增大，辐射热流所占比重逐渐增大，例如Apollo飞船再入辐射热流可占总热流载荷的40%～60%[14]，此时辐射热流不可忽略。驻点的辐射热流采用Tauber-Sutton公式[15]计算，该公式适用于地球或者火星再入。Tauber-Sutton热流公式的表达式为

(2)

在弹道上的每个点，将对流热流和辐射热流相加得到总加热热流。

1.2 热化学烧蚀计算模型

如图1所示,烧蚀材料表面能量平衡方程可以表示为

(3)

(4)

图1 烧蚀材料表面能量平衡及内剖面原理示意图

qnet=qconv(1-hw/hr)φHALφblow+qrad-

(5)

式中：hr为恢复焓；ε为辐射系数；σ为斯忒藩-玻耳兹曼常数；Tref为参考温度；φHAL为侵蚀系数，由Hove-Shih理论可知φHAL≥1[17]，对于清洁大气飞行取φHAL=1[18]；φblow是由气体进入边界层导致的，根据边界层的薄膜理论可以由下式计算[19-20]

(6)

式中：吹气系数a′取值范围为0.3～1.3[21]，对于碳-酚醛材料取a′=1.3，碳-碳材料取a′=1.5[18]。

正常质量流量的计算式为

(7)

未校正质量流量B′0的计算式为

(8)

未校正的热传递系数gh0=(qconv/hr)φHAL，热传递系数gh=gh0φblow。壁内焓hu的计算式为

(9)

式中：Cp为气体比热；C∞=2.3 J/(g·K)；D=800 K；Tref=300 K。

为了计算烧蚀速率，可以采用碳-空气表面热化学平衡计算的简化拟合公式[22-23]:

(10)

(11)

在最低烧蚀温度范围内可以采用非平衡反应速率，碳质材料的速率控制氧化有多种表达形式：

(12)

(13)

式中：R为通用气体常数。

将速率控制表达式和平衡表达式合并成总烧蚀率为

(14)

(15)

1.3 Landau变换

对于常见的钝头热防护罩外形，一维模型就能够满足要求。轴对称几何的一维热传导方程为

(16)

式中：λ为固体热导率；Cp,s为固体比热;r为径向。

Landau变换的最初发展是用于解决一维平板的相变问题，这里采用Landau坐标系[24]实现烧蚀的动网格。该坐标系重新定义了空间坐标，从而使在Landau空间内剩余区域的无量纲厚度始终为1，如图2所示并且在求解过程中某一点的Landau坐标始终为常值。图2中z坐标系的原点固定在初始烧蚀界面处，而x坐标系与瞬时烧蚀界面相关联。任意一点的Landau坐标[25]为

图2 一维平板模型坐标系示意图

(17)

式中：η为Landau坐标；L(t)为随时间变化的防热层厚度。

对于轴对称几何外形，z坐标系和r坐标系的关系可以表示为

r=Ro-z,z=Ro-r

(18)

而η坐标系和r坐标系的关系可以表示为

(19)

利用式(19)和微积分中的链式法则，热传导方程可以转换成以下形式:

(20)

(21)

式中：j为节点号；n为时间节点；Δt为时间步长。

因此，任意节点的移动速度可以表征为与表面后退率相关的公式

(22)

烧蚀问题的网格方案通常有平移网格和收缩网格2种方案，平移网格方案的节点位移关系相对来说要更简单，但是会在背壁面位置的网格间产生很大的不连续，如图3所示，可能会导致在空间离散精度上产生负数值结果。而本文采用的基于Landau变换的收缩网格方案避免了复杂的节点删除过程，并且在整个求解过程中可以保证空间离散精度，因为该方案的节点数和无量纲单元厚度Δzj(t)/L(t)是常数[26]。

图3 均布收缩网格方案对比

1.4 有限差分计算

变换后的热传导方程采用有限差分法进行隐式离散化，对时间的偏导数采用后向差分格式：

(23)

对空间坐标的偏导数全部采用中心差分格式：

(24)

(25)

(26)

对于平板模型有式(20)中的m=0，将式(23)～式(26)代入到热传导方程中，得到

(27)

外表面(η=1)为热流边界条件qnet(nΔt)，并假设内表面(η=0)给定绝热边界条件，最终平板形式热防护的Landau坐标系下的热传导方程可化为如式(28)的矩阵形式，同理亦可得到圆柱和球头形式热防护在Landau坐标系下的热传导方程。公式推导过程比较冗长，此处不再详细给出。

(28)

1.5 计算流程

计算的有限差分网格采用均匀网格，初始给定一个防护层厚度以及所需的最大节点间距Δx，并计算出能够保持节点间距略低于或等于所需最大间距的节点数量，然后将节点编号存储在一个数组内。

另外，该方法也可以利用割线法进行迭代，用来预估防热层厚度。只要给定飞行轨迹和背壁温度约束，就可以根据防热层厚度的初始预估值按照图4所示流程进行热响应计算；然后根据割线法的需要，防热层厚度比初始预估值增加0.1 cm再进行一轮计算，如果背壁温度与之前相比没有显著差异，则将初始预估值减小一半，并重复计算过程。一旦这些初始计算完成，割线法就可以按部就班地校正防热层厚度，直到计算所得的背壁温度与要求的背壁温度约束相差小于0.1 ℃。该过程通常经过10～15次迭代即可收敛。

图4 计算流程图

2 结果与讨论

为了验证本文提出的计算方法，以典型钝头体的地球再入和PICA电弧风洞模拟算例分别针对相对高密度和中低密度的材料进行了计算和对比验证。

2.1 碳-碳材料钝头体地球大气再入

高弹道系数(95 760 N/m2)的钝头再入飞行器再入地球大气，高度在30 s内从92 km降低到海平面。再入初始速度为6 000 m/s，驻点半径Rn=0.05 m。热防护材料为碳-碳材料，密度为1 900 kg/m3，材料的热物性参数如表1所示。来流的密度和速度与时间的关系为

表1 碳-碳材料物性

(29)

利用以上再入条件将本文方法与热平衡积分(Heat Balance Integral，HBI)算法的计算结果[18]进行了对比分析，图5和图6分别为壁面温度和烧蚀后退量的结果对比。2种方法计算的峰值温度基本一致，预测的烧蚀后退量相差约0.2 mm。可见本文算法与经典的热平衡积分算法计算结果吻合较好，能够满足工程需要。

图5 表面温度计算结果对比

图6 表面烧蚀后退量计算结果对比

2.2 PICA材料电弧风洞烧蚀模拟

为了测试星尘号(Stardust)飞船的热防护材料PICA的烧蚀性能，NASA开展了一系列电弧风洞试验，同时开展了一系列相应的仿真分析[28]，对实际烧蚀表面温度与试验中热量计测量结果的差异进行对比分析。试验和仿真针对4种不同工况条件，试验中气体为干燥空气加7%质量的氩气，氩气的浓度是试验条件与实际飞行条件最大的不同之处。通过测量铜制热量计的传热速率得出气体最小焓值为30 MJ/kg，而通过分析激波层辐射的频谱得出焓值为40.6±1.4 MJ/kg。

在所有工况计算中，与文献[28]保持一致，PICA材料表面温度均设为3 000 K，PICA材料参数为：密度ρs=240 kg/m3[28]，热扩散系数为 7×10-7m2/s[29]。表2给出了利用本文方法计算的烧蚀速率和烧蚀后退率与文献计算结果的对比。

表2 计算结果对比

在热流和压力相对较低的工况1和工况2条件下，本文方法预测结果与文献结果都吻合较好，烧蚀后退率和烧蚀速率的偏差分别为18%、13.8% 和9.1%、8%。而在压力更高(60.8 kPa)的工况3和工况4条件下，后退量预测出现了较大差异。本文方法预测的结果明显低于文献计算结果，烧蚀后退率和烧蚀速率的偏差分别达到了25%、25.5%和47.7%、47.4%。研究表明[30]，当前采用的烧蚀模型对于非碳化烧蚀材料具有较高精度，而对于碳化烧蚀材料则具有一定局限性，需要进一步完善。PICA材料以碳化烧蚀为主，在高压力高热流条件下，会发生大量的碳化和热解，烧蚀模型需要更精细地考虑材料碳化、热解能量吸收、热解气体通过固体时的对流以及表面化学作用的影响。另外，PICA这样的中低密度烧蚀材料的烧蚀性能对压力高度敏感[31]。因此，针对此类具有高压力高热流的中低密度烧蚀材料的烧蚀计算在实际应用之前需要经过大量的验证。