時間領域FEMの高次吸収境界条件(Higher oder ABC)
時間領域FEMに適用可能な高次吸収境界条件(Higer oder absorbing boundary condition)の構成方法が下記の文献に記されていたので試してみました。
Dan Givoli , Beny Neta & Igor Patlashenko
"Finite element analysis of time-dependent semi-infinite wave-guides with high-order boundary treatment"
http://www.math.nps.navy.mil/~bneta/higdon_fe_pap.pdf
February 27, 2003
International Journal for Numerical Methods in Engineering
Volume 58, Issue 13, p.1955- 1983,
7 December 2003
時間領域のスカラヘルムホルツ方程式
∇2 u - (1 / c0^2) d2 u / dt2 = 0 (1)
c0：真空中の光速
HigdonのABC
[ Π j = [1, J] ( d / dx + (1 / cj) d / dt ) ] u = 0 (2)
補助変数 φ j (j = 1, ...., J - 1)を導入して書き換えると、
( d / dx + (1 / c1) d / dt ) u = φ1
( d / dx + (1 / c2) d / dt) φ1 = φ2
( d / dx + (1 / cJ) d / dt) φJ - 1 = 0
となります。または、
φ0 = u ( 定義)
φJ = 0 ( 定義)
とすると
( d / dx + (1 / cj) d / dt )φ j - 1 = φ j j = 1, ... J
と書くことができます。
φjは元のスカラヘルムホルツ方程式を満たすことを利用して、d / dxの項をd / dyで表現すると、最終的には
β0 d u / dt + d u / dx = φ1 (3)
βj d φ j / dt - α j d2 φ j - 1 / dt2 - d2 φ j - 1 / dy2 = φ j + 1 (4)
α j = 1 / (c j ^2) - 1 / (c0^2)
β0 = 1 / c1
β j = 1 / c j - 1 / c j + 1
となります。
(1)、(3)、(4)が解く方程式です。
時間領域FEMでは、(1)に対する弱形式の境界積分に(3)を代入して、uとφ1の方程式を作り、(4)に対する弱形式と連立して解きます。
また文献に述べられていますが、
(1)の弱形式から φ1を既知として、uを求める
(4)の弱形式のj = 1の場合について、φ2を既知として、φ1を求める
j = 2, ... J - 1 の場合について同様にしてφ2, ..., φJ - 1を求める
この一連の処理では、φ1等を既知と仮定したので誤差が含まれます。そこで、上記処理を反復することで精度を上げます。
実際反復処理を行わないと解が発散することがありました。なお、反復はもう1回だけ追加すれば充分のようです。
時間領域FEMの基本式は、文献の式(64)、式(65)に示されています。また時間方向に離散化した式は式(81)～式(91)に示されています。
Ｈ面導波管誘電体共振器の散乱係数の計算結果
比誘電率分布
cj をすべてc0 とした場合の結果を示します。なお、励振はTE10モードをガウシアンパルス時間変化で励振しました。
1次ABC
2次ABC
3次ABC
4次ABC
5次ABC
6次ABC
ABCの波の速度をすべて光速としても、5次のABCにすると良好な結果が得られました。高次のABCにしても未知数はそれほど増えませんし、PMLのときのような時間ステップ毎のFEM行列構築もないので計算時間はそんなにかかりません。
時間領域FEMをいろいろ試してきましたが、いまのところ今回の高次ABCによる方法が一番利用できそうな方法のようです。
＜補足＞
文献では時間方向の離散化で、uの速度(時間微分)、加速度(２階の時間微分)を用いています。上で示した結果は速度、加速度を用いた式で計算した結果です。
ただ、文献に定式化が示されていなければ、uだけの式で計算していたと思います。
つまり、
d2 u / dt2 = (u n + 1 - 2 u n + u n - 1) / ⊿t^2
d u / dt = (u n + 1 - u n - 1) / (2⊿t)
u = β u n + 1 + ( 1 - 2β) u n + β u n - 1
β = 0.25
として計算していたと思います。こちらの方法で計算した結果を次に示します。
この場合、反復計算で用いるφ1等の初期値は、直前のステップの値を使用しました。
1次ABC
2次ABC
3次ABC
4次ABC
5次ABC
6次ABC
