Transcript 第2讲Maxwell方程Yee算法

第 2 讲
Maxwell方程Yee算法
• 本讲介绍K.S. Yee提出的FDTD算法，它是电磁场
FDTD分析的基础。
• Yee的独特之处是在空间为每一个电场和磁场分量
的空间取样选择一种特殊的网格—称之为Yee网格，
在时间上，采用了蛙跳算法，使得利用一阶导数
的二阶中心差分近似从Maxwell方程获得的FDTD公
式，既满足Maxwell方程的微分形式又满足其积分
形式。
• 因此，Yee的FDTD算法非常稳固，具有很广的应用
领域。
2.1 一维Maxwell方程的Yee算法（1）
一维Maxwell方程
E x
t
 
1
 0 r
H
H
y
z
y
t
 
E x
1
 0  r z
利用一阶导数的二阶中心差分近似，上面的方程变为
n 1
Ex
(k )  E x (k )
n
t
H
n 1 / 2
y
( k  12 )  H
t
 
n 1 / 2
y
H
1
n 1 / 2
y
( k  12 )  H
 0 r (k )
( k  12 )
 
n 1 / 2
y
( k  12 )
z
1
E x ( k  1)  E x ( k )
 0  r ( k  12 )
z
n
n
2.1 一维Maxwell方程的Yee算法（2）
采用归一化磁场
~
H 
0
0
H
使得电场与归一化磁场有相同的数量级，于是可以
得到FDTD迭代公式为
1
1
c t
~ n 1 / 2
~ n 1 / 2
n
n
Hy
(k  )  H y
(k  ) 
E x ( k  1)  E x ( k )
2
2
 z  r ( k  12 )

n 1
Ex
(k )  E x (k ) 
式中，c 
n
1
 0 0
c t
1
1 
~ n 1 / 2
 ~ n 1 / 2
H
(
k

)

H
(
k

)
y
y

 z r (k ) 
2
2 
为自由空间中的光速。

2.1 一维Maxwell方程的Yee算法（3）
用计算机语言表示的FDTD公式
Hy [ k ]  Hy [ k ]  ca [ k ] *  Ex [ k  1]  Ex [ k ] 
Ex [ k ]  Ex [ k ]  cb [ k ] *  Hy [ k ]  Hy [ k  1] 
式中，时间变量已隐含在迭代公式中，以及
1
~
Hy [ k ]  H y ( k  ); Ex [ k ]  E x [ k ];
2
c t
c t
ca [ k ] 
; cb [ k ] 
1
 z r (k  2 )
 z r (k )
只要给定了所有空间点上电/磁场的初值，就可以一步
一步地求出任意时刻所有空间点上的电/磁场值。
2.1 一维Maxwell方程的Yee算法（3）
Ex
n 2
H
y
n  3/2
Ex
n 1
H
y
n 1/ 2
Ex
n0
0
1
2
3
k
电场与磁场分量的空间-时间分布图
2.1 一维Maxwell方程的Yee算法（4）
Main loop in 1D FDTD C-program:
for (k=0;k<=kmax;k++)
{ Hy[k]=0;
Ex[k]=0;}
for (n=1;n<=nmax;n++)
{
Ex(0)=Source(n);
for (k=0;k<kmax;k++)
{ Hy[k]=Hy[k]-ca[k]*(Ex[k+1]-Ex[k]); }
for (k=1;k<kmax;k++)
{ Ex[k]=Ex[k]-cb[k]*(Hy[k]-Hy[k-1]); }
Ex(kmax)=Boundary;
}
2.2 三维Maxwell方程的Yee算法（1）
考虑非时变、线性、各向同性媒质填充的无源区域，
Maxwell旋度方程为
H
x
t



H
  E   H  
t
H
y
t
H
t
z

E z
1  E y


 H x 

  z
y

E x
1  E z



 H y 

  x
z

E y

1  E x


 H z 

  y
x

2.2 三维Maxwell方程的Yee算法（2）
以及
E x


  H  E  

E
t
t
E y
H y

1  H z
 

 E x 
  y
z

H z
1  H x




E
y 
t
   z
x


H x
E z
1  H y
 

 E z 
t
  x
y


2.2 三维Maxwell方程的Yee算法（3）
• Yee 首先将空间
按立方体分割，
电磁场的六个分
量在空间的取样
点分别放在立方
体的边沿和表面
中心点上，电场
与磁场分量在任
何方向始终相差
半个网格步长。
Hz
Ex
Ey
Ex
Ez
Ey
Ez
Hy
Ez
Hx
Ey
z
( i, j, k )
y
x
Ex
2.2 三维Maxwell方程的Yee算法（4）
在时间上，
Yee 把电场
分量与磁场
分量也差半
个步长取样。
Ez
Ez
Ez
Ez
t=2t
Hy
Hy
Hy
t=1.5t
Ez
Ez
Ez
Ez
t=t
Hy
Hy
Hy
t=0.5t
Ez
0
Ez
x
Ez
2x
Ez
3x
x
t=0
2.2 三维Maxwell方程的Yee算法（5）
于是，利用
一阶导数的
二阶中心差
分近似就可
以导出旋度
方程的FDTD
公式。如：
n
H
H
n
t
E y
z
E z
y
2
x

x
1
1
i , j  , k 1
2
n
Ey

1
i , j  , k 1
2
Ez

1
1
i , j  ,k 
2
2
2
x
1
i , j  ,k
2
 Ey
n
i , j 1, k 
1
 Ez
2
y

2


2


2

 O  t 
n
1
i , j  ,k
2
z
1
1
i , j  ,k 
2
2
n
1
t
1
1
i , j  ,k 
2
2
n
H
n
 O  z 
n
i , j ,k 
1
2
 O  y 
2.2 三维Maxwell方程的Yee算法（6）
将上述公式代入相应的方程，得
n
H
1
 H
2
x
i, j 
1
,k 
2
1
n
1
2
x
i, j 
2
1
,k 
2
1
2
t
n

1

i, j 
Ez

Ey
1
i, j 
, k 1
2
n
2
i, j ,k 
1
 Ez
2
1
,k
2


 O  z 
2

n
i , j  1, k 
1
2
y
i, j 
z
1
2
n
 Ey
2
{
,k 
1

 O  t 

 O  y 
2
 
i, j 
1
2
,k 
1
2
H
n
x
i, j 
1
2
,k 
1
2
}
2.2 三维Maxwell方程的Yee算法（7）
采用时间平均近似
n
H
H
n
1
1
x i , j  ,k 
2
2

1
2
x
1
1
i , j  ,k 
2
2
H
n
1
2
x
1
1
i , j  ,k 
2
2
2

 O  t 
最后，忽略高次项，得
H
x
i ,
j , k , n  1  D aHx i , j , k H
 D bHx
x
i ,
j, k , n 
 E y i , j , k  1, n   E y i , j , k , n  





z
i , j , k 
E z i , j  1, k , n   E z i , j , k , n  




y


2

2.2 三维Maxwell方程的Yee算法（8）
同理，可以得到其他2个磁场分量的FDTD方程
H y i , j , k , n  1  D aHy i , j , k H y i , j , k , n 
 E z i  1, j , k , n   E z i , j , k , n 




x
 D bHy i , j , k 

E x i , j , k  1, n   E x i , j , k , n 




z
H z i , j , k , n  1  D aHz i , j , k H z i , j , k , n 
D aHw
 E x i , j  1, k , n   E x i , j , k , n  




y
 D bHz i , j , k 

E y i  1, j , k , n   E y i , j , k , n 




x


2   t
2t

D bHw 
2     t H 所在空间位置
2     t H 所在空间位置
w
w
2.2 三维Maxwell方程的Yee算法（9）
上面公式之间有明显的规律，便于记忆，如：
• 系数Da和Db在空间的位置就是方程左边项的场分量
的空间位置；
• 右边第一项的场分量与左边的相同，但为n时间步，
而左边场分量的时间步为n+1；
• 右边第二、第三项的场分量与左边的相反（电场与磁
场）三者的坐标分量满足循环关系：x-y-z-x;
• 右边第二、第三项为空间差分形式。右边第二项分子
上场量的坐标分量与分母上空间步长的坐标分量也满
足x-y-z-x的循环关系。而右边第三项分子上场量的坐
标分量和空间步长坐标分量与第二项恰好对调。
• 第三项符号为负。
2.2 三维Maxwell方程的Yee算法（10）



利用对偶原理： H  E ， E   H ，    ,    ,
Da  C a , Db  Cb
，并注意到E与H在时间上差半个步
长，可以直接从磁场FDTD公式得到电场的FDTD公式。
如：
E x i , j , k , n  1  C aEx i , j , k E x i , j , k , n   C bEx i , j , k 
 H y i , j , k  1, n  1  H y i , j , k , n  1 





z

H z i , j  1, k , n  1  H z i , j , k , n  1 




y


C aEw 
2    t
2    t
C bEw 
E w 所在空间位置
2t
2    t
E w 所在空间位置
2.2 三维Maxwell方程的Yee算法（11）
• 给出n=0时刻电磁场的初值和媒质参数；
• 由磁场FDTD公式，根据n时间步的电场值和磁场值
求得n+1时间步空间所有点的磁场分量；
• 由电场FDTD公式，根据n时间步的电场值和n+1时间
步的磁场值求得n+1时间步空间所有点的电场分量；
• 如此迭代，可获得任何时刻空间所有点的电磁场值。
每一过程常称为蛙跳法(leapfrog)。
• 对于连续变化的媒质，FDTD法需储存的量有n时间
步和n+1时间步的六个场分量，六个D参数和六个C
参数，所以总储存量近似为24N, N为空间网格数。
但对于均匀媒质，D参数和C参数为常数，故总储
存量减少为12N。
2.2 三维Maxwell方程的Yee算法（12）
媒质参数赋值
在所有空间点给电磁场分量赋初值
求所有空间离散点上n+1时间步的磁场
求所有空间离散点上n+1时间步的电场
n=n+1
No
Yes
n>nmax
结束
2.3 以积分形式的Faraday和 Ampere 定理
解释Yee算法（1）
• 上面介绍的FDTD算法是从点的观点对Maxwell方程
微分形式中的两个旋度方程直接进行导数二阶中心
差分近似得到的。
• 这种观点对理解FDTD如何模拟波在媒质中的传播是
有用的。但是，当模拟细几何结构如导线、槽和曲
面时，点的观点对于指导为了获得适当解需要作怎
样的算法修正却帮助甚少。
• 为了解决这一问题，我们从积分形式的Ampere和
Faraday定理出发来解释Yee算法。
2.3 以积分形式的Faraday和 Ampere 定理
解释Yee算法（2）
Ey
• 仅讨论自由空间的情
况。考虑右图中实线
网格的y-z平面上包围
面积S1的矩形围线C1。
• 沿C1应用Faraday定理:

t
Ez
C1
Hx
Hx
Ez
Hz
S1
Ey
( i , j, k )
S2

 

  0 H  ds    E  dl
s1
Hz
C2
c1
( i+1/2, j+1/2, k-1/2 )
Hx
2.3 以积分形式的Faraday和 Ampere 定理
解释Yee算法（3）
假设在围线任一边上和围线所围的面积上相应的场
值不变，时间导数采用中心差分近似，则Faraday定
理近似为
H
 0  y z

 Ey


n  12
x
1
1
i , j  ,k 
2
2
n  12
x
1
1
i , j  ,k 
2
2
t
n
i, j
H
1
2
, k 1
 Ey
n
i, j
1
2
,k


y   E

 z


n
i , j 1, k 
1
2
整理后便可得到自由空间中FDTD公式。
 Ez
n
i , j ,k 
1
2

z


2.3 以积分形式的Faraday和 Ampere 定理
解释Yee算法（4）
• 以相同的方式，把Ampere 定理

t
 
  0 E  ds 


 H  dl
s2
c2
应用于图中虚线网格的x-z平面上包围面积S2 的矩形围
线C2 ，并作类似的假设，也可以得到相应的FDTD公式。
• 所以，FDTD公式既是微分形式的Maxwell旋度方程的
中 心 差 分 近 似 ， 也 自 然 满 足 积 分 形 式 的 Ampere 和
Faraday定律。
2.4 Yee算法的无散性（1）
• 对于无源区域，满足Maxwell两个旋度方程的场
也一定满足Maxwell的两个散度方程或它们的积
分形式

  D  0,

s


D  ds  0,
B  0

 
B  ds  0
s
• 下面证明对于从旋度方程近似而来的FDTD 公式也
满足两个散度方程。
2.4 Yee算法的无散性（2）
在自由空间的一个Yee网格上考虑 

t

Yee cell
 
 
B  ds  0  H
t 

 
B  ds
s
x i, j
1
2
,k 
1
H
 0 
 H y
t 
 
 0 Hz
t 
x i 1, j 
2
1
1
i  , j 1, k 
2
2
1
1
i  , j  , k 1
2
2
1
2
H
H
y
,k 
1
2
1
1
i  , j ,k 
2
2
z i  1 , j  1 ,k
2
2
，有

 y z


  x  z


 x z

利用磁场分量的FDTD公式，把与磁场分量时间导数相关
的电场空间有限差分代入上式中各项，可得

t

 
B  d s  (Term 1)  y  z  (Term 2 )  x  z  (Term 3 )  x  y  0
Yee cell
2.4 Yee算法的无散性（3）
于是



B  d s  const
Yee cell
设初始时磁场为零，则



B  d s  const  0
Yee cell
所以，对于无源区域，FDTD公式满足磁场Gauss定理，
即对于磁场是无散的。
同理可以证明，对于电场，FDTD公式也满足Gauss定
理，即电场也是无散的。
结 论 2 （1）
本讲介绍了求解矢量Maxwell方程的FDTD Yee算法，归纳
起来，Yee算法的主要特点有：
• Yee算法采用耦合的Maxwell旋度方程，同时在时间和
空间求解电场和磁场，而不是采用波动方程只求解电
场或磁场。
– 同时使用E和H信息比只使用其中一个的优点是获得的解
更稳固（robust），即算法可以适用非常广泛的电磁波
物理结构，并且电场和磁场的特性可以用更直接的方式
模拟。
– 如果同时使用电场和磁场，每一种场的独立特性，如边
沿和角处切向磁场的奇异性、细线附近磁场的奇异性以
及靠近点、边沿和细导线处径向电场的奇异性就能够独
立地模拟。
结 论 2 （2）
• Yee网格在三维空间这样安排E和H分量，使得每一个E
或H分量由四个H或E循环的分量所环绕。
– 这提供了一幅三维空间中由相互交链的Faraday定理和
Ampere定理围线阵列构成的优美而简单的图画。保证了
Yee算法同时模拟了Maxwell方程点意义上的微分形式和
宏观的积分形式。后者对于处理边界条件和奇异性是极
其有用的。
– 旋度算子中空间导数的差分公式是二阶精度的中心差分。
– 如果不同材料的交界面平行于Yee网格的一个坐标轴，在
交界面上切向E和H的连续性自然保持。
– 在Yee算法隐含地执行了两个高斯定律。所以，同时保证
了无源区域中电磁场的无散性。
结 论 2（3）
• Yee算法以蛙跳算法在时间上安排E和H分量。在某一
时刻，使用前一时刻的E数据计算所有H分量。然后，
再使用刚计算的H数据计算所有的E分量。如此循环，
直至完成时间步进过程。
– 蛙跳时间步进过程是全显式的，所以完全避免了因求
解联立方程和矩阵求逆所带来的问题。
– 旋度方程中时间导数的差分公式是二阶精度的中心差
分。
– 时间步进算法是无数值损耗的，即在网格中传播的数
值波并不产生寄生衰减。
习
题 2
2.1 证明对于电场，Yee算法也满足Gauss定理，即对于
电场Yee算法也是无散的。
2.2 试推导二维TE模（Ez=0）和TM模（Hz=0）的FDTD Yee
算法。
2.3 编制二维TM模Yee 算法的程序。假设模拟区域为自
由空间单位正方形，时间步为  t   x /( c 2 ) ，x方向
与y方向步长相等。模拟区域的边界为理想电导体。
设在区域的中心电场分量Ez随时间按高斯或正弦变
化。在外向波到达区域边界之前和之后的一些时刻
求区域内的电场和磁场分布。对于正弦激励的情况，
确定外向波振幅离开源点随径向距离的衰减特性，
并与二维解析Green函数相比较。