EM Simulation 11 · 境界条件と PML — 無限の空間を有限の箱に収める

Chapter 11

境界条件と PML — 無限の空間を有限の箱に収める

この章のゴール.

「計算領域は有限なのに、なぜ無限空間の答えが出せるのか」に答えられること。 PML の座標伸長という原理を理解し、 なぜ PML が「どんな角度・どんな周波数でも反射しない」のかを説明できること。 そして PML の設定(層数・領域までの距離)を根拠を持って決められるようになること。

この章で使う既出の用語(定義は各リンク先). FDTD(01 章 5 節)、FEM(01 章 5 節)、表面(01 章 5 節)、領域(01 章 5 節)、発散(02 章 2 節)、垂直(03 章 5 節)、磁気的対称面(03 章 4 節)、電気的対称面(03 章 4 節)、EMC(08 章 6 節)

1. 問題 — 箱の壁で波が跳ね返る

FDTD や FEM は空間を格子で刻む。格子は有限だから、どこかで打ち切る必要がある。 何もしなければ、打ち切った面は完全反射の壁になる。

問題 — 打ち切った面は完全反射の壁になるアンテナ跳ね返る何もしなければ壁で全反射吸収PML: 出ていった波は戻ってこない求められているのは「窓」ではなく「無限の彼方」——出た波が二度と戻らない条件である
ABC(吸収境界条件)は垂直入射では良いが斜め入射で悪化する。現代のソルバはほぼ PML を使う

アンテナの放射を計算しているのに、 壁で跳ね返った波が戻ってきて元の場と干渉すれば、答えは完全な嘘になる。

求められているのは「窓」ではなく「無限の彼方」である.

境界の外には無限に広い空間が続いていて、 そこへ出ていった波は二度と戻ってこない——それを有限の箱で再現しなければならない。

歴史的には 2 つのアプローチがあった。

手法中身反射
ABC(吸収境界条件)境界で「外向きに進む波」の方程式を課す(Mur、Engquist–Majda など)垂直入射では良いが、斜め入射で悪化
PML(完全整合層)境界の外側に電磁波を吸収する人工材料の層を貼る角度・周波数によらず極小

現代のソルバはほぼ PML を使う。ABC は軽いので今も併用される。

2. なぜ「ただの吸収材」ではだめなのか

素朴には「損失のある材料(\(\sigma > 0\))を境界に貼れば波が減衰する」と考えたくなる。 しかしこれは失敗する。材料が変われば、その界面で反射が起きるからである。

平面波が媒質 1 から媒質 2 へ垂直入射するときの反射係数は

\[ \Gamma = \frac{\eta_2 - \eta_1}{\eta_2 + \eta_1}, \qquad \eta = \sqrt{\frac{\mu}{\varepsilon}}\quad\text{(波動インピーダンス)} \]

吸収させるために \(\sigma\) を入れる(= \(\varepsilon\) を複素にする)と \(\eta_2 \ne \eta_1\) となり、 入り口で反射してしまう。吸収層に入る前に跳ね返されては意味がない。

3. Bérenger の発明 — 電気と磁気の損失を釣り合わせる

1994 年、Bérenger は決定的なアイデアを出した。

電気的な損失 \(\sigma\) だけでなく、磁気的な損失 \(\sigma^*\) も同時に入れ、 波動インピーダンスが変わらないように釣り合わせる。

条件は

\[ \boxed{\;\frac{\sigma}{\varepsilon} = \frac{\sigma^*}{\mu}\;} \qquad\text{(整合条件)} \]

このとき媒質のインピーダンスは

\[ \eta_{\text{PML}} = \sqrt{\frac{\mu + \sigma^*/(j\omega)}{\varepsilon + \sigma/(j\omega)}} = \sqrt{\frac{\mu\left(1 + \frac{\sigma^*}{j\omega\mu}\right)}{\varepsilon\left(1 + \frac{\sigma}{j\omega\varepsilon}\right)}} = \sqrt{\frac{\mu}{\varepsilon}} = \eta_0 \]

——括弧の中が同じになって約分され、真空と同じインピーダンスのままになる。 つまり入り口で反射せずに入り、中で減衰する材料ができた。

PML の原理 — 電気と磁気の損失を釣り合わせる素朴な吸収材(失敗)真空 η₀σ を入れた材料入り口で反射するη が変わってしまうからPML(Bérenger 1994)真空 η₀σ と σ* を両方反射せずに入り、中で減衰するσ/ε = σ*/μ (整合条件)括弧の中が約分され、η が真空と同じままになる
波動インピーダンス η = √(μ/ε) が変わらなければ界面で反射しない。それを損失を入れつつ実現するのが PML である

4. 座標伸長という見方 — 現代的な定式化

より見通しの良い定式化が座標伸長 PML(stretched-coordinate PML)である。

発想: 空間座標を複素数に「引き伸ばす」。

境界に垂直な方向(\(x\) とする)の微分を

\[ \frac{\partial}{\partial x} \;\longrightarrow\; \frac{1}{s_x}\frac{\partial}{\partial x}, \qquad s_x = \kappa_x + \frac{\sigma_x}{j\omega\varepsilon_0} \]

と置き換える。これはマクスウェル方程式を、複素数に伸びた座標系

\[ \tilde{x} = \int_0^x s_x(x')\,dx' \]

の中で解くことに相当する。

座標伸長 — なぜ角度・周波数に依存しないのか∂/∂x → (1/sₓ) ∂/∂x, sₓ = κₓ + σₓ/(jωε₀)= 空間座標を複素数に「引き伸ばす」PML斜め入射xy伸長は x 方向にだけ掛かる→ 接線方向の波数 k_y は保存される→ 界面での位相整合が完全→ どんな入射角でも反射ゼロこれが "Perfectly Matched"(完全整合)という名前の意味である
整合条件は ω を含まないので周波数にも依存しない。だから広帯域の FDTD でそのまま使える

この座標系では、\(x\) 方向に進む波 \(e^{-jk\tilde{x}}\) は

\[ e^{-jk\tilde{x}} = e^{-jk\int \kappa_x dx}\cdot \underbrace{e^{-\frac{k}{\omega\varepsilon_0}\int \sigma_x dx}}_{\text{指数減衰}} \]

となり、実部が位相の進み、虚部が減衰を与える。

なぜ角度に依存しないのか——これが PML の核心である.

座標伸長は \(x\) 方向にしか掛けていない。 斜めに入射する波は \(x\) 方向成分と \(y\) 方向成分を持つが、 \(y\) 方向は全く手を加えていないので、\(y\) 方向の波数 \(k_y\) は保存される。

界面での反射は「両側で \(k_y\) が一致するか」で決まる(スネルの法則)。 \(k_y\) が同じなら位相整合が完全で、どんな入射角でも反射がゼロになる。 ——これが "Perfectly Matched"(完全整合)という名前の意味である。

さらに周波数にも依存しない。 \(s_x\) の中の \(\sigma_x/(j\omega\varepsilon_0)\) は \(\omega\) を含むが、 これは減衰量を決めるだけで整合条件は \(\omega\) に依らない。 だから広帯域の FDTD でそのまま使える。

5. 現実の PML — 3 つの設定

理論上は完全でも、離散化すると反射が残る。原因と対策は以下の通り。

実務の設定 — 層数と、構造からの距離λ/4 以上PML 8〜12 層PML最も間違えられる設定PML は伝搬する波を吸収する装置である構造のすぐ近くには「リアクティブ近傍界」(蓄えられて放射しない場、15 章)があるこれが PML に触れると、共振周波数と入力インピーダンスが狂う距離は「解析帯域の最低周波数の波長」で判断する1〜10 GHz を見るなら λ は 30 cm(1 GHz のもの)であって 3 cm ではない
σ を滑らかに立ち上げる(多項式分布)、κ>1 でエバネッセント波も吸う、CPML で低周波の安定性を確保する
問題原因対策
界面での離散化反射\(\sigma\) が急に立ち上がると、格子上では不連続に見える\(\sigma\) を滑らかに増やす(多項式分布、通常 3〜4 次)
吸収しきれず裏で反射PML が薄い/\(\sigma\) が小さい層数を増やす/\(\sigma_{\max}\) を上げる
エバネッセント波が減衰しないエバネッセント波=伝搬せずその場で指数関数的に減衰する波(構造のごく近傍にだけ存在する)。これは \(\sigma\) で吸われにくい\(\kappa > 1\)(実部の伸長)を使う
低周波・長時間で不安定PML の実装によっては後期に発散するCFS-PML / CPML(複素周波数シフト付き)を使う

現代の標準は CPML(Convolutional PML)で、\(s_x\) に極をずらす項を加えた形

\[ s_x = \kappa_x + \frac{\sigma_x}{\alpha_x + j\omega\varepsilon_0} \]

を使う。\(\alpha_x > 0\) が低周波での安定性とエバネッセント波の吸収を改善する。

実務の設定値

項目目安
PML の層数8〜12 層(標準)。低周波・強い放射なら 16 層以上
構造から PML までの距離\(\lambda/4\) 以上(推奨 \(\lambda/2\))。最も低い周波数の波長で判断する
理論反射係数 \(R(0)\)\(10^{-5}\sim10^{-8}\) を目標に \(\sigma_{\max}\) を決める(ソルバが自動計算することが多い)

「PML までの距離」が最も間違えられる設定である.

PML はあくまで伝搬する波を吸収する装置である。 構造のすぐ近くにはリアクティブ近傍界(15 章)——蓄えられて放射しない場——が存在し、 これが PML に触れるとエネルギーを吸い取られて、共振周波数や入力インピーダンスが狂う。

だから「構造から \(\lambda/4\) 以上離す」。 そして波長は解析帯域の最低周波数で計算する—— 1〜10 GHz を見るなら \(\lambda\) は 30 cm(1 GHz のもの)であって 3 cm ではない。 低い周波数を含む解析ほど、箱を大きく取らねばならない。

PML の効きを見る

2 次元 FDTD をブラウザで走らせる。 点源から波を出し、 境界を PEC(完全反射)にした場合と PML にした場合を切り替えられる。 壁で跳ね返った波が計算領域を汚す様子と、PML がそれを消す様子を比較してほしい。

6. 境界条件の一覧 — 使い分け

境界意味使いどころ
PML / 放射境界無限に開いた空間アンテナ、散乱、EMC
PEC完全導体の壁導波管、筐体内部、電気的対称面(03 章)
PMC完全磁気導体磁気的対称面
周期境界単位セルの繰り返し周期構造(アレイアンテナ、メタマテリアル=自然界にない実効的な誘電率・透磁率を人工構造で作った材料、FSS=周波数選択表面。特定の周波数だけ通す/反射する金属パターンの繰り返し)
インピーダンス境界有限導電率の表面厚い金属の代替(03 章)

周期境界という強力な道具

無限に並ぶアレイやメタ表面は、単位セル 1 個だけを解けばよい。 境界で場が位相差 \(e^{-jk_x d}\) を持って接続するという条件(Floquet 境界条件)を課す。

1000 素子のアレイを 1 素子の計算で扱えるので、計算量が桁違いに減る。 ただし無限周期の仮定なので、実際の有限アレイの端(エッジ効果)は再現できない。 中央付近の素子の振る舞いを知る道具として使う。

7. この章のまとめ

ポイント内容
問題有限の箱の壁が完全反射する。無限空間を再現する必要
素朴な吸収材の失敗インピーダンス不整合で入り口で反射する
PML の原理電気と磁気の損失を釣り合わせる(\(\sigma/\varepsilon = \sigma^*/\mu\))/座標を複素数に伸長する
角度非依存の理由接線方向の波数 \(k_y\) が保存されるので位相整合が完全
実務設定8〜12 層、構造から \(\lambda/4\) 以上(最低周波数の波長で!)、CPML
他の境界PEC/PMC(対称面)、周期境界(アレイを 1 素子で)