高鐵鎖，董維中，丁明松，江 濤

（中國(guó)空氣動(dòng)力研究與發(fā)展中心，四川 綿陽(yáng) 621000）

0 引 言

高超聲速飛行器在大氣層中以高超聲速飛行時(shí)，由于激波和粘性作用，波后氣體溫度超過(guò)10000K，這時(shí)高溫空氣發(fā)生振動(dòng)能量激發(fā)、離解和電離等化學(xué)反應(yīng)過(guò)程，成為部分電離氣體，并在飛行器周?chē)纬筛邷氐入x子體流場(chǎng)。由于飛行器周?chē)入x子體分布直接關(guān)系到電磁波的傳輸特性［1－6］，因此研究高超聲速飛行器流場(chǎng)中等離子體的產(chǎn)生機(jī)理和分布特性，對(duì)于高超聲速飛行器的無(wú)線電通信技術(shù)的發(fā)展具有重要意義。

飛行器周?chē)鲌?chǎng)等離子體分布與飛行高度、速度、姿態(tài)和飛行器外形尺寸、表面材料催化特性等多種因素有關(guān)。采用求解N－S方程的數(shù)值方法對(duì)其進(jìn)行模擬時(shí)，除了需要解決熱化學(xué)非平衡源項(xiàng)引起的剛性問(wèn)題外，還需要選取合適的化學(xué)反應(yīng)模型，因?yàn)橛?jì)算中考慮的組分越多計(jì)算越復(fù)雜，對(duì)計(jì)算資源的要求越高，而采用含組分少的簡(jiǎn)化模型有時(shí)不能全面反映高溫流動(dòng)中等離子體的產(chǎn)生機(jī)制，影響數(shù)值模擬的精準(zhǔn)度［7］。數(shù)值模擬高溫空氣電離流動(dòng)常用的有7組分和11組分化學(xué)反應(yīng)模型［8－10］，如得到廣泛應(yīng)用的Blottner的7組分反應(yīng)模型［8］和Dunn與Kang的11組分反應(yīng)模型［9］等。由于高溫電離流動(dòng)問(wèn)題的復(fù)雜性，目前還沒(méi)有明確和統(tǒng)一的化學(xué)反應(yīng)模型選取準(zhǔn)則。另外，準(zhǔn)確模擬高溫電離流動(dòng)還需要研究熱力學(xué)非平衡效應(yīng)及壁面催化和壁面溫度等邊界條件的影響［11－14］。高溫氣體熱力學(xué)非平衡涉及分子振動(dòng)能級(jí)之間和振動(dòng)與平動(dòng)能級(jí)之間的能量交換以及振動(dòng)－離解耦合效應(yīng)等復(fù)雜物理化學(xué)過(guò)程，對(duì)熱力學(xué)非平衡效應(yīng)的數(shù)值模擬一直是高溫氣體動(dòng)力學(xué)領(lǐng)域的前緣課題之一［7，11，13－14］。本文針對(duì)復(fù)雜升力體外形，采用數(shù)值模擬方法，研究不同反應(yīng)模型和非平衡氣體模型以及壁面催化和溫度條件對(duì)流場(chǎng)等離子體數(shù)值模擬結(jié)果的影響。

1 復(fù)雜外形戰(zhàn)術(shù)彈頭非平衡流場(chǎng)的計(jì)算方法

1.1 控制方程

控制方程是三維熱化學(xué)非平衡N－S方程，其無(wú)量綱化形式如下［13］：

其中Q＝（ρi，ρEV，ρ，ρu，ρv，ρw，ρE）T，W＝（wi，wV，0，0，0，0，0）T，ρi是組分密度，EV是分子組分的總振動(dòng)能量，wi和wV是化學(xué)非平衡源項(xiàng)和振動(dòng)非平衡能量源項(xiàng)。

1.2 熱化學(xué)模型

采用11組分（O2，N2，NO，O，N，NO＋，e，O＋2，，O＋，N＋）和7組分（O2，N2，NO，O，N，NO＋，e）化學(xué)反應(yīng)模型，具體反應(yīng)式及反應(yīng)速率常數(shù)見(jiàn)文獻(xiàn)［9，13］。對(duì)于正向和逆向反應(yīng)速率常數(shù)為kfi與kbi，第i個(gè)反應(yīng)表達(dá)式（2）對(duì)應(yīng)的第j個(gè)組分的化學(xué)反應(yīng)的生成源項(xiàng)見(jiàn)式（3）：

振動(dòng)非平衡能量源項(xiàng)為兩項(xiàng)之和，一項(xiàng)是平動(dòng)能與振動(dòng)能的交換項(xiàng)，另一項(xiàng)是化學(xué)反應(yīng)中分子離解引起的振動(dòng)能量的變化項(xiàng)，其表達(dá)式如下［13］：

1.3 壁面催化條件

圖1 計(jì)算與飛行測(cè)量電子數(shù)密度分布比較Fig.1 Comparison of computation with flight experiment

1.4 數(shù)值方法

采用全隱式LU－SGS數(shù)值方法對(duì)方程（1）求解［13－15］，無(wú)粘項(xiàng)的離散采用Yee的TVD格式［15－16］，粘性項(xiàng)的離散采用中心差分方法。

2 計(jì)算結(jié)果分析

2.1 RAM－C模型繞流等離子體分布的比較

RAM－C模型為球錐體，球頭半徑Rn＝0.1542m，半錐角θ＝9°，長(zhǎng)度L＝1.295m；飛行條件：U∞＝7650.0m／s，H＝61km、71km，壁面溫度TW＝1500K。

采用7組分化學(xué)反應(yīng)模型和兩溫度熱力學(xué)非平衡模型，考慮完全非催化（NCW）和完全催化（FCW）壁面條件，主要考核電子數(shù)密度分布情況。圖1給出了兩個(gè)飛行高度的電子密度分布與文獻(xiàn)［14］的飛行數(shù)據(jù)比較，其中圖1（a）和圖1（b）對(duì)應(yīng)流場(chǎng)剖面電子數(shù)密度峰值沿軸向分布?？梢钥闯?，電子密度計(jì)算與飛行測(cè)量數(shù)據(jù)符合較好，誤差只在3倍左右。飛行速度不變，隨著高度的降低，流動(dòng)越接近熱化學(xué)平衡態(tài)，完全非催化和完全催化條件的計(jì)算結(jié)果越接近。從圖1（c）可見(jiàn)，催化條件對(duì)X／Rn＝8.1處的電子數(shù)密度徑向分布具有較大影響。總體上看，完全催化壁條件的計(jì)算結(jié)果更接近飛行測(cè)量值。

2.2 化學(xué)反應(yīng)模型和非平衡模型對(duì)流場(chǎng)等離子體分布的影響

計(jì)算外形為升力體外形，計(jì)算高度60km，飛行馬赫數(shù)為20，攻角12°。認(rèn)為壁面組分滿(mǎn)足完全非催化壁面條件，壁面溫度取600K。P1和P2分別表示靠近升力體頭部和底部的兩個(gè)剖面位置。下面主要給出升力體對(duì)稱(chēng)面上迎風(fēng)區(qū)的計(jì)算結(jié)果。

圖2給出了采用7組分與11組分反應(yīng)模型計(jì)算的壁面離子和電子數(shù)密度沿軸向分布的比較?？梢钥闯?，11組分化學(xué)反應(yīng)模型的電子數(shù)密度計(jì)算值高于7組分的計(jì)算值，二者在升力體球頭區(qū)域差別較大，最大相差一個(gè)量級(jí)，隨著流動(dòng)向下游發(fā)展，其差別越來(lái)越小。頭部流場(chǎng)溫度最高達(dá)10000K以上（見(jiàn)圖3a），此時(shí)除了考慮NO＋離子，化學(xué)反應(yīng)模型中還應(yīng)考慮O＋的作用。NO＋主要由O與N原子的締合電離反應(yīng)機(jī)制產(chǎn)生，O＋則主要由O原子與電子及其它組分的碰撞反應(yīng)機(jī)制產(chǎn)生。由于頭部流動(dòng)區(qū)域溫度較高，使得O原子電離起主導(dǎo)作用，O電離生成O＋對(duì)電子數(shù)密度的貢獻(xiàn)超過(guò)NO電離生成NO＋的貢獻(xiàn)，而在下游區(qū)域流場(chǎng)溫度逐漸降低，NO的化學(xué)電離逐漸起主導(dǎo)作用，7組分和11組分兩種反應(yīng)模型的計(jì)算結(jié)果也逐漸趨于一致。高溫流動(dòng)中化學(xué)電離過(guò)程除了與溫度密切相關(guān)外，還與流動(dòng)壓力與非平衡馳豫過(guò)程有關(guān)，壓力有抑制電離和離解的作用。

圖3和圖4給出P1和P2剖面的溫度、離子組分質(zhì)量分?jǐn)?shù)和電子數(shù)密度分布。其中圖3（a）和圖4（a）給出完全氣體模型（PG）和非平衡氣體模型的流場(chǎng)溫度的比較，可以看出，完全氣體模型的流場(chǎng)溫度和激波脫體距離的計(jì)算值明顯高于考慮真實(shí)氣體效應(yīng)的非平衡氣體模型的計(jì)算值。這是因?yàn)榉肿与x解和電離等吸熱反應(yīng)過(guò)程使得氣體流場(chǎng)溫度降低，密度增加，而化學(xué)過(guò)程對(duì)流動(dòng)壓力影響很小，故使得激波脫體距離減小。從圖3和圖4還可以看出，在升力體的頭部流場(chǎng)區(qū)域（P1剖面），由于流場(chǎng)溫度很高，兩種化學(xué)反應(yīng)模型計(jì)算的剖面電子數(shù)密度差別很大，此時(shí)電子主要來(lái)源于O電離生成O＋的化學(xué)過(guò)程，然后以次是氣體分子生成NO＋、O＋2和N＋2的化學(xué)過(guò)程。在升力體的后部（P2剖面），由于流場(chǎng)溫度降低，二者的剖面電子數(shù)密度分布差別變得很小，此時(shí)O＋的作用較頭部明顯減弱。

圖5給出P1和P2流場(chǎng)剖面的溫度與電子數(shù)密度分布云圖?？梢?jiàn)，迎風(fēng)面的等離子體厚度明顯小于背風(fēng)面的等離子體厚度，隨著流動(dòng)從升力體頭部（P1剖面）發(fā)展到后部（P2剖面），剖面溫度最大值從12144K降低到4118K，剖面最大電子數(shù)密度從5.98×1014／cm3衰減到2.49×1010／cm3，但等離子體厚度明顯增加，并且在升力體對(duì)稱(chēng)面附近等離子體相對(duì)較厚，沿展向逐漸變薄。

圖2 壁面離子及電子數(shù)密度分布Fig.2 Distribution of ions andelectron number densities on the wall

圖3 P1剖面流場(chǎng)參數(shù)分布Fig.3 Distribution of flow parameter of P1profile

圖4 P2剖面流場(chǎng)參數(shù)分布Fig.4 Distribution of flow parameters of P2profile

圖5 剖面溫度與電子數(shù)密度云圖（Ns＝11）Fig.5 Temperature andelectron number densities contours

圖6給出了采用化學(xué)非平衡（CNeq）和熱化學(xué)非平衡模型（TCNeq）計(jì)算的頭部軸線溫度和電子數(shù)密度分布，可見(jiàn)在壁面附近區(qū)域，兩種模型的計(jì)算結(jié)果一致，但在離開(kāi)壁面較遠(yuǎn)的流動(dòng)區(qū)域二者出現(xiàn)差異，化學(xué)非平衡模型計(jì)算的脫體距離較小，電子數(shù)密度衰減較快。熱化學(xué)非平衡模型中考慮了分子振動(dòng)內(nèi)能模式的非平衡馳豫過(guò)程和振動(dòng)－離解耦合效應(yīng)，目前此模型還在發(fā)展之中［7，14］，對(duì)其可靠性驗(yàn)證還存在很大困難。

圖6 非平衡模型對(duì)溫度和電子數(shù)密度分布的影響Fig.6 Effects of nonequilibrimmodels on distribution of temperature and electron number densities

2.3 壁面條件對(duì)流場(chǎng)等離子體分布的影響

圖7（a）給出了不同壁面催化條件對(duì)P2剖面電子數(shù)密度分布的影響。這里沒(méi)有考慮NO＋等離子在壁面的催化復(fù)合反應(yīng)，只是考慮了O和N原子在壁面的有限催化反應(yīng)，并認(rèn)為它們的壁面催化速率常數(shù)為5m／s。從圖7（a）可以看出，有限催化（KW＝5m／s）條件下，電子數(shù)密度計(jì)算值接近完全非催化壁面的計(jì)算結(jié)果。

圖7 壁面條件對(duì)電子數(shù)密度分布的影響Fig.7 Effects of wall condition on distribution of electron number densities

圖7（b）給出了不同壁面溫度條件對(duì)頭部軸線流場(chǎng)電子數(shù)密度分布的影響。在完全非催化條件下，壁面溫度對(duì)電子數(shù)密度的影響主要集中在壁面附近很小的區(qū)域，貼近壁面的電子數(shù)密度以及電子數(shù)密度下降的梯度隨壁面溫度增加而減小。這里給出的結(jié)果是在同樣計(jì)算網(wǎng)格下獲得的，改變網(wǎng)格的數(shù)量和分布可能對(duì)計(jì)算結(jié)果產(chǎn)生不同程度的影響。

3 結(jié)束語(yǔ)

建立了復(fù)雜外形飛行器三維流場(chǎng)等離子體數(shù)值計(jì)算手段，實(shí)現(xiàn)了多塊網(wǎng)格并行運(yùn)算功能，具備模擬復(fù)雜外形高溫流場(chǎng)電離效應(yīng)、熱化學(xué)非平衡效應(yīng)和壁面催化效應(yīng)的能力。RAM－C模型繞流等離子體分布數(shù)值計(jì)算與飛行測(cè)量結(jié)果具有較好的一致性，表明了模型和數(shù)值方法的可行性。通過(guò)研究熱化學(xué)模型與壁面條件對(duì)升力體流場(chǎng)等離子體分布的影響，可得到以下兩點(diǎn)認(rèn)識(shí)：

（1）準(zhǔn)確數(shù)值模擬高溫流場(chǎng)等離子體分布，必須選擇合適的化學(xué)反應(yīng)模型?；瘜W(xué)反應(yīng)模型的選取主要取決于流場(chǎng)等離子體的溫度和壓力，另外還與流動(dòng)的非平衡松馳過(guò)程有關(guān)。對(duì)于本文研究的升力體外形和飛行狀態(tài)，頭部流場(chǎng)溫度最高達(dá)10000K以上，此時(shí)除了考慮NO＋離子，O＋離子對(duì)電子數(shù)密度的也有重要貢獻(xiàn)，化學(xué)反應(yīng)模型中應(yīng)包括O＋。采用不同非平衡氣體模型（熱化學(xué)非平衡或化學(xué)非平衡氣體模型）對(duì)高溫流場(chǎng)等離子體分布計(jì)算結(jié)果也會(huì)產(chǎn)生一定影響。

（2）壁面條件特別是壁面催化條件對(duì)流場(chǎng)等離子體分布具有較大影響。為了準(zhǔn)確數(shù)值模擬升力體飛行器繞流等離子體分布，需要結(jié)合飛行和地面試驗(yàn)測(cè)量結(jié)果，獲得可靠的材料催化特性和壁面溫度數(shù)據(jù)，在此基礎(chǔ)上對(duì)其周?chē)入x子體環(huán)境進(jìn)行精細(xì)模擬。

［1］MUYLAERT J，CIPOLLINI F，et al.In－flight research on real gas effects using the ESA EXPERT vehicles［R］.AIAA 2003－6981，2003.

［2］MORABITO D，KORNFELD R，et al.The mars phoenix communications brownout during entry into the martian atmosphere［R］.IPN Progress Report，2009：42－179.

［3］HARTUNIAN R A，STEWART G E，et al.Implications and mitigation of radio frequency blackout during reentry of reusable launch vehicles［R］.AIAA 2007－6633，2007.

［4］鐘育民，諶明，盧滿(mǎn)宏，等.無(wú)線電波在再入等離子體中傳輸?shù)乃p模型及仿真驗(yàn)證［J］.遙測(cè)遙控，2010，31（2）：1－6.

［5］MATHER D E，PASQUAL J M，SILLENCE J P.Radio frequency（RF）blackout during hypersonic reentry［R］.AIAA 2005－3443.2005.

［6］STARKEY R P.Electromagnetic wave／magnetoactive plasma sheath interaction for hypersonic vehicle telemetry blackout analysis［R］.AIAA 2003－4167，2003.

［7］SCALABRIN L C，BOYD I D.Numerical simulation of weakly ionized hypersonic flow for reentry configurations［R］.AIAA 2006－3773，2006.

［8］BLOTTNER F G.Viscous shock layer at the stagnation point with nonequilibrium air chemistry［J］.AIAA Journal，1969，7（12）：2281－2288.

［9］DUNN M G，KANG S W.Theoretical and experimental studies of reentry plasmas［R］.NASA－CR－2232，1973.

［10］GUPTA R N，YOS J M，et al.A review of reaction rates and thermodynamic and transport properties for an 11－species air model for chemical and thermal nonequilibrium calculations to 30000K［R］.NASA－RP－1232，1990.

［11］LIN T C，SPROUL L K，et al.Reentry plasma effects on electromagnetic wave propagation［R］.AIAA 95－1942.

［12］NUSCA M J，COOPER G R.Computational simulation of EM attenuation by plasmas formed in hypervelocity atmosphereic flight［R］.AIAA－98－1573，1998.

［13］董維中.熱化學(xué)非平衡效應(yīng)對(duì)高超聲速流動(dòng)影響的數(shù)值計(jì)算與分析［D］.［博士學(xué)位論文］.北京航空航天大學(xué)，1996.

［14］CANDLER G V，MacCORMACK R W.The computation of hypersonic ionized flows in chemical and thermal nonequilibium［R］.AIAA－88－051，1988.

［15］高鐵鎖，李椿萱，董維中，等.高超聲速電離繞流的數(shù)值模擬［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2002，20（2）：184－191.

［16］YEE H C，KLOPFER G H，et al.High－resolution shock－capturing schemes for inviscid and viscous hypersonic flows［R］.NASA－TM－100097，1988.