新型分裂步長(zhǎng)時(shí)域有限差分法

林智參+班濤

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

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

中圖分類號(hào)： TN802?34； O441 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 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 引 言

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

分步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方案，對(duì)矩陣方程式（3）～式（6）的右端進(jìn)行近似，進(jìn)一步化簡(jiǎ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]的單位矩陣。四個(gè)分步內(nèi)需要求解的方程式以第（1）步為例化簡(jiǎn)如下：

[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]的三對(duì)角矩陣方程式，其具體形式如下：

[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）中，此時(shí)式（13）為關(guān)于[Hyn+14i+12， j]的一個(gè)顯式方程，可直接進(jìn)行求解。因此，在分步1內(nèi)只需要求解一個(gè)隱式方程和一個(gè)顯式方程。分步2，3和4內(nèi)也采用與分步1內(nèi)相同的方式進(jìn)行處理，而分步3，4中得到的不是三對(duì)角矩陣方程式，是顯式方程，求解變得更簡(jiǎn)單。

2 數(shù)值色散分析

利用Fourier方法[10]在第[n]個(gè)時(shí)間步內(nèi)，場(chǎng)分量在空間區(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）中，得到一個(gè)完整時(shí)間步長(zhǎng)內(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)生的電磁場(chǎng)滿足：

[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，λ代表波長(zhǎng)，][N]表示單位波長(zhǎng)元胞數(shù)，傳輸角度為θ。圖1為[S=0.5，][N=8]時(shí)歸一化數(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)，電場(chǎng)分量[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)，電場(chǎng)分量[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]場(chǎng)分量上加正弦波激勵(lì)源[sin（2πft），]其中，[f=1.5] GHz。網(wǎng)格尺寸為[Δx=Δy=1 cm，]為激勵(lì)源最高頻率對(duì)應(yīng)波長(zhǎ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]場(chǎng)分量的空間分布圖，圖4為傳統(tǒng)FDTD算法和NSS?FDTD算法比較圖。

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

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

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

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

5 結(jié) 論

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

參考文獻(xiàn)

[1] 葛德彪，閆玉波.電磁波時(shí)域有限差分方法[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.