NUMERICAL STUDY OF HIGH-TEMPERATURE NONEQUILIBRIUM AIR NOZZLE FLOW WITH STATE-TO-STATE MODEL
-
摘要: 采用态−态模型在总温T0 = 2000 K ~ 8000 K、总压p0 = 1 ~ 20 MPa范围, 开展高温空气准一维非平衡喷管流动数值模拟. 考虑5种化学组元(N2, O2, NO, N, O), 其中N2, O2, NO分别有61, 46, 48个振动能级, 将不同振动能级上的粒子视为不同组元, 共157个组元. 对未见公开文献给出速率系数的分子能级跃迁过程, 综合振动能方程松弛时间和其他类似微观过程跃迁速率系数进行折算. 计算结果表明, 喷管喉道前流动接近平衡, 喉道后出现非平衡, 喉道下游不远处发生化学组元质量分数、较低振动能级分子数和表征振动能的振动温度的冻结, N2的振动温度冻结较NO和O2早, 冻结值也更高; 对振动能级跃迁起主导作用的微观机制是平动−振动能量交换(VT)过程, 复合反应生成的分子更多位于中等振动能级; 喷管非平衡和冻结区域分子能级分布偏离振动温度下的玻尔兹曼分布, 高能级出现过分布; 提高驻室总压能够降低喷管流动非平衡程度, 推迟热化学冻结发生.Abstract: Using state-to-state model, quasi-one-dimensional nozzle flow of high-temperature nonequilibrium air is investigated numerically. The five chemical species mixture N2/O2/NO/N/O is considered with 61 bound vibrational levels for N2, 46 for O2, and 48 for NO. Each vibrational state is regarded as a pseudo species, which leads to a total of 157 species for the air mixture. The state-specific transition rate coefficients of some processes, which have no available data, are calculated based on the relaxation time and the rate coefficients of other similar processes. The flow simulation and analysis are made for reservoir temperature from 2000 to 8000 K and pressure from 1 to 20 MPa. The nozzle flow is essentially in equilibrium before the throat, but deviates from equilibrium shortly after the throat. The mass fraction of chemical species, populations of lower energy levels, and vibrational temperatures are frozen in the downstream not far away from the throat. The vibrational temperature of N2 is freezing earlier and has a higher frozen value than that of NO and O2. The process of vibration-translation (VT) energy exchange is predominant for vibrational transition, the recombination reaction generates molecules preferably to middle vibrational levels. Throughout the nonequilibrium and frozen zone, the vibrational population distributions are far from Boltzmann distribution at vibrational temperature, and feature a large overpopulation of the high-lying vibrational levels. Increasing the reservoir pressure could reduce the nonequilibrium to a certain extent and delay the flow thermochemical freezing.
-
引 言
高焓风洞是产生超高速气流的重要地面设备, 能模拟再入或高超声速飞行器周围气体发生的复杂物理化学过程[1]. 风洞驻室条件下发生了振动激发、离解或电离的试验气体在喷管中的膨胀加速伴随着复合反应和振动松弛, 由于气流速度迅速增大、温度密度迅速下降, 来不及达到当地温度压力下热化学平衡态, 喷管流动为典型的非平衡流动. 相对于激波后离解占优的非平衡流, 复合占优的喷管流动表现出不同的非平衡特征. 喷管非平衡流一直受到高温气体动力学领域学者的关注, 数值模拟方面的工作随着热化学模型和数值方法的发展也在持续前进[2-3].
非平衡流数值模拟采用的热化学模型主要有宏观多温度(multi-temperature, MT)模型[4]和精细的态−态(state-to-state, StS)模型[5]两类. 多温度模型采用不同的温度(如振动温度、电子温度)体现分子各模式内能间的热力非平衡, 在化学反应控制温度中引入振动温度和电子温度的影响体现热力−化学耦合[6]. 多温度模型应用广泛但存在一定局限[7-8], 其体现热力−化学耦合的方式隐含着采用了各内能模式的能级分布满足其对应温度下玻尔兹曼分布的假设. 态−态模型是一种直接针对粒子量子态的精细模型, 将位于不同能级的粒子当作独立组元, 构建起从具体量子态反应物到产物的化学反应速率计算模式[9]. 按照对内能能级区分的精细程度, 态−态模型可分为转振态−态(rv-StS)模型[10]、振动(v-StS)态−态模型[11]以及电子能级态−态(e-StS)模型, e-StS一般也被称为碰撞−辐射(collisional-radiative, CR)模型[12]. 由于分子振动能与平动能很容易达到平衡, 而电子能在很高的温度才有明显的激发, 因此目前态−态模型一般特指v-StS, 考虑分子在基态电子能级上不同振动能级间的跃迁, 和各振动能级上分子的化学反应.
态−态模型在20世纪90年代被引入空气动力学, 由于计算量巨大, 早期研究集中于不涉及流动的零维问题, 或是正激波后[13]、边界层[14-16]和准一维喷管[17-18]非平衡流动, 除文献[17-18]针对空气外, 大多针对纯N2或纯O2, 考虑的振动能级数较少, 略去了某些复杂的跃迁过程. 采用态−态模型的早期研究工作促进了对非平衡流动内在机理的认识, 也为改善宏观热化学模型提供了参考依据.
随着量子化学的发展, 态−态模型基元反应速率系数的计算方法和数据近十多年来也有了更新和扩充. 采用态−态模型的非平衡流研究应用了这些最新的数据, 涵盖的微观跃迁机制也更丰富. 正激波后非平衡流动的模拟工作有文献[19-23], 其中文献[19-20, 22]包含了辐射特性, Ninni等[24-25]开展了纯O2和空气的双锥双楔流动模拟.
下面重点介绍近20年来采用态−态模型开展的喷管流动模拟研究. 公开发表的研究工作都采用了准一维假设, 针对高温空气的一般采用了5种化学组元(N2, O2, NO, N和O), 还有针对纯氮(N2, N)或纯氧(O2, O)的.
针对纯氮, Guy等[26]2013年开展了膨胀比为11.6的喷管流动模拟. 考虑了原子碰撞引起的多量子跃迁, 但对分子间碰撞的振动−平动能量交换(VT)过程和振动−振动能量交换(VV)过程, 都仅考虑单量子跃迁. 在驻室T0 = 10000 K, p0 = 0.1 MPa和T0 = 5535 K, p0 = 10 MPa两个条件下的计算结果表明, 喉道下游高振动能级存在明显的“过分布”. 发现限制原子碰撞的VT过程的最大跃迁量子数对N2的摩尔分数和振动能级分布都有明显影响, 当只允许单量子跃迁时, 喷管中N复合程度较低. 作者还提出了一种简化的态−态模型, 将振动能级分为n组, 每组有自己的振动温度, 组内振动能级服从自己振动温度下的玻尔兹曼分布, 形式上与传统的多温度模型相似. T0 = 10000 K算例中, 将能级分为3组即可得到与完整态−态模型相近的振动能级分布. 徐丹等[27]2014年在驻室T0 = 3000 K ~ 10000 K, p0 = 0.1 ~ 1 MPa范围, 对膨胀比为64的喷管采用时间推进法进行了非平衡流动模拟. 考虑了与Guy等[26]相同的跃迁和反应机制, 计算结果表明, 喉道下游振动能级显著偏离振动温度下的玻尔兹曼分布, 复合占优的膨胀流中粒子偏离平衡分布的程度比离解占优的激波后流动中更强.
高温空气喷管非平衡流模拟早期的代表性工作是Colonna等[17-18]1998和1999年对膨胀比为64的喷管进行的. 在驻室T0 = 4000 K ~ 8000 K, p0 = 0.1 MPa条件下, 考虑5种化学组元(N2, O2, NO, N, O), 其中O2有33个, N2有45个振动能级, NO位于振动基态. 计算结果发现喉道下游O2高振动能级分布呈现出一个“平台”分布, 但N2振动能级分布更加复杂, 两者都严重偏离玻尔兹曼分布. 计算结果还表明喉道下游处NO的生成速率远超阿伦尼乌斯公式, 作者指出这起因于置换反应N2(v) + O ↔ NO + N的反应物N2的高振动能级分布远超平衡值, 而这又是因为由N复合反应生成的分子更多位于较高能级. Nagnibeda等[28]2018年进行的空气非平衡喷管流模拟也未考虑NO的振动能级, 同样发现喉道下游N2和O2分子振动能级显著偏离玻尔兹曼分布, 中间振动能级呈“平台”分布. 另外还研究了忽略某些微观过程对喷管出口处温度、压力等的影响.
Gu等[29]2022年针对高焓风洞拉瓦尔喷管开展的空气非平衡流动模拟中, 特别关注了收缩段流动的热化学状态, 研究了驻室总温(1500 K, 4000 K和7000 K)和总压(1, 10和100 MPa)和不同收缩段形状的影响. 结果发现提升总温和总压均能促进喉道处流动趋近热化学平衡, 总压1 MPa条件下态−态模型结果与平衡流结果差异最大: 总温1500 K时为N2的振动温度(低14%), 总温4000 K时为O的质量分数(高8%), 总温7000 K时为N的质量分数(高2%); 总压提高至10 MPa后, 总温4000 K和7000 K条件下态−态模型结果与平衡流结果已无差异, 总温1500 K条件下N2振动温度与平衡流的差别降至5%. 总压100 MPa条件下, 总温1500 K的态−态模型结果也与平衡流完全一致了. 计算结果还表明不同收缩段形状几乎不影响喉道处的热化学状态. 因此作者指出, 对于高焓风洞的典型运行条件, 可以在喷管收缩段采用平衡流假设计算得到喉道处的平衡化学组成、振动能级玻尔兹曼分布, 然后作为采用态−态模型进行扩张段非平衡流空间推进计算的起始条件, 这将使扩张段采用态−态模型进行二维流场的计算量可以接受. Gu等[30]2022年还针对典型的膨胀风洞运行条件专门研究了其扩张喷管的非平衡流动, 喷管入口速度8.0 km/s, 压力20 kPa, 温度3500 K, 视为热化学平衡, 总焓约35 MJ/kg. 研究发现将N2和O2的振动−振动−平动过程(VVT)跃迁量子数限制在3以下可以得到准确结果, 计算得到的喷管出口静压与实验测量值一致. 和双温度模型的结果对比表明, 态−态模型给出的扩张管中气流的热化学松弛速率更慢. 作者以不同热化学模型得出的喷管出口参数作为自由流参数, 进行试验段模型激波后流动模拟, 据此研究自由流热化学状态对模型流场的影响. 发现重要的影响发生在紧靠激波的区域, 当激波后振动激发和离解开始进行并趋于平衡后影响减小, 直至消失. 对自由流热化学特性敏感的激波后气流特性是振动温度、NO和O2的摩尔分数, 它们是影响辐射特性的重要物理量, 因此, 不同自由流条件(分别由态−态模型和双温度模型确定)下, 激波后NO的中红外辐射发射系数峰值差别达到20%, O2的Schumann-Runge带发射系数峰值差别达到50%. 不过将自由流分子振动能级分布设为态−态模型得到的振动温度下的玻尔兹曼分布, 和直接使用态−态模型给出的非平衡分布得到的激波后特性一致.
上述喷管流动研究中, 早期工作采用的跃迁机制和速率系数计算不够完善, 例如文献[18]对NO 参与的反应仅考虑了N2(v) + O ↔ NO + N, 且未考虑NO的振动能级, 速率系数也都来自早期的经验公式. 那之后态−态模型速率系数计算有了较大进展. Adamovich等[31]采用半经典无微扰的强迫谐振子(forced harmonic oscillator, FHO)模型计算分析了N2(v) + N2(w), O2(v) + O2(w), O2(v) + N2(w)的VT过程与VV过程引起的多量子跃迁速率系数, 并根据渐近分析推导了速率系数计算的解析公式, Silva等[32-34]基于FHO模型计算了分子振动能级上的离解反应速率系数[32], 并进一步完善了5化学组元空气中分子−分子间的VT过程[33-34].
FHO模型仍然属于半经验性模型, 近年来, 通过量子化学计算确定态−态速率系数的办法得到越来越多的关注[35-38]. 最准确的速率系数确定方法是耦合分子平动、转动和振动量子态的量子力学(quantum mechanics method, QM)[36]方法, 以及仅对转动−振动或仅对振动做量子力学处理的混合量子经典(mixed quantum classical method, MQC)[38]方法, 但是这些方法计算量很大, 如对于N2−N2的VV过程, 所有可能的跃迁过程数量达到105个, 目前还很难给出完整多量子跃迁过程的速率系数[39].
为了克服经典力学在量子效应方面的局限性和QM方法耗时长的问题, 对轨道始末内能态用量子力学描述但核运动通过Hamilton方程求解的准经典轨道(quasi classical trajectory, QCT)方法得到了国内外研究者的广泛关注, Esposito等[40-44]计算了O2(v) + O[40]和N2(v) + N[41]的VT过程及离解反应的速率系数, O2(v) + N[42-43]和N2(v) + O[44]的VT过程、离解反应及置换反应的速率系数; Fangman等[45]计算了N2(v) + N2与N2(v) + N的VT过程与离解反应速率系数.
虽然文献[29-30]基本采用了新发展的态−态模型跃迁过程和化学反应速率系数, 但是仍有不少VT过程速率系数采用了宏观振动松弛时间折合式, 特别是较为重要的N2(v) + O和O2(v) + N的VT过程, 而这两个过程已有MQC[38]或QCT[42-43]计算结果. 本文在此基础上, 整合最新的速率系数, 对少量仍未见公开文献给出速率系数的VT过程提出处理方法, 构建更加完善的高温空气态−态模型反应机制.
前述采用态−态模型的高温空气非平衡喷管流动模拟, 文献[18]主要关注NO总生成速率, 文献[29]对高焓风洞喷管重点关注收缩段和喉道处的热化学状态, 文献[30]则是针对膨胀风洞的扩张喷管. 本文基于最新的5化学组元空气态−态模型速率系数, 在驻室总温为2000 K ~ 8000 K, 总压为1 ~ 20 MPa条件下, 开展高焓风洞非平衡喷管流模拟. 详细研究流动中粒子能级分布的非平衡演化特点、起支配作用的跃迁和反应机制, 分析驻室总温总压条件对喷管热化学非平衡特性的影响规律及内在机理.
1. 热化学模型与流动控制方程
1.1 态−态模型微观过程与速率系数
对包含5种化学组元(N2, O2, NO, N, O)的高温空气, 采用Lopez等[46]通过量子化学计算得到的分子电子基态上的振动能级数及能值. N2有61个, O2有46个和NO有48个振动能级, 相应数据从STELLAR (state-to-state elementary rates)数据库[34]中获取.
表1给出了振动能级跃迁过程及其正向反应速率系数计算模型, 所有VT与VV过程均包含多量子跃迁. 表中部分过程的速率系数直接引用了文献数据, 对无公开文献数据或现有计算方法不适用于多量子跃迁的过程, 本文根据现有模型或改进方法计算, 标记为“present”. 式中v和w分别表示跃迁初态和末态振动能级. 对跃迁过程1 ~ 7, Adamovich等[31, 47]指出, 在200 K < T < 8000 K范围内, 可以将双原子分子间的VVT过程分解为VT过程和VV过程(满足v1 + v2 = w1 + w2)以简化计算, 因此本文根据FHO模型给出的解析公式计算N2与N2, O2, O2与N2, O2间的VV, VT过程速率系数, 其中FHO模型参数来自文献[34]. N2和O2各自与原子间的VT过程8 ~ 10采用Esposito等[40-43]应用QCT方法计算的速率系数: N2(v) + N[41], O2(v) + O[40]以及O2(v) + N[42-43], 过程11[N2(v) + O]则采用Hong等[38-39]应用MQC方法计算的速率系数, 其中叠加了振动−电子能交换(VE)过程的贡献. 过程14 ~ 17的速率系数采用Oblapenko等[48]结合实验数据和QCT方法修正的FHO模型计算.
表 1 振动能级跃迁过程Table 1. Vibrational processes involved in the StS modelNo. Processes Model Ref. 1 N2(v1) + N2(v2) ↔ N2(w1) + N2(w2) FHO present 2 O2(v1) + O2(v2) ↔ O2(w1) + O2(w2) FHO present 3 N2(v1) + O2(v2) ↔ N2(w1) + O2(w2) FHO present 4 N2(v) + N2↔ N2(w) + N2 FHO present 5 O2(v) + O2 ↔ O2(w) + O2 FHO present 6 N2(v) + O2↔ N2(w) + O2 FHO present 7 O2(v) + N2 ↔ O2(w) + N2 FHO present 8 N2(v) + N ↔ N2(w) + N QCT [41] 9 O2(v) + N ↔ O2(w) + N QCT [42-43] 10 O2(v) + O ↔ O2(w) + O QCT [40] 11 N2(v) + O ↔ N2(w) + O MQC+
RT_Ratio[38], present 12 N2(v) + NO ↔ N2(w) + NO RT_Ratio present 13 O2(v) + NO↔ O2(w) + NO RT_Ratio present 14 NO(v) + N2 ↔ NO(w) + N2 FHO present 15 NO(v) + O2 ↔ NO(w) + O2 FHO present 16 NO(v) + NO ↔ NO(w) + NO FHO present 17 NO(v) + N ↔ NO(w) + N FHO present 18 NO(v) + O ↔ NO(w) + O RT_Ratio present 对VT过程11 ~ 13和18, 无直接可用的速率系数数据或当前模型不适用. 常见的做法[19-20, 30, 39]是: 根据分子此类碰撞过程对应的振动能变化率方程中的松弛时间$ \tau _{{\text{AB}} - X}^{{{\mathrm{vib}}} } $, 折合出由v能级跃迁到w能级的速率系数$ k_{{\text{AB}} - X}^{{\text{VT,}}v \to w} $
$$ k_{{\text{AB}} - X}^{{\text{VT,}}v \to w} = \frac{1}{{{n_X}\tau _{{\text{AB}} - X}^{{{\mathrm{vib}}} }}}\dfrac{{\exp \left( { - \dfrac{{\varepsilon _{{\text{AB(}}w{\text{)}}}^{{{\mathrm{vib}}} }}}{{{k_{\text{B}}}T}}} \right)}}{{Q_{{\text{AB}}}^{{{{\mathrm{vib}}}}}}} $$ (1) 式中, AB表示分子, X表示参与碰撞的粒子, X ∈ (N2, O2, NO, N, O), nX为其数密度, $ \varepsilon _{{\text{AB(}}w{\text{)}}}^{{{\mathrm{vib}}} } $为w能级上一个分子的振动能, $ Q_{{\text{AB}}}^{ {\text{vib}}} $为配分函数, kB为玻尔兹曼常数. 该方法通过振动能松弛时间直接折合出跃迁速率系数, 这里将其记为RT_D (from relaxation time_directly), 它具一定合理性, 且空气5种化学组元中分子与其余各粒子间碰撞过程的振动松弛时间数据完备. 但需要注意的是, 根据各能级粒子数变化率导出振动能变化率方程时, 只考虑了单量子跃迁, 对跃迁速率系数采用了递推规律和细致平衡原理
$$ {k_{v,v - 1}} = v{k_{1,0}} $$ (2) $$ {k_{v - 1,v}} = {k_{v,v - 1}}\exp \left( { - \frac{{h\nu }}{{{k_{B} }T}}} \right) $$ (3) 式中, k1,0为1能级跃迁至基态的速率系数, h为普朗克常数, $ \nu $为分子的基振频率. 那么式(1)采用振动松弛时间折算跃迁率系数, 也隐含了只考虑单量子跃迁. 另外, 式(1)没有区分跃迁初态v = w – 1和v = w + 1的速率系数, 这不合理. 其实式(3)通过细致平衡原理体现了由高或低一能级跃迁至w能级的速率系数的差异.
在振动松弛时间基础上进一步考虑多量子跃迁, 并对具体不同跃迁初态能级详细分配速率系数比较复杂. 笔者认为可根据宏观振动松弛时间能体现微观跃迁速率的内涵, 通过两个类似跃迁过程的松弛时间比值推算其能级跃迁速率系数比值, 从而由已知的微观过程速率系数折合出另一微观过程速率系数. 谐振子、单量子跃迁假设下得出的振动能变化率方程中松弛时间$ {\tau ^{{{\mathrm{vib}}} }} $与速率系数k1,0的关系为[49]
$$ {\tau ^{{{\mathrm{vib}}} }} = \frac{1}{{{k_{1,0}}\left[ {1 - \exp \left( { - \dfrac{{h\nu }}{{{k_{B} }T}}} \right)} \right]}} $$ (4) 结合式(4)和式(2)、式(3)可知, 振动松弛时间与跃迁速率系数成反比. 对于AB(v) + X的VT过程, 若已知跃迁速率系数的X粒子有s个, 则对于未知速率系数的AB(v) + Y过程可采用下式折算
$$ k_{{\text{AB}} - Y}^{{\text{VT,}}v \to w} = \sum\limits_{i=1}^s {\left( {\dfrac{1}{s}\frac{{\tau _{{\text{AB}} - {X_i}}^{{{\mathrm{vib}}} }}}{{\tau _{{\text{AB}} - Y}^{{{\mathrm{vib}}} }}}k_{{\text{AB}} - {X_i}}^{{\text{VT,}}v \to w}} \right)} $$ (5) 式中, Xi为采用FHO或QCT等方法已得到其与AB(v)间VT过程跃迁速率系数的粒子. 将这种综合利用振动松弛时间比值和其他过程跃迁速率系数的处理方法记为RT_Ratio (from relaxation time ratio). 相较于直接应用振动松弛时间的RT_D方法(式(1)), 该方法能够利用类似的微观过程的跃迁速率系数信息, 可涵盖单量子或多量子跃迁, 并且包含了非谐振子假设等的影响.
本文采用RT_Ratio方法补全了空气5种化学组元VT过程的能级跃迁速率系数. 对于除MQC方法计算得到部分跃迁速率外的过程11和整个过程12, 由过程4, 6, 8折合; 对于过程13, 由过程5, 7, 9, 10折合; 对于VT过程18, 由过程14 ~ 17折合得到. 振动松弛时间采用Park[6]修正的Millikan和White振动松弛模型计算.
以过程9(O2(v) + N)为例检验RT_Ratio方法的合理性. 将直接采用振动松弛时间折算的RT_D方法、本文提出的RT_Ratio方法与QCT方法的计算结果进行对比, 见图1, 其中采用RT_Ratio计算时, 参与折合的VT过程有5, 7和10. RT_D方法对应单量子跃迁, RT_Ratio与QCT方法列出了O2(0, 5, 12) + N → O2(w) + N的多量子跃迁速率系数, 横坐标为跃迁末态能级w. RT_Ratio方法得到的速率系数与QCT计算得到的速率系数趋势吻合, 当初态为基态, 跃迁量子数较小时, 两者一致性较好, 随着跃迁量子数增加, 差别增大, 但此时跃迁速率系数本身已经很小, 故影响有限; 当初态能级较高(如v = 12)时, RT_Ratio方法仍能得到和QCT方法相近的多量子跃迁结果.
![]()
图 1 采用不同方法计算得到的O2 + N的多量子跃迁VT过程速率系数 (T = 5000 K)Figure 1. Multi-quantum VT transition rates for O2 + N by different methods (T = 5000 K)表2给出了离解复合和置换反应及其正向反应速率系数计算模型和来源文献. 采用DRM和DRA分别表示分子或原子参与碰撞的离解复合反应, ZE表示置换反应. 反应1 ~ 6为3种分子的DRM, 参与碰撞的分子为N2或O2, 来自Silva等[34]应用推广FHO模型计算的速率系数. 反应7 ~ 10为N2或O2的DRA, 来自Esposito等应用QCT方法的计算结果: N2(v) + N[41], O2(v) + N[42-43], N2(v) + O[44]和O2(v) + O[40]. 置换反应11和12采用了Bose等[50-51]的QCT方法计算结果. 反应1 ~ 7和10 ~ 12的速率系数取自STELLAR数据库[34], 反应8和9的速率系数采用文献提供的插值程序计算得到.
表 2 离解复合反应和置换反应机制Table 2. Dissociation/recombination and Zeldovich exchange reactions involved in the StS modelNo. Reaction Model Ref. 1 N2(v) + N2 ↔ 2 N + N2 FHO [34] 2 O2(v) + N2 ↔ 2 O + N2 FHO [34] 3 NO(v) + N2 ↔ N + O + N2 FHO [34] 4 N2(v) + O2 ↔ 2 N + O2 FHO [34] 5 O2(v) + O2 ↔ 2 O + O2 FHO [34] 6 NO(v) + O2 ↔ N + O + O2 FHO [34] 7 N2(v) + N ↔ 3 N QCT [41] 8 O2(v) + N ↔ 2 O + N QCT [42-43] 9 N2(v) + O ↔ 2 N + O QCT [44] 10 O2(v) + O ↔ 3 O QCT [40] 11 N2(v) + O ↔ NO(w) + N QCT [50] 12 O2(v) + N ↔ NO(w) + O QCT [51] 13 N2(v) + NO↔ 2 N + NO reaction 1 [6, 20] 14 O2(v) + NO↔ 2 O + NO reaction 2 [6, 20] 15 NO(v) + NO↔ N + O + NO 20 × (reaction 3) [6, 20] 16 NO(v) + N ↔ 2 N + O 20 × (reaction 3) [6, 20] 17 NO(v) + O↔ N + 2 O 20 × (reaction 3) [6, 20] 反应13 ~ 15分别是NO作为参与碰撞粒子的3种分子的DRM, 反应16和17是NO的DRA, 均无可用速率系数, 本文按照宏观化学反应模型中应用的催化物系数[6]确定, 例如令N2(v) + NO离解反应的速率系数和N2(v) + N2的相等, 令NO(v) + N离解反应的速率系数等于NO(v) + N2的20倍.
表1和表2中微观过程的逆过程速率系数由细致平衡原理得到[22, 52-53].
1.2 粒子数密度变化率方程
记位于振动能级v的分子AB(v), 其数密度nAB(v)的变化率为
$$ \begin{split} & \frac{{{\text{d}}{n_{{\text{AB}}(v)}}}}{{{\text{d}}t}} = {\left\{ {\frac{{{\text{d}}{n_{{\text{AB}}(v)}}}}{{{\text{d}}t}}} \right\}_{{\text{VT}}}} + {\left\{ {\frac{{{\text{d}}{n_{{\text{AB}}(v)}}}}{{{\text{d}}t}}} \right\}_{{\text{VV}}}}+ \\ &\qquad {\left\{ {\frac{{{\text{d}}{n_{{\text{AB}}(v)}}}}{{{\text{d}}t}}} \right\}_{{\text{DR}}}} + {\left\{ {\frac{{{\text{d}}{n_{{\text{AB}}(v)}}}}{{{\text{d}}t}}} \right\}_{{\text{ZE}}}} \end{split} $$ (6) 式中等号右端前两项分别为来源于VT过程、VV过程引起的振动能级跃迁的贡献, 后两项则为来源于离解复合反应、置换反应的贡献. 下面给出AB(v)在各过程中发生的变化率.
由VT过程引起的nAB(v)变化率为
$$ \begin{split} &{{\left\{ \frac{\text{d}{{n}_{\text{AB}(v)}}}{\text{d}t} \right\}}_{\text{VT}}}=\\ &\qquad \sum\limits_{X}{\sum\limits_{\begin{smallmatrix} w=0 \\ w\ne v \end{smallmatrix}}^{v_{\max }^{\text{AB}}}{\left( k_{\text{AB}-X}^{\text{VT,}w\to v}{{n}_{\text{AB}(w)}}{{n}_{X}} \right.}}\left. -k_{\text{AB}-X}^{\text{VT,}v\to w}{{n}_{\text{AB}(v)}}{{n}_{X}} \right) \end{split}$$ (7) 式中, $v_{\max }^{{\text{AB}}}$为最高振动能级.
对VV过程, 令v1 = v. 当碰撞对为同种分子时, w1 = v1或w2 = v1的VV过程实际上是VT过程, 在式(7)中已考虑. 故由VV过程引起的nAB(v)变化率为
$$ \begin{split} & \left\{ \frac{\text{d}{{n}_{\text{AB}(v)}}}{\text{d}t} \right\}_{\operatorname{VV}}^{{{v}_{1}}=v}=\\ &\qquad \sum\limits_{\begin{smallmatrix} {{w}_{1}}=0 \\ {{w}_{1}}\ne {{v}_{1}} \end{smallmatrix}}^{v_{\max }^{\text{AB}}}{\sum\limits_{\begin{smallmatrix} {{w}_{2}}={{w}_{1}} \\ {{w}_{2}}\ne {{v}_{1}} \end{smallmatrix}}^{v_{\max }^{\text{AB}}}{\sum\limits_{\begin{smallmatrix} {{v}_{2}}=0 \\ {{v}_{2}}\ne {{v}_{1}} \end{smallmatrix}}^{v_{\max }^{\text{AB}}}{\left( k_{{{w}_{1}}{{w}_{2}}\to {{v}_{1}}{{v}_{2}}}^{\text{VV,S}}{{n}_{\text{AB}({{w}_{1}})}}{{n}_{\text{AB}({{w}_{2}})}} \right.}}} -\\ &\qquad \left. k_{{{v}_{1}}{{v}_{2}}\to {{w}_{1}}{{w}_{2}}}^{\text{VV,S}}{{n}_{\text{AB}({{v}_{1}})}}{{n}_{\text{AB}({{v}_{2}})}} \right)+ \\ &\qquad \sum\limits_{\begin{smallmatrix} {{w}_{1}}=0 \\ {{w}_{1}}\ne {{v}_{1}} \end{smallmatrix}}^{v_{\max }^{\text{AB}}}{\sum\limits_{\begin{smallmatrix} {{w}_{2}}={{w}_{1}} \\ {{w}_{2}}\ne {{v}_{1}} \end{smallmatrix}}^{v_{\max }^{\text{AB}}}{2\left( k_{{{w}_{1}}{{w}_{2}}\to {{v}_{1}}{{v}_{1}}}^{\text{VV,S}}{{n}_{\text{AB}({{w}_{1}})}}{{n}_{\text{AB}({{w}_{2}})}} \right.}}-\\ &\qquad \left. k_{{{v}_{1}}{{v}_{1}}\to {{w}_{1}}{{w}_{2}}}^{\text{VV,S}}n_{\text{AB}({{v}_{1}})}^{2} \right)+ \\ &\qquad \sum\limits_{\begin{smallmatrix} {{w}_{1}}=0 \\ {{w}_{1}}\ne {{v}_{1}} \end{smallmatrix}}^{v_{\max }^{\text{AB}}}{\sum\limits_{\begin{smallmatrix} {{w}_{2}}=0 \\ {{w}_{2}}\ne {{v}_{2}} \end{smallmatrix}}^{v_{\max }^{\text{CD}}}{\sum\limits_{{{v}_{2}}=0}^{v_{\max }^{\text{CD}}}{\left( k_{{{w}_{1}}{{w}_{2}}\to {{v}_{1}}{{v}_{2}}}^{\text{VV,D}}{{n}_{\text{AB}({{w}_{1}})}}{{n}_{\text{CD}({{w}_{2}})}} \right.}}}-\\ &\qquad \left. k_{{{v}_{1}}{{v}_{2}}\to {{w}_{1}}{{w}_{2}}}^{\text{VV,D}}{{n}_{\text{AB}({{v}_{1}})}}{{n}_{\text{CD}({{v}_{2}})}} \right) \end{split} $$ (8) 式中CD表示与AB不同的分子, $ k_{{v_1}{v_2} \to {w_1}{w_2}}^{{\text{VV,S}}} $和$ k_{{w_1}{w_2} \to {v_1}{v_2}}^{{\text{VV,S}}} $分别为同种分子VV过程的正逆速率系数, $ k_{{v_1}{v_2} \to {w_1}{w_2}}^{{\text{VV,D}}} $和$ k_{{w_1}{w_2} \to {v_1}{v_2}}^{{\text{VV,D}}} $分别为异种分子VV过程的正逆速率系数.
由离解复合反应引起的nAB(v)变化率为
$$ {\left\{ {\frac{{{\text{d}}{n_{{\text{AB}}(v)}}}}{{{\text{d}}t}}} \right\}_{{\text{DR}}}} = \sum\limits_X {\left( {k_{{\text{AB}}(v) - X}^{\text{R}}{n_{\text{A}}}{n_{\text{B}}}{n_X}} \right.} \left. { - k_{{\text{AB}}(v) - X}^{\text{D}}{n_{{\text{AB}}(v)}}{n_X}} \right) $$ (9) 式中, $ k_{{\text{AB}}(v) - X}^{\text{R}} $为AB和X碰撞复合反应的速率系数, $ k_{{\text{AB}}(v) - X}^{\text{D}} $为离解反应速率系数.
以N2(v) + O ↔ NO(w) + N (表2中反应11)为例给出置换反应引起的分子振动能级数密度变化率. 该反应引起的NO(w)数密度变化率为
$$ \begin{split} & {\left\{ {\frac{{{\text{d}}{n_{{\text{NO}}(w)}}}}{{{\text{d}}t}}} \right\}_{{\text{ZE}}}} = \sum\limits_{v = 0}^{v_{\max }^{{{\text{N}}_{\text{2}}}}} {\left( {k_{{{\text{N}}_2}{\text{(}}v{\text{)}} - {\text{O}}}^{{\text{ZE}}}{n_{{{\text{N}}_2}{\text{(}}v{\text{)}}}}{n_{\text{O}}}} \right.}- \\ &\qquad \left. {{\text{ }} k_{{\text{NO(}}w{\text{)}} - {\text{N}}}^{{\text{ZE}}}{n_{{\text{NO(}}w{\text{)}}}}{n_{\text{N}}}} \right) \end{split} $$ (10) 式中, $ k_{{{\text{N}}_2}{\text{(}}v{\text{)}} - {\text{O}}}^{{\text{ZE}}} $和$ k_{{\text{NO(}}w{\text{)}} - {\text{N}}}^{{\text{ZE}}} $分别为正逆反应速率系数. 该反应引起的N2(v)数密度变化率为
$$ \begin{split} & {\left\{ {\frac{{{\text{d}}{n_{{{\text{N}}_2}(v)}}}}{{{\text{d}}t}}} \right\}_{{\text{ZE}}}} = \sum\limits_{w = 0}^{v_{\max }^{{\text{NO}}}} {\left( {k_{{\text{NO(}}w{\text{)}} - {\text{N}}}^{{\text{ZE}}}{n_{{\text{NO(}}w{\text{)}}}}{n_{\text{N}}}} \right.}- \\ &\qquad {\text{ }}\left. { k_{{{\text{N}}_2}{\text{(}}v{\text{)}} - {\text{O}}}^{{\text{ZE}}}{n_{{{\text{N}}_2}{\text{(}}v{\text{)}}}}{n_{\text{O}}}} \right) \end{split} $$ (11) 置换反应O2(v) + N ↔ NO(w) + O引起的O2(v)和NO(w)的变化率类似可得.
N和O的粒子数密度变化率仅涉及DR和ZE反应, 可按照式(9)和式(10)、式(11)类似得出.
1.3 流动控制方程
定常准一维喷管流动的控制方程如下.
组元连续方程
$$ \frac{{\text{d}}}{{{\text{d}}x}}\left( {\rho u{c_i}A} \right) = A{\dot \omega _i},\quad i = 1,2, \cdots ,{N_s} $$ (12) 动量方程
$$ \frac{{{\text{d}}p}}{{{\text{d}}x}} + \rho u\frac{{{\text{d}}u}}{{{\text{d}}x}} = 0 $$ (13) 能量方程
$$ p\frac{{{\text{d}}u}}{{{\text{d}}x}} + \rho u\frac{{{\text{d}}e}}{{{\text{d}}x}} + \frac{{pu}}{A}\frac{{{\text{d}}A}}{{{\text{d}}x}} = 0 $$ (14) 式中, x为喷管轴线坐标, A是喷管横截面面积, ρ为混合气体的密度, u为速度, p为压力.
态−态模型下, 位于不同振动能级的分子视为独立组元, 分子组元有61个N2(v)、46个O2(v)和48个NO(v), 原子组元为N和O, 组元总数Ns = 157, 故有157个组元连续方程. 式(12)中ci为第i组元的质量分数, $ {\dot \omega _i} $为单位时间单位体积中生成的i组元的质量, 以分子组元AB(v)为例, 其质量生成项为
$$ \dot \omega _{{\text{AB}}(v)}^{} = \frac{{{\mathcal{M}_{{\text{AB}}}}}}{{{N_{\text{A}}}}}\frac{{{\text{d}}{n_{{\text{AB}}(v)}}}}{{{\text{d}}t}} $$ (15) 式中, $ {\mathcal{M}_{{\text{AB}}}} $为摩尔质量, NA为阿伏伽德罗常数.
能量方程中e为单位质量气体的内能
$$ \begin{split} & e = \sum\limits_{i = 1}^{{N_s}} {{c_i}{e_i}} = \sum\limits_{i = {{\mathrm{mol}}} . } {\left( {\frac{5}{2}{c_i}{R_i}T + {c_i}e_i^{{{\mathrm{vib}}} }} \right)}+ \\ &\qquad {\text{ }} \sum\limits_{i = {{\mathrm{atom}}} } {\frac{3}{2}{c_i}{R_i}T} + \sum\limits_{i = 1}^{{N_s}} {{c_i}h_i^{0{\text{f}}}} \end{split} $$ (16) 式中, i = mol.代表分子组元, i = atom代表原子组元. $ {e_i} $为单位质量i组元的内能, $ h_i^{0{\text{f}}} $为单位质量i组元的生成焓, Ri为其气体常数, $ e_i^{{{\mathrm{vib}}} } $为单位质量分子组元的振动能, 以位于v能级的分子AB(v)为例, 有
$$ e_i^{{{\mathrm{vib}}} } = \frac{{\varepsilon _{{\text{AB(}}v{\text{)}}}^{{{\mathrm{vib}}} }{N_{\text{A}}}}}{{{\mathcal{M}_{{\text{AB}}}}}} $$ (17) 将式(16)中的平动、转动能采用混合气体的平动−转动比定容热容$ {\bar c_{V,tr}} $表示后, 有
$$ e = {\bar c_{V,tr}}T + \sum\limits_{i = {{\mathrm{mol}}} . } {{c_i}e_i^{{{\mathrm{vib}}} }} + \sum\limits_{i = 1}^{{N_s}} {{c_i}h_i^{0{\text{f}}}} $$ (18) 状态方程为
$$ p = \rho \bar RT $$ (19) 式中$ \bar R $为混合气体的气体常数.
$$ \bar R = \sum\limits_{i = 1}^{N{ _s}} {{c_i}{R_i}} $$ (20) 定常流动中混合气体的密度ρ根据确定的质量流率qm = ρuA由速度和截面积确定.
态−态模型中气体的振动能可直接由各分子组元能级粒子数分布得到. 为直观表征各化学组元振动能与平动能间的非平衡程度, 根据振动能折算等效振动温度
$$ \sum\limits_{v = 0}^{v_{\max }^{{\text{AB}}}} {\frac{{{n_{{\text{AB(}}v{\text{)}}}}}}{{{n_{{\text{AB}}}}}}} \varepsilon _{{\text{AB(}}v{\text{)}}}^{{{\mathrm{vib}}} } = \frac{{\displaystyle\sum\limits_{v = 0}^{v_{\max }^{{\text{AB}}}} {\varepsilon _{{\text{AB(}}v{\text{)}}}^{{{\mathrm{vib}}} }\exp \left( { - \frac{{\varepsilon _{{\text{AB(}}v{\text{)}}}^{{\text{vib}}}}}{{{k_{\text{B}}}{T_{{\text{vib,AB}}}}}}} \right)} }}{{Q_{{\text{AB}}}^{{{\mathrm{vib}}} }}} $$ (21) 式中Tvib,AB即AB分子的等效振动温度. 后面将分析喷管流动中气流振动温度的变化, 并详细考察等效振动温度下的玻尔兹曼分布与态−态模型直接求得的粒子非平衡分布间差异.
2. 数值方法及边界条件处理
2.1 数值方法
控制方程中式(12) ~ 式(14)构成了关于组元质量分数ci、速度u和内能e的常微分方程组, 内能与温度、压力与密度间的关系分别由式(18)和式(19)确定, 另由质量流率确定混合气体的密度. 参照文献[18, 29]的做法, 选定ci, u和T为求解变量, 将常微分方程组改写为
$$ \frac{{{\text{d}}{c_i}}}{{{\text{d}}x}} = \frac{{{{\dot \omega }_i}}}{{\rho u}} $$ (22) $$ \frac{{{\text{d}}u}}{{{\text{d}}x}} = \frac{{{{\bar c}_{V,tr}}T\displaystyle\sum\limits_{i = 1}^{{N_s}} {{R_i}{{\dot \omega }_i}} - \bar R\sum\limits_{i = 1}^{{N_s}} {{e_i}{{\dot \omega }_i}} - \frac{{up\left( {{{\bar c}_{V,tr}} + \bar R} \right)}}{A}\frac{{{\text{d}}A}}{{{\text{d}}x}}}}{{\left( {p - \rho {u^2}} \right){{\bar c}_{V,tr}} + p\bar R}} $$ (23) $$ \frac{{{\text{d}}T}}{{{\text{d}}x}} = \frac{{\left( {u - \dfrac{p}{{\rho u}}} \right)\displaystyle\sum\limits_{i = 1}^{{N_s}} {{e_i}{{\dot \omega }_i}} + \frac{{{u^2}p}}{A}\frac{{{\text{d}}A}}{{{\text{d}}x}} - \frac{{pT}}{{\rho u}}\sum\limits_{i = 1}^{{N_s}} {{R_i}{{\dot \omega }_i}} }}{{\left( {p - \rho {u^2}} \right){{\bar c}_{V,tr}} + p\bar R}} $$ (24) 由于各类跃迁机制的速率差别大, 导致常微分方程组的数值求解存在严重刚性, 本文使用了变系数的常微分方程求解器DVODE[54].
2.2 入口边界条件处理
从喷管入口开始空间推进求解控制方程时, 需要给出入口气流速度. 对等熵的准一维喷管流动, 可根据喉道达到声速的条件由驻室总温总压确定质量流率qm = ρuA, 进一步获得入口密度、速度等其他气流参数. 但非平衡流动不等熵, 无法直接给出质量流率, 实际求解时需要先假设流量初值, 进行迭代计算. 另外流动在喉道附近存在的鞍点奇性问题, 使得从亚声速到超声速的积分变得困难.
下面以式(24)为代表, 分析喉道附近的奇性问题, 并据此给出迭代确定流率的方法. 式(24)分母为
$$ \begin{split} & \left( {p - \rho {u^2}} \right){{\bar c}_{V,tr}} + p\bar R = \rho \left( {\frac{{{{\bar c}_{p,tr}}}}{{{{\bar c}_{V,tr}}}}\bar {R}T - {u^2}} \right){{\bar c}_{V,tr}} = \\ &\qquad \rho \left( {a_f^2 - {u^2}} \right){{\bar c}_{V,tr}} \end{split} $$ (25) 式中$ {\bar c_{p,tr}} $为平动−转动比定压热容, af为冻结声速. 将式(25)代入式(24), 整理得到
$$ \frac{{{\text{d}}T}}{{{\text{d}}x}} = \frac{{\dfrac{1}{{{q_m}}}\left[ {\left( {u{q_m} - pA} \right)\displaystyle\sum\limits_{i = 1}^{{N_s}} {{e_i}{{\dot \omega }_i}} - pTA\sum\limits_{i = 1}^{{N_s}} {{R_i}{{\dot \omega }_i}} } \right] + \dfrac{{{u^2}p}}{A}\dfrac{{{\text{d}}A}}{{{\text{d}}x}}}}{{\rho {{\bar c}_{V,tr}}\left( {a_f^2 - {u^2}} \right)}} $$ (26) 当u趋近于af时, dT/dx的分母趋近于0, 欲使dT/dx有意义, 其分子也必须趋近于0, 这就是确定质量流率qm的条件.
以文献[27]中的喷管为例说明迭代确定qm的方法. 在驻室T0 = 5000 K, p0 = 1 MPa条件下, 设定某一qm初值, 根据总温总压条件联立确定入口气流速度、温度和密度等, 开始空间推进求解. 图2给出了不同质量流率对应的喷管喉道附近温度变化情况. qm = 28.94918 kg/s是迭代得出的使dT/dx分子趋近于0的qm值, 对应图中黑实线. 入口质量流率小于该值时, u不会趋近于af, 喷管内为全亚声速流动; 而质量流率大于该值时, 则dT/dx的分母趋于零时分子不趋近于零, dT/dx趋于−∞.
![]()
图 2 不同质量流率对应的喉道附近温度分布Figure 2. Computed temperature profile around nozzle throat for different values of qm计算时, 通过判断dT/dx的值用二分法进行迭代以确定质量流率, 文献[29]采用Crank-Nicolson格式, 在3.0 GHz Intel i7 CPU平台上需迭代1 ~ 2 h. 以驻室T0 = 5000 K, p0 = 1 MPa条件为例, 本文程序在单核2.3 GHz Intel i7 CPU平台上运行时, 迭代5 min左右即可得到足够精度的质量流率, 进而用30 s左右求解完整的喷管流动. 值得一提的是, 包含多量子VV跃迁使计算耗时大幅增加, 如表3所示.
表 3 喷管流动求解与迭代确定流量耗时Table 3. Execution time for simulation nozzle flow and iteration for mass flow rateModel Iterativetime/s Flow calculation time/s full VV 298 32 single quantum VV 60 6 2.3 程序验证
为验证本文反应速率系数和喷管流动计算程序, 对文献[29]中驻室T0 = 4000 K, p0 = 1 MPa条件下的锥形喷管流动进行模拟. 喷管喉道半径为15 mm, 收缩段半锥角为10°, 扩张段半锥角为6°, 在喉道处光滑过渡, 但文献中未给出具体过渡形式, 本文采用抛物线过渡. 喷管总长1.5 m, 喉道位于x = 0 处, 半径沿轴线变化设定为
$$ \left. \begin{aligned} & {r_1} = - 0.176x + 0.011\;128,\quad - 0.4 < x \leqslant - 0.044 \\ & {r_2} = 2{x^2} + 0.015,\quad - 0.044 < x \leqslant 0.026 \\ & {r_3} = 0.105x + 0.013\;622,\quad 0.026 < x \leqslant 1.1 \end{aligned} \right\} $$ (27) 文献[29]既进行了喷管全程的态−态模型(记为full StS)计算, 又进行了入口到喉道按平衡流动、喉道后用态−态模型的计算(记为simplified StS).
图3给出了本文和文献得到的膨胀过程中温度T与O质量分数cO变化. 本文T变化曲线与文献结果基本重合, 出口处相对误差为4.9%; cO介于文献的full StS与simplified StS结果之间, 相较于full StS结果相对差别为1.4%. 可见本文计算结果与文献吻合良好.
![]()
3. 算例与结果分析
在驻室总温2000 K ~ 8000 K、总压1 ~ 20 MPa范围, 对一膨胀比为64的轴对称喷管进行非平衡流动模拟与分析. 喷管半径沿轴线变化曲线为抛物线
$$ r\left( x \right) = 3.5{x^2} - 3.5x + 1 $$ (28) 喷管总长L = 1 m, 喉道位于x = 0.5 m处.
3.1 算例设置
选取总温T0 = 2000 K, 5000 K, 8000 K和总压p0 = 1, 5, 20 MPa共9个驻室条件, 设定驻室的高温空气处于热化学平衡, 氮氧元素质量比为0.76: 0.24. 表4给出了9个条件下驻室的N, O和NO质量分数. 可见总温2000 K条件下驻室基本未发生化学反应; 总温5000 K时, 各总压条件下O2都有大量离解, 而N2基本未离解, 且明显存在NO; 总温8000 K条件下O2已离解完毕, N2在p0 = 1 MPa时约离解一半, 20 MPa时只有少量离解.
表 4 算例驻室条件Table 4. The reservoir conditions tested in this workCase T0/K p0/MPa cO,res cN,res cNO,res 1 2000 1 6.2 × 10−5 0 2.4 × 10−3 2 2000 5 2.8 × 10−5 0 2.4 × 10−3 3 2000 20 1.4 × 10−5 0 2.4 × 10−3 4 5000 1 0.203 4.3 × 10−3 0.032 5 5000 5 0.148 1.9 × 10−3 0.053 6 5000 20 0.093 9.0 × 10−4 0.069 7 8000 1 0.240 0.355 7.5 × 10−3 8 8000 5 0.232 0.177 0.021 9 8000 20 0.216 0.089 0.044 本文首先在总压1 MPa条件下, 对3个总温条件选取代表性流动热化学特性进行分析. 对总温2000 K条件, 重点分析气流膨胀加速过程中O2振动能级分布的演化; 对总温5000 K条件, 同时分析膨胀过程中的O复合反应和O2, NO振动能级跃迁, 考察热力−化学耦合特点; 对总温8000 K条件, 则重点分析N的复合反应和N2分子的振动能级跃迁及二者之间的联系和影响. 在上述分析基础上, 选取总温5000 K条件, 分析总压对流动热化学非平衡特性的影响.
3.2 T0 = 2000 K条件下的流动特征分析
该条件下几乎不发生化学反应, 驻室cO仅为10−5量级. 伴随气流膨胀加速的就是N2与O2的热力非平衡过程. 图4给出了平动温度和O2, N2的振动温度分布. 平动温度与振动温度在x/L = 0.4之前都基本保持不变, 此处气流速度仅335 m/s, 喉道处(x/L = 0.5)速度达到812 m/s, 喉道后速度迅速增大, 对应平动温度迅速降低, 到出口处速度达1937 m/s, 平动温度降为245 K. 振动温度在平动温度降低的过程中与之分离, 在喉道下游一定位置冻结(本文以达到物理量冻结值105%或者95%的位置为冻结位置). N2的热力非平衡程度远高于O2, N2振动温度在x/L = 0.45处开始与平动温度分离, 在x/L = 0.58处趋于冻结值1685 K; O2振动温度在x/L = 0.57处开始与平动温度分离, 在x/L = 0.72附近趋于冻结, 冻结值为1044 K.
![]()
图 4 T0 = 2000 K条件下膨胀过程中平动温度与振动温度变化Figure 4. Translational and vibrational temperature profiles during expansion for T0 = 2000 K驻室条件下约1/3的O2跃迁到了激发态能级. 膨胀过程中高能级粒子向低能级跃迁, 激发态能级粒子数占比逐渐下降, 到出口处约为1/10. 图5给出了各振动能级的O2质量分数变化过程, 基态O2 (v = 0)的质量分数膨胀过程中一直在增大, 而能级v = 1 ~ 10的质量分数一直降低. 能级v > 10的质量分数在x/L = 0.6后出现了明显的上升, 这是因为氧原子复合到了高振动能级(cO从驻室的6.2 × 10−5降低到喷管出口处的3.1 × 10−5), 由于高能级粒子数本身很少, 所以体现明显. 图6给出了若干x/L位置O2的振动能级分布, 在喉道处高振动能级就开始偏离振动温度下的玻尔兹曼分布, 到喉道下游后中间能级偏离得更严重.
![]()
图 5 T0 = 2000 K条件下膨胀过程中O2(v)质量分数变化Figure 5. Mass fraction profiles of O2(v) during expansion for T0 = 2000 K![]()
图 6 T0 = 2000 K条件下不同站位处O2振动能级分布Figure 6. Vibrational population distributions of O2 at different position for T0 = 2000 K由于高振动能级粒子数占比很小, 气体的振动能及对应的振动温度其实取决于较低的几个振动能级. 下面通过观察膨胀过程中N2与O2几个低能级粒子数的演化分析N2与O2振动温度差异的原因. 图7给出了膨胀过程中v = 0, 1, 3, 6能级质量分数的导数与对应分子总质量分数比值的变化. 基态N2与O2的质量分数导数峰值及趋于零的位置均存在差异, (dcN2(0)/dx)/cN2的峰值为0.33, 峰值位置为x/L = 0.49, 在 x/L = 0.6附近趋于零, 而(dcO2(0)/dx)/cO2的峰值为1.17, 峰值位置为x/L = 0.55, 在 x/L = 0.76附近趋于0. 这说明N2向基态跃迁的相对量比O2小, 跃迁速率开始急剧下降对应的温度更高: N2(0)在平动温度降至1200 K后dcN2(0)/dx就趋于0了, 而对O2(0), 这个温度值是约500 K. 图7中各能级质量分数变化主要由VT过程贡献, 且占优机制为N2(1) → N2(0)和O2(1) → O2(0). 在x/L = 0.48处, T = 1790 K, $ k_{{{\text{N}}_2} - {{\text{N}}_2}}^{{\text{VT,}}1 \to 0} $ = 1.5 × 10−15 cm3/s, $ k_{{{\text{N}}_2} - {{\text{O}}_2}}^{{\text{VT,}}1 \to 0} $ = 9.2 × 10−16 cm3/s, $ k_{{{\text{O}}_2} - {{\text{N}}_2}}^{{\text{VT,}}1 \to 0} $ = 1.1 × 10−13 cm3/s, $ k_{{{\text{O}}_2} - {{\text{O}}_2}}^{{\text{VT,}}1 \to 0} $ = 1.1 × 10−13 cm3/s, 在x/L = 0.6处, T = 1200 K, 相应的$ k_{{{\text{N}}_2} - {{\text{N}}_2}}^{{\text{VT,}}1 \to 0} $ = 8.7 × 10−17 cm3/s, $ k_{{{\text{O}}_2} - {{\text{O}}_2}}^{{\text{VT,}}1 \to 0} $ = 1.3 × 10−14 cm3/s. 可见膨胀过程中O2(1) → O2(0)的VT过程速率系数约为N2(1) → N2(0)的100倍. 这就是N2振动冻结温度比O2高, 冻结也更早的原因.
![]()
图 7 T0 = 2000 K条件下膨胀过程中N2(v)和O2(v)质量分数导数变化Figure 7. The normalized mass fraction gradient profiles of N2(v) and O2(v) during expansion for T0 = 2000 K3.3 T0 = 5000 K条件下的流动特征分析
p0 = 1 MPa和T0 = 5000 K条件下驻室中约90%的O2离解, N2离解很少, 但置换反应生成了NO, cNO = 0.032. 气流膨胀过程中伴随着多种热化学非平衡机制: O发生复合反应, 位于高振动能级的N2, NO和O2向低能级跃迁.
图8给出了膨胀过程中组元质量分数变化, 可见在喉道下游不远处, 大约x/L = 0.65后就发生了化学冻结. cN2从驻室的0.737增至0.749, cN从4.3 × 10−3降为0, 新生成的N2约2/3来自置换反应. $c_{{\mathrm{O}}_{2}} $从0.024增长到0.077, cO从0.203降为0.159, 喉道下游约2/3的O2是新生成的, 其中超过80%来源于O的复合. NO减少了一半, 质量分数从0.032减少到0.015. 总的来说该驻室条件下膨胀过程中的化学反应主要是O的复合反应和两个置换反应, 主要发生在x/L = 0.4 ~ 0.65区域.
![]()
图 8 T0 = 5000 K条件下膨胀过程中组元质量分数变化Figure 8. Mass fraction profiles during expansion for T0 = 5000 K图9给出了膨胀过程中温度变化. 在x/L = 0.65前流动基本处于热力平衡, 之后NO和N2的振动温度与平动温度分离, 但O2的振动能一直基本与平动能保持平衡, 出口处平动温度为767 K, O2的振动温度约为881 K. N2的振动温度在x/L = 0.83处冻结于1965 K, NO的振动温度则在x/L = 0.9处才发生冻结, 冻结值为1413 K. 与T0 = 2000 K条件下相比, N2振动非平衡程度有所减轻, 其振动温度冻结位置离喉道更远, 且与平动温度差值减小: T0 = 2000 K时N2振动温度冻结值比平动温度约高1500 K, T0 = 5000 K时这一差值降至约1200 K. 原因是该总温条件下气流中存在大量O, 而VTA过程促进分子实现能级跃迁上更加高效. 例如5000 K时N2(1) + O → N2(0) + O的速率系数要比N2(1) + N2 → N2(0) + N2的高3倍左右.
![]()
图 9 T0 = 5000 K条件下膨胀过程中平动温度与振动温度变化Figure 9. Translational and vibrational temperature profiles during expansion for T0 = 5000 K图10给出了膨胀过程中O2各振动能级的质量分数变化. 驻室中尚未离解的O2中有约2/3处于激发态, 而在喷管出口处, 处于激发态的仅占约8%. 膨胀过程中基态能级的质量分数单调上升, v > 38能级的质量分数则单调下降. v = 1 ~ 6激发能级的质量分数在喉道附近有所增大, 是来自复合反应生成O2的贡献; 这之后随着向低能级跃迁, 质量分数迅速下降. 中间振动能级的质量分数在x/L = 0.85附近从降转变为略升, 对应图中的重叠区域. 图11给出了若干x/L位置O2的振动能级分布. 喉道前及附近O2振动能级都服从平衡分布, 从x/L = 0.6开始, 中高振动能级明显超越振动温度下的玻尔兹曼分布, 呈现“过分布”形态, 在喷管出口处能值在3.5 ~ 4.5 eV间的振动能级分布出现了一个明显的“突起”, 对应于图10中能级质量分数的重叠部分, 文献[29]也有类似发现.
![]()
图 10 T0 = 5000 K条件下膨胀过程中O2(v)质量分数变化Figure 10. Mass fraction profiles of O2(v) during expansion for T0 = 5000 K![]()
图 11 T0 = 5000 K条件下不同站位处O2振动能级分布Figure 11. Vibrational population distributions of O2 at different position for T0 = 5000 K下面通过观察各能级生成项分析图10中不同振动能级质量分数变化特点的原因. 图12给出了膨胀过程中O2的v = 0, 1, 3和6能级质量分数导数dcO2(v)/dx变化. dcO2(0)/dx在膨胀过程中一直为正, 对v = 1, 3, 6, dcO2(v)/dx在喉道附近为正, 之后下降, 分别在x/L = 0.665, 0.593, 0.553处取0, 后继续减小为负, 然后再增大.
![]()
图 12 T0 = 5000 K条件下膨胀过程中O2(v)质量分数导数变化Figure 12. The mass fraction gradient profiles of O2(v) during expansion for T0 = 5000 K选取dcO2(0)/dx, dcO2(1)/dx均达到峰值的x/L = 0.55和cO2接近冻结的x/L = 0.7两个典型位置, 进一步考察各微观机制的贡献, 见图13. x/L = 0.55处, dcO2(0)/dx = 0.174 m−1, dcO2(1)/dx = 0.085 m−1, v < 7的dcO2(v)/dx均为正. 各机制对dcO2(v)/dx贡献为正对应的能级为: VT过程是v = 0 ~ 3和v > 40, 且是主要贡献; O复合生成的O2主要位于中等激发能级, v = 17对应峰值dcO2(17)/ dx = 1.03 × 10−2 m−1; 置换反应在v < 9时贡献为正, 即NO + O → O2 + N起支配作用, 最大值为dcO2(1)/dx = 0.02 m−1 (其贡献占各机制总和的1/4), v > 9后逆向反应占主导. VT过程和复合反应是dcO2(v)/dx的主要贡献源, 二者共同作用下, v < 7前总的dcO2(v)/dx > 0, 之后为接近于0 (10−4量级)的负值. VV过程对各激发态能级的贡献都在10−4量级或以下, v > 7后ZE过程的贡献也降至10−4量级以下.
![]()
图 13 T0 = 5000 K条件下O2(v)质量分数导数及各微观过程的贡献Figure 13. Mass fraction gradient profiles of O2(v ) and contributions from different processes for T0 = 5000 K在x/L = 0.7位置, 低振动能级基本上只受VT过程影响, 且此处VT过程对全部激发态能级的dcO2(v)/dx贡献都为负. 图10中出现的重叠区域可通过此位置v = 19和v = 20的dcO2(v)/dx及各微观过程的贡献分析原因, 见表5. 对较高能级, dcO2(v)/dx的正负由VT过程和复合反应决定, VV过程和ZE过程的贡献可以忽略. dcO2(v)/dx > 0是因为复合反应的正贡献略大, 其实VT过程和复合反应的贡献基本上抵消了, 二者对dcO2(v)/dx的贡献均为10−4 m−1量级, 差别只是10−7 m−1量级.
表 5 x/L = 0.7处v = 19和v = 20的dcO2(v)/dx及各微观过程的贡献Table 5. dcO2(v)/dx and contributions from different processes for v = 19 and v = 20 at x/L = 0.7dcO2(v)/dx O2(19) O2(20) VT −3.356 × 10−4 −3.426 × 10−4 VV −1.273 × 10−6 −2.071 × 10−6 DR 3.368 × 10−4 3.449 × 10−4 ZE −5.800 × 10−8 −6.657 × 10−8 total −1.858 × 10−7 2.028 × 10−7 p0 = 1 MPa, T0 = 5000 K条件下, 驻室cNO = 0.032, 约x/L = 0.65处降至冻结值cNO = 0.014. 膨胀过程中NO的主要热化学过程是向低能级跃迁和置换反应. 图14给出了膨胀过程中NO各振动能级的质量分数变化, 各激发态能级均是单调降低. 图15给出了若干x/L位置的NO振动能级分布, 可见在x/L = 0.7之后开始出现过分布, 出口处呈现阶梯状.
![]()
图 14 T0 = 5000 K条件下膨胀过程中NO(v)质量分数变化Figure 14. Mass fraction profiles of NO(v) during expansion for T0 = 5000 K![]()
图 15 T0 = 5000 K条件下不同站位处NO振动能级分布Figure 15. Vibrational population distributions of NO at different position for T0 = 5000 K图16给出了膨胀过程中NO的v = 0, 1, 3, 6能级的dcNO(v)/dx变化, 可见各能级均在x/L = 0.6附近达到最小值. 基态NO质量分数在x/L = 0.655由减小变为增加, 是因为这之前, NO(0) + N → N2 + O占主导, 之后高能级向基态跃迁的贡献超过置换反应. 以x/L = 0.65处为例 (见图17), 此处两个置换反都使得NO (v < 20)减少, 并且对v = 0 ~ 2主要是反应NO + N → N2 + O, 而对v = 4 ~ 10主要是反应NO + O → O2 + N.
![]()
图 16 T0 = 5000 K条件下膨胀过程中NO(v)质量分数导数变化Figure 16. The mass fraction gradient profiles of NO(v) during expansion for T0 = 5000 K![]()
图 17 T0 = 5000 K条件下NO(v)质量分数导数及各微观过程的贡献Figure 17. Mass fraction gradient profiles of NO(v) and contributions from different processes for T0 = 5000 K3.4 T0 = 8000 K条件下的流动特征分析
p0 = 1 MPa和T0 = 8000 K条件下驻室中O2已离解完毕, N2也离解了一半, 驻室cN2 = 0.398. 图18给出了膨胀过程中各化学组元质量分数的变化. 可见气流膨胀降温过程中部分N复合生成N2, 但O基本没有复合, 尽管出口处气流温度已降至1086 K, cO = 0.24几乎保持不变. 约在x/L = 0.65处出现化学冻结, 冻结值cN2 = 0.537, cN = 0.220.
![]()
图 18 T0 = 8000 K条件下膨胀过程中组元质量分数变化Figure 18. Mass fraction profiles during expansion for T0 = 8000 K驻室的N2中处于激发态的占比约2/3, 气流膨胀过程中原有的N2和复合生成的N2均发生向低振动能级的跃迁, 出口处N2中处于激发态的比例降为约1/3. 图19给出了膨胀过程中N2各振动能级的质量分数变化, cN2(v > 6)单调减少然后趋于冻结, 不产生“重叠”区域. 图20给出了若干x/L位置的N2振动能级分布, 与T0 = 5000 K, 2000 K条件下的O2类似, 也呈现出“阶梯”形, 且高振动能级呈现出“过分布”形态.
![]()
图 19 T0 = 8000 K条件下膨胀过程中N2(v)质量分数分布Figure 19. Mass fraction profiles of N2(v) during expansion for T0 = 8000 K![]()
图 20 T0 = 8000 K条件下不同站位处N2振动能级分布Figure 20. Vibrational population distributions of N2 at different position for T0 = 8000 K图21给出了膨胀过程中N2的v = 0, 1, 3, 6能级的dcN2(v)/dx变化, 图22给出了dcN2(1)/dx达到峰值的x/L = 0.6处各振动能级的dcN2(v)/dx和各微观机制的贡献. 与T0 = 5000 K下的O2 (图12和图13)类似, x/L = 0.6处, v = 0 ~ 3的dcN2(v)/dx > 0, 且由VT过程主导, dcN2(0)/dx = 0.619 m−1, dcN2(1)/dx = 0.230 m−1, v = 2时降至0.054 m−1, 后降为负, 最小值在v = 5处取得, 为dcN2(6)/dx = −0.044 m−1, 然后绝对值继续减小, 当v > 18后dcN2(v)/dx均为很小的负值(10−4量级). 置换反应生成的N2主要位于0 ~ 10能级, 对dcN2(v)/dx的贡献在0.01 m−1量级. 复合反应生成的N2分子基本处于中高振动能级(v > 10), 峰值在v = 33处, 为2.9 × 10−3 m−1, 复合反应和少量置换反应使高能级粒子数增加的效果与VT过程作用下粒子向低能级跃迁的效果基本抵消.
![]()
图 21 T0 = 8000 K条件下膨胀过程中N2(v)质量分数导数变化Figure 21. The mass fraction gradient profiles of N2(v) during expansion for T0 = 8000 K![]()
图 22 N2(v)质量分数导数及各微观过程的贡献Figure 22. Mass fraction gradient profiles of N2(v) and contributions from different processes of N2(v)从3.2 ~ 3.4节不同总温条件下喷管流动的分析可见, 3种分子在喷管出口处都呈现出高振动能级“过分布”形态. O2和N2都是因为复合反应生成的分子更多位于中等振动能级, 与VT过程综合影响下形成“过分布”形态. Hong等[55]采用态−态模型计算了钝体的驻点线非平衡流动, 发现壁面附近也存在类似的“过分布”特征, 并且同样主要是由复合反应引起. 但两类流动又有明显差异, 喷管流动中出现热化学冻结现象, 高能级的“过分布”也与此相关. 而驻点线上气流向壁面趋近的降温过程中, 高振动能级粒子数先减小(向低能级跃迁), 在紧靠壁面的区域又开始上升, 是因为低壁温非催化壁条件下的大量原子复合到了较高的振动能级.
3.5 驻室总压的影响
本节以T0 = 5000 K条件为例, 通过对比p0 = 1, 5, 20 MPa条件下喷管流场, 分析驻室总压对流动热化学非平衡特征的影响. 图23给出了3个总压条件下膨胀过程中的温度变化. 可见驻室总压升高后, 流动的热力非平衡程度有所减弱, N2和NO振动温度的冻结位置离喉道更远, 冻结值也更低. 相应地, 气流平动温度下降平缓, 出口处平动温度升高.
![]()
图 23 不同总压条件下膨胀过程中平动温度与振动温度变化(T0 = 5000 K)Figure 23. Translational and vibrational temperature profiles during expansion f or different reservoir pressure (T0 = 5000 K)图24给出了N, O2和NO的质量分数分布. 总压升高首先抑制了驻室的离解(详见表4), 然后在气流膨胀过程中又增大了原子的复合程度, 推迟了热化学冻结. 例如p0 = 1 MPa条件下驻室的$c_{{\mathrm{O}}_2} $ = 0.024, 在x/L = 0.635处冻结, 冻结值为0.077, O大约复合了1/5; 而p0 = 20 MPa条件下驻室的$c_{{\mathrm{O}}_2} $ = 0.114, 在x/L = 0.665处冻结, 冻结值为0.213, 该总压条件下O复合了约4/5.
![]()
图 24 不同驻室总压条件下组元质量分数分布(T0 = 5000 K)Figure 24. Mass fraction profiles for different reservoir pressure (T0 = 5000 K)图25给出了3个总压条件下若干x/L位置的O2与NO振动能级分布. 提高驻室总压增大了粒子间碰撞频率, 进而使分子振动能级服从玻尔兹曼分布的区域扩大.
![]()
图 25 不同驻室总压条件下O2与NO振动能级分布(T0 = 5000 K)Figure 25. Vibrational population distributions of O2 and NO at different position for different reservoir pressure (T0 = 5000 K)4. 结 论
采用振动能级精细化态−态模型(v-StS), 在总温T0 = 2000 K ~ 8000 K、总压p0 = 1 ~ 20 MPa范围, 开展高温空气准一维非平衡喷管流动数值模拟, 分析气流热化学特性演化特点和重要微观机制, 得到主要结论如下.
(1) 对NO参与碰撞的N2和O2振动能级跃迁过程, 未见公开文献给出采用QCT等高精度方法的针对能级的速率系数. 本文提出了一种综合振动能方程松弛时间和其他类似微观过程跃迁速率系数进行折算的方法, 有效利用了现有QCT等方法的结果.
(2) T0 = 2000 K条件下几乎无化学反应. 驻室中发生了振动激发的N2和O2在气流的膨胀加速过程中向低能级跃迁, 微观机制中起主导作用的是VT过程. 喉道前流动接近热力平衡, 体现分子振动能的振动温度与平动温度一致, 喉道后二者分离, 并在喉道下游不远处发生振动温度冻结. 振动温度冻结对应着较低能级粒子数的冻结, 高能级粒子数虽然仍有变化, 但因占比很小对振动能贡献小, 不影响振动温度. 在非平衡和冻结区域, 分子能级分布均偏离了振动温度下的玻尔兹曼分布, 高能级出现“过分布”. N2的VT过程速率远低于O2, 因而振动温度冻结更早, 冻结值更高, 与平动温度差别更大.
(3) T0 = 5000 K和T0 = 8000 K条件下驻室发生了强烈的振动激发和化学反应. 膨胀过程伴随着原子的复合反应和分子从高振动能级向低能级跃迁. 对分子能级跃迁起主导作用的仍是VT过程, 主要化学反应在T0 = 5000 K条件下是O的复合反应和置换反应, T0 = 8000 K条件下则是N的复合反应. 复合反应生成的分子更多位于中等振动能级. 喷管喉道前流动接近热化学平衡, 喉道后出现非平衡, 并在喉道下游不远处发生化学冻结. T0 = 5000 K条件下O2的前几个激发态向基态跃迁速率较高, 因而对振动能贡献大的低能级粒子分布更接近平衡, 膨胀过程中振动温度基本与平动温度一致; NO和N2则发生振动温度冻结, N2的冻结更早, 冻结值更高. T0 = 8000 K条件下分子只有N2, 喉道下游不远处发生振动温度冻结. 两个总温条件下喉道下游分子能级分布都偏离了振动温度下玻尔兹曼分布, 出现高能级的“过分布”.
(4) 提高驻室总压首先能抑制驻室的离解反应, 然后可增大膨胀过程中的能级跃迁和复合反应速率, 促进振动能和平动能的交换, 推迟热化学冻结. 总压升高后振动温度与平动温度差异减小, 分子能级分布偏离玻尔兹曼分布的程度降低.
数据可用性声明
支撑本研究的科学数据已在中国科学院科学数据银行(Science Data Bank) ScienceDB平台公开发布, 访问地址为https://www.doi.org/10.57760/sciencedb.j00140.00040.
-
![]()
图 1 采用不同方法计算得到的O2 + N的多量子跃迁VT过程速率系数 (T = 5000 K)
Figure 1. Multi-quantum VT transition rates for O2 + N by different methods (T = 5000 K)
![]()
图 2 不同质量流率对应的喉道附近温度分布
Figure 2. Computed temperature profile around nozzle throat for different values of qm
![]()
![]()
图 4 T0 = 2000 K条件下膨胀过程中平动温度与振动温度变化
Figure 4. Translational and vibrational temperature profiles during expansion for T0 = 2000 K
![]()
图 5 T0 = 2000 K条件下膨胀过程中O2(v)质量分数变化
Figure 5. Mass fraction profiles of O2(v) during expansion for T0 = 2000 K
![]()
图 6 T0 = 2000 K条件下不同站位处O2振动能级分布
Figure 6. Vibrational population distributions of O2 at different position for T0 = 2000 K
![]()
图 7 T0 = 2000 K条件下膨胀过程中N2(v)和O2(v)质量分数导数变化
Figure 7. The normalized mass fraction gradient profiles of N2(v) and O2(v) during expansion for T0 = 2000 K
![]()
图 8 T0 = 5000 K条件下膨胀过程中组元质量分数变化
Figure 8. Mass fraction profiles during expansion for T0 = 5000 K
![]()
图 9 T0 = 5000 K条件下膨胀过程中平动温度与振动温度变化
Figure 9. Translational and vibrational temperature profiles during expansion for T0 = 5000 K
![]()
图 10 T0 = 5000 K条件下膨胀过程中O2(v)质量分数变化
Figure 10. Mass fraction profiles of O2(v) during expansion for T0 = 5000 K
![]()
图 11 T0 = 5000 K条件下不同站位处O2振动能级分布
Figure 11. Vibrational population distributions of O2 at different position for T0 = 5000 K
![]()
图 12 T0 = 5000 K条件下膨胀过程中O2(v)质量分数导数变化
Figure 12. The mass fraction gradient profiles of O2(v) during expansion for T0 = 5000 K
![]()
图 13 T0 = 5000 K条件下O2(v)质量分数导数及各微观过程的贡献
Figure 13. Mass fraction gradient profiles of O2(v ) and contributions from different processes for T0 = 5000 K
![]()
图 14 T0 = 5000 K条件下膨胀过程中NO(v)质量分数变化
Figure 14. Mass fraction profiles of NO(v) during expansion for T0 = 5000 K
![]()
图 15 T0 = 5000 K条件下不同站位处NO振动能级分布
Figure 15. Vibrational population distributions of NO at different position for T0 = 5000 K
![]()
图 16 T0 = 5000 K条件下膨胀过程中NO(v)质量分数导数变化
Figure 16. The mass fraction gradient profiles of NO(v) during expansion for T0 = 5000 K
![]()
图 17 T0 = 5000 K条件下NO(v)质量分数导数及各微观过程的贡献
Figure 17. Mass fraction gradient profiles of NO(v) and contributions from different processes for T0 = 5000 K
![]()
图 18 T0 = 8000 K条件下膨胀过程中组元质量分数变化
Figure 18. Mass fraction profiles during expansion for T0 = 8000 K
![]()
图 19 T0 = 8000 K条件下膨胀过程中N2(v)质量分数分布
Figure 19. Mass fraction profiles of N2(v) during expansion for T0 = 8000 K
![]()
图 20 T0 = 8000 K条件下不同站位处N2振动能级分布
Figure 20. Vibrational population distributions of N2 at different position for T0 = 8000 K
![]()
图 21 T0 = 8000 K条件下膨胀过程中N2(v)质量分数导数变化
Figure 21. The mass fraction gradient profiles of N2(v) during expansion for T0 = 8000 K
![]()
图 22 N2(v)质量分数导数及各微观过程的贡献
Figure 22. Mass fraction gradient profiles of N2(v) and contributions from different processes of N2(v)
![]()
图 23 不同总压条件下膨胀过程中平动温度与振动温度变化(T0 = 5000 K)
Figure 23. Translational and vibrational temperature profiles during expansion f or different reservoir pressure (T0 = 5000 K)
![]()
图 24 不同驻室总压条件下组元质量分数分布(T0 = 5000 K)
Figure 24. Mass fraction profiles for different reservoir pressure (T0 = 5000 K)
![]()
图 25 不同驻室总压条件下O2与NO振动能级分布(T0 = 5000 K)
Figure 25. Vibrational population distributions of O2 and NO at different position for different reservoir pressure (T0 = 5000 K)
表 1 振动能级跃迁过程
Table 1 Vibrational processes involved in the StS model
No. Processes Model Ref. 1 N2(v1) + N2(v2) ↔ N2(w1) + N2(w2) FHO present 2 O2(v1) + O2(v2) ↔ O2(w1) + O2(w2) FHO present 3 N2(v1) + O2(v2) ↔ N2(w1) + O2(w2) FHO present 4 N2(v) + N2↔ N2(w) + N2 FHO present 5 O2(v) + O2 ↔ O2(w) + O2 FHO present 6 N2(v) + O2↔ N2(w) + O2 FHO present 7 O2(v) + N2 ↔ O2(w) + N2 FHO present 8 N2(v) + N ↔ N2(w) + N QCT [41] 9 O2(v) + N ↔ O2(w) + N QCT [42-43] 10 O2(v) + O ↔ O2(w) + O QCT [40] 11 N2(v) + O ↔ N2(w) + O MQC+
RT_Ratio[38], present 12 N2(v) + NO ↔ N2(w) + NO RT_Ratio present 13 O2(v) + NO↔ O2(w) + NO RT_Ratio present 14 NO(v) + N2 ↔ NO(w) + N2 FHO present 15 NO(v) + O2 ↔ NO(w) + O2 FHO present 16 NO(v) + NO ↔ NO(w) + NO FHO present 17 NO(v) + N ↔ NO(w) + N FHO present 18 NO(v) + O ↔ NO(w) + O RT_Ratio present 
下载: 导出CSV
表 2 离解复合反应和置换反应机制
Table 2 Dissociation/recombination and Zeldovich exchange reactions involved in the StS model
No. Reaction Model Ref. 1 N2(v) + N2 ↔ 2 N + N2 FHO [34] 2 O2(v) + N2 ↔ 2 O + N2 FHO [34] 3 NO(v) + N2 ↔ N + O + N2 FHO [34] 4 N2(v) + O2 ↔ 2 N + O2 FHO [34] 5 O2(v) + O2 ↔ 2 O + O2 FHO [34] 6 NO(v) + O2 ↔ N + O + O2 FHO [34] 7 N2(v) + N ↔ 3 N QCT [41] 8 O2(v) + N ↔ 2 O + N QCT [42-43] 9 N2(v) + O ↔ 2 N + O QCT [44] 10 O2(v) + O ↔ 3 O QCT [40] 11 N2(v) + O ↔ NO(w) + N QCT [50] 12 O2(v) + N ↔ NO(w) + O QCT [51] 13 N2(v) + NO↔ 2 N + NO reaction 1 [6, 20] 14 O2(v) + NO↔ 2 O + NO reaction 2 [6, 20] 15 NO(v) + NO↔ N + O + NO 20 × (reaction 3) [6, 20] 16 NO(v) + N ↔ 2 N + O 20 × (reaction 3) [6, 20] 17 NO(v) + O↔ N + 2 O 20 × (reaction 3) [6, 20]
下载: 导出CSV
表 3 喷管流动求解与迭代确定流量耗时
Table 3 Execution time for simulation nozzle flow and iteration for mass flow rate
Model Iterativetime/s Flow calculation time/s full VV 298 32 single quantum VV 60 6
下载: 导出CSV
表 4 算例驻室条件
Table 4 The reservoir conditions tested in this work
Case T0/K p0/MPa cO,res cN,res cNO,res 1 2000 1 6.2 × 10−5 0 2.4 × 10−3 2 2000 5 2.8 × 10−5 0 2.4 × 10−3 3 2000 20 1.4 × 10−5 0 2.4 × 10−3 4 5000 1 0.203 4.3 × 10−3 0.032 5 5000 5 0.148 1.9 × 10−3 0.053 6 5000 20 0.093 9.0 × 10−4 0.069 7 8000 1 0.240 0.355 7.5 × 10−3 8 8000 5 0.232 0.177 0.021 9 8000 20 0.216 0.089 0.044
下载: 导出CSV
表 5 x/L = 0.7处v = 19和v = 20的dcO2(v)/dx及各微观过程的贡献
Table 5 dcO2(v)/dx and contributions from different processes for v = 19 and v = 20 at x/L = 0.7
dcO2(v)/dx O2(19) O2(20) VT −3.356 × 10−4 −3.426 × 10−4 VV −1.273 × 10−6 −2.071 × 10−6 DR 3.368 × 10−4 3.449 × 10−4 ZE −5.800 × 10−8 −6.657 × 10−8 total −1.858 × 10−7 2.028 × 10−7
下载: 导出CSV
-
[1] 汪运鹏, 姜宗林. 高超声速喷管设计理论与方法. 力学进展, 2021, 51(2): 257-294 (Wang Yunpeng, Jiang Zonglin. A review of theories and methods for hypersonic nozzle design. Advances in Mechanics, 2021, 51(2): 257-294 (in Chinese) doi: 10.6052/1000-0992-20-002 Wang Yunpeng, Jiang Zonglin. A review of theories and methods for hypersonic nozzle design. Advances in Mechanics, 2021, 51(2): 257-294 (in Chinese) doi: 10.6052/1000-0992-20-002
[2] Teixeira O, Páscoa J. Catalytic wall effects for hypersonic nozzle flow in thermochemical non-equilibrium. Acta Astronautica, 2023, 203: 48-59 doi: 10.1016/j.actaastro.2022.11.031
[3] Gu S, Hao J, Wang Q, et al. Influence of thermochemical nonequilibrium on expansion tube air test conditions: A numerical study. Physics of Fluids, 2023, 35(3): 036106 doi: 10.1063/5.0141281
[4] Park C, Lee SH. Validation of multitemperature nozzle flow code. Journal of Thermophysics and Heat Transfer, 1995, 9(1): 9-16 doi: 10.2514/3.622
[5] Capitelli M, Armenise I, Gorse C. State-to-state approach in the kinetics of air components under re-entry conditions. Journal of Thermophysics and Heat Transfer, 1997, 11(4): 570-578 doi: 10.2514/2.6281
[6] Park C. Review of chemical-kinetic problems of future nasa missions. I-earth entries. Journal of Thermophysics and Heat Transfer, 1993, 7(3): 385-398 doi: 10.2514/3.431
[7] Candler GV. Rate effects in hypersonic flows. Annual Review of Fluid Mechanics, 2019, 51(1): 379-402 doi: 10.1146/annurev-fluid-010518-040258
[8] Wang X, Hong Q, Hu Y, et al. On the accuracy of two-temperature models for hypersonic nonequilibrium flow. Acta Mechanica Sinica, 2022, 39(2): 122193
[9] Yang X. State-to-state dynamics of elementary bimolecular reactions. Annual Review of Physical Chemistry, 2007, 58: 433-459 doi: 10.1146/annurev.physchem.58.032806.104632
[10] Pan TJ, Stephani KA. Rovibrationally state-specific collision model for the O2(sigmag-3) + O(p3) system in dsmc. The Journal of Chemical Physics, 2021, 154(10): 104306 doi: 10.1063/5.0027411
[11] Campoli L, Kunova O, Kustova E, et al. Models validation and code profiling in state-to-state simulations of shock heated air flows. Acta Astronautica, 2020, 175: 493-509 doi: 10.1016/j.actaastro.2020.06.008
[12] Du Y, Sun S, Tan M, et al. Non-equilibrium simulation of energy relaxation for earth reentry utilizing a collisional-radiative model. Acta Astronautica, 2022, 193: 521-537 doi: 10.1016/j.actaastro.2022.01.034
[13] Adamovich IV, Macheret SO, Rich JW, et al. Vibrational relaxation and dissociation behind shock waves part 2: Master equation modeling. AIAA Journal, 1995, 33(6): 1070-1075 doi: 10.2514/3.48339
[14] Capitelli M, Armenise I, Colonna G, et al. Nonequilibrium vibrational kinetics during hypersonic flow of a solid body in nitrogen and its influence on the surface heat flux. Plasma Chemistry and Plasma Processing, 1995, 15(3): 501-528 doi: 10.1007/BF03651420
[15] Armenise I, Capitelli M, Gorse C. Nitrogen nonequilibrium vibrational distributions and non-arrhenius dissociation constants in hypersonic boundary layers. Journal of Thermophysics and Heat Transfer, 1998, 12(1): 45-51 doi: 10.2514/2.6300
[16] Armenise I, Capitelli M, Colonna G, et al. Nonequilibrium vibrational kinetics in the boundary layer of re-entering bodies. Journal of Thermophysics and Heat Transfer, 1996, 10(3): 397-405 doi: 10.2514/3.803
[17] Colonna G, Tuttafesta M, Capitelli M, et al. No formation in one-dimensional nozzle air flow with state-to-state nonequilibrium vibrational kinetics//7th AIAA/ASME Joint Thermophysics and Heat Transfer Conference. Albuquerque, USA, 1998
[18] Colonna G, Tuttafesta M, Capitelli M, et al. Non-arrhenius no formation rate in one-dimensional nozzle airflow. Journal of Thermophysics and Heat Transfer, 1999, 13(3): 372-375 doi: 10.2514/2.6448
[19] Su W, Bruno D, Babou Y. State-specific modeling of vibrational relaxation and nitric oxide formation in shock-heated air. Journal of Thermophysics and Heat Transfer, 2018, 32(2): 337-352 doi: 10.2514/1.T5271
[20] Gu S, Hao J, Wen C. On the vibrational state-specific modeling of radiating normal shocks in air. AIAA Journal, 2022, 60(6): 3760-3774 doi: 10.2514/1.J061438
[21] Gimelshein SF, Wysong IJ, Fangman AJ, et al. Kinetic and continuum modeling of high-temperature air relaxation. Journal of Thermophysics and Heat Transfer, 2022, 36(4): 870-893 doi: 10.2514/1.T6462
[22] 郑伟杰. 采用态态模型的高温非平衡流场辐射特性研究. [硕士论文]. 长沙: 国防科技大学, 2019 (Zheng Weeijie. The study of radiation characteristics of thermo-chemical nonequilibrium flow field using state-to-state model. [Master Theses]. Changsha: National University of Defense Technology, 2019 (in Chinese) Zheng Weeijie. The study of radiation characteristics of thermo-chemical nonequilibrium flow field using state-to-state model. [Master Theses]. Changsha: National University of Defense Technology, 2019 (in Chinese)
[23] 洪启臻, 王小永, 孙泉华. 高温非平衡流动中的氧分子振动态精细分析. 力学学报, 2019, 51(6): 1761-1774 (Hong Qizhen, Wang Xiaoyong, Sun Quanhua. Detailed analysis of vibrational states of oxygen in high temperature non-equilibrium flows. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(6): 1761-1774 (in Chinese) doi: 10.6052/0459-1879-19-145 Hong Qizhen, Wang Xiaoyong, Sun Quanhua. Detailed analysis of vibrational states of oxygen in high temperature non-equilibrium flows. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(6): 1761-1774 (in Chinese) doi: 10.6052/0459-1879-19-145
[24] Ninni D, Bonelli F, Colonna G, et al. On the influence of non equilibrium in the free stream conditions of high enthalpy oxygen flows around a double-cone. Acta Astronautica, 2022, 201: 247-258 doi: 10.1016/j.actaastro.2022.09.017
[25] Ninni D, Bonelli F, Colonna G, et al. Unsteady behavior and thermochemical non equilibrium effects in hypersonic double-wedge flows. Acta Astronautica, 2022, 191: 178-192 doi: 10.1016/j.actaastro.2021.10.040
[26] Guy A, Bourdon A, Perrin MY. Consistent multi-internal-temperatures models for nonequilibrium nozzle flows. Chemical Physics, 2013, 420: 15-24 doi: 10.1016/j.chemphys.2013.04.018
[27] 徐丹, 曾明, 张威等. 采用态−态模型的热化学非平衡喷管流数值研究. 计算物理, 2014, 31(5): 531-538 (Xu Dan, Zeng Ming, Zhang Wei, et al. Numerical study of thermochemical nonequilibrium nozzle flow in state-to-state model. Chinese Journal of Computational Physics, 2014, 31(5): 531-538 (in Chinese) doi: 10.3969/j.issn.1001-246X.2014.05.004 Xu Dan, Zeng Ming, Zhang Wei, et al. Numerical study of thermochemical nonequilibrium nozzle flow in state-to-state model. Chinese Journal of Computational Physics, 2014, 31(5): 531-538 (in Chinese) doi: 10.3969/j.issn.1001-246X.2014.05.004
[28] Nagnibeda E, Papina K, Kunova O. State-to-state modeling of non-equilibrium air nozzle flows//The Eighth Polyakhov’S Reading: Proceedings of the International Scientific Conference on Mechanics. Saint Petersburg, Russia, 2018
[29] Gu S, Hao J, Wen C. Air thermochemistry in the converging section of de laval nozzles on hypersonic wind tunnels. AIP Advances, 2022, 12(8): 085320 doi: 10.1063/5.0106554
[30] Gu S, Hao J, Wen C. State-specific study of air in the expansion tunnel nozzle and test section. AIAA Journal, 2022, 60(7): 4024-4038 doi: 10.2514/1.J061479
[31] Adamovich IV, Macheret SO, Rich JW, et al. Vibrational energy transfer rates using a forced harmonic oscillator model. Journal of Thermophysics and Heat Transfer, 1998, 12(1): 57-65 doi: 10.2514/2.6302
[32] Silva MLD, Guerra V, Loureiro J. State-resolved dissociation rates for extremely nonequilibrium atmospheric entries. Journal of Thermophysics and Heat Transfer, 2007, 21(1): 40-49 doi: 10.2514/1.24114
[33] Silva MLD, Loureiro J, Guerra V. A multiquantum dataset for vibrational excitation and dissociation in high-temperature O2–O2 collisions. Chemical Physics Letters, 2012, 531: 28-33 doi: 10.1016/j.cplett.2012.01.074
[34] Silva MLD, Lopez B, Guerra V, et al. A multiquantum state-to-state model for the fundamental states of air: The stellar database//5th International Workshop on Radiation of High Temperature Gases in Atmospheric Entry, Barcelona, Spain, 2012
[35] Andrienko DA, Boyd ID. Vibrational energy transfer and dissociation in O2-N2 collisions at hyperthermal temperatures. The Journal of Chemical Physics, 2018, 148(8): 084309 doi: 10.1063/1.5007069
[36] Li N, Zhang H, Cheng X. Quantum mechanical investigation of N + N2 collision: State-to-state non-reaction and exchange reaction probabilities and rate constants. Journal of Physics B : Atomic, Molecular and Optical Physics, 2021, 54(22): 225202
[37] Dani R, Makri N. Quantum state-to-state rates for multistate processes from coherences. The Journal of Physical Chemistry Letters, 2022, 13(34): 8141-8149 doi: 10.1021/acs.jpclett.2c02286
[38] Hong Q, Bartolomei M, Pirani F, et al. Vibrational deactivation in O(3P) + N2 collisions: From an old problem towards its solution. Plasma Sources Science and Technology, 2022, 31(8): 084008 doi: 10.1088/1361-6595/ac86f3
[39] 洪启臻. 高温热化学非平衡流动的精细模拟研究. [博士论文]. 北京: 中国科学院力学研究所, 2022 (Hong Qizhen. Study on detailed simulation of high temperature Thermochemical nonequilibrium flow. [PhD Thesis]. Beijing: Institute of Mechanics, Chinese Academy of Sciences, 2022 (in Chinese) Hong Qizhen. Study on detailed simulation of high temperature Thermochemical nonequilibrium flow. [PhD Thesis]. Beijing: Institute of Mechanics, Chinese Academy of Sciences, 2022 (in Chinese)
[40] Esposito F, Armenise I, Capitta G, et al. O–O2 state-to-state vibrational relaxation and dissociation rates based on quasiclassical calculations. Chemical Physics, 2008, 351(1-3): 91-98 doi: 10.1016/j.chemphys.2008.04.004
[41] Esposito F, Armenise I, Capitelli M. N–N2 state to state vibrational-relaxation and dissociation rates based on quasiclassical calculations. Chemical Physics, 2006, 331(1): 1-8 doi: 10.1016/j.chemphys.2006.09.035
[42] Esposito F, Armenise I. Reactive, inelastic, and dissociation processes in collisions of atomic nitrogen with molecular oxygen. Journal of Physical Chemistry A, 2021, 125(18): 3953-3964 doi: 10.1021/acs.jpca.0c09999
[43] Armenise I, Esposito F. N + O2(v) collisions: Reactive, inelastic and dissociation rates for state-to-state vibrational kinetic models. Chemical Physics, 2021, 551: 111325 doi: 10.1016/j.chemphys.2021.111325
[44] Esposito F, Armenise I. Reactive, inelastic, and dissociation processes in collisions of atomic oxygen with molecular nitrogen. Journal of Physical Chemistry A, 2017, 121(33): 6211-6219 doi: 10.1021/acs.jpca.7b04442
[45] Fangman AJ, Andrienko DA. Vibrational-specific model of simultaneous N2−N and N2−N2 relaxation under postshock conditions. Journal of Thermophysics and Heat Transfer, 2022, 36(3): 568-583 doi: 10.2514/1.T6441
[46] Lopez B, Silva MLD. Non-boltzmann analysis of hypersonic air re-entry flows//11th AIAA/ASME Joint Thermophysics and Heat Transfer Conference, Atlanta, USA, 2014
[47] Adamovich IV, Macheret SO, Rich JW. Vibrational relaxation and dissociation behind shock waves. Part 1-Kinetic rate models. AIAA Journal, 1995, 33(6): 1064-1069
[48] Oblapenko GP. Calculation of vibrational relaxation times using a kinetic theory approach. J Phys Chem A, 2018, 122(50): 9615-9625 doi: 10.1021/acs.jpca.8b09897
[49] Anderson JD. Hypersonic and High Temperature Gas Dynamics, 2nd ed. Reston: AIAA, 2006
[50] Bose D, Candler GV. Thermal rate constants of the N2 + O → NO + N reaction using ab initio 3A'' and 3A'potential energy surfaces. The Journal of Chemical Physics, 1996, 104(8): 2825-2833 doi: 10.1063/1.471106
[51] Bose D, Candler GV. Thermal rate constants of the O2 + N → NO + O reaction based on the A2' and A4' potential energy surfaces. The Journal of Chemical Physics, 1997, 107: 6136-6145 doi: 10.1063/1.475132
[52] Nagnibeda E, Kustova E. Non-equilibrium Reacting Gas Flows. Berlin: Springer, 2009
[53] Kim JG. Expansion of the equilibrium constants for the temperature range of 300 K to 20 000 K. International Journal of Aeronautical and Space Sciences, 2016, 17(4): 455-466 doi: 10.5139/IJASS.2016.17.4.455
[54] Brown PN, Byrne GD, Hindmarsh AC. Vode: A variable-coefficient ode solver. SIAM Journal on Scientific and Statistical Computing, 1989, 10(5): 1038-1051 doi: 10.1137/0910062
[55] Hong QZ, Wang XY, Hu Y, et al. Development of a stagnation streamline model for thermochemical nonequilibrium flow. Physics of Fluids, 2020, 32(4): 046102
-
期刊类型引用(1)
1. 郭凯，郭志文，董温哲，程雨轩，于美琪，张红升. 并排小间距比塔式容器的流致振动试验研究. 压力容器. 2024(11): 1-11+32 . 
百度学术
其他类型引用(1)
计量
- 文章访问数: 436
- HTML全文浏览量: 101
- PDF下载量: 77
- 被引次数: 2
