軸対称問題の FDTD 法 (長方形メッシュ)

山本昌志^{∗}

2005 年 11 月 29 日

1 はじめに

時間領域の電磁場の解析には，時間領域差分法(Finite Difference Time Domain:FDTD) がよく使われる．これについては，いろいろな教科書で説明されている．ほとんどの教科書はカーテシアン座標系 (x, y, z) が使われているが，我々の解析対象である加速器では円柱座標系 (r, θ, z) の方が便利である．そこで，円柱座標系の FDTD 法で電磁場を解析するときの式を導くことにする．

ここでは，我々の主な解析対象は TM_{0n} モードのみを取り扱うものとする．他のモードでも同じようにして式が導くことが出きるので，必要に応じて計算すれば良い．

当然，ここに載せている式はそんなに難しいものでは無い．誰かが，とっくの昔に論文や教科書に公表しているかもしれない．だが，私はそれを持っていないので，引用文献を載せることはできない．引用文献とすべきものがあれば教えて欲しい．

2 Maxwell の方程式

2.1 解析対象

我々は，加速空洞内の電磁場の伝搬を計算している．内部は真空で，当面は電荷が無いものとして計算を進める．電荷がある場合については，次のステップで考えるつもりである．

さらに，取り扱うモードは TM_{0} モードとする．これは，加速器空洞で使われるもっとも一般的なモードである．ほとんどの加速空洞は軸対称構造なので，このモードを円柱座標系で解析することになる．一般的，電場は 3 個の成分 (E_{r}, E_{θ}, E_{z})，磁場も 3 個の成分 (H_{r}, H_{θ}, H_{z}) をもち，それぞれは位置 (r, θ, z) の関数である．しかし，ここでは TM_{0} モードが解析対象なので，

E_{θ} = 0  H_{r} = 0  H_{z} = 0  (1)

となる．さらに，残りの E_{r}, E_{z}, H_{θ} は，位置 (r, z) の関数となり， θ に依存しない． θ 方向に一様なのである．

∗国立秋田工業高等専門学校 電気情報工学科

2.2 解くべき方程式

真空中で，電荷も電流もない場合の Maxwell の方程式は，

∇ · E = 0  (2)
∇ · H = 0  (3)
∇ × E = −μ_{0} \frac{∂H}{∂t}  (4)
∇ × H = ε_{0} \frac{∂E}{∂t}  (5)

である．ここで，解くべきものは，式 (4) と (5) である．なぜ，最初の 2 つを無視してよいかはここでは述べない．時間があるときに述べることにする．

FDTD 法で計算すべき式は積分形で表す方が良い．これら，微分形の式をストークスの定理を用いて，積分形に直すと

∮ E · dℓ = −μ_{0} \frac{d}{dt} ∫_{S} H · ndS  (6)
∮ H · dℓ = ε_{0} \frac{d}{dt} ∫_{S} E · ndS  (7)

となる．真空中の誘電率と透磁率の間には，

c = \frac{1}{\sqrt{ε_{0}μ_{0}}}  (8)

の関係があることを忘れてはならない．ここで， c は光速を表す．また，これら電場の強さ E と束密度 D，磁場の強さ H と磁束密度 B には次の関係がある．

D = ε_{0}E  (9)
B = μH  (10)

3 メッシュと時間の刻み

軸対称構造の解析なので，本来は 2 次元問題である．しかし，解くべきマクスウェルの方程式 (6) と (7) は 3 次元で表現されている．そこで，ここでは 3 次元の式から出発して， 2 次元に直すことにする．

解析の対象は図 1 に示すような軸対称構造物である．対称性を考慮すると， rz 平面の上半分のみを計算すればよい．解析の対象としている TM_{0n} では，上半分と下半分は同じ電磁場となるからである．このような単極モードではなくとも多極モードでも，上半分が計算できれば下半分の計算ができるので，やはり半分のみを解析の対象とすればよい．

rz 平面での形状が長方形で形作られているが，これが最終的に計算するメッシュとなる．この長方形の面の中心に磁場 H_{θ} が，辺の中心に電場 E_{r} と E_{z} がある．このように長方形により離散化された電磁場を計算することになる．多くの長方形により計算する領域が構成されることになるが，長方形の大きさは一様である必要はない．だが一般的には計算を容易にするために，実際にはプログラムを容易にするために，一様な大きさの正方形を用いる．

電磁場のみならず， FDTD 法では時刻も離散化しなくてはならない．本間らの教科書 [1] に示しているように，電場と磁場を交互に計算するように，時刻を計算する (図 3)．長方形メッシュにより，電場と磁場を計算する位置が交互になっているように，時間領域でも交互に計算するのである．計算する時間間隔もまた一定にする必要はないが，一定間隔にする方が計算が容易である．

電場も磁場も位置 (z, r) と時刻 (t) の関数で，通常の表記だと式が長くなる．そこで，簡潔に式を記述するために，参考文献 [1] と同様，以下の表現を用いることにする．

E_{z}(z_{i+1}, r_{j+3/2}, t_{n+1/2}) = E_{z i+1,j+3/2}^{n+1/2}  (11)
E_{r}(z_{i+3/2}, r_{j+1}, t_{n+1/2}) = E_{r i+3/2,j+1}^{n+1/2}  (12)
H_{θ}(z_{i+1}, r_{j+1}, t_{n+1}) = H_{θ i+1,j+1}^{n+1}  (13)

図 1: FDTD のメッシュ． rz 平面では，長方形メッシュとなる．

図 2: FDTD 法で電磁場計算のメッシュ．

図 3: FDTD 法で電磁場を計算するときの時刻の刻み方．

4 FDTD の式

4.1 磁場の θ 方向

FDTD 法で計算すべき式は，式 (6) と式 (7) を評価することにより得られる．これは，図 1 の構造を形作っている多数の立体，これを体積素と言うことにして，それらを個別に積分することにより求められる．

この体積素を構成する平面，これを面積素ということにして，このうちの一つ rz 平面で式 (6) を評価してみよう．体積素が小さいとして，この式の左辺は電場の線績分を表し，図 4 に示すように経路では，

∮ E · dℓ ≃ [E_{r i+1/2,j}^{n+1/2} − E_{r i−1/2,j}^{n+1/2}] (r_{j+1/2} − r_{j−1/2}) + [E_{z i,j−1/2}^{n+1/2} − E_{z i,j+1/2}^{n+1/2}] (z_{i+1/2} − z_{i−1/2}) (14)

と近似ができる．右辺の方の近似は，

−μ_{0} \frac{d}{dt} ∫_{S} H · ndS ≃ −\frac{μ_{0}}{t_{n+1} − t_{n}} [H_{θ i,j}^{n+1} − H_{θ i,j}^{n}] (z_{i+1/2} − z_{i−1/2})(r_{j+1/2} − r_{j−1/2})  (15)

となる．この両辺は等しいので，

[E_{r i+1/2,j}^{n+1/2} − E_{r i−1/2,j}^{n+1/2}] (r_{j+1/2} − r_{j−1/2}) + [E_{z i,j−1/2}^{n+1/2} − E_{z i,j+1/2}^{n+1/2}] (z_{i+1/2} − z_{i−1/2})

≃ −\frac{μ_{0}}{t_{n+1} − t_{n}} [H_{θ i,j}^{n+1} − H_{θ i,j}^{n}] (z_{i+1/2} − z_{i−1/2})(r_{j+1/2} − r_{j−1/2}) (16)

となる．これから，

H_{θ i,j}^{n+1} ≃ H_{θ i,j}^{n} + \frac{t_{n+1} − t_{n}}{μ_{0}} [\frac{E_{r i−1/2,j}^{n+1/2} − E_{r i+1/2,j}^{n+1/2}}{z_{i+1/2} − z_{i−1/2}} + \frac{E_{z i,j+1/2}^{n+1/2} − E_{z i,j−1/2}^{n+1/2}}{r_{j+1/2} − r_{j−1/2}} ]  (17)

が得られる．この式の右辺は，時刻 t_{n+1/2} の電場 (E_{z}, E_{r}) と t_{n−1} の磁場 (H_{θ}) である．左辺は，時刻 t_{n+1} の磁場 H_{θ} である． t_{n+1} の磁場は，それ以前の情報があれば計算できると言っているのである．

図 4: zr 平面の積分路

4.2 電場の r 方向

次に，先ほどと同様のことを図 5 の領域について行う．ここで TM モードについて計算しているので， H_{z} が無いことに注意して積分を行う．式 (7) の左辺の積分は，

∮ H · dℓ ≃ [H_{θ i−1,j}^{n} − H_{θ i,j}^{n}] r_{j}Δθ  (18)

となる．式 (7) の右辺の積分は，

ε_{0} \frac{d}{dt} ∫_{S} E · ndS ≃ \frac{ε_{0}}{t_{n+1/2} − t_{n−1/2}} [E_{r i−1/2,j}^{n+1/2} − E_{r i−1/2,j}^{n−1/2}] (z_{i} − z_{i−1})r_{j}Δθ  (19)

となる．両辺は等しいので，

H_{θ i−1,j}^{n} − H_{θ i,j}^{n} ≃ \frac{ε_{0}}{t_{n+1/2} − t_{n−1/2}} [E_{r i−1/2,j}^{n+1/2} − E_{r i−1/2,j}^{n−1/2}] (z_{i} − z_{i−1})  (20)

である．これから，電場の r 方向成分を求める式

E_{r i−1/2,j}^{n+1/2} ≃ E_{r i−1/2,j}^{n−1/2} + \frac{t_{n+1/2} − t_{n−1/2}}{ε_{0}(z_{i} − z_{i−1})} [H_{θ i−1,j}^{n} − H_{θ i,j}^{n}]  (21)

が得られる．先ほど同様，これも t_{n+1/2} の電場 E_{r} は，それ以前の電磁場から求められると言っている．

図 5: θz 平面の積分路

4.3 電場の z 方向

最後に，いままでと同様のことを図 6 の領域について行う．ここで TM モードについて計算しているので， H_{r} が無いことに注意して積分を行う．式 (7) の左辺の積分は，

∮ H · dℓ ≃ [H_{θ i,j+1}^{n}r_{j+1} − H_{θ i,j}^{n}r_{j}] Δθ  (22)

となる．式 (7) の右辺の積分は，

ε_{0} \frac{d}{dt} ∫_{S} E · ndS ≃ \frac{ε_{0}}{t_{n+1/2} − t_{n−1/2}} [E_{z i,j+1/2}^{n+1/2} − E_{r i,j+1/2}^{n−1/2}] \frac{1}{2}(r_{j+1}^{2} − r_{j}^{2})Δθ  (23)

となる．両辺は等しいので，

H_{θ i,j+1}^{n}r_{j+1} − H_{θ i,j}^{n}r_{j} ≃ \frac{ε_{0}}{2(t_{n+1/2} − t_{n−1/2})} [E_{z i,j+1/2}^{n+1/2} − E_{z i,j+1/2}^{n−1/2}] (r_{j+1}^{2} − r_{j}^{2})  (24)

となる．これから，電場の z 方向成分を求める式

E_{z i,j+1/2}^{n+1/2} ≃ E_{z i,j+1/2}^{n−1/2} + \frac{2(t_{n+1/2} − t_{n−1/2})}{ε_{0}(r_{j+1}^{2} − r_{j}^{2})} [H_{θ i,j+1}^{n}r_{j+1} − H_{θ i,j}^{n}r_{j}]  (25)
≃ E_{z i,j+1/2}^{n−1/2} + \frac{2(t_{n+1/2} − t_{n−1/2})}{ε_{0}(r_{j+1} + r_{j})(r_{j+1} − r_{j})} [H_{θ i,j+1}^{n}r_{j+1} − H_{θ i,j}^{n}r_{j}]  (26)

が得られる．いままで同様，これも t_{n+1/2} の電場 E_{z} は，それ以前の電磁場から求められると言っている．

図 6: rθ 平面の積分路

5 まとめ

軸対称構造の TM_{0n} モードを円柱座標系で表した FDTD 法の式は，以下のようになる．

H_{θ i,j}^{n+1} ≃ H_{θ i,j}^{n} + \frac{t_{n+1} − t_{n}}{μ_{0}} [\frac{E_{r i−1/2,j}^{n+1/2} − E_{r i+1/2,j}^{n+1/2}}{z_{i+1/2} − z_{i−1/2}} + \frac{E_{z i,j+1/2}^{n+1/2} − E_{z i,j−1/2}^{n+1/2}}{r_{j+1/2} − r_{j−1/2}} ]  (27)
E_{r i−1/2,j}^{n+1/2} ≃ E_{r i−1/2,j}^{n−1/2} + \frac{t_{n+1/2} − t_{n−1/2}}{ε_{0}(z_{i} − z_{i−1})} [H_{θ i−1,j}^{n} − H_{θ i,j}^{n}]  (28)
E_{z i,j+1/2}^{n+1/2} ≃ E_{z i,j+1/2}^{n−1/2} + \frac{2(t_{n+1/2} − t_{n−1/2})}{ε_{0}(r_{j+1} + r_{j})(r_{j+1} − r_{j})} [H_{θ i,j+1}^{n}r_{j+1} − H_{θ i,j}^{n}r_{j}]  (29)

参考文献

[1] 本間利久, 五十嵐一, 川口秀樹. 数値電磁力学. 計算電気・電子工学シリーズ 3. 森北出版株式会社, 2002.