FDTD法模拟无源区域二维TE波的传播

FDTD法模拟无源区域二维TE波的传播

摘要：自1966年Yee在其论文中首次提出时域有限差分以来，时域有限差分法在电磁研究领域得到了广

泛的应用。本文首先介绍了FDTD的基本原理，推导了无源区域二维TE波的FDTD差分方程公式，接着推导了一阶，二阶Mur吸收边界条件差分公式及四个角点的迭代公式，然后用matlab模拟了TE波的传播过程，最后简要地比较Mur一阶和二阶吸收边界条件。

关键词：FDTD,TE波，Mur吸收边界条件

1. 引言

FDTD自Yee(1966年)提出以来发展迅速，获得广泛应用。FDTD法以Yee元胞为空间电

磁场离散单元，将Maxwell方程转化为差分方程，表述简明，容易理解，结合计算机技术能处理十分复杂的电磁问题。在时间轴上逐步推进地求解，有很好的稳定性和收敛性，因而在工程电磁学领域倍受重视。

2. FDTD的基本原理

根据电磁场基本方程组的微分形式，若在无源空间，其空间中的媒质是各向同性、线性

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

??H???E??E?t （2-1a）

??E????H??mH?t （2-1b）

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

在直角坐标系中，可化为如下六个标量方程：

??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在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)?x?F(x,y,z,t)?y?F(x,y,z,t)?zFn(i?12,j,k)?Fn(i?12,j,k)2????O?xx?i?x?x?? (2-5a) ? (2-5b) ? (2-5c)

Fn(i,j?12,k)?Fn(i,j?12,k)2????O?yy?j?y?y?Fn(i,j,k?12)?Fn(i,j,k?12)2????O?zz?k?z?z??F(x,y,z,t)?tFn?12(i,j,k)?Fn?12(i,j,k)2????O?tt?n?t?t?? (2-5d)

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

3. FDTD差分公式

3.1 三维的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) （3-1a）

2222n?1n?1211?Hz(i?,j?,k)?Hz2(i?1,j?1,k)?y??? 2222 ???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) (3-1b) 2222n?1n?1211?Hx(i,j?,k?)?Hx2(i,j?1,k?1)?z??? 2222 ??11?n?n?221111?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) (3-1c) 2222n?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) (3-1d)

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) (3-1e) 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????n?1n?1Hz2(i?1,j?1,k)?CP(i?1,j?1,k)?Hz2(i?1,j?1,k)?CQ(i?1,j?1,k) (3-1f)

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????式中

CA(i,j,k)?2?t2?(i,j,k)??(i,j,k)?t， (3-1g) 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， (3-1h) CQ(i,j,k)?2?(i,j,k)??m(i,j,k)?t2?(i,j,k)??m(i,j,k)?tCP(i,j,k)?(3-1)式就是FDTD的基本差分方程组。从式中可以看出，方程组中含有半个空间步和半个时间步，为了便于编程，可将（3-3）式改写成如下形式：

Ex(i,j,k)?CA(i,j,k)Ex(i,j,k)?CB(i,j,k)?HZ(i,j,k)?HZ(i,j?1,k)Hy(i,j,k)?Hy(i,j,k?1)?(3-1i)

?????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)?(3-1j)

??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)?(3-1k)

?????x?y??Hx(i,j,k)?CP(i,j,k)Hx(i,j,k)?CQ(i,j,k)?E(i,j?1,k)?EZ(i,j,k)Ey(i,j,k?1)?Ey(i,j,k)?(3-1l)

??Z???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)?3-1m)

??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)?(3-1n)

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

3.2 二维TE波的FDTD差分方程

二维问题就是各个场分量只与两个坐标轴有关，而与第三个坐标轴无关，??z?0，

??Ex1?Hz??t??[?y??Ex]???Ey1?Hz?[???Ey]TE波的偏微分方程为：? （3-2a） ?t??x???Hz1?Ex?Ey?[???mHz]???y?x??t差分公式为：

11n?n??111122H(i?,j?)?H(i?,j?)?n?1zz11n2222?Ex(i?,j)?CA(m)Ex(i?,j)?CB(m)22?y?11?n?n?11112?Hz(i?,j?)?Hz2(i?,j?)2222?En?1(i,j?1)?CA(m)En(i,j?1)?CB(m)yy?22?x?111nn?E(i?,j?)?E(i,j?)11yyn?1111?n?22222??Hz(i?2,j?2)?CP(m)Hx(i?2,j?2)?CQ(m)[?x?11n?Exn(i?,j?1)?Ex(i?,j)?22]??y?

其中CA(m),CB(m),CP(m),CQ(m)参数与三维情况完全类似，m的取值与公式左边分量节点的取值相同。

?(m)?t?2?(m)?t2?CA(m)??(m)?(m)?(m)?t?1??t22?(m)?(m)?(m)1??(m)?t?2?(m)?t2?CP(m)??(m)?(m)?(m)?t?1??t22?(m)?(m)?(m)1?4. 时域有限差分法相关技术 4.1 数值稳定性问题

?t1?(m)CB(m)???(m)?(m)?(m)?t?1??t22?(m)

?t1?(m)CQ(m)???(m)?(m)?(m)?t?1??t22?(m)上述FDTD方程是一种显式差分方程，在执行时，存在一个重要的问题：即算法的稳定性问题。这种不稳定性表现为在解显式方程时，随着时间步数的继续增加，计算结果也将

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

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

?t??s (4-2)

v3其中，v为计算空间中的电磁波的最大速度。 4.2 数值色散

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

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

??c?22?kx2?ky?kz2 (4-4)

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

最小波长。由上分析说明，数值色散在用FDTD法分析电磁场传播中的影响是不可能避免的，但我们可以尽可能的减小数值色散的影响。 4.3 离散网格的确定

无论是简单目标还是复杂目标，在进行FDTD离散时网格尺寸的确定，除了受计算资源的限制不可能取得很小外，还需要考虑以下几个因素：

1．目标离散精确度的要求。网格应当足够小以便能精确模拟目标几何形状和电磁参数。 2．FDTD方法本身的要求。主要是考虑色散误差的影响。设网格为立方体

?x??y??z??，所关心频段的频率上限为f值色散要求

max，对应波长为?min，则考虑FDTD的数

???min （4-5）

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

3．入射波的要求。入射波的上限截止频率5. 吸收边界条件

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

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

fc应包含所关心频率范围，即fc?fmax。

??2?2?21?2????x2??y2??z2?c2?t2??f?0 (5-1) ??式中，

f表示直角坐标系下任意电磁场分量。

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

一直角场分量。G. Mur对表中的吸收边界条件引入了一种简单有效的差分数值算法，即对时间和空间的偏微分取二阶中心差分近似，将单向波方程离散化，便形成了著名的G.Mur的一阶和二阶吸收边界条件，其总体虚假反射在1%～5%。图5-1给出了反射系数（反射波与入射波）与入射角的关系，由图5-1可看出，仅当入射角较小时其反射系数较小。

表5-1 三维长方体FDTD区域的一阶和二阶吸收边界条件

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

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

根据Yee氏元胞网格的特点，在FDTD截断边界面上只有电场切向分量和磁场法向分

量。以x?0界面为例，此界面仅有Hx，Ey，Ez节点。由于FDTD中Hx的计算式不涉及x?0区域，即不涉及截断边界界面外的节点。所以，吸收边界条件将不考虑Hx，而只考虑电场切向分量Ey和Ez。以Ez为例，对一阶和二阶吸收边界条件分别有

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

?????EZ???0 (5-2b)

?1?21?21??2?2?22?????y2??z2c?x?tc?t2??x?0将式(5-2b)分别在距离边界半个空间步长的辅助网格点?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)? (5-3)

c?t??式中

f代表截断边界上切向场分量。

将式(5-3)所示二阶吸收边界推广到长方体元胞各边?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??x

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

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

对于二维电磁场问题，在用FDTD计算边界处的元胞时，将涉及到截断边界外侧的节点。对于TE波只需给出边界处切向场分量Hz的吸收边界条件。表5-2为矩形截断边界面上节点的二阶Mur吸收边界条件。

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

截断边界位置 二阶Mur吸收边界条件 TE x?0 ??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 ??HZ1?HZc??Ey???y?c?t?2?x??0 ??y?0??HZ1?HZc??Ey???y?c?t?2?x??0 ??y?b在二维矩形计算区域的角点，吸收边界条件的离散式需特殊考虑。假设角点附近只有向外传播的行波，且传播方向沿角点处元胞的对角线，如图5-2所示的矩形区域左下角点处的TM元胞为例，导出适用于角点的吸收边界条件离散式。根据Courant稳定条件c?t??式，有2c?t2?2?。在点

?i,j?与对角点?i000?1,j0?1?之间取一点P，对于沿对角

线的外行波，有以下等式：

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

n?2n?2n?2n?2????i0?1,j0?1?（5-6） EZi0,j0??EZP?EZi0,j0??EZ ?2c?t2?由式（5-5）解出EZ?n?i0,j0???EZ?1??n?2?P?后代入（5-4）式得

2c?t2c?t?n?2??EZ?i0,j0?????n?2?i0?1,j0?1?EZ （5-7）

或

n?1n?i0,j0??EZ?i0?1,j0?1??EZc?t?c?t?2??EZn?1?i0?1,j0?1??EZn?i0,j0??（5-8） 2?对于矩形域的其它角点可作类似的处理。对于二维TE波情况，将式（5-6），（5-8）式中的EZ

换为HZ即可。

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

图5-2 矩形域四个角点

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

为了使计算达到稳定,通常计算所需要时间步按照电磁波往返穿越FDTD计算区对角线3～5次来估计。若FDTD计算区总元胞数为N，则对角线上元胞约为按照Courant稳定条件，设计算中心区

3N13 （三维)。

?t??(2c)，即穿越对角线一次需要时间步为

23N13。总计算时间步约需23?(6~10)N13步。对于二维情况则约为

22?(6~10)N12。或者说，计算时间步大约等于FDTD计算区对角线上元胞数目的12～

20倍。实际上，计算所需时间步还与散射体具体形状、结构有关。图5-3给出了应用FDTD分析电磁场问题时的程序流程图

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

图5-3 FDTD 程序流程图

6. MATLAB仿真模拟出图及结果

6.1 Mur一阶吸收边界条件

6.2 Mur二阶吸收边界条件

、

相对于Mur一阶吸收边界条件,Mur二阶吸收边界条件具有较高的计算精度,但需要略多一些的存储空间,实现难度也略大. 不过这里存储量和实现难度的增加其实是非常小的。

6.3 金属柱对TE波传播的影响

电磁波在金属导体中不能传播，TE波在传播的过程中遇到金属柱体后，将会绕开金属柱进行传播。

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