Scientific and Educational Reports of the Facultyof Science and Technology, Kochi University
Vol. 2 (2019), No. 1

粒子法（SPH）とその計算プログラムの実装についての解説
An introduction to SPH and its implementaion as computer programs
中込照明
Teruaki Nakagomi〒780-8520 高知市曙町 2 ー 5 ー 1
高知大学理工学部情報科学科

Department of Information Science, Faculty of Science and Technology, Kochi University,
Kochi 780-8520 Japan

概要
本論文では流体力学のシミュレーション分野で使われる SPH（Smoothed Particle

Hydrodynamics，粒子法）の原理とその計算プログラムの実装について解説する．主な特徴は以下の三点である．
（１）SPH における基礎方程式は厳密な導出過程を経て得られるものではなく，かなりの程度発見法的である．そのため，この分野になれていない人には取りつきにくい面がある．
これを改善すべく，定義，要請，仮定を整理し，基礎方程式の導出過程をできるだけ見通しの良いものにした．
（２）SPH 粒子は粒子間相互作用を行う．N 個の粒子があれば，粒子対の個数は N 2 であり，相互作用の計算にかかる計算量は計算式をそのまま実行すると O(N 2 ) となり，通常使
われるパソコンでは N がすこし大きくなるだけで実行困難になる．この点を改良して，計算量を O(N ) に落とすための手法を提供した．
（３）本論文で解説した手法にのっとって作成した SPH のシミュレーション・プログラム
をダウンロードによって提供している．このプログラムでは，２次元，３次元モデルにつき
様々な流体現象がインタラクティブにシミュレートできるようになっている．例えば，揺れ
動く容器の中の流体の動き，水柱の崩壊，容器の穴からの流体の流出，熱対流など．

1 はじめに（SPH とは）
SPH（Smoothed Particle Hydrodynamics) とは連続体の運動を粒子系の運動として近似する方法の一種であり，特に流体運動の近似に使われるものである．
連続体を剛体，弾性体，流体と分けるなら，このうち最も粒子系近似がしやすいのは弾性体である．弾性体の内部に格子状の構造を想定し各格子点に一定質量の粒子を置き，近傍粒子を結ぶ各辺にバネ状の相互作用を設定することにより，その運動をかなり良く近似するこ1
Received: November 11, 2018
Reviewed by anonymous referee(s), and accepted: December 5, 2018
Published: March 27, 2019

とができる．このばねを固くした極限が剛体となるので，剛体も同様に近似できる．この場
合はもちろん完全な剛体ではなく疑似剛体となるが，コンピュータ・シミュレーションは簡単になり，コマの動きなど容易に再現できる．格子構造の特徴は，ある粒子と相互作用する
相手はその粒子の近傍粒子だけであり，しかも，この近傍粒子の集合が決まっていて変化しないという点である．
これに対して，流体を粒子の集まりとして記述しようとするなら，相互作用は近傍粒子に限られるけれども，その近傍粒子の集合は一定ではなく，常に入れ替わり得るとしなければならない．このため格子模型とは異なる近傍粒子間の相互作用が必要になり，考え出されたのが SPH である．SPH でイメージする粒子は完全な質点あるいは剛体球のようなものではなく，その質量が粒子の位置を代表する点（中心点）の周りになめらかに分布するというものであり，その分布の範囲が粒子の相互作用領域となる．この質量の広がりを通して，粒子系と流体系とで物理量の対応関係が付けられ，流体系の Navier-Stokes 方程式から粒子系の運動方程式が導かれる．最初にこの手法が現れたのは天体物理学分野でのシミュレーショ
ンに関する Gingold と Monaghan (1977)[1] の論文あるいは Lucy (1977)[2] の論文と言われている．その後この手法は改良され様々な現象に応用されている．解説書としては [3]，[5] などがある．
SPH における基礎方程式の導出の過程は厳密なものではなく，かなりの程度発見法的（heuristic）であり，厳密志向の研究者には受け入れがたいかもしれない．しかし，得られた方程式はニュートン力学の運動方程式であり，容易にコンピュータ・シミュレーションが実行でき，その結果は非常によく流体の運動を再現している．特に，界面を持つ流体の動きや容器からの流体の流失の様子，熱対流など，従来の時空での偏微分方程式の有限要素法では扱いが困難であった流体運動が極めて簡単に扱える点が特筆に値する．
なお，SPH の運動方程式には影響半径や粒子数などの未決定パラメータがあり，それらを適切に決定すためにはシミュレーションを実行して結果を見るというのが現在のところ取り得る最良の方法の一つである．そのために，インタラクティブにシミュレーションが実行できる高速プログラムが望まれる．
以下，第２節，第３節では SPH の基本方針を詳述する．その際 SPH で使われる定義，仮定あるいは要請を明瞭にし論理的一貫性を持たせるべく工夫した．第４節では得られた粒子
系運動方程式をコンピュータプログラムに載せるための方針が示される．特に，相互作用の計算での計算量を O(N 2 ) から O(N ) に落とすための方法について解説する（N = 粒子数）．
第５節では熱伝導の方程式を組み込み，熱対流などが扱えるようにするための補足が加えられる．第６節ではこの論文で解説した運動方程式を Java で実装したインタラクティブなプログラムがダウンロードできる URL を示し，その使い方について簡単に解説している．第７節では若干の補足を述べる．

2 SPH の基本方針■密度場 体積 V ，密度 ρ0 の流体を N 個の粒子の集まりとして表現する．各粒子
i = 0, 1, . . . , N − 1 は質量 mi ，位置 r i = (xi , yi , zi )，速度 v i = (vxi , vyi , vzi ) を持ち，質点系の力学に従う．まず流体の質量とそれを記述する粒子系の質量は同じであるべきだから，

m0 + m1 + · · · + mN −1 = V ρ02各粒子 i の質量 mi は粒子の中心位置 r i の回りに一定の分布関数 w(r − r i ) で分布していると考える．すると流体の密度分布 ρ(r) は∑ρ(r) =w(r − r j )mjjで与えられる．ここに分布関数 w(r) はカーネル（kernel）関数と呼ばれ，原点周りに回転
対称，台は有界，２階連続微分可能であり，かつ∫w(r) ≥ 0,

w(r)d3 r = 1(規格化条件)R3を満たす．ここに d3 r は体積素すなわち d3 r = dxdydz を表す．また R3 は 3 次元の実空間全体を表わす．w(r) の回転対称性より，w(r) の台は球体である．この球の半径 h は影
響半径（influence radius）と呼ばれる．密度分布 ρ(r) を全空間に渡って積分すると∫
3

ρ(r)d r =

∑

R3

∫
mj
R3

j

w(r − r j )d3 r =∑mj = V ρ0jとなって流体の全質量に一致する．カーネル関数の回転対称性により w(r) = W (|r|) と表される．W (r) もカーネル関数と呼ばれる．W (r) は 2 階連続微分可能でW (r) ≥ 0,

W (−r) = W (r),∫ ∞4π

W (r)r2 dr = 1

0

を満たす．
なお，考えている流体は密度が一定値に指定されているのであるから，当然非圧縮性流体である．したがって，ρ(r) がこの流体の密度場を表すとするなら，ρ(r) = ρ0 でなければならないが，SPH 近似の立場では，これは要請しない．疑似剛体の場合と同様に，粒子運動の結果として，密度場に変化が生じそれによって局所的力の作用が系全体に広がっていくというイメージである．SPH 系が表す流体は疑似的非圧縮性となり，運動の結果として，

ρ(r) ≈ ρ0 となることが想定されている．

■速度場と加速度場と連続の式 粒子 i ごとの速度 v i が与えられたとき，それらによって生成される流体系の速度場 v(r) は運動量密度の場 ρ(r)v(r) として∑
ρ(r)v(r) =w(r − r j )mj v jjと定義される．また流体系の加速度場 v̇(r) は力密度の場 ρ(r)v̇(r) として

ρ(r)v̇(r) =∑w(r − r j )mj v̇ jjと定義される．ここで定義した速度場と密度場は連続の式を満たす．

∂t ρ(r) =∑j=−∂t w(r − r j )mj∑∇w(r − r j ) · v j mj

j= −∇ · (ρ(r)v(r))∂∂∂∂，∇ = (∂x , ∂y , ∂z ) = ( ∂x
, ∂y
, ∂z
) である．
ただし，∂t = ∂t

なお，後で示されるように加速度場は速度場の Lagrange 微分に対応する．

3

■物理量の流体系ー粒子系対応と整合性の要請 SPH 粒子は統計力学的粒子ではなく流体の小さな塊のようなものである．したがって，各粒子は質量と位置と速度以外に流体で定義されている巨視的物理量に対応した物理量を持つことになる．そこで両者の関係を決めておく必要がある．流体系で定義される物理量の場 f (r) に対応する粒子系での物理量 fi は

f (r) ≈∑w(r − r j )j

mj f jρj

ただし ρj = ρ(r j )を満たすものとする．これは粒子 i の流体場への寄与の割合は w(r − r i )mi /ρ(r i ) であることを意味する．さらに，次の整合性の要請を置く．

f (r i ) ≈ fi
この要請は結局f (r i ) ≈∑w(r i − r j )mj f (r j )ρjw(r i − r j )

mj fj
ρj

w(r i − r j )

mj
1
ρj

j

または

fi ≈

∑
j

が成り立つことを意味する．特に，

1≈

∑
j

が成立する．また，流体の密度場 ρ(r)，運動量場 ρ(r)v(r) に対して上記の対応関係と整合性の要請を使うと

ρ(r) ≈∑w(r − r j )

j

=

∑

mj ρ jρjw(r − r j )mjjρ(r)v(r) ≈∑w(r − r j )

j

=

∑

mj ρ j v jρjw(r − r j )mj v jj

最初の定義と矛盾しない結果が得られる．加速度場についても同様である．

■カーネル関数の具体例 カーネル関数 w(r) = W (|r|) としては計算のしやすさから次の形が採用されることが多い．

{
W (r) =

C(h2 − r2 )30

C は規格化条件から決定される．3 次元系の場合
C=

315
64πh9

4

if r2 < h2otherwise

となる．ここで，以下で必要になる ∇w(r), ∇ · ∇w(r) を計算して置く．r = |r|, r 2 < h2の範囲でr
= −6C(h2 − r2 )2 r
r
1
∇ · ∇w(r) = W ′′ (r) + 2W ′ (r) = 6C(h2 − r2 )(7r2 − 3h2 )r∇w(r) = W ′ (r)なお，r 2 ≥ h2 の範囲では，これらはいずれも 0 である．
以下このカーネル関数について話を展開する．

■カーネル関数の台の半径 h の決め方 与えられているパラメータは流体の体積 V ，密度 ρ0 ，SPH 粒子の個数 N である．これにさらに，相互作用可能な近傍粒子の個数すなわ
ちカーネル関数の台の中に入り得る粒子の個数 nh を指定する．まず粒子は体積 V 中に均一に分布しているとして，粒子一つ当たりの体積は V0 = V /N であり，半径 h の球体に入り得る粒子の個数は 34 πh3 /V0 である．したがって，
4πh3
= nh
3V0
が成り立つ．これを解いて，

(
h=

3nh
4π

)1
3

d0

1

ただし，d0 = V03 である．粒子を直方体として表せば，その一辺の長さが d0 である．コン
ピュータ・グラフィックでは粒子は通常直径 d0 の球として表されるが，SPH としては粒子の形に意味はない．なお，d0 を用いれば，粒子の質量はm = ρ0 d30と表される．
さて，近傍粒子の個数 nh はどのようにして決定したらよいであろうか．これには特に指針があるわけではなく，実際にシミュレーションを行って，うまくいきそうな値を採用するということになる．これを適切に実行するには，プログラムはインタラクティブに後から様々なパラメータが変更できるようになっていなければならない．ちなみに私の作ったプロ
グラムでは nh の初期値は３次元系では nh = 12，２次元系では nh = 6.7 としてある．これで大体うまくいくようである．

3 基本方程式の導出
前節の仮定とカーネル関数の具体的形に基づいて，Navier-Stokes 方程式から，SPH の基本方程式を導く．

■粒子系の運動方程式 二つの SPH 粒子 i, j は w(r j − r i ) > 0 で決まる近傍にあるとき，相互に力を及ぼしあい，通常のニュートン力学の法則にしたがう．その力を F ij とすると，重力を受ける地上での運動方程式はṙ i = v i∑mi v̇ i =F ij + mi gj

5

となる．ここで，g = (0, −g, 0) は y 軸マイナス方向の重力を表す．g = 9.8 [m/s2 ] は地
上での重力加速度である．粒子間力 F ij は作用―反作用の条件を満たす．

F ij = −F ji
この Fij を流体の Navier-Stokes 方程式に関連付けて決定する．

■Navier-Stokes 方程式 流体は粘性率 µ が時間空間的に一定で，疑似的非圧縮性を想定す
る．疑似的非圧縮性とは密度場 ρ は変化しうるが，その量は小さく ρ∇·v = ∂t ρ+v·∇ρ ≈ 0とおけるものを意味する．このとき，Navier-Stokes 方程式は運動量の保存式として∂t (ρv) + ∇ · (ρvv) = −∇p + µ∇ · ∇v + ρgで与えられる．左辺は連続の式を使って∂t (ρv) + ∇ · (ρvv) = ρ∂t v + ρv · ∇v + (∂t ρ + ∇ · (ρv)) v = ρDt vと表すこともできる．ただし Dt = ∂t + v · ∇ は Lagrange 微分である．
以下各項の SPH 表現を与える．ここでは，次の記法を使う．
［流体系の場の量］→［対応する SPH 表現］■時間微分：Navier-Stokes 方程式の左辺∑
w(r − r j )v j mj∂t (ρv) → ∂t
j

=−

∑

∇w(r − r j ) · v j v j mj +j∇ · (ρvv) → ∇ ·∑∑w(r − r j )v̇ j mjjw(r − r j )v j v j mj =∑j∇(r − r j ) · v j v j mj

jより

∂t (ρv) + ∇ · (ρvv) →∑w(r − r j )v̇ j mj = ρv̇

j

■圧力項
( )
( )
( )
∇p
p
1
p
p
=∇
− p∇
=∇
+ 2 ∇ρ
ρ
ρ
ρ
ρ
ρ
( )
∑
p
pj
∇
→
∇w(r − r j )mj 2
ρ
ρj
j
∑
∇ρ →
∇w(r − r j )mj
j

より

∑
∇p
→
∇w(r − r j )mj
ρ
j

6

(

pj
p
+ 2
2
ρj
ρ

)

なお，単純に∑∇p→∇w(r − r j )mj
ρ
j

(

pj
ρρj

)

ととることも可能であるが，粒子間力の反対称性の要請を満たさないので，これは採用し
ない．

■粘性項
∇ · ∇v →∑∇ · ∇w(r − r j )mj

j
対称化

≈

∑

vj
ρj

∇ · ∇w(r − r j )mjjvj − vρj

なお，対称化での追加項は前節で述べた整合性の要請により ≈ 0 である．すなわち∑∇ · ∇w(r − r j )mjj= v∇ · ∇∑v
ρjw(r − r j )jmj≈ v∇ · ∇1 = 0ρj■以上まとめて， 上記の Navier-Stokes 方程式は
(
)
∑
p
pj
ρv̇ = − ρ∇w(r − r j )mj
+ 2
ρ2j
ρ
j
∑
vj − v+µ∇ · ∇w(r − r j )mj
ρj
j
+ ρg
と変換される．

■r = r i として整合性の要請を使うと
(
)
∑
pj
pi
v̇ i = −∇w(r i − r j )mj
+ 2
ρ2j
ρi
j
∑
vj − vi+µ∇ · ∇w(r i − r j )mj
ρj ρi
j
+g
これが SPH の基本方程式である．粒子系の運動方程式として書けば，

mi v̇ i =∑(F pressure+ F viscosity) + mi g
ij
ij

j
viscosity
F pressureと F ijは，それぞれ，粒子間の圧力勾配力と粘性力であり，
ij
)
(
pi
pj
pressure+ 2
F ij= −mi mj ∇w(r ij )
ρ2iρj= µmi mj ∇ · ∇w(r ij )F viscosityiv ji
ρi ρ jと与えられ，反対称性の要請を満たす．ただし，r ij = r i − r j , v ij = v i − v j である．
これら粒子間力につきカーネル関数の具体例に基づき最終的表現を与えることができる．

7

■圧力勾配力) (6C pi + pj mi mj (h2 − |r ij |2 )2 r ij
ρ2i
ρ2j
F pressure
=
ij

0

(|r ij | < h)(|r ij | ≥ h)

C は前節カーネル関数の具体例で示したもの．各粒子 i の圧力 pi は状態方程式により別に
定める．

■粘性力


6Cµ v ji mi mj (h2 − |r ij |2 )(7|r ij |2 − 3h2 )viscosityρ i ρjF ij=
0(|r ij | < h)(|r ij | ≥ h)

ただし µ は粘性係数．

■圧力 各粒子の圧力 pi には状態方程式を仮定する．
{
γ(ρi − ρ0 ) if (ρi − ρ0 ) > 0
pi = k0 (ρi − ρ0 ) +0otherwisek0 ≥ 0, γ ≥ 0 は定数．k0 は気体定数に相当し，γ は剛性を表す．ここでは，疑似的非圧縮性流体を想定しているので，通常は k0 = 0 とし，γ には比較的大きな値が設定される．なお，ρi = ρ(r i ) は前節の密度場の定義式より決定される．この状態方程式は重要である．これにより，ρi ≈ ρ0 がほぼ成り立ち，各粒子は互いに大きく重なり合うことが阻止され，一定の容積のある流体が生成される．

4 運動方程式の差分近似とコンピュータプログラム■運動方程式 [粒子 i に働く力]/mi を f i とすると，運動方程式は
dr i
= vi
dt
dv i
= f i ((r · ), (v · ))dtと書ける．ただし(r · ) = (r 0 , r 1 , . . . , r N −1 )(v · ) = (v 0 , v 1 , . . . , v N −1 )■積分して
r i (t + ∆t) = r i (t) + v̄ i (t)∆t
ただし

1v̄ i (t) =∆t
これはv̄ i (t + ∆t) = v̄ i (t) +1∆t∫ t+∆tv i (s)dst∫ t+∆t∫ s+∆t
ds

t

s

を満たす．

8

ds′ f i ((r · (s′ )), (v · (s′ )))■被積分項の中間点近似 ここで，被積分項 f i ((r i (s′ )), (v i (s′ ))), t < s′ < t + 2∆t, を
f i ((r · (s′ )), (v · (s′ ))) ≈ f i ((r · (t + ∆t)), (v̄ · (t)))と近似する．結局r i (t + ∆t) = r i (t) + v̄ i (t)∆tv̄ i (t + ∆t) = v̄ i (t) + f i ((r · (t + ∆t)), (v̄ · (t))∆tを得る．ここでさらに

ui = v̄ i ∆t,f̃ i ((r · ), (u· )) = f i ((r i ), (

u·
))∆t2
∆t

とおくと，

r i (t + ∆t) = r i (t) + ui (t)ui (t + ∆t) = ui (t) + f˜i ((r · (t + ∆t)), (u· (t)))■C 言語風コンピュータプログラムの形で書くと
loop{
for(i = 1, ..., N ) r i += ui ;for(i = 1, ..., N ) f̃ i = f̃ i ((r · ), (u· ));for(i = 1, ..., N ) ui += f̃ i ;}となる．適当な初期条件の下この 3 行の計算を繰り返すことにより，コンピュータ内で粒子系の運動をシミュレートすることができる．

■f̃ i の最終形 同一質量粒子 m = mj として，f̃ i の最終形をまとめると)∑′ ( p i
pjf̃ i ((r · ), (u· )) = Cpress+(h2 − |r ij |2 )2 r ij
2
2
ρ
ρ
i
j
j
∑′ uji 2+ Cvisco(h − |r ij |2 )(7|r ij |2 − 3h2 )ρi ρjj+ (0, −Cgrav , 0)
各係数は945m∆t232πh9945mµ∆t
Cvisco = 6Cmµ∆t =32πh92Cgrav = g∆t

Cpress = 6Cm∆t2 =である．ここで和の記号∑′jは |r ij | < h を満たす j についての和

∑′
j

∑

=

j:|r ij |<h

を表す．

9

■容器壁との相互作用について 上述のプログラム方式により流体系の SPH によるシ
ミュレーションはできるのであるが，具体例に適用するには流体を入れる容器の壁との相互作用を記述する必要がある．壁の取り扱いで最も簡単なのは，壁を SPH 粒子の集まりとして記述する方法である．この場合，壁粒子は流体粒子の運動方程式の力の項に流体粒子と同じ形式で参加することになる．その際，壁粒子の質量及び粘性係数が流体粒子のそれらと異なり得るが，境界条件を後から付け加える必要はない．そして，壁粒子としてのもっと重要な点は，壁粒子の位置は固定されるあるいは外部的に制御されるという点である．もちろん，壁粒子が流体粒子からの作用を受けて集団的な運動するというモデルも作り得るが，今回は考えない．

■空間セル分割法による計算量の低減化 上で得られた運動方程式のプログラムを実行する場合，最も時間のかかるのは粒子間の力の計算である．粒子対 (i, j) は N 2 個あるので，loop の各ステップで，オーダー O(N 2 ) の計算量となる．各粒子ごとの力の計算に現れる和演算は |r ij | < h を満たす j についての和であるが，このような j の個数は総粒子数
N とは無関係な一定数に抑えられるので和演算そのものの計算量は O(1) である．問題は|r ij | < h を満たすか否かを判定する計算である．この計算をそのまま実行すると，これには粒子対 (i, j) ごとの計算が必要になり，O(N 2 ) の計算量となる．
実はこの計算量は O(N ) まで低減することができる．粒子の存在する空間を一片の長さ
h の立方体のセルに分割し，各粒子の位置が決まるごとにその粒子を粒子が所属するセルに登録する．すると粒子 i に対し |r ij | < h をみたす粒子 j は i と同じセルあるいはそれに隣接する 26 個のセルに登録されている粒子に含まれる．これらのセル内の粒子の総数は全粒子数 N とは無関係に一定数に抑えられるので，|r ij | < h の判定にかかる計算量は O(1) となる．
各セルは空間構造を保つように 3 次元整数 Z3 で番号付けすることにより，粒子の位置r = (x, y, z) に対して所属すべきセルの番号が，([x/h], [y/h], [z/h]) と簡単に決定できる．
よって粒子のセル登録に要する計算量はもちろん N に無関係で，O(1) である．
上記 loop の各ステップではこれら O(1) の計算を N 回実行することになるので，計算量は O(N ) となる．

■粒子系の初期設定 例えば，穴の開いた容器から流体が流失する様子を SPH でシミュレートする場合，初期条件として粒子を容器中に配置することになる．通常は規則的に並べて配置する．この状態は粒子間の相互作用によって得られる自然状態ではない．したがって，初期配置からいきなりシミュレーションを開始してはならない．まずは，穴を閉じた状
態でしばらく置いた後，おもむろに穴をあけるという操作が必要になる．このような操作を実行するためにもプログラムはインタラクティブでなければならない．

5 熱伝導
疑似的非圧縮性流体の近似が使える状況で単位質量あたりの比熱 c が一定であるとして，
内部エネルギー（熱）は cρT と表される．ここに T は絶対温度であるが，以下の式は温度の差のみに関係するので，基準点は氷点あるいは室温など任意に選ぶことこができる．このとき．熱伝導方程式は∂t (cρT ) + ∇ · (cρT v) = κ∇ · ∇T

10

と与えられる．ここで，熱伝導率 κ は一定とする．連続の式を使うと，

ρDt T =

κ∇ · ∇Tcとも表される．

■熱伝導方程式の SPH 表現 温度の場 T (r) の SPH 表現は∑mj
T (r) →w(r − r j )Tjρjjとなる．Ti は SPH 粒子 i の温度であり，整合性の条件 Ti ≈ T (r) が満たされているものとする．粘性項の場合と同様にして，

ρi Ṫi =κ∑Tj − Ti∇ · ∇W (rij )mjc jρjを得る．∇ · ∇W (rij ) の具体的式を代入し，さらに mj = m としてṪi = 6Cmκ ∑′ 2Tj − Ti22(h − rij)(7rij− 3h2 )c j
ρi ρj

■SPH 熱伝導方程式の差分化 コンピュータでの計算に載せるため，上記の 1 階の常微分方程式を差分方程式に直し，さらに C 言語風の代入式に書き換えると
Ti += Cthermo∑′22(h2 − rij)(7rij− 3h2 )jTj − Ti
;
ρi ρj

係数は

κ
945mκ∆t
Cthermo = 6Cm ∆t =c
32πh9 c

(C =

315
)
64πh9

この代入式を先の loop 文の中に組み込めば，系の各部の温度変化も含めたシミュレーショ
ンが可能になる．

loop{
for(i = 1, ..., N ) r i += ui ;for(i = 1, ..., N ) {f̃ i = f̃ i ((r · ), (u· ), (T· )); fiT = fiT ((r · ), (T· )); }for(i = 1, ..., N ) {ui += f̃ i ; Ti += fiT ; }}ここで，f̃ i ((r · ), (u· ), (T· )) は元の f̃ i ((r · ), (u· )) の重力項に下記の変更を加えたものであり，また fiT ((r · ), (T· )) は上記温度変化の代入式の右辺を表す．

■熱伝導の運動への効果 熱伝導のある流体においては重力環境下で熱対流が起こる．
熱対流は系全体でみると，温度による密度分布の変化により生じた位置エネルギーの不安定状態を解消させるべく引き起こされる運動である．粒子法的観点から見ると，温度によって異なる重さをもつ粒子の集まりの運動となるが，質量を温度依存にすると質量保存則に抵触する上に粒子の質量はカーネル関数の台の半径および基準密度に密接に関係しているため，安易に粒子の質量だけ変更すると粒子法の基本構造を破壊しかねない．そこで，質量を慣性

11質量と重力質量とに分けて考え，粒子が軽くなるという効果は重力質量だけに現れるとす
る．こうすると，粒子法の基本構造はそのままにして，粒子 i の運動方程式に現れる重力項

g を g/(1 + bTi ) に変更するだけでよいことになる．ここで b は適当な正の係数であり，膨張率に関係する量である．また温度 Ti は適宜設定した基準点から測ったものである．
疑似的非圧縮性流体で熱対流を扱う場合，通常は Boussinesq 近似が使われるようである．これは新たに，重力項 −bT g を付け加えるというものである．SPH の運動方程式では
重力項 g を (1 − bTi )g と変更することに相当する．bTi が小さい場合には両者ほぼ同じである．

■界面からの熱放出 容器の底を温めて流体に永続的対流を起こすには流体の界面から熱が空間中に逃げるという効果を取り入れなければならない．容器に入れられた SPH 粒子系
では界面は重力の効果により自然に出来上がる．界面として特別に境界を設ける必要はない
が，粒子が界面に近いところにいるか否かを判定する必要が出てくる．界面は近傍粒子数が平均より少ないか否かで判定できる．空間への熱放出量は近傍粒子数の少なさに比例して決定される．飛び散った粒子などには近傍粒子がないので，急速に冷やされる．

6 プログラムの説明
本論文で解説した SPH の基礎方程式を具体的に Java で実装したプログラムが以下のサ
イトからダウンロードできる．

http://www.is.kochi-u.ac.jp/~nakagomi/public.html
ここには，2 次元，3 次元の Java 実行ファイル（jar ファイル）とソースファイル（java
ファイル）が収められている．以下使い方について簡単に説明する．

■2 次元 SPH プログラム
ファイル SPH2D.jar をクリックする，あるいは，コマンドラインで

java -jar sph2D.jar
を実行すると左図のウィンドウが現れる．上部にメニューボタンがある．それらをマウスで押すことによってさまざまな操作ができる．[END] プログラムの終了，[Reset] 運動
設定の初期化，[Start/Stop] 運動の開始と停止，[Remove] 上部容器の削除，[Break] 仕
切り壁 (ダム) の崩壊，[HeatB] 底からの加熱，[HeatL] 左壁の冷却，[HeatR] 右壁の冷却，[Swing] 容器の左右運動の開始と停止，[Add] 液体の追加，[Show] データの表示 (コマンドラインから立ち上げた場合のみ)，[Control] 諸パラメータの変更を行うためのパネルが立ち上がる．またマウスでドラッグして図の移動，ホイールボタンを回して拡大縮小ができる．
以下操作ボタンを押して得られる運動の様子を示す．

12

[Start/Stop][Open]−−−−−−−→[Remove]−−−−→[Break]−−−−−→−−−−→

スタートからダムの崩壊まで．

[Add]→[Open][Remove]−−−−−−−−→[HeatB]−−−−−→−−−−→→−

液体の追加から底加熱，対流発生まで．

[Swing]

−−−−→

→
−

→
−

容器の横揺れ運動．

■3 次元 SPH プログラム
ファイル SPH3D.jar をクリックする，あるいは，コマンドラインで

java -jar sph3D.jar
を実行すると左図のウィンドウが現れる．2次元系の場合と同様に上部にメニューボタン
があり，同様な操作が可能であるが，下部にもメニューボックスがあり様々な形の容器が選択できるようになっている．同じ形でも機能（運動の種類，温度設定，変形の可能性）が異なる容器は別の名前が与えられている．機能の違いはメニューボタンおよび [Control]のパネルに反映される．また，マウス操作では，右ボタンと左ボタンで機能が分かれる．
左ボタンは表示位置の移動，右ボタンは視点の 3 次元回転を行う．

以下若干の操作例を示す（操作ボタンの表示は省略）．

スタートから液体の追加を経てダムの崩壊まで．

13

トーラス型容器に液体を注入（5 回）し，下半分を加熱して，一方向への回転が起こるまで，最後は視点の回転あり．

7 終わりに■MPS 法について 本論文では取り上げなかったが，もう一つの粒子法として，MPS 法（Moving Particle Semi-implicit method[4]）がある．これは，完全な非圧縮性流体の記述に用いられるもので，圧力の決定に状態方程式を使わず，運動方程式を解くときに密度が非圧縮性の条件を満たすように圧力を逐次的に変えていくという方法を取る．非圧縮性流体の数値解法としては優れた方法である．興味のある方は [4]，[5] などを参照していただきたい．

■計算量について 本論文で述べた空間のセル分割法により相互作用計算の計算量がO(N ) まで落とせるのであるが，力の計算は各粒子ごとに独立に行えるので，並列度 K のシステムが使えるなら，各並列分枝の計算量は O(N/K) となり，K が大きい場合大幅な計
算時間の短縮が期待される．GPU など並列計算システムを使った研究が望まれる．
■SPH プログラムについて 本論文で紹介したプログラムでは，様々な流体現象（例えば，揺れ動く容器の中の流体の動き，水柱の崩壊，容器の穴からの流体の流出，熱対流など）
がインタラクティブにシミュレートできるようになっているが，単一成分の流体を想定している．今後，多成分流体や流体中に浮遊する物体の動きをシミュレートするプログラムに拡張する予定である．また，３次元の SPH プログラムでは，様々な形状の容器が選択できるようになっているが，これら容器は Wall クラスを継承して作られている．Wall クラスを継承したクラスが jar ファイル内，あるいはカレントディレクトリ内にあれば，自動的にメ
ニューボックスに載るようになっているので，新しい容器の追加は簡単である．興味のある方は容器のソースファイル（Dam.java, Bowl.java, Bucket.java, Torus.java など）
を参照して独自容器の作成を試みられたい．

参考文献[1] R. A. Gingold, J. J. Monaghan (1977), Smoothed particle hydrodynamics: theory
and application to non-spherical stars , Mon. Not. R. Astron. Soc. 181: 375-389.
[2] L. B. Lucy (1977), A numerical approach to the testing of the fission hypothesis ,Astron. J. 82: 1013-1024.
[3] J. J. Monaghan (1988), “An Introduction to SPH”, Comput. Phys. Comm. 48: 8996.
[4] 越塚誠一，日本計算工学会（編） (2005)，「粒子法」，丸善株式会社，ISBN 978-4-62107522-7．
[5] 後藤仁志（2018），「粒子法」，森北出版，ISBN 978-4-627-92231-0．

14