FDTD(Finite Diffarence Time Domain)法のメモ書き

Qiitaアドカレ-数値解析の13日目です！初めてのアドカレ参加となります。

はじめに

以前の記事ではRCWA法の紹介を行っている。今回学ぶFDTD(Finite Difference Time Domain)法は古くYeeによって1966年に提唱された[Yee1966]。

この提唱当時では”障害物が波長に対して大きい場合…”との記述があることや、レーザー開発と同時期[maiman1962]であることからマイクロ波領域での利用が自然と考えられる。また、現代までに電磁場計算といえばこちらの計算手法がよくつかわれている。

また、PCの高性能化による計算パワーの増大などによって光領域においての利用も注目されており、メタサーフェス素子での開発でもたびたび利用されている。

このように非常に多岐にわたる応用があるため、市川さんによる解説記事[市川2006]、柴山さんによる解説記事[柴山2017]、Qiitaの平易な解説記事[Qiita link]、Pythonを使ったFDTDの教科書[岡本2019]、計算の具体例のスライド[柴山さんスライド]を代表とする日本語での解説記事がある。

このなかで本記事では、原論文となっているYeeによる論文[Yee1966]をベースにFDTDの中でも特に著名なTafloveによる教科書[taflove1995]を参考資料として、自分の理解を加えていく形でFDTDの紹介を行う。

多数の解説記事がある中で、本記事においては自分の理解とリファレンスノートを目的としている。

maxell方程式自分の理解

等方性の媒質においてMaxell方程式の微分形は以下のように表せる。(上からファラデーの法則、アンペールの法則、電場のガウスの法則、磁場のガウスの法則となっている。)

\begin{align}
\frac{\partial \boldsymbol{B}}{\partial t} + \nabla \times \boldsymbol{E} &= 0\\
\frac{\partial \boldsymbol{D}}{\partial t} – \nabla \times \boldsymbol{H} &= \boldsymbol{J}\\
\boldsymbol{B} &= \mu \boldsymbol{H}\\
\boldsymbol{D} &= \epsilon \boldsymbol{E}
\end{align}

ここで\( J, \mu, \epsilon\)はそれぞれ電流密度、透磁率、誘電率で空間と時間による関数である。直交座標系においてスカラーによる方程式に書き換えることができて、

\begin{align}
\frac{\partial B_z}{\partial t} &= \frac{\partial E_x}{\partial y} – \frac{\partial E_y}{\partial z}, \tag{1a}\\
\frac{\partial B_y}{\partial t} &= \frac{\partial E_z}{\partial x} – \frac{\partial E_x}{\partial z}, \tag{1b}\\
\frac{\partial B_x}{\partial t} &= \frac{\partial E_y}{\partial x} – \frac{\partial E_z}{\partial y}, \tag{1c}
\end{align}

\begin{align}
\frac{\partial D_x}{\partial t} &= \frac{\partial H_z}{\partial y} – \frac{\partial H_y}{\partial z} – J_x, \tag{1d}\\
\frac{\partial D_y}{\partial t} &= \frac{\partial H_x}{\partial z} – \frac{\partial H_z}{\partial x} – J_y, \tag{1e}\\
\frac{\partial D_z}{\partial t} &= \frac{\partial H_y}{\partial x} – \frac{\partial H_x}{\partial y} – J_z. \tag{1f}
\end{align}

ここで任意のベクトル場の成分は\( \boldsymbol{A} = (A_x, A_y, A_z) \)と表している。

Yeeアルゴリズム

Yeeアルゴリズムのメリット

図1. (a)Yeeアルゴリズムによる電磁場配置 (b)電磁場の時間配置

YeeアルゴリズムこそFDTDの中心である[Yee1966]。ここからは先述したmaxell方程式をもとに議論をしていくが、その前に概要とYeeアルゴリズムを採用するメリットを教科書[taflove1995]のp59-61ベースで述べていく。

1. Yeeアルゴリズムは電場と磁場について時間/空間依存で解いている点。
   – モーメント法による複合場積分法(combined-field integral equation)と同様に電場と磁場に境界条件を課すことで数値計算安定性が増す。
   – 電場と磁場の両方の情報を使うことで、片方のみの波動方程式を解く場合よりもロバストで、より広い種類の電磁構造を正確に扱える。
   – エッジや角の近くの磁場特異点のtangential成分、細いワイヤ近くの磁場特異点のazimuthal成分(azimuthには円周成分)などの各種特異点についても電場と磁場の計算をすることで、再現することが可能になる。
   
2. 図1(a)が示すように、3次元空間ではある電場成分の周りに4つの磁場成分またはある磁場成分の周りに4つの電場成分が存在する点。
   – これはとても美しく図となっており、ファラデーの法則とアンペールの法則を関連付けている。例えば、磁場ループに対応して電場の循環(電束の時間変化)になり、電場ループに関連している磁場の循環(磁束の時間変化)を捉えることができる。これはYeeアルゴリズムが微分形と積分系のmaxwell方程式を満たすことを意味している。後者は場の境界条件や特異点を指定するのに極めて有用である。
   – Yee格子での空間導関数は中心差分で離散化されるため、回転（curl）演算子は二次精度の中心差分として実装される。
   – 異種材料の界面がある基準軸に平行な場合、tangential方向の電場と磁場の連続性は自動的に満たされる。これは電磁場計算をする際に、特別な境界条件を課す必要が生まれない。実際には、各格子点に透磁率と誘電率を設定するだけでいい。長方形のYee格子にむけて、これはYee格子のセルサイズによる空間分解能に従って表面と構造の内部座標が階段近似されることによってもたらされる。
   – Yee格子内や中心差分作用による電場と磁場の場所は暗黙的に2つのガウスの法則を満たす。すなわちYee格子は初期条件では電場と磁場は勾配(divergence)フリーで、それによってモデル化される磁場/電場源がない空間では自由電荷や自由磁荷が発生しない。(補足：これは通常の有限差分法ではdivergence – free(div D = 0)を常に満たすわけではないので、数値誤差が蓄積することで偽の電荷(div D = ρ:偽の電荷)や磁荷が発生する)
   
3. 図1(b)に示すように、Yeeアルゴリズムは電場と磁場を時間方向に半ステップ(leapfrog)ずらして更新する。計算対象の3次元空間の電場(磁場)をぜっぶ計算し、それを使って次の磁場を計算している。このサイクルは交互に繰り返され、時間のステップが完了するまで繰り返される。
   – この半ステップずらし更新法(leapfrog)は完全に明示的(explicit)で、連立方程式や逆行列に関連する問題を回避できる。
   – 空間微分（curl 演算子に含まれる部分）は中心差分で離散化され、二次精度を持つ。
   – 時間ステップのアルゴリズムは非散逸(non dissipative)である；メッシュ内を伝搬する数値波が、time-steppingアルゴリズムの数値計算特有の減衰(物理に即していない減衰)が生じない。

Yeeアルゴリズムの数式的理解

Yeeアルゴリズムを考えていく、格子インデックスが\( (i,j,k) \)のとき、単位距離\( (\Delta x, \Delta y, \Delta z)\)を用いて実際の格子座標は

\begin{equation}
(i\Delta x, j\Delta y, k\Delta z)
\end{equation}

と表すことができる。また、空間と時間の任意の関数について次のように置く、

\begin{equation}
F(i\Delta x, j\Delta y, k\Delta z, n\Delta t) = F^n(i,j,k)
\end{equation}

ここで境界条件を完全導体としてスカラーで表されるMaxwell方程式(式\( 1a\))を中心差分近似(参考:鳥取大学講義資料)を用いて差分方程式に書き直すと

\begin{align}
&\frac{
B_x^{\,n+\frac{1}{2}} \left(i, j+\frac{1}{2}, k+\frac{1}{2}\right)
– B_x^{\,n-\frac{1}{2}} \left(i, j+\frac{1}{2}, k+\frac{1}{2}\right)
}{
\Delta t}\\
&= \frac{
E_y^{\,n} \left(i, j+\frac{1}{2}, k+1\right)
– E_y^{\,n} \left(i, j+\frac{1}{2}, k\right)
}{
\Delta z
}
–
\frac{
E_z^{\,n} \left(i, j+1, k+\frac{1}{2}\right)
-E_z^{\,n} \left(i, j, k+\frac{1}{2}\right)
}{
\Delta y
}.\tag{2a}
\end{align}

また、同様に磁場(式\( 1b\))に着目して式変形すると

\begin{align}
&\frac{
D_x^{n} \left(i+\frac{1}{2}, j, k\right)
– D_x^{\ n-\frac{1}{2}} \left(i+\frac{1}{2}, j, k \right)
}{
\Delta t}\\
&= \frac{
H_z^{n-\frac{1}{2}} \left(i+\frac{1}{2}, j+\frac{1}{2}, k\right)
– H_z^{n-\frac{1}{2}} \left(i+\frac{1}{2}, j-\frac{1}{2}, k\right)
}{
\Delta y
}
–
\frac{
H_y^{n-\frac{1}{2}} \left(i+\frac{1}{2}, j, k+\frac{1}{2}\right)
-H_y^{n-\frac{1}{2}} \left(i+\frac{1}{2}, j, k-\frac{1}{2}\right)
}{
\Delta z
}\\
&+ J_x^{n-\frac{1}{2}}\left(i+\frac{1}{2},j,k \right).\tag{2b}
\end{align}

となる。これはほかのMaxell方程式でも同様に成り立つ。これを視覚化したのが図1となる。まさにこれがYeeアルゴリズムと呼ばれる部分であり、FDTDの肝となる。

境界条件

境界条件についてはFDTDの元論文[Yee1966]では完全導体としている。完全導体の条件下では境界部分では接線成分の電場と法線成分の磁場が0になる。したがって、導体表面は座標軸に平行な立方体の集合として近似される。これらは閉じた空間に適用するのに向いており、完全導体境界条件は構造体の共振モードや導波路モードのような閉じた系 (ケージ、箱、導波路、金属壁キャビティ など) の解析に向いていると考えられる。

一方で、昨今のFDTDでは吸収境界条件(ABCs: Absorbing Boundary Conditions)が用いられることが多い。FDTDの教科書[taflove1995]においてもABCsでのChapter(p145, chapter.7)があるくらいには詳細に書かれている。ABCsでは無限の空間を想定したもので、電磁場が反射せずに逃げていくイメージとなる[加古2005]。ここでABCsでは具体的にどういった境界条件となるかは様々な種類があり、Perfectly Matched Layer(PML)やHuygens ABCsをはじめとしていくつか存在する。本記事では詳細に記すことはないが、ABCsの歴史をreviewした論文[Jean2015]などを参考にするとよさそうである。このABCsでは構造 (回折格子、フォトニック結晶、メタレンズなど) を含んだ開放領域のような解析に向いていると考えられる。

こちらの解説については本記事で書きれないので、上記に記した参考文献を参考にされたい。

格子サイズと計算の安定性(stabirity criterion)

論文には文間が多かったので、教科書[taflove1995]を参考にする。

\(\mu = 1, \epsilon = 1, \sigma = 0, \rho’=0, c=1 \)とするときに複素数\(j=\sqrt{-1}\)としてMaxell方程式を書き直すと

\begin{equation}
j \nabla \times ( \boldsymbol{H} + j \boldsymbol{E}) = \frac{\partial}{\partial t} ( \boldsymbol{H} + j \boldsymbol{E})
\end{equation}

これをもっとシンプルに表記すると

\begin{equation}
j \nabla \times \boldsymbol{V} = \frac{\partial \boldsymbol{V}}{\partial t}
\end{equation}

ここで\( \boldsymbol{H} + j \boldsymbol{E} = \boldsymbol{V} \)である。

数値解析からの要請

数値解析における安定性のためには次の固有値問題を考慮することでシンプルにすることができ、

\begin{align}
\left. \frac{\partial}{\partial t} \right|_{numerical} \boldsymbol{V} &= \Lambda \boldsymbol{V} \tag{3a} \\
\left. \nabla\right|_{numerical} \times \boldsymbol{V} &= \Lambda \boldsymbol{V} \tag{3b} \\
\end{align}

Yeeアルゴリズムを数値解析での時間微分関数として使うと、式(3a)を以下のように書き換えることができ

\begin{equation}
\frac{V^{n+1/2} – V^{n-1/2}}{\Delta t} = \Lambda V^{n} \tag{4}
\end{equation}

ここで、空間上ではある一点を考慮している(\(i,j,k\)は一定)。不安定因子(growth factor)を次のように定義する。

\begin{equation}
q = \frac{V^{n+1/2}}{V^{n}} = \frac{V^{n}}{V^{n-1/2}}
\end{equation}

これはすべての\( n\)で成り立つ。ここでこの不安定因子\(q\)について、\( |q| \leq 1\)のときに、すべての空間座標系において計算が収束する。数式(4)を不安定因子を用いて示すと

\begin{equation}
\frac{q V^{n} – V^{n}/q}{\Delta t} = \Lambda V^{n}
\end{equation}

共通因数\( V^n\)をくくりだすと

\begin{align}
\frac{q ^2 – 1}{q \Delta t} = \Lambda \\
q ^2 + ( \Lambda \Delta t )q -1 = 0
\end{align}

二次方程式の解の公式から

\begin{align}
q &= \frac{\Lambda \Delta t}{2} \pm \sqrt{\left( \frac{\Lambda \Delta t}{2} \right)^2 + 1}\\
& \equiv a \pm \sqrt{a^2+1}
\end{align}

ここで\( a = \frac{\Lambda \Delta t}{2} \)とした。ここで厳しい条件の\( q = 1\) (物理的には損失なしのモデル)として方程式を解くと\( a\)は純虚数となる。これを示していく。まずは\( a\)を純虚数と置くと

\begin{equation}
a = \rm{Re}(a) + j\rm{Im}(a)
\end{equation}

純虚数なので、\(\rm{Re}(a) = 0\)。また、\( -1 \leq \rm{Im}(a) \leq 1\)すると

\begin{equation}
\sqrt{a^2+1} = \sqrt{[j\rm{Im}(a)]^2 + 1} = \sqrt{1 – [\rm{Im}(a)]^2}:\rm{実数}
\end{equation}

ここで\(q\)を書き直すと

\begin{equation}
q = j\rm{Im}(a) \pm \sqrt{a^2+1}
\end{equation}

第一項は純虚数、第二項は実数となるので

\begin{align}
|q| &= \sqrt{ [\rm{Im}(a)]^2 + \sqrt{1 – [\rm{Im}(a)]^2}^2}\\
&= \sqrt{ 1 – [\rm{Im}(a)]^2 + [\rm{Im}(a)]^2} = 1　\tag{5}
\end{align}

となるので\( q = 1\)が成り立つことがわかる。

さて、\( a = \frac{\Lambda \Delta t}{2} \)と置いたのと\(-1 \leq \rm{Im}(a) \leq 1\)を利用すると

\begin{equation}
-1 \leq \frac{\Delta t}{2} \rm{Im}(\Lambda) \leq 1
\end{equation}

ここで\(\Lambda\)は純虚数(\(\Lambda = j\rm{Im}(\Lambda)\))を利用した。これを整理して

\begin{equation}
-\frac{2}{\Delta t} \leq \rm{Im}(\Lambda) \leq \frac{2}{\Delta t} \tag{6}
\end{equation}

Yeeアルゴリズムからの要請

式(3b)の

\begin{equation}
\left. \nabla\right|_{numerical} \times \boldsymbol{V}= \Lambda \boldsymbol{V}
\end{equation}

を考えていく。\(\boldsymbol{V}\)は以下のように置き換えることができ

\begin{equation}
\boldsymbol{V} = \boldsymbol{V_0} e^{j(Ik_x\Delta x + Jk_y\Delta y + K k_z\Delta z)}
\end{equation}

となる。ここで\(I,J,K\)は格子空間上での格子位置、\(k_x,k_y,k_z\)は波数ベクトルの表記は省略している。

続いて、式(2a)-式(2b)で表すようにcurl演算子を作用させて差分化する。curl の \(x\)成分は

\begin{equation} (\nabla \times \boldsymbol{V})_x = \frac{\partial V_z}{\partial y} – \frac{\partial V_y}{\partial z} \end{equation}

であるため、\(\Lambda \boldsymbol{V}\) の \(x\) 成分は (\(V_y, V_z\) ) を考える。 まず \( \partial V_z/\partial y \) の差分をとると

\begin{equation} \frac{\partial V_z}{\partial y} \approx \frac{e^{jk_y(J+1/2)\Delta y} – e^{jk_y(J-1/2)\Delta y}}{\Delta y} (\boldsymbol{V_0})_z e^{jk_x I\Delta x + jk_z K\Delta z} \end{equation}

となり，\(e^{j\theta}-e^{-j\theta} = j2\sin\theta\)を使って整理すると

\begin{equation} (\Lambda \boldsymbol{V})_{x,1} = 2 e^{j(k_xI\Delta x + k_yJ\Delta y + k_zK\Delta z)} \frac{(\boldsymbol{V_0})_z}{\Delta y} j\sin(k_y\Delta y/2) \end{equation}

同様に，\(-\partial V_y/\partial z\) は

\begin{equation} (\Lambda \boldsymbol{V})_{x,2} = -2e^{j(k_xI\Delta x + k_yJ\Delta y + k_zK\Delta z)} \frac{(\boldsymbol{V_0})_y}{\Delta z} j\sin(k_z\Delta z/2) \end{equation}

である。よって curl の (x) 成分は

\begin{equation} (\Lambda\boldsymbol{V})_x = -2\left[ \frac{\tilde{z}}{\Delta y}\sin(k_y\Delta y/2) – \frac{\tilde{y}}{\Delta z}\sin(k_z\Delta z/2) \right] \end{equation}

ここで

\begin{equation} \tilde{y}= (\boldsymbol{V_0})_y e^{j(k_xI\Delta x + k_yJ\Delta y + k_zK\Delta z)},\quad \tilde{z}= (\boldsymbol{V_0})_z e^{j(k_xI\Delta x + k_yJ\Delta y + k_zK\Delta z)} \end{equation}

と置いた。また、これらのチルダ成分は単位ベクトル\( (\boldsymbol{V_0})_y,(\boldsymbol{V_0})_z \)を含んでいることに注意したい。同様に\((y,z)\) 成分も

\begin{align} (\Lambda\boldsymbol{V})_y &= -2\left[ \frac{\tilde{x}}{\Delta z}\sin(k_z\Delta z/2) – \frac{\tilde{z}}{\Delta x}\sin(k_x\Delta x/2) \right]\\ (\Lambda\boldsymbol{V})_z &= -2\left[ \frac{\tilde{y}}{\Delta x}\sin(k_x\Delta x/2) – \frac{\tilde{x}}{\Delta y}\sin(k_y\Delta y/2) \right] \end{align}

よって式(3b)のcurl演算子を置き換えて

\begin{equation}
-2\left[ \frac{\tilde{x}}{\Delta x } \sin (k_x\Delta x/2) + \frac{\tilde{y}}{\Delta y } \sin (k_y\Delta y/2) + \frac{\tilde{z}}{\Delta z } \sin (k_z\Delta z/2) \right] \times \boldsymbol{V}= \Lambda \boldsymbol{V}
\end{equation}

となる。外積は歪対象行列であり、\( \boldsymbol{A} \times \boldsymbol{X} = a \boldsymbol{X}\)としたときこの固有値は\(0,\pm|A|\)ということが知られている(参考：3×3の歪対称行列の特異値分解)。そこで \(\Lambda^2\)について解いてみると

\begin{equation}\Lambda^2 = -4\left[ \frac{1}{(\Delta x)^2}\sin^2(\tilde{k}_x\Delta x/2) +\frac{1}{(\Delta y)^2}\sin^2(\tilde{k}_y\Delta y/2) +\frac{1}{(\Delta z)^2}\sin^2(\tilde{k}_z\Delta z/2) \right] \end{equation}

さらに歪対象行列の固有値は純虚数であるため

\begin{equation} \rm{Re}(\Lambda) = 0,\qquad \rm{Im}(\Lambda)=\pm 2 \sqrt{\frac{\sin^2(\tilde{k}_x\Delta x/2)}{(\Delta x)^2} +\frac{\sin^2(\tilde{k}_y\Delta y/2)}{(\Delta y)^2} +\frac{\sin^2(\tilde{k}_z\Delta z/2)}{(\Delta z)^2}} \end{equation}

また、\( \rm{Im}(\Lambda)\)の取りうる範囲は

\begin{equation} -2\sqrt{\frac{1}{(\Delta x)^2} +\frac{1}{(\Delta y)^2} +\frac{1}{(\Delta z)^2}} \leq \rm{Im}(\Lambda) \leq 2\sqrt{\frac{1}{(\Delta x)^2} +\frac{1}{(\Delta y)^2} +\frac{1}{(\Delta z)^2}}\tag{7} \end{equation}

となる。

数値安定性の条件

ここまでで導出してきた式をまとめていく。より計算を安定させるためには、式(7)は式(6)の条件を満たす必要がある。(空間固有値スペクトル ⊂ 時間固有値許容領域である。数値計算を安定させたうえ(式(6))でFDTD計算(式(7))を行うから)

\begin{equation} 2\sqrt{ \frac{1}{(\Delta x)^2}+\frac{1}{(\Delta y)^2}+\frac{1}{(\Delta z)^2} } \leq \frac{2}{\Delta t}\end{equation}

これを整理して

\begin{equation} \Delta t \leq \frac{1}{2\sqrt{ \frac{1}{(\Delta x)^2}+\frac{1}{(\Delta y)^2}+\frac{1}{(\Delta z)^2} }}\end{equation}

となる。これがCourant条件[岡本2019]と呼ばれる。

さいごに

本記事ではFDTDの基礎的な部分(論文の前半部分)についてをまとめた。具体的な計算手法(論文の後半)については今後の記事で時間があれば他の記事で書いていきたい。

基本の計算に関しては非常にシンプルでありつつも、そのメリットが大きいことを示せたと考えている。

また、次の記事の予定ではRCWAの具体的な計算部分をまとめていきたいと思っている。