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

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

0 引 言

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

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

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

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

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

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

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

(1)

式中i=0～46。式(1)右端分別為振動(dòng)-平動(dòng)躍遷、振動(dòng)-振動(dòng)躍遷、離解-復(fù)合化學(xué)反應(yīng)和輻射幾類微觀過程對(duì)基態(tài)電子能級(jí)的i振動(dòng)能級(jí)上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)電子能級(jí)的簡并度，εel,X和εel,B分別為基態(tài)和激發(fā)態(tài)電子能級(jí)的能值，kB為玻爾茲曼常數(shù)。激發(fā)態(tài)電子能級(jí)上各振動(dòng)能級(jí)的分布則根據(jù)振動(dòng)溫度Tvib下的玻爾茲曼分布確定。系統(tǒng)中總的O2分子數(shù)密度為:

(6)

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

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

(7)

(8)

(9)

(10)

(11)

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

(12)

式(12)右端當(dāng)i=0時(shí)僅含第1和第3項(xiàng)，i=46時(shí)僅含第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)

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

(27)

(28)

1.4 振動(dòng)溫度的確定

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

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

振動(dòng)溫度Tvib下的玻爾茲曼分布為:

(29)

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

(30)

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

(31)

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

(32)

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

(33)

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

(34)

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

(35)

(36)

(37)

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

(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ù)隨時(shí)間緩變的多組元化學(xué)反應(yīng)模擬。考慮到Euler隱式算法形式簡單、應(yīng)用方便，對(duì)本文算例從計(jì)算精度和計(jì)算消耗上可以接受，本文采用了Euler隱式算法。

(39)

其中，

(40)

(41)

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

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

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

首先根據(jù)文獻(xiàn)[35]計(jì)算結(jié)果對(duì)本文編制的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)在起始時(shí)刻突然升溫或降溫并保持等溫等容條件下的熱化學(xué)非平衡過渡過程。

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

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

(a) T0=100 K，升溫

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

圖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) 低振動(dòng)能級(jí)

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

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

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

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

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

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

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

(a) 低振動(dòng)能級(jí)

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

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

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

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的兩個(gè)算例條件，對(duì)O2/O系統(tǒng)研究高溫下耦合輻射的熱化學(xué)非平衡過渡過程，考察輻射生成項(xiàng)對(duì)粒子能級(jí)分布的貢獻(xiàn)。

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

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

(42)

圖12給出了由T0=30 000 K或T0=10 000 K降溫或升溫至T=20 000 K后三個(gè)總輻射強(qiáng)度I、Ieq、Ieq,B輻射強(qiáng)度隨時(shí)間的變化。為分析輻射與熱化學(xué)非平衡過程的相互影響特點(diǎn)，圖13、圖14還給出了相應(yīng)的O2、O數(shù)密度和振動(dòng)溫度的時(shí)間演化曲線。由圖12可見隨著時(shí)間的推進(jìn)，輻射強(qiáng)度I和平衡輻射強(qiáng)度Ieq最終將達(dá)到設(shè)定溫度下的黑體平衡輻射強(qiáng)度Ieq,B。可將輻射強(qiáng)度I達(dá)到黑體平衡輻射強(qiáng)度Ieq,B的時(shí)間定義為非平衡輻射特征時(shí)間，則由T0=10 000 K升溫至設(shè)定T=20 000 K的算例中，非平衡輻射特征時(shí)間在10-5s量級(jí)，而振動(dòng)松弛時(shí)間(參見圖14)在10-7～10-6s量級(jí)。由T0=30 000 K降溫至T=20 000 K的條件下振動(dòng)松弛時(shí)間與升溫條件接近，但非平衡輻射特征時(shí)間在10-6s量級(jí)，明顯低于升溫條件。

圖12 由T0=10 000 K或T0=30 000 K升溫或降溫至T=20 000 K后輻射總強(qiáng)度隨時(shí)間變化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ù)密度隨時(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后振動(dòng)溫度隨時(shí)間變化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。這反映出在該時(shí)間段內(nèi)大量基態(tài)電子能級(jí)上的分子吸收光子向激發(fā)態(tài)電子能級(jí)躍遷，超越了設(shè)定溫度對(duì)應(yīng)的平衡輻射躍遷數(shù)，導(dǎo)致激發(fā)態(tài)電子能級(jí)的粒子數(shù)開始超過設(shè)定溫度下的平衡分布粒子數(shù)、之后再回遷的特點(diǎn)。Ieq在10-9s～10-7s區(qū)間的下降程度更大，這與圖14體現(xiàn)出的該時(shí)間段內(nèi)振動(dòng)溫度(由非平衡的振動(dòng)能級(jí)分布確定)的大幅變化是一致的。I和Ieq的差異體現(xiàn)了輻射本身的非平衡特點(diǎn)(I不能像Ieq那樣直接由分子的能級(jí)分布確定)，Ieq與Ieq,B的差異則體現(xiàn)了真實(shí)的非平衡分子能級(jí)分布與設(shè)定溫度下玻爾茲曼分布的差異。I與Ieq在10-9s～10-7s時(shí)間段內(nèi)的明顯差異說明真實(shí)的輻射躍遷數(shù)遠(yuǎn)低于根據(jù)能級(jí)分布確定的輻射躍遷數(shù)，輻射的非平衡程度很強(qiáng)。

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

(a) t=10-10s

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

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

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

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

4 結(jié) 論

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

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

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

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

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