時間領域FEMにおける完全整合層の定式化ーTsaiらによる方法－
時間領域FEMで完全整合層(PML:Perfectly Matched Layer)を用いる場合の定式化にはいくつかあるようです。
いままで用いてきたのは、Jiaoらによる定式化です。彼らの定式化による結果は、時間領域FEMで解く場合もっとも良好な結果が得られています。
今回、下記の文献のTsaiらによる定式化を試してみました。
Hsiao-Ping Tsai, Yuanxun Wang, and Tatsuo Itoh
"An unconditionally stable finite element time domain solution of active nonlinear microwave circuits using perfectly matched layers"
http://www.ee.ucla.edu/~mwlab/publications/2001c/MTTPhoenix/hp.pdf
Microwave Symposium Digest, 2001 IEEE MTT-S International. IEEE, 2001. p. 2059-2062
上記の定式化では、係数行列が時間に依存しないのでこれはいいかもしれないと思い試してみました。
文献の(5)式は、TMzモード、x方向PMLのみの場合、
E = u z
σy = σz = 0
[I] = [ σx 0 0 ]
[ 0 0 0 ]
[ 0 0 0 ]
[J] = [ 0 0 0 ]
[ 0 σx 0 ]
[ 0 0 σx ]
[K] = 0
[L] = 0
を代入すると
d/dx (1/μ d u / dx)
+ d/dy (1/μ du/dy) + 2σx /(jωε0) d/dy(1/μ du /dy) - σx^2 / (ω^2ε0^2) d/dy(1/μ du/dy)
+ ω^2εu - jω 2σx ε / ε0 u - σx^2 ε / ε0^2 u = 0
となります。この式で、
jω → d/dt
-ω^2 → d2/dt2
1 / (jω) → ∫t
- 1 / ω^2 → ∬t
と読み替えると時間領域の支配方程式となります。時間領域の支配方程式にガラーキン法を適用して離散化すると（9)式に対応する式として、
[K] u + [Kyσ1] v + [Kyσ2] w + [M] d2 u / dt2 + [Mσ1] du/dt + [Mσ2] u = 0
ここで
[K]ij = <∇Ni, p ∇Nj>
[M]ij = ε0μ0 <Ni, q Nj>
[Kyσ1]ij = <dNi/dy, (2σx/ε0) p dNj/dy>
[Kyσ2]ij = <dNi/dy, (σx^2/ε0^2) p dNj/dy>
[Mσ1]ij = ε0μ0 <Ni, (2σx/ε0) q Nj>
[Mσ2]ij = ε0μ0 <Ni, (σx^2/ε0^2) q Nj>
p = 1 / μr
q = εr
εr = ε / ε0
μr = μ / μ0
が得られます。
v = ∫u dt
w = ∬u dt2 = ∫v dt
これらの時間積分値は台形公式( Trapezoidal rule )より
v n = v n - 1 + (⊿t / 2) (u n - 1 + u n )
w n = w n - 1 + ⊿t v n - 1 + (⊿t^2 / 4) (u n - 1 + u n )
で計算します。文献では(13)式に記されていますが、2重積分の方の式が間違っているので注意が必要です。（この文献に記されている式は結構記載ミスがあります。)
さてこれまでと同じように、σxはPML開始位置からの2次の分布とし、各々の要素内では、σxは一定とみなして計算してみました。計算に用いた条件は次の通りです。
導波管幅W = 240 mm
セルの刻み幅⊿l = W / 20
時刻刻み幅⊿t = 0.5 ⊿l / (√2 c0)
ガウシアンパルスの遅延時間T0 = 40 ⊿t
ガウシアンパルスの時間幅Tp = T0 / (4 √2)
ガウシアンパルス
Ez inc(t) = E0(y) exp( - (t - T0)^2 / (2 T0^2))
E0(y)：TE10モードの分布
H面導波管直線回路の散乱係数の計算
完全整合層を10層としたとき
σxの分布
散乱係数計算結果
全然駄目でした。
完全整合層を50層としたとき
σxの分布
散乱係数の計算結果
改善はしましたが、Jiaoらの定式化でPMLを10層としたときの計算結果より悪い結果です。
H面導波管誘電体装荷共振器の散乱係数の計算
完全整合層を50層として計算しました。
比誘電率分布
σxの分布
散乱係数の計算結果
3000ステップ後の電界Ezの分布
やはりきちんと吸収できていないように思われます。
つづく。
