[共同研究成果]
スーパーコンピュータを用いたトランジスタのスイッチング動作シミュレーション
春日貴志

秋田大学

1.

井上浩

工学資源学部

はじめにDC-DC コンバータやスイッチング電源などの電子機器は、負荷に応じて内部回路のスイッチング周波数を変化させて、高効率を実現する機器である。これら電子機器内のスイッチング回路を搭載したプリント回路基板（Printed Circuit Board: PCB）から放射する電磁ノイズが電子機器周辺で問題となっている[1]。しかし、低周波での電磁ノイズに関する具体的な対策が行われていないのが現状である。従来はノイズ放射を実験によって測定し対策する方法が主流であったが[2]、近年では電磁界数値シミュレーションによってノイズ放射を推定する方法も行われている[3]。
我々は、非線形で動作するスイッチング回路を搭載した PCB からのノイズ放射を解明することを目的として、有限差分時間領域法[4]（Finite Difference Time Domainmethod: FDTD method）を用いた数値シミュレーションによって、近傍から遠方までのノイズの放射伝搬を推定することを試みている[5, 6]。また、同時に実験結果との比較を行い、シミュレーションの妥当性も確かめながら研究の発展を計っている。
FDTD 法は、マクスウェルの電磁界方程式を差分化し、空間と時間領域で数値解析するシミュレーション法である[4]。この方法はモデル化が容易であることや、時間領域すなわち波形での解析が出来ることなどから、近年アンテナや電磁界の解析に良く用いられる解析手法である[7]。微分方程式を数値計算する方法であるため、解析空間内のセル数を細かくとれば計算精度の向上が望まれるものの、解析量が膨大となるため一般のコンピュータでの計算は困難である。精度良く広範囲の計算を行うためには、スーパーコンピュータのような大規模な計算機システムを利用することが必要である。
これまで、東北大学情報シナジーセンターのスーパーコンピュータ（SX-4）を利用して、高速計算に適した FDTD シミュレーションのアルゴリズムについて検討してきた[8]。また、トランジスタを等価回路によりモデル化し、FDTD 法に組み込んだ解析法の開発を目指してきた[9]。本稿では PCB モデルからの電磁放射解析のための高速化と、非線形解析を FDTD 法へ組み込みを行った例の結果を述べる。

2.
2.1.

FDTD シミュレーションの高速化
FDTD 法[4,7]

FDTD 法はマクスウェルの電磁界方程式を差分化することによって数値計算する方法である。マクスウェルの方程式を差分化すると、空間的な電界と磁界の配置は
図 1 に示す単位セルの配置となる。単位セルの配置に従い、z 成分について時間と空間で差分化した式を式(1a)、(1b)に示す。
n−n−∂Η x ⎞⎟∆t ⎛⎜ ∂H yεnn −1−Ez =Ez +∂y ⎟ε + σ∆t
ε + σ∆t ⎜ ∂x
12⎝Hn + 12z=Hn − 12z

12(1a)⎠n∆t ⎛⎜ ∂E y ∂E xn ⎞⎟−−∂y ⎟⎠µ ⎜⎝ ∂x(1b)x、y 成分についても同様に差分化することができ、3 次元解析空間の場合、図 1 の単位セルを積み重ねた解析空間において、電界と磁界の 6 つの式によって数値計算が行われる。FDTD 法では、解析モデルの電気的な特性であるε、σ、µをマクスウェルの電磁界方程式に与えることによってモデル化される。
図 2 に PCB とその周辺の電磁界解析モデルを示す[5]。解析空間内の数値は単位セル数である。3 次元 FDTD 法で解析する場合、電界と磁界の数式は x、y、z 座標による 3 次元ループとなる。図 2 上図の直方体（外側）は解析空間であり、内部にハッチングした PCB モデルが置かれている。図中の数値は、上右側の座標系を用いたときのセルの位置を示している。図 2 下図は PCB モデルの拡大図である。信号は、ストリップラインの信号源側の抵抗 50Ωを介して印加される。また、負荷側では 1kΩの抵抗を介して PCB の裏のアース面に接続される。
図 2 のように解析空間が非常に広い場合、3 重ループの計算量が膨大となるため汎用コンピュータでは計算が困難となる。また、解析モデルの電気的特性に応じて電界と磁界のパラメータが異なるため、計算時にはそれぞれの媒質毎の電界もしく

図1単位セル

図2

解析例[5]

は磁界計算を行う必要がある。このため、大規模計算システムすなわちスーパーコンピュータを利用する計算が必要となる。この際、スーパーコンピュータの高速化に適したアルゴリズムを適用することが重要である[8]。

2.2.

高速化の手法[10]SX-4 はベクトル・並列型のスーパーコンピュータであり、1 ノード当たりの理論最大演算性能は 64 GFLOPS、メモリは 8 GB、CPU は 32 個あり、最大 32 並列までの計算が可能であった。ベクトル化は最大 256 ループ計算を 1 度に処理する方法であり、並列化はループ計算を複数の CPU に割り当てる方法である。コンパイラには
do k=2,NZ1
do j=2,NY1
do i=1,NX1
if (idone(i,j,k)==0) then! 空間ex(i,j,k)=ex(i,j,k)+∆hz/∆y-∆hy/∆zelse
if (idone(i,j,k)==2) then! 誘電体
ex(i,j,k)=ex(i,j,k)+ ∆hz/∆y-∆hy/∆zelse
if (idone(i,j,k)==1) then! 導体
ex(i,j,k)=0.0;else! 抵抗ex(i,j,k)=ex(i,j,k)+ ∆hz/∆y-∆hy/∆z
end ifend ifend ifend doend doend do(a)IF を用いたアルゴリズム
図3

do k=2,NZ1
do j=2,NY1
do i=1,NX1! 空間
ex1(i,j,k)=ex(i,j,k)+ ∆hz/∆y-∆hy/∆z
end doend doend dodo k=KA1,KA2
do j=2,NY1
do i=1,NX1! 誘電体
ex1(i,j,k)=ex(i,j,k)+ ∆hz/∆y-∆hy/∆z
end doend doend do！導体………
！抵抗………
！電界の更新do k=2,NZ1
do j=2,NY1
do i=1,NX1ex(i,j,k)=ex1(i,j,k)
end doend doend do(b)

IF を用いないアルゴリズム

高速化アルゴリズムの工夫例

C 言語と FORTRAN90 が用意されている。高速化を図るため、(1)プログラミング言語による計算速度の違い(2)ベクトル化・並列化に適したプログラミング(3)ノイズスペクトルの推定に必要な最短信号長の 3 点について検討した。(2)の具体的な対策としては、ベクトルループ中に IF-ELSE文があると計算効率が悪くなることから、図 3 に示すような IF-ELSE 文を用いないプログラミングを行った。さらにループ長の最も長いループにはベクトル化命令を、次に長いループには並列化命令を指定した。(3)については、繰り返し連続波形を想定して、スペクトルが推定可能な最短データ長について検討した。
本解析に用いたPCBモデルを図 4 に示す。基板はFR4（εr=4.9）を用い、裏面は全面銅箔と設定した。線路終端は 1kΩで終端した。実験およびシミュレーション共に、入力信号は小型発振器からの 3 MHzの方形波を用いた。FDTD法における単位セル
サイズは∆x=2.5mm、∆y=5mm、∆z=0.35mmとし、時間ステップ∆tは、クーラントの条件を満たすように、1.15psに設定した。3MHzの 3 周期に必要な計算時間ステップ数Nは 8.7×105回である。吸収境界条件には 2 次のMurを用いた。解析空間のセル数と計算時間の関係を考察するため、2 種類の解析空間(100×100×200cell、300×300×400cell)について計算した。また、微小磁界プローブ（Shielded Loop Probe）を用いて線路上の近傍磁界を測定し、実験と計算結果とを比較した[13]。

図42.3.

高速化シミュレーションに使用した PCB モデル

最適な FDTD 計算条件の検討プログラミング言語の違いによる計算時間について検討した。C 言語とFORTRAN90 による計算プログラムは、それぞれシングル、8、16、32 並列で計算し、並列 CPU の違いについても比較した。また、時間ステップ数 N は 1,000 回とした。シングルは User Time、並列は計算で使用した CPU のうち最も利用時間の長いCPU 時間を表 1 にまとめた。100×100×200cell における解析では、8 並列の結果が最も良く、C 言語よりも FORTRAN90 の方が 6.3 倍高速計算可能であった。16、32並列の計算に時間がかかったのは、オーバーヘッドによる影響と考えられる。
また、ベクトル、並列化に適したプログラミングの改善前後の比較、および解析空間の違いによる比較も表 1 に示した。8 並列で比較すると、改善によって約 1.3倍の高速化が可能であった。IF-ELSE 文の除去により、計算回数が増加するものの、ベクトルループ中の計算効率が高まったため、計算時間の短縮につながったと考えられる。300×300×400cell の解析では、100×100×200cell よりも並列ループにかかる計算量が増加したため、16 並列で最も効率が良くなったものと考えられる。解析量の異なる二つのモデルについて、時間ステップ数に対する CPU 時間の関係を図 5に示す。二つの解析モデルは、最適な計算条件によって解析を行っている。時間ステップ数が増加するにつれて、CPU 時間も比例して増加する。解析空間の計算量によって比例係数は異なるが、時間ステップ数と CPU 時間に比例関係が成り立つため、少ない時間ステップ数から最終的な計算に要する時間が予測できる。あらかじめ少ない時間ステップ数によって最適な並列数を探ることは、計算効率の点で重要であると考えられる。
表1
CPUCPU 時間の比較（NT=1000 ステップ）

C language

Fortran 90

100×100×200cell

100×100×200 cell
Unimproved

Improved

300×300×400
cell

Single

410 s

158 s

83 s

1556 s8 parallel

196 s

31 s

23 s

217 s

16 parallel

256 s

92 s

44 s

178 s

32 parallel

384 s

148 s

137 s188 s

図5

時間ステップ数に対する計算時間

2.4.

波形計算例

給電点から 1 cm、線路上 3 mm における近傍磁界の 3 周期分の計算波形を図 6(a)に示す。波形の 3 周期と 1 周期をそれぞれ高速フーリエ変換（FFT）し、それぞれのスペクトルのエンベロープを図 6(b)に示す。1 周期と 3 周期の結果は一致した。
また、実験結果ともほぼ一致することから、1 周期の波形を計算することによってスペクトルが推定できると考えられる。FDTD 法で 1 周期の波形を求めるために要する時間は約 10 時間であることから、数 MHz 帯の低周波ノイズの FDTD シミュレーションが 6 dB 程度の工学的な誤差範囲内で十分可能であると考えられる。

(a)

磁界波形

(b)
図63.
3.1.

磁界スペクトル近傍磁界ノイズスイッチングトランジスタのシミュレーション[9]トランジスタの大振幅動作のモデル化

FDTD法を用いた解析の応用例として、スイッチングトランジスタを搭載したPCBモデルからの電磁波放射について解析を試みた[9]。トランジスタのFDTDモデリン
グとしては、エバース・モルの等価回路を用い[4]、印加電圧によって変化する拡散容量と接合容量を付加したモデル[11]を元に、電荷蓄積効果を含む低周波スイッチング解析のFDTDアルゴリズムを提案した[9]。図 7 にトランジスタの等価回路と、スピードアップコンデンサCB によって高周波特性を改善したインバータ型スイッチング回路を示す。C'-B間は拡散容量CDC、接合容量CJC、電流源とダイオードからなり、各成分を等価的に 1 つの素子と見なすことが出来る。B-E間は順バイアス状態にあるため各容量並びに電流源を省略した。C-C'間にスイッチング特性に影響の大きいコレクタ寄生抵抗rC を付加した[12]。C'-B間、B-E間に流れる電流IC'B 、IBE 、およびC'-B間容量CC'Bは次式で表される。

⎧− qVC ′B⎞ − 1⎫ − α I ⎧exp⎛ qV BE⎞ ⎫ + C dVC ′B
I C ′B = I C 0 ⎨exp⎛⎜⎜⎟⎬⎟ ⎬F E0 ⎨C ′BkTkT⎠⎭⎠ ⎭dt⎩ ⎝⎩ ⎝(2a)⎧qV⎞ − 1⎫
I BE = − I E 0 ⎨exp⎛⎜ BE⎟ ⎬kT⎠ ⎭⎩ ⎝(2b)C’-B間の容量CC'Bは接合容量CJCと拡散容量CDCの合成容量として与えられ、次式で表される。

⎛ V ⎞CC′B = C JC (0)⎜⎜1 + C′B ⎟⎟Φ0 ⎠⎝−m+qτ DC I C 0 exp⎛⎜ − qVCB kT ⎞⎟kT⎝⎠(3)

FDTD 法で回路素子を計算する場合、式(4)に示すように差分化式に素子の電流項を加えることで計算が可能である。
n +1nEx i , j , k = Ex i , j , k +∆tε0n+ 1∇ × H i , j , k2 −∆tn+ 1
I x i , j , k2
ε 0 ∆y ∆z(4)トランジスタのC’-B間、B-E間の電界は、式(4)のIxに式(2)のIC'BとIBEを代入することで計算が出来る。

図73.2.

Ebers-Moll のトランジスタ等価回路

シミュレーション例

スイッチング用トランジスタ(2SC2671)の物理パラメータ値例を表 2 に示す。I0、αF、Φ0、CJC(0)、mは測定によって求めた。拡散容量は測定することが難しいため、未知の値である τ DC はFDTDシミュレーションによって求められる出力電圧VCE の波形が実験値と一致する値を用いた。3 MHz、振幅 5Vppの方形波を印加した時の出力電圧VCEのスペクトルを図 8 に示す。RC=470Ω、RB=3.3kΩ、CB=47pFとし、図 7 の信号線路と終端負荷は開放とした。単位セルサイズは∆x=1mm、∆y=2.5mm、∆z=0.255mmとし、時間ステップは 0.82psである。FDTDとSPICEによるシミュレーション結果はよく一致し、実験値とも最大で 4 dBの誤差内で一致した。

表2トランジスタ 2SC2671 の物理パラメータ[9]ParametersValues-9

逆方向飽和電流 I 01x10 A

電流増幅率 α F0.984

拡散電位 Φ 00.8 V

動作温度 T

300 K

ボルツマン定数 k-23

1.38x10J/K

接合容量 C JC (0)3.5 pF

接合容量の乗数 m0.15

少数キャリアのベース通過時間 τ DC
コレクタ内部抵抗 r C [12]図8180.0 ns
10 Ω

トランジスタの出力電圧VCEの周波数特性[9]図 9 にスイッチングトランジスタを搭載したPCBモデルを示す。基板の裏面は全面銅箔とし、基板の比誘電率εrは 3.4 であった。信号線路の特性インピーダンスはTDR測定により 38Ωであった。信号線路の入力端はスイッチングトランジスタ(2SC2671)を接続し、終端には不整合したモデルを想定して、1kΩチップ抵抗を使用した。裏面に装着した小型発振器によって方形波（3 MHz、5Vpp）を印加した。発振器並びに電源用小型リチウム電池は銅ケースで覆った。PCB近傍磁界Hxはシールデットループプローブを用いて測定した[13]。
図 10 に給電点から 1 cmにおける信号線路上の近傍磁界Hxの波形と、その周波数スペクトルをそれぞれ示す。計算値は図 10(a)の波形をFFT計算して、周波数スペクトルを求めた。0.16µsでトランジスタはOFFからONに遷移し、これに伴う大きな磁界が発生する。図 10(b)の実験値とFDTD値 は周波数の広い帯域内で一致して

図9(a)
図 10

トランジスタを装荷した PCB モデル

磁界波形

(b)

磁界スペクトル

信号線路上の磁界波形とスペクトル（給電点より 1cm）[9]おり、ノイズ放射の FDTD 解析が可能であることが示されている。
また、この解析では最速条件である 16 並列で解析を行い、解析に使用したメモリ量が 1.2GB、CPU 時間が約 17 時間であった。計算の出力データ量によっては、計算時間が大きく増加する場合もあるが、およそ 17 時間で解析が可能であるため、スーパーコンピュータによる解析が有益であると考えられる。

4.

おわりにPCB からのノイズ放射の FDTD シミュレーションについてスーパーコンピュータを用いて実行するため、高速化アルゴリズムの検討、並びに非線形スイッチングトランジスタを組み込んだ計算アルゴリズムの開発を行った。高速アルゴリズムの検討により、低周波で駆動した PCB からのノイズ放射の FDTD シミュレーションが可能となった。また、トランジスタを組み込んだ FDTD 法による PCB からのノイズ放射の解析では、実験と良く一致するシミュレーションが可能であった。スーパーコンピュータを使用した解析術が一歩進歩したと考えられる一方、スーパーコンピュータもメモリサイズの大容量化・計算速度の高速化が進んでいるので、今後本稿で提案したシミュレーション技術を応用したノイズ解析およびその応用研究が望まれている。

謝辞
本研究は東北大学情報シナジーセンターとの共同研究、｢広帯域・広範囲電磁界解析用高速・高精度計算機シミュレーション法におけるスーパーコンピュータの有効利用に関する研究｣として行われた。

参考文献[1] 電気学会，電子機器のノイズアイソレーション技術，コロナ社，(1998)[2] 木下敏雄，EMC の基礎と実践−電磁障害とノイズ対策−，日本工業新聞社[3] 山下榮吉，電磁波問題解析の実際，電子情報通信学会，(1993)[4] A. Taflove，Computational Electrodynamics – The Finite Difference Time Domain
Method, Boston, Artech house Publishers, (1995)[5] T. Kasuga, M. Tanaka, and H. Inoue, "FDTD Simulation and Experimental Study onLine Impedance and Magnetic Near Field Noise for a Simple Printed Line Model",IEICE Trans. Commun., Vol. E83-B, No. 3, pp. 561-568, (2000-03).
[6] T. Kasuga, M. Tanaka, and H. Inoue, "Estimation of Spatial Distribution of Wideband
Electromagnetic Noise around a Printed Circuit Board", IEICE Trans. Commun.,Vol. E86-B, No. 7, pp. 2154-2161, (2003-07).
[7] 宇野亨，FDTD 法による電磁界およびアンテナ解析，コロナ社，1998.
[8] 春日貴志，井上浩，”プリント基板からの電磁波雑音放射の FDTD 法による一シミュレーション”，信学技報，EMD98-85，(1999-01)[9] 春日貴志，田中元志，井上浩，”スイッチングトランジスタを搭載した PCB モ
デルの FDTD モデリング”，信学論(C)，Vol. J85-C，No. 4，pp.304-305，（2002-04）.
[10] 東北大学大型計算機センター広報 SENAC 編集委員会，“スーパーコンピュータ特集号”，SENAC，Vol. 30，No. 4, (1997-11)[11] P. Ciampolini, P. Mezzanotte, L. Roselli, and R. Sorrentino, “Accurate and Efficient
Circuit Simulation with Lumped-Element FDTD Technique”, IEEE Trans. MTT,vol.44, no.12, pp.2207-2215, (1996-12).
[12] 玉井德迪，半導体回路設計技術，日経 BP 社, (1987).
[13] 滝田栄志，田中元志，井上浩，“近傍磁界プローブの校正と PCB 近傍におかれた平面導体における磁界分布の変化計測”，信学技報，EMCJ2000-39，(2000-07).