苗文博，黃 飛，程曉麗，俞繼軍

（中國(guó)航天空氣動(dòng)力技術(shù)研究院，北京 100074）

文章編號(hào)：1001?246X（2015）01?0027?06

再入飛行器等離子體預(yù)測(cè)與氣體組分相關(guān)性

苗文博，黃 飛?，程曉麗，俞繼軍

（中國(guó)航天空氣動(dòng)力技術(shù)研究院，北京 100074）

對(duì)一類(lèi)大鈍頭再入飛行器，理論分析電子密度預(yù)測(cè)與氣體組分的相關(guān)性，使用經(jīng)驗(yàn)證的數(shù)值模擬方法研究化學(xué)反應(yīng)模型對(duì)等離子體預(yù)測(cè)的影響.研究發(fā)現(xiàn)：馬赫數(shù)是影響等離子體預(yù)測(cè)與氣體組分相關(guān)性的重要參數(shù)，馬赫數(shù)越高，不同氣體組分模型所得電子密度差異越大.氣體組分模型對(duì)等離子體預(yù)測(cè)的影響在駐點(diǎn)強(qiáng)壓縮區(qū)域和身部位置基本一致.對(duì)于大鈍頭再入飛行器，高度H＝60 km，馬赫數(shù)大于23時(shí)應(yīng)該采用11組分化學(xué)反應(yīng)模型.

再入飛行器；等離子體；氣體組分；化學(xué)反應(yīng)模型；數(shù)值模擬

0 引言

當(dāng)飛行器以高超聲速再入大氣層時(shí)，激波后溫度快速升高，達(dá)6 000 K～10 000 K，甚至10 000 K以上，氣體大量離解并伴隨電離反應(yīng)進(jìn)行.飛行器周?chē)a(chǎn)生高密度、強(qiáng)電離的等離子層，形成“等離子體鞘套”.這種“等離子體鞘套”對(duì)通信天線的電磁波信號(hào)形成屏蔽作用，進(jìn)而產(chǎn)生“黑障”.“黑障”問(wèn)題最先在航天器再入過(guò)程中發(fā)現(xiàn)，針對(duì)高超聲速再入飛行器在進(jìn)入地球大氣時(shí)遇到的等離子體屏蔽問(wèn)題，NASA基于Apollo探月工程進(jìn)行了大量的理論分析、地面試驗(yàn)以及飛行試驗(yàn)研究，我國(guó)的探月和載人探月研究也會(huì)遇到類(lèi)似的等離子體屏蔽問(wèn)題.

等離子體由激波后氣體大量電離產(chǎn)生，實(shí)際飛行中，氣體電離程度受來(lái)流條件、壁面燒蝕引射產(chǎn)物等因素影響.Candler［1］、Ghislain［2］等人研究了不同飛行狀態(tài)下等離子體的分布規(guī)律，研究發(fā)現(xiàn)馬赫數(shù)越高，高度越低，電子密度越高.Scalabrin［3］針對(duì)高馬赫數(shù)狀態(tài)（Ma＞25）對(duì)比分析了不同的化學(xué)反應(yīng)模型對(duì)電子密度預(yù)測(cè)影響，發(fā)現(xiàn)11組分反應(yīng)動(dòng)力學(xué)模型所得電子密度高于7組分，與飛行試驗(yàn)結(jié)果更接近.程曉麗［4］在地面電弧風(fēng)洞實(shí)驗(yàn)研究了材料燒蝕對(duì)尾跡電子密度的影響，燒蝕將小幅降低尾跡電子密度.董維中［5］分析了再入體表面硅基材料燒蝕對(duì)電子密度的影響，燒蝕對(duì)頭部區(qū)域電子密度影響微弱，將降低身部及尾跡區(qū)電子密度.這些影響因素中，氣體組分及化學(xué)反應(yīng)模型的選擇對(duì)電子密度預(yù)測(cè)的影響最顯著.高鐵鎖研究發(fā)現(xiàn)當(dāng)飛行速度極高，激波后氣體溫度達(dá)10 000 K以上時(shí)，除了需要考慮NO＋離子還要考慮O＋離子對(duì)等離子體分布的影響［6］.Park［7］研究認(rèn)為，11組分化學(xué)動(dòng)力學(xué)模型可以更準(zhǔn)確地模擬流場(chǎng)的化學(xué)反應(yīng)過(guò)程，但是其同時(shí)帶來(lái)更大的計(jì)算成本.如何選擇盡可能少的氣體組分及相應(yīng)的化學(xué)動(dòng)力學(xué)模型，使得計(jì)算結(jié)果更準(zhǔn)確，計(jì)算成本更低？本文基于平衡流假設(shè)和Gibbs自由能原理分析氣體組分與電子密度預(yù)測(cè)的相關(guān)性，并采用經(jīng)驗(yàn)證的數(shù)值模擬方法，分析了一類(lèi)大鈍頭飛行器再入過(guò)程中，氣體組分及化學(xué)反應(yīng)動(dòng)力學(xué)模型對(duì)等離子體分布的影響規(guī)律和機(jī)理，獲得等離子體預(yù)測(cè)與氣體組分相關(guān)性的認(rèn)識(shí).

1 控制方程

求解的控制方程為三維化學(xué)反應(yīng)完全N?S方程［8］，并基于以下假定：流動(dòng)滿(mǎn)足化學(xué)非平衡及熱力學(xué)平衡；忽略輻射以及徹體力的影響；流動(dòng)質(zhì)量擴(kuò)散采用雙組元?dú)怏w模型假設(shè).三維化學(xué)反應(yīng)完全N?S方程如下，其中連續(xù)性方程為各組元連續(xù)性方程，同時(shí)包含化學(xué)反應(yīng)源項(xiàng).

其中

其中，ρi＝［ρ1，ρ2，…，ρns］T，ns為氣體組分的個(gè)數(shù)，ρ＝∑ρi為氣體總密度，u，v，w分別為x，y，z方向速度，P為壓強(qiáng)，e為單位質(zhì)量的總能量．Dim為氣體組分i的等效擴(kuò)散系數(shù)，為氣體組分i的化學(xué)反應(yīng)生成項(xiàng).

2 數(shù)值格式

選用具有較高流場(chǎng)分辨率的AUSM＋格式對(duì)控制方程進(jìn)行離散.AUSM＋格式［9］的主要思想認(rèn)為流場(chǎng)在傳播中存在對(duì)流影響與聲波影響，為分別考慮兩個(gè)過(guò)程，將無(wú)粘項(xiàng)分為對(duì)流項(xiàng)和壓力項(xiàng)進(jìn)行處理，基于馬赫數(shù)對(duì)兩者分別進(jìn)行特征分裂，是適用于化學(xué)反應(yīng)流求解的一類(lèi)高精度格式.下面給出其通量表達(dá)式.僅以項(xiàng)為例說(shuō)明.

其中，Ψ＝（ci，u，v，w，H），g＝（ci，nx，ny，nz，0），其中下標(biāo)i＋1表征界面參數(shù)，L和R分別表征左右臨近界面參數(shù).

3 數(shù)學(xué)物理模型

進(jìn)行流場(chǎng)以及電子密度場(chǎng)求解時(shí)分別使用了7組分7反應(yīng)模型和11組分20反應(yīng)模型［10］.7組分模型考慮了N，O，N2，O2，NO，NO＋，e－7個(gè)組分，只考慮了NO分子的電離過(guò)程，11組分模型考慮了N，O，N2，O2，NO，NO＋，，，N＋，O＋，e－11個(gè)組分，同時(shí)考慮了NO，N2，O2，N，O這5種氣體組分的電離.

選擇雙子星再入飛行器進(jìn)行電子密度預(yù)測(cè)與氣體組分的相關(guān)性分析，其具有和探月飛行器類(lèi)似的倒鐘型布局特征.圖1給出飛行器基準(zhǔn)外形幾何尺寸示意圖.

4 算法驗(yàn)證

選擇RAM?CII飛行試驗(yàn)結(jié)果對(duì)電子密度預(yù)測(cè)方法進(jìn)行驗(yàn)證［11］.RAM?CII飛行器是一個(gè)頭部半徑為0.152 4 m，半錐角為9°，長(zhǎng)度為1.295 m的球錐體.選用非催化壁壁面條件，計(jì)算高度 H＝71 km，Ma＝25.9的飛行狀態(tài)，壁面溫度取Tw＝1 500 K，分別進(jìn)行了7組分和11組分化學(xué)反應(yīng)模型.此時(shí)，飛行器繞流流場(chǎng)中氣體充分電離.圖2給出該狀態(tài)下沿軸向的電子密度峰值分布對(duì)比.預(yù)測(cè)結(jié)果與實(shí)際飛行試驗(yàn)測(cè)量結(jié)果及分布趨勢(shì)吻合較好.并且可以看到，對(duì)于該狀態(tài)，11組分預(yù)測(cè)結(jié)果比7組分預(yù)測(cè)結(jié)果大約1個(gè)量級(jí).11組分模型預(yù)測(cè)結(jié)果高估峰值電子密度，而7組分模型預(yù)測(cè)結(jié)果低估峰值電子密度.

圖1 飛行器幾何尺寸示意圖Fig.1 Outline of the vehicle

圖2 峰值電子密度沿流向?qū)Ρ菷ig.2 Peak electron density along x direction

5 結(jié)果及分析

對(duì)于鈍頭再入飛行器，頭部駐點(diǎn)區(qū)域受強(qiáng)激波壓縮，化學(xué)反應(yīng)非常劇烈，流動(dòng)可以近似看作平衡流［12］.基于Gibbs自由能原理求解可以獲得激波后氣體平衡組分濃度分布［13］.通過(guò)對(duì)比氣體組分中離子的含量，可以確定影響等離子體分布的主要組分，從而為數(shù)值模擬的氣體組分選擇提供參考.此時(shí)，激波后氣體壓力由正激波關(guān)系式給出，激波后氣體總焓依據(jù)總能守恒原理由來(lái)流總焓給出.

表1給出高度H＝60 km多個(gè)馬赫數(shù)狀態(tài)下的平衡狀態(tài)正激波后氣體電離組分摩爾分?jǐn)?shù)以及N＋，O＋，和離子之和對(duì)電子密度所做貢獻(xiàn)與NO＋所做貢獻(xiàn)之比.此處取N＋，O＋，和離子對(duì)電子密度的貢獻(xiàn)與NO＋離子貢獻(xiàn)之比?為參考量.由表1可以看到，對(duì)于高度H＝60 km，隨著馬赫數(shù)升高激波后氣體電離程度快速增加.在馬赫數(shù)M＝18時(shí)，?僅0.037 8，在馬赫數(shù)M＝25時(shí)，?已經(jīng)超過(guò)1，預(yù)示此時(shí)N＋，O＋，和離子對(duì)電子密度的貢獻(xiàn)已超過(guò)NO＋離子，7組分氣體模型已經(jīng)失效.

表1 H＝60 km處，平衡假設(shè)下激波后離子摩爾分?jǐn)?shù)Table 1 M olar fraction of ions after the shock at H＝60 km

依據(jù)Gibbs自由能原理和駐點(diǎn)區(qū)域平衡流假設(shè)僅給出頭部強(qiáng)壓縮區(qū)域的氣體組分與等離子體預(yù)測(cè)隨來(lái)流條件變化的規(guī)律，對(duì)于身部以及尾跡區(qū)域同樣的規(guī)律是否依然存在？由于身部存在明顯的膨脹以及大分離流動(dòng)，氣體平衡假設(shè)很難滿(mǎn)足，無(wú)法使用平衡流假設(shè)進(jìn)行分析.以往研究表明，等離子體主要由頭部強(qiáng)壓縮區(qū)域產(chǎn)生，隨著流動(dòng)逐漸向下游發(fā)展，上游的流動(dòng)特征將對(duì)下游產(chǎn)生一定的影響.因此設(shè)置以下算例使用數(shù)值模擬的方法進(jìn)行分析，高度H＝60 km，飛行攻角20°，飛行速度分別為10 km·s－1，7.0 km·s－1，5.5 km·s－1，對(duì)應(yīng)的馬赫數(shù)分別為32.66，22.86和17.96.對(duì)這三個(gè)狀態(tài)分別使用7組分模型和11組分模型進(jìn)行了數(shù)值模擬分析.

圖3給出3個(gè)狀態(tài)不同化學(xué)反應(yīng)模型所得電子密度對(duì)比，可以看到在Ma＝17.96時(shí)，11組分模型所得電子密度略高于7組分模型，相差在5%以?xún)?nèi)，隨著馬赫數(shù)升高，11組分模型所得電子密度與7組分模型差異逐漸明顯，Ma＝22.86的時(shí)候不同化學(xué)反應(yīng)模型所得電子密度相差30%左右，Ma＝32.66時(shí)不同化學(xué)反應(yīng)模型差異1－2個(gè)量級(jí).駐點(diǎn)區(qū)域數(shù)值預(yù)測(cè)所示規(guī)律與基于Gibbs自由能原理分析所得基本一致，在Ma＝23時(shí)，不同氣體組分所得電子密度差異已大于30%.這種差異在頭部強(qiáng)壓縮區(qū)域與身部大面積區(qū)域同時(shí)存在，不同位置的氣體組分對(duì)電子密度預(yù)測(cè)的影響規(guī)律一致，基于平衡流假設(shè)和Gibbs自由能原理可以獲得電子密度預(yù)測(cè)與氣體組分相關(guān)性較合理的認(rèn)識(shí).

圖3 兩種化學(xué)反應(yīng)模型迎風(fēng)子午面峰值電子密度x方向分布Fig.3 Peak electron density along windward generatrix in chemical reaction models

對(duì)比發(fā)現(xiàn)，數(shù)值預(yù)測(cè)所得差異略小于理論分析所得，這主要是因?yàn)槭褂?1組分化學(xué)反應(yīng)模型時(shí)，激波后氣體內(nèi)能除了用于激發(fā)NO＋電離之外還要用于激發(fā)N＋，O＋，和離子電離，而使用7組分化學(xué)反應(yīng)模型時(shí)更多的能量用于NO＋電離，NO＋離子貢獻(xiàn)的電子密度要高于11組分反應(yīng)時(shí)NO＋離子貢獻(xiàn)的電子密度.

駐點(diǎn)區(qū)域是流動(dòng)化學(xué)反應(yīng)最強(qiáng)烈的位置，對(duì)馬赫數(shù)Ma＝17.96和Ma＝32.66兩個(gè)狀態(tài)進(jìn)一步分析不同化學(xué)反應(yīng)模型的駐點(diǎn)線流場(chǎng)參數(shù).圖4給出兩個(gè)狀態(tài)駐點(diǎn)線溫度分布，Ma＝32.66時(shí)，激波后氣體溫度在10 000 K左右，考慮N＋，O＋，和離子氣體電離過(guò)程使得更多的氣體內(nèi)能轉(zhuǎn)化為化學(xué)能，降低了激波后氣體溫度，同時(shí)使得激波層更靠近壁面.而對(duì)于Ma＝17.96，激波后氣體溫度在6 000 K左右，氣體電離程度減弱，即使考慮更多的氣體電離過(guò)程仍對(duì)流場(chǎng)能量分布狀態(tài)影響微弱.

圖4 駐點(diǎn)線溫度分布（左：Ma＝17.96，右：Ma＝32.66）Fig.4 Temperature along stagnation line（left：Ma＝17.96，right：Ma＝32.66）

圖5給出2個(gè)狀態(tài)駐點(diǎn)線離子數(shù)密度分布對(duì)比.對(duì)于Ma＝32.66狀態(tài)，受高溫高壓作用，氣體電離以原子電離為主，N＋，O＋，和離子電離明顯，特別是N＋和O＋離子.對(duì)于Ma＝17.96狀態(tài)，N＋，O＋，和離子相比NO＋的電離程度微弱的多.可見(jiàn)當(dāng)流場(chǎng)中N＋，O＋，和離子的電離度與NO＋相當(dāng)或者高于NO＋離子時(shí)，必須考慮11組分化學(xué)反應(yīng)模型.此處假定N＋，O＋，和離子的電離度與NO＋的電離度之比大于20%時(shí)，必須考慮11組分化學(xué)反應(yīng)模型.由此可以認(rèn)為，對(duì)于大鈍頭再入飛行器，當(dāng)H＝60 km，馬赫數(shù)高于23的時(shí)候，應(yīng)該選擇11組分化學(xué)反應(yīng)模型進(jìn)行電子密度預(yù)測(cè).

圖5 駐點(diǎn)線離子數(shù)密度分布（左：Ma＝17.96，右：Ma＝32.66）Fig.5 Number density of ions along stagnation line（left：Ma＝17.96，right：Ma＝32.66）

6 小結(jié)

高超聲速再入飛行器再入過(guò)程往往經(jīng)歷“黑障”，要進(jìn)行飛行器通信中斷邊界分析需要獲得準(zhǔn)確的電子密度分布，本文使用數(shù)值模擬和理論分析手段研究了大鈍頭再入飛行器再入過(guò)程時(shí)氣體組分對(duì)電子密度預(yù)測(cè)的影響規(guī)律和機(jī)理.獲得以下結(jié)論：對(duì)于高馬赫數(shù)再入飛行狀態(tài)，馬赫數(shù)是影響氣體組分與電子密度預(yù)測(cè)相關(guān)性的重要參數(shù)，馬赫數(shù)越高，不同化學(xué)反應(yīng)模型所得電子密度差異越大.氣體組分對(duì)等離子體預(yù)測(cè)的影響規(guī)律在駐點(diǎn)強(qiáng)壓縮區(qū)域和身部大面積區(qū)域一致，理論分析與數(shù)值模擬所得結(jié)論基本一致，可以通過(guò)理論分析快速選擇等離子體預(yù)測(cè)時(shí)的氣體組分模型.對(duì)于大鈍頭再入飛行器，高度H＝60 km，馬赫數(shù)Ma大于23時(shí)應(yīng)考慮11組分化學(xué)反應(yīng)模型.

［1］ Candler G V.The computational of weakly ionized flow in non?equilibrium［D］.California：Stanford University，1988.

［2］ Ghislain T.Numerical prediction ofweakly ionized high enthalpy flow in thermo?chemical non?equilibrium［R］.AIAA 2004－2462，2004.

［3］ Scalabrin L C，Boyd ID.Numerical simulation ofweakly ionized hypersonic flow for reentry configurations［R］.AIAA 2006－3773，2006.

［4］ Cheng X L，Dong Y H，Li T L.Computational and experimental studies of ablation effect on electronic characteristics in vehicle wake［J］.Acta Aeronautica et Astronautica Sinica，2007，28（4）：796－780.

［5］ DongW Z，Gao T S，Ding M S，Jiang T.Numerical analysis for the effect of silicon based material ablation on the flowfield around re?entry blunt body［J］.Acta Aerodynamica Sinica，2010，28（6）：708－714.

［6］ Gao TS，DongW Z，Ding M S.The effects of physicochemicalmodels on distribution of plasma in high temperature flow field ［J］.Acta Aerodynamica Sinica，2013，31（5）：541－545.

［7］ Park C，Griffith W，Nonequilibrium hypersonic aerothermodynamics［M］.New York：Wiley，1990.

［8］ 歐陽(yáng)水吾，高溫非平衡空氣繞流［M］.北京：國(guó)防工業(yè)出版社，2001.

［9］ Liou M S.A further developmentof the AUSM＋scheme towards robustand accurate solutions for all speeds［R］.AIAA Paper 2003－4116，2003.

［10］ Gnoffo P A，Gupta R N，Shinn J L.Conservation equations and physical models for hypersonic air flows in thermal andchemical non?equilibrium［R］.NASA TP－2867，1986.

［11］ Scbexnayder C J，Huber P W，Evans J S.Calculation of electron concentration for a blunt body at orbital speed and comparison with experimental data［R］.NASA TN?D－6294，1971.

［12］ Park C，Howe JT，Jaffe R L.Review of chemical?kinetic problems of future NASAmissions，I：Earth entries［J］.Journalof Thermo?physics and Heat Transfer，1993，7（3）：383－392.

［13］ Liu X Z.Numerical simulation of hypersonic thermochemical equilibrium flows［D］.Beijing：China Academy of Aerospace Aerodynamics，1994.

Plasma Prediction of Reentry Vehicle and Gas Com ponents

MIAOWenbo，HUANG Fei，CHENG Xiaoli，YU Jijun
（China Academy ofAerospace and Aerodynamics，P．O．Box 7201?16，Beijing 100074，China）

Hypersonic flows around a bluntbody are studied with numerical simulation and theoretical analysis.Correlations between plasma prediction and gas components are obtained.It shows that it is Mach number which affects the peak number density of electrons.Numerical simulation agreeswith theoretical analysis well.Differences between two gasmodels get greater as Mach number increases.The differences follow same trend at back taper with stagnation region.A 11?species chemicalmodel should be applied to increase accuracy when reentry capsule flight at height of60 km and Mach number is over 23.

reentry vehicle；plasma；gas components；chemical reaction model；numerical simulation

V211.3

A

2013－12－30；

2014－04－11

苗文博（1980－），男，博士，高級(jí)工程師，主要從事非平衡流研究，E?mail：tingles＠126.com

?通訊作者：黃飛（1982－），高級(jí)工程師，E?mail：huang05013＠163.com

Received date： 2013－12－30；Revised date： 2014－04－11