SPH
法による流体と剛体の連成シミュレーションCoupled simulation of fluid and rigid body using SPH method
物理学専攻 本間慎吾
Homma Shingo
目的
本研究では、まず
SPH
法が得意とする非圧縮性流体の自由表面流れの数値計算プログラムを作成し、ベン チマーク問題であるダムブレークのシミュレーションを行うことによりプログラムの評価を行う。次に非圧縮 性流体と剛体の連成問題として、ダムブレークと剛体の落下を組み合わせたシミュレーションと、石の水切り 現象のシミュレーションを行ない、その挙動を計算することで、SPH
法の有用性について調べる。内容
SPH(Smoothed Particle Hydrodynamics)
法とは、流体を粒子の集まりとして近似してラグランジュ的に 計算するシミュレーション手法である。計算空間を格子で分割し、固定された格子点で物理量を離散化する格 子法と比べ、複雑に変形する自由表面流れの計算が得意であるという利点を持っている。流体の支配方程式である
NS
方程式Dv Dt = − 1
ρ ∇ p + ν ∇
2v + a (1)
を
SPH
法で解くためには、まず計算空間における物理量を計算する。空間のある座標r
での物理量ϕ(r)
は、近傍粒子が離散的に保持する物理量の足し合わせとして
ϕ(r) = ∑
j
m j
ϕ j ρ j
W (r − r j , h) (2)
と表せる。
m j
、ρ j
は粒子j
の持つ質量と密度で、W
、h
はカーネル関数と物理量の広がりに関する長さであ る。さらに物理量の勾配∇ ϕ(r)
は∇ ϕ(r) = ∑
j
m j
ϕ j
ρ j ∇ W (r − r j , h) (3)
とカーネル関数の微分量のみで記述することができる。
NS
方程式をSPH
法を用いて離散化すると次のようになる。Dv i
Dt = − ∑
j
m j (
p i
ρ
2i + p j
ρ
2j )
∂W ij
∂r i + ∑
j
m j µ i + µ j
ρ i ρ j v ij ( 1
r ij
∂W ij
∂r ij )
+ f i (4)
NS
方程式には2
階の微分項である粘性項が存在するが、SPH
法で2
階微分項を離散化するには工夫が必要 なため、上式の右辺第2
項のような人工的な粘性に置き換えている。f i
は外力項である。非圧縮性流体を記述する場合、連続の式
∇ · v = 0 (5)
1
を解かなければならないが、
SPH
法では連続の式を解かないため、代わりに密度と圧力を関連付ける状態方 程式を与える。粒子i
の密度は式(2)
よりρ i = ∑
j
m j W (r i − r j , h) (6)
となる。流体の圧力は上式で求められた密度から
p i = { c
2ρ
0γ
{( ρ
iρ
0) γ
− 1 }
ρ i > ρ
00 ρ i < ρ
0(7)
を経て決定される。
c
は流体の音速、γ
は正の定数であり、基準密度ρ
0より低下した粒子は自由表面にあると みなし圧力をゼロとする。音速c
を代表的な流速に比べて大きく設定しておけば、密度のゆらぎが小さく抑え られ実質的に非圧縮性流体を扱うことができる。剛体の計算について、
SPH
法では1
個の剛体を複数の粒子で表し、流体と相互作用させつつ運動させる。各時間刻みにおいて剛体粒子の重心変化量と回転角変化量を計算し、それらが保存するように新たな粒子の座 標を決定する。
成果
SPH
法により、流体物理の基礎的な問題を解くことによりプログラムの評価を行った。ダムブレークのシ ミュレーションではKoshizuka
ら[8]
が開発をしたMPS
法での計算結果とSPH
法での計算結果について、水平到達距離の比較を行ない、良い一致を得ることができた。図
1
は水柱の幅L
、高さ2L
としてシミュレー ションした結果を示す。両手法での誤差はNS
方程式の離散化方法や、粘性計算の違いによるものと思われる。1 1.5 2 2.5 3 3.5 4 4.5 5 5.5
0 0.5 1 1.5 2 2.5 3 3.5 4
Z/L
t(2g/L)
1/2SPH MPS
図
1
両手法による水平到達距離計算プログラムの検証の結果、良い一致を得られたため、非圧縮性流体と剛体の連成シミュレーションに移 行した。
ダムブレークと剛体の落下のシミュレーションを行なった。流体粒子のパラメータは前述のダムブレークの 場合と同じであり表面張力は計算していない。比重が流体より大きい場合と小さい場合で行なった。
流体と剛体が同時に崩壊と落下を始め、衝突する。その際、密度の違いによる剛体の運動の様子を確認する ことができた。比重が小さい場合、流体からの影響を強く受け、流れに乗って運ばれている様子が確認でき た。また剛体が流体から受ける浮力についても確認することができた。
2
石の水切りのシミュレーションを行った。粘性は水と同程度であり、流体と剛体の比重
σ = 2.7
としてい る。剛体である円板の大きさは半径r = 2.4cm
、厚さd = 0.6cm
を想定している。図
2 t = 0.00s
図3 t = 0.03s
図4 t = 0.06s
Nagahiro
ら[5]
によると水切り現象を微分方程式モデルで考えた時に、速度の反転が起こる最小速度v
minを表す式は
v
min=
√ 2gr cos(ϕ + θ)
√
ξ
∗sin ϕ + cos ϕ
2S
wetsin
2ϕ (8)
また速度の反転が起こる最大の入射角
θ
maxはθ
max= arccos
√ 2 F
(
ξ
∗sin ϕ + cos ϕ 2π sin
2ϕ
)
− ϕ (9)
と表される。ここで
S
wet= π
として近似し、ξ
∗はフィッティングパラメータで、ξ
∗= 2.0
とする。またF
はフルード数であり、F = v
20/gr
である。この解析解とシミュレーションによって得られた結果を比較することでシミュレーションの再現度を考 える。
図
5
は式(8)
においてθ = 20
として、ϕ
とv
minの関係をプロットしたもので、赤の実線が理論値、青の点 がシミュレーション結果である。同じく図5
は式(9)
においてv
0= 3.5
としてϕ
とθ
maxの関係をプロット したものである。双方ともϕ > 20
の範囲で理論値とほぼ一致する傾向がみられたが、ϕ < 20
の範囲では理論 値よりも大幅なずれが生じてしまった。これは流体表面に衝突した際に発生する圧力の疎密波が底面で反射を して、再び剛体に衝突したことによる影響が大きいためと考えられる。1 2 3 4 5 6 7
0 10 20 30 40 50 60 70
VminPHI ODE model
SPH result
図
5
入射角θ = 20
°の平面内における反発領域10 15 20 25 30 35 40 45 50 55 60
0 10 20 30 40 50 60 70
THETAmaxPHI ODE model
SPH result
図
6
初速度v
0= 3.5m/s
の平面内における反発領域3
結論
ダムブレークにおける低レイノルズ数での自由表面流れのシミュレーションでは、
SPH
法での計算結果と 他の手法での計算結果や実験値との良い一致を得た。自由表面流れの非圧縮性流体と剛体の連成シミュレー ションとして行なったダムブレークと剛体の落下の組み合わせでは、比重の違いにより剛体の運動の様子が変 化しているのが確認できた。水切りのシミュレーションでは、圧力の疎密波の影響が大きい場合を除き、理論 値と傾向がほぼ一致する結果を得た。以上の結果から、本研究で使用した
SPH
法によるシミュレーションは信頼度が高いと考えられる。しかし 高レイノルズ数流れをSPH
法を用いて計算する場合、非圧縮性がより重要になり本研究のような状態方程式 を用いた手法では適用が難しい。またシミュレーション精度を高めるために、粒子数を増加させるとともに計 算時間を短縮させる必要がある。さらに、物理法則に基づいた粘性計算、厳密な境界条件の取り扱いなどが今 後の課題である。参考文献
