等離子體電磁散射的多物理場研究

逯貴禎，洪楚雨

(中國傳媒大學(xué) 信息工程學(xué)院，北京 100024)

1 引言

等離子體是一種物質(zhì)的存在形式，就像固體、液體和氣體一樣。物理學(xué)上對等離子體的定義是，由大量帶電粒子組成的氣態(tài)物質(zhì)。與固體、液體和氣體不一樣，等離子體是由電子、離子和中性原子的混合物組成的。這些帶電粒子可以自由移動并相互作用。因此，等離子體的性質(zhì)明顯地區(qū)別于其它三種物質(zhì)形態(tài)，具有特有的行為和規(guī)律。

等離子體振蕩是一種集體運動模式，由于電子的等離子頻率遠(yuǎn)大于離子的等離子頻率，所以通常用電子等離子頻率作為等離子體的振蕩頻率。在電磁波的等離子散射研究中，當(dāng)電磁波的頻率遠(yuǎn)大于等離子頻率時，等效介電常數(shù)可以近似為1，此時等離子體對電磁散射的影響很小。然而當(dāng)電磁波的頻率接近或小于離子頻率時，介電常數(shù)就會接近于零或是負(fù)值，對電磁散射的影響就會比較大。同時，電子與離子的碰撞頻率會對電磁波的傳輸產(chǎn)生損耗，由于這些影響的存在，等離子的介電常數(shù)有一個和碰撞頻率有關(guān)的虛數(shù)部分。等離子體的性質(zhì)和許多因素有關(guān)，這些因素包括速度場、溫度場、分子之間的化學(xué)反應(yīng)性質(zhì)。

高速飛行物體的等離子鞘，是飛行器在大氣飛行過程中，由于激波壓縮和摩擦加熱所產(chǎn)生的高溫電離氣體層。它與高速飛行器的通訊、制導(dǎo)、遙測、引爆以及電子對抗等問題有著密切的關(guān)聯(lián)。例如，等離子鞘會使飛行器上電磁天線系統(tǒng)的工作受到重大干擾，甚至?xí)雇ㄓ嵧耆袛?。對等離子鞘的研究，大致從六十年代初開始，其涉及的方面相當(dāng)廣泛。從研究目的來看，包括電離機制、電磁效應(yīng)及減輕措施等。從研究深度來看，從低功率電磁波與等離子體的線性相互作用，到高功率的非線性作用，都有所涉及。從研究廣度來看，包括信號衰減、擊穿特性、湍流噪音等。等離子體的特點是有相當(dāng)大的頻率色散，也就是說其介電系數(shù)會隨電磁波的頻率改變而改變。等離子體的傳播特性是電磁波頻率的函數(shù)。其中，等離子頻率是一個重要的特征參數(shù)，它與等離子體中電子密度分布有關(guān)。當(dāng)工作頻率大于等離子頻率時，等離子體就稱作欠稠密等離子體，電波可以透過等離子體傳播;反之，當(dāng)工作頻率小于等離子頻率時，則稱作過稠密等離子體，電波按指數(shù)形式衰減;工作頻率等于等離子頻率時，則稱作臨界等離子體，電波在表面會發(fā)生全反射，一般稱這個頻率為截止頻率。由于等離子體的溫度相當(dāng)高，其壓力梯度的作用便不可忽略。當(dāng)帶電粒子在電磁波作用下發(fā)生強迫振蕩時，壓力會發(fā)生變化，引起帶電粒子密度的變化。與這種電荷累積相關(guān)聯(lián)的交變電場的傳播，便是電聲波。這種變化出現(xiàn)時，鑒于有一部分電磁能量會轉(zhuǎn)換成電聲縱波的功率，其比例取決于等離子體的狀況，因而等離子散射性能亦發(fā)生改變。此外，等離子鞘參數(shù)在空間有一定的分布，其變化的尺度可以與電磁波的波長相比較，因此要考慮非均勻效應(yīng)電磁散射特性。

在國內(nèi)的研究方面，文獻(xiàn)［1］采用分段線性電流密度遞歸卷積時域有限差分方法計算導(dǎo)電金屬半球高超聲速繞流流場電磁散射特性，分析半球等離子體包覆繞流流場雷達(dá)散射截面(RCS)隨入射電磁波頻率、雙站角、飛行馬赫數(shù)和高度的變化特性，其中繞流流場等離子體特征頻率及電子碰撞頻率利用了參考文獻(xiàn)中的計算結(jié)果。文獻(xiàn)［2］采用迎風(fēng)型矢量分裂格式求解N-S流體力學(xué)方程，計算了等離子鞘中電子密度分布，利用高頻WKB方法計算了電波在等離子層中的傳播特性。文獻(xiàn)［3］發(fā)展了一套隱式數(shù)值求解化學(xué)非平衡三維層流及湍流邊界層的高效計算方法，并成功地應(yīng)用于球鈍錐和飛船返回艙這樣的倒錐體在20°大攻角下的計算，得到了包括背風(fēng)區(qū)在內(nèi)的熱流分布和等離子體特性參數(shù)。文獻(xiàn)［4］利用真實氣體效應(yīng)情況下等離子體流場計算方法和移位算子時域有限差分電磁散射建模方法，計算和分析了以零攻角再入時的低頻電磁散射特性。文獻(xiàn)［5］以等離子體的介質(zhì)模型為基礎(chǔ)，給出了等離子體雙站RCS公式與計算實例，分析了等離子體電參量對欠稠密等離子體柱雙站RCS的影響。文獻(xiàn)［6］研究了再入目標(biāo)特性理論模型框架，建立了再入飛行器目標(biāo)特性的基本計算方法，對典型再入飛行器的紅外輻射和電磁散射特性進(jìn)行了數(shù)值模擬，其中對于等離子體包覆體的RCS特性，把層流等離子體等效為分層電介質(zhì)，求解分層電介質(zhì)包覆導(dǎo)體的散射場。采用物理光學(xué)方法求解每一面元散射場，然后對散射場矢量求和。對于尾跡的RCS的特性，關(guān)鍵是亞密湍流等離子體的電磁散射特性。采用一階畸變波Born近似方法，電子密度譜函數(shù)采用hkarofsky譜函數(shù)形式。文獻(xiàn)［7］對再入飛行器等離子體尾跡及其雷達(dá)散射特性進(jìn)行了分析、研究和大量的計算。討論了物形、流場各因素對尾跡雷達(dá)散射截面的影響。流場計算使用準(zhǔn)一維粘性尾跡方程，以修正基爾方法(多值法)求解，用一階Born近似完成亞密雷達(dá)散射截面(RCS)計算。計算中使用8組元混合空氣、14個非平衡化學(xué)反應(yīng)模型，考慮5種不同尺度的小鈍頭錐形物體，沿再入軌道取65至34公里，共13個高程的飛行條件。通過計算得到了再入體尾跡各流場參數(shù)、電子密度分布及湍流亞密尾跡的RCS。

在國外研究方面，文獻(xiàn)［8］進(jìn)行了覆蓋等離子的球體前向散射試驗，根據(jù)實驗結(jié)果推導(dǎo)了電子密度的徑向分布密度模型。文獻(xiàn)［9］研究了不同金屬表面與大氣的摩擦系數(shù)，研究了大氣壓力和摩擦之間的關(guān)系，這些研究主要是等離子體生成的機理研究。文獻(xiàn)［10］發(fā)展了一種新的關(guān)于連續(xù)方程的數(shù)值算法，進(jìn)行高溫、化學(xué)非平衡狀態(tài)下高超聲流體的高效計算;同時研究了N-S方程的改進(jìn)算法，并用于了等離子體的計算。文獻(xiàn)［11］采用熱化學(xué)和非平衡熱動力學(xué)研究軸對稱物體的粘性激波層，與一維激波管的試驗進(jìn)行比較分析，討論了等離子體的產(chǎn)生特性。文獻(xiàn)［12］分別采用粒子方法和連續(xù)介質(zhì)方法計算分析了高超聲速流場中的物體，分析了物體形狀與流場特性的關(guān)系。文獻(xiàn)［13］研究了再入飛行器的黑障問題，提出了采用電磁場降低再入飛行器激波層中電子密度的方法。文獻(xiàn)［14］使用磁流體模型靜電多流體反應(yīng)模型研究弱電離的等離子體問題，并研究了高速流體流過鈍狀飛行器的周圍等離子體分布特性，計算了電磁波在等離子體中的傳播性質(zhì)。文獻(xiàn)［15］針對等離子體研究中的化學(xué)反應(yīng)模型，對常見的 Gardiner、Moss、Dunn and Kang和Park化學(xué)反應(yīng)模型進(jìn)行了比較分析，研究了對于不同馬赫數(shù)條件下化學(xué)計算模型的可信度。文獻(xiàn)［16］研究了亞密湍流邊界層的非相干散射截面，研究了導(dǎo)體平面和柱體平面周圍等離子界面層的影響。文獻(xiàn)［17］研究了在低緯度過密與亞密氣體等離子體的電磁散射特性。研究過程中對等湍流參數(shù)相關(guān)函數(shù)類型、相關(guān)長度、湍流的各向同性、電導(dǎo)率等參數(shù)做了假設(shè)。文獻(xiàn)［18］研究了再入飛行器在軸對稱情況下的超高聲速流體問題，利用Park’s兩溫度模型模擬熱化學(xué)的非平衡狀態(tài)和弱電離效應(yīng)，并對N-S方程進(jìn)行修正，使用有限體積方法分析求解，得到了壓力、熱傳輸速率、電子密度的變化規(guī)律，并與試驗進(jìn)行了比較。文獻(xiàn)［19］研究了熱非均勻磁等離子體的參數(shù)耦合問題，研究了非線性波的相互作用，該結(jié)果可用于確定電場頻率偏移的門限值。文獻(xiàn)［20］研究了再入飛行器表面在高速流體和高溫下的化學(xué)反應(yīng)問題。綜上所述，目前對等離子鞘的研究可以總結(jié)為:1)研究高速流動條件下的等離子體電子密度分布;2)等離子鞘的電磁散射研究。

在等離子體散射的研究中，數(shù)值計算仿真是一個非常重要的方法。目前在電磁散射的分析中，主要采用FDTD、高頻方法、有限元方法。在流體分析中，主要采用有限差分方法、有限體積方法、有限元方法。對于流體分析，有限差分方法相對比較直觀，易于編程和并行計算，但是難以處理不規(guī)則區(qū)域，對區(qū)域的連續(xù)性等要求比較苛刻。有限元方法對偏微分方程的離散較容易，適合處理復(fù)雜區(qū)域，并且計算精度可靠。對于能使用偏微分方程描述的物理問題，都能使用有限元方法進(jìn)行模擬。有限體積法適于流體計算，可以應(yīng)用于不規(guī)則網(wǎng)格，但由于有限體積法的截取誤差是不定的，它的精度基本上只能是二階。因此，在實用性、適用性以及擴展性方面，有限元方法具有更大的優(yōu)勢，也是現(xiàn)在應(yīng)用最為廣泛的一種數(shù)值計算方法。一般來說，物理現(xiàn)象都不是單獨存在的。這種物理系統(tǒng)的耦合就是我們所說的多物理場。有限元法在多物理場方面的應(yīng)用，有著得天獨厚的優(yōu)勢。在以下研究中，我們采用基于有限元的多物理場分析方法研究等離子體的電磁散射問題。

2 多物理場分析方法

在多物理場分析方法中，為了得到等離子體的電磁散射特性，最重要的是要得到等離子體中介電常數(shù)的空間分布特性，而介電常數(shù)的空間分布與電子密度分布有著密切的關(guān)系，因此需要求解分子的擴散與對流方程，并要考慮離子生成的化學(xué)反應(yīng)過程，而這些過程又與溫度場的分布有密切的關(guān)系，因此需要考慮這些物理現(xiàn)象的綜合效果。

2.1 溫度場分析

溫度場的分布與熱量傳遞的方式有著密切的關(guān)系，熱量傳遞有三種基本方式:傳導(dǎo)、對流和輻射。傳導(dǎo)是分子的熱運動而產(chǎn)生的熱量傳遞，其中傳導(dǎo)系數(shù)是表征物質(zhì)性能的重要參數(shù)。對流是由于物質(zhì)中的流體宏觀運動產(chǎn)生的熱量傳遞過程。流體的速度是表征對流影響的重要參數(shù)。輻射是以通過電磁波的方式進(jìn)行熱傳遞。在我們目前所研究的問題中，輻射的影響可以忽略。研究溫度場分布的基本方程是熱傳輸方程，該方程具有如下形式:

其中T是溫度，ρ是密度，Cp是等壓條件下的熱容量，k是熱導(dǎo)率，Q是熱源，u是流體速度。采用有限元方法研究溫度場時，邊界條件是非常重要的，不同的邊界條件會對形成的有限元方程組產(chǎn)生很大的影響，而這些因素會對求解的精度產(chǎn)生影響。此外，對于輸運問題，如果對流項占據(jù)主要部分，會對解的穩(wěn)定性產(chǎn)生影響，原因是溫度場可能會出現(xiàn)急劇的變化。此外，有限元解如果出現(xiàn)振蕩現(xiàn)象，會影響解的收斂性，為了克服解的收斂性問題，可以采用加入人工擴散因子的方法。

2.2 電子密度空間分布

在非平衡條件下，由于系統(tǒng)中定向運動的動量、溫度、粒子數(shù)密度等因素均存在空間不均勻性，因而會產(chǎn)生定向動量、能量、質(zhì)量等的輸運過程。在等離子體中，由于溫度分布的不均勻性，空間不同位置粒子的密度不一樣，會產(chǎn)生粒子從密度高的地方向密度低的地方的轉(zhuǎn)移過程。粒子密度的輸運可以通過擴散和對流兩個過程來實現(xiàn)，粒子質(zhì)量輸運方程可以表示為:

其中c是粒子的濃度分布，D是擴散系數(shù)，u是流動速度，R是粒子產(chǎn)生速率。在等離子研究中，為了得到電子密度的分布，需要了解分子的電離過程。在大氣等離子體中，氣體分子的電離涉及幾十種化學(xué)反應(yīng)過程，這些反應(yīng)過程與空間的溫度分布有密切的關(guān)系。由于對于每一類粒子，都會有相應(yīng)的粒子數(shù)分布，因此完全求解各種粒子的空間分布是一個計算量非常的任務(wù)。但是對于電磁散射問題來講，研究者最關(guān)心的是電子密度的空間分布。所以我們對電離過程進(jìn)行了簡化，只考慮電子的產(chǎn)生過程。電子密度分布與溫度場的耦合關(guān)系主要體現(xiàn)在電離過程的熱分解過程中。

2.3 電磁散射分析

研究等離子體的電磁散射，最主要的是需要了解等離子體的介電常數(shù)空間分布特性。等離子體的介電常數(shù)有三種模型，它們是Lorentz-Drude模型、Drude模型和德拜模型。對于大氣等離子體，通常采用Drude模型，因為該模型反應(yīng)了氣體分子沒有約束振動的實際情況。Drude模型可以表示為:

其中ωp是等離子頻率，它和電子密度有關(guān)，并且是空間位置的函數(shù)。v是碰撞頻率，與粒子數(shù)密度、溫度有密切的關(guān)系。在很多情況下，碰撞頻率與粒子密度的比例關(guān)系大約是0.0001。碰撞頻率不為零，意味著介電常數(shù)具有虛部，物理上對應(yīng)著能量的耗散，會對電磁波傳輸產(chǎn)生損耗。電磁散射分析采用麥克斯韋方程:

在電場的微分方程中，我們特意把介電常數(shù)隨空間變化的關(guān)系表示出來。介電常數(shù)與電子密度的關(guān)系體現(xiàn)在等離子頻率的計算。

3 計算結(jié)果與討論

為了分析等離子體的散射問題，先考慮一個金屬圓柱的電磁散射問題。金屬圓柱的半徑為0.2米。電磁波的入射頻率為3GHz，波長為0.15米。在多物理場的分析中，我們采用有限元方法。對于任何一個物理場，都可以采用有限元方法對偏微分方程進(jìn)行離散化處理，這也是有限元處理多物理場問題的一個優(yōu)點。分析過程分為三個步驟，首先求得溫度場的分布;其次求電子密度的分布;最后進(jìn)行電磁散射分析。

3.1 溫度場分析

用有限元方法分析溫度場，利用公式(1)作為分析方程，密度、熱容量和熱導(dǎo)率均采用空氣介質(zhì)的材料參數(shù)。在分析過程中，首先考慮非流動的空氣情況，因此此時空氣的流速為零，溫度場的空間分布由傳導(dǎo)過程實現(xiàn)。考慮到很多情況下，物體與空氣摩擦產(chǎn)生高溫，這里假定圓柱體的溫度處于高溫條件，設(shè)為1000K;在遠(yuǎn)離柱體的位置仍然為常溫條件，設(shè)為300K。通過求解相應(yīng)的有限元方程，得到圍繞圓柱體的溫度空間分布如圖1所示。

圖1 溫度場空間分布

從圖1可以看到，溫度場是圍繞圓柱體的徑向?qū)ΨQ分布，這是由于假定只有熱傳導(dǎo)過程，沒有考慮對流過程造成的。如果考慮流體的流動速度，溫度場的分布將不再是對稱分布。

3.2 電子密度分布

電子密度分布服從公式(2)的質(zhì)量輸運方程。對于等離子體而言，電離過程涉及化學(xué)反應(yīng)。在等離子散射研究中，作為簡化計算，通常將電子密度處理為不同的空間分布，比如指數(shù)分布、拋物線分布、高斯分布等不同空間分布形式。作為實際情況的近似，電子密度可以通過流場分析的方法得到電子密度的空間分布特性。常見的有采用熱化學(xué)非平衡方法，其中11組元的Dunn-Kang空氣化學(xué)模型是一種公認(rèn)較為接近實際情況的分析模型。在本文的研究中，為了減少計算量，只考慮電子密度的質(zhì)量輸運過程，電子的電離過程采用公式(5)。

將(2)與(5)結(jié)合，通過將溫度場與電子密度分布相耦合，得到電子密度空間分布如圖2所示。圖2中采用的是每立方米摩爾分布，如果要轉(zhuǎn)化為電子密度分布，需要乘以每摩爾的電子數(shù)目。從圖2可以看到，由于溫度場是對稱分布，電子密度分布也是對稱分布。

圖2 電子密度分布(mol/m^3)

3.3 電場分布與雷達(dá)雙站散射截面

研究等離子體包覆體的RCS特性時，需要了解等離子的介電常數(shù)特性，利用前面得到的電子密度分布，可以通過Drude模型，即(3)式得到介電常數(shù)的空間分布，如圖3所示。從目前的計算中可以看到介電常數(shù)介于0和1之間，在柱體附近介電常數(shù)非常小，遠(yuǎn)離圓柱后，介電常數(shù)接近于1。

在電磁散射的研究中，平面波從右向左入射，電磁波的頻率是3GHz，散射場的分布如圖4所示。從電場分布可以看到，后向散射電場很強，正向散射場較弱。根據(jù)散射場的分布，計算的雷達(dá)雙向站散射截面計算如圖5所示。為了比較，圖5給出了圓柱體在自由空間的雷達(dá)雙站散射截面。可以看到，在背散射方向，有電離的散射截面高于自由空間的情況，在前向散射方向，有電離的雷達(dá)散射截面低于自由空間的情況。

圖3 介電常數(shù)空間分布

圖4 電磁波從右向左入射，散射電場分布

4 結(jié)論

綜上所述，本文采用有限元方法對等離子的電磁散射問題進(jìn)行了多物理場的分析，分析是對特定條件進(jìn)行的，從分析中可以得到，等離子的散射過程是一個很復(fù)雜的過程。對雷達(dá)散射截面而言，有時需要增強雷達(dá)散射截面，也有時也可以降低雷達(dá)散射截面。

圖5 雷達(dá)雙站散射截面

［1］陳偉芳，常雨.高超聲速半球繞流流場電磁散射特性分析［J］.航空動力學(xué)報，2009(3):552-557.

［2］李江挺，郭立新，方全杰.高超聲速飛行器等離子鞘套中的電波傳播［J］.系統(tǒng)工程與電子技術(shù)，2011(5):969-973.

［3］楊宏，曹文祥，潘梅林.化學(xué)非平衡三維邊界層及相關(guān)等離子體鞘套的數(shù)值模擬［C］.中國空氣動力學(xué)學(xué)會第十屆物理氣體動力學(xué)專業(yè)委員會會議論文集，2001:121-126.

［4］韋笑，彭世鏐，殷紅成.基于平衡流場的再入飛行器電磁散射特性分析［J］.系統(tǒng)工程與電子技術(shù)，2011(3):506-510.

［5］黃興中，呂善偉.欠密等離子體柱空間散射特性的數(shù)值分析［J］.航空電子技術(shù)，1997(4):36-41.

［6］張志成，高鐵鎖，董維中.再入飛行器目標(biāo)特性建模研究［J］.實驗流體力學(xué)，2007(3):7-12.

［7］牛家玉，許國斌，曹榮達(dá).再入飛行器尾跡流場及其雷達(dá)散射效應(yīng)研究［J］.空氣動力學(xué)學(xué)報，1996(12):422-429.

［8］Twardeck T G，Quinn R G.Microwave Scattering by a Plasma Covered Sphere［J］.IEEE Trans-AP，1969，17(6):823-824.

［9］Mishina H.Atmospheric Characteristics in Friction and Wear of Metals［J］.Wear，1992，152(1):99-110.

［10］Chapman D R，MacCormack R W.Improved Computational Fluid Dynamics for Continuum Hypersonic Flow［R］.Final report for a r o contract daal 03-36-K-0139，1989:1-6.

［11］Mott D R，Gally T A，Carlson L A.Viscous Normal Shock Solutions Including Chemical，Thermal and Radiative Nonequilibrium［J］.Journal of Thermophysics and Heat Transfer，1995，9(4):577-585.

［12］Wang W L，Boyd I D，Candler G V.Particle and Continuum Computations of Hypersonic Flow Over Sharp and Blunted Cones［C］.35th AIAA Thermophysics Conference，2001:45-63.

［13］Keidar M.Electromagnetic Reduction of Plasma Density During Atmospheric Reentry and Hypersonic Flights［J］.Journal of Spacecraft and Rockets，2008，45(3):445-453.

［14］Kundrapu M，Loverich J，Beckwith K.Modeling and Simulation of Weakly Ionized Plasmas Using Nautilus［C］.51st Aiaa Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition，2013.

［15］Tchuen G，Burtschell Y.Physico-Chemical Modelling in Nonequilibrium Hypersonic Flow Around Blunt Bodies［J］.Aeronautics and Astronautics，2013，9:125-158.

［16］Ruquist R D.Electromagnetic Scattering From Turbulent Boundary Layers［R］.MIT Technical Note，1972:1-13.

［17］Luzzi T.Analysis of electromagnetic scattering from turbulent low altitude rocket plumes［R］.Grumman Aerospace Corporation，1972:3-45.

［18］Scalabrin L C，Boyd I D.Numerical Simulation of Weakly Ionized Hypersonic Flow for Reentry Configurations［C］.9th AIAA/ASME Joint Thermophysics and Heat Transfer Conference，2006.

［19］T P Khan.Nonlinear Interaction of Waves in a Hot Inhomogeneous Magnetized Plasma［J］.International Journal of Mathematics and Mathematical Sciences，1979，2(2):325-336.

［20］Prakash A，Zhong X.Numerical Simulation of Planetary Reentry Aeroheating over Blunt Bodies with Non-equilibrium Reacting flow and Surface Reactions［C］.Orlando:47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition，2009.