文部科学省 先端研究基盤共用・プラットフォーム形成事業『みんなのスパコン』 TSUBAME による日本再生

東京工業大学 TSUBAME 産業利用
平成 26 年度利用終了課題 利用成果報告書

東京工業大学 学術国際情報センター
共同利用推進室
http://www.gsic.titech.ac.jp/kyoyo

平成 26 年度終了課題 利用成果報告書（利用区分別採択日順）
利用
区分

課題名

企業名

頁

古河電気工業株式会社
次世代電池研究開発センター

1

無機材料開発への第一原理計算の活用

株式会社 ニコン

5

電子写真システム設計における電磁場計算の高速化

株式会社 リコー

9株式会社 資生堂
生産技術開発センター

13

リチウムイオン二次電池正極材料の第一原理計算

ホモジナイザーを用いた撹拌混合による乳液製造のスケールアップに関する解析トライアルユース

LTE-Advanced における
株式会社パナソニックシステ
大型車両内電磁界特性に関する基礎検討 ムネットワークス開発研究所

17

大規模シミュレーションによるレーダの車両搭載時の特性把握

古河電気工業株式会社
自動車電装技術研究所

21

大容量の行動関連データを使った最適シミュレーション及び暗黙知の利用

株式会社
電通国際情報サービス

25

減衰を考慮した高周波数領域までの音響構造連成シミュレーション大規模化技術の検討

フォスター電機株式会社

29

清水建設株式会社
技術研究所

33

ルネサスシステムデザイン
株式会社

37

ワイヤレス電力伝送による株式会社パナソニックシステ
漏えい電波の環境解析技術の研究開発 ムネットワークス開発研究所

41

広域都市環境の大規模計算による検討
三次元電磁界シミュレータを用いた静電気放電イミュニティ試験に於けるPCB/Package/Chip のイミュニティ解析

利用
区分

課題名

有償利用 成果公開

企業名頁

気象イベントを考慮した建築環境の解析評価

清水建設株式会社
技術研究所

47

鋼中析出物の水素捕捉能の高精度計算

新日鐵住金株式会社51

拡張アンサンブルシミュレーションによるタンパク質とリガンドの結合構造予測の展開

武田薬品工業株式会社
医薬研究本部

55

電子デバイス材料の計算機設計

太陽誘電株式会社

59

理論計算に基づく有機半導体材料の開発

住友化学株式会社
先端材料探索研究所

63

TOTO 株式会社
技術開発センター

69

トヨタ自動車株式会社

77

衛生陶器混相流シミュレーションの商品設計および販促への展開

塗工スラリーの分子シミュレーション

ワイヤレス電力伝送による株式会社パナソニックシステ
漏えい電波の環境解析技術の研究開発 ムネットワークス開発研究所

81

大容量データ伝送用ミリ波アンテナのレドームに関する基礎検討87

スタッフ株式会社()

（様式 20）成果報告書

平成 26 年度 TSUBAME 産業利用トライアルユース 成果報告書

利用課題名 リチウムイオン二次電池正極材料の第一原理計算
英文： An ab initio study of Li-ion battery cathode material
風間 吉則、谷村 雄大、幡谷 耕二
Yoshinori Kazama, Yudai Tanimura, Koji Hataya
古河電気工業株式会社 横浜研究所 次世代電池研究開発センター
Battery Materials R&D Center, Yokohama R&D Laboratories, FURUKAWA ELECTRIC CO., LTD.
http://www.furukawa.co.jp/
リチウムイオン二次電池用の正極活物質としてケイ酸塩系正極活物質 Li2MSiO4(M: 遷移金属)が注目を集めている。Li2MSiO4 は組成式内に Li を 2 つ含んでおり、2 電子反応を起こす高容量材料として期待されている。理
論容量は 330 mAh/g であるが、Li2FeSiO4 は 1 電子反応までしか起こらず、Li2MnSiO4 はサイクル特性に課題があることが知られている。更なる特性向上を目指し、M や Si に代わる置換元素を選定するため、TSUBAME を用いた大規模計算の有用性を検討することを目的とした。本報告では LixMSiO4 (M = Mn、Fe、Co、Ni)について
格子定数の構造緩和と平均電位計算を行ったことについて報告する。
Silicate-based active materials Li2MSiO4 (M: transition metal) attracts attention as cathode active
materials for the lithium ion battery. Li2MSiO4 includes two Li per formula unit and is expected as highcapacity materials causing 2 electron reactions. The theoretical capacity is 330 mAh/g, however,
Li2FeSiO4 exhibits only 1 electron reaction, and Li2MnSiO4 has inferior cycling performance. Aiming atfurther characteristic improvement, we examined the usefulness of the large-scale calculation using
TSUBAME to choose a substituted element for M and Si. In this report, we performed the structure
optimization of the lattice constant and average electric potential calculation about LixMSiO4 (M = Mn,
Fe, Co, Ni).

Keywords: 第一原理計算、リチウムイオン二次電池、正極材料、Li2MSiO4、結晶構造

背景と目的

し、活物質として最適な組成や結晶構造を、実験的に低炭素社会の実現に向け、電気自動車や電力貯蔵

見出だしていくのは非常に時間を要する。そこで、計算用途への大容量リチウムイオン電池への要求が高まっ

科学的手法を用い、容量特性やサイクル特性に優れたている。より一層の性能向上のため、正極、負極、電解材料設計を行い、実験にフィードバックしていく。多数の質など、リチウムイオン電池の様々な部材の研究開発

元素種と原子数の取り扱いを検討するため、通常のが進められている。

PC だ け で は 計 算 資 源 に 限 界 が あ る 。 そ こ で

次世代電池研究開発センターでは、正極活物質材

TSUBAME を使用することで、大規模な計算を短期間料の研究開発を進めており、特に近年注目を集めていで行い、実験では検討しきれない部分や、実験では知る、ポリアニオン系の Li2MSiO4(M: 遷移金属)の可能

り得ない特性などを分析し、開発を促進することが目的性を検討した。Li2MSiO4 は組成式内に Li を 2 つ含んである。材料開発の発展のために、大規模計算およびでおり、2 電子反応を起こす高容量材料として期待され高速計算による設計の精密化は不可欠であり、実験的ている。理論容量は 330 mAh/g を持ち、低コストかつには観測することができない情報を計算により得ること

熱安定性に優れ、多くの大学や企業が研究開発を行っで、他の材料への更なる開発を促進するとともに、高性ている。

能・高効率な二次電池を実現していく。

正極活物質 Li2MSiO4 において、Li2FeSiO4 は 1 電子反応までしか起こらず[1]、Li2MnSiO4 はサイクル特概要

性に課題があることが知られている[2]。これらの課題

リチウムイオン二次電池用の正極活物質 Li2MSiO4を解決するにあたり、M や Si に代わる置換元素を選定(M: 遷移金属)の第一原理計算を行う。図 1 に示す空1

（様式 20）成果報告書
間群 Pmn21 の結晶構造を用い、正極材料 Li2MSiO4

効率を向上させることができ、npool はその core の分の安定結晶構造と電子状態、および凝集エネルギーを割数を表す。npool の最適値は k 点数や core 数に依シミュレーションにより検討した。また、リチウムイオン

存するので、系を変えるたびに最適化が必要である。

脱挿入による充放電での結晶構造の変化と安定状態

図 2 に示すように、標準 node 数および標準 npoolを計算し、系の持つ平均電位計算を行った。本報告で数の使用に対して、node 数を増やすことで計算時間を構築した計算手法などを用いて、将来的には大規模な大幅に短縮することができたが、一方で TSUBAME ポ

系の取り扱いを行うことにより、より実際の系に近い計

イントは若干多く必要となった。npool 数を改善すること算を行い、実験値との比較を行うことは長期的な目標

で計算時間は更に短縮され、かつ使用したである。

TSUBAME ポイントも標準条件よりは少なく抑えることができた。従って本利用課題においては、この計算条
件を基本として検討を進めた。
表 1 および表 2 に、LixFeSiO4 および LixMnSiO4の構造緩和計算結果による格子定数の値を示している。
構造緩和計算の結果、文献値とほぼ同等の値を得るこ

図 1. Li2FeSiO4 の結晶構造. (空間群: Pmn21)
計算時間, TSUBAMEポイント

120%

結果および考察
第一原理計算ソフトは Quantum Espresso を用い、500 eV のカットオフエネルギーを有する平面波基底関数を用い、k 点サンプリング数は構造最適化においては 2×2×2 を、エネルギー算出時においては 4×4×4

100%
80%

計算時間
TSUBAMEポイント

60%

40%
20%
0%
標準node数
標準npool数を使用した。交換相関ポテンシャル GGA+U では U =node数改善後
標準npool数

node数改善後
npool数改善後

5 eV を使用し、電子スピンを考慮した。シミュレーション

図 2. Li2FeSiO4 のエネルギー計算における node 数との初期構造には空間群 Pmn21 を持つ結晶構造を用い、

npool 数の関係と計算時間と TSUBAME ポイント.

格子定数の初期値は文献値[3](計算値)を使用した。
空間群 Pmn21 の単位格子内には Li が 4 個、遷移金

表 1.LixFeSiO4 の格子定数の文献値[3]と計算値[Å].

属及び Si が 2 個、O が 8 個含まれており、全部で最大

Li2FeSiO4

LiFeSiO4

16 個の原子を取り扱い、計算を行った。

文献

計算

文献

計算

文献

計算

FeSiO4

上記計算を行うに当たり、本利用課題で使用できる

a

6.32

6.27

6.08

6.05

6.07

6.05

計算資源(TSUBAME ポイント)を有効に利用しつつ、

b

5.38

5.35

5.63

5.63

5.61

5.73

計算時間をなるべく短くするため、一度の計算で利用

c

5.00

4.99

5.04

4.98

5.30

5.49

する node 数と npool 数の最適化を行った。使用した 1
node に対しては 8 core を基本として検討を行った。計

表 2.LixMnSiO4 の格子定数の文献値[3]と計算値[Å].

算対象には、Li2FeSiO4 のエネルギー計算を採用した。

Li2MnSiO4

LiMnSiO4

図 2 に標準 node 数および標準 npool 数を用いたとき

文献

計算

文献

計算

文献

計算

a

6.37

6.31

6.26

6.21

6.17

6.12

b

5.43

5.41

5.42

5.40

5.70

5.72

c

5.04

5.02

5.15

5.18

4.85

4.96

の計算時間と消費した TSUBAME ポイントを示した。
ここで Quantum Espresso で設定できるオプションの一つ「npool」について説明する。各 k 点に対して、
core を分割して割り当てた並列計算を行うことで、計算

2

MnSiO4

（様式 20）成果報告書とができ、上記計算条件で第一原理計算による物性予5.5

測に十分な値を得られることが分かった。

5.0

図 3 には、第一原理計算により平均電位を求める方
4.5
V[V]

法の概念図を示している。第一原理計算により、組成式内の Li が 2 個、1 個、0 個のときの、系の構造最適

4.0

化と凝集エネルギーを求め、更に Li 単体での構造最3.5

適化と凝集エネルギーを求める。それらの値を図 3 に計算結果 V1
計算結果 V2
実験値 V1
実験値 V2
文献値 V1
文献値 V23.0

示す Nernst の式に当てはめ平均電位を算出する。

Mn

表 3 には、LixFeSiO4 および LixMnSiO4 の平均電

Fe

Co

Ni

図 4. 平均電位計算の結果と文献値[3]との比較.

位の計算結果を示している。文献値(計算値)と比較す
ると、0.1～0.2 V 程度の差異はあるものの、概ね文献

まとめ、今後の課題
リチウムイオン二次電池用の正極活物質 Li2MSiO4

値と同様の値を得ることができ、上記計算条件にて平均電位の見積もりが可能であることが分かった。

について、第一原理計算を用いた結晶構造の構造最図 4 には、上記と同様の方法で LixCoSiO4 および

適化と、凝集エネルギーの計算による平均電位の算出

LixNiSiO4 についても平均電位計算を行った結果を、を行った。文献値を初期値として、それらの値と比較を

LixFeSiO4 および LixMnSiO4 の結果とともに示してい行ったところ、文献値とほぼ同等の平均電位を得ることる。LixCoSiO4 は V1 および V2 のどちらも文献値の値かができ、TSUBAME を用いて行った本報告における計ら約 0.3 V 低い値となった。文献[4]によっては 0.1V 程算手法で活物質の物性予測が行えることが分かった。

度の差しかないものもあるが、LixCoSiO4 については

今後はこれらの計算手法を用い、更に他の M や Si

良好な実験結果が存在しないため、この計算結果が妥に置換する元素を検討し、活物質の材料開発にフィー当な結果であったか判断は難しい。LixNiSiO4 についドバックしたい。更に、単位格子を複数個用意した超格ては文献値とほぼ同等の値を得ることができ、本報告

子の結晶構造を用い、大きなセルを用いたシミュレーシの方法で系の予測を行えることが分かった。

ョンを行うことで、実際の系に近い状況を作り、電子状
態の細かい変化などを追えるようにしたい。その中で、

V
Li

充放電に伴う Li 拡散の様子や、Li 脱挿入に伴って起G Etotal  E (Li x MSiO 4 )  E (Li x 1MSiO 4 )  E (Li)nchargechargeV1

V2

discharge

A. Li2MSiO4

discharge

B. LiMSiO4

こる結晶構造の変化などの観察、更には Li の拡散に伴う活性化障壁の計算や、Li の移動パスの予測など、実験的には観測が難しい部分のメカニズム解明に適用

C. MSiO4

していく。

図 3. 平均電位計算の概念図.

これら超格子の計算や、Li の拡散パス・挙動に関す
る計算については、大規模な計算を行う必要があり、

表 3. LixMSiO4 の平均電位[V]. 文献値は[3]の値.
Fe

かつ企業の研究開発においては特にその迅速性が求められる。これらを満足するために TSUBAME のよう

Mn

文献

計算

文献

計算

なスーパーコンピュータを使用することは非常に有効で

V1

3.2

3.42

4.1

4.18

あり、今回の検討により、その有用性を確認することが

V2

4.8

4.70

4.5

4.38

できた。今後、上記のような更なる課題の検討を進め、電池材料のみならず、これらの計算ノウハウを他の材料へも展開していきたい。

3

（様式 20）成果報告書
参考文献[1] A. Kojima, T. Kojima, T. Sakai, J. Electrochem.
Soc. 2012, 159 (5): A525-A531.
[2] T. Muraliganth, K. R. Stroukoff, A. Manthiram,
Chem. Mater., 2010, 22 (20), 5754–5761.
[3] M. E. Arroyo-de Dompable, M. Armand, J.M.
Tarascon, U. Amador, Electrochem. Comm. No.8,
2006, 1292-1298.
[4] G. Zhong, Y. Li, P. Yan, Z. Liu, M. Xie, H. Lin, J.
Phys. Chem. C 2010, 114, 3693-3700.

4

（様式第 20）成果報告書

平成 26 年度 TSUBAME 産業利用トライアルユース 成果報告書

利用課題名 無機材料開発への第一原理計算の活用
英文： First principles calculations for development of inorganic material
松成 秀一
Shuichi Matsunari(株)ニコン
Nikon

固体酸化物の第一原理的電子状態計算は、密度汎関数理論の発展により、数百原子程度のセルにおいても適用可能となってきている。しかし、一般的な交換相関相互作用近似ではバンドギャップが小さめに計算される傾向にある。本課題では、40 原子イットリア（Y2O3）のバンドギャップ計算に G0W0 近似を適用し、その近似の有効性および、並列計算の有効性を検証した。まず GW 計算による補正効果を検証した結果、一般的な汎関数(GGA)では実験値
に対し 70%程度の計算値になったのに対し、GW 計算では実験値に近い計算結果が得られた。今回計算に用いた条件では、36core 並列まで直線的なスピードアップが見られ、並列計算が十分有効であることが分かった。
In this report, we discuss band gap calculations of a 40-atom yttria (Y2O3) crystal within the G0W0
approximation and the parallelization efficiency of those calculations. The calculated band gap isconfirmed to be improved by the approximation. When four hundred bands are used, a linear decreasein the total calculation time is observed for parallelization with up to 36 cores (3 nodes). This improved
efficiency is important for larger systems.
Keywords: 電子構造, バンドギャップ, GW 近似, 第一原理計算

界を乗り越えるため、交換項に HF 近似を混合したハイ

背景と目的
第一原理電子状態計算は方法論の発展と計算機のブリッド汎関数を用いる方法など高精度な計算手法が高性能化により、欠陥の形成・ダイナミクスから光吸収開発されている。その中で、本研究では多体量子論にスペクトルなどの光学特性まで幅広く、応用へ向けた研

基 づ く 摂 動 補 正 を 加 え る GW 近 似 に 注 目 し 、究が活発に進んでいる。電子系の交換相関相互作用

TSUBAME を用いてその有用性を検討した。

に 対 し 、 局 所 密 度 近 似 (LDA) や 一 般 化 勾 配 近 似

GW 計算は一般的な GGA 計算よりも負荷が大きく、(GGA)は第一近似として広く用いられ、結晶の格子定数 10 原子規模の計算でも並列計算を有効に使う必要数・結合距離等の構造的物性や化学反応過程におけがある。本プロジェクトでは、無機材料における GW 計

るエネルギー変化、相安定性などの算出について、信算、及び GW 計算並列化の有効性を検証することを目

頼性の程度も調べられており、実験だけからはわから的とし、40 原子の酸化イットリウムを例として実施した。

ない状態の検証・予測に応用することも可能となってきている。しかし、電気伝導特性・光学特性にとって重要概要

な物性値であるバンドギャップは実験値に対し大きく過本研究では、図 1 に示さている酸化イットリウム Y2O3

小評価されるという問題がよく知られている。実験値がについて、GGA 及び G0W0 計算によりバンドギャップをわかっている系であれば、それを用いて経験的に補正調べた。バンドギャップの GGA 計算値は実験値のをすることも可能であるが、実験的な観測結果がない

70%程度であるが、GW 計算による値は実験値に近く、系や不純物や結晶性の影響で材料本来の値と異なっ

高精度計算による改善が見られた。今回用いた G0W0ている実験例もある。よって、第一原理計算を無機材

計算は無機材料におけるバンドギャップの予測に利用料開発に活用する上で、バンドギャップを高精度に予できる。

測することは重要な課題である。LDA, GGA による限

また、今回検証した条件では 36 コア並列までスピー

5

（様式第 20）成果報告書
ドは直線的に上昇した。この結果により TSUBAME の結果および考察ようなスパコンを用いれば、今回検証したより大きなサ

GGA 計算によって得られた Y2O3 (conventional

イズの結晶セル（より多くのバンド数が必要となる系）のcell)のエネルギーバンド図を図２に示す。伝導帯の底計算も可能であることがわかった。

はガンマ点にある。一方で、価電子帯の頂上は詳細に見るとガンマ点からわずかにずれている。しかし、ガンマ点付近ではかなりなだらかになっており、エネルギー
的にはほとんどガンマ点と変わらないため、本研究のGW 近似ではガンマ点でのギャップを計算し比較する。
GGA 計算によって得られたバンドギャップの値は 4.08
eV であり実験値 5.8 eV[6]に対して、70%程度過小評
価されている。

図１．酸化イットリウム(Y2O3; Bixbiyte, 希土類 C 型蛍石）の結晶構造 (conventional cell, 80 原子)。緑色の球がイットリウム、赤色が酸素を表す。

計算方法
ABINIT コードを用いて計算を行った[1,2]。交換相

図２．Y2O3 エネルギーバンド構造。

関汎関数には GGA を用い、内殻電子はノルム保存擬
ポテンシャルを用い明示的には扱わない[3]。GW 計算

GW 計算では摂動計算のための非占有軌道の数にでは DFT 計算から直接得られるグリーン関数 G0 と、対する収束性が問題となる。近年、Bruneval ら が開分極関数、誘電関数の計算を介して得られる遮蔽クー

発した手法により、バンド数に対する収束性が著しく

ロン相互作用 W0 を用いた G0W0 計算を使用した。その向上することがわかっている[4]。ABINIT ではこの手

際、G0W0 計算に使用するバンド数を少なくする手法を法を実装しており、GW 計算の実用を後押ししている。

一貫して用いた[4]。Y2O3 結晶の格子定数には GGA

図３に Y2O3 の系でバンド数に対するバンドギャップ値

で 緩 和 し た 値 10.878Å を 用 い た 。 こ れ は 実 験 値の収束性を検討した結果を示す。占有軌道のバンド数

10.60Å[5] に 対 し て 2.6% 程 度 の 相 対 誤 差 で あ る 。

は 96 であるので、4 から 5 倍程度のバンド数で収束しG0W0 計 算 で は 、 同 じ 格 子 定 数 を も つ 40 原 子 のていることがわかる。

primitive cell を用いた。G0W0 計算における逆格子空バンド数 500 での G0W0 計算によるバンドギャップ値

間 Brillouin zone 内の平均には、4×4×4 の k 点メッは 5.76 eV となり、実験値（5.8eV[6]）と近い値が得らシュを用いたが、計算時間ベンチマークでは 2×2×2れた。その他の半導体・酸化物系でも系統的に、G0W0の k 点メッシュを用いた。

計算により約 10%以内の誤差に収まることが知られており、今回検証した Y2O3 の系でも、妥当な計算が可能であることが確認できた。

6

（様式第 20）成果報告書
今回、摂動補正の最低次である G0W0 計算を用いたが、より高精度にバンドギャップを計算するためには、
セルフコンシステントに計算する方法や vertex 項を加える方法などを検討することが必要である。また、今回は価電子帯頂上と伝導帯の底の DFT 一電子エネルギ
ー（固有値）への補正を計算し、バンドギャップのみに注目したが、バンド図を計算する場合はワニエ関数を使用するなどの工夫が必要になる。

図３．Y2O3 におけるバンドギャップ値のバンド数に対する収束性。

以上の結果を踏まえ、並列計算の有効性の検証を行った。バンド数は 400 とし、2×2×2 の k 点メッシュを用いた。ベンチマーク計算を行った結果を図４に示す。
ノード数( 並列コア 数)と し ては 1(12), 2(24), 3(36),4(48), 6(72)で計算を行い、DFT 波動関数を求める部分も含めた G0W0 計算全体の計算時間をまとめた。図４(b)から、今回用いた条件では 36 コア(3 ノード)並列程度までスピードは直線的に上昇し、それ以上で並列効率が著しく下がることがわかる。12 コアより小さな並列
数では 4 日以内で計算が終わらない可能性があるため検証を行っていないが、GW 計算において並列計算が有効に利用できることがわかった。より大きなセル（バ

図４．G0W0 計算ベンチマーク結果。(a)計算コア数と全ンド数）の計算では、今回よりも優れた並列効率が期待体の計算時間。 (b)ノード数と(1 ノード計算に対する)ス

できる。

ピードアップ。

まとめ、今後の課題

参考文献

無機材料のバンドギャップ値の予測に関して、GGA
計算では不十分であるが、G0W0 計算を用いることで、[1] X. Gonze et al., Z. Kristallogr. 220, 558—562 (2005).

GGA 計算より信頼性の高い結果が得られ、40 原子程[2] X. Gonze et al., Phys. Comm. 180, 2582—2615 (2009).

度のセルでも実用的な時間内で計算が可能である。よ[3] M. Fuchs and M. Scheffler, Comput. Phys. Commun. 119,

り大きなセルでの並列効率は今回よりも優れているた67 (1999).

め、TSUBAME のような 100 コア並列も可能なスパコ[4] F. Bruneval and X. Gonze, Phys. Rev. B 78, 085125 (2008).

ンを使えば、GW 計算を材料開発へ活用することが十[5] E. A. Stern, M. Newville, B. Ravel et al. Physica B,

分可能であると考えられる。100 原子程度の GW 計算

208—209 (1995).

は荷電点欠陥の形成エネルギー計算にも適用可能で[6] J. Robertson, J. Vac. Sci. Tech. B, 18(3), 1785—1791 (2000).

ある。

7

8

（様式 20）成果報告書

平成 26 年度 TSUBAME 産業利用トライアルユース 成果報告書

電子写真システム設計における電磁場計算の高速化
The Study of Speed up for an electrostatic field simulationin electrostatic photography system design
利用課題責任者 羽山 祐子
Yuko Hayama
株式会社 リコー
RICOH Co.Ltd
http://www.ricoh.co.jp
電子写真方式を利用した画像装置内での用紙詰まりの要因となる，用紙巻き付きのメカニズムを解明する
為に，電界計算と用紙変形計算を連成した解析ツールを開発した．このツールは，用紙の搬送速度と電荷の移動速度の差が大きく，連成時に膨大な計算回数が必要となり，計算時間が実用的な解析ツールとしてのレベルに達していなかった．今回の利用課題では，並列化と数値解析ライブラリを利用し，用紙巻き付きシミ
ュレータの計算速度の向上を目指した．その結果，計算速度の向上が実現し，巻き付き解析ツールとしての，実用的な適用の可能性を確認することが出来た．
The simulation tool has been developed to figure out the mechanism of paper deformation at second transfer
processes in electrophotography system. The simulation demands enormous time steps, because the electric charge
speed is almost million times faster than the paper feeding speed. It is difficult to practically use at design stage. In use
challenge of TSUBAME, we aimed to improve computational speed of the simulation tool by using parallelization andnumerical analysis library．As a result, it could be showed the possibility of practical applying of the numericalmethods.

Keyword: paper deformation, electrophotography system, numerical calculation, parallelization

背景と目的

界計算にオームの法則による部材内の電荷移動や，ロ

用紙に文字や画像を印刷する画像形成装置では，ーラと用紙間の放電現象を考慮しており，実機内の詳

高い品質を提供し続けるために，用紙の搬送状態を細な現象まで再現するように設計された．

管理することが重要である．Fig.1 に示すような電定着

子写真方式では，電界を利用して画像形成を行うこ

電源とから，用紙が周囲の電界によって変形するため，
印加ローラ

中間転写ベルト

用紙挙動はより複雑なものとなる．例えば，中間転A

写ベルトから用紙に画像を転写するプロセスでは，

st1 転写ローラ

用紙周辺に高い電界が生じるため，適切なバイアスを用いなければ，用紙がローラに巻きつく不具合が

(Y)

(M)

(C)

(B)
nd
用紙 2 転写ローラ

生じる．このような用紙の静電巻き付きは，電子写帯電器

真方式の不安定性をもたらす要因であり，高品質を現像器達成するためには，そのメカニズムを解明すること

感光体
用紙トレイ

Fig.1 電子写真装置の概略図が重要となる．

一方でこのツールには計算時間に課題を残す．本ツこのような静電巻き付きのメカニズムを解明するため

ールは一回の計算に数日程度の計算時間を要する．よに，我々はこれまでに，電界計算と用紙の変形計算にって，用紙やローラの物性値，プロセス条件等を振ってよる連成解析ツールを開発してきた．このツールは，電9（様式 20）成果報告書計算が可能な，解析ツールとしての実用化にまでは至の電気抵抗，厚み，などによっても異なることが分かっっていない．これは，用紙の搬送がミリ秒オーダーで起ている．

こるのに対して，電荷移動がナノ秒オーダーの現象で印加電圧あることから，十分な計算精度を得るためには膨大な計（V）

算ステップが必要になることが原因である．

-1000

-2000

-3000

-4000

-5000

紙

実験結果

今回の利用課題では，TSUBAME を利用して，並列

2 次転写ローラ

印加ローラ＋ベルト

Fig.3 印加電圧による用紙変形の違い

化と数値計算ライブラリを利用したチューニングによる用紙巻き付きシミュレータの計算速度の向上を確認し，《解析モデルと手法》

静電巻き付き解析への実用化の見込みを得ることを目

静電気力による用紙の変形計算の概要と手順を説

的とする．

明する．
計算の対象は，Fig.1 内の A の転写領域であり，印

概要

加ローラ，転写ベルト，用紙，各物質間の空気層，2 次

電子写真プロセスでは，画像が形成される過程で，転写ローラの層により形成される．Fig.4 に今回対象と電場が重要な役割を果たす. その為，電場の状態を管した，解析モデルのメッシュ例を示す．このモデルにて，理することが重要である．しかし一方で，電場を直接観各ローラ，転写ベルト，用紙内の電荷移動と，パッ察，計測することが困難である．したがって，シミュレー

シェンの法則に基づく各物質間での電荷移動を考慮

ションを利用した現象解析，予測のニーズは高い．

した差分法による電界計算が実施される〔1〕
〔2〕
〔3〕
．

《解析対象となる現象》

静電気力は，計算された電界強度と用紙上に発生し解析対象である，Fig.1 内の A の領域について，Fig.2

た電荷により計算される〔4〕．

に実験機による可視化結果を示す．用紙は正しく搬送
空気され，正常に次のプロセスに送られていることが分かる.
用紙

用紙
空気

印加ローラ
+転写ベルト

進行方向

nd

2

印加
ローラ

2

nd

転写ローラ

転写ローラ

Fig.4 解析モデルメッシュ例

Fig.2 転写プロセス内の用紙搬送状態可視化結果

用紙変形は，用紙を梁とみなし，計算された用紙

しかしながら，搬送される用紙が薄く腰が弱い場合に上の静電気力により，一次元の FEM 計算で求める〔5〕。

.2 用紙進行状態の例（実験）印加なし
は，印加した電圧により発生した電界により，用紙が変形全体の計算のフローを Fig.5 に示す．初期は用紙して搬送される．Fiｇ.3 に印加電圧を，-1000V～-5000V

がローラ間のニップ部から排出した時の解析モデルの範囲で異ならせた条件での，用紙搬送の可視化結果を作成する．そのモデルに対し，オームの法則，放を示す．-1000V の印加条件では，用紙が印加ローラ側

電を考慮した電界解析を行う．電界解析では，電荷に大きく変形していることが分かる．印加電圧を下げていの移動時間を考慮して，計算のタイムステップは 10-6～

くと，変形量が小さくなり，-4000V では変形方向が 2 次転

10-8 秒に設定している．電界解析の計算を一定時間実

写ローラ側に変化する．

施し，計算が収束した後，電界内で用紙に働く静電気

このように，用紙の変形挙動は，印加する電圧によって力を計算する．次にその静電気力に対して，用紙変も変形度合いや，変形方向が異なってくる．また，用紙

形計算を行う．その変形位置に合わせて，電界計算

10（様式 20）成果報告書
直接法では 2.3 倍の計算速度向上を達成した．

の領域のメッシュを作りなおし，次の静電界解析を行う．このサイクルを一定の用紙排出量ごとに繰り返した．計算は所定の排出量において終了させる．

5.0

計算速度倍率4.0

start
メッシュモデル作成

3.0
2.0
SOR

1.0

直接法

静電界計算
0.0

オーム法則,放電の考慮

0Convergence

4

8
12
スレッド数

16

Fig.6 電界計算の計算速度

静電気力算出次に，二つのソルバを利用して計算速度向上の最適化を行った．検討の結果，Red-Black SOR 法をベースと用紙変形計算して，一定の反復回数を超えたときのみ，直接法を利用
用紙の移動するハイブリッドソルバが最も計算時間を短縮できることが分かった．SOR 法の計算時間は，その反復計算の初

用紙進行終了

期値に大きく依存する．よって，電界計算内で放電判定が生じた場合や，用紙変形計算を行った直後は電界がend

大きく変化するため，前計算結果を初期値として利用しFig.5 静電気力による用紙変形計算のフロー

た場合，収束には多くの反復が必要になる．一方で，直結果および考察

接法は常に一定の計算時間で解を得ることができる．よ本計算手法のホットスポットである，電界計算の連立って，SOR 法の反復回数に閾値を設けて，一定回数を方程式ソルバについて２種類の高速化を検討した．

超えたときのみに直接法を利用するように最適化した．

一つ目は，Red-Black SOR 法を用いた並列化である．

最後に，本シミュレータを利用して印加電圧条件を振差分法の領域を 4 分割して，OpenMP を利用した並列った際の用紙巻き付きの傾向について，実験結果との化を行った．

比較を行った．計算条件は Table.1 の通りである．各印

二つ目は，数値計算ライブラリである LAPACK[6]を加電圧の用紙巻き付きの計算結果を Fig.7 に示す．結利用した方法である．本計算モデルでは係数行列が帯果を，Fig.3 に示した実験結果と比較すると，本シミュレ

行列となるので，帯行列用の直接法ソルバを利用した．

ータにより，用紙巻き付き方向の傾向を再現できたこと

また，本ソルバは標準で OpenMP の並列化に対応してが分かる．また，本計算は先述したハイブリッドソルバをおり，並列化による高速化も検討した．

利用することで，高速化を行う前と比べて，計算時間を以上の高速化手法を静電巻き付き解析シミュレータ

8 分の 1 に短縮することが出来た．

に実装した．以下にテスト計算の結果について説明する．
まず，二つのソルバの並列化効率を確認した．Fig.6に 16 スレッドまでの計算速度倍率を示す．Fig.6 より，両ソルバで 16 スレッドまで並列化による計算速度向上が確認できる．それぞれ，16 スレッドで SOR 法は 3.8 倍，11（様式 20）成果報告書
Table 1 計算条件
-1000, -2000, -3000, -4000, -5000

印加電圧[V]

用紙特性
厚み[um]
電気抵抗[Ω m]
ヤング率[N/m2]
密度[kg/m3]
印加電圧[V]-1000

70
5×10 8
6.48×10 9
845
-2000

-3000

-4000

-5000

計算結果

Fig. 7 用紙巻き付き計算結果

まとめ，今後の課題
本プロジェクトでは，TSUBAME を利用し，静電巻き付きシミュミレータの高速化を確認した．電界計算に，Red-BlackSOR 法と，直接法のハイブリッ

ド化を利用し計算速度の向上を図り，最終的に，8倍の高速化を達成した．本結果により，静電巻き付
きシミュレータの実用化への見込みが得られたと判断する．
今後，本シミュレータをメカニズム解析や，設計開発へ展開していく．

参考文献
[1]

長尾剛次他,静電転写部における用紙変形の基礎的解析

Japan Hardcopy 2003 Fall 論 文 集

P.53-56[2]

Toyoshige Sasaki etc. Second Transfer Process
Simulation Ⅱ

Society for Imaging Science andTechnology Nip24 and Digital Fabrication 2008
p.342—345
[3]

門永雅史
ーション

中間転写ベルトの２次元電界シミュレ
日本画像学会誌

244 号

2003

p.128—135[4]

Ｈ.Ｈ.ウッドソン,J.Ｒ.メルヒャー
キュラム電気力学
材料力学

MIT コアカリ

Ⅱ 産業図書

[5]

大学課程

オーム社

中山秀太郎

[6]

Anderson E., Bai Z., Bischof C., Blackford S.,
Demmel J., Dongarra J., Du Croz J, Greenbaum A.,
Hammarling S., McKenney A., Sorensen D., LAPACK
User’s Guide 3rd edition, SIAM, 2000

12

（様式 20）成果報告書

平成 26 年度 TSUBAME 産業利用トライアルユース 成果報告書

利用課題名 ホモジナイザーを用いた撹拌混合による乳液製造のスケールアップに関する解析

CFD for an establishment of scale-up procedure in manufacturing of emulsion by using of
homogenizer
利用課題責任者 横川 佳浩
Yoshihiro Yokokawa
所属 株式会社資生堂 生産技術開発センター
Technology & Engineering Center
Shiseido Co., Ltd.
http://www.shiseido.co.jp
乳液の製造においては、撹拌混合の強さによって、得られる乳化粒子の大きさが異なり、その結果、粘度等の物性が変化し、使用感の全く異なるものが出来上がる。機械力による乳化の場合、設備のせん断力が乳化の進行に寄与していると考えられる。乳液の調製に広く使われる設備として、ホモジナイザーがある。粒子法に基づく流体解析により、非ニュートン高粘度流体が、ホモジナイザーを通過する際に受ける最大せん断力を算出した。その結果、流体の受けるせん断力は分布を持ち、回転数が大きいほど分布幅は大きいことが判った。回転数の違いによるせん断力の強さおよび分布の差は、サイズ違いのホモジナイザーでも生じる。少量スケール用
ホモジナイザーで調製される乳液を大容量用ホモジナイザーを用いて再現するためには、今回、得られた最大
せん断力分布が有効な指標になると考えられた。
In manufacturing milky lotion, that is emulsion, the size of emulsion particle is determined by thestrength of mixing equipment. For details, the shear force of the mixer is thought to contribute theprogress of emulsification. The homogenizer is widely used to produce emulsion. We found out
the strength and distribution of the shear force of liquid treated by the homogenizer in various
rotating speed. As a result, the acquired data of the shear force led us the understanding theability of emulsification in the different size of homogenizers. The analysis of the shear force data
could be useful indicator in the scale-up process.

Keywords: 流体解析、粒子法、乳化、ホモジナイザー、せん断

背景と目的るホモジナイザー【注釈 2】が生み出すせん断力を明ら

化粧品、食品を初め、広く製造業では、液体の撹拌混

かにすることを目的とした。新規処方の化粧品中味の合により目的物を調製することが多い。化粧品の製造

製造においては、まず少量スケール（1 L 程度）で確立

においては、本来、混合されにくい複数の原料を均一された調製方法を参考に大容量（50～1000 L）の製造に混合し、目的とする混合物（化粧品中味）を得るが、釜を用いた製造工程にスケールアップすることとなるが、全く同一の原料を同じ比率で使用しても、混合の仕方

少量用および大容量用の異なるサイズのホモジナイザ

によって得られる混合物の物性（使用性）は大きく異なーのせん断力を明らかにすることで、スケールアップ時る。特に乳液・クリームの調製においては、強いせん断の理論的な指標が設定可能になると考えられる。ホモ

力をかけて乳化【注釈 1】することが多いが、得られたジナイザーは、0.5 mm の狭い空間に強いせん断力を乳化粒子の大きさによって粘度等の物性が異なり、製かけて乳化を引き起こす設備である。解析においては、

品特性に大きな影響を与える。そのため、せん断力の

この微小空間を通過する流体を解析可能な精度で設制御は極めて重要となるが、実際の製造設備が液体に

備全体の解析を行う必要があるため、TSUBAME によ与えるせん断力を測定することができないため、理論る大規模計算を行うこととした。

ではなく実験結果から製造工程を設定することとなる。

本研究では、液体流れの数値解析に有効な粒子法に我々は、流体解析を用いて、乳化工程によく用いられ基づく商用プログラム 「Particleworks」を用いて解析

13

（様式 20）成果報告書を行なった。さらに、非定常計算において、粒子が過去12.880 Pa・sec(1.03 1/s)に受けたせん断力を保持し、最大値を表示するモジュールも開発した。その結果、ホモジナイザーを通過した
流体が受けたせん断力の最大値の分布を明らかにすることができ、装置のせん断力を把握することができた。
今後、得られたデータをもとにスケールアップの指標を設定することが可能と考えられる。
【注釈 1】
乳化：水と油のように混ざり合わない 2 つの液体に界面活性剤を加え、一方の液体（連続相）中に他方が微粒

図 2.高粘度流体のせん断速度・粘度データ子（分散相）として分散して存在させること。化粧品の乳液、クリームには、W/O 乳化（連続相：水、分散相：油）か O/W 乳化（連続相：油、分散相：水）があり、これに粉末が加えられるものもある。
【注釈 2】
ホモジナイザー：本研究では、プライミクス社製 ホモミ
クサーMARKII をモデルとして解析を行った。
概要
乳液やクリーム等の製品の製造においては、乳化工程が存在する。乳化方法の中でも、設備の機械力による強いせん断力で乳化を行う場合、設備のせん断力を把図 3-1.回転数・吐出量データ（水）握することで、乳化粒子の制御が可能と考えられる。
我々は、粒子法に基づく流体解析ソフト「Particleworks」を用いて、乳化物の調製のために広く用いられるホモジナイザーのせん断力の解明を試みた。乳化に用いられるホモジナイザーは、攪拌翼が周囲の壁から 0.5 mm 内側を毎分、数千回転で回転する設備で、この高速回転により強いせん断力を生み出す
とともに流体を吸引・吐出し、製造釜内の流体の混合も促進する装置である。（図１）解析の対象とする流体は、ニュ

図 3-2.回転数・吐出量データ（高粘度流体）

ートン流体である水および乳液の代表として B 型粘度計で

図1
およそ 10 Pa・s の粘度特性を有する高粘度の非ニューナイザー下部の吸引部分に流入口を設定し、得られた吐出量を与トン流体（図２）とした。ホモジナイザーの吐出量およびえて計算した。（図４）乳化の進行

通過した中味のせん断速度分布を粒子法解析によりとともに生成する乳化粒子は小さ求めようとすると計算負荷が極めて高くなるため、吐出くなるが、乳化粒子の大きさは、量については、格子法による解析および実測により求流体が受けた最大のせん断力にめた。（図３）せん断力解明のための解析では、ホモジ

図414（様式 20）成果報告書
よって決定されると考えられる。定常状態での瞬間のせん断速度分布は通常の解析で得ることができるが、ホモジナイザーを通過する流体が受けた最大のせん断力を知る必要がある。そこで、粒子法解析において、各粒子が受けたせん断力の履歴を保持し、その最大値を示すモジュールも合わせて作成した。（図 5）

瞬間せん断速度コンター図

図 5.最大せん断速度表示モジュール・イメージ図

結果および考察
ホモジナイザー処理された流体が受けるせん断力を知るため、最大せん断力分布表示モジュールを用い、簡素化したモデルで解析を行った。水および高粘度流体の流入量をそれぞれ

8.65×10-4

図 6-1.最大せん断速度分布（水）m3/s、4.92×10-4

m3/s とし、7000 rpm で回転させたホモジナイザー通
過後の流体が受けた最大せん断力はそれぞれ（図６）に示す分布を持つことが明らかになった。いずれも正規分布に近い分布を示し、水と比較して高粘度流体の方が受けるせん断力は小さいことが確認された。次にホモジナイザーの回転数を変えて解析を行った。その結果、水および高粘度流体の最大せん断速度分布はそれぞれ（図７）に示す分布となった。いずれの流体も瞬間せん断速度コンター図

回転数が高いほど強いせん断力分布となり、分布の幅
は広がることが判った。実際の乳化工程においては、処理時間が長くなるにつれて、流体は何度もホモジナ
イザーを通過することになる。そこで、繰り返しホモジナ
イザーを通過した高粘度流体が受ける最大せん断速度分布と回転数の関係を調べた。理論上は低い回転数によって生まれる弱いせん断力でも時間をかけて繰り返しホモジナイザーを通過させることで高回転数で処理した場合と同程度のせん断力を与えることができることとなる。しかし、解析から得られた最大せん断速度分
布からは、化粧品製造の現実的な時間範囲では、回転

図 6-2.最大せん断速度分布（高粘度流体）

15（様式 20）成果報告書

瞬間せん断速度コンター図

図 8.パス回数による最大せん断速度分布の推移(7000rpm）

図 7-1.回転数・最大せん断速度分布データ（水）図 9.10 パス時の最大せん断速度分布の比較
大きさのホモジナイザーを用いることが多い。大容量用
瞬間せん断速度コンター図

ホモジナイザーは、少量用とほぼ同じ機構であるが、サイズが違うため、回転によって流体に与えるせん断力
およびその分布は異なる。乳液の処方ごとに、目的とする乳化物を得るために必要なせん断力を最大せん
断力分布から求めることで、スケールアップの指標を設定できると考えられた。

まとめ、今後の課題
我々は、粒子法解析により、ホモジナイザーが流体に与えるせん断力および本設備を通過した流体が、そ
図 7-2.回転数・最大せん断速度分布データの過程で受ける最大せん断力の分布を把握することに（高粘度流体）成功した。この結果から、目的とする乳液の製造のた数の違いによるせん断力の差を補うことができないとめに必要なせん断力を算出することができるようになっ考えられた。（図８）例えば、9000 rpm で 10 回ホモジナ

た。この必要せん断力がスケールアップの指標としてイザーを通過した高粘度流体の 90 ％以上が受けるせ有効と考えられた。今後、実際の乳液製造に使用するん断力を、7000 rpm による繰り返し処理で達成するこ

大容量用ホモジナイザーのせん断力について解析を行とは困難である。（図９）乳液の製造においては、これまうとともに、そこから算出したスケールアップ指標の設で解析してきた少量用ホモジナイザーよりおよそ 5 倍の定および有用性の検証を行う。

16

（様式 20）成果報告書

平成 26 年度 TSUBAME 産業利用トライアルユース 成果報告書

利用課題名 LTE-Advanced における大型車両内電磁界特性に関する基礎検討
英文 A Study on the Electromagnetic-field Characteristic in a Large-scale Vehicle in LTE-Advanced
利用課題責任者 斎藤 裕
Yutaka Saito
所属 株式会社パナソニック システムネットワークス開発研究所

Panasonic System Networks R&D Lab. Co., Ltd.
http://panasonic.co.jp/avc/psnrd/
邦文抄録（300 字程度）
近年，自動車から得られる多種多様な情報を活用するクラウドサービスが注目されている．無線通信シ
ステムとしては LTE-Advanced や無線 LAN などが想定されており，それぞれのシステムに適したアンテナ
を自動車に搭載することが必要となる．この自動車クラウド化の流れはバスなどの業務用車両から普及することが考えられ，大型車両に搭載されるアンテナ性能の評価が求められる．しかしながら，大型車両に搭載されたアンテナの放射性能を測定するためには多大な時間と労力が必要であり，また電磁界シミュレーションでは空間メッシュ数が大規模になる課題がある．そこで本報告では，LTE-Advanced を解析対象のシステムとし，3.5GHz 帯アンテナをバスに搭載した場合のアンテナ性能について，スーパーコンピュータ
TSUBAME2.5 を用いた電磁界シミュレーションにより解析する．
英文抄録（100 words 程度）
The cloud service system which uses various data from vehicles has been attracting much attention in recent years.
The use of LTE-Advanced and Wireless LAN is proposed for its wireless communication system. Since the cloudservice is introduced in advance into the vehicle for business use, the evaluation of the antenna performance installedon the large-scale vehicle, such as the bus and the truck, is required. However, the antenna measurement withlarge-scale vehicle involves immense amount of time and effort, and the antenna electromagnetic simulation requires
the enormous analysis meshes because of the volume of the vehicle. In this paper, we adopt the LTE-Advanced as thetarget system of the cloud service, and evaluate the 3.5 GHz antenna performance installed on the bus by electromagnetic simulation using supercomputer TSUBAME 2.5.
Keywords: Electromagnetic simulation, Large-scale vehicle, Cloud service, LTE-Advanced

1. まえがき用いてアンテナ性能を解析することが望ましい[2]．また，

近年，自動車から得られる多種多様な情報（ビッグデー

例えば LTE-Advanced や無線 LAN などの数 GHz 帯におけタ）を活用するクラウドサービスに注目が集まっている

る大型車両全体を含む電磁界シミュレーションでは，解析[1]．自動車のクラウド化を実現するための無線通信シス

空間を波長の数十分の一程度の大きさで分割する必要がテムとして，LTE や LTE-Advanced, 無線 LAN などが想定あることから，大規模な空間メッシュでの解析が必要となされており，それぞれのシステムに適したアンテナを自動る．現在，電磁界シミュレーションで用いられている一般車に搭載することが必要となる．この自動車クラウド化の的なスタンドアロン型計算機では，1 億個程度のメッシュ

流れは，一般乗用車への適用に先駆けてバスやトラックな数のモデルの解析が限界であるため，並列処理によって大どの業務用車両から普及が進むことが想定され，これらの規模メッシュの解析が可能なスーパーコンピュータの利

大型車両に各無線通信システム用アンテナを搭載した場用が有効となる．

合のアンテナ性能評価が求められる．

本 利 用 課 題 で は ， 基 礎 検 討 と し て ， 3.5GHz 帯 のしかしながら，大型車両に搭載されたアンテナの性能を

LTE-Advanced を解析対象のシステムとし，大型車両のバ

車両全体の影響を含めて測定するには大掛かりな設備とスに搭載されるアンテナの性能をスーパーコンピュータ

膨大な労力が必要となるため，電磁界シミュレーションを

TSUBAME2.5[3]を用いて電磁界解析する．実際の構造を

17

（様式 20）成果報告書
モデリングしたバスモデル車両内に 3.5GHz で動作するア

図 1 に本検討で用いる大型車両のバスモデルを示す．大ンテナを配置し，バス全体の影響を含めたアンテナ放射指きさは全長 8990mm×全幅 2315mm×全高 2980mm，ボデ

向性を解析すると共に，バス車両内の電磁界分布を可視化

ィ及びホイールは金属，座席シート，ダッシュボード，窓する．また，バス車両内の構造物の座席シートなどによるガラス，バンパーは誘電体でモデリングしている．

反射・減衰特性について解析し，その影響を検証する．以

図 2 にアンテナモデルを示す．3.5GHz で動作するアン上の検討により，使用ノード数と解析時間の関係を把握し，テナとして，一般的な車両搭載アンテナである 1/4 波長モ

電磁界シミュレーションにおける TSUBAME2.5 活用の有

ノポールアンテナを採用した．

効性について検討する．

図 3 にアンテナ配置位置を示す．バス車両内へのアンテ
ナ配置として，①前方ダッシュボード上配置，②後方ダッ2. 解析モデルと解析条件

シュボード上配置，③ルーフ中央裏配置，④ルーフ後方裏

表 1 に解析条件を示す．計算機システムはスーパーコン

配置の 4 通りを検証する．アンテナの配置向きは，①②のピュータ TSUBAME2.5，電磁界解析には有限積分法を用

ダッシュボード上配置では図 2 に示す座標軸の向き（アンいた電磁界シミュレータである CST 社の MICROWAVE

テナ素子がルーフ側）であり，③④のルーフ裏配置では図STUDIO[4]を使用する．解析周波数は LTE-Advanced を想2 においてアンテナを X 軸中心に 180 回転させた向き（ア定し 3.5GHz としている．

ンテナ素子が床側）である．
以上の解析条件でモデル 1 条件あたりの解析メッシュ数

表 1 解析条件

は 24.9 億メッシュとなり，一般的なスタンドアロン型計

計算機システムTSUBAME 2.5

電磁界解析

MICROWAVE STUDIO

解析周波数

3.5 GHz

車両モデル

バス

アンテナ

1/4 波長モノポールアンテナ

算機での MICROWAVE STUDIO で解析可能な約 1 億メッ
シュを大きく超えているため，TSUBAME2.5 での解析が必要である．

2mm①前方ダッシュボード上配置
アンテナ②後方ダッシュボード上配置

配置位置③ルーフ中央裏配置

100mm

21mm

④ルーフ後方裏配置
解析メッシュ数

Port

24.9 億メッシュ

8990mmGND

図 2 解析モデル（アンテナ）

2315mm

2980mm（a） 全景（b） 断面図① 前方ダッシュボード上配置② 後方ダッシュボード上配置③ ルーフ中央裏配置④ ルーフ後方裏配置

図 3 アンテナ配置位置

図 1 解析モデル（バス）18（様式 20）成果報告書
3. 解析結果

3.2. 車両内の電界分布の比較

3.1. アンテナ放射指向性の比較

図 5 に各アンテナ配置位置における車両内の電界分布

図 4 に各アンテナ配置位置におけるアンテナ放射指向の比較を示す．それぞれ XZ 面の電界分布であり，図中の性を示す．XY 面（水平面）と XZ 面（垂直面）を記載しピンク点はアンテナ配置位置を示している．③のルーフ中ている．車両搭載条件において，アンテナ配置位置によっ

央裏配置の場合に，比較的バス車両全体に電界が分布しやて放射指向性が大きく変化することが確認できる．XZ 面

すい傾向が確認できる。一方，①②のダッシュボード上配の放射指向性より，①②のダッシュボード上配置では天頂

置や④のルーフ後方裏配置では，アンテナが車両の端に配方向，③④のルーフ裏配置では地面方向に最大放射方向が置される条件のため，アンテナが配置される車両端の反対向くことがわかる．また，車両の影響で放射指向性に多数

側における電界強度が低いことがわかる．また，座席シーのヌルが発生している．

トや運転席背面の仕切り板の遮蔽の影響によって電界強YZ度が低下する傾向も確認できる．

Eq 成分
Ef 成分

3.3. 座席シートの影響

X

X

（dBi）

（dBi）

図 6 に座席シートの有無に対するアンテナ放射指向性の比較を示す．アンテナ配置位置は③のルーフ中央裏配置とした．座席シートを取り除いても，最大放射方向などの

-30

10

-30

放射指向性の傾向は大きく変化しないことが確認できる．

10① 前方ダッシュボード上配置

YZ

ただし，座席シートなしでは座席シートによる誘電損失分
Eq 成分
Ef 成分
大

X

X

（dBi）

-30

10

（dBi）
-30

① 前方ダッシュボード上配置

10② 後方ダッシュボード上配置

Y

Z

X

小

Eq 成分
Ef 成分

X
② 後方ダッシュボード上配置

（dBi）
-30

10

（dBi）
-30

10

Z

Eq 成分
Ef 成分

③ ルーフ中央裏配置

Y③ ルーフ中央裏配置

X

X

（dBi）
-30

10

（dBi）
-30

10

④ ルーフ後方裏配置④ ルーフ後方裏配置

図 4 アンテナ放射指向性

図 5 車両内の電界分布

19

（様式 20）成果報告書
だけアンテナ効率が改善するため，放射指向性が全体的にとを把握した．また，TSUBAME2.5 の 48 ノードを用いた大きくなっている．その改善量は今回の解析モデルでは約並列解析により約 25 億メッシュのモデルを約 5 時間で解2dB である．

析でき，電磁界シミュレーションにおける TSUBAME2.5

図 7 に座席シートの有無に対する車両内の電界分布のの有用性を確認した．

比較を示す．アンテナ配置位置は③のルーフ中央裏配置と本検討はスーパーコンピュータ活用による電磁界解析した．座席シートの遮蔽によって電界が減衰する傾向や，の基礎検討の位置付けであり，さまざまなアンテナ形式の座席シートからの反射波に起因する定在波が発生するこ差，乗客などの人体の影響，更には性能改善の検討などがとが確認できる．

今後の研究課題となる．

4. TSUBAME ノード数と解析時間の関係

参考文献[1] K. Mase, “Information and Communication Technology
and Electric Vehicles －Paving the Way towards a Smart
Community,” IEICE TRANS. COMMUN., Vol. E95-B,
No. 6, Jun. 2012.

図 8 に TSUBAME ノード数と解析時間の関係を示す．
解析モデルは③のルーフ中央裏配置のモデルを用い，解析
メッシュ数は表 1 の通り 24.9 億メッシュである．図 8 よ[2] S. Horiuchi et al., “Comparisons of Simulated and Measured Electric Field Distributions in Cabin of a Simplified
Scale Car Model,” IEICE TRANS. COMMUN., Vol. E90-B, No. 9, Sep. 2007.

り，ノード数を増やすことで解析時間が低減する傾向が確認できる．今回の解析モデルでは，48 ノード程度以上で解析時間がほぼ飽和すると考えられる．

[3] 東京工業大学 学術国際情報センター TSUBAME（http://www.gsic.titech.ac.jp/tsubame）

5. まとめ[4] CST MICROWAVE STUDIO（http://www.cst.com/）

本利用課題では，スーパーコンピュータ TSUBAME2.5

大を用いて，バスに搭載される 3.5GHz 帯アンテナの特性を
電磁界シミュレーションにより解析した．バス車両内のアンテナ配置位置により放射指向性が大きく変化し，車両中
央のルーフ裏にアンテナを配置する場合に比較的バス全③ ルーフ中央裏配置，座席シートあり

体に電界が分布する傾向である．座席シート等の遮蔽によ

小って電界強度が低下するため，車両の端にアンテナが配置される場合は，その反対側における電界強度が低くなるこ

Y

Z

Eq 成分
Ef 成分

③’ ルーフ中央裏配置，座席シートなし

図 7 電界分布に対する座席シートの影響

X

X
25000

（dBi）
10

（dBi）
-30

Caluculation time （sec.）

-30

10

③ ルーフ中央裏配置，座席シートあり

Y

Z

X

X

（dBi）
-30

10

Eq 成分
Ef 成分

（dBi）
-30

23000
21000
19000
17000
15000
16

10

③’ ルーフ中央裏配置，座席シートなし

24

32
40
48
56
Nodes for caluculation

図 8 ノード数と解析時間の関係

図 6 放射指向性に対する座席シートの影響

20

64

（様式 20）成果報告書

平成 26 年度 TSUBAME 産業利用トライアルユース 成果報告書

利用課題名 大規模シミュレーションによるレーダの車両搭載時の特性把握
英文：Evaluation of electromagnetic characteristics of radar within car environment by large scale
simulation
利用課題責任者 井上 大輔
Daisuke Inoue所属 古河電気工業株式会社
Furukawa Electric Co., Ltd.
http://www.furukawa.co.jp/
邦文抄録（300 字程度）
TSUBAME を用い、準ミリ波レーダの自動車搭載を想定した大規模電磁界シミュレーションを行った。昨年度のトラ
イアルに続き、さらなる規模拡大やモデル詳細化への対応として計算の効率化の妥当性確認と、実際の車両構造に即したモデル化スキームの開発と車両搭載を想定したシミュレーションを実施した。計算の効率化として CPU マルチ
コア駆動による TSUBAME の能力の十分な活用やメッシュの改善の効果が確認できた。また車両構造のモデル化
では新たに 3D スキャナを活用した大規模解析ルーチンの可能性や、より実物に近い解析を示すことができた。
英文抄録（100 words 程度）
The large-scale electromagnetic simulation which estimates the car installing situation of thequasi-millimeter wave band radar was performed using TSUBAME. Following the last year's trial, we
newly confirmed the validity of efficiency of calculation as correspondence to further scale expansion
and model detailing and developed the modeling scheme with the actual vehicle structure and thesimulation which assumed the car installing situation. We showed the enough utilization of
TSUBAME with driving CPU multi-core and the effect of the mesh improvement. Furthermore, we
showed the possibility of the large-scale analysis routine with 3D scanner newly and the more actual
analysis by modeling of the vehicle structure.
Keywords: ５つ程度

Electromagnetic, simulation, radar, bumper, 3Dることが現状の技術的な課題となっている。

背景と目的

交通事故死傷者数を減少させるための安全技術

従来は多大な工数の試作試験を実施することに開発や今後の段階的な自動運転へ技術開発が各国より、搭載時の特性の把握や改善検討を行ってきた。

で進められている。国内では、今後高齢ドライバーの一方、シミュレーションを用いた検証は、開発の加速増加する状況に対しても、予防安全技術の適用が望への期待は高いものの、準ミリ波波長で車体を解析まれる。現状、交通事故の約半数以上が交差点にする場合、その規模の大きさが課題であった。

おいて事故が発生しているとされており [1] 、周辺を昨年度の TSUBAME トライアルユースより、電磁見渡すことができるような広角なセンサは有用であ

界シミュレーションの実験結果との整合性の確認や、る。これまでにレーザレーダ、電波式レーダ、画像セ

GPU 活用による計算の高速化が確認できている。

ンサなど様々なセンサが商品化されてきているが、実際の車両規模を想定した大規模モデルであってもその中で電波式レーダは耐候性に優れているとされシミュレーションを開発ルーチンへの組み込むことがている。また、ミリ波に比べ、準ミリ波のレーダは一できる可能性を示せていた。
しかしながら、上記はあくまで簡易モデルでの解般的に安価であり、普及拡大が現在進んでいる。
通常レーダはバンパー裏側、何らかの車体構造

析にすぎない。今後、実際に即した解析構造におけ

物上に配置される。しかしながら、レーダはこれら搭る解析のためのモデルの詳細化や解析領域の拡大

載状態によって電波特性に大きな影響をうけるとさ等に対応するため、新たなモデル化スキームの開発れている[2]。この状況でも十分な電波特性を保持すやシミュレーションの効率化検証が必要であると考21

（様式 20）成果報告書えられる。本プロジェクトでは、TSUBAME により大今回は新たに 3D スキャナの活用を導入し、実際の規模シミュレーションを実施し、これら準ミリ波レーダ

車両約 1.5m 程度の範囲の金属ボディ、バンパーその自動車搭載にかかる現実的な課題に対して、シミれぞれを 3D データとして記録し、メッシュ化やファイュレーションのさらなる有用性について示す。

ル形式の変換等を行った上で MW-studio のモデラ
上へインポートを実現した。これらによりシミュレータ
上でのモデル化だけでなく、実物にある程度即した

概要

CST 社製 MW-studio ver2013 また ver2014 によ構造の解析ができることが期待される。なお、図 2 にるシミュレーションを TSUBAME 上で実施した。モデおいて、半透明で示すものはバンパーである。

ルとしては準ミリ波 ISM バンド 24.15GHz で動作するレーダ単体を基本状態（＝①）とし、その前面にバンパーや車体をごく簡易に模擬したプラスチック板と金属板を装荷したモデル（＝②）、また実際の車両構造物を模擬すべく新たに 3D スキャナを活用し車両
ボディ

構造をデータ化したモデル（＝③）を試行した。
まず、TSUBAME 上 MW-studio のバージョンア
ップによる大規模解析の効率化について確認する。

バンパー

ver2013 においては MPI クラスタ時に限り、マトリク
図 2 ③車両構造装荷ス計算(matrix cal)プロセス時に CPU がマルチコア駆動できずシングルコア駆動となる課題があった。
TSUBAME2.5 において CPU コア 12 ヶ（Thin ノー

結果および考察ド 2CPU × 6 コ ア ） が 本 来 駆 動 可 能 で あ る が 、図 3 に計算の効率化検討に関する結果を示す。

ver2014 からは MPI クラスタ時も同プロセスでマル

図 3(a)ではメッシュ数、図 3(b)では計算時間を各計

チコア駆動が可能となった。さらに、ver2014 では従算条件において示す。②の曲面板装荷モデルにお来のメッシュ生成設定が変更されており、新たな設

いて、各バージョン、各メッシュによる比較を行った。

定を利用することでメッシュ数の低減が可能とされてver2014 へのアップグレードにより、マトリクス計算プ

いる。上記 2 点で期待される計算効率化について②

ロセスにおける CPU マルチコア駆動による高速化

モデルにおいて確認を行う。

が確認できた。当モデルではトータル計算時間として
約 3 割程度の高速化が可能となり、TSUBAME クラ
スタのポテンシャルを十分活用できるようになった。
一方メッシュ生成に関しては、従来のオブジェクトの構造特徴点全てをメッシュ起点としていた 設定

図 1 ①レーダ単体（=legacy）と、新たに導入されたオブジェクトのエッジに接する面をメッシュ起点とする設定（=new）の差異を確認した。結果、今回の設定においては CST 社推奨のとおり、メッシュ数低減が図れ、それに応じた解析時間の低減も実現することができた。
なお計算結果に関しては、バージョン違いによる

図 1 ②曲面板装荷

特性差異はなく、メッシュ設定違いによる特性差異
はごく微小であったことは確認している（結果割愛）。

次に実際の車両構造の詳細化について検討行う。

22

また、ここでの解析はいずれも S キューでの実施、①

（様式 20）成果報告書においてのみ 6 ノード、②③は 24 ノードという条件でイアルにおける簡易板と同様の結果であるが、今回ある。

新たに実際の車両構造においても確認された。
なお、③モデルは 32 億メッシュであり、当方経験
上過去最大メッシュ数となっている。②モデルと比較
すると若干の増加にとどまるが、前述のメッシュ数低減等計算効率化を図った上での数字であり、実質上昨年度の 1.5 倍程度の規模のモデルの計算を実現することができた。ここで③がメッシュ数規模に対して計算時間が長い結果となっているが、これはメッシ
ュの各ノード振り分け時に偏りが発生したことと、
図 3(a) 各条件におけるメッシュ数

GPU メモリを一部オーバし CPU 計算を行ったことが要因である。この状況の改善やさらなる規模拡大の可能性として、より多ノードの使用は試みたが、24 ノードより多いノードでの解析に課題があり、実施できていない。解析モデルによって一部ノードのメモリ使用量が異常増加する現象が原因であり、今後のソフ
トウェアバージョンアップ等での改善が望まれる。
また、3D スキャナを活用したモデル化に関しては具現化できたものの、レーダ、バンパー、ボディの相図 3(b) 各条件における計算時間

対位置関係の合わせ込み等に課題も見られた。

次に実際の車両搭載を想定した結果を示す。図 4に図 1 ①②に加え、図 2 ③の車両構造装荷モデルのシミュレーション結果として、レーダ水平方向の放射特性を示す。この結果からは、①の基本状態から
従来②簡易板モデルで発生した細かなリップルは③車両構造モデルにおいて見られないものの、一部角度にて大きな放射特性の落ち込みの発生が確認される。これらから、細かなリップルは簡易板が 2 次元的（上下方向で変化がない）形状となっており端部条件が揃うことで発生していた特殊なものであり、実際
図 4 放射特性シミュレーションの車両構造では発生の程度が小さいことが想定される。一方で実際の車両搭載において、特定角度での大幅な利得劣化、検知性能劣化の存在は示唆さ車体側れるため、注意深い搭載条件調整が必要とされる。
さらに図 5 に図 2 ③におけるレーダ断面の電界分布を示す。図 5 からは、バンパーとボディの間の広範な領域に電磁界成分の広がりが確認できることか

電波放射側

ら、車両構造の影響をレーダ電波環境の設計の考
図 5 電界分布 車両構造モデル

慮にいれる必要性が考えられる。これは昨年度トラ

23

（様式 20）成果報告書

まとめ、今後の課題

TSUBAME を用い、準ミリ波レーダの自動車搭載を模擬した大規模電磁界シミュレーションを行った。
今回は新たにさらなる規模拡大への対応として計算の効率化の妥当性確認と、実際の車両構造に即したモデル化スキームの開発と車両搭載を想定したシミュレーションを実施することができた。これらから、より実物の車両状況に近い高度なシミュレーションの実現性を示すことができた。
なお、今回は大規模性として昨年度トライアルに比べ 1.5 倍オーダ規模 32 億メッシュの解析を実現できたが、さらなる大規模化実現には多ノードの安定活用が不可欠であろうと考えられる。また TSUBAME
自体の今後のバージョンアップによるポテンシャルの伸長にも期待したい。
参考文献[1] 平成 24 年中の交通事故発生状況，警視庁[2] Adv. Radio Sci., 7, 61–65, 2009 79 GHz UWB
automotive short range radar – Spectrum
allocation and technology trends

24（様式 20）成果報告書

平成 26 年度 TSUBAME 産業利用トライアルユース 成果報告書

利用課題名 大容量の行動関連データを使った最適シミュレーション及び暗黙知の利用
英文：Tacit knowledge and use of optimal simulation using large amounts of data related behavior
佐々木 武
Takeshi Sasaki
所属
株式会社 電通国際情報サービス
Information Service International-Dentsu,Ltd.
URL http://www.isid.co.jp

我々の課題は、膨大なデータが日々蓄積される今日の消費者生活において、広告業務の媒体計画作成に携わる人向けの業務システムを提供しており、顧客が求める期日迄に媒体計画を策定する機能を実現することが困難なケースがあるということである。昨年度、シンプルな課題において高速化が図れたが、より複雑な条件における高速化を目指し、CPU レベルの並列化、GPU レベルの並列化をチューニング検証することが今回の目的である

Keywords:

Media-Mix Optimized Media Planning Unit- Vehicle

背景と目的

概要

広告キャンペーンにおいて、様々な条件 例えば、「地媒体計画立案のための条件を下記のようにした。

区」、「媒体種類」、「媒体予算」、「媒体出稿期間」、「最適化

広告ターゲット＝「４歳以上男性または女性」

指標」などがあり、その条件から媒体への接触ローデータを

広告出稿エリア＝「関東エリア」

元に、ビジネスロジックの実装をコンピュターリソースを用い

対象媒体 A 案＝「TV」、「ラジオ」て最適化された媒体計画を算出するアプリケーションを

対象媒体 B 案＝「TV」、「ラジオ」、「雑誌」、「新聞」、

我々は提供している。昨今の広告媒体の状況は取り分け「屋外・交通」、「Internet」

Internet 媒体の種類の多さ、ビジネス・ルールの複雑さな

各媒体は スポット広告も番組に見立てて集計

どがめまぐるしい速さで変わっており、媒体組合せが膨大で

出稿予算 A 案＝「１億円」（各媒体 10%刻みで構成）あり、日々増大している。複数の媒体及びその個々メニュー



出稿予算 B 案＝「１億円」（各媒体 10%刻みで構成）のうち、目標に対し最適化された媒体計画を出力するため

最適化指標＝広告接触回数「３回」に掛ける時間が膨大になりつつある。集計時間を待つ非効

最適化目標＝最適化指標到達率「３５％」率性と機械損失を含めた



予算許容指数＝１．２（倍）これらの条件を元に、最適化された媒体計画を算出する。
本プロジェクトでは、昨年度の TSUBAME 上での実行時
間「４時間５８分」に対し半分の時間「２時間３０分」を目標に、２倍の処理速度を得るために 既存のボトルネック部分
CPU の並列化を行い、GPU 演算に回る部分も大きな分岐の単位で並列化をめざす。

25

（様式 20）成果報告書
結果および考察◆新たな処理概要は下記に処理概要フローを示す。

◆現行の処理概要は下記に処理概要フローを示す。

（A 案）

（A 案）

●CPU 並列処理導入●CPU 処理の中に GPU 処理への変更【広告ターゲット（標的）処理】

条
件

サンプル
属性マスタ

【TV媒体処理】

条
件

標的：該当
サンプル

番組候補

TV

番組ｘ個人
広告効果

C
P
U
並
列
処
理

C
P
U
1

統合処理

TV番組ｘ個人
広告効果

G
P
U
1

【ラジオ媒体処理】
番組候補
条
件

ラジオ

番組ｘ個人
広告効果

C
P
U
２

統合処理

TV番組ｘ個人
広告効果

G
P
U
２

共通の処理は、従来型の GPGPU の処理内容を踏襲。
GPGPU の（B 案）の場合（A 案）においては、「TV」「ラジオ」の組合せが２媒体で固定されていたので、最終的な GPGPU の処理は、１つの組合せブロックのまま処理するだけで済む。

（B 案）の場合においては、６媒体の並列処理を行う。
この場合、組合せが少ない媒体新聞などは、他媒体より先に終了するが、全媒体の組合せ最適を行うには全部が終了することが必要になる。そのため、管理機能を新たに構築して、「媒体組合せ」、「予算比の組合せ」、「GPU 集計の制御」、「共通部分の GPU 結果共有」、「GPU 結果の保存」の機能を管理機能側が制御し、GPU の演算分と切り分ける。

26

（様式 20）成果報告書
Step5: 「従系」媒体および予算の組合わせを CPU の管理機能が組合せを作成し、「GPU に処理依頼」を実行◆管理機構概要する。その際、「主系」の結果を共有し無駄な演算をせずに済むようにする。
Step6: 「従系」媒体の集計が終了した場合、結果を管理機構に戻す。
Step7: 管理機構は「従系」で「未集計」の組合せを、新たに GPU に割り当て集計処理を実施する
Step8： 全ての割り当ての集計が終了したら、集計結果
をファイルに出力し終了する。
【各 GPU での処理】
GPGPU で処理する主要処理の概念図。

※ 複数媒体の媒体組合せは集計開始時に条件指定できる。
Step1: 全ての媒体において、単媒体の処理を集計する。
Step2: 単媒体処理が終わったもので、２媒体の組合せが可能なものを「主系」とし、アサイン可能な CPU に対して処理を指示する。
Step3： 「主系」の「予算比率」分の組合せ演算を「GPU」
で集計する。
Step4: その結果を共通で利用するため一度 CPU 側のメモリに保存する

27

（様式 20）成果報告書の効率が落ちている。この部分は、複数媒体の組合せ集計

実測度の計測：１） MPI による CPU 側処理の並列化の影響調査において「管理機能」が存在するオーバーヘッドが、若干発（同一条件での試行回数 5 回の平均値）

生している。

【A 案】２媒体
No. サーバ様式
1 NormalLinux（NonMPI）
2 DualCPU Linux（MPI)
3 DualCPU Linux（MPI)【Tuning】

経過時間
354 分 （5時間54分）
312 分 （5時間12分）
183 分 （3時間3分）

CPU 処理の高速化を図り、MPU 処理の昨年比約「1.70

まとめ、今後の課題

今期の目情である２媒体で「２時間 30 分」以内は、実現できた。課題については、（B 案）の場合において、６媒体の並列処理で、組合せが少ない媒体新聞などは、他媒体よ倍」の高速化を実現できた。従来型（一昨年前）の MPI 無し

りの場合と比較した場合には、「1.93 倍」の高速化が図れた。

先に終了するが、全媒体の組合せ最適を行うには全部が終了することが必要になる。結果一部の GPU のリソースが【B 案】６媒体「待ち」になっている。

これ以上の高速化を図る場合、一番長く処理時間が掛る媒体を更に並列処理をすることが必要である。今回の試みで媒体数を増やした【B 案】の場合においても、今回のはその部分まで踏むこむことが出来なかった。１媒体あたり

Tuning 版・MPU における高速処理の比率は【A 案】とほぼの処理を内部で細分化して、共通機構側に進捗情報を共同等であった。（Tuning 無し）MPI 比「1.71」倍。MPI 無し・有する仕組みが新たに必要でありこの部分の改善が今後

比「1.93 倍」という数字であった。

の高速化を図るチューニング・ポイントとなる。

６媒体の組合せにおいて、媒体の特性などがあり新聞や雑誌などは組合せが限定的にしか発生しないため、TV などに比べて合計処理時間が短いため、単純に媒体数の比（２：６）である「３倍」という時間にはなっていない。
（NormalLinux NonMPI で「１．９倍」）

２） GPU 並列化の影響調査（同一条件での試行回数 5 回の平均値）【A 案】２媒体

MPI+Tuning による高速化係数 ＝ 「2.17」倍

【B 案】６媒体

MPI+Tuning による高速化係数 ＝ 「2.1」倍

媒体数を増やした「B 案」においても、２媒体の「A 案」とほぼ同じで、MPI＋Tuning を行った方が約「2.1」倍高速化が実現できた。ただし、「B 案」の方が「0.07」倍ほど、高速化

28

（様式 20）成果報告書

平成 26 年度 TSUBAME 産業利用トライアルユース 成果報告書

利用課題名 減衰を考慮した高周波数領域までの音響構造連成シミュレーション大規模化技術の検討
英文：Study of large scale simulation technology of acoustic-structural couplingwith dumping at a high frequency range
利用課題責任者 笹島 学
Manabu Sasajima
所属 フォスター電機株式会社
Foster Electric Co., Ltd.
URL http://www.foster.co.jp
邦文抄録
小型の音響機器は、狭い音響伝達経路での空気の粘性による減衰効果の影響が相対的に大きくなる。しかし、一般的な音響解析手法はこの影響が少ないものとして、無視して行われることが多い。また、解析対象とする周波数を人の可聴周波数帯域とした場合、特に高い周波数での振動を表現するには、波長に依存した非常に小さなサイズの要素分割が必要となる。しかし、要素数が膨大となり、計算規模が大幅に増加してしまう問題が生じるため実用的な解析は困難であった。本件は空気粘性を考慮した、可聴域の高周波数帯域までの音響構造連成シミュレーション
実用化のための大規模解析手法を検討する。
英文抄録
In the case of small audio equipment, effect of damping of the attenuation material and the shape of
the acoustic transmission path effect are larger than the case of general audio equipment. However,acoustic analysis is performed without regard to this effect in general. The audible frequency band of
human frequency to be analyzed, the element size must be very small size that depends on the
wavelength. However, analysis became difficult by large number of elements. This study is to examine
the large-scale analysis for structural acoustic coupled simulation that considers the air viscosity.
Keywords: 音響機器、有限要素法、音響構造連成解析、空気粘性

背景と目的の最大サイズを約 3mm 以下にする必要がある。さらに、イヤホンなどの小型の音響機器では、音響伝達経路

空気要素の節点変位を未知数として定式化した場合は、の狭さによる減衰効果が音響特性に与える影響が他

スムーズに節点を変位させるため、さらに領域の要素の音響問題と比べて相対的に大きい。しかし、一般的分割数を増やし、要素サイズを小さくする必要がある。

によく行われている音響解析手法は空気の粘性による

しかし、解析モデルの要素数が膨大となり、計算規模影響が少ないものとして、この影響を無視して行われるが大幅に増加してしまう問題が生じ、一般的なワークス

ことが多い。よって、一般的な音響解析手法を小型のテーションでは実用的な解析は困難であることがわか

音響機器に適用すると、実験結果とは大きく異なったシってきた。さらに、スーパーコンピュータ TSUBAME を用

ミュレーション結果となってしまう。また、特に閉じた狭いても、その解析時間から、所望の解析を行うことは難い空間では音源となる振動板の振動が周囲の空気のしいことがわかった。本件ではスーパーコンピュータ

圧力変動による影響を受けることもわかっているが、一

TSUBAME を用いて空気の粘性を考慮した、可聴域の般的な音響解析手法は空気の圧力である音圧を未知

高周波数帯域までの音響構造連成大規模解析のシミ数として定式化されるため、音源となる振動板の振動を、

ュレーション実用化のための手法を検討する。

音響解析と連成させて解析することは簡単ではなかった。一方、解析対象とする周波数の範囲を、人の可聴概要

周波数帯域とされている 20Hz から 20kHz とした場合に（１）音響構造連成有限要素の定式化

は、特に高い周波数で振動の現象を表現するためには、

直交デカルト座標系と、四面体定ひずみ要素を考え波長に依存した要素サイズが必要となることから、要素る。要素内の任意の点でのｘ方向変位を ux、ｙ方向変位29

（様式 20）成果報告書~

を uy、ｚ方向変位を uz とすると、ひずみエネルギU は以（２）プログラムについて

下のように表せる。
プログラムの開発言語は FORTRAN、コンパイラは∂u y ∂u z  ∂u~ 1 dxdydz ･･･(1)U = E ∫∫∫  x +
+
e
2
∂y
∂z 
 ∂x
．
ここで、E は体積弾性率である。変位の時間微分を u と~表現すると、運動エネルギ T は以下のように表せる。
2

T~ 1
T = ∫∫∫ ρ {u&} {u&}dxdydze2Intel コンパイラを用いている。このプログラムを用い、
TSUBAME 上で約１００万要素の構造音響連成モデルを解析してみたところ 1 周波数あたり約９０時間かかる
ことがわかった。全可聴領域の周波数特性を把握する･･･(2)
ためには少なくとも数十回の解析が必要になることを考えると、一つの周波数特性を得るために１から２カ月

ここで、ρは要素密度、T は転置行列を表す。また、粘~かかることになり、解析時間が実用上の大きな問題と性エネルギ D は以下のように表せる。

{ }

1
~
D = ∫∫∫ ΤT Γdxdydz
e2

なっていた。そこで、行列ソルバの改良により解析時間･･･(3)の短縮をはかった。

{T }は粘性による応力ベクトル、Γはひずみベクトルで

目標は従来プログラムを用いた場合の解析時間の十分の一以下とし、実用化の数値目標としては、100 万~ある。さらに、ポテンシャルエネルギ V は以下のように要素程度のモデルを解析する場合の解析時間を、1 周
波数あたり「30 分以下」とした。この場合、20 周波数分

表せる。

~
T
T
V = ∫ Γ {u} {P }dΓ + ∫∫∫ {u} {F }dxdydz ･･(4)の計算を約 10 時間でできるため、実用的な計算時間でeあると考えた。プログラムの計算プロセスごとの時間を{ }調査したところ、解析時間の大半はスカイライン法によ

ここで、 P は表面力ベクトル、｛Ｆ｝は物体力ベクトル、

る複素連立１次方程式を解くプロセスに費やされている∫ dΓ は要素境界での積分を表す。次に、粘性減衰をΓ

ことが分った。そこで、プログラムの行列ソルバを効率の良いものに入れ替え、解析時間全体の短縮を試み考慮した音響解析モデルとして要素運動方程式の定式~た。

化を考える。全エネルギ E は次式となる。

~ ~ ~ ~ ~E = U + D − T −V･･･(5)

さらに、以下のラグランジュの方程式を用いると、要素離散化方程式が求められる。

~
~
~
~
~
d ∂T
∂T
∂U ∂V
∂D
−
+
−
+
=0
dt ∂u& ei ∂u ei ∂u ei ∂uei ∂uei･･･(6)

ここで、 {u e } 、 {u& e } は節点変位ベクトル、節点速度ベクトルである。式(1)～式(4)を用いて式(6)を整理すると、165kmodel 解析モデル

次の要素離散化方程式を得る。ただし、変位は角周波数ωの周期応答としている。
165kmodel 解析結果例− ω 2 [M e ]{u e }+ [K e ]{u e }+ jω [C e ]{u e } = { f e } ･･･(7)
図１ 使用した四面体要素と解析モデル、解析例
ここで、 [M e ] 、 [K e ] 、 [Ce ] はそれぞれ要素質量行列、要素剛性行列、要素減衰行列である。また、 { f e }今回は、直接法で比較的取り扱いが容易ということ
で、インテルのマスカーネルライブラリ（MKL）から直接は節点力ベクトルである。この方程式を変位につい法の疎行列連立一次方程式ソルバである PARDISO をて解けば、全ての節点変位が求められる。

30

（様式 20）成果報告書導入した。PARDISO は、大規模な疎行列対称連立線性を考慮した可聴域の高周波数帯域までの音響構造

形方程式を解くことができ、高速で、メモリ効率がよく、連成大規模解析のシミュレーション実用化のため、解安定しているソフトウェアとされている。また、複素連立析の効率化手法を検討した。

線形方程式も解くことが可能である。

結果および考察
図１に今回使用した四面体要素、および例として、約165,000 要素の解析モデルと解析結果を示す。また、解析に用いた全てのモデル名称と概略の要素数を表１にまとめる。各モデルの要素は全て四面体で、節点はｘ、ｙ、ｚ方向の３自由度を持っているが、境界条件や形状等は統一していないので、計算の条件としては完全に同一にはなっていない。
図２に行列ソルバにスカイライン法を用いた従来プロ
グラムと、行列ソルバを PARDISO に置き換えたプログ

図２ 各モデルの解析時間比較

ラムを用いた解析について、各解析モデルの解析時間を示す。縦軸の解析時間は対数目盛りとしてある 。
PARDISO を用いることによって、どのモデルにおいても解析時間は約１００分の１となっていることが分る。可聴領域全域の周波数特性を把握するためには少なくとも数十回の解析を繰り返すことを考えると、実用上の効
果は非常に大きいと思われる。たとえば、これまで 1 カ
月必要であった周波数特性を得るための計算が半日でできるようになった。
また、この計算時のメモリ使用量をまとめ図３に示す。
この図より、PARDISO を用いることによって、どのモデルにおいても使用メモリ量は約 4 分の１となっているこ

図３ 各モデルの使用メモリ量比較とが分る。従来プログラムと比べ大幅に使用メモリの削減ができていることが確認できた。

従来のプログラムでは、解析時間が実用上の大きな問題となっていた。そこで、行列ソルバの改良により解表 1 テストモデル名および概略要素数

析時間の短縮を図った。今回は疎行列連立一次方程式の直接法ソルバである PARDISO を導入し、いくつか

モデル名

概略要素数（個）50K model

50,000165K model

1,650,000

600K model

600,000

1M model1,000,000のモデルを計算してみたところ解析時間が約１００分の
１、使用メモリ量が約４分の１になった。よって、本検討
での目標は十分に達成できた。
今後、従来解析が困難であった全可聴帯域での、高精度の音響解析シミュレーション技術の実用化が進むと考えられる。引き続き、PARDISO の解析規模の上限
まとめ、今後の課題を確認するとともに、場合によってはさらなる大規模化

スーパーコンピュータ TSUBAME を用いた空気の粘のため、反復法の行列ソルバも試していきたい。

31

32

（様式 20）成果報告書

平成 26 年度 TSUBAME 産業利用トライアルユース 成果報告書
利用課題名 広域都市環境の大規模計算による検討
英文：Large scale simulation for the global urban environment
利用課題責任者 佐々木 澄
Sasaki Kiyoshi
所属： 清水建設（株）技術研究所
Affiliation： Institute of Technology, Shimizu Corporation
URL: http://www.shimz.co.jp
本利用課題では，「商用アプリバンドル型トライアルユース」に提供された商用アプリを活用して，TSUBAME ス
パコンの計算資源を利用するとともに，広域都市環境評価における大規模計算検討を行った．また，プリ・ポストの自動処理環境を構築し，広域市街地や複雑形状市街地を対象として最大 2 億メッシュ規模の都市環境解析を実施した．本商用アプリにおいてその大規模計算の性能特性を確認し，システムの実用性を明らかにした．
A large scale simulation for a global urban environment is investigated using commercial software
which is provided on TSUBAME2.5 supercomputer. An automatic Pre-post processing system has
been developed to handle the large-scale model instead of using the GUI system of the software.
Urban environment analysis of 200M meshes has also carried out to clarify the performance of thesoftware in the computation of such large-scale model.

Keywords: Urban Environment, Large-Scale Simulation, Building, Commercial Software

背景と目的と相当）を対象として，商用アプリのプリ・ポストよりモデ

新たに高層建物，あるいは大規模構造物を建てる際ルの作成と可視化までの一連処理を確認していた．

には，風の状況変化に従う影響は無視できない場合が本利用課題では，数千万また数億計算格子規模の多く，風環境や風騒音などの問題が生じる場合もある．

風環境解析を対象として，モデルの作成や計算結果のこうした建設による問題や障害の発生を未然に防ぐため

可視化処理とその計算の安定性等，商用アプリの検討には，風環境の変化や風騒音の発生状況を風洞実験

課題として行い，必要となる環境を整備するとともに，本や数値流体計算などの方法により予測し，事前に調査・商用アプリの大規模計算環境を確立する．

検討および対策の立案を行う必要がある．都市部では，高層ビルが密集し，周辺から受ける影響が大きくなるが，概要

実験装置の制約などにより風洞実験による予測は限界「商用アプリバンドル型トライアルユース」に提供されがあり，広域を考慮できる大規模数値流体計算の実施

た商用アプリは，風・温熱環境や騒音などの建築環境が必要となる．

問題の評価に多くの実績を持っているアプリケーション

前年度では，TSUBAME スパコンの「商用アプリバンとして知られているため，本利用課題はその応用検討を

ドル型トライアルユース」に提供された商用アプリを活用

重視せず，TSUBAME2.5 計算資源の利用による大規

し，利用計算環境の整備と，実中層市街地の建物群に模計算検討と大規模データハンドリングおよび一連自対象とした風環境解析を実施した．また，数千 CPU コア

動計算環境の構築を行った．ここで，シミュレーションの規模より本アプリの並列性能を調べており，アプリのシミ部分について商用アプリのコア部分のバッチ処理を活ュレーションのコア部分に対して大規模並列処理とその用して，プリ・ポスト部分ついて独自な環境を整備した．

計算時間の短縮化を確認していた．なお，本検討にはまた，5 千万および２億計算メッシュ規模の広域市街地

1 千万格子規模（現状の建築分野における大規模計算などの風環境解析を実施し，構築された計算環境の有33

（様式 20）成果報告書
2. 計算モデルと結果

効性を明らかにした．

本課題で構築した自動処理環境の有効性を確認する
結果および考察ためには，広域市街地また複雑形状市街地を対象とし

1. データハンドリングて，計算モデルの作成と可視化およびシミュレーション

GUI は視認性に優れ，直感的な操作ができるため，を行った．図-2 には 4km×4km 広域市街地の計算領域，

利便性の観点で多くの商用アプリが整備されている．し計算メッシュ，および計算から得られた地上 50m の風速

かしながら，数億計算メッシュ規模の扱いや大規模処理と建物壁面風圧コンターの拡大図を示す．計算メッシュを行う際には，現状の計算機の性能よりその生産性がは非構造格子であり，その全体は約 5 千万メッシュであ著しく低下し，本来のアプリの性能を十分に発揮できるる．風上側の高層建物群から発生されたウェークは，周ケースがある．特に，TSUBAME の計算環境ではデー辺に与える影響も確認できる．

タハンドリングや計算処理などが遠隔操作より行われているため，その操作性が不適である．
2km

本利用課題では，その大規模処理を行うため，図-1に示すように既存の商用アプリに提供された処理プロセ
スに対して新たな一連の自動処理環境を構築していた．
まずは，大規模モデルの作成について独自のプリ自動処理環境を構築して，メッシュデータの変換などにより本商用アプリの必要な計算モデルを作成した．また，計

4km
4km

算部分について商用アプリのバッチ処理より行われており，得られた計算結果はデータの変換より独自の自動ポスト処理より可視化される．

(a) 計算領域

これよりすべての作業は自動的に行われており，大規模モデルの取扱いに対してモデルの作成から解析
及び可視化の一連のプロセスはスムーズに行われる．
・計算モデル作成

GUI処理

・計算条件の設定
・シミュレーション

GUI・専用バッチ処理

・画像処理・可視化

GUI処理(b) 計算メッシュの一部拡大図(a) 商用アプリの既存処理方法
・計算モデル作成
・データ変換オープンソース・独自コード
バッチ処理

・計算条件の設定
・シミュレーション
・計算データの変換

対象商用アプリのバッチ処理（汎用フォーマット）

・画像処理・可視化

オープンソース・独自コード
バッチ処理(b) 本課題の大規模計算処理方法(c) 風速・風圧コンターの一部拡大図

図-1 計算データハンドリング

図-2 広域市街地

34

（様式 20）成果報告書
大規模計算を実現できており，本商用アプリのシミュレ
ーションコアの有効性を確認できた．今後は計算環境等を応じて本商用アプリを積極的に活用していきたいと考えている．

(a) 計算領域(b) 計算メッシュとその拡大図(c) 風速・風圧コンター
図-3 複雑形状市街地

一方，図-3 には 1.6km×1.6km 市街地であるが，建物形状が複雑であり，その形状を再現するために非常に高解像度の計算メッシュを利用する必要がある．図-3(b)は計算メッシュとその一部拡大図を示す．建物形状の他，ベランダなども細かく再現されることが分かる．全体は約 2 億計算メッシュである．図-3(c)は計算から得られた地上 5m の風速と建物壁面風圧コンターとその拡大図を示す．建物形状の影響や周辺による影響等は確認できる．

まとめ，今後の課題
本利用課題では，「商用アプリバンドル型トライアルユース」に提供された商用アプリの活用および大規模計
算検討を行った．TSUBAME 計算環境では本商用アプ
リ GUI よりモデル作成・可視化等，遠隔操作で大規模データハンドリングとしては困難であることを確認した．
また，本商用アプリ GUI 依存しない新たな独自のプリ・ポスト処理環境を構築し，データフォーマットの変化作
業と自動バッチ処理の活用により数億メッシュレベルの

35

36

平成 26 年度 TSUBAME 産業利用トライアルユース 成果報告書

利用課題名 三次元電磁界シミュレータを用いた静電気放電イミュニティ試験に於けるPCB/Package/Chip のイミュニティ解析 （継続課題）
英文： Immunity analysis of PCB/Package/Chip in the electrostatic discharge test

using the three-dimensional electromagnetic field simulator
賀 文寿、永野 民雄
Fumikazu Ga, Tamio Nagano
所属: ルネサスシステムデザイン株式会社
Affiliation: Renesas System Design Co., Ltd.
URL: http://www.renesas.com/
静電気放電（ESD）イミュニティの実験環境を含む大規模モデルに対して、3 次元電磁界シミュレータによる ESD
イミュニティ解析を検討した。シミュレーションと実測の相関性を確認し、実用的な時間範囲内で解析できることを確認した。スーパーコンピューターTSUBAME による大規模 3 次元電磁界解析は、ESD イミュニティ設計に有用であることを示した。
Electrostatic discharge (ESD) immunity analysis using 3-dimensional (3D) electromagnetic field
simulator for large-scale models, including ESD immunity test environment, have been studied. The
correlation between simulations and measurements is confirmed, and the analysis can be finished withina practical time. The large-scale 3D electromagnetic field analysis by the TSUBAME supercomputer isdemonstrated to be useful in designing ESD immunity.

Keywords: 電 磁 界 シ ミ ュ レ ー シ ョ ン Electromagnetic field simulation 、 静 電 気 放 電 Electrostatic
discharge（ESD)、イミュニティ Immunity、プリント回路基板 Printed Circuit Board (PCB)、
大規模集積回路 Large Scale Integration (LSI)TSUBAME を利用し、ESD イミュニティ試験の 3 次元

背景と目的
電子機器の低電圧化、高速化に伴い、静電気が起電磁界解析を行い、ノイズの伝搬経路の可視化を含む

因の誤動作が発生しやすくなってきている。静電気放ESD イミュニティの設計手法を構築することを目的とし電（ESD）による破壊や誤動作を防止するため、ESDた。

イミュニティ試験[1]が規格化され、製品に対する ESD

平成 25 年度に同課題名で、3 次元電磁界解析によイミュニティ試験が行われている。

るノイズ伝搬のメカニズムについて検討を行った。しか低コストで ESD イミュニティに強い PCB/Package/

し、実機試験結果とシミュレーションとの比較考察にはChip を設計するには、ESD イミュニティ試験でノイズ

至らず、LSI モデルや、解析条件などについて課題ががどの経路をたどってチップに伝搬し、誤動作を引き起残る結果となった[4]。

こしているかを理解することが重要である。一つのアプ

本検討では、平成 25 年度の課題に対して、 LSI モ

ローチとして、3 次元電磁界シミュレーションによるノイデル、ケーブルモデル、励起源などをより実体に合せてズの可視化がある[2]、[3]。

3 次元電磁界解析を実施し、試験と解析結果は傾向一

一般の PC 環境では、基準グラウンド面、水平結合

致という結果を得た。また、ノイズの伝搬経路の可視化

板、ESD ガン、PCB/Package/Chip 等を含む ESD イを実施し、実用的な時間範囲内で解析できることを確

ミュニティ試験の実環境を 3 次元電磁界シミュレータで認した。

解析することは、使用メモリ、解析時間の観点から困難

前回と今回の検討を通じて、TSUBAME による大規
模 3 次元電磁界解析は、ESD イミュニティ設計に有用である。
本プロジェクトでは、スーパーコンピューター

であることを示した。

37

概要

そこで、下記のモデル改良を行い、3 次元電磁界解

CST STUDIO SUITE 2014 を用いて[5]、ESD

析を実施した。

イミュニティ解析の検討を行った。

－LSI モデルの電源構造を改良する

シミュレーションモデルを図 1 に示す。実試験環境に－ケーブルモデルの長さ、設置状態をより実体に合わせて、基準グラウンド面、水平結合板、ケーブル、合わせる

ESD ガン[6]、絶縁シート、PCB、LSI 等をモデル化し－ESD ガンの出力は実測に合うように、励起信号をた。図 1 のモデルを用いて 3D 電磁界解析を行い、試調整する
図 2(b)に改良モデルのノイズ波形を示す。LSI に伝験との比較、ノイズの伝搬経路の可視化を行った。

搬するノイズは電源分離ボードの方（赤線）が大きい。
試験と解析結果は傾向一致という結果が得られた。
ノイズの伝搬経路を可視化した表面電流分布を図 3に示す。TSUBAME を使用して、図 1 に示すモデル
全体の表面電流分布が得られた。しかし、表示用 PC
マシンの性能制限で、局部の表面電流分布しか表示できず、課題があった。
(a) 全体赤：電源分離ボード
青：電源共通ボード(a) 前年度モデルのノイズ波形（b） 評価ボード周辺
図 1 シミュレーションモデル
赤：電源分離ボード結果および考察

青：電源共通ボード

PCB 上電源配線パターンが共通と分離の 2 種類の評価ボード（以降、電源共通ボードまたは電源分離ボ(b) 改良後モデルのノイズ波形ードと呼ぶ）に対して、IEC61000-4-2 に準拠した試験

図 2 シミュレーション結果

環境で、ESD ガンによる接触放電試験を実施した。
試験結果により、LSI の誤動作耐圧は、電源共通ボードより電源分離ボードの方が低い。よって、LSI に伝搬するノイズは、電源分離ボードの方が大きいと考えられる。
図 2 に Chip の電源端子でのノイズ波形のシミュレ
ーション結果を示す。図 2(a)に前年度モデルのノイズ波形を示すように、LSI に伝搬するノイズは電源共通
ボードの方（青線）が大きい。この結果は試験傾向と異図 3 表面電流分布

なる。

38

今回のシミュレーション時間を表 1 に示す。約 200M

ノイズの可視化の検討を行った。しかし、大規模モデ

メッシュ規模のモデルに対して、TSUBAME の 4 ノールのノイズ可視化データに対して、表示用 PC マシンのドでの計算時間は約 5 h であり、実用的な時間範囲内

性能が追いつかず、全体表示困難の課題があった。

で解析できることを確認した。

今後、メッシュ数を削減して、計算の高速化、表示のPC と TSUBAME の計算時間の比較を図 4 に示容易化を図る。

す。一般的な PC 環境で、基準グラウンド面、水平結合
板、ESD ガン、PCB/Package/Chip 等を含む ESD イ
ミュニティ試験の実環境を解析すると、推定の計算時

謝辞

間は 1200 h 以上となる。一方、TSUBAME での計算

本プロジェクトにおきまして、多くのご支援を頂いた東

時間は約 5 h で、240 倍以上の時間短縮効果が得ら

京工業大学学術国際情報センター共同利用推進室、れた。

CST STUDIO SUITE ソフトウェアのサポートを頂いた
表 1. シミュレーション時間

株式会社エーイーティーに深く感謝いたします。

メッシュ数

使用ノード数*

計算時間

200M

4

5h
参考文献

*1 ノード：12CPU コア、3GPU[1] IEC61000-4-2 ， Testing and measurement
techniques - Electrostatic discharge immunity
test，Edition 2.0，2008.12[2] http://www.fujitsu.com/jp/solutions/businesstechnology/tc/fields/cae/poynting/esd-pcbchassis.html[3] 賀 文寿, 秋山 雪治, 白井 淳一, “電磁界シミュレーションによる ESD イミュニティの解析”2013
年第 27 回エレクトロニクス実装学会春季講演大会 14C-11[4] 秋本 哲也, “三次元電磁界シミュレータを用いた

図 4 PC と TSUBAME の計算時間の比較

静電気放電イミュニティ試験に於けるPCB/Package/hip のイミュニティ解析” 平成 26

まとめ、今後の課題

年度 東京工業大学 TSUBAME 産業利用シン

TSUBAME を利用して、ESD イミュニティの実験環

ポジウム pp.147~150

境を含む大規模モデルに対して、3 次元電磁界解析を[5] https://www.cst.com/

実施した。

[6] S.

LSI モデル、ケーブルモデル、励起源などをより実体

Caniggia,

F.

Maradei,

Numerical

析結果は傾向一致の結果を得た。また、200M メッシュ

Discharge Generators," IEEE Transactions on

規模のモデルに対して、TSUBAME の 4 ノードでの計

Industry Applications, Vol. 42, pp. 1350-1357,

算時間は約 5 h であった。これらの結果から、

Nov. 2006.

解析の実用性を確認できた。本技術はノイズの可視化、
誤動作のメカニズム解明、イミュニティ設計・対策の応用に期待される。

39

of

and

に合せて 3 次元電磁界解析を実施した結果、試験と解TSUBAME による大規模モデルの ESD イミュニティ

Modeling"Circuit

Electrostatic

40

（様式 20）成果報告書

平成 26 年度 TSUBAME 産業利用トライアルユース 成果報告書

利用課題名 ワイヤレス電力伝送による漏えい電波の環境解析技術の研究開発
英文 Research and Development of Environment Analysis Technique of Leaked Electromagnetic Field from Wireless Power Transfer
利用課題責任者 池田 和彦
Kazuhiko Ikeda
所属 株式会社パナソニック システムネットワークス開発研究所

Panasonic System Networks R&D Lab. Co., Ltd.
http://panasonic.co.jp/avc/psnrd/
邦文抄録（300 字程度）
家庭用電子機器などの充電に用いられるワイヤレス電力伝送（WPT：Wireless Power Transfer）システム
が近年検討されている．WPT システムから漏えいする電磁界は他の様々な無線機器に影響を与える可能性があり，その干渉影響を定量的に評価する必要がある．WPT システムの設置が想定される戸建て住宅など
では，近傍に無線通信機器が存在し，それらに影響を及ぼす恐れのある周波数帯の漏えい電磁界評価が必要である．しかしながら，戸建て住宅全体の漏えい電磁界を測定するには多大な時間と労力がかかり，また電磁界シミュレーションでは空間メッシュ数が大規模になる課題がある．そこで本報告では，WPT シス
テムの高調波である 87 MHz，815 MHz，2.497 GHz の漏えい電界の戸建て住宅内特性を，スーパーコンピ
ュータ TSUBAME2.5 を用いて電磁界解析し，測定結果と比較検証する．
英文抄録（100 words 程度）
The wireless power transfer (WPT) system for electronics devices has been actively developed in recent years. Since
the electromagnetic field leaked from the WPT system may cause the false operation of other devices, it is required toquantitatively evaluate the electromagnetic interference. In the housing environment, the WPT system is locatedgenerally close to the wireless communication devices. Therefore, the quantitative evaluation of the leaked
electromagnetic field in various installation environments is indispensable. However, the measurement in the wholehouse involves immense amount of time and effort, and the electromagnetic simulation requires the enormous
analysis meshes because of the volume of the house. In this paper, we evaluate the leaked electric field at 87 MHz,815 MHz and 2.497 GHz in the house by electromagnetic simulation using supercomputer TSUBAME 2.5. The
comparison between the simulation and measurement results is discussed.
Keywords: Electromagnetic simulation, Wireless power transfer, Leaked electromagnetic field

1. まえがきら，戸建て住宅内全体の漏えい電磁界強度を測定するには近年，電気自動車や家庭用電子機器などの充電に用いら膨大な時間と労力がかかるため，電磁界シミュレーション

れるワイヤレス電力伝送（WPT：Wireless Power Transfer）

を用いて漏えい電磁界分布を解析することが望ましい．しシステムが検討されている[1],[2]．WPT システムから漏えたがって，WPT システムの設置が想定される戸建て住宅

いする電磁界は他の様々な無線機器の誤作動の要因になや商用施設における，広範囲の周波数領域の漏えい電磁界り得るため，漏えい電磁界の強度規制値が盛んに議論されの解析技術が必要である．戸建て住宅全体を含む電磁界シている．多種多様な WPT システムからの漏えい電磁界にミュレーションは大規模な空間メッシュでの解析となる

よって引き起こされる干渉問題が重要な課題となっておため，並列処理によって大規模メッシュの解析が可能なス

り，様々な設置環境における定量的な評価が WPT システ

ーパーコンピュータの利用が有効となる．

ムの実現には必要不可欠である．

本利用課題では，ラジオ受信機，携帯端末，無線 LAN

例えば，戸建て住宅では，WPT システムの近傍に様々ルータへの干渉を想定し，WPT システムの高調波であるな無線通信機器が存在し，それらに影響を及ぼす恐れのあ87 MHz，815 MHz，2.497 GHz の漏えい電界の戸建て住宅

る周波数帯の漏えい電磁界評価が必要である．しかしなが内特性をスーパーコンピュータ TSUBAME2.5[3]を用いて

41

（様式 20）成果報告書
電磁界解析する．戸建て住宅及び漏えい電界の波源をモデ

図 1 に本検討で用いる住宅モデルを示す．住宅モデルとリングして住宅モデル内の漏えい電界分布を解析し，住宅

して，パナソニック保有の戸建て実験住宅をモデリングし内における一部の場所での測定結果と比較することで，解ている．住宅モデルの大きさは 19×14×8 m である．

析モデルの妥当性を検証する．また，各周波数における使図 2 に波源モデルを示す．それぞれの周波数で高い電界用ノード数と解析時間の関係を把握し，電磁界シミュレー

強度を放射するために周波数ごとに波源モデルを設計しションにおける TSUBAME2.5 活用の有効性を検討する．

ており，87 MHz では直径 5 cm 巻き数 3T のスパイラルコ
イル、
815 MHz では 88.5×88.5 mm 素子のパッチアンテナ，

2. 解析モデルと解析条件

2.497 GHz では 28.9×28.9 mm 素子のパッチアンテナを波

表 1 に解析条件を示す．計算機システムはスーパーコン

源として用いる．

ピュータ TSUBAME2.5，電磁界解析には有限積分法を用

図 3 に受信機器モデルを示す．87 MHz では電界プロー

いた電磁界シミュレータである CST 社の MICROWAVE

ブ，815 MHz ではスマートフォン（端末下部のモノポール

STUDIO[4]を使用する．解析周波数は，ラジオ受信機，携アンテナ），2.497 GHz では無線 LAN ルータ（端末側部

帯端末，無線 LAN ルータへの干渉を想定し 87 MHz，815のダイポールアンテナ）を受信機器として用いる．

MHz，2.497 GHz としている．

図 4 に波源の配置位置及び観測区間の場所（＝受信機器の配置位置）を示す．波源は 1 階リビングの壁掛け TV へ

表 1 解析条件

の電力伝送を想定し，1 階リビング壁面の TV ラック内に計算機システムTSUBAME 2.5

電磁界解析

MICROWAVE STUDIO

解析周波数

87 MHz，815 MHz，2.497 GHz

87 MHz波源

受信機器

配置している．波源の配置向きは，図 1 及び図 2 の座標軸

Z

190mm

直径 5cm 巻き数 3T コイル

Y
X

120mm
88.5mm815 MHz

88.5×88.5 mm パッチアンテナ

2.497 GHz

28.9×28.9 mm パッチアンテナ

87 MHz

電界プローブ

815 MHzスマートフォン

2.497 GHz

無線 LAN ルータ

88.5
mm

（a） 87MHz（b） 815MHz

190mm

28.9mm
28.9
mm

120
mm

（c） 2.497GHz

図 2 解析モデル（波源）

住宅モデル

戸建て実験住宅モデル

波源の配置位置

1 階リビングの TV ラック内

受信機器の配置位置

1 階キッチン，2 階子供部屋

Z

Z
Y

X

電界
プローブ

Y

（a） 87MHzX（b） 815MHz（c） 2.497GHz

図 3 解析モデル（受信機器）
観測区間の始点観測区間の終点

8m

キッチンの
観測区間

14m

19m

子供部屋の
観測区間

波源

Y

図 1 解析モデル（戸建て実験住宅）（a） 1階（b） 2階Z

図 4 波源の配置位置と観測区間（受信機器の位置）

42X

（様式 20）成果報告書の通りであり，壁面に対してコイルまたはパッチアンテナ

電界プローブを 1 階キッチン及び 2 階子供部屋の観測区間が水平（＝地面に対して垂直）となる．また，様々な無線内を移動させた場合の電界強度の解析結果と測定結果を機器の設置が想定される 1 階キッチン及び 2 階子供部屋に記載している．電界強度は 3 軸合成値とした．図 6 より，

図 3 の受信機器を配置し，図 4 に示す観測区間内を移動さ

1 階キッチンにおいて解析結果と測定結果の電界強度のせた場合の電界強度または受信電力を観測する．

差分が大きいことがわかる．これは，初期住宅モデルに実験住宅との差異があるためと考えられ，住宅モデルの改良

3. 解析結果と測定結果の比較が必要である．

本章では，住宅内の漏えい電界について解析結果と測定

図 7 に住宅モデルの改良箇所を示す．初期住宅モデルと結果を比較検証し，解析モデルを構築する．解析周波数

実験住宅との差分として，以下の 8 項目を抽出した．

87 MHz，815 MHz，2.497 GHz の順で検討し，全ての周波数において解析結果と測定結果の差分が小さい解析モデ①

基礎のコンクリート内部の鉄筋を追加

ルの構築が目標である．

②

階層間の金属製の梁を追加③

バルコニー床の材質を変更（鋼板を模擬）④

階段手すりの金属製の支柱を追加⑤

窓ガラスの材質を変更（UV カットガラスを模擬）⑥

屋根の材質を変更（遮熱仕様鋼板を模擬）⑦

外壁の材質を変更（サイディング材を模擬）⑧

内壁と天井の材質を変更（石膏ボードを模擬）

3.1. 87 MHz の解析結果と測定結果の比較
87 MHz では，解析結果と測定結果の比較の指標として
電界プローブで観測した電界強度を用いる．
図 5 に初期住宅モデルを示す．初期住宅モデルは，住宅モデルの各構造に設定する材質を 4 種類のみとした簡易モデルであり，屋根・外壁・内壁の材質を木に，地面・石垣の材質をコンクリートに，玄関のドア・浴室・太陽光パ
ネル・屋外設備の材質を金属に，窓ガラスの材質をガラスに設定している．
図 6 に初期住宅モデルにおける電界強度の解析結果と① 基礎内部の鉄筋追加② 階層間の梁追加③ バルコニー床の金属化④ 階段手すりの支柱追加⑤ 窓ガラスの材質変更⑥ 屋根の金属化⑦ 外壁の材質変更⑧ 内壁・天井の材質変更

測定結果の比較を示す．波源の送信電力は 10 dBm とし，
屋根
材質：木

浴室
材質：金属
屋外設備
材質：金属

石垣
地面
玄関のドア 窓ガラス
材質：コンクリート 材質：コンクリート 材質：金属 材質：ガラス（a） 全景（b） 1階の断面図

図 5 初期住宅モデル
120

120

電界強度 [dBuV/m]

100

80
60
87MHz, 10dBm
初期住宅モデル1階キッチン

40

120

Meas.
Sim.

20

80
60
87MHz, 10dBm
初期住宅モデル2階子供部屋

40
20

0

1

2

3

120

Meas.
Sim.

100

電界強度 [dBuV/m]

Meas.
Sim.

100

電界強度 [dBuV/m]

図 7 住宅モデルの改良箇所

80
60
87MHz, 10dBm
改良住宅モデル1階キッチン

40
20

0

1

2

3

Meas.
Sim.

100

電界強度 [dBuV/m]

太陽光パネル
材質：金属

内壁
材質：木

外壁
材質：木

80
60
87MHz, 10dBm
改良住宅モデル2階子供部屋

40
20

0

1

2

3

0

1

2

移動距離 [m]移動距離 [m]移動距離 [m]移動距離 [m]（a） 1階キッチン（b） 2階子供部屋（a） 1階キッチン（b） 2階子供部屋

図 6 初期住宅モデルでの電界強度（87 MHz）

図 8 改良住宅モデルでの電界強度（87 MHz）

43

3

（様式 20）成果報告書
図 8 に①～⑧の項目を全て反映した改良住宅モデルに

移動区間中央値（区間長：4 波長）を算出してグラフに記おける電界強度の解析結果と測定結果の比較を示す．初期

載している．図 9 より，1 階キッチン，2 階子供部屋共に住宅モデルでは差分が大きかった 1 階キッチンで解析結解析結果と測定結果の電界強度がよく一致することが確

果の電界強度が増大し，測定結果との差分が低減すること認できる．観測区間における区間中央値の差分は 1 階キッが確認できる．改良住宅モデルを用いることで，観測区間

チンでは 0.1 dB，2 階子供部屋では 1.5 dB である．

における解析結果と測定結果の電界強度の区間平均値の3.3. 2.497 GHz の解析結果と測定結果の比較

差分を 1 階キッチンでは 2.2 dB，2 階子供部屋では 2.0 dB

2.497 GHz では，815 MHz の場合と同様に，解析結果とまで低減できている．また，①～⑧の項目をそれぞれ検証

測定結果の比較のための指標として伝搬損失を用いる．住した結果，①基礎内部の鉄筋追加，③バルコニー床の金属
宅モデルは 87 MHz 及び 815 MHz と同様の改良住宅モデ

化，⑦外壁の材質変更が 87 MHz の電界強度に与える影響

ルを使用する．

が特に大きいことを確認している．
図 10 に改良住宅モデルにおける伝搬損失の解析結果と3.2. 815 MHz の解析結果と測定結果の比較

測定結果の比較を示す．815 MHz の場合と同様，伝搬損失

815 MHz では，解析結果と測定結果の比較のための指標の移動区間中央値（区間長：4 波長）をグラフに記載してとして，受信機器（スマートフォン）で観測した受信電力

いる．電界強度の区間中央値の差分は 1 階キッチンではから波源の送信電力を減算することで算出した伝搬損失

2.3 dB，2 階子供部屋では 3.3 dB 程度である．2.497 GHzを用いる．住宅モデルとしては，87 MHz で構築した改良においても，87 MHz 及び 815 MHz と同様の改良住宅モデ

住宅モデルを使用する．

ルを用いることで解析結果と測定結果がよく一致してお図 9 に改良住宅モデルにおける伝搬損失の解析結果とり，本解析モデルは妥当であると考えられる．

測定結果の比較を示す．815 MHz では波長が観測区間より

4. 住宅内の漏えい電界分布の解析も短く，区間内にヌルが多数存在するため，伝搬損失の
-20

Meas.
Sim.

-40

伝搬損失 [dB]

伝搬損失 [dB]

周波数の住宅内における漏えい電界を解析する．

-20

Meas.
Sim.

-40
-60
-80

-100

本章では，3 章で構築した解析モデルを用いて，各解析

815MHz
改良住宅モデル1階キッチン

-80

-100

-120

大

-60

815MHz
改良住宅モデル2階子供部屋

-120
0.0

0.5

1.0

1.5

0.0

0.5

1.0

移動距離 [m]移動距離 [m]（a） 1階キッチン（b） 2階子供部屋

1.5

波源
小

図 9 改良住宅モデルでの伝搬損失（815 MHz）
-20

-20

Meas.
Sim.

Meas.
Sim.

-40

伝搬損失 [dB]

伝搬損失 [dB]

-40
-60
-80

-100

（a） 1階2.497GHz
改良住宅モデル1階キッチン

-60
-80

-100

-120

2.497GHz
改良住宅モデル2階子供部屋

-120
0.0

0.5

1.0

1.5

0.0

0.5

1.0

移動距離 [m]移動距離 [m]（a） 1階キッチン（b） 2階子供部屋

1.5（b） 2階図 11 住宅内における電界分布（87 MHz）

図 10 改良住宅モデルでの伝搬損失（2.497 GHz）

44

（様式 20）成果報告書
図 11，図 12，図 13 に 87 MHz，815 MHz，2.497 GHz

する解析モデルを構築し，住宅内における 87 MHz，815の住宅内電界分布の解析結果をそれぞれ示す．住宅モデル

MHz 及び 2.497 GHz の漏えい電界分布を解析した．

は改良住宅モデル，波源は表 1 及び図 2 で説明した通りで

TSUBAME2.5 の 56 ノード並列計算により約 52 億メッシあり，送信電力は各周波数で同様としている．図 11～図ュのモデルを約 10 時間で解析でき，電磁界シミュレーシ

13 の比較により，波源の違い（87 MHz：微小コイル，815

ョンにおける TSUBAME2.5 の有用性を確認した．

MHz 及び 2.497 GHz：パッチアンテナ）による電界分布のなお，本研究は総務省平成 26 年度電波資源拡大のため差や，周波数が高いほど漏えい電界の直進性が強まる傾向の研究開発「ワイヤレス電力伝送による漏えい電波の環境などが確認できる．

解析技術の研究開発」の一部である．

5. 解析モデルのメッシュ数と解析時間

参考文献[1] Q. Chen et al., “Antenna Characterization for WirelessPower Transmission System Using Near Field Coupling,”
IEEE Antennas and Propagation Magazine, Vol. 54, No. 4,
Aug. 2012.

表 2 に各解析周波数におけるメッシュ数，ノード数，解析時間の関係を示す．それぞれ改良住宅モデルでの値を記載している．本検討で用いた電磁界シミュレーションでは，大な解析時間が必要となる．このため，815 MHz では[2] J. Kim et al., “Electromagnetic Interference andRadiation from Wireless Power Transfer Systems,” 2014
IEEE International Symposium on Electromagnetic
Compatibility, Aug. 2014.

TSUBAME 2.5 による 16 ノード並列計算，2.497 GHz では[3] 東京工業大学 学術国際情報センター TSUBAME（http://www.gsic.titech.ac.jp/tsubame）

56 ノード並列計算を用い，解析時間を約 10 時間に抑えた．

[4] CST MICROWAVE STUDIO（http://www.cst.com/）

周波数に応じて解析メッシュ数が指数関数的に増大し，膨6. まとめ

表 2 メッシュ数・ノード数・解析時間の関係

本利用課題では，スーパーコンピュータ TSUBAME2.5

解析周波数

メッシュ数

ノード数

解析時間

を用いて，漏えい電界の住宅内特性を電磁界シミュレーシ

87 MHz

0.5 億

1

4 時間

815 MHz

11 億

16

10 時間

2.497 GHz

52 億

56

10 時間

ョンにより解析した．解析結果と測定結果を比較すること
で，住宅の基礎内部の鉄筋や外壁の材質が漏えい電界に与える影響が大きいことを把握した．測定結果と傾向が一致

大

大

波源

波源

小

小
（a） 1階（a） 1階（b） 2階（b） 2階図 13 住宅内における電界分布（2.497 GHz）

図 12 住宅内における電界分布（815 MHz）

45

46

（様式 20）成果報告書

TSUBAME 共同利用 平成 26 年度 産業利用 成果報告書

利用課題名 気象イベントを考慮した建築環境の解析評価
英文：Building Environment Analysis under considering Weather Events
利用課題責任者
PHAM VAN PHUC
所属： 清水建設（株）技術研究所
Affiliation： Institute of Technology, Shimizu Corporation
URL： http://www.shimz.co.jp
本利用課題では，TSUBAME スパコンに WRF 気象モデルを導入して，気象イベントを考慮できる建築環境の上空における大型台風や竜巻・突風等の局所的な気象現象の再現を試みる．本報では，過去に発生した竜巻イベントを対象として，初期データや解析条件及びネスティング方法等の検討により積乱雲や竜巻等の発達できるメソサイクロンの再現を確認し，そのイベントの再現性と WRF 気象モデルの有効性を明らかにした．
A meso scale meteorogical model has been used to develop the analysis system of build
environment considering the effects of extreme weather events such as large typhoons, tornadoes.
Several initial meteorogical conditions and nesting methods were investigated and it success to
reproduce a mesocyclone that needs for the developing of the cumulonimbus and tornado.

Keywords: Typhoon, Tornado, WRF, Building Environment, Weather Event
ル

背景と目的

WRF(Weather Researchand Forecasting

台風や竜巻・突風などの気象イベントによる人的，経Model)1)を採用した．WRF は圧縮性の非静力学モデ

済的ロスは極めて大きい．特に，竜巻・突風では建物やルで，物理・境界層過程と地表面過程を含んでいる．

鉄道施設に甚大な被害が発生しており，内閣府の防災

また，対象とした気象イベントは 2012 年 5 月 6 日 12基本計画にも風水害に係る事象として竜巻などの突風

時 35 分頃に茨城県常総市で発生した竜巻（フジタスケが明記され，竜巻・突風の発生予測は先端科学技術のール F3）によって，特に茨城県常総市とつくば市で建開発として重点的に取り組まれている．また，米国で築物等の多くの被害が生じた．図-1 はその竜巻による2011 年に発生した竜巻は，多数の死傷者の発生に加

建物の転倒被害である．建物の安全設計においてはこ

えて原子力発電施設の稼働に影響を与えるなど，気象

れらの転倒被害の詳細な解明は急務となる．

イベントは国際社会にも大きなインパクトを与えた．

解析としては，異なる数値予報データとネスティング

そこで，本利用課題はまず，TSUBAME スパコンの方法を用いて，そのイベントの再現性を検討していた．

計算資源を活用して，気象モデルを導入することにより，こうした局所的なイベントの再現にもかかわらず，用いた

これらの気象イベント，特に竜巻イベントの再現を試み

計算領域の影響が大きいことを確認し，広域な気象モ

た．これより，既存に開発された広域建築環境解析とのデル解析の必要性を確認できた．

連携を行い，記録結果や観測データとの比較により，その気象イベントの特徴と，地上付近建築環境等へ与える影響を調べており，気象イベントを考慮できる合理的な建築環境の構築を目指すものである．

概要
本研究では，気象モデルとして，米国大気研究センター（NCAR）等が中心となって開発したメモ気象モデ

図-1 竜巻による建物の転倒被害

47

（様式 20）成果報告書
計算概要

第1領域

計算期間は対象とした竜巻発生日（2012 年 5 月 6 日

第4領域

第2領域

0:00:00~18:00:00）とした．
初期値・境界値は，異なる気象データと計算手法を

対象イベントの発生位置

第3領域検討したが，本報では NCEP (National Centers for

第4領域

Environmental Prediction) の全球客観解析データ
FNL と気象庁の全休数値予報モデル GSM (Global
Spectral Model)データの結果を紹介する．

図-2 WRF 計算領域

計算領域は，親領域として水平解像度を 16km，格2012‐05‐06:03:00:00

2012‐05‐06:06:00:00

子数(nx=300,ny=300)，対象とした気象イベントの発生
位 置 で あ る 中 心 緯 度・ 経度 を 北 緯 36.11 度 ， 東 経
139.95 度とした．また，局所的な気象イベントを再現す

対象イベントの発生位置るためには，計算規模・計算間隔と計算領域の広さが計算結果に及ぼす影響を考慮して，4 段階 2-way ネス
ティングで計算領域を設定した．

2012‐05‐06:09:00:00

2012‐05‐06:12:00:00

図-2 には計算領域，ネスティング領域と対象イベントの発生地点の概要を示す．なお，鉛直方向の気圧準拠
座標については地表から 30 層を設定した．

結果および考察
図-3，4 には， GSM データと FNL データを用いた第4 領域の WRF の計算結果（5km 高度における上空の
風速・気温・相対湿度の分布と海面気圧分布）を示す．
それぞれは時刻 3 時，６時，９時と 12 時の結果となって

図-3 GSM データを用いた WRF 計算結果

いる．なお，対象とした竜巻イベントの発生時刻は１２時

2012‐05‐06:03:00:00

2012‐05‐06:06:00:00

前後となる．それぞれのデータの結果については，上
空の気温や相対湿度は概ね同様な傾向である．また，当日の気象観測データより上空の気温は-18 度や，地

対象イベントの発生位置

上の相対湿度は 61%に記録されており，計算の妥当性を概ねを確認できた．
上空風速分布については GSM データを用いた場合，

2012‐05‐06:09:00:00

2012‐05‐06:12:00:00

時刻９時に対象としたイベントの発生位の北側に小さな
メソサイクロンが発生していることを確認できるが，その後は北の方に移動してしまうことが分かる．
一方，FNL データを用いた場合，時刻９時には同様
なメソサイクロンが東側に発生している．時間を経過す
ると対象としたイベントの発生位置に移しており，時刻１
２時前後にその発生位置に到達している様子を確認できる．こうしたメソサイクロンは，竜巻をもたらした積乱雲
図-4 FNL データを用いた WRF 計算結果にともなうと考えられる．

48

（様式 20）成果報告書
まとめ，今後の課題
本利用課題では，WRF 気象モデルを導入して，
TSUBAME スパコンで複数計算条件や計算領域ネス
ティング方法と気象データを用いて，過去に発生した実在大型竜巻の再現計算を試みた．全球客観解析デー
タ FNL を用いた場合や 4 段階 2-way ネスティングでの計算により竜巻をもたらしなり得るようなメソサイクロンの再現を確認し，竜巻再現の可能性を明らかにした．
今後は，高解像度計算領域や高度の乱流モデルを用
いて，WRF の大規模計算を行うことにより，竜巻再現と地上付近の流れ場等の特性を明らかにする．

参考文献
1) WRF(2014): http://wrf-model.org/index.php

49

50

（様式 20）成果報告書

TSUBAME 共同利用 平成 26 年度 産業利用 成果報告書

利用課題名 鋼中析出物の水素捕捉能の高精度計算
英文：Calculation of hydrogen trapping at precipitate in steel
澤田 英明
Hideaki Sawada
新日鐵住金株式会社 先端技術研究所
Advanced Technology Research Laboratories, Nippon Steel & Sumitomo Metal Corporation
http://www.nssmc.com
鉄鋼業では、いかに鋼を高強度化するかということに注力してきたが、高強度化が進むにつれて、水素脆化という問題にも対処する必要性が高まってきた。その対処のために利用されているのが析出物であり、析出物は水素を一時的に捕捉し、後に環境の水素が少なくなった時に捕捉されていた水素を放出することでその役割を果たしている。今回、大規模第一原理電子状態計算パッケージ OpenMX での計算に必要な Fe、Nb、C の基底関数の最適化を行い、
bcc-Fe、NbC の格子定数、体積弾性率などの基本的物性値を再現することのできる基底関数の作製に成功した。
加えて、大規模計算を実施するために必要なオーダーN 法のパラメータについても網羅的に検討を実施し、計算精度と計算速度の両者を満足するパラメータの選択に成功した。
So far steel industry has mainly focused on strengthening of steel, however, hydrogen embrittlement
has become important subject as strength of steel is enhanced. Precipitates are used to inhibit thehydrogen embrittlement. Specifically hydrogen atoms are trapped at precipitates temporally, and thetrapped hydrogen atoms are released when hydrogen concentration in steel decreases. The basis
functions of Fe, Nb and C were successfully optimized so as to reproduce the basic properties, namely,
lattice constants and bulk moduli of bcc iron and NbC. Furthermore, the parameters for order-N
method were also optimized in order to carry out large scale simulation.

Keywords: first-principles calculation, order-N method, hydrogen enbrittlement, precipitate, steel
背景と目的

どこに水素が捕捉されているかは明確には分かってい鉄鋼材料に要求される特性として特に重要なのは、ないために、水素捕捉能を高める方策は不明である他、構造物を維持するための強度であり、それは自動車用

析出物種による水素捕捉能の違いは分かっていない。

鋼板のように衝突時にキャビンを守る強度から、建造

昨年度、NaCl 型の析出物である NbC と鉄母相との部物のように長期にわたって構造を保つ強度、ボイラー

分整合界面での水素捕捉位置を調べた。しかし、他析

管などのように高温での安定した強度まで多岐にわた出物や他の欠陥の水素捕捉能との比較をするためにる。そして、鉄鋼業の歴史としては、いかに鋼を高強度は計算精度を上げる必要があることが分かり、それに化するかということに注力してきた。しかし、高強度化取り組んだ。計算手法については、第一原理計算の大が進むにつれて、水素脆化という問題にも対処する必規模化可能なオーダーN 法のプログラム OpenMX（東要性が高まってきた。それは、強度が高まるほど水素

京大学尾崎教授開発）を用い、計算機資源として脆化に対する感受性が高まるためであるが、その対処

TSUBAME を利用させて頂き、計算精度向上のためのために利用されているのが析出物である。析出物との基底関数の最適化に取り組んだ。

して頻繁に使われているのは、VC などの NaCl 型の析出物であり、水素が大量に鋼中に侵入した際に、析出

概要物で一時的に捕捉し、後に環境の水素が少なくなった

最適化を行った基底関数を用い、Fe/NbC 整合界面時に捕捉されていた水素を放出している。それによって、の計算を対角化法とオーダーN 法で行い、その結果を脆性破壊につながる粒界の水素捕捉量を少なくするこ比較した。図 1 はその結果を示すが、オーダーN 法とが可能となる。しかし、NaCl 型の析出物においても、の結果は対角化法の結果を良く再現していることが51

（様式 20）成果報告書分かる。

乗に比例しており、計算規模を大きくすることは容易ではない。例えば、20 年後に計算可能な計算規模は、計

total energy (hatree)-566.64-566.65

diagonalization

算機の能力がムーアの法則に従って向上したとしても、-566.66

order-N_1

従来の第一原理計算法では 52 倍程度であり、1 万原

order-N_2

子程度といったところである。その問題を解決する 1 つ-566.67-566.68の方法としてオーダーN 法と呼ばれる手法がある。オ-566.69

ーダーN 法というのは計算量が計算規模に比例する方-566.70法という意味であり、この手法を用いれば、20 年後に-566.71は 145000 倍の規模、つまり、原子数にすれば、5000-566.72

万原子に近い原子数の計算が可能ということになる。

10.0 10.5 11.0 11.5 12.0 12.5

逆にいえば、オーダーN 法の開発を行わない限り、近c axis

い将来に計算規模を大きく増やすことは不可能である。
（図 1）Fe/NbC 整合界面の全エネルギーの c 軸長依存

オーダーN 法は、ワニア関数や密度行列の局在性を利性。

用することによって、計算規模を低減させる手法である。
現在、オーダーN 法の 1 つのプログラムとして、東京大学尾崎教授が開発している OpenMX を使用している結果および考察
材料の物性、特性は全て電子状態に起因している。

が、その中では局在基底として原子様数値基底（一電しかし、全ての物理量を電子状態から求めることは現子コーン・シャム軌道を擬原子軌道、つまり、球面調和

時点では不可能なため、原子間に働く力を経験的、若関数と動径分布関数の積、の線型結合で表現）を用い

しくは、第一原理計算によって決定されたパラメータをている。実際に計算を行うためには、この局在基底を用いる古典分子動力学法を用いたり、原子を扱わず、適用対象に対して最適なものを選ぶ必要があり、その空間的なメッシュ点での相を熱力学関数によって決定最適化を行った。

しながら組織変化を追うフェーズフィールド法なども用

0.08

いて、多種多様な物性値を計算機シミュレーションによvasp
Fe5.5H-s2p2d1
Fe5.5H-s2p2d1(order-N_1)
Fe5.5H-s2p2d1(order-N_2)

0.06

total energy (eV)って導出することが試みられている。ただ、古典分子動力学法の原子間に働く力として、高い信頼性を有しているのは単元素のポテンシャルが中心で、複数種の元素が存在する系でのポテンシャルについては、信頼性が低下するというのが現状である。そのため、複数種

0.04
0.02
0.00

の元素を含む系や欠陥を含む系の物性値を高い精度で求めようとすれば、電子状態から決定していくことが

-0.02
2.78

肝要となる。電子状態から決定するということは第一原

2.82

2.86

2.90

2.94

lattice parameter (A)理計算と呼ばれる手法を用いることに直結するが、現在の計算機を用いて一般的に行うことができる計算規（図 2）基底関数 Fe5.5H-s2p2d1 を用いた時の対角化

模は数 100 原子以下であり、複雑な構造（金属系の材法とオーダーN 法での bcc-Fe の全エネルギーの格子料を考えた場合には、粒界、異相界面などの欠陥がそ

定数依存性。

の 1 つである）の物性の導出を考えた場合には少なくと

図 2 は bcc の Fe における全エネルギーの格子定数も 1000 原子以上の大規模な電子状態計算を行う必要依存性である。Fe5.5H-s2p2d1 は、カットオフ半径ががある。しかしながら、従来から用いられている第一原

5.5a.u.で、基底関数として s 軌道を 1 つ、p 軌道を 2 つ、

理計算手法に必要な計算量は計算規模（原子数）の 3d 軌道を 1 つ使用したもので、充分な精度で VASP の52

（様式 20）成果報告書
結果を再現していることが分かる。対角化法では計算果を示した。表2から分かるように、オーダーN法

する系全体の電子状態を一度に計算することになるが、の計算によって、対角化法によって計算された格オーダーN 法ではある原子の電子状態をその周囲の子定数と体積弾性率を良く再現していることが物理的に切り出されたクラスタに対する固有値問題を分かる。

解くことによって求め、その足し合わせによって全系の（表 2）対角化法とオーダーN 法で計算された NbC

電子状態を得ることになる。そのため、これらの基底関の平衡格子定数と体積弾性率。

数を用いてオーダーN 法の計算を行うためには、どの計算手法ようなクラスタを考えるかを決める必要がある。図 3 に平 衡格子定 数

体積弾性率（MPa）（Å）

は、クラスタに入る原子数（表 1 の計算手法の（）中の対角化法

4.5014(-0.3%)

298.0(3.8%)

数字）を変化させたときの計算結果を示したが、オーダ

ｵｰﾀﾞｰN法1

4.4980(-0.4%)

301.2(4.9%)

ーN 法での計算が対角化法の結果を再現していること

（340）

が分かる。

ｵｰﾀﾞｰN法2

4.4980(-0.4%)

301.2(4.9%)（表1）対角化法とオーダーN法で計算された（256）

bcc-Feの平衡格子定数と体積弾性率。

上で最適化された基底関数を用い、Fe/NbC 整合界

計算手法
対角化法
ｵｰﾀﾞｰN法1（330）
ｵｰﾀﾞｰN法2（258）

平 衡格子定 数（Å）

面の計算を実施した。具体的には、基底関数

体積弾性率（MPa）

2.8432(0.5%)

201.3(6.8%)

2.8374(-0.3%)

216.3(14.8%)

2.8432(0.5%)

201.3(6.8%)

Fe5.5s2p2d1 、 Nb7.0-s2p2d1 、 Ti7.0-s2p2d1 、
C6.0-s2p2d1 を用い、Fe/NbC 整合界面の計算を、対角化法とオーダーN 法で行い、その結果を比較した。
図 1 には Fe/NbC の整合界面に対する結果を示すが、オーダーN 法の結果は対角化法の結果を良く再現していることが分かる。

0.08
0.06

total energy (eV)

まとめ、今後の課題

vasp
Nb7.0-s2p2d1
Nb7.0-s2p2d1(order-N_1)
Nb7.0-s2p2d1(order-N_2)

Fe 、Nb、C に対して基底 関数の最適化を行い 、
bcc-Fe、NbC の格子定数、体積弾性率などの基本的

0.04

物性値を再現することのできる基底関数の作製に成功
0.02

している。加えて、大規模計算を実施するために必要なオーダーN 法のパラメータについても網羅的に検討0.00を実施し、計算精度と計算速度の両者を満足するパラメータの選択に成功した。

-0.02
4.40

4.45

4.50

4.55

4.60

今年度作製した基底関数を用い、鋼中での析出物

lattice parameter (A)の水素捕捉能のより高精度な計算を実施していく予（図 3）基底関数 Nb7.0-s2p2d1、C6.0-s2p2d1 を用い定である。

た時の対角化法とオーダーN 法での NbC の全エネル
ギーの格子定数依存性。

図3はNbCにおける全エネルギーの格子定数依存性であるが、Nb7.0-s2p2d1、C6.0-s2p2d1が充分な精度でVASPの結果を再現していることが分かる。更に、クラスタに入る原子数（表2の計算手法の（）中の数字）を変化させたときの計算結

53

54

TSUBAME 共同利用 平成 26 年度 産業利用 成果報告書

利用課題名 拡張アンサンブルシミュレーションによるタンパク質とリガンドの結合構造予測の展開
英文：Applications of Prediction Method of Protein-ligand Binding Modes by Generalized-ensemble
Simulations
利用課題責任者 田中 稔祐
Toshimasa Tanaka
所属 武田薬品工業株式会社 医薬研究本部 化学研究所
Medicinal Chemistry Research Laboratories, Pharmaceutical Research Division,
Takeda Pharmaceutical Company Limited
http://www.takeda.co.jp
邦文抄録 二次元レプリカ交換法を膜タンパク質（GPCR の一つであるβ2 adrenergic receptor）に対して適用したところ、予測構造は実験構造とよく類似しており、結合モードを再現できることが分かった。このことは本
課題を通して開発した手法が水溶性タンパク質だけでなく膜タンパク質に対しても適用可能であることを示唆している。一方で、新たな課題として膜タンパク質の系ではリガンドが脂質二重膜の膜面にトラップされてしまう場
合があることが分かった。これはサンプリング上必ずしも好ましくないため、今後さらに改良と検証を進めて行きたい。さらに JNK3 に対して結合自由エネルギー計算を実施したところ、複数の異なる結合モードから正しい結合モードを計算により判別できた。
英文抄録 We applied two-dimensional replica-exchange method for the prediction of a membraneprotein (β2 adrenergic receptor) system, and found that our method successfully predicted thestructure close to the experimental one. This indicates that the method developed through thisresearch might be applied not only to soluble proteins but also membrane proteins. However, we
encountered a new difficulty that a ligand tends to be trapped around the membrane surface. We
also performed binding free energy calculations for the JNK3 system, and found that the correctbiding mode was distinguished computationally.

Keywords: replica-exchange umbrella sampling, molecular dynamics, protein-ligand binding,
computer-aided drug design, free energy calculation.

精度化と検証が必要とされている。本利用課題では相背景と目的タンパク質に結合する分子（リガンド）の結合構造や互作用の物理化学的背景の理解に基づいた論理的な結合の強さを予測することは、様々な生体反応を理解

医薬品の分子設計を可能とするために、昨年度までにし制御するために極めて重要であり、ドラッグデザイン

開発した二次元レプリカ交換法を膜タンパク質に対してにとっても有用性が高い。

適用し、結晶構造解析が困難な膜タンパク質に対して

そこで我々はこれまでにレプリカ交換法に基づく高精の有効性を一例で示した。図 1 は膜タンパク質の一つ

度なリガンド結合モードの予測法を開発してきた[1-4]。

であるβ2 adrenergic receptor に対して開発した手法

本課題では、特に膜タンパク質の一つであるβ2を適用した結果であり、図 1 左のようにリガンドが溶媒

adrenergic receptor に対して、開発した手法を適用し、側にいる構造から始めて、図 1 右のように実験と近い

膜タンパク質系に対する適用可能性と課題を調査する

構造を安定構造として検出できた。

と共に、結合自由エネルギー計算との連携を強化することを目的として研究を実施した。
概要
従来の手法では実験を置き換えるレベルに至っておらず、それを目指したシミュレーション手法のさらなる高55

図 1 β2 adrenergic receptor に対する結合モード予ての平均力ポテンシャルを計算し、安定な結合距離を測計算の初期構造（左）と予測結合モード（右）予測することができる。図 1 左に、予測されたリガンドの結合モードと実験の結合モードとの比較を示す。完全ではないものの、計算により予測された構造は概ね位

また、自由エネルギー計算との連携を強化に取り組んだ。図 2 はドッキングから予測されるリガンドの 4 つの置向き共に実験構造と一致していることが分かる。

異なる結合モードであり、それぞれに対して結合自由エ

次に、JNK3 の系（PDB ID: 1PMV）でドッキングに

ネルギー計算を実施したところ、実験と近い正しい結合より予測した 4 つの異なる結合モード（図 2 左）から正しモードが最安定であった。

いモードを結合自由エネルギー計算により識別できるか調べた。図 4 は複合対中での<dU/dλ>の静電項とvdW 項をそれぞれ示す。共に収束にとって必要な滑らかな曲線となっている。

図 2 JNK3 に結合するリガンドの結合モードの結合自由エネルギー計算による判別
図4
さらに主成分解析などの手法を用いた機能性ダイナ計算から予測された結合自由エネルギーは、図 2 のミクスの解析をめざし、GPU を効率的に使用した長時

左上、右上、左下、右下の結合モ ード で、そ れぞれ

間分子動力学計算を実施した。現在、構造の安定性、-11.9, -8.9, -8.9, -8.0 kcal/mol であり、左上の結合モ

長時間運動モード、構造変異のしやすさ、などについてードが安定性が高いと予測された。図 2 右に示すように実験データと比較しながら解析を鋭意進めている。

実験の結合モードと予測モードを比較したところよく一結果および考察致していた。また、予測された値-11.9 kcal/mol は実験の結合自由エネルギー-10.5 kcal/mol とも数値的に近最初に、開発した手法を膜タンパク質の一つである
β2 adrenergic receptor の系（PDB ID: 2RH1）に対しかった。

て適用した（手法の詳細については[1-4]を参照）。

まとめ、今後の課題
二次元レプリカ交換法を膜タンパク質（GPCR の一つであるβ2 adrenergic receptor）に対して適用し、実験と近い結合モードを予測できるという有望な結果を得た。このことは開発した手法が水溶性タンパク質だけでなく膜タンパク質に対しても適用可能であることを示唆している。特に膜タンパク質は結晶化が一般的に困難であり、計算によりリガンドの結合モードを高精度に予図 3 β2 adrenergic receptor の系のスナップショット

測できれば、大きなインパクトになると考えられるため、今後さらに改良と検証を進めて行きたい。また、JNK3

図 3 はシミュレーションの代表的なスナップショットの系に対して結合自由エネルギー計算を実施したとこを示しており、リガンドが溶媒側にある場合や、膜面付ろ、複数の異なる結合モードから正しい結合モードを計

近にある場合、タンパク質のへリックス間にある場合な算により判別できることが分かった。さらに、GPU を使ど様々な構造を探索していることが分かる。全てのレプ用した長時間ダイナミクスを調べる計算も開始し、いく

リカの情報を集めてタンパク質とリガンドの距離についつかのターゲットタンパク質に対して 1 マイクロ秒程度

56

の長時間計算を、比較のため異なる複数の初期構造や初期速度を用いて実施した。構造の安定性、長時間運動モード、構造変異のしやすさ、などのタンパク質による違いの知られている実験データと比較しながら解析を鋭意進めている。
参考文献

1. H. Kokubo, T. Tanaka, Y. Okamoto, Ab initio
prediction of protein-ligand binding structuresby replica-exchange umbrella sampling simulations, J. Comput. Chem. 32, 2810–2821 (2011).
2. H. Kokubo, T. Tanaka, Y. Okamoto,Two-dimensional replica-exchange method forpredicting protein–ligand binding structures, J.

Comput. Chem. 34, 2601-2614 (2013).
3. H. Kokubo, T. Tanaka, Y. Okamoto, Predictionof Protein–Ligand Binding Structures byReplica-Exchange Umbrella Sampling
Simulations: Application to Kinase Systems, J.

Chem. Theory Comput. 9, 4660-4671 (2013).
4. Y. Okamoto, H. Kokubo, T. Tanaka, Predictionof Ligand Binding Affinity by the Combination of
Replica-Exchange Method andDouble-Decoupling Method, J. Chem. Theory

Comput. 10, 3563-3569 (2014)

57

58

共同利用 成果報告書 平成 ２６ 年度 課題種別

利用課題名 電子デバイス材料の計算機設計
英文：Computational design of materials for electronic devices
渥美 照夫、 岩崎 誉志紀
Teruo Atsumi, Yoshiki Iwazaki,
太陽誘電株式会社
TAIYO YUDEN CO., LTD.
http://www.yuden.co.jp/jp/
本報告では、電子デバイス用途材料として多用される遷移金属酸化物と電極界面の酸素欠陥安定性に関する第一原理計算として、ペロブスカイト型構造酸化物のチタン酸ストロンチウムと白金(Pt)電極を例に報告する。さらに、電子デバイス用途・ 有機電解液の大規模計算に向けて、 GPGPU(General-purpose computing on graphics
processing units)を活用した古典 MD の検討結果についても報告する。
We report first-principles theoretical study of oxygen vacancy at the interface between perovskite-typeoxide (strontium titanate) and metal-electrode (platinum). Especially, we will force on the stability ofthe oxygen vacancy at the metal/insulator interface. Furthermore, we also present our study onclassical MD simulations with GPGPU aiming the clarification of properties observed in electrolytic
solutions used for energy storage devices.

Keywords: First-principles calculation, perovskite, interface, classical molecular dynamics, electrolytic
solution

背景と目的(1) 第一原理計算による誘電体/電極界面検討 1

電子デバイスに用いられる無機材料、有機材料には様々な機能・特性の発現が求められている。半導体や概要(1)酸化物等の無機材料においては、ナノレベルオーダー

誘電体/電極界面における酸素欠陥(VO)の挙動は積の電極界面構造がデバイスの信頼性や特性に大きく

層セラミックキャパシタ(MLCC)の特性と高い関連性を影響する事が知られており、これらの理解と制御が材持っており、例えば界面付近の VO 量とデバイス寿命と料開発における鍵となる。また、リチウムイオン電池やの関連性や、チタン酸バリウムと Ni 電極界面に、O がリチウムイオンキャパシタ等のエネルギーデバイス開欠乏した Ba-Ti-Ni 合金層が存在する事が報告されて発では、使用される電解液の物性予測に対し古典分子いる。しかし界面近傍の VO を、直接実験で測定する事

動力学(Molecular dynamics: MD)シミュレーションには困難であり、その挙動については不明点も多い。そよる大規模計算が有効である事が示されているが、ここで本研究では、代表的なペロブスカイト型酸化物であれ ら の 理 論 計 算 は GPGPU (General-purpose

るチタン酸ストロンチウム(SrTiO3:ST)と Pt 電極をモデ

computing on graphics processing units)の活用にルとして取り上げ、第一原理計算により酸化物/金属界より高速化や大規模化が可能である事が期待される。

面における VO の安定性を調査した。

本利用課題では、実験のみのアプローチでは理解が困難なこれらの課題に対して原子レベルでの計算機シ

第 一 原 理 計 算 は 、 密 度 汎 関 数 法 の GGA+U 法

ミュレーションを行い、そのベースとなる基礎技術開発(UTi=3eV)を用い、計算パッケージは VASP を使用した。

を目指した。本報告書では、「(1)第一原理計算による

構 造 は 、 ス ラ ブ モ デ ル を 用 い 、 真 空 中 ST 薄 膜

誘電体/電極界面検討」と「(2)古典 MD シミュレーション(ST/Vac.)と ST-Pt 電極(ST/Pt)について検討した(Fig.

による電解液検討」について報告する。

1) 。ST/Vac.は 、ST(100) 面 の両 端を TiO2-term か59

SrO-term になるように切り出し、真空層 20Å を加えた。

移動量と、EVO には相関があり、VO の安定化が Pt 電

ST/Pt は、ST(100)面に Pt(100)面を接合させたモデル

極への電荷移動によるものであると推測された。

をいくつか検討し、最安定な構造(TiO2-term の ST に

6.0

Pt を接合)を用いた。いずれも構造最適化計算により5.5

安定化させ、これに VO を l = 1-10 (l は界面からの層
5.0

EVO [eV]数)の位置に導入し、さらに構造最適化計算(格子定数は固定)を行った。
l = 1 2 3 4 5 6 7 8 9 10 9 8 7 6 5 4 3 2 1

4.5
4.0

Bulk
ST(SrO-term)/Vac.

3.5

ST(SrO-term)/Vac.

ST(TiO2-term)/Vac.
ST(TiO2-term)/Pt

3.0

ST(TiO2-term)/Vac.

2.5
ST(TiO2-term)/Pt

1

Figure 1. Structure of slab models. l mean number of layer from interface.

2

3

4

5
6
7
VO position l

8

9

10

Figure 2. Formation energy of oxygen vacancy EVO.

結果および考察(1)1Charge transfer ST to Pt [e-]電極界面からの距離に対する VO の生成エネルギー
EVO を下記の式(1)に従って求めた。
EVO =E(O defect ST slab) + 1/2E(O2 gas)
－E(Perfect ST slab)(1)

ここで、E(O defect ST slab)は酸素欠陥を含むスーパーセルの全エネルギー、E(Perfect ST slab)は完全 ST 結晶の全エ
ネルギー、E(O2 gas)は酸素分子の全エネルギーである。

0.8

0.6

0.4

O defect slab
0.2

Perfect slab

得られた結果を Fig. 2 に示す。Fig. 2 には ST: 3×3×3
サイズのセルを用いて計算されたバルク ST のEVO が0

1

黒実線で同時に示されている。ST/Vac.では、VO 位置

2

3

4

5
6
7
VO position l

8

9

10

Figure 3. Number of charge transfer from ST to Pt.

や末端構造によらず、bulk とほぼ同じ VO 生成エネル
ギー値となった一方、ST/Pt では Pt 電極に近いほどまとめ、今後の課題(1)EVO が低くなった。特に l = 1, 2 の VO は bulk におけ第一原理計算により ST/Pt 界面では、Pt 電極に近いる値と比べて、それぞれ 2.56eV, 1.32eV 安定化しておほど VO が安定化することが分かった。また、安定化要

り、Pt 電極近傍で VO が大きく安定化することが分かっ因の一つとして界面間電荷移動が関係していると推測た。Pt 電極による VO 安定化が生じる要因はいくつか想された。しかしながら、今回用いたモデルでは、セルサ定されるが、VO の生成により結晶内に放出された余剰

イズが小さいため、深さ方向への収束性を検討するに電子が Pt 電極へ移動して VO の生成が促進されているは不十分であり、面方向の影響を考慮することが出来可能性がある。そこで VO 導入前後の ST/Pt に対して電ないため、今後の課題としてより大きなセルによる検討

荷解析（Bader charge analysis）を適用し、ST から Ptが必要と考えている。

への電荷移動量 (CTSTtoPt)を求めた (Fig. 3)。VO 添
加前(破線)の CTSTtoPt は 0.04 であり、ST と Pt を接合(2) 古典 MD シミュレーションによる電解液検討しただけでは、殆ど電荷移動しない。一方、VO 添加後概要(2)の CTSTtoPt は、0.41～0.97 となっており、Pt 電極に近昨年度の検討では、純溶媒の古典 MD シミュレーシ

い程 ST から Pt へ電荷移動することが分かった。このョン(Gromacs)のベンチマークテストを行い、GPGPU

60

の有効性と 10 万原子程度のシミュレーションが実用的

実デバイスで用いられている電解液の組成は、より複に計算可能であることを確かめた。本年度は、リチウム

雑（添加物、混合溶媒）であり、これを扱うには、より多イオン電池等の開発で重要となるリチウムイオンの拡くの分子による計算が必要となると予想される。

散係数を求める事を目的として、電解質を添加したシミュレーションを行った。

(a) [BF4-+Li+]

計算条件は、溶媒 2,000 分子に電解質を 0.0, 0.5,
D (solvent) [10-5cm2/sec]

1.0, 2.0M 添加したものについて、平衡化の後に熱平衡状態のシミュレーションを行った。溶媒は、エチレンカ
ーボネート(EC), プロピレンカーボネート(PC), ガンマ
ブチロラクトン(GBL), ジメチルカーボネート(DMC), エ(b) [PF6-+Li+]

1.6

1.6

1.4

1.4

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

チルメチルカーボネート(EMC), ジエチルカーボネート

0(DEC)の 6 種、電解質は[BF4-+Li+], [PF6-+Li+]の 2 種

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

D (Li+) [10-5cm2/sec]

初期構造は、分子をランダム配置したものを用い、NPT ア ンサンブル ( 温度と 圧力はそ れぞれ 300K,
1atm)のシミュレ―ションを実施した。平衡化は時間刻
み 0.5fs を 6,000,000step(300ps)、熱平衡状態は時間

0.5

1

1.5

2

0

D (anion) [10-5cm2/sec]

テンシャルとして、溶媒は general AMBER 力場を、電解質は文献値[Journal of Molecular Liquids 148(2009) 99–108]を用いた。

0.5

1

1.5

2

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

GBL

0

0.5

1

1.5

2

0

0.5

1

1.5

2

0.5

1

1.5

2

0

0

0

0.5

1

1.5

2

[BF4-+Li+] concentration [M]

結果および考察(2)

PC

DEC

0

0

刻み 0.2fs を 20,000,000step(4ns)で行った。また、ポ

EC

EMC

0

0

とした(原子数は 20,000~39,881)。

DMC

0[PF6-+Li+] concentration [M]Figure 4. Diffusion coefficient(D) of (a) [BF4-+Li+] or

熱平衡状態のシミュレーションの平均二乗変位(b) [PF6-+Li+] added electrolytic solution.

(MSD)から求めた拡散係数 D を Fig. 4 に示す。これから、電解質濃度が高いほど、いずれも D が減少する事が示された。また、電解質の種類に対する D の変化と

1 第 32 回強誘電体応用会議

して、D(solvent)>D(anion)>D(Li+)といった、実験と対ト型酸化物/金属界面における酸素欠陥生成メカニ

応する傾向が得られた。今回の計算結果では、BF4-のズムに関する第一原理計算」

ほうが PF6+より D(anion)と D(Li+)の値が大きくなって
おり、実験とは逆の傾向を示しているが、この原因の一つとして、電解質が全て解離しているという条件で計算を行っている事などが考えられる。
まとめ、今後の課題(2)
電解質を含めた電解液の古典 MD シミュレーションを行い、リチウムイオン等の拡散係数を求めた。その結果、実験とおおよそ対応する傾向が得られた。今後の課題として、電解質による差を再現するために、解離度の考慮等を検討していく必要があると考えられる。また、

61

20-B-04「ペロブスカイ

62

TSUBAME 共同利用 平成 26 年度 産業利用 成果報告書

理論計算に基づく有機半導体材料の開発
Theoretical Design of Advanced Organic Semiconductor
栗田 靖之
Yasuyuki Kurita
住友化学株式会社・先端材料探索研究所
Advanced Materials Research Laboratory, Sumitomo Chemical Co., Ltd.
http://www.sumitomo-chem.co.jp/
第一原理バンド計算により、チオフェン系高分子の固体状態では、一部の例外を除き、スタッキング相互作用
によって孤立状態よりもヘテロ環間の捩れ角が小さいという結果を得た。この捩れ角減少に伴い、フルオレン系
高分子も含めたイオン化ポテンシャルの固体状態における実測値と高分子鎖１本に関する計算値との間の相関係数は 0.778 から 0.949 に向上した。固体状態モデルから切り出した高分子鎖３本について計算したイオン化ポ
テンシャルは、高分子鎖間の分子軌道相互作用により、高分子鎖１本を切り出した場合より小さくなるが、固体状態における実測値との相関係数は 0.953 であり、向上幅は僅かであった。
With first-principles band calculations, the torsion angle between heterocycles in the solid state
thiophene-type conducting polymer is predicted to be smaller than that in the isolated state except a
few cases. By this change in torsion angles, the correlation coefficient between the experimental
ionization potentials of solid state polymers and the calculated ones for single polymer chains is
improved from 0.778 to 0.949. The ionization potential calculated for the three polymer chainsextracted from the solid state model is lower than that calculated for the single chain because of themolecular orbital interaction between polymer chains. The correlation coefficient for the threepolymer chains model is 0.953, which is almost the same as that for the single polymer chain model.

Keywords: 有機エレクトロニクス材料／導電性高分子／イオン化ポテンシャル／密度汎関数法

体について、密度汎関数法を用いて計算した孤立状態

背景と目的
有機ＥＬ（電界発光）ディスプレーや有機トランジスタの１分子構造を固定し、複数分子間の相対配置のみをー、有機太陽電池といったデバイスを実現するための分子力場法で最適化して固体状態モデルを構築し、イ

材料として有機半導体の開発が進められている。 そのオン化ポテンシャルを計算した。その際、チオフェン系

中でも導電性高分子材料は、インクジェット法や印刷法

高分子について、アルキル側鎖の炭素原子を、それが等の安価なデバイス作製プロセスを適用可能なため、属するヘテロ環平面上に固定すると、ヘテロ環間の捩

注目を集めている。また、その柔軟性を活かし、フレキ

れ角が小さくなり、フルオレン系高分子も含めたイオン

シブルディスプレーやフレキシブル太陽電池など、新し化ポテンシャルの実測値と計算値の相関が高くなること

い形態の製品を作ることが可能である。

を見出した１）。一方、Ｘ線回折実験に基づいて、上記ア

上記有機半導体材料の重要な特性の一つに、電極

ルキル側鎖はヘテロ環平面上にないと推定した研究がからのキャリア（電子 or ホール）注入性がある。キャリアある２）。以上のことから、固体状態では、ヘテロ環間の捩注入性が良いほど、低電圧でデバイスを駆動することがれ角は孤立状態における捩れ角よりも小さいが、その原でき、消費電力を低く抑えることができる。有機半導体

因は、アルキル側鎖がヘテロ環平面上に存在することのイオン化ポテンシャルは、上記ホール注入性に関係ではなく、分子間のスタッキング相互作用にあると推定があると考えられており、有機半導体材料を設計する上した。

で重要な指標である。

そこで本年度は、高分子間の相対配置だけでなく、本プロジェクトでは、2013 年度、高分子系有機半導分子間スタッキング相互作用による高分子の構造変化

63も考慮した計算を、第一原理バンド計算法を用いて実Gaussian09 プログラム７）を用いた。イオン化ポテンシ施した。その結果、チオフェン系高分子の固体状態に

ャル実測値（固体状態、光電子分光法による測定）

おけるヘテロ環間の捩れ角は、孤立状態における捩れ

は文献を参照した８）。

角よりも小さいと推定された。また、フルオレン系高分子

チオフェン系高分子

フルオレン系高分子

C6H13も含めたイオン化ポテンシャルの実測値と計算値の相T1

関は、孤立状態の構造で計算するよりも高くなった。

H17C8
θ1

F1

θ1 S
S

6
C6H13

H17C8

C8H17
2.5

C8H17

C14H29T2

概要

θ1 S θ2

S

S

3

S

F2

θ1 S
H17C8

C14H29

θ2
S

2.3

C8H17

CO2C12H25１）高分子系有機半導体の固体状態モデルの構築

T3高分子の構造は、第一原理バンド計算プログラム
VASP ３） を利用し、周期的境界条件下、密度汎関数

T4法を用いて計算した。密度汎関数には、ファンデル

S

T5

孤立状態の構造は、大きなユニットセルを設定し、高分子間の相互作用が十分小さい状態で計算した。
計算対象とした高分子の側鎖は炭素数６以上で、高

T6

S

S
C12H25O2C
S
C16H33

ワールス力を考慮した PBE-D2４）を採用した。

θ1 S

θ2
3

θ1

F3
H17C8

C8H17

θ2 θ3
N
θ4

2.5

C16H33S
C16H33S
θ1 S
2S
C16H33S

C16H33
θ2
θ1 S
S

S

3

S
C16H33
C16H33
θ1 S

S
S

C16H33S

3.5

図１．計算した高分子系有機半導体のモデル

分子間の相互作用を小さくするには、かなり大きなユ
ニットセルを設定する必要があり、計算時間がかかる
ため、側鎖を全て n-Pr で置換した（孤立状態の場合、

結果および考察
各高分子の固体状態モデルの構造を図２～図１０にn-Pr よりアルキル側鎖を長くしてもヘテロ環間の捩示す（上下方向にスタックした高分子鎖は２本のみ示し、れ角とイオン化ポテンシャルの計算値の変化は小さそれぞれの高分子鎖を色分けした）。チオフェン系高分い）。

子では、高分子主鎖間でヘテロ環同士がπ-πスタック

固体状態の構造を計算する場合は、実際の側鎖をするが、フルオレン系高分子では、フルオレン環の上下用いた。高分子間のスタッキング相互作用に側鎖間に出たアルキル側鎖のため、π-πスタックが起こらないの相互作用も加わるため、n-Pr に置換すると高分子構造が計算された。

間の距離が変化し、その結果、イオン化ポテンシャ
ルの計算値も変化するからである。構造は複数構築した中から最も安定なものを選択した。
２）高分子系有機半導体のイオン化ポテンシャルの計算
１）で求めたユニットセルを拡張してスーパーセルを構築し、その中からチオフェン環１２個程度の長さの部分を切り出し、末端を水素原子で封止してモデ
ル高分子鎖とした（図１）。孤立状態の構造について
は、高分子鎖１本から電子を１個取り去った状態と元の状態とのエネルギー差を密度汎関数法（密度汎関
数 B3LYP５）、基底関数 6-31G*６））で計算し、イオン化ポテンシャルとした。固体状態については、高分
子鎖１本、或いは、スタッキングした高分子鎖３本をスーパーセルから切り出し、イオン化ポテンシャルを

図２．T1(図１)の固体状態モデルの構造計算した。イオン化ポテンシャルの計算には

64

図３．T2(図１)の固体状態モデルの構造

図５．T4(図１)の固体状態モデルの構造

図６．T5(図１)の固体状態モデルの構造

図４．T3(図１)の固体状態モデルの構造65

図９．F2(図１)の固体状態モデルの構造

図７．T6(図１)の固体状態モデルの構造
図６．T5(図１)の固体状態モデルの構造

図１０．F3(図１)の固体状態モデルの構造

図１１に孤立状態と固体状態の共役環間の捩れ角
計算値の差を示す。チオフェン系高分子では T4 以外
は固体状態の捩れ角の方が孤立状態の捩れ角よりも小さい。これは高分子鎖間のスタッキング相互作用が原因と考えられる。
図１２に孤立状態モデルの高分子鎖１本、図１３に固体状態モデルから切り出した高分子鎖１本、
図８．F1(図１)の固体状態モデルの構造
10Δθ (deg.)

5
0
-5
-10
-15
-20
-25
-30

θ1

θ1
θ2

θ1
θ2

θ1
θ2

θ1
θ2

θ1

θ1

θ1
θ2

θ1 θ3
θ2 θ4

T1

T2

T3

T4

T5

T6

F1

F2

F3

図１１．固体状態と孤立状態の共役環間の捩れ角の差 (θ1～4 の定義は図１参照)

66

図１４に固体状態モデルから切り出した高分子鎖３6.0

Ionization Potential (eV, Exp.)

本のイオン化ポテンシャル計算値に対し、固体状態の高分子のイオン化ポテンシャル実測値をプロットした。相関係数は、それぞれ 0.778、0.949、0.953である。孤立状態から固体状態になる際の共役環間の捩れ角の変化により、相関が大幅に向上した。高分子鎖間の分子軌道相互作用によりイオン化ポテン
シャルの計算値は小さくなるが、相関係数の向上幅は小さい。

r=0.953
F3

5.5

T6
T2 T4

5.0

F2

F1

T3
T5

T1
4.5

4.0
4.0

Ionization Potential (eV, Exp.)6.0

r=0.778
F3

5.5

T5
T4

5.0

5.0

5.5

6.0

F1

F2

図１４．固体状態における高分子のイオン化ポテ

T6 T3

ンシャル実測値と固体状態モデルから切り出した

T2

高分子鎖３本のイオン化ポテンシャル計算値との

T1

関係

4.5

4.0
4.0

4.5

5.0

5.5

Ionization Potential (eV, Calc.)

まとめ、今後の課題

6.0

第一原理バンド計算により、チオフェン系高分子の固体状態では、一部の例外を除き、スタッキング相互

図１２．固体状態における高分子のイオン化作用によって孤立状態よりもヘテロ環間の捩れ角がポテンシャル実測値と孤立状態モデルの高分小さいという結果を得た。この捩れ角減少に伴い、フ子鎖１本のイオン化ポテンシャル計算値とのルオレン系高分子も含めたイオン化ポテンシャルの関係
6.0

Ionization Potential (eV, Exp.)

4.5

Ionization Potential (eV, Calc.)

固体状態における実測値と高分子鎖１本に関する計
算値との間の相関係数は 0.778 から 0.949 に向上し
r=0.949
F3

5.5

F2

た。固体状態モデルから切り出した高分子鎖３本にF1ついて計算したイオン化ポテンシャルは、高分子鎖間の分子軌道相互作用により、高分子鎖１本を切り

T3 T6
T5 T2
T4

5.0

出した場合より小さくなるが、固体状態における実測値との相関係数は 0.953 であり、向上幅は僅かであっ

T1

た。

4.5

今後の課題は、イオン化ポテンシャルの実測値と計算値の差がより小さい計算モデルと計算方法の探

4.0
4.0

4.5

5.0

5.5

Ionization Potential (eV, Calc.)

6.0

索である。図１３においては、実測値と計算値の差は比較的小さいが、より現実に近いモデルに関する図１

図１３．固体状態における高分子のイオン化

４では差が大きい。また、別系統の高分子系有機半ポテンシャル実測値と固体状態モデルから切導体についても検討し、チオフェン系、及びフルオレ

り出した高分子鎖１本のイオン化ポテンシャル

ン系と比較する予定である。

計算値との関係

67

参照文献

L. Martin, K. Morokuma, V. G. Zakrzewski, G. A.

１）栗田靖之, 平成26年度東京工業大学TSUBAME産

Voth, P. Salvador, J. J. Dannenberg, S. Dapprich,

業利用シンポジウム, 203-208 (2014).

A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V.

(http://www.gsic.titech.ac.jp/sites/default/files/

Ortiz, J. Cioslowski, and D. J. Fox, Gaussian, Inc.,

H25_sangyo_sumitomo-chem.pdf)

Wallingford CT, 2009.

２）Y. Olivier, D. Niedzialek, V. Lemaur, W. Pisula, K.

８）a)T. Yamanari, T. Taima, J. Sakai, and K. Saito,

Müllen, U. Koldemir, J. R. Reynolds, R. Lazzaroni,

Solar Energy Materials & Solar Cells, 93,

J. Cornil, and D. Beljonne, Adv. Mater., 26, 2119-

759-761 (2009).

2136 (2014)．

b) I. McCulloch, M. Heeney, M. L. Chabinyc,

３）a) G. Kresse and J. Hafner, Phys. Rev. B, 47, 558

D. DeLongchamp, R. J. Kline, M. Cölle, W.
Duffy, D. Fischer, D. Gundlach, B. Hamadani,-561 (1993).
b) G.Kresse and J. Furthmüller, Phys. Rev. B, 54,

R. Hamilton, L. Richter, A. Salleo, M. Shkunov,
D. Sparrowe, S. Tierney, and W. Zhang, Adv.

11169-11186 (1996).
４）a) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys.

Mater. , 21, 1091-1109 (2009).

Rev. Lett., 77, 3865-3868 (1996).

c) J. Hwang, E.-G. Kim, J. Liu, J.-L. Brédas, A.

b) S. Grimme., J. Comp. Chem.. 27, 1787-1799

Duggal, and A. Kahn, J. Phys. Chem. C, 111,

(2006).

1378-1384 (2007).

５）a)A. D. Becke, J. Chem. Phys., 98, 5648-5652(1993).
b) C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B,
37, 785-789 (1988) .
６）P. C. Hariharan and J. A. Pople, Chem. Phys. Lett.,
16, 217-219 (1972).
７）Gaussian 09, Revision D.01, M. J. Frisch, G. W.
Trucks, H. B. Schlegel, G. E. Scuseria, M. A.
Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B.
Mennucci, G. A. Petersson, H. Nakatsuji, M.
Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J.
Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M.
Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M.
Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai,
T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F.
Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K.
N. Kudin, V. N. Staroverov, R. Kobayashi, J.
Normand, K. Raghavachari, A. Rendell, J. C.
Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N.
Rega, J. M. Millam, M. Klene, J. E. Knox, J. B.
Cross, V. Bakken, C. Adamo, J. Jaramillo, R.
Gomperts, R. E. Stratmann, O. Yazyev, A. J.
Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R.

68

（様式 20）成果報告書

TSUBAME 共同利用 平成 26 年度 産業利用 成果報告書

利用課題名 衛生陶器混相流シミュレーションの商品設計および販促への展開
英文：Multiphase fluid simulations of sanitary wares for designing and promotions
池端 昭夫
Akio Ikebata

吉田 慎也
Shinya Yoshida

TOTO 株式会社 生産技術本部
Production Technology Division, TOTO LTD.
http://www.toto.co.jp
邦文抄録（300 字程度）
平成 25 年 12 月～平成 26 年 3 月の有償利用により、TSUBAME2.5 においてあらゆる衛生陶器商品の薄膜陶器
表面流れをシミュレートできるようになり、現在新商品設計に活用している。次の課題として、洗浄時の水はねの評価を行うための技術開発を実施した。新しい保存型運動量解法を提案し、実機とは異なる「水はね」現象を抑制することができた。また流体解析技術を販売の第一線におけるお客様プレゼンツールとしても活用するための技術開発を併せて行った。レイトレーシング法により複雑な気液混相流を鮮明に可視化することで、展示会や新商品説明において顧客に製品の価値伝達をスムーズに伝えることに寄与した。
英文抄録（100 words 程度）
We are now applying our multiphase simulation code which was developed on our research in TSUBAME projects
from November, 2013 to March, 2014 to designing various products. Moreover, we have been developing theevaluation method of droplet motions at water flushing. For the sake of suppressing “spurious motions of droplets”
on the numerical solutions, we will propose new numerical method, “Conservative Momentum Method”. On theother hand, we continue to apply our simulation technology to supporting sales of our products as the presentation
tools. We contributed to sending “values of our products” to many customers with our visualization tools using
Ray-tracing method.
Keywords: multiphase flow, house hold equipment, sanitary ware, UTI-VSIAM3, CIP

こともまた求められている。

背景と目的
少子高齢化の進展により住宅着工戸数が伸び悩む

こういった課題に対し、昨年度まで TSUBAME 有償

中、住宅設備機器業界ではリフォーム分野に力を入れ利用で大規模多相流体解析コードを開発し、陶器表面てきており、取り換え需要喚起の近年の大きな販売施の薄膜流れの評価が可能となったが、商品開発に適用策の柱として、エコ・省エネの観点からの節水型商品のした段階で、洗浄時の水はねが新たな課題と判明した。

販売が拡大してきている。また国内のみならず海外に

これは節水型トイレの場合、規定の洗浄性能をクリア

おいても、地球温暖化対応や気候変動による砂漠化問するために流速が増大する傾向にあり、洗浄時の水同題、発展途上国の衛生問題の観点から、節水型住設士の衝突における水はねを発生しやすくなっているた

機器の輸出や現地生産が拡大している。しかし衛生陶めと推測される。すなわち洗浄性能と水はねは相反す

器では競争激化から、12L→6L→4L と急激に節水化る傾向にあり、双方のバランスを取った設計が必要となが進展し、衛生陶器表面の汚れ残りや衛生面での懸る。したがって洗浄時の水はね量も、事前評価で定量念などの声が上がってきており、お客様の懸念を払拭

的にシミュレーションできる技術が新たに求められるよ

するための評価技術の開発は急務となっている。実機うになった。

評価では測定のばらつきや条件設定、測定負荷等の問題により全てをまかなうことは不可能であり、シミュレ

概要

ーションによる洗浄性能評価の実現は積年の課題であ

平成 25 年 12 月～平成 26 年 3 月の有償利用により、る。さらに世界的な競争に生き残るため、これらの節水

TSUBAME2.5 に最適化した大規模並列 GPGPU クラス

技術を積極的にアピールし顧客からの信頼を獲得する

ター対応の混相流シミュレーションコードを開発し、あら69

（様式 20）成果報告書
ゆる衛生陶器商品において薄膜陶器表面流れを精度

ドは、一般的な計算手法と同様に、NS 運動方程式

よく再現できるようになった。現在当社での様々な衛生

u   uu  R.H.S t

陶器新商品設計に活用している。
次の課題として、洗浄時の水のはねの評価を行うた(1)を非保存型の定式化

u1
 u  u  R.H.S t

めの数値計算手法およびシミュレーションコードを開発
し、TSUBAME2.5 にて検証を実施する。また流体解析(2)技術を販売の第一線におけるお客様プレゼンツールとに変形した式に沿って流体計算することでロバストなしても活用するための高品位可視化処理手法の開発

シミュレーションを実現している。これは、水と空気のを併せて行い、TSUBAME2.5 での混相流シミュレーシ

ような密度比が 1000 倍も変わる混相流において、元ョン結果にて効果を確認する。

の NS 運動方程式の形で計算した場合、流速はu

結果および考察u

(3)

すなわち運動量を密度で割って計算されるため、気１． 水はねシミュレーション
昨年度までに開発した UTI-VSIAM3[1]に基づく混液界面近傍セルにおいて、式(3)の分母の密度の計

相流シミュレーションコードを 4 種類の衛生陶器製品

算誤差により、非常に大きな流速のジャンプを生じるに適用し、衛生陶器上部への水はね量を計算した。

可能性があり、非保存型解法ではこの問題を回避す

図１に製品Ａにおける実験結果との比較を示す。明らることができるからである。しかしながら、この非保存

かに水はね量が実験よりも過大であり、符合しないこ型方程式では、陽的に運動量保存性を保証していなとが判明した。

いため、数値計算上の誤差により、運動量保存性が次に、シミュレーションコードに LES 乱流モデルを大幅に悪化する可能性がある。今回の水はね不具合

追加し、再び 4 種類のモデルの計算を実施した。LES

は、この運動量保存性の誤差に起因していると推測

モデルは通常のスマゴリンスキーモデルを採用した。

される。

図２で明らかなように水はね量は実験値に近づけるこ

そこで、従来の非保存型運動量保存式解法に変わとができているが、図３に示す 4 種類の製品間の傾向

り、保存型運動量保存式解法を用いて流速を計算すが実験とは一部異なることが分かった。その原因を詳る新しい手法を開発した。式(3)における流速の不安

しく確認したところ、製品Ｂにおいて図４に示すように、定性を抑制するため、分母の運動量と分子の密度を実機では観測されない、吐水口での「水はね」が発生全く同じ移流計算定式化により更新することにより、していることが分かった。水はねの原因としては、図４

流速の安定化を図った。

の流速ボリュームレンダリング値を見ると、「水はね」

式(1)の運動量保存則を運動量について計算する

直前に赤色で示される高流速の塊が吐水口付近で際に、同時に質量保存則を密度について計算する。

発生しており、この塊が引き金となって吐水口付近でけで、進行方向ではない上方向に向かって液滴が上u   uu  0t   u   0t

昇する現象は考えにくい。そこで一つの仮説として、式(4)、(5)の左辺第一項により、運動量および密度が吐水口近傍の空気が、高流速の塊によって強く圧縮更新されるが、このままでは式(3)の問題は解消されされ、それが出口で解放される際に空気の破裂流速ない。そこで、各時刻ステップに流速と流体率から運が、液体部に数値計算誤差により拡散し、運動量保動量と密度を都度計算し、式(4)、(5)に代入して運動存則に見合わない水とびが発生していると推察した。

量と密度を更新して、式(3)により流速を求める。流速

破裂が生じたものと推測される。しかし単に高流速の塊が生じるだけでは、通常は吐水速度が増大するだ

昨年度までに開発した混相流シミュレーションコーを決めたら、運動量と密度は破棄する。

70

(4)

(5)

（様式 20）成果報告書
度の保存性については問題が生じないと考えられ計算精度とロバスト性を考慮し、式(4)、(5)を次元分る。

割型 CIP-CSL 法により計算する。UTI-VSIAM3 と同様に、物理量の VIA(=Volume Integrated Average)とまず一般的なダム崩壊流れテスト計算により、本手

SIA(=Surface Integrated Average)の 2 種類の値を用い法の妥当性を確認した。図６に計算結果を従来手法て高次精度計算を行う。

と比較して示す。計算の途中経過は従来手法同様極めて安定であり、また計算結果もほぼ同等の結果が本手法に基づく多相流体解析の具体的な計算手得られることが分かった。

順は次のとおりである。

次に、気体と液体の密度比が 106 倍と極めて大き
1) 流体率 F n より、密度 n の VIA および SIA を計算い条件における液滴運動のベンチマークテスト[3]を2) i=x,y,z について、次元分割型移流計算により保実施した。計算条件は、図７の左側から右側方向に
水平に 1.0 の流速を円形の液滴に与え、気体部分は存型 CIP 法を用いて計算
2-1)  の SIA-i を補正
n初期流速 0 として計算を実施した。計算結果は、図７2-2)u*   n u n  t  ui un2-3) *   n  tのとおり従来の非保存型運動量解法では液滴の変xi形が生じるのに対し、今回の保存型解法では変形がui nxi

生じることなくほぼ理論解どおりの結果が得られており、計算誤差による異なる相への流速の数値拡散が抑えられていることが分かった。

2-4) u*  u  /  **2-5)以上のベンチマークテストの結果を踏まえ、今回の保存型運動量解法を衛生陶器製品Ｂのシミュレーシ

Time Evolution Converting[2]を用いて SIA-j(i  j) の u を更新*

ョンに適用した。図８にシミュレーション結果を示す。

3) 流体率 F n+1 を UTI-VSIAM3 および STAA により流速分布を見ると従来手法とは大きな違いは見当たらないものの、従来手法ではどうしても抑えることがで移流計算きなかった上方への水はねが完全に抑制できており、

4) 粘性項、表面張力項、重力項および圧力項を通
n+1

常のプロジェクション法にて計算し、u実験と同様の吐水の傾向が得られた。この結果より、を更新

今回の保存型運動量解法は運動量の保存性を改善
様々なテスト計算の結果では、上記アルゴリズム

し、水の不自然な動きを解消する効果があることが分のとおり、運動量および密度の移流計算は次元分割かった。

型で計算し、一方水と空気の界面を表す流体率の移従来型手法における製品Ｂの計算結果で見られた

流計算は従来どおり非次元分割型解法で計算する方吐水口の「水はね」は、製品Ａおよび製品Ｃの計算結法が、計算の安定性の面で最も良い結果が得られ果ではほとんど生じていないことが判明した。したがた。

って今回の保存型運動量解法により、図３のＢの水

上記 2-1)の SIA の補正は、2-4)の分母の密度 が

はね量が他の製品の水はね量よりも相対的に少なく

急激に小さくならないようにするためのものである。具なることが期待される。今後、この保存型運動量解法

体的には、2-3)の密度の移流計算において、図５のよを他の製品にも同様に適用し、水はね量が実験と同

うな上流側が空気の気液界面が通過するときに、 n様の傾向を得られるか検証していく。

*の SIA-i の急激な密度低下が生じる可能性があるため、上流側の n -VIA が、空気に近いある一定の値以下になったときに、界面の  n -SIA を上流側の n-VIA の値に書き換える上流化を行う。この計算手法
では n -VIA については何も補正を行わないので、密71

（様式 20）成果報告書

図１．製品Ａの実験と解析の水とび量比較（実験値=1、乱流モデルなし）

図２．製品Ａの実験と解析の水とび量比較（実験値=1、乱流モデルあり）
図４．製品Ｂの吐水口の「水はね」
上図：上方から見た流路の流速の絶対値
下図：側方から見た水はねのようす

図３．製品Ａ～Ｄの実験と解析の水とび量比較
図５. 界面 j-1/2 において n -SIA が急激に低下する（製品Ａ=1、乱流モデルあり）

状況の例

72

（様式 20）成果報告書

図７. 液滴移動ベンチマークテスト結果（密度比 106）
左図：従来の非保存型運動量解法
右図：今回の保存型運動量解法

図６．ダム崩壊テスト計算結果
上図：従来の非保存型運動量解法
図８．非保存型運動量解法による製品Ｂの計算結果

下図：今回の保存型運動量解法

上図：上方から見た流路の流速の絶対値
下図：側方から見た水はねのようす73

（様式 20）成果報告書
ためには、水の反射や屈折による輝度の変化を適切

２． 混相流シミュレーションのコミュニケーションツール
への適用に表現する必要があるからである。水の流れの実験

TSUBAME により非常に細かいメッシュサイズのシを高速度カメラで動画撮影するときにストロボやレフ

ミュレーションが可能となり、実製品の流れの傾向を

板の位置を細かく調整するのと全く同様であり、可視良く表せるようになってきたことから、TSUBAME のシ化という点において両者に違いはないといえる。これ

ミュレーションを実物に変わる代替可視化ツールとしまで当社のレンダリングでは、色がない水の識別を際ても利用できるようになってきた。衛生陶器の流れの立たせるために背景に様々な物体を置いて、その背場合、水の流れは至るところで渦や乱流が発生し複景像が屈折することで流れを可視化するという一般

雑な流れ場となるが、これらの流れを分かりやすく可的な手法を用いていたが、今回の高品位レンダリング

視化することで、製品ごとの流れの特徴を第三者にでは、照明の調整を最適化することにより、背景の物

明確に伝えることができる。このためには、体を配置することなく水の動きを明確化することを目①流体シミュレーションの精度が十分高いこと指した。

②シミュレーション結果の可視化精度が十分高いこと

図９に当社の衛生陶器海外市場向け新製品におの２点の条件が必要である。これまでの検討で①はける TSUBAME シミュレーション結果を可視化レンダ

かなり達成できるようになってきたことから、今回②にリングにより作成した動画の一例を示す。本シミュレ

ついても検討を行う。

ーションは、0.8mm メッシュを使用している。本動画に流体シミュレーション結果の可視化精度とは、各時より、当社が「トルネード洗浄」とよんでいる螺旋状の刻における水や空気、物体などのシーンにおける各渦流れが、大小無数の気泡の動きにより的確に捉え

構造体に対して行う照明シミュレーションの精度に相られ、螺旋状の渦が固体を巻き込むように排出してい当する。水の場合は透過体であるため、水の表面にく過程が分かりやすく示されている。このことにより、到達した光線は屈折および反射を繰り返し、複雑な当社製品独自の流れである「トルネード洗浄」の商品

光線経路を辿る。この複雑な計算を実施する手法とし価値を分かりやすく顧客に伝えることが可能となった。

て、レイトレーシングが広く用いられている。レイトレー

図１０には、内部流れを微小のトレーサーにより可視

シングは一般的に CG の手法として認識されているが、化した例を示す。トレーサーは非常に軽いと仮定し、本来は忠実に光線追跡を実行し、実写と同様の画像

水の各時刻、各座標における流速ベクトルに沿ってを得る照明シミュレーション手法である。

並進させ、その流速ベクトルの rot u    u を計算

これまで TSUBAME では GPGPU 並列処理を用いし粒子の各軸方向の三次元角速度を計算して回転さたレイトレーシングプログラムが開発された事例[4]がせる方法を用いた。この粒子の動きを追跡することにあるが、今回は水の高品位レンダリングを目的とする

より、時々刻々変化する内部の大まかな三次元流れため汎用ツール POV-Ray Ver3.7 を利用した。Ver3.7

構造を把握することができる。

では、処理ピクセルごとに光線追跡深度に大きなばら

また図１１には、他社製品との流れの違いを、流速つきを生じる水のレンダリングをマルチスレッド並列ベクトルとともに比較したものを示す。他社製品は両処理により効率的に処理可能である。ただし側から水が流れるため中心で両者がぶつかり、エネ

POV-Ray 自体に水を忠実に表示するマテリアルは用ルギーロスを起こしているが、当社製品は均一な反意されていないため、水の屈折率は実物に近い値を時計回りの流れを生成するため、無駄なエネルギー

使用する一方、拡散反射係数、鏡面反射係数などのロスがなくしっかり陶器表面を洗い流せることが明確

各物体表面物性は、レンダリング結果を見ながら細になった。

かく調整していった。またもう一つの重要な点として照このようにして作成された可視化動画を元に映像明の強さや位置、個数の最適化も行った。これは、水用に加工されたものは、海外大規模展示会 KBIS（米、自体には本来色がないため、水の存在を可視化する

2015.1）および ISH（独、2015.3）にて展示され、来場74

（様式 20）成果報告書
者に当社製品の独自性と価値を伝達することに寄与した。

図１０．海外市場向け衛生陶器製品流れ（トレーサー微粒による内部流れの可視化）

図１１．海外市場向け衛生陶器製品流れ（他社製品との流れの比較、上：当社、下：他社）

図９．海外市場向け衛生陶器製品流れ（固体を含む内部の流れの可視化）

75

（様式 20）成果報告書
スケールへの飛翔） 成果報告書 (2010)

まとめ、今後の課題
衛生陶器の流体シミュレーションは、TSUBAME 上での高精度な流体解析により、以前の PC ベースのシミュレーションよりも画期的に適用範囲が拡大しており、これまではシミュレーション精度や可視化品質の面で活用できなかった洗浄時の水はね評価や営業支援に応用範囲を拡大している。本年度の検討では、水はね
評価について本来実機では発生しない数値計算上の「水はね」という課題が明らかになり、保存型運動量解法に基づく新しい計算手法の提案および検証を行い、課題の解決を確認した。また営業支援については、
TSUBAME の高精細なシミュレーション結果を活用して、従来のポストプロセッサと比較して内部の複雑流れの大幅に鮮明な可視化を達成した。
今後は保存型運動量解法による新しいシミュレーションコードを他の衛生陶器製品にも適用し、水はね量の傾向が実機と一致するか検証していく。また可視化については様々な製品の営業支援に展開していく予定で
ある。

参考文献
[1] A. Ikebata, Y. Muraoka, F. Xiao: Multiphase fluid
simulations on a multiple GPGPU PC using unsplit
time integration VSIAM3, Progress in NUCLEARSCIENCE and TECHNOLOGY, 2, pp.491-497(2011)[2] F. Xiao: Unified formulation for compressible
and incompressible flows by using multi-integrated
moments I: one-dimensional inviscid compressible
flow, Journal of computational physics, 195,
pp.629-654 (2004)[3] M. Raessi, H. Pitsch: Consistent mass and momentum
transport for simulating incompressible interfacial
flows with large density ratios using the level set
method, Computers & Fluids, 63, pp.70-81(2012)[4] P. V. PHUC: 建築物の室内外環境の連成解析とその高速化技術の開発, 共同利用（産業利用トライアルユース：先端研究施設共用促進事
業『みんなのスパコン』TSUBAME によるペタ

76

（様式 20）成果報告書

TSUBAME 共同利用 平成 26 年度 産業利用 成果報告書

利用課題名 塗工スラリーの分子シミュレーション
英文：Molecular Dynamic Simulation of Slurry Coating Process
諸星 圭
Kei Morohoshi
トヨタ自動車株式会社
Toyota Motor Corporation
http://www.toyota.co.jp/
スラリー塗工プロセスを分子レベルでシミュレーションする手法を開発し、TSUBAME と粗視化分子動力学法を用いた大規模計算を行うことで材料開発への活用を可能にした。溶媒の蒸散モデルと多孔質構造中の拡散評価法を開発して塗工プロセスを計算した。溶媒蒸散に伴い微粒子と高分子の構造化が進み、多孔質構造体が形成される様子を再現することができた。また、塗工に用いる分子が塗工後の特性に及ぼす影響についても調べた。
Slurry coating process has been investigated by using coarse-grained molecular dynamics
simulations. Supercomputer TSUBAME enables to deal with large-scale calculations, and that canrealize a design of molecular structure. We reproduce the process that some fillers and polymer chainsare gradually gathered during evaporation and finally form to complicated porous structure.

Keywords: electrode structure, evaporation, coarse-grained molecular dynamics, bead-spring model

背景と目的
スラリーの塗工は、燃料電池やリチウムイオン電池

材料開発としては、微粒子の表面改質や溶質高分の電極作製などで幅広く行われている工業プロセスで、子の分子構造変更により構造体としての特性向上を図微粒子と溶質を溶媒に混ぜたスラリーを基板の上に塗りたいが、そういった分子レベルの情報がどう構造形成

布した後に溶媒を蒸発させて微粒子と溶質からなる多に作用するのかを直接観測することは難しく、また個々

孔質構造を得るプロセスである（図 1）。構造体の形状の分子を表現してシミュレーションを行うには規模が大によって電極特性が変わるので、特性を高めるためのきすぎる。そこで、TSUBAME を利用することにより規改善が生産技術と材料開発の両面から行われている。

模の問題を解決し、分子モデルでの動力学シミュレー

しかしその方法は経験的・試行錯誤的であることが多く、ションを実施して塗工プロセスで起きている現象を解明

特に材料開発面では用いる材料の組成と得られる構することを目指した。ただし、TSUBAME を使用しても造体の特性との関係がよくわかっていない。

個々の原子をあらわに扱うことは困難であるため、複数の原子をまとめてひとつの粒子として扱う粗視化分
子動力学法を用いることにした。

溶媒蒸発

概要
計算には粗視化手法のひとつであるビーズ‐スプリン
グモデルを用いた[1]。スラリー中の微粒子・溶質高分
子・溶媒低分子は図 2 のようにモデル化し、微粒子を構成する粒子間には結合を設定して球形状を保つようにした。非結合粒子間の相互作用には Lennard-Jones
図 1 スラリー塗工プロセス

ポテンシャルを用いた。

緑：微粒子、赤：溶質高分子、青：溶媒低分子

77

（様式 20）成果報告書
𝜎𝑖𝑗 12𝜎𝑖𝑗 6)−() } + 𝜀LJ , 𝑟 ≤ 𝑟𝑐
𝑈LJ (𝑟) = {𝑟
𝑟
0,𝑟 > 𝑟𝑐
ここで𝑟は粒子間距離、𝜀LJ は𝑟 = 𝑟𝑐 でポテンシャルを連

電池の水素イオンなら高分子中を、リチウムイオン電

続にするための定数である。本研究では全ての粒子でを評価しなければならない。そこで、蒸散後の粒子分布

𝜎𝑖𝑗 = 𝜎、𝑟𝑐 = 2.5𝜎とした。結合粒子間には FENE ポテ

から微粒子・高分子・空孔の 3 次元メッシュを作成し（図ンシャルを追加した[1]。

4）、メッシュ毎に拡散方程式を解いて特性を評価するこ

4𝜀𝑖𝑗 {(

池のリチウムイオンなら電解液部分になる空孔中を、また電子なら微粒子中を移動するため、経路毎に拡散性とにした[3,4]。

高分子の分子構造を変えると鎖の剛直性が変わることがあり、材料開発では重要な特性のひとつである。
そのため、連続して結合している 3 粒子間に調和振動
子型の曲げポテンシャルを適用して剛直性を表現した。
𝑈bend (𝜃) =1𝑘 (𝜃 − 𝜃0 )22 𝑖

ここで𝜃は結合角度の外角、𝑘𝑖 は曲げ弾性定数、平衡
角𝜃0 は 0 とした。

図 2 微粒子・溶質高分子・溶媒低分子のモデル

溶媒の蒸散は、気相と液相が共存するときの気相領域の端から粒子を削除することでモデル化し、削除す
る時間間隔を変えることで蒸散速度の違いを表現できるようにした[2]。溶媒粒子間に働く引力が弱いと液体としての凝集を再現できないため、図 3 に示すように引力

図 4 拡散経路の抽出の大きさを変えて気液界面を保持できるか調べた結果

緑：微粒子、赤：高分子、青：空孔から、ε = 1.0を用いることとした。
モデリングには材料物性解析ソフトウェア J-OCTA を[5]、計算には並列化分子動力学エンジン VSOP を[5]、拡散評価には衝撃・構造解析ソフトウェア LS-DYNA
を[6]、それぞれ用いて行った。

結果および考察
スラリー塗工プロセスの計算結果を図 5 に示す。溶媒が蒸散していき、それに伴って微粒子と高分子の構
図 3 気液界面の形成

造化が起こり、複雑な構造体が形成される様子が再現できている。最終構造が多孔質形状となっているか調

蒸散後の構造体として重要となる輸送特性は、燃料べるために、隙間に球を詰めていくことにより空孔の分

78

（様式 20）成果報告書
布を解析した。結果を図 6 に示す。空孔が構造体全域に存在し、大きなものでは微粒子と同程度の大きさをもつ多孔質構造であることがわかった。

図 7 高分子の分布

まとめ、今後の課題
スラリーの塗工プロセスを分子レベルでシミュレーシ
図 5 スラリー塗工プロセスの計算

ョンする手法を開発し、TSUBAME と粗視化分子動力学法を用いることで初めて材料開発への活用を可能にした。溶媒の蒸散に伴い微粒子と高分子が構造化して多孔質構造が形成される様子を再現した。溶質高分子
が輸送特性に与える影響について調べたところ、剛直性が高分子の経路を細くして拡散性を低下させることがわかった。今後は微粒子の表面改質や蒸散速度の違いによる影響について検討する。

参考文献[1] K. Kremer ,G. S. Grest, J. Chem. Phys., 92,

図 6 空孔の分布5057 (1990)
材料開発の手段のひとつとして溶質高分子の分子[2] M. Tsige, T. R. Mattsson, G. S. Grest,

構造を変更し、鎖の剛直性が変わった場合の影響につMacromolecules, 37, 9132 (2004)

いて調べた。柔軟な高分子では曲げ弾性定数𝑘を 0 に、

[3] G. Dorenbos, V. A. Pomogaev, M. Takigawa,

剛直な高分子では 1.5 に設定して計算した。蒸散後のK. Morohoshi, Electrochemistry

構造体における高分子中の拡散を評価すると、柔軟なCommunications, 12, 125 (2010)[4] G. Dorenbos, K. Morohoshi, Energy Environ.

高分子では 0.46 なのに対して剛直な高分子では 0.22と約半分であった。拡散評価に用いた高分子のメッシュ

Sci., 3, 1326 (2010)

を見ると微粒子を被覆している高分子が薄くなっており[5] http://www.j-octa.com/jp/（図 7）、また空孔分布を解析して空孔率を求めると[6] http://cae.jsol.co.jp/product/struct/lsdyna/

0.53 から 0.51 に減少していることから、高分子が剛直になると微粒子への吸着量が減って一部が空孔中で凝集することがわかった。従って輸送特性を考えるときに、剛直性は高分子の経路を細くして拡散性を低下させていることが明らかになった。

79

80

（様式 20）成果報告書

TSUBAME 共同利用 平成 26 年度 産業利用 成果報告書

ワイヤレス電力伝送による漏えい電波の環境解析技術の研究開発
Research and development of environment analysis technique of leaked electromagnetic field fromwireless power transfer
池田 和彦
Kazuhiko Ikeda(株)パナソニックシステムネットワークス開発研究所
Panasonic System Networks R&D Lab. Co., Ltd.
http://panasonic.co.jp/avc/psnrd/
家庭用電子機器などの給電に用いられるワイヤレス電力伝送(WPT)システムが近年検討されている．WPT シス
テムから漏えいする電磁界は，他の電子機器の動作に影響を及ぼす可能性があるため，漏えい電磁界による干渉影響の把握が課題であり，様々な設置環境や周波数帯における定量的な評価が必要となる．WPT システムでは，
数次高調波における漏えい電磁界も発生することから，GHz 帯を使用する無線機器への干渉も懸念されるため，
GHz 帯の漏えい電磁界評価も必要である．そこで本研究開発では，これまでに筆者らが開発した Ray Launching法による電波伝搬シミュレータを用いて，戸建て実験住宅内の WPT システムからの漏えい電磁界特性を解析し，測定結果との比較検証を行う．また，TSUBAME 2.5 を用いた並列計算により，実用的な時間内に大規模な空間における漏えい電磁界特性の解析が可能となることを示す．

The wireless power transfer (WPT) system has been actively developed in recent years. Since theelectromagnetic field leaked from the WPT system may interfere with the other devices that use theGHz band, it is required to quantitatively evaluate the electromagnetic interferences over a wide
range of frequencies in various installation conditions. In this paper, we propose an approach foranalysis the electromagnetic field of the high-harmonic component or the microwave WPT system inthe case-study house by simulation using the Ray Launching method. The validity of the proposed
simulation method is confirmed by comparing the simulation and the measurement results. Moreover,
we also show that the electromagnetic field properties in a huge space can be analyzed within apractical time using the supercomputer TSUBAME 2.5.

Keywords: 無線電力伝送, 漏えい電磁界, Ray Launching

することが必要となる．しかしながら，数十 kHz～数

背景と目的
近年，電波の利用方法としては，通信だけではなく，GHz の周波数範囲において漏えい電波の様々な特電力伝送技術が注目されている．2015 年以降の実用性をシミュレーションにより評価する技術はこれまで確化が予定されているワイヤレス電力伝送（WPT）システ

立 さ れ て い な い ． 漏 え い 電 波 の 特 性 をムとして，電気自動車用の大電力（数 kW）を伝送する

FDTD(Finite-Difference Time-Domain)法などの電ものや，家電機器用の中電力（数百 W）を伝送するも磁界シミュレーション技術により解析する場合は，計算のなど，多種多様な方式開発が進められている．

規模は波長に依存する．例えば 3.5GHz 帯や 5.3GHz

WPT システムは，一般家庭や集合住宅への普及が帯などの周波数帯では，住宅内またはそれ以上の空見込まれており，使用範囲として住宅内に限らず屋外間で計算規模が膨大となり，実用的な時間内に解析すの駐車場も含めた広い空間が想定されている．上記のることが困難となる．

ような環境において，今後普及が進む WPT システム

本研究開発では，WPT システムをはじめとする各と他の機器との共存を図るためには，WPT システム

種電子機器が密集して設置された環境における漏えいから発生する数十 kHz～数 GHz の広い周波数範囲に電波の特性を解析できるシミュレーション環境の構築をわたる漏えい電波が近接する機器に与える影響を解析目的として，3.5GHz 帯と 5.3GHz 帯の漏えい電磁界を81（様式 20）成果報告書
対象に次に示す検証を行う．幾何光学的解析手法に基づく Ray Launching 法電波伝搬シミュレータを開発し，戸建て実験住宅内における WPT システム漏えい電磁界の測定結果とシミュレーション結果の差分±5dB 以内を目標に比較検証を行うことで，本シミュレーションの妥当性を示す．また，本シミュレータに対し，スーパーコ
ンピュータ TSUBAME 2.5 を用いた並列計算の適用
図 2 TSUBAME 2.5 を用いた並列計算方法

により，戸建て実験住宅周辺の構造物を含む大規模な空間において，実用時間内に漏えい電磁界特性の解

戸建て実験住宅におけるシミュレーション条件

析が可能となることを示す．

図 3 及び図 4 に実験住宅モデルを示す．また，図 4
Ray Launching 法の概要に示すように，住宅周辺の建物からの反射等により住Ray Launching 法を用いた幾何光学伝搬解析手法

宅内の漏えい電磁界電界特性が影響を受けるため，は，送信点から放射される電波を図 1 に示すように電150m 四方の構造物もモデル化した．WPT システムの波を光線(Ray)とみなし，送信点から離散的に放射した設置位置（送信位置）は実験住宅内の一階リビングとし，Ray が受信点に到達するまでの軌跡を求めることで送伝送距離が数 m 程度のマイクロ波帯を利用するシステ

受信間の電波伝搬特性を推定する解析手法である．こムを想定して図 5 に示す放射パターンを有する指向性の解析手法は，計算量が送信点から放射される Ray

アンテナを用いる．受信機器には被干渉機器 として間の離散角度に依存する．そのため，大規模な空間に

3.5GHz 帯は携帯電話端末，5.3GHz 帯は無線 LAN

おいても，他の解析手法(例えば FDTD 法)と比較してルータを想定し，図 6 に示す放射パターンを有するア

計算量を抑えられる利点がある．

ンテナを用いる．測定結果とシミュレーション結果の比Ray Launching 法の特性として，放射した各 Ray が較は，見通し内外を含む 4 箇所について，漏えい電磁

受信点に到達するまでの軌跡計算を独立に行えること

界の受信電力の中央値を比較することで行う．また計が挙げられる．本特性を利用し，各 Ray について並列

算条件は反射回数 10，透過回数 10，回折回数 2，放計算を行うことにより，計算時間を短縮することができ

射離散角度 1 度とする．

る．図 2 に，本シミュレーションで採用した並列計算方法を示す．TSUBAME 2.5 は複数の計算ノードが高速通信で接続された構成であるため，送信点から放射される Ray を方向毎に分割し，各計算ノードに配分することで複数の計算ノードを用いた並列計算を行う．

図 3 実験住宅モデル
図 1 Ray Launching 法82

（様式 20）成果報告書
シミュレーション結果と測定結果の比較検証
3.5GHz 帯及び 5.3GHz 帯について，二階子供部屋
における測定結果とシミュレーション結果の 30cm 区間
毎の中央値を比較した結果を図 7 に示す．図 7 中のSim.model(a)は，参考文献[3]に基づいて電気定数を
実験住宅内の各壁面に適用した結果であり，Sim.model(b)は外壁の電気定数を変更した結果である．Sim.model(a)は，一部区間で測定結果との差分が大きく，測定結果と差異があることがわかる．

【考察】
図 4 実験住宅周辺の構造物

Sim.model(a)における差分要因として，住宅内の部
材のうち，外壁を経由して到来する伝搬路に着目する．
外壁経由の伝搬路は図 8 に示すように屋内における外壁反射波と屋外からの外壁透過波が到来することをシミュレーション結果より解析した．外壁は内部に断熱材や内壁側の石膏ボード等の複数材質によって構成され，構造が複雑である．従って，外壁の電気定数が測定環境と異なり，透過時または反射時の減衰量の差異
が差分要因と推測した．
そこで外壁の電気定数を測定環境と同等に調整するため，実験住宅の屋外において WPT システムからの漏えい電波を測定し，そのパスロスから外壁における透過波の減衰量を推定した．表 1 は外壁透過波の減衰量から推定した電気定数の変更前後を比較したものである．電気定数変更後の Sim.model(b)は透過時の

図 5 WPT システムの放射パターン

減衰量が電気定数変更前に対し約 5dB 低下し，同様に反射時の減衰量は電気定数変更前に対し約 8dB 低下することがわかる．

図 7 二階子供部屋における測定結果とシミュレーショ
図 6 受信機器の放射パターン

ン結果の比較

83

（様式 20）成果報告書

図 8 着目した二階子供部屋に到来する伝搬路

表 1 外壁の電気定数変更前後の比較
外壁の

比誘

導電率

反射時の

透過時の

電気定数

電率

[S/m]

減衰量

減衰量

Sim.model(a)2.8

5×10-29.4dB

12.2dB

Sim.model(b)10.9

1×10-2

1.5dB

7.5dB

【結果】
前記考察と同様に実験住宅内の各壁面(内壁, ドア,

図 10 各測定箇所における測定結果とシミュレーショ

窓等)について電気定数を変更した結果を図 9 及び図

ン結果の比較 (5.3GHz 帯)

10 に示す．図 9 に示す 3.5GHz 帯及び図 10 に示す
計算速度向上効果検証

5.3GHz 帯の共に測定結果とシミュレーション結果の差分が±5dB 以内となることを確認できる．以上より，壁

前節で述べたように，本シミュレーションではスーパ

面の電気定数を適切に設定することで測定結果とほぼーコンピュータ TSUBAME 2.5 について複数の計算ノ一致するシミュレーション環境を構築した．

ードを用いた並列計算が可能である．並列計算による計算速度の向上効果を検証するため，表 2 に示す計
測条件に従って計算時間を計測した． 表 3 は計算時間の計測結果であり，図 11 は計算ノード数に対する計算速度向上度を示したものである．なお，計算速度向上度は 1 計算ノード時の計算時間に対する計算時間の短縮割合を示し，下記の式より計算する．
計算速度向上度 =計算ノード数 1 の計算時間
計算ノード数 n の計算時間(式 1)

表 3 に示すように，複数計算ノードを用いた並列計
算により約 20 倍の計算速度向上効果が得られ，1 時間以内に住宅モデルの漏えい電磁界解析が可能となった．一方，図 11 が示すように，計算ノード数 30 以上において計算速度向上度が飽和することがわかった．

【考察】

図 9 各測定箇所における測定結果とシミュレーション

所定以上の計算ノード数において計算速度向上度

結果の比較 (3.5GHz 帯)が飽和する要因として，計算ノード毎の計算量に偏りがあることが考えられる．前節で述べたように本シミュレーションの並列計算は方向毎に分割した放射 Ray を各

84

（様式 20）成果報告書
計算ノードに配分する方法を採用している．例えば，特定の方向に遮蔽物となる壁面が多い場合は，該当方向に配分された計算ノードの計算量が増加し，計算ノード間で計算量に偏りが生じることが考えられる．このような場合，特定の計算ノードの計算量がボトルネックとなり，計算ノード数を増加した場合においても所望の計算速度向上効果を得られない．
計算ノード間の計算量の偏りを改善する方法として，各計算ノードに配分した放射 Ray を定期的に再配分する方法が考えられる．改善案では，放射 Ray の再配分

図 11 計算ノードに対する計算速度向上度の際に，計算ノード間に放射 Ray に関する情報のデー
タ転送が発生するが，TSUBAME 2.5 の計算ノード間

まとめは表 2 中に示すように最大 80Gbps の高速通信が可

Ray Launching 法シミュレーションにより戸建て実

能な Infiniband を用いて接続されているため，全体の験住宅内の WPT 漏えい電磁界特性を解析し，住宅内計算時間に対するデータ転送時間の占める割合は少部材のうち壁面の電気定数を適切に設定することで測ないと考えられ，計算ノード間の計算量偏りの解消によ定結果とシミュレーション結果の差分が±5dB 以内となり計算速度の向上が期待できる．

り，シミュレーションの妥当性を示した．
また，TSUBAME 2.5 を用い，Ray Launching 法に

表 2 計算時間の計測条件
値

並列計算を行うことにより，実験住宅内の WPT システ

Intel Xeon X5670 2.93GHz,

ム漏えい電磁界特性解析を 1 時間以内で計算が可能

項目
CPU

おける放射 Ray 毎の軌跡計算を複数の計算ノード間でなシミュレーションシステムを構築した．

6cores x 2slot, with Hyper Threading
54 GByte

今後は，戸建て実験住宅以外の住宅モデルにおけ計算ノード間接続

QDR Infiniband x2 (80Gbps)

る検証を通してのシミュレーション精度の向上と，大型計算ノードあたりの

1 プロセス

店舗等の大規模モデル解析に向けて計算ノード間の計

搭載メモリ

算量の偏り解消による計算速度改善に取り組む．

MPI プロセス数

なお，本研究は総務省平成 26 年度電波資源拡大の24 スレッド

計算ノードあたりのための研究開発「ワイヤレス電力伝送による漏えい電

マルチスレッド数
最大使用メモリ

48 GByte

計算ノード数

1～100 ノード

周波数帯1) 3.5GHz 帯,

波の環境解析技術の研究開発」の一部である．

参考文献[1] 篠原 真毅, 電気学会論文誌Ｂ, Vol.J96-B ,No.9,2) 5.3GHz 帯
MPI バージョン

pp.881-893, “無線電力伝送の送電距離に対する理

1.1

論と技術”

表 3 計算時間の計測結果
周波数帯[2] 高 木 他 , 信 学 技 報 ,AP2003-159,RCS2003-165,

計算時間
ノード数 1

ノード数 30

3.5GHz 帯

21.36 hour0.97 hour

5.3GHz 帯

11.54 hour0.68 hour"FDTD 法と Ray-Launching 法による携帯機アンテナの実効利得評価"[3] 近藤孝亮他,信学技報,EMCJ2005-109 (2005-11),“室内における建築材料電気特性モデルベース測定法の測定条件についての実験的検討”

85

86

（様式 20）成果報告書

TSUBAME 共同利用 平成 26 年度 産業利用 成果報告書

利用課題名 大容量データ伝送用ミリ波アンテナのレドームに関する基礎検討
英文：A Study on radome for millimeter-wave antenna
利用課題責任者 千葉 修二
Shuji Chiba所属 スタッフ株式会社
STAF corporation．
URL http://www.staf.co.jp/
邦文抄録（300 字程度）
屋外で使用される大容量データ伝送用ミリ波帯アンテナに装着されるレドームの厚みに注目し検討を行った。レドームの厚みは電気的特性、機構的強度を考慮し 2mm~3mm に設定した。TSUBAME を用いた電磁界解析で、レドームの厚みを細かく可変させた時のアンテナの指向性ならびにアンテナ利得を示している。その結果をも元にレドームの厚みの最適値について示した。
英文抄録（100 words 程度）
We studied the thickness of the radome in the millimeter-wave band antenna for outdoor. The
thickness of the radome, in consideration of electrical characteristics and mechanical strength, was set
from 3 mm to 2 mm. Shows an antenna directivity and antenna gain obtained when varied thethickness of the finely radome in electromagnetic analysis using TSUBAME. We shows the optimum
value of the thickness of the radome the result based on.

Keywords: radome, electromagnetic simulation, horn antenna,

背景と目的している。レドームのパラメータとしては、厚み、誘電率、

近年モバイル端末の普及、およびコンテンツ情報量

形状など様々あるが、本課題においては、レドームのの増大に伴い、通信ネットワークトラフィックが爆発的に厚みに対する指向性の変化と利得の変化に的を絞り、増大しているなか、ミリ波帯 E-Band（71GHz～86GH

基礎的なデータの取得を行った。

ｚ）を使用した大容量のデータ伝送が可能な無線装置が注目されている。使用されるアンテナにおいては、通概要

信距離の確保、通信の安定化などを実現する為にアン

本課題は、基礎的なデータの取得を目的としている

テナに求められる性能は、高利得化、低サイドローブ化

為、検討アンテナはホーンアンテナを採用している。そなど極めて高く、また屋外で使用するため、気象の悪条の構成を図 1 に示す。

件も考慮したものが必要となる。気象の悪条件からアレドームの影響を確認する為に電磁界解析実施モデ

ンテナを保護する役割でレドームを装着するが、レドー

ルは、レドームなしモデル図 1 を基本状態とし、次にレムに求められる性能は電波の透過率が高く、機械的にドームありモデル図 2 としている。

強く、耐候性に優れていることが求められる。ミリ波帯

レドームは天面の厚みを機構的な強度を加味しアンテナにおいてレドームの影響は大きいことが予想さ

2mm~3mm とし 0.1mm 刻みで可変（図 3）、電磁界解れ、高いアンテナ性能を維持しつつ優れた耐候特性を析を実施し、指向性と利得をレドームなし状態とレド確保するためには、レドームの最適化は極めて重要でームの厚みｔを可変したものと比較する事で、ミリ波あり大きな課題となる。

E-Band 帯の 71GHz と 76GHz におけるレドームの影響を確認する。

上記課題に対し電磁界解析を用いてミリ波帯アンテナにおけるレドームの影響を明らかにすることを目的と

87

（様式 20）成果報告書
開口部程度悪くなり、2.4mm と 2.8mm でレドームなし状態と同等の利得となっている。
71GHz と 76GHz では利得の変化傾向の違いが見られた。これは 71GHz の波長λ＝4.23mm に対し76GHz では波長λ＝3.95mm となりこの差が影響していると考える。
71GHｚの利得とレドーム厚みの関係

給電部

レドームなし19.0

図 1 ホーンアンテナ：角すいホーン

レドーム厚み可変
17.0/レドームなしモデル基本状態

利得【dBi】

15.0

レドーム：比誘電率 4.35
210mm

13.0

11.0
9.0

図４（a）

7.0
5.0
1.9

2

2.1

2.2

2.3

2.4

2.5

2.6

2.7

2.8

2.9

3

3.1

レドーム天面の厚み【mm】

図４（a）
図 2 レドームありモデル

レドームなし76GHｚの利得とレドーム厚みの関係

レドーム厚み可変19.0
17.0

利得【dBi】

15.0
13.0
11.0
9.0
7.0

16.7mm

レドーム天面厚みｔ

5.0
1.9

2

2.1

2.2

2.0mm~3.0mm

2.3

2.4

2.5

2.6

2.7

2.8

2.9

3

3.1

レドーム天面の厚み【mm】

0.1mm 刻みで可変

図４（b）
次に図 5（a）（b）に 71GHｚの指向性とレドーム厚みの関係を示す。グラフにはレドームなし状態図 1 と、利図 3 レドーム天面厚みｔ可変

得と厚みの関係から利得の変化が大きかった 2.1mm、2.8mm と変化が少なかった 2mm、2.4mm、3mm を比結果および考察較したものとなっている。H 面の結果からは、レドーム

図 4（a）（b）に図 1 と図 3 の電磁界解析の結果としての厚み 2.4mm、2.8mm において指向性の変化が大き

71GHz、76GHｚの利得とレドーム厚みの関係を示す。

く、一方 E 面は、2.4mm において 0°～10°付近で指厚みを可変する事で利得が変化している事が解る。

向性の変化が見られたが、それ以外の所では厚みに71GHｚの利得は、レドームなし状態を基準とすると関係なく変化量は少ない結果となった。今回の条件に2.1mm で 1dB 悪く、2.8mm で 1.5dB 程度良くなり、

おいて 71GHｚではレドームの厚みを 2mm または

2mm、2.4mm、3mm でレドームなし状態と同等の利3mm にする事でレドームによる影響を最小限に抑える得となっている。76GHｚの利得はレドームなし状態をことが可能と考える。

基準とすると 2.1mm で 2dB 程度良く、2.6mm で 2dB

88

（様式 20）成果報告書

レドームの厚みとH面の指向性の変化

76GHz
20

15

15

10

10

5

5

利得【dBi】

利得【dBi】

71GHz
20

0

レドーム厚み2.0mm-5

レドーム厚み2.1mm-10

レドーム厚み2.4mm

-15

0
-5

レドーム厚み2.0mm-10

レドーム厚み2.1mm
レドーム厚み2.4mm-15

レドーム厚み2.8mm

レドームの厚みとH面の指向性の変化

レドーム厚い2.6mm-20

レドーム厚み3.0mm-20

レドーム厚み2.8mm-25

レドームなし-25

レドームなし

-30

-30

0

10

20

30

40
50
角度【deg】

60

70

80

90

0

10

20

30

76GHz
20

15

15

10

10

5

5
利得【dBi】

利得【dBi】

レドームの厚みとE面の指向性の変化

200
レドーム厚み2.0mm-5

レドーム厚み2.1mm-10

レドーム厚み2.4mm

-15

60

70

80

90

80

90

図 6（a）

図 5（a）
71GHz

40
50
角度【deg】

レドームの厚みとE面の指向性の変化

0-5

レドーム厚み2.0mm-10

レドーム厚み2.1mm

レドーム厚み2.4mm

レドーム厚み2.8mm

-15

-20

レドーム厚み3.0mm-20

レドーム厚み2.8mm-25

レドームなし-25

レドームなし-30

レドーム厚み2.6mm

-30

0

10

20

30

40
50
角度【deg】

60

70

80

90

0

10

20

30

40
50
角度【deg】

6070

図 6（d）

図 5（b）
まとめ、今後の課題
次に図 6（a）（b）に 76GHz の指向性とレドーム厚み

ミリ波アンテナにおけるレドームの厚みの影響を電の関係を示す。グラフにはレドームなし状態図 1 と、利磁界解析にて確認した。レドームの厚みを可変すること

得と厚みの関係から利得の変化の大きかった 2.0mmで、利得、指向性が変化する事を確認できた。また、こと 2.1mm、2.6mm、利得の変化が少なかった 2.4mmの結果から影響の少ないレドームの厚みが確認でき、と 2.8mm、を比較したものとなっている。H 面の結果は今回の条件においてはレドームの厚みを 2mm にする

厚み 2.0mm が最も指向性の変化が少なく、2.4mm で事でレドームの影響を最小限に抑えられることが解っ

指向性の変化が大きい結果となっている。一方 E 面のた。

結果はどの厚みにおいても指向性の変化が少ない結今回は 71GHz と 76GHz に周波数を絞ってレドーム果となっている。今回の条件において 76GHz では、レの影響を確認した が、 今 後は 周波数の 範囲を広 げ

ドームの厚みを 2.0mm にする事でレドームの影響を最（81GHz~86GHz）レドームの影響を確認し、厚み以外小限に抑える事が可能と考える。

のパラメータ（比誘電率、レドーム形状など）の影響も確認する予定である。

89文部科学省 研究開発施設共用等促進費補助金
先端研究基盤共用・プラットフォーム形成事業

東京工業大学 TSUBAME 産業利用
平成 26 年度利用終了課題 利用成果報告書

平成 28 年 2 月 29 日発行
国立大学法人 東京工業大学
学術国際情報センター 共同利用推進室〒152-8550 東京都目黒区大岡山 2-12-1 E2-6
URL http://www.gsic.titech.ac.jp/
E-mail kyoyo@gsic.titech.ac.jp
本書に記載の記事・写真等の二次使用を禁じます。これらの情報は著作権法上認められた「私的利用」または「引用」の条件をみたした場合を除いて、著作権者に無断で転載、複製、放送、公衆送信、翻訳、販売、貸与等の利用を禁じます。