8マルチGPU 計 算による衛 生 陶 器のための
気 液 二 相 流シミュレーション
Numerical simulation of air/water multiphase flows
for ceramic sanitary ware design by multiple GPUs

Liグラファイト層 間 化 合 物におけるLi 拡 散に対するハイブリッド量 子 古 典シミュレーション
A Hybrid Quantum-Classical Simulation study on the Li Diffusion
in Li-Graphite Intercalation Compounds

建 築 分 野の数 値 流 体 解 析における大 規 模 計 算
A Large-scale Simulation on CFD in Construction Industry

大 規 模 地 震における強 震 動 評 価と屋 内 収 容 物の被 害 評 価Estimation of strong ground motion and damage of indoor objects during giant earthquake

マルチGPU計算による衛生陶器のための気液二相流シミュレーション
池端 昭夫 *

吉田 慎也 ** ＴＯＴＯ株式会社・生産技術センター

肖 鋒 **** 東京工業大学・大学院総合理工学研究科

衛 生 陶 器など住 宅 水まわり商 品の開 発において、様々な設 計 課 題に対 応するため気 液 二 相 流 体 解 析コードを自 社 開 発
している 。空 気と水の複 雑な界 面を適 切な精 度で再 現するためには、信 頼 性の高い計 算 手 法を用いた大 規 模 計 算が必 要であり、大きな課 題となっている 。当 社では計 算 手 法の改 良を進めるとともに、実 製 品への適 用において大 規 模 計 算を加 速するため GPUシステムに適したコード開 発を行っている 。
今 回、TSUBAME 2.0 において並 列 GPGPUアーキテクチャを効 率 的に利 用できる大 規 模シミュレーションコードを実 装することで、飛 躍 的に計 算 性 能とシミュレーション精 度を向 上させることができた 。

1

はじめに

当社では、衛生陶器や浴室、キッチンなどのあらゆる水まわり住宅
設備機器を開発し、国内外で販売している。増大する省エネルギー

本稿では、多相流における数値計算手法の改良およびGPGPUハー
ドウェアへのコード実装における最近の我々の取り組みを示す。

移流方程式の高精度化

や節 水へのニーズに直 面しており、設 計や製 造プロセスにおいて2高い品 質 基 準をクリアする必 要がある。そのため、CAE 設 計ツー
ルとして高精度な流体解析技術を開発し効果的に利用することは以下の移流方程式について考える。

大変重要となっている。
Navier-Stokes 方 程 式に基づいて、気 泡 混 入や薄 膜 流れのような空 気と水が混 在する流れを計 算する数 値 計 算モデルにおいて、移 流 項の計 算 手 法は計 算 安 定 性や精 度において重 要な役 割を果たす。MUSCL 法 [10] は最も代 表 的な移 流 計 算 手 法とされており、

φは移 流される物 理 量を表し、u 、v 、wは三 次 元 流 速 場の各 成 分

市 販の流 体 解 析ソフトウェアに広く採 用されている。MUSCL はを表す。

基 本 的に一 次 式による内 挿関数を用い、また風上 化やリミッタの

一 次 元 の 保 存 型 CIP 法 で あ る CIP-CSL 法 [5] で は、計 算 す る 変

利 用により安 定 化および単調性を保持している。しかし後述する数として少なくとも二 種 類のモーメントを用いる。多 次 元 化す

とおり、MUSCLは拡 散 的であり、流れ場において急 激な値の変 化

る定 式 化は種々考えられているが、メモリ容 量や計 算 量が比 較や不 連 続 性を解くことは難しいと考えられる。さらに例えば液 柱

的 小 さ い 方 法 と し て、図 1 に 示 す よ う に 体 積 積 分 平 均 値 VIA（が壁に垂 直に衝 突した際に、現 象が等 方 的にも関わらずグリッド

Volume Integrated Average）および面積積分平均値 SIA（Surface

依 存 性を示して薄 膜が不 自 然に広がり、妥 当な結 果を再 現できな

Integrated Average）の二種類のモーメントのみを用いる方法がいことがある。

提 案された [6] 。この方 法では、各 軸 方 向に分 割して順に計 算する[1-3]がロバストかつよりグリッド依 存 性

次 元 分 割 法を用いる。各 軸の計 算では、一 次 元の保 存 型 CIP 法をの小さい高 精 度な移 流 計 算 手 法として提 案されている。CIPでは用いる。一つのセルにおいて二 次 関 数もしくは有 理 関 数による補マルチモーメントの概 念を用いており、物 理 量の点 値のみならず間多項式を構成することで、気液二相流れのような不連続性を持

勾配値も更新する。CIPは三次精度を持ち、さらにセミラグランジェ

つ複 雑な流れ場に適 用できる。本 手 法では各 軸 方 向にVIAとそれ法と組み合わせて効 率 的に計 算することが可 能である。しかしなをはさむ二つの SIA で直ちに補 間 多 項 式を求めることができる。

がらCIPは保存性を保証しないため、多相流における小さい水滴やSIAはセミラグランジェ法で更新し、VIAは一般的な有限体積法同気 泡、薄 膜 流れの計 算など、質 量もしくは体 積 保 存 性が必 要 条 件様のフラックスの収支計算で更新するため、保存性が保証される。

となる多くの実例において問題を生じることが分かっている。

このマルチモーメント法はVSIAM3（ Volume/Surface Integrated

他の手 法として、CIP 法

実際の商品設計に適用できる数値シミュレーションを目的とし

Average based Multi-Moment Method）とよばれ、移流方程式のて、保存型 CIP法の一種であるUTI-VSIAM3法 [4] を提案するとともに、

みならず一般的な流体解析ソルバー全体の定式化として提案され二相流シミュレーションコードを実装した。このコードはGPGPUている [6] 。

システムにも移 植され、さらに大 規 模 高 精 細シミュレーションのためにTSUBAME 2.0スーパーコンピュータに最適化された。

02

図 1 VIA および SIA の定 義

図 2 SIA の非 次 元 分 割 計 算

もとのVSIAM3 では、多 次 元 計 算において一 次 元 移 流 方 程 式 計

移 流 計 算 手 法の精 度を検 証する。初 期 値 分 布として、TOTOの算を繰り返す簡易な方法を用いている。別の手法として、次元分割"T"の字の内 部を 1 、外 部を 0 と設 定し、反 時 計 回りの剛 体 回 転 速することなく多次元計算を行う以下の手法を提案する。

度 場とする。図 3 に一 回 転 後の数 値 計 算 結 果を示す。VSIAM3 と図 2 に示すようにSIA-Xの移流計算について考える。移流方程式

UTI-VSIAM3 はMUSCL 法よりも明らかに高精度な結果である。本を、以下のように二つのステップに分割する。

テストではVSIAM3 とUTI-VSIAM3 に明 確な差は現れていないが、ステップ1：UTI-VSIAM3 の次 元 対 称 性は実 際の適 用 問 題において良 好な結果を期待できる。

ステップ2：

ステップ1 では、式（ 2）から面の法線方向、すなわちX 軸方向について一次元保存型 CIP 法を用いてSIA-Xの中間値を計算する。次にステップ2 では、式（3）から面内方向、すなわちY 、Z 軸方向について中心差分を用いてSIA-Xの中間値を更新し、次の時刻ステップのSIA-Xとする。SIA-Y 、SIA-Zについても同様の手順で求められる。
VIAについては、数値積分にて各流入面のフラックスを計算し、有限体積における収支計算によって更新する。数値積分点の値は上記 SIAの計算方法を流用できる。この手法はUTI-VSIAM3（ Unsplit
Time Integration VSIAM3）とよばれ、本質的に各次元方向について完全対称な定式化である [4] 。

図 3 移 流 方 程 式 計 算ベンチマーク

03

マルチGPU 計 算による衛 生 陶 器のための
気 液 二 相 流シミュレーション

気液界面の追跡解法

3

気 液 二 相 流れにおいて、オイラー格 子 法では、各 格 子の流 体 体 積
占 有 率、すなわち流 体 率はもっとも気 液 界 面を認 識するのに適した物 理 量である。VOF 法では、気 液 界 面の追 跡のためにこの流 体率を直 接 計 算する。この方 法では、他の粒 子 法ベースの界 面 追 跡
法と異なり、気 泡や液 滴の分 裂・合 体など界 面 形 状トポロジーの変化に柔軟に対応でき、様々な気液二相流れに汎用的に利用可能である。またVOF 法では流 体 率の保 存 性が保 証されているのも実際の適用問題において大きな利点であり、気液界面追跡手法として最も広く用いられてきた。従 来 型のVOF 法では、流 体 率を更 新する際に用いるフラックス計 算のために界 面の形 状を表す関 数を再構築する必要がある。近年ではフラックスの計算精度を上げる
ため、PLIC 型の界面再構築がよく用いられている [9] 。しかしながら
PLICは一般にアルゴリズムが複雑であり、実装は容易ではない。
一 方 我々は気 液 界 面 追 跡をより高 精 度かつ効 率 的に行うため、
STAA 法（ Surface Tracking by Artificial Anti-diffusion）を提案し

図 4 気 液 界 面 移 流 方 程 式 計 算ベンチマークている [7] 。STAA 法は保 存 型 CIP 法 等と組み合わされ、流 体 率の移流 計 算 後に、数 値 拡 散 補 正のための後 処 理を行う。保 存 型 CIP 法
は高い移流計算精度を持つものの、流体率は気液界面近傍で急激に値が変動するため、他の高次精度解法と同様に数値拡散が生じ、このままでは気 液 界 面の追 跡には利 用できない。そこで界 面の不連 続 性を終 始 一 定に保つため、移 流 計 算の後に、人 工 的な反 数 値

並列 GPGPU システムでの流体解析

4

拡散フラックスを発生し流体率に補正を加える。その際に流体率からレベルセット関 数を生 成し、界 面 法 線ベクトルを計 算する。

Navier-Stokes 流 体 計 算は、フラクショナルステップ法により実補正量はその法線ベクトルに沿って界面方向に強制的に移流され行する。すなわち移 流 方 程 式を解いた後、非 移 流 項である粘 性 項る。これにより、界 面 形 状の忠 実 性を損なうことなく界 面のコンや重 力 項、表 面 張 力 項を計 算し、最 後に圧 力ポアソン方 程 式を解

パクトさを回復することができる。このようにSTAA 法は極めて効

いて次の時間ステップの流速に更新する。圧力ポアソン方程式を率 的かつ実 装が容 易な手 法である。図 4 に STAA 法での前 述の移除き、いずれのステップの計 算も陽 的なベクトル計 算で構 成され流方程式計算結果を示す。界面形状の忠実性において、PLIC 法よているため、並 列 計 算システムにおいてコードを実 行することはり優れた計算結果であることが明らかである。界面近傍の遷移部容易である。今回、GPGPUアクセラレータで構成された分散メモ

分の厚みがPLIC 法よりも厚いが、しきい値の設 定で制 御 可 能であリ型 並 列 計 算システムに適 応させるため、MPI および OpenCLラ

り、設定値よりも広がっていくことはない。STAA 法は実際の問題

イブラリを用いてコードを実 装した。またOpenCLを実 行するデへの適 用においても計 算が破 綻することはなく、極めて実 用 的なバイスは当 面 nVIDIA 社の CUDAアーキテクチャを想 定している方法である。

が、AMDやIntel 等の他のデバイスや将来のシステムにも対応する
上で、OpenCLライブラリを利 用することは柔 軟 性において利 点がある。ただし現 在のCUDAアーキテクチャでは、原 因は不 明であるがOpenCLコードはCUDAコードに対し十～数 十％程 度 計 算
速度が遅いことが分かっており、コンパイル時にマクロによって、OpenCLライブラリをCUDA Driver APIに変換することでネイティ
ブ動 作を実 現している。さらに複 雑なOpenCL 関 数を簡 易 化し流
体 計 算コードをコンパクトかつ見やすく保つため、独 自の簡 易な
ラッパー関 数を定 義し、通 常はラッパー関 数をコールすることで透過的にOpenCL 関数を実行する方法を採用している。

04

図５に、圧力ポアソン方程式解法のCPU 計算とGPU 計算との比較を示している。計 算の収 束 性を考 慮し、全て倍 精 度 演 算を使 用している。GPGPU 化により、コア単体比較で約 8 倍、4コア比較でTSUBAME2.0 への実装

10 倍弱の計算速度向上が得られた。

5

衛 生 陶 器のシミュレーション結 果の一 例を図 6 に示す。複 雑な境界形状に沿った流れが、浮遊物体の相互作用とともに適切に再

近 年、住 宅 水まわり設 備 機 器は、ほぼ例 外なく節 水のニーズが高現されている。浮 遊 物体は、気液界面の識別と同様に各格 子の物まっている。このためそれら機器において少ない水をいかに有効に

体 体 積 占 有 率で表している [8] 。レベルセット関 数を生 成して物 体利用して機能を発揮するかが大きな課題となっている。それに伴い、同 士の面の距 離を算 出し、ばねモデルを利 用して接 触、反 発の計

薄膜流れや微小気泡など計算精度を特に要求する事例が増加して算を行っている。これらの計算も全て並列 GPGPU 化しているため、いる。例えば 2ミリ程度の薄膜流れを再現するためには、膜厚方向高 速に計 算 可 能である。このような従来であればスパコンレベルに数～十数程度のメッシュ分割が必要であり、1mm 未満の格子幅の高性能計算システムが必要であった固気液三相流解析が、GPUが必須となる。このような現象が機器の広い範囲で同時多発的にボードを挿した 1 台の低価格 PCで短時間で手軽に計算できるよう

発 生するような場 合、全 体で数 千 万メッシュ以 上の計 算 規模となになり、これらのPCを複数台利用することで、商品設計時に必要り、従 来のPCベースシステムでは対 応が困 難である。そこで東京な多くのケーストライを短時間で実施可能となった。

工業大学のGPGPUスパコンTSUBAME2.0を用い、大容量高精細メッ
シュでの混相流シミュレーション技術を開発している。
大 規 模 並 列 GPGPU 計 算での問 題として、GPGPU 計 算は通 常のCPU 計算よりも 1プロセスあたりの計算速度が飛躍的に速いため、分割領域間の通信のオーバーヘッドが見かけ上大きくなってしま
い、並 列 化 効 率が上がりにくくなることが挙げられる。一 方 実 製品への適 用では、単なる矩 形 領 域ではない複 雑な流 路をもっており、流路に沿った要素分割では、理想的なスケーラビリティを達成することは難しい。しかし流路を無視して全体を矩形領域に要素分 割する方 法では、流 路のみを要 素 分 割する場 合と比 較して計算
要素数が数倍から十数倍となる場合が多く、実用的とはいえない。
そこで、さまざまな制 約の中で実 用 的な計 算 時 間で実 製 品のシ
ミュレーションを可能にするという観点より、図 7 に示す境界不整

図 5 PCG 法による圧 力 方 程 式 解 法の CPU および
GPU 計 算 速 度 比 較

合な三次元領域分割の計算手法を提案する。本手法は流路格子を
Z 軸→ Y 軸→ X 軸方向に順に、各々の分割領域の要素数がなるべく等しくなるように自動分割する。したがって領域間の隣接部は不整合になり、通信手順が極めて複雑となる。そこで図8のようにオーバーラップ通信と組み合わせ、領域間の通信を同期させずに、あらかじめ決めておいた通信マップにしたがって非同期送信、受信コマ
ンドを一斉に発行する。そしてMPI_Waitallにより、全ての通信が完了するまで待つ。どの順番で実際の通信が行われ、完了するかは
全てシステム側に任せる。
図 9 にTSUBAME2.0 上での実 製 品シミュレーションにおけるポアソン方 程 式のベンチマーク結 果を示す。約 1.5 億の計 算 要 素 数の実製品モデルにおいて、100GPUで 80 倍以上、150GPUで 100 倍
以 上の高 速 化を達 成し、提 案 手 法の妥 当 性を示すことができた。
流体計算においては、大半の計算時間が圧力ポアソン方程式の求解に費やされるため、今回のベンチマーク結果から、高効率な大規模二相流体シミュレーションをTSUBAME2.0 上で実現できたといえる。
図 10 に当社製品のシミュレーション結果を示す。いずれの画像
もTSUBAME2.0 上でPOV-Rayを利 用して出 力した。トイレブース

図 6 衛 生 陶 器シミュレーション結 果

全 体の画 像で示すようにタンク式の衛 生 陶 器であり、計 算 条 件と05

マルチGPU 計 算による衛 生 陶 器のための
気 液 二 相 流シミュレーション

してはタンクに規 定 水 量を貯めて、固 体モデル化したタンク内バ
ルブを規 定 速 度で開くことでタンクの水が落 下し、陶 器 内に流れ込むようにした。計算結果では、タンクからの給水の過程で水路の空気を巻き込みながら無数の気泡を発生させ陶器内を上昇していくようすが捕捉できている。TSUBAME 2.0 上の計算でメッシュを細分化したことで、ミリ単位の気泡、薄膜、飛沫の再現性が飛躍的に向上した。これにより衛生陶器の洗浄性能をより精度よく評価できるようになった。

図 9 PCG 法による圧 力 方 程 式 解 法の複 数 GPU 実 行 性 能

図 7 三 次 元 境 界 不 整 合 領 域 分 割（概 念 図）

図 8 オーバーラップ通 信

06

おわりに

6

衛 生 陶 器 設 計で必 要な精 度および計 算 速 度の気 液 二 相 流 体 解 析を実 現するために、UTI-VSIAM3 およびSTAA 法に基づく計 算 手 法を開 発した。さらに境 界 不 整 合 三 次 元 領 域 分 割により大 規 模 並
列 GPGPU 計算を高効率に実行できるようにした。その結果、実製品の適用においてTSUBAME2.0 上で数ミリオーダーの気泡、飛沫、
薄 膜 流れを捕 捉できるようになった。今 後はさらなる高 速 化のために、圧力方程式解法の収束性を向上させていくとともに、濡れ性を始めとする細部の流体計算モデルを精度向上させ、実製品の設計への適用を加速させていきたい。
参考文献
[1]

T. Yabe, T. Aoki: A Universal Solver for Hyperbolic-Equations
by Cubic-Polynomial Interpolation 1. One -Dimensional
Solver, Comp. Phys. Comm., Vol. 66, pp.219 (1991)[2] 		 T. Yabe, T. Ishikawa, P. Y. Wang, T. Aoki, Y. Kadota, F. Ikeda:
A Univer s al Solver for Hy p erb olic- Equations by CubicPolynomial Interpolation 2. 2-Dimensional and 3-Dimensional
Solvers，Comp. Phys. Comm., Vol. 66, pp233 (1991)[3]		 T. Aoki: Multi-dimensional Advection of CIP (Cubic-Interpolated
Propagation) Scheme, CFD Journal, Vol. 4, pp.279 (1995)[4] 		 A. Ikebata, Y. Muraoka, F. Xiao: Multiphase Fluid Simulations
on a Multiple GPGPU PC Using Unsplit T ime Integration
VSIAM3, Progress in Nuclear Science and Technology, Vol.2,
pp.491 (2011)[5] 		 T. Yabe, R. Tanaka, T Nakamura, F. Xiao: Exactly Conservative
Semi-Lagrangian Scheme (CIP-CSL) in One Dimension, Mon.
Wea. Rev., Vol. 129, pp.332 (2001)[6] 		 F. Xiao, A. Ikebata, T. Hasegawa: Numerical Simulations of
Free-interface Fluids by a Multi-integrated Moment Method,
Computers and Structures, Vol. 83, pp.409 (2005)[7] 		 A . Ikebata, F. Xiao: Development of Conser vative Front
Capturing Scheme and Applications to Multi-f luid
Simulations, Proc. Annual Conference of Japan Society of
Mechanical Engineering, 3, pp.301 (2002) (in Japanese).
[8] 		 F. Xiao: A Computational Model for Suspended Large Rigid
Body in 3D Unsteady Viscous Flow, J. Comp. Phys., Vol. 155,
pp.348 (1999)[9] 		 W.J. Rider, D.B. Kothe: Reconstructing Volume Tracking, J.
Comp. Phys., Vol. 141, pp.112 (1998)[10] B.V. Leer: Towards the Ultimate Conser vative Dif ference
図 10 TSUBAME2.0 での衛 生 陶 器シミュレーション能Scheme Part V : A Second Order Sequel to Godunov's Method, J.
Comp. Phys., Vol. 32, pp.101 (1979)07

Liグラファイト層間化合物におけるLi拡散に対するハイブリッド量子古典シミュレーション
大庭 伸子 *

尾形 修司 ***（株）豊田中央研究所** 名古屋工業大学

河野 貴久 ***

旭 良司 *** 東京大学・物性研

Li-グラファイト層 間 化 合 物における Li 拡 散 過 程を、ハイブリッド量 子 古 典シミュレーション法により調べた 。本シミュ
レーションでは、挿 入 Li とその周 辺グラファイトの炭 素 原 子を電 子 状 態 計 算（= 量 子）領 域とし、他の炭 素 原 子は経 験 的
原 子 間 相 互 作 用モデルを用いる古 典 領 域とする 。量 子 領 域は Li の移 動に伴って再 選 択する 。一 定 温 度でグラファイトの平 均 層 間 距 離を様々に設 定した一 連のシミュレーションにより、ゼロ圧 力あるいは少し伸 張された場 合には、層 間が圧縮された場 合に比べ拡 散が速くなり、ホッピング的な拡 散に加え弾 道 的に移 動することもあることを確 認した 。さらに、
グラファイト層 間の中 央 付 近に Li が位 置する瞬 間に、Li が速く動くことも分かった 。そこで、電 場を用いて層 間におけるLi位置を制御することを検討した結果、適切な電場強度と周波数によってLiの面内熱拡散が促進される可能性を見出した。

はじめに

1

（iii）古典 MDで必須となるLiとC 原子間のポテンシャルモデルを作成する必要がない。（後で述べるが、）Li 拡散の活性化エネ
ルギーはかなり小さく、高精度のLi-C 間ポテンシャルモデル
を作成することは容易ではない。
さらに、層間におけるLi位置と拡散挙動についての詳細を解析し、グラファイトに代表される層状物質はその層間にイオンなど他の外的遥動、たとえば、電場を用いて層間におけるLi 位置を制御する

物質を取り込み、層間化合物を形成する。特に、層間にリチウム（Li）ことによる拡散の促進・抑制可能性について検討した。

が挿 入されたLi-グラファイト層 間 化 合 物（ Li-GIC）はリチウムイ
オン電 池の負 極にも実用化されている材料である。リチウムイオ
ン電池における基本原理は、電解質相を介して正負極間でLiが移動することであり、電池性能の観点からは、グラファイト中 Liの拡
散速度は電池から取り出すことができる出力に直接対応する。Li-

計算方法

GICではLi 挿入に伴い層間距離がおよそ 10 ％伸張することが確認

2

されており、グラファイトの構造変化により生成する応力場は、Li
拡散に影響を及ぼすと考えられる。このような系をシミュレーショ
ンで取り扱う際には、長距離まで作用する応力場を考慮することが重要となる。

境界部分のモデル化

ハイブリッド量 子 古 典シミュレーションにおいて、量 子 力学計 算

そこで本研究では、ハイブリッド量子古典シミュレーション法をが行われる領 域と古 典 MD 計 算のみで取り扱われる領 域の境 界で用いて、まず、グラファイトの平均層間距離とLi 拡散の関係についは、原 子 間の結 合が切 断され境 界 表 面にあたる原 子にダングリンて調べた。ハイブリッド量子古典シミュレーション法では、電子状

グボンドが形 成される。本 来、このダングリングボンドは存 在し態の変化が重要となる領域には並列性能が良い実空間表現での密ないはずのものであるので、領域内部の電子状態や結合距離等に度汎関数法（DFT）などの量子力学計算法を適用し、その他には経対するその影響を取り除くことが必要となる。しかし、ダングリン

験的相互作用モデルに基づく古典分子動力学（MD）法を適用する。

グボンドを単に連 結 原 子で終 端するだけでは、表 面 原 子の緩 和にこれにより、着目する化学反応に対し、周囲の歪緩和を考慮した系

伴う表面再構成の影響が全系の原子にも及び、本来保持すべき構の形状および大きさを取り入れた現実的な大規模原子系のシミュ

造から大きなずれを生じる場 合がある。そこで本 研 究では、量 子

レーションがTSUBAME2.0 のような多 数コアのスパコンで可 能と領域と古典領域の境界部分のモデル化に関して、Ogata[1] によってなる。特に、Li-GICに本手法を適用する利点は以下の通りである。

考案されたbuffered cluster 法を用いた。Buffered cluster 法では、

（i）ハイブリッド量子古典シミュレーションでは、リーズナブルな領域の境界表面で切断された結合を終端するために付加的な原子計算コストで、まわりの多数のC 原子の影響を考慮したLi 拡散（buffer atom）を置く。その位置は古典 MD 計算によってのみ決定過程を調べることができる。古典 MDでは、LiとC 原子の間のされ、量子力学計算を行う際に緩和されない。さまざまな対象系に電荷移動（化学反応）を高精度で取り扱うことができないし、対して表面再構成によるひずみが生じないことが実証されており、第一原理 MDで大規模計算を行うには計算コストがかかり過計算コストも低い、高精度で汎用的な結合モデルであるといえる。

ぎる。
（ii）ハイブリッド法では、古典系の層間ポテンシャルとしてグラファ

08

２．１

２．２

量子力学計算イト層間に働く弱い分子間力（分散力）を容易に考慮すること

量 子 力 学 計 算には、計 算の大 規 模 化に対 応しやすく、さらに境 界ができる。通常のDFTでは分散力を再現することができない。

条 件を柔 軟に設 定できることから、実 空 間 密 度 汎 関 数 法を用いた。計 算プログラムの詳 細は文 献 [ 2 、3 、4 ] を参 照されたい。電子と原 子 核との相 互 作用には、Troullier-Martins 型ノルム保存擬
ポテンシャル [5] を用いて価 電 子のみを考 慮した。交 換 相 関エネ
[6]

２．４

計算モデル

Li 挿入の初期過程を想定し、3072 個のC 原子から成るAB 型構造のルギー項には、部 分 内 殻 補 正 を取り入れ、一 般 化 密 度 勾 配 近 似

グラファイトにLiを 1 個挿入した（図 1）。層間方向をz 、層間に垂（ Generalized Gradient Approximation:GGA）を適用し、Perdew,

直な面をx-y 面とし、温度 T=0 Kにおけるハイブリッド計算より、グ[7]Burke 及びErnzerhofらによって導入された表式（PBE）
を用いた。

ラファイトの平 衡 距 離としてL x =34.44 Å, L y =29.83 Å, L z =26.97 Åを

Kohn-Sham 軌道とHartreeポテンシャルを、実空間直交座標系の得た。グラファイトの面 内 方 向の熱 膨 張は小さいので、有 限 温 度

格子点上の値で離散化 [8] し、二階微分は四次の有限差分法によりにおいても面 内 方 向（ L x , L y ）の平 衡 距 離は変 化しないとする。層近似した。

間（ z ）方 向に関しては、古 典 MD 計 算から得られた層 間 方 向の線-5

空 間 格 子を計 算ノードごとの小 領 域に分 割した並 列 化を行い、膨 張 係 数 6.7 × 10 を乗じ、T=423 K における平 衡 距 離として L z =計算の高速化のため格子点上のデータを計算ノード毎に保存した。

27.74 Åとする。これを基準として、
平均層間距離 L z を 0.6,2.9 %伸張、波動関数の空間分布における長波長成分の収束性加速にはマルチ

3.9,1.7 %圧縮した場合のLi 拡散過程を調べる。

グリッド法 [9] を使用している。空間分割のメッシュサイズとして、
h=0.45 a.u.（1 a.u. ～ 0.5292 Å）を用いた。これは、波動関数の平面波
2展開におけるカットオフエネルギーに換算すると（π /h） ～ 49 Ry（ 1 Ry = 13.6 eV）に相当する。なお、原子核近傍では擬ポテンシャ
ルの為にさらに細密なメッシュ（ h/3 ）を用いた。
外 部 電 場 Eが印 加された場 合、密 度 汎 関 数 理 論におけるKohnSham 方程式は以下となる。

ただし、、 はそれぞれ、固有軌道および固有エネルギーであり、

v eff は外 部 電 場が無い時の有 効ポテンシャルである。外 部 電 場 項によりＩ 番目の原子に働く力はF I = Z I E（ Z I は価電子数）となる。
２．３

古典分子動力学計算

古 典 MD 計 算では、周 期 境 界 条 件を課し、ニュートンの運 動 方 程式の数 値 積 分に速 度ベルレの方 法を用い、時 間 刻みを 0.97 fsとした。また、C 原子間ポテンシャルとしてブレンナ型ポテンシャル [10]図 1 y 方 向から見た計 算モデル図 。
大、中、小 球はそれぞれ Li 、量 子 力 学 計 算の対 象とした量子領域（QM region）のC 、それ以外のC 原子（古 典 領 域 :CL region ）を表す 。

を適 用した。本 計 算では、カーボン系 材 料の結 合エネルギー等をできるだけ再 現するように最 適 化されたブレンナ型ポテンシャルのパラメータ（文献 [10] : TABLE1 を参照）を用いている。ただし、平衡距離については、量子力学計算による平衡原子間距離と等しくなるよう 1.33116 Åとした。ハイブリッドシミュレーションでは、
量子 - 古典の境界において人工的な応力を発生させないために、そのような微調整が必要となる。
ブレンナ型ポテンシャルでは、原 子 間 相 互 作 用のカットオフ距
離が 2.0 Åに設 定されており、グラファイトに対しては面 内の最 近接原子間の相互作用のみが考慮され、面間（層間）に働くファンデ
ルワールス相 互 作 用が取り込まれない。そこで、レナードジョー
ンズ型を基 本として設 定した原 子 間ポテンシャルV vdw [11] を組み込
み、層間相互作用を補正する。

09

Liグラファイト層 間 化 合 物におけるLi 拡 散に対するハイブリッド量 子 古 典シミュレーション

層間距離と Li 拡散

3

距 離 L z 依 存 性について、以 下の特 徴を知ることができる。層 間を0.6 % 、2.9 %伸張した場合の平均二乗変位は、圧縮した場合よりも大きくLiの拡散が速い。ただし、2.9 %伸張した場合には、長時間の拡散はやや遅くなる傾向がある。これは、Liが上下どちらかのグラ
ファイト層にトラップされるためであり、その詳細を 4 章に述べる。

３．１

グラファイトの構造変化とLi 拡散圧縮の場合に長時間の拡散が低下しているのは、先に述べたcage
効果によるものと考えられる。

10 ピコ秒（ ps）間のLiの軌跡を図 2 に示す。開始（最終）のLi 位置

を白抜き矢印（ 大きい赤球 ）で表している。Liは、（ x 、y ）座標が同じ値で位 置するC 原 子を避けるような経 路でホッピング的に拡 散している。グラファイトの層 間 距 離が長くなるとLiは広い範 囲にわたってよく拡散する。さらに、ホッピング的な拡散に加えて、瞬間的に開いた上下のC 原子間を弾道的に進む動き（ 図 2 黒矢印）がみられた。弾 道 的な動きが起こる際には、層 間 距 離が 3.85 Å 以 上
広がっていることが必 要である。一 方、層 間 距 離が短くなるとLiの動く範 囲は限 定 的となっており、これはグラファイトによりLiが閉じ込められたcage 効果によるものと考えられる。

図 3 平 均 層 間 距 離 L z を伸 張（0.6 % , 2.9 %）および
圧 縮（-1.7 % , -3.9 %）した際の、Li の平 均 二 乗 変 位の時間変化。

図 2 z 方 向から見た 10 ピコ秒 間の Li の軌 跡 。
開 始の Li 位 置を白 抜き矢 印、最 終の Li 位 置を大きい
赤 球で表している 。緑 球は QM 領 域における C 原 子 。
黒 矢 印は x , y 座 標が同じである上 下のＣ 原 子の間をLi が通り抜けた場 所を示す 。
Liの平 均 二 乗 変 位を図 3 に示す。本シミュレーションでは系にLiイオンが 1 個しか存在していないため、より統計性の良いデータを得るため、一つの長いシミュレーションから 0.01 ps 毎に時 間の原点をとった時系列データを切り出して平均二乗変位を計算した。
したがって、図 3 に示したデータは拡 散 係 数を見 積もることが可能なほどには十 分な直 線 部が得られていないが、拡 散の平 均 層 間

10

図 4 グラファイト中 Li の基 底 状 態（ground state ）と遷 移 状 態（transition state ）の位 置 。

３．２

活性化エネルギー

Li 拡 散の詳 細を調べるため、T=0 Kにおけるホッピング的な拡 散の
活 性 化エネルギーを計 算した。図 4 は拡 散の基 底 状 態（ groundstate）と遷移状態（transition state）におけるLiイオンの位置を示している。活性化エネルギーは、遷移状態と基底状態のエネルギー差から見積もった。
図 5 に各 平 均 層 間 距 離 L z における活 性 化エネルギーを示す。全エネルギーから見 積もられた値は、量 子（ QM）領 域と古 典（ CL）
領 域の両 方の寄 与を含んでいる。活 性 化エネルギーはいずれも0.07 eV 以下であり、とても小さい。従って、グラファイト中のLiイ

オンは熱エネルギーにより容 易に拡 散する。また、層 間 距 離が短くなると活 性 化エネルギーは増 加する。これは、圧 縮の際に拡 散

表 1 VASPコードを用いた基 底 状 態と遷 移 状 態のエネル
400K
400K
400K
400K
ギーの計算結果。E ele , E vib , S vib , G = E ele + E vib
400K+pV-TS vib は、それぞれ T=0 K における電 子 状 態のエネルギー、T=400 K におけるゼロ点振動を含む振動のエネルギー、振動のエントロピー、ギブスの自由エ
ネルギーである。電子状態のエネルギーは原子の標準状態を基準として算出された。

が遅くなることと対 応する。Li 周りの局 所 的な効 果を含んでいると考えられるQM 領 域のエネルギーから見 積もられた活 性 化エネ
ルギーのL z に対する変化は小さい。一方、CL 領域（Liから離れた領域）の寄与は、特に層間距離が圧縮された場合に大きくなり、Liから離れたC 原 子の変 位が活 性 化エネルギーに影 響している。活 性

外部電場による Li 拡散の促進

化エネルギーに関する計算精度確認のため、電子状態計算ソフト

4

VASP [12] を用いても計算を行った。AB 型積層構造のグラファイト
に対して 3 × 3 × 2 のスーパーセルを用意し（ 72 個のC 原子）、Li 原子を一つ挿入した。基底状態（GS）と遷移状態（TS）における計算

４．１

電場 E = 0 の場合のLi 位置と拡散

結 果を表 1 に示す。電 子のエネルギーから見 積もられた活 性 化エ

はじめに、外 部 電 場 E=0 の場 合における、層 間でのLi 位 置と拡 散

ネルギーは 0.09 eVと、ハイブリッド法で見 積もられたものとオー

挙動について考察する。Liが良く拡散する状態で解析を行うため、ダーが同じであり、両者の計算モデルの違いを考慮すると、良く一層間距離が基準から 2.9 % 伸張した場合（ L z = 28.55 Å）を解析する。

致しているといえる。また、MedeA のphononパッケージ [13] を用

図 6 に、Liイオンと上下グラファイト層の間の距離 d u , d l の時間発展

いて、有限温度における自由エネルギーを計 算し、400 K におけるを示す。より詳細にはd u（ d l ）は、上（下）グラファイト層の量子力

振 動エントロピーの寄 与として TΔS = -0.14 eV を得た。したがっ学計算の対象としたC 原子群の重心とLiとの距離である。Liのt=0て、ギブスの自 由エネルギーから見 積もられた活 性 化エネルギーにおける位 置を基 準としたx,y,z 方 向の変 位の時 間 変 化も図 6 の下は 0.04 eVとなり、とても小さいことに変わりない。

d l はd u より短くなっており、パネルに示す。時間 t = 0 から 1 psでは、Liが下グラファイト層にトラップされている状態であることを意味する。さらに、同じ時間帯のLi 変位から、Liはほとんど動いていな
いことも分かる。一方、d l とd u がほとんど同じ値となっている時間
t=1 から 1.2 psの領域では、Liの面内（ x , y 方向）の変位が大きくLi
d u とLiの面内変位との関係は、が良く拡散している。このようなd l 、
t=1.6 と 4.2 ps 付近でも見られる。したがって、上下グラファイト層の真ん中あたりにLiが位置しているときには、拡散の抵抗力が弱くなり、Liがよく拡散するといえる。

図 5 温 度 T=0 K における活 性 化エネルギーの平 均 層 間
距離 L z 依存性。黒丸は全エネルギーから見積もられた値であり、白丸はQM 領域のエネルギーのみから見 積もられた値 。エラーバーは実 空 間 密 度 汎 関
数 法における収 束 誤 差を示している 。

11

Liグラファイト層 間 化 合 物におけるLi 拡 散に対するハイブリッド量 子 古 典シミュレーション

ギーと関 係があると考えられる。本 計 算で設 定した層 間 距 離はL z=28.55 Åであるので、拡散の活性化エネルギーは図 5 より 0.04 eV

程 度と見 積もられる。また、遷 移 状 態と基 底 状 態の間の層 間（ z）
方 向の距 離はおよそ 0.2 Åである。したがって、熱エネルギーによる助け無しでこの活 性 化 障 壁を越えるのに必 要な電 場の閾 値は、0.04/0.2 = 0.2 V/Åと見 積もられる。この閾 値を超える電 場 強 度になると、Liはグラファイト層どちらかにトラップされる傾向が強くなり、拡散が遅くなる。閾値以下では、それほどの影響は無い。

（a ） 電 場 E=0.4 V/Å, 周 波 数 f=0.4 THz

図 6 電 場 E = 0 の時の、外部電場強度 E（上パネル）、z
Liイオンと上下グラファイト層間の距離 d u 、d（中l|パネル）、Li 変 位の絶 対 値 |(R Li(t) － R Li(0) ) x,y,z（下パネル）の時 間 変 化 。

４．２

電場印加時における拡散

前 節の結 果より、たとえば、外 部 電 場を層 間に印 加することにより、グラファイト層の間で Li の位 置（存 在 確 率 ）を制 御することで、Li 拡散が促進・抑制できることが予想される。そこで、強度が0.14,0.43 V/Å 、周波数 f=0.4, 0.8,1.7 THz の電場をグラファイトの層間方向に印加しハイブリッド量子古典シミュレーションを行った。

（b ） 電 場 E=0.4 V/Å, 周 波 数 f=1.7 THz

CL 領 域の炭 素 原 子はイオン化されておらず外 部 電 場に対する応答はほとんど無いと考えられるので、QM 領 域の原 子 群に対して、電場を印加した。

図 7 さまざまな強 度 E amp z および周 波 数 f の外 部 電 場を印 加した際の、Li の平 均 二 乗 変 位の時 間 変 化 。

12z
図 8 電 場 E = 0.4 V/Å 印 加 時の、外 部 電 場 強 度 E（上パネル）、Liイオンと上 下グラファイト層 間の距d（中パネル）l、Li 変位の絶対値 |(RLi(t) − RLi(0)) x.y,z|
離d u 、（下パネル）の時間変化。中パネル横のグラフは、
5ピコ秒間の統計をとったd u 、d l の存在確率分布を示す 。

温 度 T=423 KにおけるLiの平 均 二 乗 変 位を図 7 に示す。電 場の周

電 場 強 度を 0.43 V/Åに固 定し周 波 数を変 化させた場 合には、周波 数を 0.4 THzに固 定したとき、電 場 強 度が大きくなるにつれLiの

0.4 THzと 0.8 THzでは平均二乗変位にあまり違いが無いが、波数が、

平 均 二 乗 変 位が小さくなっている。これは層 間 方 向に働く電 場

1.7 THz としたときには Li 拡 散が速くなった。周 波 数 0.4 およびからLiが感じる力のために、Liは層 間の上 下の安 定サイトのどち

1.7 THz の時の電場、Liと上下層間距離 d u 、d l 、変位の絶対値 （R| Li(t)らかに誘 導されていることが原 因であり、拡 散の活 性 化エネル

− RLi(0) ）x.y,z |の時間変化を図 8 に示す。図 8（ a）より周波数が 0.4 THz

では、上（下）向き電場が印加されるとLiは上（下）層との距離が短れ、その発生と検出が難しく課題も多い。しかし、高出力テラヘル

くなり、上（下）層にトラップされる傾向がある。故に、先に述べたツ光源であるジャイロトロンの開発 [15] なども進んでおり、本シミュ様に、そのような比較的低い周波数ではLi は上下の安定サイトのレーションで設定した電場効果を検証する実験も将来的には可能どちらかに留まることになり、ほとんど拡散しなくなる。一方図 8となると考えられる。

（ b)より、周 波 数が 1.7 THz では、外 部 電 場により生じる層 間 方 向の力を交互に受け層間におけるLi の位置は無秩序に変化している。
言い換えると、電 場の方 向 転 換があまりに速いためLi はその動きについていくことができず、もう一 方の安 定サイトに到 達する前に電 場から働く力によりその方 向を変えてしまっているようであ

5

おわりに

る。その結果、上下層の中央付近にLiが存在する確率が高くなり、Li 拡散が促進される。実際、図 8（b）の中パネル横に示したとおり、
5ピコ秒間の統計をとったd u 、d l の存在確率分布はどちらも層間のハイブリッド量 子 古 典シミュレーション法により、Li-GICのグラ

中央付近にピークを持っている。

ファイトの層 間 距 離とLi 拡 散との関 係を調べた。層 間が圧縮され強 度が 0.43 V/Åである電 場を静 的に印 加したときの電 子 密 度 分

ると拡 散は遅くなり、伸 張されると速くなる。層 間が伸 張された布を図 9 に示す。電 場 強 度が大きい場 合には、電 子 雲が原 子 核 位

場 合、ホッピング的な拡 散に加えて、炭 素 原 子の間を弾 道 的に拡

置から大きくずれることがあり得る。本 研 究で検 討した電 場 強 度散する様 子が見られた。圧 縮された場 合はホッピングによる拡 散0.43 V/Å はかなり強い電 場ではあるが、電 子 分 布は量 子 領 域の炭のみが見られ、cage 効 果により長 時 間 限られた範 囲にとどまる。

素原子の周りに留まっている。

このような層 間 距 離とLi 拡 散の関 係が、ステージ混 在 状 態における拡散係数の低下の一因となっていると考えられる。
グラファイトの層 間 方 向に電 場を印 加し、Liの熱 拡 散が促 進される条件として、電場強度 0.43 V/Å , 周波数 1.7 THzを見出した。これは、上 下 層 間の真ん中あたりに位 置するLiの存 在 確 率が高くなるためであると考えられる。
本研究では、量子力学計算として並列化に適した実空間密度汎関数法による電子状態計算コード（RGDFT）を適用しているが、系3に含まれる原 子 数 N に対して計 算 量がN で増 大する為に、計 算コ
ストの観点から多くの原子を取り扱うことが難しい。そこで、量子
領 域を分 割 領 域の重ね合わせで表 現することによりオーダー N 化を可能にしたDivide-and-Conquer（ DC）型 -RGDFTコードの開発を別途進めている [16] 。開発したDC-RGDFTコードの計算時間に関する TSUBAME2.0 上でのベンチマークテストの結 果を図 10 に示す。計 算は原 子 数 N=48,384,1296,3072,6000 の Alクラスタをそれ
図 9 電場強度 0.43 V/Åを印加した際の（010）面の価電子
密 度 分 布 。黒、灰 色 球はそれぞれ Li ,C 原 子を表す 。
赤くなるほど電 子 密 度が高い 。

ぞれ、1,8,27,64,125 個の領域に分割して行い、電子状態計算の繰り返し回数（SCF）を 25 回とした。計算機 1コアあたりに分配される
原 子 数が 4 原 子 程 度となるよう計 算コア数を設 定したところ、従3来の計算手法では全計算時間は原子数 N に対してN で増大するが、-1

テラヘルツ領 域は波 数に換 算すると 3 THz ≈ 100 cm であるの

DC-RGDFTでは原 子 数 N に比 例することを確 認した。今 後は、古でちょうど格 子 振 動 数 程 度に相 当する。しかし、Li-GICでは、1 ∼典 Ｍ Ｄ 計 算プログラムとのハイブリッド化を行い、グラファイト

2 THz 領 域の格 子 振 動 数ピークが観 測されているとの報 告 例は無中に複 数 Liが存 在する場 合に適 用し、Li 拡 散におけるLi-Li 相 互 作

い。一 方、グラファイトの二つの層が平 行に動く振 動モードE 2g(1)用の影響などを調べる。

として 42 cm -1 が観 測されている [14] 。これは周 波 数に換 算すると、1.26 THz に相 当する。もし、外 部 電 場の周 波 数をこの振 動モードに同調するようにすると、一方の層に張り付いているLiを剥がして上 下層間の真ん中辺りのLiの存 在 確 率を大きくし、Li 拡 散を増 長する可能性が考えられる。しかし、今回のシミュレーションではそのような振 動モードとの共 鳴は観 測されなかった。テラヘルツ波
は光と電波の境界領域の周波数帯であり、長年未開拓電磁波とさ

13

Liグラファイト層 間 化 合 物におけるLi 拡 散に対するハイブリッド量 子 古 典シミュレーション

based on K. Parlinski, phonon 3.11, 2002.
[14] R. Kostic, M. Miric, T. Radic, M. Radovic, R. Gajic, and Z. V.
Popovic: Acta Phys. Pol. A 116, 718 (2009).
[15] T. Saito, T. Notake, Y. Tatematsu, A. Fujii, S. Ogasawara, L.
Agusu, T. Idehara, S. Kubo, T. Shimozuma, K. Tanaka, M.
Nishiura, K. Kawahata, and V. N. Manuilov: J. Phys. Conf. Ser.
227, 012013 (2010).
[16] N. Ohba, S. Ogata, T. Kouno, T. Tamura and R. Kobayashi,
Comput. Phys. Commun. 183, 1664 (2012).

図 10 全計算時間と原子数 Nの関係。

謝 辞
本計算の一部は、文部科学省の先端研究施設共用促進事業の補助をいただき、東京工業大学のスパコンTSUBAME2.0 を用いて行った。
また、科学研究費補助金・基盤研究（B）課題番号 233100074「電場中でのナノ液固界面の化学反応に関するハイブリッド量子古典計算法の開発」から支援を頂いた。記して謝意を表す。
参考文献
[1]

S. Ogata, Phys. Rev. B 72, 045348 (2005).

[2]

S. Ogata, E. Lidorikis, F. Shimojo, A. Nakano, P. Vashishta, and
R. K. Kalia, Comput. Phys. Commun. 138, 143 (2001).

[3]

S. Ogata, F. Shimojo, R. K. Kalia, A. Nakano, and P. Vashishta,
Comput. Phys. Commun. 149, 30 (2002).

[4]

F. Shimojo, T. J. Campbell, R. K. Kalia, A. Nakano, S. Ogata, P.
Vashishta, and K. Tsuruta, Future Gen. Comp. Sys. 17, 279 (2000).

[5]N. Troullier and J. L. Martins, Phy. Rev. B 43, 1993 (1991).

[6]

S. G. Louie, S. Froyen, and M. L. Cohen, Phy. Rev. B 26, 1738
(1982).

[7] J. P. Perdew, K. Burke, and M. Ernzerhof , Phy. Rev. Lett. 77,
3865 (1996).
[8]

J. R. Chelikowsky, N. Troullier, and Y. Saad, Phys. Rev. Lett. 72,
1240 (1994).

[9]

A. Brandt, Math. Comput. 31, 333 (1977).

[10] D.W. Brenner, Phys. Rev. B 42, 9458 (1990); 46, 1948(E) (1992).
[11] N. Ohba, S. Ogata, T. Tamura, S. Yamakawa and R. Asahi,
CMES-Comp. Model. Eng. 75, 247 (2011).
[12] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
[13] K. Parlinski, Medea-phonon, materials design, inc., 2003,

14

建築分野の数値流体解析における大規模計算
ファム バン フック*野津 剛 *

菊池 浩利 * 日比 一喜 *** 清水建設株式会社 技術研究所・総合解析技術センター** 株式会社 数値フローデザイン

実 応 用 問 題や計 算 手 法に応じて計 算 機の性 能を予 測することは重 要である 。
本 稿では、建 築 分 野の流 体 問 題を対 象に、計 算ノード本 体の性 能および計 算 機のネットワークと計 算 規 模による性 能 特
性を調べ、大 規 模モデルの大 規 模 並 列 計 算によりその性 能 向 上が可 能であることを確 認した 。また、高 解 像 度の利 用による計 算 精 度の向 上 化や解 析モデルの精 緻 化、及び 3 次 元 広 域 津 波 解 析などのような新たな解 析 手 法の実 用 化も確 立されつつ、建 築 分 野での大 規 模 計 算の重 要 性を明らかにした 。

はじめに

1

数値流体問題における計算性能の動向

210 年 間でスパコンの計 算 性 能を 1000 倍に向 上させること [1] が、

数 値 流 体 問 題における計 算 機の性 能 特 性を把 握することは実 応

複数計算ノードと高速なネットワークおよびGPUなどの最新デバ

用 促 進において重 要である。本 章では、計 算ノード本 体の性 能おイスの活 用により着々と実 現されている。一 方で、産 業などのもよび計 算 機のネットワークと計 算 規 模による性 能を調べ、その計のづくりの現場において計 算 機の日 常 利 用としてどのように変 化

算性能の動向を明らかにする。なお、本検討では数値流体解析オーしているのか、またはどのように変化していくべきか、あまり議論

プンソースOpenFOAM[2] を用いた。

されていないのが現状である。
建築分野では、周辺環境の配慮や安全・安心かつ快適な建築環

２．１ 単体計算ノード内の性能特性 [3]

境の提 供が重 要な課 題の一つである。この課 題を解 決するために

図 1 には非 構 造 格 子の 1 千 万メッシュ規 模の風 環 境 問 題を対 象

は、たとえば風荷重、風環境、温熱環境、光環境、さらには津波荷重にして、年 代 CPUによる計 算 向 上を示す。図 中 縦 軸は 2006 年 時などの様々な評 価が必要となってくる。経験 的な評価または実 験

点の CPUコアまたは計 算ノードの計 算 速 度を 1 とした加 速 率をによる評価は限界となりつつあり、数値シミュレーションによる評示す。利 用した CPU は 2006 年の TSUBAME1.0 で採 用した CPU

価は有効な手法の一つとして見られている。また、近年の建築分野

AMD Opteron と 2010 年 の TSUBAME2.0 で 採 用 し た IntelCPU

ではBIM（BuildingInformation Modeling）などの大規模 3 次元デー

Westmere-EP 、及 び 市 販 計 算 機 の 2008 年 の Intel CPU Xeon とタの活用や多様なニーズの対応により、数値シミュレーションでは2012 年のIntel CPUSandyBridgeである。2008 年 以 降では、CPU解 析モデル自 体が複 雑となり、計 算 量も膨 大になる。計 算 機の実のクロック周波数が上がっていないが，浮動素数演算能力やキャッ行性能の把握及びその能力を最大化することは必要かつ緊急な課

シュ効率性能などの向上によりCPUコアにおける加速率が少しず題となってくる。

つ伸びている。ただし、その加速率が減速しているため、今後の単

本稿では、TSUBAMEスパコンの活用を通じて、ユーザーという立

体 CPUコアを用いたシリアル問題ではその計算性能を向上しにく場で建築分野の流体問題における計算機の実行能力の見極めを試いとみられる。一 方、ノード内で複 数 CPUコアを用いた並 列 処 理

みるために、各種の建築環境解析を実施するとともに、建築分野でではコア数の増加に伴い、その加速率は 2 年間で 1.5 倍～ 2 倍で順の大規模計算の重要性を示す。

調に伸びているため、ノード内の並 列 処 理は今 後もますます重 要である。ただし、10 年間で 32 倍以上を性能向上させることは現状で難しいと分かる。
また、図 2 には，数 値 流 体 計 算で最も計 算 時 間を要する部 分であるPoisson 方 程 式のそれぞれの計 算 手 法に対してGPUによる加速率を示す [4] 。6 倍以上を得る手法もあるが、GPU 等の新デバイスの活 用と、それに適 切する新たな計 算 手 法の開 発とチューニング
は計算性能の向上において今後も重要な要素となってくる。

15

建 築 分 野の数 値 流 体 解 析における大規模計算

図 1 年 代 CPU による計 算 性 能 向 上

図 3 計 算 規 模による加 速 率の変 化

大規模数値流体計算の可能性

3

前 章では、計 算 機ハードウェアに着 目し、数 値 流 体 解 析の計 算 性
能を議 論しており、その能 力を最 大 化するためには大 規 模モデル
での大 規 模 並 列 計 算が重 要であることを示した。本 章では、数 値
流体計算の応用面に着目して、大規模計算による可能性を示すことが目的である。
図 2 GPU による計 算 性 能の向 上

数値流体解析を大規模に計算することは、計算格子の高解像度による高精度化、解析モデルの精緻化及び解析領域の広域化などを行い、より実 物と近い形 状や実 現 象と近いスケールを再 現する

２．２ 大規模並列計算の性能特性 [3]ことである。ここでは、風 洞 気 流を高 精 度に再 現するための高 解
像度の計算格子利用、複雑な 3 次元建物形状を再現するために小近 年の計 算 機ではネットワークの性 能が良くなっており、複 数

さいメッシュ幅の利用によるモデルの作成と、従来では困難であっ計 算ノードを用いた大 規 模 並 列 処 理による計 算 性 能の向 上がたVOF 法による広域 3 次元津波解析について紹介する。

期 待されている。図 3 には 1 千 万 格 子 規 模 及び 4 億 格 子 規 模の風 環 境 問 題 を 対 象 に し て、2006 年 の TSUBAME1.0 で 採 用 し て

３．１ 風洞実験の再現 [5]いたInfinibandSDR（ 2 本 SDR 、バンド幅：20GBps）と 2010 年の

図 4 はLESによる大型境界層風洞の助走路及び測定部分を忠実に

TSUBAME2.0 で採用していたInfinibandQDR（ 2 本 QDR 、バンド幅：再現した一例である。計算では 100 万格子（A）と 1000 万格子（B）

80GBps）における単ノードに対する加速率を示す。同じ 1 千万格

及び 1 億 格 子（ C）モデルを用いた。図 5 は風 洞 内で発 達してきた子 規 模においても 2006 年から 2010 年までの 4 年 間でその加 速 率

境界層気流の風速変動のパワースペクトル密度の計算値と実験結は 6 倍以上に伸びており、4 億格子の大規模モデルではその性能は

果との比較である。図内では格子幅と格子規模が記載されている。

さらに向上する。前節のノード内の年間加速率の結果と合わせて、計算結果は実験結果と比べると無次元周波数 nL x /U H =1 以上の領

大 規 模 並 列 処 理 及び大 規 模 計 算モデルの利 用により 10 年 間で数域において格 子 Aのパワースペクトル密 度は急 激に減 少している

値流体解析の計算性能を 1000 倍に向上できると分かる。

が、より高解像度の格子 Bまたは格子 Cの利用によりその高周波数成 分が実 験に漸 近している。ちなみに、実 験 結 果を完 全に再 現す
るためには結果の外挿により約 1 千億計算格子が必要となる。

16

図 4 風 洞 実 験を再 現した計 算モデル

図 5 風 速のパワースペクトル密 度

図 6 建 物に作 用する外 圧の平 均 値（左）と負 側ピーク値（右）図 6 は実 験および解 析から得られた建 物の耐 風 設 計に必 要な平均 外圧係数および負側のピーク値の比較一例である。3 千 万 格 子
規模および高解像度の 1.2 億格子規模のいずれの計算結果も実験
結 果と類 似した結 果となっているが、入隅部付近の側面での局 部
風圧の再現については高解像度の格子を用いる必要がある。
また、図 7 には建物近傍の平均流線を示す。3 千万格子では、地上付近の流れは低層屋根面の前縁角で剥離して、高層部からの吹き降しにより低層屋根の入隅部近傍で馬締形渦が生じており、その渦が巻き込みながら側面に流れていく。一方、1.2 億格子では、入隅部近傍で生じた馬締形渦は側面に移流し、高層屋根面に巻き込まれる様子が見られており、建物近傍では新たな円錐渦の 3 次元的な流れ場を確認できた。数値流体解析の大規模化は計算精度の向上
や新たな流れ場のメカニズムの解明にも重要な役割を持っている。

図 7 建 物 近 傍の平 均 流 線 流れ場

17

建 築 分 野の数 値 流 体 解 析における
大規模計算

３．
２ 複雑な 3 次元建物形状の再現 [3]３．
３ VOF 法による広域 3 次元津波解析 [6]

最近の建築分野では高度な 3 次元データの利用により、ニーズに応

従来の津波解析では平面 2 次元浅水長波モデルを用いて、津波の遡じて建物の形状も次第に複雑となっている。数百万格子規模の解上様子や広域の津波の浸水範囲を予測している。しかし、建物レベ析モデルを得るためには数 mから数 十 cm 程 度の格 子 幅を利 用すルにおける周辺の影響と建物本体の津波波力などの詳細な予測は困る必 要がある。従 来の計 算 格 子の作 成 方 法では、3 次 元 形 状デー難である。本章では前章の結果を踏まえて、数値計算では、VOF 法 [7]タをそのまま利用できないため、数値流体解析用にデータのクリーによる気液二相流シミュレーションを用い、複雑地形を考慮できる

ニングや形 状の単 純 化などを行う必要がある。形状の詳細に伴い

広域 3 次元津波解析を実施し、その実用性を確認した。対象は沿岸

その作業も困難となる。

の街区全体とした 6km × 6kmの計算領域である。全体の計算格子対 象としたものより小さなメッシュサイズを利 用することにより簡単にその形状を再現できる。図 8 に示す複雑な 3 次元データ形

数は約 5 千万であり、実現象での時間 7 分間で、計算時間は 198CPU
で約 8 日間であった。

状は数 cm 程 度のメッシュ幅の採 用により室 内 屋 外を含む計 算モ

図 -10 には、津波が汀線付近に到達し、街区への遡上が開始する直デルを簡 単に作 成できた（ 図 9）。得られる解 析モデルは約 2 億 格

前の図を示す。右上部分の山が海岸線に迫り急峻な地形になってい子の大 規 模モデルになっているが、建 物の手 摺りや室 内の設 備なる部分では既に津波の遡上が開始しており、局所的な波の集中や複ども精緻に再現されることも分かる。

雑な反射波の発生が確認できる。図 -11 と 12 には、その後の街区での津波の遡上の様子を示すが、大きな建物のある部分では津波遡上が阻止され、路地や空地等に沿って進行する様子が見られる。
以上のように、3次元VOF法による津波遡上シミュレーションでは、地形や街区形状を詳細に考慮することで、陸域における複雑な津波の挙動を再現することが可能となる。

図 8 三 次 元 CAD データ

図 10 計 算 領 域

図 11 市 街 地の津 波 遡 上 様 子
図 9 計算格子作成例18

図 12 時 刻の変 化による建 物 近 傍の津 波 挙 動

おわりに

4

建築分野の実応用流体問題を対象とした場合、計算機の計算性能を最 大 化するためには大 規 模 並 列 処 理による大 規 模 計 算の実 施が必 要となる。また、大 規 模 計 算を行うことで、計 算 格 子の高 解
像 度 化による計 算 精 度の向 上や計 算モデルの精 緻 化の重 要 性を確 認し、さらにはこれまで困 難であったVOF 法による広 域 3 次 元
津波解析も実施できた。計算機の発展とともに新たな計算手法が
確 立されつつあり、建 築 分 野においてもその手 法に応じた新たな評価方法の確立が今後も急務な課題となっている。
謝辞
本 研 究は平 成 20 年 度と 21 年 度に東 京 工 業 大 学の共 同 利 用（産業 利 用）および先 端 研 究 施 設 共 用 促 進 事 業『みんなのスパコン』
TSUBAMEによるペタスケールへの飛 翔 トライアルユースにより実施した [6,7] 。ここに、記して謝意を表する。
参考文献[1] Top500: http://www.top500.org/[2] OpenFOAM: http://www.openfoam.com/[3] ファム バン フック：オープンソースを用いた数値流体解析，
日本流体力学会， Vol.31,pp.277-284 (2012)[4] ファム バン フック，野津 剛，菊池浩利，日比一喜：数値流
体 解 析の大 規 模 並 列 計 算 及びGPU 化の試み，ハイパフォー
マンスコンピューティングと計 算 科 学シンポジウム論 文 集，pp.75(2011)[5] ファム バン フック，野津剛，菊池浩利，日比一喜，田村幸雄：
LESによるセットバックした建 物の局 部 風 圧の検 討，風 工 学
シンポジウム論文集， Vol.22， pp.359-364 (2012)[6] ファム バン フック，長谷部 雅伸，髙橋 郁夫 : VOF 法を用いた3 次元津波解析に関する研究，土木学会論文集 B2（海岸工学），
Vol.68, pp.71-75 (2012)[ 7 ] C．W．Hirt， et al.: SOLA-VOF: A solution algorithm fortransient fluid flow with multiple free boundaries , Los
Alamos Scientific Laboratory， pp.1-32 (1980)

19

大規模地震における強震動評価と屋内収容物の被害評価
正月 俊 行株式会社 構造計画研究所

地 震による被 害 軽 減 策を考える上で、基 礎 的なデータとなる地 震の揺れの強さや建 物 内の収 納 物の挙 動を想 定するためのプログラムを TSUBAME2.0 に移 植し、大 規 模な計 算の実 行が可 能な環 境を構 築した 。地 震の揺れの強さを想 定する波 動 伝 播 解 析では、プログラムの GPU 化を行ない、CPU のみの場 合の約 10 倍の計 算 速 度が出ることを確 認した 。
また、巨 大な計 算 領 域を必 要とする東 海・東 南 海・南 海 地 震のシミュレーションを行ない、多 数の GPU を用いることで、実 用 的な時 間で計 算が出 来ることを確 認した 。建 物 内の収 納 物の挙 動 解 析については、倉 庫 内に積み上げられた積荷が崩れる様子をシミュレーションし、従来の約 2 倍の速度が計算を実行できることを確認した。

はじめに

1

重 要なファクターであることは自 明である。地 震の強さの想 定 方
法としては様々な手法が存在し、現状では、実務上の手軽さや時間、費用の制約から、簡易な手法 例えば[1] が用いられることが多い。トラ
イアルユースにて取り組んだ波 動 伝 播シミュレーションも地 震の強さの想 定 手 法の一つである。この手 法を用いると、簡 易な方 法

2011 年 3 月に発 生した東 北 地 方 太平洋沖 地 震では、非 常に広い範では考 慮 出 来ない情 報を加 味することで、より本 物の地 震 波に近囲で甚大な被害が発 生した。今後、発 生確率が高いと言われているい波（ 主に、長周期帯域 ）が生成出来る。しかし、膨大な計算量と東 海・東南海・南海 地 震が起こった場 合も同様に、都市部を含めた

計 算 時 間が必 要なため、一 部の重 要 構 造 物を除いては、ほとんど

広い範 囲で被 害が発 生すると予測されている。このような被 害を軽用いられていないのが現 状である。ただ、このような手 法を積 極減するためには、事前に、地 震による揺れや被害を予測し、対 策を立的に取り入れていくべきという指針が学会から出されており[2] 、高てておくことが重 要である。

度な手法の必要性は今後増していくものと思われる。

本 稿では、まず、安 全な建 物の設 計という観 点から、地 震に起 因

次に、地 震 時の屋 内 被 害 評 価に関する実 務の現 状について述べする被 害 評 価に関 連する実 務の現 状と課 題について述べる。続いる。建物の構造設計とは異なり、室内の家具の転倒防止や収納物て、その課 題に対して、平成 23 年度のTSUBAME2.0 の産 業利用トラの落 下 防 止などに関 連する厳 密な取り決めは、現 状では存 在しな

イアルユースにて、実務への適 用が困難な地 震 動 評 価手法と屋内被い。ただ、近年、特に、東北地方太平洋沖地震以降、このような屋

害評 価 手法の実 用性向上に取り組んだ結果について述べる。

内 被 害に対する意 識が高まってきており、屋 内 被 害に関する数 多トライアルユースでは、まず、地 震 動の波 動 伝 播 解 析プログラム

くの被害調査報告が出され例えば[3], [4] 、対策に関する検討会 例えば[5] 等と屋 内 被 害のシミュレーションプログラムの 2 種 類のプログラムをも開かれている。屋 内 被 害のシミュレーションを用いると、詳 細TSUBAME2.0 に移植した。そして、移植したプログラムを用いて、東海・な情 報を反 映して、リアリティのある被 害 評 価を行なうことが出

東 南海・南海 地 震の波 動 伝 播シミュレーション、および、屋内 被 害来る。しかし、精 度のよい計 算をしようとすると膨 大な時 間が必評 価の一 例として、倉庫内の積荷が崩れる様 子のシミュレーションを要になるため、実務で用いられることは少ないのが現状である。計行なうことで、TSUBAME2.0 を用いた実 用的な計算が実 行出来るこ

算 時 間を大 幅に短くできれば、屋 内 被 害 対 策を立てるための検 討との確 認を行なった。

をより手軽に行なうことが出来るようになると思われる。
構造計画研究所では、波動伝播解析の実用性向上に向けて、平成 21 年度から取り組みを行なっている。初年度は、地球シミュレー

実務における地震被害評価の現状

2

タのトライアルユースを利 用し、プログラムの移 植を行なった。
その次の年 度では、プログラムの最 適 化や機 能 追 加を重ね、大 幅
な計 算 速 度の向 上を実 現するとともに、より短い周 期 帯 域に対する計算精度向上を目的とした検討も行なった。そして、平成 23 年
度は、屋内被害シミュレーションを加え、さらなる高速化・実用性

本題に入る前に、安全な建物の設計という観点から、地震に起因す向上を目指して、プログラムをGPU 化してTSUBAME2.0 に移植をる被害評価（ここでは、地震の強さの想定と屋内被害評価に限る）

行なった。

に関連する実務の現状と課題について述べる。その後、平成 23 年
度 TSUBAMEトライアルユースを行なうに至った背 景と経 緯についても述べる。
地 震の強さの想 定が、建 物の耐 震 安 全 性を考える上で、非 常に

20波動伝播解析3

考に行なった。その結果、GPUを用いると、CPUのみの場合（-fastsse
オプションを使用）と比較して、10 倍程度計算速度が向上（約 31
億メッシュ規模、100ノード、300GPUの場合）することが確認できた。図 1 にプログラムの弱スケーリング性能を示す。並列数を 600にしても、ある程度性能がスケールしていることが分かる。現状でこれより、平成 23 年度 TSUBAMEトライアルユースの結果について

は、MPIの通信時間がボトルネックになっているため、さらなる高述べる。

速化を図るためには、通信量を小さくするか、通信時間を短くする

まず、地 震動の波動伝播解析プログラムのTSUBAME2.0 への移ことが必要である。また、GPUのグロバールメモリがさらに大きく

植・GPU 化について述べる。このプログラムは、3 次元差分法を用なれば、同じ規模の計算をより少ない並列数で実行できるので、効いて波 動 伝 播を計 算する。TSUBAME2.0 への移 植にあたっては、率が上がると思われる。なお、今回計算を実行した最大のメッシュ

既にMPI 化により複数ノードでの並列化がなされていたため、主に

規 模は約 120 億メッシュ規 模で、ステップ数 77600（継 続 時 間 約

各プロセスにおける計 算のGPU 化を行った。なおコンパイラは、

232 秒）の計算をGPU 810 基（270ノード）を用いて、4 時間程度でFortranからCUDAを用いるため、PGI Fortranコンパイラを用いた。

解くことが出来ることを確認している。

GPU 化は、同一の計算のGPU 化を実施している竹中（2011）[6] を参続いて、改造したプログラムを用いて行なった東海・東南海・南※ 周 期 2 秒より短い周 期をフィルターでカットしている

図 1 波 動 伝 播 解 析プログラムの弱スケーリング性 能

図 3 名 古 屋における速 度 波 形

図 2 東 海・東 南 海・南 海 地 震の震 源 域と
計算領域

21

大 規 模 地 震における強 震 動 評 価と屋 内 収 容 物の被 害 評 価

海 地 震の 3 連 動の地 震の波 動 伝 播シミュレーションについて述べ

揺れはじめ、それ以降、300 秒後も揺れがほとんど収まっていないこる。図 2 に想 定 地 震の震 源 域と計 算 領 域を示す。震 源は、中 央 防とが分かる。図 4 には波動伝播の様子のスナップショットを示す。

災会議（2003）の震源情報 [7] 、および、山本・吉村（2012）[8] を参考

地表面の速度の大きさをコンターで示している。断層の破壊開始に設定した。計算規模を 31 億メッシュ規模（550x750x60km）、
メッから約 75 秒で大阪平野や濃尾平野に波が到達し、その後、平野でシュ間隔を 0.2x0.2x0.2km 、計算時間刻みを 0.006 秒とし、計算ス

揺れが長く続いており、この部分からも平野部で強い揺れが長く続

テップ数は 50000（継続時間 300 秒）とした。計算可能周期は、2 秒

く現象をシミュレーションで再現できていることが分かる。今後は、以上である。図 3 に名古屋における速度波形を示す。50 秒前後で計算結果の妥当性について検証を行なっていく。

図 4 波 動 伝 播 解 析の様 子（地 表の速 度をコンターで表 示）

22屋内被害評価

4

屋 内 被 害 評 価については、剛 体の 3 次 元 物 理 運 動シミュレーショ
ンプログラムを用いて、倉 庫 内の積 荷の挙 動を計 算した。このプ
ログラムを用いると、地震の揺れにより、倉庫の収納ラックから物が滑り落ちる様子などがシミュレーション出 来る（図 5 参 照）。剛体同士が衝突した際の接触力については、拘束条件を運動方程式

図 5 自 動ラックからの収 納 物の落 下 解 析

と連立させて解く解析法 [9],[10],[11] を用いて計算を行なっている。大まかな計算フローを以下に示す。
1. 衝突判定：剛体の位置・速度・形状から衝突判定を行なう。

2. 拘束条件の更新：接触判定の結果に基づいて、接触による拘束を追加する。
3. 拘束力の計算：拘束条件が与える拘束力を計算する
4. 運動方程式に基づき、速度・位置を更新するプログラムをTSUBAME2.0 に移植するにあたって、まず、プログラ
ムを大 規 模なデータに対 応させる改 造を行なった。改 造 前のプロ
グラムは、剛体の数が 1000 個程度以上となると、計算がストップしてしまうような仕様であったが、改造後は、10000 個程度以上でも計 算が実 行できるようになった。プログラムの並 列 化は、IntelのTBBを用いて行ない、複 数のGPUを用いて計 算が実 施 出 来るよ
うに改 造を行なった。計 算 時 間の大きな部 分を占めていた剛 体の図 6 倉 庫の積 荷の地 震 時 挙 動
シミュレーションの計 算モデル

衝突判定の処理部分をメインにGPU 化を行ない、この部分については、CPUに比べて 3 ～ 4 倍の計算スピードとなることを確認した。
全体の計算速度としては、改造前の 2 倍程度となった。
倉庫内の積荷の地震時挙動のシミュレーションの計算モデルを
図 6 に示す。約 1m 角のパレット（ 図 中 赤 線で表 示 ）の上に、麻 袋を想定した剛体（ 図中白線で表示 ）48 個を、方向を変えながら積みおわりに上げたモデルを一塊として、さらにその塊を2列分積み上げている。

5

片方の列は、壁際に設置している。全体の剛体数は約 1400 個である。今回の課題では、最 大で剛 体 数 10000 個 強のケースが実 行 可

本 稿 で は、平 成 23 年 度 トライアルユ ー ス に て、地 震 動 の 波 動

能であることを確認しているが、数十秒間のシミュレーションに 2

伝 播 解 析 プログラム と 屋 内 の 収 容 物 の 挙 動 解 析 プログラム を

週 間 以 上かかる上、剛 体の数が増えると計 算が不 安 定になりやす

TSUBAME2.0 に移 植し、大 規 模な計 算を従 来よりも短い時 間で行いことが分かったため、剛体数を少なめに設定した。

なえる環 境を構 築した結 果について述べた。今 後は、これらのプ

図 7 に、図 3 で示した波形を計算モデルの床に入力した場合のシ
ミュレーション結 果を示す。図を見ると、揺れが強くなる 75 秒を

ログラムを実 務に生かしていくことを念 頭に、計 算 結 果の工 学 的
妥当性の検証を行なっていく予定である。

過ぎたあたりで、積み上げていた積荷が崩れ始め、崩れた壁際の積荷が隣の単 独 配 置の列に覆いかぶさっている様 子が見て取れる。

謝

今後は、計算結果の工学的妥当性について検証を行なっていく必

本 稿に記 載した計 算は、文 部 科 学 省の先 端 研 究 施 設 共 用 促 進 事要がある。

業の補助をいただき、東京工業大学 学術国際情報センターのスパ

辞

コンTSUBAME2.0 を用いて行なった。記して謝意を表す。
参考文献[1] 翠川 三郎、小林 啓美 : 地震断層を考慮した地震動スペクトルの推定 , 日本建築学会論文報告集 , Vol. 282, pp.71-81, 1979.

23

大 規 模 地 震における強 震 動 評 価と屋 内 収 容 物の被 害 評 価[2] 日本建築学会 : 最新の地盤震動研究を活かした強震波形の作
成法 , 2009.
[3] 東 京 消 防 庁 : 平 成 23 年（ 2011 年）東 北 地 方 太 平 洋 沖 地 震 調
査報告書 , 2012.
[4] 日 本 建 築 学 会 : 長 周 期 地 震 動 対 策に関する公 開 研 究 集 会 ,2012.
[5] 気 象 庁 地 震 火 山 部 : 長 周 期 地 震 動に関する情 報のあり方 検
討会 , 2012.
[6] 竹 中 博 士 : 並 列 GPUを用いた大 規 模 地 震 波 伝 播シミュレー
ション, 学際大規模情報基盤共同利用・共同研究拠点
平成 23 年度共同研究中間報告書 , 2011, 12.
[7] 内閣府中央防災会議 :「東南海 , 南海地震等に関する専門調査会」（第 16 回）資料 3 図表集 2, 2003, 12.
[8] 山本 優・吉村 智昭 : 3 次元大規模 FEMによる東海・東南海・
南海連動地震の長周期地震動シミュレーション, 日本建築学
会構造系論文集 , Vol. 7, No. 677, p.1055-1064, 2012.
[9] D. Baraﬀ: Analytical Methods for Dynamic Simulation of
Non-Penetrating Rigid Bodies, Computer Graphics Proceedings(SIGGRAPH 89) Vol. 23, pp. 223-232, 1989.
[10] Kenny Erleben: Multibody Dynamics Animation Lecture 12,
h t t p : / / w w w. d i k u . d k / O L D / u n d e r v i s n i n g / 2 0 0 5 f / 1 0 1 /lecture12.pdf .
[11] 長谷川晶一・田崎勇一・三武裕玄・他 : “Springhead2”,
http://springhead.info/wiki/index.php .

図 6 倉 庫の積 荷の地 震 時 挙 動シミュレーションの結 果

24

● TSUBAME e-Science Journal No.8
2013 年 3 月12 日 東京工業大学 学術国際情報センター発行 ©
ISSN 2185-6028
デザイン・レイアウト：キックアンドパンチ
編集： TSUBAME e-Science Journal 編集室
青木尊之 ピパットポンサー・ティラポン
渡邊寿雄 佐々木淳 仲川愛理
住所： 〒 152-8550 東京都目黒区大岡山 2-12-1-E2-1
電話： 03-5734-2087 FAX：03-5734-3198
E-mail： tsubame_j@sim.gsic.titech.ac.jp
URL： http://www.gsic.titech.ac.jp/

25

vol.

TSUBAME 共同利用サービ ス『みんなのスパコン』TSUBAME 共同利用サービスは、
ピーク性能 2.4PFlops 、18000CPUコア、4300GPU 搭載
世界トップクラスの東工大のスパコンTSUBAMEを
東工大以外の皆さまにご利用いただくための枠 組みです。

課 題 公募する利 用 区 分とカテゴリ
共同利用サービスには、「学術利用」「産業利用」、「社会貢献利用」、の3つの利用区分があり、さらに「成 果公開」と「成 果非公開」のカテゴリがあります。
ご利用をご検 討の際には、下記までお問い合わせください。

他大学や公的研究機関の研究者の 学術利用［有償利用］

TSUBAME

民間企業の方の 産業利用［有償・無償利用］

共同利用とは…

その他の組織による社会的貢献のための 社会貢献利用［有償利用］

共同利用にて提 供する計算資源
共同 利 用サービスの利 用区 分・カテゴリ別の利 用課 金 表を下 記に示します。TSUBAME
における計算 機 資 源の割 振りは口 数を単 位としており、1 口は標 準 1ノード
（12 CPUコア,
3GPU, 55.82GBメモリ搭載）の 3000 時間分（≒約 4 ヵ月）相当の計算 機資 源です。
1000 CPUコアを1.5日利用する使い方や、100 GPUを3.75日利用する使い方も可能です。
利用区分

利用者

制度や利用規定等

カテゴリ

利用課金（※）

学術利用

他大学または研究機関等共同利用の利用規定に基づく

成果公開

1口：100,000 円

産業利用

民間企業を中心としたグループ「先端研究基盤共用・プラットフォーム形成事業」に基づく

非営利団体、公共団体等共同利用の利用規定に基づく

社会貢献利用

成果公開

トライアルユース（無償利用）
1口：100,000 円

成果非公開

1口：400,000 円

成果公開

1口：100,000 円

成果非公開

1口：400,000 円（※） 利用課金は平成 24 年度の価格で、平成 25 年度より値上されます。平成 25 年度の価格は、下記の URL をご参照ください。
http://www.gsic.titech.ac.jp/node/58

産業利用トライアルユース制度（先端研究基盤共用・プラットフォーム形成事業）
東 工 大のスパコンTSUBAMEを、より多くの企 業の皆さまにご利 用いただくため、初めて
TSUBAMEをご利 用いただく際に、無 償にてご試 用いただける制 度です。
（文部 科 学省 先 端 研 究 基 盤 共 用・プラットフォーム形成事業による助成）

詳しくは、下 記までお問い合わせください。

お問い合わせ●東京工業大学 学術国際情報センター 共同利用推進室●e-mail kyoyo@gsic.titech.ac.jp Tel. 03-5734-2085 Fax. 03-5734-3198
詳しくは http://www.gsic.titech.ac.jp/tsubame/をご覧ください。

8