Q-K模型在氮氧离解复合反应中的评估

陈 浩, 李林颖, 张 斌,2,*, 刘 洪(. 上海交通大学 航空航天学院, 上海 200240; 2. 上海交通大学 数学系, 上海 200240)

0 引 言

高超声速再入飞行通常发生在稀薄气体环境中，由于激波层和高空低密度的特征，在激波下游会产生显著的热化学非平衡松弛现象[1]。流场中的化学反应对飞行器气动性能和表面传热有着显著影响，多原子振动激发和离解反应的发生对高超声速飞行器来说是很重要的特征，低温下的振动激发和有限的氧气离解仍然很重要[2]。

高速稀薄气体环境下，直接模拟蒙特卡洛(DSMC)方法[3]在模拟热化学非平衡方面具有很好的优势，该方法是目前从自由分子流区到稀薄过渡流区用于数值分析气体流动问题的主要模拟工具。DSMC方法的化学反应模型发展有着近40年的历史，总能碰撞(TCE, total collision energy)模型是最为广泛使用的一种[4]。该模型基于碰撞理论，将标准Arrhenius形式的连续速率方程转变为反应碰撞截面，因此需要每一个化学反应的Arrhenius速率系数[3-4]。虽然DSMC化学反应模型发展比较缓慢并且种类不多，但在精细化学反应模拟方面的研究很多。这一类模型直接采用QCT(Quasi classical trajectory)计算结果[5]，或者采用态-态碰撞截面实验数据[6]来描述碰撞截面。Bird[7]于2011年提出了一种全新的化学反应模型：量子-动理学(Q-K, quantum-kinetic)模型。该方法基于量子振动模型，几乎所有的化学反应概率都基于碰撞分子的振动能态而描述，不需要采用任何的宏观速率方程，也不存在气体平衡态的假设[2,4,7]。

TCE模型不能细致的描述内能(振动能、转动能)对反应过程的贡献，针对不同反应采用不同的碰撞总能来满足对连续反应速率的准确模拟[2,4,8-9]。而Q-K模型将分子振动能态与化学反应直接相关联，根据平衡原理得到精准的模拟结果。本文关注的便是两种模型在这一类细致模拟的问题上会有怎样的不同，因此本文旨在评估Q-K模型模拟平衡态、非平衡离解反应速率的能力，完成Q-K模型和TCE模型在离解度、非平衡反应速率和振动非平衡效应方面的对比研究工作。

下一节简略地介绍了Q-K模型的特点，以及离解、复合反应的在DSMC程序中的实施策略。结果与讨论部分分别在平衡态离解度和非平衡离解问题上对Q-K模型进行了数值评估，同时对比分析了Q-K和TCE模型在模拟非平衡离解速率和离解速率中的温度非平衡效应两方面的差异，并与经典Park双温度模型进行了对比分析。

1 Q-K化学反应模型

振动能态在Q-K模型中具有重要的作用，离解和置换反应的发生都是基于碰撞分子对的振动能态而描述的，采用量子化碰撞温度能够实现振动能平衡均分的结果[10]。Q-K模型采用量子振动和简谐振子模型，处理振动能松弛时采用“量子碰撞温度”求解松弛碰撞数。量子碰撞温度可写为：

(1)

式中，imax为碰撞分子对所具有的最大能级，ω为气体黏性系数中的温度指数。采用量子碰撞温度Tcoll定义振动松弛碰撞数如下[11]：

(2)

其中，Zref为参考温度Tref下的振动松弛碰撞数，该松弛碰撞数可通过TCE模型中求解振动碰撞数的方式获得[2]：

(3)

式(2)和(3)中的参考温度Tref通常用振动特征温度θv表示[2]。

本文所考虑的化学反应类型包括离解反应和复合反应两种，Q-K模型中两种化学反应的处理方式简要表达如下小节所示，详细的Q-K理论可查阅文献[7]。

1.1 离解反应

离解反应的发生条件可描述为：一次碰撞的能量足够大以至于能够达到离解阀值所对应的振动能级。考虑一个离解反应AB+T=A+B+T，AB为离解对象，T为碰撞对象，A和B为分子AB的离解产物。碰撞后的最大振动能级imax由方程(4)确定：

(4)

imax>θd/θv

(5)

其中，θd为分子AB的特征离解温度。

1.2 复合反应

A和B发生复合反应的条件为：有另外一个分子或者原子进入到碰撞体积Vcoll中，这里的体积为“三元碰撞体积”。因此，复合反应A+B+T=AB+T发生的概率可描述为：

Prec=nVcoll=naTbVref

(6)

式中，n为分子数密度，Vref为参考碰撞体积，其值为半径等于三个碰撞分子(A、B和T)半径之和的球体的体积；系数aTb中的温度一般采用宏观温度，则T=Tcell/θv，Tcell为单位网格内的宏观温度；系数aTb中的参数a、b根据理论平衡常数进行调整，推荐值如表1所示。

2 结果与讨论

2.1 基于离解度的Q-K模型评估和对比分析

文献[7]对Q-K理论直接推导的反应速率进行了详细的验证，Body教授[2]通过零维热浴和圆柱扰流数值模拟将Q-K模型的DSMC采样结果与实验数据[12-15]、传统TCE模型[14]进行了详尽的对比，证明了Q-K模型在DSMC程序中的可行性，因此本文不再赘述平衡态反应速率的对比结果。工程上关注更多的是化学反应模型在模拟气体组分和反应后流场温度上的精准度，因此我们对标准分子数密度状态下(n=2.687×1025m-3)的空气进行热浴实验模拟，对比Q-K和TCE化学反应模型在模拟双原子离解度上的精准性。

表1 三元碰撞体积参数Table 1 Parameters in ternary collision volume

空气组分为78%的氮气和22%的氧气，采用1×104个模拟分子进行热浴模拟。注意的是，该模拟中只考虑氮气和氧气的离解、复合反应，因此不会有NO分子的生成。两种化学反应模型的对比结果如图1所示，图中每一个点代表当前平衡温度下所对应的离解度。圆形和三角形标记分别为Q-K和TCE模型的计算结果，线条为离解度理论分析值[16]，红色表示氮气，蓝色表示氧气。我们发现，Q-K和TCE模型所得到的结果几乎完全一致。相比于氮气，氧气表现出极易容易离解的特征。当温度约为5000 K左右时，氮气开始离解，而氧气离解度则高达90%。此外，如图1椭圆形区域所示，Q-K模型表现出比TCE模型更精确的仿真能力，TCE模型得到的离解度过低，当平衡温度约为6000 K时，离解度相对误差高达75%。

图1 不同温度下的氮氧离解度Fig.1 The degree of dissociation for nitrogen and oxygen

根据两种化学反应模型的理论和实现过程，可理性地得到TCE模型在低离解度时产生大误差的原因。TCE离解反应的发生是建立在振动松弛条件之下的，即离解反应的触发是由振动能唯一控制的。然而空气中的化学反应并非由振动能唯一控制[16]，在分子振动能还没有高度激发的状态下，平动能和转动能对化学反应的影响将变得不可忽略。Q-K模型能得到精准离解度的原因在于，Q-K离解反应的触发是独立于振动松弛而设置的，并且由方程(4)可看出，离解反应发生条件中的最大振动能级imax考虑到了相对平动能对反应的贡献。

2.2 基于非平衡离解的Q-K模型评估和对比分析

通常来说，非平衡效应有两种。一种是速度分布不再满足Maxwell分布，另一种是内能松弛引起的温度非平衡。本节讨论的非平衡便是由平动温度、转动温度和振动温度所体现的温度非平衡效应。

考虑三种非平衡的情况：1、平动-转动温度保持在6000 K，振动温度逐渐上升；2、平动-转动温度保持在10000 K，振动温度逐渐上升；3、平动-转动温度保持在20000 K，振动温度逐渐升高。在DSMC程序中很容易实现这种温度非平衡：振动能在碰撞过程中不再与平动能、转动能进行能量交换，而是按照初始给定的振动温度给予振动能Boltzmann分布。作为比照的参考，本文采用Park[12]发展的双温度模型从理论上描述前述三种非平衡条件下的非平衡反应速率。在Park的双温度模型中，采用一个平动-转动温度和振动温度耦合的几何平均温度来替代平衡态Arrhenius化学反应速率方程中的宏观温度，以此来表征振动温度对化学反应速率的影响，该几何平均温度Ta定义如下：

(7)

式中，Tv为振动温度，T为平动-转动温度。针对氮气离解反应，Sharma[17]建议取值q=0.3。

针对上述的人为的给定的三种非平衡离解反应过程，对Q-K模型、TCE模型和Park双温度模型得到的非平衡离解反应速率进行了对比。比较结果如图2所示，线条表示Park非平衡离解反应速率、标记块为Q-K模型和TCE模型模拟所得到的非平衡速率，具体标记区别如图2中图例所示。

如图2所示，TCE模型所描述的非平衡离解速率更接近Park理论，而Q-K模型的结果比Park非平衡速率要大近1～2个量级左右。还可以发现，随着平动-转动温度的增加，非平衡反应速率对振动温度的敏感性变得很弱，其中TCE模型的这种趋势更为明显。当平动-转动温度降低时，振动温度的升高能有效的提高化学反应速率，其中Q-K模型的这种趋势更为明显；当振动温度很低时，两种化学反应模型DSMC模拟结果与Park双温度模型之间的误差很大，其原因在于：Park双温度模型中的非平衡反应速率对几何平均温度Ta的依赖导致了反应速率的快速下降，以及DSMC程序中采用了相对平动能来满足反应所需的激活能。

图2 Q-K模型和TCE模型的非平衡离解速率Fig.2 The nonequilibrium dissociation rate using Q-K and TCE model

另外，本文还对Q-K模型和TCE模型描述温度非平衡效应的能力进行了对比研究。专著[18]中建议采用比值Kf(T,Tvib)/Kf(T,T)来考察振动非平衡效应对离解反应速率的影响，本文也采用该比值进行对比分析。随着振动温度与平动-转动温度之间的温差越大，这种温度非平衡效应也就越强。为得到较大程度的振动非平衡效应，本文取定振动温度为Tvib=6000 K，而平动-转动温度Ttran-rot从6000 K逐渐上升至20000 K。同样，采用Park双温度模型作为理论参照，并使用两种q值，一种为Park推荐值[12]：q=0.5，另一种为Sharma推荐值[17]：q=0.3。

图3 温度非平衡对离解速率的影响Fig.3 The effect of vibrational nonequilibrium on the dissociation rate

Q-K模型、TCE模型以及Park双温度模型的比较结果如图3所示，纵坐标为比值Kf(T,Tvib)/Kf(T,T)的对数结果，横坐标为平动-转动温度。首先，由图左上角所示，四条曲线都相交于一点(6,0)，该点意味着：当平动-转动温度与振动温度均等于6000 K时，比值Kf(T,Tvib)/Kf(T,T)为1，即此时并没有非平衡效应。随着平动-转动温度的升高，曲线逐渐下降并逐渐趋于稳定，这表明振动非平衡效应先增强，然后会趋于一定的程度，该趋势与图2中所描述的“当平动-转动温度很高时，振动温度对非平衡反应速率的影响变得很小”这一特征相吻合。更重要的是，图3中四条曲线的对比表明：在描述离解速率中的温度非平衡效应时，Q-K模型与Park双温度模型(q=0.5)一致吻合，而TCE模型却过大地预估了这种非平衡效应。

3 结 论

基于工程上比较关心的大气高超声速流动中的氮氧离解反应，本文对近期提出的量子-动理学(Q-K, quantum-kinetic)模型[7]在模拟平衡态离解度、非平衡离解反应速率以及温度非平衡方面进行了评估，并基于广泛使用的TCE模型进行详细的对比研究，得到如下结论：

1) 在模拟平衡态离解度方面，Q-K模型与理论值精准吻合，而TCE模型在低离解度时具有较大误差。该结果表明，Q-K模型中将离解判断和振动松弛独立而设，并采用能级判断的方式是满足物理真实的；

2) 在模拟非平衡离解速率方面，两模型均与理论描述大致吻合，相比而言TCE结果与Park双温度模型更为接近，而Q-K非平衡离解速率则比Park理论大1～2个量级；

3) 在模拟离解速率中的温度非平衡效应方面，Q-K模型与Park双温度模型高度吻合，而TCE模型却过大地预估了这种非平衡效应。

[1]Park C. Review of chemical-kinetic problems future NASA missions, I: Earth entries[J]. Journal of thermophysics and heat transfer, 1993, 7(3): 385-398

[2]Scanlon T J, Boyd I D, et al. Open-source direct simulation Monte Carlo chemistry modeling for hypersonic flows[J]. AIAA Journal, 2015, 53(6): 1670-1680

[3]Bird G A. Molecular gas dynamics and the direct simulation of gas flows[M]. Oxford: Oxford Univ. Press, 1994

[4]Wysong I, Gimelshein S, et al. Comparison of direct simulation Monte Carlo chemistry and vibrational models applied to oxygen shock measurements[J]. Physics of Fluids, 2014, 26(043101): 1-15

[5]Levin D A, Gimelshein S F. Modeling of OH vibrational distributions using molecular dynamics with the direct simulation Monte Carlo method[C]//34thThermophysics Conference, Doi: 10.2514/6.2000-2432

[6]Kim J G, Boyd I D. Monte Carlo simulation of nitrogen dissociation based on state-resolved cross sections[J]. Physics of Fluids, 2014, 26(012006): 1-20

[7]Bird G A. The Q-K model for gas-phase chemical reaction rates[J]. Physics of Fluids, 2011, 23(106101): 1-13

[8]梁杰, 李志辉, 杜波强, 等. 探月返回器稀薄气体热化学非平衡特征性数值模拟[J]. 载人航天, 2015, 21(3): 295-302

[9]黄飞, 赵波, 程晓丽, 等. 返回器高空稀薄气动特性的真实气体效应研究[J]. 宇航学报, 2014, 35(3): 283-290

[10]Bergemann F, Boyd I D. New discrete vibrational energy model for the direct simulation Monte Carlo method[J]. Progress in Astronautics and Aeronautics, 1994, 158: 174-183

[11]Bird G A. A comparison of collision energy-based and temperature-based procedures in DSMC[C]//Rarefied Gas Dynamics: 26thInternational Symposium, 2008, 1084: 245-250

[12]Park C. Nonequilibrium hypersonic aerothermodynamics[M]. New York: A wiley-interscience publication, 1989

[13]Gupta R N, Thompson R A, et al. A review of reaction rates

and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30000 K[R]. Hampton, NASA Langley Research Center, 1990

[14]Carlson A B, Bird G A. Implementation of a vibrationally linked chemical reaction model for DSMC[R]. Hampton, NASA Langley Research Center, 1994

[15]Gallis M A, Bond R B. Assessment of collision-energy-based models for atmospheric species reactions in hypersonic flows[J]. Journal of the Thermophysics and Heat Transfer, 2010, 24(2): 241-253

[16]杨浩森, 张斌, 刘洪, 等. DSMC热化学模拟分子数对氮氧离解反应计算影响的分析[J]. 空气动力学学报, 2014, 32(4): 511-517

[17]Sharma S P, Huo W M, Park C. Rate parameters for coupled vibrational-dissociation in a generalized SSH approximation[J]. Journal of Thermophysics and Heat Transfer, 1988, 6: 9-21

[18]吴其芬, 陈伟芳. 高温稀薄气体热化学非平衡流动的DSMC方法[M]. 长沙: 国防科技大学出版社, 1999.