ＦＤＴＤ解析法（Matlab 版、２次元）プログラム解説 v2.01
１．概要
ここでは、２次元ＦＤＴＤ法の定式化とプログラミングを１次元からの拡張として説明する。さ
らに、時空間分割に関する制約と、安定な数値解析のための条件について述べる。
２．２次元空間におけるFDTD法の定式化
マクスウェル方程式
G
G
⎧
∂H
⎪⎪∇ × E = − μr μ0
∂t
G
⎨
⎪∇ × HG = ε ε ∂E + σ EG + JG
r 0
⎪⎩
∂t
(1)
において、伝搬方向をｘｚ面内にとり、ｙ軸方向には一様な空間を考えると
∂
= 0 であるから、
∂y
２次元面内の電磁波の伝搬はＴＭ波（ Ex , Ez , H y ）およびＴＥ波（ H x , H z , E y ）の独立した場合
に区別することができる。すなわち、それらの直角座標系の成
分についての方程式は、
ｙ
TM
ＴＭ波
∂H y
⎧ ∂Ex ∂Ez
−
= − μ r μ0
⎪
∂x
∂t
⎪ ∂z
∂E
⎪ ∂H y
= ε r ε 0 x + σ Ex + J x
⎨−
∂t
⎪ ∂z
⎪ ∂H y
∂E
= ε r ε 0 z + σ Ez + J z
⎪
∂t
⎩ ∂x
Hｙ
Ex
Ez
(2)
x
伝搬
z
ｙ
TE
Eｙ
ＴＥ波
∂E
⎧ ∂H x ∂H z
−
= ε rε 0 y + σ Ey + J y
⎪
∂x
∂t
⎪ ∂z
∂H x
⎪ ∂E y
= − μ r μ0
⎨−
∂t
⎪ ∂z
⎪ ∂E y
∂H z
= − μ r μ0
⎪
∂t
⎩ ∂x
Hx
Hz
(3)
z
x
伝搬
となる。まず、２次元 Yee 格子を下記のように定義し、これに基づいて方程式を差分近似する。
ｘｚ平面内の
ＴＭ波
x(i)
Ezn,i +1,k +1/2
i+1
Exn,i +1/2,k
i+1/2
Exn,i +1/2, k +1
H yn,+i 1/2
+1/2, k +1/2
H yn,+i 1/2
+1/2, k −1/2
i
Ezn,i ,k +1/2
i-1/2
o
ｘｚ平面内の
ＴＥ波
H yn,+i 1/2
−1/2, k +1/2
k-1/2
k
k+1/2
k+1
z(k)
x(i)
E yn,i +1,k
i+1
H zn,+i +1/2
1/2, k
i+1/2
H xn,+i ,1/2
k −1/2
E
i
E yn,i ,k +1
H zn,+i −1/2
1/2, k
i-1/2
o
H xn,+i ,1/2
k +1/2
n
y ,i , k
k-1/2
k
k+1/2
k+1
z(k)
差分近似式はそれぞれ、次のようになる。
ＴＭ波
n −1/2
⎧ E xn,i +1/2,k +1 − Exn,i +1/2,k Ezn,i +1,k +1/2 − Ezn,i ,k +1/2
H yn,+i 1/2
+1/2, k +1/2 − H y ,i +1/2, k +1/2
−
= − μ r μ0
⎪
Δz
Δx
Δt
⎪
n +1/2
n +1/2
n +1
n
n +1
⎪⎪ H y ,i +1/2,k +1/2 − H y ,i +1/2,k −1/2
Ex ,i +1/2, k − E x ,i +1/2,k
Ex ,i +1/2,k + Exn,i +1/2,k
+ J xn,+i +1/2
= ε rε 0
+σ
⎨−
1/ 2, k
2
Δz
Δt
⎪
n +1/2
⎪ H yn,+i 1/2
E zn,+i ,1k +1/2 − Ezn,i ,k +1/2
Ezn,+i ,1k +1/2 + Ezn,i ,k +1/2
+1/2, k +1/2 − H y ,i −1/2, k +1/2
+ J zn,+i ,1/2
= ε rε 0
+σ
⎪
k +1/2
2
Δx
Δt
⎪⎩
(4)
ＴＥ波
n +1/2
n +1/2
⎧ H xn,+i ,1/2
E yn,+i 1,k − E yn,i ,k
E yn,+i 1,k + E yn,i ,k
H zn,+i +1/2
k +1/2 − H x ,i , k −1/2
1/2, k − H z ,i −1/2, k
+ J yn,+i 1/2
−
= ε rε 0
+σ
⎪
,k
2
Δ
z
Δ
x
Δ
t
⎪
n +1/2
⎪⎪ E yn,i ,k +1 − E yn,i ,k
H xn,+i ,1/2
k +1/2 − H x ,i , k +1/2
= − μ r μ0
⎨−
Δz
Δt
⎪
n
n
n +1/2
⎪ E y ,i +1,k − E y ,i ,k
H z ,i +1/2,k − H zn,+i +1/2
1/ 2, k
= − μ r μ0
⎪
Δx
Δt
⎪⎩
(5)
ここで Δx および Δz は Yee セルの大きさ、すなわち空間ステップであり、Δt は時間ステップ（繰
り返し計算における微小時間間隔）を表す。i, k はそれぞれ x, z を離散化した時の空間の位置を
表す指標（空間インデックス）、すなわち E や H のサンプル点の位置を表し、ｎは現時点の時間
が何番目であるか（すなわち、時間インデックス）を表す。これらの差分方程式を最新の時間の
n
電界 E や磁界 H
n +1/2
について計算する。すなわち、
ＴＭ波
n
n
n
n
⎧ n +1/2
Δt ⎧⎪ Ex ,i +1/2,k +1 − Ex ,i +1/2,k Ez ,i +1,k +1/2 − Ez ,i ,k +1/2 ⎫⎪
n −1/2
−
⎪ H y ,i +1/2,k +1/2 = H y ,i +1/2,k +1/2 −
⎨
⎬
μr μ0 ⎩⎪
Δz
Δx
⎪
⎭⎪
⎪
n +1/2
n +1/2
⎛ H y ,i +1/2,k +1/2 − H y ,i +1/2,k −1/2
⎞
2ε r ε 0 − σΔt n
2Δt
⎪ n +1
Ex ,i +1/2,k +
− J xn,+i +1/2
⎜⎜ −
⎨ Ex ,i +1/2,k =
1/2, k ⎟
⎟
2ε r ε 0 + σΔt
2ε r ε 0 + σΔt ⎝
Δz
⎪
⎠
⎪
n +1/2
n +1/2
⎛ H y ,i +1/2,k +1/2 − H y ,i −1/2, k +1/2
⎞
2ε ε − σΔt n
2Δt
⎪ E n +1
Ez ,i ,k +1/2 +
= r 0
− J zn,+i ,1/2
⎜
⎟
,
,
1/2
1/2
z
i
k
k
+
+
⎪
⎟
2ε r ε 0 + σΔt
2ε r ε 0 + σΔt ⎜⎝
Δx
⎠
⎩
(6)
ＴＥ波
n +1/2
n +1/2
⎧ n +1 2ε ε − σΔt n
⎛ H xn,+i ,1/2
⎞
H zn,+i +1/2
2Δt
k +1/2 − H x ,i , k −1/2
1/2, k − H z ,i −1/2, k
r 0
E y ,i , k +
−
− J xn,+i ,1/2
⎪ E y ,i , k =
⎜⎜
k ⎟
⎟
2ε r ε 0 + σΔt
2ε r ε 0 + σΔt ⎝
Δz
Δx
⎪
⎠
⎪
n
n
Δt ⎧⎪ E y ,i ,k +1 − E y ,i ,k ⎪⎫
⎪ n +1/2
n −1/2
(7)
⎨ H x ,i ,k +1/2 = H x ,i ,k +1/2 −
⎨−
⎬
μr μ0 ⎪⎩
Δz
⎪⎭
⎪
⎪
n
n
⎪ H n +1/2 = H n −1/2 − Δt ⎧⎪ E y ,i +1,k − E y ,i ,k ⎫⎪
⎨
⎬
z ,i +1/2, k
⎪ z ,i +1/2,k
μr μ0 ⎩⎪
Δx
⎭⎪
⎩
となる。
実際の計算においては、係数をまとめて、
⎧
2ε r ε 0 − σΔt
⎪C0 =
2ε r ε 0 + σΔt
⎪
⎪
2Δt
⎨Ce =
2ε r ε 0 + σΔt
⎪
⎪
Δt
⎪Cm =
μr μ0
⎩
(8)
としておき、それぞれの方程式はこれらの係数を用いて
ＴＭ波
n
n
n
n
⎧ n +1/2
⎪⎧ Ex ,i +1/2,k +1 − Ex ,i +1/2,k Ez ,i +1,k +1/2 − Ez ,i ,k +1/2 ⎪⎫
−
−
C
⎪ H y ,i +1/2,k +1/2 = H yn,−i1/2
⎬
+1/2, k +1/2
m ⎨
Δz
Δx
⎪
⎩⎪
⎭⎪
⎪
n +1/2
⎞
⎛ H yn,+i 1/2
⎪ n +1
+1/2, k +1/2 − H y ,i +1/2, k −1/2
n
− St J xn,+i +1/2
⎨ Ex ,i +1/2,k = C0 Ex ,i +1/2,k + Ce ⎜⎜ −
1/2, k ⎟
⎟
Δz
⎪
⎠
⎝
⎪
n +1/2
⎛ H yn,+i 1/2
⎞
+1/2, k +1/2 − H y ,i −1/2, k +1/2
n
⎪ E n +1
= C0 Ez ,i ,k +1/2 + Ce ⎜
− St J zn,+i ,1/2
+
+
,
,
1/2
1/2
z
i
k
k
⎪
⎜
⎟⎟
Δ
x
⎝
⎠
⎩
ＴＥ波
(9)
n +1/2
n +1/2
⎧ n +1
⎛ H xn,+i ,1/2
⎞
H zn,+i +1/2
k +1/2 − H x ,i , k −1/2
1/2, k − H z ,i −1/2, k
n
−
− St J xn,+i ,1/2
⎪ E y ,i ,k = C0 E y ,i ,k + Ce ⎜⎜
k ⎟
⎟
Δz
Δx
⎪
⎝
⎠
⎪
n
n
⎧⎪ E y ,i ,k +1 − E y ,i ,k ⎫⎪
⎪ n +1/2
n −1/2
⎨ H x ,i ,k +1/2 = H x ,i ,k +1/2 + Cm ⎨
⎬
Δz
⎪⎩
⎪⎭
⎪
⎪
n
n
⎪ H n +1/2 = H n −1/2 − C ⎧⎪ E y ,i +1,k − E y ,i ,k ⎫⎪
⎬
z ,i +1/2, k
m ⎨
⎪ z ,i +1/2,k
Δx
⎩⎪
⎭⎪
⎩
(10)
として繰り返し計算する。ここで S t は励起パルスの時間変化を設定するための係数である。
３．時間および空間分割に関する条件
３－ａ．空間分割と数値分散誤差
電磁界は本来連続な物理量であるがＦＤＴＤ法などの数値解析においては離散的な数値の配列
によってこれを表現している。この空間の変化を精度良く表現するためには、まず標本定理によ
る最も粗い分割に対する制約がある。すなわち、電磁波の波長を λ 、空間分割を Δz とすると、標
本定理より
Δz ≤
λ
(11)
2
でなければ連続な波を離散的な数値により表現できない。
次に、ＦＤＴＤ法において、十分な精度で解析を行うためには一般に Δz ≤ λ 10 ないし λ 20 程
度の空間分割が必要となるが、このことを以下に説明する。簡単のため、１次元のマクスウェル
方程式を考える。すなわち、
アンペールの法則：
∂Ex
1 ∂Hy
=−
ε ∂z
∂t
(12)
∂Hy 1 ∂Ex
=
μ ∂z
∂t
(13)
ファラデーの法則： −
において、差分近似は空間分割 Δz および時間分割 Δt を用いて、
n +1
アンペールの法則：
Ex i − Ex i
Δt
n
n +1/2
n +1/2
1 Hy i +1/2 − Hy i −1/2
=−
Δz
ε
(14)
等となる。１次元のＦＤＴＤ解析によって得られるＴＥＭ波（平面波）は、
電界： Ex i = Eo exp{ j[ω n Δt − k i Δz ]}
n
(15)
n+
磁界： Hy i + 1 2 = H o exp{ j[ω ( n + 1 ) Δt − k (i + 1 ) Δz ]}
1
2
2
と表される。ただし、 k =
2π
λ
2
(16)
は波数である。いま、(15)および(16)式を(14)の差分式に代入する
と、
exp{ j[ω (n + 1)Δt − kiΔz ]} − exp{ j[ω nΔt − kiΔz ]}
=
Δt
(17)
1
1
1
1
H o exp{ j[ω (n + 2 )Δt − k (i + 2 )Δz ]} − exp{ j[ω (n + 2 )Δt − k (i − 2 )Δz ]}
−
ε
Δz
Eo
となり、これをオイラーの公式によりまとめると、(12)と(13)式に相当する離散表現の式は
アンペールの法則：
Eo
H
ω Δt
k Δz
sin(
) = − o sin(
)
2
2
Δt
ε Δz
ファラデーの法則： −
Ho
E
ω Δt
k Δz
sin(
) = o sin(
)
2
2
Δt
μ Δz
(18)
(19)
となる。これらの両辺をそれぞれ乗算し、 εμ = 1 / c の関係を利用すると、
2
ωΔt ⎫ ⎧ 1
k Δz ⎫
⎧ 1
sin(
) ⎬ = ⎨ sin(
)⎬
⎨
2 ⎭ ⎩ Δz
2 ⎭
⎩ cΔt
2
2
(20)
が得られる。すなわち、この式は時空間分割 Δt と Δz が与えられた時の角周波数 ω と波数 k の関
係を表し、 ω = ck からの差が離散的な数値解析における誤差（数値分散誤差）を表す。概ね、
Δz ≤ λ 20 において数値分散誤差は１％程度となる。
π
規格化週波数
ω = ck
⎛ ω Δ t = k Δz ⎞
⎜
⎜⎜
⎝
ω Δt
Δt
c=
Δz
⎟
⎟⎟
⎠
数値分散誤差
k Δz ⎫
ωΔt ⎫ ⎧ 1
⎧ 1
sin(
) ⎬ = ⎨ sin(
)⎬
⎨
2 ⎭ ⎩ Δz
2 ⎭
⎩ cΔt
2
]
d
a
r
[
４分割
８分割
0
2
２分割/波長
π
k Δz
規格化波数 [rad]
分割： Δz
≤
λ
2
に相当
以上をまとめると、各空間次元数により、数値分散特性は次の様に与えられる。
１次元ＦＤＴＤ法： ω = ck （ただし、 cΔt = Δz を満たすとき）
ωΔt ⎫ ⎧ 1
k Δx ⎫
⎧ 1
sin(
) ⎬ = ⎨ sin(
) ⎬ （任意の Δt と Δx のとき）
１次元ＦＤＴＤ法： ⎨
2 ⎭ ⎩ Δx
2 ⎭
⎩ cΔt
(21)
2
2
k y Δy ⎫
k x Δx ⎫ ⎧ 1
ωΔt ⎫ ⎧ 1
⎧ 1
２次元ＦＤＴＤ法： ⎨
sin(
) ⎬ = ⎨ sin(
) ⎬ + ⎨ sin(
)⎬
2 ⎭
2 ⎭ ⎩ Δx
2 ⎭ ⎩ Δy
⎩ cΔt
(22)
2
2
2
2
2
2
k y Δy ⎫ ⎧ 1
k x Δx ⎫ ⎧ 1
k z Δz ⎫
ωΔt ⎫ ⎧ 1
⎧ 1
３次元ＦＤＴＤ法： ⎨
sin(
) ⎬ = ⎨ sin(
) ⎬ + ⎨ sin(
) ⎬ + ⎨ sin(
)⎬
2 ⎭
2 ⎭ ⎩ Δx
2 ⎭ ⎩ Δy
2 ⎭ ⎩ Δz
⎩ c Δt
2
(23)
３－ｂ．時間分割とCoulantの安定条件
必要な数値解析精度から空間分割が決定されると、次は、安定に解析できる時間分割を決定し
なければならない。この条件は Coulant の安定条件と呼ばれており次のように求めることが出来
る。２次元の場合を例に取ると、その数値分散式は(22)である。ここで、最も粗い分割である(11)
の条件では
sin(
k x Δx
) =1
2
および、 sin(
k y Δy
2
) =1
(24)
であるから、(22)式から
ωΔt ⎫ ⎧ 1 ⎫ ⎧ 1 ⎫
⎧ 1
sin(
)⎬ = ⎨ ⎬ + ⎨ ⎬
⎨
2 ⎭ ⎩ Δx ⎭ ⎩ Δy ⎭
⎩ cΔt
2
2
2
(25)
が得られる。ここで周波数 ω が実数であるためには、
2
ωΔt ⎫ ⎧ cΔt ⎫ ⎧ cΔt ⎫
⎧
0 ≤ ⎨sin(
)⎬ = ⎨
⎬ +⎨
⎬ ≤1
2 ⎭ ⎩ Δx ⎭ ⎩ Δy ⎭
⎩
2
2
(26)
でなければならないから、右辺の関係において Δt で整理すると、
2
2
1 ⎡⎛ 1 ⎞ ⎛ 1 ⎞ ⎤
２次元ＦＤＴＤ法： Δt ≤ ⎢⎜
⎟ ⎥
⎟ +⎜
c ⎢ ⎝ Δ x ⎠ ⎝ Δy ⎠ ⎥
⎣
⎦
−1
2
(27)
が得られる。これが２次元ＦＤＴＤ解析における時間分割の選び方を規定し、これを満たすよう
な Δt を用いると安定に解析が出来る。同様に
1 ⎡⎛ 1 ⎞
１次元ＦＤＴＤ法： Δt ≤ ⎢⎜
⎟
c ⎢⎣⎝ Δx ⎠
2
⎤
⎥
⎥⎦
−1
2
=
Δx
c
1 ⎡⎛ 1 ⎞ ⎛ 1 ⎞ ⎛ 1 ⎞
３次元ＦＤＴＤ法： Δt ≤ ⎢⎜
⎟ +⎜ ⎟
⎟ +⎜
c ⎢ ⎝ Δ x ⎠ ⎝ Δ y ⎠ ⎝ Δz ⎠
⎣
2
2
(28)
2
⎤
⎥
⎥⎦
−1
2
(29)
が得られる。特に各方向の分割が等しいとき、
２次元（ Δx = Δy ）： Δt ≤
Δx
c 2
３次元（ Δx = Δy = Δz ）： Δt ≤
(30)
Δx
c 3
(31)
となり、これらの結果から、次元数が大きいと時間分割は小さくしなければならないことがわか
る。
© Copyright 2024 ExpyDoc
