電磁場解析SUPERFISH二次元軸対称TMモード

解いている偏微分方程式

本プログラムは二次元軸対称（r,z）座標系での時間調和電磁場（$ e^{-i\omega t}$）のTMモードを扱う．円筒座標系において，軸方向成分$E_z(r,z)$を主変数として扱うと，真空中（比例定数$\mu,\epsilon$）での支配方程式はスカラーのHelmholtz方程式に帰着する：

$$ \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial E_z}{\partial r}\right)+\frac{\partial^2 E_z}{\partial z^2}+k^2 E_z=0,\qquad k^2=\omega^2\mu\epsilon. $$

この方程式は導出上，Maxwell方程式（無電荷・無電流）から始まり，$\partial/\partial\phi=0$（軸対称）かつTMモード（磁場が周方向成分を持ち，電場の主要成分は$E_z$）を仮定して得られる．導出の概略は次の通りである：

1. Maxwell方程式より，時間ハーモニック場での回転と発散の関係を用いる．
2. 軸対称かつTMモードの仮定の下で，$E_z$の周方向微分が消え，他の成分は$E_z$から微分で表される．
3. これらを代入して整理すると，上のHelmholtz方程式が得られる．

境界条件：導体表面では伸長電場が消えるため

$$ E_z=0\quad\text{(導体表面)} $$

開放境界や対称軸($r=0$)では適切な対称条件や吸収境界条件（PML 等）を用いる．

有限要素法による式

Helmholtz方程式を有限要素法で近似するため，まず重み付き残差法（変分法）により弱形式を導く．試験関数を$v(r,z)$とすると，円筒座標の体積要素は$r\,dr\,dz$であるため，弱形式は次の形になる：

$$ \int_\Omega \left[\left(\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial E_z}{\partial r}\right)+\frac{\partial^2 E_z}{\partial z^2}\right)v + k^2 E_z v\right] r\,dr\,dz =0. $$

部分積分を行い，境界での自然境界項を取り除くと（Dirichlet境界で$E_z=0$として），弱形式はより標準的な形になる：

$$ \int_\Omega \left(\frac{\partial v}{\partial r}\frac{\partial E_z}{\partial r}+\frac{\partial v}{\partial z}\frac{\partial E_z}{\partial z}\right) r\,dr\,dz -\int_{\partial\Omega}\! v\,\frac{\partial E_z}{\partial n} r\,ds + k^2\int_\Omega v E_z r\,dr\,dz =0. $$

Dirichlet 境界（導体）では第2項は消える．これを有限要素分解で近似する：局所基底関数$N_i$の線形結合で$E_z\approx\sum_j N_j c_j$，試験関数にそれらを取ると離散化により一般化固有値問題が得られる：

$$ K\mathbf{c}=\omega^2 M\mathbf{c}, $$

ここで剛性行列$K$と質量行列$M$は要素毎の積分で組み立てられる：

$$ K_e=\int_{\Omega_e} (B^T B)\, r\,dA,\qquad M_e=\int_{\Omega_e} (N^T N)\, r\,dA, $$

ただし，$B$は基底関数の勾配（$\partial/\partial r,\partial/\partial z$を並べた行列），$N$は基底関数ベクトルである．数値積分（ガウス積分）を用いて各要素で評価し，グローバル行列へ組み立てる．

要点：

• 円筒座標由来の重み$r$を積分に含めること（体積要素に注意）。
• 要素レベルでの線形形関数（3角形・4節点）や高次形関数を用いれば精度向上。
• 境界条件は強制（Dirichlet）あるいは自然（Neumann）として処理。

固有値探索と固有ベクトルの計算

有限要素法の離散化により得られる一般化固有値問題は

$$ K\mathbf{c}=\lambda M\mathbf{c},\qquad \lambda=\omega^2 $$

である．固有周波数は$\omega=\sqrt{\lambda}$，周波数$f=\omega/(2\pi)$で与えられる．数値的には以下の方法が用いられる：

1. シフト＆逆反復（shift-and-invert）：特定の周波数近傍の固有値を効率よく求める．$(K-\sigma M)^{-1}M$に対する反復で，シフト$\sigma$を適切に選ぶ．
2. レイリー商反復（Rayleigh quotient iteration）：収束が高速だが初期値に依存する．
3. サブスペース反復法・Lanczos法：複数の最初の固有値固有ベクトルを同時に求めるのに有効．

反復法で得られた固有ベクトル$\mathbf{c}$（離散化された$E_z$の係数）は正規化して物理量に結び付ける．例えば電場エネルギーで正規化する場合：

$$ W=\frac{1}{4}\int_\Omega \epsilon |\mathbf{E}|^2\,dV = \frac{1}{4}\mathbf{c}^T M_\epsilon \mathbf{c}, $$

ここで$M_\epsilon$は誘電率を重みとする質量行列であり，必要に応じて$\mathbf{c}$を$W=1$にスケーリングする．またモードの直交性は質量行列に関して成り立つ：

$$ \mathbf{c}_i^T M \mathbf{c}_j =0\quad (i\neq j). $$

数値実装上の注意：

• 剛性行列$K$はしばしば対称正定値（境界条件による零固有値に注意）．
• 数値安定化（行列のスパース性の保持，適切な前処理器の使用）が計算効率に重要．
• 得られた$E_z$から他の場成分（$E_r,H_\phi$等）を復元して，散逸・損失の評価やQ値計算を行う．

参考：Superfish の入門・物理解説は[references/Docs/SFINTRO.DOC](references/Docs/SFINTRO.DOC#L1)および[references/Docs/SFPHYS1.DOC](references/Docs/SFPHYS1.DOC#L1)等に詳しい。

• びっくりマーク'!'があると，それ以降はコメント文になる．

Poissonも同様であるが，Superfish の計算結果の図で余白が大きすぎ見苦しいことがある．この場合，以下のようにして調整できる．もっと，良い方法があるかもしれませんが，私の方法です．

1. *.T35ファイルをダブルクリックして，計算結果の図を表示させる．メニューの HardCopy の Driver を SVG(Scalable Vector Graphics)に設定する．再度，メニューの HardCopy の Start を選択する．すると，「ファイル名.SVG」という図のファイルができあがる．
2. このファイルを inkscapeのようなSVGフォーマットの図を編集できるソフトウェアーで編集する．図がグループ化されている場合，それを解除してから編集をする．
3. 不要な余白を取り除いたら，eps形式で保存する．私はTeXを使うのでepsで保存しますが，自分に敵したフォーマットで保存すれば良いでしょう．

参考資料

更新履歴