学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

13-NA24

防災・環境のための大規模流体シミュレーションとその可視化
樫山和男（代表者：中央大），牛島省（京都大），松本純一（産総研），藤岡奨（水資源機構），京谷孝史（東北大），寺田賢二郎（東北大），中畑和之（愛媛大），浅井光輝（九州大），永井学志（岐阜大），高瀬慎介（東北大），車谷麻緒（茨城大），加藤準治（東北大），小山田耕二（京都大），岩下武史（京都大）
概要 防災・環境シミュレーションにおける津波氾濫解析として，東日本の海域（沖合部），沿岸部といった広範囲な領域から内陸部（都市部）の局所的領域を一体的に計算可能な大規模並列計算技術の検討・開発を行う．また，環境流体解析として，沿岸部付近の湖沼をとりあげ，気液二相流れと密度変化を考慮した 3 次元流体解析技術の検討・開発を行う．津波氾濫解析，環境流体解析ともに実スケール地形での大規模並列計算を実施し，大規模データを有効利用した高効率かつ効果的な可視化の検討を行う．

1. 研究の目的と意義

本研究では，防災・環境シミュレーションにお東日本大震災の大津波や日本でも発生し得る巨ける津波氾濫解析として，東日本の海域（沖合部），大台風などの大規模自然災害により，広範囲な領沿岸部といった広範囲な領域から内陸部（都市部）域を対象とした流れ現象，かつ，沿岸部や都市部の局所的領域を一体的に計算可能な大規模並列計の局所的領域における影響の把握や予測などを目

算技術，環境流体解析として，沿岸部付近の湖沼的とした防災・環境シミュレーションへの重要性をとりあげ，気液二相流れと密度変化を考慮した 3が近年高まっている．ここで，日本全体あるいは次元熱流体解析技術，津波および環境シミュレー

東日本や西日本の大部分の領域を計算対象としたションの共通課題である流体中における複数物体シミュレーションにおいては，膨大な計算量を要の輸送解析技術の開発を目的とする．また，津波する大規模解析となり，スーパーコンピュータな氾濫解析，環境流体解析ともに実スケール地形でどの超高速計算機の利用が不可欠となる．日本全の大規模並列計算を実施し，大規模データを有効体，東日本，西日本などの広域な計算を行いなが利用した高効率かつ効果的な可視化の検討を行う．

ら，都市部など局所的領域での詳細計算を実施す

近年，深刻かつ重大な問題となっている大津波

るうえで，広域では 2 次元（準 3 次元）解析，局や巨大台風による洪水などの大規模自然災害に関所的領域では 3 次解析を用いた合理的な数値解法する研究課題であり，今後の国民の安全で安心なが開発できれば，実用的な時間で計算結果を得る社会環境のための基盤研究として意義がある．

ための非常に有効な手段となり得る．また，局所的領域での 3 次元詳細計算においては，温度によ2. 当拠点公募型共同研究として実施した意義

る密度変化の考慮が現象の再現性に重要な役割を(1) 共同研究を実施した拠点名および役割分担果たす場合や流れが構造物に作用する力の影響の本研究では「学際大規模情報基盤共同利用・共把握が重要となる場合も多い．一方，近年の可視同研究拠点」のネットワーク型拠点のうち，「京都化技術の進歩により，非常に写実的な可視化や立大学」を共同研究先としており，計算機環境とし体視表現が可能になりつつある．近年の可視化技て AMD Opteron 6000 シリーズのアーキテクチャに術と計算結果（数値データ）の有効利用により，基づく CRAY XE6 を利用した．研究の実施内容は次

数値シミュレーションで得られた現象を疑似的にの 6 つに大別される．

体験するといった試みが行われている．
1

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

1) 実スケール地形での Boussinesq 方程式による(3) 当公募型共同研究ならではという事項など

2 次元浅水流れにおける並列化，大規模解析

日本全体あるいは東日本や西日本の大部分の領

2) 実スケール地形での Navier-Stokes 方程式に域を計算対象とした広範囲なシミュレーションによる 3 次元気液二相流れにおける並列化，大規

おいては，上記「(1) 1)～3)」に示した計算手法

模解析の開発が必要であり，沿岸部，内陸部の計算では，

3) Boussinesq 方程式，Navier-Stokes 方程式（3構造物に作用する力の評価，密度変化や物体の輸

次元気液二相流）による結合手法の検討および送などの「(1) 4), 5)」における 3 次元の詳細な並列計算法の開発計算が必要となる．
「1)～5)」における実用的なシ

4) Navier-Stokes 方程式（3 次元気液二相流）にミュレーションを実現するには，大規模並列計算基づく密度変化を考慮した計算手法による並機および計算の高速化が，得られた大規模かつ複列化，大規模解析

雑な計算結果を数値計算が専門でない一般の方も5) Navier-Stokes 方程式（3 次元気液二相流）にわかるような分かり易い表現を実現するためには，基づく流体中における複数物体の輸送計算手効率的な可視化技術と効果的な可視化が必要であ法による並列化，大規模解析る．以上に示した研究を目的とする本課題では，

6) 沖合部，沿岸部，内陸部における効率的な可視流体解析と構造解析分野の研究者，計算機科学と化および VR 技術を用いた効果的な可視化の検可視化分野における研究者との体制および連携は討不可欠である．

「1), 2), 6)」については主に樫山（中央大），「3)」3. 研究成果の詳細については主に松本（産総研），「4), 5)」については主に牛島（京都大），藤岡（水資源機構）が担

本研究では，図-1 に示すような広範囲な領域を当する．「2), 3), 5)」における解析には，構造物対象とした流れ現象，かつ，沿岸部や都市部の局（物体）に作用する力の評価，流体－構造連成問所的領域における影響の把握や予測に関する防題，物体の移動（輸送）などの構造解析に関する災・環境のための大規模流体シミュレーションを知見が必要となるため，構造解析に関連する評価想定している．

や手法の妥当性について，参加研究者である京谷・寺田・加藤・高瀬（東北大），中畑（愛媛大），浅井（九州大），永井（岐阜大），車谷（茨城大）との連携をはかっている．また，「1)～5)」の連立一次方程式の並列化に関する高速化や妥当性については岩下（京都大），「6)」の効率的な可視化および VR 技術を用いた効果的な可視化の検討については小山田（京都大）との連携をはかっている．
参加者全員のメーリングリストを作成し，当該メーリングリストを用いて，研究の進捗における情報共有を行い，定期的に情報交換会を開催し，

図-1: 防災・環境のための大規模流体シミュレーション研究の推進および連携を進めている．
「2.(1)」に示した実施内容に関する研究テーマと(2) 共同研究分野

して，下記(1)～(4)の内容について報告する．(1)

本課題の共同研究分野は「超大規模数値計算系は「2.(1) 4)」，(2)は「2.(1) 5)」に関する研究応用分野」である．

テーマである．(3)は「2.(1) 2)」に関する研究テ
2

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

ーマであり，実スケール地形計算で用いた計算機（CRAY XE6）のアーキテクチャ(AMD Opteron 6000シリーズ)に関する並列効率について検討し，大規模計算を実施している．また，「2.(1) 6)」に関する可視化例が示されている．(4)は「2.(1) 3)」に関する研究テーマであり，AMD Opteron 6000 シリーズを用いた並列計算による計算例である．
(1) 汽水湖の塩水交換に関する 3 次元解析（池田・牛島）
沿岸に位置する汽水湖の塩水交換に関する 3 次元解析を行った．計算領域の水平方向の大きさは，
約 5.7km x 3.1km であり，一方汽水湖の水深は 5m

図-3: 計算開始から 494 秒後の塩分水平分布程度であるため，今回の計算では水平方向のみ領域分割を施して試計算を行った．並列数は，320である．
汽水湖の地形は，主要な位置の標高データを用意して，それを計算セル位置で空間内挿して定めた．地中内に位置する計算セルは，固体セルとして扱い，セル内に部分的に陸地部分を含む場合には，固液あるいは固気セルとして，混相場のモデルと同様の扱いとしている．
図-2 に，計算開始から 494 秒後の塩分鉛直分布を示す．汽水湖の初期塩分濃度は 0 としており，海域の底面付近から塩水がわずかに浸入し始めて

図-4: 計算開始から 494 秒後の塩分水平分布いることがわかる．
図-3 から図-6 に，水平方向の塩分濃度の時間的な変化を示す．汽水湖の開口部を通じて，塩水・淡水の交換が進行している状況が計算された．

図-5: 計算開始から 700 秒後の塩分水平分布

図-2: 計算開始から 494 秒後の塩分鉛直分布3学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

されている．図-11 は，各プロセスが保有する物体数の最大値 Npmax の時間的な変化を示したもの
で，動的領域分割法の方が平均して 73％ほど最大物体数が減少することが確認された．両手法の物
体運動に要した時間(動的負荷分散に伴う通信時間を含む)の合計は，動的領域分割法では 1355.4秒，静的領域分割法では 1674.7 秒となり，動的領域分割では約 80.9％計算時間が短縮された．
また，本計算手法が物体数 100 万を超える大規模計算にも適用できることを示すため，128 ノードで計算可能な最大物体数を求めたところ，物体数 3 億，流体セル数 42 億(物体数の約 14 倍と設定)

図-6: 計算開始から 1,660 秒後の塩分水平分布の計算が可能であることを確認した．計算結果の(2) 3 次元自由表面流れによる多数物体の大規模出力には並列 IO を利用しているが，このような大輸送解析（丸山・牛島）規模計算ではファイルへの出力時間が増大するた

3 次元自由水面流れにより，流体と同じ密度をめ，今後は大規模データの出力方法を検討する必有する多数の物体が輸送される計算を行った．流要があると考えられる．

体と物体間の力学的な相互作用は，圧力項と粘性項を体積積分する手法により精度良く評価されている．
並列計算は領域分割法に基づいて行われるが，流体計算を行うサブドメインと，物体の衝突および並進・回転運動を計算するサブドメインは独立に設定される．物体計算のサブドメインの大きさを計算中に動的に変化させる方法を「動的領域分割法」，また初期の領域分割のまま変化させない方法を「静的領域分割法」と呼ぶ．今回の大規模計

図-7: 初期の水面形と物体の分布算では，両手法の物体運動計算時間を比較するこ（黄色線は物体計算の領域を表す）とで動的負荷分散手法の有効性を確認した．
比較対象は水柱崩壊の影響を受ける 100 万の物
体輸送計算である．物体数は 100 万，流体セル数は物体数の約 14 倍にあたる 1,415 万 5,776 とし，流体には水と空気の物性値を設定した．また，１つの物体形状は 121 個の四面体で構成される回転楕円体とし，水と同密度とする． 計算プロセス数は 4,096 である．
4.8 秒間の計算結果を図-7 から図-10 に示す．
図-7 から図-9 に示された黄色の分割線は，動的領図-8: 2.4 秒後の水面形と物体の分布域分割法におけるサブドメインを表しており，物（黄色線は物体計算の領域を表す）

体数の偏りを解消するように計算中に領域が調整
4

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

方程式を，離散化手法としては SUPG/PSPG 法に基づく安定化有限要素法を用いた大規模津波遡上計算を行うことを目的とし，並列計算手法として図-12 に示す 3 つの手法を導入して並列化効率の観点からそれらの比較を行った.なお，使用したコンピュータは京都大学の CREY XE6 である．

図-9: 4.8 秒後の水面形と物体の分布（黄色線は物体計算の領域を表す）

図-12: 並列計算手法

図-13 に解析モデルを示す．計算条件として，総節点数は 6,738,732，総要素数は 37,560,556 と図-10: 4.8 秒後の計算結果の拡大図

し，最小メッシュ幅は地表面近傍で 0.50[m]とする．気体および液体は空気と水を仮定し，密度はそれぞれ 1.293[kg/m3]，1000[kg/m3]，粘性係数はそれぞれ 1.8×10-5[Pa・s]，1.0×10-3[Pa・s]とする．微小時間増分量は 0.05[s]とし，境界条件として壁面で Slip 条件を与える．

図-11: 各プロセスに含まれる物体数の最大値 Npmax の時間的な変化

図-13: 解析モデル(3) 実モデル地形における 3 次元気液二相流大規

図-14，15 に使用コア数とスピードアップおよ

模解析（樫山）び並列化効率の関係を示す．これらの図から，Flat

本研究では，基礎方程式として非圧縮性粘性流MPI， Hybrid-B，Hybrid-A の順に高い効率が得ら体を仮定した VOF 法に基づく 3 次元 Navier-Stokesれていることが確認できる．この原因の考察は十5

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

分でないが，プログラム中の境界条件処理を行う部分で OpenMP によるスレッド並列化が行えていないことや，ハイブリッド並列時にプロセスと CPUコアのバインドの仕方を指定していないことが大きな要因と考えられる．

図-16: 解析モデル

総 節 点 数 は 23,233,249 ， 総 要 素 数 は125,280,464 とし，最小メッシュ幅は地表面近傍で 0.25[m]，微小時間増分量は 0.05[s]，境界条件は Slip 条件を与える．並列数は 2048 とし，64 ノードを使用して計算を行う．
図-17 に可視化ソフト（Lumion）を用いて高品質な可視化を行った例を示す．津波が都市部に流入する様子が分かる．

図-14: 使用コア数とスピードアップの関係

図-17: より高品質な可視化（38 秒後）(4) 2 次元浅水長波流れ，3 次元気液二相流れにおける連成解析（松本）
2 次元非線形浅水長波方程式による広範囲な領域での計算を目的とした 2 次元（準 3 次元）解析，3 次 元 Navier-Stokes 方 程 式 に 基 づ い た

図-15: 使用コア数と並列化効率の関係

Phase-Field モデルによる局所的な領域での計算
以上の結果を踏まえ，Flat MPI を用いてさらなを目的とした 3 次元気液二相流解析の双方における大規模計算を行う．図-16 に解析モデルを示す．

る並列化を実施した．AMD Opteron 6000 シリーズを用いた並列計算の検証を行い，2 次元計算と 3次元計算の連成解析手法の検討を行った．図-18,19 に弱連成による計算結果を示す．
6

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月(a)Time = 0.0 s(a)Time = 0.0 s(b)Time = 0.2 s(b)Time = 0.3 s(c)Time = 0.4 s(c)Time = 0.6 s(d)Time = 0.6 s(d)Time = 0.9 s(e)Time = 0.8 s(e)Time = 1.2 s

図-18: ダムの崩壊問題【2D-3D 弱連成解析】

図-19: ダムの崩壊問題【2D-3D 弱連成解析】（3D から 2D への受け渡し，マニングの粗度係数=0.0）（2D から 3D への受け渡し，マニングの粗度係数=0.0）

7

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

図-18 の弱連成は，3 次元解析で得られた物理量(a)Time = 0.0 s(b)Time = 0.3 s(c)Time = 0.8 s(d)Time = 1.6 s(e)Time = 2.6 s（流速，気液界面）を 2 次元解析へ受け渡す結合
手法，図-19 の弱連成は，2 次元解析で得られた物
理量（平均流速，水位変動量）を 3 次元解析へ受け渡す結合手法により，時々刻々と変化する非定常計算を行った結果である．本結合手法は，連続の条件（流れの連続性）を満足するように 2 次元と 3 次元の物理量を受け渡している．緑と青の水面は 2 次元と 3 次元の計算を示す．
図-18 は，3 次元計算で求められる波が 2 次元計算で伝搬するかの検証を，図-19 は，2 次元計算で求められる波が 3 次元計算で伝搬するかの検証を目的として行った計算である．図-18 では，2 次元の計算領域において 3 次元計算からの波が，図-19
では，3 次元の計算領域において 2 次元計算からの波が伝搬している様子が確認できる．
図-20 は，2 次元解析で得られた物理量（平均流速，水位変動量）を 3 次元解析へ受け渡す結合手法を用いて，実際の広範囲な領域（沖合部）から局所的な領域（沿岸部や内陸部）を想定し，海底摩擦を考慮した 2 次元計算（沖合部）から 3 次元計算に移行した後に，傾斜のある底面（沿岸部），さらに構造物（内陸部）を設けた計算である．図-20 では，傾斜のある底面において伝搬されてきた波の上昇により，水際が一時的に沖合部にひき，その後，上昇した波が内陸部に伝わり，構造物に達している様子が確認できる．

4. これまでの進捗状況と今後の展望(1) 汽水湖の塩水交換に関する 3 次元解析
本中間報告では，比較的少ない並列数で，汽水湖と海域の間で行われる塩分交換の解析を試みた．
今後は，大規模計算により，より詳細かつ長時間にわたる現象の再現を行う予定である．
(2) 3 次元自由表面流れによる多数物体の大規模輸送解析
動的負荷分散手法を導入したことで，扱える計算対象の幅が広がったと考えられる．今後は実際の現象に合わせた計算を行い，漂流物の挙動や衝図-20: ダムの崩壊問題【2D-3D 弱連成解析】

突の影響について考察していきたい．大規模計算（傾斜と構造物を考慮，マニングの粗度係数=0.025）8

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月に着目した今後の展望としては，より効果的な動5. 研究成果リスト

的負荷分散手法の開発，計算結果のファイルへの(1) 学術論文書き出しの高速化があげられる．また現状，物体1) 牛島省, 福谷彰, 牧野統師 : 3 次元自由水面数が 1 億を超えるような場合に計算結果をそのま流中の接触を伴う任意形状物体運動に対するま可視化することが困難なため改善が必要である．

数値解法, 土木学会論文集, B, Vol.64, No.2,(3) 実モデル地形における 3 次元気液二相流大pp.128-138, 2008.
(2) 国際会議プロシーディングス

規模解析
現在は図-21 に示すような沖合からの 2 次元津

1) T. Fumuro, S. Tanaka and K. Kashiyama :波波動伝播解析と都市部の 3 次元自由表面流れ解

Large-Scale Tsunami Simulation Based on

析の連成解析が実現できていない．今後は，連成Three-Dimensional Parallel SUPG-VOF Method,解析の実現を行うとともに，都市部の解析領域を

Proc. of APCOM & ISCM 2013.
(3) 国際会議発表

より高領域にした大規模解析を行う予定である.

1) K. Kashiyama : Modeling and Simulation for
Safety and Environmental Problems Using
Virtual Reality Technique, APCOM & ISCM 2013
(Plenary Lecture).
2) T. Fumuro, S. Tanaka and K. Kashiyama :Large-Scale Tsunami Simulation Based onThree-Dimensional Parallel SUPG-VOF Method,
APCOM & ISCM 2013.
3) J. Matsumoto and K. Kashiyama : A Coupling
Method using FEM of 2D Shallow Water Flow and3D Gas-Liquid Two-Phase Flow, APCOM & ISCM
2013.
(4) 国内会議発表

図-21: 2D-3D 連成解析1) 丸山紀尚, 牛島省 : 粘性流れによる物体群輸(4) 2 次元浅水長波流れ，3 次元気液二相流れにお送の動的負荷分散に基づく並列計算法, 第 27

ける連成解析

回数値流体力学シンポジウム, 2013.

3 次元解析で得られた物理量を 2 次元解析へ受

2) 松本純一, 高田尚樹, 松本壮平 : 非構造格子

け渡す結合手法および 2 次元解析で得られた物理を用いた Phase-Field モデル自由表面問題に量を 3 次元解析へ受け渡す結合手法の検討を行っ

おける Allen-Cahn 方程式の計算法の検討, 第た．現在，上記(3)と月 1 回程度のミーティングを18 回計算工学講演会, 2013.

実施しながら 2 次元解析において非線形浅水長波3) 松本純一, 高田尚樹 : 非構造格子を用いた方程式から Boussinesq 方程式への拡張に着手しPhase-Field モデル自由表面問題における化ている．今後は，上記(3)と連携して，Boussinesq

学ポテンシャルとモビリティーの効果の検討,方程式への拡張および 2 次元計算と 3 次元計算の第 26 回計算力学講演会, 2013.

双方向での物理量の受け渡しが可能な結合手法を開発し，実験結果との比較を行う予定である．

9