3次元のMaxwell方程式（周波数定常）を有限要素法で解きたくて調べたメモ

3次元のMaxwell方程式（周波数定常）を有限要素法で解きたくて調べたメモ

弱形式
∫( (α∇×E)・(∇×F) - βE・F)dV - ∫(α∇×E)・n dS＝∫ f・F dV
但し、α＝1/μ、β=εω^2、nは領域外周の外向きのベクトル、fはソース項、Eは電界ベクトル、Fは試行関数ベクトル

離散化
・分割された要素内において、電界ベクトルEを形状関数Niと自由度eiを用いて
E = (N1, N2, ...)・(e1, e2, ...)^T と近似する。

形状関数
・Nedelecエッジ要素 (curl-conforming element) を用いる
・一般の形状を考えるのは大変なので、まず簡単な形状のリファレンス要素で考える。その後、座標変換して望みの形状の要素に対する形状関数を得る

Nedelecエッジ要素の具体的な形
・1次要素は、1辺に1自由度を持つ。簡単な四面体のリファレンス要素（頂点がp0=(0,0,0), p1=(1,0,0), p2=(0,1,0), p3=(0,0,1))のとき、形状関数をηijとすると
η01 = (1-y-z, x, x), η02 = (y, 1-x-z, y), η03 = (z, z, 1-x-y)
η12 = (-y, x, 0), η23 = (0, -z, y), η31 = (z, 0, -x)
確かに全てdivを取ると0になる。また、体積座標λiを用いて、ηij = λi∇λj - λj∇λi とも書ける。

・2次以上の高次要素は簡単には書き下せない。代わりに、別の基底ベクトルの組{φ_j}を考え、各辺・面・領域に対応する自由度に対する条件を満たすよう係数u_jを決め、形状関数η_i=Σu_ij φ_jを求める。四面体に対する2次要素の自由度数は、各辺2自由度、各面2自由度の計20。詳細は文献[1]など参照だが、調べ中。

座標変換
・リファレンス要素における3次元座標をξ、望みの四面体要素における3次元座標をxとするとき、リファレンス要素から望みの四面体要素へのアフィン変換Fは
x = F(ξ) = Bk ξ + bk, 但しBkは3x3の行列、bkは3次元ベクトル
となる。四面体要素の頂点の3次元座標をa0, a1, a2, a3とすると、Bk = (a1-a0, a2-a0, a3-a0)、bk=a0 となる。
・形状関数はベクトルなので、リファレンス要素から座標変換するときにはPiola変換が必要になる。リファレンス要素における形状関数をη(ξ)、望みの四面体要素の形状関数をN(x)とすると
N(x) = (Bk^T)^-1 η(F^-1(x))
となる。
・curlに関する座標変換は
∇×N(x) = 1/det(Bk) Bk ∇×η(F^-1(x))

係数行列の作成
・質量行列 Mij = ∫Ni・Nj dx
・剛性行列 Kij = ∫(∇×Ni)・(∇×Nj) dx
・エッジの向きがあるため、自由度eiにも向きがある。リファレンス要素のエッジの向きと自由度として保持するエッジの向きが逆のときは負号をつける

積分関連
・∫dx = |det(Bk)|∫dξ
・ガウス積分で求める。正確に積分される多項式の次数を2とすると、
三角形：頂点(0,0), (1,0), (0,1)のとき、積分点は(1/2, 0), (0, 1/2), (1/2, 1/2)の3つで、重みは全て1/3×三角形の面積1/2。
四面体：頂点(0,0,0), (1,0,0), (0,1,0), (0,0,1)のとき、積分点は(α, β, β, β), (β, α, β, β), (β, β, α, β), (β, β, β, α)の4つで、重みは全て1/4×四面体の体積1/6。但し、α=0.58541020, β=0.13819660

参考文献
[1] M. Olm, S. Badia and A. F. Martin, Adv. Eng. Softw. 132, 74-91 (2019)
[2] I. Anjam and J. Valdman, Comput. Appl. Math. 267, 252-263 (2015)
[3] Peter Monk, "Finite Element Methods for Maxwell's Equations", Oxford Science Publications, 2003
[4] 小柴正則, 光・波動のための有限要素法の基礎, 森北出版, 1990

間違いなどあればぜひ教えてほしいです。

5月7日23時57分
fartrip