数値電磁界解析のサージ解析への適用

作成,更新

マクスウェルの方程式を直接解く数値電磁界解析をサージ解析へ適用します。本稿では,適用に必要な理論式を丁寧に記載し,MATLAB や Python で書いた FDTD のプログラム(コード)を掲載しています。

目標

  • 数値過渡電磁界解析手法の概念を理解する
  • FDTD 法を用いたサージ解析ができる

目次

数値電磁界解析

マクスウェル (Maxwell) の方程式から直接に数値的に解を求める数値電磁界解析手法について説明する。

有限差分時間領域法

マクスウェルの方程式をディジタルコンピュータを用いて時間軸で直接に解いていくために,方程式が空間と時間について離散され,結果として電磁界は空間において方程式を満たすように,各離散点に電磁界変数が配置された構造となる。

1960 年代半ば以降になると,汎用大型計算機が出現し,数値解析手法の具体化が各方面で進んだ。1966 年に,Yee が各電磁界変数を空間の離散点と 1 対 1 に対応づけ,マクスウェルの方程式におけるファラデー (Fraday) の方程式とアンペア (Ampere) の方程式を交互に用いた逐次計算による実用的な数値解析手法(Yee アルゴリズム,Yee's algorithm)を提案した。さらに,実用性を示すために長方形プレートの 2 次元散乱問題を解析した。この方法が有限差分時間領域法(FDTD : Finite Difference Time Domain 法)であり,電磁パルスの散乱,電磁波吸収問題などに適用され,発展してきた。

ここで,有限差分とは微分を差に置き換える,時間領域とは時間空間について問う,という意味である。

有限差分時間領域法(FDTD 法)の特徴

有限差分時間領域法(FDTD 法)の特徴の特徴は以下のとおりである。

  • マクスウェルの方程式を直接時間領域で解く方法である。
  • マクスウェルの方程式を満足するように電磁界変数が各空間離散点に 1 対 1 で対応づけられている。具体的には,電界の時間微分が磁界の回転に等しいアンペアの法則,磁界の時間微分が電界の回転に等しいファラデーの法則を満たす。
  • 計算時間ステップを $\Delta t$ として,電界と磁界を $\Delta t/2$ ごとに交互に計算するリープフロッグ演算を行う。
  • 時間領域で計算結果が得られ,直感的に理解しやすい。
  • 計算アルゴリズムは比較的簡単で理解しやすく,プログラミングが容易,ベクトル計算向き,かつ,並列計算向きのアルゴリズムであり,高速化が容易である。
  • 短所として,かなりの計算機資源(メモリの容量,CPU の処理時間)が必要となる。

リープフロッグ演算

計算時間ステップを $\Delta t$ として,電界 $\boldsymbol{E}$ と磁界 $\boldsymbol{H}$ を $\Delta t/2$ ごとに交互に計算するリープフロッグ演算のイメージを示す。電界 $\boldsymbol{E}$ と磁界 $\boldsymbol{H}$ は $\Delta t/2$ 時間ステップずれた時刻に配置される。これにより一時差分の中で最も精度がよい中心差分 (central differences) を用いて,マクスウェルの方程式の近似表現を与えることができる。これを Yee アルゴリズム (Yee's algorithm)という。

FDTD 法における電界と磁界の時間軸上での配置
図 FDTD 法における電界と磁界の時間軸上での配置

吸収境界

解析領域の境界では Yee アルゴリズムを使えないため,境界上の電磁界は別の方法で決めなければならない。解析領域が完全導体壁や完全磁気導体壁で囲まれている場合,壁面上の電界あるいは磁界の接線成分は 0 である。これに対して,電磁波の散乱問題やアンテナ解析などの開放領域の問題を扱う場合には,解析領域の反射が起こらないような仮想的な境界で閉じておく必要がある。この境界に課せられる無反射条件を吸収境界条件といい,これまでに数多くの手法が提案されてきた。

Mur の 1 次吸収境界条件 (Mur's first order absorbing boundary condition)

図のように $z=z_0$ を解析領域の外壁,すなわち吸収境界の一つとする。この境界 ABL に向かって $x$ 成分を持つ平面波が垂直に入射するものとすると,平面波の電界 $E_x = E_x (z+vt)$ と表され,次の微分方程式を満足する。

\[ \frac{\partial E_x}{\partial z} - \frac{1}{v}\frac{\partial E_x}{\partial t}=0 \]
1 次吸収境界条件
図 1 次吸収境界条件

吸収境界で反射がないなら,電界は形を保ったまま伝搬し,$z = z_0$ の吸収境界でも微分方程式が成り立つ。そこで,Yee アルゴリズムと同様に $z = z_0 + \Delta z/2$,$t = (n - 1/2)\Delta t$ で中心差分する。

\[ \frac{{E_x}^{n-1/2}(1)-{E_x}^{n-1/2}(0)}{\Delta z} - \frac{1}{v}\frac{{E_x}^n (1/2) - {E_x}^{n-1} (1/2)}{\Delta t}=0 \]

ここで,次式のように平均で近似する。

\[ {E_x}^{n-1/2} = \frac{{E_x}^n + {E_x}^{n-1}}{2} \] \[ {E_x}^n (1/2) = \frac{{E_x}^n (0) + {E_x}^n (1)}{2} \]

近似式を代入し,${E_x}^n (0)$ についてまとめる。

\[ {E_x}^n (0) = {E_x}^{n-1}(1)+\frac{v\Delta t - \Delta z}{v\Delta t + \Delta z}\{{E_x}^n (1) - {E_x}^{n-1}(0)\} \]

ここで,$\displaystyle c_z = \frac{v\Delta t - \Delta z}{v\Delta t + \Delta z}$ とおく。

\[ {E_x}^n (0) = {E_x}^{n-1}(1)+ c_z \{{E_x}^n (1) - {E_x}^{n-1}(0)\} \]

同様にして,$z = NZ$ では次式が成り立つ。

\[ {E_x}^n (NZ) = {E_x}^{n-1}(NZ-1)+ c_z \{{E_x}^n (NZ-1) - {E_x}^{n-1}(NZ)\} \]

マクスウェルの方程式

マクスウェルの 4 つの方程式を示す。ここで,$\boldsymbol{D}$ は電束密度 (electric flux density),$\boldsymbol{E}$ は電界 (electric field),$\boldsymbol{B}$ は磁束密度 (magnetic flux density),$\boldsymbol{H}$ は磁界 (magnetic field),$\rho$ は電荷密度 (charge density),$\boldsymbol{j}$ は電流密度 (current density) である。

\[ \boldsymbol{\nabla} \cdot \boldsymbol{D} = \rho \] \[ \boldsymbol{\nabla} \cdot \boldsymbol{B} = 0 \] \[ \boldsymbol{\nabla} \times \boldsymbol{E} + \frac{\partial \boldsymbol{B}}{\partial t} = 0 \] \[ \boldsymbol{\nabla} \times \boldsymbol{H} - \frac{\partial \boldsymbol{D}}{\partial t} = \boldsymbol{j} \]

電束密度と電界の関係,磁束密度と磁界の関係は,次の近似式で表される。ここで,$\varepsilon$ は誘電率 (permittivity),$\mu$ は透磁率 (permeability) である。

\[ \boldsymbol{D} \approx \varepsilon \boldsymbol{E} \] \[ \boldsymbol{B} \approx \mu \boldsymbol{H} \]

真空の場合,真空の誘電率 $\varepsilon_0$ = 8.854 × 10-12 [F/m],真空の透磁率 $\mu_0$ = 1.257 × 10-6 [H/m] を用いると,次式が成り立つ。

\[ \boldsymbol{D} = \varepsilon_0 \boldsymbol{E} \] \[ \boldsymbol{B} = \mu_0 \boldsymbol{H} \]

さらに,FDTD 法では,電界と電流を関係づける以下の式も利用する。ここで $\sigma$ は導電率 (conductivity),もしくは,電気伝導率と呼ばれる。

\[ \boldsymbol{j} = \sigma \boldsymbol{E} \]

(参考)ベクトル微分演算子 $\boldsymbol{\nabla}$

マクスウェルの方程式に用いた $\boldsymbol{\nabla}$ (ナブラ,nabla)は,スカラー関数に作用してベクトルを作るベクトル微分演算子であり,ハミルトン (Hamilton) の記号とも呼ばれる。

\[ \boldsymbol{\nabla} \equiv \frac{\partial}{\partial x}\boldsymbol{i} + \frac{\partial}{\partial y}\boldsymbol{j} + \frac{\partial}{\partial z}\boldsymbol{k} \]

ベクトル演算子を用いた式は,次のように展開できる。

\[ \boldsymbol{\nabla} \times \boldsymbol{E} = (\frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z})\boldsymbol{i} + (\frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x})\boldsymbol{j} + (\frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y})\boldsymbol{k} \] \[ \boldsymbol{\nabla} \times \boldsymbol{H} = (\frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z})\boldsymbol{i} + (\frac{\partial H_x}{\partial z} - \frac{\partial H_z}{\partial x})\boldsymbol{j} + (\frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y})\boldsymbol{k} \]

$z$ 方向へ進む 1 次元の電磁波

マクスウェル方程式の離散化

簡単のため,真空中を $z$ 方向に進む 1 次元の電磁波を考える。すなわち,$E_z=0$,$H_z=0$,$\displaystyle \frac{\partial}{\partial x}=0$,$\displaystyle \frac{\partial}{\partial y}=0$ をマクスウェルの方程式に代入する。

\[ -\frac{\partial E_y}{\partial z}\boldsymbol{i}+\frac{\partial E_x}{\partial z}\boldsymbol{j} =-\mu_0(\frac{\partial H_x}{\partial t}\boldsymbol{i}+\frac{\partial H_y}{\partial t}\boldsymbol{j}) \] \[ -\frac{\partial H_y}{\partial z}\boldsymbol{i}+\frac{\partial H_x}{\partial z}\boldsymbol{j} =\varepsilon_0(\frac{\partial E_x}{\partial t}\boldsymbol{i} + \frac{\partial E_y}{\partial t}\boldsymbol{j}) \]

上の 2 式の両辺を比較すると,以下の 4 式が得られる。

\[ \frac{\partial E_y}{\partial z} = \mu_0 \frac{\partial H_x}{\partial t} \] \[ \frac{\partial E_x}{\partial z} = -\mu_0 \frac{\partial H_y}{\partial t} \] \[ \frac{\partial H_y}{\partial z} = -\varepsilon_0 \frac{\partial E_x}{\partial t} \] \[ \frac{\partial H_x}{\partial z} = \varepsilon_0 \frac{\partial E_y}{\partial t} \]

上の 4 式を離散化し,リープフロッグ演算を適用できるようにする。

\[ {E_x}^{n+1}(i)={E_x}^{n}(i)-\frac{\Delta t}{\varepsilon_0 \Delta z}\{{H_y}^{n-\frac{1}{2}}(i+\frac{1}{2})-{H_y}^{n-\frac{1}{2}}(i-\frac{1}{2})\} \] \[ {E_y}^{n+1}(i)={E_y}^{n}(i)+\frac{\Delta t}{\varepsilon_0 \Delta z}\{{H_x}^{n-\frac{1}{2}}(i+\frac{1}{2})-{H_x}^{n-\frac{1}{2}}(i-\frac{1}{2})\} \] \[ {H_x}^{n+\frac{1}{2}}(i+\frac{1}{2})={H_x}^{n-\frac{1}{2}}(i+\frac{1}{2})+\frac{\Delta t}{\mu_0 \Delta z}\{{E_y}^{n}(i+1)-{E_y}^{n}(i)\} \] \[ {H_y}^{n+\frac{1}{2}}(i+\frac{1}{2})={H_y}^{n-\frac{1}{2}}(i+\frac{1}{2})-\frac{\Delta t}{\mu_0 \Delta z}\{{E_x}^{n}(i+1)-{E_x}^{n}(i)\} \]

解析ステップ

$z$ 方向へ進む 1 次元の電磁波を FDTD 法で解析するステップを示す。

  1. 電界,磁界の初期条件を与える。便宜上,時刻 0 の電界,時刻 $\Delta t/2$ の磁界とする。
  2. 時刻 $\Delta t$ の電界を,時刻 0 の電界,時刻 $\Delta t/2$ の磁界より求める。
  3. 吸収境界条件を電界に適用する。
  4. 時刻 $\Delta t/2$ の磁界を,時刻 $\Delta t$ の電界,時刻 $\Delta t/2$ の磁界より求める。
  5. 吸収境界条件を磁界に適用する。
  6. 以降,2. から 5. を繰り返す。

上記の FDTD 法で解析するステップをフローチャートで示す。

FDTD 法のフローチャート(1 次元解析)
図 FDTD 法のフローチャート(1 次元解析)

ガウシアンパルスの伝搬

初期条件が,次式で示すガウシアンパルスの電界 $E_x$ [V/m] で与えられる電磁波の伝搬を解析する。

\[ E_x = \exp(\frac{z-z_0}{\sigma^2}) \]

FDTD 法によりガウシアンパルスの伝搬を解析した結果を示す。横軸は $z$ [mm],縦軸は $E_x$ [V/m] または $H_y$ [A/m] として描画している。

FDTD 法によるガウシアンパルスの伝搬解析(Ex)
図 FDTD 法によるガウシアンパルスの伝搬解析($E_x$ [V/m])
FDTD 法によるガウシアンパルスの伝搬解析(Hx)
図 FDTD 法によるガウシアンパルスの伝搬解析($H_y$ [A/m])

解析コード

TEz モード

マクスウェルの方程式(アンペアの法則,ファラデーの法則)を示す。ここで,簡単のため電流源は 0 とした。

アンペアの法則(電界の時間微分が磁界の回転に等しい)
\[ \boldsymbol{\nabla} \times \boldsymbol{H} = \varepsilon \frac{\partial \boldsymbol{E}}{\partial t} \]
ファラデーの法則(磁界の時間微分が電界の回転に等しい)
\[ \boldsymbol{\nabla} \times \boldsymbol{E} = -\mu \frac{\partial \boldsymbol{H}}{\partial t} \]

ここで,$\boldsymbol{H} = H_z \boldsymbol{z}$,$\displaystyle \frac{\partial}{\partial z} = 0$ の TEz モード (transverse electric mode) を考える。

\[ \frac{\partial H_z}{\partial t} = -\frac{1}{\mu}\frac{\partial E_y}{\partial x} - \frac{1}{\mu}\frac{\partial(-E_x)}{\partial y} \] \[ \frac{\partial E_y}{\partial t} = -\frac{1}{\varepsilon}\frac{\partial H_z}{\partial x} \] \[ \frac{\partial (-E_x)}{\partial t} = -\frac{1}{\varepsilon}\frac{\partial H_z}{\partial y} \]

有限差分を適用するため,TEz モードにおいて,電界 $\boldsymbol{E}$ と磁界 $\boldsymbol{H}$ を半セルだけずらして配置したイメージを示す。これにより,空間の偏微分は中心差分が適用できる。

図 TEz モード (transverse electric mode)

有限差分の適用

マクスウェルの方程式に有限差分法を適用する。添え字 $t + \Delta t$ をつけたもの以外は,時刻 $t$ の電界,磁界を表す。

アンペアの法則
\[ \frac{{H_z}^{t+\Delta t} - H_z}{\Delta t} = -\frac{1}{\mu}\frac{E_y (i+1, j + \frac{1}{2}) - E_y (i,j+\frac{1}{2})}{\Delta y} -\frac{1}{\mu}\frac{-E_x (i+\frac{1}{2}, j + 1) + E_x (i + \frac{1}{2},j)}{\Delta x} \]

上式を ${H_z}^{t+\Delta t}$ について整理する。時刻 $t + \Delta t$ の磁界 ${H_z}^{t+\Delta t}$ は,時刻 $t$ の電界 $E_x$,$E_y$,磁界 $H_z$ より求めることができる。

\[ {H_z}^{t+\Delta t} = H_z -\frac{\Delta t}{\mu\Delta y}\{E_y (i+1, j + \frac{1}{2}) - E_y (i,j+\frac{1}{2})\} -\frac{\Delta t}{\mu\Delta x}\{-E_x (i+\frac{1}{2}, j + 1) + E_x (i + \frac{1}{2},j)\} \]
ファラデーの法則
\[ \frac{{E_y}^{t + \Delta t} - E_y}{\Delta t} = -\frac{1}{\varepsilon}\frac{H_z(i+\frac{1}{2}, j+\frac{1}{2})-H_z(i-\frac{1}{2},j+\frac{1}{2})}{\Delta x} \] \[ \frac{-{E_x}^{t + \Delta t} + E_x}{\Delta t} = -\frac{1}{\varepsilon}\frac{H_z(i+\frac{1}{2}, j+\frac{1}{2})-H_z(i+\frac{1}{2},j-\frac{1}{2})}{\Delta y} \]

上式を ${E_y}^{t + \Delta t}$,${E_x}^{t + \Delta t}$ について整理する。時刻 $t + \Delta t$ の電界 ${E_y}^{t+\Delta t}$,${E_x}^{t+\Delta t}$ は,時刻 $t$ の電界 $E_x$,$E_y$,磁界 $H_z$ より求めることができる。

\[ {E_y}^{t + \Delta t} = E_y -\frac{\Delta t}{\varepsilon \Delta x}\{H_z(i+\frac{1}{2}, j+\frac{1}{2})-H_z(i-\frac{1}{2},j+\frac{1}{2})\} \] \[ {E_x}^{t + \Delta t} = -E_x +\frac{\Delta t}{\varepsilon \Delta y}\{H_z(i+\frac{1}{2}, j+\frac{1}{2})-H_z(i+\frac{1}{2},j-\frac{1}{2})\} \]

参考文献

  1. サージ現象に関する数値電磁界解析手法調査専門委員会 編,「数値過渡電磁界解析手法 ――サージ現象への適用――」,オーム社,2008年3月25日(初版 第 1 刷)
  2. 宇野 亨 編著,何一偉・有馬卓司 著「数値電磁界解析のためのFDTD法-基礎と実践-」,コロナ社,2016年5月25日(初版 第 1 刷)
  3. 過渡数値電磁界解析手法の応用調査専門委員会 編,「FDTD 法を中心とした数値電磁界解析手法によるサージ解析の最新動向」,電気学会技術報告 第 1216 号,2011年3月
  4. 一般社団法人電気学会 電力・エネルギー部門 用語解説「第8回テーマ : FDTD 法

Qiita の FDTD に関する良質な記事

inserted by FC2 system