PML吸收边界条件下TM波的时域有限差分法分析

1. 引 言
目前，电磁场的时域计算方法越来越引人注目。时域有限差分（Finite Difference Time Domain，FDTD）法作为一种主要的电磁场时域计算方法，最早是在 1966 年由 K. S. Yee 提出的。这种方法通过将 Maxwell 旋度方程转化为有限差分式而直接在时域求解，通过建立 时间离散的递进序列，在相互交织的网格空间中交替计算电场和磁场。经过三十多年的发展， 这种方法已经广泛应用到各种电磁问题的分析之中[1]。
http://www.paper.edu.cn
PML 吸收边界条件下 TM 波的时域有限差分法分析
马晓晖 1，刘丽卓 2，赵维 2
1.北京邮电大学光通信中心，北京（100876） 2.陕西渭南师范学院计算机科学系，陕西渭南（714000）
E-mail：mxhui123@gmail.com
摘 要：介绍了时域有限差分法的基本原理，通过场区的 FDTD 公式分析了 PML 边界及其 相关参数设置的要求。最后以自由空间中 TM 波传播为实例，利用 C 语言编写了 FDTD 程 序。并对结果进行了详细分析。 关键词：完全匹配层，TM波，时域有限差分法 中图分类号：TN252
2. 场区 FDTD 公式的讨论
∂E z = − µ ∂H x − σ mH x
∂y
∂t
⎫ ⎪ ⎪
∂E z = µ ∂H y + σ mH y
∂x
∂t
⎪ ⎬
TM 波
⎪
(1)
∂H y ∂x
−
∂H x ∂y
=
ε
∂E z ∂t
+
σ
E
⎪ z⎪
⎭
-1-
http://www.paper.edu.cn
对于 TM 波， Ex = Ey = Hz =0，FDTD 公式为：
电导率分布阶数的改变不会影响计算量，却会影响吸收边界效果。选择常数电导率分布
会带来较强的数值反射，选取线性分布可以改善结果。通常，选取 n=2 或 3，可进一步改进 吸收效果。
4. TM 波实例分析
4. 1 问题提出
二 维 TM 波 在 自 由 空 间 传 播 ， 采 用 PML 吸 收 边 界 ， 激 励 源 使 用 脉 冲 源
以立方指数形式增大。将(14)代入(13)得：
∫ R(θ ) = exp(− 2 cosθ
ε 0c0
d 0
σ
( ρ )dρ )
=
exp⎨⎧− ⎩
2dσ max cosθ (n + 1)ε 0c0
⎫ ⎬ ⎭
(15)
当θ =0(垂直入射)时，
R (0)=
exp⎨⎧− ⎩
2dσ max (n + 1)ε 0c0
4. 3 仿真分析
1)T=20 步的情况
-4-
http://www.paper.edu.cn
2)T=40 步的情况
图 4 T=20 步时场传播情况
3)T=80 步的情况
图 5 T=40 步时场传播情况
图 6 T=80 步时场传播情况
图 6 看上去出现了反射，但是我们应该看到，这些反射都是出现在 PML 区域（共有八 层），而在场区没有出现任何反射，右侧图可以证明这一点。 3) T=120 步的情况
-5-
http://www.paper.edu.cn
图 7(a) T=120 步时场传播情况(有吸收边界)
图 7(b) T=120 步时场传播情况(无吸收边界)
由脉冲源的性质可知，T=120 步时，场区的左下部分应该已经没有电场，即电场值为零， 而在左上部分仍可以看到场在传播，这是脉冲源没放在中心处的原因。图 7(a)正反应了这一 情况。但在图 7(b)中，左下部分出现了负的电场值，即产生了反射波，这正是因为没有吸收 边界导致的。
固有误差，通常选取 R (0)=1%。
(2)当 PML 层衰减值一定时，理论上可以将其厚度取得尽可能小(由(13)PML 表面处的反
射系数是乘积σ d 的函数可知)，但电导率跃变太大会带来数值反射；增加 PML 厚度，可以
使局部及总体误差都单调地减少，但 PML 层数过多会使计算量剧增。折衷考虑吸收边界效 果与计算量，通常选取 N 为 4-8 层。
(
i
,
j)
(3)
E
n+1 z
(i,
j)
=
CA(m)
⋅
E
n z
(i,
j)
+ CB
(m)
⋅
⎡ ⎢ ⎢⎣
H
n +1/ y
2
(i
+
1/
2
,
j) − ∆x
H
n +1/ y
2
(i
−
1/
2,
j)
-
H
n x
+1
/
2
(i,
j
+ 1/
2)
−
H
n +1/ x
2
(i,
∆y
j
−1/ 2) ⎤ ⎥ ⎦
(4)
CA(m)
=
ε (m)
传播到理想导体边界处会反射回来，重新回到 FDTD 区域。这样，PML 的反射系数不再等
于零。图 3 中侧边上的介质层参数为(σx ,σmx ,0,0,)和(0,0，σy ,σmy )，在距离内侧界面 ρ 处，
外行波的振幅为：
Ψ(ρ)
=
Ψ
0
exp(−
σ cosθ ε 0c0
ρ)
(12)
其中θ 为外行波相对于分界面的入射角；σ 表示σx 或σy 。穿过 PML 媒质后，波将被
理想导电壁反射，然后再次穿过媒质分界面进入自由空间，
因此，如果 PML 介质的厚度为 d,则 PML 介质内表面的反射系数为[6]：
R(θ ) = exp(− 2σ cosθ d )
(13)
ε 0c0
-3-
http://www.paper.edu.cn
由上式可知，当入射方向很接近分界面(即θ 趋近于 π/2)时，无论σ 为多少，R 都接近
由 FDTD 差分形式可知，边界网格点上的计算需要截断边界外场的信息，故也需要适 用于截断边界网格点处计算的算法。为了让实际的计算空间与实际的无限的物理空间等效， 须对有限空间的周围边界做特殊处理，使得向边界面行进的波在边界处保持“向外行进”特 点，无反射现象，即好像被吸收了一样，并且不会使内部空间的场产生畸变，具有这种功能 的边界条件称为吸收边界条件，吸收边界条件处理成为了 FDTD 的一大重点和发展方向， PML 是近年来发展起来的吸收边界条件[4－5]。
∆t ε (m)
− σ (m) 2
+ σ (m)
=
1− 1+
σ (m)∆t 2ε (m) σ (m)∆t
(5)
∆t 2
2ε (m)
∆t
CB(m)
=
ε(m)
1 +σ(m)
=
ε(m) 1+σ(m)∆t
(6)
∆t 2
2ε(m)
CP(m)
=
µ(m)
∆t µ(m)
− σm(m) 2
+ σm(m)
1− σm(m)∆t
本文将理想匹配层(PML ) 吸收边界条件用于 TM 波自由空间传播的时域有限差分
-6-
http://www.paper.edu.cn
(FDTD) 法分析中，给出了二维情况下的 PML 吸收边界条件,并对场区 FDTD 中心差分式进 行了推导， 最后结合二维 TM 波在自由空间传播的实例，编制程序，仿真分析了不同时间步数情况下， 空间电场的传播分布情况，并对结果做了详细的分析，所得结果与理论值非常一致, 证明了 PML 吸收边界条件应用于 TM 波自由空间传播分析的有效性，与传统的 Mur 吸收边界条件 相比, 可以大大提高计算的正确性和精确性。
Ei(t) =
exp(− 1 (t − t0 2 τ2
)2 ) ，其中 t0
=
20s ,τ
=
6s
。计算
Ez 值的空间分布，并分析场的情
况。
4. 2 模型建立
空间步长 ddx=0.01m，时间步长 dt=ddx/6e8=16.67ps,满足 courrant 条件。空间网格数 100*100，即为一个 1m*1m 区域中的脉冲源的传播情况。其中脉冲源置于网格坐标(30,30) 处。
因此，当 PML 介质采用以上参数时，通过 PML 的平面波以光速传播，且穿过边界和
角顶区域 PML 的分界面时均无反射。
在实际计算中，PML 介质层不可能延伸到半空间，只能是有限厚度。PML 层的外侧通
常采用理想导体(Perfect Electric Condition, PEC，其中电场分量为零)截断。透入 PML 层的波
H
n x
+1
/
2
(i
,
j
+ 1/
2)
=
CP (m)
⋅
H
n x
−1
/
2
(
i
,
j
+ 1/
2)
−
CQ (m)
⋅
E
n z
(i
,
j
+ 1) − ∆y
E
n z
(i
,
j)
(2)
H
n y
+1
/
2
源自文库
(i
+
1/
2,
j)
=
CP
(m)
⋅
H
n y
−1
/
2
(
i
+ 1/
2,
j)
+
CQ
(m)
⋅
E
n z
(
i
+ 1,
j) − ∆x
E
n z
⎫ ⎬ ⎭
(16)
我们使用 PML 时，首先要选定 3 个参数：PML 层数，电导率分布阶数 n, PML 内表面
反射系数 R (0)。大量的数值试验表明[7]：
(1)当 PML 厚度固定时，减少 R (0)，即增加 PML 的衰减，可以使局部及总体误差都单
调地减少。然而，当 R (0)小于 10-5 时，这种现象不再出现，原因是存在由数值网格引起的
ϖ 计算 t2=t1+ ∆ t/2 时刻空间各处 H 的值
ϖ 计算 t1=t2+ ∆ t/2 时刻空间各处 E 的值
图 2 二维情形 TM 波 FDTD 计算电磁场的逐步推进步骤
一般情况下， ∆ t 与 ∆ x, ∆ y 按下面的规则选取[3]： 1) 在每个最小波长中至少要包含 10 个网格。如： ∆ x<0.1 λ 。
2) ∆ t 可由 Courant 稳定性条件得到。即：
对于二维情形:
∆t≤
1
(9)
c0
1+1 (∆x)2 (∆y)2
若 ∆ x= ∆ y= δ ，上式化为:
∆t≤ δ
(10)
c0 2
在实际计算中，我们一般取:
∆ t= δ
(11)
2c0
-2-
http://www.paper.edu.cn
3. PML 吸收边界条件
外行波在四边的 PML 界面上是无反射的。
图 3 完全匹配层的设置
在四个角顶区域，介质参数为 PML( σx , σmx , σy , σmy ),而于之相邻的边上的参数为 PML(σx ,σmx ,0,0,)和 PML(0,0，σy ,σmy )。根据上面的讨论，满足匹配条件时，在角顶边和
侧边的界面上，如图 3 中的 BB1 和 BB2，平面波同样可以无反射地传播。
5. 总结
时域有限差分法(FDTD)作为研究电磁场传播的有效工具，在微波电路的时域分析、天 线辐射特性、电磁场散射及生物医学工程等领域有着广泛的应用。FDTD 法的一个关键问 题是吸收边界的处理。吸收边界的效果直接关系到 FDTD 计算的正确性和精确性, 是影响 计算品质的决定因素。J. P.Berenger 于 1994 年首先提出了二维理想匹配层(PML)吸收边界 条件的概念, 与传统的二阶 Mur 吸收边界条件相比, PML 吸收边界条件可提高精度 40dB 左 右, 是目前最好的吸收边界条件[8]。
时域有限差分法的主要思想是把 Maxwell 方程在空间、时间上离散化，用差分方程代替 一阶偏微分方程，求解差分方程组，从而得出各网格单元的场值。FDTD 空间网格单元上电 场和磁场各分量的分布如图 1 所示。
图 1 FDTD 网格中电磁场分量的分布
电场和磁场被交叉放置，电场分量位于网格单元每条棱的中心，磁场分量位于网格单元 每个面的中心，每个磁场（电场）分量都有 4 个电场（磁场）分量环绕。这样不仅保证了介 质分界面上切向场分量的连续性条件得到自然满足，而且还允许旋度方程在空间上进行中心 差分运算，同时也满足了法拉第电磁感应定律和安培环路积分定律，也可以很恰当地模拟电 磁波的实际传播过程[2]。
=
1+
2µ(m) σm(m)∆t
(7)
∆t 2
2µ(m)
∆t
CQ(m)
=
µ(m)
1 + σm(m)
=
1+
µ(m) σm(m)∆t
(8)
∆t 2
2µ(m)
上式中 CA 、 CB 、 CP 、 CQ 中标号 m 的取值与(2),(3)和(4)式左端场分量节点的空间位置
相同[3]。
由此分析可以得知：
ϖ
ϖ
已知 t1=t0=n ∆ t 时刻空间各处 E 分布及(n-1/2) ∆ t 空间各处 H 分布
于 1，这并不奇怪，如入射方向很接近左右侧分界面的波将最终几乎垂直地传到上或下侧边 界，在那里被吸收掉。
计算中 PML 层电导率变化通常采用下面形式：
σ (ρ ) = σ max( ρ )n n=0，1，2，3，为分布阶数
(14)
d
n=0 时为常数电导率；当 n=1 时，为线性变化；n=2 时，σ 以抛物线形式增大；n=3 时，
二维情况下 PML 吸收层区分为图 3 所示的边和角顶两种区域。设靠近 PML 层的 FDTD
区为真空，则 AB 与 CD 的 PML 界面垂直于 x 轴，其参数应为 PML(σx ,σmx ,0,0)。同样， BC 和 DA 的 PML 界面垂直于 y 轴，其参数应为(0,0,σy ,σmy )。这样，满足阻抗匹配条件时