含電離化學(xué)反應(yīng)DSMC模擬的稀有組分權(quán)重因子方法

方 明, 杜波強(qiáng), 李中華, 李丹楊

(1. 中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所, 四川 綿陽 621000； 2. 國家計(jì)算流體力學(xué)實(shí)驗(yàn)室, 北京 100191; 3. 北京大學(xué) 物理學(xué)院, 北京 1000871)

0 引 言

稀薄氣體電離現(xiàn)象是航天器再入速度持續(xù)增大面臨的新問題，直接導(dǎo)致通信黑障等傳統(tǒng)難題向稀薄區(qū)域大幅延伸。傳統(tǒng)的再入物理[1]通過求解含化學(xué)反應(yīng)的N-S方程組[2]得到再入體繞流電子數(shù)密度分布，但是該方程組對于稀薄氣體流動失效。近年來發(fā)展的稀薄氣體數(shù)值模擬方法中，基于Boltzmann方程的分析和計(jì)算[3-4]尚無法包含電離過程，只有基于粒子隨機(jī)統(tǒng)計(jì)模擬的DSMC(Direct Simulation Monte Carlo)方法[5-9]是最有可能解決稀薄氣體電離數(shù)值預(yù)測難題的手段。

第一宇宙速度下的稀薄氣體電離度極低，通常在10-5量級，相對于N2、O2等來流氣體組分以及反應(yīng)生成的NO、N和O，電子和離子屬于稀有組分。使用傳統(tǒng)DSMC方法模擬計(jì)算網(wǎng)格中的一個電子，將需要105個其它粒子，計(jì)算代價過高而不可接受?；谶@一原因，早期稀薄氣體電離過程的DSMC模擬[10-11]主要針對10 km/s以上的低維再入流動，此時電離度在1%量級，一維駐點(diǎn)線和軸對稱計(jì)算的稀有組分統(tǒng)計(jì)漲落可以通過增大仿真分子數(shù)目控制在可以接受的范圍內(nèi)。文獻(xiàn)[12]建立了三維復(fù)雜外形航天器再入稀薄電離特性分析的DSMC方法和并行計(jì)算平臺，但是文中結(jié)果表明，RAM-C II的電子數(shù)密度分布存在極為顯著的統(tǒng)計(jì)漲落。

電子等稀有組分粒子數(shù)密度的統(tǒng)計(jì)波動帶來若干不利影響。首先，真實(shí)的電子數(shù)密度分布應(yīng)該較為光滑，存在較大波動的電子數(shù)密度分布在物理上是不真實(shí)的，這種統(tǒng)計(jì)波動是數(shù)值原因造成的，應(yīng)該可以用數(shù)值的方式予以抑制。第二，電子和離子在網(wǎng)格中統(tǒng)計(jì)數(shù)目的波動，給網(wǎng)格中分子碰撞對數(shù)目的確定帶來了一定程度的不確定性，而分子碰撞對數(shù)目的正確計(jì)算是含電離化學(xué)反應(yīng)過程正確模擬的基礎(chǔ)，亦即稀有組分粒子數(shù)密度的統(tǒng)計(jì)波動可能導(dǎo)致物理上不真實(shí)的化學(xué)反應(yīng)計(jì)算。第三，從應(yīng)用角度說，比如電子數(shù)密度對通信設(shè)計(jì)的預(yù)測分析，如果電子數(shù)密度分布本身波動較大，在應(yīng)用上的指導(dǎo)意義也大打折扣。

現(xiàn)有稀薄氣體DSMC模擬的稀有組分處理方法可以籠統(tǒng)地分為三種類型：間接方法、權(quán)重格式和權(quán)重格式的改進(jìn)格式。間接方法不直接模擬稀有組分粒子，如Bird的工作[13]并不真正計(jì)算NO+和電子，只是在N和O碰撞時標(biāo)識其發(fā)生電離的比例，根據(jù)這一比例和N、O數(shù)密度計(jì)算電子數(shù)密度分布，然而這一處理方式僅考慮N和O的聯(lián)合電離，導(dǎo)致其預(yù)測的電子數(shù)密度低估了2到10倍[14]，而且這一方法無法處理離子的電荷交換過程；Shevyrin的工作[15]不將電子作為一種粒子模擬，其數(shù)密度由網(wǎng)格中的離子之和得出，記錄電子能量用以計(jì)算電子溫度，對聯(lián)合電離逆反應(yīng)單獨(dú)處理，這一方法顯然未包括電子-原子碰撞導(dǎo)致的直接電離。

權(quán)重格式最早由Bird提出[16]，用每個仿真粒子代表的真實(shí)粒子數(shù)FNi作為粒子權(quán)重，權(quán)重較小的粒子對應(yīng)稀有組分。兩個不同權(quán)重的粒子發(fā)生碰撞時，將權(quán)重較大的粒子一分為二，其中一個與稀有組分具有相同權(quán)重；讓兩個具有相同權(quán)重的粒子二體碰撞，然后將碰后非稀有組分粒子與分裂得到的另一個粒子融合。Boyd指出[17]，上述格式不能保證每次碰撞動能守恒，進(jìn)而提出了守恒格式，即將上述每次碰撞的動能損失記錄下來，添加到下一個時間步長的常規(guī)組分碰前相對平動能中。隨后，Boyd等將其守恒權(quán)重格式發(fā)展用于模擬OH紫外輻射的變權(quán)重形式[18]。受當(dāng)時的研究局限，上述權(quán)重格式均未拓展到含電離的化學(xué)反應(yīng)過程。

權(quán)重格式的另一種形式是對所有組分使用相同的FN，但是對每個組分設(shè)置不同的權(quán)重Wi；常規(guī)組分的權(quán)重為1，稀有組分的權(quán)重小于1. 這一方法在稀薄氣體電離DSMC模擬中的應(yīng)用可以追溯到Bartel等人的工作[19]，但他們并未考察不同權(quán)重粒子之間的碰撞細(xì)節(jié)。Boyd[14]在Bird[16]的基礎(chǔ)上，人為增大電離化學(xué)反應(yīng)概率R倍，同時復(fù)制電離反應(yīng)產(chǎn)物S次，即賦予帶電粒子權(quán)重為1/(R×S)；這一權(quán)重用于確定正確的碰撞對分子數(shù)目，兩個不同權(quán)重的粒子碰撞過程仍采用Bird的分裂思想。

權(quán)重格式的改進(jìn)包括Boyd等的重疊方法[20]，而樊菁等的TSS算法[21]可以看成是權(quán)重格式的最新發(fā)展。重疊方法是一種兩步流場計(jì)算方法，第一步(基礎(chǔ)模擬)使用較大權(quán)重，不考慮稀有組分；第二步(重疊模擬)使用較小權(quán)重，常規(guī)組分的產(chǎn)生受到抑制。在重疊模擬中，常規(guī)-稀有組分碰撞中的常規(guī)組分粒子由根據(jù)網(wǎng)格中常規(guī)組分屬性抽樣產(chǎn)生的虛擬粒子代替，常規(guī)-常規(guī)組分碰撞生成稀有組分粒子數(shù)目由網(wǎng)格中的反應(yīng)物粒子數(shù)密度等參數(shù)計(jì)算得到。重疊方法主要用于輻射模擬，后來又發(fā)展處理表面過程[22]，但尚未見到處理稀薄氣體電離過程的報導(dǎo)。TSS方法與常規(guī)DSMC的主要區(qū)別是稀有組分的產(chǎn)生不依賴碰撞過程，而根據(jù)網(wǎng)格中反應(yīng)物的粒子數(shù)目計(jì)算，稀有組分的內(nèi)能狀態(tài)根據(jù)當(dāng)?shù)睾暧^參數(shù)抽樣得到。

本文工作基于稀有組分權(quán)重格式處理稀薄氣體含電離化學(xué)反應(yīng)過程。與前述方法不同的是，本文在前期工作[12]的基礎(chǔ)上，對不同權(quán)重粒子之間的碰撞依概率改變內(nèi)能狀態(tài)，對不同權(quán)重粒子之間的化學(xué)反應(yīng)依概率刪除或保留反應(yīng)物粒子，同時對生成物中的稀有組分粒子進(jìn)行復(fù)制。該方法在對成熟DSMC模擬框架和程序作極小化修正的前提下，可顯著提高電子等稀有組分?jǐn)?shù)密度等值線光滑性和計(jì)算精度。

1 稀薄氣體電離反應(yīng)的DSMC方法

1.1 化學(xué)反應(yīng)模型

再入過程的純大氣組分二體碰撞化學(xué)反應(yīng)，通?？梢詫懗扇缦滦问?/p>

A+B→C+D+E(1)

涉及電離的化學(xué)反應(yīng)包括五種：雙原子分子的離解反應(yīng)、置換反應(yīng)、聯(lián)合電離反應(yīng)及其逆反應(yīng)、直接碰撞電離反應(yīng)、離子交換反應(yīng)。對于離解反應(yīng)，E代表原子；對于直接電離反應(yīng)，E代表電子；而對于其它反應(yīng)，E為空。

化學(xué)反應(yīng)速率通常表示為如下Arrhenius形式的方程[23]

k(T)=ΛTηexp(-Ea/κT)(2)

其中Λ和η為常數(shù)，Ea為反應(yīng)活化能，κ為Boltzmann常數(shù)，T為溫度。當(dāng)碰撞總能量Ec>Ea時，上述常數(shù)Λ和η通過下式與化學(xué)反應(yīng)概率P相關(guān)聯(lián)[24]

1.2 弱等離子體效應(yīng)簡化處理

電子與離子及常規(guī)組分的巨大質(zhì)量差異是稀薄氣體電離過程DSMC模擬的本質(zhì)困難，電子真實(shí)質(zhì)量的使用和運(yùn)動軌跡的捕捉要求極短的時間步長，進(jìn)而導(dǎo)致不可接受的計(jì)算代價。受Boyd工作[14]啟發(fā)，本文采用增大電子質(zhì)量三個數(shù)量級的方式對上述等離子體效應(yīng)簡化處理。電子的真實(shí)質(zhì)量為9.11×10-31kg，增大三個數(shù)量級至9.11×10-28kg。以氧原子O為例，其質(zhì)量為2.656×10-26kg，調(diào)整后的電子質(zhì)量約為氧原子質(zhì)量的3.4%。在溫度相同的情況下，調(diào)整后電子的隨機(jī)熱運(yùn)動速度均值比氧原子大不到一個數(shù)量級。同時，為保證化學(xué)反應(yīng)過程質(zhì)量守恒，增大電子質(zhì)量須同步減小離子質(zhì)量。以氧離子O+為例，真實(shí)質(zhì)量為2.656×10-26kg ，調(diào)整后的質(zhì)量為2.5649×10-26kg ，質(zhì)量變化約為3.4%，亦不會對其熱運(yùn)動速度造成明顯影響。

表1 涉及電子的相關(guān)化學(xué)反應(yīng)常數(shù)Table 1 Chemical reaction rate coefficients involving electron

2 含電離化學(xué)反應(yīng)的稀有組分權(quán)重因子方法

2.1 權(quán)重因子對粒子碰撞的影響

本文采用權(quán)重格式的第二種形式，即對每個組分設(shè)置不同的權(quán)重Wi，N2、O2、NO、N、O等常規(guī)組分的權(quán)重因子設(shè)置為1，離子和電子等稀有組分的權(quán)重因子小于1。

如Boyd[16]指出，粒子權(quán)重在兩個方面影響DSMC的碰撞計(jì)算：碰撞對數(shù)目和二體碰撞機(jī)制?；贜TC(Non-Time-Counter)方法[26]，文獻(xiàn)[12]提出將電子作為一個組P，而將除電子以外的其它組分作為另一個組Q；不考慮P組內(nèi)電子的相互碰撞，于是碰撞可以分為兩類，即Q組內(nèi)粒子之間的相互碰撞以及P組與Q組粒子之間的相互碰撞。同時，統(tǒng)計(jì)P組的電子數(shù)目NP和Q組的粒子數(shù)目NQ，于是網(wǎng)格i中Q組內(nèi)的粒子碰撞對數(shù)為

P組與Q組粒子之間的分子碰撞對數(shù)為

Wj為組分j粒子對應(yīng)的權(quán)重。式(5)中的NP亦須包含電子的權(quán)重，即

粒子權(quán)重對DSMC碰撞計(jì)算的第二個方面是二體碰撞機(jī)制。當(dāng)兩個粒子權(quán)重相同時，碰后的粒子速度為

Um為碰撞粒子對的質(zhì)心速度，g為碰撞粒子對的相對速度。當(dāng)兩個粒子權(quán)重不同時，碰撞后較低權(quán)重粒子的速度由式(8a)或式(8b)確定，較高權(quán)重粒子則僅有一部分速度發(fā)生變化(由式(8a)或式(8b)確定)；較高權(quán)重粒子速度發(fā)生變化的概率為

粒子的碰撞后的內(nèi)能狀態(tài)類似確定。

上述方法并不保證每一步的動量和能量守恒，但是在足夠大量的碰撞時可以保持動量和能量近似平均守恒。

2.2 權(quán)重因子對化學(xué)反應(yīng)的影響

采用上述權(quán)重因子法的根本目的是提高稀有組分的仿真分子數(shù)目，進(jìn)而抑制其統(tǒng)計(jì)漲落，得到稀有組分較為光滑的分布等值線。依據(jù)上述處理方式，稀有組分的仿真分子數(shù)目擴(kuò)大到1/Wr倍。

3 典型再入的稀有組分權(quán)重因子效果驗(yàn)證

3.1 RAM-C II飛行試驗(yàn)稀薄電離過程計(jì)算驗(yàn)證

RAM-C II飛行試驗(yàn)是目前可查到的僅有飛行試驗(yàn)電子數(shù)密度測量公開案例，常作為檢驗(yàn)電離計(jì)算模型和程序的經(jīng)典算例。文獻(xiàn)[12]考察了不使用稀有組分權(quán)重因子方法的電離過程模擬，結(jié)果顯示盡管在測點(diǎn)方向的最大電子數(shù)密度與試驗(yàn)測量值吻合較好，球錐體周圍流場的電子數(shù)密度分布存在很大統(tǒng)計(jì)波動。采用與文獻(xiàn)[12]完全相同的計(jì)算條件，取稀有組分權(quán)重因子Wr為0.01。圖1給出了是否采用稀有組分權(quán)重因子方法得到的宏觀流場參數(shù)比較情況。很明顯，使用稀有組分權(quán)重因子方法與否，得到的流場宏觀參數(shù)幾乎不存在差異，意味著電離導(dǎo)致的稀有組分對宏觀流場參數(shù)幾乎不存在影響，也從一個方面印證了第一宇宙速度下常規(guī)氣動特性計(jì)算通常不考慮電離過程在工程意義上的合理性。

然而，若考察電子數(shù)密度分布等再入物理領(lǐng)域極為關(guān)心的核心問題，權(quán)重格式使用與否，對于能否得到光滑的流場電子數(shù)密度分布，有著極為重要的影響。圖2給出了使用稀有組分權(quán)重因子方法與否得到的電子數(shù)密度分布云圖。不使用權(quán)重因子方法時，在電子密度較高的頭部附近，電子數(shù)密度等值線較為光滑，但是電子數(shù)密度較低的身部和尾部區(qū)域，巨大的統(tǒng)計(jì)漲落直接淹沒了電子數(shù)密度等值線。使用權(quán)重因子方法時，上述情形得到了極大改善，電子數(shù)密度較高和較低的區(qū)域均能得到光滑的等值線。

(a) 粒子數(shù)密度

(b) 全局溫度

(c) 壓強(qiáng)

(d) 馬赫數(shù)

圖3給出了使用稀有組分權(quán)重因子方法與否得到的沿測點(diǎn)方向最大電子數(shù)密度比較情況。該圖表明，在不使用稀有組分權(quán)重方法時，盡管總體上DSMC計(jì)算值與試驗(yàn)測量值較好吻合，在中部兩個測點(diǎn)仍存在一些偏差；使用稀有組分權(quán)重方法以后，上述偏離能得以較好矯正，計(jì)算值與試驗(yàn)測量值明顯更為靠近。甚至可以認(rèn)為，不使用稀有組分權(quán)重因子時的偏差是因?yàn)榻y(tǒng)計(jì)漲落造成的，事實(shí)上圖2上半部分也確實(shí)存在極其明顯的統(tǒng)計(jì)漲落，這一漲落在使用權(quán)重因子方法時可以得到有效抑制，進(jìn)而得到與試驗(yàn)測量值更為靠近的計(jì)算結(jié)果。

圖2 使用與不使用稀有組分權(quán)重因子方法得到的RAM-C II 81km電子數(shù)密度分布比較Fig.2 Comparison of electron number density distributions between using rare species weighting scheme or not for RAM-C II at 81 km

圖3 使用與不使用稀有組分權(quán)重因子方法得到的沿RAM-C II 81 km測點(diǎn)方向最大電子數(shù)密度比較Fig.3 Comparison of maximum electron number density along the probe direction between using rare species weighting scheme or not for RAM-C II at 81 km

3.2 彗星探測器Stardust再入稀薄段電離過程計(jì)算驗(yàn)證

彗星探測器Stardust[27]實(shí)現(xiàn)了迄今人造飛行器的最高速度, 它以12.8 km/s的速度再入大氣層，必將引起強(qiáng)烈的化學(xué)反應(yīng)和電離過程。采用與文獻(xiàn)[12]完全相同的計(jì)算條件，取稀有組分權(quán)重因子Wr為0.1。圖4給出是否采用稀有組分權(quán)重因子方法得到的宏觀流場參數(shù)比較情況，與RAM-C II算例的情況一樣，稀有組分權(quán)重因子的使用不會給宏觀流場參數(shù)分布帶來影響。

圖5給出了核心關(guān)注的Stardust在80 km高度的電子數(shù)密度分布情況。該圖表明，對于12.8 km/s量級的再入速度，不使用權(quán)重因子方法時，DSMC方法計(jì)算得到的電子數(shù)密度分布等值線在流場大部分區(qū)域已經(jīng)較為光滑，但在尾流的低密度區(qū)仍存在明顯的統(tǒng)計(jì)波動。權(quán)重因子方法的使用，使得上述波動被進(jìn)一步抑制，得到的等值線光滑性更好。事實(shí)上，相對于文獻(xiàn)[12]，使用權(quán)重因子方法得到的電子數(shù)密度與文獻(xiàn)[27]相比，數(shù)據(jù)吻合性更好，如圖6所示。如引言所述，電子等稀有組分?jǐn)?shù)密度的統(tǒng)計(jì)波動可能導(dǎo)致不真實(shí)的化學(xué)反應(yīng)計(jì)算，其原因是統(tǒng)計(jì)波動給網(wǎng)格中碰撞分子對計(jì)算帶來一定程度的不確定性，這也是圖5中壁面附近電子數(shù)密度顯示出一定差異的可能原因。

(a) 粒子數(shù)密度

(b) 全局溫度

(c) 壓強(qiáng)

(d) 馬赫數(shù)

圖5 使用與不使用稀有組分權(quán)重因子方法得到的Stardust 80 km電子密度分布比較Fig.5 Comparison of electron number density distributions between using rare species weighting scheme or not for Stardust at 80 km

圖6 80 km處Stardust繞流電子密度與參考文獻(xiàn)對比Fig.6 Comparison of electron density for Stardust at 80 km

4 結(jié) 論

在前期工作[12]的基礎(chǔ)上，提出了一種處理含電離化學(xué)反應(yīng)DSMC模擬的稀有組分權(quán)重因子方法。該方法對DSMC二體碰撞的影響包括碰撞分子對數(shù)的確定和碰撞機(jī)制兩個方面，前者的修正需考慮權(quán)重因子，后者則根據(jù)權(quán)重因子之比確定的概率進(jìn)行狀態(tài)確定。對于空氣11組元的含電離化學(xué)反應(yīng)，歸類為四種情況分別處理，基本思想是根據(jù)權(quán)重因子對生成的稀有組分粒子進(jìn)行復(fù)制、對反應(yīng)物中的常規(guī)組分按概率保留或刪除。數(shù)值測試表明，在諸如RAM-C II飛行條件的極弱電離情況下，權(quán)重因子方法的使用使得電子數(shù)密度等值線光滑性大幅改善；在Stardust飛行條件的強(qiáng)電離情況下，權(quán)重因子方法的使用也可以改善低電子數(shù)密度區(qū)域等值線的光滑性，得到的電子數(shù)密度分布與參考文獻(xiàn)值吻合較好。同時，權(quán)重因子方法的引入不會給宏觀流場參數(shù)計(jì)算帶來影響。