充气式再入与减速技术(inflatable reentry and descent technology，IRDT)是一种新型的大气再入减速技术^[1]，该技术特点在于其减速结构是由耐高温柔性材料制成的充气式结构。与传统返回舱式减速技术相比，充气式减速技术具有重量轻、体积小、成本低等优点，拥有广泛的应用前景。充气减速结构作为IRDT系统的核心组成部分，有多种不同的外形设计方案^[2-4]，其中典型的充气减速结构为多个充气圆环组成的锥形结构^[5-6]。

在IRDT系统减速的过程中，会在弓形激波与锥形减速结构表面形成高温激波层，对外部热环境进行准确的预测是IRDT系统热防护设计的前提。以氧化铝纤维为代表的充气减速结构热防护材料的耐热性能要低于传统的返回舱热防护材料^[7]，这也对充气式减速结构热环境预测提出了更高的要求。对于航天器再入大气层的热环境问题，试验手段包括高空飞行试验^[8-9]和高焓风洞试验^[10-11]，高空飞行试验成本高且风险大，而高焓风洞试验数据的分析、外推和使用与常规风洞试验相比更加困难和复杂^[12]。目前高超声速流场气动热的计算方法主要有工程算法、工程数值耦合方法及数值仿真方法。工程算法^[13-15]借助经验公式求解驻点热流密度，计算效率高且可以保证一定的精度；工程数值耦合方法^[16]对边界层外无黏流场进行数值求解，对边界层内流场用工程算法求解，可应用于复杂外形；数值仿真方法主要包括求解连续流动的计算流体动力学(computational fluid dynamics, CFD)方法^[17-19]和求解稀薄流动的直接模拟Monte Carlo(direct simulation of Monte Carlo，DSMC)方法^[20]，虽然可以获得复杂流场的流动特征，但是计算量大，特别是在计算包含复杂的热化学反应时。

在高超音速流场中激波层温度非常高，分子内的振动能得到激发，同时会引起边界层内的气体出现离解甚至电离^[21]，因此想要得到正确的计算结果，数值仿真计算必须考虑激波层内的热化学反应。本文基于有限体积方法，采用双温非平衡热化学反应模型，对锥形减速结构的外部热环境进行仿真计算，分析了沿驻点线的温度变化与化学组元浓度变化，以及结构表面的热流及压强分布；同时，对4种不同半锥角的锥形减速结构进行仿真计算，对减速结构外部热流及压强受半锥角变化的影响进行研究。

1 计算方法 1.1 非平衡热化学反应N-S方程

气体的内能模式包含分子的平动能、转动能、振动能及电子能^[21]。在低温下，只有分子的平动能和转动能处于激发状态并达到平衡，气体温度才可以用单一平动温度T进行描述。在高温气体中，振动能及电子能被激发，其能量传递需要用振动温度T_v及电子温度T_e进行描述。在双温模型^[22]中，原子振动能及电子能处于平衡状态，其能量传递可以统一用T_v进行描述。基于双温模型，包含N种不同化学组元的气体在非平衡状态下的N-S(Navier-Stokes)方程积分形式可表示为

$ \frac{\partial}{\partial t} \int_V \boldsymbol{Q} \mathrm{d} \sigma+\oint_{\partial V}\left(\boldsymbol{F}_{\mathrm{c}}-\boldsymbol{F}_{\mathrm{v}}\right) \mathrm{d} s=\int_V \boldsymbol{S} \mathrm{d} \sigma . $

(1)

其中：F_c和F_v分别为对流通量及黏性通量，S为源项，t为时间，V和ƏV分别为控制体体积及表面积，dσ和ds分别为体积微元和面积微元，Q为流场守恒变量。Q可表示如下：

$ \boldsymbol{Q}=\left[\begin{array}{c} \rho \\ \rho_1 \\ \vdots \\ \rho_{N-1} \\ \rho u \\ \rho v \\ \rho w \\ E \\ E_{\mathrm{v}} \end{array}\right], $

(2)

$ E=\sum\limits_{i=1}^N \rho_i h_i+\frac{1}{2} \rho|\boldsymbol{V}|^2, $

(3)

$ h_i=\left[\frac{H_i^0(T)-R}{M_i}\right]+e_{\mathrm{v}, i}, $

(4)

$ E_{\mathrm{v}}=\sum\limits_{i=1}^N \rho_i e_{\mathrm{v}, i}, $

(5)

$ e_{\mathrm{v}, i}=\frac{R}{M_i} \frac{\theta_i}{e^{\theta_i / T_{\mathrm{v}}}-1} . $

(6)

其中：ρ为气体密度，ρ_i为第i种气体组元密度，i=1，2，…，N-1；速度矢量V=(u, v, w)，u、v和w为速度分量；E为气体的总内能，E_v为气体的振动内能；h_i为单位质量气体组元的焓；H_i⁰(T)为每摩尔气体组元在温度T下的焓值，可通过文[23-24]获得；R为理想气体常数；M_i为气体组元相对分子质量；e_{v, i}为单位质量气体组元的振动能；θ_i为气体组元振动特征温度。

对流通量F_c可表示如下：

$ \boldsymbol{F}_{\mathrm{c}}=\left[\begin{array}{c} \rho V_n \\ \rho_1 V_n \\ \vdots \\ \rho_{N-1} V_n \\ \rho u V_n+p n_x \\ \rho v V_n+p n_x \\ \rho w V_n+p n_y \\ (E+p) V_n \\ E_{\mathrm{v}} V_n \end{array}\right] . $

(7)

其中：P为气体压强；n(n_x, n_y, n_z)为控制体表面ƏV法向量，x、y、z分别为3个坐标轴方向；V_n=V·n为流体速度在n方向投影，n为法向量n的方向。

扩散通量F_v可表示如下：

$ \boldsymbol{F}_{\mathrm{v}}=\left[\begin{array}{c} 0 \\ D_1 \frac{\partial \rho_1}{\partial n} \\ \vdots \\ D_{N-1} \frac{\partial \rho_{N-1}}{\partial n} \\ \tau_{n x} \\ \tau_{n y} \\ \tau_{n z} \\ u \tau_{n x}+v \tau_{n y}+w \tau_{n z}+ \\ q_n+q_{\mathrm{v} n}+\sum\limits_{i=1}^N D_i h_i \frac{\partial \rho_i}{\partial n} \\ q_{\mathrm{v} n}+\sum\limits_{i=1}^M D_i e_{\mathrm{v}, i} \frac{\partial \rho_i}{\partial n} \end{array}\right], $

(8)

$ \boldsymbol{\tau}_n=\boldsymbol{\tau} \cdot \boldsymbol{n}=2 \mu\left[{\mathit{\pmb{ε}}}-\frac{1}{3}(\nabla \cdot \boldsymbol{V}) \boldsymbol{I}\right] \cdot \boldsymbol{n}, $

(9)

$ {\mathit{\pmb{ε}}}=\frac{1}{2}\left[(\nabla \cdot \boldsymbol{V})+(\nabla \cdot \boldsymbol{V})^{\mathrm{T}}\right], $

(10)

$ q_n=k \frac{\partial T}{\partial n}, $

(11)

$ q_{\mathrm{v} n}=k_{\mathrm{v}} \frac{\partial T_{\mathrm{v}}}{\partial n}. $

(12)

其中：D_i为气体组元的质量扩散系数，假设各气体组元的二元扩散系数^[25]相同，则D_i=μ/ρSc，其中，Sc为Schimidt数，对于中性粒子组元可取0.5；τ为切应力张量, τ_n为τ在n方向投影，τ_n=(τ_nx, τ_ny, τ_nz)；μ为气体黏性系数；I为单位张量；ε为应变率张量；q_n和q_vn分别为气体平动和振动热流；k和k_v分别为气体平动导热系数和振动导热系数；对于混合气体，k、k_v和黏性系数μ可由Wilke^[26]半经验方法通过各组元的黏性系数和导热系数确定。

源项S表示如下：

$ \boldsymbol{S}=\left[\begin{array}{c} 0 \\ \dot{\omega}_1 \\ \vdots \\ \dot{\omega}_i \\ 0 \\ 0 \\ 0 \\ 0 \\ \sum\limits_{i=1}^N\left(\rho_i q_{\mathrm{T}-\mathrm{v}, i}+\dot{\omega}_i e_{\mathrm{v}, i}\right) \end{array}\right], $

(13)

$ \dot{\omega}_i=M_i\left(\nu_i^{\prime \prime}-\nu_i^{\prime}\right)\left(k_{\mathrm{f}} \prod\limits_{i=1}^N\left(\frac{\rho_i}{M_i}\right)^{\nu_i^{\prime}}-k_{\mathrm{b}} \prod\limits_{i=1}^N\left(\frac{\rho_i}{M_i}\right)^{\nu_i^{\prime \prime}}\right) . $

(14)

其中：$\dot{\omega}$_i为组元生成速率；q_{T-V, i}为平动能-振动能转化速率；ν'_i和ν″_i分别为反应物和生成物的化学反应式计量系数；k_f和k_b分别为正向及逆向反应速率，可由Arrhenius公式^[27]表示，两者形式相同，以k_f为例，

$ k_{\mathrm{f}}=A T^B \exp (-C / T) . $

(15)

其中：A、B和C为反应常数，可通过实验获得。k_f和k_b满足关系

$ \frac{k_{\mathrm{f}}}{k_{\mathrm{b}}}=K_{\mathrm{eq}}=\exp \left(-\frac{\Delta_{\mathrm{f}} G^0}{R T}\right) . $

(16)

其中：K_eq为反应平衡常数，Δ_fG⁰为正向反应过程中Gibbs自由能的变化量^[21]。当k_f为已知时，可通过式(16)计算k_b。

q_{T－V, i}可通过Landau-Teller公式^[28]表示，即

$ q_{\mathrm{T}-\mathrm{v}, i}=\frac{e_{\mathrm{v}, i}^*-e_{\mathrm{v}, i}}{\tau_{\mathrm{v}, i}} . $

(17)

其中：e_{v, i}^*为平衡振动能，τ_{v, i}为振动松弛时间。

本文采用有限体积法对式(1)进行离散求解，迭代采用LU-SGS方法，对流通量采用Roe通量分裂格式，通过MUSCL方法结合Van Leer限制器^[29]进行变量重构获得二阶精度，黏性通量采用中心格式。

1.2 热化学反应模型

对于空气中不同组元热化学反应式的建立及式(15)中各化学反应常数的确定，有多种不同的模型供选择。常用的热化学反应模型包括Park85^[30]、Park93^[31]及Dunn-Kang^[32]模型等。其中，Park85和Park93模型只给出了正向反应的参数值，逆向反应需通过反应平衡常数求得，而Dunn-Kang同时给出了正向及逆向反应的参数值。在实际应用中，为提高计算效率，一般对上述化学反应模型组元进行简化。在工程计算满足精度需求的情况下，5组元Park85模型可以兼顾计算的准确度和效率^[33-34]。本文计算采用基于Park85模型的5组元模型，相应的反应式及反应参数如表 1所示。

热化学反应壁面边界条件可分为非催化壁及完全催化壁^[21]，本文采用非催化壁，即壁面化学反应速率为0。

2 算法验证

为了验证计算方法的准确性，本文首先对钝体标准模型ELECTRE进行仿真计算，并将结果与文[11]飞行试验数据及计算结果进行对比。ELECTRE模型为钝锥形旋成体，总长度为2 m，球形头部半径R为0.175 m，半锥角为4.6°，如图 1所示。流场环境参数与文献中飞行试验一致，即高度为53.3 km，来流速度为4 230 m/s，温度为265 K，压强为53 Pa，Ma为13，密度为6.94×10^-4 kg/m³，壁面温度为343 K。计算网格如图 2所示，基于1/4外形划分结构网格，网格总数为36万，其中壁面法向网格节点数为80。为验证计算结果的网格相关性，壁面网格法向尺寸分别设置为1.7×10^-5 m和1.7×10^-4 m，对应网格Reynolds数Re_c分别为3和30。

部分变量的归一化残差迭代收敛曲线如图 3所示，各变量残差值均收敛于稳定结果。壁面热流计算结果与飞行试验数据及文[11]计算值对比如图 4所示，X即图 1中x方向坐标值, 其中文[11]使用模型为5组分Dunn-Kang模型。两套网格计算结果整体相差在3%以内，可认为计算结果满足网格无关性。在驻点附近本文结果与试验及文[11]结果一致，在流场下游3者变化趋势相同，本文计算结果低于文[11]及试验值。造成结果偏差的主要原因为两者选取的化学反应模型不同。在飞行试验中壁面催化效应及高温气体的热辐射效应都会导致壁面热流增加，但在计算中没有考虑这两方面因素。

沿驻点线平动温度和振动温度计算结果与文[11]对比，如图 5所示，X即图 1中x方向坐标值。定义组元浓度Y_i=ρ_i/ρ，沿驻点线组元O及NO浓度与文[11]对比如图 6所示，X为图 1中x方向坐标值。本文激波脱体距离与文[11]结果一致，激波后平动温度T略低于文[11]结果。T与T_v变化趋势及组元浓度变化趋势一致，因选取的化学反应模型不同，本文激波后组元离解速度更高，导致温度计算结果低于文[11]值。本文与文[11]计算结果对比及误差如表 2所示。

3 锥形减速结构流场仿真

本文计算所用锥形减速结构如图 7所示，R为0.48 m，半锥角为60°，底边半径为2.3 m。流场环境为高度70 km，来流速度3 861 m/s，压强5.2 Pa，温度220 K, Ma为13，密度8.28×10^-5 kg/m³，壁面温度1 600 K。锥形减速结构计算网格如图 8所示，基于1/4外形划分结构网格，网格总数为21.6万，其中壁面法向网格节点数为80，壁面网格法向尺寸为6.0×10^-4 m，网格Re_c为3。

沿驻点线平动温度及振动温度变化曲线如图 9所示, X即图 7中x方向坐标值。激波距离前缘约0.05 m，激波后温度约7 000 K，激波后空气振动能被激活，振动温度逐渐升高，在前缘附近平动温度与振动温度达到平衡。沿驻点线组元O、NO及N浓度变化曲线如图 10所示, X即图 7中x方向坐标值。激波后空气中N₂及O₂发生离解反应生成O及N，进而通过置换反应生成NO。因N₂离解反应活化温度(即表 1中系数C)较高，N组元浓度较低且在驻点处几乎完全复合为N₂。O及N组元浓度持续升高，到驻点附近分别达到最大值。

在IRDT方案设计中，柔性结构热响应及受力变形分析需要以表面热流及压强分布作为边界输入条件。锥体表面热流及压强沿展向变化曲线如图 11所示, Y即图 7中y方向坐标值。在驻点处热流及压强最大值分别为84 kW/m²和1 140 Pa，在驻点附近热流及压强沿展向快速降低，在距离驻点展向距离0.5 m以后，热流呈线性缓慢下降，压强近似为常量。在实际应用中，锥形结构中心部位设计为刚性防热结构，并与外部柔性结构相连接，根据上述计算结果可知，当柔性结构内缘尺寸满足条件时，其表面热流及压强分布可以认为是线性，柔性结构热响应及受力变形分析边界输入条件可以得到简化。

为了分析锥形减速结构半锥角变化对表面热流及压强分布的影响，本文比较了50°、55°、60°及65°共4种半锥角仿真结果，不同半锥角下沿驻点线平动温度变化对比如图 12所示。50°、55°及60°半锥角激波位置、激波后最高温度及温度变化基本相同；相比前3种半锥角，65°半锥角激波脱体距离增加了0.027 m，激波后最高温度略低，约为6 870 K。

不同半锥角沿展向表面热流对比如图 13所示。50°、55°及60° 3种半锥角表面热流结果相近，当半锥角增加到65°以后，表面热流出现整体下降，幅度约为10 kW/m²。激波脱离距离对表面热流产生影响，一方面因为激波层厚度增加使温度梯度降低，另一方面激波层内反应时间增加导致空气离解更加充分。不同半锥角沿驻点线组元浓度对比如图 14所示，65°半锥角的2种离解气体组元浓度均高于60°半锥角。

不同半锥角沿展向表面压强对比如图 15所示。除65°半锥角略高外，其他3种半锥角驻点压强基本相等。在锥体外围表面压强随半锥角增大而升高，分别为747、856、962和1 079 Pa，表面压强与半锥角变化近似为线性关系。

4 结论

本文对IRDT方案中使用的锥形减速结构流场进行了化学非平衡气动热仿真计算。计算模型基于有限体积法对N-S方程进行求解，利用5组元非平衡双温度模型计算流场热化学反应，壁面化学反应条件为无催化条件。对钝体标准模型ELECTRE进了仿真计算，仿真结果与飞行试验及文献计算结果基本相符。

对半锥角为60°的锥形减速结构在70 km高度、来流Ma为13工况下的流场进行了仿真计算，结果显示：激波后气体平动温度约为7 000 K，气体振动能被激活后沿驻点线平动能向振动能转化，振动温度逐渐升高，平动温度和振动温度在驻点附近达到平衡并降低至壁面温度。激波层内N组元浓度较低，并在驻点处重新降低为0；O及NO组元沿驻点线逐渐升高，在驻点处达到最大值；表面热流及压强在驻点处为最大值，在驻点附近沿径向快速降低，随后热流呈线性下降，压强近似为常量。

对50°、55°、60°及65°共4种半锥角锥形减速结构仿真结果进行对比，结果表明：50°、55°和60°半锥角激波位置及表面热流基本相同；相比前3种半锥角，65°半锥角激波距离前缘点更远，同时表面热流更低；4种半锥角驻点压强基本相同，外围压强随半锥角增加呈线性增加。