三維周期結(jié)構(gòu)弱無條件穩(wěn)定FDTD算法研究

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

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

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

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

1 算法的差分公式

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

其中：

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

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

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

顯然，（15）式中含有未知場量Hz，對其進行重新排列和替代，可以得到：

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

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

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

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

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

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

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

3 數(shù)值檢驗

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

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

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

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

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

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

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

圖3 本文所提方法取不同時間步長的反射系數(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)的電磁散射問題，從數(shù)值計算結(jié)果中可以得出，用本文所提的方法，當沿軸方向設(shè)置16層UPML層時，可將反射誤差降到-100 dB以下，即使設(shè)置4層UPML，反射誤差也在-50 dB以下。因此，本文所提方法在解決周期結(jié)構(gòu)電磁散射具有一定的優(yōu)勢。

[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.