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

14-NA16

防災・環境のための大規模流体シミュレーションとその可視化
樫山和男（中央大学）
防災・環境シミュレーションにおける津波氾濫解析として，東日本の海域（沖合部），沿岸部といった広範囲な領域から内陸部（都市部）の局所的領域を一体的に計算可能な大規模並列計算技術，環境流体解析として，沿岸部での格納容器の微細なクラックから流出するガスの流動における低マッハ数非等温圧縮性流体の 3 次元解析技術の開発を目的とする．また，津波氾濫解析，環境流体解析ともに実スケール地形での大規模並列計算を実施し，大規模データを有効利用した高効率かつ効果的な可視化の検討を行う．

1. 共同研究に関する情報

が近年高まっている．日本全体あるいは東日本や(1) 共同研究を実施した拠点名

西日本の大部分の領域を計算対象としたシミュレ

京都大学ーションにおいては，膨大な計算量を要する大規(2) 共同研究分野

模解析となり，スーパーコンピュータなどの超高(3)■超大規模数値計算系応用分野

速計算機の利用が不可欠となる．日本全体，東日

超大規模データ処理系応用分野

本，西日本などの広域な計算を行いながら，都市

超大容量ネットワーク技術分野

部など局所的領域での詳細計算を実施するうえで，

超大規模情報システム関連研究分野

広域では 2 次元（準 3 次元）解析，局所的領域では 3 次解析を用いた合理的な数値解法が開発でき

参加研究者の役割分担
参加研究者

れば，実用的な時間で計算結果を得るための非常樫山和男（中央大学：代表者）に有効な手段となり得る．また，沿岸部における牛島省（京都大学）

局所的領域での 3 次元詳細計算においては，防災

松本純一（（独）産業技術総合研究所）に関する問題とともに環境問題と密接な関係が生

浅井光輝（九州大学）

じる場合も考えられる．例えば，原子力使用済み

小山田耕二（京都大学）

燃料格納容器には，2 気圧程度の高圧状態で放射性岩下武史

物質が設置されており，格納容器内の高圧ヘリウ（北海道大学）

役割分担（担当分野等）

ムガスの漏洩は，環境に対する重大な問題を引き

津波氾濫解析技術開発と並列計算：

起こす可能性がある．

樫山，松本，浅井環境流体解析技術開発と並列計算：

一方，近年の可視化技術の進歩により，非常に牛島

写実的な可視化や立体視表現が可能になりつつあ小山田

る．近年の可視化技術と計算結果（数値データ）

連立一次方程式の高速解法に係る分野：岩下の有効利用により，数値シミュレーションで得ら

計算結果の可視化に係る分野：

れた現象を疑似的に体験するといった試みが行わ

2. 研究の目的と意義れている．

東日本大震災の大津波や日本でも発生し得る巨

本研究は，防災・環境シミュレーションにおけ

大台風などの大規模自然災害により，広範囲な領る津波氾濫解析として，東日本の海域（沖合部），

域を対象とした流れ現象，かつ，沿岸部や都市部

沿岸部といった広範囲な領域から内陸部（都市部）の局所的領域における影響の把握や予測などを目の局所的領域を一体的に計算可能な大規模並列計的とした防災・環境シミュレーションへの重要性

算技術，環境流体解析として，沿岸部での格納容1

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

器の微細なクラックから流出するガスの流動にお3. 当拠点公募型共同研究として実施した意義

ける低マッハ数非等温圧縮性流体の 3 次元解析技

本研究では「学際大規模情報基盤共同利用・共術の開発を目的とする．また，津波氾濫解析，環同研究拠点」のネットワーク型拠点のうち，「京都

境流体解析ともに実スケール地形での大規模並列大学」を共同研究先としており，計算機環境とし計算を実施し，大規模データを有効利用した高効て AMD Opteron 6000 シリーズのアーキテクチャ

率かつ効果的な可視化の検討を行う．

に基づく CRAY XE6 を利用した．研究の実施内容
は次の 6 つに大別される．

日本全体あるいは東日本や西日本の大部分の領域を計算対象とした広範囲なシミュレーションに
おいては，広域では 2 次元（準 3 次元）解析，局

1) 実スケール地形での Boussinesq 方程式による 2

所的領域では 3 次解析を用いた合理的な数値解法

次元浅水流れ大規模並列計算および 2D-3D 連の開発が必要であり，沿岸部においては，格納容成解析のためのモデリング

器の微細なクラックから流出するガス流動におけ2) 実スケール地形での Navier-Stokes 方程式による低マッハ数非等温圧縮性流体の 3 次元詳細計算

る 3 次元気液二相流れ大規模並列計算およびが必要となる．これらの実用的なシミュレーショ

2D-3D 連成解析のためのモデリング

ンを実現するには，大規模並列計算機および計算

3) 有 限 要 素 法 を 用 い た 実 ス ケ ー ル 地 形 で のの高速化が，また，得られた大規模かつ複雑な計

Boussinesq 方程式，Navier-Stokes 方程式（3 次

算結果を数値計算が専門でない一般の方もわかる

元気液二相流）による連成手法の検討および並ような分かり易い表現を実現するためには，効率列計算法の開発
4) 3 次元自由表面流れによる多数物体の大規模輸的な可視化技術と効果的な可視化が必要である．
以上に示した研究を目的とする本研究課題では，

送解析
5) 粒子法を用いた浅水長波方程式，Navier-Stokes

流体解析の研究者，計算機科学・可視化の研究者との体制および連携は不可欠である．

方程式（3 次元気液二相流）による連成手法の本防災解析は，近年，深刻かつ重大な問題とな検討および並列計算法の開発っている大津波や巨大台風による洪水などの大規

6) 沖合部，沿岸部，内陸部における効率的な可視模自然災害に関する研究課題であり，今後の国民化および VR 技術を用いた効果的な可視化の検の安全で安心な社会環境のための基盤研究として討意義がある．
本環境解析では，格納容器の微細なクラックか「1), 2), 6)」については主に樫山（中央大），「3)」ら流出するガスの流動を大規模並列計算環境で予については主に松本（産総研），「4)」については測し，外界への影響評価を行う．低マッハ数の非

主に牛島（京都大），「5)」については主に浅井（九等温圧縮性流体を計算する計算アルゴリズムは，

州大）が担当する．
「1)～5)」の連立一次方程式の当グループで開発・検証中であり，これをメニー

並列化に関する高速化や妥当性については岩下

コアシステムを基盤とする大規模計算環境で適切（北海道大），「6)」の効率的な可視化および VRに動作するものとし，従来困難であった環境影響技術を用いた効果的な可視化の検討については小評価を試みる．当該研究課題で得られた成果は，

山田（京都大）との連携をはかっている．

地熱発電や圧縮気体を利用する波力発電などの新参加者全員のメーリングリストを作成し，当該

エネルギー開発にも資すると考えられ，グリーン

メーリングリストを用いて，研究の進捗におけるイノベーションのための基盤研究として意義があ

情報共有を行い，定期的に情報交換会を開催し，る．

研究の推進および連携を進めている．
2

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

4. 前年度までに得られた研究成果の概要

5. 今年度の研究成果の詳細

前年度までに得られた研究成果のテーマ(1)～

本研究では，図-3 に示すような広範囲な領域を対(4)およびテーマ(2), (4)の計算結果（図-1, 2）を以象とした流れ現象，かつ，沿岸部や都市部の局所下に示す．

的領域における影響の把握や予測に関する防災・環境のための大規模流体シミュレーションを想定(1) 汽水湖の塩水交換に関する 3 次元解析

している．

(2) 3 次元自由表面流れによる多数物体の大規模輸
送解析(3) 実モデル地形における 3 次元気液二相流大規
模解析(4) 2 次元浅水長波流れ，3 次元気液二相流れにおける連成解析

図-3: 防災・環境のための大規模流体シミュレーション「3.」に示した実施内容に関する研究テーマとして，下記(1)～(4)の内容について報告する．(1)は「3.
3)」，(2)は「3. 4)」，(3)は「3. 5)」，(4)は「3. 6)」に

図-1: 3 次元自由表面流れによる関する研究テーマである．

多数物体の大規模輸送解析(1) Boussinesq 方程式，Navier-Stokes 方程式（3 次Green surface

元気液二相流）による連成手法の検討および

2D Shallow Water Flow

並列計算法の開発（松本）Blue surface(1)-1 連成手法の検討

3D Two-Phase Flow

昨年度は，2 次元非線形浅水長波流れ解析，3 次
元気液二相流れ解析について，3 次元解析で得ら
図-2: 2 次元浅水長波流れ，3 次元気液二相流れれた物理量（流速，気液界面）を 2 次元解析へ受における連成解析

け渡す結合手法，および，2 次元解析で得られた物理量（平均流速，水位変動量）を 3 次元解析へ

テーマ(2)は，今年度課題の実施内容「5. (2)」へ拡受け渡す結合手法を開発し並列計算を行った．今張されている．テーマ(4)は，今年度の実施内容「5.

年度は，2 次元から 3 次元，3 次元から 2 次元への(1)」
（実スケール地形への計算を目的とした 2D-3D

双方向の物理量を同時に受け渡す連成手法の開発

連成解析）へと展開され，前年度の採択課題では，および並列計算を実施した．連続の条件（流れの参加した研究機関間にて連名での成果発信が行わ連続性）を満足するように 2 次元から 3 次元およれており，当該共同研究を実施したことによる有び 3 次元から 2 次元の物理量を受け渡している．

意義な研究が進められている．
3

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

AMD Opteron 6000 シリーズを用いた並列計算の検証を行い，開発した連成手法の検討を行った．
図-4 に弱連成による計算結果を示す．緑と青の水
面は 2 次元と 3 次元の計算を示す．
計算モデルは，昨年度実施した実際の広範囲な領域（沖合部）から局所的な領域（沿岸部や内陸部）を想定し，海底摩擦を考慮した 2 次元計算（沖合部）から 3 次元計算に移行した後に，傾斜のある底面（沿岸部），さらに構造物（内陸部）を設けたものを用いた．
今回の計算は，2 次元から 3 次元と 3 次元から 2次元への双方向の物理量を同時に受け渡す連成手
法の検証を行うために，水深を浅く設定し，押しよせた波が傾斜を遡った後に，沖合に戻って行くといった計算である．図-4 では，押しよせた波が傾斜を遡った後に，波が沖合に戻って行くといった押し波，引き波の様子が確認できる．

図-4: 双方向の物理量を同時に受け渡す連成解析(1)-2 Boussinesq 方程式「(1)-1」の連成解析は，分散性が考慮されていない 2 次元非線形浅水長波方程式を用いている．
ここでは，分散性を考慮した浅水長波流れの計算が可能な次式に示す Boussinesq 方程式を用いた計算手法の開発を行う．
運動方程式

(1)
連続式
(2)
式(1)の青色の項が分散項である．これらの分散項
は空間の二階導関数，かつ，時間の一階導関数からなっており，非線形浅水長波方程式に比べて時間方向の離散化において，陽解法を採用するのは4

学際大規模情報基盤共同利用・共同研究拠点 平成 26 年度共同研究 最終報告書 2015 年 5 月(1)-3-1 流体の基礎方程式

難しい方程式である．本研究では，非構造格子に

非圧縮粘性流れにおける Navier-Stokes 方程式のおける陰解法に基づいた大規模並列計算手法の開発を行ってきているため，これまでの開発手法を運動方程式と連続式は次式によって表される．

拡張して Boussinesq 方程式においても並列計算が可能である．

(3)

検証問題として孤立波の伝搬解析を行う．図-5(4)に初期条件を，図-6 に非線形浅水長波方程式と， ， は流速，圧力，外力であり， は粘性係

Boussinesq 方程式の計算結果を示す．

数である．
時間方向の離散化には，安定性に優れ時間増分を大きくとれる陰的解法を適用し，連続式は完全陰的に扱い，計算効率の良い分離型解法を用いた．
Navier-Stokes 方程式における空間方向の離散化には，数値安定性を考慮した MINI 要素（直交基底
気泡関数要素安定化法）を採用した．
(1)-3-2 気液界面の基礎方程式

図-5: 孤立波の伝搬解析における初期条件

Phase-Field モデルによる自由表面流れの界面位置を決定する界面関数は，流速を ，界面関数をとして次式の Cahn-Hilliard 方程式により表される．
(5)界面関数は，液体で，気体での値をとり，

を気液界面と定義する関数である．
とした場合，各節点上での密度 ，粘性係
数 は，界面関数 と液体，気体の密度 ,(a) 非線形浅水長波方程式

係数

,

，粘性

を用いて式(6)によって計算され，式(3),(4)の Navier-Stokes 方程式にわたされる．
(6)
(7)
式(5)の

はモビリティー（易動度）である．

は化学ポテンシャルを示し，式(8)よって計算される．

(b) Boussinesq 方程式(8)

図-6: 孤立波の伝搬解析における計算結果

非線形浅水長波方程式の計算結果は，波の非線は曲率である．Phase-Field モデルで用いられる界形性が表れ進行方向に対して波がたってしまうよ

面関数 は液相と気相の界面近傍で連続的に変化

うな形状になっている．一方，分散性を考慮した

し により界面上の滑らかさを制御する．式(8)のBoussinesq 方程式の計算結果は，波の進行方向に

化学ポテンシャルにおいて，対して波は立たず，後方に波の分散性からくる小される．

は次式よって表さな波形も確認できる．
(9)(1)-3 Phase-Field モデルを用いた気液二相流解析
内陸部の計算は，Phase-Field モデルを用いた 3(10)

次元気液二相流れによる並列計算手法を開発する．
5

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

並列効率100％の場合
400x400x240分割を規格化するために導入したパラメータであり，計算速度（倍）は8と （要素毎に定義）は次式のように

表される．
(12)
, ,

,

6

4

2

はそれぞれ界面モビリティー，界面エ

ネルギー，連続的に変化する気液界面幅（1 分割の場合は 1.0, 2 分割の場合は 2.0），要素毎に定義されるメッシュ分割幅である．実際の計算では

0

8

16

24

32

ノード数

,

を指定する．

図–7: ノード数と計算速度

時間方向の離散化に陰的解法を採用し，非線形方程式の解法に Newton-Raphson 法に基づいた反表–1，図–7 は 4 ノードの計算速度を 1 倍とした場

復計算を用いた．Cahn-Hilliard 方程式における空合の各ノートにおける計算速度を示している．ど間方向の離散化には，数値安定性を考慮した MINIのノードの計算効率もほぼ 100%となっている．特要素（直交基底気泡関数要素安定化法）を用いたに，4 ノードの 8 倍となる 32 ノードでは，計算速た．

度が 8.58 倍（計算効率 107.23%）となっており，(1)-4 並列効率の考察

計算効率が 100%を超えている．

これまでに開発した並列計算手法に基づき，乱領域分割法を用いた並列計算では，分割された流モデルを導入した気液二相流解析に拡張した．

各領域どうしで，データの送受信をおこなうため，並列効率の評価のため，以前に実施した一億自由

一般的には，ノード数の増加に伴い並列効率は下度解析の計算モデル（ミルククラウン）を用いてがるが，表–1，図–7 は，それとは異なった結果と検討を行った．本解析は，200GB 程度のメモリ容なっている．この原因を考察するため，非定常拡量を必要とする．京都大学の計算機（CRAY XE6）

散方程式における有限要素解析を行った．空間方は 1 ノード当たりのメモリ容量が 64GB のため，

向には四面体一次要素，時間方向には 4 次精度ル200GB 以上のメモリ容量が確保可能な 4 ノードか

ンゲクッタ法による陽解法を採用した．非定常拡

らの計算を実施した．

散方程式における陽解法は，省メモリのため１コ

表–1，図–7 に 4～32 ノードまでの計算時間（時アからの計算が可能となる．また，本研究で用い間方向に 3 ステップの計算を行い，3 回計測したている陰解法に基づく分散メモリ型の並列化は，値の平均値）と並列効率の結果を示す．Flat MPI

主に行列・ベクトル積の部分に関する並列計算がによる領域分割法を用いたプログラムの並列計算支配的であるため，非定常拡散方程式における陽を実行した．

解法の並列性能の結果と類似の傾向となることが予想される．
表–1: 計算時間の並列効率

表–2，図–8 に 1～128 ノードまで（ノード間）

ノード数
計算時間（秒） 計算速度（倍） 計算効率（％）
4( 128コア)
428.1801.00
100.00214.369
1.99
99.87
8( 256コア)
105.9194.04
101.0616( 512コア)
49.916
8.58
107.2332(1024コア)の計算時間（時間方向に 100 ステップの計算を行い，3 回計測した値の平均値）と並列効率の結果を示す．Flat MPI による領域分割法を用いたプログラムの並列計算を実行した．
6

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

表–2: 計算時間の並列効率

表–3: 計算時間の並列効率
コア数
1
2
4
8
16
32

ノード数
計算時間（秒） 計算速度（倍） 計算効率（％）
1( 32コア)
145.4281.00
100.0073.892
1.97
98.41
2( 64コア)
36.763
3.96
98.90
4( 128コア)
17.9918.08
101.048( 256コア)
8.853
16.43
102.6716( 512コア)
4.456
32.64
101.9932(1024コア)
2.290
63.51
99.23
64(2048コア)
1.174
123.87
96.78
128(4096コア)

32

並列効率100％の場合
300x300x300分割
計算速度（倍）128

並列効率100％の場合
300x300x300分割
計算速度（倍）計算時間（秒） 計算速度（倍） 計算効率（％）
1589.830
1.00
100.001241.912
1.28
64.01
771.708
2.06
51.50
686.841
2.31
28.93
290.177
5.48
34.24
145.428
10.93
34.16

96

64

24

16

8

32

0

8

16

24

32

コア数
0

32

64

96

図–9: コア数と計算速度

128

ノード数

ており，その場合の 1 コア当たりの節点数はそれ

図–8: ノード数と計算速度

ぞれ 58877，29767 であった．このことから 1 コア

表–2，図–8 をみると，どのノードの計算効率も当たりの計算規模が 3 万節点から 6 万節点の間で100%に近い値となっている．本計算結果では，8計算速度が向上していることが考えられ，その場

ノード，16 ノード，32 ノードの計算において計算

合には，1 ノード内の計算速度が（1 コアからみて効率が 100%を超えている．表–3，図–9 に 1～32

32 コアで約 11 倍から 12 倍弱に）向上するためでコアまで（ノード内）の計算時間（時間方向に 100はないかと考察される．

ステップの計算を行い，3 回計測した値の平均値）(1)-5 実スケール地形による 2D-3D 連成解析と並列効率の結果を示す．Flat MPI による領域分

実スケール地形での計算例として，女川港をモ

割法を用いたプログラムの並列計算を実行した．

デル化した連成解を行う．3 次元領域は内陸部を1 コアを 1 倍とした場合に，32 コアで約 11 倍の実地形データにより，モデリングし 400m×250m

計算速度となっている．CRAY XE6 で採用されてで 任 意 分 割 を 行 っ た 節 点 数 3,957,882 要 素 数

いる CPU のアーキテクチャ(AMD Opteron 6000 シ

23,155,890 のメッシュを，2 次元領域は 200m×リーズ)は，1 モジュールの中に 2 つの整数演算コ250m で節点数 59,605 要素数 118,236 のメッシュをアを有し，1 ノード内で 16 モジュール（32 コア）用いた．図-10 は 2 次元と 3 次元の接合条件である．

となっている．本計算は主に実数演算であるので，
1 コアからみた 32 コアの計算速度が半分以下になっているのは，そのためであると推測される．
気液二相流解析では，32 ノード（1024 コア）の際に，計算効率が 107.23%となっており，その場
合の 1 コア当たりの節点数は 41967 である．また，
非定常拡散方程式での計算結果では，16 ノード，32 ノードで計算効率が 102.67%，101.99%となっ

図-10: 2 次元と 3 次元の接合条件

7

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

2 次元計算と 3 次元計算を連成した解析結果を図流れにより，流体と異なる密度を有する多数の物-11 に示す．図-11 を見ると，2 次元から 3 次元に体が，静止構造物との衝突を伴いながら輸送され水面（波）が伝搬し，3 次元計算の内陸部へ向かる計算を行った．流体と物体間の力学的な相互作って水面が進んで行く様子が確認できる．

用は，圧力項と粘性項を体積積分する手法により精度よく評価されている．
並列計算は Flat MPI を用いて，領域分割法に基づいて行われる．また，流体計算を行うサブドメ
インと，物体の衝突および並進・回転運動を計算するサブドメインは独立に設定されている．物体
計算において，大規模分散メモリ並列環境で効率的な剛体運動の計算を行うため，剛体のデータ構造を分離している．
計算対象は水柱崩壊の影響を受ける 18,032 個の物体と，20 本の静止した円柱群との衝突・輸送計算である．物体数は合計 18,212，流体セル数は
32,768,000 とし，流体には水と空気の物性値を設定した．また，本計算にはテトラポット，円柱構
成物体，回転楕円体，車の模型の４種類の物体を用いており，構成四面体数はそれぞれ 353，1,320，
121，1,382，密度はそれぞれ 1.1×103 [kg/m3]，2.4×103[kg/m3]，1.1×103 [kg/m3]，0.1×103 [kg/m3]である．
計算プロセス数は 4,096 である．
実時間 2.525 秒間の計算結果を図-12(a)-(j)に示す．時間が進行するにしたがって，自由水面流れによって物体群が相互に衝突しながら輸送される様子が確認できる．ウィークスケーリングに基づいた並列化効率は，並列数 512 を基準として並列
数 4,096 で約 0.91 であり，高い効率が得られていることを確認した．

図-11: 実スケール地形による 2D-3D 連成解析の結果(2) 3 次元自由表面流れによる多数物体の大規模(a) t = 0[s]輸送解析（牛島）
３次元多相場モデルを用いて，３次元自由水面
8

学際大規模情報基盤共同利用・共同研究拠点 平成 26 年度共同研究 最終報告書 2015 年 5 月(b) t = 0.2 [s](g) t = 1.5 [s](c) t = 0.45 [s](h) t = 1.8 [s](d) t = 0.7 [s](i) t = 2.0 [s](f) t = 1.0 [s](j) t = 2.525 [s]図-12: 3 次元多相場モデルを用いた計算結果

9学際大規模情報基盤共同利用・共同研究拠点 平成 26 年度共同研究 最終報告書 2015 年 5 月(3) 粒子法を用いた流体剛性連成解析（浅井）

0.48m

8.3m
8.3m
1.54m

0.8m

0.49m

SPH (Smoothed Particle Hydrodynamics) 法 や

0.69m

0.2m

MPS2)法などの粒子法は，完全なメッシュレス法(a)ポストテンション式(b)プレテンション式

図-14: 歌津大橋上部構造断面図

であり，自由表面の変動の激しい流体問題だけでなく破壊を考慮した固体等での利用が注目されて8.0m子を用いたまま剛体までを扱うことができ，流体

0.8m

1.5m

4.8m

0.8m

(a)ポストテンション式モデルている．本研究では，固体間の接触を含む問題に

0.8m

6.0m

0.2m

0.8m

剛体連成解析手法の応用例がすでに複数提案され

8.0m

0.5m

0.5m

いる．さらに，DEM の概念を導入すれば，同じ粒(b)プレテンション式モデル

図-15: 橋梁モデル断面図

おいては，流体はナビエ・ストークス方程式，剛体は運動方程式を解き，両者を境界面での力によ

表-4: 解析条件

り連成する方法が妥当とされる定式化を採用し，

粒子
間隔
0.2m

著者らがこれまでに提案してきた高精度の境界処理法を組み合わせることで，より適切な境界条件

総粒子
数
約 106
万

時間増
分
0.001sec

実時間

解析時間

10sec

約 10 時間を有した流体剛体連成解析手法へと発展させる．
その後，この手法を東日本大震災による宮城県南歌津大橋は今回の震災により，第 3 径間から第三陸町の橋梁の流失被害の再現解析へと適用した．

10 径間の計 8 径間が流失した(図-13)．震災後の現さらに，同地区での津波遡上解析を行い，流体剛

地調査により，計 8 径間において橋梁上部構造の体連成解析を行う際の流入条件への適用方法を考

一部は裏返しに流失しており，第 3 径間から第 7察した．

径間にかけては転覆せずにそのまま押し流される
ように流失していた．なお，裏返しに流失してい(3)-1 橋梁上部構造の流失シミュレーション

たのは桁高の大きいポストテンション式の上部構(3)-1-1 解析モデル

造，水平に流失した第 3 から第 7 径間の上部構造「(3)」で述べた流体剛体連成解析手法を用いて，
図-13 に示す宮城県南三陸町歌津大橋を模した橋

はすべて桁高の小さいプレテンション式となっていることに注目した．本研究では，この 2 種の上梁上部構造の流失シミュレーションを行う．

部構造の形状より，桁高と総幅員，地覆の張り出し長を簡易的に再現し(図-14,15)，上部構造の断面形状の違いによる流失挙動の変化を見ることとした．

プレテンション式(G3-G7)

図-16 に解析モデルの概要を示す．津波を模擬する水の形状は非段波状の波を簡易的に再現した．
津波高は 8m を想定し，流入条件として 10m/secの初速を水に与え，貯水部の先端位置を通過する

ポストテンション式(G8-G10)

水に継続的に 10m/sec を与えることで津波高を保ったまま模型に衝突させることとした．

図-13: 歌津大橋の橋梁流失状況

10

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

津波衝突時に鉛直方向に受ける流体力が大きくな

＜平面図＞
11m

1.0m
1.0m

75.0m

が遅れることにより，より回転するようなモーメ
ントが作用するためと考えられる．ポストテンシ＜側面図＞

ョン式上部構造を模した橋梁モデルがより回転す

10m/s
30.0m

るといった現象は，歌津大橋の実際の被害報告と

8.0m

20.0m

ること，また桁高の高さにより橋脚と床版の接触

10.0m

6.7m
25.0m

類似したものとなった．本解析では歌津大橋上部構造の断面形状を正確に再現したものではなく，

8.0m

支承やサイドブロック，橋脚と上部構造間の摩擦

28.0m

20.0m

20.0m

力，また津波の流入条件も実際のものとは一致し

図-16: 解析モデル概要

ないため，定性的な評価にとどまるが，紹介した
流体剛体連成解析により，断面形状の違いが異なTime=1.70s

Time=2.75s

る流失挙動を与えるといった著者らの予想を裏付ける結果となった．ただし，まだ支承部のモデル化が不十分なことに加えて，流入条件が実際の被害とはことなるため，この結果をもって妥当性確(a)ポストテンション式モデル
Time=1.20s

認が実施できたとは断言はできない．

Time=2.75s

解析範囲
歌津大橋

(b)プレテンション式モデル

図-17: 各モデルの流出挙動の比較

モデルの後方には後流側の壁面に衝突した波が反流入領域射して戻ってこないようにタンクを設けている．
橋脚は完全に固定し，水の密度は 1000kg/m3 上部構造の密度は鉄筋コンクリートを想定し
津波領域

2450kg/m3 とした．上部構造と橋脚のポアソン比は鉄筋コンクリートの一般的な値 0.28 を採用し，弾性定数は 200N/m2 を用いた．解析条件の詳細は表

図-18: 歌津地区解析領域および津波流入領域-4 に示す．
表-5: 解析条件(3)-1-2 解析結果

粒子
間隔
3m

1.20s~2.75s 時の解析結果を図-17 に示す．なお
上部構造は緑，橋脚は灰，水は青で色付けしている．図-17 よりプレテンション式上部構造モデルと比べポストテンション式上部構造モデルはより回転を伴う流出挙動が確認された．これは，ポストテンション式上部構造は地覆張り出し長が大きく，
11

総粒子
数
約 330
万

時間増
分
0.005sec

実時間

解析時間

250sec

約 42 時間

学際大規模情報基盤共同利用・共同研究拠点 平成 26 年度共同研究 最終報告書 2015 年 5 月(3)-2 津波遡上解析ため，同図上向きの流れが支配的となり，左方向

前章の結論を受けて，本章では橋梁流失シミュへの浸水が若干遅れているが，実際の被害状況とレーションにおける津波の流入条件を見直すため

おおむね一致しており，地形モデルについては妥に，同地区での津波遡上解析を行った事例を報告

当なモデルが作成できたと考えている．

する．
Time=15s

Time=15s

Time=120s

Time=110s

Time=135s

Time=135s(3)-2-1 解析モデル
解析の対象地は，前章で橋梁流失シミュレーシ
ョンのモデルとして用いた歌津大橋が位置する宮城県南三陸町歌津地区と伊里前湾とした．図-18に解析領域と歌津大橋の位置関係を示す．ここでは予め，解析範囲のほかに津波領域を縦約 1km×
横約 900m，津波高 10m，水深 5m で用意しておき，流入領域を通過する水に継続的に 10m/s を与える
条件としている．今後，浅水長波方程式を解くことで，震源から津波伝搬現象を再現したデータをそのまま入力波として用いる計画であるが，今回は地形モデルの実際よりも短時間かつ速い流れを入力することで地形モデルの再現性についてなどa)G3～G7周辺領域b)G8～G10周辺領域を考察するに留める．また解析範囲の周囲に壁を[m/s]

設け，津波流入領域以外からの解析範囲への流入，
図-20: 津波衝突時の橋梁周りの流速分布

流出が無いようにしている．解析条件の詳細は表-5 に示す．

次に，図-20 に歌津大橋の位置に津波が到達した
際の，橋梁周りの水の流速を示す．ここでも，解析結果上に灰色で表示した橋梁モデルを表示している．なお，今回の解析結果において上部構造の高さに到達した時点での，上部構造に対して津波
流入方向の流速を評価すると，3~4m/s 程度となっており，「(3)-1」で示した解析における流入条件と図-19: 浸水被害報告と解析結果の比較

比べて小さな値をしめした．また回転を伴う流失が見られたポストテンション式上部構造の周囲は，(3)-2-2 解析結果
図-19 は歌津地区における津波の浸水範囲と本

プレテンション式上部構造周囲とくらべ，解析時

解析結果を比較したものを示す．なお浸水被害報

間全体を通して流速が大きくなり，流速の最大値

告図において灰色が浸水領域を示し，解析結果では約 2 倍の違いが現れた．このことからも再度，

は水は青，地形はオレンジで色付けしている．ま適切な流入条件のもと，この津波遡上解析を実施

た，解析モデル作成の際に橋梁はモデル化していし直し，橋梁モデルに作用させる際の流入条件をないが，解析結果上に緑色で表示することで橋梁再検討する必要がある．

の場所を示すことにした．前述のように実際よりは短時間かつ流れの速い津波モデルを流入させた
12

学際大規模情報基盤共同利用・共同研究拠点 平成 26 年度共同研究 最終報告書 2015 年 5 月(4) 津波シミュレーションの可視化（樫山）

凹凸を表現する指数および気泡量の調整等を行う

東北地方太平洋沖地震津波を契機に，我が国でことで図-22 に示す可視化映像を作成した．

は津波に対して，災害規模に応じた対応を取ることとなった．すなわち，設計外力以下の津波に対しては，対策構造物で浸水等を防ぎ，設計超過津波に対しては，早期の警報発令や避難誘導といった情報による対処で人的被害を防ぐことを目指すというものである．しかしながら，ハザードマッ
プを示すだけでは，地域住民が事前に襲来する津波をイメージすることは難しく，避難情報としては不十分と思われる．そこで，本研究では，地域防災教育の一助となり得る可視化について検討した．

図-22: 可視化例

本解析では，検討地域に対して，差分法を用いた VOF 法に基づく CADMAS-SURF/3D1)による津

6. 今年度の進捗状況と今後の展望

波浸水解析を実施した．計算格子間隔は水平方向(1) Boussinesq 方程式，Navier-Stokes 方程式（3 次は 5m で一定としたが，鉛直方向は津波の遡上へ

元気液二相流）による連成手法の検討および並の影響が大きいので，解像度を上げて 1.0m とした．

列計算法の開発とその可視化解析に用いた建物レベルデータをそのまま可視

中間報告では，2 次元非線形浅水長波流れ解析，化に使用した場合，建物形状は，図-21 右上に示す

3 次元気液二相流れ解析について，2 次元から 3 次

通り，側面が台形形状となる．ここを同図右下に

元と 3 次元から 2 次元への双方向の物理量を同時示すように矩形へ修正して陸上標高データと同化に受け渡す連成手法の開発および並列計算を実施させた後，航空写真画像をテクスチャーとして貼した．分散性を考慮した浅水長波流れの計算が可り付けた．

能な Boussinesq 方程式を用いた計算手法の開発を実施した． Phase-Field モデルを用いた 3 次元気液

全領域

二相流れによる並列計算手法を開発し，京都大学

元の計算格子の計算機（CRAY XE6）による並列効率の検討を行った．最終報告では，「5. (4)」と協力して，2 次
x

元浅水流れと実スケール地形での 3 次元気液二相

新節点により作成したモデル

流れにおける 2D-3D 連成解析のためのモデリングzyおよび並列計算を行った．
(2) 3 次元自由表面流れによる多数物体の大規模輸

z

x
x
y

送解析

y

本中間報告では，流体と密度の異なる多数の物

図-21: 建物形状モデルの修正例

体が，静止構造物との衝突を伴いながら，自由表
津波の水表面は，CADMAS-SURF/3D における

面流れにより輸送される解析を試みた．今後は，VOF 関数計算値 0.5 の位置を水表面とし，現実感

実際の市街地を模して地形や静止構造物を計算領の向上を目的として水表面のポリゴン勾配に応じ域内に設定し，多数の任意形状物体の輸送解析をて白波が立った水流のＣＧ効果の付加や，水面の行いたい．また，流体や輸送される物体が静止構
13

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

造物に与える流体力，衝突力の影響について考察

SPH 法による流体剛性連成解析とその橋梁流していきたい．

失被害予測への応用, 土木学会論文集 A2（応用力学）, 70, No.2, p.I_317-I_328, 2014.

(3) 粒子法を用いた流体剛性連成解析4) 浅井光輝, 合田哲朗, 小國健二, 磯部大吾郎,

橋梁上部構造の流失被害を再現することで数値

樫山和男, 一色正晴 : 安定化 ISPH 法を用いた解析手法の実スケール問題の妥当性を確認するこ

津波避難ビルに作用する流体力評価, 土木学とを目的として，粒子法を用いた流体剛体連成解会論文集 A2（応用力学）, 70, No.2, p.I_649-I_658,

析手法を提案した．その後に，東日本大震災時に2014.

流失被害が発生した南三陸町にある歌津大橋をと

5) 森本敏弘, 浅井光輝, 笠間清伸, 藤澤和謙, 井

りあげ，被害調査報告を基に構造形式と照合させ

元佑介 : 安定化 ISPH 法による拡張ダルシー

ることで被害形態を 2 種類に分類したうえで，こ則とナビエ・ストークス方程式の統一解法, 土れらを簡易にモデル化した橋梁上部構造モデルを

木 学 会 論 文 集 A2 （ 応 用 力 学 ） , 70, No.2,

用いて，橋梁の流失挙動予測を行った．桁高の大p.I_213-I_221, 2014.

きいポストテンション形式の上部構造はプレテン

6) A. M. Aly and M. Asai : Modelling of non-Darcy

ション形式の上部構造と比べて回転量が大きくな

Flows through porous media using an extended

る傾向が確認でき，提案した流体剛体連成解析手

Incompressible Smoothed Particle Hydrodynamics,法の妥当性の一端は示すことができた．ただ，まNumerical Heat Transfer, Part B:Fundamental, 67,

だ完全に転覆現象を再現するには至っておらず，

pp. 255-279, 2014.

今後は橋梁に作用させる津波の流入条件の見直し， 7) 浅井光輝, 市村強, M.L.L.Wijerathne, 堀宗朗 :サイドブロックや支承による抵抗力のモデル化を

スパコンを用いた自然災害被害予測シミュレ

実施する計画である．最後に，津波の流入条件のーションの現状, 第 7 回土砂災害に関するシン見直しのために実施している同地区における津波

ポジウム論文集【招待論文】, pp.199-204, 2014.

遡上解析例を紹介した．遡上解析モデルの作成は

8) N. Takada, J. Matsumoto and S. Matsumoto : A

無事に実施できていることが確認でき，既に実施

Diffuse-interface

Tracking

Methodforthe

した震源からの津波伝搬解析結果とこの遡上解析

Numerical Simulation of Motions of a Two-phaseとのシームレスに連動させることで橋梁周辺の流

Fluid on a Solid Surface, The Journal of

れをより詳細に評価できるものと考える．

Computational Multiphase Flows, 6(3), pp.283-298,
2014.

7.

9) 高橋佑典, 桜庭雅明, 樫山和男 : 津波シミュレ

研究成果リスト(1) 学術論文

ーションのための CIVA-安定化有限要素法, 土1) 田邊将一, 浅井光輝, 中尾尚史, 伊津野和行 : 3

木学会論文集 A2（応用力学）, 70, No. 2, p.
I_349-I_356, 2014.

次元粒子法による橋桁に作用する津波外力評
価 と そ の 精 度 検 証 , 構 造 工 学 論 文 集 , 60A,

10) 岩塚雄大, 古牧大樹, 西畑剛, 川辺赳史, 樫山

293-302, 2014.

和男 : 地域防災教育のための 3 次元津波浸水

2) A. M. Aly and M. Asai : Incompressible Smoothed
Particle

Hydrodynamics

Simulations

解析とその可視化に関する研究, 土木学会論of

文集 F3（土木情報学）, 70(2), I_152-I_159, 2014.

Fluid-Structure Interaction on Free Surface Flows,11) 川辺赳史, 樫山和男, 宮地英生, 岩塚雄大, 古International Journal of Fluid Mechanics Research,

牧大樹, 西畑剛 : 可聴化技術を用いた津波疑41, pp.471-484, 2015.

似体験システムの構築, 土木学会論文集 F3（土木情報学）, 70(2), I_235-I_242, 2014.

3) 田邊将一, 浅井光輝, 宮川欣也, 一色正晴 :14

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

5) S. Tanabe, M. Asai, K. Terada, K. Kashiyama, S.

(2) 国際会議発表
1) J. Matsumoto and K. Kashiyama : An Interaction

Moriguchi and M. Kurumatani : Numerical

Method using Implicit FEM of Shallow Water

Evaluation of Impact Tsunami Force Acted on

Flow

Bridge Girders by Using a Particle Method,and

Gas-Liquid

Two-Phase

Flow,

COMPSAFE, 2014.

COMPSAFE, 2014.

2) J. Matsumoto and K. Kashiyama : A coupling(3) 国内会議発表

Method of Free Surface Flow using FEM for
Boussinesq

Equationsand1) 松本純一, 樫山和男 : 有限要素法を用いた 2D

Navier-Stokes

浅水流れと 3D 気液二相流れにおける結合法のEquations, WCCM XI, 2014.

検討, 第 19 回計算工学講演会, 2014.

3) T. Morimoto, M. Asai and K. Kasama: Fundamental

2) 岩塚雄大, 古牧大樹, 西畑剛, 川辺赳史, 樫山Study for Seawall Collapse Simulation during

和男 ： 地域防災教育のための３次元津波浸水

Tsunami by Using a Particle Method, COMPSAFE,

解析とその可視化に関する研究, 土木学会土

2014

木情報利用技術講演集, 39, pp.181-184, 2014.

4) A. M. Aly and M. Asai : Simulation of3) 青木一真, 牛島省, 鳥生大祐 : 自由水面流れに

Fluid-Structure Interaction in Flood Disastersより輸送される多数の浮遊物体と静止構造物

Using ISPH Method, pp.590-593, COMPSAFE,との衝突の数値計算, 第 28 回数値流体力学シ2014.

ンポジウム, A06-4, 2014.

15