狠狠撸

aokomoriuta (2013-10-12): 計算シミュレーション勉強会#1

粒子法の復習
（陽解法と陰解法の比較から）

青子守歌


はじめに

-2-

＜はじめに＞


おことわり

時間の都合等により
不正確な表現を
多数含んでいます
-3-


自己紹介

-4-

＜自己紹介＞


名前とアバター

青子守歌 / aokomoriuta

※call me “あお”
-5-

＜自己紹介＞


いろんなとこのアカウント

twitter: @aokomorita
github: aokomoriuta
web: j.mp/ao_komoriuta
-6-

＜自己紹介＞


専門分野

粒子法（MPS法）を用いた計算力学
＠土木工学（主に水工学分野）
分類：
土木工学
＞水理学
＞計算力学
-7-


＜自己紹介＞

計算力学とは

海岸工学?水工学分野において
File:Gileston20.jpg
by Mick Lobbi : CC-BY-SA 2.0 Generic

File:Coastal defences, Beesands - geograph.org.uk - 69707.jpg
by Nigel Chadwick : CC-BY-SA 2.0 Generic

File:Dolos.jpg
by Adam Brinki : CC-BY-SA 3.0 Unported

現実空間
実験空間

計算空間

代替となりうるか？

-8-

＜自己紹介＞


なので

精度がとても大事！
（超重要）
-9-


陽解法vs陰解法の結果

- 10 -


＜結果＞

いきなり比較1

陽解法

→resultExplicit1.mp4

陰解法

- 11 -

→resultStandard1.mp4


＜結果＞

いきなり比較2

陽解法

陰解法

- 12 →resultExplicit2.mp4



＜結果＞

いきなり比較3

陽解法

陰解法




陽解法？陰解法？

- 14 -

＜陽解法？陰解法？＞


陽関数と陰関数

陽関数：
陰関数：

y= f ( x)

f ( x , y)=0
- 15 -



陽関数と陰関数の例

陽関数：
陰関数：

y= √ 1? x
2

2

x + y =1

2

- 16 -



陽解法と陰解法

陽解法 x= f ( x )
：
陰解法 f ( x , t )=0
：
- 17 -

x ：未知
t ：既知



陽解法と陰解法の例

dx
dy
=sin ( y ) , =cos( x)
　　　　　　　　　　に対して
dt
dt

= x +sin ( y ) Δ t

陽解法 x
：

k +1

陰解法 x
：

k +1

k

k

= x +sin ( y
k

k

k

k +1

)Δt

x , y：k番目の値
- 18 -



？？？

つまり

- 19 -



陽解法と陰解法とは

陽解法：既知の値だけで
どうにか解こう
陰解法：未知の値は
未知の値のまま解こう
- 20 -


粒子法における
陽解法と陰解法
- 21 -

＜粒子法における陽と陰＞


粒子法とは（狭義）

例えば：非圧縮性ニュートン流体のナビエ?ストークス方程式：
2

? ui
? ui
? ui
1 ?p
+u j
= f i? ρ
+ν
?t
?xj
? xi
? xj? xj
のような

と

? ui
=0
? xi

連続体を計算点（粒子）に
離散化して
ラグランジュ的に解く方法
主にSPH法とMPS法
- 22 -



弱圧縮SPH法

Smoothed Particle Hydraulics
非圧縮性流れ用

弱圧縮SPH

(Lucy, 1977)

(Monagahan, 1994)

1. 重力?粘性項を計算する
2.弱圧縮で圧力を計算

- 23 -



（通常の）MPS法

Moving Particle Semi-implicit
(Koshizuka&Oka,1996)

1. 重力?粘性項を計算する
2.非圧縮で圧力を計算

- 24 -



粒子法の詳細

また機会があれば
そのうち
末尾の参考文献もどうぞ
- 25 -



第一段階

1. 重力?粘性項によって、
粒子を移動させる

圧縮！？
- 26 -



第ニ段階

弱圧縮

2. 圧力を計算

少しなら

圧力

非圧縮
絶対に×
- 27 -



弱圧縮の場合

(1)密度

ρ i =∑ m j W ( r ij ) 圧力が
j≠i

陽的に
(2)圧力
γ
ρi
解ける！
pi =B ρ ?1
0

[( ) ]
- 28 -



非圧縮の場合

(1)密度

※非圧縮＝密度が変化しない

ρ i =∑ m j W ( r ij ) 圧力が
j≠i

(2)圧力

pi =?? ? ? ? ?
- 29 -

陽的に
解けない！



非圧縮の場合どうする？(1)

(1)密度を戻す

密度の修正量←速度の修正量
Δρ
+ ρ ? Δ u=0
Δt
Dρ
+ ρ ? u=0
Dt
連続式　　　　　　　より
- 30 -




(2)圧力勾配で戻す

速度の修正量←圧力勾配
Δu
1
=? ρ ? p
Δt
Du
1
2
= f ? ρ ? p+ν ? u
Dt
NS式　　　　　　　　より
- 31 -




(3) 1と2を混ぜる

(1)に代入←(2)の発散
Δρ
+ ρ ? Δ u=0
Δt

?Δu
1 2
=? ρ ? p
Δt

ρ Δρ
? p=?
Δt Δt
2

- 32 -




(4)ラプラシアンモデル
2

? p=?

ρ Δρ
Δt Δt

2d
? pi =
∑ ( p j ? pi ) W ( r ij )
λ n0 j≠i
標準MPS法の例
- 33 -




(4)圧力のポアソン方程式
標準MPS法の例

ρ n ' ?n0
2d
k +1
k +1
∑ ( p j ? pi ) W ( r ij )=? Δ t 2 n
λ n0 j≠i
0

陰解法なら解ける！
- 34 （ただの連立一次方程式）
-



粒子法における陽と陰は

弱圧縮

陽解法

非圧縮

陰解法
- 35 -



SPHとMPS

陽解法
（弱圧縮）
SPH WCSPH

Monagahan, 1994

MPS

E-MPS
Oochi&Koshizuka, 2010

- 36 -

陰解法
（非圧縮）
ISPH

Shao&Lo, 2003

（オリジナル）MPS
Koshizuka&Oka, 1996


陽解法と陰解法の比較

- 37 -

＜陽と陰の比較＞


比較1

陽解法

→resultExplicit1.mp4

陰解法

- 38 -




比較1から

陽解法

ふつう

陰解法

ふつう
- 39 -



比較2

陽解法

陰解法





比較2から

陽解法

ぽよぽよ

陰解法

ふつう

弱圧縮のせい？

- 41 -



比較3

陽解法

陰解法





比較3から

陽解法

沈んだ！

陰解法

ふつう

弱圧縮のせい！

- 43 -



1ステップの計算時間

陽解法

超速い

陰解法

遅い
- 44 -

連立一次方程式
の解法があるから



1[s]の計算時間

陽解法

速い？

時間刻み小さい

陰解法

遅い？
時間刻み大きい
- 45 -


まとめ

- 46 -

＜まとめ＞


粒子法の陽解法?陰解法

[陽解法]
[陰解法]

?弱圧縮性
　→精度が落ちる
?1ステップの計算時間は速い
?時間刻み小さい
　→トータルで速い???？
?非圧縮性
　→精度はそれなり
?1ステップの計算時間は遅い
?時間刻みは大きい
　→トータルで???？
- 47 -

＜まとめ＞


参考文献

?越塚誠一(2006)：『計算力学レクチャーシリーズ5 粒子法』、丸善。
　→粒子法/MPS法をやるならまずコレ
?後藤仁志(2004)：『数値流砂水理学 -粒子法による混相流と粒状体の計算
力学-』、森北出版。
→格子法?粒子法の分類などに詳しい、DEMについての情報もあり
?土木学会編(2012)：『数値波動水槽 -砕波波浪計算の深化と耐波設計の革
新を目指して-』、「8.粒子法による数値波動水槽」、丸善。
→高精度化などの最新情報はココにまとまっている
?W.G.Hoover(2006)：『粒子法による力学 -連続体シミュレーションへの展
開-』、志田晃一郎訳、森北出版。
→SPH法についての本、やや難読感あり
- 48 -


return 0;
Otherwise noted, all text and images are available
under the Creative Commons Attribution-Share Alike or Attribution-Noncommercial 3.0 Unported.
一部で特別な指定があるものを除き、全ての文章と画像は
クリエイティブ?コモンズ表示-継承または表示-非営利 3.0 非移植で利用可能です。
CC-BY-SA: http://creativecommons.org/licenses/by-sa/3.0

CC-BY-NC: http://creativecommons.org/licenses/by-nc/3.0

狠狠撸

[計算シミュレーション勉強会#1] 粒子法の復習（陽解法と陰解法の比較から）

Recommended

More Related Content

What's hot (20)

Viewers also liked (20)

Recently uploaded (8)

[計算シミュレーション勉強会#1] 粒子法の復習（陽解法と陰解法の比較から）