摘要
采用空气五组分模型对充气式返回舱的气动特性进行了化学非平衡数值模拟研究，考察了返回舱外形变化对流场温度和压力的分布以及壁面热流密度和压强的影响，分析了流场中组分分布情况。研究结果表明，外形变化在总体上对返回舱流场特性影响较小，会使舱体壁面处的热流密度有所增加；对流场组分分布的研究发现，由于N2比O2更难解离，在整个流场中N的摩尔分数远低于O的摩尔分数；在目前的计算条件下，沿轴线和壁面的氮氧比与来流基本保持一致。仿真结果在总体上与实验结果符合较好，验证了仿真模型的可靠性。
Abstract
The aerodynamic characteristics of the inflatable reentry decelerator were studied with chemical nonequilibrium numerical simulation using five-species air model. The effects of reentry decelerator′s shape deformation on distributions of flow temperature and pressure , heat flux and pressure along the wall were investigated, and the mole fractions of each species in the flow filed were presented. The modeling results show that the shape deformation of inflatable reentry decelerator has little influence on flow characteristic, and leads an increase of heat flux on the capsule surface. The study of species distribution shows that since N2 is more difficult to dissociate than O2, the mole fraction of N is much lower than the mole fraction of O in the flow filed. Under the current calculation conditions, the ration of nitrogen and oxygen along the wall and axis is almost consistent with the incoming flow. The modeling study results are in good agreement with the experiment results, which verify the reliability of the model.
随着航天技术的不断发展，传统的刚性返回舱由于质量大、有效载荷小等缺点已经很难满足未来空间任务的需求。充气式返回舱作为一种新型的再入飞行器，具有质量小、收拢体积小、气动阻力面大等优点[1]，具有广泛的应用前景，各航天大国均开展了广泛的研究[2-8]。不同于传统的刚性返回舱，充气式返回舱在高速气流的冲击下会发生结构变形甚至被破坏。因此，充气式返回舱的结构和防热的设计十分重要，而获得准确的气动特性参数是进行充气式返回舱结构和防热设计的重要前提。
为了能够缩短研制周期，降低研制成本，数值模拟是学者们进行充气式返回舱气动特性研究的常用手段[9-14]。王帅等基于Workbench仿真平台对气动热作用下的高超声速充气减速结构的性能进行了研究，比较了气动力和气动热两种气动载荷对蒙皮结构的影响[9]。赵晓瞬等以美国研制的充气式返回舱IRV-3为研究对象，研究了来流气体迎角和速度变化对返回舱气动特性和阻力特性的影响[10]。Guo等采用SST湍流模型对高超声速充气式空气动力减速器（hypersonic inflatable aerodynamic decelerator， HIAD）进行了研究，分析了结构变形对流场结构、蒙皮结构表面压力以及热流密度的影响[12]。充气式返回舱高速再入时，头部形成强激波，波后空气温度急剧升高，并伴随着复杂的解离以及电离过程。目前充气式返回舱大气再入气动特性仿真大多基于理想气体假设，没有考虑空气化学反应对充气式返回舱气动特性的影响。
在高速再入飞行器气动特性的仿真研究中，空气的化学反应过程是不可忽略的[15-16]。例如，在Apollo飞行高度为53 km、飞行马赫数为32.5的数值模拟中，采用理想气体模型得到头部驻点温度约为60 000 K，而实际飞行中头部驻点温度约为10 000 K[17]。吕俊明等对Apollo返回舱的研究表明化学非平衡效应是影响气动力变化的主要原因之一[15]。杨建龙等对高超声速绕流的研究指出热化学非平衡效应使激波距离变小，并且能够更好地体现实际热载荷特性[18]。由此推断，为获得准确的充气式返回舱高速再入的气动特性，需要采用化学非平衡仿真模型来进行返回舱气动特性的研究。
本文采用化学非平衡模型对日本研制的充气式返回舱实验模型HWT-MAAC进行数值模拟研究，研究了充气式返回舱外形变化对流场温度、压力分布以及壁面热流密度和压强的影响，考察了流场中各组分的分布特性。
1 计算模型
本文在N-S方程的基础上添加组分方程来进行充气式返回舱气动特性的研究，计算中采用的二维轴对称可压缩流动的控制方程如式（1）所示。
(1)
式中，U=[ρ， ρu， ρv， ρe， ρfi]T， E=[ρu， ρuu+p， ρuv，（ρe+p）u， ρfiu]T， F=[ρv， ρuv， ρvv+p，（ρe+p）v， ρfiv]T， H= [ρv， ρuv， ρvv+p，（ρe+p）v， ρfiv]T， Ev=[0， τzz， τzr， uτzz-vτzr-qz， Jiz]T， Fv=[0， τrz， τrr， uτrz+vτrr-qr， Jir]T， Hv= [0， τrz， τrr-τθθ， uτrz+vτrr-qr， Jir]T， Sc=[0， 0， 0，-Ec， Ci]T，Ec为化学能损失， Ci为化学反应引起的粒子变化率， fi、Jiz、Jir分别为组分i的质量分数、z方向的扩散通量和r方向的扩散通量。
控制方程中变量的表达式如式（2）~（8）所示。
(2)
(3)
(4)
(5)
(6)
(7)
(8)
其中，u、v分别为速度在z、r方向的分量，ρ、p、e、T分别为气体的密度、压强、内能和温度，μ、k、γ为气体的黏性、热导率和比热比，τ为切应力，hi为组分i的比焓。扩散通量的计算采用Fick定律，表达式如式（9）所示。
(9)
式中，Ji为组分i的扩散通量，Di为组分i的扩散系数。本文采用的湍流模型为SST模型，具体方程参考文献[19]。
本文采用Park的五组分化学反应模型，包含的组分有氮分子（N2）、氧分子（O2）、氮原子（N）、氧原子（O）和氮氧分子（NO），包含的化学反应如表1所示，模型具体参数可参考文献[20]。
表1化学反应
Tab.1 Chemical reactions
本文对控制方程的对流项采用Roe格式离散，对耗散项采用中心差分格式离散。本文以充气式返回舱实验模型HWT-MAAC为研究对象，该模型为飞行演示验证实验中采用的充气式返回舱SMAAC的缩比模型，其结构示意如图1所示，具体的结构尺寸可参见文献[21]。
图2给出了本文的计算域，在流场入口给定来流速度、来流压强和来流温度，出口参数外推，壁面采用无滑移边界条件，壁面温度设为300 K。
图1HWT-MAAC模型结构示意图
Fig.1Model structure sketch of HWT-MAAC
图2计算域
Fig.2Computational domain
2 计算结果
2.1 网格无关性检查
本文首先采用了粗（100×150）、中（200×300）、细（400×600）三种不同数量的网格进行了计算，同时对近壁面网格进行加密保证无量纲壁面距离小于1。计算工况为来流马赫数10、来流温度45 K、来流压强60 Pa，图3给出了采用这三种网格计算获得的沿返回舱壁面热流密度。从图中可以看出粗网格的计算结果显著高于中网格和细网格的计算结果，中网格的计算结果与细网格的计算结果十分接近。因此，为了能够保证计算的准确性，同时也不大幅增加计算量，本文采用中网格来进行计算研究，即网格数为200×300，图4给出了网格划分情况。
图3沿返回舱壁面热流密度
Fig.3Heat flux along the reentry decelerator wall
图4网格划分
Fig.4Grids partitioning
2.2 模型验证
为了验证本文计算模型的可靠性，将上述工况的计算结果与实验相同工况的测量结果进行比较，图5给出了本文计算的沿返回舱壁面热流密度分布与文献[21]实验结果的比较。实验中采用高超声速风洞对充气式返回舱实验模型HWT-MAAC进行了研究，采用辐射平衡假设对表面热流密度进行测量。从图中可以看出在蒙皮处，本文计算得到的热流密度比实验值略高，在其他区域符合较好，计算得到的热流密度变化规律与实验结果一致，总体上计算结果与实验结果符合较好。
图5壁面热流密度与实验结果比较（HWT-MAAC）
Fig.5Comparison of heat flux along the wall with experiment result (HWT-MAAC)
为进一步验证本文模型的可靠性，本文对钝头体模型ELECTRE的高速流场进行了仿真研究。模型总长度为2 m，球形头部半径为0.175 m，半锥角为4.6°，计算工况高度为53.3 km，来流速度为4 230 m/s，温度为265 K，压强为53 Pa[22]。仿真结果与文献[22]实验结果的比较如图6所示。从图中可以看出本文计算结果在总体上与实验结果符合较好。
图6壁面热流密度与实验结果比较（ELECTRE）
Fig.6Comparison of heat flux along the wall with experiment results (ELECTRE)
2.3 外形变化影响
在高速气流的冲击下，充气式返回舱会发生结构变形，Guo等研究了充气式返回舱HIAD外形变化对表面热流密度分布等的影响[12]。由于HWT-MAAC与HIAD的构型不同，经高速气流冲击后外形变化情况也不同。文献[13]采用虚拟粒子方法对充气式返回舱SMAAC在高度为10 km和43 km的变形情况进行了分析，结果表明在10 km高度时，蒙皮弯曲，而充气环几乎没有变形，在43 km高度时，蒙皮弯曲，充气环也在高速气流的冲击下被压瘪[23]。基于以上研究，本文对两种外形变化情况进行研究，图7给出了充气式返回舱迎风面外形变化的示意图。
图8和图9分别给出了充气式返回舱未形变和两种形变情况的流场压强和温度分布情况，计算工况均为来流马赫数10、来流温度45 K、来流压强60 Pa。从图中可以看出高速气流在返回舱前端附近形成弓形激波，气体经过激波后被迅速压缩，气体压强和温度急剧上升，激波前后存在非常大的温度梯度和压力梯度。随着气体向两侧充气环方向流动，气体压强和温度逐渐下降，激波层内沿流场纵向也存在较大的温度梯度和压力梯度。比较发现形变1流场中存在两个高压区，一个在舱体附近区域，另一个在蒙皮凹陷处，未形变时只在舱体附近区域存在高压区，如图8（a）所示。可以推测，当蒙皮凹陷加剧时，凹陷处的压强会继续上升，最终导致蒙皮结构被高压破坏。从图8（b）中可以看出，形变2流场中压强分布规律和未形变时分布规律基本一致。从图9中可以看出两种形变流场温度分布与未形变流场温度分布的差异并不大。
图7迎风面外形变化示意图
Fig.7Sketch of windward shape deformation
图8流场压强分布比较
Fig.8Comparison of pressure distribution in flow fields
图9流场温度分布比较
Fig.9Comparison of temperature distribution in flow fields
在高速气流的冲击下，充气式返回舱迎风面将承受剧烈的气动热作用，迎风面的气动热特性是充气式返回舱防热设计的重要参考。图10给出了外形变化对迎风面热流密度的影响，从图中可以看出未形变时，热流密度在舱体和充气环壁面处均存在峰值。当外形变化时，在蒙皮和舱体连接处也产生一个峰值，且舱体壁面处的热流密度比未形变时热流密度更高，因此在外形变化后对舱体的防热设计提出了更高要求。在蒙皮处，未形变和形变2的热流密度分布类似，向圆环方向逐渐下降，随后又逐渐上升，在充气环处达到峰值。而形变1的热流密度向圆环方向逐渐上升，在充气环处达到峰值。高超声速飞行器前缘驻点热流密度与其驻点曲率半径的平方根成反比，外形变化对舱体和圆环附近的曲率半径产生了较大影响。形变1使得舱体和圆环附近的曲率半径减小，热流密度增加。而形变2使充气环处的曲率半径增大，导致充气环附近的热流密度下降；使舱体处曲率半径减小，导致舱体附近的热流密度增加。
图10外形变化对迎风面热流密度的影响
Fig.10Effect of shape deformation on windward heat flux
返回舱迎风面的气动力主要来源于表面压强的作用，图11给出了外形变化对迎风面压强的影响，从图中可以看出未形变时压强在舱体和充气环处存在峰值，从蒙皮向充气环方向，压强逐渐下降。比较发现，外形变化对舱体壁面的压强影响并不大，在舱体和蒙皮连接处存在一个压强谷值，随后压强逐渐上升。形变1在蒙皮处的压强高于在充气环处的压强，形变2在充气环处存在压强峰值，充气环处的压强高于蒙皮处的压强。
从仿真结果来看，充气式返回舱外形的变化给其自身受热受力环境带来了巨大的影响，因此防热和结构的设计应该兼顾原始外形和不同形变情况来统筹考虑。
2.4 组分分布
化学非平衡模型的一个重要特点在于能够描述流场中的组分演化过程，而流场中的组分直接决定了气体的热力学和输运物性，进一步影响飞行器的气动特性。深入地了解流场中的组分分布特性有助于增进对高速再入流场特性的认识，为仿真计算模型的优化提供思路。
图11外形变化对迎风面压强的影响
Fig.11Effect of shape deformation on windward pressure
组分分布与来流参数存在一定关系，为细致地研究流场中的组分分布特性，本文考察了不同来流参数下流场中的组分分布特性，如表2所示。
表2计算工况
Tab.2 Calculation condition
以工况2、3为例，介绍流场中组分分布特性，图12~13给出了流场中各组分的摩尔分数分布，两个工况流场中的最高温度分别为4 163 K和6 830 K。从图12~13中可以看出由于激波后气体温度急剧上升，N2和O2开始解离，N2和O2摩尔分数下降，N2的最小摩尔分数分别为0.77和0.65，O2的最小摩尔分数分别为0.17和0.02，相比于来流N2摩尔分数分别减少了3%和18%，O2摩尔分数分别减少了16%和92%。产生这种差异的主要原因是N2的解离能为9.759 eV，O2的解离能为5.115 eV[24]，N2比O2更难解离。这也导致在整个流场中N的摩尔分数远低于O的摩尔分数，它们的分布规律也存在较大差异，N的最大摩尔分数在稍微远离壁面附近，而O的最大摩尔分数在充气环附近。对比图12（e）和图13（e）可知，NO的摩尔分数分布存在较大差异，工况2最高摩尔分数在充气环附近，工况3最高摩尔分数在激波和壁面之间的位置。由表1可知，NO的来源有反应3、5的逆反应和反应4的正反应。由图12（c）可知N在充气环处的含量十分低，由此推断工况2中NO在充气环处的主要来源是由N2和O的反应获得，即反应4的正反应。而工况3中N的摩尔分数相比于工况2已有极大上升，前述的三种途径均为NO的重要来源。这也反映出来流参数不同，流场中的化学动力学过程也存在较大差异。
图12各组分摩尔分数分布（工况2）
Fig.12Distribution of species mole fraction (case 2)
图13各组分摩尔分数分布（工况3）
Fig.13Distribution of species mole fraction (case 3)
图14、15给出了表2中各工况沿轴线和壁面的氮氧比分布。其中氮包括N2、N、NO，氧包括O2、O、NO，氮氧比计算如式（10）所示。
(10)
式中，rN/O表示氮氧比，yi表示各组分的摩尔分数。
从图14、15可以看出激波前气体主要成分是N2和O2，氮氧比保持稳定。激波后，不同工况沿轴线的氮氧比分布变化存在一定差异。但是，从总体来看不论沿轴线，还是沿壁面，各工况氮氧比的变化都非常小，在3.82左右，与来流基本保持一致。
图14沿轴线氮氧比分布
Fig.14Nitrogen-oxygen ratio distribution along the axis
图15沿壁面氮氧比分布
Fig.15Nitrogen-oxygen ratio distribution along the wall
在大多数多组分气体的流动中，由于流场内各组分扩散速度不同，流场内混合气体的比值与计算域入口混合气体的比值是不相等的，如电弧加热发动机、自由燃烧电弧、等离子体发生器以及自由射流等的流场[25-29]。而在本文的研究中，激波前，气体温度较低，化学反应较弱，气体组分主要为N2和O2且参数梯度较低，组分扩散作用较弱，氮氧比保持恒定。激波后，随着气体温度的上升，化学反应开始变得剧烈，气体成分变得复杂，组分参数梯度增加，但氮氧比仍然与来流十分接近。由此推断，组分扩散的作用对组分分布的影响相对较弱，组分分布主要还是由各组分的对流和化学反应主导。
组分扩散的计算由于涉及组分间两两之间的扩散以及流场中多种扩散驱动力的耦合作用，计算量较大且比较复杂，并且随着考虑组分的增多，计算量和复杂程度增加十分显著。为此，学者们提出了一些简化的计算方法，比如本文采用的Fick定律以及Murphy提出的组合扩散系数法等[30-31]。而本文的研究发现组分扩散的作用相对较弱，是否可以舍去组分扩散以简化计算模型值得进一步研究。当然，本文目前的研究还可进一步改进。本文采用的扩散计算方法为Fick定律，也是目前学者们常用的扩散计算方法，该方法考虑了组分浓度梯度对流场的影响。但高超声速飞行器再入流场中除存在浓度梯度外，还存在压力和温度梯度。此外，Fick定律还存在扩散通量不守恒的问题[32-33]。本文采用的是空气五组分模型，若考虑更多组分，是否能够得出一致的结论也值得进一步研究。种种迹象表明，对于高超声速飞行器再入流场组分扩散值得进一步研究，研究结果可为模型的合理优化提供理论参考。
3 结论
本文采用化学非平衡模型对充气式返回舱进行了数值模拟，研究了外形变化对流场温度和压力、壁面热流密度和压强分布的影响，考察了流场中各组分的分布特性。研究结果表明：
1）形变1使得舱体和圆环附近的曲率半径减小，热流密度增加。而形变2使充气环处的曲率半径增大，导致充气环附近的热流密度下降；使舱体处曲率半径减小，导致舱体附近的热流密度增加。
2）外形变化对舱体壁面的压强影响不大，形变1在蒙皮处的压强高于在充气环处的压强，而形变2在充气环处的压强高于蒙皮处的压强。
3）由于N2比O2更难解离，O2相比于来流减少得更多，在整个流场中N的摩尔分数远低于O的摩尔分数，N的最大摩尔分数在稍微远离壁面附近，而O的最大摩尔分数在充气环附近。
4）在目前的计算条件下，沿轴线和壁面的氮氧比与来流基本保持一致，由此推断组分扩散的作用较弱，组分分布特性由各组分的对流和化学反应主导。
