カール要素とは (そしてなぜそれが必要なのか) ?

2019年 12月 30日

カール要素はエッジ要素またはベクトル要素と呼ばれることもあり, 電磁気の問題を求解するために有限要素法で広く使用されています. このブログでは, COMSOL Multiphysics® ソフトウェアで使用される理由と方法など, このタイプの要素を包括的に紹介しています. ただし, カール要素が使用される理由と方法を理解するのは簡単ではありません. したがって, 最初にマックスウェル方程式と有限要素法に関する背景情報を確認します.

2つの電磁気学で解くべき方程式の基本形

電磁気学モデリングの目的は,特定の境界条件に従うマックスウェル方程式を解くことです. 微分形式のマックスウェル方程式は次のように与えられます:

(1)

\nabla \cdot \textbf{D} = \rho

 

(2)

\nabla \times \textbf{E}=-\frac {\partial \textbf{B}}{\partial t}

 

(3)

\nabla \cdot \textbf{B} = 0

 

(4)

\nabla \times \textbf{H} = \textbf{J} + \frac {\partial \textbf{D}}{\partial t}

 

ここで, \textbf{E}\textbf{H} はそれぞれ電場と磁場です; \textbf{D}\textbf{B} は, それぞれ電束密度と磁束密度です; {\rho}\textbf{J} は, それぞれ電荷密度と電気伝導電流密度です.

系を閉じるために, マックスウェル方程式には, 媒体の巨視的特性を表す構成関係が含まれています. COMSOL Multiphysics を使用すると, 非線形および異方性材料を含む, さまざまな媒体をモデル化できます. ただし, 簡単にするために, 残留効果を無視し, 媒体が線形で等方性であると仮定します. その結果, 次のような単純な構成関係が得られます:

(5)

\textbf{D} = \epsilon \textbf{E}

\textbf{B} = \mu \textbf{H}

 

ここで, \epsilonは誘電率, \muは媒体の透磁率です.

(1)(4) は, 連続媒体中の点で成り立ちますが, 媒体1と2の間の不連続界面では, 境界条件 は次のように表されます.

(6)

\textbf{n}_2 \cdot (\textbf{D}_1 – \textbf{D}_2) = \rho_s

 

(7)

\textbf{n}_2 \times (\textbf{E}_1 – \textbf{E}_2) = \textbf{0}

 

(8)

\textbf{n}_2 \cdot (\textbf{B}_1 – \textbf{B}_2) = 0

 

(9)

\textbf{n}_2 \times (\textbf{H}_1 – \textbf{H}_2) = \textbf{J}_s

 

ここで, \bf{n_2}は媒体2からの外向きの法線であり, {\rho_s}\textbf{J}_sはそれぞれ表面電荷密度と表面電流密度です.

(7)(8) は, 接線方向の電場と磁束場の法線成分が境界で連続であることを意味します. 式 (6)(9) は, 電束場と接線磁場の法線成分が不連続になる可能性があることを意味します.

特定の問題については, マックスウェル方程式から導かれる別の形の電気, 磁気,または電磁気の定式化を使うのが便利です. これらの定式化に基づいて, COMSOL Multiphysics には, いくつかの電磁気モジュール (AC/DC モジュール, RF モジュール, 波動光学モジュールなど) と, モジュールごとにいくつかの異なるフィジックスインターフェースがあります.

たとえば, AC/DC モジュールの静電インターフェースは, 静電荷のみが存在する静電問題をモデル化します. この状況では, 方程式 (1)(2) を解くだけで済みます. 電位Vを導入して\textbf{E}=-\nabla Vを定義することにより, 方程式 (1)(2) は, \nabla \times \nabla \textit{V} = 0が常に成り立つため, 単一の方程式であるポアソン方程式に還元されます. ポアソン方程式は次のように表されます:

(10)

\nabla \cdot (\epsilon \nabla \textit{V} )= – \rho

COMSOL Multiphysics では, すべての電位スカラーポテンシャルインターフェースがポアソン方程式のバージョンとして定式化されています. 電流のない磁場の場合, 式(4)は静電場と同様の形になります. 磁気スカラーポテンシャルを導入することにより, 問題をポアソン方程式に定式化することもできます. その他の場合は, ベクトルの定式化を処理する必要があります. たとえば, 静電流のみが存在する静磁場の状況では, 方程式 (3)(4) を解くだけで済みます. 磁気ベクトルポテンシャル\textbf{A}を導入し, \textbf{B}=\nabla \times \textbf{A}を定義することにより, 式 (3)(4) は, \nabla \cdot (\nabla \times \textbf{A} )= 0が常に成り立つため, 以下の1つの式に還元されます.

(11)

\nabla \times (\frac {1}{\mu} \nabla \times \textbf{A} )= \bf{J}

(11) は, \nabla \times (\nabla \times .)演算子を含むベクトル式の形式で記述されます. 実際, 電磁波などの他の状況もこの形式で定式化されます.

議論を容易にするために, 例として時間調和電場を取り上げます. 周波数領域では, 式 (2) を式 (4) に挿入すると次のようになります.

(12)

\nabla \times (\frac {1}{\mu} \nabla \times \textbf{E} ) – \omega^2 \epsilon \textbf{E}= – j \omega \textbf{J}

ここで, \omegaは周波数, jは虚数単位です.

形状関数とラグランジュ要素

COMSOL Multiphysics では, 偏微分方程式 (PDE) を解くために有限要素法 (FEM) が一般的に使用され, マックスウェル方程式も例外ではありません. 有限要素法は, 次のようないくつかのステップで偏微分方程式を解くために使用されます:

  1. ドメインをメッシュ要素と呼ばれる多くの小さな重複しない要素に分割
  2. 解を局所形状関数または基底関数で各要素上で近似
  3. PDEを弱形式で記述し, 各メッシュ要素を離散化してローカル行列を取得
  4. ローカル行列をグローバル行列にアセンブルして求解

FEM がどのように機能するかを示す例として, ポアソン方程式 (10) を取り上げましょう. 使用する形状関数に応じて, さまざまな要素を使用できます. 簡単にするために, ドメインは 2D であると見なし, 線形三角形のラグランジュ要素を使用します. 頂点i,j,kを持つ要素の電位Vは, 線形関数N_{i,j,k}(x,y)で次のように近似できます.

(13)

V(x,y) =

\begin{bmatrix}

V_i, V_j, V_k

\end{bmatrix}

\begin{bmatrix}

N_i(x,y) \\

N_j(x,y) \\

N_k(x,y)

\end{bmatrix}

= V_i N_i(x,y)+V_j N_j(x,y) + V_k N_k(x,y)

ここから, 各頂点が自由度 (DOF) を追加し, 対応する形状関数が頂点で1に等しく, 他のすべての頂点では0であることがわかります. 高次のラグランジュ要素は, 頂点だけでなく, 要素内にある他のノードにも DOF を追加します. したがって, ラグランジュ要素は節点要素とも呼ばれます. 1次三角形のラグランジュ要素の形状関数を以下に可視化します.

1次三角形のラグランジュ要素の形状関数の図.

次の図に示すように, ラグランジュ要素は各要素で連続しているだけでなく, 境界を越えて連続しています.

境界を越えて連続しているラグランジュ要素の図.

媒体1と媒体2は同じ境界kiを共有します. 境界kiに近い電場は, 媒体1と2にそれぞれ青と赤の矢印でプロットされます. 電位Vは境界を越えて連続的であるため, Vの接線方向の勾配 (つまり, 接線方向の電場) は連続的です. 言い換えれば, 境界条件式 (7) は自動的に満たされます. したがって, 式を考慮するだけで済みます. 式(6) 表面電荷が存在する場合. 次に, これを FEM で自然に処理する方法を示します.

ポアソン方程式(10)の弱形式の導出は難しくありません. テスト関数\phi を式(10)の両側に乗算すると, 定義域\Omegaの積分により次のようになります.

(14)

\int_{\Omega} \nabla \cdot (\epsilon \nabla \textit{V} )\phi= \int_{\Omega} – \rho \phi

左側の部分積分を適用すると, 次のようになります.

(15)

-\int_{\Omega} \epsilon \nabla \textit{V} \cdot \nabla \phi + \int_{\partial \Omega} \epsilon \nabla \textit{V} \cdot \textbf{n} \phi= \int_{\Omega} – \rho \phi

\partial \Omegaはドメインの境界であり, \textbf{n}はドメインから離れる方向に向いた法線ベクトルです.

このステップは非常に重要であり, 2つの主な利点があります:

  1. 空間導関数の最大次数を減らします. これにより, 線形 (1次) 要素を使用して2次偏微分方程式を解くことができます
  2. 方程式の自然境界条件 (考慮しない場合は自動的に満たされる) が何であるかを明確にする

左側の2番目の積分は, \nabla \textit{V}の法線成分が境界で消えると消えます. それがどのように機能するかを明確に理解するには, 式 (15) を次のように書き直すと便利です.

(16)

\int_{\Omega} \textbf{D} \cdot \nabla \phi – \int_{\partial \Omega} \textbf{D} \cdot \textbf{n} \phi= \int_{\Omega} – \rho \phi

ここで, 自然境界条件が\textbf{D} \cdot \textbf{n}=0であることが明らかです.

AC/DC モジュールの静電インターフェースでは, この自然境界条件がデフォルトとして設定されています. 境界条件が電位ではなく電荷密度で与えられる場合, 式 (16) に境界条件式(6) を組み込むのは簡単です.

カール要素

ラグランジュ要素ベースの FEM を使用して, ベクトル方程式 (12) をさまざまな成分\textbf{E}_{x,y,z}に分割することにより, ベクトル方程式 (12) を解くこともできます. ただし, 特に不規則な境界の場合, 境界条件の実装は複雑です. それでも, ラグランジュ要素は各成分\textbf{E}_{x,y,z}が境界を越えて連続するように強制するため, 結果は誤っている可能性があります. これは, 電場の通常の成分が材料の境界面, 特に鋭い角からの境界で不連続になる可能性があるという事実と反します.

例として, アプリケーションライブラリのSierpinski フラクタルモノポールアンテナモデルを取り上げましょう. 1.6GHz でのフラクタルラジエーターの表面電場を以下に示します. \textbf{E}_{x,y,z}がエッジ間で大きく変化することは明らかです.

スカラー場を近似するためにスカラー形状関数を使用した, COMSOL Multiphysics の ラグランジュ要素の図.

前に示したように , ラグランジュ要素はスカラー形状関数を使用してスカラー場を近似します. ベクトル場を近似するためにベクトル形状関数を使用するのは自然なことです. たとえば, 式 (12) の電場は次のように表すことができます.

(17)

\textbf{E} = E_i \textbf{W}_i+E_j \textbf{W}_j + E_k \textbf{W}_k

(17) が示すように, ベクトル場にも方向があるため, ノードに DOF を追加することは適切ではありません. エッジに沿ったベクトル場成分 (接線成分) はスカラーであるため, スカラー値を取得するには, 各エッジに DOF を追加できます. 次に, 式 (17) を次のように書き直します.

(18)

\textbf{E} = E_{ij} \textbf{W}_{ij}+E_{jk} \textbf{W}_{jk} + E_{ki} \textbf{W}_{ki}

ラグランジュ形状関数と同様に, 1つのエッジに沿ったベクトル形状関数の接線成分はゼロ以外の定数であり, 他のエッジに沿ってはゼロに等しくなります. 上記の特性を満たす形状関数の1つのタイプは, Whitney 1形式の基底関数 (参照1, 参照2) であり, 次のように表されます.

(19)

\textbf{W}_{ij} = N_i \nabla N_j – N_j \nabla N_i

\textbf{W}_{jk} = N_j \nabla N_k – N_k \nabla N_j

\textbf{W}_{ki} = N_k \nabla N_i – N_i \nabla N_k

1次三角エッジ要素の形状関数を以下にプロットします.

1次三角エッジ要素の形状関数のプロットの図.

数学の観点から, エッジ要素はH(curl)空間に適合しているため (参照3), カール要素とも呼ばれます. COMSOL Multiphysics では, より一般的な名前カール要素を使用します. これは, 高次の要素がエッジだけでなく要素の内部にも DOF を追加するためです. 2つの隣接するラグランジュ要素に対してプロットしたものと同様に, 同じ境界kiを共有する2つのカール要素の基底関数を以下に示します.

同じ境界を共有する2つのカール要素の関数のプロットの図.

境界を横切る接線方向の電場が連続していることがわかります. これは, ラグランジュ要素を使用してポアソン方程式を解く場合と非常によく似ています. したがって, カール要素を使用することにより, 境界条件式は次のようになります. 式(7) は自動的に実行されます. 次に, 境界方程式がどのようになっているのかを示しましょう. 式(9) は FEM によって処理されます.

ポアソン方程式\textbf{W}を解くのと同様に, テスト関数を方程式 (12) の両側に乗算すると, 定義域の積分により次のようになります.

(20)

\int_{ \Omega} \nabla \times (\frac {1}{\mu} \nabla \times \textbf{E} ) \textbf{W}=\int_{ \Omega} (\omega^2 \epsilon \textbf{E} – j \omega \textbf{J}) \textbf{W}

左側の部分積分を適用すると, 次のようになります.

(21)

\int_{ \Omega} (\frac {1}{\mu} \nabla \times \textbf{E} ) (\nabla \times \textbf{W}) + \int_{\partial \Omega} \textbf{n} \times (\frac {1}{\mu} \nabla \times \textbf{E} )\textbf{W} =\int_{ \Omega} (\omega^2 \epsilon \textbf{E} – j \omega \textbf{J}) \textbf{W}

左側の2番目の部分を書き換えると, 次のようになります.

(21)

\int_{ \Omega} (\frac {1}{\mu} \nabla \times \textbf{E} ) (\nabla \times \textbf{W}) + \int_{\partial \Omega} \textbf{n} \times (\frac {1}{\mu} \nabla \times \textbf{E} )\textbf{W} =\int_{ \Omega} (\omega^2 \epsilon \textbf{E} – j \omega \textbf{J}) \textbf{W}

ここで, 境界条件式(9) 弱式に簡単に組み込むことができます.

カール要素の長所と短所

カール要素を使用することにより, マックスウェル方程式の境界条件を自然に処理できることを示しました. カール要素は, 接線方向の電場が連続するように強制し, 通常の電場が境界を越えてジャンプできるようにします. さらに, 他の境界条件の制約も容易になります (参照4).

たとえば, 電磁波 (周波数領域)インターフェースでは, 金属などの導電性表面をモデル化するために, 完全電気導体境界機能がデフォルトとして設定されています. 境界条件により, 接線方向の電場がゼロになります. つまり, \textbf{n} \times \textbf{E} = \textbf{0} となります. 磁場インターフェースでは, デフォルトの境界条件は磁気絶縁です. これは, 接線方向の磁気ベクトルポテンシャルをゼロに制限します. つまり, \textbf{n} \times \textbf{E} = \textbf{0} です. 未知数は境界上の接線場であるため, これらの境界条件は, 点ごとの制約を使用して FEM で簡単に考慮することができます. 特に電磁波の問題を求解するために, スプリアス解を排除するなど, カール要素を使用することには他にも利点があります (参照5). ただし, 境界条件の自然な処理が支配的である必要があります.

カール要素を使用することにはいくつかの欠点があります. たとえば, 線形カール要素には, 次数O(h)の局所近似誤差があります. ここでhは要素サイズですが, 1次ラグランジュ要素を使用すると, 局所誤差はO(h^2)に収束します (参照6). これは, カール要素が混合要素であり, 形状関数の順序がさまざまな方向で変化するためです. たとえば, エッジijおよびエッジijに平行な方向に沿った形状関数\textbf{W}_{ij}の接線成分は定数ですが, 他の方向では線形関数です. 混合特性は成分を確認することによっても表示できます. たとえば,

(23)

\textbf{W}_{ij} = N_i \nabla N_j – N_j \nabla N_i

=(a_i+b_i x+c_i y)(b_j, c_j)-(a_j+b_j x+c_j y)(b_i,c_i)

=(d_i + d_{ij} y, d_j – d_{ij} x)

ここで, a,b,c,dは要素の座標のみに依存する定数です.

\textbf{W}_{ij}の成分は, x軸に沿って一定であり, y軸に沿って線形です. 各成分の空間導関数の精度は大幅に異なります. このため, カール要素をポスト処理する場合, 場の高次の空間導関数は使用できません. 式 (23) は, 線形カール要素の形状関数が特定の方向に沿って一定になり得ることも示しています. これにより, 局所誤差の収束が1次ラグランジュ要素を使用する場合よりも遅くなります. 一方では, この欠点は高次のカール要素またはより細かいメッシュを使用することで解消できます. 他方, 境界条件を自然に処理することの利点と比較して, 式 (12) を解くためにラグランジュ要素を使用することによって遭遇する困難は, 高階形状関数を使用したりメッシュを改良したりすることによって補うことができないため, 欠点ははるかに重要ではなくなります.

概要

このブログでは, マックスウェル方程式とその境界条件から始めて, 電磁気学モデリングでよく見られる2つの基本的なタイプの方程式を紹介しました. つまり, スカラーポアソン方程式と\nabla \times (\nabla \times .)演算子を使用したベクトル方程式です. 古典的なラグランジュ要素を使用してポアソン方程式を解くことにより, 連続接線電場の条件が自然に満たされることを示しました. ただし, 関連する境界条件の処理が難しいため, ラグランジュ要素を使用してベクトル方程式を解くのは適切ではありません.

次に, 連続的な接線電場の条件を満たすことができるカール要素を導入しました. 同時に, 弱式を詳細に導出することにより, 他の境界条件が FEM にどのように組み込まれるかを示しました. 最後に, カール要素を使用することのいくつかの欠点について説明しました. これは, 利点よりもはるかに重要ではありません

参考文献

  1. H. Whitney, Geometric integration theory, Princeton UP, Princeton, 1957.
  2. A. Nentchev, Numerical analysis and simulation in microelectronics by vector finite elements, 2008.
  3. J.C. Nédélec, Mixed finite elements in ℝ3. Numerische Mathematik, 35(3), pp. 315–341, 1980.
  4. J.P. Webb, Edge elements and what they can do for you. IEEE Transactions on Magnetics, 29(2), pp. 1460–1465, 1993.
  5. J.M. Jin, Theory and computation of electromagnetic fields. John Wiley & Sons, 2011.
  6. G. Mur, Edge elements, their advantages and their disadvantages. IEEE transactions on magnetics, 30(5), pp. 3552–3557, 1994.

コメント (0)

コメントを残す
ログイン | 登録
Loading...
COMSOL ブログを探索