今回から電磁界シミュレーションの作成に入ります。
前回の記事はこちら↓
電磁界解析の手法には大きく分けて FDTD法 と 有限要素法 の二つがあります。このブログでは、FDTD法の解析方法を使用します。
FDTD法とは
FDTD法とは、Finite Difference Time Domainの略で、「空間と時間を離散化してそれぞれの領域においてマクスウェルの方程式を解く」解析手法のことを言います。
具体的には下記のマクスウェルの方程式を空間と時間で離散化して解いていきます。
\begin{align}
\boldsymbol{\nabla} \times \boldsymbol{E} &= -\frac{\partial \boldsymbol{B}}{\partial t}\tag{1}\\
\boldsymbol{\nabla} \times \boldsymbol{H} &= \frac{\partial \boldsymbol{D}}{\partial t}+\boldsymbol{J} \tag{2}
\end{align}
時間・空間における離散化の考え方
電界$\boldsymbol{E}$、磁界$\boldsymbol{H}$は、位置と時間の変数を持つため, 3次元空間を対象とした場合、$\boldsymbol{E}$, $\boldsymbol{H}$はそれぞれ$\boldsymbol{E}(x,y,z,t)$, $\boldsymbol{H}(x,y,z,t)$とかけます。
時間の離散化
電界と磁界それぞれの時間項は$\Delta t$の離散ステップを持つものとし、一般項nを用いて$n\Delta t$とかきます。
ここで、式更新時にメモリの節約をするために、電界と磁界は$\frac{1}{2}$ステップ時間をずらして表記することとします。イメージとしては下図のとおりです。どうしてメモリが節約できるかは後ほど解説します。
なお、今後の表記は時間項を変数の右上に書くようにします。具体的には下図の通りです。
\begin{align}
\boldsymbol{E}(x,y,z,n\Delta t)&=\boldsymbol{E}^n (x,y,z)\\
\boldsymbol{H}(x,y,z,(n+\frac{1}{2}) \Delta t)&=\boldsymbol{H}^{n+\frac{1}{2}} (x,y,z)
\end{align}
空間の離散化
次に、空間の離散化について説明します。
こちらも少し特殊な工夫が入ります。マクスウェルの方程式を解くにあたって、解析空間を$\Delta d$の立方体で分割していきます。一つの立方体（2次元であれば正方形）を解析空間におけるユニットセルと呼び、セルにおける電界と磁界の配置を下図に示します。
$i$,$j$,$k$は原点からどの立方体の頂点を指しているかを示しています。原点からx方向に$i$、y方向に$j$、z方向に$k$離れていたとき、座標は$(i\Delta d, j\Delta d, k\Delta d)$で表しますが、$\Delta d$を何度も書くのは大変なため、これを$(i,j,k)$のように書くこととします。
電界と磁界成分の配置がしれっと書かれていますが、ここに工夫がなされています。これはYeeセルと呼ばれる電磁界のセル配置方法で、電界と磁界が各方向に半セルずれて配置されています。こうすることで回転の計算がユニットセル内で簡潔に計算できます。
時間の離散化
時間、空間の離散化について概要を示しましたので、具体的にマクスウェルの方程式を離散化していきましょう。まずは時間の離散化からです。
式(1)、(2)を離散化してみます。まずは式(1)からです。磁束密度$\boldsymbol{B}$は透磁率と磁界の積ですので、式変形をすると
\begin{align}
\frac{\partial \boldsymbol{H}}{\partial t}=-\frac{1}{\mu}\boldsymbol{\nabla} \times \boldsymbol{E}\tag{3}\\
\end{align}
次に時間を離散化していきます。式(3)で時間の項は左辺の偏微分のみです。
偏微分を微小時間$\Delta t$の差分に変換します。ここで、左辺の差分は$n+\frac{1}{2}$と$n-\frac{1}{2}$とすると、右辺の時間はその間である$n$となるので、
\begin{align}
\frac{\boldsymbol{H}^{n+\frac{1}{2}}-\boldsymbol{H}^{n-\frac{1}{2}}}{\Delta t}=-\frac{1}{\mu}\boldsymbol{\nabla} \times \boldsymbol{E}^n\tag{4}\\
\end{align}
\begin{align}
\boldsymbol{H}^{n+\frac{1}{2}}=\boldsymbol{H}^{n-\frac{1}{2}}-\frac{\Delta t}{\mu}\boldsymbol{\nabla} \times \boldsymbol{E}^n\tag{5}\\
\end{align}
これで、式(1)の時間の離散化は完了です。
続いて式(2)について時間の離散化をします。電束密度$D$は誘電率$\epsilon$と電界$E$の積で表せるため、式変形すると
\begin{align}
\frac{\partial \boldsymbol{E}}{\partial t}=\frac{1}{\varepsilon}\boldsymbol{\nabla} \times \boldsymbol{H}-\frac{1}{\varepsilon}\boldsymbol{J} \tag{6}\\
\end{align}
ここで、左辺の差分は$n+1$と$n$とすると、右辺の時間はその間である$n+\frac{1}{2}$となり、$J$は印加電流密度$\boldsymbol{J_{0}}$とオームの法則からなる$\sigma \boldsymbol{E}$からなるので
\begin{align}
\frac{\boldsymbol{E}^{n+1}-\boldsymbol{E}^n}{\partial t}=\frac{1}{\varepsilon}\boldsymbol{\nabla} \times \boldsymbol{H}^{n+\frac{1}{2}}-\frac{1}{\varepsilon}\boldsymbol{J_{0}}-\frac{1}{\varepsilon}\sigma \boldsymbol{E}^{n+\frac{1}{2}}\tag{7}
\end{align}
ここで、$\boldsymbol{E}^{n+\frac{1}{2}}$は時間のルール上存在しないため差分で置換します
$$\boldsymbol{E}^{n+\frac{1}{2}}=\frac{\boldsymbol{E}^{n+1}-\boldsymbol{E}^{n}}{2}\tag{8}$$
式(7)に式(8)を入れて整理すると、
$$\boldsymbol{E}^{n}=\frac{1-\frac{\Delta t}{2 \varepsilon}}{1+\frac{\Delta t}{2 \varepsilon}}\boldsymbol{E}^{n-1}+\frac{\frac{\Delta t}{\varepsilon}}{1+\frac{\sigma \Delta t}{2\varepsilon}}\boldsymbol{\nabla} \times \boldsymbol{H}^{n-\frac{1}{2}}-\frac{\frac{\Delta t}{\varepsilon}}{1+\frac{\sigma \Delta t}{2\varepsilon}}\boldsymbol{J_{0}}\tag{9}$$
かなり複雑になってしまいましたが、これで式(2)の時間の離散化が完了しました。
空間の離散化
時間の離散化が終わりましたので、次は式(5)、(9)に対して、空間の離散化を行います。
磁界の離散化
式(5)について、$\boldsymbol{H}^{n+\frac{1}{2}}$をx, y, z成分に分解します。その際、右辺の回転の式についても各成分に分解し、ユニットセルの座標を用いて表記することとします。
\begin{align}
H_{x}^{n+\frac{1}{2}}\Biggl( i,j+\frac{1}{2},k+\frac{1}{2}\Biggr)=H_{x}^{n-\frac{1}{2}}\Biggl(i,j+\frac{1}{2},k+\frac{1}{2}\Biggr)-\frac{\Delta t}{\mu}\Biggl(\frac{\partial E_{z}^n\Bigl(i,j+\frac{1}{2},k+\frac{1}{2}\Bigr)}{\partial y}-\frac{\partial E_{y}^n\Bigl(i,j+\frac{1}{2},k+\frac{1}{2}\Bigr)}{\partial z}\Biggr) \tag{10} \\
H_{y}^{n+\frac{1}{2}}\Biggl( i+\frac{1}{2},j,k+\frac{1}{2}\Biggr)=H_{y}^{n-\frac{1}{2}}\Biggl(i+\frac{1}{2},j,k+\frac{1}{2}\Biggr)-\frac{\Delta t}{\mu}\Biggl(\frac{\partial E_{x}^n\Bigl(i+\frac{1}{2},j,k+\frac{1}{2}\Bigr)}{\partial z}-\frac{\partial E_{z}^n\Bigl(i+\frac{1}{2},j,k+\frac{1}{2}\Bigr)}{\partial x}\Biggr) \tag{11}\\
H_{z}^{n+\frac{1}{2}}\Biggl(i+\frac{1}{2},j+\frac{1}{2},k\Biggr)=H_{z}^{n-\frac{1}{2}}\Biggl(i+\frac{1}{2},j+\frac{1}{2},k\Biggr)-\frac{\Delta t}{\mu}\Biggl(\frac{\partial E_{y}^n\Bigl(i+\frac{1}{2},j+\frac{1}{2},k\Bigr)}{\partial x}-\frac{\partial E_{x}^n\Bigl(i+\frac{1}{2},j+\frac{1}{2},k\Bigr)}{\partial y}\Biggr) \tag{12}\\
\end{align}
右辺に偏微分の項があるので、これをユニットセル内における差分形式に変形します。
例えば、式(10)の偏微分に対して、差分形式に変換させると下式のようになります。
$$\frac{\partial E_{z}^n\Bigl(i,j+\frac{1}{2},k+\frac{1}{2}\Bigr)}{\partial y}=\frac{E_{z}^n\Bigl(i,j+1,k+\frac{1}{2}\Bigr)-E_{z}^n\Bigl(i,j,k+\frac{1}{2}\Bigr)}{\Delta d}$$
これを式(10),(11),(12)に適用すると、
\begin{align}
H_{x}^{n+\frac{1}{2}}\Biggl( i,j+\frac{1}{2},k+\frac{1}{2}\Biggr)=H_{x}^{n-\frac{1}{2}}\Biggl(i,j+\frac{1}{2},k+\frac{1}{2}\Biggr)-\frac{\Delta t}{\mu \Delta d}\Biggl(E_{z}^n\Bigl(i,j+1,k+\frac{1}{2}\Bigr)-E_{z}^n\Bigl(i,j,k+\frac{1}{2}\Bigr)-E_{y}^n\Bigl(i,j+\frac{1}{2},k+1\Bigr)+E_{y}^n\Bigl(i,j+\frac{1}{2},k\Bigr)\Biggr) \tag{13} \\
H_{y}^{n+\frac{1}{2}}\Biggl( i+\frac{1}{2},j,k+\frac{1}{2}\Biggr)=H_{y}^{n-\frac{1}{2}}\Biggl(i+\frac{1}{2},j,k+\frac{1}{2}\Biggr)-\frac{\Delta t}{\mu \Delta d}\Biggl(E_{x}^n\Bigl(i+\frac{1}{2},j,k+1\Bigr)-E_{x}^n\Bigl(i+\frac{1}{2},j,k\Bigr)-E_{z}^n\Bigl(i+1,j,k+\frac{1}{2}\Bigr)+E_{z}^n\Bigl(i,j,k+\frac{1}{2}\Bigr)\Biggr) \tag{14}\\
H_{z}^{n+\frac{1}{2}}\Biggl( i+\frac{1}{2},j+\frac{1}{2},k\Biggr)=H_{z}^{n-\frac{1}{2}}\Biggl(i+\frac{1}{2},j+\frac{1}{2},k\Biggr)-\frac{\Delta t}{\mu \Delta d}\Biggl(E_{y}^n\Bigl(i+1,j+\frac{1}{2},k\Bigr)-E_{y}^n\Bigl(i,j+\frac{1}{2},k\Bigr)-E_{x}^n\Bigl(i+\frac{1}{2},j+1,k\Bigr)+E_{x}^n\Bigl(i+\frac{1}{2},j,k\Bigr)\Biggr) \tag{15}\\
\end{align}
電界の離散化
式(9)に対しても同様にx,y,z成分に分解して式を立式します。係数が複雑なため、下記のように代数を定義します。
\begin{align}
\boldsymbol{E}^{n}&=A\boldsymbol{E}^{n-1}+B\boldsymbol{J_{0}}+C\boldsymbol{\nabla} \times \boldsymbol{H}^{n-\frac{1}{2}}\tag{16} \\
A&=\frac{1-\frac{\sigma \Delta t}{2\varepsilon}}{1+\frac{\sigma \Delta t}{2\varepsilon}} \\
B&=\frac{\frac{\Delta t}{\varepsilon}}{1+\frac{\sigma \Delta t}{2\varepsilon}} \\
C&=\frac{\frac{\Delta t}{\varepsilon}}{1+\frac{\sigma \Delta t}{2\varepsilon}}
\end{align}
磁界の離散化と同様に式(16)をx,y,z成分に分解すると
\begin{align}
E_{x}^{n+1}\Biggl( i+\frac{1}{2},j,k\Biggr)=AE_{x}^{n}\Biggl(i+\frac{1}{2},j,k\Biggr)-BJ_{0x}\Biggl( i+\frac{1}{2},j,k\Biggr)+\frac{C}{\Delta d}\Biggl(H_{z}^{n+\frac{1}{2}}\Bigl(i+\frac{1}{2},j+\frac{1}{2},k\Bigr)-H_{z}^{n+\frac{1}{2}}\Bigl(i+\frac{1}{2},j-\frac{1}{2},k\Bigr)-H_{y}^{n+\frac{1}{2}}\Bigl(i+\frac{1}{2},j,k+\frac{1}{2}\Bigr)+H_{y}^{n+\frac{1}{2}}\Bigl(i+\frac{1}{2},j,k+\frac{1}{2}\Bigr)\Biggr) \tag{17} \\
E_{y}^{n+1}\Biggl( i,j+\frac{1}{2},k\Biggr)=AE_{y}^{n}\Biggl(i,j+\frac{1}{2},k\Biggr)-BJ_{0y}\Biggl( i,j+\frac{1}{2},k\Biggr)+\frac{C }{\Delta d}\Biggl(H_{x}^{n+\frac{1}{2}}\Bigl(i,j+\frac{1}{2},k+\frac{1}{2}\Bigr)-H_{x}^{n+\frac{1}{2}}\Bigl(i,j+\frac{1}{2},k-\frac{1}{2}\Bigr)-H_{z}^{n+\frac{1}{2}}\Bigl(i+\frac{1}{2},j,k+\frac{1}{2}\Bigr)+H_{z}^{n+\frac{1}{2}}\Bigl(i-\frac{1}{2},j+\frac{1}{2},k\Bigr)\Biggr) \tag{18}\\
E_{z}^{n+1}\Biggl( i,j,k+\frac{1}{2}\Biggr)=AE_{z}^{n}\Biggl( i,j,k+\frac{1}{2}\Biggr)-BJ_{0z}\Biggl( i,j,k+\frac{1}{2}\Biggr)+\frac{C}{ \Delta d}\Biggl(H_{y}^{n+\frac{1}{2}}\Bigl(i+\frac{1}{2},j,k+\frac{1}{2}\Bigr)-H_{y}^{n+\frac{1}{2}}\Bigl(i+\frac{1}{2},j,k+\frac{1}{2}\Bigr)-H_{x}^{n+\frac{1}{2}}\Bigl(i,j+\frac{1}{2},k+\frac{1}{2}\Bigr)+H_{x}^{n+\frac{1}{2}}\Bigl(i,j-\frac{1}{2},k+\frac{1}{2}\Bigr)\Biggr) \tag{19}\\
\end{align}
のようになります（間違っていたらごめんなさい）
まとめ
今回はマクスウェルの方程式のうち、磁界と電界が伝播する方程式について空間と時間の離散化をしてみました。離散化するにあたりYeeセルや電界と磁界の周期を1/2だけズラすことによる計算の工夫が入り、複雑になってしまいますが、これでFDTD法によるシミュレーションをするための準備が整いました。
次回はいよいよシミュレーションを実践してみようと思います。
コメントを残す