スーパーコンピュータを利用した大規模電磁界解析法の提案
春日貴志

秋田大学

1井上浩

工学資源学部

はじめに電子機器から放射される電磁ノイズは、システムの高速化や様々な高性能回路の利用によって、近年非常に増加している。電子機器からの電磁ノイズ放射を推定する方法としては、測定とシミュレーション法による解析がある。測定法の一つとして、産業・民生機器から放射する電磁ノイズ量の規格に対する適合性を認証する方法が確立されている。この方法は、放射源から 3m、10m、30mでの電界を測定し、CISPR などで定めた妨害放射強度量の許容値と比較する方法である [1, 2, 3]。また、電子機器内の電磁ノイズ放射源を特定する方法としては、近傍磁界を測定する方法がある。
一方、シミュレーション法は、産業用や民生用機器の設計支援として用いられる電磁ノイズ解析や、電磁ノイズ放射や対策の基礎研究用に幅広く用いられている。産業用シミュレーション法の利点としては、予めシステム設計の段階から電磁ノイズ放射量を予測することで、製品を試作すること無しに EMI 対策ができるため、システム設計時間と費用を短縮する利点がある。また基礎研究用としては、測定が難しい場所での電磁ノイズを予測するための手法や、電磁ノイズ放射の理論解析の手法として、EMI 対策の理論的背景を築くことに用いられる。
しかしながら、従来から行われている測定とシミュレーション法は、近傍もしくは 3m 以上の遠方界における電磁界を求める方法であり、電子機器から 3mまでの電磁界を推定する方法は無い。このため、電子機器周辺での電磁妨害に対する電磁ノイズ推定法や対策法は確立されていない。
筆者らは、これまでスーパーコンピュータ（東北大学情報シナジーセンター及び秋田大学総合情報処理センター）を用いて、FDTD 解析の高速化 [4] と広範囲の電磁界解析法 [5] について検討を行ってきた。
本研究では、これまで電磁界シミュレーション法として用いられてきた有限差分時間領域法（Finite Diﬀerence Time Domain method: FDTD 法）を応用することで、従来よりも広範囲の電磁界解析が可能な FDTD-Multiple Analysis

Space 法（FDTD-MAS 法）を開発した [6]。従来の FDTD 法によって、例えば電子機器内のプリント回路基板（Printed Circuit Board: PCB）からの電磁放射を解析する場合、解析における単位セルサイズは基板の構造（基板厚、線路幅）によって決定する必要があり、単位セルサイズは非常に微小となる。このため、PCB からの電磁放射を広範囲で求める場合には、計算量が膨大になってしまう。著者らが開発した FDTD-MAS 法は、PCB と極近傍の電磁界解析はPCB の構造によって決定したセルサイズによって構成した解析空間（InternalAnalysis Space: IAS）で解析し、IAS よりも遠方は IAS のセルサイズの整数倍のセルサイズによって構成した解析空間（External Analysis Space: EAS）によって解析を行う方法である。この手法により、空間のセル数と時間ステップ数を大幅に減らすことが出来るため、これまでよりも広範囲かつ広帯域の大規模解析が実現できる可能性が拡大した。
本稿では、FDTD-MAS 法の原理と、ダイポールアンテナを放射源とした時の 6 × 6 m2 における電磁界分布を、スーパーコンピュータ SX-7 を用いて計算した例を述べ、大規模解析の実証例を示す。個々で示す解析方法は、SX-7 のような高速かつ大容量のスーパーコンピュータを用いることにより実現可能なものであり、その有用性を示すことが出来る。

2FDTD-MAS 法
FDTD-MAS 法の特徴は、全く独立した複数の解析空間を組み合わせることで、広範囲の電磁界解析を実現できることである。例えば PCB からの電磁放射を解析する場合、PCB とその周辺の空間における電磁界を細かいセルサイズにより構成した解析空間で解析し、さらに遠方での電磁界を粗いセルサイズで構成した解析空間で解析する。図 1 に FDTD-MAS 法の解析モデル図の概念を示す。内部の細かいセルの解析空間を Internal Analysis Space (IAS)、外部の粗いセルの解析空間を External Analysis Space (EAS) と名付けた。内部にダイポールがおかれている例である。IAS と EAS の単位セルサイズと時間ステップ数はそれぞれ独立して与える。IAS と EAS のセルサイズをそれぞれ、∆IAS と∆EAS としたとき、セルサイズ比 Ra は式 (1) となる。
Ra = ∆EAS /∆IAS(1)

IAS の時間ステップ ∆tIAS と EAS の時間ステップ ∆tEAS の比も、式 (2) となる。
Ra = ∆tEAS /∆tIAS(2)

EAS の時間ステップ数は、IAS よりも 1/Ra に減らすことが出来る。IAS と EASの解析データを連結させる方法として、IAS 内にある放射源を取り囲むデータ

(iI2,jI2,kI2)

S
θ
(IA,JA,KA)

ϕr

∆IAS

(iI1,jI1,kI1)-0.15 0.15m
16 layer PML

Szyx-0.15㨪0.15m-3.5㨪3.5m(a)
図1∆EAS=Ra∆IAS(b)

FDTD-MAS 法の解析空間の基本概念図、(a) 内部解析空間 (IAS)、(b) 外部解析空間 (EAS)

変換面上の電磁界を平均し、EAS へ与える。図 2 は IAS から EAS へのデータ変換面（x − z 平面）におけるの詳細図である。IAS の計算データは EAS のセルサイズ ∆xIAS と ∆zEAS の面積の範囲内で平均化される。すなわち、セルサイズ比で平均化面積を表すと −Ra/2 から Ra/2 (もしくは Ra/2 − 1) の範囲となる。EAS における x − z データ変換面上の電界 Ex と磁界 Hx は式 (3)、(4) となる。

ExEAS (iE + 1/2, j, kE ) =1(Rax − 1)Raz·x −1
iE + Ra2

zkE + Ra2

x k =j − Raz
iI =IE − RaIE22HxEAS (iE , j + 1/2, kE + 1/2) =·x
iE + Ra2

z −1kE + Ra2x k =k − Raz
iI =IE − RaIE
22HxIAS (iI , j +ExIAS (iI + Ra2 x , j, kI ) (3)1Rax (Raz − 1)
RayRaz2 , kI + 2 )(4)

ここで、ExEAS と HxEAS はそれぞれ EAS における電界と磁界、ExIAS と HxIASはそれぞれ IAS における電界と磁界である。x − z 面であるため、式 (3)、(4)の j は j = je1 もしくは j = je2 となる. iE は EAS、iI は IAS の座標である. 式(3)、(4) と同様に y 、z 成分や x − y 、y − z 平面も同様に定式化が可能である。

HzEAS(ie2,je1+1,ke2)
EyEAS
ExEAS
HyEAS

HxEAS
EzEASS(ie1,je1,ke1)∆ze∆xeexternal analysis space
ExIAS∆ziz ∆xiyinternal analysis spacex

図2

IAS から EAS へのデータ変換面の詳細

図 3 に、Ra=6 の時の IAS と EAS における電界と磁界のタイミングチャー
トを示す。IAS と EAS の時間ステップは異なるため、IAS と EAS の解析には、時間的同期を取る必要がある。図 3 より、IAS と EAS の電界の時間間隔は、式(5) で一致し、このとき IAS から EAS への電界値を変換する。
(nEAS )∆tEAS = (nIAS )∆tIAS ,(5)ここで、nEAS と nIAS はそれぞれ EAS と IAS の時間ステップ数である。
一方、磁界の保存タイミングは、Ra が奇数と偶数の時で異なる。Ra が奇数の場合、IAS と EAS の磁界の時間間隔は(nEAS + 1/2)∆tEAS = (nIAS + 1/2)∆tIAS(6)の時保存される。しかし、Ra が偶数の場合、IAS の磁界の時間が (nIAS +1/2)∆tIASの時、EAS における同じ時間は電界の解析となるため、磁界に関する式 (6) の条件を満たさない。そこで、(n − 1/2)∆tIAS と (n + 1/2)∆tIAS の磁界の平均値を与える。

EIAS∆tIASt

HIAS
EEAS

∆tEAS

t=n∆tEAS
HEAS

図3t=(n+1)∆tEAS
t=(n+1/2)∆tEASt

IAS と EAS における電磁界のタイミングチャート（Ra=6）

Initialize andsetup parameter

Initialize andsetup parameter
Calculation ofmagnetic fieldCalcultion ofelectric fieldt=(n+1/2)∆tEAS

Excitation ofelectric fieldt=n∆tEAS

yes
outputmagnetic fieldnoCalculation ofmagnetic fieldCalcultion ofelectric field

Input
magnetic fieldyes
Input electricfield

H-PML
E-PML
t=t+∆tIAS/2

H-PML

yes
outputelectric fieldt=t+∆tIAS/2

t=t+∆tIAS/2

E-PML

t=t+∆tIAS/2

IAS
図4EAS

FDTD-MAS 法のフローチャート

図 4 に FDTD-MAS 法のフローチャートを示す。左が IAS、右が EAS である。

IAS と EAS の解析を別々に行う場合、IAS の変換データは一度ハードディスクに保存され、EAS の解析時にハードディスクから読み出す。一方、コンピュー
タメモリを大量に保有しているコンピュータを用いる場合、IAS と EAS の解析を同時に行うことが出来る。スーパーコンピュータ SX-7 は 256GB のメモリを有するコンピュータであるため、大規模な計算を行うことが可能である。このため、ハードディスクとの I/O にかかる時間が必要ないため、計算の高速化が実現できる。

3ダイポールアンテナ放射源による FDTD-MAS 法の評価
FDTD-MAS 法による広範囲の電磁界解析の可能性を検討するため、1 GHzの半波長ダイポールアンテナを放射源として、3 m までの電界分布を計算して理論式と比較する。半波長ダイポールは、アンテナ構造が簡単である上、放射する電磁界理論がすでに確立されている。すなわち、ダイポールの理論式により求めた電界分布と比較することで、新しいシミュレーション法を評価が可能であると考えられ、半波長ダイポールアンテナを用いた。図 1 では、1 GHz の半波長ダイポールアンテナを解析空間内に含んで示した。ダイポールアンテナ

Ra

表1計算条件

∆IAS

∆EAS

∆tIAS

∆tEAS

I0mmmmps

psmA

62.5000

1.443

27.3781.8750

1.604

28.21101.5000

1.804

28.90

121.2500

141.0714

2.405

29.35

160.9375

2.886

29.48

180.8333

3.608

29.59

200.7500

4.811

29.62

152.062

28.86

29.12の励振源を IAS の中心におき、座標を (x, y, z) = (IA, JA, KA) と定義する。ダイポールアンテナの入力抵抗は 73.13 Ω とした [7]。ダイポールアンテナ長 L を考慮した電界分布と指向性パターンの理論式は式 (10) で与えられる [8]。

Ey (r)|ϕ=90degree =−jI0 Z22π r 2 + L4exp

− jβr2 +j60I0 ejω{t−(r/c)} cos {(βL cos θ)/2}Eθ =rsin θL24

(10a)

(10b)

ここで、r はダイポールアンテナからの距離、L はダイポール長、Z は自由空間における固有インピーダンス、I0 は時間における電流の最大値である。式 (10a)中の ϕ は 90 度と仮定し、式 (10b) 中の r はアンテナ中心からの距離である。

IAS と EAS のセルサイズ比を 6 から 20 まで変えた時の電界分布を計算した。
EAS のセルサイズは 1 GHz の波長の 1/20 にするため、15 mm とし、セルサイズ比に応じて IAS のセルサイズを変えた。各セルサイズ比における計算条件を表 1 に示す。ダイポール長 L は 0.15 m で、アンテナ半径 φ は 0.15 mm とした。
アンテナ半径は IAS のセルサイズよりも小さいため、サブセル法 [9, 10] を適用しアンテナ半径を考慮した解析を行った。吸収境界条件としては Mur よりも高精度な Perfect Matched Layer (PML) を用いた。PML の層数は 16 層、反射係数は 1 × 10−8 、導電率の乗数を 3.1 とした。
図 5 に Ra = 10 の時のダイポールアンテナ上の電流分布の計算結果を示す。
ダイポール上の電流分布は時間によって変化するため、一定時間における電流

Current [mA]

3020100-0.1
図50Distance [m]0.1

ダイポールアンテナ上の電流分布（Ra = 10）の最大値をプロットした。典型的な半波長の分布となっている。各セルサイズ比における最大電流値 I0 を表 1 に併せて示す。ダイポール半径を 0.15 mm と一定にしたものの、IAS のセルサイズを変えたため、セルサイズ比により電流の最大値も変化する。Hockanson や渡辺らはセルサイズに対するダイポールアンテナの入力インピーダンスの変化について検討している [11, 12]。最も影響が大きいのは、セルサイズを変化させたことによってダイポールアンテナの給電点のギャップ長も変化することである。しかし、今回の検討は IAS から EAS へのデータの連続性と広範囲の電磁界解析の可能性を検討することであるので、セルサイズによる電流もしくは放射電界の変化分については、Ra = 10 の最大電流値で正規化することにする。
図 6 に Ra = 20 における x − y 平面の電界 |E| 分布を示す。図 6(a) は 6 × 6 m2の電界分布、(b) は IAS で計算した 0.3 × 0.3 m2 の電界分布を示している。図6(a) 中央部の破線のかこみは IAS から EAS へのデータ変換面を示しており、破線内部は IAS の解析領域である。ダイポールアンテナは y 方向に配置されている。ダイポールアンテナの直上方向の電界は非常に弱く、典型的な半波長ダイポールからの放射パターンとなっている。図 7(a)〜(c) はダイポールアンテ
ナからの電界分布（0.1〜3m）と、|Eθ |、|Ey (ϕ)|r=3 m の指向性特性である。指向性はダイポールアンテナの励振源から 3 m における電界強度分布である。実線は FDTD-MAS 法による計算結果であり、|Eθ | と |Ey (ϕ)|r=3 m はそれぞれ式(11a)、(11b) によって算出した。
|Eθ (x, y, KA)| =Ex2 (x, y, KA) + Ey2 (x, y, KA)|Ey (ϕ)| = |Ey (x, JA, z)|.

(11a)
(11b)

3
180

2y [m]

IAS

140

0120-1|E| [dBµV/m]

160

1100-280-3-3-2-10x [m]123(a)
0.15
1800.1160
1400120-0.05|E| [dBµV/m]

y [m]

0.05

100-0.1
80-0.15
-0.15-0.1-0.05

0 0.05x [m]0.10.15(b)図6ダイポールからの電界分布の解析結果, (a) 6 × 6 m2 範囲の分布, (b) IAS 内の分布

破線は式 (10) による理論式より求めた値である。図 7(a) において、横軸は波数 1/k(= λ/2π) により正規化した。FDTD-MAS 法によって計算した電界分布において、Ra を 6 から 20 まで変化させても電界分布の差は非常に小さい。また、IAS から EAS への境界面での分布は連続的に求められており、広範囲で連続的な電界分布の計算が行われている。0 度と 180 度方向の |Eθ | は電界が弱く、|Ey (ϕ)|r=3 m は全ての角度において円状分布となっている。距離特性と指向性特性の FDTD-MAS 法の計算結果と理論値は 1 dB 以内で一致することから、FDTD-MAS 法を用いて広い空間範囲内の電磁界解析が可能であることが示された。
図 7(a) に全解析空間を EAS のセルサイズ 15 mm で構成した解析モデル(Coarse grid) で計算した結果をあわせて示す。ダイポールアンテナの最大電Electric field strength[dBµV/m]150
IAS EAS
Line: FDTD (Ra=6 to 20)

140

Chain: Coarse Grid

130

Broken: Theoretical Value

12011015 1050
Normalized Distance (2πr/λ)
(a)

0150

1001005090180(b)図715050270

FDTD(Ra=6 to 20)
Theoretical value
Coarse grid

090270

FDTD(Ra=6 to 20)
Theoretical value
Coarse grid180(c)

電界の距離特性と指向性特性, (a) |Ey (r)|ϕ=90degree の距離特性、(b) |Eθ | の指向性、(c) |Ey (ϕ)|r=3 m の指向性

流値 I0 は 18.8 mA であったので、Ra = 10 の電流値との差分を求め、差分に比例して電界値を正規化した。粗いセルサイズによって計算した電界値と理論値との差は 4 dB であった。粗いセルで解析空間をモデル化することによって、ダイポールを構成するセル数も大幅に減少する。このため、放射電界の誤差が大きくなったと考えられる。FDTD-MAS 法の計算精度は、セルサイズを粗くした計算よりも高精度であることが明らかとなった。

4まとめ
広範囲の電磁界解析を目的とした FDTD-MAS 法を開発し、ダイポールアンテナを用いて有効性の検討を行った。その結果、3m までの電界分布を連続的に解析できることができた。また、粗いセルサイズで構成した場合よりも計算精度が良いことが明らかとなった。東北大学情報シナジーセンターのスーパーコ
ンピュータ SX-7 は 256GB のメモリを有するために、広範囲の大規模解析には有効であると考えられる。今後は、電磁界問題への適用について検討したい。

参考文献[1] 電気学会・電磁波の散乱・吸収計測技術調査専門委員会編，電磁波の散乱・吸収計測と建築電磁環境，コロナ社, 2000.

[2] 井上正弘，“家庭用電気機器および照明装置からの妨害とイミュニティ”,
EMC, no. 133, pp.38-43, May 1999.
[3] 大澤千春，井上博史，“パワーエレクトロニクス関連 EMC 規格の動向”,
EMC, no. 184, pp.28-36, Aug. 2003.
[4] 春日貴志，井上浩，“スーパーコンピュータを用いた PCB 周辺における電磁界分布のシミュレーション”，東北大学情報シナジーセンター大規模科学
計算システム広報，SENAC vol.37, no.3, pp.29-36, 2004.

[5] 春日貴志，井上浩，“スーパーコンピュータを用いたトランジスタのスイッチング動作シミュレーション”，東北大学情報シナジーセンター大規模科学
計算システム広報，SENAC vol.36, no.4, pp.21-30, 2003.

[6] T. Kasuga, and H. Inoue, “Novel FDTD Simulation Method using
Multiple-Analysis- Space for Electromagnetic Far Field”, IEEE Trans.
on Electromagnetic Compatibility, （印刷中）[7] J. D. Kraus, Antennas, McGraw-Hill Inc., 1988.
[8] 長谷部望, 電波工学, コロナ社, 1995.
[9] A. Taﬂove, Computational Electrodynamics –The Finite-Diﬀerence TimeDomain Method–, Artech House Publishers, Boston, 1995.
[10] 宇野享, FDTD 法による電磁界およびアンテナ解析, コロナ社, 1998.
[11] D. M. Hockanson, J. L. Drewniak, T. H. Hubing, and T. P. Van Doren,“FDTD modeling of common-mode radiation from cables”, IEEE trans.
on Electromagnetic Compatibility, vol. 38, pp.376-387, Aug. 1996.
[12] 渡辺聡一、多気昌生, “FDTD 法におけるアンテナ給電点のモデル化に関する検討”, 信学技報, EMCJ96-118, Mar. 1997.