第三章 时域有限差分法

第三章時域有限差分法
此章節將介紹時域有限差分法(Finite Difference Time Domain Method)的發展起源、理論、精確度與穩定條件。

3.1 FDTD發展歷史
FDTD的概念是Yee[18]在1966年所提出，當初只用來解決電磁場在空間中的傳播情形。

1975年Taflove和Brodwin修正Yee的演算法並且推導出FDTD穩定收斂的條件[19][20]。

1977年Holland[21]、Kunz 和Lee[22]應用在解析電磁脈衝(Electromagnetic Pulse, EMP)的問題上。

1981年Mur[23]提出了吸收邊界條件(Absorbing Boundary Condition, ABC)的觀念，當電磁波傳播至模擬邊界時可以減小或消除反射的電磁波，也稱為非反射邊界條件(Non-reflecting Boundary Condition)或輻射邊界條件(Radiating Boundary Condition)，在有限的空間裡模擬電磁波在邊界一去不復返的效果。

1994年Berenger[24]提出一個更有效更完美的吸收邊界，完全匹配層(Perfect Matched Layer, PML)；並與Katz等人[25]將完全匹配層延伸至三維空間上。

3.2 馬克斯威爾方程式
F D T D 的基礎計算是由電磁波的馬克斯威爾方程式 (Maxwell ’s equation)中與時間相關的式子而來的，一個為法拉第定律(Faraday ’s Law )另一個為安培定律(Ampere ’s Law)。

其式子如下: 法拉第定律:
E t
B
⨯-∇=∂∂ (3.1) 安培定律:
J H t
D
-⨯∇=∂∂ (3.2) 假設介質為線性(linear)、等向性(isotropic)、非色散材料(nondispersive materials )，我們可以用簡單的式子來表示D 與E 和B 與H 的關係:
H B μ= (3.3)
E D ε= (3.4) E J σ= (3.5)
其中0μμμr =為導磁係數(magnetic permeability )，0εεε
r =為介電係
數(electrical permittivity )，σ為導電係數(electric conductivity)。

r μ為相對導磁係數(relative permeability );
0μ為在真空中的導磁係數
(free space permeability, 7
104-⨯πH/m);
r ε為相對的介電係數(relative
permittivity);
0ε為真空中的介電係數(free space permittivity,
1210854.8-⨯F/m )。

將(3.3)、 (3.4) 、(3.5)代入(3.1)、(3.2)可以得到:
t
H
E ∂∂-=⨯∇μ (3.6)
t
E
E H ∂∂=-⨯∇εσ (3.7) (3.6)與(3.7)兩個旋度方程式可以推出以直角坐標系統(x, y, z)表示的六個純量方程式，:
⎪⎪⎭⎫
⎝⎛∂∂-∂∂=∂∂z y t y z x E E 1H μ (3.8.1) ⎪⎭
⎫
⎝⎛∂∂-∂∂=∂∂x z t z x y E E 1H μ (3.8.2) ⎪⎪⎭⎫ ⎝⎛∂∂-∂∂=∂∂y x t x y z E E 1H μ (3.8.3) ⎪⎪⎭⎫ ⎝⎛-∂∂-∂∂=∂∂x y z x z y t E H H 1E σε (3.8.4) ⎪⎭⎫ ⎝⎛-∂∂-∂∂=∂∂y z x y
x z t E H H 1E σε (3.8.5) ⎪⎪⎭
⎫ ⎝⎛-∂∂-∂∂=∂∂z x
y z y x t E H H 1E σε (3.8.6)
3.3 FDTD 基本架構與運算
FDTD 的原理就是利用電磁場的交互作用，在時間與空間中以差分法分割成離散的座標來表示。

差分法可分為前差分、中間差分、後差分，如圖3-1。

各種差分法的表示如下: (1) 前向差分法
h
x f h x f x f dt d )
()()(-+=
(3.9)
(2) 中間差分法
h
h x f h x f x f dt d )
()()(--+=
(3.10)
(3) 後向差分法
h
h x f x f x f dt d )
()()(--=
(3.11)
圖3-1 )(a 前向差分法 )(b 中間差分法 )(c 後向差分法
(x f x
(x f x
(x f x
)(a )
(b )
(c
為了方便FDTD 計算電磁場的問題，Yee 將空間細分為許多大小一樣的均勻網格，使時間與空間座標從連續變成離散的，所以每個網格點的座標可以標示成
) , ,( ) , ,(z k y j x i k j i ∆∆∆= (3.12)
其中x ∆、y ∆、z ∆是代表在x 、y 、z 直角坐標系統裡三個方向的間距變化量，即每個網格的三個邊長；而i 、j 、k 為整數。

在此我們先定義一個空間與時間的方程式為
n , ,) , , ,(k j i v t n z k y j x i v =∆∆∆∆ (3.13)
其中t ∆為時間間距(time step)，n 為整數。

圖3-2 FDTD 空間中電場與磁場的配置圖
圖3-2是空間中電場與磁場的配置圖，其中電場與磁場互相交錯並間距為半個網格，每一個電場周圍都被磁場包圍著，同樣，每一個磁場都被電場包圍著；此配置必須符合安培定律與法拉第定律。

在時間中電場與磁場也是互相交錯的，如圖3-3，其計算的方式稱為跳步式(leapfrog)計算，即若要計算當前的電場，必須利用前一時刻所被儲存的磁場與電場來計算;相同的，若要計算當前磁場，必須利用前一時刻所儲存的電場與磁場來計算。

圖3-3 FDTD時間中電場與磁場的配置圖
現在我們將方程式(3.13)分別對空間與時間利用中央差分法來表示，分別可以得到 (1)對空間
x
v v t n z k y j x i x v
n k j i n k j i ∆-=∆∆∆∆∂∂-+,,21,,21) , , ,( (3.14) (2)對時間
t
v v t n z k y j x i t v n k
j i n k j i ∆-=∆∆∆∆∂∂-+2
1,,21,,) , , ,( (3.15) 接下來利用(3.14)、(3.15)兩式將之前以直角坐標系統(x, y, z)表示的六個純量方程式(3.8.1) ~ (3.8.6)以離散的形式來表示。

以x H 方向來說明計算過程，如(3.8.1)式。

從圖3-2我們可以發現x H 的周圍環繞著y E 與
z E ，將此平面放大即為圖3-4。

y
E y
E z
E E x
H 2
1
,,(+k j i )
2
1
,1,++k j i )
,2
1,(k j i +)
1,2
1
,(++k j
i y
z
圖3-4 FDTD 網格放大圖
從圖3-4可以知道y E 的座標為 ⎪⎭⎫ ⎝⎛+k j i ,21, 與 ⎪⎭⎫
⎝⎛++1,21,k j i ，
z E 的座標為⎪⎭⎫ ⎝
⎛
+21,,k j i 與⎪⎭⎫ ⎝⎛++21,1,k j i 。

將(3.8.1)利用(3.14)、(3.15)
兩式展開後得到
y
k j i k j i y k j i k j i t k j i k j i k j i ∆++-+-∆+-++=⎥⎥⎥⎥
⎦⎤⎢⎢⎢⎢⎣⎡∆++-++++-+)1 ,21 ,(E ) ,21 ,(E )21 , ,(E )21 ,1 ,(E )21 ,21 ,(H )21 ,21 ,(H )21 ,21 ,(n
y
n y n z n z 21n x 21n x μ
(3.16)
將上式整理過後為
⎥⎥⎥⎦⎤⎢
⎢⎢⎣⎡∆++-+-∆+-++++=∆++-++-+z k j i k j i y k j i k j i k j i t
k j i k j
i )1 ,21 ,(E ) ,21 ,(E )21 , ,(E )21 ,1 ,(E )21 ,21 ,(1 )
2
1
,21 ,(H )21 ,21 ,(H n
y n y n z n z 21
n x 21
n x μ (3.17)
同理，我們也可推出其他方向分量的式子
⎥⎥
⎥⎦⎤⎢⎢⎢⎣⎡∆++--++-∆+--++-++-=∆++--++--+x k j i k j i z k j i k j i k j i t
k j i k j i
)21 1 ,1(E )21 ,1 ,(E ) ,1 ,21(E )1 ,1 ,21(E )
2 ,1 ,2(1 )
2
1
,1 ,21(H )21 ,1 ,21(H n
z n z n x n x 21
n y 21
n y μ
(3.18)
⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡∆+--++--∆++--++++-=∆++--++-
-+
y k j i k j i x k j i k j i k j i t
k j i k j i )1 , ,21(E )1 ,1 ,21(E )1 ,21 ,1(E )1 ,21 ,(E )1 ,2 ,2(1 )
1 ,21
,21(H )1 ,21 ,21(H n
x n x n y n y 21
n z 21
n z μ
(3.19)
)1 ,1 ,2
1
(E )
1 ,1 ,2
(1
)21, 1 ,21(H )23, 1 ,21(H )1 ,21 ,21(H )1 ,23 ,21(H )
1 ,1 ,2
1(E )1 ,1 ,21(E 21
-n x
21
n y
21n y 21-n z 21-n z 1
n x n x ++-
-++-⨯
⎥⎥
⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∆++--++--
∆++--++-=∆++--++----k j i k j i z k j i k j i y k j i k j i t
k j i k j i σε(3.20)
)
1 ,21 ,(E )
1 ,2
1 ,(1 )1, 21
,21(H )1, 21 ,21(H )21 ,21 ,(H )23 ,21 ,(H )1 ,2
1 ,(E )1 ,21 ,(E 21
-n y 21
n z 21n z 21-n x 21-n x 1
n y
n y ++-++⨯⎥⎥
⎥⎥
⎦
⎤⎢⎢⎢
⎢⎣⎡∆++--+++-∆++-++=∆++-++---k j i k j i x k j i k j i z k j i k j i t
k j i k j i σε (3.21)
)
2
1 ,1 ,(E )
2
1 ,1 ,(1 )21, 21 ,(H )21, 23 ,(H )21 ,1 ,,21(H )21 ,1 ,,21(H )
21 ,1 ,(E )21 ,1 ,(E 21
-n z
21
n x 21n x 21-n y 21-n y 1
n z n z ++-++⨯
⎥⎥
⎥⎥
⎦⎤⎢⎢⎢
⎢⎣⎡∆++-++-∆+++-++-=∆++-++---k j
i k j i y k j i k j i x k j i k j i t
k j i k j i σε (3.22)
3.4 精確度與穩定條件
為了能讓模擬的準確度更精確，網格大小的選擇就變得很重要，切得太大模擬時就容易產生誤差，切得太小會使得電腦的記憶體負擔過大且模擬的運算時間也會拖得太長。

一般來說，網格大小的選擇必須小於所欲模擬的波長的十分之一(λ/10)。

時間間距t ∆的選擇也會影響到模擬的精確度，為了獲得一個較穩定的模擬，在空間與時間的大小其必須遵守Courant 條件:
⎪⎪⎭
⎫ ⎝⎛∆+∆+∆≤
∆2221111
z y x t c (3.23)
其中c 為光速，x ∆、y ∆、z ∆是代表在x 、y 、z 直角坐標系統裡三個方向的間距變化量。

假設∆=∆=∆=∆z y x ，則 (1)一維空間
c
t ∆
≤
∆ (3.23.1) (2) 一維空間
2
c t ∆
≤
∆ (3.23.2) (3) 一維空間
3
c t ∆
≤
∆ 3.23.3)。