FDTD法の基礎

0. はじめに

0.1 まえがき

今までは偏微分方程式の数値解析を研究していたのだが、突然FDTD法を利用する必要が出てきた。
これまではどちらかというと静電界で現れる方程式の研究だったため、FDTD法のような(動)電磁界解析は触れていなかった。

これまでの研究のおかげで、数値解析の部分はすんなりと理解できるのだが、電磁気学の知識が若干乏しいのもあって、結構詰まったりもする。
そこで、今回は自分のための知識整理、という意味で、ここにFDTD法の基礎を書き留めようと思う。

0.2 アウトライン

この記事では、1次元FDTD法の実装を目標とする。
2次元や3次元のほうがグラフで可視化したときにグラフィカルで面白いかもしれないが、
1次元のほうが、可視化したときにどんな電界や磁界の動きが見やすくてイメージしやすいため、1次元を選んだ。

といっても、実際に利用するときは3次元がほとんどだと思うため、3次元の定式化について簡単に説明したあと、プログラムは1次元FDTDで実装する、という流れにしようと思う。
なお、プログラムはpythonで作る。

また、吸収境界条件やフーリエ変換など、FDTD法を学習するにあたって欠かせないテクニックも簡単に解説していく。

電磁界解析のイメージが沸かない人は、3.4節の解析結果の図を軽く眺めて見てから読んだほうが良いかもしれない。

実装したプログラムは以下に置いている。
atily17/FDTD1d(github)

1. FDTD法のアルゴリズム

1.1 FDTD法とは

FDTD法は、Finite Difference Time Domain Methodの略で、日本語化すると、有限差分時間領域法である。
「有限差分」は簡単に言うと差分法であることを意味する。微分を差に置き換えるという意味だ。
「時間領域」は時間空間について解くことを意味する。対義語は「周波数領域」だ。

「有限差分」を「有限要素」に変えたFETD(Finite Element Time Domain)法や、
「時間領域」を「周波数領域」に変えたFDFD(Finite Difference Frequency Domain)法などもある。ただ、FDTD法と比べるとマイナー。

FDTD法は有限差分法を基礎においているため、定式化が結構簡単である。
ただ、収束条件などいろんな注意事項があるため、必要な知識はできる限り説明していきたい。

1.2 マクスウェル方程式

まず、簡単に電磁気学の基礎方程式を解説しよう。

電磁気学では、マクスウェル方程式がすべてを表しているため、どんな現象もマクスウェル方程式で説明できる。
マクスウェル方程式は以下の四式だ。

\begin{alignat}{2}
&\displaystyle\nabla \cdot \boldsymbol{D} & =& \rho \\\\
&\displaystyle\nabla \cdot \boldsymbol{B} & =& 0 \\\\
&\displaystyle\nabla \times \boldsymbol{E} + \frac{\partial \boldsymbol{B}}{\partial t} & = & 0 \\\\
&\displaystyle\nabla \times \boldsymbol{H} - \frac{\partial \boldsymbol{D}}{\partial t} & = & \boldsymbol{j}
\end{alignat}

ここで、$\boldsymbol{D}$は電束密度、$\boldsymbol{E}$は電界、$\boldsymbol{B}$は磁束密度、$\boldsymbol{H}$は磁界を表し、$\rho$は電荷密度、$\boldsymbol{j}$は電流密度である。¹

ただ、実はこれだけでは足りなく、$\boldsymbol{D}$と$\boldsymbol{E}$の関係と、$\boldsymbol{B}$と$\boldsymbol{H}$の関係の二式も必要になる。
材料それぞれの物性値である、誘電率$\varepsilon$と、透磁率$\mu$の比例関係が成り立つと仮定される。これは、実験で与えられた近似値であることを覚えておこう。

\begin{align}
\boldsymbol{D} & \approx \varepsilon \boldsymbol{E} \\\\
\boldsymbol{B} & \approx \mu \boldsymbol{H}
\end{align}

あえて、ここではイコールではなく、$\approx$で書くことにした。
材料によってはテンソルになったり、単純な比例関係じゃなかったりもするが、今回は、この比例関係が成り立つことを前提に進めていこう。

ただし、真空の場合は真空の誘電率$\varepsilon_0$と真空の透磁率$\mu_0$から

\begin{align}
\boldsymbol{D} & = \varepsilon_0 \boldsymbol{E} \\\\
\boldsymbol{B} & = \mu_0 \boldsymbol{H}
\end{align}

が厳密に成り立つことも覚えておこう。

さらに、FDTD法では、電界と電流を関係付ける式である、
$$
\boldsymbol{j} = \sigma \boldsymbol{E}
$$
という式も利用する。この$\sigma$は導電率、もしくは電気伝導率と呼ばれる。
これは、オームの法則を微分形で表現したものである。²

これで、FDTD法の定式化で必要な式は揃った。

1.3 FDTD法の定式化

1.3.1 マクスウェル方程式の変形

FDTD法で使うマクスウェル方程式は以下の二つだ。

\begin{align}
\displaystyle\nabla \times \boldsymbol{E} + \frac{\partial \boldsymbol{B}}{\partial t} = 0 \\\\
\displaystyle\nabla \times \boldsymbol{H} - \frac{\partial \boldsymbol{D}}{\partial t} = \boldsymbol{j}
\end{align}

ここに、誘電率、透磁率、導電率の関係式である、

\begin{align}
\boldsymbol{D} &= \varepsilon \boldsymbol{E} \\\\
\boldsymbol{B} &= \mu \boldsymbol{H} \\\\
\boldsymbol{j} &= \sigma \boldsymbol{E}
\end{align}

を代入しよう。

式(1) マクスウェル方程式の変形1

$$
\displaystyle\nabla \times \boldsymbol{E} = - \mu \frac{\partial \boldsymbol{H}}{\partial t} \tag{1}
$$

式(2) マクスウェル方程式の変形2

$$
\displaystyle\nabla \times \boldsymbol{H} = \varepsilon \frac{\partial \boldsymbol{E}}{\partial t} + \sigma \boldsymbol{E} \tag{2}
$$

この、$\mu$と$\epsilon$と$\sigma$は材料によって決まる既知の値である。
よって、これでたくさんあった変数がたった二つになった。
あとはこれを差分式にするだけだ！　簡単！

1.3.2 Yee格子

さて差分近似！　と行きたいところだが、そのまえにYee格子について説明する必要がある。

FDTD法は差分法のため、計算領域を格子分割する必要がある。
Yee格子は簡単に言うと、$\boldsymbol{E}$と$\boldsymbol{H}$で格子のとり方をずらすことで、計算しやすく、そして、高精度にできるというものだ。
まあ、そんなに難しくないため、簡単に説明しよう。

$E$の格子は、下のように、直方体の格子をとする。

ここで、直方体の大きさは、$\Delta x \times \Delta y \times \Delta z$とする。

直方体の頂点の座標は$(i\Delta x, j\Delta y, k \Delta z)$である。
毎度毎度、$\Delta x$など書くのは面倒なので、$[i,j,k]$と表現することにする。
例えば、上の図の、濃い緑の直方体の手前左下の頂点は、$[3,2,0]$となる。

そして、この表現を使って、上の図の$E_x$, $E_y$, $E_z$は、
$$
\begin{matrix}
E_x:& [&3, & \displaystyle3-\frac{1}{2}, &\displaystyle 1-\frac{1}{2} &] \\
E_y:& [&\displaystyle 3-\frac{1}{2}, &3, &\displaystyle 1-\frac{1}{2} &] \\
E_z:& [&\displaystyle 3-\frac{1}{2}, &\displaystyle 3-\frac{1}{2}, &1 &]
\end{matrix}
$$
の位置にあることになる。
つまり、それぞれの電界は、
$$
\begin{matrix}
E_x:& [&i, &\displaystyle j-\frac{1}{2}, &\displaystyle k-\frac{1}{2} &] \\
E_y:& [&\displaystyle i-\frac{1}{2}, &j, &\displaystyle k-\frac{1}{2} &] \\
E_z:& [&\displaystyle i-\frac{1}{2}, &\displaystyle j-\frac{1}{2}, &k &]
\end{matrix}
$$
の位置にあることになる。
この格子に関しては、特になんの工夫もない。

次に、磁界$\boldsymbol{H}$の格子について考えよう。
これも、特に違いはない。

直方体の大きさは、電界と同じ$\Delta x \times \Delta y \times \Delta z$に設定する。
そして、ここがYee格子の工夫点で、$\boldsymbol{H}$の格子の頂点が、$\boldsymbol{E}$の格子のそれぞれの直方体の中心に来るようにずらして置くのである。
以下に図を書いた³。

このようにずらした格子のため、磁界の位置は、
$$
\begin{matrix}
H_x:& [&\displaystyle i-\frac{1}{2}, &j, & k &] \\
H_y:& [& i, & \displaystyle j -\frac{1}{2}, &k &] \\
H_z:& [& i, & j, &\displaystyle k -\frac{1}{2}&]
\end{matrix}
$$
となる。

説明の都合上、EとHは以上の座標のとり方になったが、原点を$(\Delta x/2, \Delta y/2, \Delta z/2)$ずらせば、以下のようになる。

式(3): Eの位置

$$
\begin{matrix}
E_x:& [&\displaystyle i+\frac{1}{2}, &j, & k &] \\
E_y:& [& i, & \displaystyle j +\frac{1}{2}, &k &] \\
E_z:& [& i, & j, &\displaystyle k +\frac{1}{2}&]
\end{matrix} \tag{3}
$$

式(4): Hの位置

$$
\begin{matrix}
H_x:& [&i, &\displaystyle j+\frac{1}{2}, &\displaystyle k+\frac{1}{2} &] \\
H_y:& [&\displaystyle i+\frac{1}{2}, &j, &\displaystyle k+\frac{1}{2} &] \\
H_z:& [&\displaystyle i+\frac{1}{2}, &\displaystyle j+\frac{1}{2}, &k &]
\end{matrix}\tag{4}
$$

多くの解説書では、このように設定している場合が多いため、今後はこっちで説明することにする。

2次元の場合は、3次元で適切な$yz$平面を切り出すことで、

となり、
1次元の場合は、2次元格子のある直線を切り出すことで、

となる。

また、マクスウェル方程式には時間微分も現れるため、時間も離散化する必要がある。
時間刻み幅を$\Delta t$とすると、$t = n\Delta t$と表すことができる。
Yee格子では、空間と同様に時間も$\boldsymbol{E}$の計算と$\boldsymbol{H}$の計算を$\Delta t/2$ずらす。
つまり、$\boldsymbol{E}$の値は、$t = n\Delta$の時間に、$\boldsymbol{H}$の値は、$t=(n+1/2)\Delta t$の時間のものを求める($n$は自然数)。

時間$t=n\Delta t$における$E_{x}[i,j-1/2,k-1/2]$を$E_{x}^{n}[i,j-1/2,k-1/2]$と、時間$t=(n+1/2)\Delta t$における$H_{x}[i-1/2,j,k]$を$H_{x}^{n+1/2}[i-1/2, j, k]$と表現することにしよう。

まあ、こんなふうに、Yee格子は$\boldsymbol{E}$の格子から空間や時間を半分ずらしたものを$\boldsymbol{H}$とする、ってだけだ。

1.3.3 差分式

さて、ついに差分近似だ。

今更説明することもないかもしれないが、差分近似とは、
(偏)微分の定義
$$
\frac{d f}{d x} = \lim_{\Delta x\rightarrow 0} \frac{f(x+\Delta x) - f(x)}{\Delta x}
$$
において、$\Delta x$が十分に小さいとすると、
$$
\frac{df}{dx} \approx \frac{f(x+\Delta x) - f(x)}{\Delta x}
$$
のように近似できる、というものだ。
$x \rightarrow x-\Delta x/2$と平行移動すると、上の式は、
$$
\frac{df}{dx} \approx \frac{\displaystyle f(x+\frac{\Delta x}{2}) - f(x-\frac{\Delta x}{2})}{\Delta x}
$$
とできる。これは、中心差分と呼ばれ、このほうが一般に精度がよくなる。⁴

さて、もう一度、(変形した)マクスウェル方程式を書き出しておこう。

$$
\displaystyle\nabla \times \boldsymbol{E} = - \mu \frac{\partial \boldsymbol{H}}{\partial t} \tag{1}
$$
$$
\displaystyle\nabla \times \boldsymbol{H} = \varepsilon \frac{\partial \boldsymbol{E}}{\partial t} + \sigma \boldsymbol{E} \tag{2}
$$
まず、式(1)の$x$成分

$$
\frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z} = - \mu \frac{\partial H_x}{\partial t}
$$
のみについて考えてみよう。

(3)と(4)を思い出すと、
$$
\begin{matrix}
E_y:& [& i, &\displaystyle j+\frac{1}{2}, &k &] \\
E_z:& [& i, & j, &\displaystyle k +\frac{1}{2} &]
\end{matrix} \tag{3}
$$
$$
\begin{matrix}
H_x:& [& i, &\displaystyle j+\frac{1}{2}, &\displaystyle k+\frac{1}{2} &]
\end{matrix} \tag{4}
$$
だった。

整数と半整数(整数-1/2)がどの座標かに注意して、
$$
\begin{array}{l}
\text{時間} : n\Delta t,\\
\text{空間} : [i,j+\frac{1}{2}, k+\frac{1}{2}]
\end{array}
$$
について、全て中心差分をとると、
$$
\begin{array}{ll}
\text{左辺第一項}:&\displaystyle \frac{\partial E_z}{\partial y} \approx \frac{E_z^n [i, j+1, k+\displaystyle\frac{1}{2}] - E_z^n [i, j, k+\displaystyle\frac{1}{2}]}{\Delta y} \\
\text{左辺第二項}:&\displaystyle \frac{\partial E_y}{\partial z} \approx \frac{E_y^n [i, j+\displaystyle\frac{1}{2}, k+1] - E_y^n [i, j+\displaystyle\frac{1}{2}, k]}{\Delta z} \\
\text{右辺第一項}:&\displaystyle \frac{\partial H_x}{\partial t} \approx \frac{H_x^{n+1/2}[i,j+\displaystyle\frac{1}{2}, k+\frac{1}{2}] - H_x^{n-1/2}[i,j+\displaystyle\frac{1}{2}, k+\frac{1}{2}]}{\Delta t}
\end{array}
$$
と近似できるため、これを代入することで、

\begin{multline}
\frac{E_z^n [i, j+1, k+\displaystyle\frac{1}{2}] - E_z^n [i, j, k+\displaystyle\frac{1}{2}]}{\Delta y} - 
\frac{E_y^n [i, j+\displaystyle\frac{1}{2}, k+1] - E_y^n [i, j+\displaystyle\frac{1}{2}, k]}{\Delta z} \\\\
= - \mu \frac{H_x^{n+1/2}[i,j+\displaystyle\frac{1}{2}, k+\frac{1}{2}] - H_x^{n-1/2}[i,j+\displaystyle\frac{1}{2}, k+\frac{1}{2}]}{\Delta t}
\end{multline}

という式ができる。式が長く複雑に見えるが、全然難しいことはしていない。
$y$成分、$z$成分に関しても同様に導出できる。
ぜひ自分で計算してみよう。

次に、式(2)の$x$成分も差分化してみよう。
$$
\frac{\partial H_z}{\partial y}-\frac{\partial H_y}{\partial z} = \varepsilon\frac{\partial E_x}{\partial t} + \sigma E_x
$$

式(3)と式(4)から
$$
\begin{matrix}
E_x:& [& i+\displaystyle\frac{1}{2}, &j, &k &] \\
\end{matrix} \tag{3}
$$
$$
\begin{matrix}
H_y:& [&\displaystyle i+\frac{1}{2}, &j, &\displaystyle k+\frac{1}{2} &] \\
H_z:& [&\displaystyle i+\frac{1}{2}, &\displaystyle j+\frac{1}{2}, &k &]
\end{matrix} \tag{4}
$$
だった。

式(1)のときと半分だけズラした、
$$
\begin{array}{l}
\text{時間} : \displaystyle(n+\frac{1}{2})\Delta t,\\
\text{空間} : [i+\frac{1}{2},j, k]
\end{array}
$$
について中心差分を取ると、次のようになる。
$$
\begin{array}{ll}
\text{左辺第一項}:&\displaystyle \frac{\partial H_z}{\partial y} \approx \frac{H_z^{n+1/2} [\displaystyle i+\frac{1}{2}, j+\frac{1}{2}, k] - H_z^{n+1/2} [i + \frac{1}{2}, j - \frac{1}{2}, k]}{\Delta y} \\
\text{左辺第二項}:&\displaystyle \frac{\partial H_y}{\partial z} \approx \frac{H_y^{n+1/2} [i + \displaystyle\frac{1}{2}, j, k+\frac{1}{2}] - H_y^{n+1/2} [i+ \displaystyle\frac{1}{2}, j, k-\frac{1}{2}]}{\Delta z} \\
\text{右辺第一項}:&\displaystyle \frac{\partial E_x}{\partial t} \approx \frac{E_x^{n+1}[i+\displaystyle\frac{1}{2},j, k] - E_x^{n}[i+\displaystyle\frac{1}{2},j, k]}{\Delta t} \\
\text{右辺第二項}:&\displaystyle E_x = E_x^{n+1/2}[i+\frac{1}{2},j, k] \approx \frac{E_x^{n+1}[i+\displaystyle\frac{1}{2},j, k] + E_x^{n}[i+\displaystyle\frac{1}{2},j, k]}{2}
\end{array}
$$
特筆すべきは右辺第二項だ。
これは、もともと微分じゃないので、差分近似する必要はない。
しかし、Yee格子では電界の値は時間が$n\Delta t$のときしか使わないと決めているのに対して、$(n+\frac{1}{2})\Delta t$となってしまっているため、前後の時間($n\Delta t$と$(n+1)\Delta t$)の二つの値の平均をとってしまおう！　という発想である。⁵

長くなるので、わざわざ式は書き下さないが、これらを式(2)に代入すれば、離散式の完成だ。こちらも$y$, $z$成分も同様に求めることができる。
これで、式(1)と合わせて計6個の式が得られることになる。

1.3.4 1次元FDTD法の定式化

式が煩雑になるので、ここからは実装しようと考えている1次元で考えてみよう。

$\boldsymbol{E}$, $\boldsymbol{H}$の電磁波があった時、進行方向は、
$$
\boldsymbol{E} \times \boldsymbol{H}
$$
の方向になる。
つまり、電磁波が$z$方向に進行していたとして、電界は$E_x$成分のみとすると、磁界は$H_y$成分のみになるのである。

つまり、一次元FDTDでは、$E_x$, $H_y$以外は0と考えればよい。

例えば、式(1)の$y$成分は、
$$
\frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x} = - \mu \frac{\partial H_y}{\partial t}
$$
だが、$E_z$は0のため、
$$
\frac{\partial E_x}{\partial z} = - \mu \frac{\partial H_y}{\partial t}
$$

となる。⁶
面倒なので、今後は、$E_x=E$, $H_y=H$と書くことにし、$x$座標、$y$座標はもう必要ないため、$E(z) =E(i\Delta z)$のとき、$E[i]$のように表現することにしよう。
1次元のときのYee格子は

• $E$: 時間 $n\Delta t$, 空間 $[i]$
• $H$: 時間 $(n+1/2)\Delta t$, 空間 $[i+1/2]$

の値が与られる。

上の式(1)の$y$成分は

式(5): 1次元マクスウェル方程式(1)

$$
\frac{\partial {E}}{\partial z} = - \mu \frac{\partial H}{\partial t} \tag{5}
$$

であり、式(2)の$x$成分を考えると、

式(6): 1次元マクスウェル方程式(2)

$$
-\frac{\partial {H}}{\partial z} = \varepsilon \frac{\partial E}{\partial t} + \sigma E \tag{6}
$$

である。
この式(5)と式(6)の二つを利用して、1次元FDTD法を定式化する。

定式化も3次元と同様で、以下のようになる。導出は省略しよう。
$$
\frac{E^{n}[i+1] - E^{n}[i]}{\Delta z} = -\mu \frac{H^{n+1/2}[i+\frac{1}{2}] + H^{n-1/2}[i+\frac{1}{2}]}{\Delta t} \\
-\frac{H^{n+1/2}[i+1/2] - H^{n+1/2}[1-1/2]}{\Delta z} = \varepsilon \frac{E^{n+1}[i] E^n[i]}{\Delta t} + \sigma\frac{E^{n+1}[i] + E^n[i]}{2}
$$

アルゴリズムにしやすいように整理すると、次のようになる。

式(7): 1次元FDTD(1)

$$
H^{n+1/2}[i+\frac{1}{2}] = H^{n-1/2}[i+\frac{1}{2}]-\frac{1}{\mu}\frac{\Delta t}{\Delta z}(E^{n}[i+1] - E^{n}[i]) \tag{7}
$$

式(8): 1次元FDTD(2)

$$
E^{n+1}[i] = \frac{2\varepsilon - \sigma \Delta t}{2\varepsilon + \sigma \Delta t}E^{n}[i]-\frac{2}{2\varepsilon + \sigma \Delta t}\frac{\Delta t}{\Delta z}(H^{n+1/2}[i+\frac{1}{2}] - H^{n+1/2}[i-\frac{1}{2}]) \tag{8}
$$

これで、1次元FDTD法に必要な式は完成だ。

1.3.5 アルゴリズム

FDTD法のアルゴリズムは以下のようになる。

1. $E^{0}[i]$と$H^{-1/2}[i+1/2]$を初期条件として与える。
2. 式(8)で、$E^{1}[i]$を求める。
3. 式(7)で、$H^{1+1/2}[i]$を求める。
4. 2,3を繰り返して、$E^{n}[i]$と$H^{n+1/2}[i]$を交互に求めていく。
5. 適当なところで打ち切って終了

基本的に以上の流れで計算できる。
ただし、FDTD法は、格子分割するため、領域は無限に大きくできないことは注意しよう。
例えば、$[0]$から$[I]$までの有限領域で計算しようと思っても、$E^{n+1}[0]$を求めようとした場合、領域外の値$H^{n+1/2}[-1/2]$が必要になってしまう。
そのため、境界条件というものを与える必要がある。

例えば、完全導体(電界が0になる)で囲まれた空間内で計算しようと思った場合は、$E^{n}[i]$を計算する段階で、$E[0]=0$, $E[I]=0$とすればいいだけだ。

ただ、無限に大きい空間を計算しようとした場合は、ちょっと面倒で、吸収境界条件というものを使う必要がある。
これについては後々話そう。

3次元の場合も、上と同様の流れで計算していけばいい。

2. FDTD法で必要な技術・知識

上記まででFDTD法のアルゴリズムは実装できる。
しかし、FDTD法や電磁界解析には、その他いろんな技術や知識が必要になる。

ここからは、ちょっと雑多な情報になるかもしれないが、全て必須とも言える内容だ。
全部しっかりと理解しておこう。

2.1. フーリエ変換

3次元FDTDを解くと、時間$t$, 位置$x,y,z$, 電界$E_x, E_y, E_z$, 磁界$H_x, H_y, H_z$の10個のパラメータの関係が得られる。めちゃくちゃ多い…。
1次元だとしても、時間$t$, 位置$z$, 電界$E(=E_x)$, 磁界$H(=H_y)$の4つのパラメータの関係が得られることになる。

これから紹介するフーリエ変換をすると、更にそこに周波数$f$とその周波数成分の大きさ(スペクトル)$\mathcal{F}$が出てくる。
しかも、$\mathcal{F}$は複素数のため、より複雑だ。
こんなに多いと、混乱しそうだ。ここらへんを整理しながら、フーリエ変換について話していこう。

フーリエ変換はFDTD法で得られた結果を時間領域から周波数領域に変換する手法である。
これはどちらかというと、FDTD法というよりは、FDTD法で得られた結果から、その性質を見たりする、段階で使うものだ。

ところで、時間領域とか周波数領域とか何なんだ、と思う人もいると思う。
ここでは、あまり深いところには触れずに、簡潔に説明する。より深く知りたい人は、色んな所で詳しく解説されているのでフーリエ変換、離散フーリエ変換、高速フーリエ変換あたりのキーワードでググろう。

まず、以下のグラフを見てほしい。

(gif動画のため、動かない場合はクリック)

上図の上のグラフは$z$方向に進む電界である。単純化のため、速度1で進んでいるとしよう。
ここで、観測点(青い点線)の電界の時間変化を見てみよう。これが、下のグラフである。
この下のグラフが時間領域のグラフである。

さて、次に周波数領域について見てみよう。
時間領域のグラフに注目して、(観測点上の)電界$E$は時間$t$の関数のため、$E(t)$と表す。

これを、以下の式でフーリエ変換しよう。⁷
$$
\mathcal{F}(f) = \int_{-\infty}^\infty E(t)e^{-j t (2\pi f)}dt
$$
$j$は虚数単位である。ここで、$2\pi f=\omega$とすれば、

式(9): フーリエ変換

$$
\mathcal{F}(\omega) = \int_{-\infty}^\infty E(t)e^{-j t \omega}dt
$$
となる。この$\omega$は角周波数と呼ばれる。
わりと、角周波数もまとめて周波数と表現している資料も多いが、ここでは厳密に言葉を使い分けていこう。つまり、周波数と言ったら$f$だし、各周波数と行ったら$\omega$である。

式(9)は、無限領域で積分しているが、例えば、上のグラフでは、$0.5\le t \le 0.6$以外では0になっている。
実際にフーリエ変換が必要になるときも、大抵は十分な時間$T$がたてば0になると考えることができて、
$$
\mathcal{F}(\omega) = \int_{0}^{T} E(t)e^{-j t \omega}dt
$$
とすることができる。
これを以下のように数値積分することで、コンピュータでも簡単に計算できるようになる。

式(10): 離散フーリエ変換(DFT)

$$
\mathcal{F}(\omega) = \sum_{n=0}^{N} E(n\Delta t)e^{-j n\Delta t \omega}\Delta t \tag{10}
$$
ここで、$t=n\Delta t, T=N\Delta t$で、$N$は分割数である。
この式を、離散フーリエ変換(Discrete Fourier Transform; DFT)という。⁸⁹

それでは、実際にフーリエ変換してみよう。

下のグラフがフーリエ変換した結果だ。
え？　こんなの見ても全然意味がわからない？　その意味をこれから書いていこう。

まず、電界(上のグラフ)が

$$
E(t) = A_0\cos{(0\times 2\pi)} + A_1\cos{(1\times 2\pi)} + A_2 \cos{(2\times 2\pi)} + \cdots + A_N(\cos{N \times 2\pi}) \\
+B_1\sin{(1\times 2\pi)} + B_2\sin{(2\times 2\pi)} + \cdots + B_N\sin{(N\times 2\pi)}
$$

という感じで、三角関数の和で表すことができると仮定しよう。
すると、この$A_n, B_n$(の絶対値)が大きければ大きいほど、「その周波数の波が強い」ことがわかる。
たとえば、上の図だったら、なんとなく周波数$f=8$くらいまでは結構強くて、それ以降は徐々に小さくなっている、という事が見て取れる。つまり、なんとなく低い周波数の寄与が大きいのかな～って事がわかるだろう。
また、よくみると、ある周波数の周期で($15, 25, 35,\cdots$)で強くなっていることがわかる。
これも重要な特性の一つだ。

前述の通り、スペクトルは複素数の値を取る。
このとき、実部が$\cos$成分の大きさ$A_n$を表し、虚部が$\sin$成分の大きさ$B_n$に対応しているのである。

たとえば、周波数$f=1$のとき、スペクトルはだいたい、$-0.19+6.02j$である。
つまり、$A_1 = -0.19$, $B=6.02$ということがわかるのである。
以下のアニメーションを見てほしい。もとの形状に戻っていくことがわかるだろうか？

(gif動画のため、動かない場合はクリック)

その周波数帯の強さ(振幅)は、$\sqrt{(\text{Re}(\mathcal{F}))^2+(\text{im}(\mathcal{F}))^2}$で表される。$\text{Re}$と$\text{Im}$は実数部と虚数部を表す。
これは、三角関数の合成公式
$$
A_n\cos{(2\pi ft)}+B_n\sin{(2\pi ft)} =\sqrt{A_n^2 + B_n^2} \sin{(2\pi t + \alpha)}
$$
から、理解できる。また、$\alpha$に対応するところは、位相差と呼ばれたりする。
普通、グラフ化するときは、実数部・虚数部よりも、こちらの強さと位相差のほうで表す。

で、こんな事がわかって何が嬉しいか？　という話だが、電磁波というものは周波数帯で全然違う振る舞いをする。
例えば、電磁波が鉄などの導体にぶつかった時、周波数が高ければ高いほど、急激に減衰する表皮効果というものがある。
もし、とても薄い鉄の板に、低い周波数から高い周波数まで一様に分布している電磁波があったとすると、高い周波数成分は減衰して消えてしまい、低い周波数成分のみが残るのである。
つまり、外部に電磁波を漏らさないように、鉄で囲っていたとして、もし低い周波数成分が強い！　と言うことがわかると、外に漏れてしまうことが解る。
逆に、高い周波数成分が大きくて、低い周波数成分がそうでもないなら、鉄の板をもっと薄くして軽量化できそう、などいろいろなことが考えられるのだ。

それ以外にも、周波数が重要なのはアンテナだ。
例えば、bluetoothなんかは、2.4GHz帯の電磁波で通信しており、この2.4GHz帯に対応する周波数帯はどこまで届いているだろうか？　ということを計測できるのである。
最近流行りの5Gなんかも同様だ。

ところで、FDTD法は電界$E^t[i]$や磁界$H^t[i]$を順次解いていく方法のため、式(10)のDFTをFDTDをときながら計算していくことが可能である。
ここに関しては、実装のところで詳しく話そう。

2.2 波源モデル(初期値)

FDTD法を解くには、初期値を与える必要がある。これは、実際の波源に対応する。
初期値として強さ1の電界を与えるとなると、どうやって与えればいいだろうか？

最も簡単なのは、一つの適当な要素に1という値を与えて、それ以外を0にする、つまり、デルタ関数として与える手法だ。
実際には、格子分割して有限領域で求めるため、上のフーリエ変換で説明したような、矩形パルスになる。

しかし、このように設定してしまうと、分割数が大きくなるにつれて高周波成分が大きくなってしまうし、上のパルスを例に取ると、上述したように$f=15, 25, 35,\cdots$と、周期的に強い成分ができているし、逆に、$f=20, 30, 40,\cdots$あたりはほぼ0になっていることがわかる。
このようなスペクトルのパルスを波源とすると、観測したスペクトルが波源によるものなのか、その他の要素のものなのかがわからず、不便な時がある。

そこで、よく使われるガウスパルスについて説明しよう。

ガウスパルスは、ガウス関数
$$
p(\tau) = \exp{\left(-\frac{\tau^2}{2w^2}\right)}
$$
ここで、$w$はパルス幅である¹⁰。
$\tau$は光の速度を$c$として、$\tau = t + z/c$で計算できる。$z/c$の項は、速度$c$で進んでいる電磁波が位置$z$のとき、$z=0$よりどれくらい進んでいるかを表す。

$\tau=0$の時、ガウスパルスは最大値となり$p(\tau)=1$になる。

ガウスパルスの便利なところは、フーリエ変換してもガウス関数のままである。

上の図は、ガウスパルスとそのスペクトル特性だ。
スペクトル特性は、ガウス関数の右半分になっている。
$w$を小さくすればするほど、ガウスパルスはより鋭く、スペクトル特性は、よりなだらかに、高い周波成に伸びていく。
逆に、$w$を小さくすると、ガウスパルスはなだらかになり、スペクトル特性は、低い周波帯のみになっていく。

例えば、1GHzオーダーの解析をしたいのであれば、ガウスパルスの幅を狭く鋭くして、高周波成分にまで残せばいいし、逆に1MHzオーダーだったら、ガウスパルスの幅はなだらかにして、余計な高周波成分を含めないようにすればいい。

ガウスパルスの他にも、ガウス関数を微分した微分ガウスパルスや、sin波で変調したガウスパルスなどもよく用いられる。これらは、ガウスパルスと違い、平均が0になるため、直流成分をなくしたい場合によく用いられる。

なお、ガウス関数は、$\tau$をどれだけ大きくしても、0にならない。
そこで、以下のようにある程度の時刻で打ち切るのが普通である。

$$
p(\tau)=\left\{
\begin{array}{ll}
\exp{\left(\displaystyle-\frac{\tau^2}{2w^2}\right)} & \text{if }(-N_t w < \tau < N_tw)\\
0 & \text{else}
\end{array}
\right.
$$

ここで、$N_t=2$のときは、$p(\tau)\simeq 0.135=(13.5 \%)$, N_t=3のときは、$p(\tau)\simeq 0.0111(=1.11 \% )$となる。
適切なところで打ち切ればいい。

2.3 吸収境界条件

次に、吸収境界条件について解説する。

今後は、簡単のため、以下の図の1次元Yee格子をもとに説明していく。

前述したとおり、FDTD法は、有限領域を解析するため、領域の境界には境界条件を設定する必要がある。

境界条件には、以下のように、境界をある値と固定するというものだ。
$$
E^{n}[0]=A, E^{n}[I] =B
$$
このような条件をディリクレ境界条件という。
$A=B=0$のとき、これは完全導体に囲まれていることに値する。

しかし、電磁界解析では、無限に広いと仮定して解析したいことも多い。
このとき、使うのが吸収境界条件だ。

一言でいうと、吸収境界条件は「電磁界が反射しないような境界」である。
この吸収境界条件はいろいろな手法が提案されており、Murの吸収境界とBerengerのPML吸収境界が有名だ。
Murの吸収境界は実装が簡単だが精度に劣る。一方、BerengerのPML吸収境界は高精度だが、実装が難しく、また少しメモリも必要になる。
他にも、LiaoやHigdonの吸収境界などがあるが、上2つと比べるとあまり有名じゃないのかもしれない。詳しく知りたい人は、参考文献の「FDTD法による電磁界およびアンテナ解析」を読もう。

今回は、Murの吸収境界条件だけ説明する。
一次元FDTDではこれで十分だ。

2.3.1 Murの吸収境界条件

電磁波というのは、光速$c$の速度で進んでいる。
それならば、境界の外側へそのまま光速$c$で進んでいるという式を立てる、というのが、Murの境界条件のアイデアだ。

ところで、FDTD法の定式化をしたマクスウェル方程式に速度$c$が出てこないのに、突然$c$が出てくるのに不安を覚える人もいるかも知れない。これについて簡単に解説しよう。
マクスウェル方程式をいじると、$1/\sqrt{\varepsilon_0 \mu_0}$の速度で電磁波が伝搬しているというのがわかる。
そして、この値が光速と同じだ！　ってことから、
$$
c = \frac{1}{\sqrt{\varepsilon_0 \mu_0}}
$$
と言うことが解ったという歴史がある。
そのため、直接は見えないが、真空の誘電率と真空の透磁率の値で関連付けられるのである。

まず、上の一次元Yee格子の、$E^n[I]$側について考えよう。
もし、$E$が(形状を変えずに)速度$c$で$z$方向に進んでいるとすれば、
$$
E(z,t)=E(z-ct)
$$
と表すことができる。
なぜか？　まず、位置が$z$, 時間が$t_0$の電界は$E(z,t_0)$と表せる。
それから$t$たったとき、$E(z,t_0)$は$z$方向に$ct$進んでるのだから、
$$
E(z+ct,t_0+t) = E(z,t_0)
$$
と表せるだろう。さらに、$z+ct\rightarrow z$, $t_0+ t \rightarrow t$と置き直そう。すると、
$$
E(z,t) = E(z - ct, t_0)
$$
となるが、$t_0$は定数のため、結局、$E(z,t) = E(z-ct)$となる。

つまり、$E(z,t) = E(z-ct)$を境界条件としてしまえば、形状を変えず、反射もせず電磁波は境界を通り抜けていくことになる。

上の式の両辺を$z$で偏微分してみよう。
$$
\begin{array}{ll}
\displaystyle\frac{\partial E(z,t)}{\partial z} & =\displaystyle\frac{\partial E(z-ct)}{\partial z}\\
&\displaystyle= \frac{\partial (z-ct)}{\partial z} \frac{\partial E(z-ct)}{\partial (z-ct)} \\
&\displaystyle= \frac{\partial E(z-ct)}{\partial (z-ct)}
\end{array}
$$
次に、$t$でも偏微分する。
$$
\begin{array}{ll}
\displaystyle\frac{\partial E(z,t)}{\partial t} & =\displaystyle\frac{\partial E(z-ct)}{\partial t}\\
&\displaystyle= \frac{\partial (z-ct)}{\partial t} \frac{\partial E(z-ct)}{\partial (z-ct)} \\
&\displaystyle= -c\frac{\partial E(z-ct)}{\partial (z-ct)}
\end{array}
$$
二つの式の右辺が似ていることに注目して、代入すると、
$$
\displaystyle\frac{\partial E(z,t)}{\partial t} = -c\displaystyle\frac{\partial E(z,t)}{\partial z}
\iff \displaystyle\frac{\partial E(z,t)}{\partial z} +\frac{1}{c}\displaystyle\frac{\partial E(z,t)}{\partial t} = 0
$$

という式が得られる。

あとは差分近似して、
$$
\frac{E^{n-\frac{1}{2}}[I] - E^{n-\frac{1}{2}}[I-1]}{\Delta z} + \frac{1}{c}\frac{ \displaystyle E^{n}[I-\frac{1}{2}] - \displaystyle E^{n-1}[I-\frac{1}{2}] }{\Delta t} = 0
$$
となり、FDTD法の導出のときと同じく、Yee格子に合わせるために、平均を取り、

\begin{multline}
\frac{1}{\Delta z} \left( \frac{E^{n}[I] + E^{n-1}[I]}{2} - \frac{E^{n}[I-1] + E^{n-1}[I-1]}{2} \right) \\\\
+ \frac{1}{c\Delta t}\left( \frac{E^n[I]+E^n[I-1]}{2} - \frac{E^{n-1}[I]+E^{n-1}[I-1]}{2} \right) = 0
\end{multline}

となり、整理すると、

式(11): Murの吸収境界(+z側)

$$
E^n[I] = E^{n-1}[I-1] + \frac{c\Delta t - \Delta z}{c\Delta t + \Delta z}\left( E^n[I-1] - E^{n-1}[I] \right) \tag{11}
$$
となる。これがMurの吸収境界である。
$E^{n}[0]$の方は、電磁波の進行方向が逆のため、
$$
E(z,t)=E(z+ct)
$$
として計算すれば良い。
計算過程は省略するが(自分で計算してみよう)、結果は、式(11)の$I\rightarrow 0$, $I-1 \rightarrow 1$と書き換えた、

式(12): Murの吸収境界(-z側)

$$
E^n[0] = E^{n-1}[1] + \frac{c\Delta t - \Delta z}{c\Delta t + \Delta z}\left( E^n[1] - E^{n-1}[0] \right) \tag{12}
$$
となる。

Murの吸収境界は、一次元の問題の時、(差分近似の部分を除けば)反射係数は0になり、理想的な吸収境界となる。これが、上で一次元の場合はMurの境界条件で十分といった理由だ。
しかし、二次元以上の問題のときは、境界に対して斜めに入射してくる電磁波も無数に存在する。
このような問題では、$E=E(z-ct)$は成り立たず、精度が落ちてしまう。
その誤差が無視できない場合には、PML吸収境界など、他の吸収境界を使った方が良い。

2.4 Courantの安定条件(収束条件)

FDTD法は$\Delta z$を小さくすると、$\Delta t$も小さくしないと発散して正しくない解になる可能性がある。
結論から言おう。1次元FDTDのCourantの安定条件は、

式(13): Courantの安定条件(1次元)

$$
c\Delta t \le \Delta z \tag{13}
$$
である。ここで$c$は速度を表す。

この式の意味を簡単に説明しよう。

電磁波が速度$c$で進んでいるとする。
すると、$t=n \Delta t$のとき、$z=0$の電界$E^n[0]$は、$t=(n+1)\Delta t$は$z=c\Delta t$まで進んでいることになる。

ここで、例えば、Courantの安定条件を満たさない、$c \Delta t = 2\Delta z$と設定したとする。
すると、速度$c$は
$$
c=(2\Delta z)/\Delta t
$$
となる。

以下にイメージしやすいように図を用意した。

Yee格子で表すと、$\Delta t$秒で、$z=[0]$から$z=[2]$まで進むことになる(図の白抜きの赤丸の部分)。

しかし、1次元FDTDの式(7)と式(8)を思い出すと、$E^{n}[0]$を利用して求めることができる値は、$H^{n+1/2}[1/2]$(図の$t=(n+1/2)\Delta t$の行の赤丸)である。
さらに、$H^{n+1/2}[1/2]$で求めることができる値は、$E^{n}[1]$(図の$t=(n+1)\Delta t$の行の赤丸)である。

つまり、どう頑張っても、式(7)と式(8)では、$E^{n}[0]$にある値を$E^{n}[2]$まで送ることができないのだ。

そのため、Courant安定条件$c \le \Delta z/\Delta t$を満たさなければおかしな値が得られることになるのだ。
以上は感覚的な説明で、もう少し厳密な証明もあるが、ここでは説明しない。

なお、1次元FDTDではなく、$d$次元FDTDでは

式(14): Courantの安定条件(d次元)

$$
c\Delta t \le \frac{\Delta}{\sqrt{d}} \tag{14}
$$
である。ただし、上式は立方(正方)格子の場合で、$\Delta=\Delta x = \Delta y = \Delta z$である。

3. 実装

説明がとても長くなってしまったが、実装に入ろう。

解く問題は何度も行っているように一次元の電磁界だ。
今回は、$z$方向に進む平面波($xy$平面に一様な波)がどんなふうに動くか？
というシミュレーションをしてみよう。

ソースコードは以下のURLにおいた。

必要なパッケージはnumpy, scipy, matplotlibがあれば動くと思う。

• atily17/FDTD1d(github)

以下では、このソースの中から抜粋して重要なところを解説していく。

3.1 問題設定

上述の通り、平面波のシミュレーションをする。
パルスは、わかりやすくガウスパルスとしよう。
以下のようなグラフの電界と磁界を右に移動させていくイメージだ。

ガウスパルスは大きさ、パルス幅などを設定できるようにし、
境界条件も、ディリクレ境界条件と吸収境界条件の二つを与えられるようにしよう。

また、媒質の位置やパラメータも与えられるようにする。
これで結構面白い波の波形が得られる。

3.2 Yee格子

まず、Yee格子を作っていこう。

Yee格子の前に、媒質(Medium)クラスを作ることにした。

class Medium:
    def __init__(self):
        self.setMediumInit()
        self.left = 0
        self.right = 0

    def setRegion(self, left, right):
        self.left=left
        self.right=right

    def setMediumInit(self):
        self.eps = constants.epsilon_0
        self.mu = constants.mu_0
        self.sgm = 0

    def setMedium(self, reps, rmu, sgm):
        self.setMediumInit()
        self.eps = reps*self.eps
        self.mu = rmu*self.mu
        self.sgm = sgm

これらで、媒質の各パラメータを設定できる。
簡単に説明すると、self.epsは真空の誘電率を初期値として入れている。透磁率も同様だ。
そして、setMedium関数に、reps(比誘電率$\varepsilon_r$), rmu(比透磁率$\mu$), sgm(導電率)を引数にとって、誘電率$\varepsilon=\varepsilon_r \varepsilon_0$, 透磁率$\mu = \mu_r \mu_0$と、導電率$\sigma$の値を設定できる。

そして、YeeGridクラスを作ろう。
まず、メンバ変数は以下の通りだ。

class YeeGrid:
    def __init__(self, directory="output"):
        self.domain = 0 # 領域長さ
        self.n = 0      # 分割数
        self.dz = 0     # 空間分割幅
        self.dt = 0     # 時間分割幅
        self.mediums = []   # 媒質
        self.e_element = np.array([])   # E要素の座標(numpy)
        self.h_element = np.array([])   # H要素の座標(numpy)

まず、解析する幅だが、左の境界の幅は$z=0$で固定して、右の境界だけを設定するようにした。
これは、初期値なのでとりあえず全部に0を入れている。
媒質はリストで、Mediumオブジェクトを追加していく。

次に、グリッドを作る関数は以下の通りだ。

class YeeGrid:
    def setGrid(self, domain, N, ratio=0.99):
        self.domain = domain
        self.n = N
        self.dz = self.domain/self.n
        self.e_element=np.linspace(0,self.domain, N+1)
        self.h_element=np.linspace(0+self.dz/2,self.domain-self.dz/2, self.n)
        c = constants.c
        self.dt = self.dz/c*ratio

引数のdomainに領域の範囲で、Nは分割数である。
これで、self.e_elementにはN+1個の要素が、self.h_elementにはN個の要素ができる。

setGrid関数でself.e_elementとself.h_elementに各格子の要素の座標を入れる。
self.e_element=np.linspace(0,self.domain, N+1)と
self.h_element=np.linspace(0+dz/2,self.domain-dz/2, N)の部分がそれだ。
磁界$H$の格子は$E$から半分ずらすため、格子間隔dzの半分をずらしている。
なお、self.e_elementの要素数がself.h_elementの要素数より1多いのも、Yee格子の図を見れば容易にわかる。

たとえば、setGrid(1, 10)なら、

self.e_element = [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
self.h_element = [0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95]

となる。
また、ここでCourantの安定条件から、$\Delta t$も決定した。
式は、以下の通りだ。
$$
\Delta t = r\frac{\Delta z}{c}
$$
理屈上、$\Delta t = \Delta z / c$でもいいのだが、丸め誤差により発散する可能性もあるため、$r=0.99$として、それよりも小さくした。

以下のaddMedium関数は、媒質を追加する関数だ。

class YeeGrid:
    def addMedium(self, left, right, reps, rmu, sgm):
        m = Medium()
        m.setRegion(left, right)
        m.setMedium(reps, rmu, sgm)
        self.mediums.append(m)

特に説明することはないと思う。

3.3 FDTDソルバー

それでは、FDTDのメイン箇所を書いていこう。

3.3.1 初期条件

まず、初期条件を作っていこう。
Pulseクラスを作った。メンバ変数は以下の通りだ。

class Pulse:
    def __init__(self):
        self.w = 100e-12    # パルス幅
        self.threshold = 4  # 打ち切る位置
        self.E = 1          # パルスの大きさ

self.thresholdはパルスの式
$$
p(\tau)=
\left\{
\begin{array}{ll}
\exp{\left(\displaystyle-\frac{\tau^2}{2w^2}\right)} & \text{if }(-N_t w < \tau < N_tw)\\
0 & \text{else}
\end{array}
\right.
$$
の$N_t$に値する。

これの計算は、YeeGridオブジェクトを引数とした。

class Pulse:
    def calc(self,grid):
        c = constants.c
        mu0 = constants.mu_0
        eps0 = constants.epsilon_0
        impedance = np.sqrt(mu0/eps0)
        z0 = c*self.w*self.threshold+2*grid.dz
        taues = (grid.e_element-z0)/c
        tauhs = grid.dt/2 + (grid.h_element-z0)/c
        return {"E": self.E*self._gaussian(taues, self.w), 
                "H":self.E
/impedance*self._gaussian(tauhs, self.w)}

ここで、注意すべきは、$p(\tau)$の$\tau$だ。$\tau = t + z/c$と書いたが、まず、$z$は左境界を0にしたため、そのままでは、ガウス関数の右半分しか計算できない。
そのため、$p(\tau)=0$となるように、$z$を$N_tw$だけずらす必要がある。
そこで、$z \rightarrow z-N_tw$と設定する。
上のプログラムでは、さらに$2\Delta z$分ずらしているが、ぴったしだと丸め誤差などから、意図通りにならない可能性があるからだ。

self._gaussian関数は、上の式を再現しているだけで、ここでは載せない。

3.3.2 1次元FDTDの式(7), (8)を実装する。

次に、FDTDを実装しよう。

以下の流れになっている。

class Solver:
# __init__などは省略
    def calc(self):
        self.calcInit() # 初期値
        self.save(0)    # そのときの電界、磁界をセーブ
        self.time = 0   # 今の時間
        self.calcFourier()  # フーリエ変換の計算
        for t in range(1, self.NTime+1):
            self.nextstep(t)    # 次のステップへ
            self.save(t)        # そのときの電界、磁界をセーブ
            if t % 10 == 0:
                print("iter="+str(t)+", ", end="")
                if self.judgeConvergence(bPrint=True):  # 収束判定
                    break

文字で書くと次のとおりだ。

1. まず、初期値を決める。
2. FDTD法の計算式で次のステップへ
3. 電界、磁界の値のセーブ
4. ループ10回につき収束判定をする。収束しているなら、ループから抜ける。

フーリエ変換とセーブの部分は後々説明する。

次に、2.のところを細かく見てみよう。

class Solver:
    def nextstep(self, t):
        preE = self.E.copy()    # 電界Eのコピー
        self.calcH()            # Hを式(7)から計算
        self.time += self.grid.dt/2   # dt/2 時間足す
        self.htimes.append(self.time) # 計算に使った時間を保管
        self.calcFourier("H")   # Hのフーリエ変換
        self.calcE()            # Eを式(8)から計算
        self.calcBorder(preE)   # Eの境界条件を計算
        self.time += self.grid.dt/2   # dt/2 時間足す
        self.etimes.append(self.time) # 計算に使った時間を保管
        self.calcFourier("E")   # Eのフーリエ変換

self.Eとself.Hが電界、磁界の値が入ったnumpy配列である。
重要なのは、calcH()とcalcE(), そしてcalcBorder()で、これが式(7)と式(8)と境界条件に値する。
最初に電界Eをコピーしているのは、Murの吸収境界条件、式(11), 式(12)で、「一つ前の時間の電界」の値が必要だからだ。

calcH()は以下のようになる。

class Solver:
    def calcH(self):
        mu = self.properties["mu"]
        dt = self.grid.dt
        dz = self.grid.dz
        c1 = 1/mu*dt/dz
        self.H = self.H - c1 * (self.E[1:] - self.E[:len(self.E)-1]) # 式(7)

ここで、mu = self.properties["mu"]の、self.propertiesには各要素の媒質の透磁率や誘電率、導電率などの値がnumpy配列で入っている。
高速化のため、forループではなく、numpyで一括計算している。式(7)と見比べてみよう。

calcEは以下の様になる。

class Solver:
    def calcE(self):
        I=len(self.E)
        eps = self.properties["eps"][1:I-1]
        sgm = self.properties["sgm"][1:I-1]
        dt = self.grid.dt
        dz = self.grid.dz
        c1 = (2*eps-sgm*dt)/(2*eps+sgm*dt)
        c2 = 2*dt/((2*eps+sgm*dt)*dz)
        self.E[1:I-1] = c1*self.E[1:I-1] - c2*(self.H[1:]-self.H[:len(self.H)-1]) # 式(8)

これも、式(8)と見比べておこう。

そして、最後に境界条件を求めよう。

class Solver:
    def calcBorder(self, preE):
        c = constants.c
        dz = self.grid.dz
        dt = self.grid.dt
        if self.borderL == BoundaryCond.ABSORB:
            c1 = (c*dt - dz)/(c*dt+dz)
            self.E[0] = preE[1] +  c1 * (self.E[1] - preE[0]) # 式(12)
        elif self.borderL == BoundaryCond.DIRICHLET:
            self.E[0] = 0

        I = len(self.E) - 1
        if self.borderR == BoundaryCond.ABSORB:
            c1 = (c*dt - dz)/(c*dt+dz)
            self.E[I] = preE[I-1] +  c1 * (self.E[I-1] - preE[I]) # 式(11)
        elif self.borderR == BoundaryCond.DIRICHLET:
            self.E[I] = 0

ここで、引数のpreEは一つ前の電界だ。

このプログラムは、吸収境界条件と、ディリクレ境界条件を表している。
self.borderLは左境界、self.borderRは右境界で、EnumクラスBoundaryCondに、ABSORB, DIRICHLETから判断できる。

まず、BoundaryCond.ABSORBは吸収境界条件である。
これは、式(11), 式(12)を再現している。

そして、BoundaryCond.DIRICHLETはディリクレ境界条件で、右辺を0にしている。
これは、Eの値を0にしているだけだ。

3.3.3 収束判定する

吸収境界条件でFDTDを回し続けていると、すべての要素でE,Hの値が0になる。
このあとも、FDTDを回し続けるのは無意味なので、ある程度収束したら、打ち切るようにしよう。

このために、以下のようなjudgeConvergence関数を作った。

class Solver:
    def judgeConvergence(self, cc = 1e-3, bPrint=False):
        emax = self.pulse.E
        eave = np.linalg.norm(self.E, 1)/len(self.E)
        conv = eave/emax
        bConvergence = conv < cc
        if bPrint:
            print("time="+str(self.time)+", ave/max="+str(conv))
        return bConvergence

簡単にいうと、最大値と平均の比がある程度以下になったら打ち切るというものだ。
$$
\frac{\text{max}({E^0}[k])}{\text{average}(E^t[k])} < \epsilon
$$
最大値は、$t=0$のときの最大値、平均は今の時間の平均値だ。

このへんの収束判定条件に関しては自分の好きなように設計しても良いとおもう。

3.3.4. フーリエ離散変換

フーリエ離散変換の式を再掲しよう。
$$
\mathcal{F}(\omega) = \sum_{n=0}^{N} E(n\Delta t)e^{-j n\Delta t \omega}\Delta t \tag{10}
$$

class Solver:
    def calcFourier(self, eh=""):
        j = 1j
        f = np.exp(-j*self.time*2*np.pi*self.freqs)*self.grid.dt # 式(10)
        for observe in self.observes:
            if eh != "H":
                self.fourierE[observe["name"]] += self.E[observe["E"]]*f
            if eh != "E":
                self.fourierH[observe["name"]] += self.H[observe["H"]]*f

ここで、self.freqsは計算する周波数(numpy)、self.observesは観測点(リスト)だ。
今回は、全要素で計算せずに、観測点で設定した部分だけをフーリエ変換することにした。

各ステップごとに、4行目の式(10)を計算して、どんどん足して最終的に勝手に計算
引数ehは、Eに設定すれば電界のフーリエ変換を、Hに設定すれば磁界のフーリエ変換をする。

他にも細かいことは色々あるが、ここでは省略しよう。
詳しくはソースコードを見てほしい

3.4 解析結果

計算結果が得られたら、次は可視化だ。
これについては、matplotlibが基本になるので解説は省略しよう。

ここでは、解析結果の例をあげる。

3.4.1 媒質なし、両側とも吸収境界条件

解析領域は、0~1、分割数は500で計算した。
両方の境界条件は吸収境界条件にし、媒質は置かない。

N=500

grid = YeeGrid.YeeGrid()
grid.setGrid(1, N)

solver = Solver.Solver(grid)
solver.setInit(1,100.0e-12)
solver.setBorder("absorb", "absorb")
solver.calc()

graph = Graph.SpaceGraph(solver)
graph.animate(intervl=50)

graph = Graph.FourierGraph(solver)
graph.plotE(bSave=True)

下のGraphがグラフを出力するモジュールだ。

結果を見てみよう

• アニメーション

(gif動画のため、動かない場合はクリック)

左側の軸が電界の値、右側の軸が磁界の値だ。

いい感じに流れているのが見て取れる。右境界で吸収されているのもわかるだろう。

• 周波数数特性(観測点: 0.5)

こちらもいい感じだ。
前に説明したように、ガウス関数のフーリエ変換はガウス関数だ。
ログスケールのためちょっと見にくいかもしれないが、滑らかになっているのがわかる。

3.4.2 媒質なし、右側境界をディリクレ境界条件

(gif動画のため、動かない場合はクリック)

上手く反射しているのがわかる。
上のグラフより速度が早いように見えるが、Qiitaの制限である1MB以下にしようとしたらこうなってしまった。
同じ速度であることに注意しよう。

3.4.3 媒質あり、両側境界を吸収境界条件

比誘電率$\varepsilon_r=3$の媒質を中心に入れてみたものだ。
ガウスパルスが誘電体に入る動きが面白い

(gif動画のため、動かない場合はクリック)

4. まとめ

今回はFDTD法の基礎について解説した。
FDTD法のアルゴリズムと、フーリエ変換、波源モデル、そしてMurの吸収境界条件を解説し、そしてこれらを用いてpythonで実装した。

まだまだ2次元以上のFDTD法の実装やPML吸収境界条件、遠方界の解析など、書ききれないことが多いが、ひとまずここまでにしておこう。予定より遥かに長くなってしまった。

筆者もまだ学んで数ヶ月のため、まだまだ勉強中だ。
なにか指摘があれば、コメントしてほしい。

参考文献

• FDTD法による電磁界およびアンテナ解析
  • 式展開などをかなり細かく(ちょっと冗長に感じるほど)してくれるので、これ一冊読めば、FDTD法はだいたい理解したと思っていいと思う。
• 株式会社EEM
  • OpenFDTDの説明や、EEMドキュメントの理論説明書, セミナー資料はめちゃくちゃ参考にさせてもらった。
  • 株式会社と書いているが、社員たった一人らしい。それで４つの製品と４つのOSS(フリーソフト)を出してるのはすごい。OpenFDTDはフリーのFDTDソフトとしてかなり有名なソフトらしい。

1. $\nabla \cdot$や$\nabla \times$の説明もしようと思ったが、冗長になると思ったのでやめた。 ↩

2. 手元にある電磁気学の本にあまり詳しい説明がなかったのだが、おそらくこれも誘電率と同様に近似式だと思う。 ↩

3. この図を作るのめっちゃ頑張った(どうでもいい) ↩

4. 中心差分は二次精度と呼ばれ、$\Delta x$が10倍小さくなると、誤差が$10^2$倍小さくなる。これは厳密にはテイラー展開で導出できる。 ↩

5. これも二次精度であることをテイラー展開で導出できる。 ↩

6. マクスウェル方程式で$x$,$y$方向に一様と考えてもいい。電磁気学の知識を使わず、式のみで考えたいときは、こっちで考えたほうがわかりやすい。 ↩

7. なぜ式(9)がこんな形になるかは説明しないが、色んな所で解説されてるはずだ。 ↩

8. このとき、$\Delta t$が省略されている流儀が多いが、今回は直感的なわかりやすさを重視して、$\Delta t$を利用しよう。全体のスケールが変わるである。 ↩

9. さらに、DFTを高速にした高速フーリエ変換(FastFourierTransform;FFT)もある。簡単に概要を説明すると、計算の順番をうまいこと工夫して、DFTの計算量を減らしたものだ。 ↩

10. いろいろ調べたが、パルス幅というのはいろんな定義があるようだ。こちらのサイトを参考に、RMSwidthを採用した。 ↩