高溫空氣熱化學物理計算模型研究

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

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

1 引 言

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

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

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

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

2 理論模型

2.1 平衡態組分

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

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

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

?

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

以O2和O為例，其內配分函數的表達式為：

O原子內部電子激發對原子內配分函數的貢獻可以表示為：

O2分子的轉動配分函數可以表示為：

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

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

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

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

O2分子電子部分對內配分函數的貢獻可以表示為：

式中，ε2i為O2分子內部電子能級，g2i為對應的簡并度.

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

2.2 空氣熱物性模型

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

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

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

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

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

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

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

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

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

在飛行器高速飛行推出強激波作用下，其駐點處的高溫空氣熱力學狀態滿足Hugoniot能量約束關系：

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

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

3 數值結果與分析

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

圖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

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

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

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

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

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

圖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.

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

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

4 結 論

本文基于局域熱動平衡和局域化學平衡假設，針對N2-O2-Ar空氣體系，建立了考慮18種化學反應產生22組元混合體系的熱化學物理計算模型，模型能夠依據溫度及密度的不同自適應地實現從初始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)高溫稀薄空氣組分及熱物性參數的計算，其可靠性通過與國際流行多組分模型的對比得到了驗證.基于該模型，研究了高溫空氣化學反應、組分、以及焓、熱容和比熱比等熱物性參數隨溫度和海拔高度的變化趨勢，分析了高溫空氣化學反應發生的溫度條件(即反應特征溫度)及其隨空氣稀薄程度不同的變化規律和高溫引起化學反應導致的熱容和比熱比波動變化現象，并對飛行器高速飛行時駐點處的高溫空氣組分、壓縮比、溫度、熱容以及比熱比等參數隨飛行速度的變化給出了預測.所取得的結果及規律性認識，有助于我們更好地認識高速飛行器周圍稀薄空氣因強激波壓縮產生的高溫誘導發生的物理化學反應、現象和作用機理，并為飛行器熱防護材料等設計提供科學依據.