《计算电磁学》第五讲

y 2
1 c2
2 t 2
Dx2
Dy2
1 c2
Dt2
(5-2)
于是方程(5-1)可写为
LU 0
(5-3)
2020年6月11日1时57分
算子L还可以通过因式分解写为
L LU 0
其中，
L
Dx
Dt c
1 s2
L
Dx
Dt c
1 s2
s Dy Dt / c
2020年6月11日1时57分
(5-4)
(5-5) (5-6) (5-7)
2U zt
1 c
2U t 2
c 2
2U x2
c 2
2U y 2
0，z=h边界
(5-28)
对于矢量Maxwell方程的FDTD仿真，近似吸收边界条件(5-23)-(5-28) 中的U表示位于网格边界上的E和H的各个切向分量。
2020年6月11日1时57分
Mur差Βιβλιοθήκη Baidu格式
对于上述一阶、二阶近似解析吸收边界条件，Mur提出了一种简单有效 的差分数值算法。利用它们来截断FDTD仿真区域，总体虚假反射在1%-5%。
(1, j)
(0, j)
(1, j )
(0, j)
2t x
x
(5-34)
2020年6月11日1时57分
将该式对时间t的偏导数写为相邻两点(0, j)和(1, j)处对时间偏导的平均， 其表达式为
2W n
1
2W
n
2W n
t 2 ( 1, j) 2
2 t 2 (0, j)
t 2
(1, j)
2U xt
1 c
2U t 2
c 2
2U y 2
c 2
2U z 2
0
(5-23)
当s很小时，(5-23)是准确吸收边界条件LU 的很好近似。
对于其他网格边界相应的二阶近似解析吸收边界条件为
2U
1 2U
c
2U
c
2U
0
，x=a边界
xt c t 2 2 y2 2 z2
(5-24)
2020年6月11日1时57分
2020年6月11日1时57分
Reflection wave
x
o
a
相似地，算子L 作用于波函数，将构成x=a处的准确吸收边界条件。
对式(5-5)和(5-6)中根式的处理，可以用Taylor级数展开。
将 1 s2 在s=0附近展开为Taylor级数，
1 s2 1 s2 ... 2
当s很小时，只取一项，
W
n (1, j,k 1)
2W
n (1, j,k )
W
n (1, j,k 1)
2
(z)2
(z)2
(5-38)
将这些差分表达式代入式(5-23)，解出W
n1 0, j,k
就得到三维情形下W分量在
x=0网格边界出的二阶吸收边界条件
2020年6月11日1时57分
W n1 W n1 ct x W n1 W n1
2U 1 2U c 2U c 2U 0 ，y=0边界 yt c t 2 2 x2 2 z2
(5-25)
2U yt
1 c
2U t 2
c 2
2U x2
c 2
2U z 2
0，y=b边界
(5-26)
2U 1 2U c 2U c 2U
zt c
t 2
2
x2
2
y 2
0 ，z=0边界
(5-27)
2W y 2
n (1, j) 2
1
2W
2 y2
n (0, j)
2W y 2
n
(1, j)
1
W
2
n (0, j1)
2W n (0, j) (y)2
W
n (0, j1)
W
n (1, j1)
2W n (1, j ) (y)2
W
n (1, j1)
(4)
Substituting
(5-34)-(5-36)
© Mei-Fang超吸收边界条件(1992) ；
© PML完全匹配层
2020年6月11日1时57分
Engquist-Majda吸收边界条件 考虑二维情形时的齐次波动方程，
2U x2
2U y 2
1 c2
2U t 2
0
(5-1)
其中U为标量场分量，c为波的相速度。定义偏微分算子
L
2 x2
2
以一阶情形，x=0边界为例。在 x 1 x 处、t (n 1/ 2)t 时刻
2
用中点差分来代替式(5-12)中的偏微分，得
U
1 2
x, (n x
1 )t 2
1 x
n1 U 2
n1 (1) U 2 (0)
(5-29)
2020年6月11日1时57分
1 c
U
1 2
x, (n t
1 2
)t
从式(5-23)出发，此时前面的网格图2位于z kz 网格平面，
在距离边界半个步长的辅助网格点(1 , j, k) 处，对分量 W n (1 , j, k)
2
2
用中点差分代替(5-23)中的偏导运算。
偏导 2W
xt
、2W
t 2
和
2W y 2
与式(5-34)-(5-36)相同。只是在 z
kz
平面进行计算。
U
n1 (0)
U
n
(1)
ct ct
x x
U
n 1 (1)
U
n
(0)
(5-33)
这就是Mur一阶吸收边界条件(Mur’s 1st order ABC)。 下面推导Mur二阶吸收边界条件。仍旧以x=0边界为例，如图2所示。
2020年6月11日1时57分
W (0, j 1)
W (1, j 1)
W (0, j)
in
(5-23),
one
can
obtain
W
n1 (0, j)
.
Thus, the 2nd ABC at x=0 boundary can be got, which is
(5-36)
2020年6月11日1时57分
W n1 W n1 ct x W n1 W n1
(0, j,k )
(1, j) ct x (1, j)
2020年6月11日1时57分
有哪些吸收边界条件？
© 基于Sommerfeld辐射条件的Bayliss-Turkel吸收边界条件； © 基于单向波动方程的Engquist-Majda吸收边界条件(1977年)； © Mur吸收边界条件(1981) ； © Trefethen-Halpern近似展开(1986) © Higdon算子；
L
2 x2
2 y 2
2 z 2
1 c2
2 t 2
Dx2
Dy2
Dz2
1 c2
Dt2
(5-20)
将L分解为L和 L ,得到与(5)和(6)相同的准确吸收边界条件算子。
不同的是，s为
s
Dy Dt /
c
2
Dz Dt /
c
2
(5-21)
2020年6月11日1时57分
L 算子作用于波函数U，将在网格左边界x=0处准确地吸收以任意角度 入射到边界的平面波。
利用Taylor级数近似展开式(5-10)，可得到x=0处的 一阶吸收边界条件，其形式与(5-12)相同。
利用Taylor级数近似展开式(5-13)，可得到x=0处的二阶吸收边界条件。 其表达式为
Dx
Dt c
cDy2 2Dt
cDz2 2Dt
U
0
(5-22)
2020年6月11日1时57分
两边同乘以Dt ，得
1 ct
U
n1( 1 ) 2
U
n
(
1 2
)
(5-30)
其中半网格点和半时间步长时刻的值，可用下列二阶精度的平均公式计算
U
n1 2
(m)
1 2
U
n 1 (m)
U
n
(m)
U
n
(m
1) 2
1 2
U
n
(m
1)
U
n
(m)
(5-31) (5-32)
2020年6月11日1时57分
Substituting (5-31), (5-32) in (5-29), (5-30) and (5-12), after the manipulation, gets
具体为，在距离网格边界半步长的辅助网格点(1 , j) ，对W n (1 , j) 分量,
2
2
将(5-23)中的偏微分用中心差分来代替。
将式中关于x、t的偏导用中心差分格式写为
2W n
1
W
n1
W n1
xt ( 1, j) 2
2t
x
(1, j) 2
x
(
1, 2
j)
1
W n1 W n1 W n1 W n1
1 s2 1 s2 2
Substitution (5-13) in (5-5) yields
L
Dx
Dt c
1
1 2
s2
Dx
Dt c
1
1 cDy
2
Dt
2
Dx
Dt c
cDy2 2Dt
(5-13) (5-14)
2020年6月11日1时57分
Substituting (5-14) in (5-8) has
1 2U c t 2
c 2U 2 x2
0
，y=0边界（下边界）
(5-17)
2U yt
1 c
2U t 2
c 2
2U x2
0
，y=b边界（上边界）
(5-18)
考虑三维情形时的齐次波动方程，
2U 2U 2U 1 2U
0
x2 y2 z2 c2 t 2
(5-19)
2020年6月11日1时57分
此时，偏微分算子
(0, j)
2x W n W n
ct x (1, j)
(0, j)
(ct)2 x 2(y)2 (ct x)
W
n (0, j1)
2W
n (0, j)
W
n (0, j 1)
W
n (1, j1)
2W
n (1, j )
W
n (1, j1)
(5-37)
2020年6月11日1时57分
对于三维情形，考虑x=0边界。
在网格边界，如x=0处，将算子L作用于波函数将完全吸收以任意角度 入射
到边界的平面波，即将
LU 0
(5-8)
用于图1中的边界x=0，可构成一个准确的解析吸收边界条件。 它将吸收来自区域Ω内的波。
y
b
Incident wave
i
图1 二维吸收边界条件
Fig.1. Two-dimensional ABC.
2020年6月11日1时57分
© 利用差值技术的廖氏吸收边界条件
廖氏吸收边界条件比Mur二阶吸收边界条件在网格外边界引起的 反射要小一个数量级(20dB)，对外向波传播角度或数值色散均不敏感， 并且在矩形计算区域的角点处也易于实现。
M. Moghaddam, R. L. Wagner和W. C. Chew (周永祖) 曾指出，如果 采用单精度计算，可能导致使用廖氏吸收边界条件的FDTD算法不稳定， 而采用双精度则可改善稳定性。
2U 1 2U c 2U
xt c
t 2
2
y 2
0
(5-15)
这就是所分析区域左侧边界x=0的二阶近似吸收边界条件。 对于图1中的其他边界，相应的二阶近似解析吸收边界条件为
2U 1 2U c 2U
xt c
t 2
2
y 2
0
，x=a边界（右边界）
(5-16)
2020年6月11日1时57分
2U yt
2020年6月11日1时57分
偏导 2W z 2
可表达为在(0,
j,
k)和(1,
j,
k)处对z的偏导的平均值，其为
2W n
1
2W
n
2W n
z 2 (1, j,k) 2
2 z2 (0, j,k )
z 2
(1, j,k )
1
W
n (0, j,k 1)
2W
n (0, j,k)
W
n (0, j,k 1)
第五讲 吸收边界条件
Dr. Ping Du(杜平) E-mail: pdu@hfut.edu.cn
HeFei University of Technology
2020年6月11日1时57分
为何要用吸收边界条件？
由于计算机内存空间有限，让我们的分析区域不能无限大，必须在某处 截断。另外，即使内存空间很大，也必须在某处截断。因为大的分析空间就 意味着巨大的计算时间。从节省计算时间的角度考虑，也必须截断。
W (1 , j)
2
X
W (1, j)
W (0, j 1)
W (1, j 1)
y ( j 1)y y jy y ( j 1)y
x 0 x x
图2 Mur吸收边界条件
Fig. 2. Mur’s 2nd ABC.
图中W n (0, j) 表示位于x=0网格边界的E或H的切向分量。
2020年6月11日1时57分
(0, j,k )
(1, j,k ) ct x
(1, j ,k )
(0, j,k )
2x W n W n
ct x
(1, j ,k )
(0, j,k )
(ct)2 x 2(y)2 (ct x)
W
n (0, j 1,k )
2W
n (0, j,k )
W
n (0, j 1,k )
W
n (1, j 1,k )
2W
n (1, j,k )
W
n
(1, j 1,k )
(ct)2 x 2(z)2 (ct x)
1
W
2
n1 (0, j)
2W n (0, j) (t ) 2
W
n1 (0, j)
W
n1 (1, j)
2W n (1, j) (t ) 2
W
n1 (1, j)
(5-35)
2020年6月11日1时57分
(3) 将该式对y的偏导数写为相邻两点(0, j)和(1, j)处对y偏导的平均， 其表达式为
1 s2 1
(5-9)
(5-10)
2020年6月11日1时57分
将(5-10)代入(5-5)中，有
L
Dx
Dt c
将其代入(5-8)，可得
U 1 U 0 x c t
(5-11)
(5-12)
这就是所分析区域左侧边界x=0的一阶近似吸收边界条件。
2020年6月11日1时57分
将(5-9)中的级数取两项，有