EM Simulation 07 · FEM — 弱形式にして、任意形状の要素で解く

Chapter 07

FEM — 弱形式にして、任意形状の要素で解く

この章のゴール.

「弱形式」とは何を弱めたのかを説明でき、部分積分で境界条件が自然に入る仕組みを追えること。 なぜ電磁界では節点ではなく辺に未知数を置くのか(辺要素)を理解すること。 そして適応メッシュという FEM 最大の武器が何をしているかを知ること。

この章で使う既出の用語(定義は各リンク先). FDTD(01 章 5 節)、FEM(01 章 5 節)、MoM(01 章 5 節)、次元(01 章 5 節)、領域(01 章 5 節)、PEC(03 章 4 節)、異方性(03 章 6 節)、非線形(03 章 6 節)、階段近似(05 章 6 節)、周波数依存(06 章 3 節)

1. 出発点 — ベクトル波動方程式

FEM(Finite Element Method、有限要素法)は周波数領域で解く。 02 章の 2 式から \(\mathbf{H}\) を消去すると、\(\mathbf{E}\) だけの方程式になる:

\[ \nabla\times\left(\frac{1}{\mu_r}\nabla\times\mathbf{E}\right) - k_0^2\,\varepsilon_r\,\mathbf{E} = -j\omega\mu_0\mathbf{J} \]

ここで \(k_0 = \omega\sqrt{\mu_0\varepsilon_0} = \omega/c\) は自由空間の波数である。 これをベクトル波動方程式またはベクトルヘルムホルツ方程式と呼ぶ。

この式をそのまま解こうとすると、解に2 階微分可能性が要求される。 しかし現実の解は、材料の境界で微分が不連続になる。そのままでは扱えない。

2. 弱形式 — 何を「弱める」のか

そこで発想を変える。「方程式が各点で厳密に成り立つ」ことを要求するのをやめ、 「任意の試験関数と掛けて積分したときに成り立つ」ことだけを要求する。

残差を \(\mathbf{R}\) として、任意の試験関数 \(\mathbf{W}\) に対して

\[ \int_V \mathbf{W}\cdot\mathbf{R}\ dV = 0 \]

を課す。これが重み付き残差法である。試験関数を「どれだけたくさん」取るかで精度が決まる。

弱形式 — 何を「弱める」のか強形式各点で方程式が厳密に成り立つことを要求→ 解に 2 階微分可能性が要る弱める弱形式任意の試験関数と掛けて積分したときに成立∫ W·R dV = 0部分積分(グリーンの定理)を適用する① 微分の階数が下がるE と W が 1 階ずつ微分されるだけになる→ 1 次多項式のような単純な基底が使える② 境界条件が式に現れる面積分に n̂×∇×E ∝ n̂×H が出る→ ポート励振・吸収境界をここに書き込める
「各点で厳密」を「積分で平均的に」に緩める。その副産物として境界条件が自然に入るのが FEM の設計思想である

なぜこれで十分なのか.

「すべての試験関数に対して積分がゼロ」なら、実質的に「\(\mathbf{R}\) はいたるところゼロ」と同じである。 しかし有限個の試験関数に限れば、有限個の条件になる——これが離散化になる。 「各点で厳密」という強い要求を「積分で平均的に」という弱い要求に緩めたので弱形式という。

部分積分が生む 2 つの利益

\(\mathbf{W}\cdot\nabla\times(\frac{1}{\mu_r}\nabla\times\mathbf{E})\) の項にベクトル版の部分積分 (グリーンの定理)を適用すると

\[ \int_V \frac{1}{\mu_r}(\nabla\times\mathbf{W})\cdot(\nabla\times\mathbf{E})\,dV - k_0^2\int_V \varepsilon_r\,\mathbf{W}\cdot\mathbf{E}\,dV = -\oint_S \mathbf{W}\cdot\left(\hat{n}\times\frac{1}{\mu_r}\nabla\times\mathbf{E}\right)dS + \cdots \]

2 つの重要な効果がある。

効果中身
微分の階数が下がる2 階微分が消え、\(\mathbf{E}\) と \(\mathbf{W}\) がそれぞれ 1 階微分されるだけになる。だから 1 次多項式のような単純な関数を基底に使える
境界条件が式に現れる右辺の面積分に \(\hat{n}\times\nabla\times\mathbf{E} \propto \hat{n}\times\mathbf{H}\) が出る。ポート励振や吸収境界を、この面積分に直接書き込める(自然境界条件)

これが FEM の設計思想の核である.

「解きにくい強形式」を「解きやすい弱形式」に翻訳し、 その翻訳の副産物として境界条件が自然に入ってくる。 PEC のように解自体に課す条件(\(\mathbf{E}_{\text{接線}}=0\))は必須境界条件として基底関数側で処理し、 ポートや放射境界は自然境界条件として面積分に入れる——役割分担が明確である。

3. ガラーキン法と行列方程式

未知の場を有限個の基底関数 \(\mathbf{N}_n\) で展開する:

\[ \mathbf{E}(\mathbf{r}) \approx \sum_{n=1}^{N} x_n\,\mathbf{N}_n(\mathbf{r}) \]

そして試験関数に基底関数と同じものを使う(\(\mathbf{W} = \mathbf{N}_m\))。 これをガラーキン法と呼ぶ。弱形式に代入すると、\(N\) 個の未知数に対する \(N\) 本の連立 1 次方程式になる:

\[ \sum_n \left[\underbrace{\int_V \frac{1}{\mu_r}(\nabla\times\mathbf{N}_m)\cdot(\nabla\times\mathbf{N}_n)dV}_{\text{剛性行列 } S_{mn}} - k_0^2\underbrace{\int_V \varepsilon_r\,\mathbf{N}_m\cdot\mathbf{N}_n\,dV}_{\text{質量行列 } T_{mn}}\right] x_n = b_m \]

すなわち

\[ \boxed{\;(\mathbf{S} - k_0^2\mathbf{T})\,\mathbf{x} = \mathbf{b}\;} \]
ガラーキン法 — 疎行列の連立方程式へE(r) ≈ Σ xₙ Nₙ(r)基底関数で展開し、係数を未知数にする試験関数 = 基底関数(ガラーキン)(S − k₀² T) x = b疎行列行列の性質が実務を決める疎(基底は局所的)→ メモリ O(N) / 対称(相反媒質なら)k₀² を含む → 周波数ごとに解き直し。だから高速掃引(有理関数補間)がある
FEM の最大の弱点は「1 周波数 1 回」。共振が鋭い帯域では補間が外れるので、疑わしければ離散点で検算する

行列の性質が FEM の実務を決める。

性質理由帰結
疎(スパース)基底関数は 1 つの要素の周辺だけで非ゼロ。離れた要素同士の積分はゼロメモリが \(O(N)\) で済む。反復解法・スパース直接解法が使える
対称(相反媒質なら)\(S_{mn}=S_{nm}\)記憶量が半分、解法が高速
\(k_0^2\) に依存周波数が変わると行列が変わる周波数ごとに解き直す必要がある(掃引が高価)

FEM の最大の弱点は「1 周波数 1 回」であること.

広帯域の S パラメータが欲しいなら、周波数点の数だけ行列を解く。 だから商用ソルバは高速掃引(fast sweep)を持っている—— 数点で厳密に解き、その結果から有理関数で内挿する(Padé 近似・AWE=漸近波形評価法・モデル次数低減)。 「離散掃引」と「補間掃引」の違いはここにある。 共振が鋭い帯域では補間が外れることがあるので、疑わしければ離散点を増やして検算する。

4. 辺要素 — 電磁界で節点を使ってはいけない理由

構造力学の FEM では、未知数を節点(要素の頂点)に置く。 電磁界で同じことをすると、スプリアス解——物理的に存在しない偽の解——が大量に出る。

原因は 2 つある。

問題中身
発散条件が課せない節点基底は \(\nabla\cdot\mathbf{E}=0\) を自動的には満たさない。満たさない偽の解が答えに混ざる
境界の連続性が違う03 章の通り、材料境界で連続なのは接線成分だけ。節点基底は 3 成分すべてを連続にしてしまう

解決策が辺要素(edge element、Whitney 要素、ベクトル要素)である。

辺要素 — 電磁界で節点を使ってはいけない理由節点要素(使うと失敗する)未知数を頂点に置く辺要素(Whitney 要素)未知数は「辺に沿った電界の線積分」節点だと: 発散条件が課せず偽の解(スプリアス)が出る/3 成分すべてを連続にしてしまう辺要素なら: 接線成分だけ連続(材料境界を正しく表現)・∇·N = 0 が構造的に成立
05 章の Yee 格子も電界を辺に置いていた。「電界は辺に住む」という原理は手法を超えて共通している

未知数を「節点の値」ではなく「辺に沿った電界の線積分」にする。

四面体の 6 本の辺それぞれに 1 つの未知数を割り当て、基底関数を

\[ \mathbf{N}_{ij} = \lambda_i\nabla\lambda_j - \lambda_j\nabla\lambda_i \]

(\(\lambda_i\) は節点 \(i\) の重心座標)と定義する。この関数には美しい性質がある。

05 章の Yee 格子と同じ思想である.

Yee 格子も電界を「辺」に置いていた。 辺要素はそれを任意形状の四面体に一般化したものと言ってよい。 「電界は辺に住む」という原理は、手法を超えて共通している—— 微分幾何でいう微分形式の構造が背後にある。

5. 適応メッシュ — FEM 最大の武器

FDTD では人間がメッシュを設計する。FEM ではソルバが自動で最適なメッシュを作れる。

適応メッシュ — FEM 最大の武器粗いメッシュで解く誤差推定どこが怪しいか怪しい要素だけ細分化解き直すΔS < 閾値?満たせば終了満たさなければ繰り返す04 章の「一律細分化は最も高い」への直接の回答 —— 場が急変する場所だけが自動で細かくなる停止条件(ΔS < 0.02 を 2 回連続など)が、そのままメッシュ収束の判定になる(13 章)ただし「収束した」は「メッシュに対して安定」であって、モデルが正しい保証ではない
FDTD では人間がメッシュを設計するが、FEM ではソルバが自動で最適化できる。これが大きな実務差になる
段処理
1粗いメッシュで解く
2誤差推定: 要素間で場が不連続になっている量、残差の大きさなどから「どこが怪しいか」を推定する
3怪しい要素だけを分割して細かくする(\(h\)-refinement)/次数を上げる(\(p\)-refinement)
4解き直す。前回の解との差(たとえば S パラメータの変化 \(\Delta S\))が閾値以下なら終了
5でなければ 2 へ

04 章で見た「一律細分化は最も高い」という問題への直接の回答である。 場が急変する場所(導体のエッジ、給電点、共振部)だけが自動的に細かくなる。

適応メッシュの停止条件が「収束」の意味を持つ.

多くの商用 FEM ソルバは「\(\Delta S < 0.02\) を 2 回連続で満たしたら終了」のような条件を持つ。 これは 13 章で扱うメッシュ収束そのものであり、 FEM では収束判定が自動化されているのが FDTD との大きな違いである。 ただし「収束した」は「メッシュに対して答えが安定した」だけであって、 モデルが正しいことは保証しない——ここを混同しないこと。

弱形式と基底関数を見る

要素数を変えて、1 次元の問題を基底関数の重ね合わせで近似する様子を見られる。 要素を増やすと真の解に近づく——有限要素法が何をしているかの本質がここにある。

6. FEM の得手不得手

場面適性理由
複雑形状・曲面◎四面体でそのまま貼れる。階段近似がない
高 Q 共振・フィルタ◎周波数領域なので減衰待ちがない
不均質・異方性材料◎要素ごとに材料テンソルを持てる
狭帯域で高精度◎適応メッシュで精度を追い込める
広帯域(多数の周波数点)△1 点ずつ解き直し。高速掃引で緩和
電気的に非常に大きい構造△〜×未知数が爆発。メモリが厳しい
開放領域(放射問題)○吸収境界(11 章)が要る。MoM より不利なことも
時間応答・非線形×定常解を解く枠組み

7. この章のまとめ

ポイント内容
弱形式「各点で厳密」を「積分で成立」に弱める。部分積分で微分階数が下がり、境界条件が自然に入る
ガラーキン試験関数=基底関数。結果は \((\mathbf{S}-k_0^2\mathbf{T})\mathbf{x}=\mathbf{b}\) の疎行列
周波数依存行列が \(k_0^2\) を含む → 周波数ごとに解き直し。高速掃引で緩和
辺要素未知数は辺の線積分。接線連続・発散ゼロでスプリアス解を防ぐ
適応メッシュ誤差推定 → 局所細分化の自動ループ。FEM 最大の武器
適性曲面・共振・不均質に強く、広帯域・超大規模に弱い