李洪宇，王發(fā)年，劉宗信，2，李林，李軍

（1.解放軍95972部隊(duì)甘肅酒泉735018；2.解放軍理工大學(xué)江蘇南京210007）

在電子元器件中，許多結(jié)構(gòu)都具有一維或二維周期性，如光柵、相控天線陣、光電子帶隙（PBG）及頻率選擇表面（FSS）等等。時(shí)域有限差分法由于自身的一些優(yōu)點(diǎn)，如引入周期邊界條件（PEC）就可以通過(guò)研究一個(gè)單元的電磁散射特性來(lái)獲取整個(gè)結(jié)構(gòu)的電磁散射特性，在解決此類問(wèn)題時(shí)具有很大的優(yōu)勢(shì)，被廣泛用來(lái)模擬這些器件的電磁散射特性。但是，傳統(tǒng)時(shí)域有限差分法受Courant-Friedrich-Lewy（CFL）穩(wěn)定性條件的限制，在處理一些細(xì)小結(jié)構(gòu)如薄層介質(zhì)、孔、縫等上有局限性。為了消除CFL穩(wěn)定性條件的限制，人們努力尋找各種無(wú)條件穩(wěn)定算法或弱條件穩(wěn)定算法。1999年，T.Namiki首先將交變隱式差分方法（Alternating-Direction Implicit Method，ADI）應(yīng)用到FDTD中，提出了無(wú)條件穩(wěn)定ADIFDTD方法[1]，該算法在分析細(xì)微結(jié)構(gòu)電磁問(wèn)題時(shí)，縮短了計(jì)算時(shí)間，提高了計(jì)算效率。由于ADI-FDTD突破了穩(wěn)定性條件的限制，時(shí)間步長(zhǎng)可以適當(dāng)增加，使計(jì)算時(shí)間大幅減少，因此得到了較為廣泛的應(yīng)用。但是ADI-FDTD法也存在著不足，當(dāng)時(shí)間步長(zhǎng)增大時(shí)，將會(huì)導(dǎo)致較大的數(shù)值誤差，降低了計(jì)算精度。2003年，Binke Huang等人提出了一種弱條件穩(wěn)定（Weakly Conditionally Stable）的FDTD算法[2]，該算法的的時(shí)間步長(zhǎng)只由兩個(gè)方向上的空間步長(zhǎng)決定，相對(duì)于傳統(tǒng)FDTD方法，其穩(wěn)定性條件減弱，該方法非常適用于在一個(gè)方向上具有精細(xì)結(jié)構(gòu)的電磁問(wèn)題的模擬。2007年，Juan Chen等人在雙向弱條件穩(wěn)定FDTD方法的基礎(chǔ)上，又提出了一種新的單向弱條件穩(wěn)定LOD-FDTD方法[3]，該方法的時(shí)間步長(zhǎng)只由一個(gè)方向的空間步長(zhǎng)決定，相對(duì)于傳統(tǒng)FDTD方法其穩(wěn)定性條件進(jìn)一步減弱，一經(jīng)提出，就得到了廣泛應(yīng)用。

當(dāng)將LOD-FDTD方法用于解決周期結(jié)構(gòu)的電磁散射問(wèn)題時(shí)，跟傳統(tǒng)FDTD方法相同，只需在兩個(gè)方向設(shè)置吸收邊界，因此UPML要比PML具有優(yōu)勢(shì)[4]，因此，本文擬將UPML引入到LOD-FDTD方法中用于解決周期結(jié)構(gòu)的電磁散射問(wèn)題。

1 算法的差分公式

不失一般性，在電導(dǎo)率為σ，導(dǎo)磁率為σ*的介質(zhì)中，頻域內(nèi)的Maxwel方程可以寫做：

其中：

在無(wú)損條件下，假設(shè)吸收邊界設(shè)置在z軸方向，即sx=sy=1，這里定義兩個(gè)輔助變量：

將（5）和（6）式帶入（1）和（2）式中，可以得到：

假設(shè)窄邊沿著x軸方向，應(yīng)用文獻(xiàn)[7]所介紹的弱無(wú)條件穩(wěn)定時(shí)域有限差分算法，從方程（7）到（14）可以得到場(chǎng)量Ey、Hz以及輔助變量Qz的下述迭代方程，用同樣地方法也可以得到其余場(chǎng)量的相關(guān)迭代方程。

顯然，（15）式中含有未知場(chǎng)量Hz，對(duì)其進(jìn)行重新排列和替代，可以得到：

場(chǎng)量可以通過(guò)求解一特殊的非三角矩陣獲得，求得場(chǎng)量后，場(chǎng)量和輔助變量可以通過(guò)（16）式和（17）式的直接迭代求得。同理，其余各場(chǎng)量可以按照同樣地方法求得。

2 周期結(jié)構(gòu)中的應(yīng)用

考慮垂直入射波，PBC可以寫成下列形式：

其中，Nx1，Nx4代表元胞邊界上電場(chǎng)所在的節(jié)點(diǎn)位置。將（23）式帶入（18）式可以得到：

系數(shù)矩陣[M]并非三對(duì)角矩陣，因此并不能用前向消去和后向迭代法直接求解，應(yīng)用Sherman Morrison公式，按照文獻(xiàn)[6]所介紹的方法可將其化作如下的兩個(gè)線性方程進(jìn)行求解：

因此，（24）式的解可以通過(guò)下述方程求得：

顯然，通過(guò)觀察可以發(fā)現(xiàn)（25）式中的矩陣[M]與矩陣[N]相關(guān)，矩陣[N]是三對(duì)角矩陣，輔助線性方程可以利用文獻(xiàn)[8]介紹的前向消去和后向迭代有效解決。

3 數(shù)值檢驗(yàn)

為了驗(yàn)證所提算法的精度和效率，引入一調(diào)制過(guò)的高斯脈沖作為入射波，所討論局域的前端垂直于軸，入射場(chǎng)可寫作：

其中：f0=10 GHz，t0=1.125×10-10，τ=1.0×10-10，將它應(yīng)用于計(jì)算周期陣列的反射系數(shù)，周期陣列在x方向上帶有細(xì)縫，如圖1所示。參考文獻(xiàn)[9]，設(shè)置大區(qū)域和小區(qū)域，稱為ΩB和ΩS，大區(qū)域計(jì)算空間為100×20×318，小區(qū)域計(jì)算空間為100×20×78，細(xì)縫的寬度為0.1 mm，單元尺寸取為Δy=Δz=5Δx=Δ=0.5 mm，計(jì)算用16層UPML在方向截?cái)?，反射誤差表示為：

用傳統(tǒng)FDTD進(jìn)行計(jì)算時(shí)，時(shí)間步長(zhǎng)必須滿足穩(wěn)定性條件，本文中為用本文所提方法進(jìn)行計(jì)算時(shí)，時(shí)間步長(zhǎng)僅由Δy、Δz決定，為，本文中分別選用Δt=

圖1 帶有細(xì)縫的金屬薄片周期陣列Fig.1 Periodic array of metallic patches with thin slots

為了便于比較，圖2給出了用本文所提方法將UPML層設(shè)置為4層、16層以及用傳統(tǒng)的FDTD方法將UPML設(shè)置為16層時(shí)的數(shù)值反射誤差。當(dāng)用16層UPML在方向截?cái)鄷r(shí)，本文所提方法與傳統(tǒng)FDTD算法的數(shù)值反射誤差相差很小，即使用本文所提方法，設(shè)置4層UPML用于截?cái)嘤?jì)算局域時(shí)，反射誤差也在-50 dB以下。

圖3為使用不同的時(shí)間步長(zhǎng)計(jì)算時(shí)的周期結(jié)構(gòu)的反射系數(shù)隨頻率的變化關(guān)系，從圖中可以看出，既是時(shí)間步長(zhǎng)選為時(shí)，本文所提方法的計(jì)算結(jié)果也與傳統(tǒng)FDTD方法的計(jì)算結(jié)果保持高度一致。但其計(jì)算效率卻比傳統(tǒng)FDTD算法的計(jì)算效率高出好多，使用本文所提方法僅需要169.2 s，而傳統(tǒng)FDTD算法卻要耗時(shí)385.9 s。

圖2 本文所提方法的反射誤差與傳統(tǒng)FDTD方法反射誤差的比較Fig.2 Comparison of the reflection error of conventional method and proposed method

圖3 本文所提方法取不同時(shí)間步長(zhǎng)的反射系數(shù)與傳統(tǒng)TDTD（Δt=Δ/6c）反射系數(shù)的比較Fig.3 Comparison of reflection coefficient for conventional method（Δt=Δ/6c）and proposed method

4 結(jié)論

本文中，將UPML應(yīng)用于LOD-FDTD算法中用于解決周期結(jié)構(gòu)的電磁散射問(wèn)題，從數(shù)值計(jì)算結(jié)果中可以得出，用本文所提的方法，當(dāng)沿軸方向設(shè)置16層UPML層時(shí)，可將反射誤差降到-100 dB以下，即使設(shè)置4層UPML，反射誤差也在-50 dB以下。因此，本文所提方法在解決周期結(jié)構(gòu)電磁散射具有一定的優(yōu)勢(shì)。

[1] T.Namiki.A new FDTD algorithm based on alternatingdirection implicit method[J].IEEE Trans.Microwave Theory Tech.，1999，47（10）：2003-2007.

[2] Binke Huang，Gang Wang，Yansheng Jiang.A Hybrid implicit-explicit FDTD scheme with weakly conditional stability[J].Microwave and Optical Tech.Lett.，2003，39（2）:97-101.

[3] Chen J，Wang J.3-D FDTD method with weakly conditional stability[J].Electronics Letters，2007，43（1）：2-3.

[4] Berenger J P.A perfectly matched layer for the absorption of electromagnetic waves[J].J.Comput.Phys.，1994（114）：185-200.

[5]Taflove A，Hagness S C.Computational Electrodynamics:The Finite-Difference Time-Domain Method[M].second edition，Boston，MA:Artech House，2000.

[6] Thomas J W.Numerical Partial Differential Equations:Finite Difference Methods[M].Berlin，Germany:Springer Verlag，1995.

[7] Huang B K，Wang G，Jiang Y S，et al.A hybrid implicit explicit FDTD scheme with weakly conditional stability[J].Microw.Opt.Tech.Lett.，2003（39）：97-101.

[8] Zhao A P.Two special notes on the implementation of the unconditionally stable ADI FDTD method[J].Microw.Opt.Technol Lett，2002，33（4）：273-277.

[9] Thomas G M，Jeffrey G B，Taflove A，et al.Theory and application of radiation boundary operators[J].IEEE Trans.Antennas Propagat.，1988，36（12）:1797-1812.