時間領域FEMによるH面導波管誘電体共振器の散乱係数の計算

FDTD、WE-FDTDと計算してきて時間領域の計算に慣れてきたので、今度はFEMで時間領域の計算をしてみました。周波数領域のFEMと異なり固有モード展開による入出力ポート境界条件が使えなさそうなので、なんらかの無反射条件を課す必要があります。ここでは最も簡単な１次の吸収境界条件(1st order ABC,)を使用しました。

　時間領域のFEM方程式は、H面導波管(TMzモード)の場合、

　[K] Ez + [M] d^2 Ez / d^2 = - (1/vp) [Qb] dEz / dt

　　vp = ω / βx

　　 βx = √{ k0^2 - (π / W)^2 }

　　[K]ij = <dNi/ dx, (μr^-1) dNj / dx>Ω + <dNi / dy,  (μr-1) dNj / dy> Ω

　　[M]ij = <Ni, (ε0 μ0 εr) Nj> Ω

　　[Qb]i,= <Ni, (μr^-1) Nj> Γ

となります。 Ω:解析領域、Γ:入出力境界(1st order ABC適用)です。また

　　<a, b>Ω = ∬a b dΩ

　　<a, b>Γ = ∫ a b dΓ

です。

時間に関してはNewmark Betaのγ = 0.5、β = 0.25を適用して離散化しました。離散化すると

　( [M] / ⊿t^2 + β[K] + (1 / vp)(1 / (2⊿t)) [Qb]  } Ezn+1

　　= { (2 / ⊿t^2) [M] - (1 - 2β) [K] } Ezn + { - [M]  / ⊿t^2 - β[K] + (1 / vp)(1 / (2⊿t)) [Qb] } Ezn-1

となります。ここで

　d^2u/dt^2 = (u n+1 - 2 u n + u n-1) / ⊿t^2

　du / dt = (u n+1 - u n-1) / (2⊿t)

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

を用いました。β = 1/4のとき無条件で安定(unconditionally stable)となるようです。

入射(励振)項は、

　(- 2 / vp) <Ni, μr-1 Nj> d Ez inc / dt 

を右辺に追加して計算しました。

H面導波管誘電体共振器は、FDTDやWE-FDTDで扱ったのと同じ寸法のものを扱います。領域分割は、FDTDの場合、導波管幅を80分割しましたが、今回のFEMでは20分割とします。

励振はガウシアンパルス（時間変化)、TE10モード分布とします。

　分割幅dl = 3.0 x 10 -3 m

　導波管幅W = 80 x ⊿l

　時間刻み⊿t = 0.5 x ⊿l / √2

　Ezinc = - E0 exp ( - (t- t0)^2 / (2 Tp^2) )

　　t0 = 30 ⊿t 

　　Tp = t0 / (4√2)

　　E0: TE10モード分布

　　ポート1境界から5⊿l の位置で励振

　入出力導波管の長さ：20⊿l

　観測点

　　ポート1：励振源から5⊿l  (境界から10⊿l)

　　ポート2：境界から - 5⊿l

比誘電率分布

散乱係数の計算結果は次の通りです。

周波数特性の傾向は周波数領域FEMと同じ物が得られましたが、あまり精度は良くないようです。ちなみにFDTDでもMurの吸収境界条件(今回の1st order ABCと等価だと思います)を用いた場合はよく似た結果が得られたので、恐らく吸収境界条件の問題だと思われます。

なお、散乱係数の周波数特性を計算するだけなら、周波数領域で解いた方が短時間で精度の良い周波数特性が得られるので、時間領域で解くメリットはあまりないと思われますが、FDTDのようにガウシアンパルス + FFTを使って散乱係数を計算をしてみたかったので試してみました。

【追記 (2013-5-19)】

上記記事では、入射項を境界積分<>Γで計算しましたが、下記のように領域積分<>Ωで励振する方法を試したので記します。

　入射項：

　　　[K] { Ez inc } + [M] d^2 { Ez inc } / dt^2

ここで、{ Ez inc }は領域全体の入射電界分布ですが、境界Γ上の入射電界のみを指定しました。結果は下記のとおりです。