ベクトル直交基底関数を用いた時間領域FEMによるH面導波管の散乱係数計算

磁界ベクトルを用いた時間領域FEMの場合、質量行列を対角化する方法として直交基底関数を用いる方法があるようです。

引用元はこちら。

D. A. White

"Orthogonal vector basis functions for time domain finite element solution of the vector wave equation"

https://e-reports-ext.llnl.gov/pdf/233577.pdf

Lawrence Livermore National Laboratory, Center for Applied Scientific Computing,

UCRL-JC-129188, December 1997

IEEE Trans. Magn., vol.35, pp. 1458- 1461, May 1999

スカラ場の場合、質量行列を対角化する方法としては集中質量近似がありますが、ベクトル場では適用できないのか、関連文献を見つけることができませんでした。

通常、磁界ベクトルの補間関数には辺要素の補間関数が用いられますが、この補間関数は直交基底ではないので質量行列は非対角項を含みます。もし、補間関数を直交基底にすれば<Ni, Nj> = 0 (i != j) となるので、質量行列は対角項のみとなります。

このような補間関数の構成方法が上記文献に記されていました。

文献では、異なる節点の補間関数の内積のガウス数値積分した結果が０となるように補間関数を構成しています。

つまり、ノルムを

　<Ni, Nj> = Σk =[1,3] α k Ni ( mk ) Nj ( mk ) 

　mk：辺kの中点

と近似すれば、辺の中点だけで直交関係を満足するような補間関数を構成すればよいことになります。

上の数値積分は、２次元問題のハンマーらの積分公式の次数２、積分点３の式に相当するので、

　α k = (1 / 2) (1 / 3)  |J|

　(|J| = 2Ae:ヤコビアン)

です。

さて、直交基底ベクトル補間関数は、辺方向成分に対する３つの補間関数と、辺の法線方向成分に対する３つの補間関数からなります。

　Ni = Zi         i = [1, 3]

　        Bi- 3    i = [4, 6]

　Zi = Wi - Σ k=[1,3]  ( <Wi, Bk> / <Bk, Bk> ) Bk 

　Bi = Lj Lk ni

ここで

　Wi = Lj ∇Lk - Lk∇Lj

　Li:節点iの面積座標

Ziは要素境界で接線方向成分が連続ですが、Biは要素間で連続性は課さないので内部節点の扱いです。

時間領域FEM + 1次ABCにベクトル直交基底関数適用した結果を示します。解析対象はいつも通りH面導波管共振器です。

メモリ不足でいつもの領域はとれなかったので、入力導波管長さを10⊿l、出力導波管長さを5⊿lにしました。励振位置は 1⊿l で、観測点は励振位置から5⊿l、ポート2端から1⊿lとしました。

微分ガウシアンパルスを励振したときの散乱係数の計算結果は次の通りです。

上記は、Newmark β法で。γ=0.5、β = 0として、係数行列が質量行列のみになるようにして計算しました。したがって、無条件安定でなく、時刻刻みにFDTDのような制限があります。このモデルでは0.38 ⊿l / c0 で計算した結果です。（これより時間刻みを緩和すると発散しました)そのため、計算ステップ数はいつもの倍のt6000回としています、しかしながら、方程式を解く必要がないため計算時間は、Newmark β = 0.25、計算ステップ回数3000回の場合と比べると短くて済みました。

なお、Newmark β法で、β = 0.25としたときの計算結果は次の通りです。この場合は、無条件安定ですが、方程式を解く必要があります。時間刻み幅 0.5⊿l / c0、計算ステップ数 3000回としています。