石の水切りのSPHシミュレーション

仙台高専機械システム工学科
永弘進一郎

2010年12月9日木曜日

様々な衝突の研究
弾性球と弾性壁の斜め衝突
H. Kuninaka and H. Hayakawa
Phys. Rev. Lett. 93, 154301 (2004)
M. Y. Louge & E. Adams
Phys. Rev. E 65, 021303 (2002)

弾性ゲルと壁の衝突
Y.Tanaka et. al.
Europhys. Lett. 63, 146 (2003)

小さな液滴と壁の衝突

R0+X

K. Okumura et. al.
Europhys. Lett. 62, 237 (2003)

J. W. Glasheen and T. A. McMahon,Nature, 380 (1996)
C. Clanet, F. Hersen and L. Bocquet.
Nature, 427 (2004)2010年12月9日木曜日

hh(a)流体表面と物体の衝突

R0+X(b)

物体と水面の斜め衝突：実験
Weight[Cwt.]

distance[Yards]Number ofricochet

8428501665290032

反発の臨界角：√θc = 18/ σ

cWater2010年12月9日木曜日

σ=

物体の密度
流体の密度

G. Birkhoff et. al. (1944)
W. Johnson & S. R. Reid (1975)
I. M. Hutchings (1976)

円柱衝突の理論
Birkhoﬀモデル
1

2

3

5

6y
Oydv
m=dt

w

0!

pn dS

x4

E. G. Richardson, Proc. Phys. Soc.
London, Sect. A 61 (1948)

surface

π sin αpn =ρv 24 + π sin α(Rayleighの式)√θc = 17.8/ σ
2010年12月9日木曜日

Flow direction

物体と水面の斜め衝突：実験（２）フルード数：v2F =∼1g!

v

円板（石）に十分大きい回転を加えて水面へ衝突させる→ φはほとんど変化しない。
初期条件:

入射角度θ
円板の水面に対する迎え角度φ
衝突初期速度 v0

円板の半径：5 [cm], 入射速度 〜 3 [m/sec.]
C. Clanet, F. Hersen and L. Bocquet. Nature, 427 (2004)2010年12月9日木曜日

660

550

入射角度（θ）[deg.]

Umin[m/sec]

石の水切りの実験
Skipping Stone

432
100204060

円盤の傾き角度（φ）[deg.]v = 3.5[m/sec.]

4030
Skipping stone
2010002040

円盤の傾き角度（φ）[deg.]“Magic Angle” φm〜 20°60

θφv

C. Clanet, F. Hersen and L. Bocquet. Nature, 427 (2004)2010年12月9日木曜日

発表の項目
１, 数値計算手法 ­SPH法
２, 水面と円柱の２次元衝突シミュレーション
円柱・流体表面の衝突ショミュレーション。反発のベキ則と比較を行い，SPH法の衝突問題への妥当性と問題点を検証。
３, 水面と円板の３次元衝突とモデル解析
円板と水面の衝突の３次元ショミュレーション。
実験結果の再現。数値解析より得られた結果を元に，円板が受ける力の性質を調べ，衝突を単純化したモデルを解析する。
４, 結論、今後の課題

2010年12月9日木曜日

１, 数値計算手法 ­SPH法
計算例：液滴と水面の衝突

2010年12月9日木曜日

自由表面を扱う数値計算手法
Lagrange描像に基づく手法
・SPH( Smoothed Particle Hydrodynamics)法
・MPS( Moving Particle Semi-implicit ) 法
非圧縮性を保証した粒子法。
Euler描像に基づく手法
・CIP( Cubic Interpolated Propagation)法
・Level-set法
衝突を扱う都合上，運動量やエネルギーの保存が保証されるスキームが望ましい
SPH法

2010年12月9日木曜日

SPH法 (1)W(r j):
Gaussian karnel.

流体の密度ρの補間＜ρ＞：!ρ"(r) =!

vjmW (r − r j )j

これを用いて，任意の物理量Aを! m!A"(r) =Aj W (r − r j )!ρ"jj! m!∇A#(r) =Aj ∇W (r − r j )!ρ#jj!ρ"j ≡ !ρ"(r j )%&!ρ"i2

cρ
−1

0

ρ0

圧力を決定する状態方程式： pi = = 0
2010年12月9日木曜日!ρ"i ! ρ0!ρ"i < ρ0SPH法 (2)Navier-Stokes方程式∂u1+ (u · ∇)u = − ∇p + ν∇2 u +∂tρ!

"ζ+ ν/3 ∇∇ · uρ

! m!A"(r) =Aj W (r − r j )!ρ"jjに基づく補間の手続き
SPHの運動方程式! mbd!u"i = −dt
ρi ρjj"uij · r ij
pi + pj − νξ
2rij#∇i Wij斥力相互作用する粒子系の運動方程式

2010年12月9日木曜日

SPH法 (3) ­ 境界条件
固体の表面に粘性によって固定された流体粒子を仮定。
固定された粒子の受ける力の総和から壁の受ける力を計算できる。

vi
Length ofsmoothing kernel
Fluid particles
fixed on solid wall (Wall particles)

2010年12月9日木曜日

2, 水面と円柱の２次元衝突シミュレーション
円柱と水面の斜め衝突r=1zv

gxLz

Lx

2010年12月9日木曜日

円柱と水面の衝突
ー衝突の様子と速度プロファイル

上図a-d; E. G. Richardson,
Proc. Phys. Soc. London, Sect. A 61 (1948)より引用
2010年12月9日木曜日

斜め衝突の間に円柱が受ける力・粒子の初期配置の影響：粒子の配置を正方格子とランダム配置に選んだ二つの結果を比較
force
v02 r

1.50.6force
0.5v02 r

Square Lattice
Random

1

Square Lattice
Random

0.4

fz

0.3

fz

0.5

0.2
0.1

T=r/cで00-0.1-0.5-1fx

時間平均処理 -0.2

fx

-0.3
-0.4

01234t5670812345t

6782r/v02r/v0・境界(底面)の影響：深さの異なる水槽用いたシミュレーションの結果

0.4

force
0.3
v02 r

fz

Lz=7r
Lz=4r

fz

fx

0.2fz

0.1

fx

0-0.1

fx

-0.2-0.301

23

45t2r/v0
2010年12月9日木曜日

678

反発の臨界角度の再現性√

θc = 18/ σ
100tank size lx=30r, lz=4r
lx=42r, lz=8r

Critical Angle [deg.]

θc-0.5
10

10.1

1

10

Specific gravity
2010年12月9日木曜日

100

3, 水面と円板の３次元衝突とモデル解析

円板と水面の斜め衝突v

2010年12月9日木曜日

水面と円板の衝突の様子１ー円板が回転しながら衝突する場合

5cm

θ = 12.5◦ , φ = 20◦ , ω = 60[rounds/sec.]
2010年12月9日木曜日

水面と円板の衝突の様子２ー円板に回転が無い場合

5cm
θ = 12.5◦ , φ = 20◦ , ω = 0[rounds/sec.]
2010年12月9日木曜日

反発条件の実験との比較
v

7Experiment
SPH5[deg.]

6Stone Skip domain

4Incident Angle

Minimum Velocity [m/sec.]

60

3
2Experiment
SPH

100

2010年12月9日木曜日

10203040
Tilt Angle [deg.]

5060

5040
30
20

Stone Skip domain

1000102030405060入射角・迎え角のTransfer Function
ffv

φ0 = 16.0

θ0 = 8.030

40352530θf

202515201015
1050

502010年12月9日木曜日

10θ203000

10203040水面と円板の衝突のモデル解析
仮定：・回転速度は十分大きい → φ = const.
・衝突による水面の変形を無視。
・粘性を無視して、円盤にはp ∼ ρ(v · n)2で表される圧力のみが働く。
・空気の効果は無視（その密度は水に比べて0.001倍と軽いため。）
円板にかかる力： f =pndS

fxRv0

2010年12月9日木曜日

x'Water!

CD ρ(v · n)2 n/2dS

disk! s face== const.
d=disk! s facegz'z!

1CD ρ(v · n)2 Sn2

衝突のODEモデル（１）
無次元化した運動方程式ẍz̈（xz座標系）
1= − sin φF1! 2= κS(z )ż − cos φFCD Rρκ=2πdρ!
(x,z)：円盤下側の角の位置v02
F =gR：フロード数gz'zfxR= const.

v0d#!
"!
""2!
!
!
!
zzπz
− 1+
S(z ! ) = − arcsin 1 +1− 1+2
sin φ
sin φ
sin φ

流体にひたっている面積。

2010年12月9日木曜日

x'WaterODEモデルとSPHシミュレーションの比較
定数κの決定κ= 0.94,(Fitting parameter)

円板の受ける力

円板の軌道
0.6

0.35

fz

0.3
0.25

0.4

fx

0.3

0.15

0.2

z/r

0.2

0.1

0.1

0.05

00-0.1fx-0.05

2010年12月9日木曜日

00.5-0.2

SPH
ODE-0.1-0.15

11.52time
2r/v02.5ODE
SPH

0.5fz

33.5-0.34-0.40

1234x/r

5678

ODEモデルとSPHシミュレーション及び実験との比較

7Experiment
SPH
Theory65[deg.]Stone-skip domain
Criterion A4Incident Angle

Minimum Velocity vmin [m/sec.]

60

32Experiment
SPH
Theory100

2010年12月9日木曜日

10203040
Tilt Angle [deg.]

50504030
20

Stone-skip domain

106000102030Tilt Angle

4050[deg.]60

衝突のODEモデル̶解析解”反発”の定義…通常は｢水面の高さを基準｣にとる。位置条件
ここでは，「重力方向の速度を基準にとる。速度条件「円盤の重力方向の速度が反転したら，反発が 起こったものと，見なす。」

ẍ
(a)z̈
(b)

(c)
(d)

2010年12月9日木曜日=1− sin φF=1
κS(z )ż − cos φF!

2

変曲点の存在条件をしらべる。

球と水面の衝突の実験で得られた衝突後の球の軌道。
E. G. Richardson, Proc. Phys. Soc.
London, Sect. A 61 (1948)

速度条件を用いた時の反発条件：
解析解と実験・シミュレーションの比較√!
" 12σd cos φ2gRx∗ sin φ +vmin =cos(θ + φ)
CD R sin2 φθmax = arccos!

2F"σd cos φx∗ sin φ +CD R sin2 φ#−φ7Experiment
SPH
TheoryStone-skip domain(速度条件)

5432Experiment
SPHTheory10

2010年12月9日木曜日[deg.]

6Incident Angle

Minimum Velocity vmin [m/sec.]

60010203040Tilt Angle [deg.]

5050403020

Stone-skip domain(速度条件)

106000102030
Tilt Angle

4050[deg.]60

"Magic angle"と入射角度θの関係
35

σ＝2.0Angle[deg.]

30
25201510Experiment
SPH
Theory50

2010年12月9日木曜日

0102030
4050
Incident Angle [deg.]60

反発回数の見積もり

v＝25[m/sec.]

θ＝２０°n>n2002>

10203040Angle of incidence [deg.]

2010年12月9日木曜日

n>10051030n > 3815252050n>

52020n>5103010n0>140n>

15Angle of incidence [deg.]n>

2530n>

Incident speed [m/sec.]

38303035

Tilt angle [deg.]4045

結論とまとめ
SPHによる衝突の直接シミュレーション
・2次元，３次元のSPHシミュレーションは，水面と円柱及び円板の衝突の反発条件を良く再現した。
・円板と水面の衝突について，衝突過程で円板が受ける力積をシミュレーションによって調べ，最適角度φ=20 に対して直感的な理解を得た。
問題点として：平均流速が大きい場合に，粒子法の特性に起因する余分な粘性が表れる。
衝突速度が大きい場合，この粘性の効果によって，現象が正しく再現されない場合がある。
流体表面と物体の衝突のモデル化
・フルード数が１程度の領域に於ける円板と水面の衝突の実験事実が，単純なモデルによって十分に説明できた。
物体と流体間の衝突が極めて単純な力学過程として理解できる事を示した

2010年12月9日木曜日

シミュレーション法の応用分野と今後の課題

・”濡れ”を伴う粉体系へのシミュレーション法の開発→ 地盤の液状化現象への応用
R0+X

・非常に大きな変形をともなう粘弾性体のシミュレーション−ゲルの衝突−水滴の衝突

R0+X
hh(a)

(b)

K. Okumura et. al. Europhys. Lett. 62, 237 (2003)−高圧の物理学の実験
・自由表面を持つ様々な現象への数値的なアプローチ−円形跳水の問題
Clive Ellegaard et. al., Nonlinearity, 12, 1 (1999)
Nature, 392, 23 (1998)2010年12月9日木曜日