時間領域FEMの完全整合層(PML)の定式化（まとめ）
スカラヘルムホルツ方程式に限るとこれまでに次の定式化を試しました。
(1) Jiaoらの定式化
[M] d2 u / dt2 + [Mσ] d u / dt + [K] u + [Kx] Φ + [Ky] v = 0
u = Ez
p = 1/μr
q = εr
Φ = - σx / ε0 exp( - σx t / ε0 ) ★ u
v = σx / ε0 ∫ u dt
[M] = ε0μ0 <Ni, q Nj>
[Mσ] = ε0μ0 <Ni, (σx / ε0) q Nj>
[K] = <∇Ni, p ∇Nj>
[Kx] = <dNi/dx, p dNj/dx>
[Ky] = <dNi/dy, p dNj/dy>
時間領域FEMによるH面導波管誘電体共振器の散乱係数の計算（２） －完全整合層(PML)の場合－
http://ryujimiya.hatenablog.com/entry/2013/05/19/214317
(2) Tsaiらの定式化
[M] d2 u / dt2 + [Mσ1] du / dt + [Mσ2] u + [K] u + [Kyσ1] v + [Kyσ2] w = 0
u = Ez
p = 1/μr
q = εr
v = ∫u dt
w = ∬u dt = ∫v dt
[M] = ε0μ0 <Ni, q Nj>
[Mσ1] = ε0μ0 <Ni, (2σx / ε0) q Nj>
[Mσ2] = ε0μ0 <Ni, (σx^2 / ε0^2) q Nj>
[K] = <∇Ni, p ∇Nj>
[Kyσ1] = <dNi/dy, (2σx / ε0) p dNj/dy>
[Kyσ2] = <dNi/dy, (σx^2 / ε0^2) p dNj/dy>
時間領域FEMにおける完全整合層の定式化ーTsaiらによる方法－
http://ryujimiya.hatenablog.com/entry/2013/06/30/234959
Jiaoの定式化ではPMLの層数は10で充分でしたが、Tsaiらの定式化では、Jiaoの5倍にしても結果はあまりよくありませんでした。
そこで、今回はTsaiらの定式化を検証してみます。
Tsaiらの定式化 の検証
Tsaiらの定式化では支配方程式は、周波数領域で
d/dx (1/μ du / 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/(jωε0) εu - σx^2 / (ω^2ε0^2) ε u = 0
でした。これは、sx = 1 + σx / (jωε0) を使って書き換えると
d/dx (1/μ du / dx) + sx^2 (1/μ du / dy) + sx^2 ω^2ε u = 0 (1)
となります。
一方、Uniaxial PMLの基本式は
d/dx (1/sx 1/μ du / dx) + d/dy (sx 1/μ du /dy) + ω^2 sx ε u = 0 (2)
です。
sxは、y方向には一定なので(2)の第2項のsxはyに関する偏微分の外に出すことができます。
一方第1項はsxがxについて一定とみなせる場合はxに関する偏微分の外に出すことができるので
1/sx d/dx (1/μ du / dx) + sx d/dy (1/μ du / dy) + ω^2 sx ε u = 0
両辺にsxを掛けると
d/dx (1/μ du / dx) + sx^2 (1/μ du / dy) + sx^2 ω^2ε u = 0
となり、(1)の式と一致します。
あれ？ sxはx方向には一定ではないですね。
sx = 1 + σx / (jωε0)
σx = σmax |x - x0|^2 / L^2
L：PMLの長さ
x0：PMLの開始位置x座標
sxがx方向に変化することを考慮して(2)の第1項を展開すると
d/dx (1/sx 1/μ du / dx) = - 1/sx^2 (d sx / dx) (1/μ du / dx) + 1 / sx d/dx(1/μ du / dx) (3)
となり、dsx/dxなる項が生じます。
したがって、(1)の式を修正すると
- 1 / sx (d sx / dx)(1/μ du / dx) + d/dx(1/μ du / dx) + sx^2 (1/μ du / dy) + sx^2 ω^2ε u = 0 (4)
となります。
Tsaiらの定式化で計算したときPMLの層数を多くとらないと波の吸収が悪かったのは、定式化でd sx /dxの項を無視していることに原因がありそうです。
Tsaiらの定式化の修正
(4)式を弱形式にすると周波数領域では
∬ 1/sx (dsx/dx) p u (du/dx) + p du/dx du/dx + p du/dy du/dy + 2σx / (jωε0) p du/dy du/dy - σx^2/(ω^2 ε0^2) p du/dy du/dy
- ω^2 ε0μ0 q u u + jω (2σx / ε0) ε0μ0 q u u + (σx^2/ ε0^2) ε0μ0 q u u dΓ = 0 (5)
よって、Tsaiらの定式化に
∬ 1/sx (dsx/dx) p u (du/dx) dΓ = ∬(1 / jω) (1 / sx) (1 / ε0) (dσx /dx) p u (du/dx)dΓ (6)
なる項を追加することになります。ここで、
dsx /dx = (1 / jωε0) dσx /dx (7)
を用いました。
この項を時間領域に変換すると
∬(1 / ε0) (dσx /dx) p u { F-1[(1 / jω) (1 / sx)] ★ (du/dx) } dΓ (8)
F-1：フーリエ逆変換
{ g ★ f } = ∫τ=[0 t] g(τ)f(t - τ)dτ：畳み込み積分
ここで、
F-1[(1 / jω) (1 / sx)] = exp(-σx t / ε0) u(t)
u(t)：単位ステップ関数
より、
∬(1 / ε0) (dσx /dx) p u { F-1[(1 / jω) (1 / sx)] ★ (du/dx) } dΓ
= ∬(1 / ε0) (dσx /dx) { exp(-σx t / ε0) u(t) ★ (du/dx) } dΓ
= ∬(- 1 / σx) (dσx /dx) p u Φ (9)
ただし、
Φ = ( - σx/ε0) exp(-σx t / ε0) ★ (du/dx) (10)
です。
(8)式をガウスの求積法で数値積分することにすると(10)式は、ガウスの積分点について
Φn = - (1 - exp(-σx ⊿t / ε0)) (du /dx) + exp(-σx ⊿t / ε0) Φn - 1 (11)
として評価します。
また、今回の場合、修正項にdσx / dxが含まれているのでσxの分布は要素内では最低でも1次補間する必要があると思われます。
（Jiaoの定式化や修正前のTsaiの定式化ではσxの要素内分布は一定として計算していました)
Tsaiらの定式化に修正項を追加したときの計算結果
PMLの層数を10としたときのH面導波管直線回路の計算結果は次の通りです。
σxの分布
散乱係数の計算結果
前回の修正項なしの場合と比べると、修正項なしでPML層数50の計算結果より、今回の修正項ありでPML層数10の計算結果の方が良好な結果となりました。
この結果はJiaoらの定式化と同等の精度です。
以上よりTsaiらの定式化には修正項が必要であることが分かりました。
H面導波管誘電体装荷共振器の散乱係数の計算
PMLの層数を10としたときの計算結果を示します。
比誘電率分布
σxの分布
散乱係数の計算結果
こちらは、低周波側では依然計算結果に劣化が見られますが、PMLの層数は前回の修正前の定式化の1/5に低減できました。
ただ、Jiaoらの定式化と比べると計算精度は良くないようです。
