电磁场数值分析 第16讲

第16讲 边界条件
为了描述方便，令
L− = Dx − Dt 1− S 2 c
；L+ = Dx + Dt
c
1− S 2 ∂ ∂t
S=
Dy Dt / c
；Dx =
∂ ∂x
；Dy =
∂ ∂y
(1-15)
； Dt =
首先考察 L− 中的根式。当S很小时，保留一项的泰勒 级数为 1 − S 2 ≅ 1 , S很小意味着同外行波相对于c规一化 的时间偏导相比，其 y 方向偏导数可以忽略不计。换言 ˆ 方向朝x=0边界传播（近轴 之，外向波以非常接近于 − x 近似），对 L− 中的根式采用一阶近似可得：
x =0
=0
(1-12)
根据时频域对应关系：
jk → 1∂ c ∂t
jk y →
∂ ∂y
第16讲 边界条件
则上式变换到时域为:
∂ ∂ 2 − j k2 − ky = − ∂x ∂x
2 → j 2k 2 − j 2k y
∂ ∂ 1 ∂ − 2 2 − 2 ∂x ∂y c ∂t =0
(1-13) (1-14a)
第16讲 边界条件
波在边界处保持“外向行进”的特征，以模拟电磁波无反 射地通过截断边界，向无限远处传播。理想的ABC是难以 实现的，通常只能采用近似的ABC。对于近似ABC，要求： (1) 能够模拟向外传播的波；(2) 引入的反射应足够小， 对计算结果的影响可忽略； (3) 保证算法稳定。 我们可以将大多数的吸收边界分为两大类：一种是由微 分方程推导出的吸收边界，另一种是由吸收媒质构成的吸 收边界。Mur提出的吸收边界属于第一种吸收边界，1994 年 Berenger 提出的理想匹配层 (PML) 与 1996 年 Gedney 提出 的UPML属于第二种吸收边界。PML吸收边界目前在应用 中较之其他的所有吸收边界具有很高的精度。
2 ∂ 1 ky ( − jk − )f ∂x 2 jk
(1-25) (1-26)
x =0
= 0 ⇒ ( jk
∂ 1 2 + k 2 − ky )f ∂x 2
x =0
由于 L− f = 0 ，可得
L− f −
和
L− f +
： (1-27)
2 k 1 ∂ ∂ 1 ∂2 ∂2 2 ( − 2 2+ ) f − = ( −k k 2 − k y + k 2 − y ) f− 2 c ∂x ∂t c ∂t 2∂y 2
第16讲 边界条件
FDTD方法中三种最常用的边界条件是PEC、PMC以 及ABC （吸收边界条件）。PEC或者 PMC边界条件比较 容易实现，PEC简单的在边界处把切向电场设置为0，而 PMC则把切向磁场设置为0即可实现。 时域有限差分最重要也是研究最多的问题之一是如何 截断开域问题的计算区域，即吸收边界条件（Absorbing Boundary Condition ABC）。由于FDTD计算时，每个单元 网格上的六个场分量均需在任一时间步上贮存起来供下一 步时间计算之用。因此所取的问题空间愈大，要求存储量 也愈大，很难想象计算机的存储量是无限的。因此问题空 间是有限的，要求它能将被研究的模型“装入”，并实施 FDTD的运算过程。为了让这种有限空间与无限空间等效， 需对有限空间的周围边界进行处理，使得向边界面行进的
第16讲 边界条件
表1-1 二维截断边界条件的近似
边界 一阶近似 二阶近似
x=0
∂ 1 ∂ ( − )f ∂x c ∂t
∂ 1 ∂ ( + )f ∂x c ∂t
x =0
2 2 2 1 f f c f ∂ ∂ ∂ =0 ( ) =0 + − 2 2 x =0 2 ∂y ∂x∂t c ∂t
x=a
x=a
− +
(1-32)
将 L 分解为 L ，L 得到与（1-15）类似的准确吸收边界 算子，不同的S由下式给出
S= ( Dy Dt / c )2 + ( Dz 2 ) Dt / c
(1-33)
第16讲 边界条件
利用泰勒级数展开，可得泰勒一阶展开x=0处与二维结果 相同。 在x=0处二阶展开
∂2 f 1 ∂2 f c ∂2 f =0 ( ) =0 − + 2 2 x =a 2 ∂y ∂x∂t c ∂t
∂2 f 1 ∂2 f c ∂2 f ) =0 + − =0 ( ∂y∂t c ∂t 2 2 ∂x 2 y =0 ∂2 f 1 ∂2 f c ∂2 f =0 ( ) =0 − + ∂y∂t c ∂t 2 2 ∂x 2 y =b
L+ = (
记:
∂ 2 ) − j k 2 − ky ∂x
∂ 2 + j k2 − ky ) ∂x
将(1-5)式中左向行波 f − 和右向行波 f 分别与算子 作用，可得: L− f − = 0 (1-10) + +
L f = 0
第16讲 边界条件
+ − L L 因此称 为左行波算子， 为右行波算子 − 若将左行波算子 L 作用在平面波(1-4)上，可得:
| f+ f− |
0
0
5
0.0019
15
30
45
85
一阶近似
0.0173
0.0718
0.172
0.84
二阶近似
0
3.6 × 10 −6
3.0 × 10 −4
0.0052
0.029
0.71
第16讲 边界条件
1.0 0.8 0.6 0.4 0.2 0.0 一阶 近似 二阶 近似
f+
f− |
第16讲 边界条件
本章将讨论目前应用最广泛的 Mur 吸收边界条件， Berenger完全匹配层以及单轴各向异性完全匹配层吸收 边界条件。
1. Engquist-Majda吸收边界条件
如果一个偏微分方程允许波沿一定的方向传播，则称 它为单向波方程。当在FDTD网格外边界处应用单向波 方程时，可以吸收外向散射或辐射波。 Engquist和Majda导出了适合直角坐标FDTD网格吸 收边界条件的单向波方程，他们的单向波方程可以用偏 微分算子分解因式得到，以直角坐标系中的二维波动方 程为例： ∂2 f ∂2 f 1 ∂2 f + 2 − 2 2 =0 (1-1) 2
x =0 2 = ( j k2 − ky −j
ω
c
2 −j ) f − + (− j k 2 − k y
ω
c
) f+
x =0
=0
即 令
k y = k sin θ
(1-18)
f+ f−
=
x =0
2 k 2 − ky −k 2 k 2 − ky +k
(1-19)
θ 为入射角，则 ，
f+ f− =
x =0
第16讲 边界条件
2 ky ∂2 1 ∂ ∂ 1 ∂2 2 2 2 − 2 2 + ( ) f + = (k k − k y + k − ) f + 2 c ∂x ∂t c ∂t 2∂y 2
(1-28) (1-29)
f+ f−
=−
x =0
2 − k k 2 − ky + k2 −
2 ky
2 2
k k −k +k −
∂x ∂y c ∂t
第16讲 边界条件
f是一标量场分量，波的相速度为c的平面波解为： j(ωt - k x x - k y y ) (1-2) f(x, y, t) = Ae
k = k +k =
2 2 x 2 y
ω2
c2
(1-3)
y
f−
入射波
θ θ
x=0 f +
反射波
x
图1-1 截断边界
第16讲 边界条件
cos θ − 1 cos θ + 1
(1-20)
第16讲 边界条件
这就是一阶近似式作用在 x ＝ 0 界面后所残留的反射波 与入射波之比（反射系数）。进一步考察泰勒级数中保 留两项的情形，此时
1− S 2 ≅ 1− 1 2 S 2
(1-21)
显然，在相同精度要求下，这一近似与一阶近似相比
ˆ 可取稍大的S，或者说允许向网格边界x=0传播的波同 − x Dy 方向的夹角稍大一些，因为 的比值可以更大一些 ,总
2 2 y 2
2 ky
令
k y = k sin θ
，则：
f+ f− sin 2 θ − Fra Baidu bibliotekos θ + 1 − 2 =− sin 2 θ cos θ + 1 − 2
(1-30)
x =0
这就是二阶近似吸收边界条件作用后，边界处残留 的反射波与入射波之比。
第16讲 边界条件
我们将两种近似情况下所得 | f + f − | 列在表1-2中， 从表1-2可见，二阶近似比一阶近似性能有所改善。图 1-2中进一步给出了该比值与入射角的关系，可见，入 射角越小反射系数越小。 表1-2 边界处反射波与入射波之比
+
设 0<k
y
<
ω
c
, x≥0
, 将上式第一项记为 f ，第二项记为 f ，即: j ( wt + k − k x + k y ) (1-5a) f − = A− e
2 2 y y
f = A+ e
+
j ( wt − k 2 − k 2 y x +k y y )
(1-5b)
−
对所讨论的左侧界面x=0处而言，位于x>0的 f 为左向 + 行波，代表入射波；而 f 为右向行波，代表反射波。 将 (1-4) 代入 (1-1) ，因为考察 x 方向的截断面，故保留 对x的导数，可得
体看两者 可以更接近，这样可以得到
2 Dt cD y + L− ≅ D x − c 2 Dt
Dt
(1-22)
第16讲 边界条件
通过
L− f = 0
可得到偏微分方程
∂2 f 1 ∂2 f c ∂2 f − + =0 ∂x∂t c ∂t 2 2 ∂y 2
(1-23)
它可以在左边界 x=0 处数值实现，对于小 S 值，泰 勒级数展开形式 (1-21)式成立，式(1-23)是准确吸收边 界条件式 L− f = 0 的很好近似，即对于比较小入射角 传向网格外边界x=0的数值平面波，式(1-23)代表一种 近于无反射的网格截断。
可得:
1 ∂2 ∂ ∂2 ( − 2 2 − 2 )f ∂x ∂y c ∂t
x =0
同理，在右侧边界相应公式为:
∂ ∂2 1 ∂2 ( + 2 2 − 2 )f c ∂t ∂x ∂y
x =b
=0
(1-14b)
这就是Engquist和Majda给出的吸收边界条件。 由于 L− , L+ 中含根式运算这种特性限制其直接用于数值 实现吸收边界条件。为此可采取许多种算子的近似形式。
y=0
∂ 1 ∂ ( − )f ∂y c ∂t
∂ 1 ∂ ( + )f ∂y c ∂t
y =0
y=b
y =b
第16讲 边界条件
变换到频域可以得到：
ky ∂ − L = − jk 1 − k ∂x
2
(1-24)
]
=0
2
近似后：
1 ky ∂ − − jk[1 − L = 2 ∂x k
∂2 f 2 2 k k + ( − y)f = 0 2 ∂x
频域
(1-6)
第16讲 边界条件
∂2 2 2 定义 L = ∂x 2 + (k − k y ) ，则上式可写为
Lf = 0
(1-7) (1-8) (1-9a) (1-9b)
+
将 L 进行因式分解:
L=( ∂ ∂ 2 2 )( + j k 2 − k y ) − j k2 − ky ∂x ∂x L− = (
0
10
20
30
40
50
60
70
80
90
入射角度( Deg)
图1-2 一阶以及二阶近似时反射波与入射波之比
第16讲 边界条件
对于三维波动方程为：
∂2 f ∂2 f ∂2 f 1 ∂2 f + + − =0 ∂x 2 ∂y 2 ∂z 2 c 2 ∂t 2
(1-31)
相应的偏微分算子为：
∂2 ∂2 ∂2 1 ∂2 1 2 L = 2 + 2 + 2 − 2 2 = D x2 + D y + D z2 − 2 Dt2 ∂x ∂y ∂z c ∂t c
设 x ＝ 0 平面为截断边界，如图 1-1 所示，图中 x>0 区 域同时存在有入射波与反射波，因此在此区域中有
f ( x, y , t ) = A− e
2 j ( wt + k 2 − k y x +k y y )
+ A+ e
−
2 j ( wt − k 2 − k y x +k y y )
(1-4)
L− ≅ D x − Dt c
(1-16)
第16讲 边界条件
将其应用于左边界的波，可得一阶近似解析吸收边 界条件为： ∂ 1 ∂ )f =0 L− f = ( − (1-17) ∂x c ∂t x =0 同理可得右边界的一阶近似吸收边界条件。下面我们检 查（1-17）式的近似程度
( ∂ 1 ∂ − )f ∂x c ∂t
L− f = L− f − + L− f
+
= L− f
+
(1-11)
其结果带有与右向行波相关联的部分。因而，若在x=0 边界处设置条件 L− f x =0 = 0 就相当于使截断边界面处右向 − 行波，即反射波分量等于零，将算子 L 具体表达式代 入上式可得
( ∂ 2 − j k2 − ky )f ∂x