はじめに
数値計算は楽しいので、数値計算のプログラムを書きます。
今回は、主として電磁界などの解析に用いられる、時間領域有限差分（Finite Difference Time Domain, FDTD）法を試したいと思います。
弁明
本記事は、何となく高く見える数値計算・物理シミュレーションの簡単なご紹介と出来ればと思ってますので、厳密な理論・精度は追求しません。
理論・語句用法の間違いへのご指摘は、優しい言葉で沢山いただけるとありがたいです。
また、本記事の投稿にあたっては、全編通して、『数値電磁界解析のためのFDTD法- 基礎と実践』[1]を参考にしております。詳しい内容が知りたい方は、当該書籍をご参照いただければと思います。
FDTD法とは
FDTD法は、時間領域の陽解法を用いて、計算空間内の（主に）波動方程式を解く手法です。
もう少し詳しく説明すると、
(1)計算空間を以下のように有限に小さい領域で離散化し、
(2)波動方程式を時間に対する漸化式とすることで、
各時間ステップごとに任意の位置の数値解を求める方法です。
まさしく数値計算らしい数値計算手法で、物理シミュレーションというと、こういったものを想像されるのではないでしょうか。実際に、時間領域での陽解法は、他領域でも幅広く使用されています。
一般的に、波動方程式等は解析解を求めることが非常に難しいことが多いため、こうした数値計算手法は有効です。
数値計算の準備
楽しい数値計算に入る前には、少し辛い作業を済ませなければなりません。
それは、数値計算に用いる数式の準備です。
しかしながら、多くの場合はインターネットに情報が転がっているのでとりあえずコピペして実行しちゃいましょう。
ここでは、勉強のため、簡単に式展開の概要だけ説明いたします。厳密な話は他の文献・書籍を参照いただければと思います。読み飛ばしていただいて問題ありません。
Maxwell方程式の展開
この世の電磁界理論は、全てMaxwell式で説明されますので、FDTD法でも漏れなく以下のMaxwell方程式$(1), (2)$をガチャガチャして使います。
\begin{aligned}
\nabla \times \boldsymbol{E}(t, x) + \mu_0 \frac{\partial \boldsymbol{H}(t, x)}{\partial t} &= 0 &…(1)\\
\nabla \times \boldsymbol{H}(t, x) - \epsilon_0 \frac{\partial \boldsymbol{E}(t, x)}{\partial t} &= 0 &…(2)
\end{aligned}
ここで$\boldsymbol{E}$は電界、$\boldsymbol{H}$は磁界、その他記号は電気定数を表すもので、$\epsilon$は誘電率、$\mu$は透磁率、$\sigma$は導電率です。$\nabla$はベクトル表記の微分演算子ですね。
$\nabla$および$\boldsymbol{E}, \boldsymbol{H}$について展開すると以下の式$(3) - (8)$を得ます。ここでは簡単のため、$y, z$方向には一様であるとして、$x$方向のみ考えましょう。
\begin{aligned}
\mu_0 \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z} - \frac{\partial E_z}{\partial y} &…(3)\\
\mu_0 \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z} &…(4)\\
\mu_0 \frac{\partial H_z}{\partial t} = \frac{\partial E_x}{\partial y} - \frac{\partial E_y}{\partial x} &…(5)\\
\epsilon_0 \frac{\partial E_x}{\partial t} = \frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z} &…(6)\\
\epsilon_0 \frac{\partial E_y}{\partial t} = \frac{\partial H_x}{\partial z} - \frac{\partial H_z}{\partial x} &…(7)\\
\epsilon_0 \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} &…(8)
\end{aligned}
頭がいたいですね。
今回は一次元($x$成分のみ)で考えるので$y, z$成分の変化は$0$として以下の式$(9) - (14)$を使います。
\begin{aligned}
\mu_0 \frac{\partial H_x}{\partial t} &= 0 &…(9)\\
\mu_0 \frac{\partial H_y}{\partial t} &= \frac{\partial E_z}{\partial x} &…(10)\\
-\mu_0 \frac{\partial H_z}{\partial t} &= \frac{\partial E_y}{\partial x} &…(11)\\
\epsilon_0 \frac{\partial E_x}{\partial t} &= 0 &…(12)\\
-\epsilon_0 \frac{\partial E_y}{\partial t} &= \frac{\partial H_z}{\partial x} &…(13)\\
\epsilon_0 \frac{\partial E_z}{\partial t} &= \frac{\partial H_y}{\partial x} &…(14)
\end{aligned}
一次元で考えれば、少し見通しが良くなります。
電磁界の$x$成分は時間領域での変化は生じないことがわかりましたので、初期条件が保存され続けるようです。
また、電界の時間差分は磁界の空間差分で、磁界の時間差分は空間の時間差分で、それぞれ相互に表されており、時間変化に対して電界→磁界→電界→磁界→…と解いていけばうまいこと行きそうな気がしてきました。
離散化
では前章で求めた$(3) - (8)$式をFDTD法で用いるため、時間に対する漸化式として変形しましょう。
有限な時間差分長を$\Delta t$として、また有限な空間差分長を$\Delta x$表現すると、偏微分は以下のように展開できます。
\begin{aligned}
\frac{\partial A}{\partial t} & \approx \frac{A(t) - A(t-\Delta t)}{\Delta t}\\
\frac{\partial A}{\partial x} & \approx \frac{A(x) - A(x-\Delta x)}{\Delta x}
\end{aligned}
上記の展開を用いて、$(9) - (14)$式を展開しましょう。
ここで求めたいのは、ある時間$t$における、電磁界の$x$成分です。
- 式$(4)$について
\begin{aligned}
(左辺) &= \mu_0 \frac{\partial H_y}{\partial t} \approx \mu_0 \frac{H_y(t) - H_y(t-1)}{\Delta t}\\
(右辺) & = \frac{\partial E_z}{\partial x} \approx \frac{E_z(x) - E_z(x-1)}{\Delta x}
\end{aligned}
より、
\begin{aligned}
H_y(t) = H_y(t-\Delta t) &+ \frac{\Delta t}{\mu_0\Delta x}\Bigl(E_z(x) - E_z(x-\Delta x)\Bigr)
\end{aligned}
- 式$(7)$について
\begin{aligned}
(左辺) & = -\epsilon_0 \frac{\partial E_y}{\partial t} \approx -\epsilon_0 \frac{E_y(t) - E_y(t-\Delta t)}{\Delta t}\\
(右辺) &= \frac{\partial H_z}{\partial x} \approx \frac{H_z(x) - H_z(x-\Delta x)}{\Delta x}
\end{aligned}
より、
\begin{aligned}
E_y(t) = E_y(t-\Delta t) &+ \frac{\Delta t}{\epsilon_0\Delta x}\Bigl(H_z(x) - H_z(x-\Delta x)\Bigr)
\end{aligned}
式$(5), (8)$についても同様なので省略して、以下6つの時間発展に関する式$(15) - (20)$を得ます。
\begin{aligned}
H_x(t) &= H_x(0) &…(15)\\
H_y(t) &= H_y(t-\Delta t) + \frac{\Delta t}{\mu_0\Delta x}\Bigl(E_z(x) - E_z(x-\Delta x)\Bigr) &…(16)\\
H_z(t) &= H_z(t-\Delta t) - \frac{\Delta t}{\mu_0\Delta x}\Bigl(E_y(x) - E_y(x-\Delta x)\Bigr) &…(17)\\
E_x(t) &= E_x(0) &…(18)\\
E_y(t) &= E_y(t-\Delta t) - \frac{\Delta t}{\epsilon_0\Delta x}\Bigl(H_z(x) - H_z(x-\Delta x)\Bigr) &…(19)\\
E_z(t) &= E_z(t-\Delta t) + \frac{\Delta t}{\epsilon_0\Delta x}\Bigl(H_y(x) - H_y(x-\Delta x)\Bigr) &…(20)\\
\end{aligned}
お疲れさまでした。
以上で、任意の時間 $t$ における電磁界を、一つ前の時間ステップの電磁界と、電磁界の空間差分で表すことができました。
実践
辛い立式が終わったところで、実際に上式$(15) - (20)$を用いて、以下の条件で簡単な数値計算をはじめましょう。ようやく楽しくなります。
自由空間
でははじめに、以下図に示すような、自由空間（何もない空間）における電界の伝搬を可視化してみましょう。
$x, y, z$方向に無限に広がっているイメージです。
以下に使用するコードを示しました。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from tqdm import tqdm
# 定数の初期化
dx = 0.002 # 空間差分間隔[m]
c = 299792458.0 # 光速[m/s]
dt = dx/c * 0.99 # 時間差分間隔[s]
f = 3.5e9 # 周波数[Hz]
nx = 10000 # 計算点数
nt = 500 # 計算ステップ数
# 変数初期化
t = 0.0
eps = np.full(nx, 8.854187817e-12)
mu = np.full(nx, 1.2566370614e-6)
E_y = np.zeros(nx)
E_z = np.zeros(nx)
H_y = np.zeros(nx)
H_z = np.zeros(nx)
fig = plt.figure()
image_list = []
for _ in tqdm(range(nt)):
# 電界のy成分を励振
E_y[nx//2] = np.sin(2.0 * np.pi * f * t)
t += dt/2
E_y += -dt/(eps*dx) * (H_z - np.roll(H_z, shift=1)) # 式(16)
E_z += dt/(eps*dx) * (H_y - np.roll(H_y, shift=1)) # 式(17)
# 電界のy成分を励振
E_y[nx//2] = np.sin(2.0 * np.pi * f * t)
t += dt/2
H_y += -dt/(mu*dx) * (E_z - np.roll(E_z, shift=-1)) # 式(19)
H_z += dt/(mu*dx) * (E_y - np.roll(E_y, shift=-1)) # 式(20)
img = plt.plot(E_y[nx//2 - 200 : nx//2 + 200], color="blue")
plt.ylim(-3.0, 3.0)
image_list.append(img)
ani = animation.ArtistAnimation(fig, image_list, interval=50)
ani.save("freeSpace.gif", writer="pillow")
電界・磁界は交互に求める必要があるため、時間ステップを$1/2$にして与えています。
上記のコードで模擬した実行結果は以下図です。
図の中心（$x$=200）において励振された$\sin$波が隣接するセルからセルへどんどん伝搬していく様子が確認できるかと思います。
また、時間的波長（周期）と、空間的波長は、与えた周波数によって決定されており、いずれもただしい値になっているはずです（未確認）。
誘電体
では次に、以下図に示すような、誘電体を配置した空間における電界の伝搬を可視化してみましょう。
$y, z$方向に無限に大きく、$x$方向に半無限に大きい誘電体があるイメージです。（$x$=250のあたりから右方が誘電体）
以下に使用するコードを示しました。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from tqdm import tqdm
# 定数の初期化
dx = 0.002 # 空間差分間隔[m]
c = 299792458.0 # 光速[m/s]
dt = dx/c * 0.99 # 時間差分間隔[s]
f = 3.5e9 # 周波数[Hz]
nx = 10000 # 計算点数
nt = 500 # 計算ステップ数
# 変数初期化
t = 0.0
eps = np.full(nx, 8.854187817e-12)
mu = np.full(nx, 1.2566370614e-6)
eps[nx//2 + 50] = 100 * 8.854187817e-12 # 比誘電率100の誘電体を設置
E_y = np.zeros(nx)
E_z = np.zeros(nx)
H_y = np.zeros(nx)
H_z = np.zeros(nx)
fig = plt.figure()
image_list = []
for _ in tqdm(range(nt)):
# 電界のy成分を励振
E_y[nx//2] = np.sin(2.0 * np.pi * f * t)
t += dt/2
E_y += -dt/(eps*dx) * (H_z - np.roll(H_z, shift=1)) # 式(16)
E_z += dt/(eps*dx) * (H_y - np.roll(H_y, shift=1)) # 式(17)
# 電界のy成分を励振
E_y[nx//2] = np.sin(2.0 * np.pi * f * t)
t += dt/2
H_y += -dt/(mu*dx) * (E_z - np.roll(E_z, shift=-1)) # 式(19)
H_z += dt/(mu*dx) * (E_y - np.roll(E_y, shift=-1)) # 式(20)
img = plt.plot(E_y[nx//2 - 200 : nx//2 + 200], color="blue")
plt.ylim(-3.0, 3.0)
image_list.append(img)
ani = animation.ArtistAnimation(fig, image_list, interval=50)
ani.save("freeSpace.gif", writer="pillow")
実は、以下のたった一行を追加するだけで、誘電体を配置することができます。
eps[nx//2 + 50] = 100 * 8.854187817e-12 # 比誘電率100の誘電体を設置
上記のコードで模擬した実行結果は以下図です。誘電体を配置した位置近辺で、電界がなにやら怪しい動きをしています。また、誘電体内部では電界の振れ幅が何となく小さく見えますね。
誘電体が何かはよくわかりませんが、多分そういうものなんでしょう。あとで本で勉強します。
まとめ（ポエム）
概して、個人的には数値計算に以下の楽しさがあると思っています。
・よくわからない理論を可視化して理解できる楽しさ
・計算速度を追求するアスリート的なアツさ
・少し手法の改善をするだけで論文が出せる容易さ（これは怒られるので嘘）
中でも、数値計算の楽しいところは、【理論の理解を後回しにして、現象論的な理解をすすめることができるところ】です。
最初の立式は面倒な部分が多いですが、最初は兎にも角にもインターネットからコードを拾ってきて、実行してしまえば取りあえずは現象を”見る”ことができます。
百聞は一見に如かずとは言いますが、実際に目で見て得られる情報量は、理解の足がかりとして非常に大きなものですので、ぜひご活用いただければと思います。
次回は、今回のFDTD法のプログラムを3次元に拡張して、大規模な数値計算にチャレンジします。
蛇足
金属に入射する電磁界の動き
せっかくなので、以下図のように金属体を配置した場合についても計算してみましょう。
さっきの誘電体をそのまま金属に置き換えた場合です。（$x$=250のあたりから右方が金属）
以下に使用するコードを示します（手抜きではないです）。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from tqdm import tqdm
# 定数の初期化
dx = 0.002 # 空間差分間隔[m]
c = 299792458.0 # 光速[m/s]
dt = dx/c * 0.99 # 時間差分間隔[s]
f = 3.5e9 # 周波数[Hz]
nx = 10000 # 計算点数
nt = 500 # 計算ステップ数
# 変数初期化
t = 0.0
eps = np.full(nx, 8.854187817e-12)
mu = np.full(nx, 1.2566370614e-6)
E_y = np.zeros(nx)
E_z = np.zeros(nx)
H_y = np.zeros(nx)
H_z = np.zeros(nx)
fig = plt.figure()
image_list = []
for _ in tqdm(range(nt)):
# 電界のy成分を励振
E_y[nx//2] = np.sin(2.0 * np.pi * f * t)
t += dt/2
E_y += -dt/(eps*dx) * (H_z - np.roll(H_z, shift=1)) # 式(16)
E_z += dt/(eps*dx) * (H_y - np.roll(H_y, shift=1)) # 式(17)
E_y[nx//2 + 50] = 0.0 # 金属体に対する接線成分を0とする
E_z[nx//2 + 50] = 0.0 # 金属体に対する接線成分を0とする
# 電界のy成分を励振
E_y[nx//2] = np.sin(2.0 * np.pi * f * t)
t += dt/2
H_y += -dt/(mu*dx) * (E_z - np.roll(E_z, shift=-1)) # 式(19)
H_z += dt/(mu*dx) * (E_y - np.roll(E_y, shift=-1)) # 式(20)
img = plt.plot(E_y[nx//2 - 200 : nx//2 + 200], color="blue")
plt.ylim(-3.0, 3.0)
image_list.append(img)
ani = animation.ArtistAnimation(fig, image_list, interval=50)
ani.save("freeSpace.gif", writer="pillow")
計算結果は以下です。確かに金属を配置した位置より右方には電波が伝搬していませんね。（結果抜けてたので追記しました[2019/10/13]）
これも本当にかんたんなことです。
以下のたった二行を追加するだけで、金属を配置することができます。
E_y[nx//2 + 50] = 0.0 # 金属体に対する接線成分を0とする
E_z[nx//2 + 50] = 0.0 # 金属体に対する接線成分を0とする
なぜ、接線成分を消せば金属体を表現できるか、については電磁気学の教科書を参照ください。
数値計算・物理シミュレーションによって、現象を記述し、逆算的に理論の理解を深める助けとして活用できるのも、楽しさですね。
[次回]数値計算に親しむ ～FDTD法による電磁界解析編 (2)3次元への拡張とアンテナ放射解析～
https://qiita.com/sandshiP/items/27bc6ee6f079a88b0f08
参考文献
[1]宇野亨, 何一偉, 有馬卓司, "数値電磁界解析のためのFDTD法- 基礎と実践", コロナ社, 2016
本記事は、全編通して上記書籍を参考にしております。
理論的な部分から実装まで、詳細に記載されており、FDTD法の実装をする上では欠かせない名著です。
電磁界解析のみならず数値計算の入門者にとってもおすすめの一冊です。