李江挺 郭立新 金莎莎 方全杰

(西安電子科技大學(xué)理學(xué)院,陜西西安710071)

1.引 言

飛行器在空間高速飛行時,無線電信號會呈現(xiàn)出一定程度的衰減。如果飛行器的速度繼續(xù)提高達到十幾馬赫或者幾十馬赫時,飛行器與大氣強烈作用,在頭部形成弓形脫體激波,波后氣體溫度、壓強急劇升高,使大氣離解、電離,在飛行器周圍形成等離子體鞘套[1-2]。等離子體鞘套電磁參數(shù)是流場參數(shù)的函數(shù),其分布特性與高超聲速流場特性緊密相關(guān),同時與飛行器外形、飛行馬赫數(shù)及飛行高度有關(guān)。等離子體鞘套的存在,使電磁波產(chǎn)生反射、折射及散射[3],同時吸收電磁波能量,使地面站與飛行器間通信受擾,如果中斷的時間太長,導(dǎo)致目標可能消失。為了有效地識別和跟蹤飛行器,必須研究飛行器表面等離子體鞘套對電磁波傳播的影響。

隨著計算機技術(shù)和計算流體力學(xué)的飛速發(fā)展,數(shù)值模擬已在工程實際中發(fā)揮越來越大的作用,而差分格式則是計算流體力學(xué)中最為核心的因素之一。本文采用雜交通量分裂格式中的AUSMPW+格式。這種格式是在迎風分裂格式(AUSM)的基礎(chǔ)上引入壓力權(quán)函數(shù)修正等技術(shù),構(gòu)造簡單,無矩陣運算,捕捉激波能力強且穩(wěn)定性好[4-5]。對飛行器表面流場的求解,采用了熱化學(xué)非平衡流Navier-Stokes控制方程組[6]、二溫度模型以及11組元的姜鄧恩(Dunn-Kang)空氣化學(xué)模型計算不同飛行速度下的流場參數(shù)。計算結(jié)果與通量分裂格式(Roe格式)、直接模擬蒙特卡羅方法(DSMC方法)[7]以及無線電信號衰減測量(RAM-C)實驗結(jié)果[8]進行了比較,相對其他兩種方法本文計算結(jié)果與實驗數(shù)據(jù)符合較好。同時,在研究電磁波在非均勻等離子體中傳播特性時,本文在傳統(tǒng)WKB方法的基礎(chǔ)上,應(yīng)用差分傳輸矩陣方法[9-10],將WKB解由一階解提高至二階解,分析電磁波在等離子鞘套中的衰減特性。計算時,首先通過流場模擬得到不同馬赫數(shù)下的電子密度分布,再根據(jù)改進的WKB方法計算了等離子鞘套的透射系數(shù)。同時在此基礎(chǔ)上,分析了電磁波在鞘套中傳輸?shù)南囝l特性,并給出了群時延與飛行器飛行速度及入射電波頻率的關(guān)系。

目前國內(nèi)外研究飛行器包覆流場的電磁特性[11-12]時,通常將電子密度函數(shù)近似為指數(shù)分布、拋物線分布或者在其基礎(chǔ)上做出一些修正,這種近似與實際電子密度剖面存在較大差異,而模擬實際流場下的電波傳播計算則鮮見報道。本文通過建立的飛行器模型模擬等離子體流場,求解空氣流動控制方程組,得出電子密度剖面,其更接近流場中的真實情況,其結(jié)果在飛行器載入通信、制導(dǎo)與測控等領(lǐng)域,具有重要的理論意義及應(yīng)用價值。

2.等離子鞘套模型

2.1 AUSMPW+格式的Navier-Stokes方程組

模擬高速飛行器表面繞流流場采用的軸對稱熱化學(xué)非平衡流動的無量綱化控制方程組為

式中:Q為守恒變量組成的矢量;F、G分別為x、r方向的對流通量矢量;FV、HV分別為x、r方向的粘性項矢量;H、H V為有粘部分和無粘部分的源項矢量;W為化學(xué)反應(yīng)和振動能量源項矢量Re為流動的雷諾數(shù)。求解控制方程(1),將其離散為

其中

Φn為無粘數(shù)值通量項,Ωi,j為粘性項與源項之和。為了提高激波捕捉效率和計算精度,這里采用AUSMPW+格式[4],將數(shù)值通量項寫為

式中馬赫數(shù)分裂函數(shù)M±和壓力項分裂函數(shù)p±分別為

式中ML,R=u L,R/a1/2,u L,R為相鄰左右兩分界面處的標量速度,界面處聲速a1/2定義為

Hmormal為激波法線方向的總焓,γ為比熱容比。

2.2 狀態(tài)方程及化學(xué)反應(yīng)繞流模型

控制方程中,守恒變量組成的矢量Q′為

其中,兩溫度非平衡氣體內(nèi)能E為

e為各組分的平功能、束縛電子激發(fā)能及組分生成能之和。u、ν分別為速度在x、r方向的分量。而組分的振動能量EV為

ci、Mi、eνi分別為 i 組分的質(zhì)量分數(shù) 、分子量和分子振動能量。ce、Me、TV分別為振動項的質(zhì)量分數(shù)、分子量和振動溫度。

因此,可得兩溫度非平衡氣體狀態(tài)方程為

其中ρ為繞流組分密度。

繞流計算的化學(xué)模型采用11組分26個化學(xué)反應(yīng)的Dunn-Kang空氣化學(xué)模型第 j個組分的化學(xué)反應(yīng)生成源項為

其中,γij為平衡常數(shù),χi為第i個化學(xué)反應(yīng)的生成源項為

上式中正、逆反應(yīng)速率常數(shù)取為

A、B、C分別為指前因子、指數(shù)因子的活化能。

流場計算中的繞流場的粘性系數(shù)μ及導(dǎo)熱系數(shù)ktr、kν和ke采用維爾克(Wilke)的半徑驗公避開[13]

式中

角標表示各組分的輸運系數(shù)。

2.3 化學(xué)反應(yīng)繞流場計算

算例中飛行器模型選取球頭半徑Rn=0.15 m、長為1.3 m、半錐角為9°的鈍頭錐模型。飛行器模型本身具有旋轉(zhuǎn)對稱性,且沿飛行器軸線方向飛行,因此流場也具有旋轉(zhuǎn)對稱性質(zhì),本問題簡化為2D問題。算例中飛行器所處空間背景壓強設(shè)為 50 Pa,背景溫度為270 K,飛行器速度為22馬赫。對繞流場的模擬采取了AUSMPW+格式,模擬結(jié)果如圖1所示。從圖1可以看出,飛行器高速飛行時在飛行器表面形成了高溫氣體,由于邊界層的影響,壁面附近溫度峰值可達4800 K左右。

圖1 飛行器表面溫度分布

圖2 為飛行器表面垂直軸線方向氧與一氧化氮的質(zhì)量百分數(shù)分布。由圖2可知,AUSMPW+格式與Roe格式計算結(jié)果符合較好。圖3、圖4給出采用AUSMPW+格式計算的電子密度分布。從圖3可以看出在飛行器頭部區(qū)域電子密度可達1×1018m-3左右,且沿著飛行器軸線方向逐漸減小。圖4表明x=8Rn位置處,電子密度沿垂直于飛行器軸線方向緩慢增加。計算結(jié)果與RAM-C實驗數(shù)據(jù)、DSMC模擬結(jié)果以及Roe格式模擬結(jié)果比較[7-8]表明,AUSMPW+格式計算的結(jié)果相比其他兩種算法更加符合RAM-C實驗數(shù)據(jù)。

3.基于修正WKB方法的電波傳播模型

對于非均勻等離子體中的電波傳播問題,本文從WKB解出發(fā),采用傳輸矩陣方法修正WKB解[14],研究非均勻等離子鞘套中TE波與TM波的傳播特性。非均勻等離子體中的TE(電矢量在y軸方向)波動方程為

考慮到傳播方向上反射波的耦合效應(yīng),方程(18)的解具有如下形式

式中E±(x)的上標“±”表示電磁波沿軸的正向與反向傳播,且滿足如下方程

求解方程(20)時,根據(jù)傳輸矩陣法可知

Q0→x為從鞘套邊界傳播至x位置處的傳輸矩陣。

式中

若入射波為TM波,推導(dǎo)方法類似,其中

4.等離子鞘套中的電波傳播特性

在飛行器高速飛行時,飛行器附近高溫區(qū)域內(nèi)空氣中分子產(chǎn)生強烈的振動、離解和電離,形成等離子體鞘套。電磁波在等離子鞘套中傳播時,將發(fā)生相移、時滯、色散、反射、折射和吸收等效應(yīng)。由于存在碰撞,等離子鞘套的相對介電常數(shù)是一復(fù)數(shù)

式中:ωp為等離子體頻率;υ為電子的碰撞頻率[15]。

式中:ne為電子密度;e為基元電荷;me為電子質(zhì)量;ε0為真空中的介電常數(shù);κ為波爾滋曼常數(shù);T為溫度。

研究電磁波在等離子鞘套中的幅頻特性,定義透射系數(shù)為

由如前所述的等離子體鞘套模型,且空間背景壓強設(shè)為 50 Pa,背景溫度為270 K,采用AUSMPW+格式分別模擬飛行速度為 8 Mach、15 Mach、22 Mach時的等離子鞘套。根據(jù)式(19)和(31)計算不同飛行速度下垂直入射電波在鞘套中的透射系數(shù),如圖5所示。

圖5 透射系數(shù)隨入射電波頻率的變化

由圖5可知,當電磁波入射等離子鞘套時,考慮到等離子體的碰撞,電磁波在其中傳播時的衰減隨入射電波頻率增加而減小。當入射電磁波頻率較低時,等離子體對電磁波的碰撞衰減較小,而飛行器飛行速度增加時,飛行器表面溫度迅速升高,等離子體碰撞頻率也隨之增大,所以對入射電磁波的衰減迅速增大。當入射電磁波頻率較高時,相應(yīng)的等離子體碰撞頻率也較大,因而提高飛行速度對電磁波衰減的影響不如較低頻率時明顯。

電磁波在等離子體中傳輸產(chǎn)生的相移定義為

式中d為等離子體厚度,且

將流場仿真計算的等離子體鞘套參數(shù)帶入式(33),計算傳輸相移。圖6給出不同飛行速度下入射電波在鞘套中產(chǎn)生的相移隨入射頻率的變化。由圖6可以看出在等離子體頻率附近相移達到最大值,隨入射波頻率增大相移逐漸減小,且在較高頻率范圍等離子鞘套產(chǎn)生的相移與飛行速度成正比,飛行速度越大,產(chǎn)生的相移越大。

描述系統(tǒng)相移特性的另一種方法是用群時延特性來表示。群時延定義為相移特性的微分。圖7分別給出了不同馬赫數(shù)等離子鞘套產(chǎn)生的群時延,可以看出在等離子體頻率附近群時延有最大值,且飛行速度越高,群時延越大,22 Mach時最大可達0.6 ns左右。

5.結(jié) 論

本文通過建立的飛行器模型,采用兩溫度模型、11組元的Dunn-Kang空氣化學(xué)模型,由AUSMPW+格式求解NS方程,模擬等離子鞘套內(nèi)的流場分布,得出電子密度剖面,計算鞘套影響下的電波幅頻特性、相頻特性,得到如下結(jié)論:1)飛行器速度越高,鞘套內(nèi)電子密度越大,從高頻-L頻段的電波接近等離子體頻率,衰減明顯,而S頻段以上的較高頻率范圍,頻率越高衰減越小。2)電磁波在等離子鞘套中傳輸?shù)南囝l特性受等離子體頻率以及碰撞頻率的影響,而這兩者取決于鞘套的流場特性,飛行速度越高,鞘套產(chǎn)生的相移越大。3)相移隨著入射波頻率增大而增大,接近等離子體頻率時到達極大值,而后緩慢減小,相應(yīng)的群時延在等離子體頻率附近最大。本文算例中的電子密度剖面源于流場模擬的結(jié)果,相較于以往文獻中的各種近似剖面,其更接近流場中的真實情況,研究成果為航天活動中超高聲速飛行器測控、通信與制導(dǎo)等方面的研究,提供了有效的分析手段。

[1] MATHER D E,PASQUAL J M,SILLENCE J P.Radio Frequency(RF)Blackout During Hypersonic Reentry[R].AIAA/CIRA 13th International Space Planes and Hypersonic Systems and Technologies,2005.

[2] SUKHAREVSKY O I,ZALEVSKY G S,NECHITAYLO S V,et al.Simulation of scattering characteristics of aerial resonant-size objects in the VHF band[J].Radioelectronics and Communications Systems,2010,53(4):213-218.

[3] 莫錦軍,劉少斌,袁乃昌.非均勻等離子體覆蓋目標隱身研究[J].電波科學(xué)學(xué)報,2002,17(1):69-73.

MO Jinjun,LIU Shaobin,YUAN Naichang.On the stealth effect of non-uniform plasma covered radar targets[J].Chinese Journal of Radio Science,2002,17(1):69-73.(in Chinese)

[4] KIM K H,KIM C,RHO O H.Methods for the accurate computations of hypersonic flows I.AUSMPW+scheme[J].Joumal of Computational Physics,2001,174(1):38-80.

[5] KIM K H,KIM C.Accurate efficient and monotonic numerical methods formulti-dimensional compressible flows Part II:Multi-dimensional limiting process[J].Jourmal of Computational Physics,2005,208(2):570-615.

[6] LIU Jun.Experimental and numerical research on thermo chemical nonequilibrium flow with radiation phenomenon[D].Changsha:National University of Defense Technology,2004.

[7] PADILLA J F,BOYD I D.Assessment of Rarefied Hypersonic Aerodynamics Modeling and Windtunnel Data[R].AIAA Paper 2006-3390,2006.

[8] BOYD ID.M odeling of associativeionization reactions in hypersonic rarefied flows[J].Physics of Fluids,2007,19(9):3-14.

[9] EGHLIDI M H,MEHRANY K,RASHIDIAN B.Modified WKB method for solution of wave propagation in inhomogeneous structures with arbitrary permittivity and permeability profiles.Proceedings of the 37th European Microwave Conference[C]//Munich:Causal Productions Pty Ltd,2007:1373-1376.

[10] EGHLIDI M H,MEHRANY K,RASHIDIAN B.Modified differential transfer matrix method for solution of one dimensional linear inhomogeneous optical structures[J].JOpt Soc Amer B,Opt Phys,2005,22(7):1521-1528.

[11] 閻玉波,董 慧,李清亮.等離子體涂覆三維目標散射特性的PLRC-FDTD分析[J].電波科學(xué)學(xué)報,2007,22(4):563-566.

YAN Yubo,DONG Hui,LI Qingliang.Analysis for scattering of 3-dimentional target coated with plasma by PLRC-FDTD technique[J].Chinese Journal of Radio Science,2007,22(4):563-566.(in Chinese)

[12] 王 飛,葛德彪,魏 兵.SO-FDTD法計算磁化等離子體層的反射透射系數(shù)[J].電波科學(xué)學(xué)報,2008,23(4):704-708.

WANG Fei,GE Debiao,WEI Bing.SO-FDTD method for computation of reflection and transmission coefficients for magnetized plasma layer[J].Chinese Journal of Radio Science,2008,23(4):704-708.(in Chinese)

[13] GUPTA R N,YOS J M,THOMPSON R A.A review of reaction rates and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30000K[M].NASA-TM-101528,1990.

[14] KHORASANI S,MEHRANY K.Differential transfermatrix method for solution of one dimensional linear nonhomogeneous optical structures[J].J Opt Soc Amer B,Opt Phys,2003,20(1):91-96.

[15] PETRIN A B.Transmission of Microwaves Through Magnetoactive Plasma[J].IEEE Trans on Plasma Sci,2000,29(3):471-478.