Go to Qiita Advent Calendar 2024 Top
LoginSignup
8
12

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@sandshiP

数値計算に親しむ　～FDTD法による電磁界解析編 (2)3次元への拡張とアンテナ放射解析～

Last updated at Posted at 2019-10-13

はじめに

前編に引き続き、楽しい数値計算をします。

[前編]数値計算に親しむ　～FDTD法による電磁界解析編 (1)概要と基礎実践～
https://qiita.com/sandshiP/items/f07f6c1443e024cb775a

今回も、時間領域有限差分（Finite Difference Time Domain, FDTD）法を試します。
前回は、1次元FDTD法を使った立式・基礎実践まで行いましたので、今回からは3次元FDTD法を用いた実用的な電磁界解析に入っていきます。立式パートは前回の拡張なのでほとんど飛ばします。

立式（飛ばし読み可）

今回も、以下のMaxwell方程式$(1) - (2)$から、電磁界の各成分に分解した式$(3) - (8)$を離散化していきます。式$(2)$については、前回と異なり、電流要素$\boldsymbol{j}$が追加されています。理由は本章最後に記載します。

\begin{aligned}
\nabla \times \boldsymbol{E}(t, \boldsymbol{r}) + \mu \frac{\partial \boldsymbol{H}(t, \boldsymbol{r})}{\partial t} = 0 　&…(1)\\
\nabla \times \boldsymbol{H}(t, \boldsymbol{r}) - \epsilon \frac{\partial \boldsymbol{E}(t, \boldsymbol{r})}{\partial t} = \boldsymbol{j}(t, \boldsymbol{r}) 　&…(2)\\

\mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z} - \frac{\partial E_z}{\partial y}　&…(3)\\
\mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}　&…(4)\\
\mu \frac{\partial H_z}{\partial t} = \frac{\partial E_x}{\partial y} - \frac{\partial E_y}{\partial x}　&…(5)\\
\epsilon \frac{\partial E_x}{\partial t} + j_x = \frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z}　&…(6)\\
\epsilon \frac{\partial E_y}{\partial t} + j_y = \frac{\partial H_x}{\partial z} - \frac{\partial H_z}{\partial x}　&…(7)\\
\epsilon \frac{\partial E_z}{\partial t} + j_z = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}　&…(8)
\end{aligned}

以下の近似を用いてガシガシ展開し、最新の時刻における$A(t)$について解いていきます。

\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}

解説のため、磁界・電界の更新式としてそれぞれ、式$(3), (6)$についてのみ、過程を示します。

  • 式$(3)$について
\begin{aligned}
(左辺)=\mu\frac{\partial H_x}{\partial t}  &\approx \mu\frac{H_x(t) - H_x(t-\Delta t)}{\Delta t}\\
(右辺)=\frac{\partial E_y}{\partial z} - \frac{\partial E_z}{\partial y} &\approx \frac{1}{\Delta z}\Bigl(E_y(z) - E_y(z-\Delta z)\Bigr) \\ &- \frac{1}{\Delta y}\Bigl(E_z(y) - E_z(y-\Delta y)\Bigr)
\end{aligned}

したがって

\begin{aligned}
\mu\frac{H_x(t) - H_x(t-\Delta t)}{\Delta t} &= \frac{1}{\Delta z}\Bigl(E_y(z) - E_y(z-\Delta z)\Bigr)\\ 
&- \frac{1}{\Delta y}\Bigl(E_z(y) - E_z(y-\Delta y)\Bigr)
\end{aligned}

最新時刻の磁界$H$の$x$成分、$H_x(t)$について解くと、

\begin{aligned}
H_x(t) = H_x(t-\Delta t) &+ \frac{\Delta t}{\mu\Delta z}\Bigl(E_y(z) - E_y(z-\Delta z)\Bigr)\\ 
&- \frac{\Delta t}{\mu\Delta y}\Bigl(E_z(y) - E_z(y-\Delta y)\Bigr)
\end{aligned}
  • 式$(6)$について
    媒質の導電率を$\sigma$とすると、$j_x=\sigma E_x$ より、
\begin{aligned}
(左辺)=\epsilon\frac{\partial E_x}{\partial t} + \sigma E_x &\approx \epsilon\frac{E_x(t) - (1-\Delta t\sigma)E_x(t-\Delta t)}{\Delta t}\\
(右辺)=\frac{\partial H_y}{\partial z} - \frac{\partial H_z}{\partial y} &\approx \frac{1}{\Delta z}\Bigl(H_y(z) - H_y(z-\Delta z)\Bigr) \\ &- \frac{1}{\Delta y}\Bigl(H_z(y) - H_z(y-\Delta y)\Bigr)
\end{aligned}

したがって

\begin{aligned}
\epsilon\frac{E_x(t) - \Bigl(1-\frac{\Delta t\sigma}{\epsilon}\Bigr) E_x(t-\Delta t)}{\Delta t} &= \frac{1}{\Delta z}\Bigl(H_y(z) - H_y(z-\Delta z)\Bigr) \\ &- \frac{1}{\Delta y}\Bigl(H_z(y) - H_z(y-\Delta y)\Bigr)
\end{aligned}

最新時刻の電界$E$の$x$成分、$E_x(t)$について解くと、

\begin{aligned}
E_x(t) = \Bigl(1-\frac{\Delta t\sigma}{\epsilon}\Bigr)E_x(t-\Delta t) &+ \frac{\Delta t}{\epsilon\Delta z}\Bigl(H_y(z) - H_y(z-\Delta z)\Bigr)\\ 
&- \frac{\Delta t}{\epsilon\Delta y}\Bigl(H_z(y) - H_z(y-\Delta y)\Bigr)
\end{aligned}

電流要素$\boldsymbol{j}$の追加によって、右辺第一項に係数が現れたかと思います。

右辺第一項の係数$\Bigl(1-\frac{\Delta t\sigma}{\epsilon}\Bigr)$は、書籍によっては異なる形で記載されている場合がありますが、差分の取り方によって異なる形式で表示されるためです。本記事においては後進差分と呼ばれる一次精度の差分を用いておりますが、実際は精度の問題から二次精度の中心差分を使用するほうが望ましいです。（詳細は他書籍等ご参照ください。）

さて、電流要素$\boldsymbol{j}$が追加された理由ですが、前編においては真空中（非導体：無損失媒質）で計算する前提があったため、暗黙的に$\sigma=0, \boldsymbol{j}=\boldsymbol{0}$としていましたが、今回は損失（導体：電流を通す媒質）を含む媒体も含めて考慮するため、本項が現れています。

実践

今回も、前回と同じくらい面倒な立式パートでしたが、ここからは楽しい実践パートです。
立式パートで求めた方法と同様の手法で、式$(3) - (8)$をすべて差分展開すると、以下の式$(9) - (14)$を得ます。

\begin{aligned}
H_x(t) = H_x(t-\Delta t) &+ \frac{\Delta t}{\mu\Delta y}\Bigl(E_z(y) - E_z(y-\Delta y)\Bigr)\\ 
&- \frac{\Delta t}{\mu\Delta z}\Bigl(E_y(z) - E_y(z-\Delta z)\Bigr)&...(9)\\

H_y(t) = H_y(t-\Delta t) &+ \frac{\Delta t}{\mu\Delta z}\Bigl(E_x(z) - E_x(z-\Delta z)\Bigr)\\ 
&- \frac{\Delta t}{\mu\Delta x}\Bigl(E_z(x) - E_x(x-\Delta x)\Bigr)&...(10)\\

H_z(t) = H_z(t-\Delta t) &+ \frac{\Delta t}{\mu\Delta x}\Bigl(E_y(x) - E_y(x-\Delta x)\Bigr)\\ 
&- \frac{\Delta t}{\mu\Delta y}\Bigl(E_x(y) - E_x(y-\Delta y)\Bigr)&...(11)\\

E_x(t) = \Bigl(1-\frac{\Delta t\sigma}{\epsilon}\Bigr)E_x(t-\Delta t) &+ \frac{\Delta t}{\epsilon\Delta y}\Bigl(H_z(y) - H_z(y-\Delta y)\Bigr)\\ 
&- \frac{\Delta t}{\epsilon\Delta z}\Bigl(H_y(z) - H_y(z-\Delta z)\Bigr)&...(12)\\

E_y(t) = \Bigl(1-\frac{\Delta t\sigma}{\epsilon}\Bigr)E_y(t-\Delta t) &+ \frac{\Delta t}{\epsilon\Delta z}\Bigl(H_x(z) - H_x(z-\Delta z)\Bigr)\\ 
&- \frac{\Delta t}{\epsilon\Delta x}\Bigl(H_z(x) - H_z(x-\Delta x)\Bigr)&...(13)\\

E_z(t) = \Bigl(1-\frac{\Delta t\sigma}{\epsilon}\Bigr)E_z(t-\Delta t) &+ \frac{\Delta t}{\epsilon\Delta x}\Bigl(H_y(x) - H_y(x-\Delta x)\Bigr)\\ 
&- \frac{\Delta t}{\epsilon\Delta y}\Bigl(H_x(y) - H_x(y-\Delta y)\Bigr)&...(14)
\end{aligned}

自由空間

コード

では実際にコードを実装していきましょう。
前編同様にnumpyの行列計算を用いて、サクサク動かしていきましょう。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

dx = 0.01 # x方向空間差分間隔[m]
dy = 0.01 # y方向空間差分間隔[m]
dz = 0.01 # z方向空間差分間隔[m]

c = 2.99792458e8 # 光速[m/s]

dt = 0.99/(c * np.sqrt((1.0/dx ** 2 + 1.0/dy ** 2 + 1.0/dz ** 2))) # 時間差分間隔[s] c.f.Courantの安定条件
f = 1.0e9 # 周波数[Hz]

nx = 100 # x方向計算点数
ny = 100 # y方向計算点数
nz = 100 # z方向計算点数

nt = 250 # 計算ステップ数

# 電気定数初期化と更新係数の計算
eps = np.full((nx, ny, nz), 8.854187817e-12)
mu = np.full((nx, ny, nz), 1.2566370614e-6)
sigma = np.full((nx, ny, nz), 0.0)

dhx = dt / (mu * dx) # 式(9) - (11)の右辺係数
dhy = dt / (mu * dy) # 式(9) - (11)の右辺係数
dhz = dt / (mu * dz) # 式(9) - (11)の右辺係数

ce = 1.0 - dt * sigma / eps # 式(12) - (14)の右辺第一項係数
dex = dt / (eps * dx) # 式(12) - (14)の右辺第二項係数
dey = dt / (eps * dy) # 式(12) - (14)の右辺第二項係数
dez = dt / (eps * dz) # 式(12) - (14)の右辺第二項係数

# 電磁界初期化
t = 0.0

E_x = np.zeros(shape=(nx, ny, nz))
E_y = np.zeros(shape=(nx, ny, nz))
E_z = np.zeros(shape=(nx, ny, nz))
H_x = np.zeros(shape=(nx, ny, nz))
H_y = np.zeros(shape=(nx, ny, nz))
H_z = np.zeros(shape=(nx, ny, nz))

image_list = []
for _ in range(nt):
    
    # 電界のz成分を励振
    E_z[nx//2, ny//2, nz//2] = np.sin(2.0 * np.pi * f * t)
    
    t += dt/2
    
    # 電界各成分計算
    E_x = ce * E_x + dez * (H_z - np.roll(H_z, shift=1, axis=1))\
                   - dey * (H_y - np.roll(H_y, shift=1, axis=2))

    E_y = ce * E_y + dey * (H_x - np.roll(H_x, shift=1, axis=2))\
                   - dez * (H_z - np.roll(H_z, shift=1, axis=0))

    E_z = ce * E_z + dez * (H_y - np.roll(H_y, shift=1, axis=0))\
                   - dex * (H_x - np.roll(H_x, shift=1, axis=1))
    
    E_amp = np.sqrt(E_x ** 2 + E_y ** 2 + E_z ** 2)
    
    # 電界のz成分を励振
    E_z[nx//2, ny//2, nz//2] = np.sin(2.0 * np.pi * f * t)
    
    t += dt/2
    
    # 磁界各成分計算
    H_x = H_x + dhz * (E_z - np.roll(E_z, shift=-1, axis=1))\
              - dhy * (E_y - np.roll(E_y, shift=-1, axis=2))

    H_y = H_y + dhy * (E_x - np.roll(E_x, shift=-1, axis=2))\
              - dhz * (E_z - np.roll(E_z, shift=-1, axis=0))

    H_z = H_z + dhz * (E_y - np.roll(E_y, shift=-1, axis=0))\
              - dhx * (E_x - np.roll(E_x, shift=-1, axis=1))
    
    # 可視化
    fig = plt.figure()
    sns.heatmap(E_z[nx//2, :, :], cmap="Reds", vmax=0.1, vmin=0.0)
    fig.savefig("output/free/YZ/freeYZ_{}_map.png".format(str(_).zfill(6)))
    
    fig = plt.figure()
    sns.heatmap(E_z[:, :, nz//2], cmap="Reds", vmax=0.1, vmin=0.0)
    fig.savefig("output/free/XY/freeXY_{}_map.png".format(str(_).zfill(6)))
    
    plt.close('all')

計算結果

今回は上記コードで出力された複数のPNG形式の画像のうち、よさそうな時点を切り出してます。
動画っぽくする場合は上記コード実行の上、フリーソフト等でGIFに変換できます。

  • $XY$平面
    freeXY_000113_map.png

  • $YZ$平面
    freeYZ_000142_map.png

地味すぎて計算できているかよくわかりませんね。
ただ、発散とかはしていないので、多分大丈夫でしょう。
対数軸で見れば、多分比較的きれいな放射が確認できると思いますが、今回の本意ではないのでやりません。

半波長ダイポールアンテナ

コード

先程の計算結果はあまりにも地味で見栄えしないので、給電点周りに金属をおいてあげて、放射を手助けしてあげましょう。ここでは、半波長のダイポールアンテナを設置してみます。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

dx = 0.01 # x方向空間差分間隔[m]
dy = 0.01 # y方向空間差分間隔[m]
dz = 0.01 # z方向空間差分間隔[m]

c = 2.99792458e8 # 光速[m/s]

dt = 0.99/(c * np.sqrt((1.0/dx ** 2 + 1.0/dy ** 2 + 1.0/dz ** 2))) # 時間差分間隔[s] c.f.Courantの安定条件
f = 1.0e9 # 周波数[Hz]

nx = 100 # x方向計算点数
ny = 100 # y方向計算点数
nz = 100 # z方向計算点数

nt = 250 # 計算ステップ数

# 電気定数初期化と更新係数の計算
eps = np.full((nx, ny, nz), 8.854187817e-12)
mu = np.full((nx, ny, nz), 1.2566370614e-6)
sigma = np.full((nx, ny, nz), 0.0)

# 半波長ポールアンテナ配置
half_wavelength = int((c/f)//(2.0 * dz)) # アンテナz方向長さ(半波長)
eps[nx//2, ny//2, (nz - half_wavelength)//2 : (nz + half_wavelength)//2] = 156 * 8.854187817e-12
sigma[nx//2, ny//2, (nz - half_wavelength)//2 : (nz + half_wavelength)//2] = 5.76e7

dhx = dt / (mu * dx) # 式(9) - (11)の右辺係数
dhy = dt / (mu * dy) # 式(9) - (11)の右辺係数
dhz = dt / (mu * dz) # 式(9) - (11)の右辺係数

ce = (2.0 * eps - sigma * dt)/(2.0 * eps + sigma * dt) # 式(12) - (14)の右辺第一項係数
dex = 2.0 * dt /((2.0 * eps * dx) + (sigma * dt * dx)) # 式(12) - (14)の右辺第二項係数
dey = 2.0 * dt /((2.0 * eps * dy) + (sigma * dt * dy)) # 式(12) - (14)の右辺第二項係数
dez = 2.0 * dt /((2.0 * eps * dz) + (sigma * dt * dz)) # 式(12) - (14)の右辺第二項係数


# 電磁界初期化
t = 0.0

E_x = np.zeros(shape=(nx, ny, nz))
E_y = np.zeros(shape=(nx, ny, nz))
E_z = np.zeros(shape=(nx, ny, nz))
H_x = np.zeros(shape=(nx, ny, nz))
H_y = np.zeros(shape=(nx, ny, nz))
H_z = np.zeros(shape=(nx, ny, nz))

image_list = []
for _ in range(nt):
    
    # 電界のz成分を励振
    E_z[nx//2, ny//2, nz//2] = np.sin(2.0 * np.pi * f * t)
    t += dt/2
    
    # 電界各成分計算
    E_x = ce * E_x + dez * (H_z - np.roll(H_z, shift=1, axis=1))\
                   - dey * (H_y - np.roll(H_y, shift=1, axis=2))

    E_y = ce * E_y + dey * (H_x - np.roll(H_x, shift=1, axis=2))\
                   - dez * (H_z - np.roll(H_z, shift=1, axis=0))

    E_z = ce * E_z + dez * (H_y - np.roll(H_y, shift=1, axis=0))\
                   - dex * (H_x - np.roll(H_x, shift=1, axis=1))
    
    E_amp = np.sqrt(E_x ** 2 + E_y ** 2 + E_z ** 2)
    
    # 電界のz成分を励振
    E_z[nx//2, ny//2, nz//2] = np.sin(2.0 * np.pi * f * t)
    t += dt/2
    
    # 磁界各成分計算
    H_x = H_x + dhz * (E_z - np.roll(E_z, shift=-1, axis=1))\
              - dhy * (E_y - np.roll(E_y, shift=-1, axis=2))

    H_y = H_y + dhy * (E_x - np.roll(E_x, shift=-1, axis=2))\
              - dhz * (E_z - np.roll(E_z, shift=-1, axis=0))

    H_z = H_z + dhz * (E_y - np.roll(E_y, shift=-1, axis=0))\
              - dhx * (E_x - np.roll(E_x, shift=-1, axis=1))
    
    
    # 可視化
    fig = plt.figure()
    sns.heatmap(E_z[nx//2, :, :], cmap="Reds", vmax=0.1, vmin=0.0)
    fig.savefig("output/dipole/YZ/dipoleYZ_{}_map.png".format(str(_).zfill(6)))
    plt.close('all')
    
    fig = plt.figure()
    sns.heatmap(E_z[:, :, nz//2], cmap="Reds", vmax=0.1, vmin=0.0)
    fig.savefig("output/dipole/XY/dipoleXY_{}_map.png".format(str(_).zfill(6)))
    plt.close('all')

コードの変更点は、以下二点です。

  • アンテナ設置

アンテナを模擬するために高い導電率・比誘電率の導体を$z$軸方向に配置します。

# 半波長ポールアンテナ配置
half_wavelength = int((c/f)//(2.0 * dz)) # アンテナz方向長さ(半波長)
eps[nx//2, ny//2, (nz - half_wavelength)//2 : (nz + half_wavelength)//2] = 156 * 8.854187817e-12
sigma[nx//2, ny//2, (nz - half_wavelength)//2 : (nz + half_wavelength)//2] = 5.76e7
  • 中心差分化

後進差分で実施した際には、精度の問題からか発散してしまったため、中心差分形式に変更しました。

ce = (2.0 * eps - sigma * dt)/(2.0 * eps + sigma * dt) # 式(12) - (14)の右辺第一項係数
dex = 2.0 * dt /((2.0 * eps * dx) + (sigma * dt * dx)) # 式(12) - (14)の右辺第二項係数
dey = 2.0 * dt /((2.0 * eps * dy) + (sigma * dt * dy)) # 式(12) - (14)の右辺第二項係数
dez = 2.0 * dt /((2.0 * eps * dz) + (sigma * dt * dz)) # 式(12) - (14)の右辺第二項係数

計算結果

  • $XY$平面
    dipoleXY_000125_map.png

  • $YZ$平面
    dipoleYZ_000233_map.png

先程に比べて、高いレベルで電界が伝搬していますね。
なんだか形もいい感じで、アンテナが放射器としてよく働いていることが確認できました。

$YZ$平面の動画も示します。（重たいので一枚だけ）
dipole_YZ.gif

端面で怪しい動きをしていますが、境界条件をつけていないため、電界が全反射しており、あたかも金属壁が配置されているように計算されてしまっているためです。

まとめ

ひとまずこれで、現実の問題を解く上で必要な3次元FDTD法の実装ができました。

本来であれば、計算領域の端面に吸収境界を配置したり、もっと色々機能を追加していく上での説明の必要がありますが、長くなってしまうので割愛します。

次回以降は、実際にこちらの3次元FDTD法を用いて、現実的な以下の問題①, ②と、数値計算の醍醐味である③に挑んでいきたいと思います。

① FDTD法を用いたアンテナ性能解析
② 誘電体による電界の散乱界の逆問題
③ 解析精度・速度の向上

ようやく現実的な話ができる準備が揃ったので、次回以降はもう少し面白い話にできるのではないかと思っています。

参考文献（いつもの）

[1]宇野亨, 何一偉, 有馬卓司, "数値電磁界解析のためのFDTD法- 基礎と実践", コロナ社, 2016

本記事は、全編通して上記書籍を参考にしております。
理論的な部分から実装まで、詳細に記載されており、FDTD法の実装をする上では欠かせない名著です。
電磁界解析のみならず数値計算の入門者にとってもおすすめの一冊です。

8
12
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
8
12

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?