新型分裂步长时域有限差分法

林智参+班涛

摘 要： 提出一种新型的分裂步长时域有限差分（NSS?FDTD）法，并对其数值色散进行分析。该方法基于Split?Step方案和Crank?Nicolson方案，采用新的矩阵分解形式，与传统的FDTD算法、SS?FDTD算法相比，减少了计算复杂度。新型算法的推导程序简单，且具有良好的数值色散特性，还加入了一阶Mur吸收边界条件，给出一阶Mur吸收边界差分方程。将数值实验的结果和传统FDTD方法及理论值进行比较，数值结果一致性较好。

关键词： 时域有限差分法； 分裂步长； Split?Step方案； 数值色散

中图分类号： TN802?34； O441 文献标识码： 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，FDTD）是一种简单直观的全波分析时域算法[1?3]，该方法以Yee氏立体网格作为电磁场离散单元，将麦克斯韦方程转化为差分方程，能够方便有效地结合计算机技术处理复杂的电磁场问题，目前已经在电磁学的各个领域中得到了广泛应用。然而，应用传统的FDTD方法也明显体现出其不足之处，为减小差分近似带来的数值色散[4]，空间网格尺寸必须远小于波长， 这样反而增加运算负担，因而，出现了多钟FDTD的改进方法[5?7]。本文以TM波为例，提出一种基于Split?Step方案[8]和Crank?Nicolson方案[9]的时域有限差分法，所提出的算法采用新的矩阵分解方式，简化计算复杂度、减少差分近似所带来的数值色散。

1 NSS?FDTD算法理论推导

考虑空间一无源区域，在均匀无耗、各向同性介质中，介电常数为[ε，]磁导率为[μ，]二维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）的右端进行近似，进一步化简得到以下形式：

[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]的单位矩阵。四个分步内需要求解的方程式以第（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内，只需要求解方程式（11）和（13），并且方程式（11）和（13）是互为耦合的方程式。将式（13）代入到式（11）中，消去[Hyn+14i+12， j，]得到关于[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）为关于[Hyn+14i+12， j]的一个显式方程，可直接进行求解。因此，在分步1内只需要求解一个隐式方程和一个显式方程。分步2，3和4内也采用与分步1内相同的方式进行处理，而分步3，4中得到的不是三对角矩阵方程式，是显式方程，求解变得更简单。

2 数值色散分析

利用Fourier方法[10]在第[n]个时间步内，场分量在空间区域内的表达形式如下：

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

将式（15）代入到式（7）～（10）中，并将式（7）～（10）进一步整理为如下矩阵形式：

[Un+14=A2dUn] （16）

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

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

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

将式（16）～（18）代入到式（19）中，得到一个完整时间步长内的矩阵方程式：

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

可以求得[Ω]的特征值，其结果如下：

[λ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]，假定一角频率为[ω]的电磁波产生的电磁场满足：

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

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

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

其中，[Un]与初始值[U0]相关，其具体关系式如下：

[Un=U0ejωΔtn] （24）

为了使式（23）中 [Un]有非零解，[Un]的系数行列式的值应为零，即：

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

由式（21）得[Ω]的特征值，可以得到 NSS?FDTD 算法的数值色散表达式，其形式如下：

[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的数值色散特性，作如下定义：[S=cΔtΔx，Δx=Δy，N=λΔx，λ代表波长，][N]表示单位波长元胞数，传输角度为θ。图1为[S=0.5，][N=8]时归一化数值相位速度随传输角度θ的变化曲线。

图1 归一化数值相位速度随传输角度[θ]的变化曲线

从图1可以看出，NSS?FDTD算法的归一化相位速度大于传统的FDTD算法和传统SS?FDTD算法的归一化数值相位速度，并接近1，且曲线变化值在0.98～0.995范围内，比较平缓，归一化数值相位速度各向异性误差也比较小。

3 一阶Mur吸收边界条件

在分步1内，电场分量[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内，电场分量[Ez]在[i=1，][i=Imax，][j=1]和[j=Jmax]上的一阶Mur吸收边界差分方程式与分步1内的差分方程式类似，此处不再一一进行展开说明。

4 计算结果

将NSS?FDTD算法运算于尺寸101 cm×101 cm的自由空间，以一阶Mur为边界条件，以二维TM波为例，在中心区域[Ez]场分量上加正弦波激励源[sin（2πft），]其中，[f=1.5] GHz。网格尺寸为[Δx=Δy=1 cm，]为激励源最高频率对应波长的[120，]计算网格数为101×101。观察点在中心区域和吸收边界之间的中心位置。NSS?FDTD数值计算结果仿真如图2～图4所示。其中图2，图3两图为NSS?FDTD数值计算过程中[Ez]场分量的空间分布图，图4为传统FDTD算法和NSS?FDTD算法比较图。

图2 运行50步时[Ez]的3D图

图3 运行200步时[Ez]的3D图

图4 两种FDTD算法在观察点处的电场值[Ez]

由图可见，NSS?FDTD算法的计算结果与传统的FDTD算法的计算结果吻合的很好，且符合电磁场理论，从而证实新型分裂步长时域有限差分法的可行性。

5 结 论

本文基于Split?Step和Crank?Nicolson方案提出了一种新型的二维FDTD算法，改进算法采用新的分解形式，与传统的SS?FDTD算法相比较，减少了计算复杂度，优化了计算公式，使推导过程更简单。结合算例使用Matlab对NSS?FDTD算法进行编程分析，结果表明，该算法具有良好的预期效果。

参考文献

[1] 葛德彪，闫玉波.电磁波时域有限差分方法[M].西安：西安电子科技大学出版社，2002.

[2] 王秉中.计算电磁学[M].北京：科学出版社，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算法的数值色散分析[J].电子与信息学报，2005，27（10）：1662?1665.

[6] 党涛，郑宏兴.关于二维ADI?FDTD方法的数值色散分析[J].中国民航学院学报，2004，22（2）：42?46.

[7] 夏冬，党涛，郑宏兴.一维ADI?FDTD方法的数值色散分析[J].中国民航学院学报，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.