时域有限差分法(FDTD算法)的基本原理及仿真

时域有限差分法(FDTD算法)的基本原理及仿真
时域有限差分法（FDTD 算法）
时域有限差分法是1966年K.S.Yee 发表在AP 上的一篇论文建立起来的，后被称为Yee 网格空间离散方式。

这种方法通过将Maxwell 旋度方程转化为有限差分式而直接在时域求解, 通过建立时间离散的递进序列, 在相互交织的网格空间中交替计算电场和磁场。

FDTD 算法的基本思想是把带时间变量的Maxwell 旋度方程转化为差分形式，模拟出电子脉冲和理想导体作用的时域响应。

需要考虑的三点是差分格式、解的稳定性、吸收边界条件。

有限差分通常采用的步骤是：采用一定的网格划分方式离散化场域；对场内的偏微分方程及各种边界条件进行差分离散化处理，建立差分格式，得到差分方程组；结合选定的代数方程组的解法，编制程序，求边值问题的数值解。

1.FDTD 的基本原理
FDTD 方法由Maxwell 旋度方程的微分形式出发，利用二阶精度的中心差分近似，直接将微分运算转换为差分运算，这样达到了在一定体积内和一段时间上对连续电磁场数据的抽样压缩。

Maxwell 方程的旋度方程组为：
E E H σε
+∂∂=⨯∇t H H
E m t
σμ-∂∂-=⨯∇ （1） 在直角坐标系中，（1）式可化为如下六个标量方程：
⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫+∂∂=∂∂-∂∂+∂∂=∂∂-∂∂+∂∂=∂∂-∂∂z z x y y y z x x x y
z E t E y H x H E t E x H z H E t E z H y H σεσεσε，⎪⎪⎪⎪⎭
⎪⎪⎪
⎪⎬⎫
-∂∂-=∂∂-∂∂-∂∂-=∂∂-∂∂-∂∂-=∂∂-∂∂z m z
x y y m y z x x m x y z H t H y E x E H t H x E z E H t H z E y E σμσμσμ （2）
上面的六个偏微分方程是FDTD 算法的基础。

Yee 首先在空间上建立矩形差分网格，在时刻t n ∆时刻，F(x,y,z)可以写成
),,(),,,(),,,(k j i F t n z k y j x i F t z y x F n =∆∆∆∆= （3）
用中心差分取二阶精度： 对空间离散：
()[]
2
),,21(),,21()
,,,(x O x
k j i F k j i F x t z y x F n n x
i x ∆+∆--+≈∂∂∆= ()[]
2
),21,(),21,()
,,,(y O y
k j i F k j i F y t z y x F n n y
j y ∆+∆--+≈∂∂∆= ()[]
2
)21,,()21,,()
,,,(z O z
k j i F k j i F z t z y x F n n z
k z ∆+∆--+≈∂∂∆=
对时间离散：
()[]
2
2121),,(),,()
,,,(t O t
k j i F k j i F t t z y x F n n t
n t ∆+∆-≈∂∂-+∆= （4） Yee 把空间任一网格上的E 和H 的六个分量，如下图放置:
o
y
x
z
Ey
Hz
Ex
Ez
Hx
Ey
Ey
Ez
Ex Hy
Ez
Ex
图1 Yee 氏网格及其电磁场分量分布
在FDTD 中，空间上连续分布的电磁场物理量离散的空间排布如图所示。

由图可见，电场和磁场分量在空间交叉放置，各分量的空间相对位置也适合于Maxwell 方程的差分计算，能够恰当地描述电磁场的传播特性。

同时，电场和磁场在时间上交替抽样，抽样时间间隔相差半个时间步，使Maxwell 旋度方程离散以后构成显式差分方程，从而可以在时间上迭代求解，而不需要进行矩阵求逆运算。

因此，由给定相应电磁问题的初始条件，FDTD 就可以逐步推进地求得以后各个时刻空间电磁场的分布。

根据这一原则可以写出六个差分方程:
]
)
2/1,,2/1()2/1,,2/1(),2/1,2/1(),,2/1([.)
,,2/1(2),,2/1(11.),,2/1()
,,2/1(.)
,,2/1(2),,2/1(1),,2/1(2),,2/1(1),,1(2/12
/12/12/11z
k j i H k j i H y
k j i H k j i H k j i t
k j i k j i t k j i E k j i t k j i k j i t k j i k j i E n y n y n z n z n
x n x ∆++--++
∆-+-++∆+++∆+-++∆++
+∆+-
=++++++εσεεσεσ
其余的也如法可以写出，每个网格点上的个场分两的新值依赖于该点在前一时间步长时刻的值机该点周围的临近点上另一场量在早半个时间步长时的值。

因此任一时刻可一次算出一个点，并行算法可计算出多个点。

通过这些运算可以交替算出电场磁场在各个时间步的值。

（5）
根据上述FDTD 差分方程组可得出计算电磁场的时域推进计算方法，如图2所示。

循环n 次
图2 FDTD 在时域的交叉半步逐步推进计算
2.数值稳定性条件
时间步长t ∆，空间步长x ∆，y ∆，z ∆必须满足一定的关系，否则就使得数值表现不稳定，表现为：随着计算步数的增加，计算场量的数值会无限的增大，这种增大不是由于误差积累造成的，而是由于电磁波的传播关系被破坏造成的。

所以t ∆，x ∆，y ∆，z ∆必须满足一定的关系以保证稳定性。

Taflove 等在1975年对Yee 氏差分格式的稳定性进行了讨论，并导出了对时间步长的限制条件。

数值解是否稳定主要取决于时间步长t ∆与空间步长x ∆、y ∆、z ∆的关系。

对于非均匀媒质构成的计算空间选用如下的稳定性条件：
2
22)1()1()1(1
z
y x v t ∆+∆+∆≤
∆ （6） 若采用均匀立方体网格：s z y x ∆=∆=∆=∆， 3
v s
t ∆≤∆ （7） 而一般取：c
x
t 2∆=
∆，c 为光速。

当x ∆，y ∆，z ∆不相等时，c
z y x t 2)
,,min(∆∆∆=∆ （8）
3.数值色散
FDTD 网格中，会导致数字波模在网格中发生改变，这种改变是由于计算网格本身引起的，而非物理因素，所以必须考虑。

即在FDTD 网格中，电磁波的相速与频率有关，电磁波的相速度随波长、传播方向及变量离散化的情况不同而改变。

色散将导致非物理因素引起的脉冲波形畸变、人为的各向异性和虚假
已知t n t t ∆==01 0时刻空间各处的电磁场初始计算2/12t t t ∆+= 时刻空间各处的磁场值
计算2/21t t t ∆+= 时刻空间各处的电场值
折射等现象。

显然，色散与空间、时间的离散间隔有关，如下式所示：
()()()()()⎪⎭
⎫ ⎝⎛∆∆+⎪⎭⎫ ⎝⎛∆∆+⎪⎪⎭⎫ ⎝⎛∆∆+⎪⎭⎫ ⎝⎛∆∆=⎪⎭⎫ ⎝⎛∆∆2sin 12sin 12sin 12sin 12sin 12222222222z k z z k z y k y x k x t t c z z y x ω （9） 与数值色散关系相对应，在无耗介质中的单色平面波，色散解析关系是： ()
2222
z y x k k k c ++=ω （10）
由式（9）可知，当式（9）中的t ∆、x ∆、y ∆、z ∆均趋于零时，它就趋于式（10）。

也就是说数值色散是由于用近似差分替代连续微分而引起的，而且
在理论上可以减小到任意程度，只要此时时间步长和空间步长都足够小。

为获得理想的色散关系，问题空间分割应按照小于正常网格的原则进行。

一般选取的最大空间步长为20m in m ax λ=∆，min λ为所研究范围内电磁波的最小波长。

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

现在适当选取时间和空间步长，传播方向，可以得到理想情况，如下所示： 3-D 方形网格：（数值稳定的极限状态，可得理想色散关系）
取波沿对角线传播3/k k k k z y x ===，3
,με
δδ=
∆=∆=∆=∆t z y x （11）
2-D 方形网格:也是沿对角线传播2/k k k k z y x ===，2
με
δ=∆t （12）
1-D 网格： με
δ=∆t （13）
4.吸收边界条件
在电磁场的辐射和散射问题中，边界总是开放的，电磁场占据无限大空间，而计算机内存是有限的，所以只能模拟有限空间。

即：时域有限差分网格将在某处被截断。

这要求在网格截断处不能引起波的明显反射，因而对向外传播的波而言，就像在无限大的空间传播一样，一种行之有效的方法是在截断处设置一种吸收边界条件。

使传播到截断出的波被边界吸收而不产生反射。

下面只给出Engquist-Majda 吸收边界条件，采用Mur 差分格式，其总体虚假反射在1%~5%之间。

一维一阶近似情形, x=0边界：)]0()1([)1()0(1
1u u x
t c x t c u u n n
n -∆+∆∆-∆+
=++ （14）
二维二阶近似情形, x=0边界：
)]
1,1(),1(2)1,1()1,0(),0(2)1,0([.)
()(2)()],1(),0([2)]
,0(),1([),1(),0(2211
11-+-++-+-+∆+∆∆∆∆++∆+∆∆+-∆+∆∆-∆+
-=-+-+j W j W j W j W j W j W x t c y x t c j W j W x t c x j W j W x
t c x t c j W j W n n n n n n n n n n n n
（
三维二阶近似情形, x=0边界：
)]1,,1(),,1(2)1,,1()1,,0(),,0(2)1,,0([)
()(2)()]
,1,1(),,1(2),1,1(),1,0(),,0(2),1,0([.)
()(2)()],,1(),,0([2)]
,,0(),,1([),,1(),,0(2
22211
11-+-++-+-+∆+∆∆∆∆+-+-++-+-+∆+∆∆∆∆++∆+∆∆+-∆+∆∆-∆+
-=-+-+k j W k j W k j W k j W k j W k j W x t c z x t c k j W k j W k j W k j W k j W k j W x t c y x t c k j W k j W x t c x k j W k j W x
t c x t c k j W k j W n n n n n n n n n n n n n n n n n n 5.仿真
m 文件见附件，程序表现的是使用二维FDTD 算法对TE 波的仿真。

（1。