何新 江濤 張振福 楊俊波

1) (國(guó)防科技大學(xué)文理學(xué)院,長(zhǎng)沙 410073)

2) (中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所,綿陽(yáng) 621000)

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

1 引言

在載人航天、太空探測(cè)等應(yīng)用中,高超聲速飛行器返回時(shí)被大強(qiáng)度激波包裹,激波層中的空氣將變?yōu)榈入x子體[1].當(dāng)飛行器速度很高時(shí),其所承受熱負(fù)荷的很大一部分源于激波層的輻射加熱[2].為了評(píng)估輻射加熱、指導(dǎo)熱防護(hù)設(shè)計(jì),必須研究空氣等離子體的輻射特性[3,4].

分析原子激發(fā)是這其中的一個(gè)十分重要的方面,因?yàn)樵诟叱曀偌げ▽涌諝獾入x子體中,原子的輻射將占據(jù)主導(dǎo)地位[5,6].例如阿波羅飛行器返回時(shí),約90%的激波層輻射來(lái)自于原子發(fā)射譜線(xiàn)[7].

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

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

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

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

2 計(jì)算方法

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

根據(jù)束縛態(tài)特征溫度法[20],若設(shè)某種原子的電離能為I,則表征該種原子能級(jí)布居的特征參數(shù)Tb可由下式計(jì)算:

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

其中Ei和gi分別是原子第i能級(jí)的能量和簡(jiǎn)并度,分別是一價(jià)離子第j能級(jí)的能量和簡(jiǎn)并度.

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

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

3 典型空氣等離子體算例及結(jié)果

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

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

由(1)—(3)式可知,計(jì)算能級(jí)布居還需要原子及其離子的能級(jí)參數(shù).本文中,N 和O 原子能級(jí)數(shù)據(jù)來(lái)自文獻(xiàn)[22],N 和O 原子電離能采用美國(guó)國(guó)家標(biāo)準(zhǔn)與技術(shù)研究院(NIST)數(shù)據(jù)[23],N+和 O+離子能級(jí)數(shù)據(jù)來(lái)自文獻(xiàn)[1].

圖1 為狀態(tài)點(diǎn)1634-25 的氮原子非簡(jiǎn)并能級(jí)布居.標(biāo)有“Johnston”的數(shù)據(jù)由文獻(xiàn)[7]中方法計(jì)算而得;標(biāo)有“Boltzmann”的數(shù)據(jù)由Te下的Boltzmann 分布,即(3)式中的Tb換成Te計(jì)算而得(下同);標(biāo)有“Saha”的數(shù)據(jù)由另一形式的Saha 方程計(jì)算而得(下同):

眾所周知,若等離子體處于熱平衡態(tài),則Te下的Boltzmann 分布與Saha 方程的結(jié)果必然是相等的.圖1 中,標(biāo)有“Boltzmann”和“Saha”的數(shù)據(jù)只存在輕微偏離,說(shuō)明該狀態(tài)點(diǎn)空氣等離子體處于近平衡態(tài).根據(jù)圖1,對(duì)于N 的低能級(jí)和高能級(jí)布居,計(jì)算結(jié)果與CR 模型一致;對(duì)于N 的中間一些能級(jí)布居,計(jì)算結(jié)果介于CR 模型和Johnston 方法.

圖1 狀態(tài)點(diǎn)1634-25 的氮原子能級(jí)布居Fig.1.Energy level populations for N of Case 1634-25.

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

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

CR 模型是更為精細(xì)的物理模型,準(zhǔn)確度更高.圖1 和圖2 的結(jié)果表明,若以CR 模型為參照討論計(jì)算精度,則本文計(jì)算精度與Johnston 方法接近,優(yōu)于Spradian 模塊.

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

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

圖3 狀態(tài)點(diǎn)1634-7 的氮原子能級(jí)布居Fig.3.Energy level populations for N of Case 1634-7.

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

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

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

圖5 為狀態(tài)點(diǎn)1643-5 的氧原子非簡(jiǎn)并能級(jí)布居.此狀態(tài)點(diǎn)對(duì)應(yīng)Fire II 的飛行高度低,激波層中粒子能量松弛充分,因此空氣等離子體處于近平衡態(tài).由圖5 可知,計(jì)算結(jié)果與CR 模型一致.

圖5 狀態(tài)點(diǎn)1643-5 的氧原子能級(jí)布居Fig.5.Energy level populations for O of Case 1643-5.

上述5 個(gè)不同狀態(tài)點(diǎn)的結(jié)果表明,本文計(jì)算得到的空氣等離子體原子能級(jí)布居與CR 模型基本一致,是合理有效的.在計(jì)算精度上,本文與QSS Abba和Johnston 簡(jiǎn)化方法接近,優(yōu)于QSS Park 和Spradian 簡(jiǎn)化方法.

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

4 結(jié)論

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