先駆的科学計算に関するフォーラム2015～先端的計算科学研究プロジェクト成果報告～

沿岸構造物におけるFEM-SPH法による大規模シミュレーション
荻野正雄 (名大・情報基盤センター),
浅井光輝 (九大・工学研究院社会基盤部門),
田上大助 (九大・マスフォアインダストリ研究所),
渡部善隆 (九大・情報基盤研究開発センター)
2015年4月24日(金)津波による複合被害
構造
洗掘＋液状化による構造物の不安定化

流体力を受ける構造

剥落？
転倒？

地盤

水理
胸壁の破壊

洗堀

基礎地盤の洗掘

衝撃的な波力研究成果の概要• 2012年度(1年目)の成果– [FEM] 地盤-橋梁モデルにおける反復法の収束性改善– [SPH] 緩和パラメータ・乱流モデル導入等による安定化– 構造物・海底を含む高精度地形モデル構築• 2013年度(2年目)の成果– [SPH] 津波荷重の実験値との比較による妥当性確認– [SPH] 流体-構造(剛体)連成による橋桁流出シミュレーション– 防災・減災向け映像教材の作成• 2014年度(3年目)の成果– [SPH] 圧力境界条件処理方法の改善– [SPH] 剛体移動処理方法の改善– [SPH-FEM] 流体-構造(弾性体)連成による道路橋モデル解析– 津波防災教育向けに高知市モデルの作成と解析

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

流体構造連成解析
レベル1解析

レベル2解析

解像度:1m～2m

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

解像度:10cm～20cm

津波遡上現象

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

シミュレーションベース・ハード防災～橋梁・堤防の安全性評価～

粒子法における基礎式(Lagrange記述)
移流項を無視∂vi∂vi 1 ∂σ ij 1+ vj−= bi∂t∂x j ρ ∂x j
ρσij =− pδij + 2µDij
圧力項

粘性項1  ∂vi ∂v j +Dij =2  ∂x j ∂xi µ 1Dv
1= − ∇p + ∇ ⋅  ∇v  + bDtρρ ρStabilized ISPH (projection method)
Our formulation

Original ISPH
Mass conservation law

Mass conservation law

Dρ+ ρ∇ ⋅ u = 0Dt

Dρ+ ρ∇ ⋅ u = 0Dt

×Pressure Poisson eq.

∇ 2 pin+1

Pressure Poisson eq.

ρ=
∇⋅u∆t

∇ 2 pin+1*0iρ − ρρ
=∇⋅u +∆t∆t
0*in +1nii

2Stabilized ISPH
Proposed PPE∇ 2 pin +1 =ρ0∆t∇ ⋅ ui* + α

Relaxation coefficient：ρ 0 − ρ in∆t 2α (0 ≤ α ≤ 1)

The second term in our formulation should be vanished in the incompressible condition as theoriginal formulation. Then, we put a small α.

陰的なNeumann境界処理法の提案
陰的ノイマン境界処理法の基本原理∇ Pi2n +1

ρ i − ρ in C r (ε i )  ∇ ⋅ v i∗ ρ i=+ αε i2εi ∆t∆t


A

Pwater=Pfreesurface
Pwall

b=0
ディリクレ
境界条件

水粒子だけの行列に縮約化

A’
壁面上に配置したマーカーを介してノイマン条件を満たすように係数行列を修正．

PwaterPfreesurface
Pwall

==0
ディリクレ
境界条件

b'1b2陰的なNeumann境界処理法の提案
解析結果(理論値との比較)
2000

圧力(Pa)

1800

理論値

1600

陰的ノイマン法
提案手法

1400

従来法

1200
1000
800
600
400
200
0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
高さ(cm)

シミュレーションベース・ハード防災
研究テーマ②～橋梁の流失被害予測～
レベル0解析
解像度:10m～100m

レベル1解析

レベル2解析

解像度:1m～2m

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

解像度:10cm～20cm

津波遡上現象

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

解析モデル
橋桁

津波
橋脚
初速度 10m/s

単位 : m

12.2

10.0m

平面図

113.3

40.0

20.0

24.0

30.0

10.32

15.0

18.0

正面図

25.0

[MN]

上部構造固定時の流体力
8水平力

流体力

6上揚力

4200123[sec]

時間24橋桁流失解析(流体剛体連成)

解析条件

粒子間距離

総粒子数時間増分

実時間

解析時間(500node=8000コア)

6cm

55 million

0.001s3s

7h25

SPH-FEM片方向連成解析の流れ

粒子

粒子生成

SPH法解析

境界上節点CAD
データ

境界上節点抽出
プログラム

モデル形状
CAD
データ

: 粒子法 preprocess: FEM preprocess: 開発ツール

粒子圧力

SPH-FEM
カプラー
境界上要素面への圧力

メッシュ情報メッシュ生成境界条件付加メッシュFEM解析

境界条件
付加メッシュ26

構造物境界面の設定
 粒子モデルは構造物の内部を充填するように生成 粒子モデルにメッシュの境界面にある節点の位置を中心とするダミー粒子を追加
 SPH法の補間計算方法に基づいてダミー粒子上の物理量を計算

メッシュの境界面上の節点を抽出
改善前
: 水粒子
: 構造物粒子ツールを開発改善後: ダミー粒子（境界面上節点）: 節点

メッシュ境界面

実際の境界面27荷重境界条件の設定
 境界面上の節点に、対応するダミー粒子の持つ物理量を持たせる
 境界面上の要素面を構成する節点の圧力平均を要素面に加える
ツールを開発X-Y-Z図Y-Z図
𝒑𝒑𝒕𝒕𝟏𝟏YZ
𝑘𝑘

1
𝑝𝑝𝑗𝑗𝑡𝑡 = � 𝑝𝑝𝑖𝑖𝑡𝑡
𝑘𝑘𝑖𝑖=1𝑝𝑝2𝑡𝑡

𝑝𝑝3𝑡𝑡

𝑝𝑝4𝑡𝑡

𝒑𝒑𝒕𝒕𝟏𝟏

𝒑𝒑𝒕𝒕𝟐𝟐

𝒑𝒑𝒕𝒕𝟑𝟑

𝑝𝑝1𝑡𝑡

YZ
X

𝑝𝑝𝑖𝑖𝑡𝑡 ∶ ある時刻の要素面𝑗𝑗を構成する節点𝑖𝑖に対応する粒子の保持する圧力
𝑘𝑘 ∶ 要素面𝑗𝑗を構成する節点の数

28

津波伝播問題設定

 歌津大橋の震災で流失した第8径間の断面を元に作成
 津波高 : 15𝑚𝑚
 𝑥𝑥 = 20 𝑚𝑚 以下の水粒子に初速度10 𝑚𝑚⁄𝑠𝑠𝑠𝑠𝑠𝑠 を常に与える
 先端の波には初期条件に流速10 𝑚𝑚⁄𝑠𝑠𝑠𝑠𝑠𝑠 を与える
自由落下
総粒子数時間増分

実時間6𝑐𝑐𝑐𝑐

約5,500万0.001 𝑠𝑠𝑠𝑠𝑠𝑠

2.55 𝑠𝑠𝑠𝑠𝑠𝑠

正面図113.3m

Z30.0m

YX40.0m

非段波状X

20.0m

24.0m

15.0m

18.0m

10.0m

12.2m

初速度 10m/s

10.32m平面図

粒子間隔25.0m29

構造解析用メッシュの作成
 歌津大橋の震災で流失した第8径間の断面を元に作成
 橋桁と橋脚は接触部で剛結合
 4面体1次要素

198[𝑐𝑐𝑐𝑐]

834[𝑐𝑐𝑐𝑐]

24830[𝑐𝑐𝑐𝑐]

節点数

要素数

338,835

889,250

1000[𝑐𝑐𝑐𝑐]

181032[𝑐𝑐𝑐𝑐]

𝑌𝑌

𝑌𝑌
𝑋𝑋𝑍𝑍

1220[𝑐𝑐𝑐𝑐]30

構造解析用境界条件

物性値設定（鉄筋コンクリート） ヤング率 = 5.0 × 1010 𝑁𝑁⁄𝑚𝑚2
 ポアソン比 = 0.2
 質量密度= 2.45 × 104 𝑁𝑁⁄𝑚𝑚31. 橋脚底面(𝑦𝑦 = 0.0)を地面で固定されていると仮定して、 x, y, 𝑧𝑧方向変位を拘束

2. 橋桁は𝑧𝑧軸方向に十分長いと仮定して、面 (z = 1.1, 𝑧𝑧 = 11.1) の𝑧𝑧方向変位を拘束
1.

2.

境界上要素面数

Pressure値 有り

Pressure値 無し

時間増分

実時間

662,934

284,335378,559

0.001 𝑠𝑠𝑠𝑠𝑠𝑠

0.97 𝑠𝑠𝑠𝑠𝑠𝑠

31

SPH法解析結果

SPH法による解析 : 圧力コンター図1201201209090906060603030300 [kPa]

0 [kPa]0 [kPa]𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 0.10[𝑠𝑠𝑠𝑠𝑠𝑠]

𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 0.30[𝑠𝑠𝑠𝑠𝑠𝑠]

𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 0.50[𝑠𝑠𝑠𝑠𝑠𝑠]

12012012090909060606030300 [kPa]

𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 0.70[𝑠𝑠𝑠𝑠𝑠𝑠]

300 [kPa]

0 [kPa]𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 0.90[𝑠𝑠𝑠𝑠𝑠𝑠]

𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 1.10[𝑠𝑠𝑠𝑠𝑠𝑠]

32

FEM解析結果
FEM解析 : 節点応力( Nodal Equivalent Stress )コンター図

𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 0.00[𝑠𝑠𝑠𝑠𝑠𝑠]

𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 0.10[𝑠𝑠𝑠𝑠𝑠𝑠]

𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 0.40[𝑠𝑠𝑠𝑠𝑠𝑠]

𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 0.55[𝑠𝑠𝑠𝑠𝑠𝑠]

𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 0.70[𝑠𝑠𝑠𝑠𝑠𝑠]

𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 0.972[𝑠𝑠𝑠𝑠𝑠𝑠]

33

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

レベル1解析

レベル2解析

解像度:1m～2m

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

解像度:10cm～20cm

津波遡上現象

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

解析場所の選定
南海トラフ地震
将来M8～9の発生確率30年以内に約60～70%(地震調査研究推進本部)

南海トラフでの最大クラスの地震による高知県の主な津波被害予想
死者数 4万9千人 避難者数 51万人
被害額 10兆6千億円
地震発生から約3分で津波高1mの津波が到達
最大津波高 34m(全国1位)(内閣府より)南海トラフ(文部科学省 地震調査研究推進本部)

最も人災が予想される
高知市の津波遡上解析

解析モデル(一部領域のみ)
約2250m

津波
高さ 5m
初速度 5m/s

約1800m(高知県南海トラフ地震対策課より)

津波 約1000m

解析条件

九州大学スーパーコンピューターの8node使用

総粒子数

粒子間隔

時間増分

総ステップ数

実時間

解析時間

約794万

2m

0.005sec

33000step

165sec

約89時間まとめと今後の課題• SPH法による自由表面流れ解析の精度改善– 圧力境界条件の処理見直し– 浮体等を剛体で取り扱う場合の処理見直し• SPH-FEM連成による流体構造(弾性体)連成解析の実装– 微小変形を前提とした片方向連成– 既存のSPHコードとFEMコードに対するカプラー開発– 歌津大橋道路橋モデルにおける津波荷重による応力解析を実現• 防災教育への貢献– 高知市モデルの作成と解析– 教育向けCG映像の作成• 今後の課題– 流体(SPH)-構造(FEM)の双方向連成解析