鄭偉杰, 曾 明, 王東方, 謝 易

(國防科技大學(xué) 空天科學(xué)學(xué)院, 長沙 410073)

0 引 言

在高超聲速飛行條件下，流場中的高溫導(dǎo)致氣體各個內(nèi)能模式激發(fā)并且發(fā)生化學(xué)反應(yīng)，而高空條件氣體密度低又使流場處于熱力和化學(xué)非平衡狀態(tài)。當(dāng)溫度足夠高時，氣體進(jìn)一步發(fā)生輻射躍遷。各類微觀反應(yīng)和躍遷機(jī)制影響到飛行器的宏觀氣動熱和氣動物理特性。更深入地理解和更精確地模擬熱化學(xué)非平衡和輻射現(xiàn)象在學(xué)術(shù)上和工程應(yīng)用上都具有重要意義。

有兩類描述流動中振動非平衡和振動-化學(xué)耦合的方法。一類是多溫度模型法，以Park在1985年提出的雙溫度模型[1]為代表，引入不同于平動溫度的振動溫度Tvib表征非平衡振動能，在化學(xué)反應(yīng)控制溫度中引入振動溫度反映熱力-化學(xué)耦合。另一類是態(tài)-態(tài)模型法[2-3]，直接描述各振動能級間的分子躍遷過程和各能級分子的化學(xué)反應(yīng)，求解各能級粒子數(shù)變化率方程。

雙溫度模型及后來拓展的多溫度模型，計(jì)算消耗小，且一定程度上反映了振動非平衡與振動-化學(xué)耦合，獲得了廣泛應(yīng)用。許多氣動熱力學(xué)程序采用了該模型，包括Gnoff等發(fā)展的LAURA (Langely Aerothermodynamic Upwind Relaxation Algorithm)程序[4-5]和GASP(General Aerodynamic Simulation Program)程序[6-7]。采用該模型的數(shù)值模擬工作還有很多，例如：Tissera[8]等進(jìn)行的高超聲速繞鈍錐-柱-裙的流動模擬，Sorensen等[9]進(jìn)行的有限催化壁條件下的高超聲速圓柱繞流模擬，美國和歐洲研究者們獨(dú)立進(jìn)行的為評估CFD方法預(yù)測高超聲速條件下激波相互作用能力的一系列數(shù)值計(jì)算[10]，汪球[11]、周凱[12]等進(jìn)行的高焓風(fēng)洞噴管非平衡流場計(jì)算，等。雙溫度模型在許多條件下獲得了良好的數(shù)值結(jié)果，但也存在局限[13]。在溫度很高條件下以及降溫復(fù)合反應(yīng)過程中，可能會產(chǎn)生較大誤差，甚至出現(xiàn)模型失真[14]。雙溫度模型隱含著振動能級分布滿足振動溫度下玻爾茲曼(Boltzmann)分布的假設(shè)[15]，而態(tài)-態(tài)模型方法的結(jié)果表明該假設(shè)在許多條件下并不成立。

態(tài)-態(tài)模型法將粒子在能級間的躍遷與化學(xué)反應(yīng)視為同一性質(zhì)的過程研究，描述振動非平衡及其與化學(xué)反應(yīng)耦合的方法則更為精細(xì)和直接[16]。它直接研究各能級上粒子的化學(xué)反應(yīng)，避免了多溫度模型中經(jīng)驗(yàn)性的熱力-化學(xué)耦合方法。隨著分子動力學(xué)的研究和計(jì)算機(jī)計(jì)算能力的發(fā)展，態(tài)-態(tài)模型作為一種高保真度模型逐漸進(jìn)入人們的視野并在近十年來快速發(fā)展起來。采用態(tài)-態(tài)模型的研究由不涉及流動的零維問題[17]，逐步發(fā)展到正激波后、邊界層[18-21]和準(zhǔn)一維噴管[22]非平衡流動問題。徐丹對N2/N系統(tǒng)采用態(tài)-態(tài)模型開展了零維問題研究[23]，還數(shù)值求解了耦合態(tài)-態(tài)模型的N2/N準(zhǔn)一維噴管非平衡流方程[16]。Druguet[24]采用態(tài)-態(tài)模型(考慮N2基態(tài)電子能級上68個振動能級)求解了N2/N的球頭非平衡流場，考察態(tài)-態(tài)模型與流動耦合的可行性。應(yīng)用CFD程序PINENS(并行隱式非平衡N-S方程求解程序)求解了包含69個組元連續(xù)方程的軸對稱非平衡流N-S方程，結(jié)果表明僅在邊界層外振動能級分布滿足振動溫度下的玻爾茲曼分布。

本文對封閉靜止的O2/O氣體系統(tǒng)，采用態(tài)-態(tài)模型研究熱化學(xué)非平衡和輻射過程。考察不同初始條件和設(shè)定溫度條件下(溫度由初始值瞬時改變?yōu)樵O(shè)定值后保持系統(tǒng)等溫等容)不同振動能級粒子分布隨時間的演化特點(diǎn)，分析比較各類躍遷或化學(xué)反應(yīng)過程以及輻射躍遷過程對系統(tǒng)化學(xué)組成和分子振動能級分布變化的貢獻(xiàn)。為更好地理解O2/O和高溫空氣熱化學(xué)反應(yīng)機(jī)理提供新的思路，也為今后采用態(tài)-態(tài)模型研究高溫空氣的非平衡過程并進(jìn)一步與流動耦合奠定基礎(chǔ)。

1 態(tài)-態(tài)模型基本方程

1.1 采用態(tài)-態(tài)模型的粒子數(shù)密度變化率方程

(1)

式中i=0～46。式(1)右端分別為振動-平動躍遷、振動-振動躍遷、離解-復(fù)合化學(xué)反應(yīng)和輻射幾類微觀過程對基態(tài)電子能級的i振動能級上O2分子數(shù)密度變化的貢獻(xiàn)。熱化學(xué)反應(yīng)過程將在1.2節(jié)詳細(xì)描述，輻射躍遷過程在1.3節(jié)詳細(xì)描述。

O原子數(shù)密度nO的變化僅起因于離解-復(fù)合化學(xué)反應(yīng)：

(2)

(3)

則,

(4)

(5)

gX和gB分別為基態(tài)和激發(fā)態(tài)電子能級的簡并度，εel,X和εel,B分別為基態(tài)和激發(fā)態(tài)電子能級的能值，kB為玻爾茲曼常數(shù)。激發(fā)態(tài)電子能級上各振動能級的分布則根據(jù)振動溫度Tvib下的玻爾茲曼分布確定。系統(tǒng)中總的O2分子數(shù)密度為:

(6)

1.2 不涉及輻射的能級躍遷與化學(xué)反應(yīng)過程

不涉及輻射的微觀過程包括：分子間的振動-平動能量交換過程VTm，分子與原子間的振動-平動能量交換過程VTa，振動-振動能量交換過程VV，原子作為碰撞參與者的離解-復(fù)合反應(yīng)DRa，分子作為碰撞參與者的過程DRm。具體如下：

(7)

(8)

(9)

(10)

(11)

對VTm過程只考慮單量子數(shù)躍遷，對VTa過程考慮多量子數(shù)躍遷，Δv為躍遷量子數(shù)。VV過程中，i和j代表O2基態(tài)電子能級上的不同振動能級。

(12)

式(12)右端當(dāng)i=0時僅含第1和第3項(xiàng)，i=46時僅含第2和第4項(xiàng)。

(13)

(14)

(15)

(16)

O原子數(shù)密度變化方程式(2)右端的兩項(xiàng)生成項(xiàng)為:

(17)

(18)

1.3 輻射躍遷過程

輻射躍遷包括三類——自發(fā)發(fā)射、誘導(dǎo)發(fā)射、輻射吸收。具體表示如下：

(19)

(20)

(21)

(22)

式(22)中c為光速。Aα′i′,αi、Bα′i′,αi、Bαi,α′i′分別為自發(fā)發(fā)射、誘導(dǎo)發(fā)射和輻射吸收的愛因斯坦系數(shù)。根據(jù)細(xì)致平衡原理，它們之間滿足如下關(guān)系：

(23)

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

(24)

(25)

(26)

本文對三種輻射躍遷機(jī)制的譜線均采用多普勒增寬的線型分布，譜線互不重疊。將式(22)對頻率積分后得：

(27)

(28)

1.4 振動溫度的確定

1.4.1 根據(jù)非平衡振動能級分布確定振動溫度

一種方法是通過假設(shè)非平衡的各振動能級分布盡可能接近振動溫度下的玻爾茲曼分布導(dǎo)出一個統(tǒng)一的振動溫度。需要采用線性回歸的方法。

振動溫度Tvib下的玻爾茲曼分布為:

(29)

式(29)中εvib,i為位于振動能級i的一個分子具有的振動能。Qvib為振動配分函數(shù),

(30)

對式(29)兩邊取對數(shù)，得:

(31)

近似將上式中的振動配分函數(shù)Qvib視為常數(shù)，則i能級粒子數(shù)密度的對數(shù)便為該能級能值的線性函數(shù)，記作:

(32)

其中K1為該線性函數(shù)的斜率,

(33)

根據(jù)已知求得的各振動能級非平衡粒子數(shù)密度，采用線性回歸方法確定式(32)的斜率，即得到對應(yīng)的振動溫度Tvib。

(34)

式(34)中wi為權(quán)重因子，根據(jù)該能級的粒子占比確定：

(35)

(36)

(37)

1.4.2 根據(jù)平均的非平衡振動能確定振動溫度

(38)

2 計(jì)算方法

由于各類躍遷機(jī)制的速率差別大，導(dǎo)致常微分方程組的數(shù)值求解存在嚴(yán)重剛性問題，可以采用Euler隱式方法求解。另外，VODE(變系數(shù)的常微分方程求解器)[37]是變系數(shù)的BDF(向后分化公式)方法，可以高效地解決具有剛性的問題。α-QSS(擬穩(wěn)態(tài)逼近方法)[38]是具有二階精度的預(yù)測-校驗(yàn)方法，常用于反應(yīng)速率系數(shù)隨時間緩變的多組元化學(xué)反應(yīng)模擬。考慮到Euler隱式算法形式簡單、應(yīng)用方便，對本文算例從計(jì)算精度和計(jì)算消耗上可以接受，本文采用了Euler隱式算法。

(39)

其中，

(40)

(41)

由于系統(tǒng)在初始瞬間改變溫度后保持等溫等容，各生成項(xiàng)中的速率常數(shù)在起始時刻根據(jù)系統(tǒng)設(shè)定溫度計(jì)算一次即可。時間推進(jìn)的步長視具體溫度和壓力條件，可在10-11～10-8s之間取值。本文涉及的不耦合輻射的算例條件(設(shè)定溫度分別為3000 K、10 000 K、20 000 K)下，時間步長均取為10-10s。包含輻射機(jī)制時，由于輻射強(qiáng)度隨時間變化率較大，時間推進(jìn)步長需進(jìn)一步減小至10-14～10-13s，本文含輻射的算例條件(設(shè)定溫度20 000 K)時間步長均取為10-13s。

3 算例及計(jì)算結(jié)果分析

3.1 態(tài)-態(tài)模型熱化學(xué)非平衡程序校驗(yàn)

首先根據(jù)文獻(xiàn)[35]計(jì)算結(jié)果對本文編制的O2/O態(tài)-態(tài)模型熱化學(xué)非平衡程序進(jìn)行校驗(yàn)。文獻(xiàn)[35]給出了兩組算例，初始的O2/O數(shù)密度均為nO2,0=1×1017/cm3、nO,0=9×1017/cm3，設(shè)定的系統(tǒng)溫度分別為T=3000 K和T=10 000 K。兩組算例分別設(shè)定不同的初始溫度：T=3000 K組的初始溫度分別為T0=100 K和T0=10 000 K;T=10 000 K組的初始溫度分別為T0=100 K和T0=20 000 K。模擬了兩組系統(tǒng)在起始時刻突然升溫或降溫并保持等溫等容條件下的熱化學(xué)非平衡過渡過程。

圖1、圖2分別給出了兩組算例系統(tǒng)的振動溫度隨時間的變化，振動溫度為應(yīng)用式(38)根據(jù)基態(tài)電子能級上振動能級分布導(dǎo)出的Tvib。圖中“dv”代表VTa過程中的振動躍遷量子數(shù)Δv，all代表允許Δv取任意可能值?？梢姳疚挠?jì)算結(jié)果與文獻(xiàn)[35]的計(jì)算結(jié)果一致性良好，驗(yàn)證了本文方法和程序的正確性。

圖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均可發(fā)現(xiàn)限制躍遷量子數(shù)Δv的影響，可見VTa過程中的多量子數(shù)躍遷的重要性。不過Δv<6與不限制Δv的結(jié)果差異就不大了。無論是升溫還是降溫至T=3000 K條件，隨著時間推進(jìn)，振動溫度都是近似以指數(shù)方式趨近于平動溫度，約在10-6s時與平動溫度相等，基于振動溫度的定義，可以說振動能和平動/轉(zhuǎn)動能此時基本達(dá)到了平衡。但設(shè)定溫度T=10 000 K條件的情況中振動溫度隨時間的變化特點(diǎn)不同，見圖2。初始T0=100 K條件下，振動溫度在10-7s升至6000 K附近后有一段平臺區(qū)，到10-6s后又開始迅速爬升，達(dá)到10 000 K，整體呈現(xiàn)了臺階形狀。初始T0=20 000 K條件下，振動溫度在10-7s前一直下降，達(dá)到6000 K左右，之后和T0=100 K條件的變化曲線幾乎重合，經(jīng)歷平臺區(qū)后再迅速爬升至T=10 000 K。振動溫度在10-7～10-6s期間的平臺區(qū)與此時開始的化學(xué)反應(yīng)有關(guān)，后面再進(jìn)一步分析。

3.2 O2/O系統(tǒng)的熱化學(xué)非平衡過程分析

圖3分別給出了由T0=100 K或T0=10 000 K升溫或降溫至T=3000 K的兩個算例O2分子歸一化振動能級分布隨時間的變化曲線，圖中不同能級采用了不同顏色的曲線，符號標(biāo)識對VTa過程中振動躍遷量子數(shù)Δv的限制。這里和后面的振動能級分布均指基態(tài)電子能級上的振動能級分布。初始時刻振動能級分布為初始溫度下的平衡分布，因此由T0=100 K升溫的算例(見圖3a)中，隨著時間推進(jìn)，低能級分子數(shù)減小，高能級分子數(shù)增大；低能級(i<10)分布在10-7s時就接近平衡，但i>20能級上的O2分子，在10-6～10-4s區(qū)間存在一個粒子數(shù)密度基本不變的平臺區(qū)，到10-2s后才接近平衡。由T0=10 000 K降溫的算例，隨時間推進(jìn)，高、低能級分子數(shù)變化的趨勢與由T0=100 K升溫的算例相反，但分布達(dá)到平衡的時間和過渡過程出現(xiàn)的平臺區(qū)與初始溫度為T0=100 K的情況類似。這是因?yàn)檎駝榆S遷和化學(xué)反應(yīng)速率均取決于設(shè)定的系統(tǒng)溫度(T=3000 K)。兩個算例中都體現(xiàn)了限制躍遷量子數(shù)的影響，特別是對高振動能級，多量子數(shù)的躍遷更加重要。與圖1給出的振動溫度隨時間變化曲線比較，可知高振動能級分布達(dá)到平衡的松弛時間比振動溫度的長得多，這正是高振動能級粒子分布不滿足振動溫度下的玻爾茲曼分布的體現(xiàn)。不過，由于高能級分子占比很小，其對總體振動能的影響也很小。

(a) T0=100 K，升溫

振動能級分布在非平衡過渡過程中的變化與化學(xué)反應(yīng)相關(guān)。這里首先以初始T0=10 000 K降溫至T=3000 K的算例為代表進(jìn)行分析。圖4給出了O2與O粒子數(shù)密度隨時間的變化，圖5對比了包含所有振動躍遷、化學(xué)反應(yīng)過程(VTm、VTa、VV、DRa、DRm)和去除化學(xué)反應(yīng)計(jì)算得到的振動能級分布隨時間的變化，圖中含所有振動躍遷和化學(xué)反應(yīng)過程的用實(shí)心符號，去除化學(xué)反應(yīng)的則用空心符號。圖6對比了包含或去除化學(xué)反應(yīng)計(jì)算得到的振動溫度變化曲線。圖7則分別給出了包含或去除化學(xué)反應(yīng)時振動能級分布與對應(yīng)振動溫度下的玻爾茲曼分布的對比。

圖4 由T0=10 000 K降溫至T=3000 K后O2和O數(shù)密度變化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后振動溫度隨時間變化(含或不含化學(xué)反應(yīng)的對比)Fig.6 Time evolution of vibrational temperature (with or without DR processes, T0=10 000 K, T=3000 K)

(a) 去除化學(xué)反應(yīng)

相對于算例設(shè)定的溫度T=3000 K下的平衡化學(xué)組成，初始時刻的nO,0/nO2,0偏大，非平衡過渡過程中將發(fā)生O原子的復(fù)合反應(yīng)?；瘜W(xué)反應(yīng)有一定的“孵化”時間，由圖4可見該設(shè)定溫度T=3000 K條件下大約10-6s復(fù)合反應(yīng)才起步，到10-5s后才有明顯的反應(yīng)，大約10-2s反應(yīng)達(dá)到平衡。結(jié)合圖5給出的有無化學(xué)反應(yīng)的振動能級分布對比，可以看到二者的高振動能級分布的明顯差異正是化學(xué)反應(yīng)開始的時間。可以分析得出這里的復(fù)合反應(yīng)生成的O2分子主要是在高振動能級。高振動能級的分子分布受到化學(xué)反應(yīng)的明顯影響，而低能級分布幾乎不受影響。10-6s～10-4s之間出現(xiàn)的高振動能級分子數(shù)密度的平臺區(qū)，是該區(qū)間復(fù)合生成高能級O2和高能級O2通過振動松弛向低能級躍遷互相平衡的結(jié)果。

由圖6可知算例條件下，是否包含化學(xué)反應(yīng)對振動溫度沒有明顯影響，這和高能級分子占比小因而對振動能貢獻(xiàn)小的特點(diǎn)是一致的。另外可看出以振動溫度衡量的振動松弛時間大約為10-6s，遠(yuǎn)小于以高振動能級分布衡量的松弛時間和離解-復(fù)合反應(yīng)達(dá)到平衡所需要的時間(約10-2s)。

從圖7可見，包含或去除化學(xué)反應(yīng)時振動能級分布均偏離了對應(yīng)振動溫度下的玻爾茲曼分布，包含化學(xué)反應(yīng)時偏離程度更大，特別是高振動能級的分布。從隨時間的變化來看，由圖7(b)可見t<10-7s時，化學(xué)反應(yīng)尚未發(fā)生，振動松弛占主導(dǎo)地位，位于高能級的O2分子數(shù)低于振動溫度下的玻爾茲曼分布。而t>10-6s后，復(fù)合反應(yīng)開始，生成大量的高能級O2分子，使得位于高能級的O2分子數(shù)密度高于此時振動溫度下的玻爾茲曼分布。

接下來對初始T0=100 K升溫至T=10 000 K的算例，分析化學(xué)組成隨時間的變化(見圖8)，振動能級分布和振動溫度的演化特點(diǎn)及其受化學(xué)反應(yīng)的影響(見圖9、圖10)，不同時刻振動能級分布與振動溫度下玻爾茲曼分布的差異(見圖11)。相對于算例設(shè)定的溫度T=10 000 K下的平衡化學(xué)組成，初始時刻的nO,0/nO2,0偏小，非平衡過渡過程中將發(fā)生O2分子的離解反應(yīng)。在設(shè)定溫度下，振動躍遷和化學(xué)反應(yīng)的速率都較T=3000 K算例條件下更快，松弛時間也大大減小。

圖8 由T0=100 K升溫至T=10 000 K后O2和O數(shù)密度變化Fig.8 Time evolution of number densities of O2 and O(T0=100 K, T=10 000 K)

由圖8～圖10可知，化學(xué)反應(yīng)達(dá)到平衡的特征時間、以振動溫度或高振動能級分布表征的松弛時間，都在10-7s～10-6s之間?；瘜W(xué)反應(yīng)對振動能級分布的影響，比T=3000 K算例條件下更大。離解反應(yīng)主要發(fā)生在10-8s～10-6s時間段內(nèi)，且存在高振動能級優(yōu)先離解的特點(diǎn)。由圖11(a)可見，不包含化學(xué)反應(yīng)時，高振動能級的O2分子數(shù)密度高于振動溫度下的玻爾茲曼分布。而包含化學(xué)反應(yīng)后，較高能級的分子優(yōu)先離解，導(dǎo)致高能級分子數(shù)低于振動溫度下的玻爾茲曼分布值，見圖11(b)。

(a) 低振動能級

圖10 由T0=100 K升溫至T=10 000 K后振動溫度隨時間變化(含或不含化學(xué)反應(yīng)的對比)Fig.10 Time evolution of vibrational temperature (with or without DR processes, T0=100 K, T=10 000 K)

(a) 去除化學(xué)反應(yīng)

對比分析T=3000 K和T=10 000 K算例可知，較低的設(shè)定溫度下，化學(xué)反應(yīng)所需時間比振動松弛時間大得多，兩個過程相對獨(dú)立，化學(xué)反應(yīng)對振動松弛產(chǎn)生的影響較小。而高設(shè)定溫度下，化學(xué)反應(yīng)特征時間和振動松弛時間可以相比擬，彼此耦合程度更高，互相影響更大。

3.3 O2/O系統(tǒng)的熱化學(xué)非平衡與輻射過程分析

選取設(shè)定溫度為T=20 000 K，初始的O2和O數(shù)密度均為nO2,0=1×1017/cm3、nO,0=9×1017/cm3，初始溫度分別為T0=10 000 K和T0=30 000 K的兩個算例條件，對O2/O系統(tǒng)研究高溫下耦合輻射的熱化學(xué)非平衡過渡過程，考察輻射生成項(xiàng)對粒子能級分布的貢獻(xiàn)。

首先研究總的輻射強(qiáng)度。將式(27)給出的針對躍遷α′i′→αi的輻射帶強(qiáng)度對所有的躍遷求和，即得到總的輻射強(qiáng)度隨時間的變化率方程。本文考慮的輻射躍遷包括由激發(fā)態(tài)電子能級上振動能級i′=0～14至基態(tài)電子能級上振動能級i=0～21的共計(jì)330個輻射譜帶躍遷?？偟妮椛鋸?qiáng)度記為I。

令式(27)為0，得到平衡輻射譜帶強(qiáng)度：

(42)

圖12給出了由T0=30 000 K或T0=10 000 K降溫或升溫至T=20 000 K后三個總輻射強(qiáng)度I、Ieq、Ieq,B輻射強(qiáng)度隨時間的變化。為分析輻射與熱化學(xué)非平衡過程的相互影響特點(diǎn)，圖13、圖14還給出了相應(yīng)的O2、O數(shù)密度和振動溫度的時間演化曲線。由圖12可見隨著時間的推進(jìn)，輻射強(qiáng)度I和平衡輻射強(qiáng)度Ieq最終將達(dá)到設(shè)定溫度下的黑體平衡輻射強(qiáng)度Ieq,B?？蓪⑤椛鋸?qiáng)度I達(dá)到黑體平衡輻射強(qiáng)度Ieq,B的時間定義為非平衡輻射特征時間，則由T0=10 000 K升溫至設(shè)定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后輻射總強(qiáng)度隨時間變化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數(shù)密度隨時間變化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都不是單調(diào)趨近Ieq,B的。以T0=30 000 K的降溫條件為例，在10-9s～10-7s區(qū)間，I降至低于Ieq,B，之后迅速回升，再漸近趨近Ieq,B。這反映出在該時間段內(nèi)大量基態(tài)電子能級上的分子吸收光子向激發(fā)態(tài)電子能級躍遷，超越了設(shè)定溫度對應(yīng)的平衡輻射躍遷數(shù)，導(dǎo)致激發(fā)態(tài)電子能級的粒子數(shù)開始超過設(shè)定溫度下的平衡分布粒子數(shù)、之后再回遷的特點(diǎn)。Ieq在10-9s～10-7s區(qū)間的下降程度更大，這與圖14體現(xiàn)出的該時間段內(nèi)振動溫度(由非平衡的振動能級分布確定)的大幅變化是一致的。I和Ieq的差異體現(xiàn)了輻射本身的非平衡特點(diǎn)(I不能像Ieq那樣直接由分子的能級分布確定)，Ieq與Ieq,B的差異則體現(xiàn)了真實(shí)的非平衡分子能級分布與設(shè)定溫度下玻爾茲曼分布的差異。I與Ieq在10-9s～10-7s時間段內(nèi)的明顯差異說明真實(shí)的輻射躍遷數(shù)遠(yuǎn)低于根據(jù)能級分布確定的輻射躍遷數(shù)，輻射的非平衡程度很強(qiáng)。

圖15展示了由T0=30 000 K降溫至T=20 000 K條件下，在不同時刻各個躍遷或反應(yīng)機(jī)制以及輻射對各振動能級數(shù)密度變化的貢獻(xiàn)，即式(1)右端的各生成項(xiàng)。圖中給出各項(xiàng)的絕對值，其中DRa和DRm的貢獻(xiàn)為負(fù)(算例條件下是發(fā)生離解反應(yīng)，參見圖13)，其他各項(xiàng)采用“+”、“-”表示正負(fù)。由于初始輻射強(qiáng)度設(shè)定為T0=30 000 K下的黑體平衡輻射強(qiáng)度，所以降溫條件下系統(tǒng)是吸收輻射，基態(tài)電子能級各振動能級上的粒子數(shù)是減小的，rad過程的貢獻(xiàn)也為負(fù)，到10-7s后rad過程貢獻(xiàn)變?yōu)檎?。由圖15可見在t<10-9s時，離解反應(yīng)尤其是DRa過程起主導(dǎo)作用。到了10-8s時，DR反應(yīng)的生成項(xiàng)減小，VT過程生成項(xiàng)的貢獻(xiàn)逐漸凸顯。圖14中t=10-8s～10-7s時間段內(nèi)振動溫度的迅速回升就主要是VT過程貢獻(xiàn)。到10-7s時VT和DR過程的生成項(xiàng)迅速下降，比10-8s時低4個數(shù)量級，由圖14也可知此時振動溫度基本達(dá)到設(shè)定溫度。rad過程生成項(xiàng)在10-10s時比VT過程低一個量級，而到10-9～10-8s時則比VT過程低了2～3個數(shù)量級，但其后rad生成項(xiàng)量級變化不大，達(dá)到了與DRa、VTa項(xiàng)相同的數(shù)量級。該算例條件下，VV過程的生成項(xiàng)一直遠(yuǎn)低于其他機(jī)制。

(a) t=10-10s

3.4 計(jì)算耗時統(tǒng)計(jì)

表1 各算例計(jì)算耗時統(tǒng)計(jì)Table 1 Calculation time cost of each case

由于反應(yīng)速率隨設(shè)定溫度升高而升高，高溫下系統(tǒng)較快達(dá)到平衡，而溫度降低后達(dá)到平衡所需時間明顯變長。本文含輻射算例(T=20 000 K)的τeq在10-6s～10-5s量級，不含輻射的算例中，T=20 000 K的τeq在10-6s量級，而T=3000 K的τeq則達(dá)到了10-1s量級?？赡苡捎贓uler隱式方法解決數(shù)值計(jì)算剛性問題的效果不夠好，即使對T=3000 K算例，起始的時間步長也只能取為Δt=10-10s，在推進(jìn)至10-2s后再增大時間步長以減小計(jì)算消耗。后續(xù)研究應(yīng)嘗試更有效的VODE方法等。

不包含輻射時，T=3000 K和T=10 000 K算例一個時間步的計(jì)算時間約0.5 ms，但T=20 000 K算例的計(jì)算時間僅0.2 ms左右，這與不同溫度下剛性問題程度不同因而需要的Euler隱式內(nèi)迭代次數(shù)不同有關(guān)。T=20 000 K條件包含輻射后，由于求解方程增多，每個時間步的計(jì)算時間增大至0.6 ms左右?？梢姂B(tài)-態(tài)模型的計(jì)算消耗確實(shí)很大，但未來隨著計(jì)算機(jī)性能的進(jìn)一步提高和發(fā)展更先進(jìn)的計(jì)算技術(shù)，將態(tài)-態(tài)模型用到二維非平衡流動模擬(如文獻(xiàn)[24])，或者是三維流場中非平衡程度特別強(qiáng)的局部區(qū)域采用態(tài)-態(tài)模型，也是具有工程可行性的。

4 結(jié) 論

本文對封閉靜止的O2/O系統(tǒng)采用態(tài)-態(tài)模型模擬包含輻射的熱化學(xué)非平衡過渡過程，得到以下結(jié)論：

1) 根據(jù)非平衡振動能級分布導(dǎo)出的振動溫度隨時間變化情況與雙溫度模型所描述的按指數(shù)律逼近平衡值的特點(diǎn)不同，有些條件下還出現(xiàn)了振動溫度隨時間并非單調(diào)變化的現(xiàn)象。各算例條件下的非平衡振動能級分布均不滿足當(dāng)時振動溫度下的玻爾茲曼分布。

2) 當(dāng)溫度較低時，化學(xué)反應(yīng)特征時間遠(yuǎn)大于振動松弛時間，振動非平衡過程對化學(xué)反應(yīng)的影響不強(qiáng)。而當(dāng)溫度較高(T>10 000 K)時，化學(xué)反應(yīng)和振動松弛過程同時發(fā)生，同時對振動能級分布演化產(chǎn)生影響。這時雙溫度模型不能很好反映粒子的非平衡振動能級分布，也就不能很好處理化學(xué)反應(yīng)和振動松弛過程的耦合。

3) 本文算例顯示出O2/O系統(tǒng)在高溫條件下的非平衡輻射特征，不過輻射躍遷對熱化學(xué)非平衡過程的影響不明顯。非平衡過程初期，粒子能級分布變化率中輻射生成項(xiàng)的貢獻(xiàn)遠(yuǎn)低于平動-振動躍遷。非平衡輻射的特征時間也遠(yuǎn)大于振動松弛時間，非平衡輻射特征時間在由T0=10 000 K升溫至設(shè)定T=20 000 K的算例中為10-5s量級，由T0=30 000 K降溫至T=20 000 K的條件下為10-6s量級；而振動松弛時間在升溫和降溫條件下均在10-7～10-6s量級。

4) 包含各類躍遷機(jī)制的非平衡過程數(shù)值計(jì)算中剛性問題嚴(yán)重，本文采用Euler隱式算法求解時要求極小的時間推進(jìn)步長，限制了計(jì)算效率。未來研究中需選用對剛性問題更加高效的VODE方法或α-QSS方法。