3Dで学ぶ電磁気学

G-1 有限要素法の基礎

有限要素法は、偏微分方程式をいきなり巨大な行列へ置き換える手法ではありません。解析領域を小さな要素へ分け、各要素で場を補間し、弱形式から得た局所的な関係を全体へ組み立て、境界条件を入れて解く方法です。このページでは、高橋則雄先生の有限要素法の体系に沿って、磁気ベクトルポテンシャル $A_z$ を未知量とする静磁界問題を入口にします。

このページでつなぐ四つの段階

1 分割と近似連続な場を三角形ごとの一次関数で近似し、要素を細かくしたときの誤差を観察します。
2 補間関数面積座標 $N_1,N_2,N_3$ が節点値から要素内部の値を作る様子を動かします。
3 組立と境界要素係数行列を全体行列へ加算し、固定境界と自然境界の違いを確かめます。
4 未知量の置き場所四面体の節点自由度と辺自由度を比べ、電磁界で辺要素が必要になる理由へ進みます。

可視化 G-1-1 メッシュを細かくすると何が改善するか

単位正方形上の既知の場

\[ A_z(x,y)=\sin(\pi x)\sin(\pi y) \]

を三角形一次要素で補間します。この場が一定の磁気抵抗率 $\nu$ のPoisson型方程式を満たすとき、右辺は

\[ -\nabla\cdot(\nu\nabla A_z)=J_z, \qquad J_z=2\pi^2\nu\sin(\pi x)\sin(\pi y) \]

です。ここでは解法誤差を混ぜず、節点上では正しい値を与え、要素内部を一次補間したときの誤差だけを見ます。分割数を増やし、折れた近似面が滑らかな場へ近づくことと、計算量が増えることを同時に確認してください。

図を動かして動かしてください。$h=1/n$ を半分にすると、滑らかな解に対する一次補間の $L^2$ 誤差はおおむね $O(h^2)$、勾配の誤差はおおむね $O(h)$ で減ります。

\[ \lVert A_z-I_hA_z\rVert_{L^2(\Omega)}=O(h^2), \qquad \lVert \nabla(A_z-I_hA_z)\rVert_{L^2(\Omega)}=O(h) \]

可視化 G-1-2 一次三角形要素の中を面積座標で補間する

三角形の節点を $P_1,P_2,P_3$、節点値を $a_1,a_2,a_3$ とすると、要素内部の近似値は

\[ A_h(x,y)=N_1(x,y)a_1+N_2(x,y)a_2+N_3(x,y)a_3 \]

で与えます。一次三角形では $N_i$ は面積座標であり、

\[ N_i(P_j)=\delta_{ij}, \qquad N_1+N_2+N_3=1 \]

を満たします。評価点 $P$ をドラッグするか、$x,y$ を動かしてください。色付きの面積比が $N_1,N_2,N_3$ となり、節点値の重み付き平均が $A_h(P)$ になります。

三角形内部で $\nabla N_i$ は一定です。したがって一次要素の $A_h$ は連続でも、その勾配から求める磁束密度

\[ \vec B_h=\nabla\times(A_h\vec e_z) =\dfrac{\partial A_h}{\partial y}\vec e_x -\dfrac{\partial A_h}{\partial x}\vec e_y \]

は要素ごとに一定となり、要素境界で跳びます。細分によって磁束密度分布が滑らかになる理由はここにあります。

可視化 G-1-3 要素係数行列を全体へ組み立てる

磁気抵抗率 $\nu$ が要素内で一定のとき、一次三角形要素の係数行列は

\[ K_{ij}^{(e)} =\int_{\Omega_e}\nu\nabla N_i\cdot\nabla N_j\,d\Omega =\dfrac{\nu}{4\Delta_e}(b_ib_j+c_ic_j) \]

です。局所節点番号を全体節点番号へ写すBoolean行列 $\boldsymbol L_e$ を使えば、全体行列は

\[ \boldsymbol K =\sum_e\boldsymbol L_e^{\mathsf T}\boldsymbol K^{(e)}\boldsymbol L_e, \qquad \boldsymbol K\vec a=\vec f \]

となります。スライダーで要素を一枚ずつ組み立て、同じ全体節点へ複数要素の寄与が足し込まれる様子を見てください。4要素を組み立てた自然境界の行列では、各行の和が0となり、定数ベクトルが零モードです。

自然境界だけでは定数を加えても勾配が変わらないため、行列に零空間が残ります。固定境界を入れると、既知の境界値 $a_1$―$a_4$ は自由な未知量から除かれます。この例で残る自由系は中心節点だけです。全要素を組み立て、中心へ単位荷重 $f_5=1$ を与えると、

\[ \boldsymbol K_{ff}=[K_{55}]=[4], \qquad a_5=K_{55}^{-1}f_5=\dfrac14 \]

となります。確認を確定した後、自然境界と固定境界を切り替え、零モードを持つ5自由度系が正則な1自由度系へ変わることを行列の表示で確かめてください。

\[ \boldsymbol K\vec 1=\vec 0 \quad\text{(自然境界のみ)}, \qquad A_z\vert_{\Gamma_D}=0 \quad\text{(固定境界)} \]

可視化 G-1-4 四面体では節点と辺のどちらへ未知量を置くか

スカラーポテンシャルには節点要素が自然です。一方、磁気ベクトルポテンシャル $\vec A$ や電界 $\vec E$ のように接線連続性を持つ1-formは、辺に沿う積分値を自由度とする辺要素で表します。四面体の重心座標を $\lambda_i$ とすると、辺 $i\to j$ の一次Nédélec基底は

\[ \vec W_{ij}=\lambda_i\nabla\lambda_j-\lambda_j\nabla\lambda_i \]

です。この基底は選んだ辺に沿う線積分だけを1にします。

\[ \int_{e_{kl}}\vec W_{ij}\cdot d\vec\ell =\delta_{ik}\delta_{jl}-\delta_{il}\delta_{jk} \]

表示を節点要素と辺要素で切り替え、辺要素では選択した辺を変えてください。矢印は選択した辺基底 $\vec W_{ij}$ の要素内部での向きです。

なぜ電磁界解析では辺要素がよいのか

静磁界の磁気ベクトルポテンシャル法では、$\vec B=\nabla\times\vec A$ を用いて、代表的には次の弱形式を解きます。

\[ \int_{\Omega}\nu (\nabla\times\vec A)\cdot(\nabla\times\vec v)\,dV =\int_{\Omega}\vec J\cdot\vec v\,dV, \qquad \vec A,\vec v\in H(\operatorname{curl};\Omega) \]

この問題が要求するのは、要素境界をまたぐ接線成分の整合です。三成分をすべて節点値として連続にするベクトル節点要素は、必要以上に強い連続性を課します。Nédélec辺要素は接線成分を連続に保ちながら、法線成分には材料境界やゲージに応じた不連続を許せるため、curlを含むMaxwell問題の関数空間に適合します。

辺要素の未知量は、点でのベクトル値ではなく、向きを付けた辺に沿う循環です。

\[ a_e=\int_e\vec A\cdot d\vec\ell, \qquad \vec A_h=\sum_{e}a_e\vec W_e \]

したがって、隣り合う四面体は共有辺の自由度を同じ一個の未知量として使えます。辺の向きを逆にすれば $a_e$ の符号も反転するため、Stokesの定理で面を回る循環を足し合わせる操作が、行列の符号としてそのまま組み立てられます。

\[ \int_f(\nabla\times\vec A_h)\cdot\vec n\,dS =\oint_{\partial f}\vec A_h\cdot d\vec\ell =\sum_{e\subset\partial f}C_{fe}a_e \]

さらに、節点から辺、辺から面へ進む離散微分は、連続場の恒等式を保ちます。

\[ H^1\xrightarrow{\nabla}H(\operatorname{curl}) \xrightarrow{\nabla\times}H(\operatorname{div}), \qquad \boldsymbol C\boldsymbol G=\boldsymbol0, \qquad \boldsymbol D\boldsymbol C=\boldsymbol0 \]

このため、勾配場をcurlした成分や、curlを発散した成分が離散化だけを原因として現れにくくなります。特にMaxwell固有値問題では、適合するメッシュ・空間・組立てと組み合わせることで、物理的でない偽モードを抑える重要な条件になります。

| 比べる点 | ベクトル節点要素 | Nédélec辺要素 |

|---|---|---|

| 自由度 | 節点での各成分 | 向き付き辺の線積分 |

| 要素間で連続にする量 | ベクトルの全成分 | 接線成分 |

| 適合する代表空間 | $(H^1)^3$ | $H(\operatorname{curl})$ |

| 得意な未知量 | 変位など全成分を連続にしたいベクトル | $\vec E$、$\vec A$ などcurlを使う場 |

| Maxwell構造 | 後から整合を工夫する必要がある | Stokes則と離散複体を直接組み込める |

ただし、辺要素を選ぶだけで解析が自動的に正しくなるわけではありません。磁気ベクトルポテンシャルには $\vec A$ と $\vec A+\nabla\phi$ が同じ $\vec B$ を与えるゲージ自由度が残るため、ゲージ条件、木・余木分解、あるいは適切な正則化が必要です。穴のある領域では切断面やcohomology自由度も別に扱います。逆に、未知量がスカラーポテンシャルのような0-formなら節点要素、磁束のような2-formなら面要素が自然です。大切なのは「辺要素が常に優秀」ではなく、Maxwell方程式が積分する幾何学的対象と自由度の置き場所を一致させることです。

次に進む電磁界数値解析

弱形式部分積分、境界項、Galerkin法から有限要素方程式を作ります。 離散Hodgeと材料則接続行列と材料行列を分け、辺・面自由度を組み立てます。 電磁力と節点力Maxwell応力、エネルギー微分、節点力を比較します。 位相・切断面・ゲージ穴を回る大域自由度と零空間を扱います。 保証付き誤差評価相補場からエネルギー誤差の上下界を作ります。 高次基底と適応解析h細分・p次数上昇・hp適応を選びます。 開境界とDtN作用素無限外部を境界作用素へ縮約し、弱形式へ戻します。

主な参照資料