不同尺度飛行器周圍等離子體分布及電磁波傳輸效應

高鐵鎖,江濤,傅楊奧驍,丁明松,劉慶宗,董維中,許勇,李鵬

(中國空氣動力研究與發展中心 計算空氣動力研究所,四川 綿陽 621000)

0 引言

高超聲速技術是未來先進飛行器發展的核心技術之一。在對高超聲速飛行器優化設計時,除了分析評估飛行器的氣動特性和飛行性能,還需要評估飛行器周圍等離子體流動對目標特性和通信性能的影響[1-3]。高超聲速飛行器在大氣中飛行時,與周圍氣體發生強烈作用產生強激波效應,使得波后氣體溫度急劇升高,氣體分子發生振動激發、分解和電離等復雜氣動物理化學現象,在飛行器周圍形成高溫電離氣體層俗稱等離子體鞘套。電磁波通過等離子鞘套傳播時,被等離子體反射、折射和吸收,發生強度衰減、傳播方向偏折和相位畸變等效應,導致電磁波作用距離縮短和信噪比下降等電磁性能退化現象,情況嚴重時,等離子體使電磁波傳輸中斷,出現“黑障”現象,這些現象一般統稱為氣動電磁波傳輸效應[4]。氣動電磁波傳輸效應特別是“黑障”現象嚴重影響飛行器與地面之間的微波通信性能,對飛行器實時控制和飛行安全構成嚴重威脅[5-6]。因此,研究氣動電磁波傳輸效應對“黑障”現象的分析評估及高超聲速飛行器通信系統的設計都具有重要意義。

隨著臨近空間飛行器飛行速度的不斷提高,氣動電磁波傳輸效應對電磁通信的影響問題愈加突出。考慮到飛行試驗模擬的高成本及地面試驗對實際飛行條件模擬能力限制等因素,對氣動電磁波傳輸效應問題的研究目前仍以數值模擬為主要手段。1995年,Lin等[7]總結分析了飛行器等離子體鞘套的計算和試驗研究狀況,給出了均勻等離子體中各種電磁波傳輸效應的解析表述;1998年,Nusca等[8]采用7組分化學非平衡N-S方程和波動方程,計算分析了平面電磁波在無線電衰減測量C(Radio Attenuation Measurement C,RAM-C)再入等離子體鞘套中衰減效應;2002年,Funaki等[9]結合固體火箭發動機地面試驗,預估分析了發動機羽流等離子體對不同波段電磁波的衰減特性;2003年,Starkey等[10]采用Park化學模型和解析方法,分析了航天飛機等多種飛行器等離子體鞘套對電磁波通信的影響;2006年,White等[11]基于高階數值方法,分析了RAM-C再入等離子體鞘套的電磁波傳輸效應;2009年,Kim等[12]采用非平衡等離子體數值模擬方法及傳輸效應解析方法,研究了高超聲速飛行器軌道再進入試驗(Obiter Re-entry Experiment,OREX)繞流等離子體對電磁波傳輸的影響,探索外加磁場降低通信窗口等離子體傳輸效應的方法;2018年,龔旻等[1]對臨近空間飛行器“黑障”現象的數值和地面模擬方法進行了回顧總結;2019年,左光等[13]采用數值及解析方法,分析了大鈍頭返回艙和類X-37B升力體飛行器再入等離子體鞘套對電磁波傳輸的影響特征。上述研究工作的共性特點是把流體力學和電磁學結合起來,分析高超聲速飛行器等離子體鞘套中電磁波傳輸效應,但對不同特征幾何尺度高超聲速飛行器等離子體鞘套及電磁波傳輸效應的產生機制和規律分析較少。

本文基于求解非平衡流場N-S控制方程和求解電磁場波動方程的模擬方法,采用自主研發的氣動物理流場數值模擬軟件和再入黑障預測分析軟件,研究高超聲速再入體周圍等離子體的產生機制及氣動電磁波傳輸效應,重點分析不同頻率電磁波在再入等離子體鞘套中傳播的衰減效應,研究不同球頭半徑尺寸再入體的周圍等離子體的分布規律及其對電磁波衰減的影響特征。

1 高超聲速等離子體流場數值模擬研究

1.1 控制方程

高超聲速飛行器從外層空間再入大氣層過程中,由于熱化學過程與流動過程的特征時間尺度效應,飛行器繞流等離子體一般要經歷熱力學和化學非平衡過程。文獻[14]的研究發現,熱力學非平衡效應對高超聲速飛行器非平衡繞流中等離子體分布的影響很小,此時可以不單獨考慮氣體分子的振動能量模態,而是基于單溫度氣體模型,通過求解含化學反應源項的化學非平衡流動N-S方程,數值模擬高超聲速行器周圍等離子體鞘套,三維化學非平衡流動的N-S方程的無量綱化形式[14]如下：

(1)

式中：Q為守恒變量,Q=(ρi,ρ,ρu,ρv,ρw,ρE)T,ρi和ρ是組分的分密度和混合氣體的總密度,u、v、w為直角坐標下3個方向的速度,E為混合氣體的總能;Re是來流雷諾數;F、G、H和FV、GV、HV分別對應3個方向的無黏和黏性通量項;W為化學非平衡源項,W=(wi,0,0,0,0,0)T,wi為i組分生成源項,

(2)

(3)

1.2 化學模型

高超聲速飛行器與來流空氣發生作用形成脫體激波,波后壓縮空氣溫度急劇增高,使得來流壓縮空氣不僅發生離解反應和置換反應,而且發生各種電離反應,在不計空氣中微量組分的情況下,認為來流空氣由O2和N2組成,此時高溫空氣中的電離反應主要包括締合電離和碰撞電離反應機制[14]：

2 數值計算方法

2.1 等離子體流場數值模擬方法

基于有限差分方法對流動控制方程式(1)進行數值求解：對于方程中的無黏通量項,采用上下對稱的高斯-賽德爾(Lower-Upper Symmetric Gauss Seidel,LU-SGS)隱式處理分法[18],以解決無黏通量雅可比矩陣直接求逆帶來的運算量大的問題;為了解決流動和化學過程特征時間尺度效應導致的剛性問題,采用全隱式耦合方法對方程中化學反應源項進行處理[19],即用同一時間尺度同時求解流動方程和化學反應方程以保證計算的穩定性和收斂性;為了準確模擬流動中激波結構和邊界層特性,同時保證計算的魯棒性,采用高分辨壓力權函數修正的迎風型矢通量分裂(Advection Upstream Splitting Method by Pressure-based Weight functions ,AUSMPW+)格式[20]離散無黏通量,采用中心差分格式離散項黏性通量項;為了模擬高空稀薄效應的影響,引入壁面參數滑移修正模型[21];為了模擬多組分混合氣體輸運特性,混合氣體的黏性系數和熱傳導系數用Wilke半經驗公式計算,各組分的輸運系數基于Blotter曲線擬合公式和Eucken關系式計算,擴散系數采用等效二元擴散模型計算,具體計算方法詳見文獻[18-19]。

2.2 等離子體中電磁波傳輸特性預測方法

考慮一維情況,即等離子體參數僅在z軸方向非均勻分布,且平面電磁波沿z軸正向傳播,此時電場只是隨著z軸發生變化,設電場平行于y軸,此時波動方程[22-24]為

(4)

式中：Ey為電場強度;k為波數。對于緩變介質,即介質的電磁參數在z軸方向變化較小,該方程的Wentzel-Kramers-Brillouin(WKB)解[24]為

(5)

式中：Ey0為z=0 m處的電場強度。設電磁波從z=0 m處垂直入射到等離子體內部,并在z=d界面處透射出來,電磁波在此處的能量衰減為

(6)

在各向同性非磁化等離子體介質中波數k為復數,其表達式為

(7)

式中：ω與ωp分別表示電磁波與等離子體的角頻率;c為光速;r為等離子體介質的相對介電常數;ν為等離子體的碰撞頻率,碰撞頻率基于工程經驗關系進行計算[4]。

3 數值計算結果

3.1 不同天線位置等離子體分布及其電磁傳輸特性

研究狀態為RAM-C再入飛行狀態[25]：再入體為球錐對稱體外形,再入高度H=71 km,飛行速度V=7.65 km/s,壁面催化條件為非催化壁面(Non Catalytic Wall,NCW),壁面溫度Tw=1 500 K,飛行攻角近似為0°。在零攻角飛行狀態下等離子體流動為軸對稱流動,等離子體分布沿軸向變化是主要的,因此通信天線的軸向安裝位置將嚴重影響電磁波的傳輸效果。

圖1 不同軸向位置法向峰值與臨界等離子體分布Fig.1 Peak and critical plasma parameters distribution at different axial positions

實際上,再入體周圍等離子體鞘套處于空間(定常情況下)非均勻分布狀態,繞流等離子體厚度和大小分布均在變化。為了定量分析再入體周圍等離子體鞘套對電磁波衰減的影響,以再入等離子體參數為基礎,認為電磁波沿物面法向傳播,采用上述WKB方法預測電磁波通過等離子體的衰減量,如圖2所示。從圖2(a)可以看出,對于C頻段和X頻段高頻電磁波,其在等離子體中的衰減沿軸向總體上變小,而低頻VHF電磁波的衰減量在球頭附近區域出現極小值和極大值,使得電磁波衰減沿軸向呈現非線性變化。盡管越靠近頭部峰值電子數密度越高,但頭部等離子體鞘套的厚度較小,且等離子體參數沿物面法向和流向非均勻分布,電磁波衰減隨等離子體參數(電子數密度和碰撞頻率)非線性變化,導致不同頻率電磁波衰減沿軸向變化出現不同變化規律。總之,不管是低頻(VHF頻段)還是高頻(C、X頻段)電磁波,其通過飛行器后部等離子體鞘套的衰減量均小于通過頭部等離子體的衰減量,而且在同一軸向位置處,電磁波頻率越高其衰減量越小,因此把天線位置安裝在靠近飛行器后身部,或適當選擇頻率較高的電磁頻段,有利于減少飛行器等離子體鞘套對電磁通信的影響。

圖2 不同軸向位置的電磁波能量衰減和峰值碰撞頻率Fig.2 Electromagnetic energy attenuation and peak collision frequency at different axial positions

為了分析等離子體碰撞頻率的影響,圖2(a)中還給出了理想無碰撞等離子體(假定碰撞頻率ν=0 Hz)的預測結果。從圖2(a)可見,碰撞頻率對 C頻段和X頻段電磁波的衰減的幾乎沒有影響,而對低頻VHF電磁波在球頭附近區域等離子體中的衰減有較大影響。圖2(b)給出了物面法向等離子體剖面的峰值碰撞頻率隨軸向的變化情況。由圖2(b)可知：C頻段和X頻段的電磁波頻率顯著高于等離子體的碰撞頻率,而VHF頻段電磁波頻率低于再入體球頭附近等離子體的碰撞頻率;碰撞頻率對電磁波衰減的影響主要與碰撞頻率與電磁波頻率的相對大小有關[22],即當電磁波頻率接近或小于等離子體碰撞頻率時,碰撞頻率對電磁波衰減才會表現出來,二者差異越大,等離子體碰撞頻率對電磁波衰減的影響越大。

3.2 不同尺寸飛行器等離子體分布及電磁傳輸效應

以球頭半徑Rn=15.24 cm的RAM-C再入等離子體鞘套[25]為基礎,改變球頭半徑尺寸,保持長度不變(L=1.295 m),NCW條件,壁面溫度Tw=1 500 K,飛行攻角近似為0°。分析不同球頭半徑(Rn為15.24 cm、7.5 cm、2.5 cm)對再入體周圍等離子體分布的影響。

圖3給出不同球頭半徑飛行器等離子體分布,圖中每條等值線上的數字表示對應的電子數密度。可見隨著球頭半徑減小,等離子體鞘套層的厚度明顯減小,鞘套層的外邊界越貼近壁面。圖4進一步給出不同尺寸再入體頭部駐點線的流場溫度和電子數密度分布情況,圖中R為壁面法向距離。隨著球頭半徑增大,沿駐點線激波層的厚度也隨著增大,激波電離層電子數密度的峰值也隨之增大,如圖4(b)所示。大尺寸球頭激波層氣體電離度的增強意味著較強電離效應,而氣體離解和電離等吸熱反應使得大尺寸球頭來流氣體過激波后沿駐點線溫度下降得更快,如圖4(a)所示。實際上,在同一飛行條件下,飛行器頭部尺寸越大,脫體激波面的曲率半徑越大。曲面激波可看作是多個微元段平面斜激波段的組合,曲率半徑越大,對應每個微元面的激波角越大,頭部接近正激波面元段就越多,激波對來流氣體的壓縮效應就越強,波后壓縮氣體溫升使得化學電離效應就越強。此時高速來流氣體通過激波滯止過程中動能減少,一部分動能轉化為氣體粒子的熱運動能,另一部則轉化為離解和電離產物的零點能。

圖3 不同球頭半徑再入體周圍等離子體鞘套特征Fig.3 Plasma sheaths around the reentry body with different sphere radii

圖4 不同球頭半徑再入體壁面法向等離子體分布(x=0 m)Fig.4 Wall-normal plasma distribution (x=0 m)

圖5給出再入體身部1.234 m處(通信天線附近)壁面法向剖面的NO+質量分數和電子數密度分布,比較圖4(b)和圖5(b)不難看出,從頭部駐點區域到尾部天線附近1.234 m處,對于不同尺寸的再入體,等離子體鞘套剖面內平均電子數密度降低均超過兩個量級,而且頭部半徑越小,等離子體中電子數密度衰減越快,這是因為沿再入體長度方向等離子體特征電子數密度與球頭半徑的冪次方呈正比關系[6,22],球頭半徑越大,從頭部開始等離子體層厚度相對于上游增加的速度越快,等離子體電子數密度相對于上游衰減越慢。從圖5還可以發現,在天線附近1.234 m處,等離子體剖面內平均電子數密度和NO+質量分數隨頭部半徑的變化規律基本一致,二者均隨著頭部半徑的減小而減小,表明此時NO的電離效應對等離子體電子數密度分布起主導作用,這可以從圖6中各種離子組分的分布特征得到進一步印證。從圖6可以看出,對于不同頭部半徑尺寸的再入體,在此尾部天線附近位置處,沿物面法向的所有離子組分分布中,NO+質量分數最大,是最主要的電離組分,說明有關NO電離反應機制對等離子體中電子數密度分布的貢獻最大。

圖5 再入體壁面法向等離子體分布(x=1.234 m)Fig.5 Wall-normal plasma distribution (x=1.234 m)

圖6 再入體壁面法向離子組分質量分數分布(x=1.234 m)Fig.6 Ion mass fraction along wall normal (x=1.234 m)

圖7給出球頭半徑對等離子體中電磁波傳輸特性影響,可見,隨著球頭半徑增加,由于峰值電子數密度隨之增加,等離子體對電磁波傳輸過程中產生的衰減量也隨之增加。實際上,隨著球頭半徑增加,等離子體的厚度也隨之增加,即球頭半徑對等離子體參數大小與厚度分布均產生影響,從而影響電磁波在其中的傳輸特性。從圖7中還可以看出,再入體的球頭半徑改變時,等離子體對不同頻率電磁波衰減的影響的程度有所不同,但電磁波衰減量沿軸向的變化規律基本一致。

圖7 球頭半徑對電磁波衰減特性的影響Fig.7 Effect of sphere radius on energy attenuation of electromagnetic wave

3.3 OREX通信中斷預測分析

在對飛行器OREX通信中斷預測分析之前,先對本文計算模型方法及軟件的可靠性進行驗證分析。沿飛行彈道的高超聲速飛行器周圍等離子體及通信中斷數據可以通過飛行測量獲得,但目前公開發表飛行測量數據很少,而地面試驗還不能完全模擬實際飛行條件,且缺乏天地換算相似準則。RAM-C飛行器再入試驗提供了相對全面包括等離子體分布和通信中斷的測量數據[12,25],該數據是目前少有的能夠同時驗證等離子體流場及其對電磁波通信影響預測方法的測量數據,這里基于此測量數據及彈道條件開展驗證分析。數值模擬等離子體流場采用7組分Park化學反應模型[16],壁面采用完全催化壁面(Full Catalytic Wall,FCW),壁面溫度為1 500 K,通信天線位置在x=6.4Rn處,飛行器以零攻角再入。圖8(a)給出典型再入高度條件下剖面峰值電子數密度沿軸向分布,可見數值模擬結果與飛行測量結果符合很好。在沿彈道對RAM-C飛行器再入流場等離子體分布數值模擬的基礎上,計算了對應飛行彈道條件下電磁波通過天線附近等離子體的能量衰減,計算中認為電磁波沿天線壁面位置的法向傳播。一般認為電磁波衰減30 dB以上時發生通信中斷[23],這樣就可以通過預測的衰減量判斷是否發生通信中斷,獲得沿彈道的通信中斷區間,如圖8(b)所示,其中測量數據包括VHF波段通信中斷區間,中斷時間持續約半分鐘[6],而飛行測量只獲得C與X波段通信中斷區間,因此圖中只比較了這兩個波段通信中斷的起始高度。從圖(8)可見,RAM C-II沿彈道再入過程中,VHF、C和X波段電磁波通信中斷預測與測量數據具有較好的一致性,預測的通信中斷起始或結束高度最大誤差在5 km以內。

圖8 RAM-C飛行器再入預測結果與飛行測量對比Fig.8 Comparison between predictional and experimental data for RAM-C vehicle reentry

下面針對OREX條件[12,26]開展計算分析,其外形尺寸、等離子體探針及通信天線布置情況見圖9。在OREX飛行器以零攻角再入[12]飛行條件下,飛行器繞流對稱,探針附近等離子體分布可以反映通信天線附近的等離子體特征。圖10給出了FCW和NCW條件下電子數密度的數值模擬結果與電子探針測量結果的比較,其中橫坐標表示探針分布方向,離散點表示飛行器肩部分布探針的測量結果。由圖10可以看出,數值模擬結果在總體上和飛行測量結果具有較好的一致性,高度84.01 km的符合程度更好,高度79.9 km的計算結果相對測量結果偏高,但計算誤差在3倍左右,這與之前對RAM-C飛行器再入等離子體預測精度相當[14],表明了數值模擬等離子體模型算法及軟件的可行性。至于高度79.9 km的數值模擬結果偏高的原因,很可能是高溫空氣反應特別是電離反應速率常數數據的不確定性所致,特別是高溫條件下,不同化學動力學模型對等離子體分布的數值模擬結果可存在量級上差異[27],而高度79.9 km下流場中氣體電離效應相對84.01 km更強,有關電離反應化學動力學數據對數值結果的影響更大。圖11給出OREX飛行器在高度59.6 km條件下繞流等離子體及頭部駐點線上離子組分的分布情況。從圖11(a)可見,飛行器頭部激波層內等離子體效應很強,從頭部到飛行器接近尾部的肩部位置,等離子體濃度沿流向雖然有所減弱,但衰減速率比較緩慢,這與RAM-C飛行器的等離子體沿流向的變化特征存在明顯不同(參見圖1和圖3),RAM-C飛行器周圍等離子體濃度沿流向的衰減效應相比OREX飛行器更強,這是由于二者在飛行條件接近的情況下,OREX飛行器屬于大鈍頭體外形,其球頭半徑和錐角都明顯大于RAM-C飛行器球頭半徑和錐角,從3.2節的分析可知,飛行器頭部球頭尺寸及錐角的大小造成了飛行器周圍等離子體分布及沿流向衰減程度的不同。從圖11(b)中電離組分的分布可見,頭部激波層等離子體流場中,NO的電離機制是最主要的,NO+離子對等離子體濃度的貢獻最大,其次是O的電離效應,其他組分電離效應較弱。圖12給出了OREX飛行器天線附近不同高度和不同催化壁條件下的壁面法向等離子體的分布情況,圖中標識H84.0V7416表示再入高度84.0 km與對應再入速度7 416 m/s,其他類同。從圖12不難看出,FCW和NCW的計算結果只是在靠近壁面附近有所差別,在離開壁面的大部分流場區域二者非常一致,沿再入彈道發展過程中,等離子體強度在高度 67 km 左右達到最高。圖13是以此等離子體流場數據為基礎獲得的不同頻段電磁波能量衰減沿再入彈道高度的變化情況,計算時認為電磁波的傳播方向沿壁面安裝天線位置的法線方向,不難發現,沿再入彈道等離子體對電磁波衰減沿再入高度降低逐漸增大并在高度67 km左右達到峰值然后開始減小。這是因為對于特定飛行器外形和飛行姿態,來流大氣密度和再入速度對飛行器周圍等離子體的分布其主導作用[12,22],而等離子體強度的沿高度分布規律決定了通信電磁波的衰減規律,在電磁波衰減出現峰值以前,再入速度沿高度降低減速較慢,大氣密度沿高度降低增大,對等離子體分布和衰減起主導作用,而在電磁波衰減出現峰值以后,再入飛行速度沿高度快速降低成為等離子體衰減的主導因素。從圖13還可以看出,FCW與NCW的衰減計算值相差不大,在認為電磁衰減30 dB以上發生通信中斷情況下[22-23],X波段電磁波在高度84 km左右開始出現中斷現象,而L和C波段出現中斷的高度更高,由于再入段末端飛行減速很快,3個波段電磁波在高度50 km左右時衰減迅速降低而恢復通信,可見隨著通信電磁波頻率提高,等離子體對電磁波衰減及通信中斷的影響減弱,發生通信中斷的高度區間縮小。

圖9 OREX飛行器外形尺寸、靜電探針及天線位置Fig.9 Geometry with electrostatic probes and antenna of OREX vehicle

圖10 OREX等離子體數值模擬與飛行測量值的比較Fig.10 Comparison between computational and experimental data of OREX

圖11 OREX等離子體組分質量分數(H=59.6 km)Fig.11 Plasma and ion mass fraction for OREX (H=59.6 km)

圖12 天線附近壁面法向電子數密度分布Fig.12 Electron number density near antenna

圖13 OREX等離子體對電磁波衰減影響Fig.13 Electromagnetic wave attenuation due to OREX

4 結論

本文針對RAM-C和OREX鈍頭體高超聲速試驗飛行器,開展不同特征尺度下飛行器周圍等離子體分布特性及氣動電磁波傳輸效應的數值模擬和分析。得出以下主要結論：

1)在同一再入條件下,對于鈍頭體飛行器,保持長度不變而改變其球頭尺寸,隨著球頭尺寸的增加,球頭脫體激波內壓縮氣體的化學效應更強,來流氣體動能更多地轉化為波后氣體離解和電離的化學能,導致飛行器周圍等離子體鞘套厚度、等離子體的濃度以及對傳播電磁波的衰減均隨之增大,NO電離反應對飛行器天線附近等離子體分布特性起主導作用。

2)再入飛行器等離子體鞘套空間分布決定了對電磁波傳輸效應的不同特征：飛行器頭部等離子體對電磁波衰減效應更強,隨著流動向下游發展等離子體效應減弱;而在同樣等離子體分布條件下,隨著電磁波頻率的提高,等離子體對電磁波的衰減越弱;可以通過合理選擇天線安裝位置和適當提高通信電磁波頻率的方法減緩或消除等離子體鞘套對微波通信的影響,天線安裝最好選擇飛行器的后身部位。

3)大鈍頭高超聲速飛行器(OREX飛行器)周圍等離子體鞘套厚度大且電離度強,沿流動方向衰減慢,在沿再入彈道飛行過程中,隨著飛行高度降低,等離子體鞘套對電磁波衰減逐漸增強出現峰值,然后隨之逐漸減弱;隨著通信電磁波工作頻率提高,等離子體鞘套對電磁波通信中斷的影響區間縮小。

4)針對自主發展的計算模型方法和軟件開展驗證分析,典型飛行狀態下等離子體分布的數值模擬結果與飛行測量數據符合較好;典型飛行器沿彈道再入過程中,通信中斷預測與測量結果較為一致。該計算分析工具可為高超聲速飛行器的通信系統設計提供技術支持。