時間領域FEMによるH面導波管誘電体共振器の散乱係数の計算（２） －完全整合層(PML)の場合－
前回、時間領域FEMでH面導波管の散乱係数の計算をしましたが、1st order ABCではあまり良い結果が得られませんでした。今回は、完全整合層を用いて計算したので記します。
完全整合層は、Uniaxial PMLを用いますが、split typeのPMLでも最終的に解く方程式は同じものになります。
Σe = [1, Ne] ∬{ dNi / dx, (1/μr)(1/sx) dEz / dx + dNi / dy, (1/μr)(sx) dEz / dy
- ω^2 ε0 μ0 Ni (εr sx) Ez } dΩ = - ∫ { (1/μr)(1/sx) Ni dEz / dx } dΓ
sx = 1 + σx/ (jωε)
Ne:要素数
です。これを時間領域の汎関数に変換します。
(1 / sx) Ez ⇒ Ez + Φ
sx (dEz / dy) ⇒ dEz /dy + d v / dy
- ω^2 sx Ez ⇒ d^2 Ez / dt^2 + (σx / ε) d Ez / dt
ここで、
Φ = ( - σx / ε) exp ( - σx t / ε ) u(t) ★ Ez
v = ( σx / ε ) ∫ [0, t] Ez d τ
★：畳み込み積分(convolution)
と変換できるので、時間領域のFEMの方程式は
Σe = [1, Ne] { [M] d^2 Ez / d^2 t + (σx / ε ) [M] d Ez / dt + [K] Ez
+ [Kx] Φ + [Ky] v } = 0
となります。ここで
[Kx] = (μr^-1) < dNi / dx, dNj / dx >
[Ky] = (μr^-1) < dNi / dy, dNj / dy >
[K] = [Kx] + [Ky]
[M] = ε0 μ0 μr < Ni, Nj >
です。要素内で媒質は一定値を取るものと仮定しています。
Φの畳み込み積分は、
Φn = { exp( - σx ⊿t / ε) - 1 } Ezn + exp( - σx ⊿t / ε) Φn - 1
の関係を使って計算します。
また、vは
dv / dt = (σx / ε) Ez
の関係を後退差分で離散化して、
v n = σx ⊿t / ε + v n - 1
として計算します。
ここで、注意する点は畳み込み積分のΦを要素単位で計算する必要があるため、1st order ABCのときのようにあらかじめ全ての要素を足し合わせた行列を計算することができず、時間ステップ毎に要素の足し合わせを行う必要があるということです。
なので、かなり時間がかかります。
なお、上記の式は下記文献の（29)式に対応します。文献では３次元問題として定式化していますが、(17)式に一例として示されている式の変形はかなり参考になりました。
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
H面導波管誘電体共振器の散乱係数の計算
比誘電率の分布
領域内部の分割は前回と同じで、分割幅は
⊿l = 12,0 x 10^-3
です。１つのセルは２つの三角形要素から構成しています。
時間ステップ幅も同じで、
⊿t = 0.5 ⊿ l / √2
です。
PMLは入出力導波管にNpml = 10層設けました。これは以前FDTDで計算したときのPMLの長さ(3.0 x 10^-3 x 20)の２倍です。
PMLの媒質パラメータは以下の通りです。
σx = σx max (|x - x0| / Lpml)^2
Lpml ：PMLの長さ
x0：PMLの始点のx座標
σx max = - 3 ε0 c0 ln( R(0) ) / (2 Lpml)
Lpml = ⊿l x N pml
R(0) = 10^-8
σxの分布
入射波は時間変化ガウシアンパルス、TE10モード分布の電界としました。
Ezinc = - E0 exp { - (t - t0)^2 / (2 Tp^2) }
t0 = 30 ⊿ t
Tp = t0 / (4√2)
入射項を境界積分で評価したときの散乱係数の周波数特性の計算結果は次のようになりました。
上側パネルを参照するとわかりますが、port1(Inc)がかなり振動しています、その影響か3000回の時間ステップでは全体の特性を得ることができませんでしたが、得られた2W/λ = 1.0～2.5の周波数特性はかなり良好な結果になっていると思われます。
やはりFEMでも時間領域で解析する場合は、完全整合層を用いた方が良いようです。
入射項を領域積分で評価した場合の結果は次の通りです。
この場合は、ガウシアンパルス のパラメータt0を
t0 = 40 ⊿t
に変更しています(入射波のスペクトルがt0 = 30 ⊿tだと低周波で確保できなかったので）。
入射項を領域積分で評価した場合は、振動はほとんどなく良好な結果と言えます。
どこか間違っているのかな？
