Smoothed Particle
Hydrodynamics
岡本 崇 (筑波大学計算科学研究センター)

SPH法の特徴• 粒子法である• 流体を粒子の重ね合わせで表現• Lagrangian 的方法• 粒子と共に動く座標系• 高密度領域で高解像度• 宇宙分野で広く使われている• 天体は密度揺らぎが重力崩壊することにより形成されるから

長所と短所
短所• 長所•
密度の高い所で自動的
密度の低い所の解像度が••に解像度も高くなる• ガリレイ不変• 多次元の実装が簡単• 安定低い• 低精度（空間0次精度）• 不連続面が苦手• 安定 (変なことをしていても止まらない)Basics• ある領域での物理量 f(x) の平均値を有限で球対称な
kernel W(x;h)! を用いて次のように表す!f "(x) = W (x − x! ; h)f (x! )dx!
!
whereW (x)dx = 1 and W (x) = W (−x)!
• 微分は部分積分を用いて!∇f #(x) =W (x − x! ; h)∇x! f (x! )dx!
!
= − {∇x! W (x − x! )}f (x! )dx!
!
={∇W (x − x! )}f (x! )dx!

精度
• f(xʼ) を x の周りで展開 f (x)(2)•f (x! ) ! f (x) + f (1) (x)(x − x! ) +(x − x! )2 + · · ·2!
!f "(x) = W (x − x! ; h)f (x! )dx! に代入
h2eff 2!f "(x) = f (x) +∇ f (x) + O(h4eff )! 4

where h2eff ≡ 2

x2 W (x; h)dx

空間2次精度離散化• 積分の離散化•miと体積要素を置き換え,dxi =ρi!
mjρjj!
mj"= 1W(x−x;h)j
一般にはρjj!f "(x) #•W (x − xj ; h)fj!f "(x) # f (x)•!
j

W (x − xj ; h)なのでmj+ ···
ρjとなり、空間０次の項がすでに誤差を含む
これが SPH 法最大の問題である(ただし、Inutsuka 2002 の定式化ではこの問題は解消されている)

定式化(1)• 物理量 f として密度ρを考える! m!
ρi ! "ρ#(xi ) !

ρjjj

ρj

W (rij ; h) =mj W (rij ; h)j

• Lagrange 形式の流体の方程式
ρ̇v̇u̇
P= −ρ∇ · v,
1= − ∇P,
ρP= − ∇ · v,
ρ= (γ − 1)ρu

定式化(2)• 連続の式は不要• 運動方程式は• 第一項:• 第二項:v̇ i = −!

Pi
ρi"

∇Pi
=∇
ρi!

# Pj mjPi
ρi"+Pi∇ρi
2ρiを用いて∇W (rij , h)ρρjjj
PiPi !
∇ρi = 2mj ∇W (rij , h)
ρ2i
ρi j

∇

∇Pi=−
ρi=!
jmj

"PjPi+
ρ2j
ρ2i#∇W (rij , h)

• i と j に対して反対称であり、作用反作用の法則を満たす

定式化(3)• エネルギー方程式
ρi ∇ i · v i= ∇i · (ρi v i ) − v i · ∇i ρi!
! mj(ρj v j )∇i W (rij , h) − v i
mj ∇i W (rij , h)=ρjjj!
=mj (v j − v i )∇i W (rij , h)jPi !
Pi∴ u̇ = − ∇i · v i =mj (v i − v j ) · ∇W (rij , h)ρiρi j

• 次の対称化された形よりも性質が良い"#!
1u̇i =2mjjPi
Pj+ 2
ρ2i
ρj(v i − v j )∇i W (rij , h)Artificialviscosity(1)• 衝撃波を扱うためには人工粘性を用いる必要がある• 人工粘性項をQ として運動方程式, エネルij

ギー方程式はそれぞれ次のようになる"#v̇ iu̇i!

Pj
Pi= −mj
2 + ρ2 + Qij ∇W (rij ; h)ρjij$%!
Pi1=mj+ Qij v ij · ∇W (rij ; h)ρ2ij

• 何処にどれだけの粘性が必要かは計算をやってみるまで分からないArtificialviscosity(2)(1997) のものを紹介• Monaghan!
Qij =αvsig eij
ρ̄ij

r ij · v ij < 00otherwise

1
1
where ρ̄ij = (ρi + ρj ) and vsig = (ci + cj − 3v ij · eij )22• また、shear|∇ で人工粘性が働かないよう·v |fi =ii|∇i · v i | + |∇ × v i | + !ci /hiを導入し,1
Q̃ij = Qij (fi + fj ) を用いることが多い (Balsara ʼ95)
2

SPH Kernel• Gaussian と spline kernel が良く使われる
spline kernel• 今回用いる
30 ≤ u ≤ 21 ,12 ≤ u ≤ 1,otherwise,

spline

1.2
1W(u)


1 − 6u + 6u
σ W (r; h) = D 2(1 − u)3h 

02

1.4

r
where u = ,h4 40
8and σ = ,, andfor D = 1, 2, and 3
3 7ππ

0.8
0.6
0.4
0.2
0-1 -0.8 -0.6 -0.4 -0.2

0u

0.2 0.4 0.6 0.8

1Smoothing length• 近傍の平均粒子間隔の2倍程度以上離れた粒子からの寄与は考えない
h∼
2!

miρi

" D1

h 内の粒子数を Nngbに固定• h は粒子ごとに異なる• W(r ; h)→W(r , h ), W(r , h)→ W(r ; h )• 運動方程式は i, j に対して反対称でなくなるijiji∇W (rij ; h) → ∇W̄ij , where W̄ij =ij

ij

i

1
[W (rij , hi ) + W (rij ; hj )]2

最小作用の原理からの導出• Springel & Hernquist (2003)• 内部エネルギー, u, ではなく、エントロピー•

変数 A(s) = (γ − 1)uρ1−γ (∴ P = A(s)ργ ) を用いる

NN!
1!
1L(q, q̇) =mi ṙ 2 −mi Ai ργ−1i
2 i=1
γ − 1 i=1

where q = (r 1 , . . . , r N , h1 , . . . , hN )•拘束条件

φi (q) = αhDi ρi − mp Nngb = 04
where α = 2, π, and π, for D = 1, 2, and 3, respectively.
3
N!
d ∂L
∂L
∂φj−=λjdt ∂ q̇i∂qi∂qij=1

導出(続き)!
"#$−1
mi Pihi∂ρiλi =1+2Dαh ρiDρi ∂hi• h での偏微分から• これを用いて方程式の最初の半分からi"#$%−1Pjhj∂ρjmi v̇ i =mj 2 1 +∇ i ρjρDρ∂hjjjj=1N!

∇i ρj = mi ∇i W (rij ; hj ) + δijv̇ i = −fi =!

N!
j=1

mj"N!

k=1

mk ∇i W (rik ; hi )に注意し,#Pi
Pj
fi 2 ∇i W (rij ; hi ) + fj 2 ∇i W (rij ; hj )ρiρjhi ∂ρi
1+Dρi ∂hi"−1

Entropy• 断熱では dA/dt = 0• 衝撃波によるエントロピー生成だけを考えれば良いNdAi
1γ−1!
=mj Qij v ij · ∇i W̄ij
γ−1dt
2 ρij近傍粒子探査• 全ての粒子について例えば密度を!
Nρi =mj W (rij ; hi )jのように求めると、その手間はN2となる• 今の場合, r > h で W = 0 なので h 内に含まれる粒子だけ効率良く見つけてくれば良い.

• 重力計算が Tree 法の場合は Tree が, P M の場合
3は chain-list がそのまま使える時間積分• 一般的は Leap-frog 法が良く使われる• あまり精度の良い時間積分を用いても仕方がない• 自己重力系と異なり、加速度を求めるのに位置だけではなく速度も必要なため、一次精度で外挿して得られた速度を用いる• 内部エネルギー（またはエントロピー) も普通にやると時間一次精度になる.

実習: Shock tube
• １次元問題• 領域の左側と右側に密度・圧力の異なる静止した流体を置き、時刻 0 に間の壁を取り除く→ 衝撃波の発生
ρL = 1
PL = 1vL = 0

ρR = 0.25
PR = 0.1vR = 0Shock-tube• Makeﬁle#OPT += -DPERIODIC=1.0
OPT += -DONE̲D̲SHOCKTUBE#OPT += -DTWO̲D̲KH#OPT += -DTHREE̲D̲SEDOV#OPT += -DREDUCEVISC#OPT += -DACのように編集してコンパイル• 実行%cd shocktube と移動して、%../sph shocktube.paramと parameter ﬁle を与えると実行できる.

Parameter files• shocktube.param
OutputFileBase shocktube 結果の出力用ファイル名
Ngrid
400 粒子数 (１次元あたり)
Gamma
1.6666666666667 γ
ArtBulkViscConst 1.0 人工粘性の強度
CourantFac
0.25
時間刻みを調節BegTime0
計算開始時刻MaxTime0.1
計算終了時刻DtOutput
0.01
結果出力間隔Nngb
5SPH 計算に用いる近傍粒子の数

結果の確認• “solution.txt” に t = 0.1 での解析解が与えてある (x, v, ρ, P の順)• t = 0.1 の結果は, “shocktube̲10.txt”という ﬁle に出力されている.
# x, v, ρ, u, P

• シミュレーション結果と比較してみよ

結果
衝撃波面
密度
膨張波

接触不連続面
接触不連続面で圧力が不連続に！
圧力

人工粘性の改良• αが定数ではなく衝撃波の外側では減衰するように
αα̇ = −

τ+ αc max(0, −∇ · v) みたいな感じに• Makeﬁle の -DREDUCEVISC の前の#を取って再コンパイル• %cd shocktube• %../sph shock̲variable̲alpha.param
この parameter ﬁle の中では新たに “MinAlpha”,“ViscDecayLength” が与えられている。MinAlpha はαの下限値を、ViscDecayLength は smoothing length の何倍の距離で MinAlpha まで減衰するかを与えている.

• 実行して結果を前と比較せよ

人工熱伝導率の導入• 接触不連続面で圧力の高い方から低い方へ熱を流してやり, 圧力を連続的にする!
N"
dui !!
mj uvsig αc (ui − uj )r̂ ij · ∇i W̄ij
=dt !cond
ρ̄
ijj

• 問題点

• どうやって接触不連続面を見つけるか？
• 熱を流し過ぎると負のエントロピーを持つものが出てくる可能性が...

とりあえず実行• MakeﬁleOPT += -DONE̲D̲SHOCKTUBE
OPT += -DREDUCEVISC
OPT += -DACとしてコンパイル• 実行%cd shocktube%../sph shock̲ac.param• “shock̲ac.param” 内の “ArtBulkCondConst” が人工熱伝導率を決めている• 実行して結果を他の場合と比較せよ. とくに接触不連続面での圧力の振る舞いに注目.

Kelvin-Helmholtz
instability• 2次元問題• 速度の異なる流れの境界面上に発生する不安定性

ρH=2PH = 2.5vH = -0.5

今回の実習のセットアップ

ρL=1
PL = 2.5vL = 0.5

KHI in SPH• Okamoto et al. 2003 や Agertz et al.
2007 によると, SPH では KelvinHelmholtz 不安定性が解けない

メッシュ法

SPH

実際に動かす• Makeﬁle#OPT += -DONE̲D̲SHOCKTUBE
OPT += -DTWO̲D̲KH#OPT += -DTHREE̲D̲SEDOV#OPT += -DREDUCEVISC#OPT+= -DAC

• OpenMP を用いる場合は SYSTEM の部分をSYSTYPE = "with̲omp"#SYSTYPE = "without̲omp”
また、OPT += -DOMP=2 のようにしてコア数を指定する• kh に移動し, %../sph kh.param のようにして実行• tatara でやる場合は SYSTEM=FUJITSU ただし、今回の code はtataraでは絶望的に時間がかかるため、お勧めしない.

結果• output ﬁles
• kh̲*.txt という名前 (* は数字)
• x, y, vx, vy, ρ, u, P, の順
圧力
接触不連続面で圧力が大変なことに！

不安定性が発生しているがmixing が生じていないy人工熱伝導率を使う• Makeﬁle#OPT += -DONE̲D̲SHOCKTUBE
OPT += -DTWO̲D̲KH#OPT += -DTHREE̲D̲SEDOVOPT += -DREDUCEVISC
OPT+= -DAC

• kh3.param を使う

結果

２相が混ざっている！

熱伝導により接触不連続面での圧力の非物理的ななジャンプがなくなっている！

点源爆発• 圧力の無視できる一様な気体の中での爆発• 球対称なので本来１次元問題なのだが、今回は３次元で解く• 自己相似解が知られている (Sedov-Taylersolution)

実行• Makeﬁle#OPT += -DONE̲D̲SHOCKTUBE#OPT += -DTWO̲D̲KH
OPT += -DTHREE̲D̲SEDOV#OPT += -DREDUCEVISC#OPT+= -DAC

• %cd sedov として, %../sph sedov.param• “sedov̲raw” に t = 0.04 での解が与えてあるので sedov̲4.txtと比較せよ• output ﬁle (x, y, z, v , v , v , ρ, u, P, r, v )• sedov̲raw (r, v , ρ, P)
xryzr

結果
• みt
圧力

密度

半径
衝撃波が少しなまっているが解析解と良く一致

半径大惨事. vr も比較してみること.
論文では密度しか見せない人も多いので騙されないように.

時間変動する人工粘性• -DREDUCEVISC の前の#をとってコンパイル• sedov2.param を使う• sedov2̲*.txt という名前で結果が出力される• 前とあまり結果が変わらないことを確認せよ

人工熱伝導率も使う• -DAC の前の#もとってコンパイル• sedov3.param を使う• sedov3̲*.txt という名前で結果が出力される• 結果を解析解, また人工熱伝導率無しのシミュレーションと比較せよ.

結果
密度圧力

•

半径
半径
密度はノイズが減ってきれいになったが、衝撃波面が実際よりも前進している• 圧力も改善されたが変な構造が見える• シェルの内側の熱が外側へ輸送されてしまうのが原因• 適切なパラメータを決めるのは困難

まとめ• SPH は, 比較的歴史の浅い流体法なのでまだ発展途上• 適切に人工粘性を用いれば, post-shockの物理状態はおおむね大丈夫• 衝撃波がなまっては駄目な問題, 乱流が重要になるような問題には適さない

参考文献
••
•
•
•
•
•

Monaghan, 1992, Ann. Rev. Astron. Astrophys, 30, 543 (レビュー)
Springel & Hernquist, 2002, MNRAS, 333, 649 (今回用いた SPH)
Monaghan, 1997, J. Compt. Phys., 136, 298 (人工粘性)
Balsara, 1995, J. Compt. Phys., 121, 357 (shear-reduced artificialviscosity)
Morris & Monaghan, 1997, J. Compt. Phys.,136, 41 (time-dependent
artificial viscosity)
Price, 2008, J. Compt. Phys., 227, 10040, (人工熱伝導と KHI)
Inutsuka, 2002, J. Compt. Phys., 179, 238, (Riemann solver を用いたSPH. 人工粘性いらず！)••

Okamoto et al. 2003, MNRAS, 345, 429 (SPH では KHI が解けない)
Agertz et al. 2007, MNRAS, 380, 963 (SPH と AMR の比較)