基于多區域面積分方程的電大非均勻等離子體電磁特性高效計算方法

張慧雯 陳宇浩 黃曉偉 盛新慶

（北京理工大學 集成電路與電子學院射頻技術與軟件研究所, 北京 100081）

0 引 言

當彈道導彈、宇宙飛船等高超音速飛行器在臨近空間飛行時，飛行器前端會形成很強的激波，飛行器與空氣之間的劇烈摩擦使得飛行器周圍空氣發生離解與電離，形成一層高溫等離子體鞘套[1]。等離子體鞘套會對電磁波產生衰減、反射和散射等影響[2-3]，從而干擾飛行器的通信和雷達探測，甚至導致通信中斷，產生“黑障”現象。在對“黑障”現象進行機理分析的過程中，電磁計算可以靈活高效地對電磁問題進行建模與數值求解，成為分析電磁波與等離子體相互作用的有力工具[4-6]。

近年來，國內外涌現出一批求解電磁波與等離子體耦合作用的數值算法[4]。這些算法主要可以分為兩大類，一類是近似方法，主要包括WKB (Wentzel-Kramers-Brillouin)方法[7]、散射矩陣(scattering matrix, SM)法[8-9]、物理光學 (physical optics, PO) 法等[5-6,10]，這些方法在電磁計算過程中使用了不同程度的近似條件，例如模型近似或物理過程近似，在計算復雜的非均勻等離子體時精度嚴重受限；另一類是全波數值計算方法，例如時域有限差分 (finitedifference time-domain, FDTD) 法[11-12]、有限元法(finite element method, FEM)[13]、體積分方程 (volume integral equation, VIE) 法[14]和體面積分方程 (volume-surface integral equation, VSIE) 法[15-16]，這些算法可以計算更精細的等離子體模型，比如通過計算流體動力學(computational fluid dynamics, CFD)仿真得到的三維非均勻流場模型。然而，上述方法均需要在等離子體區域進行體剖分，導致很高的計算資源消耗，尤其是求解大尺度目標時表現乏力。針對天線搭載局部位置的等離子體電磁仿真，或者弱梯度等離子體鞘套包覆目標電磁仿真，三維非均勻等離子體可以化簡為多區域的分區均勻模型，因此只需要在各區域交界面劃分面網格，并使用面積分方程(surface integral equation, SIE)矩量法計算[17]，極大減少了未知數和仿真復雜度，具備計算電大等離子體的潛力。

多區域面積分積分方程(multi-region surface integral equation, MR-SIE)法計算等離子體的一個難點在于區域交界面電磁流連續性的處理。SIE通常使用RWG(Rao-Wilton-Glisson) 基函數進行離散[17]，但 RWG 基函數的連續性使其很難處理多區域目標的區域結合邊界[18]。通常的解決方法是在邊界處增加區域連接模型 (connect region model, CRM) 以保證邊界的連續性條件，但這需要共形區域剖分[19]。然而當目標具有精細幾何結構時，共形剖分十分繁瑣且耗費時間[20]。為解決上述問題，前人提出了區域分解技術 (domain decomposition method, DDM)，以在不同區域中使用非共形且相互獨立的幾何剖分[21]。SIE 方程的區域分解可以分為兩大類：封閉式區域分解 (體分解, volume-based SIE-DDM) 和非封閉式區域分解 (面分解, surface-based SIE-DDM)[22]。封閉式區域分解以體積為基本區域單元，在子區域的連接處需要施加 Robin傳輸條件以保證子區域電流的連續性[23-24]，而這將帶來更多的未知數。非封閉式區域分解方法中，間斷伽遼金 (discontinuous Galerkin,DG) 方法的研究最為深入[25-26]。在不同面區域的邊界處，DG 方法使用半 RWG 基函數進行剖分，同時將電流的連續性條件直接加在半 RWG 的阻抗矩陣元素上。DG 方法無需增加未知數，因此該方法的非共形剖分可以更容易地運用在復雜等離子體多區域目標中[27-34]。

面積分仿真等離子體還面臨高損耗帶來的數值不穩定性。高電子密度的等離子體通常具備很大的介電常數虛部，尤其是電磁波頻率遠低于等離子體頻率時，等離子體的介電常數實部變為負數且虛部數值較大，表現出極強的損耗特性[35]。當多層快速多極子算法(multi-level fast multipole algorithm,MLFMA)運用在SIE方程時，負介電常數目標的計算一直是一項具有挑戰性的工作[36-38]。先前的研究工作在方程的建立[39-42]、迭代器的選擇[43]、預處理的使用[44]等方面均針對負介電常數目標的計算給出了更有效的方案，但依然存在著計算不準確和收斂性較差等問題。經過大量仿真測試，我們發現負介電常數的計算問題還出現在由于數值精度引發的數值不穩定性上。在MLFMA的遠場計算中，一些特殊函數，如貝塞爾函數[45]和移置函數的特性將因為負介電常數的存在而發生改變。對于無損耗介質，這些函數原本是隨距離周期相位變化的，但在負介電常數的損耗介質中，這些函數具有了隨距離強烈衰減或劇烈增長的特性。當這些特性變化的矩陣疊加、連乘或插值時，最終的計算誤差將會很大，甚至引發計算過程中的數值不穩定性。雖然提高計算機的數值計算精度可以在一定程度上緩解問題，但當負介電常數過大時，計算不準確的問題依然會出現，阻礙了SIE法用于電大等離子體電磁計算的實施。

針對上述問題，本文提出了一種基于MR-SIE的電大、非均勻、負介電常數等離子體電磁計算方法。首先，推導了等離子體目標的MR-SIE，并引入金屬與介質混合的DG方法提高對復雜等離子體目標建模和剖分的靈活性。隨后，針對現有MLFMA求解高負介電常數目標面臨迭代變慢、精度變差的問題，使用了一種靈活有效的內問題方程截斷策略。在高負介電常數的介質中，當指數衰減的電磁波傳播距離高于一定閾值時，取消其基于MLFMA架構的傳播計算，避免了數值不穩定性。閾值的選擇和介電常數、矩陣元素的距離以及SIE方程有關，這種策略可以看作是一種等離子體到金屬的過渡，極大提高了MLFMA的魯棒性。該截斷策略可以用于多種介質SIE方程，如結合切向方程(combined tangential formulation, CTF)[46]和電磁流結合場積分方程(electric and magnetic current combined-field integral equation,JMCFIE)[47]。最后，通過不同頻率、不同分布的等離子體電磁輻射和散射問題的數值計算，驗證了本文提出方法的準確性、高效性與其廣泛的應用價值。

1 理論與公式

1.1 等離子體的介電常數

參照洛倫茲-德魯德(Lorentz-Drude)模型[35]，等離子體可以被看作是一種色散損耗介質，隨著入射頻率由低到高，等離子體介電常數從負無窮趨近于1，復數形式的介電常數和磁導率可以表示為

式中：ω,ωp,υeff(單位:rad/s)分別為入射頻率、等離子體角頻率和碰撞角頻率；ε0,μ0分別為自由空間中的介電常數和磁導率。一般使用電子碰撞頻率υe=υeff/(2π)(單位:Hz)來描述等離子體。等離子體角頻率通常與電子密度ne相關，可以表示為

式中：e為電子電荷量；me為電子質量。當ω2p>(ω2+υ2eff)時，相對介電常數的實部為負，和普通介質的介電常數有所不同。

1.2 MR-SIE

首先，考慮一個任意幾何構型的均勻目標區域?i，其介電常數為εi、磁導率為μi。使用等效原理，電場和磁場的總場可以用等效電流磁流J,M表示為：

式中：ηi=為阻抗系數；分別為入射電場、磁場；積分微分算子

Gi為均勻空間?i內的格林函數，表達式為

式中，ki=ω為波數。

之后，考慮多區域結構的SIE。若表面sk為兩個均勻區域?k1,?k2的交界面，則在交界面sk上考慮式(3)的邊界條件建立的切向與法向場積分方程(tangential and normal forms of the electric-field integral equation, T-EFIE and N-EFIE)可以表示為

同理，使用公式(4)建立的切向與法向磁場積分方程(tangential and normal forms of the magnetic-field integral equation, T-MFIE and N-MFIE)可以表示為

結合T-EFIE、N-EFIE、T-MFIE和N-MFIE，可以得到交界面sk上的兩組面積分方程：

式中，aki,bki,cki,dki為不同方程的權重系數。不同的權重組成的方程有不同的數值特性。CTF、JMCFIE這些方程對應的系數[47]為

這些切向和法向方程的系數對于求解的精度和效率是至關重要的。簡單來說，當目標存在金屬片等結構時(以輻射問題為主)，只能使用切向方程來求解，此時使用CTF以確保求解的正確性[48]；當目標不存在金屬片結構，則JMCFIE可以準確計算，同時收斂效果明顯優于CTF，因此在求解以散射為主的問題時一般使用JMCFIE[47]。使用RWG基函數對等效電磁流進行離散，并使用伽遼金匹配，可以得到介質表面的離散矩陣方程：

式中：

若交界面sk為金屬，則表面不存在等效磁流，其構建的方程只與介質側的相互作用有關，因此金屬表面的離散矩陣方程可以簡化為

最后，將多個不同交界面的方程進行結合，構建一個有所有表面未知數的矩陣，即可求解MR-SIE。

1.3 DG方法

DG DDM將原目標的外表面分解為不閉合的子區域，在子區域交界上引入電流連續性條件保證解的等效性。通過對公式(18)進一步展開，可以得到：

式中：C和C′為RWG基函數的邊；p.v.為主值積分。從式(20)～(23)可以看出：不會因為基函數是RWG還是半RWG對K算子的相關計算產生影響；而公式(20)中的線積分項會有因為半RWG的出現無法抵消的情況，其中右端項的第二和第三項的線面積分并不困難，而第四項的線線積分會產生奇異點計算，因此需要特殊處理。當表面為兩個區域的交界面時，邊界元素的處理方式參照文獻[23-24,27]，運用在多區域的交界邊界上的處理方法如下。

圖1所示為一個多區域相交邊界的二維示意圖，l為多區域相交處的交界線，此圖中為中間的一點，Xij為交界面Sij區域中在l邊上的等效電流或者磁流，為?j區域內部由Sij面指向Sjk面的切向單位向量，同時垂直于交界線l。邊界上根據電流的連續性條件有：

圖1 多區域交界處示意圖Fig.1 The junction of multiple regions

根據上述電流連續性條件，可以充分但不必要地推出以下式子[25]：

將式(25)乘上系數后強加到式(20)和(21)，不僅可以消除式(20)和(21)中奇異的線線積分計算，還可從物理上對區域交界處的電磁流連續性進行約束，確保了DG方法與連續伽遼金方法解的一致性。

1.4 MLFMA

為了計算大尺度等離子體結構，使用MLFMA對迭代求解矩量法問題進行加速。在MLFMA中，矩陣向量乘(matrix-vector multiplication, MVM)操作可以分為近相互作用計算與遠相互作用計算：

對于近相互作用，(a,b=1,2) 可以直接使用式(17)進行計算。對于遠相互作用，則需要進行聚集、轉移、發散三個步驟，其具體的表達式為

式中：Kki為球面積分的總積分點；矩陣Dnp,ki,Tp,ki，分別為發散矩陣、轉移矩陣和聚集矩陣，

式中：為單位球面上取樣點的單位矢量；ωp為積分點權重；o與o′為測試場盒子和源場盒子的中心，則rno′,roo′rom分別為源場點到源場盒子中心、測試場盒子中心到源場盒子中心、測試場盒子中心到測試場點的距離矢量，為對應的單位矢量，rno′,roo′,rom為對應的距離模值；為第二類球面漢克爾函數；Pl為勒讓德多項式。公式(29)的精度通過累加上限lki決定，一般取lki=kkiD+2ln(kkiD+π)，D為該層立方盒子的對角線長度。最后，公式(30)的右端項可以表示為

除了同層的聚集轉移發散操作，MLFMA還涉及低層盒子到高層盒子的聚集和高層盒子到底層盒子的發散。假設MLFMA一共有Lmax層盒子，并以公式(30)中的聚集矩陣為例，底層盒子聚集到高層盒子的操作還分為兩步：一是底層盒子的平面波到高層盒子中心的移置操作，二是底層平面波積分點到高層平面波積分點的插值操作。具體的表達式為

式中：為第l層盒子中心到第(l+1)層盒子中心的距離矢量；Nl為第l層的總平面波數；為插值矩陣。從高層盒子到底層盒子發散的計算與聚集計算類似，也需要一個移置系數與逆插值操作，在此不再贅述。

對于MLFMA的遠場計算，本質上是對格林函數的梯度進行近似才得到的公式(31)，因此在L算子的相關計算中無需使用矢量恒等式將面面積分轉化為線面積分和線線積分，所以DG方法將無需在遠場計算時強加電流的連續性條件，唯一的計算區別就是邊界的半RWG只需要進行一個三角形的聚集和發散計算。

1.5 高負介電常數的截斷策略

MLFMA對高負介電常數目標的計算過程中，產生數值不穩定性的地方主要在于第二類球面漢克爾函數與移置系數[49]。當波數是純實數時，二者的特性均表現為隨距離增大的周期相位變化，當波數是非常大的虛部時，二者的特性表現為隨距離快速的增大或衰減。通過大宗量計算公式發現[50]，求解漢克爾函數過程中的兩個貝塞爾函數會隨距離指數增長，而漢克爾函數本身是隨距離指數減小的，導致在有限精度的計算機中使用貝塞爾函數求解漢克爾函數的不準確性。而對于移置矩陣來說，當移置方向和球面波方向相同時，該平面波在移置后的幅值將迅速增大，而當移置方向與平面波方向相反時，該平面波在移置后的幅值將迅速減小，增大了插值計算的誤差。通過對計算不準確的函數進一步疊加或累乘，使遠場計算的結果遠遠偏離正確的計算結果，甚至影響到整個求解矩陣的性態，導致迭代過程不收斂或計算結果不準確。

綜上所述，截斷策略的意義在于避免因介電常數過大、計算距離過大的遠場計算產生的數值不穩定問題。文獻[49]給出了截斷閾值距離d的具體推導過程，其最終結果可以表示為

式中：ε為允許產生的矩陣元素精度相對誤差上限；γ為一個與公式有關的常數，

近場計算中新的近場矩陣元素可以表示為

式中，Rc為兩個RWG中心間的距離。當進行MLFMA的遠場計算時，矩陣元素的計算公式可以表示為：

公式(35)和(36)本質上都是在幾乎不改變計算精度的前提下，約束了高負介電常數區域遠相互作用的計算。一般來說，ε=10-2即可以使計算的最終RCS結果的精度幾乎不發生改變。

2 數值算例

本章節將運用多組數值算例來展現面積分仿真技術的計算效率與計算精度。程序利用廣義最小殘差 (generalized minimum residual, GMRES) 求解器來進行仿真，Krylov 子空間的大小為 30，迭代殘差為10-3。預處理求逆使用MUMPS 5.4，計算平臺為工作站，配備有2塊Intel Xeon E7-4850 CPU和 1 TB RAM。進行對比的為合元極方法(finite-elementboundary integral, FE-BI)[51-52]，對應的計算模型為體網格剖分非均勻模型，合元極方法使用有限元-吸收邊界條件(finite-element method-absorbing boundary condition, FEM-ABC)，邊界元采用結合場積分方程(combine field integral equation, CFIE)，同樣使用MLFMA加速。

首先，考慮具有解析解的多層等離子體球算例，以驗證MR-SIE方法的準確性。如圖2所示，算例為兩層介質區域的球體，內部球體半徑0.5 m，相對介電常數為-2-0.5j，分為不連續的兩個區域；外部球體半徑1 m，相對介電常數為4-0.2j，分為不連續的4個區域。入射波為沿z方向入射x方向極化，頻率為300 MHz的平面波。外層表面使用非共形剖分將目標分成4部分，使用0.05和0.075個空氣波長進行剖分；內層表面使用非共形剖分將目標分為2部分，使用0.05和0.03個空氣波長進行剖分。利用CTF和JMCFIE兩種MR-SIE計算VV 極化方向的雙站雷達散射截面 (radar cross section, RCS)并和Mie解析解進行對比，計算結果如圖3所示。可以看到，兩種方程的計算結果與理論計算Mie解析解均吻合得很好，說明本文仿真技術可以很好地應對區域分解負介電常數介質目標的計算。

圖2 雙層介質球的區域分解形式Fig.2 Domain partitioning for the double-layer dielectric sphere

圖3 雙層介質球下三種方法得到的VV 極化的雙站RCSFig.3 The bistatic VV-polarization RCS of the double-layer dielectric sphere by 3 methods

隨后，針對有等離子體鞘套包覆的再入飛行器模型進行電磁特性分析。由美國宇航局(National Aeronautics and Space Administration, NASA)主導的無線電衰減測量計劃因其數據的可公開性，計劃中使用的再入飛行器的無線電衰減測量(Radio Attenuation Measurement, RAM)實驗模型也被大量理論研究運用。如圖4所示，RAM是一個鈍頭錐目標，頭部半徑0.152 4 m，身長1.29 m，側方半角9°。圖4還給出了通過流場計算得到的RAM模型電子密度分布，RAM飛行高度為76 km，溫度195.5 K，飛行速度為27.3馬赫?？梢钥吹剑w行器周圍的高密度等離子體主要集中在頭部，且隨著等離子體的擴散，飛行器中部到尾部的等離子體密度逐漸降低，整體呈現等離子體密度跨度大、分布非常不均勻的特點。

圖4 非均勻等離子體鞘套包覆的二維RAM飛行器模型Fig.4 2D RAM aircraft model with inhomogeneous plasma sheath

考慮單層等離子體鞘套的情況。如圖5所示，飛行器的等離子體鞘套底部厚度約為35 cm，飛行器近似為完美電導體(perfect electric conductor, PEC)，從頭部入射1 GHz的V極化平面波。不同電磁參數下等離子體鞘套的平均電子密度、碰撞頻率以及所等效的相對介電常數如表1所示。剖分尺寸為0.05個自由空間波長，面未知數總數102 030。使用FE-BI的計算結果進行對比，其面未知數只分布于鞘套表面的為31 578個，但等離子體區域的體剖分未知數多達1 519 430。對比兩者計算得到的VV極化雙站RCS如圖6所示，其中MR-SIE方法使用了JMCFIE和CTF, 而FE-BI使用CFIE?？梢钥吹絻煞N數值方法有良好的計算一致性，同時隨著電子密度的增長，RAM頭部的散射逐漸增強，RAM尾部的散射改變不明顯。再對比表1的計算性能，可以看到面積分算法的峰值內存、總計算時間均遠小于FEBI, CTF的迭代次數明顯多于JMCFIE。隨介電常數的提高，MR-SIE算法的迭代次數呈現增長趨勢，而FE-BI算法的迭代次數呈現下降趨勢。綜上，單層等離子體鞘套的散射計算體現了算法的準確性和高效性。由于CTF的收斂性較差，在之后進行更復雜的等離子體模型的計算時均采用JMCFIE以兼具準確性與高效性。

表1 單層等離子體模型不同電磁參數下的數值性能對比Tab.1 Comparison of numerical performance with different methods and parameters for the mono-layer plasma sheath

圖5 單層等離子體鞘套包覆的三維RAM飛行器模型Fig.5 RAM 3D aircraft model with mono-layer plasma sheath

圖6 單層等離子體模型下不同相對介電常數三種算法得到的VV 極化的雙站RCSFig.6 Bistatic VV-polarization RCS of the mono-layer plasma sheath model with different dielectric parameters by 3 methods

進一步考慮單層等離子體鞘套的輻射問題，如圖7和圖8所示，在飛行器尾部開槽并增加工作頻率為1.5 GHz的貼片天線，其介質板相對介電常數εr=4.4，S11如圖9所示。使用0.1個自由空間波長(20 mm)進行剖分，但由于貼片天線尺寸過薄，在側面無法使用這種粗剖分；若使用共形剖分，則貼片天線側面周圍均需要5 mm的剖分才不會出現報錯或是壞三角形；而若使用非共形剖分，天線側面使用5 mm、表面可以很方便地直接使用10 mm剖分，不僅方便剖分設計，且一定程度上減少了未知數的數量。同時可以看到，在等離子體剖分與金屬飛行器的交界處(紫色與橙色剖分交界處)、 飛行器表面與開槽的交界處(橙色剖分區域內)均可以使用非共形剖分，這種方式極大程度上提高了剖分的靈活性。最終，共形剖分的未知數數量為104 868個，非共形剖分的未知數數量為99 606個，未知數減少雖然不多，但邊界處的處理明顯變得簡單。

圖7 在尾部開槽并放置貼片天線的單層等離子體鞘套包覆的三維RAM飛行器模型Fig.7 A mono-layer plasma sheath covered RAM 3D aircraft model with a slot in the tail and a patch antenna

圖8 RAM飛行器模型尾部開槽貼片天線剖分示意圖Fig.8 The meshment at the tail of the RAM aircraft model with the patch antenna

圖9 貼片天線的S11Fig.9 The S11 parameter of the patch antenna

等離子體層的電子密度分別選取2×1010cm-3,5×1010cm-3，2×1011cm-3，其中2×1010cm-3對應的等離子體體頻率略小于入射頻率，5×1010cm-3對應的等離子體體頻率高于入射頻率。碰撞頻率統一為0.5 GHz，則三種情況下等離子體相對介電常數分別為0.3565-0.2144j，-0.609-0.536j，-5.43-2.145j。由于貼片天線存在金屬片結構，無法使用法向方程進行計算，因此計算發射的MR-SIE程序使用CTF，三種情況下等離子體鞘套包覆飛行器模型的E面和H面輻射方向圖如圖10所示。可以看出：當介電常數為負時，FEKO的計算難以收斂，而本文MRSIE程序可以正常收斂；對比不同等離子體狀態下的輻射結果，當等離子體頻率低于入射頻率時，E面的輻射功率大于0 dB，但隨著等離子體密度的增加，輻射效率明顯衰減。對比無等離子體和正介電常數等離子體的FEKO商業軟件的計算結果，表明了MRSIE計算輻射問題的準確性。

圖10 單層等離子體鞘套包覆的RAM飛行器模型的輻射方向圖Fig.10 The radiation pattern of a RAM aircraft model with mono-layer plasma sheath

考慮多層復雜等離子體鞘套的建模與仿真。建立一個五分區的等離子體鞘套模型，如圖11所示，飛行條件同圖4。從左到右每個分區在1 GHz入射頻率下的平均相對介電常數值以及厚度如表2所示，該模型的水平劃分考慮了等離子體隨x方向距離的梯度變化，水平方向的等離子體變化趨勢更為明顯，垂直方向相對介電常數的計算是通過區域內的電子密度和碰撞頻率先取平均再計算得到的。由于頭部的負介電常數過大，為簡化運算使用金屬對頭部的負介電常數等離子體進行替代。入射平面波為x方向傳播z方向極化，頻率為1 GHz，剖分尺寸為0.1個空氣波長。

表2 水平分五區域等離子體模型參數Tab.2 The model parameters of the 5-region plasma sheath split horizontally

圖11 水平分五區域等離子體三維模型示意圖Fig.11 A RAM 3D model with 5-region plasma sheath split horizontally

圖12為使用五區模型的MR-SIE和FE-BI計算結果與體剖分非均勻模型的FE-BI計算結果的對比。可以看出：MR-SIE與FE-BI計算結果吻合較好； MR-SIE與體剖分非均勻模型的FE-BI計算結果整體存在較大的誤差。因此，將模型簡單地進行5層的劃分還不足以還原非均勻流場散射的效果，且該模型還忽略了垂直于飛行器表面方向的等離子體變化。

圖12 水平分五區域等離子體模型與體剖分非均勻等離子體模型VV極化下雙站RCSFig.12 Bistatic VV-RCS of the 5-region plasma sheath model with different methods and the bistatic VV-RCS of the inhomogeneous plasma sheath model

進一步優化五區等離子體模型的分區方式，使用曲面將圖4中相近的電子密度分布范圍劃分到一起，同樣進行五個區域的剖分(如圖13所示)，從頭部到尾部五區域平均等離子體電子密度如表3所示。使用 MR-SIE計算該等離子體模型和使用FE-BI計算體剖分非均勻等離子體模型的雙站RCS對比結果如圖14所示。可以看出，相較于水平分五區的結果，使用曲面分割模型計算結果更貼近于體剖分非均勻模型的結果，但在散射值較低(<30 dB)的-20°～0°區域仍存在一定誤差。原因一方面是因為兩種算法計算精度上的誤差，RCS值越低誤差越明顯；另一方面是多區域均勻模型相較體非均勻模型存在一定誤差。計算結果顯示，水平分五區的模型相較于非均勻模型絕對平均誤差為4.59 dB，曲面分割模型相較于非均勻模型絕對平均誤差僅為1.60 dB，明顯減小。

表3 自適應切分五區域等離子體模型參數Tab.3 The model parameters of the 5-region plasma sheath model split adaptively

圖13 自適應五區域等離子體三維模型示意圖Fig.13 A RAM 3-D model with 5-region plasma sheath split adaptively

圖14 多區域均勻等離子體模型與體剖分非均勻等離子體模型VV極化下的雙站RCSFig.14 Bistatic VV-RCS of the multi-region homogeneous plasma sheath model and the inhomogeneous plasma sheath model

表4給出了不同算法不同模型的數值計算效率。從總未知數上來看：MR-SIE的未知數數量相較于FE-BI方法的未知數有量級上的差距；兩種MRSIE算法都為五分區，未知數數量相近，計算時間相近，且遠小于FE-BI的計算時間。曲面剖分雖然頭部近相互作用元素較多，占用內存較大，但其內存也遠小于FE-BI，同時在未知數和計算時間上具有更大的優勢。綜上所述，MR-SIE算法相較于體剖分計算的FE-BI在計算效率上有明顯優勢，且隨著剖分方式的細化，MR-SIE計算非均勻模型的精度有了明顯提高。

表4 不同等離子體模型仿真計算數值性能對比Tab.4 Comparison of the numerical performance with different methods and different plasma sheath models

調整自適應分區的區域數，對比計算效率與計算精度的變化。如圖15所示，分別使用自適應方法進行了1~9個區域數的分區，入射波x方向頭部入射，觀察θ=-90°,φ=0°～360°時的RCS，將FE-BI非均勻模型對應的散射結果作為參考，計算平均絕對誤差隨區域數的變化。可以看到，隨著區域數的增加，RCS的誤差明顯下降，并于5個區域時達到拐點。同時，圖15還展示了剖分未知數隨區域數的變化，發現隨著區域數量的增加，未知數是近乎線性增長的，其對應的計算效率均有所下降。因此，對于該模型來說，使用5個自適應區域進行剖分是兼顧計算效率與計算精度的選擇。

圖15 不同分區數對應的RCS絕對誤差與未知數變化Fig.15 The absolute error of RCS and unknowns numbers corresponding to different partition numbers

最后，提高入射頻率至10 GHz，計算電尺寸達50個電波長的等離子體鞘套。頻率提升后相對介電常數整體變化沒有低頻時劇烈。因此，構建由兩層不同介電常數組成的等離子體鞘套模型來研究高頻問題，如圖16所示。測試了3組不同的電子密度、碰撞頻率和介電常數(如表5所示)，模型的外層等離子體使用0.15個空氣波長進行剖分，內部使用0.1個空氣波長進行剖分，整體未知數為3 111 447。入射波x方向頭部入射，觀察角為θ=-90°,φ=0°～360° 。VV極化下的雙站RCS的計算結果如圖17所示。可以看出，等離子體存在時尾部的雙站散射明顯加強，在頭部的小范圍區域明顯減弱。從表5可以看到，每一個算例均有較好的收斂性和較高的計算效率，說明使用MR-SIE算法計算大尺度非均勻等離子體鞘套的可行性。

表5 雙層等離子體模型下不同電磁參數模型的數值性能對比Tab.5 Comparison of the numerical performance with different methods and parameters for the double-layer plasma sheath model

圖16 雙層等離子體三維模型示意圖Fig.16 A RAM 3D model with double-layer plasma sheath

圖17 雙層等離子體模型不同電磁參數下VV極化雙站RCSFig.17 Bistatic VV-RCS of the double-layer plasma sheath model with different parameters

3 結 論

本文針對非均勻、存在高負介電常數的等離子體鞘套運用基于MR-SIE的仿真策略。當等離子體的密度變化梯度不是很大時，MR-SIE可以顯著提高計算效率，可為等離子體鞘套的電磁特性進行快速評估。隨著計算區域數的提高，基于DG方法的DDM可以極大化簡網格剖分的難度，且不增加計算復雜度。通過與非均勻等離子體模型使用FE-BI算法的計算結果進行比較，MR-SIE算法的多區域等離子體模型算例結果體現了其計算精度與計算效率。在處理電大且存在高負介電常數的等離子體鞘套時，使用一種截斷策略并運用于MR-SIE算法中，本文方法可以靈活避免數值不穩定性同時確保計算精度。最終，50個波長、300萬未知數的算例展現了現有算法可觀的計算能力。從算例的應用角度出發，散射問題與輻射問題均可以使用MR-SIE進行準確快速的評估。