2014/8/4

津波被害におけるマルチフィジックス

構造

粒子法シミュレーションの津波防災・減災への応用

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

流体力を受ける構造

剥落？

地盤

衝撃的な波力水理
胸壁の破壊

九州大学大学院工学研究院
社会基盤部門 浅井光輝
基礎地盤の洗掘
Structural Analysis Lab.

Structural Analysis Lab. ,Kyushu University

1

マルチスケール津波解析の必要性

代表的な粒子法（SPH vs MPS）

研究テーマ① シミュレーションベース・ソフト防災～3D ・CGによる津波防災教育～
レベル0解析
解像度:10m～100m

レベル1解析

レベル2解析

解像度:1m～2m

Smoothed P HMoving P S
articleydridynamics

has been proposed by Lucy, Monaghan and in 1977

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

article

解像度:10cm～20cm

津波遡上現象

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

emi-implicit method

研究テーマ② シミュレーションベース・ハード防災～橋梁・堤防の安全性評価～
Structural Analysis Lab. ,Kyushu University

3粒子法における基礎式（ラグランジュ記述）

MPS

Structural Analysis Lab. ,Kyushu University

SPH法における粒子離散化
Smoothed Particle Hydrodynamics(SPH)

Lagrange記述へ

物体を有限個の粒子に離散化し、その任意の点の物理量を近傍粒子の物理量の重みつき平均として近似

vi
v 1  ij 1
 vj i  bi
tx j  x j f ( x)   f ( x J )W ( x  x J , h)dV

ij   pij  2Dij D  12  xv  vx ij 1Dv1  p     v   bDt SPHSPHStructural Analysis Lab. ,Kyushu Universityijji

粒子間距離Nf ( x)  
対象粒子(i)の物理量

J 1mJf ( x J )W ( x  x J , h)Jf(x)I

近傍粒子(j)の物理量粒子近似

f(x J)
J
影響半径

重み関数

MPSW

解析対象

＜主な計算手順＞① 近傍粒子検索（各粒子Iに対するI-Jのペアとその総数N）② 重み関数の形状
SPH ③ 独立変数の更新方法
Structural Analysis Lab. ,Kyushu University

1

2014/8/4

重み関数（カーネル関数の役割）

SPH法における重み関数の条件

fj
fi

SPH＝FEM的な内挿近似

差分法的な微分の考え方＝MPS

J JNN①関数f (x)の近似
mmJf ((xx))
(x ( xxJ , xh)J , h)
 JN Jf (fx(j x)WJ )WJ J11 m
f ( x i )   j f ( x j ) W重み付き(r )j 1 W ( x  x ', h)dx '  1

Delta function property

Nj 1mjj re  1 (0  r  re )W (r )   r
(r  re ) 0コーン型

つりがね型②勾配・発散・ラプラシアンを導出f ( x )  

自由境界付近では精度が悪い
1.0  3 2 3 3c1  2 r  4 r  c2  r 2W (r )  40Normalization condition( or unity condition)

f ( x i ) h 0

粒子は均等に分布することが望ましい

Compact support condition
計算効率のため（近傍粒子内での相関を考えれば良い）W ( x  x ', h)  0,①勾配・発散・ラプラシアンの近似

f ( x j ) W (r )limW ( x  x ', h)  ( x  x ')where x  x '  hd N f j fi( x j  x i ) W (r )n0 j 1 x j  x i 2

Structural Analysis Lab. ,Kyushu University1.0

近傍粒子の差分の重み付き平均

独立変数の更新（古典的なSPH法）

具体的な関数形により変化する定数

疎 密スプライン関数がよく用いられる．

（hが決まれば，関数値が決定→関数形は固定している）

SPH

Structural Analysis Lab. ,Kyushu University

独立変数の更新（射影法の導入→非圧縮性流体）
非圧縮性NS方程式

START

 1
vin 1  vin  t   pin   2vin  bi  ①



初期条件の入力ni  m jW ( xi  x j )速度と位置の更新（ﾅﾋﾞｴｽﾄｰｸｽ方程式)

次の時間ステップへno  n  i  1
pin  B  0 


圧力評価（状態方程式）vi  vin  t  2vin  bi②

 1
vin 1  vin  t   pin   2vin  bi  

yes② pressure evaluation

 2 p n1 速度の境界条件（slip orv'  Mv ,

水粒子

境界面

仮想マーカー

Dv
nDt

－(1)

Dv
Dt

境界面
壁粒子

非すべり条件Dv
n  0Dt⇒点対称な速度V’を壁粒子に付与

Vv  Rv , Rij  δij
'

境界面V’

壁粒子v'   Mv  1   Rv
すべり条件

－(1)

水粒子の漏れを防ぐ⇒速度勾配の法線方向成分はゼロ－(2)

式(1)を満足するようにナビア・スト－クス式を参照し、非一様的圧力ノイマン条件を満足最適な条件(すべり条件と非すべり条件の間)

壁領域

SPH
11MPS

圧力の境界条件(ノイマン条件)

no-slip）

M ij  δij  2 ni n jv  0

Structural Analysis Lab. ,Kyushu University⇒鏡映対称な速度V’を壁粒子に付与

仮想マーカーV’

SPH

すべり条件
水粒子

質量保存則

1vin1  vi*  t  p n1 

水粒子

仮想マーカー
・水領域に配置・壁粒子と一対一に対応

水領域

  v

③ corrector stepENDMPSt

ポアソン方程式から圧力の導出③Structural Analysis Lab. ,Kyushu University

圧力勾配項を除いたNS式から仮の速度を導出

final time stepSPH

×

① predictor stepj

0    1

非すべり条件Dv1
  0 p    2 v  gDt

f－(3)

MPS

SPH
12p
  f nn－(2)

MPS

2

2014/8/4

シミュレーションベース・ソフト防災～3D ・CGによる津波防災教育～

高精度な地形モデルの復元<手順>①航空測量データ②平面データ作成③粒子発生

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

岩手県県土整備部河川課株)防災技術コンサルタント よりデータ提供

航空レーザ測量 + 深浅測量
地上のデータ（約1m間隔）

海底のデータ航空機から地上にレーザを照射し、反射波を観測

レベル2解析

レベル1解析
解像度:1m～2m

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

解像度:10cm～20cm

津波遡上現象

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

Structural Analysis Lab. ,Kyushu University

Structural Analysis Lab. ,Kyushu University

13

数値地図モデルと高精度モデル

粒子解析モデルへの変換<手順>①航空測量データ②スムージング③粒子発生

解析のための粒子を発生

これまでのデータ改善した地形モデル

津波

カラーコンター：標高10mの解像度
Structural Analysis Lab. ,Kyushu University

1mの解像度StructuralAnalysis Lab. ,Kyushu University
120m
-10m

田老地区での災害の再現例

津波遡上解析結果

岩手県・(株)防災技術コンサルタント→航空測量データを提供

4m間隔で粒子を発生＝1000万粒子

流入領域

V=10m/s
H=3m
Structural Analysis Lab. ,Kyushu University

Structural Analysis Lab. ,Kyushu University

3

2014/8/4

被害調査結果との比較

九州大学周辺での仮想津波被害
disaster investigation report

Simulation result

RED：tsunami inundated zone
BLUE：serious damage zone
Interferential crinkle
Structural Analysis Lab. ,Kyushu University

1920

シミュレーションベース・ハード防災
研究テーマ②～津波避難ビルの安全性評価～
レベル0解析
解像度:10m～100m

レベル2解析

レベル1解析
解像度:1m～2m

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

解像度:10cm～20cm

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

津波遡上現象Structural Analysis Lab. ,Kyushu University

Structural Analysis Lab. ,Kyushu University

橋梁および堤防の崩壊

22

津波被害におけるマルチフィジックス

構造
洗掘＋液状化による構造物の不安定化流体力を受ける構造剥落？

地盤

衝撃的な波力水理
胸壁の破壊

基礎地盤の洗掘
Structural Analysis Lab. ,Kyushu University

23

Structural Analysis Lab. ,Kyushu University

4

2014/8/4

津波波力算定式 (国総研：津波避難ビル等の構造上の要件の解説，2012)

解析モデル・解析条件

朝倉式による津波波圧の発生を仮定し，次式にて津波波力を定義する．

Qz  g z 1 (h  z ) Bdzz

●解析条件

2

粒子間隔

津波波力(Qz)

0.1m

時間増分

実時間

約3400万 0.001sec

総粒子数3.7sec

津波10.010.8 津波避難ビル
7m/s

50.0

5.0

17.2
3.0103.8

設計用浸水深構造物
αh

津波単位 : m(A)段波z2hz1αρgh

20.032.012.0<側面図>

初期条件：全ての水粒子に7m/sの速度を付与．
境界条件：津波避難ビルから10m離れた位置以降に存在する水粒子に継続的に7m/sの速度を付与．

10.0
51.0

71.0

10.0
50.0-52.0<正面図>

α は水深係数（1.5~3.0の値をとる）と呼ばれ，地形等の条件により軽減される相似の比率である．

解析結果-開口部有(背面)解析結果-開口部有(前面)014000

28000

42000

56000

70000 [Pa]

解析結果（構造全体に作用する力）

津波避難ビルモデル(a)開口部無単位 : m
2.80
2.85
2.85
2.85
2.85
3.00

●参考文献：国土交通省国土技術
政策総合研究所 一般社団法人
建築性能基準推進協会 協力 独立行政法人 建築研究所：津波避難ビル等の構造上の要件の解説

0.492
0.231(b)開口部有

0.506
0.245

津波波力[MN]

51.0
6階
5階
4階
3階
2階
1階

100000.873(c)ピロティ※赤字は幅に対して開口部が占める割合を示す．

(a)開口部無-未加工データ(a)開口部無-フィルタ処理後(b)開口部有-未加工データ(b)開口部有-フィルタ処理後(c)ピロティ-未加工データ(c)ピロティ-フィルタ処理後12時間[sec]

3

未加工データ：サンプリング周波数1000Hzでデータをプロット
フィルタ処理後：100Hzのローパスフィルター用いて処理したデータをプロット

5

2014/8/4

解析結果（フロア毎の波力）

津波波力算定式との比較●津波波力の平均値(フィルタ処理後)●フロア毎の津波波力（最大重複波力，衝撃津波波圧）
最大重複波力(平均値)
津波波力算定式

6030

00

1

2

3

1階
2階
3階
4階

60

30

00

1

2

390

津波波圧[kPa]

1階
2階
3階
4階

90

津波波圧[kPa]

津波波圧[kPa]

90

1階
2階
3階
4階

60

最大重複波力(平均値)津波波力算定式

14.54
42.45
3210Qz / gmax B(a)開口部無43210Qz / gmax B(c)ピロティ30※衝撃津波波力

最大重複波力(平均値)
津波波力算定式00

12

時間[sec]時間[sec]時間[sec](a)開口部無(b)開口部有(c)ピロティ

無次元衝撃 無次元衝撃 無次元津波 無次元津波

34.25
432

10(b)開口部有Qz / g max B

津波波力

津波波力

(1階)

(2階)

波力算定式 波力算定式
(1階)

(2階)

開口部有4.25

2.69

1.96

1.55

開口部無14.54

7.04

2.70

2.12

ピロティ

2.45

2.19

0.88

1.55

シミュレーションベース・ハード防災
研究テーマ②～橋梁の流失被害予測～結論
波頭形状■最大重複波圧は実スケールにおいても朝倉式にて評価することは可能

レベル0解析
解像度:10m～100mである．
レベル2解析

レベル1解析

■段波による衝撃的な津波波圧は朝倉式が対応するものではなく，別途解像度:1m～2m衝撃的な圧力の評価が必要である．

開口形状■津波波力算定式による津波波力の見積もりは最大重複波力に対す

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

解像度:10cm～20cmるものであり，衝撃的な波力はそれを上回る．
■開口部やピロティ構造には衝撃及び重複波力の両者を大きく低減

津波遡上現象する効果があることを確認した．また，その効果は国総研による算定式よりも期待できる可能性がある．

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

Structural Analysis Lab. ,Kyushu University

解析対象

34

補足資料：歌津地域の損傷状況
宮城県南三陸町歌津大橋

山側で支圧破壊

東日本大震災の津波被害により計8径間(第3径間～第10径間)の橋桁が流出

プレテンション式

ポストテンション式

ポストテンション式

プレテンション式

南三陸町

橋脚部の主鉄筋が破断
歌津大橋(震災後)※「Craft MAP-日本・世界の白地図」を基に作成

35

写真：東北地方整備局震災伝承館参照

写真：九州工業大学 幸左らによる津波被害分析.pdf(H.25.2.20)より62014/8/4

橋桁詳細モデル

補足資料：歌津地域の損傷状況
歌津大橋における流失した橋桁の移動距離

ポストテンション式
プレテンション式

流出後

歌津大橋モデル図
830cm
24cm

198cm

18cm

画像：九州工業大学 幸左らによる津波被害分析.pdf(H.25.2.20)より

剛体移動モデルの定式化

断面図（歌津大橋第8径間の断面形状を再現）
38

剛体の扱い方 (1/3)
・剛体移動モデルのアルゴリズム• 剛体移動の扱い方①

剛体粒子も水粒子と同様に有限個の粒子に離散化する

流体運動に制約を与えて擬似的に剛体を表す手法→速度ベースの定式化

剛体水粒子離散化
剛体粒子に作用する力を計算

剛体粒子に作用する力を流体力と外力項から計算• 剛体移動の扱い方②V

基本物理に従い境界面での力の伝達として流体剛体問題を解く手法→外力ベースの定式化

剛体速度の更新

ωv kl 1  V  ω  rk
並進速度

疑似的剛体モデル概念図

回転速度

dV
M Mg  F f  Fedt

Idω
 M f  Medt

39

剛体の扱い方(2/3)
並進運動の運動方程式

dv
M Mg  F f  Fedt
Ff 

剛体の扱い方(3/3)
水粒子

回転運動の運動方程式

水粒子

剛体粒子

dω
I
 M f  Medt

剛体粒子

流体力を受ける剛体粒子
固体間接触力を受ける剛体粒子
壁粒子

on the surface P S nii

流体力を受ける壁粒子ii

Mf 

仮想マーカー（剛体粒子と対応）
壁粒子の受圧方向＝壁の内向き法線
固体間接触力の作用方向

仮想マーカー（壁粒子と対応）

固体間接触力を受ける剛体粒子
壁粒子

 r  r  P S n

on the surface

icii

流体力を受ける壁粒子

ii

固体間接触力を受ける壁粒子

Fe  k 3 / 2

流体力を受ける剛体粒子

Me 

固体間接触力を受ける壁粒子 r  r  F

on the surface

icei

仮想マーカー（剛体粒子と対応）
壁粒子の受圧方向＝壁の内向き法線
固体間接触力の作用方向

仮想マーカー（壁粒子と対応）

ヘルツの弾性体接触理論

  r lk
並進速度の更新

ucn 1

回転速度の更新Ei E j4 r
223 1  i E j  1  j Ei

FF  uc  t  g  f  e M Mnωcn 1 ωc  I 1t M f  M e n

最終的な剛体速度の更新

M f : 剛体表面に作用する流体力によるモーメント力M e : 剛体表面に作用する他の外力によるモーメント力
I : 慣性テンソル
ωc : 回転速度
ui : 剛体速度72014/8/4

剛体流失シミュレーション

解析モデル＜箱桁断面モデル＞

11m

10.0m

＜平面図＞

速度ベースの定式化

75.0m

＜側面図＞

外力ベースの定式化(a) 箱桁モデル

10m/s

5.0m

37.0m

1.5m

8.0m

22.0m

※剛体(橋桁)密度2450kg/m3

橋脚固定5.0m 3.0m 10.0m

35.0m(b) Π桁モデル
20.0m

解析条件
粒子間隔

総粒子数時間増分

実時間

解析時間(8node)25cm

約80万0.001sec

2.4sec

約2時間
44

剛体流失シミュレーション

解析モデル
橋桁＜Π桁断面モデル＞
外力ベース定式化(正面)

外力ベース定式化(鳥瞰図)津波

橋脚
初速度 10m/s

単位 : m

12.210.0m

平面図

113.3

40.0

20.0

24.0

30.0

10.32

15.018.0

正面図

25.0
45

46

橋桁流失動画[MN]

上部構造固定時の流体力

解析条件

8水平力

流体力

6上揚力

粒子間距離

総粒子数

時間増分

実時間

解析時間(1024core)

6cm

55 million

0.001s3s

140h

420012時間

Structural Analysis Lab. ,Kyushu University

3[sec]

47488

2014/8/4

Damage mechanisms of breakwater

マルチフィジックス津波解析

1. horizontal force due tothe water level difference
between the front and rear
breakwater.
Scouring

Water pressureOverflow

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

流体力を受ける構造

剥落？

地盤Collapse!
Seepage flow

衝撃的な波力水理

2. Soil scour and/or erosion
behind the seawall during
overflow.
3. piping destruction due todecline of the bearing
capacity by seepage flow.

胸壁の破壊

Fluid-Seepage flowinteraction problem.
基礎地盤の洗掘
Structural Analysis Lab. ,Kyushu University

50

Governing equation

Mixed flow theory
Surface flow
Navier-Stokes eq.
Dv1  0 P   w   T  2v  gDtTsunamiDarcy-Brinkman model1)
Where we solve this equation with changing porosity ε, we can solve Navier-Stoke’s
equation in the pure fluid and Darcy’s law in the porous matrixv
v  v      0 P   E  2v  g v
tK   


？？

？？

ε (porosity) ： liquid volume fractionAkbari model2)
In porous region, we can solve extended Darcy’s law with changing porosity.

Cr   Dv1
  0 P   E   2v  a v  b v v  gDt

Seepage flow
Darcy’s law or extended Darcy’s law 0Cv  P   0 gv
v v  P   0 g
KKK

References :1) M. LE BARS et al. : Interfacial conditions between a pure fluid and a porous medium: implications for binary alloy
solidification, J. of Fluid Mechanics, 2006.
2) Akbari, H. : Modified moving particle method for modeling wave interaction with multi layered porous structures.
Coast. Eng., 2014.

→There need to be a unified formula between these two flows.
51

52

Validation test

Analysis result (movie)

Analysis parameters

Contour : Pressure [Pa]Number of Particles ： 637028; Initial PDTime step
： 0.0005sec ; Relaxation coef.
Boundary condition ： non-slip

Calculate viscosity
using Re number.

E v D50, Re  325
25 Re： 0.5cm： α=0.005
Scale: [cm]Navier-Stokes eqn.

Extended Darcy’s law

Navier-Stokes eqn.

29
28
37
water

2

Soil
230

References :Liu et al. : Numerical modeling of wave interaction with porous structures., J. Waterw. Port Coast. Ocean Eng., 1999.
53

4454

9

2014/8/4

Analysis result (outer surface of waer)
t=0.4st=0.8s

今後の展望
レベル0解析(2D差分法)
解像度:10m～100m

開発項目1
レベル0-レベル1の連成

レベル1解析(3D粒子法)

開発項目2.
ズーミング型粒子法の開発

解像度:1m～2m

レベル2解析(3D粒子法)
解像度:10cm～20cm

t=1.2st=1.6s開発項目3.
0300
600
900
1200[Pa]
土木三力連成解析の改良と検証

３D particle method

２D-FEM or FDM

Global tsunami

Source55

10m~100m

Local tsunami
Bay

1m~5m

Fluid structure
structure
10cm~50cm

10