殷 雄 趙振維 張 厚 孫樹計 姜聿燾 許志永

（1.空軍工程大學防空反導學院，陜西 西安710051；2.中國電波傳播研究所，電波環(huán)境特性及?；夹g重點實驗室，山東 青島266107）

引 言

由于等離子體在各個現(xiàn)代化工業(yè)領域具有廣闊的應用前景，研究這類色散媒質(zhì)的電磁特性成為當前非常重要的內(nèi)容.等離子體作為人們發(fā)現(xiàn)的第四種物態(tài)，當其與電磁波相互作用時表現(xiàn)出很多不同于一般物質(zhì)的性質(zhì)，已在民用和軍事等領域得到廣泛應用，例如高速飛行器探測［1-2］、等離子體推進［3-4］、等離子體隱身技術［5-6］等等.對等離子體電磁特性的運用技術還有待進一步提高，對它的研發(fā)一直備受各國政府重視.因此，研究等離子體這類色散介質(zhì)的電磁特性，具有非常重要的意義.

時域有限差分（Finite-Difference Time-Domain，F(xiàn)DTD）技術被擴展用于等離子體的電磁模擬，涌現(xiàn)了很多算法，主要有：分段線性遞歸卷積時域有限差分（Piecewise Linear Recursive Convolution-FDTD，PLRC-FDTD）法［7］、分段線性電流密度遞歸卷積時域有限差分（Piecewise Linear Current Density Recursive Convolution-FDTD，PLCDRC-FDTD）方法［8］、Z變換-時域有限差分（ZTransformation-FDTD，ZT-FDTD）方法［9］、基于電流密度Laplace變換的FDTD方法［10］以及移位算子時域有限差分（Shift Operator-FDTD，SO-FDTD）方法［11-12］等等.這些方法的時間步長都受到Courant穩(wěn)定性條件［13］的限制，在仿真電大尺寸或高頻的等離子體時速度很慢，效率低.而無條件穩(wěn)定的交替方向隱式時域有限差分（Alternate Direction Implicit-FDTD，ADI-FDTD）算法因其計算時間步長不受Courant穩(wěn)定性條件的限制，時間步長可以取得較大而使仿真所用的時間步數(shù)大大減少，相對傳統(tǒng)FDTD算法提高了計算效率，得到了長足的發(fā)展和應用，迄今仍吸引了國際及國內(nèi)學者的廣泛關注［14-16］.一些研究者將ADI-FDTD算法擴展到色散媒質(zhì)的電磁模擬中［17-20］.文獻［19］將關于極化強度P的輔助方程引入Maxwell方程組，提出一種計算非磁化等離子體的ADI-FDTD方法.該方法在整個計算區(qū)域采用分裂場表示所有的場分量，故數(shù)據(jù)處理復雜，計算機內(nèi)存占用量很大.文獻［20］提出了一種適于計算非磁化等離子體的分段線性遞歸卷積-交替隱式-時域有限差分（Piecewise Linear JE Recursive Convolution-ADI-FDTD，PLJERC-ADI-FDTD）方法，雖然該方法精度較高，但由于涉及耗時的卷積運算，其計算效率并不可觀.

本文在Maxwell方程組的基礎上，運用中心差分格式離散非磁化等離子體所滿足的極化電流密度方程，得到一組輔助差分方程，結(jié)合交替方向隱式技術和基于輔助差分方程的完全匹配吸收層，提出了一組改進的計算非磁化等離子體散射特性的時域方法即輔助差分方程-交替隱式-時域有限差分（Auxiliary Differential Equation-ADI-FDTD，ADE-ADIFDTD）方法.該公式推導過程不涉及復雜的卷積運算，也不需將場分量進行分裂處理，因而不僅概念簡單，易于編程，而且有效地節(jié)約了內(nèi)存占用量，提高了計算效率.另外，僅需對輔助方程的相關系數(shù)進行修改，本文提出的算法可以很容易地推廣到其他普通媒質(zhì)或有耗色散媒質(zhì)的電磁模擬中，可擴展性強.通過與其他FDTD方法對比，運用算例驗證了本文方法的正確性和高效性.在此基礎上，采用本文算法分析了非磁化等離子覆蓋導體圓柱的散射特性，得出了一些有關等離子體隱身的結(jié)論.

1 改進型ADE-ADI-FDTD算法的推導

在碰撞非磁化等離子介質(zhì)中，Maxwell方程及其極化電流密度輔助方程為［10］：

式（3）中ωp，νen分別為等離子體頻率和等離子體中電子的平均碰撞頻率.

ADI-FDTD方法中的電磁場分量在空間網(wǎng)格的分布方式與傳統(tǒng)FDTD方法一致，但需要將原來一個時間步的計算分成兩次來進行.分別以兩個步驟中的Ex、Hy為例，給出基于式（1）～（3）的三維非磁化等離子體問題的改進型ADI-FDTD迭代公式.

依據(jù)公式（1）和（2），從t＝nΔt時刻到t＝（n＋1/2）Δt時刻電磁場旋度方程的離散公式可寫為（這里僅考慮Ex和Hy）

第二分步即t＝（n＋1/2）Δt時刻到t＝（n＋1）Δt時刻的電磁場差分迭代公式與第一分步類似，只是相關方程右邊場量的顯隱格式進行了交換，例如Ex的離散公式為

觀察式（4）和式（6）可知，極化電流密度J需要作特殊處理.對輔助方程（3）在t＝nΔt時刻到t＝（n＋1/2）Δt時刻采用具有二階時間和空間精度的中心差分格式進行離散，得到極化電流密度的輔助差分方程為

式中i＝x，y，z.這樣，極化電流密度的更新可直接由電場得到.從t＝（n＋1/2）Δt時刻到t＝（n＋1）Δt時刻的極化電流密度輔助差分方程具有與式（7）相同的形式，不再寫出.將式（7）代入式（4），經(jīng)適當整理得

式中：

由于式（8）兩邊具有同一時刻的場分量，不能直接用于數(shù)值計算，因此，還需進一步推導適于數(shù)值計算的迭代過程.聯(lián)立式（5）和式（8）消去 Hn＋1/2y，化簡得

式中：Ch0＝Δt/（2μ0）；Jx｜ni＋1/2，j，k的時 間 步 進 由 式（7）完成.當k沿著z方向從小到大變化時，由式（9）得到一個聯(lián)立的線性方程組，其系數(shù)矩陣是三對角帶狀矩陣，可用追趕法求解.其余電場分量可用同樣的方法得到，且均滿足形式如式（9）的方程.求出（n＋1/2）時間步的Ex后，直接解方程（5）即可求得同一時刻的磁場分量Hy.類似的，可以完成（n＋1/2）時間步全部場量的計算.采用同樣的方法，可以計算后半個時間步長的全部電磁場分量，在此不再贅述.

對于二維和一維的ADE-ADI-FDTD迭代式，只需在三維迭代式基礎上省略相應的無關場量和空間網(wǎng)格標識符k即可.以二維橫電波（Transverse Electrical Vave，TE波）第一步迭代為例，省略Ez、Hx、Hy分量和空間網(wǎng)格標識符k后得：

方程（10）是顯式的，聯(lián)合方程（11）和（12），采用前述方法可求得En＋1/2y和Hn＋1/2z，Jnx和Jny的時間步進由式（7）完成.

關于吸收邊界的配置，本文采用文獻［21］提出的輔助差分方程-完全匹配吸收層（Auxiliary Differential Equation-Perfectly Matched Layer，ADEPML），該匹配層不需要在時域分裂場量，而僅僅在PML區(qū)為每個場量添加兩個輔助變量，可以和本文提出的場量迭代公式完美結(jié)合.本文提出的ADEADI-FDTD算法在整個問題空間采用完整場變量進行迭代計算，而文獻［19］為了和PML層中分裂場迭代保持一致，在全部區(qū)域?qū)⒚總€場變量進行分裂，生成兩個子變量，因而相對文獻［19］提出的基于全分裂場的ADI技術而言，本文的算法不僅公式推導簡單，操作容易，而且節(jié)約近一半的內(nèi)存占用量.

值得注意的是，本文算法可適用于非均勻媒質(zhì)，即若ωp和νen隨空間位置而變化，那么Ca、Cb及Cc也將是空間坐標（iΔx，jΔy，kΔz）的函數(shù).從推導過程不難看出，改變極化電流密度輔助方程的相關系數(shù)，僅對Ca、Cb及Cc的取值產(chǎn)生影響，因而，本文算法可以很容易地推廣到其他普通或有耗色散介質(zhì).例如，令Ca＝1，Cb＝1，Cc＝1，Ji＝0（i＝x，y，z），本文提出的ADE-ADI-FDTD迭代公式即還原為真空條件下的常規(guī)ADI-FDTD迭代式.可見，本文提出的算法具有很好的可擴展性.方法（ZT-FDTD和SO-FDTD）中，令α＝1以滿足Courant穩(wěn)定性條件，本文提出的算法及PLJERCADI-FDTD方法中的α設為3.入射波源為微分高斯脈沖：Einc＝（t-t0）/τ×exp［-4π （t-t0）2/τ2］，其中t0＝70Δt，τ＝140Δt；整個空間分為800個網(wǎng)格，等離子體平板占據(jù)中間300～500個網(wǎng)格，其余為空氣，兩端設為7層PML吸收邊界.

圖1給出了上述幾種FDTD方法計算的電磁波通過非磁化等離子體平板的反射系數(shù)和透射系數(shù)的收斂數(shù)值解.其中，ZT-FDTD方法與SO-FDTD方法各需運行11 000時間步才能達到較好的收斂狀態(tài)；如果運行的時間步數(shù)太少，入射脈沖的時域響應將達不到穩(wěn)態(tài)收斂，以致傅里葉變換后的頻域結(jié)果不正確.為證明這一點，圖2給出了SO-FDTD方法和ZT-FDTD方法各自運行4 000時間步的透射系數(shù)計算結(jié)果與解析解的對比，由圖可知，在0～20GHz頻段內(nèi)兩種FDTD方法的數(shù)值計算結(jié)果嚴

2 算例驗證和性能比較

假設入射波沿z軸正向入射，電磁模型為“空氣＋等離子體平板＋空氣”，等離子體板厚度為1.5 cm.利用本文提出的改進型ADE-ADI-FDTD算法，文獻［20］中的PLJERC-ADI-FDTD方法，文獻［9］中的ZT-FDTD方法，文獻［11］中的SO-FDTD方法以及文獻［12］提到的解析方法對該模型進行計算.

仿真時各種參數(shù)選擇如下：νen＝20GHz，ωp＝2π×28.7×109rad/s，空間步長Δs取75μm，仿真時間步長Δt＝α·Δs/（2c），c是真空光速，常規(guī)FDTD重偏離解析解.而PLJERC-ADI-FDTD方法及本文算法由于選擇的時間步長是常規(guī)FDTD方法的3倍，各自運行4 000時間步即達良好的收斂穩(wěn)態(tài)，如圖1所示.由圖1還可知：改進型ADE-ADI-FDTD方法與PLJERC-ADI-FDTD方法的計算結(jié)果相當一致，且與解析解吻合得非常好.雖然ZT-FDTD方法與SO-FDTD方法的計算結(jié)果也與解析解吻合得較好，但其計算精度要遜于改進型ADE-ADIFDTD方法和PLJERC-ADI-FDTD方法，這從圖1（c）、（d）不難看出.

表1列出了上述幾種FDTD方法獲得穩(wěn)定收斂解的平均計算時間.表中的計算時間為各種FDTD方法在配置為Intel（R）Core（TM）i5-2400 CPU＠3.10GHz的聯(lián)想電腦上運行6次，取其算術平均值得到.其中，ZT-FDTD方法與SO-FDTD方法各運行11 000時間步；PLJERC-ADI-FDTD方法及本文算法各運行4 000時間步.由表可知，兩種ADI-FDTD方法的計算時間要小于SO-FDTD方法及ZT-FDTD方法，這是因為ADI-FDTD方法擺脫了Courant穩(wěn)定性條件的束縛，時間步長可以是常規(guī)FDTD方法的數(shù)倍，因此，仿真的時間步數(shù)可以大大降低，計算效率隨之提高.數(shù)值實驗表明，如果進一步增大α值，ADI-FDTD方法的仿真步數(shù)將進一步減少，計算效率將進一步提高.但是，時間步長過大時會存在較大的數(shù)值色散誤差和各向異性誤差，導致S參數(shù)出現(xiàn)波動.本文實驗表明，α的取值一般不應超過6.另外，本文提出的ADE-ADI-FDTD算法由于避免了復雜的卷積運算，其計算效率要高于PLJERC-ADI-FDTD方法.

表1 不同F(xiàn)DTD方法計算電磁波通過等離子體平板的平均計算時間

3 非磁化等離子體覆蓋導體圓柱的散射特性

研究電磁波在等離子體中的電磁散射特性，能夠給我們合理利用等離子體隱身技術提供理論依據(jù).利用本文提出的算法研究了非磁化等離子體覆蓋導體圓柱對入射橫磁波（Transverse Magnetic Wave，TM波）（僅有Ez，Hx和Hy分量）和TE波（僅有Hz，Ex和Ey分量）的散射特性.等離子體覆蓋導體圓柱的FDTD計算區(qū)域及各種邊界劃分如圖3所示.假設入射波頻率為15GHz，導體圓柱的半徑r為入射波波長λ0的一半即1cm，等離子體涂覆層厚度h分別為0.5cm和1cm，等離子體頻率為12GHz，碰撞頻率νen分別取50GHz和100 GHz.仿真時選擇時間步長控制參數(shù)α＝3，空間網(wǎng)格步長Δx＝Δy＝λ0/40，整個問題空間劃分為150×150個網(wǎng)格.利用改進的ADE-ADI-FDTD方法計算得到的雙站雷達散射截面（Radar Cross Section，RCS）如圖4所示（圖中導體曲線是指半徑為1 cm的導體圓柱對入射波的散射結(jié)果，厚度是指等離子體涂覆層厚度）.

圖3 等離子體覆蓋導體圓柱的FDTD模型

觀察圖4可知，對于等離子體圓柱，TE波與TM波散射規(guī)律基本一致：與純導體圓柱的散射相比，非磁化等離子體涂覆層能夠較大地減小后向和大散射角范圍內(nèi)的RCS，但對前向散射和小散射角范圍的散射卻有小幅增強的作用.對圓柱型非磁化等離子體包層來說，等離子體對電磁波的作用主要有折射效應和吸收衰減作用，這是導致后向RCS減小的原因.對于前向附近RCS較大的解釋是，就規(guī)則目標來說，前向RCS與目標的投影面積成正比關系，若把等離子體包層和導體看作一個整體，其投影面積確實增大了［22］；另外，由于入射波頻率高于等離子體頻率，電磁波可以以行波的形式到達圓柱的陰影區(qū)［22-23］.

由圖4還可以看出，當?shù)入x子體涂覆層厚度增大時，等離子體覆蓋導體的RCS在前向附近增強了，這是由于其整體投影面積增大的緣故；而當散射角較大時，其RCS的下降非常明顯，反映出等離子體的折射和吸收作用.在相同的厚度下，等離子體碰撞頻率較大時，其RCS總體來看是降低了，但在TE波散射的個別角度范圍這個現(xiàn)象是相反的.

4 結(jié) 論

提出了一種改進的計算非磁化等離子體散射特性的無條件穩(wěn)定時域方法即輔助差分方程-交替隱式-時域有限差分（ADE-ADI-FDTD）方法.在此方法中，運用中心差分格式離散極化電流密度和等離子體特征參數(shù)所滿足的關系方程，得到一組輔助差分方程，避免了復雜、耗時的卷積運算（與文獻［20］的PLJERC-ADI-FDTD相比），提高了計算效率.采用ADE-PML，使得本文算法在整個電磁問題空間都不要分裂場變量，相對文獻［19］提出的ADI方法而言，不僅操作簡單，而且節(jié)約了近一半的內(nèi)存空間.算例的對比計算表明，本文提出的算法不論在精度上還是在效率上都高于普通的Z變換FDTD方法和SO-FDTD方法，在和文獻［20］提出的PLJERCADI-FDTD方法保持相同精度的同時取得了計算速度上的優(yōu)勢.通過對輔助方程的關鍵系數(shù)進行修改，本文提出的改進型無條件穩(wěn)定ADI-FDTD算法可以很容易地推廣到其他非均勻有耗色散媒質(zhì)的電磁模擬中，具有很好的普適性.運用本文算法對非磁化等離子體涂覆導體圓柱的散射特性的分析表明：在相當寬的散射角范圍內(nèi)，等離子體涂覆導體圓柱的RCS與導體圓柱的RCS相比要小很多，特別是大散射角時，這種現(xiàn)象更明顯；保持碰撞頻率不變時，增加非磁化等離子體涂覆層厚度能夠較大地減小后向和大散射角范圍內(nèi)的RCS；而在相同的涂覆層厚度下，等離子體碰撞頻率變大時，其RCS總體上呈降低趨勢，但在TE波散射的個別角度范圍這個變化是相反的.

［1］KIM M，KEIDAR M，BOVD I D.Two-dimensional model of an electromagnetic layer for the mitigation of communications blackout［C］//47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition.Florida，2009：10-17.

［2］MANNING R M.Analysis of Electromagnetic Wave Propagation in a Magnetized Re-Entry Plasma Sheath via the Kinetic Equation［R］.Cleveland：National Aeronautics and Space Administration Glenn Research Center，2009：1-33.

［3］SHANG J S.Plasma injection for hypersonic bluntbody drag reduction［J］.AIAA Journal，2002，40（6）：1178-1186.

［4］羅金玲，徐 敏，戴梧葉，等.高速飛行器等離子體減阻的數(shù)值模擬研究［J］.宇航學報，2009，30（1）：119-122.LUO Jinling，XU Ming，DAI Wuye，et al.Numerical simulation investigation on plasma injection for drag reduction of hypersonic vehicle［J］.Journal of Astronautics，2009，30（1）：119-122.（in Chinese）

［5］BESKAR，CHRISTOPHER.R.Cold Plasma Cavity Active Stealth Technology［R］.USA：STAVATTI White Paper，2004.

［6］潘文俊，童創(chuàng)明，周 明.等離子體與等離子體隱身技術［J］.電訊技術，2009，49（8）：108-112.PAN Wenjun，TONG Chuangming，ZHOU Ming.Plasma and plasma stealth technology［J］.Telecommunication Engineering，2009，49（8）：108-112.（in Chinese）

［7］AI X，HAN Y P，LI C Y.Analysis of dispersion relation of piecewise linear recursive convolution FDTD method for space-varying plasma［J］.Progress in Electromagnetics Research Letters，2011，22：83-93.

［8］LIU S，LIU M，HONG W.Modified piecewise linear current density recursive convolution finite-difference time-domain method for anisotropic magnetized plasma［J］.IET Microw Antennas Propag，2008，2（7）：677-685.

［9］晏 明，許 金，余錫文.不均勻非磁化等離子體柱的ZTFDTD分析［J］.海軍工程大學學報，2009，21（3）：18-22.YANG Ming，XU Jin，YU Xiwen.ZTFDTD analysis of inhomogeneous unmagnetized plasma cylinder［J］.Journal of Naval University of Engineering，2009，21（3）：18-22.（in Chinese）

［10］YANG Lixia，XIE Yingtao，YU Pingping.Study of bandgap characteristics of 2Dmagnetoplasma photonic crystal by using M-FDTD method［J］.Microwave and Optical Technology Letters，2011，53（8）：1778-1784.

［11］YANG Hongwei.A FDTD analysis on magnetized plasma of Epstein distribution and reflection calculation［J］.Computer Physics Communications，2009，180（1）：55-60.

［12］YIN Xiong，ZHANG Hou，XU Haiyang and ZENG Xianfeng.Improved shift-operator FDTD method for anisotropic magnetized plasma with arbitrary magnetic diclination［J］.Progress in Electromagnetics Research B，2012，38：39-56.

［13］TAFLOVE A，HAGNESS S C.Computational Electrodynamics：the Finite-Difference Time-Domain Method［M］.Norwood：Artech House，2005：302-313.

［14］FU W，TAN E L.Stability and dispersion analysis for ADI-FDTD method in lossy media［J］.IEEE Transactions on Antennas and Propagation，2007，55（4）：1095-1102.

［15］RAMADAN O.An Implicit 4-stage ADI waveequation PML algorithm for 2-D FDTD simulations［J］.IEEE Antennas and Wireless Propagation Letters，2009，8：391-393.

［16］王全民，陳 彬，郭 剛，等.超寬帶沖激無線電引線地面回波仿真算法［J］.系統(tǒng)仿真學報，2011，23（3）：469-473.WANG Quanmin，CHEN Bin，GUO Gang，et al.Ground echo simulation algorithm for ultra-wideband impulse-radio fuze［J］.Journal of System Simulation，2011，23（3）：469-473.（in Chinese）

［17］GARCIA S G，RUBIO R G，BRETONES A R，et al.Extension of the ADI-FDTD method to Debye media［J］.IEEE Transactions on Antennas and Propagation，2003，51（11）：3183-3186.

［18］PEREDA J A，GONZALEZ O，GRANDE A，et al.An alternating-direction implicit FDTD modeling of dispersive media without constitude relation splitting［J］.IEEE Microwave and Wireless Components Letters，2008，18（11）：719-721.

［19］湯 煒，胡茂兵.輔助方程-雙向隱式差分法的電磁散射研究［J］.電波科學學報，2011，26（5）：904-909.TANG Wei，HU Maobing.Electromagnetic scattering using ADEs-ADI-FDTD method［J］.Chinese Journal of Radio Science，2011，26（5）：904-909.（in Chinese）

［20］常 雨.超聲速/高超聲速等離子體流場數(shù)值模擬及其電磁特性研究［D］.長沙：國防科學技術大學，2009：82-93.CHANG Yu.Numerical Research on Supersonic/Hypersonic Plasma Flow and its Electromagnetic Characteristics［D］.Changsha：National University of Defense Technology，2009：82-93.（in Chinese）

［21］李 毅.雷達隱身目標電磁散射計算與實驗研究［D］.長沙：國防科學技術大學，2007：43-47.LI Yi.Study on Electromagnetic Scattering Computation and Experimentation of Radar Stealth Targets［D］.Changsha：National University of Defense Technology，2007：43-47.（in Chinese）

［22］錢志華.等離子體天線的輻射與散射特性分析［D］.南京：南京理工大學，2006：69-101.QIAN Zhihua.Analysis of Radiation and Scattering Characteristics of Plasma Antenna［D］.Nanjing：Nanjing University of Science ＆ Technology，2006：69-101.（in Chinese）

［23］葛德彪，吳躍麗，朱湘琴.等離子體散射FDTD分析的移位算子方法［J］.電波科學學報，2003，18（4）：359-362.GE Debiao，WU Yueli，ZHU Xiangqin.Shift operator method applied for dispersive medium in FDTD analysis［J］.Chinese Journal of Radio Science，2003，18（4）：359-362.（in Chinese）