位相速度の周波数特性を考慮した時間領域FEMの1次吸収境界条件
時間領域FEMで吸収境界条件(ABC)を用いる場合、高次吸収境界条件を用いないと計算結果が劣化しました。
時間領域FEMの高次吸収境界条件(Higher oder ABC) - ryujimiyaの日記
これはABCに用いる波の速度vpを周波数に関係なく一定としていることに起因しているように思われます。
実際、周波数領域で計算したとき、ABCで用いる波の速度を周波数に合わせて変化させると計算結果はモード展開を適用したときと同等の結果が得られました。
1次吸収境界条件(first-order ABC)と周波数領域FEMによるH面導波管誘電体共振器の散乱パラメータの計算 - ryujimiyaの日記
今回は、時間領域FEMで1次ABCで用いる波の速度に周波数特性を考慮してみました。
一般的な導波路の波の速度(ここでは、位相速度)は解析解が得られないので、ABCに速度の周波数特性を考慮するのは難しいですが、H面中空導波管に限定すると位相速度vp(ω)は
vp(ω) = ω / β(ω) (1)
β(ω) = √{ k0^2 - (π/W)^2 } (2)
なので何とかなりそうです。
位相速度の周波数特性を考慮した1次ABC
1次ABCでは、残差の境界積分は
-(1/vp) ∫ (1/μr) u du/dt dΓ
で評価しました。これを周波数領域に変換すると
-jω/vp ∫ (1/μr)u u dΓ = -jβ ∫ (1/μr)u u dΓ
です、この式に(2)を代入すると
-j√{ k0^2 - kc^2 } ∫ (1/μr) u u dΓ = √{ kc^2 - k0^2 } ∫ (1/μr) u u dΓ (3)
ただし、
kc = π / W : TE10モードの遮断波数
となります、したがって、√{ kc^2 + k0^2 }の逆フーリエ変換が分かれは、1次ABCを
∫ (1/μr) u (F-1[√{ kc^2 + k0^2 } ] ★ u] dΓ (4)
★:convolution
とすれば位相速度の周波数依存を考慮できると思います。
√{ kc^2 - k0^2 }のフーリエ逆変換は、次のラプラス変換の公式を利用して求めることができます。
L-1[ √{ s^2 + 1 } - s ] = J1(t) / t (5)
Ji(t):i次第一種ベッセル関数
この公式に√{ kc^2 - k0^2 }を当てはめると、jω を sと置くと
√{ kc^2 - k0^2 }
= √ { kc^2 + (s / c0)^2 }
= kc √{ 1 + (s / (kc c0))^2 }
= kc [ (s / (kc c0)) + √{ 1 + (s / (kc c0))^2 } - (s / (kc c0)) ]
= s / c0 + kc [√{ 1 + (s / (kc c0))^2 } - (s / (kc c0)) ] (6)
よって逆ラプラス変換を施すと
L-1[√{ kc^2 - k0^2 }]
= (1/c0) d/dt + kc (kc c0) J1(kc c0 t) / (kc c0 t)
= (1/c0) d/dt + kc J1(kc c0 t) / t (7)
となります。
したがって位相速度の周波数特性を考慮した1次ABCは、
- 1/c0 ∫ (1/μr) u du/dt dΓ - ∫ (1/μr) u [kc J1(kc c0 t) / t] ★ u dΓ (8)
となります。
領域を三角形要素で離散化すると
- 1/ c0 [Qb] du/dt - [Qb] Φ (9)
Φ = g ★ u
g = kc J1(kc c0 t) / t
[Qb]ij = <Ni, 1/μr Nj>
(9)式を残差に加算し、u, ΦをNewmark β法で時間方向に離散化すればFEM行列方程式の完成です。
なお、畳み込み積分は、漸化式がないので
Φ(n⊿t) = Σm = [0, n] g(m⊿t) u((n - m)⊿t) ⊿t
で評価しました。
直線導波管の散乱係数の計算結果
誘電体装荷共振器の散乱係数の計算結果
比誘電率分布
散乱係数の計算結果
どちらも従来の1次ABCを適用した場合より大幅に計算結果が改善されました。
まあ、この方法は中空導波管しか適用できないのであまり利用価値はないと思いますが、うまくいきそうだったので試してみました。
