極高超聲速稀薄氣體原子輻射效應的p-DSMC 方法

楊全順 方 明 ,1) 李埌全 粟斯堯 楊彥廣

* (中國空氣動力研究與發展中心超高速空氣動力研究所,四川綿陽 621000)

? (中國空氣動力研究與發展中心,四川綿陽 621000)

引言

氣動加熱是高超聲速飛行/再入不可回避的問題之一.飛行器在第一宇宙速度以上的極高超聲速再入/飛行過程中,頭部激波層氣體會變成對飛行器加熱的強輻射體,輻射的光子會顯著增強氣動加熱效應,輻射加熱隨著速度的增加逐漸變得重要.研究表明,飛行速度超過10 km/s 時[1-5],輻射加熱占對流加熱的30%左右,此時,飛行器的熱防護設計不僅需要考慮對流加熱、催化加熱等非輻射加熱[6-8],還需要考慮輻射加熱[9-12].

為準確計算非平衡流動中的輻射,NASA 很早開發了非平衡空氣輻射(NEQAIR)程序[13].在NEQAIR 程序中,準穩態假設用于確定非平衡氣體組分的激發態布居,被認為是非平衡流動的最佳逐線法(LBL)解算器.然而,NEQAIR-LBL 程序包含一個一維正切平板輻射輸運解算器,其成本非常昂貴,不適合二維或三維的輻射計算,此外,高分辨率的逐線方法計算需要數十萬個輻射輸運方程(RTE)解,也是計算二維或三維輻射成本太高的另一個因素.在這種情況下,光子蒙特卡羅方法(PMC)[14]應運而生,在此方法中,輻射傳輸通過允許每個網格以光子束的形式向隨機方向發射其能量來建模,對非灰輻射場更精確、更有效.之后,Ozawa 等[15]使用直接蒙特卡羅方法(DSMC)和PMC 松耦合模擬了原子在超高速再入流動中的輻射.我國學者在高溫氣體熱化學非平衡流動,開展了大量關于非平衡流動、傳熱問題的理論研究[16-19],但對非平衡輻射特性理論建模研究很少.

在實驗測試方面,Brandis 等[20]在地球大氣環境中,對激波速度8.0～11.5 km/s 的真空紫外至近紅外的空間、光譜分辨的輻射強度進行測量,與數值結果對比發現,在不同的激波速度下,LAURA(蘭利氣動熱力學迎風松弛算法)/HARA(高溫氣動熱動力學輻射)輻射熱流計算結果都有不同的高估和低估,試驗數據與理論模擬的偏差表明計算模型仍需完善.之后文獻[21]又以激波速度4.7～8.0 km/s 的混合氣體及純N2進行了輻射強度的測量,發現新的實驗數據不僅高于以往的實驗值,而且個別的新實驗值甚至高出了一個量級,這表明,假設以現有的實驗數據進行數值模擬,結果也是不準確的.2017 年,Cruden 和Brandis[22]測量了光譜范圍190～1450 nm的空氣非平衡輻射后發現,不同反應速率的DPLR/NEQAIR 對實驗數據的復現結果較差.同年,Brandis等[23]研究發現,激波速度8.0～11.5 km/s 的大氣再入環境下,NEQAIR 和HARA 的預測結果低于電弧激波管(EAST)測試的結果.綜上,近年來的輻射實驗和理論模擬研究表明,目前的流動和輻射計算模型不確定性仍然較大,并且部分實驗數據有待進一步完善.

基于以上的分析考慮,本文基于優化的原子輻射模型,提出p-DSMC 方法,研究了有無激發輻射效應影響下的壁面壓力和熱流以及沿駐點線變化的平動、振動和轉動溫度.文章的組織如下:第二部分簡要介紹計算模型和激發、輻射模式,第三部分是模型的實現與驗證,第四部分是對全文的總結.

1 輻射效應的DSMC 建模

1.1 DSMC 方法

目前,Bird[1]發展的DSMC 方法是求解稀薄流域高超聲速飛行器非平衡氣體繞流最可靠的計算工具.對于稀薄流域,由于連續介質模型失效,N-S 方程的有效性存疑,而求解玻爾茲曼方程對于巨大粒子量的真實氣體又過于復雜,尤其是對于碰撞項的積分求解.DSMC 方法用一個仿真分子代表一定量的真實氣體分子,并將分子運動與碰撞解耦,借助現象學模型和統計方法,實現流場參數和物體壁面屬性的直接計算.DSMC 的詳細陳述可見參考文獻[24],其有兩個基本要求:一是網格尺寸須小于分子平均自由程,二是時間步長須小于分子平均碰撞時間.

1.2 稀薄氣體原子輻射模型

1.2.1 原子能級模型

原子輻射模型最理想的情況是在原子輻射過程中涵蓋所有的原子能級.然而,輻射態的數量如此之多,以致于反應的數量會是一個大的未知因素.另外一個問題是大量輻射可能來自原子布居稀少的能級,而對于實際數量的模擬分子來說,這些躍遷很少發生,以至于輻射的統計散射會大得令人無法接受[1].早期的文獻[25]提到,上述問題可以通過現象學模型解決.

原子輻射呈現出一個更困難的問題是電子態的數量和輻射躍遷的數量大于分子的相應數量.計算上是可行的,但就計算機內存需求而言,這是一個重大的災難.Park[13]和Bird[1]分別提出了原子能級群組的思想,原子能級使用數字來代替光譜項表示,這樣有利于簡化能級.還有一個需要注意的問題是,能級的數量和計算時間是密切相關的,因為我們在輻射部分提到,亞時間步長必須要小于選取能級中最小的能級輻射壽命.Taylor 等[4]也在文獻中明確提出,在原有DSMC 代碼中加入激發輻射模塊后,計算時間要更久.但是,這不是最主要的問題,在激發輻射中最關心的問題應該是選取的激發輻射路徑對輻射加熱的影響,這是因為,原子的電子能級眾多,可選則的能級群組方式也多,但能級之間的躍遷必須要滿足選擇定則,此外,能級之間的躍遷即使不是禁戒的,也有強弱之分.Bird[1]和Park[13]在文獻中也對N,O 原子能級的選擇給出了建議,因此,在綜合考慮后,我們選擇了7 能級群組,這也與飛行試驗觀察到的典型譜線特征吻合.本文使用“7”能級結構如表1 所示.

表1 N,O 原子群組能級的能量和簡并度Table 1 Energy and degeneracy of energy levels of N and O atomic groups

1.2.2 原子的激發

(1)原子-原子碰撞激發的物理模型

原子-原子碰撞激發的過程以原子分子反應靜力學[26]為基礎.Gallis 和Harvey[27]沒有明確提出這種說法,但他們提出的原子-原子碰撞時形成偽分子的思想卻與原子分子反應靜力學理論不謀而合.比如N 原子,其基態為4S,運用原子分子反應靜力學原理,形成N2分子基態的過程為

假如其中的一個為激發態或兩個都是激發態,則形成N2分子的低/高激發態,例如

式中的數字代表自旋多重度,S,P 表示原子的電子態,Σ,Π,Δ 表示分子的電子態.

(2)原子激發路徑

原子能級運用群組的概念將導致能級不再使用光譜項表示能級,致使能級之間的躍遷不再嚴格地滿足各種選擇定則.各能級之間還是存在允許躍遷和禁戒躍遷,但目前的7 能級群組不存在這種情況,即

因此,碰撞對在滿足:①激發條件、②激發概率、③原子布居概率的情況下,原子可以從低能級向上躍遷到其他各個能級群組.

(3)原子的激發條件和概率

對于DSMC 方法來說,任何一個事件都是以一定概率發生的,原子與原子的碰撞激發也不例外.

對于選中的能夠碰撞的原子對,應該滿足下面的激發條件

其中,Etrans是碰撞對的相對平動能,Eu-El是躍遷能級的能級差.

由于在原子的激發過程中,標記每個原子作為一個單獨的態,所以,激發概率可以表示為[1]

式中,σex為原子-原子碰撞的激發截面,σcoll是彈性碰撞交叉截面.鑒于激發截面數據的不完整性和不可靠性,Bird[1]使用定性的手段,即電子與原子的彈性交叉截面(σcoll)量級是10-15cm2,而電子碰撞的激發交叉截面(σex)量級為10-16cm2,電子-原子碰撞的激發碰撞截面和總碰撞截面的比值為0.1,并使用概率常數作為激發概率.本文借用Bird 的思想,也選用概率常數,選取原子-原子碰撞截面比值為200-1.

(4)原子的能級布居

滿足原子激發條件,并且激發概率大于隨機數的原子激發到哪個能級,也是必須考慮的問題.分子或原子的第j電子態能級能量高于低能態的平衡分數由玻爾茲曼分布給出

式中,gj是高能級的簡并度,N為低于Ej的原子布居數,T為有效溫度并滿足下式

處于Ei能級上的原子被激發,假如隨機數小于該分數,則原子將被激發到Ej能級.

1.2.3 輻射模型

對于輻射,本文只涉及到了束縛-束縛態之間的自發輻射.在10 km/s 左右的來流條件下,雖然能量已經很高了,駐點區的溫度達到了2.0×104K,但相較于等離子體環境,能量還是偏低,因此,束縛-束縛電子態之間的躍遷產生的輻射是最明顯的.即使在當前的情況下有其他的兩種躍遷,但其躍遷強度較之束縛-束縛的躍遷強度都較弱,對其模型壁面的輻射加熱可忽略.之前文獻[5]也和本文的處理方式一致;當然,本文是我們前期計算輻射熱流的一種方式,對于更高再入速度下的輻射特征,其他兩種躍遷可能是重要的,這將作為我們后續的工作之一.

處于激發態上的原子擁有壽命,當模擬的時間步長大于其自發輻射壽命時,原子將以一定的概率自發輻射回到低能級態.

(1)輻射路徑

本文采用N,O 原子7 能級群組進行自發輻射躍遷.涉及到的躍遷能級已羅列在表2 中,路徑呈現于圖1.需要注意的是,Bird[1]已經證實,只有表2中涉及的能級躍遷才擁有較強的光譜,但本文在此添加了N 原子低激發態2 和3;O 原子低激發態4,3 和2,希望可以獲得更精確的輻射加熱數據,相應的信息也一并列于表2.

圖1 N(上)、O(下)原子群組能級躍遷Fig.1 Energy level transition of N (upper)and O (lower)atomic groups

表2 N,O 原子群組激發態的輻射壽命和躍遷輻射幾率Table 2 Radiation lifetime and transition radiation probability of excited states of N and O atomic groups

(2)輻射條件和概率

在原子的激發模型中已經提到,由于此時的能級不是使用光譜項而是數字表示,所以,能級之間的躍遷不再嚴格遵守選擇定則.但本文所選取能級的標準是能級之間不存在禁戒躍遷,處于高能級的原子可以向任意低能級進行自發輻射躍遷.需要注意的是,DSMC 選擇的時間步長遠大于能級的自發輻射壽命,假如選擇粒子的運動時間步長,則會導致大量的光譜線消失.因此,在DSMC 中添加輻射模塊時必須使得光子的運動時間步長小于7 能級中最小的輻射壽命.光子的時間步長大于能級的輻射壽命時,釋放出光子,即輻射條件為

正是由于每個能級表示的是能級群組,低能態躍遷是以一定的份數躍遷,本文將此份數作為自發輻射躍遷的概率(P).涉及到的躍遷路徑,哪個那級的自發輻射路徑的概率大于隨機數,原子將自發輻射回到該低能級.

1.2.4 光子追蹤

在目前的計算工況下,光子速度大于流場中粒子速度四個數量級,假如光子以粒子運動的時間尺度向前運動,追蹤光子是及其困難的.因此,本文采用亞時間步長的思想,在一個亞時間步長內,將光子的運動距離控制在一個網格寬度內.這樣做的優勢在于,不僅可以追蹤到光子的軌跡,而且可以減少內存分配.

光子追蹤的另一問題是跨進程并行,如圖2 所示.采用亞時間步長解決了光子追蹤問題,但光子運動的最終時間是一個時間步長.在一個時間步長內,光子肯定跨越好幾個進程進行運動.本文中已經使用MPI 的光子跨進程并行程序,最大核數測試達到500 個核,并行效果良好.

圖2 光子追蹤示意圖Fig.2 Photon tracing diagram

2 p-DSMC 方法驗證及輻射效應影響驗證

本文使用列于表3 中的初始流場條件(包括數密度N0、溫度T0、馬赫數Ma、壁溫Twall和密度ρ),相當于大氣高度86 km,壁面為全漫反射,涉及19 化學反應,運用p-DSMC 方法得到了二維圓柱繞流的流場分布.在p-DSMC 模擬中,流場涉及5 組元,每個網格中初始分布20 個模擬粒子,按照沈青[28]建議,網格寬度取大約1/3 的平均自由程,時間步長為1×10-7s.

表3 流場的初始參數Table 3 Initial parameters of flow field

2.1 不含輻射效應的DSMC 代碼驗證

使用已建立的DSMC 代碼MONACO[29]和dsmcFoma[30]與目前的p-DSMC 代碼進行比較,MONACO 和p-DSMC 代碼均采用TCE 化學模型[31],而dsmcFoma 使用D-K 化學模型[32].在所有的例子中,轉動和振動碰撞數分別設置為5 和50.在考慮化學反應且速度為6.7 km/s 時,基于表3 中的來流初始參數,在不包含N,O 原子的激發和輻射時,獲得了二維圓柱壁面的壓力和熱流,并與之前的模擬結果[29]進行了比較,二者符合的很好.壁面的壓力和熱流誤差均在5%以內,尤其是在駐點位置,誤差在1%以內;另外,獲得的平動、振動以及轉動溫度均與文獻結果符合較好,這些信息也可以從圖3中看出.此外,從圖3 中可以看出,由于MONACO與本文使用相同的化學反應模型(TCE),因此,壁面壓強和熱流都更接近MONACO 獲得的值.

圖3 馬赫數24.58 時,(a)沿壁面的熱流和壓力以及(b)沿駐點線變化的平動、轉動和振動溫度Fig.3 (a)Heating flux and pressure on wall and (b)translational,rotational and vibrational temperatures along the stagnation line at Ma=24.58

2.2 輻射對流場性質的影響

在流場條件不變的情況下,圖4 和圖5 分別展示了考慮原子輻射效應后對流場全局溫度、平動溫度和內部模態溫度的影響.從圖4 中可以看出,隨著馬赫數的增大,輻射對流場的影響越大,馬赫數越大時,熱流最大區越靠近壁面.這種效應與圖6 所示的向右速度偏移所表明的差異一致.對于非平衡氣體,全局溫度定義為平動溫度和內部模態溫度的加權平均,因此,圖5 呈現了平動、轉動和振動溫度在馬赫數增大時輻射對其產生的影響.從圖5 可以看出,輻射對平動、轉動和振動溫度的最大值幾乎沒有影響,但隨著馬赫數的增大,平動溫度和轉動溫度都有靠近駐點的趨勢,這與整個流場展示的全局溫度相對應.圖5 呈現的現象與文獻2 給出的平動、振動和轉動溫度的變化規律是一致的,即在考慮了輻射效應后,平動和內部模態溫度均沒有發生很大改變.

圖4 原子輻射效應對流場的影響Fig.4 Influence of atomic radiation effect on the flow field

圖5 原子輻射效應對(a)平動、(b)轉動和(c)振動溫度的影響Fig.5 Influence of atomic radiation effect on the (a)translational,(b)rotational and (c)vibrational temperature

圖6 原子輻射效應對速度的影響Fig.6 Influence of atomic radiation effect on the velocity

2.3 輻射效應對壁面熱流的影響

圖7 是考慮了N,O 原子激發和輻射后壁面的輻射加熱和對流加熱沿流場中心線的變化規律.從圖7 中可以看出,輻射熱和對流熱有相似的變化趨勢,即隨著遠離駐點區,輻射加熱變得越來越小,并且輻射熱流減小的速率要大于對流加熱.這是因為在背風區,溫度降低,粒子能量減小,碰撞導致分子的解離減少,致使原子濃度變少,圓柱壁面輻射加熱相應減小.這說明流場內原子濃度是影響壁面輻射熱流的一個重要因素.圖7 還展示出,速度為6.74 km/s 時,輻射加熱不明顯,大約只有5 kW/m2,只占到對流加熱的7%左右,考慮輻射效應后,對流加熱幾乎沒有變化.隨著馬赫數的增大,激發的原子數目增多,壁面的輻射加熱和對流加熱的變化變得明顯,速度為10.14 km/s 時,輻射加熱占對流加熱的30%左右,這與Bird[1]的結果一致,即飛行速度超過10 km/s 時,輻射加熱變得明顯.此外,選中的原子碰撞對碰撞后,部分相對平動能轉化為原子的激發能,致使相對平動能減小,原子速度降低,當來流速度越大時,這種現象就越明顯.因此,隨著馬赫數的增大,考慮輻射效應后,對對流加熱的影響就越明顯.可以從圖7 中看出,當來流速度增加到10 km/s時,對流加熱改變量為40 kW/m2左右,改變量是來流速度為6.47 km/s 改變量的(2.5 kW/m2左右)16 倍.

圖7 原子輻射效應對壁面的輻射加熱和對流加熱Fig.7 Influence of atomic radiation effect on the radiative heating and convective heating on the wall

3 總結

基于以上的激發輻射模型和光子追蹤技術,運用p-DSMC 方法,獲得了超高聲速二維圓柱繞流原子氣體對壁面的輻射加熱,從得到的數據中可以得出以下三條結論.

(1)來流速度低于10 km/s 時,輻射加熱不明顯,在駐點區域,輻射加熱占對流加熱比重在7%左右;來流速度大于10 km/s 時,在駐點區域,輻射加熱占對流加熱比重將超過30%.

(2)考慮了原子的輻射效應,對非平衡區的平動、轉動和振動溫度的最大值影響不大.

(3)在考慮了原子的輻射效應后,對于二維圓柱繞流,激波要更加貼體.

(4)流場中原子的濃度是影響壁面輻射熱流大小的一個重要因素之一.

在本研究的基礎上,將會持續開展高溫空氣組分的輻射效應建模研究,開發計算代碼,并加強對比驗證,為目標飛行器熱防護系統提供可靠的熱載荷數據.