前の記事で紹介した電磁界シミュレータの、中身の話です。モーメント法のソルバは、絵が出ること自体は簡単に作れます。難しいのは、出てきた数字が正しいと言えるところまで持っていくことです。この記事では、そのために確かめたことと、途中で見つかった誤りを順に書きます。
何を解いているか
導体表面の電流を未知数とする混合ポテンシャル積分方程式（MPIE）を、隣り合う 2 つのセルにまたがる屋根形の基底関数（rooftop 基底）で離散化しています。周波数 ω での行列は
Z(ω) = j(ωL − Ψ/ω) + Zs(ω)·Gの形になります。L は部分インダクタンス、Ψ は電荷の電位係数、Zs は導体の面インピーダンス（表皮効果を含む）です。L と Ψ は周波数に依存しないので、一度だけ組み立てれば全周波数で使い回せます。この性質は、3 本目の記事の高速化でも鍵になります。
定性テストは定量の誤りを見逃す
最初の版は「直線は通る」「ギャップは遮断する」「ビアで短絡する」「スタブで落ち込む」といった定性的なテストをすべて通っていました。ところが、別の方法で一つずつ数値を突き合わせると、次の誤りが見つかりました。
| 問題 | 内容 |
|---|---|
| 容量が最大 40 % 小さい | 基板の誘電体を「空気と誘電体の平均」で近似していた。薄い基板ほど外れ、全ての共振周波数が約 25 % ずれていた |
| 数値積分の破綻 | 6×6 点の積分では、薄い基板（h = 10 µm）でグリーン関数の鋭いピークを捉えられず、容量の誤差が 83 % |
| 表示の誤り | 複素電流の絶対値を足していたため、矢印が常に同じ向きを指し、定在波の反転が見えなかった |
どれも、S パラメータの形だけを見ていると気づきません。共振が 25 % ずれていても、直線は通り、ギャップは遮断するからです。
直したのは次の 2 点です。
- グリーン関数を厳密に：接地した誘電体基板の静電グリーン関数を、画像の無限級数として厳密に評価する。薄い基板の極限（平行平板の容量）と、界面の極限の両方を同時に満たす
- セル積分を閉形式に：セル面上の積分を、数値積分ではなく角の座標の差分で表す閉形式に置き換える。厳密になったうえに、評価点も 36 点から 4 点に減った
細かく切るべきは導体の縁だけ
導体の表面電流は、縁に近づくほど 1/√d（d は縁からの距離）で大きくなります。この鋭い変化を捉えるには、縁の近くを細かく分割する必要があります。
最初はパターンの全セルを一律に M×M に分割していました。ところがこの特異性が効くのは縁のごく近くだけです。そこで、内部は粗いまま、導体の縁にだけ幾何級数的に細くなる層（e, 3e, 7e…）を敷く方式に変えました。ADS Momentum の Edge Mesh や Sonnet の縁の細分化と同じ考え方です。
1 セル幅の線路の特性インピーダンスを Hammerstad-Jensen の閉形式と比べると、次のようになります。
| 方式 | 未知数 N | Z₀ の誤差 |
|---|---|---|
| 細分化なし | 18 | 16.5 % |
| 一律 3×3 | 254 | 5.2 % |
| 一律 4×4 | 474 | 4.0 % |
| 縁だけ 2 層 | 186 | 0.2 % |
未知数が少ないのに精度は 20 倍以上良くなりました。分割の数が面積ではなく縁の本数で決まるので、大きな構造ほど効果が大きく、17×17 の全面導体では未知数が 5102 から 842 に減りました。「精度を上げるには全体を細かく」という前提そのものが誤りでした。
収束試験で分かったこと
それでも、メッシュの設定（内部分割 M、縁の層の数 L、一番内側の層の幅 e）をどう選ぶかは残ります。そこで、十分に細かい参照メッシュ（M3・L4・e = 0.02Δ）を基準に、設定ごとの S パラメータの最大誤差を測りました。
| 構造 | 設定 | N | 最大誤差 |
|---|---|---|---|
| 既定基板 h = 0.8 mm の直線 | 旧既定 M1・L2・e = 0.12Δ | 195 | 3.8e-2 |
| 縁だけ細かく M1・L4・e = 0.02Δ | 418 | 7.0e-2 | |
| M2・L4・e = 0.02Δ | 809 | 1.3e-2 | |
| M3・L3・e = 0.03Δ | 1013 | 5.2e-3 | |
| 自動 M2・L3・e = 0.03Δ | 624 | 1.2e-2 | |
| 薄い基板 h = 0.2 mm の直線 | 旧既定 | 195 | 8.2e-2 |
| 自動 M2・L2・e = 0.02Δ | 425 | 2.2e-2 |
分かったことは 3 つです。
- 内部分割 M が効く：セル寸法は設計周波数で 20 セル/波長にしてありますが、最高周波数（6 GHz）では約 11 セル/波長しかありません。屋根形の基底では、これで位相の誤差が出ます（M を 1 から 2 にすると約 7 % 変わる）。「M は普通 1 でよい」という以前の説明は誤りでした
- 縁を細かくするだけでは足りない：M1 のまま縁を細かくした設定（N = 418）は、むしろ誤差が増えています。縁の細さ e への収束は √e 程度とゆっくりで、内部の粗さが残っていると縁だけ詰めても意味がありません
- 表皮厚はメッシュの尺度ではない：表皮厚は導体の厚さ方向の電流分布で、面インピーダンスに閉形式で入っています。面内で電流が縁に集まる幅を決めるのは、基板の厚さ h・セル寸法 Δ・導体の厚さ t です
ここから、メッシュを自動で決める規則を作りました（既定でオン）。
- 内部分割
M = ⌈20Δ / λg,min⌉（最高周波数で 20 セル/波長） - 縁で電流が集まる幅
ℓ = min(h, Δ/2)、一番内側の層の幅e = max(ℓ/12, t/2, 0.02Δ)（シート導体の近似なので金属厚の半分より細かくしない） - 層の数 L は
e·(2ᴸ − 1) ≥ ℓ/2となる最小値（最大 4） - 未知数が上限を超えるときは、L → M → e の順に粗くして、その旨を表示する
この規則で、旧既定より誤差が約 3 分の 1 になりました。
表示していた電流が、実動作の電流ではなかった
S パラメータを求めるには、ポート 1 に電圧をかけてポート 2 を短絡した解（Y パラメータの抽出）を使います。最初の版は、この解をそのまま表面電流として描いていました。
ところが、端を短絡した線路は半波長ごとに短絡共振します。そのため特定の周波数で電流が 56〜79 倍に跳ね上がって見えていました。この共振は、50 Ω で終端した S パラメータにはまったく現れません。表示用の電流だけが、実際には起きない状態を描いていたわけです。
周波数を変えて電流を見比べる機能を足したときに、この異常が見えるようになりました。メッシュを変えても共振周波数がほとんど動かない（2.80→2.85 GHz）ことからメッシュの問題ではないと切り分け、励振の条件そのものが原因だと特定しました。
今は 2 つのポートの解を組み合わせて、ポート 1 から入力し、ポート 2 を 50 Ω で終端した電流を描いています。直線では帯域全体で電流の大きさがほぼ一定（最大と最小の比 1.2 倍）になりました。この電流だけから S11 と S21 を計算し直すと、S パラメータと 1e-14 の精度で一致します。
確かめていること
| 項目 | 結果 |
|---|---|
| 電位のグリーン関数 vs スペクトル領域の厳密解 | 0.0〜0.5 % |
| 実効誘電率 vs Hammerstad-Jensen（εr = 2.2〜10） | 1〜3 % |
| 1 セル幅線路の Z₀ vs Hammerstad-Jensen | 0.2 % |
| 相反性（S21 = S12） | 4e-16 |
| 受動性（反射と透過の電力の和が 1 以下） | 違反なし |
| 極端な値（t = 5 mm、σ = 0.01 MS/m、h = 10 µm など 8 通り） | 計算が破綻する例なし |
| 縦の線路と横の線路（90° 回転） | 3e-12 dB で一致 |
既知の答えと突き合わせられる量を、一つずつ定量的に確かめる。地味ですが、これを飛ばすと、それらしい絵が出るだけの道具になります。
次の記事では、この計算をブラウザの中でどう速くしたかを書きます。→ ブラウザで密行列を速く解く — 適応周波数掃引と WebGPU