时域有限差分法对平面TE波的MATLAB仿真

时域有限差分法对平面TE波的

MATLAB仿真

摘 要

时域有限差分法是由有限差分法发展出来的数值计算方法。自1966年Yee

在其论文中首次提出时域有限差分以来，时域有限差分法在电磁研究领域得到了广泛的应用。主要有分析辐射条线、微波器件和导行波结构的研究、散射和雷达截面计算、分析周期结构、电子封装和电磁兼容的分析、核电磁脉冲的传播和散射以及在地面的反射及对电缆传输线的干扰、微光学元器件中光的传播和衍射特性等等。

由于电磁场是以场的形态存在的物质，具有独特的研究方法，采取重叠的研究方法是其重要的特点，即只有理论分析、测量、计算机模拟的结果相互佐证，才可以认为是获得了正确可信的结论。时域有限差分法就是实现直接对电磁工程问题进行计算机模拟的基本方法。在近年的研究电磁问题中，许多学者对时域脉冲源的传播和响应进行了大量的研究，主要是描述物体在瞬态电磁源作用下的理论。另外，对于物体的电特性，理论上具有几乎所有的频率成分，但实际上，只有有限的频带内的频率成分在区主要作用。

文中主要谈到了关于高斯制下完全匹配层的差分公式的问题，通过MATLAB程序对TE波进行了仿真，模拟了高斯制下完全匹配层中磁场分量瞬态分布。得到了相应的磁场幅值效果图。

关键词：时域有限差分 完全匹配层 MATLAB 磁场幅值效果图

目 录

摘 要 .............................................................................................................................................. 1 目 录 .............................................................................................................................................. 3 第一章 绪 论 .............................................................................................................................. 4 1.1 课题背景与意义 ................................................................................................................. 4 1.2 时域有限差分法的发展与应用 ......................................................................................... 4 2.1 Maxwell方程和Yee氏算法 .............................................................................................. 7 2.2 FDTD的基本差分方程 ........................................................................................................ 9 2.3 时域有限差分法相关技术 ............................................................................................... 11 2.3.1 数值稳定性问题 ........................................................................................................ 11 2.3.2 数值色散 .................................................................................................................... 12 2.3.3 离散网格的确定 ........................................................................................................ 13 2.4 吸收边界条件 ................................................................................................................... 13 2.4.1 一阶和二阶近似吸收边界条件 ................................................................................ 14 2.4.2 二维棱边及角顶点的处理 ........................................................................................ 17 2.4.3 完全匹配层 ................................................................................................................ 19 2.5 FDTD计算所需时间步的估计 .......................................................................................... 23 第三章 MATLAB的仿真的程序及模拟 .................................................................................. 25 3.1 MATLAB程序及相应说明 .................................................................................................. 25 3.2 出图及结果 ....................................................................................................................... 28 3.2.1程序部分 ..................................................................................................................... 28 3.2.2 所出的效果图 ............................................................................................................ 29 第四章 结 论 ............................................................................................................................ 31 参考文献 ........................................................................................................................................ 32

第一章 绪 论

1.1 课题背景与意义

20世纪60年代以来，随着计算机技术的发展，一些电磁场的数值计算方法逐步发展起来，并得到广泛应用，其中主要有：属于频域技术的有限元法（FEM）、矩量法(MM)和单矩法等；属于时域技术方面的时域有限差分法(FDTD)、传输线矩阵法(TLM)和时域积分方程法等。此外，还有属于高频技术的几何衍射理论(GTD)和衍射物理理论(PLD)等。各种方法都具有自己的特点和局限性，在实际中经常把它们相互配合而形成各种混合方法[1~2]。其中FDTD是一种已经获得广泛应用并且有很大发展前景的时域数值计算方法。时域有限差分（FDTD）方法于1966年由K.S.Yee[3] 提出并迅速发展，且获得广泛应用。K.S.Yee用后来被称作Yee氏网格的空间离散方式，把含时间变量的Maxwell旋度方程转化为差分方程，并成功地模拟了电磁脉冲与理想导体作用的时域响应。但是由于当时理论的不成熟和计算机软硬件条件的限制，该方法并未得到相应的发展。20世纪80年代中期以后，随着上述两个条件限制的逐步解除，FDTD便凭借其特有的优势得以迅速发展。它能方便、精确地预测实际工程中的大量复杂电磁问题，应用范围几乎涉及所有电磁领域，成为电磁工程界和理论界研究的一个热点。目前，FDTD日趋成熟，并成为分析大部分实际电磁问题的首选方法。

另外，利用矩量法求解电磁场问题时，要用到并失Green函数。对于某些问题，可以找到其解析形式的并失Green函数；而对于复杂的问题，很难找到其解析形式的并失Green函数，这样就使得问题无法解决。作为时域分析中的一个重要数值方法，FDTD不存在这样的问题。

1.2 时域有限差分法的发展与应用

经过四十多年的发展，FDTD已发展成为一种成熟的数值计算方法。在发展过程中，几乎都是围绕几个重要问题展开的，即数值稳定性、计算精度、数值色散、激励源技术以及开域电磁问题的吸收边界条件等。

数值稳定和计算精度对任何一种数值计算方法都是至关重要的。A.Taylor和M.E.Brodwin[4]利用本征值方法给出了直角坐标系下FDTD的空间步长与时间步长之间的关系。X.Min等[5]研究了存在边界条件时FDTD的稳定性问题。对于数值色散，与实际的物理色散不同，它是由电磁场量在空间和时间上的对波动方程作差分近似处理造成的。这种色散引起的误差造成在计算区域内传播的电磁波逐渐畸变[6~7]。K. L. Shlager 等[8]比较了二维和三维空间中几种正交网格算法的色散误差。当采用其他变形或非正交网格时，必须重新分析其数值稳定性和色散特性[9~11]，P.Monk 和 E.Suli[12]分析了不均匀长方体网格算法的稳定性。

激励源的设计和引入也是FDTD的一个重要任务。目前，应用最广泛的激励源引入技术是总场/散射场体系[12]。对于散射问题，通常在FDTD计算空间中引入连接边界，它将整个计算空间划分为内部的总场区和外部的散射场区，如图1-1。利用Huygens原理，可以在连接边界处引入入射场，使入射场的加入变得简单易行。

II 散射场区I 总场区散射体吸收边界连接边界入射波

图1-1

开域电磁问题中，为了在有限的计算空间内模拟无限空间中的电磁问题，必须在计算空间的截断边界处设置吸收边界条件。吸收边界条件从开始简单的插值边界，已经发展了多种吸收边界条件。在早期得到广泛应用的是G.Mur[13]的一阶和二阶吸收边界条件，它是基于B.Engquist和A.Majda[14]的单向波方程而提出的差分格式，在FDTD仿真区域外边界具有0.5%到5%的反射系数。目前应用最广泛的是J.P.Berenger[15-17]的分裂式完全匹配层，以及Z.S.Sacks等[18]和S.D.Gedney[20]的各向异性介质的完全匹配层，它们可使FDTD模拟的最大动态范围达到80dB。

另一方面，为了更好的拟合研究对象的形状，克服台阶逼近带来的误差，D.E.Merewether[19]提出了柱坐标系下的网格剖分方法，R.Holland[20]提出了球坐标系下的网格剖分方法，P.Monk和E.Suli[12]提出了变网格步长方法，S.S.Zivanovic等[21]和P.Thoma等[22]提出了亚网格技术（即在一般区域采用粗网格，在电磁场快变区域采用精细网格）。利用这些技术，可以更精确地模拟各种复杂的结构，适应各种复杂的介质，提高了复杂介质中数值计算的精度。

时域模拟一般获得的是近场电磁信息，为了得到诸如天线方向图或散射体雷达散射截面之类的远场信息，必须获得计算区域以外的频域场或瞬态场。多位学者在这方面做了许多工作，发展了一种高效的时域近远场变换方法[23-26]。借助这种方法，可以实现由计算区域内近场数据到计算区域外远场数据的外推。目前，粗糙面散射的FDTD，传递函数在FDTD中的应用，周期介质、各向异性介质、色散介质和含有集中元件的FDTD，以及网络并行FDTD技术等方面也取得了很大进展。

FDTD在迅速发展的同时，也获得了非常广泛的应用。目前，它几乎被应用到了电磁场工程中的各个方面，例如：电磁散射、生物电磁计量学、辐射天线的分析、微波器件和导行波结构的研究、散射和雷达截面的计算、周期结构的分析、电子封装和电磁兼容的分析、核电磁脉冲传播和散射的分析、以及微光学元器件中光的传播和衍射特性的分析等。随着新技术的不断提出，其应用范围和成效正在迅速地扩大和提高。

第二章 时域有限差分法的基本原理

Maxwell方程是描述宏观电磁现象的一组基本方程。这组方程即可以写成微分形式，又可以写成积分形式。FDTD方法由Maxwell旋度方程的微分形式出发，利用二阶精度的中心差分近似，直接将微分运算转换为差分运算，这样达到了在一定体积内和一段时间上对连续电磁场数据的抽样压缩。

2.1 Maxwell方程和Yee氏算法

根据

[27]

中电磁场基本方程组的微分形式，若在无源空间，其空间中的媒质是

各向同性、线性和均匀的，即媒质的参数不随时间变化且各向同性，则Maxwell旋度方程可写成：

??H???E??E （2-1a） ?t?H??mH （2-1b） ?t??E???式中，E是电场强度，单位为伏/米（V/m）；H是磁场强度，单位为安/米（A/m）；?表示介质介电系数，单位为法拉/米（F/m）； ?表示磁导系数，单位为亨利/米（H/m）；?表示介质电导率，单位为西门子/米（S/m）；?m表示导磁率，单位为欧姆/米（?/m）。

在直角坐标系中，（2-1）式可化为如下六个标量方程：

??Ex?Hz?Hy?????Ex??y?z?t???Ey?Hx?Hz??????Ey? （2-2） ?z?x?t??Hy?Hx??Ez?????Ez??x?y?t????Hx?Ez?Ey??????mHx??y?z?t???Hy?Ex?Ez???????mHy? （2-3） ?z?x?t??Ey?Ex??Hz??????mHz??x?y?t??这六个偏微分方程是FDTD算法的基础。

K.S.Yee[3]在1966年建立了如图2-1所示的空间网格，这就是著名的Yee氏元胞网格。

zExHzExEyEzHxEyoExyEzHyEzEyx

图2-1 Yee氏网格及其电磁场分量分布

并引入如下的差分近似方法对（2-2）、(2-3)式中的六个偏微分方程进行了差分离散。令F(x,y,z,t)代表E或H在直角坐标系中某一分量，在时间和空间域中的离散可记为

F(x,y,z,t)?F(i?x,j?y,k?z,n?t)?Fn(i,j,k) (2-4)

式中，?x、?y和?z分别是长方体网格沿x、y、z方向的空间步长，?t是时间步长，i、j、

k分别是沿x、y、z方向的网格编号，n是时间步数。对F(x,y,z,t)关于时间和空间的一阶偏导数取中心差分近似，具有二阶精度，即

?F(x,y,z,t)?xFn(i?12,j,k)?Fn(i?12,j,k)2????O?x (2-5a) x?i?x?x???F(x,y,z,t)?y?F(x,y,z,t)?z?F(x,y,z,t)?tFn(i,j?12,k)?Fn(i,j?12,k)2????O?y (2-5b) y?j?y?y??Fn(i,j,k?12)?Fn(i,j,k?12)2????O?z (2-5c) z?k?z?z??Fn?12(i,j,k)?Fn?12(i,j,k)2????O?t (2-5d) t?n?t?t?? 在FDTD中，空间上连续分布的电磁场物理量离散的空间排布如图2-1所示。由图可见，电场和磁场分量在空间交叉放置，使得在每个坐标平面上每个电场分量被磁场环绕，每个磁场分量也被电场环绕。这种电磁场的空间结构与电磁感应和电磁波传播的规律相符，在每一个网格单元都能满足法拉第感应定律和安培环流定律。各分量的空间相对位置也适合于Maxwell方程的差分计算，能够恰当地描述电磁场的传播特性。同时，电场和磁场在时间上交替抽样，抽样时间间隔相差半个时间步，使Maxwell旋度方程离散以后构成显式差分方程，从而可以在时间上迭代求解，而不需要进行矩阵求逆运算。因此，由给定相应电磁问题的初始条件，FDTD就可以逐步推进地求得以后各个时刻空间电磁场的分布。

2.2 FDTD的基本差分方程

根据上述原则，可将（2-2）、（2-3）式离散为如下的差分方程形式：

Exn?1(i?1,j,k)?CA(i?1,j,k)?Exn(i?1,j,k)?CB(i?1,j,k)

2222（2-6a）

n?1n?12211?H(i?,j?,k)?H(i?1,j?1,k)?y??? zz2222 ???n?1n?1221111?H(i?,j,k?)?H(i?,j,k?)?z??yy??2222????Eyn?1(i,j?1,k)?CA(i,j?1,k)?Eyn(i,j?1,k)?CB(i,j?1,k) 2222(2-6b)

n?1n?12211?H(i,j?,k?)?H(i,j?1,k?1)?z??? xx2222 ???n?1n?1221111?H(i?,j?,k)?H(i?,j?,k)?x??zz??2222????Ezn?1(i,j,k?1)?CA(i,j,k?1)?Ezn(i,j,k?1)?CB(i,j,k?1) 2222(2-6c)

n?1n?12211?H(i,?j,k?)?H(i?1,j,k?1)?x??? yy ??222211?n?n?221111?H(i,j?,k?)?H(i,j?,k?)?y??xx??2222????Hx2(i,j?1,k?1)?CP(i,j?1,k?1)?Hx2(i,j?1,k?1)?CQ(i,j?1,k?1) (2-6d)

22222222n?1n?1?Ezn(i,j?1,k?1)?Ezn(i,j,k?1)?y??? 22 ???nn11?E(i,j?,k?1)?E(i,j?,k)?z??yy??22????Hn?12yn?11111(i?,j,k?)?CP(i?,j,k?)?Hy2(i?1,j,k?1)?CQ(i?1,j,k?1) (2-6e) 22222222nn??Ex(i?12,j,k?1)?Ex(i?12,j,k)?z?? ??? nn11?Ez(i?1,j,k?)?Ez(i,j,k?)?x????22????Hzn?12n?11111(i?,j?,k)?CP(i?,j?,k)?Hz2(i?1,j?1,k)?CQ(i?1,j?1,k) 22222222nn??Ey(i?1,j?12,k)?Ey(i,j?12,k)?x?? ??? nn11?Ez(i?,j?1,k)?Ez(i?,j,k)?y????22?(2-6f)

???式中

CA(i,j,k)?2?t2?(i,j,k)??(i,j,k)?t， (2-7a) CB(i,j,k)?2?(i,j,k)??(i,j,k)?t2?(i,j,k)??(i,j,k)?t2?t2?(i,j,k)??m(i,j,k)?t， (2-7b) CQ(i,j,k)?CP(i,j,k)?2?(i,j,k)??m(i,j,k)?t2?(i,j,k)??m(i,j,k)?t(2-6)式就是FDTD的基本差分方程组。从式中可以看出，方程组中含有半个空间步和半个时间步，为了便于编程，可将（2-6）式改写成如下形式[28]：

Ex(i,j,k)?CA(i,j,k)Ex(i,j,k)?CB(i,j,k)?H(i,j,k)?HZ(i,j?1,k)Hy(i,j,k)?Hy(i,j,k?1)?(2-8a)

??Z???y?z??Ey(i,j,k)?CA(i,j,k)Ex(i,j,k)?CB(i,j,k)?H(i,j,k)?Hx(i,j,k?1)HZ(i,j,k)?HZ(i?1,j,k)?(2-8b)

??x???z?x??EZ(i,j,k)?CA(i,j,k)EZ(i,j,k)?CB(i,j,k)?Hy(i,j,k)?Hy(i?1,j,k)Hx(i,j,k)?Hx(i,j?1,k)?(2-8c)

?????x?y??

Hx(i,j,k)?CP(i,j,k)Hx(i,j,k)?CQ(i,j,k)?EZ(i,j?1,k)?EZ(i,j,k)Ey(i,j,k?1)?Ey(i,j,k)?(2-8d)

?????y?z??Hy(i,j,k)?CP(i,j,k)Hx(i,j,k)?CQ(i,j,k)?E(i,j,k?1)?Ex(i,j,k)EZ(i?1,j,k)?EZ(i,j,k)?(2-8e)

??x???z?x??HZ(i,j,k)?CP(i,j,k)HZ(i,j,k)?CQ(i,j,k)?Ey(i?1,j,k)?Ey(i,j,k)Ex(i,j?1,k)?Ex(i,j,k)?(2-8f)

?????x?y??根据上述FDTD差分方程组可得出计算电磁场的时域推进计算方法，如图2-2所示。

图2-2 FDTD在时域的交叉半步逐步推进计算

已知t1=t0=n?t 0时刻空间各处的电磁场初始值 计算t2=t1+?t2 时刻空间各处的磁场值 循环n次 计算t1=t2+?t2 时刻空间各处的电场值 式(2-8a)~(2-8c)的等号左边的电场值是第n次循环的电场值，等号右边的电场值是第n-1次循环存储在内存中的电场值，磁场值是本次循环计算得到的磁场值；式(2-8d)~(2-8f) 等号左边的磁场值是第n次循环的磁场值，等号右边的磁场值和电场值都是第n-1次循环存储在内存中的场值。这样，就解决了半个时间步在程序中无法表示的问题，而且也没有破坏电磁场在时间上逐步推进的逻辑关系。

2.3 时域有限差分法相关技术 2.3.1 数值稳定性问题

上述FDTD方程是一种显式差分方程，在执行时，存在一个重要的问题：即算法的稳定性问题。这种不稳定性表现为在解显式方程时，随着时间步数的继续

增加，计算结果也将无限制地增加。Taflove等[4]于1975年对Yee氏差分格式的稳定性进行了讨论，并导出了对时间步长的限制条件。数值解是否稳定主要取决于时间步长?t与空间步长?x、?y、?z的关系。对于非均匀媒质构成的计算空间选用如下的稳定性条件：

?t?1 (2-9)

111v()2?()2?()2?x?y?z(2-9)式是空间和时间离散之间应当满足的关系，又称为Courant稳定性条件。若采用均匀立方体网格：?x??y??z??s

?s?t? (2-10)

v3其中，v为计算空间中的电磁波的最大速度。

2.3.2 数值色散

FDTD方程组是对Maxwell旋度方程进行差分近似，在进行数值计算时，将

会在计算网格中引起数字波模的色散，即在FDTD网格中，电磁波的相速与频率有关，电磁波的相速度随波长、传播方向及变量离散化的情况不同而改变。这种关系由非物理因素引起，且色散将导致非物理因素引起的脉冲波形畸变、人为的各向异性和虚假折射等现象。显然，色散与空间、时间的离散间隔有关，如下式所示：

k?y?111112???t?2?kx?x?2?y2?kz?z?2?kz?z? (2-11) ??sin?sin?sin?sin?sin????????2???z?2222?c?t?2?2???x?2?2???y?2????z??????式(2-11)是三维情况下在FDTD方法中的单色平面波数值色散关系的一般形式，它表明FDTD计算中波的传播速度与传播方向有关。式中kx、ky、kz分别是波矢量沿x、y、z方向的分量，?是角频率，c是被模拟的均匀介质中的光速。与数值色散关系相对应，在无耗介质中的单色平面波，色散解析关系是：

??c?2?kx2?ky2?kz2 (2-12)

由式(2-11)可知，当式(2-11)中的?t、?x、?y、?z均趋于零时，它就趋于式(2-12)。也就是说数值色散是由于用近似差分替代连续微分而引起的，而且在理论上可以减小到任意程度，只要此时时间步长和空间步长都足够小，但这将大大增加所需的计算机存储空间和计算时间，并使累积误差增加。因此，在实际计算中要根据问题的性质和计算机的软硬件条件来选择合适的时间步长和空间步长。为获得理想的色散关系，问题空间分割应按照小于正常网格的原则进行。一般选取的最大空间步长为?max??min20，?min为所研究范围内电磁波的最小波长。由上分析说明，数值色散在用FDTD法分析电磁场传播中的影响是不可能避免的，但我们可以尽可能的减小数值色散的影响。

2.3.3 离散网格的确定

无论是简单目标还是复杂目标，在进行FDTD离散时网格尺寸的确定，除了

受计算资源的限制不可能取得很小外，还需要考虑以下几个因素：

1．目标离散精确度的要求。网格应当足够小以便能精确模拟目标几何形状和电磁参数。

2．FDTD方法本身的要求。主要是考虑色散误差的影响。设网格为立方体

?x??y??z??，所关心频段的频率上限为fmax，对应波长为?min，则考虑

FDTD的数值色散要求

???minN （2-13）

通常N?10。上式是根据已知所关心频率上限情况下来确定FDTD网格尺寸?的；反之，若给定?，则FDTD计算结果可用的上限频率也随之确定。

3．入射波的要求。入射波的上限截止频率fc应包含所关心频率范围，即

fc?fmax。

2.4 吸收边界条件

由时域有限差分法的基本原理可知，在利用时域有限差分法研究电磁场时，

需在全部问题空间建立Yee氏网格空间，并存储每个单元网格上任一时间步的六个场分量用于下一时间步的计算。而在对于辐射、散射这类开放系统的实际研究

中，不可能有无限大的存储空间。因此，必须在某处将网格空间截断，且在截断边界网格点处运用特殊的场分量计算方法，使得向边界面行进的波在边界处保持“外向行进”特性、无明显的反射现象，并且不会使内部空间的场产生畸变，从而用有限网格空间模拟电磁波在无界空间中传播的情况。具有这种功能的边界条件称之为吸收边界条件，或辐射边界条件，或网格截断条件[29~31]，如图2-3所示。

截断边界散射体入射波散射波

图 2-3 附加截断边界使计算区域变为有限域

从FDTD的基本差分方程组可以看出，在截断边界面上切向场分量的计算需要利用计算空间以外的电磁场分量，因此FDTD基本差分方程对这些截断边界面上的场分量失效。如何处理截断边界上的场分量，使之与需要考虑的无限空间有尽量小的差异，是FDTD中必须很好解决的一个重要问题。实际上，这是要求在误差可容忍的范围内，计算空间中的外向波能够顺利通过截断边界面而不引起波的明显反射，使有限计算空间的数值模拟与实际情况趋于一致，对外向波而言，就像在无限大空间中传播一样。所以，需要一种截断边界网格处的特殊计算方法，它不仅要保证边界场计算的必要精度，而且还要大大消除非物理因素引起的波反射，使得用有限的网格空间就能模拟电磁波在无限空间中的传播。但是如果处理不当，截断边界面可能造成较大反射，构成数值模拟误差的一部分，甚至可能造成算法不稳定。

加于截断边界场分量符合上述要求的算法就称为吸收边界条件（Absorbing Boundary Conditions）。

2.4.1 一阶和二阶近似吸收边界条件

在截断边界附近通常没有激励源。考虑齐次波动方程

??2?2?21?2????x2??y2??z2?c2?t2??f?0 (2-14) ??式中，f表示直角坐标系下任意电磁场分量。

B.Engquist和A.Majda[15]利用偏微分算子对式(2-13)作因式分解，并分别取其Taylor级数展开式中的第一项和前两项近似，导出了适合直角坐标系下FDTD吸收边界条件的单向波动方程，这就是Engquist-Majda吸收边界条件。设三维长方体FDTD区域0表2-1 三维长方体FDTD区域的一阶和二阶吸收边界条件

一阶近似 二阶近似 x?0 ??1?????f??xc?t?x?0?0 ?0 ?1?21?21??2?2???22???2???f2??2??y?z???c?x?tc?t?1?21?21??2?2???22???2???f2??2??y?z???c?x?tc?t?1?21?21??2?2???c?y?t?c2?t2?2??z2??x2??f????x?0?0 ?0 ?0 x?a ??1?????f??xc?t?x?ax?ay?0 ??1?????y?c?t??f??y?b ??1?????y?c?t??f??z?0 ??1?????f??zc?t?z?d ??1?????f??xc?t?

y?0?0 ?0 ?0 ?0 y?0y?b?1?21?21??2?2???c?y?t?c2?t2?2??z2??x2??f?????1?21?21??2?2???22???2???f2??2??y?x???c?z?tc?ty?b?0 ?0 z?0z?0z?d?1?21?21??2?2???22???2???f2??2??y?x???c?z?tc?tz?d?0

图2-4 近似吸收边界条件作用后残留的反射波与入射波之比

根据Yee氏元胞网格的特点，在FDTD截断边界面上只有电场切向分量和磁场法向分量。以x?0界面为例，此界面仅有Hx，Ey，Ez节点。由于FDTD中

Hx的计算式不涉及x?0区域，即不涉及截断边界界面外的节点。所以，吸收

边界条件将不考虑Hx，而只考虑电场切向分量Ey和Ez。以Ez为例，对一阶和二阶吸收边界条件分别有

??1?????EZ??xc?t?x?0?0 (2-15a)

?1?21?21??2?2???22???2???EZ2??2??y?z???c?x?tc?tx?0?0 (2-15b)

将式(2-15b)分别在距离边界半个空间步长的辅助网格点?i?12,j,k?12?处及

n?t时刻离散，此时Ez的位置在?i,j,k?12?，并对各项做差分近似，同时对

辅助网格点处的场值应用线性插值，假设?x??y??z??，可得到一阶吸收边界条件在三维情况的形式：

fn?1(i,j,k?12)?fn(i?1,j,k?12)?c?t??n?1?f(i?1,j,k?12)?fn(i,j,k?12)? (2-16) c?t??式中f代表截断边界上切向场分量。

将式(2-16)所示二阶吸收边界推广到长方体元胞各边?x、?y和?z不相等情形。设边界为x?0，对于Ez分量有

n?1n?1EZ(0,j,k?12)??EZ(1,j,k?12)?c?t??xn?1?EZ(1,j,k?12)?EZn?1(0,j,k?12)?

c?t??x2?x?EZn(0,j,k?12)?EZn(1,j,k?12)??c?t??xnn2E(0,j?1,k?12)?2E??ZZ(0,j,k?12)?x?c?t?nn??EZ(0,j?1,k?12)?EZ(1,j?1,k?12)??22??y??c?t??x???2EZn(1,j,k?12)?EZn(1,j?1,k?12)?

??nn2?EZ(0,j,k?32)?2EZ(0,j,k?12)??x?c?t?nn???E(0,j,k?12)?E(1,j,k?32)?ZZ22??z??c?t??x???2EZn(1,j,k?12)?EZn(1,j,k?12)???(2-17)

同理，可以得到所有截断边界面上切向场分量的一阶和二阶差分式。需要注意的是，用二阶近似条件计算界面上与棱边相邻的一列节点时会涉及棱边上的场值。因此，要想避免用到棱边上的场值，只需对截断边界面上切向场分量的计算按以下两种情况区别对待：（1）截断边界面上与棱边相邻的一列场分量采用一阶差分式；（2）截断边界面上其它场分量采用二阶差分式。这样，就完全不必考虑棱边上的场分量，避免了计算棱边上场分量所带来的误差。实际计算表明，这样做提高了吸收边界条件的精度和FDTD计算的稳定性。但是，就目前FDTD的发展来看，G.Mur的一阶和二阶吸收边界条件已不能满足目前高精度计算的要求。

2.4.2 二维棱边及角顶点的处理

对于二维电磁场问题，在用FDTD计算边界处的元胞时，将涉及到截断边界外侧的节点。如对于二维TM情况，电磁场分量有Ez,Hx,Hy，由图2-5可见，在用FDTD计算边界处的TM元胞

Hxn?12(i,j?12)和Hyn?12(i?12,j)时并不

n?1涉及截断边界以外E或H的节点。只有Hz(i,j)涉及截断边界外侧的H节点。

因此，只需给出边界处切向场分量Ez的吸收边界条件。同样，对于TE波只需给出边界处切向场分量Hz的吸收边界条件。表2-2为矩形截断边界面上节点的二阶Mur吸收边界条件。

EZ?i,j?12?HxHx?i?1,j?12??HyyZEZ?i,j?(i?12,j)?i?1,j?

x图2-5 二维TM左截断边界元胞

表2-2 矩形截断边界四边上的二阶Mur吸收边界条件

截断边界位置 二阶Mur吸收边界条件 TM TE x?0 ??EZ1?EZc??Hx???x?c?t?2?y??0 ??x?0??EZ1?EZc??Hx???x?c?t?2?y??0 ??x?a??HZ1?HZc??Ex???x?c?t?2?y??0 ??x?0??HZ1?HZc??Ex???x?c?t?2?y??0 ??x?ax?a y?0 y?b ??EZ1?EZc??Hy???y?c?t?2?x??0 ??y?0??EZ1?EZc??Hy???y?c?t?2?x??0 ??y?b??HZ1?HZc??Ey???y?c?t?2?x??0 ??y?0??HZ1?HZc??Ey???y?c?t?2?x??0 ??y?b在二维矩形计算区域的角点，吸收边界条件的离散式需特殊考虑。假设角点附近只有向外传播的行波，且传播方向沿角点处元胞的对角线，如图2-6所示的矩形区域左下角点处的TM元胞为例，导出适用于角点的吸收边界条件离散式。

根据Courant稳定条件c?t??2式，有2c?t?2?。在点?i0,j0?与对角点

?i0?1,j0?1?之间取一点P，对于沿对角线的外行波，有以下等式：

n?i0,j0??EZn?2?P? （2-18） EZ由于传播距离2c?t很小，上式中略去了振幅的衰减。在利用线性插值

n?2n?2n?2n?2????EZi0,j0??EZP?EZi0,j0??EZi0?1,j0?1? （2-19） ?2c?t2?由式（2-18）解出EZnZn?2?P?后代入（2-17）式得

?2c?t?n?22c?tn?2??E?i0,j0???1?EZ?i0,j0??EZ?i0?1,j0?1? （2-20） ?????或

n?1?i0,j0??EZn?i0?1,j0?1??EZc?t?2?n?1?EZ?i0?1,j0?1??EZn?i0,j0??c?t?2?（2-21）

对于矩形域的其它角点可作类似的处理。对于二维TE波情况，将式（2-19）（2-21），式中的EZ换为HZ即可。

?i,jb?0y右上角点?ia,jb?左上角点?i0?1,jb?1??ia?1,jb?1??i2?0?1,j0?1?对角点?ia?1,j0?1?EZ左下角点?i0,j0?x(ia,j0)右下角点

图2-6 矩形域四个角点

2.4.3 完全匹配层

完全匹配层[16~18]（Perfectly Matched Layer, PML）是1994年由J.P.Berenger首先提出，并将其设置在FDTD计算区域截断边界处，用来吸收外向电磁波。Berenger假设将电磁场分量在PML介质中分裂，并分别对各个分裂的场分量赋以不同的损耗。这就相当于在FDTD区域截断边界外设置了一种特殊的非物理的吸收介质层，该层介质的波阻抗与相邻介质的波阻抗完全匹配，因而外向波将无反射地穿过分界面进入PML。同时，由于PML为有耗介质，而且不依赖于外向波的入射角和波阻抗，即使为有限厚度，外向波在其中也会迅速衰减。在实际计算中，PML是目前一种很常用的吸收边界条件，有很好的吸收效果，其总的网格噪声能量是使用普通吸收边界条件时的1/107，可使FDTD模拟的最大动态范围达到80dB。

在PML介质中，其网格剖分方式与常规FDTD网格完全一致，每个场分量在网格中的位置也不变，只是都被分裂为两个子分量。这样，式(2-1)可写成

??Exy?t??yExy??(Hzx?Hzy)?y (2-22a)

?(Hyz?Hyx)?Exz???zExz?? (2-22b) ?t?z??Eyz?t?Eyx?t??zEyz??(Hxy?Hxz)?z (2-22c)

???xEyx???(Hzx?Hzy)?x (2-22d)

??(Hyz?Hyx)?Ezx??xEzx? (2-22e) ?t?x??Ezy?t?Hxy?t??yEzy???(Hxy?Hxz)?y (2-22f)

???myHxy???(Ezx?Ezy)?y (2-22g)

??(Eyz?Eyx)?Hxz??mzHxz? (2-22h) ?t?z?Hyz?t??mzHyz???(Exy?Exz)?z (2-22i)

?

??Hyx?t??mxHyx??(Ezx?Ezy)?x (2-22j)

??(Eyz?Eyx)?Hzx??mxHzx?? (2-17k) ?t?x??Hzy?t??myHzy??(Exy?Exz)?y (2-22l)

式中?x，?y，?z，?mx，?my，?mz为电导率和磁导率，描述了PML介质的各向异性。当?x??y??z??且?mx??my??mz??m时， PML介质退化为普通有耗介质;当??0且?m?0时，退化为自由空间中的Maxwell方程。同时，要满足PML介质的重要基本条件即阻抗匹配条件，如下式所示：

??m? (2-23) ?0?0式中?0和?0分别为真空的介电常数和真空的磁导率。

对式(2-18)做差分处理，可得到PML介质中的FDTD差分方程表达式：

Exy(i?12,j,k)?CAy(j)Exy(i?12,j,k)?CBy(j) ?HZn?12n?1n(i?12,j?12,k)?HZ?yn?12(i?12,j?12,k)(2-23a)

n?1nExz(i?12,j,k)?CAz(k)Eyz(i,j?12,k)?CBz(k)n?12Hx(i,j?12,k?12)?Hxn?12(i,j?12,k?12)(2-23b)

??zn?1nEyz(i,j?12,k)?CAz(k)Exy(i,j?12,k)?CBz(k)n?12Hxn?12(i,j?12,k?12)?HZ(i,j?12,k?12)(2-23c)

??zn?1nEyx(i,j?12,k)?CAx(i)Eyx(i,j?12,k)?CBx(i)n?12n?12Hz(i?12,j?12,k)?HZ(i?12,j?12,k)(2-23d)

??xn?1nEzx(i,j,k?12)?CAx(i)Ezx(i,j,k?12)?CBx(i)n?12Hy(i?12,j,k?12)?Hyn?12(i?12,j,k?12)(2-23e)

??xn?1nEzy(i,j,k?12)?CAy(j)Ezy(i,j,k?12)?CBx(j)n?12Hx(i,j?12,k?12)?Hxn?12(i,j?12,k?12) (2-23f)

??yn?12n?12Hxy(i,j?12,k?12)?DAy(j?12)Hxy(i,j?12,k?12)Ezn(i,j?1,k?12)?Ezn(i,j,k?12) (2-23g)

?DBy(j?12)??yn?12n?12Hxz(i,j?12,k?12)?DAz(k?12)Hxz(i,j?12,k?12)nEy(i,j?12,k?1)?Eyn(i,j?12,k) (2-23h)

?DBz(k?12)??zn?12n?12Hyz(i?12,j,k?12)?DAz(k?12)Hyz(i?12,j,k?12)Exn(i?12,j,k?1)?Exn(i?12,j,k) (2-23i)

?DBz(k?12)??zn?12n?12Hyx(i?12,j,k?12)?DAx(i?12)Hyx(i?12,j,k?12)nEz(i?1,j,k?12)?Ezn(i,j,k?12) (2-23j)

?DBx(i?12)??xHzxn?12(i?12,j?12,k)?DAx(i?12)Hzxnn?12(i?12,j?12,k)n ?DBx(i?12)?Ey(i?1,j?12,k)?Ey(i,j?12,k) (2-23k)

?xn?12n?12Hzy(i?12,j?12,k)?DAy(j?12)Hzy(i?12,j?12,k)nEx(i?12,j?1,k)?Exn(i?12,j,k) (2-23l)

?DBy(j?12)??y式中

CAq(?)?2???q(?)?t2???q(?)?tt，CBq(?)?2?t (2-24) 2???q(?)?ttDAq(??12)?2???m(??12)?t，2?t (2-25) DBq(??12)?2???m(??12)?t2???m(??12)?t其中，q?x,y,z；??i,j,k

在实际计算中，在FDTD计算中PML的设置不可能延伸到半无限空间，只能是有限厚度，因此PML的外侧边界需要特殊处理，通常情况下采用理想导体边界截断。这样，透入PML中的外向波到达理想导体边界处会反射回来，重新进入计算区域，PML的反射系数不再等于零。PML介质内沿q方向的电导率分布通常采取以下函数形式

?r? ?(r)???qqma?x???m （2-26）

式中，?为PML层的厚度，r为PML层靠近FDTD分界面的距离，m是电导率分布阶数，表示PML中电导率变化的程度，取整数。式（2-26）说明了实际计算中PML介质层厚度必须取若干个空间步长，使电导率从FDTD-PML分界面的0渐变到PML最外面的?qmax，避免电导率跃变太大，尽量消除数值反射。这样，如果相对于FDTD-PML分界面定义的外向波入射角为?，则PML内侧表面反射系数为

?2??qmaxcos?R????exp??m?1??0?? （2-27） ?0?使用PML吸收边界条件时，首先要选定三个参数：PML层数N，电导率分布阶数m，PML表面反射系数R(?)。大量的数值实验表明：（1）当层数N固定时，减小R(?)，即增加PML的衰减，可以使局部及总体误差都单调地减小。然而，当R(?)减小到一定程度后，这种现象不再出现，原因是存在由空间网格引起的固有误差。（2）增加PML层数N，可以使局部及总体误差都单调地减小，但是N过大又会使计算量剧增，需折衷考虑吸收效果和计算量。（3）电导率分布阶数m的改变不会影响计算量，却会影响PML的吸收效果，一般情况下，m越大，吸收效果越好。因此，在实际计算中参数的选择随所处理问题的不同而不同，需综合考虑，并由数值实验寻找最佳值。

2.5 FDTD计算所需时间步的估计

为了使计算达到稳定,通常计算所需要时间步按照电磁波往返穿越FDTD计算

区对角线3～5次来估计。若FDTD计算区总元胞数为N，则对角线上元胞约为

3N13 （三维)。按照Courant稳定条件，设计算中心区?t??(2c)，即穿

越对角线一次需要时间步为213总计算时间步约需23??6~10?N步。3N13。

12对于二维情况则约为22??6~10?N。或者说，计算时间步大约等于FDTD

计算区对角线上元胞数目的12～20倍。实际上，计算所需时间步还与散射体具体形状、结构有关。

图2-9给出了应用FDTD分析电磁场问题时的程序流程图

变量定义变量赋值 初始化波导建模加入激励源计算磁场HzStep=step+1边界和角点的处理计算电场Ex,EyStep=Max step?YES输出结果NO

图2-9 FDTD 程序流程图

第三章 MATLAB的仿真的程序及模拟

3.1 MATLAB程序及相应说明

% 二维FDTD TE波仿真 clear all; % 定义常数 pi=3.1415;

c=3.0e10; %高斯制下光速 f=1.0e15; %频率 lambda=c/f; %波长 nmax=400; %时间步数

del_s=lambda/20; %每最小波长20个采样点 del_t=0.5*del_s/c; %迭代时间步长 n=182; %真空区域网格数 np=9; %pml层数 N1=n+2*np; %总网格数 N=N1+1; %采样点数 M=4; %导电率渐变指数 sigma_max=(M+1)/1.50/pi/del_s; %最大导电率

% TE波的分量初始化 tic; figure(1);

axis([0 N 0 N -0.5 0.5]);

Ex=zeros(N1,N); %x方向为横向，采样点为网格的横向边，故行数+1 Ey=zeros(N,N1); %y方向为纵向，采样点为网格的纵向边，故列数+1 Bz=zeros(N1,N1); %矩阵行为纵向网格数，矩阵列为横向网格数,循环中用j表示行数，i表示列数 Bzx=zeros(N1,N1);

Bzy=zeros(N1,N1); Bzxx=zeros(nmax,2); %进入电磁场迭代计算 for tt=1:nmax for i=1:N1

if i>=np+1&&i<=N1-np di=0; elseif i<=np di=np-i+0.5; elseif i>=N1-np+1 di=np+i-N1-0.5;

end %di是采样点横向距PML内边界的距离 sigma_mx=sigma_max*(di/np)^M; for j=1:N1

if j>=np+1&&j<=N1-np dj=0; elseif j<=np dj=np-j+0.5; elseif j>=N1-np+1 dj=np+j-N1-0.5;

end %dj是采样点纵向距PML内边界的距离 sigma_my=sigma_max*(dj/np)^M;

if sigma_mx+sigma_my==0 %真空区 if j==100&&i==100

t=30; %可选择高斯脉冲 term=(tt-t);

% pulse=exp(-4*pi*term^2/20^2); pulse=sin(2*pi*tt/40); %可选正弦时谐源 c_miu=c*del_t/del_s;

Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j));

Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j));

Bz(i,j)=Bz(i,j)+Eterm1-Eterm2+pulse;%加入脉冲源 else

c_miu=c*del_t/del_s;

Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j)); Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j)); Bz(i,j)=Bz(i,j)+Eterm1-Eterm2; end

else %PML区

cpm=(1-2*c*sigma_mx*del_t)/(1+2*c*sigma_mx*del_t); cqm=c*del_t/(1+2*c*sigma_mx*del_t)/del_s; Bzx(i,j)=cpm*Bzx(i,j)-cqm*(Ey(i+1,j)-Ey(i,j));

cpm=(1-2*c*sigma_my*del_t)/(1+2*c*sigma_my*del_t); cqm=c*del_t/(1+2*c*sigma_my*del_t)/del_s; Bzy(i,j)=cpm*Bzy(i,j)+cqm*(Ex(i,j+1)-Ex(i,j)); Bz(i,j)=Bzx(i,j)+Bzy(i,j); end end end for i=2:N1

if i>=np+1&&i<=N1-np di=0; elseif i<=np di=np-i+1; elseif i>=N1-np+1 di=np+i-N1-1;

end %di是采样点横向距PML内边界的距离 sigma_ex=sigma_max*(di/np)^M; for j=1:N1

cam=(1-2*c*sigma_ex*del_t)/(1+2*c*sigma_ex*del_t);

cbm=c*del_t/(1+2*c*sigma_ex*del_t)/del_s; Ey(i,j)=cam*Ey(i,j)-cbm*(Bz(i,j)-Bz(i-1,j)); end end for i=1:N1 for j=2:N1

if j>=np+1&&j<=N1-np dj=0; elseif j<=np dj=np-j+1; elseif j>=N1-np+1 dj=np+j-N1-1;

end %dj是采样点纵向距PML内边界的距离 sigma_ey=sigma_max*(dj/np)^M;

cam=(1-2*c*sigma_ey*del_t)/(1+2*c*sigma_ey*del_t); cbm=c*del_t/(1+2*c*sigma_ey*del_t)/del_s; Ex(i,j)=cam*Ex(i,j)+cbm*(Bz(i,j)-Bz(i,j-1)); end end

Bzxx(tt,1)=Bz(90,50); %对靠近边界的中央磁场点采样 Bzxx(tt,2)=Bz(50,90);

3.2 出图及结果 3.2.1程序部分

figure(1); %可视化处理 clf;

mesh(Bz); %磁场的幅值 axis([0 N 0 N -0.5 0.5]); xlabel('i') ylabel('j')

drawnow; end figure(2); plot(Bzxx); figure(3);

title('磁场幅值分布图'); surface(Bz); shading interp; axis square; toc

3.2.2 所出的效果图

正弦时谐场源辐射效果图

高斯脉冲辐射效果图

第四章 结 论

通过对时域差分法的学习而进行的关于二维FDTD的TE波仿真，巩固了迭代式的形式与相应的循环编程技巧。但也存在一定问题，由于在吸收上我的编程思路不是很规范，所以在吸收效果上与理论上的完全匹配层吸收效果还有差距，但是已经基本达到了我的预期目的。可能是我的分区方法不对的缘故，原因就是出在sigma_max的设置上面，在高斯制下长度单位变为cm，所以同样的空间步长在高斯制下数量级就扩大了100倍，而sigma_max的计算式里面分母含有空间步长，所以如果不加以修正，sigma_max在高斯制下就变成了国际制下的1%，太大的导电率和太小的导电率都不能实现较好的吸收效果，下阶段应将其做适当修正并且写粒子仿真时用高斯制Maxwell方程可以参考角点反射大，多设置几层PML试一试。

参考文献

[1]J.M.Jin, M.Zunoubi, K.C.Donepudi, W.C.Chew, Frequency-Domain and Time-Domain Finite-Element Solution of Maxwell’s Equations Using Spectral Lanczos Decomposition Method. Comput. Methods Appl. Mech. Engrg. 1999,169:279~296

[2]K.Dongsoo,H.B.Lee. Hybird Full-Wave Analysis of Via-hole Grounds Using Finit-Difference and Finite-Element Time-Domain Methods. IEEE Transactions on Microwave Theory and Techniques. 1997,45(12):2217~2222 [3] Kane S. Yee Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media[J].IEEE Trans Antennas Propagation, 1966,AP-14(3):302~307 [4] A.Taylor and M.E.Brodwin. Numerical solution of steady-state EM scattering problems using the time-Maxwell’s equations. [J] IEEE Trans.Microwave Theroy Tech.,Aug. 1975，MTT-23:623~630

[5] X.Min,W.Sun and K.M.Chen. Stability analysis of the finite difference time domain method applied to unbounded electromagnetic problem [J]. IEEE Antennas and Propagation Society International Symposium,Dallas:May 1990,(4):1640~1643

[6]A.Taflove. Review of the formulation and applications of the finite-difference time-domain method for numerical modeling of electromagnetic wave interactions with arbitrary structures [J]. Wave Motion,June 1988,10(6): 547~582

[7] I.S.Kim and W.J.R.Hoefer. Numerical dispersion characteristics and stability factor for the TD-FD method [J].Electronics Letters, July 1990, 26(7): 485~487

[8] K.L.shlager,J.G.Msloney, S.L.Ray,and etc.. Relative accuracy of several finite-difference time-domain methods in two and three dimensions [J]. IEEE Transactions on Antennas and Propagation,Dec.1993,AP-41(12)：1732~1737

[9] S.L.Ray. Numerical dispersion and stability characteristics of finite-difference time-domain methods on nonorthogonal meshes [J]. IEEE Transactions on Antennas and Propagation,Feb.1993,AP-41(2)：233~235 [10] H.Shi and J.L.Drewriak. Dispersive comparison for DSI- and tenor-based nonorthogonal FDTD [J]. IEEE Microwave and Guided Wave Letters, May 1996,6(5):193~195

[11] F.Xiao and H.Yabe. Numerical dispersion relation for FDTD method in general curvilinear coordinates [J]. IEEE Microwave and Guided Wave Letters, Feb.1997,7(2):48~50

[12] R.Holland and J.Williams．Total field versus scattered field finite difference codes:a comparative assessment[J].IEEE Transactions on Nuclear Science，Dec.1983，NS-30(6)：4583～4588

[13]G.Mur．Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic field equations[J]．IEEE Transactions on Electromagnetic Compatibility，Nov.1981，EMC-23(4)：377～382

[14]B.Engquist and A.Majda. Absorbing boundary conditions for the numerical simulation of waves [J].Mathematics of the Computation，July 1977，31(139)：629～651

[15]J.P.Berenger．A perfectly matched layer for the absorption of electromagnetic waves[J]．Journal of Computational Physics，Oct.1994，114(2)：185～200

[16]J.P.Berenger．Perfectly matched layer for the FDTD solution of wave-structure interaction problem[J]．IEEE Transactions on Antennas Propagation，Jan.1996，AP-44(1)：110～117

[17]J.P.Berenger．Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[J]．Journal of Computational Physics，Sept.1996，127(2)：363～379

[18]Z.S.Sacks，D.M.Kingsland，R.Lee and J.F.Lee．A perfectly matched anisotropic absorber for use as an absorbing boundary condition[J]．IEEE

Transactions on Antennas and Propagation，Dec.1995，AP-43(12)：1460～1463

[19]D.E.Merewether.Transient currents on a body of revolution by an electromagnetic pulse [J]. IEEE Transactions on Electromagnetic Compatibility，May 1971，EMC-13(2)：41～44

[20]R.Holland．THREDS: a finite-difference time-domain EMP code in 3D spherical coordinates [J]．IEEE Transactions on Nuclear Science，Dec.1983，NS-30(6)：4592～4595

[21]S.S.Zivanovic，K.S.Yee and K.K.Mei．A subgridding method for the time-domain finite-difference method to solve Maxwell’s equations [J]．IEEE Transactions on Microwave Theory and Techniques，Mar.1991，MTT-39(3)：471～479

[22]P.Thoma and T.Weiland．A consistent subgridding scheme for the finite difference time domain method[J]．International Journal of Numerical Modeling: Electronic Network Device and Fields，1996，9(5)：359～374 [23]K.S.Yee，D.Ingham and K.Shlager．Time-domain extrapolation to the far field based on FDTD calculations[J]．IEEE Transactions on Antennas and Propagation，Mar.1991，AP-39(3)：410～413

[24]R.J.Luebbers，K.S.Kunz，M.Schneider and etc.．A finite-difference time-domain near zone to far zone transformation[J]．IEEE Transactions on Antennas and Propagation，Apr.1991，AP-39(4)：429～433

[25]R.J.Luebbers，D.Ryan and J.Beggs．A two dimensional time-domain near-zone to far-zone transformation[J]．IEEE Transactions on Antennas and Propagation，July 1992，AP-40(7)：848～851

[26]M.Kragalalott，M.S.Kluskens and W.P.Pala．Time-domain fields exterior to a two dimensional FDTD space[J]．IEEE Transactions on Antennas and Propagation，Nov.1997，AP-45(11)：1655～1663

[27]冯慈璋，马西奎. 工程电磁场导论 [M]. 高等教育出版社 2003 [28] D.M.Sullivan．Electromagnetic Simulation Using the FDTD Method[M]．New York：IEEE Press，2000

[29]B.J.Zook. A Dispersive Outer Radiation Boundary Condition for FDTD Calculation. Proceedings of the 10th Annual Review of Progress in Applied Computational Electromagnetics. 1994,1:240～247

[30]Engquist B, Majda A. Absorbing Boundary Conditions for the Numerical Simulation of Waves. Mathematics of Computation. 1977,31(7):629～651 [31]T.G.Moore,J.G.Blaschak,A.Taflobe,G.A.Kriegsmann. Theory and Application of Radiation Boundary Operators.IEEE Trans. AP.1988,36(12)1797～1812

本文来源：https://www.bwwdw.com/article/cu63.html