輔助方程－雙向隱式差分法的電磁散射研究

湯 煒 胡茂兵

（華僑大學信息學院，福建 廈門361021）

1.引 言

眾所周知，時域有限差分法（FDTD）［1］是電磁學領域一種非常有效的數值計算方法，已廣泛應用于電磁散射、天線、電磁兼容、生物電磁場［1］以及電波傳播等問題的計算與模擬。但由于FDTD是顯式差分方法，其時空步長要滿足穩定性條件要求，嚴格的時間空間步長導致計算效率較低。隨后學者們提出偽譜時域差分（PSTD），時域多分辨小波（MRTD）用以改進傳統方法。這兩種方法對空間步長的選擇有一定的改善效率但低下的問題并沒有得到完全解決。

日本學者 Namiki［2－3］將交替隱式差分格式應用到FDTD技術中，提出了ADI－FDTD算法，這一方法的提出迄今仍吸引了國際及國內學者的廣泛關注［4－6］。國內 方 面，高 本 慶［7］，湯 煒［8－9］，呂 善 偉［10］，葛德彪［11］等人都對該方法進行了一定的研究。但是遺憾的是，由于方程的復雜性使得采用ADIFDTD方法進行三維目標計算論文相對較少。文獻［12］提出在整個計算區域都采用分裂場使ADIFDTD方法得到一定的簡化，并用該方法計算了三維目標的散射問題。而本文則是在原有的基礎上，進一步完成色散媒質的電磁仿真。

關于色散媒質的研究無論在軍用還是民用都具有重要意義，如電磁屏蔽，信道衰減等。FDTD技術中色散媒質的電磁模擬大體有以下幾種方法，依次為分段線性遞推卷積方法 （PLRC）［13］，Z 變 換 方法［14］，輔助差分方程方法（ADEs）［15］，電流密度卷積［16］。幾種方法中 ADEs直接從最基本的 Maxwell方程得出，相對簡單而且物理意義明確直觀。本文以該方法為基礎，推導了ADEs－ADI－FDTD方法的迭代方程，并利用算例對這一方法進行了數值驗證及時間上的對比。

2.ADEs方法簡述［15］

為了后面闡述的方便，在此對ADEs方法進行簡單描述。一階Debye媒質的極化強度P與電場E的關系可以描述為

式中：εs和ε∞分別為0頻及∞頻時媒質相對介電常數；τ為媒質分子的馳豫時間。將上式寫為時域形式，有

利用D（ω）＝ε0E（ω）＋P（ω），并將其帶入磁場旋度方程，可得

同時根據電場旋度方程，即

方程（1）～（3）即構成了ADEs方法的核心方程，相比傳統的FDTD方法，增加極化強度P的迭代方程，故稱之為輔助方程方法。

輔助方程中極化強度P的時間和空間抽樣與電場相同，兩者具有空間位置與時刻相同的特點，故利用方程（1）即可完成對P的時間步進。

3.ADE－ADI－FDTD迭代方程的推導

3.1 全分裂場的ADI－FDTD簡述

一維和二維ADI－FDTD方法由于場量較少故相對簡單，而三維問題相對復雜，同時散射問題中存在吸收邊界與連接邊界。圖1簡單顯示該方法中迭代所涉及場量及其空間分布，其中包含四個元胞中的11個場量。當元胞跨越吸收邊界、連接邊界時，需對各場量的跨越作出大量的判斷，實現上出現極大困難。

圖1 傳統ADI－FDTD方法Ex迭代涉及場量及空間分布

本文采用分裂場方法［12］以避免出現上述情況，也即全部區域都如PML層中，將常規場量進行分裂。以Ex及Px為例，其分裂場及方程分別為

進而對微分方程進行ADI技術處理。PML層中虛擬的是各向異性媒質，而各向同性媒質可作特殊各向異性媒質處理，故在空間預先設定好各元胞電磁參數，即可消除PML吸收層與散射場區邊界處的跨越問題。

相比傳統FDTD方法，本文邊界條件是利用兩個分裂場量實現的。假設計算區域的總場區域限定在范圍內，FDTD方法中根據迭代方程所涉及場量分布，在跨越以下散射場－總場邊界時需要考慮連接邊界條件：

而分裂場中，方程（5）涉及場量為沿著y軸分布的三個場量，僅在邊界1），2）時需要考慮邊界條件；與之對應方程（6）涉及場量為沿著z軸分布的三個場量，僅在邊界3），4）時需要考慮邊界條件。綜合兩者可以看出分裂場方法與傳統方法在邊界條件上是一致的。

3.2 ADEs－ADI－FDTD方法

3.2.1 極化強度場的迭代

在 ADEs－ADI－FDTD 中，極 化 強 度 也 進 行分裂，可得

方程（1）中場量均只出現對時間的導數，可假設極化強度與電場在時間空間的抽樣完全相同，進一步可得電場與極化強度P的步進表達式。ADI－FDTD方法中場量時間步進是以0.5時間步迭代進行的，以n→n＋0.5為例可得

極化強度n＋0.5→n＋1的迭代也有相同的系數和類似的表達式。

3.2.2 電場與磁場在第一時間分步的迭代

限于篇幅，本文僅僅以相關分量在第一時間分步迭代為例來推導該方法的步進方程，其他場量均可類推。提取方程（2）中x分量并將其改寫為分裂場形式。

按照第一時間分步的差分特點，方程（9）和（10）可分別寫成如下步進方程

為便于突出重點，以（＊）縮寫形式表示與方程左邊場量空間下標相同，下同。式中S，T為關于媒質、時空步長有關的參數。

考察方程（12），可以注意到方程的右邊均為時刻場量，該方程的表達形式與傳統FDTD方法類似，可以直接進行場量更新。在ADEs－ADI－FDTD方法中，能夠直接更新的場量有第一時間分步中的Exy，Eyz，Ezx，Hxz，Hyx，Hzy及第二時間分步中的Exz，Eyx，Ezy，Hxy，Hyz，Hzx，這些場量均可以在各自時間分步中提前完成并處理好連接邊界條件。

而方程（13）右邊磁場分量均為（n＋0.5）時刻，不能直接進行迭代。但是幸運的是分裂場Hyx可以預先完成更新，因此它雖與方程左邊為相同時刻場量，但為已知量不需考慮處理。故僅剩下Hyz需要處理，按照ADI－FDTD方法，可寫出這一時間分步中Ηyz迭代方程：

與前面類似，以（＃）表示與方程左邊場量的空間位置在x，y方向相同，下同。方程（13）右邊Exy為直接迭代場量，故未知量僅為Exz.根據式（14）可以寫出的方程，并與式（13）代入（15）整理可得：

式中：

雖然方程（15）的左邊為三個后時刻的電場量，但該方程可以通過追趕法進行求解，從而完成一列電場的步進。

第一與第二時間分步中間接迭代場量都可以按照類似的方法完成場量的步進，進而完成所有場量的更新。

3.2.3 關于時間步長的選擇

本文方法中計算空間采用分裂場形式進行電磁仿真，在一定程度上增加了內存的占用，但是另一方面由于時間步長的無約束性（在實際計算過程中時間步長不能設置太大。ADI技術僅保證迭代過程中不產生數值發散現象，并不能保證計算正確。對于時間的抽樣仍需要滿足對入射波信號的抽樣定律，最小抽樣頻率為所關心最高頻率的2倍，另抽樣間隔較大，會導致計算結果出現Gibbs現象，降低計算精度），時間步長可以達到傳統FDTD技術的10倍以下，時間步長的增加也即意味著仿真過程中時間步數減少，提高了FDTD的計算效率。

4.色散媒質的電磁仿真

為了驗證上述公式推導的正確性，體現本文方法的優勢，本文通過三維算例對電磁波與色散媒質的作用過程進行模擬。在以下的計算中時空步長的選擇依據如下原則

式中：λr，λi分別為最高頻率對應實波長與虛波長；α為時空約束常數，FDTD計算時為0.5，對該參數的設定可得不同的時間步長進行計算。

4.1 水泡的后向散射截面

水作為一種常見媒質，測試結果顯示零頻介電常數εs＝81，無窮頻率介電常數ε∞＝1.8，馳豫時間τ＝9.4×10－12，本文所示的雙層球為一個水泡。水泡外徑和內徑分別為0.158mm和0.105mm.受到實波數和虛波數的限制，空間步長Δs＝4μm，其后向散射系數的頻響特性如圖2所示。

圖2 水泡后向散射的頻響特性

對該模型先后采用解析方法，ADEs－FDTD和ADEs－ADI－FDTD等三種方法進行了模擬。其中ADEs－ADI－FDTD時間步長根據取值，所得結果與理論結果吻合很好，驗證了色散媒質ADI－FDTD方法的計算精度和有效性。從算例尺寸來看，由于在最高頻率目標的電尺寸仍然處于瑞利區，其單站雷達散射截面呈單調上升區間。針對上面算例，本文對FDTD與ADI－FDTD的時間步長、時間步數以及程序運行時間進行比較，其結果如表1所示。

表1 傳統方法與本文方法計算時間比較

4.2 Von Karman型金屬體覆蓋非磁性等離子體時的后向散射［17］

各向同性等離子體屬于常見的色散媒質，通過參數之間的變換可以將其轉變為具有導電損耗Debye媒質的表達形式。本文在此計算Von Karman型金屬體覆蓋等離子體時的截面。計算模型如圖3所示。其中D1＝0.1m，D2＝0.2m，L1＝0.2m，L2＝0.4m.電磁波迎頭入射，FDTD計算時取離散間隔δ＝3.33mm.等離子體的濃度N0＝5.0×1017m－3，電子碰撞頻率ωeff＝2.0×1010Hz.圖4是后向雷達散射截面，分別給出了未覆蓋等離子體時和有等離子體覆蓋時的計算結果。由圖可見等離子體覆蓋層在3～7GHz頻段里減小了后向RCS.

需要說明的是，算例由于受到頻譜較寬和計算機內存的限制，空間步長沒有取得足夠細，因此當約束參數為4和5時，計算結果精度不是很好，論文在此設置該參數為2，計算結果與FDTD計算結果吻合較好。可以說明該方法在計算較為復雜的目標時，仍然是可行及可靠的。

5.結 論

用ADEs－ADI－FDTD方程對色散媒質與電磁波的作用進行了仿真，計算結果表明，該方法能夠有效模擬這一過程，計算精度比較高。時間步長的增加減少了迭代次數，使計算效率有極大提高。三維算例表明ADEs－ADI－FDTD方法的耗時僅僅為傳統FDTD的1／6，然而需要說明的是，ADI系列方法是時間無條件穩定性方法，也即是說，無論時間步長設置多大，模擬的結果都不會發散，但是并不能保證正確。實際上時間對波形的采樣需要滿足抽樣定理，即時間步長必須要小于關心頻段中的最小周期的一半。考慮到計算精度的影響，固定空間步長的情況下，約束條件參數α的取值介于2～5之間較為合適。在存在強反射（散射）的條件下，參數的取值不宜取得太大，而對于弱反射（散射）時，參數的取值則較為寬松，計算效率也將更高。

［1］TAFLOVE A.Computational Electrodynamics：The Finite Difference Time Domain method［M］.Boston：Artech House，2000.

［2］NAMIKI T.A new FDTD algorithm based on alternating－direction implicit method［J］.IEEE Trans.on MTT.，1999，47（10）：2003－2007.

［3］NAMIKI T.3－D ADI－FDTD method unconditionally stable time－domain algorithm for solving full vector Maxwell＇s equations［J］.IEEE Trans.on AP.，2000，48（10）：1743－1748.

［4］ROUF H K，COSTEN F，and GARCIA S G.Reduction of numberical errors in frequency dependent ADIFDTD［J］.Electronics Letters，2010，46（7）：489－490.

［5］SINGH V，AJEET A，KWATRA N，et al.Computation of induced current densities in the human body at low frequencies due to contact electrodes using the ADI－FDTD method［J］.IEEE Trans.EMC，2010，52（3）：537－544.

［6］Pereda J A，Grande A，Gonza？lez O，et al.The ADI－FDTD method for transverse magnetic waves in conductive materials［J］.IEEE Trans.AP.，2010，58（8）：2790－2793.

［7］劉 波，高本慶，薛正輝，等.無條件穩定的交替方向隱式FDTD算法［J］.電波科學學報，2002，17（2）：437－440.LIU Bo，GAO Benqing，XUE Zhenghui，et al.Numerical simulation of unc0nditi0nany stable ADI＿FDTD algorithm［J］.Chinese Journal of Radio Science，2002，17（2）：437－440.（in Chinese）

［8］湯 煒，焦培南，李清亮，等.ADI－FDTD方法在一維PBG結構中的應用［J］.電波科學學報，2003，18（3）：281－285 TANG Wei，JIAO Peinan，LI Qingliang，et al.Application ADI－FDTD in one－dimension PBG structure［J］.Chinese Journal of Radio Science，2003，18（3）：281－285.（in Chinese）

［9］湯 煒，李清亮，焦培南，等.ADI－FDTD在二維散射問題中的應用［J］.電波科學學報，2003，18（6）：620－624.TANG Wei，LI Qingliang，JIAO Peinan，et al.Twodimension scattering analysis using ADI－FDTD method［J］.Chinese Journal of Radio Science，2003，18（6）：620－624.（in Chinese）

［10］張 巖，呂善偉.ADI－FDTD＋GRT在波導電路分析中的應用［J］.電子學報，2005，33（9）：1688－1690.ZHANG Yan，LV Shanwei.Analysis of waveguide circuits using the ADI－FDTD＋GRT method［J］.Acta Electronica Sinica，2005，33（9）：1688－1690.（in Chinese）

［11］鄭奎松，葛德彪.總場－散射場方法在二維 ADIFDTD中的實現［J］.西安電子科技大學學報，2006，33（2）：201－210.ZHENG Kuisong，GE Debiao.Imaplementation of the total－field／scattered－field technique in the 2D ADI－FDTD［J］.Journal Xidian University，2006，33（2）：201－210（in Chinese）

［12］湯 煒，閆玉波，李清亮，等.一種新時域交替隱式差分算法在散射問題中的應用［J］.物理學報，2004，53（12）：4173－4179.TANG Wei，YAN Yubo，LI Qingliang，et al.A new ADI－FDTD analysis in electromagnetic scattering［J］.Acta Physica Sinica，2004，53（12）：4173－4179.（in Chinese）

［13］KELLEY D F and LUEBBERS R J.Piecewise linear recursive convolution for dispersive media using FDTD［J］.IEEE Trans.on AP，1993，41（11）：1249－1257.

［14］SULLIVAN D M.Frequency－dependent FDTD methods using Z transform［J］.IEEE Trans.on AP，1990，40（10）：1223－1230.

［15］OKONIEWSKI M，MROZOWSKI M and STUDHLY M A.Simple treatment of multi－term dispersion in FDTD［J］.IEEE Microwave Guided Wave Lett.，1997，7（5）：121－123.

［16］KELLEY D F，LUEBBERS R J.Piecewise linear recursive convolution for dispersive media using FDTD［J］.IEEE.Trans.AP，1996，44（6）：792－797.

［17］魏 兵.各向異性介質電磁散射及參數反演研究［D］.西安電子科技大學，2003.WEI Bin.Study of Electromagnetic Scattering and Reconstruction for Anisotropic Media［D］.Xidian U－niversity，2003.（in Chinese）