三維周期結構弱無條件穩定FDTD算法研究

李洪宇，王發年，劉宗信，2，李林，李軍

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

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

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

1 算法的差分公式

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

其中：

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

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

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

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

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

2 周期結構中的應用

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

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

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

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

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

3 數值檢驗

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

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

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

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

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

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

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

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

4 結論

本文中，將UPML應用于LOD-FDTD算法中用于解決周期結構的電磁散射問題，從數值計算結果中可以得出，用本文所提的方法，當沿軸方向設置16層UPML層時，可將反射誤差降到-100 dB以下，即使設置4層UPML，反射誤差也在-50 dB以下。因此，本文所提方法在解決周期結構電磁散射具有一定的優勢。

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