電磁界シミュレーションの概要と基礎原理－簡単な1次元問題による説明

電磁界シミュレーションの概要と基礎原理
－簡単な１次元問題による説明－
Basic algorithm of electromagnetic-field simulators
－Fundamental understanding using one-dimensional cases－
平野拓一
Takuichi HIRANO
東京工業大学大学院理工学研究科〒152-8552 東京都目黒区大岡山 2-12-1-S3-19
Graduate School of Science and Engineering, Tokyo Institute of Technology
2-12-1-S3-19 Ookayama, Meguro-ku, Tokyo, 152-8552 Japan
E-mail: [email protected]
Abstract
Many electromagnetic (EM) simulators, whose purpose is to solve Maxwell’s equations,
are useful these days thanks to the progress of computational resources. It is not necessary for
engineers to understand EM analysis algorithms. However, it is desirable for engineers to
understand the underlying EM analysis algorithms in order to use EM simulators effectively.
This tutorial lecture presents various EM analysis algorithms, such as the method of
moments (MoM), the finite-difference time-domain (FDTD) method, the finite-difference
frequency-domain (FDFD) method, and the finite element method (FEM), using a
one-dimensional problem for young and/or new engineers. The formulation is very easy in
one-dimensional problem. It is also possible to choose suitable EM simulator for a problem to
be solved if one has knowledge about the EM simulation algorithm. The objective of this
paper is helping engineers to understand difference of EM analysis algorithms.
本基礎講座では、電磁界分野の初学者のために 1
次元問題[3]を用いて各種電磁界解析手法、特にモー
メント法 (MoM; Method of Moments) [4], FDTD
(Finite-Difference Time-Domain) 法 [5], FDFD
B
(ファラデーの法則)
(1)
E  
(Finite-Difference Frequency-Domain)法[6], 有限要素
t
法 (FEM; Finite Element Method)[7]について紹介する。
D
1 次元問題は定式化が簡単なため、自分でプログラム
(アンペアの法則)
(2)
H  i 
t
を作成して確認することが容易である。文章で説明
を確立して以来、多くの規範問題が解かれて電磁
できないことを補うため、そして実際にプログラム
界現象が明らかとなった。しかしながら、現実に近
を動かして勉強するために Mathematica Ver.8 による
い問題を解くためには数値的手法に頼らざるを得ず、解析プログラム例を付録に掲載した。電磁界分野の
種々の電磁界シミュレーションアルゴリズムが提案
初学者の一助となれば幸いである。
された。さらに、その後のコンピュータの性能の大
きな飛躍によって電磁界シミュレーション技術は実
2. 解析モデル
用的なものになり、多くの電磁界シミュレータが市
図 1に例題として取り扱う１次元問題の解析モデ
販されるようになった。
ルを示す。電磁界は z 方向にのみ変化し、x, y 方向に
電磁界シミュレータが発展して便利になり、アル
一様である。 z1  z  z 2 の領域のみ比誘電率は
ゴリズムを知らなくても解析できるようになったが、
 r  4 (SiO2 のガラスを想定)であり、それ以外は真空
電磁界シミュレータのアルゴリズムを知らなければ
である。また z  z s には x 方向を向いた x, y 方向には
上手く使いこなせないのが現状である。電磁界シミ
ュレータの目的は、マクスウェルの方程式（電磁界
無限に広がった面電流源 J がある(電界は x 方向、磁
現象の支配方程式である微分方程式）を指定した境
界は y 方向となる)。
界条件の下で解くこと（境界値問題）である[2]。解
析アルゴリズムによって数値的扱いが異なり、それ
ぞれ問題によって長所と短所がある。
1. はじめに
マクスウェルがアンペアの法則に変位電流を導入
して電磁波の基礎方程式[1]
ˆ
 r  1 J  z ( z  z s )


 0
y z
r  4
r  1
x
r  1
 0
y
z2
z1
3.1. モーメント法 (MoM)
II
として放射電磁界を計算すればよい。領域 II では
 J x1 ,M y1 , J x 2 , M y 2 がガラス中にあるとして放射電
 J x2
J x1
磁界を計算すればよい。領域 III では  J x 2 ,M y 2 およ
J x2
M y1
J
III
 J x1
び真電流 J が真空にあるとして放射電磁界を計算す
ればよい。
 M y2
z
 M y1 M y 2
y
J  xˆJ ( z)
図 2 モーメント法の領域分割
x
II
I
 J x1
J x1
M y1
J x2
J
III
 J x2
J
 M y2
z
x
x
I
全部真空のモデル
M  yˆM y ( z )
H
H
E
E
z
y
z
y
x
E
H
E
 J x2
 M y2
x
H
x
III
 M y1 M y 2
y
(6)
うにして重み係数を求めるための連立一次方程式を
解く問題に帰着させている（これは３次元の解析で
も同様である）。電界および磁界の接線成分が等しい
（境界条件）こと、および境界の反対側では等価定
理よりの法線ベクトルが逆となるので異符号の電磁
流が流れる。このように媒質の境界面に電磁流を置
くと、各領域の電磁界の解析は真電流 J および、媒
質境界に仮定した等価電磁流からの一様媒質内での
放射で計算することができる（界等価定理より）。図
3に示すように、領域 I では J x1 , M y1 が真空中にある
3. 種々の解析法のアルゴリズムの説明
本節では例として前節で説明した１次元問題を取
り上げ、種々の解析法のアルゴリズムについて説明
する。
I
M2  M y 2 yˆ ( z  z2 )
(5)
り、 J x1 がその基底関数に対する重みである。このよ
図 1 1 次元問題の解析モデル
x
(4)
例えば、 J1 においては xˆ ( z  z1 ) が基底関数であ
z
zs
M1  M y1 yˆ ( z  z1 )
J 2  J x 2 xˆ ( z  z2 )
(a) 面電流
z
y
(b) 面磁流
図 4 面電磁流からの放射
II
 J x1
J x1
M y1
次に、この計算に必要になる面電磁流からの放射
を計算しておく。図 4 (a)に示す面電流からの放射界
は次式となる(付録 A.1 参照)。
J x2
z
y
全部真空のモデル
y
 M y1 M y 2
z
全部ガラスのモデル
図 3 各領域の電磁界の計算モデル
モーメント法(Method of Moments; MoM)[4]は周波
数領域の解析法であり、界等価定理(Field Equivalence
Thorem)[2]を基本原理とする。界等価定理によると異
なる媒質の境界（異ならなくてもよいが）に等価電
磁流を仮定すると、各領域は一様な媒質で満たされ
ていると考えることができる。図 1の問題では図 2
のように媒質境界に J1 , M1 , J 2 , M 2 の４つの電流お
よび磁流を仮定する。それらは、既知の基底関数（ま
たは展開関数）に対する重み J x1 , M y1 , J x 2 , M y 2 で次の
ように表現する。
J1  J x1 xˆ ( z  z1 )
(3)

 xˆ
E
 xˆ

J jkz
e
2
J  jkz
e
2
 J jkz
 yˆ e
, H 2
J
 yˆ e  jkz
( z  0)
2

( z  0)
( z  0)
(7)
( z  0)
また、図 4 (b)に示す面磁流からの放射界は次式と
なる。
 M jkz
e
 xˆ
E 2
M  jkz
 xˆ
e
2

M jkz

 yˆ 2 e
, H
M  jkz
 yˆ
( z  0)
e

2
( z  0)
( z  0)
(8)
( z  0)
一般に、電磁流から電界および磁界を計算すると
きは畳み込みの表現となるが、そのときに電磁流に
掛けて（ベクトル的には内積を取って）畳み込みを
行う関数は電磁流のインパルス応答に相当し、グリ
ーン関数とよばれる。
４つの等価電磁流および励振面電流源が放射する
電磁界を各領域で計算し、 z  z1 , z 2 の境界両側にお
いて電界および磁界の接線成分が等しいという方程
式を立てると次式が得られる。
Et z  z   Et z  z 
1
1

H t z  z   H t z  z 
1
1
(9)

E

E
t z  z2
 t z  z2

H t z  z2  H t z  z y
ここで、 Et , H t の下添え字 t は境界の接線

1
3.2. FDTD 法
E n1 H
は４つなので解くことができる。１次元では積分が
簡単であったが、２次元、３次元の問題では被積分
関数にある未知電磁流を求める積分方程式を行列方
程式に帰着させて解く手法である。モーメント法で
は FDTD 法や有限要素法とは違って放射条件は吸収
境界条件などの特殊な処理を用いる必要がない。
En
H
n
1
2
E n1 H
n
3
2
time
t
x
E(i  1)
E(i  1)
E(i)
1
H(i  )
2
y
1
H(i  )
2
(tangential)成分を表している。 z , z はそれぞれ
モデルで計算する。
１次元問題の場合には式(9)の連立方程式から解が
得られるが、２次元、３次元の問題では式(7)を成立
させる観測点は境界上の全ての場所としなければな
らいため、境界面上で定義された任意の関数で重み
付けして境界面で積分し、これを重み付け残差法と
言う。また、重み関数を基底関数と同一に選ぶ場合、
ガラーキン法と言う。１次元の場合は式(3)～(6)の基
底関数で重み付けすると、次式のように x,y 成分を抽
出する演算となる。
Ex z  Ex z
1
 1
H y   H y 
z1
z1

(10)

 E x z2  E x z2

 H y z2  H y z y

未知数の数は J x1 , M y1 , J x 2 , M y 2 の４つ、方程式の数
1
2
t / 2

1
z  z1 のすぐ左側と右側を表している。同様に z 2 , z 2
はそれぞれ z  z2 のすぐ左側と右側を表している。
z1 の電磁界は領域 I のモデルで、 z1 および z 2 の電

磁界は領域 II のモデルで、 z 2 の電磁界は領域 III の
n
z / 2
z
3
H(i  )
2
z
図 5 FDTD 法の時間および空間の離散化
FDTD(Finite-Difference Time-Domain)法[5][8-11]は
時間領域の解析法である。計算機で微分を扱えるよ
うにマクスウェルの方程式を時間および空間の両方
に対して差分化（差分は微分とは異なり、微小だが
有限な範囲での値の変化）して陽的（行列方程式を
解かないで漸化式を逐次計算する）に解く手法であ
る。図 5に時間および空間の離散化を示す。電界と
磁界の更新は半ステップずらして交互に行われる。
また、電界と磁界は回転の差分計算のために空間的
に半セルずれて配置されている。
式(1)-(2)において i  E  i e ( E は導電電流、i e は
電流源による励振電流)を代入すると次式が得られ
る。
H
0
t
E
H 
 E  i e
t
E  
(11)
(12)
空間、時間の差分化はどちらが先でも構わないが、
先ず図 5を見て次のように時間の差分化を行う。式
(11), (12)を時間微分の項について整理すると次式と
なる。
H
1
(13)
  E
t

E 1

1
(14)
   H  E  ie
t 


さらに、時間間隔 t で時間を差分化し、時間ステッ
プ nt のときの式 (13) および時間ステップ
(n  1/ 2)t のときの式(14)を書くと次式となる。
H
n
1
2
H
t
n
1
2

1

  En
(15)
n
En  En 1
 n 1
1 n 1
  E 2    H 2  i e 2 (16)
t



1
1
ここで、次の中央差分による時間差分化が行われた。
1
n
2
E
t

En  En 1
t
H n H

t
n
1
2
H
t
(17)
n
1
2
(18)
電磁界の上添え字は時間ステップ nt の t を除い
た部分を記したものである。式(16)において電界
E
E
n
1
2
1
n
2
のような時間配置をしていないので

En  En 1
で近似する。
2
n
E n  E n1
 E n  E n1 1
1 n 1

   H 2  i e 2 (19)
t

2


1
n
1
2
n
式(15)を H について、式(19)を E について解く
と次式が得られる。
t
1
1
2 E n1  t /    H n 2  t /  i n 12 (20)
En 
t
t
t e
1
1
1
2
2
2
1
1
n
n

t
(21)
H 2  H 2    En

式(20)からは最新時間の電界は 1 ステップ前の電
界および半ステップ前の磁界の回転と電流から求ま
ることがわかる。式(21)からは最新時間の磁界は１ス
テップ前の磁界と半ステップ前の電界の回転から求
まることがわかる。この演算を場所ごとに行えばよ
いのであるが、空間に非常に密に電磁界を配置する
わけにはいかないので、空間も図 5のように電界と
磁界の配置を半セルずらして差分化する。任意のベ
クトル A の回転の差分は次式のようになる。
xˆ
yˆ
zˆ
  A   / x  / y  / z
Ax
Ay
Az
 A Ay   Az Ax   Ay Ax 
  yˆ 

 xˆ z 


  zˆ
(22)
z   x
z   x
y 
 y
 A (i, j  1, k )  Az (i, j , k ) Ay (i, j , k  1)  Ay (i, j , k ) 
 xˆ  z


y
z


 A (i  1, j , k )  Az (i, j , k ) Ax (i, j , k  1)  Ax (i, j , k ) 
 yˆ  z


x
z


A
(
i

1
,
j
,
k
)

A
(
i
,
j
,
k
)
 y
A (i, j  1, k )  Ax (i, j , k ) 
y
 zˆ 
 x

x
y


式(20), (21)に空間差分化を考慮し、導電率を 0 とお
くと次の１次元問題の電磁界の更新式が得られる。
t
2 E n1 (i )
E xn (i ) 
t x
1
2
1
1
n
n
t /  
1
1 

 H y 2 (i  )  H y 2 (i  )

t 
2
2 
1
2
1
t / 
n

iex 2 (i )
t
1
2
1
1
n
n
1
1
H y 2 (i  )  H y 2 (i  )
2
2
t

E n (i  1)  Exn (i )
z x
1


(23)
(24)
全位置 i について式(23)で電界を更新し、次に時間
を半ステップ進めて式(24)で磁界を更新する。更に時
間を半ステップ進めて同様に式(23)で電界を更新
し・・・と繰り返すことになる。行列方程式を解く
必要が無いのでプログラムは簡単であるが、注意と
しては時間ステップには
(25)
vt  z
の制約があり、時間ステップは任意に大きく取るこ
とができない。ここで、 v は媒質中の電磁波の速度
であり、いろいろな媒質があるときは一番速い媒質
に制約を受ける。式(25)の条件を満たさなければ時間
ステップを進めると数値的に発散してしまう。この
時間ステップの制約条件は提唱者の名前を取って
CFL (Courant-Friedrichs-Lewy)条件とよばれる。CFL
条件を満たす最大時間ステップで規格化した時間ス
テップは CFL 数といわれ、工夫していない FDTD 法
では時間ステップは 1 CFL 数以下にしなければなら
ない。FDTD 法はこのように任意の時間波形を入射
させることができるので、パルスを印加して散乱波
の時間波形を得て、それをフーリエ変換することで
広い帯域の周波数特性を得ることができる。単一周
波数を印加して解析を行うこともできるが、定常状
態に達するまで時間ステップを進めなければならな
い。
次に、境界条件の処理について説明する。電気壁
(PEC)の場合には電界の接線成分が 0 になるので、電
界を更新した直後に電界の接線成分を 0 と書き換え
る。
磁気壁(PMC)の場合は磁界の接線成分が 0 になる
ので、磁界を更新した直後に磁界の接線成分を 0 に
書き換える。放射境界条件（左が吸収境界の場合）
は進行波の条件を満たす波動方程式 E x  E x ( z  ct ) ,
Ex 1 Ex

 0 から導出した次式を用いる方法があ
z c t
り、これは Mur１次の吸収境界条件といわれる。
ct  z n1
Ex (i  1)  Exn (i) (26)
ct  z
磁界も同様な形の式となるが、電界と磁界は独立
ではないので、電界あるいは磁界のどちらかに対し
て境界条件の処理を適用すればよい。
FDTD 法解析プログラムの流れを図 6に示す。
Exn1 (i)  Exn (i  1) 


・構造、媒質の設定
・メッシュ生成
・時間ステップの決定
・更新式の係数の計算
Loop {
・電界の更新
・電界の境界条件の処理
・磁界の更新
・磁界の境界条件の処理
・吸収境界条件の処理
}
#i
i
z 2i 1 z1
z1i 1
E(i  1)
#(i+1)
z 2i z1i 1
z
z 2i 1
z  z1i
z 2i  z
i
e

2
z2i  z1i
z 2i  z1i
e1i 
zi
i-th
zi
1
element
2
z
図 7 FEM の空間離散化と基底関数
 jzEx (i)  H y (i)  H y (i  1)  ziex (i) (28)
次に、吸収境界条件の方程式を導出する。  z 方向
Ex 1 E x
に速度 c で進む波動は

 0 を満たす。こ
z c t
れを差分化すると
(29)
となる。磁界も同様なので
有限要素法法(FEM)[7]は周波数領域の解析法であ
る。まずマクスウェルの方程式の 2 式から電界のみ
に対する波動方程式を導出する（磁界の波動方程式
でもよいが、電界がよく用いられる）。そして変分原
理に基づいて波動方程式の汎関数を求め、汎関数の
極値探索問題を数値的に解く手法である。ところが、
汎関数は一般には求まらないため、次式のように既
知の重み関数 W で重み付けした弱形式が提案され
た。



E
 M 
  k 0 2  r E  jk 0 0 J      dv  0
W     

 r 
 r 

(32)
1

 M  
2
   r   W     E  k 0  r W  E  W   jk0 0 J      r   dv



1
nˆ  W     EdS  0
 

(30)
となる。  z 方向に進む波を吸収したい場合は式
(29),(30)において c  c と置き換えればよい。式
(27)-(30)を用いて空間に配置した電磁界の値すべて
を未知数とする行列を作成すると次式を得る。
式 (27),式 (29)  E x    0 




式 (28),式 (30)  H y   iex 
E(i)
#(i-1)
y
1
3.3. FDFD 法
FDFD 法[6]では FDTD 法と同様の空間メッシュ分
割を用いるが、周波数領域で行列方程式を解く手法
である。式(11), (12)において  / t  j とおいて、空
間は FDTD 法と同様に差分化すると次式を得る。
(27)
Ex (i  1)  Ex (i)  jzH y (i)  0
 jz  2cH y (nz  1 / 2)
  jz  2c H y (nz  3 / 2)  0
E(i  1)
x
e*
図 6 FDTD 法解析プログラムの流れ
 jz  2c E x (i)   jz  2c E x (i  1)  0
3.4. 有限要素法 (FEM)
(31)
上式を解くと周波数領域で空間の電磁界が求まる。
行列サイズは非常に大きいが、非常に疎な行列でも
あり、疎行列に特化したソルバ―を使うことができ
る。
r
(33)
W は任意の関数だが、特に電界の基底関数と同一
のものを用いる場合はガラーキン法といわれる。１
次元問題の場合には弱形式は次のようになる。
1 Wx E x
2
 k0  rWx E x
z
r z
 [


1
 Wx  jk0 0 J x 
r

 M z M y  

 ]dz

z  
 y
(34)
z2
1
E 
  Wx x   0
z  z  z
 r
1
図
7に示すように要素 i 内の電界は
E ( z )  A1i e1i ( z )  A2i e i2 ( z ) で表される。全空間の電
i
界は
Ne
Ne
i 1
i 1
E( z )   Ei ( z )   A1ie1i ( z )  A2i ei2 ( z ) 
(35)
で表される。電界の接線成分の境界条件より、
A2i  A1i 1 となる。行列方程式を作成した際にこの条
件を無条件で満たし、必要最小限の未知数の数とな
るように図 8に示すような屋根形の基底関数(ルー
フトップ基底関数)を定義する。
N basis
Ne
n 1
i 1
 Anb n ( z )   Ai ei21 ( z )  e1i ( z )
E( z ) 
A2b 2 (b 2  e12  e12 )
(36)
Anbn (bn  en2 1  e1n )
A1b1 (b1  e11 )
Am b m
x
e11
y
1
Element
e n2 1 e1n
e12 e12 e 22
1
2
3
2
m
3
z
n
m
n
図 8 屋根形基底関数
式(36)を式(34)に代入し、Wx  bmx とすると次式を
得る。
Nbasis
1
An 

n 1
 r

b ( z ) bnx ( z )
2
dz  k0  r  bmx ( z )bnx ( z )dz 
 mxz

z

bNbasis, x ( z Ne 1 )
1
[
ANbasisbNbasis, x ( z Ne 1 )
 r ( z Ne 1 )
z

b ( z )
1
A1b1x ( z1 ) 1x 1 ]
 r ( z1 )
z
[boundary]

1
   bmx ( z ) jk0 0 J x 

r

 M y  
 
 dz
 z  
[source]
(37)
ここで、 b n  e n21  e1n より、
 b ( z )b ( z )dz   (e
  e ( z )e ( z )  e

mx

nx
m 1
2x

n 1
2x
m1
2x
m1
2x
( z )  e1mx ( z ))(e2nx1 ( z )  e1nx ( z ))dz
n 1
2x

( z )e ( z )  e ( z )e ( z )  e ( z )e ( z ) dz
n
1x
m
1x
m
1x
n
1x
  mn E   m1,n E   m ,n1 E   mn E
m
22
m
21
m
12
m
11
(38)
b ( z ) bnx ( z )
e m1 ( z ) e1mx ( z ) e2nx1 ( z ) e1nx ( z )
dz   ( 2 x

)(

)dz
 mxz

z
z
z
z
z
m1
n 1
m1
n
m
n1
 e ( z ) e2 x ( z ) e2 x ( z ) e1x ( z ) e1x ( z ) e2 x ( z ) e1mx ( z ) e1nx ( z ) 
   2x



dz

z
z
z
z
z
z
z 
 z
  mn F22m   m1,n F21m   m,n1 F12m   mn F11m
(39)
e
ここで、  mn はクロネッカーのデルタであり、 E ij ,
Fije は次の公式で計算できる。
z 2e
( z e  z e ) / 3 (i  j  1,2)
Eije   e eixe ( z )eejx ( z )dz   2e 1e
z1
( z2  z1 ) / 6 (i  j )
(40)
z 2e
Fije   e
z1
1/( z  z ) (i  j  1,2)
e ( z ) e ( z )
dz  
z
z
1/( z  z ) (i  j )
e
ix
e
jx
e
2
e
1
e
1
e
2
(41)
また、これらの積分は同一要素の中の基底関数 e 同
士の場合にしか値を持たない。そのため、プログラ
ムを作成する場合には要素のループで、要素内の行
列方程式を作成して、それらの要素を系行列の要素
に足してもよい。全体の行列方程式は次のようにな
る。
 Ne 2  1 e

2
e
   Eij  k0  r Fij   [boundary] { An }  {[source]}

 e1 i , j 1  r

(42)
これを解くと式(36)の係数が求まり、電界分布が求ま
る。有限要素法では重み関数と一致する場所の基底
関数としか行列要素の値を生成しないために疎行列
となる。モーメント法とは異なり、空間全体にメッ
シュを切る必要があるが、行列が疎になるので疎行
列に特化したソルバを使うと高速なので行列のサイ
ズだけでは優劣を比較できない。磁界はファラデー
の法則 H  
E
より計算できる。
j
次に、境界条件について説明する。電気壁(PEC)
の場合は電界の接線成分は 0 で既知の値となるので、
その部分の基底関数の重み係数を 0 で既知とすれば
よい。磁気壁(PMC)の場合は電界の接線成分の境界の
法線方向の偏微分が 0 なので、式(34)より境界成分の
寄与が 0 となるだけで、何も処理する必要がない(式
(37)の[boundary]の項を入れる必要がない)。次に、
Ex
 Ex  
z
(43)
で定義されるインピーダンス境界の場合、式(37)中の
[boundary]の項は次のようになる。
1
[boundary ]   

 r

1
r (Ne )
ANbasis

e2Nxe ( z N e 1 )  1  e11x ( z1 ) 
 AN basis

 A1

z
z 

  r 

1
Ne
1
 N e   N e e2 x ( z N e 1 ) 
A1 1   1e1x ( z1 )
 r (1)



(44)

次に、吸収境界条件について説明する。  z 方向に
速度 c (波数 k )で進む波動(進行波, 素波)は
Ex
 jkEx  0
z
(45)
となるので、インピーダンス境界の式(44)において
   jk ,   0 と置いたものに等しい。
4. 計算例
各種解析方法で計算した図 1の構造の電界分布お
よび磁界分布をそれぞれ図 9、図 10に示す。周波数
2.45 GHz、自由空間波長を z1  0.250 , z 2  0.75 0 ,
z s  1.5 0 , FDTD, FDFD, FEM のメッシュサイズは
z  0 / 40 とした。
350
180
MoM
FDTD
300
120
FDFD
FEM
250
Arg(H) (deg)
60
|E| (V/m)
200
150
0
-60
MoM
100
FDTD
-120
50
Source
 r=4
0
0
0.02 0.04 0.06 0.08
FDFD
0.1
0.12 0.14 0.16 0.18
0.2
0.22 0.24
0
0.02 0.04 0.06 0.08
Position (m)
(a) Amplitude
0.12 0.14 0.16 0.18
0.2
0.22 0.24
(b) Phase
図 10 磁界分布
MoM
FDTD
120
5. まとめ
FDFD
FEM
1 次元問題を用いて電磁界解析手法であるモーメ
ント法(MoM), FDTD 法, FDFD 法, 有限要素法(FEM)
の概要と基礎原理について紹介した。表 1に各種解
析手法の特徴をまとめた。全解析手法において結果
は良好に一致した。また、勉強の一助となるように
参考として Mathematica Ver.8 による解析プログラム
を付録に掲載した。
60
Arg(E) (deg)
0.1
Position (m)
180
0
-60
-120
Source
 r=4
-180
0
0.02 0.04 0.06 0.08
0.1
0.12 0.14 0.16 0.18
0.2
0.22 0.24
Position (m)
表 1 各種電磁界解析手法の特徴
(b) Phase
図 9 電界分布
解析法
モーメント法
(MoM)
全解析手法の結果は良好に一致した。本解析モデ
ルは 1 次元問題なので、MoM の解は厳密解と一致す
る。振幅分布より、z1  z  z s の範囲で定在波が生じ
FDTD 法
ているのがわかる。また、磁界の定在波分布では電
界とは山と谷が逆になる。FDTD 法では正弦波入射
させ、100 周期後に各場所において１周期間の波形か
ら複素フーリエ係数を計算して複素表現を得ている。
1
有限要素法(FEM)
MoM
0.9
FDTD
0.8
FDFD
FEM
0.7
|H| (A/m)
FEM
Source
 r=4
-180
0.6
特徴
 周波数領域
 メッシュ：物体の境界、メッシュ形状は三角
パッチ（柔軟性は高い）
 陰解法（密行列）
 放射、散乱問題が得意
 時間領域
 メッシュ：空間全体、直方体格子（柔軟性は
高くない）
 陽解法（前のステップの電磁界の値を用いて
値を更新）
 人体など、複雑な組成・構造の解析が得意
 放射問題のためには吸収境界条件が必要
 周波数領域
 メッシュ：空間全体、四面体要素（柔軟性は
高い）
 陰解法（疎行列）
 放射問題のためには吸収境界条件が必要
0.5
文
0.4
0.3
0.2
0.1
Source
 r=4
0
0
0.02 0.04 0.06 0.08
0.1
0.12 0.14 0.16 0.18
Position (m)
(a) Amplitude
0.2
0.22 0.24
献
[1] J.C. Maxwell, “A Dynamical Theory of the Electromagnetic
Field,” Philosophical Transactions of the Royal Society of
London, vol.155, pp.459-512, 1865.
[2] J.A. Stratton: Electromagnetic Theory, Mcgraw Hill, 1941.
[3] 平野拓一, 広川二郎, 安藤真, “1 次元問題を用いたモー
メント法、FDTD 法、有限要素法の電磁界解析教育用プ
ログラム,” 電子情報通信学会技術研究報告, vol. 112, no.
157, EST2012-31, pp. 143-148, 2012 年 7 月.
[4] R.F. Harrington: Field Computation by Moment Methods,
[5]
[6]
[7]
[8]
[9]
[10]
[11]
A.
Macmillan, 1968.
K.S. Yee, “Numerical Solution of Initial Boundary Value
Problems Involving Maxwell’s Equation in Isotropic Media,”
IEEE Trans. Antennas Propagat., Vol.AP-14, No.3,
pp.302-307, April 1966.
V. Demir, E. Alkan, A.Z. Elsherbeni, E. Arvas, “An
Algorithm for Efficient Solution of Finite-Difference
Frequency-Domain (FDFD) Methods [EM Programmer's
Notebook],” IEEE Antennas and Propagation Magazine,
vol.51, no.6, pp.143,150, Dec. 2009
J.L. Volakis, A. Chatterjee, L.C. Kempel, Finite Element
Method Electromagnetics: Antennas, Microwave Circuits,
and Scattering Applications, Wiley, 1998.
K.S. Kunz and R.J. Luebbers: The Finite Difference Time
Domain Method for Electromagnetics, CRC Press, Tokyo,
1993.
A.
Taflove
and
S.C.
Hagness:
Computational
Electrodynamics: The Finite-Difference Time-Domain
Method, Third Edition, Artech House Publishers; 3 edition,
Tokyo, 2005.
橋本修, 阿部琢美: FDTD 時間領域差分法入門, 森北
出版, 1996.
宇野亨: 「FDTD 法による電磁界およびアンテナ解析」,
コロナ社, 1998.
付録
A.1. 面電流からの平面波の放射
E
E
x
y


 0
y z
図 11 面電流からの平面波の放射
図 11に示す面電流のベクトルポテンシャルを計
算する。 A 

V
Je  jkr
dvs で計算したくなるかも
4r
しれないが、問題によっては微分方程式（ヘルムホ
ルツの方程式） A  k A   J から計算した方が
楽であり、この問題もその一例である。
2
2
 2 A  k 2 A   J
 2 Az
 k 2 Az  J ( x)
2
x
斉次解
斉次方程式:
 2 Az
 k 2 Az  0
2
x
Az  e x と仮定して代入すると、


 k 2 e x  0
    jk
2
H
1

 A
 ˆ jk jkx
zˆ
 yA e
1 Az 

0   yˆ

 x  yˆB jk e  jkx
Az


xˆ
yˆ
  / x 0

0
0
1
E

1
j
1
( x  0)
H
xˆ
yˆ
 / x 0
j
0
Hy

 zˆA

 zˆB

( x  0)
k2
e jkx
j
k 2  jkx
e
j
zˆ
1 H y
0  zˆ
j x
0
( x  0)
 zˆAje jkx

 zˆBje  jkx
( x  0) 
( x  0)
( x  0)
を決定する。
H
H
x  0 にしか波源が無いので、
 Ae jkx ( x  0)
Az    jkx
( x  0)
Be
x  0 での E, H の接線成分の境界条件から A, B
z
J  zˆJ ( x)
jkx
 jkx
斉次解: Az  Ae  Be
 xˆ  (E   E  )  0

 xˆ  (H   H  )  zˆJ
 yˆ j ( B  A)  0
J

⇒ A Bj
 jk
2k
 zˆ  ( B  A)  zˆJ

従って、
J

 yˆ e jkx ( x  0)

2
, H
J  jkx
 yˆ e
( x  0)
( x  0)
 2
A.2. 解析プログラム
Mathematica Ver.8 による解析プログラム例を示す。
A.2.1. モーメント法

 zˆ
E
 zˆ

J jkx
e
2
J  jkx
e
2
( x  0)
A.2.3. FDFD 法
A.2.2. FDTD 法
A.2.4. 有限要素法(FEM)

Download Report