基于二维TE波常用时域有限差分算法的分析

何光峰　迟洁茹

摘要： 针对传统的时域有限差分法受Courant稳定条件的限制，且存在交替方向隐式时域有限差分法数值色散较大的问题，本文以TE波为例，研究了CrankNicoloson差分方式的近似去耦时域有限差分法基本原理，并对其稳定性进行分析，证明该方法是无条件稳定。通过数值仿真，从运行时间和吸收效果方面与传统的时域有限差分法和交替方向隐式时域有限差分法进行对比。仿真结果表明，近似去耦时域有限差分法比交替方向时域有限差分法的吸收效果好，但比传统的时域有限差分法吸收效果差；近似去耦时域有限差分法比交替方向时域有限差分法运行时间长，但比传统时域有限差分法运行时间短，说明近似去耦时域有限差分法突破了Courant稳定条件的限制，且在吸收效果方面比交替方向时域有限差分法好。该研究具有广阔的应用前景。

关键词： 时域有限差分法； CrankNicoloson差分格式方案； 无条件稳定

中图分类号： O241.8；TB5文献标识码： A

1966年，Yee提出传统时域有限差分（finitedifference timedomain，FDTD）方法，该方法得到广泛应用与发展[13]。时域有限差分法是用古典显式差分的方法对麦克斯韦方程进行差分，进而对麦克斯韦方程进行微分求解。运用显式差分的FDTD方法必须满足Courant稳定条件，即时间步数的选取受到空间步长（离散网格大小）的限制。为降低数值色散误差，在空间步数选取时必须满足远远小于其采样频率的波长的条件，这就使其在解决实际问题的过程中，为得到更加真实的电磁波特性，不得不采用非常小的空间步数，进而导致计算量过大而无法实现。为克服传统FDTD方法的缺陷，T.Namiki等人[47]提出了交替方向隐式（alternative direction implicit，ADI）FDTD方法。由于ADIFDTD方法对麦克斯韦方程组的所有场分量进行交替差分，从而实现了无条件稳定，克服了传统FDTD方法必须满足Courant稳定条件的限制，但ADIFDTD方法在一个时间步数上需要进行两次迭代，增加了每步计算的时间量和存储空间，且由于每步差分方程时间不同步，对计算誤差产生了非常大的影响。2003年，Sun G L等人[811]利用CrankNicoloson半隐式差分方式，对麦克斯韦方程组进行差分，使用近似去耦（approximately decoupling，AD）法而提出的CrankNicoloson差分方式的近似去耦（cranknicoloson approximatelydecoupling，CNAD）FDTD方法。CNADFDTD方法在使用CrankNicoloson半隐式差分方式对二维麦克斯韦方程求解过程中，需要求解一个大型稀疏矩阵方程组，从而占用计算机大量内存和运行时间，因此采用近似方法求解过稀疏矩阵，将E在n时刻的值代替n+1时刻的值，从而把对稀疏矩阵求解变成对三对角矩阵求解，使计算过程简单化，并减小所需的计算机内存和运行时间。基于此，本文以横电（transverseelectric，TE）波为例，介绍了CNADFDTD方法，并对FDTD、ADIFDTD和CNADFDTD三种方法分别进行仿真实验。实验结果表明，在相同条件下，CNADFDTD比ADIFDTD的吸收效果要好，但与FDTD有一定的差距，CNADFDTD的运行时间比ADIFDTD的运行时间长。

1CNADFDTD方法

1.1CNADFDTD理论

二维TE波在x，y和z三个方向上存在Ex，Ey和Hz分量，写成CN差分方程[911]分别为

Ex|n+1i+1/2，j=Ex|ni+1/2，j+a1ΓyHz|i+1/2，j+1/2/Δy （1a）

Ey|n+1i，j+1/2=Ey|ni，j+1/2-a1ΓxHz|i+1/2，j+1/2/Δx（1b）

Hz|n+1i+1/2，j+1/2=Hz|n+1i+1/2，j+1/2+a2ΓyEx|i+1/2，j+1/Δy-a2ΓxEy|i+1，j+1/2/Δx（1c）

其中

ΓyHz|i+1/2，j+1/2=Hz|n+1i+1/2，j+1/2-Hz|n+1i+1/2，j-1/2+Hz|ni+1/2，j+1/2-Hz|ni+1/2，j-1/2（2）

式中，a1=Δt/2ε，a2=Δt/2μ；Δt为时间步长，Δx，Δy分别为空间步长；ε为介质介电常数，μ为磁导系数。其迭代过程如下：

1）将式（1c）代入式（1a）和式（1b），消去n+1时刻磁场值Hn+1，可得n+1时刻电场值En+1，实现对电场值迭代。

2）将求得的电场值En+1代入式（1c），更新磁场值Hn+1。其中，在对电场值进行迭代中，将Eny代替En+1y，可将对大型非有限带宽稀疏矩阵的求解简化为对三对角矩阵的求解，从而将式（1a）化简为

1+2b2xEx|n+1i+1/2，j-b2xEx|n+1i+1/2，j-1+Ex|n+1i+1/2，j+1=1-2b2xEx|ni+1/2，j+b2x（Ex|ni+1/2，j-1+Ex|ni+1/2，j+1）+2a1Hz|ni+1，j+1/2-Hz|ni+1，j-1/2-2bxbyEy|ni+1，j+1/2-Ey|ni，j+1/2-Ey|ni+1，j-1/2+Ey|ni，j-1/2（3）

其中，by=cΔt/2Δx，bx=cΔt/2Δy。

最终迭代过程为：对式（3）进行迭代；对式（1b）迭代求解En+1y；对式（3）和式（1b）更新；对式（1c）的磁场值迭代更新。

1.2CNADFDTD的稳定性分析

根据V.Neumann稳定性分析，二维情况下的电磁场分量Ex，Ey和Hz[12]分别为

Ex=ψAεnexp[j（kxΔx+kyΔy）]， Ey=ψBεnexp[j（kxΔx+kyΔy）]， Hz=ψCεnexp[j（kxΔx+kyΔy）]（4）

式中，ψA，ψB，ψC分別为各电磁分量初始系数；ε是增长因子；kx和ky是波常数；Δx，Δy分别代表沿x方向和沿y方向网格的大小。

将式（4）代入式（1a）～（1c）中，消去系数得关于增长因子的方程为

ε-12=-a（ε+1）2（5）

其中，a=Δt2εμΔx2sin2kxΔx/2+Δt2εμΔy2sin2（kyΔy/2）。将a（a>0）代入式（5），可求得ε为

ε=1-a±2aj1+a（6）

由于式（6）中的增长因子满足ε=1，因此时间步数无论取何值，CNADFDTD都是稳定的。

1.3CNADFDTD在完全匹配层的公式

以有限空间模拟无限空间时，以完全匹配层[1315]（perfectly matched layer，PML）作为吸收边界，则二维TE波在PML吸收介质中，按照CNADFDTD的方法进行差分[1617]，得

Qex（m）Ex|n+1i+1/2，j=Dex（m）Ex|ni+1/2，j+ΓyHz|i+1/2，j+1/2Qey（m）Ey|n+1i，j+1/2=Dey（m）Ey|ni，j+1/2+ΓxHz|i+1/2，j+1/2Hzx|n+1i+1/2，j+1/2=Qhx（m）Hzx|n+1i+1/2，j+1/2-Dhx（m）ΓxEy|i+1，j+1/2Hzy|n+1i+1/2，j+1/2=Qhy（m）Hzy|n+1i+1/2，j+1/2+Dhy（m）ΓyEx|i+1/2，j+1（7）

其中，m取值与左端场分量节点的空间位置相同；且。

Qex（m）=（2ε0+σy（m）Δt）， Dhy（m）=12μ0+σmy（m）ΔtΔtΔy（8）

将磁场分量代入迭代的点场分量求值中，利用CNADFDTD，先求出En+1x和En+1y，再对磁场值进行迭代。

2数值分析

为验证CNADFDTD方法的正确性，以二维平面上的一个点源的辐射场为例，分别对传统FDTD、ADIFDTD及CNADFDTD方法使用C语言编写程序，在Visual Studio 2015中进行仿真，对得到的数据利用Tecplot 360软件进行成像处理。

仿真条件：仿真空间为20 mm×20 mm均匀网格，磁场激励源为正弦波，位于仿真空间中心位置，频率f=1×1012 Hz，c=30×108 m/s，波长λ=3 mm，空间分辨率Δx=Δy=01 mm，时间步数Δt=0155 ps，吸收边界厚度为1 mm。PML中电导率的分布为

σz（z）=σmaxz-z0mdm（9）

式中，d是PML的厚度；z0为PML层靠近FDTD区的界面位置；m为整数。为了得到更好的吸收效果，σmax的最佳取值[1820]为

σmax=-（m+1）ln[R（θ）]2ηd（10）

其中，R（θ）是理论反射系数值，R（θ）=10-6。

通过仿真实验，得到距离中心点源为5个网格处的磁场分量，3种方法的数值比较结果如图1所示。由图1可以看出，3种方法的结果基本相同，验证了CNADFDTD方法的正确性；通过仿真计算，得到坐标为（189，189）处的反射误差如图2所示。由图2可以看出，在相同条件下，传统FDTD的吸收效果最好，CNADFDTD的吸收效果比ADIFDTD的吸收效果要好，主要由于ADIFDTD差分等式两边的时间不对称所导致，而CNADFDTD由于只需要一步，且时间上对称，故其效果比ADIFDTD要好。

3结束语

本文主要讨论无条件稳定的CNADFDTD，由于其不受稳定条件的限制，所以可通过改变时间步的大小，缩短CPU的运行时间。通过与传统的FDTD和ADIFDTD两种方法的对比可以得出，CNADFDTD的吸收效果比ADIFDTD的吸收要好，但运行时间比ADIFDTD的多，且吸收效果相比于FDTD还有一段差距。因此，CNADFDTD是一种值得深入研究和推广应用的时域算法。

参考文献：

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

[2]Yee K S. Numerical Solution of Initial Boundary Value Problems Involving MaxwellS Equation in Isotropic Media[J]. IEEE Transactions on Antennas and Propagation， 1966， 14（3）： 302307.

[3]Schneider J B. A Selective Survey of the FiniteDifference TimeDomain Literature[J]. IEEE Antennas and propagation， 1995， 37（4）： 3957.

[4]Namiki T. A New FDTD Algorithm Based on AlternatingDirection Implicit Method [J]. IEEE Transaction on Microwave Theory and Techniques， 1999， 47（10）： 20032007.

[5]Namiki T. 3D ADIFDTD MethodUnconditionally Stable TimeDomain Algorithm for Solving Full Vector Maxwells Equations [J]IEEE Transaction on Microwave Theory and Techniques， 2000， 48（10）： 17431748.

[6]Zheng F， Chen Z， Zhang J. Toward the Development of a ThreeDimensional Unconditionally Stable FiniteDifference TimeDomain Method[J]. IEEE Transaction on Microwave Theory and Techniques， 2000， 48（9）： 15501558.

[7]Zheng Hongxing， Leung K W. A Nonorthogonal ADIFDTD Algorithm for Solving Two Dimensional Scattering Problems[J]. IEEE Transaction on Antennas and Propagation， 2009， 57（12）： 38913902.

[8]Sun G L， Trueman C W. Approximate CrankNicolson Schemes for the FiniteDifference TimeDomain Method for TEz Waves[J]. IEEE Transaction on Antennas and Propagation， 2004， 52（11）： 29632972.

[9]Sun G L， Trueman C W. Unconditionally Stable CrankNicolson Scheme for Solving the TwoDimensional Maxwells Equations[J]. Electronics Letters， 2003， 39（7）： 595597.

[10]Xie X， Pan G， Hall S. A CrankNicholsonBase Unconditionally Stable TimeDomain Algorithm for 2D and 3D problem[J]. Microwave and Optical Technology Letters， 2007， 49（2）： 261265.

[11]Li J X， Jiang H L， Zhao X M， et al. Effective CNADand ADEBased CFSPML Formulations for Truncating the Dispersive FDTD Domains[J]. IEEE Antennas and Wireless Propagation Letters， 2015， 14： 12671270.

[12]王秉中. 計算电磁学[M]. 北京： 科学出版社， 2002.

[13]Berenger J P. A Perfectly Matched Layer for the Absorption of Electromagnetic Waves[J]. Journal of Computational Physics， 1994， 114（9）： 185200.

[14]Berenger J P. Perfectly Matched Layer for the FDTD Solution of WaveStructure Interaction Problem[J]. IEEE Transaction on Antennas and Propagation， 1996， 44（1）： 110117.

[15]Berenger J P. ThreeDimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves[J]. Journal of Computational Physics， 1996， 127（2）： 363379.

[16]Stanislav O， Pan G W. An Updated Review of General Dispersion Relation for Conditionally and Unconditionally Stable FDTD Algorithms[J]. IEEE Transaction on Antennas and Propagation， 2008， 56（8）： 25722583.

[17]Kusaf M， Oztoprak A Y. An Unco Nditionally Stable SplitStep FDTD Method for Low Anisotropy[J]. IEEE Microwave and Wireless Components Letters， 2008， 18（4）： 224226.

[18]Lu Yijun， Shen C Y. A Domain Decomposition FiniteDifference Method for Parallel Numerical Implementation of TimeDependent Maxwells Equations[J]. IEEE Transactions on Antennas and Propagation， 1997， 45（3）： 556562.

[19]Park J H， Strikwerda J C. The Domain Decomposition Method for Maxwells Equations in Time Domain Simulations with Dispersive Metallic Media[J]. Society for Industrial and Applied Mathematics， 2014， 32（2）： 684702.

[20]Taflove A， ME Brodwin. Computation of the Electromagnetic Fields and Induced Temperatures Within a Model of the MicrowaveIrradiated Human Eye[J]. IEEE Transactions on Microwave Theory & Techniques， 1975， 23（11）： 888896.