電磁気学共振空胴 TM0nモードマクスウェルの方程式から弱形式

図1と図2は，加速器でしばしば使われる軸対称の共振空胴です．通常はTM01モードが使われますが，TM02モードなどの高次モードを調べたい場合もあります．ここでは，軸対称共振空洞の共振周波数(固有値)と電磁場分布(固有関数)を有限要素法で計算するために必要な弱形式を説明します．

Maxwell方程式から解くべき微分方程式

内部が真空の軸対称構造の加速空洞内の電磁場を表す偏微分方程式は，

\begin{align} & \nabla^2\boldsymbol{E}+\left(\frac{\omega}{c}\right)^2\boldsymbol{E}=0 & & \nabla^2\boldsymbol{H}+\left(\frac{\omega}{c}\right)^2\boldsymbol{H}=0 & \end{align}

です．TM0nモードの場合，これらのベクトルの方程式は，成分ごとの偏微分方程式:

\begin{align} \def\pd#1#2{\frac{\partial #1}{\partial #2}} \def\pdd#1#2{\frac{\partial^2 #1}{\partial #2^2}} & \pdd{E_r}{r}+\frac{1}{r}\pd{E_r}{r}+\pdd{E_r}{z}-\frac{E_r}{r^2}+\left(\frac{\omega}{c}\right)^2E_r=0 \\[6pt] & \pdd{E_z}{r}+\frac{1}{r}\pd{E_z}{r}+\pdd{E_z}{z}+\left(\frac{\omega}{c}\right)^2E_z=0 \\[6pt] & \pdd{H_\phi}{r}+\frac{1}{r}\pd{H_\phi}{r}+\pdd{H_\phi}{z}-\frac{H_\phi}{r^2}+\left(\frac{\omega}{c}\right)^2H_\phi=0& \label{eq:solve_Hphi} \end{align}

になります．これらの三つの方程式は境界条件や式の形が異なっています．しかし，物理的に正しい境界条件を設定すれば，これら三つの偏微分方程式 (PDE) の固有値 $\omega/c$ は同一になります．これら三つの偏微分方程式のいずれを解いても，同じ共振モードを得ることができます．ここでは，境界条件を扱いやすい式\eqref{eq:solve_Hphi}を使います．

未知関数 $u=rH_\phi$

式\eqref{eq:solve_Hphi}をそのまま有限要素法で解くと，軸 $r=0$ の近くに現れる $1/r$ や $1/r^2$ の項を丁寧に扱う必要があります．そこで，磁場そのもの $H_\phi$ ではなく，

\begin{align} u(r,z)=rH_\phi(r,z) \label{eq:def_u} \end{align}

を未知関数にします．軸対称のTM0nモードでは，軸 $r=0$ の上で方位角方向は定義できないため，物理的に正則な場では $H_\phi=0$ となります．さらに軸近傍では $H_\phi\propto r$ なので，$u=rH_\phi\propto r^2$ となり，$u$ もその微分も軸上で自然に小さくなります．ただし，$u$ を使う主な理由は，単に軸上でゼロになるからではありません．$H_\phi$ の方程式に現れる $1/r$ や $1/r^2$ の特異的な項を，弱形式で扱いやすい形に整理できることが重要です．このため，$u=rH_\phi$ は軸を含む有限要素計算に適した未知関数になります．

式\eqref{eq:def_u}を式\eqref{eq:solve_Hphi}に代入すると，真空中では $u$ が満たす方程式は，

\begin{align} \pdd{u}{r}+\pdd{u}{z}-\frac{1}{r}\pd{u}{r} +\left(\frac{\omega}{c}\right)^2u=0 \label{eq:pde_u_vacuum} \end{align}

です．ここでは空洞内部を真空として扱います．数学で固有値を表すときによく使われる記号に合わせて，固有値を $\lambda=\omega$ と書き，真空中の光速を $c$ とします．式\eqref{eq:pde_u_vacuum}は，有限要素法で使いやすいように次の形にも書けます．

\begin{align} \pd{}{r}\left(r\pd{u}{r}\right) +\pd{}{z}\left(r\pd{u}{z}\right) -2\pd{u}{r} +r\left(\frac{\lambda}{c}\right)^2u=0 \label{eq:pde_u} \end{align}

式\eqref{eq:pde_u}を $r$ で割ると式\eqref{eq:pde_u_vacuum}に戻ります．求めたい共振周波数は，

\begin{align} f=\frac{\lambda}{2\pi} \end{align}

から得られます．

弱形式

有限要素法では，微分方程式を各点で直接満たす代わりに，任意の試験関数 $v$ を掛けて領域全体で積分します．ここでは，$u_r=\partial u/\partial r$，$u_z=\partial u/\partial z$，$v=\mathrm{test}(u)$ と書きます．式\eqref{eq:pde_u}に $v$ を掛けると，

\begin{align} \int_\Omega \left\{ v\pd{}{r}\left(ru_r\right) +v\pd{}{z}\left(ru_z\right) -2u_rv +r\left(\frac{\lambda}{c}\right)^2uv \right\}\diff\Omega=0 \end{align}

となります．第1項と第2項には $u$ の2階微分が含まれています．有限要素法では1階微分までに下げたいので，この二つの項だけを部分積分します．例えば $r$ 方向の項は，

\begin{align} \int_\Omega v\pd{}{r}\left(ru_r\right)\diff\Omega =-\int_\Omega ru_r\pd{v}{r}\diff\Omega +\int_\Gamma vru_r n_r \diff\Gamma \end{align}

です[証明]．$z$ 方向も同様に部分積分すると，二つの境界積分はまとめて，

\begin{align} \int_\Gamma vr\left(u_r n_r+u_z n_z\right)\diff\Gamma =\int_\Gamma vr\pd{u}{n}\diff\Gamma \end{align}

と書けます．ここで $\partial u/\partial n$ は境界の外向き法線方向微分です．完全導体壁では接線方向の電場がゼロです．電場の復元式から，接線方向の電場は $\partial u/\partial n$ に比例するため，完全導体壁上では $\partial u/\partial n=0$ となります[証明]．また，対称軸 $r=0$ では境界積分の係数に $r$ が掛かっているため，この項はゼロになります．したがって境界積分は残らず，領域内に残る弱形式は，

\begin{align} \int_\Omega \left\{ -r\left(u_r\pd{v}{r}+u_z\pd{v}{z}\right) -2u_rv +r\left(\frac{\lambda}{c}\right)^2uv \right\}\diff\Omega=0 \label{eq:weak_form} \end{align}

となります．式\eqref{eq:weak_form}の第1項は $u$ の勾配どうしの積，第2項は $u=rH_\phi$ へ変数変換したことで現れる補正項，第3項は固有値 $\lambda$ を含む項です．この三つの寄与が全領域で釣り合うような $u$ と $\lambda$ が，空洞のTM0nモードになります．

境界条件は次のように考えます．完全導体壁では接線方向の電場がゼロになるため，$u$ に対しては自然境界条件として扱えます．一方，対称軸 $r=0$ では $u=rH_\phi$ なので $u=0$ を課します．この置き換えにより，軸上の特異性を避けながらTM0nモードを計算できます．

有限要素法の計算手順

有限要素法では，空洞の断面を小さな三角形や四角形の要素に分割し，未知関数 $u$ を節点値の和で近似します．基底関数を $N_i$，未知の節点値を $a_i$ とすると，

\begin{align} u(r,z)\simeq\sum_i a_iN_i(r,z) \end{align}

と書けます．この式を弱形式\eqref{eq:weak_form}に代入すると，連続的な固有値問題は，次の行列固有値問題になります．

\begin{align} \boldsymbol{A}\boldsymbol{a}=\left(\frac{\lambda}{c}\right)^2\boldsymbol{B}\boldsymbol{a} \label{eq:matrix_eigen} \end{align}

ここで，試験関数を $N_i$，未知関数側を $N_j$ とすると，行列の各成分は，

\begin{align} A_{ij}&=\int_\Omega \left\{ r\left(\pd{N_j}{r}\pd{N_i}{r}+\pd{N_j}{z}\pd{N_i}{z}\right) +2\pd{N_j}{r}N_i \right\}\diff\Omega,\\ B_{ij}&=\int_\Omega rN_jN_i\,\diff\Omega \end{align}

です．数値計算では，それぞれの要素ごとに小さな行列を作り，それを全体行列へ足し込んでから，式\eqref{eq:matrix_eigen}を解きます．得られた固有値 $(\lambda/c)^2$ の小さい方から順に，TM01，TM02，... のような軸対称モードに対応します．

計算の流れをまとめると，次のようになります．

1. 共振空胴の $r$-$z$ 断面を作る．
2. 断面を有限要素メッシュに分割する．
3. 各要素で $A_{ij}$ と $B_{ij}$ の積分を計算する．
4. 全体の $\boldsymbol{A}$ と $\boldsymbol{B}$ を組み立てる．
5. 軸 $r=0$ の節点に $u=0$ を課す．
6. $\boldsymbol{A}\boldsymbol{a}=(\lambda/c)^2\boldsymbol{B}\boldsymbol{a}$ を解く．
7. $f=\lambda/(2\pi)$ により共振周波数を求める．

電磁場の復元

固有ベクトル $\boldsymbol{a}$ から $u$ が得られたら，磁場は

\begin{align} H_\phi=\frac{u}{r} \end{align}

で求まります．また，Maxwell方程式から真空中の電場成分も次のように計算できます．

\begin{align} E_r&=-\frac{1}{i\omega\varepsilon_0 r}\pd{u}{z},\\ E_z&= \frac{1}{i\omega\varepsilon_0 r}\pd{u}{r} \end{align}

したがって，有限要素法で直接求めるのは $u=rH_\phi$ ですが，そこから $H_\phi$，$E_r$，$E_z$ を復元できます．このため，$u$ はTM0nモードの電磁場分布を安定に求めるために便利な未知関数です．

これまでは，空洞内部が真空の場合を扱いました．ここでは，図3のように内部に誘電体が含まれる場合を扱います．この場合も，未知関数 $u=rH_\phi$ を使うと，軸上の特異性を避けながら計算できます．磁性体を含む場合も同様に扱えますが，係数に比透磁率 $\mu_r$ が現れるため，材料定数を含むMaxwell方程式から改めて式を導く必要があります．

Maxwell方程式から解くべき微分方程式

まず，比誘電率 $\varepsilon_r$ の誘電体を考えます．ここでは非磁性体を仮定し，$\mu=\mu_0$ とします．誘電率は $\varepsilon=\varepsilon_0\varepsilon_r$ です．時間依存性を $\exp(i\omega t)$ とすると，Maxwell方程式は，

\begin{align} \nabla\times\boldsymbol{H}&=i\omega\varepsilon_0\varepsilon_r\boldsymbol{E},\\[6pt] \nabla\times\boldsymbol{E}&=-i\omega\mu_0\boldsymbol{H} \end{align}

です．軸対称TM0nモードでは，磁場は方位角方向成分 $H_\phi$ だけを持つので，

\begin{align} \boldsymbol{H}=H_\phi\boldsymbol{e}_\phi,\qquad u=rH_\phi \end{align}

と書けます．このとき，$\nabla\times\boldsymbol{H}=i\omega\varepsilon_0\varepsilon_r\boldsymbol{E}$ の $r$ 成分と $z$ 成分から，

\begin{align} E_r&=-\frac{1}{i\omega\varepsilon_0\varepsilon_r r}\pd{u}{z}\label{eq:media:Er}\\ E_z&= \frac{1}{i\omega\varepsilon_0\varepsilon_r r}\pd{u}{r}\label{eq:media:Ez} \end{align}

を得ます．次に，$\nabla\times\boldsymbol{E}=-i\omega\mu_0\boldsymbol{H}$ の $\phi$ 成分を使います．円筒座標では，

\begin{align} \pd{E_r}{z}-\pd{E_z}{r}=-i\omega\mu_0H_\phi \end{align}

です．ここに式\eqref{eq:media:Er}の $E_r$ と式\eqref{eq:media:Ez}の $E_z$ を代入し，$\mu_0\varepsilon_0=1/c^2$ を使って整理すると，

\begin{align} \pd{}{z}\left(\frac{1}{\varepsilon_r}\pd{H_\phi}{z}\right) +\pd{}{r}\left[\frac{1}{\varepsilon_r r}\pd{}{r}\left(rH_\phi\right)\right] +\left(\frac{\omega}{c}\right)^2H_\phi=0 \label{eq:media:pde_Hphi} \end{align}

が得られます．これが誘電体を含む場合に解くべき偏微分方程式です．

変数変換と弱形式

この方程式を未知関数 $u=rH_\phi$ と固有値 $\lambda=\omega$ を使って書き換えます．式\eqref{eq:media:pde_Hphi}は，

\begin{align} \pd{}{z}\left(\frac{1}{\varepsilon_r}\pd{u}{z}\right) +\pd{}{r}\left(\frac{1}{\varepsilon_r}\pd{u}{r}\right) -\frac{1}{\varepsilon_r r}\pd{u}{r} +\left(\frac{\lambda}{c}\right)^2u=0 \end{align}

となります．最後に，有限要素法の弱形式で扱いやすいように全体に $r$ を掛け，$r$ 方向の項を積の微分として変形すると，

\begin{align} \pd{}{r}\left(\frac{r}{\varepsilon_r}\pd{u}{r}\right) +\pd{}{z}\left(\frac{r}{\varepsilon_r}\pd{u}{z}\right) -\frac{2}{\varepsilon_r}\pd{u}{r} +r\left(\frac{\lambda}{c}\right)^2u=0 \end{align}

が得られます．これを弱形式にすると，

\begin{align} \int_\Omega \left\{ \frac{r}{\varepsilon_r}\left(u_r\pd{v}{r}+u_z\pd{v}{z}\right) +\frac{2}{\varepsilon_r}u_rv -r\left(\frac{\lambda}{c}\right)^2uv \right\}\diff\Omega=0 \end{align}

となります．この式の第1項は $u$ の勾配と試験関数の勾配の積，第2項は変数変換による補正項，第3項は固有値を含む項です．この三つの寄与が全領域で釣り合うような $u$ と $\lambda$ が，誘電体を含む空洞のTM0nモードになります．

参考資料

更新履歴