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

湯 煒 胡茂兵

（華僑大學(xué)信息學(xué)院，福建 廈門361021）

1.引 言

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

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

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

2.ADEs方法簡述［15］

為了后面闡述的方便，在此對(duì)ADEs方法進(jìn)行簡單描述。一階Debye媒質(zhì)的極化強(qiáng)度P與電場(chǎng)E的關(guān)系可以描述為

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

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

同時(shí)根據(jù)電場(chǎng)旋度方程，即

方程（1）～（3）即構(gòu)成了ADEs方法的核心方程，相比傳統(tǒng)的FDTD方法，增加極化強(qiáng)度P的迭代方程，故稱之為輔助方程方法。

輔助方程中極化強(qiáng)度P的時(shí)間和空間抽樣與電場(chǎng)相同，兩者具有空間位置與時(shí)刻相同的特點(diǎn)，故利用方程（1）即可完成對(duì)P的時(shí)間步進(jìn)。

3.ADE－ADI－FDTD迭代方程的推導(dǎo)

3.1 全分裂場(chǎng)的ADI－FDTD簡述

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

圖1 傳統(tǒng)ADI－FDTD方法Ex迭代涉及場(chǎng)量及空間分布

本文采用分裂場(chǎng)方法［12］以避免出現(xiàn)上述情況，也即全部區(qū)域都如PML層中，將常規(guī)場(chǎng)量進(jìn)行分裂。以Ex及Px為例，其分裂場(chǎng)及方程分別為

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

相比傳統(tǒng)FDTD方法，本文邊界條件是利用兩個(gè)分裂場(chǎng)量實(shí)現(xiàn)的。假設(shè)計(jì)算區(qū)域的總場(chǎng)區(qū)域限定在范圍內(nèi)，F(xiàn)DTD方法中根據(jù)迭代方程所涉及場(chǎng)量分布，在跨越以下散射場(chǎng)－總場(chǎng)邊界時(shí)需要考慮連接邊界條件：

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

3.2 ADEs－ADI－FDTD方法

3.2.1 極化強(qiáng)度場(chǎng)的迭代

在 ADEs－ADI－FDTD 中，極 化 強(qiáng) 度 也 進(jìn) 行分裂，可得

方程（1）中場(chǎng)量均只出現(xiàn)對(duì)時(shí)間的導(dǎo)數(shù)，可假設(shè)極化強(qiáng)度與電場(chǎng)在時(shí)間空間的抽樣完全相同，進(jìn)一步可得電場(chǎng)與極化強(qiáng)度P的步進(jìn)表達(dá)式。ADI－FDTD方法中場(chǎng)量時(shí)間步進(jìn)是以0.5時(shí)間步迭代進(jìn)行的，以n→n＋0.5為例可得

極化強(qiáng)度n＋0.5→n＋1的迭代也有相同的系數(shù)和類似的表達(dá)式。

3.2.2 電場(chǎng)與磁場(chǎng)在第一時(shí)間分步的迭代

限于篇幅，本文僅僅以相關(guān)分量在第一時(shí)間分步迭代為例來推導(dǎo)該方法的步進(jìn)方程，其他場(chǎng)量均可類推。提取方程（2）中x分量并將其改寫為分裂場(chǎng)形式。

按照第一時(shí)間分步的差分特點(diǎn)，方程（9）和（10）可分別寫成如下步進(jìn)方程

為便于突出重點(diǎn)，以（＊）縮寫形式表示與方程左邊場(chǎng)量空間下標(biāo)相同，下同。式中S，T為關(guān)于媒質(zhì)、時(shí)空步長有關(guān)的參數(shù)。

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

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

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

式中：

雖然方程（15）的左邊為三個(gè)后時(shí)刻的電場(chǎng)量，但該方程可以通過追趕法進(jìn)行求解，從而完成一列電場(chǎng)的步進(jìn)。

第一與第二時(shí)間分步中間接迭代場(chǎng)量都可以按照類似的方法完成場(chǎng)量的步進(jìn)，進(jìn)而完成所有場(chǎng)量的更新。

3.2.3 關(guān)于時(shí)間步長的選擇

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

4.色散媒質(zhì)的電磁仿真

為了驗(yàn)證上述公式推導(dǎo)的正確性，體現(xiàn)本文方法的優(yōu)勢(shì)，本文通過三維算例對(duì)電磁波與色散媒質(zhì)的作用過程進(jìn)行模擬。在以下的計(jì)算中時(shí)空步長的選擇依據(jù)如下原則

式中：λr，λi分別為最高頻率對(duì)應(yīng)實(shí)波長與虛波長；α為時(shí)空約束常數(shù)，F(xiàn)DTD計(jì)算時(shí)為0.5，對(duì)該參數(shù)的設(shè)定可得不同的時(shí)間步長進(jìn)行計(jì)算。

4.1 水泡的后向散射截面

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

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

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

表1 傳統(tǒng)方法與本文方法計(jì)算時(shí)間比較

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

各向同性等離子體屬于常見的色散媒質(zhì)，通過參數(shù)之間的變換可以將其轉(zhuǎn)變?yōu)榫哂袑?dǎo)電損耗Debye媒質(zhì)的表達(dá)形式。本文在此計(jì)算Von Karman型金屬體覆蓋等離子體時(shí)的截面。計(jì)算模型如圖3所示。其中D1＝0.1m，D2＝0.2m，L1＝0.2m，L2＝0.4m.電磁波迎頭入射，F(xiàn)DTD計(jì)算時(shí)取離散間隔δ＝3.33mm.等離子體的濃度N0＝5.0×1017m－3，電子碰撞頻率ωeff＝2.0×1010Hz.圖4是后向雷達(dá)散射截面，分別給出了未覆蓋等離子體時(shí)和有等離子體覆蓋時(shí)的計(jì)算結(jié)果。由圖可見等離子體覆蓋層在3～7GHz頻段里減小了后向RCS.

需要說明的是，算例由于受到頻譜較寬和計(jì)算機(jī)內(nèi)存的限制，空間步長沒有取得足夠細(xì)，因此當(dāng)約束參數(shù)為4和5時(shí)，計(jì)算結(jié)果精度不是很好，論文在此設(shè)置該參數(shù)為2，計(jì)算結(jié)果與FDTD計(jì)算結(jié)果吻合較好?？梢哉f明該方法在計(jì)算較為復(fù)雜的目標(biāo)時(shí)，仍然是可行及可靠的。

5.結(jié) 論

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

［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］劉 波，高本慶，薛正輝，等.無條件穩(wěn)定的交替方向隱式FDTD算法［J］.電波科學(xué)學(xué)報(bào)，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結(jié)構(gòu)中的應(yīng)用［J］.電波科學(xué)學(xué)報(bào)，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在二維散射問題中的應(yīng)用［J］.電波科學(xué)學(xué)報(bào)，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在波導(dǎo)電路分析中的應(yīng)用［J］.電子學(xué)報(bào)，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］鄭奎松，葛德彪.總場(chǎng)－散射場(chǎng)方法在二維 ADIFDTD中的實(shí)現(xiàn)［J］.西安電子科技大學(xué)學(xué)報(bào)，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］湯 煒，閆玉波，李清亮，等.一種新時(shí)域交替隱式差分算法在散射問題中的應(yīng)用［J］.物理學(xué)報(bào)，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］魏 兵.各向異性介質(zhì)電磁散射及參數(shù)反演研究［D］.西安電子科技大學(xué)，2003.WEI Bin.Study of Electromagnetic Scattering and Reconstruction for Anisotropic Media［D］.Xidian U－niversity，2003.（in Chinese）