高溫空氣熱化學(xué)物理計算模型研究

楊 光，檀妹靜，聶春生，李 宇，周 禹，李潞寧，于明星，顧云軍，陳其峰

(1.中國運載火箭技術(shù)研究院空間物理重點實驗室，北京 100076；2.中國工程物理研究院流體物理研究所沖擊波與爆轟波物理重點實驗室，綿陽 621999)

1 引 言

導(dǎo)彈、火箭、飛船等各類飛行器在大氣層內(nèi)以極高的速度飛行時，飛行器周圍的氣體受到劇烈壓縮而溫度急劇升高，高溫誘導(dǎo)的氣體分子振動、轉(zhuǎn)動、氧化、電荷交換、離解、電離、復(fù)合甚至多級電離等效應(yīng)將對空氣組分、熱力學(xué)和輸運性質(zhì)產(chǎn)生重要影響[1]，而空氣電離形成的等離子體鞘套也可能導(dǎo)致通信“黒障”的出現(xiàn).因此，準(zhǔn)確地計算高溫空氣的熱物性參數(shù)對于此類飛行器的設(shè)計有著至關(guān)重要的意義.

自20世紀(jì)50年代起至今，國外的研究者已在高溫氣體熱力學(xué)性質(zhì)研究方面進(jìn)行了大量的工作[2-9]，一些代表性工作如下：Drellishak等[2]計算了氬等離子體和N2-O2混合等離子體的配分函數(shù)和熱力學(xué)性質(zhì)，其配分函數(shù)的計算方法被廣泛應(yīng)用；Bacri等[3]計算了壓力為1～200 atm(1atm為1個大氣壓)及溫度為1000～30000 K范圍內(nèi)空氣等離子體的內(nèi)配分函數(shù)、組分和熱力學(xué)函數(shù)；Gupta等[4]利用11組分模型討論了空氣的反應(yīng)速率和熱力學(xué)性質(zhì).

國內(nèi)目前對于高溫空氣熱化學(xué)模型相關(guān)的工作以應(yīng)用為主，主要是基于平衡氣體模型和非平衡氣體模型進(jìn)行了高馬赫數(shù)流場和氣動熱環(huán)境的數(shù)值模擬[10，11]，而對于高溫空氣熱化學(xué)模型基礎(chǔ)理論和建模工作相對開展較少，僅有少量研究工作報道.王海興等人[12]開展了N2-O2空氣12組元雙溫度輸運性質(zhì)研究，陳其峰等人[13]計算了2000～30000 K溫度范圍內(nèi)N2-O2空氣的離解、電離和熱力學(xué)性質(zhì).另外，陳其峰等人也對空氣主要成分氮和氧在部分離解區(qū)的熱物性和離解化學(xué)反應(yīng)開展了理論建模研究[14，15]，對稀有氣體氦、氖、氬、氪和氙在部分電離區(qū)的熱物性及輸運性質(zhì)開展了系列性理論建模及強(qiáng)激波壓縮實驗研究[16-27]，這些工作主要是針對稠密等離子體物性的研究，沒有涉及高速飛行器周圍空氣的高溫稀薄條件，其關(guān)注點與本文工作不同.

盡管國內(nèi)外研究者已經(jīng)對空氣熱物性研究了幾十年，但在足夠高的溫度下，空氣可能發(fā)生的化學(xué)反應(yīng)有多少種?不同粒子的含量及其存在形式又是如何?應(yīng)該用怎樣的勢模型描述體系中相互作用機(jī)制?到目前為止，這些問題仍然沒有統(tǒng)一的結(jié)論，需要開展研究逐步解決.本文針對初始3組元(N2-O2-Ar)空氣，建立了從初始3組元到考慮18個化學(xué)反應(yīng)22組元混合體系熱動力學(xué)模型以及相應(yīng)的數(shù)值計算方法，對高溫空氣的組分濃度以及焓和比熱等熱力學(xué)性質(zhì)進(jìn)行了數(shù)值計算，對高溫空氣的離解電離等化學(xué)反應(yīng)以及焓和熱容等物性參數(shù)隨溫度和密度的演化規(guī)律進(jìn)行了研究，并預(yù)測了飛行器高速飛行時駐點處的空氣組分、壓縮比、溫度、熱容以及比熱等參數(shù)隨飛行速度的變化.

2 理論模型

2.1 平衡態(tài)組分

根據(jù)標(biāo)準(zhǔn)大氣模型，海平面附近空氣中N2、O2、Ar的體積分?jǐn)?shù)分別為78.084%、20.9476%、0.934%，三者相加大于99.9%，因此可以忽略其它組分的影響，按此體積分?jǐn)?shù)比例的N2-O2-Ar空氣作為計算初始條件.

在飛行器高速飛行推出強(qiáng)激波作用下，其周圍空氣主要發(fā)生高溫誘導(dǎo)的離解、電離和復(fù)合等化學(xué)反應(yīng)，這些化學(xué)反應(yīng)的發(fā)生以及相應(yīng)的空氣組分因密度和溫度差異而不同，例如，在中等溫度區(qū)間(2000～6000 K)，容易發(fā)生離解和復(fù)合反應(yīng)，形成NO、N、O、NO2、N2O等組分；而在高溫區(qū)(6000 K以上)，氮氧除了離解及復(fù)合反應(yīng)外，也將發(fā)生電離反應(yīng)，引入分子離子(N2+，O2+)、原子離子(N+，O+)、甚至因多級電離形成的高價離子等新組分.為了盡可能反映強(qiáng)激波作用下高溫空氣粒子組分隨溫度及密度的變化，在局域熱動平衡(LTE)和局域化學(xué)平衡(LCE)假定下，我們考慮了18個化學(xué)反應(yīng)以及由此產(chǎn)生的22組分混合體系進(jìn)行理論建模，其化學(xué)反應(yīng)通道(即反應(yīng)式)以及所考慮的組分見表1所列.

表1 22組分模型的化學(xué)反應(yīng)通道和平衡方程Table 1 Chemical reaction channels and equilibrium equations for the 22-species model

?

表中，λ為粒子/電子熱德布羅意波長，Z1(int)～Z21(int)為編號1～21組分的內(nèi)配分函數(shù).這些內(nèi)配分函數(shù)由組分的內(nèi)部自由度決定：對于分子和分子離子，其內(nèi)配分函數(shù)的計算需要考慮電子以及振動和轉(zhuǎn)動自由度的貢獻(xiàn)；對于原子和原子離子，則只需要考慮電子部分的貢獻(xiàn).

以O(shè)2和O為例，其內(nèi)配分函數(shù)的表達(dá)式為：

O原子內(nèi)部電子激發(fā)對原子內(nèi)配分函數(shù)的貢獻(xiàn)可以表示為：

O2分子的轉(zhuǎn)動配分函數(shù)可以表示為：

式中，Θr為轉(zhuǎn)動特征溫度，I為轉(zhuǎn)動慣量，σ為對稱數(shù)(對于同核分子，取值為2；對于異核分子，取值為1).

高溫近似下(即溫度T遠(yuǎn)大于Θr)，可以將(3)式中的求和用積分代替，即：

O2分子的振動配分函數(shù)，在諧振子近似下，可以通過求和得到：

式中，Θv為振動特征溫度，ω為振動圓頻率(由分子勢函數(shù)決定).

O2分子電子部分對內(nèi)配分函數(shù)的貢獻(xiàn)可以表示為：

式中，ε2i為O2分子內(nèi)部電子能級，g2i為對應(yīng)的簡并度.

對于當(dāng)前所建立22組分模型，利用化學(xué)平衡條件可以得到18個離解、電離及復(fù)合平衡方程，再結(jié)合約束方程(粒子數(shù)守恒方程和電荷守恒方程)，就可求解任意溫度和密度狀態(tài)下各組分的粒子數(shù)密度，獲得平衡態(tài)下空氣組分.

2.2 空氣熱物性模型

與配分函數(shù)的計算類似，從形式上看，可以將自由能分為原子(原子離子)和分子(分子離子)兩種類型.因此，我們選擇以O(shè)2+、O+為例對自由能的構(gòu)建方法進(jìn)行說明.

根據(jù)單組分體系自由能的形式可知，在粒子數(shù)密度、溫度已知的條件下，單位質(zhì)量的22組元混合體系中，O2+自由能函數(shù)的形式為：

式中，ρ為22組元混合體系的密度，可以通過各組分的粒子數(shù)密度ni和粒子質(zhì)量mi給出；E3，f為O2+的形成能.

同理，單位質(zhì)量的2組元混合體系中，O+自由能函數(shù)的形式為：

相比于O2+，由于缺少振動和轉(zhuǎn)動自由度，O+的自由能函數(shù)得到簡化.仿照以上公式，可以得出表中9種分子及分子離子(O2-、O2+、N2+、O2、N2、NO、NO+、N2O、NO2)和12種原子及原子離子(O-、O、O+、O2+、O3+、N、N+、N2+、N3+、Ar、Ar+、Ar2+)的自由能函數(shù).對以上21種組分自由能函數(shù)求和，再加上自由電子的自由能函數(shù)，即可得出22組分體系的自由能：

我們知道，自由能函數(shù)是關(guān)于溫度、體積(密度)和粒子數(shù)的特性函數(shù).因此，從理論角度，得到體系的自由能函數(shù)以后，根據(jù)自由能的偏微分就可以得到所有的熱力學(xué)函數(shù).以下列出了主要的熱力學(xué)函數(shù)(例如壓力p、熵S、內(nèi)能E、焓H、定容比熱容CV、定壓比熱容Cp)與自由能之間的關(guān)系：

除了選擇以溫度和密度為獨立變量進(jìn)行建模的方法外，還可以選擇以溫度和壓力作為獨立變量建立Gibbs自由能特性函數(shù)，通過對Gibbs自由能函數(shù)做偏微分得到熱力學(xué)函數(shù).根據(jù)熱力學(xué)基本理論，Gibbs自由能和Hemholtz自由能函數(shù)滿足以下關(guān)系：

因此，在計算得到Hemholtz自由能后，利用上式即可得到Gibbs自由能函數(shù)，然后再結(jié)合熱力學(xué)函數(shù)和Gibbs自由能函數(shù)之間的微分關(guān)系及麥克斯韋關(guān)系即可得到多組分體系的熱力學(xué)量.

需要說明的是，當(dāng)高溫電離出現(xiàn)時，為了得到真實氣體的物性參數(shù)，需要考慮帶電粒子間由于Coulomb長程作用而引起的附加自由能貢獻(xiàn)，即在總自由能中引入Debye-Hückel修正，其自由能附加項為：

在飛行器高速飛行推出強(qiáng)激波作用下，其駐點處的高溫空氣熱力學(xué)狀態(tài)滿足Hugoniot能量約束關(guān)系：

式中，下標(biāo)0和H分別代表飛行器駐點處空氣在激波前和激波后的狀態(tài)，激波前的狀態(tài)由初始條件決定.而激波前后的狀態(tài)又與飛行器的飛行速度Up相聯(lián)系，其可以通過下面公式表示：

飛行器駐點處空氣壓縮比(ρ/ρ0=V0/VH)是激波前后的密度(或比容)之比，這樣，根據(jù)方程4～12，再結(jié)合由空氣初始條件確定的激波前狀態(tài)，就可以得到駐點處激波后高溫空氣的熱物性參數(shù)和空氣組分隨飛行器飛行速度變化關(guān)系.

3 數(shù)值結(jié)果與分析

空氣組分是計算熱物性的首要條件，其與高溫誘導(dǎo)的空氣離解、電離、復(fù)合等化學(xué)反應(yīng)相聯(lián)系.為了了解高空高速飛行器周圍強(qiáng)激波壓縮空氣因溫度急劇升高引起的空氣化學(xué)反應(yīng)及組分變化，我們給出兩個典型數(shù)值模擬結(jié)果，圖1和圖2分別展示了以海拔高度1 km和60 km狀態(tài)空氣為初始條件，得到的中性粒子、分子離子、原子離子和電子濃度隨溫度的變化情況.從圖1可以看到，對于海拔高度為1.0 km(地表附近)，當(dāng)溫度達(dá)到約2000 K時，氧氣開始離解，且在約5000 K溫度時，氧原子濃度達(dá)到最大值；當(dāng)溫度達(dá)到3000 K時，氮氣開始離解，當(dāng)溫度上升到約10000 K時，氮原子濃度達(dá)到最大值；當(dāng)溫度為6500 K左右時，氮和氧的一級電離發(fā)生，N+、O+開始出現(xiàn)，并且其濃度隨著溫度的增大而增大；當(dāng)溫度為約8000 K時，氬的一級電離發(fā)生，Ar+開始出現(xiàn)；當(dāng)溫度約18000 K時，氮的二級電離發(fā)生，N2+開始出現(xiàn)；當(dāng)溫度大于20000 K時，氬和氧的二級電離發(fā)生，Ar2+、O2+開始出現(xiàn).由此可見，隨著溫度不斷升高，空氣離解、電離甚至多級電離等化學(xué)反應(yīng)逐漸出現(xiàn)，新組分逐步生成，空氣組分不斷變化，最終形成多元混合等離子體.需要說明的是，在我們考慮的溫度密度范圍，O2-離子的濃度很小，基本可以忽略.O2+、N2+、NO+在較大的溫度區(qū)間內(nèi)均存在，但是其粒子濃度的最大值均小于1.0×10-3.從圖2可以看出，海拔為60 km時，空氣組分隨溫度的變化也表現(xiàn)出與海拔1 km時類似的趨勢，其化學(xué)反應(yīng)通道以及對應(yīng)的組分也是隨溫度升高逐步出現(xiàn)的.所不同的是，相對于海拔1 km而言，海拔60 km對應(yīng)更為稀薄的初始狀態(tài)，空氣離解電離的特征溫度(即離解電離出現(xiàn)時所對應(yīng)溫度)均有所降低，氧離解特征溫度降到2000 K以下，氧原子濃度在約3000 K時達(dá)到最大值；氮離解特征溫度降到3000 K以下，氮原子濃度在溫度6000～9000 K達(dá)到最大值并出現(xiàn)“平臺區(qū)”；氮氧一級電離特征溫度降到5000 K左右，氮離子和氧離子的濃度在15000-20000 K達(dá)到最大值并出現(xiàn)“平臺區(qū)”；氮、氧和氬的二級電離特征溫度降到15000 K左右；在溫度25000 K左右，氮和氧的三級電離也出現(xiàn).上述結(jié)果表明，相對于海拔1 km來說，隨著海拔高度的增加，空氣稀薄程度也隨之增加，氣體密度減小，高空稀薄空氣中的各類化學(xué)反應(yīng)能夠在相對更低的溫度下發(fā)生，也就是說，空氣稀薄條件的增加有利于化學(xué)反應(yīng)的發(fā)生.

圖1 海拔1 km高度空氣組分隨溫度變化Fig.1 Air components vary with temperature at altitude of?=1 km

圖2 海拔60 km高度空氣組分隨溫度變化Fig.2 Air components vary with temperature at altitude of?=60 km

空氣組分一旦確定，就可以根據(jù)2.2節(jié)描述的熱力學(xué)量計算方法，計算空氣的熱物性參數(shù)，如焓、定容熱容、定壓熱容以及比熱比等.為了與現(xiàn)有國際流行Gupta模型比較，圖3展示壓力為10-4～102atm、溫度為500～30000 K范圍內(nèi)我們模型計算的空氣混合體系焓值及定壓熱容與Gupta等人[4]11組元模型(NASA-RP-1260，NASA是美國國家航空和宇宙航行局National Aeronautics and Space Administration的簡稱)計算結(jié)果的比較.可以看到，本文模型在三種不同壓力下(10-4、1和100 atm)計算的焓值與NASA11組元模型所得結(jié)果較好符合，在兩種不同壓力下(1和100 atm)計算的定壓熱容也與NASA11組元模型結(jié)果基本一致.但對于壓力10-4atm，本文計算的定壓熱容在約23000 K以上卻與NASA11組元結(jié)果出現(xiàn)明顯偏差.這種偏差出現(xiàn)的原因在于：10-4atm這樣的極低壓力有利于電離反應(yīng)的發(fā)生，在23000 K以上，三級電離已經(jīng)發(fā)生，而NASA11組元模型僅僅考慮到一級電離，盡管二、三級電離的出現(xiàn)還不足以對焓產(chǎn)生大的影響，但對由焓求偏導(dǎo)數(shù)得到的定壓熱容的影響卻已經(jīng)不能忽略.

圖3 與NASA 11組分模型所得焓值與定壓熱容的比較Fig.3 Comparison of enthalpy and constant pressure heat capacity with NASA 11-species model

定壓與定容熱容是表征系統(tǒng)的熱量交換的物理量，也在一定程度上與空氣的離解電離化學(xué)反應(yīng)相聯(lián)系.為了研究空氣稀薄程度對這兩個熱物性參數(shù)的影響，圖4給出了我們模型計算海拔1 km和60 km高度下大氣定容和定壓熱容隨溫度的變化.可以看到，這兩個熱容均隨溫度增加表現(xiàn)出波動變化，經(jīng)分析發(fā)現(xiàn)，這種波動變化與空氣發(fā)生離解電離等化學(xué)反應(yīng)以及新組分的形成有關(guān).對于海拔1.0 km高度，定容比熱容出現(xiàn)的3個波峰分別對應(yīng)于氧分子的離解、氮分子的離解以及氮和氧原子的一級電離，由于N+、O+濃度隨溫度變化趨勢基本一致，因此，在定容熱容上只對應(yīng)一個波峰；對于海拔60 km高度，定容熱容的波動變化更為劇烈，出現(xiàn)了四個波峰，它們分別對應(yīng)于氧分子的離解、氮分子的離解、氮和氧原子的一級電離、以及氮原子的二級電離.對于定壓熱容來說，其隨溫度和海拔高度的變化表現(xiàn)出與定容熱容相類似的變化趨勢，但是相同溫度下的定壓熱容總是大于定容熱容.通過對比研究不同海拔高度稀薄空氣熱容變化趨勢可以看出，隨著海拔高度的增加，當(dāng)空氣越稀薄，高溫?zé)崛莶▌幼兓驮絼×?，這充分反映了能量轉(zhuǎn)化變化過程.

圖4 不同海拔高度下定容熱容與定壓熱容隨溫度的變化Fig.4 The constant volume heat capacities and constant pressure heat capacity versus temperature at different altitudes.

眾所周知，比熱比(定壓熱容與定容熱容之比)能夠更加細(xì)致地反映熱量交換波動過程.與不考慮化學(xué)反應(yīng)的理想氣體不同，高溫真實氣體的比熱比不僅是溫度的函數(shù)，同時也與密度有關(guān).為了研究不同海拔高度下比熱比隨溫度的演化規(guī)律，圖5分別給出了海拔高度1 km和60 km典型比熱比隨溫度變化的模擬結(jié)果.可以看到，常溫下，氮分子和氧分子振動自由度被凍結(jié)，平動和轉(zhuǎn)動自由度完全激發(fā)，此時，比熱比的值為1.40；隨著溫度的升高，振動自由度逐漸被激發(fā)，比熱比呈現(xiàn)下降趨勢；當(dāng)溫度達(dá)到約2000 K時，海拔高度為60 km和1 km處大氣比熱比隨溫度變化曲線開始分離，這是由于密度不同時，氮分子和氧分子離解特征溫度和離解速率不同所導(dǎo)致的；當(dāng)溫度進(jìn)一步升高時，比熱比隨因化學(xué)反應(yīng)導(dǎo)致新組分的產(chǎn)生而不斷波動變化，但是，相比于熱容的波動，比熱比的波動規(guī)律更加復(fù)雜.

圖5 不同海拔高度下比熱比隨溫度的變化Fig.5 The specific heat ratios versus temperature at different altitudes.

圖7 60 km海拔高度下飛行器駐點處空氣組分隨飛行速度的變化.Fig.7 Air components at aircraft's stagnation point with flight speed at altitude of?=60 km.

圖8 不同海拔下飛行器駐點處空氣密度壓縮比隨飛行速度的變化.Fig.8 The density compression ratios of air with flying speed at different altitudes.

圖9 不同海拔下飛行器駐點處空氣的溫度隨飛行速度的變化.Fig.9 The air temperatures at the aircraft's stagnation point with flying speed at different altitudes.

圖10 不同海拔下飛行器駐點處空氣定壓熱容隨飛行速度的變化.Fig.10 The air's constant pressure heat capacities at the aircraft's stagnation point with flying speed at different altitudes.

圖11 不同海拔下飛行器駐點處空氣比熱比隨飛行速度的變化Fig.11 The air's specific heat ratios at the aircraft's stagnation point with flying speed at different altitudes.

需要說明的是，上述高溫空氣組分及熱物性參數(shù)計算主要針對等容過程，即其結(jié)果是在確定的密度條件下得到的.當(dāng)飛行器高速飛行時，其將向周圍空氣推出激波，飛行器周圍空氣在激波作用下形成高溫等離子體，這種高溫等離子體的壓力、密度和溫度狀態(tài)、相應(yīng)的熱容和比熱比等參數(shù)、以及所包含的空氣組分隨飛行器飛行速度和所處的海拔高度而變化.其中，在激波作用區(qū)域，飛行器駐點處空氣經(jīng)激波作用所產(chǎn)生的高溫等離子體熱物性參數(shù)、所發(fā)生的物理化學(xué)反應(yīng)以及空氣組分隨飛行速度的變化趨勢對于工程設(shè)計具有重要的實際意義.為此，基于所建立的模型，我們對飛行器在不同海拔下(海拔高度H=1.8，16，32，40，55，65 km激波前分別對應(yīng)密度ρ0=1.0，0.1，10-2，5×10-3，10-3，10-4kg/m3)飛行時其駐點處高溫空氣的組分、密度壓縮比(ρ/ρ0)、駐點溫度、熱容及比熱比等熱物性參數(shù)隨飛行速度的變化進(jìn)行了預(yù)測，其結(jié)果見圖6～11所示，這些預(yù)測結(jié)果不僅有助于認(rèn)識高速飛行器周圍高溫稀薄空氣因激波加熱產(chǎn)生高溫誘發(fā)的物理化學(xué)現(xiàn)象，也可以為飛行器熱防護(hù)材料等設(shè)計提供科學(xué)依據(jù).

圖6 1km海拔高度下飛行器駐點處空氣組分隨飛行速度的變化Fig.6 Air components at aircraft's stagnation point with flight speed at altitude of?=1km.

4 結(jié) 論

本文基于局域熱動平衡和局域化學(xué)平衡假設(shè)，針對N2-O2-Ar空氣體系，建立了考慮18種化學(xué)反應(yīng)產(chǎn)生22組元混合體系的熱化學(xué)物理計算模型，模型能夠依據(jù)溫度及密度的不同自適應(yīng)地實現(xiàn)從初始3組分(N2、O2、Ar)到最多22組分(N2、N2+、O2、O2-、O2+、NO、NO+、N2O、NO2、N、N+、N2+、N3+、O、O-、O+、O2+、O3+、Ar、Ar+、Ar2+、e)高溫稀薄空氣組分及熱物性參數(shù)的計算，其可靠性通過與國際流行多組分模型的對比得到了驗證.基于該模型，研究了高溫空氣化學(xué)反應(yīng)、組分、以及焓、熱容和比熱比等熱物性參數(shù)隨溫度和海拔高度的變化趨勢，分析了高溫空氣化學(xué)反應(yīng)發(fā)生的溫度條件(即反應(yīng)特征溫度)及其隨空氣稀薄程度不同的變化規(guī)律和高溫引起化學(xué)反應(yīng)導(dǎo)致的熱容和比熱比波動變化現(xiàn)象，并對飛行器高速飛行時駐點處的高溫空氣組分、壓縮比、溫度、熱容以及比熱比等參數(shù)隨飛行速度的變化給出了預(yù)測.所取得的結(jié)果及規(guī)律性認(rèn)識，有助于我們更好地認(rèn)識高速飛行器周圍稀薄空氣因強(qiáng)激波壓縮產(chǎn)生的高溫誘導(dǎo)發(fā)生的物理化學(xué)反應(yīng)、現(xiàn)象和作用機(jī)理，并為飛行器熱防護(hù)材料等設(shè)計提供科學(xué)依據(jù).