态-态模型下的O2/O系统热化学非平衡与辐射过程

郑伟杰, 曾 明, 王东方, 谢 易

(国防科技大学 空天科学学院, 长沙 410073)

0 引 言

在高超声速飞行条件下，流场中的高温导致气体各个内能模式激发并且发生化学反应，而高空条件气体密度低又使流场处于热力和化学非平衡状态。当温度足够高时，气体进一步发生辐射跃迁。各类微观反应和跃迁机制影响到飞行器的宏观气动热和气动物理特性。更深入地理解和更精确地模拟热化学非平衡和辐射现象在学术上和工程应用上都具有重要意义。

有两类描述流动中振动非平衡和振动-化学耦合的方法。一类是多温度模型法，以Park在1985年提出的双温度模型[1]为代表，引入不同于平动温度的振动温度Tvib表征非平衡振动能，在化学反应控制温度中引入振动温度反映热力-化学耦合。另一类是态-态模型法[2-3]，直接描述各振动能级间的分子跃迁过程和各能级分子的化学反应，求解各能级粒子数变化率方程。

双温度模型及后来拓展的多温度模型，计算消耗小，且一定程度上反映了振动非平衡与振动-化学耦合，获得了广泛应用。许多气动热力学程序采用了该模型，包括Gnoff等发展的LAURA (Langely Aerothermodynamic Upwind Relaxation Algorithm)程序[4-5]和GASP(General Aerodynamic Simulation Program)程序[6-7]。采用该模型的数值模拟工作还有很多，例如：Tissera[8]等进行的高超声速绕钝锥-柱-裙的流动模拟，Sorensen等[9]进行的有限催化壁条件下的高超声速圆柱绕流模拟，美国和欧洲研究者们独立进行的为评估CFD方法预测高超声速条件下激波相互作用能力的一系列数值计算[10]，汪球[11]、周凯[12]等进行的高焓风洞喷管非平衡流场计算，等。双温度模型在许多条件下获得了良好的数值结果，但也存在局限[13]。在温度很高条件下以及降温复合反应过程中，可能会产生较大误差，甚至出现模型失真[14]。双温度模型隐含着振动能级分布满足振动温度下玻尔兹曼(Boltzmann)分布的假设[15]，而态-态模型方法的结果表明该假设在许多条件下并不成立。

态-态模型法将粒子在能级间的跃迁与化学反应视为同一性质的过程研究，描述振动非平衡及其与化学反应耦合的方法则更为精细和直接[16]。它直接研究各能级上粒子的化学反应，避免了多温度模型中经验性的热力-化学耦合方法。随着分子动力学的研究和计算机计算能力的发展，态-态模型作为一种高保真度模型逐渐进入人们的视野并在近十年来快速发展起来。采用态-态模型的研究由不涉及流动的零维问题[17]，逐步发展到正激波后、边界层[18-21]和准一维喷管[22]非平衡流动问题。徐丹对N2/N系统采用态-态模型开展了零维问题研究[23]，还数值求解了耦合态-态模型的N2/N准一维喷管非平衡流方程[16]。Druguet[24]采用态-态模型(考虑N2基态电子能级上68个振动能级)求解了N2/N的球头非平衡流场，考察态-态模型与流动耦合的可行性。应用CFD程序PINENS(并行隐式非平衡N-S方程求解程序)求解了包含69个组元连续方程的轴对称非平衡流N-S方程，结果表明仅在边界层外振动能级分布满足振动温度下的玻尔兹曼分布。

本文对封闭静止的O2/O气体系统，采用态-态模型研究热化学非平衡和辐射过程。考察不同初始条件和设定温度条件下(温度由初始值瞬时改变为设定值后保持系统等温等容)不同振动能级粒子分布随时间的演化特点，分析比较各类跃迁或化学反应过程以及辐射跃迁过程对系统化学组成和分子振动能级分布变化的贡献。为更好地理解O2/O和高温空气热化学反应机理提供新的思路，也为今后采用态-态模型研究高温空气的非平衡过程并进一步与流动耦合奠定基础。

1 态-态模型基本方程

1.1 采用态-态模型的粒子数密度变化率方程

(1)

式中i=0～46。式(1)右端分别为振动-平动跃迁、振动-振动跃迁、离解-复合化学反应和辐射几类微观过程对基态电子能级的i振动能级上O2分子数密度变化的贡献。热化学反应过程将在1.2节详细描述，辐射跃迁过程在1.3节详细描述。

O原子数密度nO的变化仅起因于离解-复合化学反应：

(2)

(3)

则,

(4)

(5)

gX和gB分别为基态和激发态电子能级的简并度，εel,X和εel,B分别为基态和激发态电子能级的能值，kB为玻尔兹曼常数。激发态电子能级上各振动能级的分布则根据振动温度Tvib下的玻尔兹曼分布确定。系统中总的O2分子数密度为:

(6)

1.2 不涉及辐射的能级跃迁与化学反应过程

不涉及辐射的微观过程包括：分子间的振动-平动能量交换过程VTm，分子与原子间的振动-平动能量交换过程VTa，振动-振动能量交换过程VV，原子作为碰撞参与者的离解-复合反应DRa，分子作为碰撞参与者的过程DRm。具体如下：

(7)

(8)

(9)

(10)

(11)

对VTm过程只考虑单量子数跃迁，对VTa过程考虑多量子数跃迁，Δv为跃迁量子数。VV过程中，i和j代表O2基态电子能级上的不同振动能级。

(12)

式(12)右端当i=0时仅含第1和第3项，i=46时仅含第2和第4项。

(13)

(14)

(15)

(16)

O原子数密度变化方程式(2)右端的两项生成项为:

(17)

(18)

1.3 辐射跃迁过程

辐射跃迁包括三类——自发发射、诱导发射、辐射吸收。具体表示如下：

(19)

(20)

(21)

(22)

式(22)中c为光速。Aα′i′,αi、Bα′i′,αi、Bαi,α′i′分别为自发发射、诱导发射和辐射吸收的爱因斯坦系数。根据细致平衡原理，它们之间满足如下关系：

(23)

gα′i′Bα′i′,αi=gαiBαi,α′i′

(24)

(25)

(26)

本文对三种辐射跃迁机制的谱线均采用多普勒增宽的线型分布，谱线互不重叠。将式(22)对频率积分后得：

(27)

(28)

1.4 振动温度的确定

1.4.1 根据非平衡振动能级分布确定振动温度

一种方法是通过假设非平衡的各振动能级分布尽可能接近振动温度下的玻尔兹曼分布导出一个统一的振动温度。需要采用线性回归的方法。

振动温度Tvib下的玻尔兹曼分布为:

(29)

式(29)中εvib,i为位于振动能级i的一个分子具有的振动能。Qvib为振动配分函数,

(30)

对式(29)两边取对数，得:

(31)

近似将上式中的振动配分函数Qvib视为常数，则i能级粒子数密度的对数便为该能级能值的线性函数，记作:

(32)

其中K1为该线性函数的斜率,

(33)

根据已知求得的各振动能级非平衡粒子数密度，采用线性回归方法确定式(32)的斜率，即得到对应的振动温度Tvib。

(34)

式(34)中wi为权重因子，根据该能级的粒子占比确定：

(35)

(36)

(37)

1.4.2 根据平均的非平衡振动能确定振动温度

(38)

2 计算方法

由于各类跃迁机制的速率差别大，导致常微分方程组的数值求解存在严重刚性问题，可以采用Euler隐式方法求解。另外，VODE(变系数的常微分方程求解器)[37]是变系数的BDF(向后分化公式)方法，可以高效地解决具有刚性的问题。α-QSS(拟稳态逼近方法)[38]是具有二阶精度的预测-校验方法，常用于反应速率系数随时间缓变的多组元化学反应模拟。考虑到Euler隐式算法形式简单、应用方便，对本文算例从计算精度和计算消耗上可以接受，本文采用了Euler隐式算法。

(39)

其中，

(40)

(41)

由于系统在初始瞬间改变温度后保持等温等容，各生成项中的速率常数在起始时刻根据系统设定温度计算一次即可。时间推进的步长视具体温度和压力条件，可在10-11～10-8s之间取值。本文涉及的不耦合辐射的算例条件(设定温度分别为3000 K、10 000 K、20 000 K)下，时间步长均取为10-10s。包含辐射机制时，由于辐射强度随时间变化率较大，时间推进步长需进一步减小至10-14～10-13s，本文含辐射的算例条件(设定温度20 000 K)时间步长均取为10-13s。

3 算例及计算结果分析

3.1 态-态模型热化学非平衡程序校验

首先根据文献[35]计算结果对本文编制的O2/O态-态模型热化学非平衡程序进行校验。文献[35]给出了两组算例，初始的O2/O数密度均为nO2,0=1×1017/cm3、nO,0=9×1017/cm3，设定的系统温度分别为T=3000 K和T=10 000 K。两组算例分别设定不同的初始温度：T=3000 K组的初始温度分别为T0=100 K和T0=10 000 K;T=10 000 K组的初始温度分别为T0=100 K和T0=20 000 K。模拟了两组系统在起始时刻突然升温或降温并保持等温等容条件下的热化学非平衡过渡过程。

图1、图2分别给出了两组算例系统的振动温度随时间的变化，振动温度为应用式(38)根据基态电子能级上振动能级分布导出的Tvib。图中“dv”代表VTa过程中的振动跃迁量子数Δv，all代表允许Δv取任意可能值。可见本文计算结果与文献[35]的计算结果一致性良好，验证了本文方法和程序的正确性。

图1 由T0=100 K或T0=10 000 K升温或降温至T=3000 K后振动温度随时间变化Fig.1 Time evolution of vibrational temperature (T0=100 K or T0=10 000 K,T=3000 K)

图2 由T0=100 K或T0=20 000 K升温或降温至T=10 000 K后振动温度随时间变化Fig.2 Time evolution of vibrational temperature(T0=100 K or T0=20 000 K,T=10 000 K)

从图1和图2均可发现限制跃迁量子数Δv的影响，可见VTa过程中的多量子数跃迁的重要性。不过Δv<6与不限制Δv的结果差异就不大了。无论是升温还是降温至T=3000 K条件，随着时间推进，振动温度都是近似以指数方式趋近于平动温度，约在10-6s时与平动温度相等，基于振动温度的定义，可以说振动能和平动/转动能此时基本达到了平衡。但设定温度T=10 000 K条件的情况中振动温度随时间的变化特点不同，见图2。初始T0=100 K条件下，振动温度在10-7s升至6000 K附近后有一段平台区，到10-6s后又开始迅速爬升，达到10 000 K，整体呈现了台阶形状。初始T0=20 000 K条件下，振动温度在10-7s前一直下降，达到6000 K左右，之后和T0=100 K条件的变化曲线几乎重合，经历平台区后再迅速爬升至T=10 000 K。振动温度在10-7～10-6s期间的平台区与此时开始的化学反应有关，后面再进一步分析。

3.2 O2/O系统的热化学非平衡过程分析

图3分别给出了由T0=100 K或T0=10 000 K升温或降温至T=3000 K的两个算例O2分子归一化振动能级分布随时间的变化曲线，图中不同能级采用了不同颜色的曲线，符号标识对VTa过程中振动跃迁量子数Δv的限制。这里和后面的振动能级分布均指基态电子能级上的振动能级分布。初始时刻振动能级分布为初始温度下的平衡分布，因此由T0=100 K升温的算例(见图3a)中，随着时间推进，低能级分子数减小，高能级分子数增大；低能级(i<10)分布在10-7s时就接近平衡，但i>20能级上的O2分子，在10-6～10-4s区间存在一个粒子数密度基本不变的平台区，到10-2s后才接近平衡。由T0=10 000 K降温的算例，随时间推进，高、低能级分子数变化的趋势与由T0=100 K升温的算例相反，但分布达到平衡的时间和过渡过程出现的平台区与初始温度为T0=100 K的情况类似。这是因为振动跃迁和化学反应速率均取决于设定的系统温度(T=3000 K)。两个算例中都体现了限制跃迁量子数的影响，特别是对高振动能级，多量子数的跃迁更加重要。与图1给出的振动温度随时间变化曲线比较，可知高振动能级分布达到平衡的松弛时间比振动温度的长得多，这正是高振动能级粒子分布不满足振动温度下的玻尔兹曼分布的体现。不过，由于高能级分子占比很小，其对总体振动能的影响也很小。

(a) T0=100 K，升温

振动能级分布在非平衡过渡过程中的变化与化学反应相关。这里首先以初始T0=10 000 K降温至T=3000 K的算例为代表进行分析。图4给出了O2与O粒子数密度随时间的变化，图5对比了包含所有振动跃迁、化学反应过程(VTm、VTa、VV、DRa、DRm)和去除化学反应计算得到的振动能级分布随时间的变化，图中含所有振动跃迁和化学反应过程的用实心符号，去除化学反应的则用空心符号。图6对比了包含或去除化学反应计算得到的振动温度变化曲线。图7则分别给出了包含或去除化学反应时振动能级分布与对应振动温度下的玻尔兹曼分布的对比。

图4 由T0=10 000 K降温至T=3000 K后O2和O数密度变化Fig.4 Time evolution of number densities of O2 and O(T0=10 000 K, T=3000 K)

(a) 低振动能级

图6 由T0=10 000 K降温至T=3000 K后振动温度随时间变化(含或不含化学反应的对比)Fig.6 Time evolution of vibrational temperature (with or without DR processes, T0=10 000 K, T=3000 K)

(a) 去除化学反应

相对于算例设定的温度T=3000 K下的平衡化学组成，初始时刻的nO,0/nO2,0偏大，非平衡过渡过程中将发生O原子的复合反应。化学反应有一定的“孵化”时间，由图4可见该设定温度T=3000 K条件下大约10-6s复合反应才起步，到10-5s后才有明显的反应，大约10-2s反应达到平衡。结合图5给出的有无化学反应的振动能级分布对比，可以看到二者的高振动能级分布的明显差异正是化学反应开始的时间。可以分析得出这里的复合反应生成的O2分子主要是在高振动能级。高振动能级的分子分布受到化学反应的明显影响，而低能级分布几乎不受影响。10-6s～10-4s之间出现的高振动能级分子数密度的平台区，是该区间复合生成高能级O2和高能级O2通过振动松弛向低能级跃迁互相平衡的结果。

由图6可知算例条件下，是否包含化学反应对振动温度没有明显影响，这和高能级分子占比小因而对振动能贡献小的特点是一致的。另外可看出以振动温度衡量的振动松弛时间大约为10-6s，远小于以高振动能级分布衡量的松弛时间和离解-复合反应达到平衡所需要的时间(约10-2s)。

从图7可见，包含或去除化学反应时振动能级分布均偏离了对应振动温度下的玻尔兹曼分布，包含化学反应时偏离程度更大，特别是高振动能级的分布。从随时间的变化来看，由图7(b)可见t<10-7s时，化学反应尚未发生，振动松弛占主导地位，位于高能级的O2分子数低于振动温度下的玻尔兹曼分布。而t>10-6s后，复合反应开始，生成大量的高能级O2分子，使得位于高能级的O2分子数密度高于此时振动温度下的玻尔兹曼分布。

接下来对初始T0=100 K升温至T=10 000 K的算例，分析化学组成随时间的变化(见图8)，振动能级分布和振动温度的演化特点及其受化学反应的影响(见图9、图10)，不同时刻振动能级分布与振动温度下玻尔兹曼分布的差异(见图11)。相对于算例设定的温度T=10 000 K下的平衡化学组成，初始时刻的nO,0/nO2,0偏小，非平衡过渡过程中将发生O2分子的离解反应。在设定温度下，振动跃迁和化学反应的速率都较T=3000 K算例条件下更快，松弛时间也大大减小。

图8 由T0=100 K升温至T=10 000 K后O2和O数密度变化Fig.8 Time evolution of number densities of O2 and O(T0=100 K, T=10 000 K)

由图8～图10可知，化学反应达到平衡的特征时间、以振动温度或高振动能级分布表征的松弛时间，都在10-7s～10-6s之间。化学反应对振动能级分布的影响，比T=3000 K算例条件下更大。离解反应主要发生在10-8s～10-6s时间段内，且存在高振动能级优先离解的特点。由图11(a)可见，不包含化学反应时，高振动能级的O2分子数密度高于振动温度下的玻尔兹曼分布。而包含化学反应后，较高能级的分子优先离解，导致高能级分子数低于振动温度下的玻尔兹曼分布值，见图11(b)。

(a) 低振动能级

图10 由T0=100 K升温至T=10 000 K后振动温度随时间变化(含或不含化学反应的对比)Fig.10 Time evolution of vibrational temperature (with or without DR processes, T0=100 K, T=10 000 K)

(a) 去除化学反应

对比分析T=3000 K和T=10 000 K算例可知，较低的设定温度下，化学反应所需时间比振动松弛时间大得多，两个过程相对独立，化学反应对振动松弛产生的影响较小。而高设定温度下，化学反应特征时间和振动松弛时间可以相比拟，彼此耦合程度更高，互相影响更大。

3.3 O2/O系统的热化学非平衡与辐射过程分析

选取设定温度为T=20 000 K，初始的O2和O数密度均为nO2,0=1×1017/cm3、nO,0=9×1017/cm3，初始温度分别为T0=10 000 K和T0=30 000 K的两个算例条件，对O2/O系统研究高温下耦合辐射的热化学非平衡过渡过程，考察辐射生成项对粒子能级分布的贡献。

首先研究总的辐射强度。将式(27)给出的针对跃迁α′i′→αi的辐射带强度对所有的跃迁求和，即得到总的辐射强度随时间的变化率方程。本文考虑的辐射跃迁包括由激发态电子能级上振动能级i′=0～14至基态电子能级上振动能级i=0～21的共计330个辐射谱带跃迁。总的辐射强度记为I。

令式(27)为0，得到平衡辐射谱带强度：

(42)

图12给出了由T0=30 000 K或T0=10 000 K降温或升温至T=20 000 K后三个总辐射强度I、Ieq、Ieq,B辐射强度随时间的变化。为分析辐射与热化学非平衡过程的相互影响特点，图13、图14还给出了相应的O2、O数密度和振动温度的时间演化曲线。由图12可见随着时间的推进，辐射强度I和平衡辐射强度Ieq最终将达到设定温度下的黑体平衡辐射强度Ieq,B。可将辐射强度I达到黑体平衡辐射强度Ieq,B的时间定义为非平衡辐射特征时间，则由T0=10 000 K升温至设定T=20 000 K的算例中，非平衡辐射特征时间在10-5s量级，而振动松弛时间(参见图14)在10-7～10-6s量级。由T0=30 000 K降温至T=20 000 K的条件下振动松弛时间与升温条件接近，但非平衡辐射特征时间在10-6s量级，明显低于升温条件。

图12 由T0=10 000 K或T0=30 000 K升温或降温至T=20 000 K后辐射总强度随时间变化Fig.12 Time evolution of total radiative intensity(T0=10 000 K or T0=30 000 K, T=20 000 K)

图13 由T0=10 000 K或T0=30 000 K升温或降温至T=20 000 K后O2和O数密度随时间变化Fig.13 Time evolution of number densities of O2 and O(T0=10 000 K or T0=30 000 K, T=20 000 K)

图14 由T0=10 000 K或T0=30 000 K升温或降温至T=20 000 K后振动温度随时间变化Fig.14 Time evolution of vibrational temperature (T0=10 000 K or T0=30 000 K, T=20 000 K)

从图12看出，无论升温条件还是降温条件，I和Ieq都不是单调趋近Ieq,B的。以T0=30 000 K的降温条件为例，在10-9s～10-7s区间，I降至低于Ieq,B，之后迅速回升，再渐近趋近Ieq,B。这反映出在该时间段内大量基态电子能级上的分子吸收光子向激发态电子能级跃迁，超越了设定温度对应的平衡辐射跃迁数，导致激发态电子能级的粒子数开始超过设定温度下的平衡分布粒子数、之后再回迁的特点。Ieq在10-9s～10-7s区间的下降程度更大，这与图14体现出的该时间段内振动温度(由非平衡的振动能级分布确定)的大幅变化是一致的。I和Ieq的差异体现了辐射本身的非平衡特点(I不能像Ieq那样直接由分子的能级分布确定)，Ieq与Ieq,B的差异则体现了真实的非平衡分子能级分布与设定温度下玻尔兹曼分布的差异。I与Ieq在10-9s～10-7s时间段内的明显差异说明真实的辐射跃迁数远低于根据能级分布确定的辐射跃迁数，辐射的非平衡程度很强。

图15展示了由T0=30 000 K降温至T=20 000 K条件下，在不同时刻各个跃迁或反应机制以及辐射对各振动能级数密度变化的贡献，即式(1)右端的各生成项。图中给出各项的绝对值，其中DRa和DRm的贡献为负(算例条件下是发生离解反应，参见图13)，其他各项采用“+”、“-”表示正负。由于初始辐射强度设定为T0=30 000 K下的黑体平衡辐射强度，所以降温条件下系统是吸收辐射，基态电子能级各振动能级上的粒子数是减小的，rad过程的贡献也为负，到10-7s后rad过程贡献变为正。由图15可见在t<10-9s时，离解反应尤其是DRa过程起主导作用。到了10-8s时，DR反应的生成项减小，VT过程生成项的贡献逐渐凸显。图14中t=10-8s～10-7s时间段内振动温度的迅速回升就主要是VT过程贡献。到10-7s时VT和DR过程的生成项迅速下降，比10-8s时低4个数量级，由图14也可知此时振动温度基本达到设定温度。rad过程生成项在10-10s时比VT过程低一个量级，而到10-9～10-8s时则比VT过程低了2～3个数量级，但其后rad生成项量级变化不大，达到了与DRa、VTa项相同的数量级。该算例条件下，VV过程的生成项一直远低于其他机制。

(a) t=10-10s

3.4 计算耗时统计

表1 各算例计算耗时统计Table 1 Calculation time cost of each case

由于反应速率随设定温度升高而升高，高温下系统较快达到平衡，而温度降低后达到平衡所需时间明显变长。本文含辐射算例(T=20 000 K)的τeq在10-6s～10-5s量级，不含辐射的算例中，T=20 000 K的τeq在10-6s量级，而T=3000 K的τeq则达到了10-1s量级。可能由于Euler隐式方法解决数值计算刚性问题的效果不够好，即使对T=3000 K算例，起始的时间步长也只能取为Δt=10-10s，在推进至10-2s后再增大时间步长以减小计算消耗。后续研究应尝试更有效的VODE方法等。

不包含辐射时，T=3000 K和T=10 000 K算例一个时间步的计算时间约0.5 ms，但T=20 000 K算例的计算时间仅0.2 ms左右，这与不同温度下刚性问题程度不同因而需要的Euler隐式内迭代次数不同有关。T=20 000 K条件包含辐射后，由于求解方程增多，每个时间步的计算时间增大至0.6 ms左右。可见态-态模型的计算消耗确实很大，但未来随着计算机性能的进一步提高和发展更先进的计算技术，将态-态模型用到二维非平衡流动模拟(如文献[24])，或者是三维流场中非平衡程度特别强的局部区域采用态-态模型，也是具有工程可行性的。

4 结 论

本文对封闭静止的O2/O系统采用态-态模型模拟包含辐射的热化学非平衡过渡过程，得到以下结论：

1) 根据非平衡振动能级分布导出的振动温度随时间变化情况与双温度模型所描述的按指数律逼近平衡值的特点不同，有些条件下还出现了振动温度随时间并非单调变化的现象。各算例条件下的非平衡振动能级分布均不满足当时振动温度下的玻尔兹曼分布。

2) 当温度较低时，化学反应特征时间远大于振动松弛时间，振动非平衡过程对化学反应的影响不强。而当温度较高(T>10 000 K)时，化学反应和振动松弛过程同时发生，同时对振动能级分布演化产生影响。这时双温度模型不能很好反映粒子的非平衡振动能级分布，也就不能很好处理化学反应和振动松弛过程的耦合。

3) 本文算例显示出O2/O系统在高温条件下的非平衡辐射特征，不过辐射跃迁对热化学非平衡过程的影响不明显。非平衡过程初期，粒子能级分布变化率中辐射生成项的贡献远低于平动-振动跃迁。非平衡辐射的特征时间也远大于振动松弛时间，非平衡辐射特征时间在由T0=10 000 K升温至设定T=20 000 K的算例中为10-5s量级，由T0=30 000 K降温至T=20 000 K的条件下为10-6s量级；而振动松弛时间在升温和降温条件下均在10-7～10-6s量级。

4) 包含各类跃迁机制的非平衡过程数值计算中刚性问题严重，本文采用Euler隐式算法求解时要求极小的时间推进步长，限制了计算效率。未来研究中需选用对刚性问题更加高效的VODE方法或α-QSS方法。