名古屋大学 不老ユーザ会，8/30, 2021

災害被害予測に向けた粒子法によるマルチフィジックスシミュレーション
九州大学
浅井 光輝,モリカワ ダニエル（D3）,辻 勲平（D2）
Structural Analysis Lab.

1

2011 Tohoku Tsunami

2017 heavy rain disaster in Kyushu

2018 heavy rain disast2

MScPHy(Multi-Scale and PHysics)ﾏﾙﾁｽｹｰﾙ粒子法
スケール間をシームレスにつなぐブリッジング技術は構築済み

レベル0解析
解像度:10m～100m

レベル1解析

レベル2解析

解像度:1m～2m

従来の津波解析（津波伝搬解析）トップダウン型マルチスケール
津波遡上現象

解像度:10cm～20cm

橋梁等の土木構造物の被害予測（橋梁の流失被害）

都市部津波遡上解析の立体視例（防災教育等を通してソフト防災に貢献）

Structural Analysis Lab. ,Kyushu University

3

MScPHy(Multi-Scale and PHysics)ﾏﾙﾁﾌｨｼﾞｯｸｽ粒子法
複合災害予測し、ハード防災に貢献

流体力を受ける構造 流体-剛体連成SPH法(数億粒子まで対応可)

構造
洗掘＋液状化による構造物の不安定化
剥落？

地盤水理
胸壁の破壊

浸透流・洗掘解析用SPH法
基礎地盤の洗掘 (実験との比較で精度確認中)
Structural Analysis Lab. ,Kyushu University

衝撃的な波力MScPHy(Multi-Scale and PHysics)ﾏﾙﾁｽｹｰﾙ粒子法
スケール間をシームレスにつなぐブリッジング技術は構築済み

レベル0解析
解像度:10m～100m

レベル1解析

レベル2解析

解像度:1m～2m

従来の津波解析（津波伝搬解析）

解像度:10cm～20cm

津波遡上現象

橋梁等の土木構造物の被害予測（橋梁の流失被害）

都市部津波遡上解析の立体視例（防災教育等を通してソフト防災に貢献）

Structural Analysis Lab. ,Kyushu University

5

マルチスケール津波解析例
5.0

Wave height on A-B
BA

2.5

Spatial scale for 3D sim.

54 [min]3D simulation period

Time after tsunami initiation in 2D simulation

43 [min]66[min]110[min]-5.0

5Velocity[m/s]

Wave velocity[m/s]in S-N direction

Velocity in main flow direction A
from 2D analysis (level0)
6-2.5

100

50-5-10-50.0

B

Temporalscale

MScPHy(Multi-Scale and PHysics)ﾏﾙﾁﾌｨｼﾞｯｸｽ粒子法
複合災害予測し、ハード防災に貢献

流体力を受ける構造 流体-剛体連成SPH法(数億粒子まで対応可)

構造
洗掘＋液状化による構造物の不安定化
剥落？

地盤

水理
胸壁の破壊

浸透流・洗掘解析用SPH法
基礎地盤の洗掘 (実験との比較で精度確認中)
Structural Analysis Lab. ,Kyushu University

衝撃的な波力固体力学
力のつり合いσ n  t
外力

u0

変位の拘束応力評価

Lagrange微分

D u 1  σ  b2Dt  x2

 ij  Dijkl kl

初期状態からの変形量（ひずみ）から応力を算出

Lagrange記述に基づく
有限要素法

数値流体力学
運動量保存式
Euler微分ナビエ・ストークス方程式
移流項

v
1 σ 1
 v v     bt x 

変位から速度へ移流項が質量保存式

v  0
応力評価

 ij   p ij   ij

運動の記述方法は？

 vi v j  ij   
x
xji



瞬間のひずみ速度から応力を算出

Euler記述（or Lagrange記述）

空間の離散化方法は？

有限要素法，有限体積法，粒子法．．．．

運動記述方法の特徴

記述方法

メリット

デメリット

Euler記述（界面捕捉法）・計算ロバスト性が高い

・細かなメッシュが必要

・メッシュによる破綻を気にしなく良い

・高性能な移流スキームが重要

Lagrange記述（界面追跡法）・粗いメッシュでも比較的高精度・計算ロバスト性が低い

・計算時間が短い・メッシュの再構築（リメッシュ）が必要

メッシュに起因するデメリット

Lagrange記述のメッシュレス法があれば，両者の恩恵が受けられる？

固体力学 vs 流体力学Dv σ  b
Dt
(力のつり合い式・運動量保存式)

 ij  Dijkl kl
固体力学アプローチ

流体力学アプローチ

変位を独立自由度に弾性体

 ij   p ij   ij
速度を独立自由度に弾塑性体

非ニュートン流体応力

せん断応力

   so   τ y , ij  K kk  ij  2G ijYield stress
K : bulk modulus G : shear modulus

ひずみニュートン流体 vi v j  ij    x j xi y
せん断ひずみ速度流体と礫(固体)の連成解析におけるマクロとミクロモデル
マクロモデルミクロモデル

流れの特徴

流体の平均的な流れを求める

流体の動きを詳細に解く
流体力

圧力から数値的に求める流速から経験則により求める
流体解像度

固体に対して 1Τ8 程度

固体粒子径と同等の大きさ

12ミクロモデルによる剛体と水の連成解析
剛体計算 → ケーソンブロックの挙動– 捨石マウンドとの接触力→ DEM– 津波により壁面に作用する流体力 → ISPH
ケーソン壁面粒子：DEMの接触判定球＋ISPHの壁粒子［並進運動］［回転運動］𝑀

𝑑𝑽= 𝑀𝒈 + 𝑭𝑐 + 𝑭𝑓
𝑑𝑡𝑑 𝑰𝜴
= 𝑴𝑐 + 𝑴𝑓
𝑑𝑡

𝑀：ケーソン質量，𝑽：速度， 𝜴：角速度，
𝑭𝑐 , 𝑴𝑐 ：接触力/トルク，𝑭𝑓 , 𝑭𝑓 ：流体力/トルク※回転運動は四元数を用いて計算．

𝒏𝑖［接触力/トルク］
𝑭𝑐 =

෍𝒇𝑐_𝑖

𝑴𝑐 =𝑖∈caisson

෍𝒇𝑐_𝑖 × 𝒙𝑖 − 𝑿𝑔

津波𝑖∈caisson［流体力/トルク］
𝑭𝑓 =෍(−𝑃𝑖 𝑎𝑖 )𝒏𝑖𝑖∈caisson

𝑴𝑓 =

෍

(−𝑃𝑖 𝑎𝑖 )𝒏𝑖 × 𝒙𝑖 − 𝑿𝑔𝑖∈caisson

𝑃：ISPHにより求まる壁面の圧力， 𝑎：粒子の投影面積（粒径2） ， 𝒏：外向き法線ベクトル，※初期値はCADから入手→𝑰：慣性テンソル， 𝑿𝑔 ：重心ベクトル

捨石マウンド

𝑃𝑖

ケーソン
ブロック

剛体と水の連成解析（ミクロモデルの解析例）

水を除いて表示

剛体と水を同時に表示

14

ISPH法による流体解析（マクロモデルへの準備）
間隙率 𝜀 を考慮したDarcy-Brinkman型支配方程式

𝐷50 = 1𝑐𝑚𝐷50 = 0.5𝑐𝑚

𝐷50 = 0.3𝑐𝑚：遷移流（地表流 ⇔ 浸透流）を計算可能
・質量保存則：ഥ𝑓𝐷 𝜌ҧ𝑓𝒗+ 𝜌𝑓ҧ 𝛻 ∙=0𝐷𝑡𝜀・運動方程式：𝒗𝑓
𝐶𝑟 (𝜀) 𝐷ഥ
1
ഥ𝑓 + 𝑭𝑟 𝜀= − 𝛻𝑃 + 𝒈 + 𝜈𝐸 𝜀 𝛻 2 𝒗
𝜀 𝐷𝑡
𝜌𝑓ഥ𝑓 = 0𝛻∙𝒗

𝜀：間隙率，𝜌ത𝑓 = 𝜀𝜌𝑓 ：見かけの流体密度，𝑭𝑟 𝜀 ：土粒子から受ける抵抗力，
ഥ𝑓 = 𝜀𝒗𝑓 ：Darcy流速， 𝜈𝐸 (𝜀) = (𝜈𝑤 + 𝜈𝑇 )/𝜀 ：有効動粘性係数, 𝜈𝑤 ：流体動粘性係数， 𝜈𝑇 ：渦粘性係数𝒗

地表流混相流（𝜀 = 1.0）（𝜀 < 1.0）

Navier-Stokes 式

拡張Darcy則𝐷𝜌𝑓+ 𝜌𝑓 𝛻 ∙ 𝒗𝑓 = 0𝐷𝑡
𝐷𝒗𝑓
1= − 𝛻𝑃 + 𝒈 + 𝜈𝐸 𝛻 2 𝒗𝑓
𝐷𝑡
𝜌𝑓

𝜀 < 0.8：浸透状態

水塊の粒子充填層透過（𝜀 ≅ 0.43）

土粒子から受ける抵抗力

𝜈𝑤 1 − 𝜀 2(1 − 𝜀) 2𝑭𝑟 𝜀 = 𝛼𝑐𝜀𝒗𝑟 + 𝛽𝑐𝜀 𝒗𝑟 𝒗𝑟
22𝜀𝐷50𝜀 𝐷50𝛼𝑐 = 150; 𝛽𝑐 = 1.75：Ergunの式（球形粒子群）
𝛼𝑐 = 200; 𝛽𝑐 = 1.10：H.Akbari(2014)の式（捨石マウンド等）

𝜀 ≥ 0.8：浮遊状態

3 𝜀 1−𝜀𝑭𝑟 𝜀 = − 𝐶𝑑 2.7𝒗 𝒗
4𝜀 𝐷50 𝑟 𝑟 （C. Wen. and Y. Yu. ,1966）※浸透/浮遊状態で抵抗力は閾値（ 𝜀 < 0.8）で不連続→ 連続的に両状態を接続する処理を実施

土粒子と水の連成解析（マクロモデルの解析例）
DEM : Discrete Element Method（個別要素法）(Cundall and Struck, 1979)–––

粒状体などの不連続体のラグランジュ型解析手法
各種外力を考慮した運動方程式を各粒子，各時間で解く．
ケーソン－粒子，粒子間の接触力：バネ‐ダッシュポットモデル

6000［並進運動］
浮力𝑚𝑠

接触力

抗力（流体の相互作用）𝜀 ≥ 0.8

１DEM粒子あたりの力に換算※浸透/浮遊状態で抵抗力は𝜀 < 0.8で不連続→連続的に両状態を接続する処理を実施

contact

0［回転運動］𝑑𝝎𝒔
𝐼= ෍ 𝒎𝑐 − 𝒎𝑟𝑑𝑡

𝑚𝑠 ：粒子質量，𝒗𝐬 ：粒子速度，𝒈：重力加速度，𝒎𝑐 ：接触力によるトルク，𝒎𝑟 ：転がり摩擦

粒子の集団沈降解析

Pressure [Pa]𝒇𝑑 𝜀 =

𝜀 < 0.82000

𝜈𝑤 1 − 𝜀 21−𝜀 2𝑉𝑠𝜀𝒗𝑟 + 𝛽𝑐𝜀 𝒗𝑟 𝒗𝑟
22
𝜀𝐷501−𝜀𝜀 𝐷503 𝜀 1−𝜀𝑉𝑠− 𝐶𝑑 2.7
𝒗𝑟 𝒗𝑟
4
𝜀 𝐷50
1−𝜀𝛼𝑐

4000𝑑𝒗𝐬= 𝑚s 𝒈 − 𝛻𝑃𝑉𝑠 + ෍ 𝒇c + 𝒇𝑑 (𝜀)𝑑𝑡流体‐地盤‐剛体連成解析手法
マクロモデル ミクロモデル

津波，流体

ケーソンブロック

捨石マウンドマクロモデル
浸透流

接触力

抗力

ケーソン壁面に作用する流体力

ミクロモデル

蛇篭によるマウンド保護工のシミュレーション
蛇篭による捨石マウンド保護工の洗掘低減効果の確認
蛇篭：金網に砕石やコンクリート片を詰めたもの．
通常，河川の護岸や斜面補強などに利用．
ケーソンブロックと異なり，透水性を有する．

Reference: Japan Jakago Association

近年，捨石マウンドを被覆する洗掘対策工法として注目．
ハザリカら（2016）：実験から蛇篭の持つ洗掘低減効果を確認［解析モデル］
上面から

正面図

実験用の蛇篭

40

固定
壁境界70

水柱(SPH)15

ケーソン（固定壁境界） 越流

5

蛇篭 (DEM)

捨石マウンド(DEM)

10408050154010浸透流

基礎地盤（多孔質領域）[cm]

構成粒子の抗力を外力とし，剛体計算を行う．

220*Hemanta Hazarika, Kengo Nishimura, Babloo Chaudhary, Model testing on resilient solution for breakwater protection against tsunami, Japanese Geotechnical Society Special Publication, Vol.3, Issue 2, pp.40-44, 2016.

蛇篭によるマウンド保護工のシミュレーション
蛇篭による捨石マウンド保護工の洗掘低減効果の確認

洗掘

巻上げ

①無補強 (1.5 s)②蛇篭補強 (3.5s)

洗掘を防止
流失

③蛇篭，重量２倍 (5.0s)

蛇篭：洗掘の直撃を防ぐ洗掘低減効果（目地に越流水が打込まれると流失の危険性有）

実務利用に向けて・施工時に必要最低限の砕石の充填率の目安を提示
・蛇篭の流失を防ぐ方法の検討②蛇篭補強（流失時の動画）

固相から液相への相変化
高粘性流体とし、斜面崩壊を表現

弾塑性体として，斜面崩壊する場所を特定
塑性ひずみにより固相から液相へと相変化
Structural Analysis Lab. ,Kyushu University

20スパコン(TYPE-I)からGPU(TYPE-II)へ
解析コードの特徴
・ラグランジュ型粒子法（運動と伴に粒子を移動させる）
・完全陽解法(ニュートン流体、3000万粒子/1GPU程度まで)から完全陰解法(非ニュートン流体、500万粒子/1GPU程度まで)までを統合的に
・１GPU vs 10 CPU 程度のパフォーマンス
・1次元領域分割(現時点では4GPUまでP2Pで通信)
現在、2次元スライスグリッド＋MPI機能を鋭意GPU解析コードの特徴
解析コードの特徴
・ラグランジュ型粒子法（運動と伴に粒子を移動させる）
・完全陽解法(低粘性ニュートン流体、3000万粒子/1GPU程度まで)から完全陰解法(高粘性・非ニュートン流体、500万粒子/1GPU程度まで)までを統合的に
・1次元領域分割(現時点では4GPUまでP2Pで通信)
MPI版も動き始めた(性能は不十分、まだまだ改良の余地あり)
Projection法を使った予測子計算(圧力勾配以外の速度予測)＋修正子計算(圧力ポアソン)
陽解法 3×粒子数nの成分の更新（速度について陽的に）陰解法 3×粒子数nの連立方程式（速度について陰的に）粒子数nの成分の更新（圧力について陽的に）粒子数nの連立方程式（圧力について陰的に）共役勾配法(CRSで保存した行列に対するスパースソルバー)
GPUで実行するために…
1．cuSPARCEによる行列ベクトル積の演算
2．ベクトルの内積は内製コード…
3．Thrustライブラリーによるソート(粒子番号ﾘﾅﾝﾊﾞﾘﾝｸﾞ)
4. CUDAFortran使用

これまでCPU(京，Intel系)ベースでチューニングしたコードがTYPE-Iではそのまま高速に！
TYPE-IIでは十分なGPU(V100)が搭載されていて高速な演算ができる！