時間領域FEMによるH面導波管誘電体共振器の散乱係数の計算（５） －磁界ベクトルとUniaxial PMLを用いた場合－

前回、Split-TypeのPMLと磁界ベクトルを用いた時間領域FEMの計算を行いましたが、あまり良い結果ではありませんでした。

今回はUniaxial PMLを使った場合の計算結果を記します。

引用文献はこちら。

Dan Jiao, Jian-Ming Jin, Eric Michielssen, and Douglas J. Riley

"Time-domain finite-element simulation of three-dimensional scattering and radiation problems using perfectly matched layers."

https://engineering.purdue.edu/~djiao/publications/DanJiao_pml3d.pdf

IEEE Transactions on Antennas and Propagation

vol. 51, no. 2, February 2003

文献では３次元電界ベクトルとPMLを用いた定式化が記されていますが、ここでは、H面導波管回路を扱うので２次元磁界ベクトルで定式化します。

X方向PML内の支配方程式は

　∇x ε-1Λ-1 ∇x H + μΛd2 H / d2 t = 0

   H = [ Hx, Hy, 0]T

　Λ = [ 1 / sx    0    0   ]

　      [    0      sx      0  ]

          [   0        0      sx ]

ガラーキン法を適用すると

　μ<Ni, Λ★d2 H / dt2>Ω + <∇xNi, ε-1Λ-1★∇x H>Ω = 0

ここで、 

　1/sx ⇒ δ(t) - (σx / ε) exp( - σx t / ε) u(t)

　　δ(t):デルタ関数

　　u(t):ステップ関数

　-ω2 sx u ⇒ d2 u / dt2 + (σx / ε) d u / dt

を用い、HをNjで

　　　H = Σ Nj uj

と離散化すると、要素内の方程式として

　[M] d2 u / d2 t> + [Mσ] d u / dt - [M1] d2 ψ / d2 t + [Kσ] ( u - ψ) = 0

　　[M] = ε0μ0 <Ni, μr Nj>

　　[Mσ] = ε0μ0 (σx / ε)<Ni, μr yy・Nj>

　　[M1] = ε0μ0 <Ni, μr xx・Nj>

　　[Kσ] = <∇x Ni, εr-1 zz・∇x Nj>

　　x, y, z：単位ベクトル

　　ψ = (σx / ε) exp ( - σx t / ε) ★ u

が得られます。この式は文献の(29)式に対応します。(補足：TMzのとき、[K] = <∇xNi, εr-1∇xNj>と[Kσ]は一致)

時間に関する離散化は、Newmark Beta法を用いました。ψは、uの関数ですがここでは独立変数として扱い、

　ψ = β ψ n +1 + (1 - 2β) ψ n + βψ n - 1 

　d2 ψ / d2 t = (1 / ⊿t^2) ( ψ n+1 - 2 ψ n + ψ n - 1 )

と離散化しました。

微分ガウシアンパルスで励振したとき

　PMLの層数は10としました。

比誘電率分布

σxの分布

散乱係数の計算結果(微分ガウシアンパルス)

ガウシアンパルスを励振したとき

　こちらは、励振位置から観測点を25⊿lだけ離したときの結果を示します。（微分ガウシアンパルスの場合は、励振位置と観測点の間の距離は5⊿lなので、5倍の距離を取っていることになります。

磁界ベクトル＋Uniaxial PMLの組み合わせで解析した場合は良好な結果が得られたようです。

これまでの計算結果からTMzに関して

●FDTDの場合は、Split-Type、UniaxialどちらのPMLでも良好な結果が得られるようです。

●FEMの場合は、Ezで計算した場合、Split-Type PMLとUniaxial PMLは等価で、良好な結果が得られます。磁界ベクトルH = [Hx, Hy]Tで計算した場合は、Uniaxial PMLでないと良好な結果が得られないようです。