WE-FDTD(波動方程式FDTD)によるH面誘電体装荷共振器の散乱係数の計算（２）－Uniaxial PMLを用いた場合－

前回、WE-FDTDをMurの吸収境界条件を用いて実装しましたが、今回は完全整合層(PML)を用いた場合を扱います。

WE-FDTDでPMLを適用する方法は、Zhou et alによって定式化されているらしいのですが、Web上で参照できなかったので下記文献を参考にしました。

(1) Yotka Rickard

Improved Absorbing Boundary Conditions for Time-domain Methods in Electromagnetics

http://digitalcommons.mcmaster.ca/cgi/viewcontent.cgi?article=2258&context=opendissertations

McMaster University Dept. of Electrical and Computer Engineering

McMaster University, 2002

5.2.1. Derivation of the IPML equations

(2) Zhou, D., Huang, W. P., Xu, C. L., Fang, D. G., & Chen, B.

"The perfectly matched layer boundary condition for scalar finite-difference time-domain method."

Photonics Technology Letters, IEEE 13.5 (2001), p.454-456.

　原文は確認していませんが、オリジナルはこちらのようです。

PMLを用いた場合補助変数の導入が必要となり、さらにその補助変数を割り当てる節点は電界格子の節点とは異なる位置になります。

WE-FDTD (TMz) Uniaxial PML:

　D.Zhou et al.

　　jωPzx = (1 / sx) d Ez / dx

　　jωFzx = (1 / sx) d (jω Pzx) / dx

　　　ここでsx = κx + σx / (jω ε)

　　Pzy, Fzyも同様に定義します。

　ヘルムホルツ方程式より

　　(jω)^2 Ez = {1/ (εμ)} { (jω)Fzx + (jω)Fzy }

　(jω) → d/dtに変換し、離散化すると、

　(1) Ez(i, j)n+1 = 2 Ez(i, j)n - Ez(i, j)n-1

　　+ (⊿t / (εμ)) { Fzx(i, j)n - Fzx(i, j)n-1 + Fzy(i, j)n - Fzy(i, j)n-1 }

　(2) Fzx(i, j)n+1 = c1fzx(i, j) Fzx(i, j)n

　　+ c2fzx(i, j) { Pzx(i + 0.5, j)n+1 - Pzx(i - 0.5, j)n+1

　　　　　　　　- Pzx(i + 0.5, j)n + Pzx(i - 0.5, j)n }

　　ここで、

　　c1fzx(i, j) = (2 ε κx - σx ⊿t ) / (2 ε κx + σx ⊿t )

　　c2fzx(i, j) = 2 ε / { (2 ε κx + σx ⊿t ) ⊿x }

　　同様にFzy(i, j)も導出できます。

　(3) Pzx(i + 0.5, j)n+1 = c1pzx(i + 0.5, j) Pzx(i + 0.5, j)n

　　+ c2pzx(i + 0.5, j) { Ez(i + 1, j)n+1 - Ez(i, j)n+1 }

　　ここで、

　　c1pzx(i, j) = (2 ε κx - σx ⊿t ) / (2 ε κx + σx ⊿t )

　　c2pzx(i, j) = (2 ε ⊿t) / { (2 ε κx + σx ⊿t ) ⊿x }

　　同様にPzy(i, j+ 0,5)も導出できます。

計算は(3)→(2)→(1)の順で行います。

PMLを用いた場合のWE-FDTDでは、Ez（現在値、１つ前の値、２つ前の値）、補助変数Fzx, Fzy, Pzx, Pzy(現在値と１つ前の値）が必要となり、計11変数を扱うことになります。標準FDTDの場合、Ez、Hx、Hy(現在値、1つ前兼用)、及びDz、Bx、By(現在値と１つ前の値)の計9変数が必要だったのに比べると逆に変数が増える結果になっています。また、PMLを用いた場合のWE-FDTDでは、Ez、Fzx、Fzyは同じ電界の格子節点、Pzx、Pzyはそれぞれx方向、y方向に0.5ずらした位置の節点となり、磁界の格子が必要となります。

H面誘電体装荷共振器の散乱係数の計算

　導波管入出力部にPMLを装荷して計算した結果を示します。

比誘電率分布

f:id:ryujimiya:20130426024013p:plain