（様式第 20）成果報告書

平成 25 年度 産業利用トライアルユース：先端研究基盤共用・プラットフォーム形成事業『みんなのスパコン』TSUBAME による日本再生
成果報告書

利用課題名 衛生陶器設計のための並列 GPGPU 気液二相流シミュレーション
英文：Numerical simulation of air-water two-phase flow conducted on GPGPU cluster
for designing sanitary ware products
池端 昭夫
Akio Ikebata

吉田 慎也
Shinya Yoshida

TOTO 株式会社 生産技術本部
Production Technology Division, TOTO LTD.
http://www.toto.co.jp
邦文抄録（300 字程度）
水まわりの住宅設備機器において、節水化のニーズの高まりに対応するため、効率的かつ適切に設計を行うための混相流シミュレーションプログラムを開発している。衛生陶器においては節水化に伴う流れの薄膜化に対応するため、TSUBAME を用いた大規模細分化メッシュにおいて陶器表面および流路内部のシミュレーションを試みた。独自の解析手法 UTI-VSIAM3 において、OpenCL と MPI を併用した大規模並列 GPGPU 流体計算コードを開発した。通常の PC や WS の 800 倍以上の計算速度を達成し、ミリ単位の薄膜流れや微小気泡流れを再現することに成功した。本技術は既に衛生陶器において節水商品の洗浄性能向上に適用されている。
英文抄録（100 words 程度）
We have been devoting our efforts to development of simulation tools on multiphase flow for
designing water-saving house hold equipment. In order to get reliable simulation results for complex
multiphase flows, such as thin layer flow on ceramic sanitary ware, we have carried out
high-resolution computations with large number of mesh elements on the TSUBAME supercomputer.
Large-scale parallel GPGPU simulation code based on the UTI-VSIAM3 Navier-Stokes solver has beendeveloped and ported with OpenCL and MPI libraries. The present code achieves over 800 times
speed-up compared with a normal PC or a workstation, which enables us to numerically reproduce thethin layer flows and bubbles of millimeter order. With the innovative performance of this state of theart simulation tool, the design process of the water-saving sanitary ware products will be greatly
accelerated in real-case applications.

Keywords: multiphase flow, house hold equipment, sanitary ware, UTI-VSIAM3, CIP

背景と目的模クラスターに適用することで、実用的な計算時間およ

当社では TOTO GREEN CHALLENGE のスロー

び計算精度で衛生陶器流れのシミュレーションを実施

ガンのもと、衛生陶器や水栓、シャワーなどの住宅水まできるようになった。本技術は既に多くの商品開発に適わり機器の節水商品開発に力を入れている。

用しており、試作評価の削減に貢献している。

衛生陶器においては、汚物の洗浄性を確保しつつ大本課題では、より一層のメッシュ細分により、社内シ幅な節水化を達成するための流路や陶器形状、水のステムでも再現できなかったミリ単位の薄膜流れや微供給条件を高度に適正化する必要があり、試作評価の小気泡を動きを再現し、試作評価の代替が可能なレベ

みならず、多相流体解析の積極的な活用を進めている。

ルの精度を実現するため、TSUBAME2 を用いてメモ

気液二相流体解析の現状の問題点としては、薄膜流リ分散並列+GPGPU による大規模流体シミュレーショ

れや微小気泡などが時々刻々大きく変動するため、計ン技術を開発した。

算メッシュ数やステップ数が膨大となることである。そこで我々は計算精度および効率に優れる独自流体計算手法を開発し、社内の GPGPU を搭載した PC や小規概要
気液二相流の流体方程式は以下のように表される。

（様式第 20）成果報告書

u11
 u  u  P  u  Tt


(1)

P2
 u  P  CS   u
t(2)る。各時刻ステップごとに図 2 の手順にて、VIA、SIAの時間発展を行う。まず直交格子の 6 面に流入するフラックスを数値積分により計算する。数値積分の方法
は図 2 ではシンプソン積分を用いているが、ガウス積分

F(3) u  F
tu, P, F は流速、圧力、流体率を表す。右辺第一項はでも同等精度が得られる。次に SIA をセミラグランジェ全て同じ形であり、移流項とよばれる。この項を精度良る こ と か ら 、 非 次 元 分 割 手 法 （ UTI=Unsplit Time

く計算することが流体解析における大きな課題である。

Integration）と呼んでおり、次元対称性を持つ。の手法で更新する。最後にフラックスの収支から VIAを更新する。これらは X、Y、Z 軸方向に同時に計算す

そこで UTI-VSIAM3 ソルバー[1]を提案している。
UTI-VSIAM3 では、VIA(=体積積分平均値）と SIA（=面積積分平均値）という 2 種類の物理量変数を定義する。通常の有限体積法では VIA のみを定義する。図１のように、一つの直交格子メッシュにおいて、X、Y および Z 方向の SIA を各々X-SIA、Y-SIA、Z-SIA と定義する。これにより、X 方向は 2 つの X-SIA と VIA で二次補間関数を構築でき、他の一般的な手法とは異なり一つの格子で高次精度化できるのが利点である。本
手法は一次元保存型 CIP 法とよばれる[2]。Y 方向、Z方向についても同様である。
一次元保存型 CIP 法の多次元化は様々な手法が提図１．VIA および SIA の定義

案されているが、UTI-VSIAM3 では、保存性および三次精度を確保しつつ、省メモリおよび簡易性を利点とす

図２．UTI-VSIAM3 移流方程式解法

（様式第 20）成果報告書
これらの計算を実行するためには、任意位置の SIA、

図４に、UTI-VSIAM3 ソルバーの計算手順を示す。

すなわちを求める必要がある。

移流項、粘性項、生成項、音響項の順番にフラクショナ

我々は、以下の 2 ステップ分割法を提案している。例とルステップ法で計算する。これらの計算は全てして SIA-X について表記するが、SIA-Y、SIA-Z につ

GPGPU 化が容易に可能なベクトル計算の集合体で構いても同様である。

成されている。流体率 F の気液界面の数値拡散を抑制[STEP1] 各面法線方向の一次元保存型 CIP 法するため、移流計算後に STAA 法[3]により界面法線

 x1   xtt x
 u '
xt(4)方向への補正移流を行う。

[STEP2] 各面面内方向の中心差分計算

 x 2   x1t

 (v'  x )1  ( w'  x )1yz v' w'   x1  y z (5)STEP2 の右辺のフラックス値は、STEP1 の保存型
CIP 法計算時に中間値として得られるフラックスを流用でき、高次精度で計算できることが利点となる。 図 3に、移流方程式のテスト計算結果を示す。φ の初期分
布で TOTO の T の字の内側を 1、外側を 0 とおき、剛体回転流速場でφ を一回転させた。市販ソフトウェアなどで安定的かつ高精度な保存型解法として広く用いられている MUSCL 法と比較し、明らかに高精度であることが分かる。

図４.UTI-VSIAM3 ソルバーの計算手順

ダム崩壊問題における気液シミュレーションの実施結果を図５に示す。飛散液滴や巻き込み気泡が高い解像度で再現できている。また STAA 法の効果により、気液界面は終始一定厚みに保たれ、数値拡散は生じていない。またミルククラウン[4]のシミュレーション結果を図
図３．移流テスト計算結果６に示す。着水時の液滴落下速度を 2.2m/s と設定し、牛乳の粘性係数および表面張力係数を用いている。こ（様式第 20）成果報告書
れまでは表面張力を考慮して王冠形状をシミュレーションで再現することは難しいといわれていたが、本解析手法では比較的良好な結果が得られた。

図７．PCG ポアソン方程式解法の GPGPU 化の効果
TSUBAME2 への実装では、並列化効率の向上が大きなポイントとなる。そこで、図８に示す三次元領域分割法および図９のオーバーラップ通信手法を用いた。
一次元領域分割法と比較して、計算メッシュ数に対するデータ通信メッシュ数の割合が大幅に小さくなるため、並列化効率が向上する。また実製品の複雑形状に対応するため、隣接分割領域を不整合とする手法を採用した。オーバーラップ通信は、計算とデータ通信を並行
図５．ダム崩壊シミュレーション結果して行うが、データ通信は非同期かつ一斉に発行し、MPI_Waitall にて通信の終了を待つことで、不整合境界での通信を簡易に実現した。

図６．ミルククラウンシミュレーション結果
UTI-VSIAM3 ソルバーの GPGPU 化を行った。コー

図８.領域分割法

ドの長寿命化を狙い、標準規格の OpenCL を用いてコーディングした。NVidia 社の CUDA デバイス上で高速に実行できるようにするため、OpenCL 関数を CUDA
Driver API に変換するマクロを作成した。また複雑なOpenCL カーネル制御関数を簡略化するため、TOTO
CL ライブラリとよんでいる独自のラッパー関数群を定義して、流体コードでは TOTO CL ライブラリを通してカーネル関数を呼び出すようにした。GPGPU 化による効果を図７に示す。CPU との比較では、1 コアあたり約 8～10 倍の高速化を達成している。

図９.オーバーラップ通信手法

（様式第 20）成果報告書
混相流シミュレーションでは、計算時間の大半を圧力結果および考察ポアソン方程式の求解で占める。TSUBAME2 のようTSUBAME2.0 上に PCG 法の圧力方程式解法を実な GPGPU と分散メモリ並列の併用アーキテクチャで装したときの計算速度結果を図１１に示す。1 メッシュあは、マルチスレッド化およびマルチプロセス化が容易なたり 0.5mm 間隔の衛生陶器モデル（1 億 8260 万メッ

PCG 法が有効である。しかし PCG 法はスケーラビリテシュ）において、約 150GPU で約 100 倍の計算速度がィに優れているものの、収束性が悪く反復回数が増大得られており、実用的なスケーラビリティを達成した。図

し、計算時間を要する。そこで CG 法の前処理にマルチ７で GPU は CPU の約 8～10 倍の計算速度であるこグリッド法を用いる MGCG 法についても検討した。マとから、150GPU で 1CPU の 800 倍以上の高速化を

ルチグリッド法は並列計算に適した Red Black SOR を得られたといえる。TSUBAME2 は 1 ノードで 3GPU 構用いた。圧力項は下式のように体積積分形で表す。

成であり、わずか 50 ノードの規模で PC の 1000 倍近1S   P  n dS  VrdV

い計算性能を実製品シミュレーションで実現している。

(6)

一般的な CPU クラスタ型スーパーコンピュータであれば数百台の計算ノードが必要と推測され、TSUBAME2 はコスト、消費電力、安定性、メンテナン

圧力の連立方程式を下のように置く。

ス性全ての面において、企業利用では圧倒的に優位で6

C7 Pi , j ,k   Cnb Pnb  R(7)あるといえる。

nb 1nb は隣接メッシュをあらわす。体積積分形であるため、Fine Grid の係数 CF を加算することで、Coarse Gridの係数 CC を容易に導出できる。例えば図１０に示す上面の係数は以下のように計算する。

CiC, j ,k 1 FC2i1, 2 j 1, 2 k  C2Fi , 2 j 1, 2 k  C2Fi , 2 j 1, 2 k  C2Fi , 2 j , 2 k2

(8)
右辺の 1/2 は、式(6)の左辺の被積分関数に∇、すなわ
ちメッシュ間隔の逆数が含まれており、Coarse Grid のメッシュ間隔は Fine Grid のメッシュ間隔の 2 倍であることをあらわしている。この定式化では、Fine Grid から Coarse Grid への制限補間では、Fine Grid の残差

図１１．TSUBAME2.0 における PCG 法の速度向上を加算するだけで済み、また Coarse Grid から Fine（1 億 8260 万メッシュ)

Grid への延長補間は、Coarse Grid の修正ベクトルをそのまま代入するだけでよい。

MGCG 法を GPGPU 搭載 PC に実装し、PCG 法と比較した結果を表１に示す。1 メッシュあたり 1.5mm 間隔の小規模な衛生陶器モデル（5,697,794 メッシュ）にて実施した。MGCG 法は反復回数が 2,485 回から 24回と 1/100 程度に削減し、計算時間も 1/5～1/6 となっている。したがって MGCG 法は衛生陶器の混相流シミュレーションにおいても、計算時間短縮に極めて有効であることが分かった。
また 1CPU と 4GPU との比較では、1.9/57.3=約 1/30

図１０．Fine Grid からの Coarse Grid 生成となっており、同コア数で換算すると 30/4=約 7.5 倍の（様式第 20）成果報告書計算速度となっている。図７の結果と比較しても概ね妥

衛生陶器流れのシミュレーション結果を図１３に示す。

当といえる。すなわち GPGPU 計算は MGCG 法にお微小気泡や薄膜流れが高解像度で再現できていることいても高い性能を発揮できることが分かった。

が分かる。本解析手法では、気泡の巻き込みや運動、分裂などの過程は一切モデル化することなく、式(1)か

表１.MGCG 法と PCG 法との比較（5,697,794 メッシュ）

PCG
MGCG
反復回数 2,485 反復回数 24
1CPU(Xeon3.6GHz)
2GPU(Geforce 690)
4GPU
(Geforce 690)

342.9 sec

57.3 sec

18.5 sec

3.9 sec

9.4 sec

1.9 sec

ら(3)の Navier-Stokes 方程式をメッシュ解像度限界の中で可能な限り精度良く解くことで、自然に一連の現象を捉えようとしているため、様々な製品の気液二相流シミュレーションに適していると考える。

本プログラムを TSUBAME2.5 にインストールし、スケーラビリティを調べた結果を図１２に示す。PCG 法と比較して MGCG 法は 30GPU 以上でスケーラビリティが低下している。これは、Coarse Grid での Red BlackSOR 法の計算において、メッシュ数が尐ないため、相対的に通信時間の割合が高くなったためと考えられる。

図１３．衛生陶器のシミュレーション結果

しかしながら MGCG 法の 70GPU で PCG 法 1GPUの 250 倍以上、また 150GPU で 300 倍以上の計算速衛生陶器表面の水の流れのようすを、従来の PC で度向上率が得られており、極めて実用的な解法であるのシミュレーション結果と比較した。図１４で明らかなよことが分かった。

うに、水の旋回のフローパターンが大幅に実機に近づいたことが分かる。

図１４.衛生陶器表面流れの比較結果

解析精度の定量的検証として、ふちの内側に生じる乗り上がり（水流の盛り上がり）の高さと位置を実機と比較した。図１５のように、旋回曲率半径が小さいふち前部右側において、水が遠心力で陶器横壁に強く押し付けられ、水が上部に盛り上がる。検証結果を図１６～１７に示す。実測値もばらつきが大きいことを考慮する
図１２．TSUBAME2.5 での圧力方程式解法結果と、TSUBAME2 での 0.5mm メッシュの結果は従来の（152,347,494 メッシュ）1.5mm メッシュと比較して実機に極めて近い結果といえる。図１７の乗り上がり水平方向位置もほぼ妥当な結（様式第 20）成果報告書
果である。

確に表現できないため、薄膜部位における解析信頼性が低く、実験での確認が必要な場合が多かったが、今回の 0.5mm メッシュシミュレーションにより実機再現性が飛躍的に高まったため、洗浄性と節水性を両立できる設計案を試作評価前に導出することが可能となった。

図１５.「乗り上がり」現象

図１８.衛生陶器表面流れ検討への適用のようすまとめ、今後の課題
図１６.乗り上がり高さの実験との比較

衛生陶器の表層流れや気泡混入流れはミリレベル以下の細分格子が必要であり、膨大な計算量となる。
また非定常流れであるため、時間方向の計算サイクル数も極めて多く、一般的な PC ベースの流体解析技術
では信頼性の高い解を得ることは困難である。計算時間の大幅な短縮のためスパコンを利用することは容易に着想できるが、投資対効果の観点より従来のスパコン の 利 用 は 有 効 と は い え な い 。 TSUBAME2 はGPGPU ベ ー ス の 大 規 模 並 列 ス パ コ ン で あ り 、GPGPU 自体が CPU から見るとスパコンレベルの演算性能を持っているため、その能力をフルに発揮でき

図１７.乗り上がりの水平方向発生位置の比較れば、投資対効果の問題を乗り越えられる切り札になると考えられる。

図１８に、衛生陶器表層流れシミュレーションの適用

我々は GPGPU 並列計算システムの性能をフルにのようすを示す。流速の方向および値を図のように可発揮するための GPGPU 化・並列化手法、さらに細分視化することで、陶器表面の局所的かつ全体的な洗浄格子において最大限の計算精度を発揮するための混状態を明確にした。陶器内部の流路構造やタンク・ポン相流解法の開発に尽力してきた。その結果、プなどからの給水条件を種々変更したときの流速分布

UTI-VSIAM3 により最新の計算手法と同等以上の計の 変 化 を ケ ー ス ス タ デ ィ で き る 。 従 来 の PC で の算精度を実証し、また OpenCL、三次元領域分割およ1.5mm メッシュの解析では、数ミリレベルの薄膜は正びオーバーラップ通信を併用した計算手法により実用

（様式第 20）成果報告書的な計算効率、スケーラビリティを達成することができた。一般的な CPU クラスタベースのスパコンと比較して、おおよそ 1/10 の規模で同程度の計算性能を発揮できていると推測され、スパコン技術が住宅設備機器メーカーの商品開発の現場に適用できる投資対効果の水準に初めて到達したと考えられる。
今後は、衛生陶器設計における様々な製品仕様を迅速に事前評価できるようにするための機能追加や、より一層の計算精度、計算速度の向上、さらには排水管も含めたより大きな系での大規模シミュレーションやその他の水まわり住宅設備商品設計への適用拡大を図っていきたい。

参考文献[1] A.

Ikebata,

MultiphaseY.
Fluid

Muraoka,F.

Simulations

Xiao:ona

Multiple GPGPU PC Using Unsplit Time
Integration

VSIAM3,ProgressinNUCLEAR SCIENCE and TECHNOLOGY,
2, pp.491-497 (2011)[2] T. Yabe, R. Tanaka, T. Nakamura and F.
Xiao,"Exactlyconservative

semi-Lagrangian scheme (CIP-CSL) in one
dimension," Monthly Weather Review, 129,332 (2001)[3] 池端昭夫, 肖鋒: 保存型自由表面捕獲スキー
ムと固液気 3 相流への適用, 日本機械学会
2002 年度年次大会, 3, pp.301-302 (2002)[4] 郡司、石井、斉藤、酒井、“ミルククラウンに関する研究” 日本流体力学会ウェブサイト「ながれマルチメディア」(2003)