新型分裂步長時域有限差分法

林智參+班濤

摘 要： 提出一種新型的分裂步長時域有限差分（NSS?FDTD）法，并對其數(shù)值色散進(jìn)行分析。該方法基于Split?Step方案和Crank?Nicolson方案，采用新的矩陣分解形式，與傳統(tǒng)的FDTD算法、SS?FDTD算法相比，減少了計(jì)算復(fù)雜度。新型算法的推導(dǎo)程序簡單，且具有良好的數(shù)值色散特性，還加入了一階Mur吸收邊界條件，給出一階Mur吸收邊界差分方程。將數(shù)值實(shí)驗(yàn)的結(jié)果和傳統(tǒng)FDTD方法及理論值進(jìn)行比較，數(shù)值結(jié)果一致性較好。

關(guān)鍵詞： 時域有限差分法； 分裂步長； Split?Step方案； 數(shù)值色散

中圖分類號： TN802?34； O441 文獻(xiàn)標(biāo)識碼： A 文章編號： 1004?373X（2015）15?0117?03

New method of finite difference time domain for split step

LIN Zhican， BAN Tao

（South China Normal University， Guangzhou 510006， China）

Abstract： A new split step finite difference time domain （NSS?FDTD） algorithm is presented， and its numerical dispersion is analyzed. The method is based on the schemes of Split?Step and Crank?Nicolson， adopted new matrix decomposition form. Compared with traditional algorithms of FDTD and SS?FDTD， the proposed algorithm can reduce computational complexity， and has simple deduction procedure and better numerical dispersion characteristic. The first?order Mur absorbing boundary condition is added in this paper， and its difference equation is presented. The numerical experiment results were compared with traditional FDTD method and theoretical values. The consistence of numerical results is better.

Keywords： FDTD method； split step； Split?Step scheme； numerical dispersion

0 引 言

時域有限差分法（Finite Difference Time Domain，F(xiàn)DTD）是一種簡單直觀的全波分析時域算法[1?3]，該方法以Yee氏立體網(wǎng)格作為電磁場離散單元，將麥克斯韋方程轉(zhuǎn)化為差分方程，能夠方便有效地結(jié)合計(jì)算機(jī)技術(shù)處理復(fù)雜的電磁場問題，目前已經(jīng)在電磁學(xué)的各個領(lǐng)域中得到了廣泛應(yīng)用。然而，應(yīng)用傳統(tǒng)的FDTD方法也明顯體現(xiàn)出其不足之處，為減小差分近似帶來的數(shù)值色散[4]，空間網(wǎng)格尺寸必須遠(yuǎn)小于波長， 這樣反而增加運(yùn)算負(fù)擔(dān)，因而，出現(xiàn)了多鐘FDTD的改進(jìn)方法[5?7]。本文以TM波為例，提出一種基于Split?Step方案[8]和Crank?Nicolson方案[9]的時域有限差分法，所提出的算法采用新的矩陣分解方式，簡化計(jì)算復(fù)雜度、減少差分近似所帶來的數(shù)值色散。

1 NSS?FDTD算法理論推導(dǎo)

考慮空間一無源區(qū)域，在均勻無耗、各向同性介質(zhì)中，介電常數(shù)為[ε，]磁導(dǎo)率為[μ，]二維TM波Maxwell微分方程組的矩陣形式如下：

[?u?t=Mu] （1）

式中：[u=[Ez，Hx，Hy]T， M=0-1ε??y1ε??x-1μ??y001μ??x00。]

將矩陣[M]分解成四個子矩陣，分別記為[A2，][B2，][C2，][D2，]矩陣形式分別如下：

[A=001ε??x0001μ??x00 B=0-1ε??y0-1μ??y00000]

[C=0-1ε??y00001μ??x00 D=001ε??x-1μ??y00000]

公式（1）可以寫為：

[?u?t=A2u+B2u+C2u+D2u] （2）

利用Split?Step方案，將式（2）分解成四個子方程式來求解。同時，從[n～][n+1]時間步，將一個時間步長等間隔分成四個子時間分步，即[n→n+14，][n+14→][n+24，][n+24→n+34]和[n+34→n+1，]得到以下四個子矩陣方程式：

分步1：

[?u?t=4?A2u n→n+14] （3）

分步2：

[?u?t=4?B2u n+14→n+24] （4）

分步3：

[?u?t=4?C2u n+24→n+34] （5）

分步4：

[?u?t=4?D2u n+34→n+1] （6）

利用Crank?Nicolson方案，對矩陣方程式（3）～式（6）的右端進(jìn)行近似，進(jìn)一步化簡得到以下形式：

[I-Δt4Aun+14=I+Δt4Aun] （7）

[I-Δt4Bun+24=I+Δt4Bun+14] （8）

[I-Δt4Cun+34=I+Δt4Cun+24] （9）

[I-Δt4Dun+1=I+Δt4Dun+34] （10）

式中：[I]為[3×3]的單位矩陣。四個分步內(nèi)需要求解的方程式以第（1）步為例化簡如下：

[Ezn+14i， j-Δt4εΔxHyn+14i+12， j-Hyn+14i-12， j= Ezni， j+Δt4εΔxHyni+12， j-Hyni-12， j] （11）

[Hxn+14i， j+12=Hxni， j+12] （12）

[Hyn+14i+12， j-Δt4μΔxEzn+14i+1， j-Ezn+14i， j= Hyni+12， j+Δt4μΔxEzni+1， j-Ezni， j] （13）

在分步1內(nèi)，只需要求解方程式（11）和（13），并且方程式（11）和（13）是互為耦合的方程式。將式（13）代入到式（11）中，消去[Hyn+14i+12， j，]得到關(guān)于[Ezn+14i， j]的三對角矩陣方程式，其具體形式如下：

[1+Δt28μεΔx2Ezn+14i， j-Δt216μεΔx2Ezn+14i+1， j+Ezn+14i-1， j= 1-Δt28μεΔx2Ezni， j+Δt216μεΔx2Ezni+1， j+Ezni-1， j+ Δt2εΔxHyni+12， j-Hyni-12， j]

由式（14）求出[Ezn+14i， j，]然后代入到式（13）中，此時式（13）為關(guān)于[Hyn+14i+12， j]的一個顯式方程，可直接進(jìn)行求解。因此，在分步1內(nèi)只需要求解一個隱式方程和一個顯式方程。分步2，3和4內(nèi)也采用與分步1內(nèi)相同的方式進(jìn)行處理，而分步3，4中得到的不是三對角矩陣方程式，是顯式方程，求解變得更簡單。

2 數(shù)值色散分析

利用Fourier方法[10]在第[n]個時間步內(nèi)，場分量在空間區(qū)域內(nèi)的表達(dá)形式如下：

[UnI， J=Une-j（kxIΔx+kyJΔy）] （15）

將式（15）代入到式（7）～（10）中，并將式（7）～（10）進(jìn)一步整理為如下矩陣形式：

[Un+14=A2dUn] （16）

[Un+24=B2dUn+14] （17）

[Un+34=C2dUn+24] （18）

[Un+1=D2dUn+34] （19）

將式（16）～（18）代入到式（19）中，得到一個完整時間步長內(nèi)的矩陣方程式：

[Un+1=D2dC2dB2dA2dUn=ΩUn] （20）

可以求得[Ω]的特征值，其結(jié)果如下：

[λ1=1，λ2=λ3*=P+jQ2-P2Q] （21）

其中：

[P=-P4xP4yb4d4+16（P4x+P4y）b2d2+104P2xP2yb2d2-224bd（P2x+P2y）+128]

[Q=8P2xP2yb2d2+32bd（P2x+P2y）+128Pα=-2ΔαsinkαΔα2，α=x，y；b=Δt2ε；d=Δt2μ ]

利用von Neumann方法[11]，假定一角頻率為[ω]的電磁波產(chǎn)生的電磁場滿足：

[Enz=EzejωΔtn，Hnα=HαejωΔtn，α=x，y] （22）

將式（22）代入到式（20）中，得：

[ejωΔtI-ΩUn=0] （23）

其中，[Un]與初始值[U0]相關(guān)，其具體關(guān)系式如下：

[Un=U0ejωΔtn] （24）

為了使式（23）中 [Un]有非零解，[Un]的系數(shù)行列式的值應(yīng)為零，即：

[det（ejωΔtI-Ω）=0] （25）

由式（21）得[Ω]的特征值，可以得到 NSS?FDTD 算法的數(shù)值色散表達(dá)式，其形式如下：

[tan2ωΔt2=P-QP+Q=bd（P4xP4yb3d3-16P4xbd-96P2xP2ybd-16P4ybd+256P2x+256P2y）P4xP4yb4d4-16（P4x+P4y）b2d2-112P2xP2yb2d2+192（P2x+P2y）bd-256]

為了方便分析NSS?FDTD的數(shù)值色散特性，作如下定義：[S=cΔtΔx，Δx=Δy，N=λΔx，λ代表波長，][N]表示單位波長元胞數(shù)，傳輸角度為θ。圖1為[S=0.5，][N=8]時歸一化數(shù)值相位速度隨傳輸角度θ的變化曲線。

圖1 歸一化數(shù)值相位速度隨傳輸角度[θ]的變化曲線

從圖1可以看出，NSS?FDTD算法的歸一化相位速度大于傳統(tǒng)的FDTD算法和傳統(tǒng)SS?FDTD算法的歸一化數(shù)值相位速度，并接近1，且曲線變化值在0.98～0.995范圍內(nèi)，比較平緩，歸一化數(shù)值相位速度各向異性誤差也比較小。

3 一階Mur吸收邊界條件

在分步1內(nèi)，電場分量[Ez]在[i=1，][i=Imax，][j=1]和[j=Jmax]上的一階Mur吸收邊界差分方程式為：

[Ezn+141， j=Ezn2， j+vΔt-4ΔxvΔt+4ΔxEzn+142， j-Ezn1， j] （26）

[Ezn+14Imax， j=EznImax-1， j+vΔt-4ΔxvΔt+4ΔxEzn+14Imax-1， j-EznImax， j]（27）[Ezn+14i，1=Ezni，2+vΔt-4ΔxvΔt+4ΔxEzn+14i， 2-Ezni，1] （28）

[Ezn+14i， Jmax=Ezni， Jmax-1+vΔt-4ΔxvΔt+4ΔxEzn+14i， Jmax-1-Ezni， Jmax] （29）

在分步2，3和4內(nèi)，電場分量[Ez]在[i=1，][i=Imax，][j=1]和[j=Jmax]上的一階Mur吸收邊界差分方程式與分步1內(nèi)的差分方程式類似，此處不再一一進(jìn)行展開說明。

4 計(jì)算結(jié)果

將NSS?FDTD算法運(yùn)算于尺寸101 cm×101 cm的自由空間，以一階Mur為邊界條件，以二維TM波為例，在中心區(qū)域[Ez]場分量上加正弦波激勵源[sin（2πft），]其中，[f=1.5] GHz。網(wǎng)格尺寸為[Δx=Δy=1 cm，]為激勵源最高頻率對應(yīng)波長的[120，]計(jì)算網(wǎng)格數(shù)為101×101。觀察點(diǎn)在中心區(qū)域和吸收邊界之間的中心位置。NSS?FDTD數(shù)值計(jì)算結(jié)果仿真如圖2～圖4所示。其中圖2，圖3兩圖為NSS?FDTD數(shù)值計(jì)算過程中[Ez]場分量的空間分布圖，圖4為傳統(tǒng)FDTD算法和NSS?FDTD算法比較圖。

圖2 運(yùn)行50步時[Ez]的3D圖

圖3 運(yùn)行200步時[Ez]的3D圖

圖4 兩種FDTD算法在觀察點(diǎn)處的電場值[Ez]

由圖可見，NSS?FDTD算法的計(jì)算結(jié)果與傳統(tǒng)的FDTD算法的計(jì)算結(jié)果吻合的很好，且符合電磁場理論，從而證實(shí)新型分裂步長時域有限差分法的可行性。

5 結(jié) 論

本文基于Split?Step和Crank?Nicolson方案提出了一種新型的二維FDTD算法，改進(jìn)算法采用新的分解形式，與傳統(tǒng)的SS?FDTD算法相比較，減少了計(jì)算復(fù)雜度，優(yōu)化了計(jì)算公式，使推導(dǎo)過程更簡單。結(jié)合算例使用Matlab對NSS?FDTD算法進(jìn)行編程分析，結(jié)果表明，該算法具有良好的預(yù)期效果。

參考文獻(xiàn)

[1] 葛德彪，閆玉波.電磁波時域有限差分方法[M].西安：西安電子科技大學(xué)出版社，2002.

[2] 王秉中.計(jì)算電磁學(xué)[M].北京：科學(xué)出版社，2002.

[3] YEE K S. Numerical solution of initial boundary value problems involving Maxwell′s equations in isotropic media [J]. IEEE Transactions on Antennas and Propagation， 1966， 14 （3）： 302?307.

[4] TAFLOVE A， HAGNESS S C. Computational electrodynamics： the finite?difference time?domain method [M]. 2nd ed. Boston： Artech House， 2000.

[5] 徐利軍，衰乃昌.高階ADI?FDTD算法的數(shù)值色散分析[J].電子與信息學(xué)報(bào)，2005，27（10）：1662?1665.

[6] 黨濤，鄭宏興.關(guān)于二維ADI?FDTD方法的數(shù)值色散分析[J].中國民航學(xué)院學(xué)報(bào)，2004，22（2）：42?46.

[7] 夏冬，黨濤，鄭宏興.一維ADI?FDTD方法的數(shù)值色散分析[J].中國民航學(xué)院學(xué)報(bào)，2005，23（2）：38?41.

[8] LEE J， FORNBERG B. A split step approach for the 3?D Maxwell′s equations [J]. Journal of Computational and Applied Mathematics， 2003， 158（2）： 485?505.

[9] SMITH G D. Numerical solution of partial differential equations： finite difference methods [M]. 3rd ed. Oxford： Oxford University Press， 1986.

[10] SMITH G D. Numerical solution of partial differential equations [M]. Oxford： Oxford University Press， 1978.

[11] PEREDA J A， VIELVA L A， VEGAS A， et al. Analyzing the stability of the FDTD technique by combining the Von Neumann method with the Routh?Hurwitz criterion [J]. IEEE Transactions on Microwave Theory and Techniques， 2001， 49（2）： 377?381.