球錐鈍頭體再入稀薄氣體電離過(guò)程三維DSMC模擬與驗(yàn)證

方 明, 李志輝,,*, 李中華, 田 穎

(1. 北京航空航天大學(xué) 國(guó)家計(jì)算流體力學(xué)實(shí)驗(yàn)室, 北京 100191;2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 超高速空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000)

球錐鈍頭體再入稀薄氣體電離過(guò)程三維DSMC模擬與驗(yàn)證

方 明1, 李志輝1,2,*, 李中華2, 田 穎2

(1. 北京航空航天大學(xué) 國(guó)家計(jì)算流體力學(xué)實(shí)驗(yàn)室, 北京 100191;2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 超高速空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000)

采用增大電子質(zhì)量三個(gè)數(shù)量級(jí)并相應(yīng)調(diào)整離子質(zhì)量的方法，拓展化學(xué)反應(yīng)的DSMC仿真方法處理稀薄氣體電離過(guò)程；采用單溫度模型處理全部化學(xué)反應(yīng)，修正涉及電子的反應(yīng)速率常數(shù)以保證真實(shí)化學(xué)反應(yīng)速率；以直角/非結(jié)構(gòu)網(wǎng)格相結(jié)合，運(yùn)用碰撞網(wǎng)格自適應(yīng)技術(shù)，基于MPI并行環(huán)境，開(kāi)發(fā)適用于真實(shí)復(fù)雜外形的三維稀薄氣體電離DSMC計(jì)算程序。對(duì)RAM-C II外形的再入繞流稀薄氣體電子密度進(jìn)行模擬驗(yàn)證，所得結(jié)果與飛行試驗(yàn)測(cè)量值吻合較好；對(duì)Stardust外形的再入稀薄段電離特性數(shù)值仿真分析，電子密度等值線云圖與參考文獻(xiàn)結(jié)果一致。相較于稀薄氣體不含電離反應(yīng)的DSMC方法，本文發(fā)展的模型和程序不會(huì)導(dǎo)致計(jì)算量的顯著增大，可直接應(yīng)用于三維復(fù)雜外形體極高速再入條件下的稀薄氣體電離計(jì)算，為工程設(shè)計(jì)提供技術(shù)支持和指導(dǎo)作用。計(jì)算結(jié)果表明，極高速再入條件下傳統(tǒng)稀薄流區(qū)的電子數(shù)密度足以引起通信黑障，需在通信設(shè)計(jì)上給予高度關(guān)注。

航天器再入；稀薄氣體；化學(xué)反應(yīng)；電離效應(yīng)；DSMC方法

0 引 言

隨著載人航天的深化發(fā)展、探月工程的逐步推進(jìn)和深空探測(cè)的日臻實(shí)踐，航天器再入速度持續(xù)增大，燒蝕、通信黑障等超高速典型現(xiàn)象愈加顯著，為其氣動(dòng)力/熱和等離子體特性預(yù)測(cè)提出了嚴(yán)峻挑戰(zhàn)。神舟飛船等近地軌道航天器以第一宇宙速度7.9 km/s再入大氣層，氣體發(fā)生強(qiáng)烈的非平衡化學(xué)反應(yīng)和微弱電離，通信黑障主要發(fā)生在高空近連續(xù)流區(qū)域；2014年秋發(fā)射并成功返回的探月飛行試驗(yàn)器返回器以11.2 km/s的速度蛙跳式再入大氣層，非平衡化學(xué)反應(yīng)和電離現(xiàn)象更為嚴(yán)重，通信黑障大幅向稀薄流區(qū)延伸；彗星探測(cè)器Stardust[1]實(shí)現(xiàn)了迄今人造飛行器的最高速度，它以12.8 km/s的速度再入大氣層，非平衡化學(xué)反應(yīng)和電離極為嚴(yán)重，通信黑障覆蓋寬廣的稀薄區(qū)域。

傳統(tǒng)的空氣動(dòng)力學(xué)[2]研究重點(diǎn)關(guān)注連續(xù)流區(qū)的氣動(dòng)力/熱問(wèn)題，百余年來(lái)極大地促進(jìn)了航空航天事業(yè)的發(fā)展。但是，隨著稀薄間段粒子效應(yīng)的增強(qiáng)，連續(xù)介質(zhì)假設(shè)失效，理論分析和數(shù)值計(jì)算難以實(shí)現(xiàn)上述極高速再入過(guò)程、尤其是稀薄段電離特性的有效預(yù)測(cè)；同時(shí)，受設(shè)備能力和測(cè)試手段的限制，地面風(fēng)洞和模型飛行試驗(yàn)也無(wú)法完成上述高溫高速非平衡條件的真實(shí)再現(xiàn)。近年來(lái)針對(duì)稀薄氣體發(fā)展的多種方法中，基于Boltzmann方程的物理分析和建模計(jì)算[3-4]目前還不能包含電離過(guò)程，高超聲速低密度風(fēng)洞雖可以實(shí)現(xiàn)稀薄卻無(wú)法達(dá)到真實(shí)的高溫高速環(huán)境，只有基于粒子隨機(jī)統(tǒng)計(jì)模擬的DSMC(Direct Simulation Monte Carlo)方法[5-9]是目前最有可能從數(shù)值上研究解決上述問(wèn)題的一種手段。

1987年，G.A. Bird[10]首次將其提出的DSMC方法用于稀薄氣體的電離模擬，這一工作開(kāi)創(chuàng)了稀薄氣體電離過(guò)程DSMC模擬研究的先河。他用41個(gè)化學(xué)反應(yīng)處理11組元空氣涉及到的電離，奠定了此后稀薄氣體電離模擬的基本框架。在這一工作中，基于流場(chǎng)電中性假設(shè)，Bird將電子與離子進(jìn)行強(qiáng)制關(guān)聯(lián)，這一做法避免了精確追蹤電子運(yùn)動(dòng)的繁瑣計(jì)算，但是電子運(yùn)動(dòng)嚴(yán)格受限的處理在物理上是不真實(shí)的。限于當(dāng)時(shí)的計(jì)算能力，Bird僅用一維駐點(diǎn)線程序模擬AFE的飛行條件。Carlson和Hassan[11]的工作也是基于一維駐點(diǎn)線流動(dòng)。由于一維駐點(diǎn)線模擬的本質(zhì)缺陷，這一方法無(wú)法考慮飛行器的幾何外形，需要指定邊界層的邊緣條件。也正是這一缺陷，Bird沒(méi)有成功實(shí)現(xiàn)與試驗(yàn)結(jié)果的比較[12]。事實(shí)上，激波脫體距離與飛行器的幾何外形直接強(qiáng)烈關(guān)聯(lián)，即使所有其它計(jì)算參數(shù)完全相同，不同幾何外形飛行器的流場(chǎng)結(jié)果存在顯著差異，駐點(diǎn)線上的流動(dòng)也明顯不同，真實(shí)計(jì)算必須考慮飛行器的幾何外形。

近年來(lái)，稀薄氣體電離過(guò)程的二維/軸對(duì)稱模擬結(jié)果陸續(xù)發(fā)表，但是鮮有三維復(fù)雜外形的公開(kāi)報(bào)道。Ozawa[1]采用軸對(duì)稱的Stardust外形，比較了Stardust極高速再入狀態(tài)下電離的DSMC和CFD計(jì)算結(jié)果，發(fā)現(xiàn)兩種方法差異很大。Boyd[12-13]基于軸對(duì)稱程序，采用稀有組分追蹤算法和電離與離子強(qiáng)制關(guān)聯(lián)的方法模擬分析了RAM-C II的再入電離情況。樊菁[14]采用稀有組分分離方法、基于軸對(duì)稱程序模擬了RAM-C II的再入電離過(guò)程。Shevyrin[15]不將電子作為一種粒子模擬，而是根據(jù)等離子體電中性假設(shè)將電子密度等于離子密度之和，在此基礎(chǔ)上比較了不同模型對(duì)于電離反應(yīng)的區(qū)別，其工作也是基于RAM-C II外形。Morsa[16]基于DS2V程序，針對(duì)Orin外形，比較了不同化學(xué)反應(yīng)模型涉及電離計(jì)算的結(jié)果，推薦在防熱設(shè)計(jì)中采用Gupta模型[17]。上述工作均基于二維或者軸對(duì)稱程序，僅適用于簡(jiǎn)單軸對(duì)稱外形。同時(shí)，它們均基于電子-離子強(qiáng)制關(guān)聯(lián)假設(shè)，而這一假設(shè)在高維情況下帶來(lái)計(jì)算量的顯著增大。

上述方法的另一個(gè)缺陷是化學(xué)反應(yīng)的多溫度模型，亦即在不同化學(xué)反應(yīng)的實(shí)現(xiàn)過(guò)程中使用平動(dòng)、振動(dòng)和電子溫度等多種溫度模態(tài)或者它們之間的組合。這一處理在不涉及電離的反應(yīng)時(shí)是可行的，但是由于電子是稀有組分，很多時(shí)候一個(gè)網(wǎng)格中甚至只有一個(gè)或者完全沒(méi)有電子，電子溫度計(jì)算必然遭遇強(qiáng)烈的統(tǒng)計(jì)漲落。這也是公開(kāi)報(bào)道結(jié)果鮮有三維的根本原因，即三維情況每個(gè)碰撞網(wǎng)格的模擬分子數(shù)目必然受到強(qiáng)烈限制，不可能像二維那樣達(dá)到數(shù)十甚至上百。理想的情況是，對(duì)于全部化學(xué)反應(yīng)使用單一溫度模型，最好是氣體全局溫度，這正是本文的工作之一。

從工程角度而言，一個(gè)有啟發(fā)意義的工作來(lái)自Boyd[18]。它將電子質(zhì)量人為增大了三個(gè)數(shù)量級(jí)，把電子當(dāng)作一種普通分子對(duì)待。在這種條件下，電子的熱運(yùn)動(dòng)速度與其它分子(原子、分子和離子)不存在量級(jí)上的差異，亦即電子將與分子整體上同步運(yùn)動(dòng)。同時(shí)，由于電離度較低，這一調(diào)整不會(huì)對(duì)流場(chǎng)結(jié)構(gòu)產(chǎn)生明顯影響。

基于上述總結(jié)和啟發(fā)，本文嘗試拓展傳統(tǒng)DSMC方法以處理極高速再入條件下的稀薄氣體電離過(guò)程，發(fā)展適用于復(fù)雜外形的三維稀薄氣體電離模擬程序，為有效預(yù)測(cè)復(fù)雜外形航天器的稀薄電離特性提供一種可能。

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

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

對(duì)于一般的二體碰撞化學(xué)反應(yīng)，通常可以寫(xiě)成如下形式：

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

化學(xué)反應(yīng)速率常數(shù)通常表示為如下Arrhenius形式的方程：

其中Λ和η為常數(shù)，Ea為反應(yīng)活化能，κ為Boltzmann常數(shù)，T為溫度。如前所述，不同反應(yīng)的Τ選取為不同溫度模態(tài)或者它們的組合，而電子溫度計(jì)算必然遭遇強(qiáng)烈統(tǒng)計(jì)漲落。本文提出：對(duì)全部化學(xué)和電離反應(yīng)，采用單一溫度模型，即氣體當(dāng)?shù)厝譁囟?。如此，?dāng)碰撞總能量Ec>Ea時(shí)，上述常數(shù)Λ和η通過(guò)下式與化學(xué)反應(yīng)概率P相關(guān)聯(lián)：

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

1.2 能量分配模型

廣泛應(yīng)用于DSMC方法中描述分子平動(dòng)能和內(nèi)能分配的模型是Larsen-Borgnakke (L-B)模型[5]。在這個(gè)模型中，假設(shè)轉(zhuǎn)動(dòng)模式的能譜是連續(xù)的，碰撞后分子的能量根據(jù)Maxwell平衡態(tài)分布函數(shù)在平動(dòng)模式和內(nèi)能模式之間重新分配。

由于涉及電離的反應(yīng)包括離解反應(yīng)、置換反應(yīng)、聯(lián)合電離反應(yīng)及其逆反應(yīng)、直接碰撞電離反應(yīng)、離子交換反應(yīng)等多種類型，有必要對(duì)不同類型的碰后內(nèi)能分配作簡(jiǎn)單說(shuō)明。

對(duì)于離解反應(yīng)，以氮?dú)鉃槔?/p>

分解為兩個(gè)過(guò)程：

和

反應(yīng)熱在式(5)中予以考慮，同時(shí)按照上述L-B能量分配規(guī)則將總碰撞能在相對(duì)平動(dòng)能、轉(zhuǎn)動(dòng)能和振動(dòng)能之間分配。對(duì)于式(6)，將N2(post)的內(nèi)能轉(zhuǎn)化為兩個(gè)N原子的相對(duì)平動(dòng)能。

對(duì)于直接電離反應(yīng)，以O(shè)為例：

也類似地分解為兩個(gè)過(guò)程：

和

需要指出的是，式(8)和式(9)均不涉及內(nèi)能的重新分配，因?yàn)檠踉雍碗娮泳痪邆鋬?nèi)自由度，將其作為兩個(gè)步驟處理是便于實(shí)施二體碰撞。

聯(lián)合電離反應(yīng)是兩個(gè)原子(不具備內(nèi)自由度)生成一個(gè)離子(有內(nèi)自由度)和一個(gè)電子(不具備內(nèi)自由度)，理論上來(lái)說(shuō)可以用L-B方法在生成物之間分配平動(dòng)能和內(nèi)能，但是實(shí)際操作中可以直接設(shè)置離子的轉(zhuǎn)動(dòng)能和振動(dòng)能為0，將全部能量指定為相對(duì)平動(dòng)能。聯(lián)合電離逆反應(yīng)則有由離子和電子反應(yīng)生成兩個(gè)原子，全部碰撞總能量指定為兩個(gè)原子的相對(duì)平動(dòng)能。

對(duì)于置換和離子交換反應(yīng)，過(guò)程較為簡(jiǎn)單，直接根據(jù)L-B方法對(duì)生成物進(jìn)行能量分配即可。

1.3 碰撞分子對(duì)抽樣方法

DSMC方法需要正確確定參與碰撞的分子對(duì)數(shù)目。在均質(zhì)氣體中，分子之間的碰撞概率與相對(duì)速度cr和總碰撞截面σT的乘積成正比。記網(wǎng)格體積為Vc，網(wǎng)格內(nèi)每個(gè)模擬分子代表FN個(gè)真實(shí)分子，時(shí)間步長(zhǎng)為Δt。

通常，兩種不同方式用于計(jì)算網(wǎng)格內(nèi)的碰撞分子對(duì)數(shù)。第一種方式是籠統(tǒng)地計(jì)算網(wǎng)格內(nèi)各種組分分子的碰撞對(duì)總和，即碰撞分子對(duì)數(shù)目為：

該方式的最大的問(wèn)題是Ni和Nj都不會(huì)很大而且漲落很大，得到的碰撞分子對(duì)數(shù)目與真實(shí)情況可能差之甚遠(yuǎn)。

本文采用分組處理的方式：將電子作為一個(gè)組P，而將除電子以外的其它組分作為另一個(gè)組Q；電子之間沒(méi)有化學(xué)反應(yīng)發(fā)生，同時(shí)電子濃度很低，可以不考慮P組內(nèi)電子的相互碰撞；于是碰撞可以分為兩類，即Q組內(nèi)分子之間的相互碰撞以及P組與Q組分子之間的相互碰撞。同時(shí)，統(tǒng)計(jì)P組的分子數(shù)目NP和Q組的分子數(shù)目NQ，于是Q組內(nèi)的分子碰撞對(duì)數(shù)為：

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

1.4 等離子體效應(yīng)

等離子體效應(yīng)是稀薄氣體電離過(guò)程DSMC模擬的一個(gè)核心困難，本文擬采用增大電子質(zhì)量三個(gè)數(shù)量級(jí)的方式對(duì)這一問(wèn)題進(jìn)行簡(jiǎn)化處理。電子的真實(shí)質(zhì)量為9.11×10-31kg，增大三個(gè)數(shù)量級(jí)至9.11×10-28kg。以氧原子O為例，其質(zhì)量為2.656×10-26kg，調(diào)整后的電子質(zhì)量約為氧原子質(zhì)量的3.4%。在溫度相同的情況下，調(diào)整后電子的隨機(jī)熱運(yùn)動(dòng)速度均值比氧原子大不到一個(gè)數(shù)量級(jí)。同時(shí)，為保證化學(xué)反應(yīng)過(guò)程質(zhì)量守恒，增大電子質(zhì)量須同步減小離子質(zhì)量。以氧離子O+為例，真實(shí)質(zhì)量為2.656×10-26kg ，調(diào)整后的質(zhì)量為2.5649×10-26kg ，質(zhì)量變化約為3.4%，亦不會(huì)對(duì)其熱運(yùn)動(dòng)速度造成明顯影響。調(diào)整后各離子的質(zhì)量及其因調(diào)整而改變的比例見(jiàn)表2.該表表明，離子質(zhì)量改變率的最大值為N+的3.9%，對(duì)其熱運(yùn)動(dòng)速度的改變亦在可以忽略的范圍內(nèi)。

表2 調(diào)整后的離子質(zhì)量及其變化情況Table 2 Modified ion mass and corresponding vibration rate

經(jīng)過(guò)上述處理之后，由于電子和離子的運(yùn)動(dòng)速度基本處于一個(gè)量級(jí)，整體上它們將一起作宏觀運(yùn)動(dòng)，因此既無(wú)需將電子和離子進(jìn)行強(qiáng)制關(guān)聯(lián)，也無(wú)需擔(dān)心電子因運(yùn)動(dòng)速度過(guò)快而疾速逃離計(jì)算區(qū)域，計(jì)算區(qū)域可維持近似電中性。

2 復(fù)雜外形航天器再入稀薄電離過(guò)程的DSMC程序設(shè)計(jì)

2.1 直角/非結(jié)構(gòu)網(wǎng)格技術(shù)

在DSMC方法中，網(wǎng)格的作用是選擇可能的碰撞分子對(duì)以及對(duì)宏觀參數(shù)取樣。從計(jì)算的角度而言，流場(chǎng)采用大小均勻的直角坐標(biāo)網(wǎng)格，追蹤分子的效率非常高，計(jì)算區(qū)域內(nèi)的模擬分子可以直接根據(jù)分子的位置坐標(biāo)來(lái)確定分子所屬網(wǎng)格，不必跟蹤分子從一個(gè)網(wǎng)格到另一個(gè)網(wǎng)格的運(yùn)動(dòng)，但是其缺點(diǎn)是無(wú)法精確地描述物面邊界條件。非結(jié)構(gòu)網(wǎng)格在表征復(fù)雜外形方面具有高度適用性，并且容易實(shí)現(xiàn)正確的邊界條件，但是追蹤分子的算法比較繁瑣和費(fèi)時(shí)。本文結(jié)合直角坐標(biāo)計(jì)算的高效率以及三角形非結(jié)構(gòu)網(wǎng)格對(duì)航天器幾何外形的精細(xì)描述，在描述物面幾何形狀的非結(jié)構(gòu)網(wǎng)格建立以后，直接將其嵌入到直角網(wǎng)格的流場(chǎng)中，使DSMC計(jì)算對(duì)流場(chǎng)網(wǎng)格的依賴程度大大降低。

直角/非結(jié)構(gòu)網(wǎng)格的詳細(xì)陳述見(jiàn)文獻(xiàn)[19]。

2.2 網(wǎng)格自適應(yīng)方法

DSMC方法的內(nèi)在要求之一就是網(wǎng)格尺度須小于氣體分子平均自由程。流場(chǎng)初始化時(shí)的網(wǎng)格設(shè)置是其尺度小于來(lái)流平均自由程，但是在超高速飛行的航天器頭部，激波及波后密度顯著增大，當(dāng)?shù)仄骄杂沙套冃?，初始化的網(wǎng)格不能繼續(xù)滿足要求，需要對(duì)網(wǎng)格進(jìn)行細(xì)化。本文采用兩級(jí)直角坐標(biāo)網(wǎng)格的策略，一級(jí)網(wǎng)格為初始化時(shí)形成的背景網(wǎng)格，二級(jí)網(wǎng)格為根據(jù)密度變化將背景網(wǎng)格加密之后形成的碰撞網(wǎng)格。

網(wǎng)格自適應(yīng)參數(shù)通常有不同選取方式，較為普遍的是根據(jù)流場(chǎng)密度梯度變化確定，但是這種方式涉及流場(chǎng)密度的梯度計(jì)算。本文采用一種簡(jiǎn)化處理方式：初始時(shí)背景網(wǎng)格即為碰撞網(wǎng)格，設(shè)置碰撞網(wǎng)格分子數(shù)上限；隨著流場(chǎng)演化，當(dāng)碰撞網(wǎng)格中分子數(shù)大于上限時(shí)，將該網(wǎng)格在每個(gè)方向一分為二，即將該網(wǎng)格細(xì)分為八個(gè)網(wǎng)格；根據(jù)流場(chǎng)演化情況，不斷重復(fù)上述過(guò)程，直到流場(chǎng)達(dá)到穩(wěn)定狀態(tài)。

2.3 DSMC并行算法

考慮到極高速航天器再入稀薄和近連續(xù)過(guò)渡流電離過(guò)程的DSMC模擬計(jì)算量巨大，單一計(jì)算機(jī)計(jì)算速度慢且無(wú)法滿足大量模擬分子導(dǎo)致的巨大內(nèi)存需求，必須發(fā)展DSMC并行算法。本文采用區(qū)域分解的策略，根據(jù)計(jì)算的處理器數(shù)目將計(jì)算區(qū)域劃分為等量的子區(qū)域，每個(gè)處理器在其分配的子區(qū)域內(nèi)獨(dú)立地計(jì)算模擬分子的碰撞和遷移，離開(kāi)子區(qū)域的模擬分子把攜帶的信息傳遞給對(duì)應(yīng)子區(qū)域的處理器。

如果模擬分子在運(yùn)動(dòng)一個(gè)時(shí)間步長(zhǎng)之后，離開(kāi)原先所在的網(wǎng)格遷移進(jìn)入相鄰的網(wǎng)格，并且改變后的網(wǎng)格隸屬于其它的處理器，這樣的分子及其攜帶的所有信息將被記錄下來(lái)并放入緩沖數(shù)組。當(dāng)所有的模擬分子被確定以后，放入緩沖數(shù)組的分子就在相應(yīng)的處理器之間進(jìn)行信息交換。

本文基于MPI并行程序設(shè)計(jì)函數(shù)，采用單一方向區(qū)域剖分，實(shí)現(xiàn)三維復(fù)雜外形航天器極高速再入稀薄氣體電離過(guò)程的DSMC模擬。三維復(fù)雜外形航天器極高速再入稀薄氣體電離過(guò)程的DSMC模擬，依據(jù)航天器尺寸和來(lái)流條件，通常需要數(shù)十至數(shù)百個(gè)計(jì)算核心。

3 航天器再入稀薄電離過(guò)程的DSMC模擬

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

RAM-C II飛行試驗(yàn)是目前可查到的僅有飛行試驗(yàn)電子密度測(cè)量公開(kāi)案例，常作為檢驗(yàn)電離計(jì)算模型和程序的經(jīng)典算例。它于1960年代后期發(fā)射，以第一宇宙速度再入大氣層，其外形如圖1所示。飛行試驗(yàn)中，沿著4個(gè)測(cè)點(diǎn)R1、R2、R3和R4測(cè)量這些方向上的最大電子密度值，沿著L1測(cè)量該方向上的電子密度分布情況。

受網(wǎng)格尺度限制，本文計(jì)算狀態(tài)采用81 km高度，此處來(lái)流條件為：粒子數(shù)密度為3.275×1020/m3，其中21%為O2、79%為N2，壓強(qiáng)為0.889 Pa；來(lái)流溫度為196.7 K，對(duì)應(yīng)聲速為281.15 m/s；來(lái)流速度為7800 m/s，對(duì)應(yīng)馬赫數(shù)28.1；氣體分子平均自由程為5.16×10-3m。壁面采用漫反射模型，溫度取為1000 K。取背景網(wǎng)格寬度5×10-3m，計(jì)算區(qū)域?yàn)閇-0.25 m,1.55 m]×[-0.5 m,0.5 m]×[0,0.5 m]，于是計(jì)算區(qū)域的網(wǎng)格數(shù)為360×200×100。采用60個(gè)CPU核進(jìn)行并行計(jì)算，取時(shí)間步長(zhǎng)為1×10-7s。

相對(duì)于RAM-C II氣動(dòng)特性和流場(chǎng)特征，本文最為關(guān)注的是RAM-C II周圍的電子密度分布情況。RAM-C II繞流流場(chǎng)的電子密度分布如圖2所示。該圖表明，盡管經(jīng)過(guò)較大樣本的統(tǒng)計(jì)抽樣，電子密度分布依然存在較大統(tǒng)計(jì)漲落，尤其是在尾流區(qū)域。造成這一現(xiàn)象主要有兩個(gè)原因：一是DSMC方法不可回避的統(tǒng)計(jì)漲落，由于電子屬于稀有組分，解決這一問(wèn)題的根本方法是發(fā)展稀有組分粒子追蹤方法[20-21]；另一個(gè)原因是本算例的特殊性，81 km高度屬于稀薄過(guò)渡流區(qū)，粒子密度相對(duì)較小，同時(shí)RAM-C II以第一宇宙速度再入，電離度很低，電子密度本身就很小。計(jì)算經(jīng)驗(yàn)

表明，隨著再入速度的增大，統(tǒng)計(jì)波動(dòng)將被明顯抑制，這也是RAM-C II算例經(jīng)常被作為電離相關(guān)計(jì)算方法和程序考核標(biāo)準(zhǔn)的原因之一。

四個(gè)測(cè)點(diǎn)R1、R2、R3和R4方向的電子密度最大值如圖3所示。一般認(rèn)為，電子密度計(jì)算值誤差在3倍以內(nèi)都是可以接受的。該圖表明，盡管沒(méi)有采用稀薄組分粒子追蹤方法，本文發(fā)展的模型和計(jì)算方法得到的電子密度與飛行試驗(yàn)測(cè)量值仍能吻合較好。

3.2 彗星探測(cè)器Stardust再入稀薄段電離特性計(jì)算分析

如前所述，彗星探測(cè)器Stardust實(shí)現(xiàn)了人造飛行器的最高飛行速度，它以12.8 km/s的速度再入大氣層，必將引起強(qiáng)烈的化學(xué)反應(yīng)和電離過(guò)程。Stardust的外形如圖4所示，其球頭半徑為0.22 m，頭部半錐角為59.5°，肩部半徑為0.02 m，前后長(zhǎng)度為0.5 m，肩部半徑為尾部半徑為0.25 m。為與參考文獻(xiàn)對(duì)比，本文計(jì)算80 km的高度，對(duì)應(yīng)來(lái)流馬赫數(shù)為45.5。同上例采用壁面漫反射模型，壁溫取1000 K。計(jì)算區(qū)域?yàn)閇-0.3 m,0.6 m]×[-0.6 m,0.6 m]×[0,0.6 m]。采用36個(gè)核進(jìn)行并行計(jì)算。

本算例關(guān)心的核心問(wèn)題依然是繞流流場(chǎng)的電子密度分布，如圖5所示。該圖表明，本文計(jì)算所得的電子密度分布與參考文獻(xiàn)[1]吻合很好。盡管在等值線分布上存在一些細(xì)微差異，兩者在電子密度分布的脫體距離、最大值上均完全吻合。該圖還表明，對(duì)于Stardust這種以極高速度再入的航天器，其在80 km高度處的電子密度高達(dá)5×1020/ m3，較上例顯著增大。這意味著此類飛行器的通信黑障區(qū)域?qū)⒋蠓蛳”^(qū)域延伸，必須引起通信設(shè)計(jì)的高度關(guān)注。

4 結(jié) 論

本文采用增大電子質(zhì)量三個(gè)數(shù)量級(jí)并相應(yīng)調(diào)整離子質(zhì)量的方法，把電子當(dāng)作一種普通分子對(duì)待，拓展DSMC模擬化學(xué)反應(yīng)方法處理稀薄氣體電離過(guò)程；采用單溫度模型處理全部化學(xué)反應(yīng)，對(duì)涉及電子的反應(yīng)速率常數(shù)進(jìn)行修正；以直角/非結(jié)構(gòu)網(wǎng)格相結(jié)合，運(yùn)用碰撞網(wǎng)格自適應(yīng)方法，基于MPI并行環(huán)境，開(kāi)發(fā)適于真實(shí)復(fù)雜外形的三維稀薄氣體電離DSMC計(jì)算程序。通過(guò)計(jì)算RAM-C II外形再入飛行高度81 km繞流流場(chǎng)的電子密度，得到計(jì)算結(jié)果與飛行試驗(yàn)測(cè)量值吻合較好；計(jì)算Stardust外形以類第二宇宙速度再入飛行高度80 km繞流流場(chǎng)的電子密度等值線云圖，得到了與參考文獻(xiàn)計(jì)算值一致的結(jié)果。相對(duì)于稀薄氣體化學(xué)反應(yīng)的DSMC方法，本文發(fā)展的模型和程序不會(huì)導(dǎo)致計(jì)算量的顯著增大，可直接應(yīng)用于三維復(fù)雜外形體極高速再入條件下的稀薄氣體電離計(jì)算。

計(jì)算經(jīng)驗(yàn)表明，RAM-C II飛行試驗(yàn)是檢驗(yàn)電離模型和算法的典型算例。彗星探測(cè)器Stardust在80 km高度以12.8 km/s速度再入時(shí)，流場(chǎng)電子密度高達(dá)5×1020/m3，這意味著此類飛行器的通信黑障區(qū)域?qū)⒋蠓蛳”^(qū)域延伸，需要引起通信設(shè)計(jì)的高度關(guān)注，也為本文進(jìn)一步深化研究指明了方向。

[1]Ozawa T, Zhong J, Levin D A, et al. Modeling of stardust reentry flows with ionization in DSMC[R]. AIAA 2007-611.

[2]Anderson J D. Fundamentals of aerodynamics (5thedition)[M]. McGraw-Hill Education, 2010.

[3]Chapman S, Cowling T G. The mathematiccal theory of no-uniform gases (3rdedition)[M]. Cambridge University Press, 1991.

[4]Li Z H, Zhang H X. Study on gas kinetic unified algorithm for flows from rarefied transition to continuum[J]. Journal of Computational Physics, 2003, 193(2):708-738.

[5]Bird G A. Molecular gas dynamics and the direct simulation of gas flows[M]. Clarendon Press, Oxford, 1994.

[6]Wu Q F, Chen W F. DSMC method for high temperature chemical nonequilibrium rarefied flows[M]. Beijing: National University of Defense Technology Press, 1999. (in Chinses)吳其芬, 陳偉芳. 高溫稀薄氣體熱化學(xué)非平衡的DSMC方法[M]. 北京: 國(guó)防科技大學(xué)出版社, 1999.

[7]Alexander F J, Garcia A L. The direct simulation Monte Carlo method[J]. Computers in Physics, 1997, 11(6):588-593.

[8]Fan J, Shen Q. The Monte Carlo direct simulation of the hypersonic nonequilibrium flow past a circular cylinder in transition flow regimes[J]. Acta Aerodynamica Sinica, 1995, 13(4): 405-413. (in Chinese)樊菁, 沈青. 過(guò)渡領(lǐng)域高超聲速圓柱繞流的直接模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 1995, 13(4): 405-413.

[9]Li Z H, Wu Z Y. DSMC simulation on rarefied aerodynamics of Apollo-CM[J]. Acta Aerodynamica Sinica, 1996, 14(2): 230-233. (in Chinese)李志輝, 吳振宇. 阿波羅指令艙稀薄氣體動(dòng)力學(xué)特征的蒙特卡羅數(shù)值模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 1996, 14(2):230-233.

[10]Bird G A. Nonequilibrium radiation during re-entry at 10 km/s[R]. AIAA 87-1543.

[11]Carlson A B, Hassan H A. Direct simulation of re-entry flows with ionization[J]. Journal of Thermophysics & Heat Transfer, 2015, 6(3):400-404.

[12]Boyd I D. Modeling of associative ionization reactions in hypersonic rarefied flows[R]. Center for Turbulence Research, Annual Research Briefs 2007.

[13]Boyd I D. Modeling of plasma formulation in rarefied hypersonic entry flows[R]. AIAA 2007-206.

[14]Fan J, Zhang Y H, Jiang J Z. Monte Carlo modeling of electron density in hypersonic rarefied gas flows[C]//29thInternational Symposium on Rarefied Gas Dynamics, 2014.

[15]Shevyrin A A, Vashchenkov P V, Bonder Y A, et al. Validation of DSMC results for chemically nonequilibrium air flows against measurements of the electron number density in RAM-C II flight experiment[C]//29thInternational Symposium on Rarefied Gas Dynamics, 2014.

[16]Morsa L, Zuppardi G. Influence of ionization on the Gupta and on the Park chemical models[C]//29thInternational Symposium on Rarefied Gas Dynamics, 2014.

[17]Gupta R N, Li K P, Thompson R A, et al. Calculations and curve fits of thermodynamic and transport properties for equilibrium air to 30000K[R]. NASA RP-1260, 1991.

[18]Boyd I D. Modeling of associative ionization reactions in hypersonic rarefied flows[J]. Physics of Fluids, 2007, 19(9):1828-1837.

[19]Liang J. DSMC method on hybrid of cartesian and unstructured grid and its application[D]. Beihang University, 2014. (in Chinese).梁杰. 直角/非結(jié)構(gòu)混合網(wǎng)格的DSMC方法及應(yīng)用研究[D]. 北京航空航天大學(xué)博士畢業(yè)論文, 2014.

[20]Boyd I D. Conservative species weight scheme for Direct Simulation Monte Carlo method[J]. Journal of Thermophysics & Heat Transfer, 2015, 10(4):579-587.

[21]Fan J, Zhang Y H, Jiang J Z. An algorithm for trace species in hypersonic rarefied gas flows[C]//CSTAM 2012-B03-0417, 2012.

Three dimensional DSMC simulation and validation of rarefied air ionization process for sphere-cone blunt body reentry

Fang Ming1, Li Zhihui1,2,*, Li Zhonghua2, Tian Ying2

(1.BeihangUniversity,NationalLaboratoryforComputationalFluidDynamics,Beijing100191,China；2.HypervelocityAerodynamicsInstituteofChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China)

Based on the assumption that the electron mass is increased by three orders of magnitude and the ion mass is adjusted correspondingly, the traditional DSMC method for chemical reactions is extended for rarefied air ionization. Single temperature model is used for full reactions, and the reaction rate coefficients involving electron are modified to guarantee the real chemical reaction rate. Combined with rectangle/unstructured hybrid grid, three dimensional rarefied air ionization DSMC program has been developed for real complex spacecraft configurations based on collision cell self-adaption and MPI parallel environment. Numerical test indicates that, the electron distribution obtained for RAM-C II case agrees well with flight experiment data, and the electron density contours completely match with those in references. Compared with the DSMC method for rarefied gas without ionization reaction, the developed model and program hardly bring significant increment of computing cost, it can be directly applied to rarefied air ionization simulation for the three dimensional complex geometries at extremely high reentry speed and provide technical support for spacecraft design. The computational results show that the electron number density during extremely high speed reentry process in traditional rarefied flow regime can be high enough to introduce communication blackout and some of attention should be paid in the communication design stage.

spacecraft reentry; rarefied air; chemical reaction; ionization effect; DSMC method

0258-1825(2017)01-0039-07

2015-07-14;

2016-01-31

973計(jì)劃項(xiàng)目(2014CB744100)；國(guó)家自然基金項(xiàng)目(11325212、11602288)；中國(guó)博士后科學(xué)基金項(xiàng)目(2014M560870)

方明(1983-),男,湖北紅安人,出站博士后,副研究員,研究方向:稀薄氣體動(dòng)力學(xué). E-mail:fangm@pku.edu.cn

李志輝(1968-),)男,出站博士后,研究員,研究方向:稀薄氣體動(dòng)力學(xué). E-mail:zhli0097@x263.net

方明, 李志輝, 李中華, 等. 球錐鈍頭體再入稀薄氣體電離過(guò)程三維DSMC模擬與驗(yàn)證[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2017, 35(1): 39-45.

10.7638/kqdlxxb-2015.0086 Fang M, Li Z H, Li Z H, et al. Three dimensional DSMC simulation and validation of rarefied air ionization process for sphere-cone blunt body reentry[J]. Acta Aerodynamica Sinica, 2017, 35(1): 39-45.

V211.3

A doi: 10.7638/kqdlxxb-2015.0086