束縛態特征溫度方法及應用*

何新 江濤 張振福 楊俊波

1) (國防科技大學文理學院,長沙 410073)

2) (中國空氣動力研究與發展中心計算空氣動力研究所,綿陽 621000)

隨著高超聲速飛行器速度增大,激波層空氣等離子體中的原子發射譜線成為輻射加熱主要來源,因此研究原子激發非常重要.考慮到處于熱非平衡態的空氣等離子體,平衡態統計理論不適用.精細物理模型(如碰撞輻射模型)雖然可以處理熱非平衡問題且準確度高,但計算量太大,難于工程應用.本文采用束縛態特征溫度法,結合FIRE II 激波管實驗中的非平衡空氣等離子體,對原子激發進行了分析.計算得到的原子能級布居與碰撞輻射模型符合,說明簡化計算是合理的,計算效率提高了2000 倍以上,且能夠保證一定的精度.

1 引言

在載人航天、太空探測等應用中,高超聲速飛行器返回時被大強度激波包裹,激波層中的空氣將變為等離子體[1].當飛行器速度很高時,其所承受熱負荷的很大一部分源于激波層的輻射加熱[2].為了評估輻射加熱、指導熱防護設計,必須研究空氣等離子體的輻射特性[3,4].

分析原子激發是這其中的一個十分重要的方面,因為在高超聲速激波層空氣等離子體中,原子的輻射將占據主導地位[5,6].例如阿波羅飛行器返回時,約90%的激波層輻射來自于原子發射譜線[7].

在高超聲速激波層中,部分空氣等離子體處于熱平衡態,部分處于熱非平衡態[8].對于熱平衡空氣等離子體可以方便地利用Boltzmann 分布或Saha方程得到原子能級布居[9].然而,對于熱非平衡空氣等離子體,計算其中的原子能級布居是一項挑戰.

特別是對于高超聲速航天器大尺度三維激波層計算,其中常常包含超大量的熱力學狀態點(有限元),在這種情況下,已發展的、較常用的碰撞輻射(collisional-radiative,CR)模型[3,4,8?11],雖然能夠處理熱非平衡問題,且準確度高,但計算耗費超大,甚至無法實現[12].

為了既能處理熱非平衡問題,又在誤差可接受的前提下降低計算成本,研究者開發了一些簡化計算方法.例如,采用準穩態(quasi-steady-state,QSS)近似、多溫度Boltzmann 分布等來計算得到空氣等離子體的原子激發數據[13?17],以便與大型飛行器流場計算耦合應用.然而,目前這些簡化方法普遍基于精細CR 模型的思想,需要用到許多微觀粒子相互作用速率系數,這些系數難以保證準確,從而導致這些簡化方法之間、它們與CR 模型之間存在較明顯的偏差[18,19],直接影響到飛行器輻射加熱評估和熱防護策略選擇.可以說,針對高超聲速大尺度三維等離子體計算需求,能夠保證一定精度且快速計算分析原子能級布居,是研究者非常關注的課題和不斷追求的目標.

在前期工作中,所提出的束縛態特征溫度法無需用到微觀粒子相互作用系數,在激光等離子體相互作用類似情形的計算中得到驗證.本文在該方法基礎上,對高超聲速空氣等離子中的原子激發進行研究,嘗試發展適用于高超聲速激波層等離子體的解決途徑.選取高超聲速實驗中典型的熱非平衡和熱平衡空氣等離子體作為研究對象,計算其中氮和氧原子能級布居,并與CR 模型、其他簡化模型的結果進行對比,分析計算準確度和計算效率.

2 計算方法

本文關注廣泛存在于載人航天器返回、太空探測器再入等航天應用領域的空氣等離子體狀態,其自由電子溫度一般不高(30000 K 以下),且為弱電離.通常在這樣的空氣等離子體中,只需考慮原子(N,O)及其一價離子(N+,O+)的存在[15].

根據束縛態特征溫度法[20],若設某種原子的電離能為I,則表征該種原子能級布居的特征參數Tb可由下式計算:

其中Tb為該種原子的束縛態特征溫度,Te為自由電子溫度,n,n+和ne分別為原子、一價離子和自由電子的數密度,me為自由電子質量,k和h分別為Boltzmann 常數和Planck 常數,Q(Tb)和Q+(Te)分別為原子和一價離子的配分函數:

其中Ei和gi分別是原子第i能級的能量和簡并度,分別是一價離子第j能級的能量和簡并度.

若已知n,n+,ne和Te,可根據(1)式計算出Tb,從而原子第i能級的非簡并布居為：

實際上,(1)式可稱為修正的Saha 方程.這樣寫的好處是對熱平衡和熱非平衡空氣等離子體都適用.對于熱平衡空氣等離子體,必然存在Tb=Te,則(1)式就是眾所周知的Saha 方程[1];對于熱非平衡空氣等離子體,所求解出的Tb將與Te不同,二者之間的差別反映了能級布居偏離熱平衡分布的程度.

3 典型空氣等離子體算例及結果

選取針對Fire II 工程的地面激波管實驗空氣等離子體為研究對象[21],其中的輻射主要來源于原子譜線,對原子能級布居(尤其是較高能級布居)計算非常重要.表1 給出了所選取的空氣等離子體狀態參數(包括自由電子溫度、粒子數密度),它們分別對應Fire II 飛行過程中的1634,1636 和1643 s 時間點.

表1 空氣等離子體參數Table 1.Parameters of air plasmas.

由(1)—(3)式可知,計算能級布居還需要原子及其離子的能級參數.本文中,N 和O 原子能級數據來自文獻[22],N 和O 原子電離能采用美國國家標準與技術研究院(NIST)數據[23],N+和 O+離子能級數據來自文獻[1].

圖1 為狀態點1634-25 的氮原子非簡并能級布居.標有“Johnston”的數據由文獻[7]中方法計算而得;標有“Boltzmann”的數據由Te下的Boltzmann 分布,即(3)式中的Tb換成Te計算而得(下同);標有“Saha”的數據由另一形式的Saha 方程計算而得(下同):

眾所周知,若等離子體處于熱平衡態,則Te下的Boltzmann 分布與Saha 方程的結果必然是相等的.圖1 中,標有“Boltzmann”和“Saha”的數據只存在輕微偏離,說明該狀態點空氣等離子體處于近平衡態.根據圖1,對于N 的低能級和高能級布居,計算結果與CR 模型一致;對于N 的中間一些能級布居,計算結果介于CR 模型和Johnston 方法.

圖1 狀態點1634-25 的氮原子能級布居Fig.1.Energy level populations for N of Case 1634-25.

圖2 給出了狀態點1634-10 的氮原子和氧原子非簡并能級布居.“Boltzmann”和“Saha”數據之間存在明顯偏離,說明此狀態點空氣等離子體為非平衡.這是由于距離激波面較近,粒子之間能量松弛不充分造成的.

如圖2(a)所示,除N 原子的第2,3 能級外,本文計算得到的其他能級布居與CR 模型一致;而此條件下采用Spradian 模塊得到的結果明顯大于CR 模型.如圖2(b)所示,除O 原子的第2—4能級外,本文計算得到的其他能級布居與CR 模型符合;而Spradian 模塊的結果明顯大于CR模型.

CR 模型是更為精細的物理模型,準確度更高.圖1 和圖2 的結果表明,若以CR 模型為參照討論計算精度,則本文計算精度與Johnston 方法接近,優于Spradian 模塊.

圖3 為狀態點1634-7 的氮原子非簡并能級布居.顯然,此狀態點為非平衡空氣等離子體.如圖3(a),計算得到的第2,3 能級占據數略小于CR 模型,但其他能級布居與CR 模型一致,這與圖2 中的現象類似.圖3(b)中還與其他簡化模型的結果進行了對比.可以看出,計算結果與Johnston 方法符合,在一些中間能級與QSS Abba 方法符合,但整體上明顯低于QSS Park 方法.如果仍然參照CR 模型討論計算精度,圖3 的結果表明,在此狀態條件下,本文計算精度與Johnston 和QSS Abba 方法接近,優于QSS Park 方法.

圖2 狀態點1634-10 的氮原子和氧原子能級布居 (a) N;(b) OFig.2.Energy level populations for N and O of Case 1634-10:(a) N;(b) O.

圖3 狀態點1634-7 的氮原子能級布居Fig.3.Energy level populations for N of Case 1634-7.

圖4 為狀態點1636-5 的氮原子和氧原子非簡并能級布居.此狀態點空氣等離子體雖然也是非平衡的,但由于粒子數密度較大,粒子之間能量松弛較圖3 所示情形充分,因此“Boltzmann”和“Saha”數據之間的偏離減小.對于N 的中間一些能級,計算得到的占據數略大于CR 模型.但總體來說,在此狀態點,計算結果與CR 模型符合.

圖4 狀態點1636-5 的氮原子和氧原子能級布居 (a) N;(b) OFig.4.Energy level populations for N and O of Case 1636-5:(a) N;(b) O.

圖2—圖4 中,本文所計算N 原子第2—3 能級、O 原子第2—4 能級的占據數比CR 模型偏小,對輻射加熱評估的影響較小.這是因為N 和O 原子的發射譜線并不包含上述能級的自發輻射躍遷[21].另外,對于N 和O 原子幾個中間能級的占據數,本文雖然與CR 模型存在一定偏差,但作為一種簡化計算結果,較其他簡化方法在計算精度上持平或有明顯提高.

圖5 為狀態點1643-5 的氧原子非簡并能級布居.此狀態點對應Fire II 的飛行高度低,激波層中粒子能量松弛充分,因此空氣等離子體處于近平衡態.由圖5 可知,計算結果與CR 模型一致.

圖5 狀態點1643-5 的氧原子能級布居Fig.5.Energy level populations for O of Case 1643-5.

上述5 個不同狀態點的結果表明,本文計算得到的空氣等離子體原子能級布居與CR 模型基本一致,是合理有效的.在計算精度上,本文與QSS Abba和Johnston 簡化方法接近,優于QSS Park 和Spradian 簡化方法.

在工程應用中,計算效率是必須考慮的問題.借助普通筆記本電腦(CPU:2×2.60 GHz,Matlab程序),計算得到表1 所有狀態點的能級布居數據約需7 s;而CR 模型借助IBM 服務器(CPU:6×2.53 GHz,Fortran 程序)約需4.5 h.即使忽略計算平臺性能、程序語言效率、算法流程設計方面的差別,計算速度也比CR 模型提高了2000 倍以上.本文計算效率的提升主要源于不需要求解能級布居方程組.CR 模型盡可能多地考慮了影響粒子能級布居的大量微觀過程,需要求解大規模的能級布居速率方程組(方程數量取決于所考慮能級數);而文中其他簡化方法雖然減少了所考慮的微觀過程,甚至進一步采用QSS 近似,但仍然需要求解能級布居方程組.因此,本文計算可與高超聲速飛行器流場計算耦合,極大降低計算成本.

4 結論

采用束縛態特征溫度法研究了高超聲速激波層空氣等離子體中的原子激發.以針對Fire II 工程的激波管實驗空氣等離子體為算例,對N 和O 原子能級布居進行了計算,并與CR 模型、其他簡化模型的結果進行了對比.所研究的空氣等離子體熱力學狀態包括近平衡態、非平衡態.結果表明,本文計算得到的原子能級布居與CR 模型基本一致,計算精度與其他簡化方法接近甚至有一定提高.在計算效率上,本文比CR 模型提高了2000 倍以上,在工程應用(如高速飛行器輻射加熱評估)中可大大節約計算成本.