1. No category

ＦＤＴＤ解析法（Matlab 版、２次元ＰＭＬ）プログラム解説 v2.11
１．概要
ＦＤＴＤ解析における吸収境界である完全整合層（Perfectly Matched Layer, PML）の定式化
とプログラミングを２次元ＴＥ波について説明する。
２．２次元ＴＥ波におけるＰＭＬの定式化
２次元ＴＥ波についてのマクスウェル方程式
∂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
TE 波
ｙ
Eｙ
(1)
Hx
Hz
z
x
伝搬方向
に対応するＰＭＬの定式化を行う。
まず、３次元ＰＭＬのためのマクスウェル方程式が、異方性ＰＭＬのテンソル s を用いて次式
で与えられていることを前提とする。（この原理を理解するためには、Berenger の電磁界分離に
よる PML 定式化等を参照のこと。）すなわち、
Ampere の法則： ∇ × H = ε r ε 0 s
∂E
+σ E + J
∂t
Faraday の法則： ∇ × E = − μ r μ0 s
∂H
∂t
(2a)
(2b)
ただし、テンソル
⎡ s y sz
⎢
⎢ sx
⎢
s =ε = μ = ⎢ 0
⎢
⎢
⎢ 0
⎢⎣
0
sz sx
sy
0
⎤
0 ⎥
⎥
⎥
0 ⎥
⎥
⎥
sx s y ⎥
sz ⎥⎦
(3)
であり、
σξ
σξ*
sξ = κξ +
= κξ +
, ξ = x, y , z
jωε 0
jωμ0
(4)
は、異方性の損失を表す。また、周波数分散がゼロの条件を満たすため
σ σ*
=
ε 0 μ0
(5)
としている。さらに、ＰＭＬ媒質の誘電体としての構成方程式を
⎧
s
⎪ Dx = ε 0ε r z Ex
sx
⎪
⎪⎪
sx
⎨ D y = ε 0ε r E y
sy
⎪
⎪
s
⎪ Dz = ε 0ε r y Ez
sz
⎪⎩
(6)
および、磁性体としての構成方程式を
⎧
s
⎪ Bx = μ0 μr z H x
sx
⎪
⎪⎪
sx
⎨ B y = μ0 μ r H y
sy
⎪
⎪
s
⎪ Bz = μ0 μr y H z
sz
⎪⎩
(7)
と置く。以上の関係から３次元空間における定式化を行い、最終的にそれらの方程式から２次元
ＴＥ波の方程式 (1) すなわち、電磁界成分（ E y , H x , H z ）に相当する方程式を抽出することとす
る。２次元ＴＭ波の方程式も同様に求めることができるので各自で試してみるとよい。
Ampere の法則(2a)を直角座標系において書き下すと、
⎡ s y sz ⎤
⎡ ∂H z ∂H y ⎤
Ex ⎥
⎢
−
⎢
⎥
s
∂
∂
y
z
⎢
⎥ J
x
⎢
⎥
⎡ x⎤
⎢
⎢ ∂H x ∂H z ⎥
sz sx ⎥ ⎢ ⎥
−
Ey ⎥ + ⎢ J y ⎥
⎢
⎥ = jωε ⎢
∂
∂
z
x
s
⎢
⎥
y
⎢
⎥
⎢
⎥ ⎢⎣ J z ⎥⎦
∂
H
⎢ y ∂H x ⎥
s
s
x
y
⎢
⎢ ∂x − ∂y ⎥
E ⎥
⎢⎣ sz z ⎥⎦
⎣
⎦
⎡ s y Dx ⎤ ⎡ J x ⎤
= jωε ⎢⎢ sz Dy ⎥⎥ + ⎢⎢ J y ⎥⎥
⎢⎣ sx Dz ⎥⎦ ⎢⎣ J z ⎥⎦
⎡ k y Dx ⎤
1
= jω ⎢⎢ k z Dy ⎥⎥ +
ε
⎢⎣ k x Dz ⎥⎦ 0
(8)
⎡σ y Dx ⎤ ⎡ J x ⎤
⎢σ D ⎥ + ⎢ J ⎥
⎢ z y⎥ ⎢ y⎥
⎢⎣σ x Dz ⎥⎦ ⎢⎣ J z ⎥⎦
となるが、この式をフーリエ逆変換によって時間に関する微分方程式に変換する。時間微分演算
子 ∂ ∂t はフーリエ変換により複素角周波数 jω となることから、(8)式は
⎡ ∂H z ∂H y ⎤
−
⎢
⎥
∂z ⎥
⎢ ∂y
⎡ k y Dx ⎤
⎢ ∂H x ∂H z ⎥ ∂ ⎢
⎥ 1
−
⎢
⎥ = ⎢ k z Dy ⎥ +
ε
∂x ⎥ ∂t
⎢ ∂z
⎢⎣ k x Dz ⎥⎦ 0
⎢ ∂H y ∂H x ⎥
⎢ ∂x − ∂y ⎥
⎣
⎦
⎡σ y Dx ⎤ ⎡ J x ⎤
⎢σ D ⎥ + ⎢ J ⎥
⎢ z y⎥ ⎢ y⎥
⎢⎣σ x Dz ⎥⎦ ⎢⎣ J z ⎥⎦
(9)
となる。
PML の誘電率に関する構成方程式(6)式からは
⎧ sx Dx = ε 0ε r sz Ex
⎪
⎨ s y D y = ε 0ε r s x E y
⎪s D = ε ε s E
0 r y z
⎩ z z
(10)
が得られ、テンソル成分を用いて書き下し、
⎧
σx
σz
⎪ ( jωκ x + ) Dx = ε 0ε r ( jωκ z + ) Ex
ε0
ε0
⎪
⎪⎪
σy
σx
⎨( jωκ y + ) Dy = ε 0ε r ( jωκ x + ) E y
ε0
ε0
⎪
⎪
σ
σ
⎪ ( jωκ z + z ) Dz = ε 0ε r ( jωκ y + y ) Ez
ε0
ε0
⎪⎩
(11)
さらに逆フーリエ変換によって時間に関する微分方程式
⎧∂
⎡∂
⎤
σx
σ
Dx = ε 0ε r ⎢ (κ z Ex ) + z Ex ⎥
⎪ (κ x Dx ) +
ε0
ε0 ⎦
⎣ ∂t
⎪ ∂t
⎪
σy
⎡∂
⎤
σ
⎪∂
Dy = ε 0ε r ⎢ (κ x E y ) + x E y ⎥
⎨ (κ y Dy ) +
ε0
ε0 ⎦
⎣ ∂t
⎪ ∂t
⎪
σ
⎪ ∂ (κ z Dz ) + σ z Dz = ε 0ε r ⎡⎢ ∂ (κ y Ez ) + y Ez ⎤⎥
⎪⎩ ∂t
ε0
ε0 ⎦
⎣ ∂t
(12)
を得る。
同様に Faraday の法則(2b)についても、(5)の関係に注意し、
⎡ s y sz
⎤
⎡ ∂Ez ∂E y ⎤
Hx ⎥
⎢
−
⎢
⎥
∂z ⎥
⎢ sx
⎥
⎢ ∂y
⎢
⎥
⎢ ∂Ex ∂Ez ⎥
sz sx
⎢
⎥
j
H
ωμ
−
=
−
⎢
⎥
y
z
x
s
∂
∂
⎢
⎥
y
⎢
⎥
⎢
⎥
⎢ ∂E y ∂Ex ⎥
sx s y
⎢
−
⎢
Hz ⎥
∂y ⎦⎥
⎢⎣ sz
⎥⎦
⎣ ∂x
⎡ s y Bx ⎤
= − jωμ ⎢⎢ sz By ⎥⎥
⎢⎣ sx Bz ⎥⎦
⎡ k y Bx ⎤
1
= − jω ⎢⎢ k z By ⎥⎥ −
μ0
⎢⎣ k x Bz ⎥⎦
⎡σ *y Bx ⎤
⎢ * ⎥
⎢σ z By ⎥
⎢σ x* Bz ⎥
⎣
⎦
⎡ k y Bx ⎤
1
= − jω ⎢⎢ k z By ⎥⎥ −
ε0
⎣⎢ k x Bz ⎥⎦
⎡σ y Bx ⎤
⎢
⎥
⎢σ z By ⎥
⎢σ x Bz ⎥
⎣
⎦
と書け、これを逆フーリエ変換し
(13)
⎡ ∂Ez ∂E y ⎤
−
⎢
⎥
∂
y
∂z ⎥
⎢
⎡ k y Bx ⎤
⎢ ∂Ex ∂Ez ⎥
∂ ⎢
⎥ 1
−
⎢
⎥ = − ⎢ k z By ⎥ −
μ0
∂x ⎥
∂t
⎢ ∂z
⎢⎣ k x Bz ⎥⎦
⎢ ∂E y ∂Ex ⎥
⎢ ∂x − ∂y ⎥
⎣
⎦
⎡σ *y Bx ⎤
⎢ * ⎥
⎢σ z By ⎥
⎢σ x* Bz ⎥
⎣
⎦
⎡ k y Bx ⎤
∂ ⎢
1
= − ⎢ k z By ⎥⎥ −
ε
∂t
⎢⎣ k x Bz ⎥⎦ 0
⎡σ y Bx ⎤
⎢
⎥
⎢σ z By ⎥
⎢σ x Bz ⎥
⎣
⎦
(14)
を得る。
PML の透磁率に関する構成方程式(7)式からは
⎧ sx Bx = μ0 μ r sz H x
⎪
⎨ s y By = μ0 μ r sx H y
⎪s B = μ μ s H
0 r y
z
⎩ z z
(15)
が得られ、対応する時間に依存する微分方程式は、(12)式を参考に
⎧∂
⎡∂
σ*
σ* ⎤
⎪ (κ x Bx ) + x Bx = μ0 μr ⎢ (κ z H x ) + z H x ⎥
μ0
μ0
⎪ ∂t
⎣ ∂t
⎦
⎪
*
*
σy
⎡∂
⎤
σ
⎪∂
By = μ0 μr ⎢ (κ x H y ) + x H y ⎥
⎨ (κ y By ) +
μ0
μ0
⎣ ∂t
⎦
⎪ ∂t
⎪
*
*
⎪ ∂ (κ B ) + σ z B = μ μ ⎡⎢ ∂ (κ H ) + σ y H ⎤⎥
0 r
y
z
⎪ ∂t z z μ0 z
μ0 z ⎥⎦
⎢⎣ ∂t
⎩
(16)
となる。
以上の式より２次元 TE 波に関する方程式を抽出すると、(9)式第２式より
σ z Dy
∂H x ∂H z ∂
−
= ( k z Dy ) +
+ Jy
ε0
∂z
∂x
∂t
(17)
(12)式第２式より
σy
⎡∂
⎤
σ
∂
(κ y Dy ) +
Dy = ε 0ε r ⎢ (κ x E y ) + x E y ⎥
ε0
ε0 ⎦
∂t
⎣ ∂t
(18)
(14)式の第１式および第３式で電磁界成分（ E y , Dy , H x , Bx , H z , Bz ）に関する項のみ考慮し
−
∂E y
∂z
∂E y
∂x
=−
=−
σ y Bx
∂
(k y Bx ) −
ε0
∂t
σ B
∂
(k x Bz ) − x z
∂t
ε0
(19)
(20)
(16)式の第１式および第３式と(5)の関係を用いて
⎡∂
⎤
σ
σ
∂
(κ x Bx ) + x Bx = μ0 μr ⎢ (κ z H x ) + z H x ⎥
ε0
ε0
∂t
⎣ ∂t
⎦
(21)
⎡∂
⎤
σy
σ
∂
Hz ⎥
(κ z Bz ) + z Bz = μ0 μ r ⎢ (κ y H z ) +
∂t
ε0
ε0
⎥⎦
⎣⎢ ∂t
(22)
を得る。
３．２次元ＴＥ波におけるＰＭＬの差分式の導出
前節(17)式から(22)式を差分化し FDTD 法による時間発展方程式を導出する。差分化には２次
元 TE 波で用いたものと同様の Yee 格子を用い、それぞれ E と D、H と B を同じ格子点にとる。
ｘｚ平面内の
ＴＥ波
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
H xn,+i ,1/2
k +1/2
E yn,i ,k
i
E yn,i ,k +1
H zn,+i −1/2
1/2, k
i-1/2
o
k-1/2
k
k+1/2
k+1
z(k)
(17)式を差分化し、
n +1/ 2
H xn,+i ,1/2
k +1/ 2 − H x ,i , k −1/ 2
Δz
−
2
H zn,+i +1/1/22,k − H zn,+i −1/1/2,
k
Δx
= κz
Dyn,+i 1, k − Dyn,i ,k
Δt
n +1
n
σ z D y ,i , k + D y , i , k
+ J yn,+i 1/,k 2
+
2
ε0
(23)
n +1
となり、これを Dy について解くと
Dyn,+i 1,k =
n +1/2
⎛ H xn,+i ,1/k +21/ 2 − H xn,+i ,1/k −21/2 H zn,+i +1/2
⎞
2ε 0κ z − σ z Δt n
2ε 0 Δt
1/ 2, k − H z ,i −1/ 2, k
− J xn,+i ,1/k 2 ⎟
−
D y ,i , k +
⎜⎜
⎟
Δz
Δx
2ε 0κ z + σ z Δt
2ε 0κ z + σ z Δt ⎝
⎠
(24)
を得る。
(18)式を差分化し、
n +1
n
n +1
n
E yn,+i 1, k − E yn,i ,k σ x E yn,+i 1,k + E yn,i , k
1 ⎡ D y , i , k − D y , i , k σ y D y , i , k + D y ,i , k ⎤
+
=
+
κ
κ
⎢
⎥
x
Δt
Δt
ε r ε 0 ⎢⎣ y
ε0
ε0
2
2
⎥⎦
n +1
となり、これを E y について解き
(25)
E yn,+i 1,k =
2ε 0κ x − σ x Δt n
1
⎡( 2ε 0κ y + σ y Δt ) Dyn,+i 1,k − ( 2ε 0κ y + σ y Δt ) Dyn,i ,k ⎤
E y ,i , k +
⎦
2ε 0κ x + σ x Δt
ε r ε 0 ( 2ε 0κ x + σ x Δt ) ⎣
(26)
を得る。
(19)式の差分式は、
E yn,i ,k +1 − E yn,i ,k
−
= −κ y
Δz
n +1/ 2
となり、これを Bx
n +1/2
x ,i , k +1/2
B
=
n −1/2
Bxn,+i ,1/2
k +1/2 − Bx ,i , k +1/2
Δt
−
n −1/2
σ *y Bxn,+i ,1/2
k +1/2 + Bx ,i , k +1/2
2
μ0
(27)
について解き
2 μ0κ y − σ *y Δt
2 μ0κ y + σ *y Δt
n −1/2
x ,i , k +1/2
B
⎧⎪ E yn,i ,k +1 − E yn,i , k ⎫⎪
2 μ0 Δt
−
⎨−
⎬
Δz
2μ0κ y + σ *y Δt ⎩⎪
⎭⎪
(28)
を得る。
同様に(20)式の差分式は、
E yn,i +1,k − E yn,i ,k
Δx
n +1/ 2
となり、これを Bz
n +1/2
z ,i +1/2, k
B
= −κ x
n −1/2
Bzn,+i +1/2
1/2, k − Bz ,i +1/2, k
Δt
−
n +1/2
n −1/2
σ x* Bz ,i +1/2,k + Bz ,i +1/2,k
2
μ0
(29)
について解き
⎧⎪ E yn,i +1,k − E yn,i ,k ⎫⎪
2 μ0κ x − σ x*Δt n −1/2
2 μ0 Δt
Bz ,i +1/2,k −
=
⎨
⎬
Δx
2 μ0κ x + σ x*Δt
2 μ0κ x + σ x*Δt ⎩⎪
⎭⎪
(30)
を得る。
(21) (22)式の差分化も同様に行う。導出は各自の課題とする。
４．ＰＭＬにおける材料パラメータの決定
ＰＭＬの厚みをｄとし、その領域内で電磁波が完全に吸収されるようにするため、誘電率と導
電率に相当するパラメータ κξ , σ ξ （ただし、 ξ = x, y , z ）を次のように決定する。
ｙ
TE 波
ｙ方向：
sy = κ y +
σy
= 1, κ y = 1, σ y = 0
jωε 0
ξ
ｘ方向：
σ
sx = κ x + x
jωε 0
ξ
ｚ方向：
σ
sz = κ z + z
jωε 0
x
解析領域
ξ
PML 領域（厚さ d）
z
ξ
図中の PML 各領域において内側から外側への座標を ξ とし、
⎛ξ ⎞
σ (ξ ) = ⎜ ⎟ σ max ,
⎝d ⎠
m
σ max = −
(m + 1) ln{R(0)}
m +1
≈ σ opt =
2ηε r d
150π ε r Δx
ただし、η は空間のインピーダンス、 R (0) は ξ = 0 の境界において達成したい反射率、および
⎛ξ ⎞
κ (ξ ) = 1 + (κ max − 1) ⎜ ⎟ ,
⎝d ⎠
m
κ max = 1 ∼ 20
のように決定して電磁波を減衰させ、内側の解析領域内では
σ = 0, κ = 1
とすることにより通常の空間を実現する。