时域有限差分法

2024-06-16

时域有限差分法（精选7篇）

时域有限差分法篇1

(一) 引言

目前, 电磁场的时域计算方法越来越引人注目。时域有限差分法 (Finite Difference Time Domain, FDTD) 作为一种主要的电磁场时域计算方法, 最早是在1966年由K.S.Yee提出的。这种方法通过将Maxwell旋度方程转化为有限差分式而直接在时域求解, 通过建立时间离散的递进序列, 在相互交织的网格空间中交替计算电场和磁场。经过三十多年的发展, 这种方法已经广泛应用到各种电磁问题的分析之中。

Matlab作为一种工程仿真工具得到了广泛应用。用于时域有限差分法, 可以简化编程, 使研究者的研究重心放在FDTD法本身上, 而不必在编程上花费过多的时间。

下面将采用FDTD法, 结合Matlab的图形功能来模拟自由空间中目标的散射近场, 说明了将二者结合起来的优越性。

(二) FDTD方法的基本原理

时域有限差分法的主要思想是把Maxwell方程在空间、时间上离散化, 用差分方程代替一阶偏微分方程, 求解差分方程组, 从而得出各网格单元的场值。FDTD空间网格单元上电场和磁场各分量的分布如图1所示。

电场和磁场被交叉放置, 电场分量位于网格单元每条棱的中心, 磁场分量位于网格单元每个面的中心, 每个磁场 (电场) 分量都有4个电场 (磁场) 分量环绕。这样不仅保证了介质分界面上切向场分量的连续性条件得到自然满足, 而且还允许旋度方程在空间上进行中心差分运算, 同时也满足了法拉第电磁感应定律和安培环路积分定律, 也可以很恰当地模拟电磁波的实际传播过程。

1. Maxwell方程的差分形式

Maxwell旋度方程为:

令f (x, y, z, t) 代表E或H在直角坐标系中某一分量, 在时间和空间域中的离散形式取以下符号表示:

其中:i, j, k和n为整数。对于二维问题, 设所有物理量均为与z坐标无关, 即∂/∂z=0, 于是将 (1) 式利用二阶精度的中心有限差分式来表示函数对空间和时间的偏导数, 下面以TE波为例, 即可得到如下FDTD基本差分式:

上式中的参数Ca、bC、aD、bD的定义如下:

其中:ε, μ, σ, σm是各点处的电磁参数。在 (3) 式的推导过程中假设时间步进Δx=Δy=δ。由于计算机存储时数组下标是整数, 所以对 (3) 式进行了修改, 修改后的更新方程如下:

利用TM波和TE波之间的对偶关系, 即可编写统一适用于TE波和TM波情况的二维FDTD计算程序。

2. 数值色散及稳定性条件

为了减小数值色散, 在选取空间网格尺寸δ时, 应满足λmin≥10δ, δ=min (Δx, Δy, Δz) , 其中λmin是被研究媒质空间的最小波长值。由此可以看出:减小网格尺寸可以减小数值色散, 但是会引起计算存储量的增大, 因此需综合考虑, 权衡处理。

为了使数值计算稳定, 时间步长的选择应满足:

(三) 自由空间中目标散射近场的可视化模拟

1. 自由空间中未加散射体时的散射近场

FDTD模拟参数为:计算区域大小150×150 (网格数) , 正中的130×130为总场区, 其外是散射场区。取空间网格大小δ=λ/40, 时间步进Δt=δ/2c。平面波从左边入射, 入射角度为0°。图2是自由空间中未加散射体时的散射近场。

从图中我们得到:在自由空间内没有散射体的情况下, 整个空间就没有散射场。入射波只存在于总场区内, 散射场区域内没有电磁波。即入射波完全由总场散射场边界产生, 又由总场散射场边界吸收, 且在总场区域内也保持了平面波传播的特性, 这与事实相符。由此可见, 该方法的有效性。

2. 自由空间中金属圆柱的散射近场

FDTD模拟参数为:计算区域大小150×150 (网格数) , 正中的130×130为总场区, 其外是散射场区。圆形金属散射体的圆心位于网格正中心, 圆柱半径r=λ, λ=1×10-2m, 取空间网格大小δ=λ/40, 时间步进Δt=δ/2c, 可见圆柱半径为40个网格。平面波从左边入射到目标, 入射角度为0°。图3是位于自由空间中金属圆柱在不同时刻的散射近场可视化模拟结果。

3. 自由空间中金属方柱的散射近场

FDTD模拟参数与上述第1点相同, 只是将散射体从金属圆柱换成金属方柱。图4是位于自由空间中金属方柱在不同时刻的散射近场可视化模拟结果。

由于散射体是金属, 其散射性很强, 从图3、图4中可以看出圆柱和方柱散射体内部均没有场, 遵循了金属散射体的电磁规律。图中总场边界 (入射波在此边界引入) 处场为不连续, 这是因为在总场边界以内为总场, 以外为散射场, 目标位于总场边界内部。通过这两个算例可以看出FDTD在可视化近场时非常实用。

(四) 自由空间中目标的远场RCS

1. 自由空间中金属圆柱的远场RCS

2. 自由空间中金属方柱的远场RCS

图5和图6给出了不同形状散射体的双站RCS (用波长归一化) , 包括TE波和TM波情况。所得结果与参考文献[1]作了对比, 可以看出一致性很好。

(五) 结语

以上结合FDTD和Matlab对自由空间中目标的近场散射做了仿真分析, 所编Matlab程序简洁明了, 运行效率也较高。FDTD法在电磁场数值分析方面有很大的优越性, 而Matlab具有强大的数据处理和图形处理功能, 可以快速地编出高效高质量的程序。将二者的优势有效地结合起来, 可以将算法迅速程序化, 并获得很好的数据处理结果, 使研究者可以集中精力在FDTD方法和研究对象本身上, 而只需花费少量的时间在程序的实现上。

参考文献

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

[2]Yee K S.Numerical solution of initial boundary value problems involving Maxwell equations in isotropic media[J].IEEE Trans, Antennas Propagat, May1996, AP-14 (3) :302-307.

[3] (美) Duane Hanselman, Bruce Littlefield著.精通Matlab7[M].朱仁峰, 译.北京:清华大学出版社, 2006.5.

[4]杨臻颖.三维散射的FDTD模拟及平面入射波补偿若干问题的研究[J].2006.

时域有限差分法篇2

一类四阶半线性抛物型方程的有限差分方法

研究了一类四阶半线性抛物方程,对其提出有限差分格式,并进行收敛性分析,得到L2范数下的误差估计.

作者：纪维强杨青 JI Wei-qing YANG Qing 作者单位：山东师范大学数学科学学院,济南,250014刊名：科学技术与工程 ISTIC英文刊名：SCIENCE TECHNOLOGY AND ENGINEERING年，卷(期)：8(23)分类号：O241.82关键词：四阶抛物方程有限差分法误差估计

时域有限差分法篇3

1 天线结构与理论计算

同轴线馈电的微带天线建模需要考虑同轴线的建模问题。对于此问题的处理可以采用阶梯近似法和一维传输线模型法[4—8]。阶梯近似法首先将同轴线分解成内、外导体及中间的介质三部分。然后对各个部件分别建立几何外形及尺寸的描述文件并进行FDTD剖分。因此接地板下方需要增加同轴线的模型空间, 同轴线的建模无疑增加了计算量。一维传输线模型法虽然不用对同轴线进行建模, 但也需要同轴线进行处理。需要将同轴线划分为接地板下的同轴线部分和馈电探针部分。同轴线部分采用一维电压V和电流I进行迭代馈电提供激励源[9]。接地板上方的馈电探针部分采用细导线模型。但是, 在同轴线与接地板的连接处需要考虑同轴线与接地板口径处的耦合问题[10]。一维传输线模型法虽然减少了计算空间, 但需要对部分空间进行特殊处理。

本文针对一维传输线模型法进行改进, 将馈电探针与贴片相连接, 探针下方与接地板之间留有一个网格作为激励源的位置, 采用点源作为激励。这样接地板与探针之间就不必采用耦合方式馈电, 于是省略了一维传输线模型的迭代。新方法建模如图1所示。

1.1 同轴探针处理

馈电探针处的处理采用细导线模型, 设细导线在 (ia, ja, k) 处。图2给出了在FDTD网格中半径为a (a<0.5Δx) 轴线与Ez分量重合的细导线。对于细导线附近的场量Hy、Ex, Hx、Ey在细导线附近有1/r的依赖关系, 以Hy、Ex为例:

对环绕细导线区域应用法拉第电磁感应定律:, 将上式带入, 得Hy的表达式为

同样的方法可以求出其他场分量。由于模型中简化了同轴线, 所以不必考虑探针处同轴线端口场量的耦合。

1.2 激励源

为了得到天线的宽频带信息激励源Vs选择高斯脉冲。激励源在细导线与接地板之间的一个网格处以电压源形式加入。电压源内阻R为50Ω, 设介质无损耗, 则激励源处的迭代式为:

1.3 电压电流取样

为了求得同轴线内电压、电流的波形, 可在探针处进行电压与电流的取样。设探针与接地板之间的电压V, 电流I;

式 (4) 中, k1, k2为接地板到探针之间的k坐标。因为探针与接地板间只有一个网格, 则上式简化为

对磁场应用安培环路定理:

值得注意的是以上的取样信息是在时间上相差半个时间步。通过得到的V、I的时域波形可以进一步求得微带天线的反射系数、输入阻抗等频域信息。

2 算例

采用以上方法编程分析了图1所示同轴探针馈电的微带天线。天线结构参数如下:贴片lx×ly=36.6 mm×26 mm, 介质厚度h=1.58 mm, 相对介电常数εr=2.17, 同轴探针半径r=0.24 mm, 馈电点位置在中线上lx, ds=6.5 mm, FDTD网格尺寸Δx=1.83 mm, Δy=1.625 mm, Δz=0.79 mm, 时间步长Δt=Δz/2c, 其中c为光速。吸收边界采用10层卷积完全匹配层 (CPML) 。图3为该微带天线的S11曲线。从图中可以看出计算结果与参考文献结果吻合很好, 说明此方法简化分析同轴馈电微带天线的可行性。图4、图5分别为微带天线的输入阻抗和方向图, 进一步说明在对于同轴线的简化问题上并不影响天线的其他参数的计算。

3 结论

简化同轴馈电天线模型同时结合同轴线的细导线分析方法分析了探针馈电的微带天线。数值结果和参考文献吻合很好, 说明这种同轴线的简化建模方式是可行的。建模过程中由于省略了对同轴线一维模型或阶梯近似方法的引入, 使得分析过程更简洁, 编程更容易。而且简化同轴线的思想并不影响天线的阻抗及辐射特性, 使得此方法更具实用性。

摘要：提出了一种简化的同轴馈电模型用于分析同轴探针馈电的微带天线。在时域有限差分方法 (FDTD) 中对于同轴馈电模型一般采用一维电压电流递推方式引入入射波形;同时结合细导线模型模拟同轴线馈电。需要对接地板同轴线接口处的电场值进行特别处理, 过程较复杂。简化方法仅采用细导线模型作为馈电端口, 在接地板与细导线探针间加电压源激励产生入射波形。省略了同轴线的建模及接口处电场值的处理过程, 使得模型更加紧凑。数值仿真验证了该算法的有效性和准确性。

关键词：时域有限差分法 (FDTD) ,同轴馈电,细导线,微带天线

参考文献

[1] 李增瑞, 王均宏, 蒋开波, 等.采用FDTD法研究微带天线雷达散射截面的减缩.电子学报, 2007;35 (6) :1065—1068

[2] 何芒, 徐晓文.探针馈电圆柱共形微带天线的FDTD法分析.微波学报, 2004;20 (3) :1—5

[3] Wong K L, Design of nonplanar microstrip antennas and transmission line.John Wiley&Sons, Inc., 1999

[4] 陈文俊, 黎滨洪, 谢涛.FDTD法分析探针加载微带天线的简便方法.电波科学学报, 2006;21 (1) :146—149

[5] 许锋, 洪伟, 冯祖伟.同轴馈电耦合微带贴片天线的时域有限差分法分析.电波科学学报, 2001;16 (1) :34—38

[6] 袁峰, 卢春兰.基于FDTD法分析计算同轴馈电贴片天线.军事通信技术, 2007;28 (2) :12—15

[7] 姜永金, 潘谊春, 傅文斌.同轴馈电耦合微带天线的MPSTD—FDTD分析.微波学报, 2011;27 (2) :19—23

[8] Supriyo D, Mitrrra R.A locally conformal finite-difference time-domain (FDTD) algorithm for modeling three-dimensional perfectly conducting objects.IEEE Trans AP, 1997;45 (7) :273—276

[9] Maloney J G, Shlager K L, Smith G S.A simple FDTD model for transient excitation of antennas by transmission lines.IEEE Trans on AP, 1994;42 (2) :289—292

时域有限差分法篇4

波动方程时域有限差分法(thewave-equationfinite-differencetime-domainmethod,WEFDTD)是电磁仿真领域的一种重要的计算方法[1],但目前WEFDTD的激励源却主要是“强迫源”,强迫源具有以下不足:(1)存在虚假反射,影响了计算结果的正确性,虽然人们采取了一些措施来减小虚假反射,但效果不很理想;(2)不能模拟散射场,而在很多问题中必须考虑散射场,此外由于吸收边界条件(Absorbingboundarycondition,ABC)只吸收外向的散射波,不能模拟散射场将导致很多情况下吸收边界条件不能使用。

时域有限差分法(thefinite-differencetime-domainmethod,FDTD)是电磁仿真领域的另一重要的计算方法[2,3,4,5,6,7,8],这一方法的“总场-散射场源”不存在“强迫源”的不足,但总场-散射场源能否用于WEFDTD,如果能又如何应用,目前还未见文献报道。

基于这种情况,首先利用等效原理,从理论上论证了总场-散射场源可以应用于WEFDTD,然后又研究了总场-散射场源在WEFDTD中的具体实现方法,为WEFDTD找到了一种新激励源,有效地克服了WEFDTD已有激励源技术的不足。

1 理论基础

根据叠加原理,空间场可以写成入射场和散射场之和,即

式(1)中E in、Hin为入射场,E sc、Hsc为散射场。如图1所示,通常把仿真区划分为总场区和散射场区,在总场区只存储总场值,在散射场区只存储散射场值。

设波如图2所示入射。为了使入射波限制在总场区,根据等效原理,在两区域的分界面A上设置面电磁流,并设A面外的场为零,如图3所示,因而A面上的等效电磁流为

式(2)中→en为面A的外法向单位矢量。所以在两区域的分界面上设置入射波的切向分量便可将入射波只引入到总场区。

2 程序实现

为简洁起见,以二维情形为例来说明总场-散射场源在WEFDTD中的程序实现方法。二维情形的波动方程为

取WEFDTD的差分网格如图4所示,把式(3)在点(i,j)和时刻n作中心差分离散,得

从中可解得

式(5)即WEFDTD的差分格式[9],适用于总场区和散射场区,但在不同的区域,U分别代表总场和散射场。

由于电磁波的波动方程是线性方程,对总场和散射场都适用,所以由电磁波的波动方程离散而得的式(5)也是对总场和散射场都适用,但在不同的区域,U分别代表总场和散射场。

但式(5)并不适用于两区域的分界面,这是因为用式(5)刷新某一点的场值时,须用到其周围四个点的场值,而且该点和其周围四个点上存储的必须同为总场值,或同为散射场值,对位于两区域分界面上的点,不可能做到这一点。这也是在WEFDTD中应用总场-散射场源的难点所在。

根据总场-散射场原理,如果能在两区域的边界上设置入射场,且对式(5)做一些修改,则不但可使式(5)适用于两区域的分界面,而且可把入射场只引入到总场区,从而实现总场-散射场源。作这些修改的基本考虑是:刷新某一总场值时,如果环绕它的有散射场值,则令该散射场值加上入射场(这样就等于总场了)。刷新某一散射场值时,如果环绕它的有总场值,则令该总场值减去入射场(这样就等于散射场了)。这样处理就可以满足式(5)只对总场或散射场适用的要求了。以下就一种常见情形来说明这一方案。

设总场区是范围为(iL≤i≤iR,jB≤j≤jA)的矩形,则入设场需设置在该矩形以及散射场区最邻近该矩形的一个矩形上。根据通常的做法,规定该矩形上存储总场值,则在该矩形的(i=iL,jB+1≤j≤jA-1)的矩形边上,式(5)需修改为

式(6)中为入射场在点(i-1,j)处的场值。在散射场区最邻近该矩形的(i=iL-1,jB≤j≤jA)的矩形边上,式(5)需修改为

式(7)中为入射场在点(i+1,j)处的场值。

两矩形其它各边或顶点上的差分格式可类似写出,为节省篇幅,此处不再给出。

此外还应注意,入射波沿着其传播方向,各点起振时刻依次要落后一些,因此入射波在各点的表达式一般并不相同。例如,设入射波沿x轴正方向传播,其在(i=iL-1)处的表达式为U=f(t),则在(i=iL)处,其表达式应为

式(8)中Δx为x方向的空间步长,为波从i=iL-1处传播到i=iL处所用的时间。

3 数值实验

为了验证以上总场-散射场源的有效性,编程进行了数值实验。程序中取空间步长Δx=Δy=0.05m,时间步长,仿真区域的范围为(0≤i≤200,0≤j≤150),总场区的范围为(50≤i≤150,20≤j≤130)。吸收边界条件采用WEFDTD的完全匹配层(perfectlymatchedlayer,PML)吸收边界条件[10],吸收层的厚度为6个网格,电导率按抛物线规律增长。

简谐波和高斯脉冲是电磁仿真中使用最多的两种信号源。程序中使用的即是这两种信号源,并设入射波沿x轴正方向行进,在i=49界面上,信号源的表达式分别取为

和

式(9)及式(10)中λ=1.0m,n0=50,ndecay=20。

图5和图6给出了总场区中一点(100,75)在分别使用两种信号源时的模拟结果,可以看出,计算值和理论值符合得很好。此外根据总场-散射场原理,此时散射场区的场值应为零,我们记录了散射场区一点(25,75)在各时刻的场值,发现其绝对值不超过0.02,非常接近于零。这表明总场-散射场源在WEFDTD中是有效的。

4 结论

总场-散射场源可以应用于WETDTD,但其具体的程序实现和在FDTD中不同。本文研究并解决了这些问题,为在WETDTD中应用总场-散射场源奠定了理论基础。此外由于所有波的波动方程都是一样的,本文的研究成果也可推广用于其它波如声波、地震波等的数值仿真。

摘要：研究了“总场-散射场源”在波动方程时域有限差分法(WEFDTD)中的应用方法,这一问题的难点是如何把基本差分格式应用于总场区和散射场区的分界面。具体做法是将激励源设置在两区域的分界面上,并基于总场-散射场原理对基本差分格式加以修正,使得用其刷新某一网格点的总场值或散射场值时,该点周围四个网格点能同时提供总场值或散射场值。这一技术有效地克服了WEFDTD的已有激励源的不足。使用这一技术,数值模拟了简谐波和高斯脉冲的传播过程,发现计算结果和解析解符合得很好。

关键词：波动方程时域有限差分法(WEFDTD),激励源,总场,散射场

参考文献

[1]喻志远,林为干.波动方程二阶差分方法的研究及在矩形波导中的应用.电子科技大学学报,1999;28(1):37—40

[2]李晓波,陈曦.基于APML-FDTD方法的指挥所UWB信道分析.科学技术与工程,2009;9(12):3 283—3 287

[3]杨勇,车小花,张菲,等.用三维时域有限差分法研究随钻声波测井仪器隔声体的设计.科学技术与工程,2009;9(3):565—567

[4] Yee K S.numerical solution of initial boundary value problems invol-ving maxwell’s equations in isotropic media.IEEE Trans AP,1966;14(3):302—307

[5] Taflove A,Hagness S C.Computational Electrodynamics:The Finite-Difference Time-Domain Method.Boston.London:Artech House,2000:194—212

[6] Namiki T.3-D ADI-FDTD method:unconditionally stable time-do-main algorithm for solving full vector Maxwell’s equations.IEEETrans Microwave Theory Tech,2000;48(10):1743—1748

[7] Wang H,OuYang Z B,Han Y L,et al.An FDTD model of photoniccrystal laser with four-energy-lever Gain Media.Chinese Journal ofComputational Physics,2007;24(4):452—456

[8]王秉中.计算电磁学.北京:科学出版社,2002:20

[9]王建永,李庆武,刘国高.一般正交曲线坐标系中的WEFDTD.计算机仿真,2007;24(2):69—71

时域有限差分法篇5

一、一维FDTD差分方程表达形式

一维波动方程, 在一维无限长介质中, 其解为u=f (x-vt) , 是以速度v向x轴正向传播的行波。用差分形式, 把波动方程进行改写:

二、初始、边界条件的确定以及编程计算的结果

边界条件:u (t, 1) =0, u (t, Nx) =0, 对t=1到Nt。

初始条件:由于在差分形式的递推公式 (1) 中, 要计算t=3的值, 需要知道t=1, 2时的值。本文规定初始条件为u (t, i) =0, 对t=1, 2.

下面是两组不同的参数, 用Matlab编程计算的结果:

第1组:Lx=1×10-4, Δx=1×10-6, Nx=100, Nt=100, v=c=3×108,

信号源, 观察点在i=80, 也就是x=80Δx处。

第2组:Lx=1×10-4, Δx=1×10-6, Nx=100, Nt=1000, v=c=3×108,

信号源, 观察点在i=80, 也就是x=80Δx处。参数中时间长度Nt与第1组不同, 其他参数相同。

第3组:Lx=1×10-4, Δx=1×10-7, Nx=1000, Nt=3000, v=c=3×108,

信号源, 观察点在i=500, 也就是x=500Δx处。

从图中可以看出:

t=1到600, 在观察点无信号, 信号还没有到达;

t=600到900, 是过渡阶段, 与数值化的计算有关, 在观察点信号尚不稳定;

t=900到2400, 观察点信号与用连续方程得到的结果一致;

t=2400以后, 观察点的信号的幅度比源的信号还大, 是因为文中取一维两端边界为零, 反射波到达观察点, 与正向波形成干涉。

用快速傅里叶变换得出的频谱透射率, 由于已经包含了开始的过渡阶段和以后的反射波的干涉的影响, 在0.8左右, 与波的无衰减传播结果接近。

三、结论

在一维的FDTD计算中, 直接使用整数时刻和整数空间坐标, 而不需要使用半整数时间和半整数空间坐标, 简化了坐标的标记形式;时间步长取, 时间总长Nt的取值, 要求在这个时间信号源的信号可以到达观察点, 但又没有界面反射波到达观察点, 这样的参数比较理想, 可以得到与连续方程相同的结果。

摘要：在用一维时域有限差分方法计算波的传播时, 为了防止出现发散的结果, 要求时间步长Δt<Δx/v, 但具体取多大的值, 没有定论。通过编程试验, 时间步长取Δt<0.75Δx/c, 而时间步数Nt大于波从波源传到观察点所需要的时间, 又小于反射波到达观察点的时间。而且计算的时间格点和位置格点全部都可以取整数。

关键词：FDTD,有限时域差分方法,参数

参考文献

[1]葛德彪, 阎玉波.电磁波时域有限差分方法[M].第2版.西安:西安电子科技大学出版社, 2005.

[2]曾辉, 杨亚培.FDTD法与平面波展开法在光子禁带计算中的差异分析[J].电子科技大学学报, 2005, 34 (6) :901-904.

[3]张国华, 袁乃昌, 付云起.FDTD方法分析光子带隙能带结构[J].微波学报, 2001, 17 (4) :14-17.

[4]张亮, 寇晓艳.FDTD的Matlab语言的实现[J].延安大学学报:自然科学版, 2009, 28 (2) :57-59.

时域有限差分法篇6

时域有限差分方法[finite-difference time-domain (FDTD) method]自1966年由Yee[1]提出来后, 一直是计算电磁学常用的方法之一, 并广泛应用于电磁散射、天线的分析与设计、雷达截口的计算等电子工业与国防工业。但是, 由于FDTD方法是条件稳定的, 即在二维情形下, 时间和空间步长分别为Δt, Δx, Δy;必须满足Courant-Friedrichs-Lewy (CFL) 稳定条件: cΔt≤[1/ (Δx) 2+1/ (Δy) 2]-1/2。其中c是在介质中的光速, ε, μ是介电常数和磁导率。为了解决上述问题, 由文献[2,3]提出了交替方向隐式时域有限差分方法 (ADI-FDTD) 方法, 并用Fourier方法证明了这种格式是无条件稳定的。近年来[4], 将分裂算子方法与FDTD方法结合提出新颖而且简单的分裂算子的时域有限差分方法 (S-FDTD) , 与二维ADI-FDTD相比都是无条件稳定的二阶格式, 但是S-FDTD格式采用算子分裂技术, 所以格式简单, 计算时间短, 并且在模拟一种散射问题时, S-FDTD比ADI-FDTD更精确。麦克斯韦方程的高阶时域有限差分方法已有较多的研究工作[4,5,6,7,8,9], 但是这些方法都没有使用分裂算子技术。

本文将高阶差分与S-FDTD方法相结合提出一种新的方法, 即在分裂后的方程基础上对空间采取四阶中心差分, 从而得到分裂 (splitting) 的高阶 (high-order) 的时域有限差分 (FDTD) 格式SHO-FDTDⅠ;并在此格式的基础上通过增加扰动项减小分裂误差给出它的修正格式SHO-FDTDⅡ, 给出了具体的计算步骤, 然后用Fourier方法证明了这两种格式是无条件稳定的, 并给出这两种格式的数值弥散关系式和数值弥散误差的计算。与条件稳定的高阶FDTD相比, 新方法是无条件稳定的;与高阶的ADI-FDTD相比, 新格式更简单, 易于编程和实现。

1 二维麦克斯韦方程的SHO-FDTD格式

考虑以下二维横向电波:

$\frac{\partial E_{x}}{\partial t} = \frac{1}{ε} \frac{\partial Η_{z}}{\partial y}$ (1.1)

$\frac{\partial E_{y}}{\partial t} = - \frac{1}{ε} \frac{\partial Η_{z}}{\partial x}$ (1.2)

$\frac{\partial Η_{z}}{\partial z} = \frac{1}{μ} (\frac{\partial E_{x}}{\partial y} - \frac{\partial E_{y}}{\partial x})$ (1.3)

其中 $\vec{E} = (E_{x} (x, y, t), E_{y} (x, y, t))$ 表示电场, $Η_{z} = Η_{z} (x, y, t)$ 表示磁场, ε和μ分别是介电常数和磁导率, Ω=[0, a]×[0, b], t ∈ (0, T]。其边界满足理想导体边界条件: $E_{x} (x, 0, t) = E_{x} (x, b, t) = E_{y} (0, y, t) = E_{y} (a, y, t) = 0$ 。初始条件: Ex0=Ex0 (x, y) , Ey0=Ey0 (x, y) , Hz0=Hz0 (x, y) 。

为了书写的方便, 我们仅讨论常系数的情形, 这里所用的方法容易推广到变系数的情形。对空间区域Ψ以及时间区域[0, T]采取如下剖分:

其中Δx和Δy分别是沿x轴方向和沿y轴方向的空间离散步长, Δt是时间步长, I, J, N为整数。

对于任意一个网格函数F (t, x, y) , 引入以下记号:

F $_{α, β}^{m}$ =F (mΔt, αΔx, βΔy) , $δ_{t} F_{α, β}^{m} = \frac{F_{α, β}^{m + \frac{1}{2}} - F_{α, β}^{m - \frac{1}{2}}}{Δ t},$

$δ_{x}^{2} F_{α, β}^{m} = \frac{1}{Δ x} (\frac{1}{24} F_{α - \frac{3}{2}, β}^{m} - \frac{9}{8} F_{α - \frac{1}{2}, β}^{m} + \frac{9}{8} F_{α + \frac{1}{2}, β}^{m} - \frac{1}{24} F_{α + \frac{3}{2}, β}^{m})$ ,

$δ_{y}^{2} F_{α, β}^{m} = \frac{1}{Δ y} (\frac{1}{24} F_{α, β - \frac{3}{2}}^{m} - \frac{9}{8} F_{α, β - \frac{1}{2}}^{m} + \frac{9}{8} F_{α, β + \frac{1}{2}}^{m} - \frac{1}{24} F_{α, β + \frac{3}{2}}^{m})$ 。

用E $_{v_{α, β}}^{m}$ (v=x, y) 表示电场Ev (tm, xα, yβ) 的近似值。H $_{z_{α, β}}^{m}$ 表示磁场Hz (tm, xα, yβ) 的近似值。将分裂算子的时域有限差分方法[9]与高阶差分方法[4]相结合, 提出如下格式:

第一步:

$\frac{E_{y_{i, j + \frac{1}{2}}}^{n + 1} - E_{y_{i, j + \frac{1}{2}}}^{n}}{Δ t} = \frac{1}{2 ε} δ_{x}^{2} = (Η_{z_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} + Η_{z_{i, j + \frac{1}{2}}}^{n})$ (1.4a)

$\frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}}{Δ t} = - \frac{1}{2 μ} δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.4b)

第二步:

$\frac{E_{x_{i + \frac{1}{2}, j}}^{n + 1} - E_{x_{i + \frac{1}{2}, j}}^{n}}{Δ t} = \frac{1}{2 ε} δ_{x}^{2} = (Η_{z_{i + \frac{1}{2}, j}}^{n + \frac{1}{2}} + Η_{z_{i + \frac{1}{2}, j}}^{n})$ (1.5a)

$\frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}}{Δ t} = \frac{1}{2 μ} δ_{y}^{2} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.5b)

此格式称为分裂算子 (splitting) 的高阶 (high-order) 时域有限差分方法, 记为SHO-FDTDⅠ。其中, 边界条件为: $E_{x_{i + \frac{1}{2}, 0}}^{n} = E_{x_{i + \frac{1}{2}, j}}^{n} = E_{y_{0, j + \frac{1}{2}}}^{n} = E_{y_{i, j + \frac{1}{2}}}^{n} = 0$ 。初始值:

求解步骤为:先解第一步, 由式 (1.4b) 得:

$Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} = Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n} - \frac{Δ t}{2 μ} δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.6)

将式 (1.6) 代入式 (1.4a) 整理得:

显然这是一个七对角的方程组, 系数矩阵元素都是常数而且在每一个时间层上都不变。因此可由一些求解线性方程组的方法如:超松弛迭代法, 共轭梯度法等, 求出 ${E_{y_{i, j + \frac{1}{2}}}^{n + 1}}$ 。然后, 将 $E_{y_{i, j + \frac{1}{2}}}^{n + 1}$ 代入式 (1.6) 显式求出 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ 。式 (1.5) 的计算与第一步相似, 但要用到第一步算出的 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ 。方程 (1.7) 中的第一个方程和最后一个方程的系数需要调整, 可根据文献[5]中的方法进行, 这里略去。

检查格式SHO-FDTDⅠ的截断误差发现关于时间它的精度不高。为了提高精度, 引入一个扰动项, 得到修正格式SHO-FDTDⅡ如下:

第一步:

第二步与SHO-FDTDⅠ中的第二步相同.

这种格式的初边值条件与式 (1.4) 、式 (1.5) 的初边值条件完全相同, 求解方法也与前一格式相同。为了了解格式SHO-FDTDⅡ的逼近精度, 由式 (1.5a) 、式 (1.5b) 、式 (1.8a) 、式 (1.8b) 消去中间项 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ , 可得到SHO-FDTDⅡ的等价格式:

$δ_{t} E_{x_{i + \frac{1}{2}, j}}^{n + \frac{1}{2}} = \frac{1}{2 ε} δ_{y}^{2} (Η_{z_{i + \frac{1}{2}, j}}^{n + 1} + Η_{z_{i + \frac{1}{2}, j}}^{n}) (1.9 a)$

$\begin{array}{l} δ_{t} E_{y_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} = - \frac{1}{2 ε} δ_{x}^{2} (Η_{z_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} + Η_{z_{i, j + \frac{1}{2}}}^{n}) + \\ \frac{Δ t}{4 μ ε} δ_{x}^{2} δ_{y}^{2} (E_{x_{i, j + \frac{1}{2}}}^{n + 1} - E_{x_{i, j + \frac{1}{2}}}^{n}) (1.9 b) \end{array}$

$\begin{array}{l} δ_{t} Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} = \frac{1}{2 μ} [δ_{y}^{2} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}) - \\ δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})] (1.10) \end{array}$

格式SHO-FDTDⅠ的等价格式与式 (1.9a) 、式 (1.9b) 、式 (1.10) 相似, 只要将式 (1.9) 最右端后一项中的“-”改为“+”。由此可以看出, SHO-FDTDⅡ的摄动误差 (二阶) 比SHO-FDTDⅠ的摄动误差 (二阶) 要高一阶, SHO-FDTDⅡ是关于时间2阶, 关于空间是4阶的式 (2.4) 格式, SHO-FDTDⅠ则是式 (1.4) 格式。

2 稳定性分析与数值弥散分析

用Fourier方法分析这两种格式的稳定性。假定格式的差分解具有下列形式。

$\vec{E}_{α, β}^{n} = \vec{E}_{0} ξ^{n} e^{- i (k_{x} α Δ x + k_{y} β Δ y)}$ ;

H $_{z_{α, β}}^{n}$ =Hz0ξne-i (kxαΔx+kyβΔy) 。

其中 $i = \sqrt{- 1}$ , $\vec{E}_{0} = (E_{x_{0}}, E_{y_{0}})^{Τ}$ 和Hz0是振幅, $\vec{k} = (k_{x}, k_{y})$ 是波矢量, $| k | = \sqrt{k_{x}^{2} + k_{y}^{2}}$ , ξ是增长因子。下面将E $_{α, β}^{n}$ , H $_{z_{α, β}}^{n}$ 代入式 (1.5a) 、式 (1.9) 、式 (1.10) , 化简得:

$(ξ - 1) E_{x_{0}} - i \frac{Δ t}{ε} (ξ + 1) b_{y} Η_{z_{0}} = 0$ (2.1)

$\begin{array}{l} (ξ - 1) E_{y_{0}} + \frac{(Δ t)^{2}}{ε μ} (ξ - 1) a_{x} b_{y} E_{x_{0}} + \\ i \frac{Δ t}{ε} (ξ + 1) a_{x} Η_{z_{0}} = 0 (2.2) \end{array}$

$\begin{array}{l} (ξ - 1) Η_{z_{0}} - i \frac{Δ t}{μ} (ξ + 1) b_{y} E_{x_{0}} + \\ i \frac{Δ t}{μ} (ξ + 1) a_{x} E_{y_{0}} = 0 (2.3) \end{array}$

其中 $a_{x} = \frac{(\frac{1}{12} \sin \frac{3}{2} k_{x} Δ x - \frac{9}{4} \sin \frac{1}{2} k_{x} Δ x)}{2 Δ x}, b_{y} = \frac{(\frac{1}{12} \sin \frac{3}{2} k_{y} Δ y - \frac{9}{4} \sin \frac{1}{2} k_{y} Δ y)}{2 Δ y}$ 。由于向量 (Ex0, Ey0, Hz0) 是非零向量, 所以关于Ex0, Ey0, Hz0的方程组的系数矩阵的行列式的值为零。整理得:

(ξ-1) (d0ξ2+2d1ξ+d0) 2=0 (2.4)

式 (2.4) 中

$d_{0} = 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2};$

$d_{1} = - 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}$ 。

解方程式 (2.4) 得: $ξ_{1} = 1, ξ_{2} = d_{0}^{- 1} (- d_{1} + i \sqrt{d_{0}^{2} - d_{1}^{2}}) 。 ξ_{3} = d_{0}^{- 1} (- d_{1} - i \sqrt{d_{0}^{2} - d_{1}^{2}})$ 。显然方程三个根的模都是1, 所以格式SHO-FDTDⅡ是无条件稳定的。对于格式SHO-FDTDⅠ, 同理可得:

c3ξ3+c2ξ2+c1ξ+c0=0 (2.5)

式 (2.5) 中

$\begin{array}{l} c_{3} = 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}, \\ c_{2} = - 3 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + 3 \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}; \end{array}$

$\begin{array}{l} c_{1} = 3 - \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + 3 \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}, \\ c_{0} = - 1 - \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2} 。 \end{array}$

方程 (2.5) 的根的表达式非常复杂 (为了简单, 这里省略) , 通过对根的表达式的研究发现 $| ξ | = 1 + Ο (Δ t)$ , 因此格式SHO-FDTDⅠ是耗散的无条件稳定的。为了直观的感受 $| ξ |$ , 下面给出它在不同情况下的变化曲线图。设kx=kcos φ。ky=ksin φ。其中k, φ。是向量k的圆柱坐标, 再设 $Ν_{λ} = \frac{λ}{h}, ω = c k, ω$ 是频率, λ是波长, Δx=Δy=h, 则Nλ是一个波长内的节点个数。

令 $S = \frac{c Δ t}{h}$ , 由CFL的定义可知

$c Δ t \sqrt{\frac{1}{Δ x^{2}} + \frac{1}{Δ y^{2}}} = \frac{c Δ t}{h} \sqrt{2} = \sqrt{2} S$ 。因此, 不妨称S为CFL数。由方程 (2.4) 或式 (2.5) 及其系数的表达式可知, 增长因子ξ是关于S, φ, Nλ的函数。例如, 可推出 $\sin \frac{3}{2} k_{x} Δ x = \sin (\frac{3 π}{Ν_{λ}} \cos φ)$ 。方程 (2.5) 的根的表达式比较复杂。下面我们用Matlab求出它的近似解, 并画出不同情况下它的根的变化曲线。

图1是在S=0.35, Nλ=10情况下给出了根 $| ξ |$ 随角度φ的变化曲线。上边的曲线是方程复数根的模, 下边的曲线是实数根的绝对值。很容易看出复数根的模比1大, 但是 $| ξ | = 1 + Ο (Δ t)$ , 其中, $Δ t = S \times \frac{1}{Ν_{λ}} = 0.035$ 。从图1中显然可以看出曲线位于: y=1±0.035两条直线之间, 这就表明SHO-FDTDⅠ格式是无条件稳定的。

图2和图3分别是在S=1.5, φ=35°和Nλ=10, φ=65°情况下, 方程 (2.5) 的根的模分别随着Nλ和S的变化情况, 从这两个图中同样可以看出SHO-FDTD格式是无条件稳定的。

2.2 数值弥散关系

假设麦克斯韦方程的差分解为: $\vec{E}_{α, β}^{n} =$ $\vec{E}_{0} e^{i (k_{x} α Δ x + k_{y} β Δ y - ω n Δ t)}$ , H $_{z_{α, β}}^{n}$ =Hz0ei (kxαΔx+kyβΔy-ωnΔt) , 其中, ω, kx, ky, $\vec{E}_{0}$ 和Hz0与2.1节的符号相同。将E $_{α, β}^{n}$ , H $_{z_{α, β}}^{n}$ 代入SH-FDTDⅡ的等价格式 (1.5a) 、式 (1.9) 、式 (1.10) , 类似于2.1节对增长因子ξ的推导, 可得:

$\begin{array}{l} \sin^{2} (\frac{1}{2} ω Δ t) = (c Δ t)^{2} [a_{x}^{2} + b_{y}^{2} + (c Δ t)^{2} a_{x}^{2} b_{y}^{2}] \\ \cos^{2} (\frac{1}{2} ω Δ t) (2.6) \end{array}$

式 (2.6) 中 $c^{2} = \frac{1}{ε μ}$ , 式 (2.6) 就是SHO-FDTDⅡ的数值弥散关系式。

同理可得SHO-FDTDⅠ格式的数值弥散关系式:

$\begin{array}{l} \sin^{2} (\frac{1}{2} ω Δ t) = (c Δ t)^{2} [a_{x}^{2} + b_{y}^{2} + (c Δ t)^{2} a_{x}^{2} b_{y}^{2} \times \\ \cos (\frac{1}{2} ω Δ t) \sin^{- 1} (\frac{1}{2} ω Δ t)] \times \\ \cos^{2} (\frac{1}{2} ω Δ t) (2.7) \end{array}$

由于 $\lim_{x \to 0} a_{x} = \lim_{x \to 0} a_{x}^{2} = \frac{k_{x}^{2}}{4}, \lim_{x \to 0} b_{y} = \lim_{x \to 0} b_{y}^{2} = \frac{k_{y}^{2}}{4}$ ;所以当Δx, Δy, Δt趋于0时, 数值弥散关系式 (2.6) 、式 (2.7) 趋向于理想的数值弥散关系: ω2=c2[k $_{x}^{2}$ +k $_{y}^{2}$ ]。此外, 比较式 (2.6) 、式 (2.7) 容易看出SHO-FDTDⅠ和SHO-FDTDⅡ的格式在时间上的精度分别是一阶的和二阶的, 并且SHO-FDTDⅡ比SHO-FDTDⅠ的数值弥散误差要小得多。

2.3 数值弥散误差

下面通过试验求出格式SHO-FDTDⅠ和格式SHO-FDTDⅡ的数值弥散误差, 并与理论结果进行对比分析。令ξ=eiωΔt是方程式 (2.1) 、式 (2.4) 的根。ω=ωR+iωI, 则ξ=e-ωIΔt[cos (ωRΔt) +isin (ωRΔt) ]从而, $\tan (ω_{R} Δ t) = \frac{Ι_{m} (ξ)}{R_{e} (ξ)}$ ; 其中Im (ξ) , Re (ξ) 分别表示ξ的虚部与实部。波的数值相速vp与光速c的比的数值弥散误差: $\frac{v_{p}}{c} = \frac{ω_{R}}{k c} = \frac{1}{c k Δ t} \arctan \frac{Ι_{m} (ξ)}{R_{e} (ξ)} = \frac{Ν_{λ}}{2 π S} \arctan \frac{Ι_{m} (ξ)}{R_{e} (ξ)}$ , 其中S是CFL数, Nλ是一个波长内的节点数, k是波长值。

图4—图6给出了Vp/c在不同情况下的变化曲线。图4给出了SHO-FDTDⅠ与SHO-FDTDⅡ的Vp/c在S=0.35, Nλ=10的条件下随着角度φ的变化图, 从图4中可以看出SHO-FDTDⅡ格式的Vp/c更接近于1。图5给出了两种格式在S=2.4, φ=120°的条件下Vp/c随着Nλ的变化图, 从图5中可以看出SHO-FDTDⅡ格式的数值弥散误差小于SHO-FDTDⅠ的误差。图6是两种格式在φ=65°, Nλ=10的条件下Vp/c随着S的变化图, 从图6中可以看出随着S的增大, 数值弥散误差变得越来越大, 但是在所有的情况下, SHO-FDTDⅡ格式的数值弥散误差比SHO-FDTDⅠ格式的误差要小得多。

参考文献

[1] Yee K S.Numerical solution of initial boundary value problems invol-ving Maxwell’s equations in isotropic media.IEEE Trans Antennasand Propagation, 1966;14:302—307

[2] Zheng F, Chen Z, Zhang J.Toward the development of a three-dimen-sional unconditionally stable finite—difference time—domain method.IEEE Trans Microwave Theory Tech, 2000;48:1550—1558

[3] Namiki T.AnewFDTD algorithmbased on alternatingdirection implicitmethod.IEEE Trans Microwave Theory Tech, 1999;47:2003—2007

[4] Turkel E, Yefet A.Fourth order compact method for the Maxwell e-quations with discontinuous coefficients.Applied Numerical Mathe-matics, 2000;33:125—134

[5] Shang J S.Higher-order compact difference schemes for time depend-ent Maxwell’s equations.Journal of Computational Physics, 2003;153:312—333

[6] Xu L J, Yuan N C.高阶ADI-FDTD算法的数值色散分析.电子信息学报, 2005;27 (3) :1662—1665

[7] Xie Z Q, Zhang B, Chan C H.An explicit fourth-order staggered finitedifference time-domain method for Maxwell’s equations.Journal ofComputational and Applied Mathematics, 2002;147:75—98

[8] Manry C W, Broschat S L, Scneider J B.Higher-order FDTD meth-ods for large problems.J Applied Computational Electromagnetics So-ciety, 1995;10:17—29

时域有限差分法篇7

时域有限差分方法 (finite-difference time-domain (FDTD) method) 是由Yee在1966年[1]提出, 在计算电磁学领域中是一种非常有效的数值计算方法之一, 已广泛应用于电磁散射、天线、电磁兼容、生物电磁场以及电波传播等问题的计算与模拟中。然而这种方法是条件稳定的, 要受到Courant-Friedrichs-Lewy (CFL) 稳定性条件的限制 (即: $c Δ t \leq [\frac{1}{(Δ x)^{2}} + \frac{1}{(Δ y)^{2}}]^{- \frac{1}{2}}$ , 其中:c是光的速度) [2]。因此当要模拟的物体具有细微结构时, 为了准确模拟其电磁特性, 空间步长必须足够小。这时为了保证解的稳定性, 时间步长也要相应地取得很小, 将使计算的时间大大延长, 有时甚至不可实现。

为了克服这一限制条件, 无条件稳定的交替方向的时域有限差分方法 (ADI-FDTD) [3,4]被提出来, 并得到电磁学界的高度关注和广泛研究[5]。最近, 文献[6]运用分裂算子的方法提出了电导率为零的麦克斯韦方程的分裂时域有限差分方法 (S-FDTD) , 该方法比二维ADI-FDTD更简单、运算时间更短, 在模拟一种电磁散射问题时更精确。

考虑电导率非零的麦克斯韦方程的分裂时域有限差分方法。提出了S-FDTDI和S-FDTDII两种格式, 分析了它们的截断误差, 得出逼近精度, 其中S-FDTDI、S-FDTDII关于时间分别为一阶、二阶, 关于空间都是二阶的差分方法。用S-FDTDI、S-FDTDII及ADI-FDTD方法给出一种矩形波导问题的计算结果和分析。数值试验表明, 文献[2]中提及的优点在解决非零电导率的电磁问题中仍然存在。

1 麦克斯韦方程、分裂格式及求解过程

考虑如下的Maxwell方程

$\frac{\partial E_{x}}{\partial t} = \frac{1}{ε} (\frac{\partial Η_{z}}{\partial y} - σ E_{x})$ (1)

$\frac{\partial E_{y}}{\partial t} = \frac{1}{ε} (- \frac{\partial Η_{z}}{\partial x} - σ E_{y})$ (2)

$\frac{\partial Η_{z}}{\partial t} = \frac{1}{μ} (\frac{\partial E_{x}}{\partial y} - \frac{\partial E_{y}}{\partial x})$ (3)

其中E= (Ex (x, y, t) , Ey (x, y, t) ) 表示电场, Hz=Hz (x, y, t) 表示磁场, (x, y) ∈Ω=[0, a]×[0, b], t∈ (0, T]。∂Ω为Ω的边界, $\overset{⇀}{n}$ 为∂Ω的外法向量, ε为介质介电常数, μ为磁导系数, σ为电导率。并且满足如下的良导体边界条件:Ex (x, 0, t) =Ex (x, a, t) =0, Ey (0, y, t) =Ey (b, y, t) =0和初始条件:E0= (Ex0 (x, y) , Ey0 (x, y) , Hz0=Hz0 (x, y) ) 。

为了简单起见, 只考虑ε, μ, σ为常数时的情况, 此种方法可以推广到变系数的情形。

对Ω采用与Yee相同的离散方法称为交错网格剖分[1]。设 $Ι = \frac{a}{Δ x}, J = \frac{b}{Δ y}$ 为正整数, 令

设Δt为时间步长, $Ν = \frac{Τ}{Δ t}$ 为正整数, 对[0, T]进行等距剖分, 令

$\begin{array}{l} t^{n} = n Δ t, t^{n + \frac{1}{2}} = (n + \frac{1}{2}) Δ t, n = 0, 1, \dots, Ν - 1, \\ t^{Ν} = Ν Δ t = Τ 。 \end{array}$

为书写简单, 对函数F (x, y, t) , 令F $_{α, β}^{m}$ =F (mΔt, αΔx, βΔy) , 现定义一些差分算子

设 ${E_{x_{i + \frac{1}{2}, j}}^{m}}, {E_{y_{i, j + \frac{1}{2}}}^{m}}, {Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{m}}$ 分别表示真解Ex, Ey, Hz在tm时刻和离散点 $(i + \frac{1}{2}, j), (i, j + \frac{1}{2}), (i + \frac{1}{2}, j + \frac{1}{2})$ 上的近似值。

根据文献[6]中的方法, 对方程 (1) —方程 (3) , 定义如下的格式S-FDTDI:

I-Stage 1:

${\begin{cases} \frac{E_{y_{i, j + \frac{1}{2}}}^{n + 1} - E_{y_{i, j + \frac{1}{2}}}^{n}}{Δ t} = - \frac{1}{2 ε} δ_{x} (Η_{z_{i, j + \frac{1}{2}}}^{*} + Η_{z_{i, j + \frac{1}{2}}}^{n}) - \\ \frac{σ}{2 ε} (E_{y_{i, j + \frac{1}{2}}}^{n + 1} + E_{y_{i, j + \frac{1}{2}}}^{n}) \\ \frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{*} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}}{Δ t} = \frac{1}{2 μ} δ_{x} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}) \end{cases} (5)$

I-Stage 2:

${\begin{cases} \frac{E_{x_{i, j + \frac{1}{2}}}^{n + 1} - E_{x_{i + \frac{1}{2}, j}}^{n}}{Δ t} = - \frac{1}{2 ε} δ_{y} (Η_{z_{i + \frac{1}{2}, j}}^{n + 1} + Η_{z_{i + \frac{1}{2}, j}}^{n}) - \\ \frac{σ}{2 ε} (E_{x_{i + \frac{1}{2}, j}}^{n + 1} + E_{x_{i + \frac{1}{2}, j}}^{n}) \\ \frac{Η_{z_{i + \frac{1}{2} j + \frac{1}{2}}}^{n + 1} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{*}}{Δ t} = \frac{1}{2 μ} δ_{y} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}) \end{cases} (6)$

如果去掉中间层上的值 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{*}$ , 可以发现式 (5) 和式 (6) 等价于带有一阶摄动项的Crank-Nicolson格式。因此, 为了提高精度, 在I-Stage 1中加入扰动项, 得到修正格式, 记为S-FDTDII:

II-Stage 1:

II-Stage 2:

从式 (5) —式 (8) 可以看到, S-FDTDI和S-FDTDII的每一个Stage包含两个方程, 要比每一个Stage包含三个方程的ADI-FDTD[4,5]方法简单, 更利于程序计算。下面以格式S-FDTDII为例说明计算步骤。

对式 (7) 中的第二式进行变形得到

$\begin{array}{l} Η_{z_{i + \frac{1}{2}, j \frac{1}{2}}}^{*} = Η_{z_{i} + \frac{1}{2}, j + \frac{1}{2}}^{n} - \frac{△ t}{2 μ △ x} (E_{y_{i + 1, j + \frac{1}{2}}}^{n^{+} 1} - \\ E_{y_{i, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + 1, j + \frac{1}{2}}}^{n} - E_{y_{i, j + \frac{1}{2}}}^{n}), (9) \end{array}$

然后将此表达式代入式 (7) 中的第一式, 变形可以得到:

$\begin{array}{l} - \frac{△ t^{2}}{4 μ ε △ x^{2}} E_{y_{i - 1, j + \frac{1}{2}}}^{n + 1} + (1 + \frac{σ △ t}{2 ε} + \frac{△ t^{2}}{2 μ ε △ x^{2}}) E_{y_{i, j + \frac{1}{2}}}^{n + 1} - \\ \frac{△ t^{2}}{4 μ ε △ x^{2}} E_{y_{i + 1, j + \frac{1}{2}}}^{n + 1} = \frac{△ t^{2}}{4 μ ε △ x^{2}} E_{y_{i - 1, j + \frac{1}{2}}}^{n} + \\ (1 - \frac{σ △ t}{2 ε} - \frac{△ t^{2}}{2 μ ε △ x^{2}}) E_{y_{i, j + \frac{1}{2}}}^{n} + \frac{△ t^{2}}{4 μ ε △ x^{2}} \times E_{y_{i + 1, j + \frac{1}{2}}}^{n} - \\ \frac{△ t^{2}}{2 μ ε △ x^{2}} (E_{x_{i + \frac{1}{2}, j + 1}}^{n} - E_{x_{i - \frac{1}{2}, j + 1}}^{n} - \\ E_{x_{i + \frac{1}{2}, j}}^{n} - E_{x_{i - \frac{1}{2}, j}}^{n}) - \frac{△ t}{ε △ x} (Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n} - Η_{z_{i - \frac{1}{2}, j + \frac{1}{2}}}^{n}) (10) \end{array}$

这是一个关于 $E_{y_{i, j + \frac{1}{2}}}, (i = 1, 2, \dots, Ι - 1)$ 的线性代数方程组, 未知量沿x方向排列, 可认为是沿x轴方向计算。其系数矩阵为三对角阵, 可用追赶法直接求出 ${E_{y_{i, j + \frac{1}{2}}}^{n + 1}}$ , 然后代入H*z的表达式式 (9) , 显式求出 ${Η_{z_{i + \frac{1}{1}, j + \frac{1}{2}}}^{*}}$ 。将H*z的值代入 II—Stage 2中, 采用与II—Stage 1相同的方法, 求出 ${E_{x_{i + \frac{1}{2}, j}}^{n + 1}}$ 和 ${Η_{z_{i + \frac{1}{1}, j + \frac{1}{2}}}^{n + 1}}$ 。

2 误差分析

式 (8) 中的第二个方程减去式 (7) 中的第二个方程得出 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{*}$ 的表达式, 代入式 (7) 和式 (8) 的第一个方程得到S-FDTDII的等价格式:

将式 (11) 中第二个方程最右端的一项改为 $\frac{△ t}{4 μ ε} δ_{x} δ_{y} (E_{x_{i, j + \frac{1}{2}}}^{n + 1} + E_{x_{i, j + \frac{1}{2}}}^{n})$ , 就得到S-FDTDI的等价格式。从上述方程可以推出S-FDTDII的截断误差。令 $ζ_{i + \frac{1}{2}, j}^{n + \frac{1}{2}}$ , $ξ_{i, j + \frac{1}{2}}^{n + \frac{1}{2}}$ , $η_{i + \frac{1}{2}, j + \frac{1}{2}}^{n + \frac{1}{2}}$ 分别表示上述等价格式中每个方程的截断误差, 根据泰勒公式可以导出它们的表达式为

$\begin{array}{l} ζ_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} = Δ t^{2} [\frac{1}{24} \frac{\partial^{3} E_{x}}{\partial t^{3}} (τ_{11}, x_{i + \frac{1}{2}}, y_{j}) + \frac{1}{8 ε} \frac{\partial^{3} Η_{z}}{\partial t^{2} \partial x} \times \\ (τ_{12}, x_{11}, y_{j}) - \frac{σ}{8 ε} \frac{\partial^{2} E_{x}}{\partial t^{2}} (τ_{13}, x_{i + \frac{1}{2}}, y_{j})] + \\ \frac{Δ y^{2}}{24 ε} \frac{\partial^{3} Η_{z}}{\partial y^{3}} (t^{n + \frac{1}{2}}, x_{i + \frac{1}{2}}, y_{11}); \end{array}$

$\begin{array}{l} ξ_{i, j + \frac{1}{2}}^{n + \frac{1}{2}} = Δ t^{2} [\frac{1}{24} \frac{\partial^{3} E_{y}}{\partial t^{3}} (τ_{21}, x_{i}, y_{j + \frac{1}{2}}) + \frac{1}{8 ε} \frac{\partial^{3} Η_{z}}{\partial t^{2} \partial x} \times \\ (τ_{22}, x_{12}, y_{j + \frac{1}{2}}) + \frac{σ}{8 ε} \frac{\partial^{2} E_{y}}{\partial t^{2}} (τ_{23}, x_{i}, y_{j + \frac{1}{2}})] - \\ \frac{Δ x^{2}}{24 ε} \frac{\partial^{3} Η_{z}}{\partial x^{3}} (t^{n + \frac{1}{2}}, x_{13}, y_{j + \frac{1}{2}}); \end{array}$

$\begin{array}{l} η_{i + \frac{1}{2}, j + \frac{1}{2}}^{n + \frac{1}{2}} = Δ t^{2} [- \frac{1}{24} \frac{\partial^{3} Η_{z}}{\partial t^{3}} (τ_{31}, x_{i + \frac{1}{2}}, y_{j + \frac{1}{2}}) + \\ \frac{1}{8 μ} \frac{\partial^{3} E_{x}}{\partial t^{2} \partial y} (τ_{32}, x_{i + \frac{1}{2}}, y_{12}) - \\ \frac{1}{8 μ} \frac{\partial^{3} E_{x}}{\partial t^{2} \partial x} (τ_{33}, x_{31}, y_{j + \frac{1}{2}})] + \\ \frac{1}{24 μ} [Δ y^{2} \frac{\partial^{3} E_{x}}{\partial y^{3}} (t^{n + \frac{1}{2}}, x_{i + \frac{1}{2}}, y_{32}) + \\ Δ x^{2} \frac{\partial^{3} E_{y}}{\partial x^{3}} (t^{n + \frac{1}{2}}, x_{32}, y_{j + \frac{1}{2}})] (12) \end{array}$

式 (12) 中tn≤τ1k, τ2l, τ3m≤tn+1; $x_{i - \frac{1}{2}} \leq x_{2 k}, x_{3 l} \leq x_{i + \frac{1}{2}}$ ; $y_{j - \frac{1}{2}} \leq y_{1 k}, y_{2 l}, y_{3 m} \leq y_{i + \frac{1}{2}}$ ;k, l, m=1, 2, 3。因此只要真解Ex, Ey, Hz足够光滑, 则有

$\begin{array}{l} | ζ_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} | + | ξ_{i, j + \frac{1}{2}}^{n + \frac{1}{2}} | + | η_{i + \frac{1}{2}, j + \frac{1}{2}}^{n + \frac{1}{2}} | \leq \\ C_{μ ε} Μ (Δ x^{2} + Δ y^{2} + Δ t^{2}) 。 \end{array}$

其中: $C_{μ ε} = \frac{1}{24} + \frac{1}{8 ε} + \frac{σ}{8 ε} + \frac{1}{8 μ} + \frac{1}{24 ε} + \frac{1}{24 μ}$ , M为常数, 与Ex, Ey, Hz的导数有关。类似地, 可推出格式S-FDTDI的截断误差界为CμεM′ (Δx2+Δy2+Δt) 。因此S-FDTDI关于Δt是一阶的, 修正格式S-FDTDII关于Δt是二阶的.所以修正格式S-FDTDII的精度比S-FDTDI高一阶。

3 数值验证

为了验证所提出的格式的有效性和与理论分析的一致性, 现用这两种格式来求解一个具体问题方程 (1) — (3) , 其中初始条件为Ex0=cosπxsinπy, Ey0=-sinπxcosπy, Hz0=-2cosπxcosπy, ε=1, μ=1, σ=3π可以验证这种具体问题的真解为

Ex=e-πtcosπxsin πy,

Ey=-e-πtsinπxcosπy,

Hz=-2e-πtcosπxcosπy。

设 $E_{x} (t^{n}, •) = E_{x} (t^{n}, x_{i + \frac{1}{2}}, y_{j}), E_{y} (t^{n}, •) = E_{y} (t^{n}, x_{i}, y_{j + \frac{1}{2}}), Η_{z} (t^{n}, •) = Η_{z} (t^{n}, x_{i + \frac{1}{2}}, y_{j + \frac{1}{2}})$ 表示真解, 定义差分解 ${E_{x_{i + \frac{1}{2}, j}}^{n}}, {E_{y_{i, j + \frac{1}{2}}}^{n}} {Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}}$ 误差:

在下列表格中, 令Err-E=‖E (tn) -En‖, Err-H=‖Hz (tn) -H $_{z}^{n}$ ‖表示差分解的绝对误差, Ra-E、Ra-H分别表示:En、H $_{z}^{n}$ 的收敛阶, CPU表示运算时间 (以秒为单位) , h=Δx=Δy表示空间步长, Δt表示时间步长。

表1—表3给出了一些计算结果:

表1给出了T=1, Δt=0.001时, S-FDTDI和S-FDTDII的绝对误差、收敛阶和运算时间, 由此看出两种格式在空间上收敛阶都能达到二阶。当Δt比较小时, 格式S-FDTDII的结果更好一些。

表2给出了T=1, h=0.001时, 格式S-FDTDI和S-FDTDII在不同时间步长下的绝对误差和收敛阶, 由此可以看出S-FDTDI在时间上的收敛阶为一阶, 而S-FDTDII在时间上的收敛阶为二阶。而且格式S-FDTDII的误差比S-FDTDI的误差小的多。

从表3可以看出, 在处理σ=0的简谐问题时, 两种格式的误差相同, 运算所耗费的时间, S-FDTDII比ADI-FDTD要少很多;而且当处理非0电导率的波导问题时, S-FDTDII方法所产生的绝对误差比ADI-FDTD小, 运算时间快了38.7%, 因此S-FDTDII比ADI-FDTD好。

参考文献

[1] Yee K S.Numerical solution of initial boundary value problems invol-ving Maxwell’s equa-tions in isotropic media.IEEE Trans Antennasand Propagation, 1996;14:302—307

[2] Courant R, Friedrichs K, Lewy H.On the partial difference equationsof mathematical physics.IBMJournal, 1967;11:215—234

[3] Namiki T.A new FDTD algorithm based on alternating direction im-plicit method.IEEE Trans Microwave Theory Tech, 1999;47:2003—2007

[4] Zheng Fenghua, Chen Zhizhang.Numerical dispersion analysis of theunconditionally stable 3-D ADI-FDTD Method.IEEE Trans on Micro-wave Theory Tech, 2001;49 (5) :1006—1009

[5] Taflove A, Hagness S.Computational electrodynamics:the finite-differ-ence time-domain method, (second ed) .Boston, MA:ArtechHouse, 2000

【时域有限差分法】推荐阅读：

时域有限差分算法09-17

时域有限差分方法09-23

时域有限元方法07-02

有限差分08-15

时域方法06-11

时域积分方程方法06-19

稳定性、时域分析07-08

控制系统的时域分析09-26

有限切开05-12

有限发展05-25