態-態模型下的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方法。