多段格子エルミート法を用いた混合粒径砂礫粒子の沈降シミュレー
ションと空隙率評価
Sedimentary Simulations of Polydisperse Granular Particles with Muti-Size-Mesh Hermite
Scheme and Porosity Evaluation

曽田康秀∗ ・渡邊明英∗∗Yasuhide SOTA and Akihide WATANABE∗ 理博∗∗ 正会員株式会社 東京建設コンサルタント（ 〒 170-0004 東京都豊島区北大塚 1 丁目 15 番 6 号）
工博 株式会社 東京建設コンサルタント（ 〒 170-0004 東京都豊島区北大塚 1 丁目 15 番 6 号）

We propose the Multi-Size-Mesh Hermite Scheme to perform the N-body simulations with shortrange interactions and apply it to the sedimentary simulations of granular materials in fluid. The
Multi-Size-Mesh Hermite Scheme is eﬀective to the simulations with a wide range of scales both
spatially and temporarily, because both the time steps and the mesh sizes of each particle are
discretized as the power of 2. As for the sedimentary simulations of monodisperse spheres, wefound that the porosity increases as the deposition intensity is enhanced. As for the sedimentary
simulations of polydisperse spheres, on the other hand, the porosity depends not only on theintensity but also on the grain-size distributions. In fact, the porosity decreases as the fraction ofsmaller particles increases for the simulations with strong deposition intensity, while it becomesalmost constant for those with weak deposition intensity. Especially, in the latter, segregation
was found to occur in the deposit, i.e., the smaller particles sink from the top to the bottom in
deposit, which suggests that the segregation plays a key role in determining the packing statesin the sedimentary simulations of polydisperse spheres.
Key Words : 粒状体, 離散要素モデル、多段格子、エルミート法

1.

程に依存すること6) が明らかにされている. すなわち,河床材料の空隙率は前述した粒度分布変化以外にも, 沈
降粒子の濃度や沈降速度などの物理パラメータに依存して変化する.

はじめに

河床変動計算では空隙率は一定の値が用いられていることが多いが, 現実の河川では掃流砂の運動状態や粒度分布の変化などによって時間的空間的に変化する. したがってこのような粒度分布の変化と空隙率の関係を明らかにしておくことが重要であると考えられる. 実
際,Sulaiman らは混合粒径球状粒子をランダムに付着させていく粒子付着法を用いて, 河床面における空隙率と粒度分布の関係について解析している1),2) .

一方, 混合粒径分布のランダム充填と空隙率の関係に関しては, 粒子付着法や充填した粒子を圧縮していく粒子圧縮法8) などが適用されているだけで, 粒子沈降の数
値解析はほとんどなされていない. これは, 単一粒度分
布に比べて, 連続粒度分布の数値解析では, 粒子分布の広がりによる粒子数の増加などの点で, 数値計算のコストが非常に大きくなるためであると考えられる. このような点を踏まえて, 本研究では, 連続粒度分布に従う粒子の沈降シミュレーションを高速に行うための計算ス
キームとして, 多段格子エルミート法を提案した9) . このスキームは分子動力学の分野で適用されているセル
インデックス法10),11) を多段格子に拡張し, それと天体力学の分野で適用されているエルミート法12),13),14) とを組み合わせたものである.

単一粒径分布のランダム充填と空隙率の関係については, 粉体工学の分野で精力的に研究が進んでおり, 粒子付着法や粒子間接触力をバネ・ダッシュポット系として取り扱うことによって多体粒子系の運動を同時に解く DE(Discrete Element) 法3) によって堆積層の空隙率が調べられている4),5),6) . 粒子付着法では粒子沈降の物理過程を考慮せずに粒子が積み上げられるのに対し,DE 法では多体粒子間接触や流体場からの影響, 回転の効果などの物理過程を考慮した堆積層生成過程が実現できる. これらの解析によって, 粉粒体の堆積層がブリッジのでき方や接触角の違いによって Random Close
packing(RCP) と Random Loose packing(RLP) に分類されること4),5),7) や, 粒度分布が固定されたとしても空隙率は一意的に決まらず, 粒子沈降および堆積の物理過

これまで河川工学の分野において, DE 法はアーマーコートなどの分級過程の解析15) や粒子法を用いた固液
混相流の解析にも用いられている 16) . しかしながら,固相が混合粒径砂礫から成るような系では, マルチス
ケールリンクなどの工夫が必要であることが示されている 17) . 一方我々の提案した多段格子エルミート法で

1

は, 個々の粒子ごとに異なる可変固有時間が与えられるため, 落下粒子と堆積層に落ちついた静止粒子とで異なる時間刻みを与えることができる. また各粒子の粒径ごとに異なるサイズの格子を用意し, 粒子間接触の判定を効率良く行うことを可能にしている. したがって多段格子エルミート法は, 河川中の混合粒径砂礫の沈降過程のように, 粒子の落下時間と粒子接触時間の大きさや各粒子の粒径の間に大きな隔たりがあるような場合には特に有効であると考えられる. 本研究では、多段格子エルミート法を静止流体中での混合粒径砂礫の沈降シミュ
レーションに適用し, 沈降過程によって得られる堆積層の空隙率と粒度分布の間のマップについて考察する.

度ω~ i で回転する粒子では, これらの力に加えて
1
2(4)
f~Li = CL ρSi |~ui | ~eL ,2で与えられる揚力 f~Li が働く. ただし ~eL は (~ωi × ~ui ) 方
向の単位ベクトル, CL は揚力係数である21) .
これらの係数のうち, 慣性係数 CM は静止流体中の球状粒子に対しては 0.5 とした22) . CD と CL は νは流体の動粘性係数として, 一般に粒子レイノルズ数Re ≡ |~ui |di /ν の関数となる. 非回転球に関しては,Re  103 で CD は 0.4 に Re  1 で 24/Re に収束す
る. Re < 103 で実験結果に従う CD の関数形としては,
)
24 (
CD =
1 + 0.15Re0.687(5)Reが知られている23) .

静止流体中での砂礫粒子沈降シミュレーションの基礎方程式

2.

角速度 ω
~ i で回転する球に対しては, CD と CL は Reおよび Ω∗ ≡ di |~ωi |/(2 |~ui |). の関数となる.
Kuroseと Komori は 1 ≤ Re ≤ 500 ,0 ≤ Ω∗ ≤ 0.25 の範囲内
で,CD と CL が以下の関数形で近似できることを数値的に示している21) .
)(4
CD (Re , Ω∗ ) = 1 + η1 (Re ) Ω0.5
∗ + η2 (Re ) Ω∗ CD (Re , 0)

一般に流体中での多体粒子系の運動に対しては、流体は粒子系から相互作用を受ける. この粒子-流体間相互作用は、流体の基礎式を局所平均化して導かれ、空間の空隙率が小さい場合は Ergun の式18) で記述されることが知られている19) . 本研究では、流体中の多体粒子計算の第一段階として、静止流体中での粒子の運動のみを考え、流体は粒子運動の影響を受けない状態を仮定している.

CL (Re , Ω∗ ) = η3 (Re ) Ω∗ .
以下の数値計算では,CD (Re , 0) に対して 式 (5),0 <Ω∗ < 0.5 の範囲で式 (6) ,Ω∗ ≥ 0.5 の範囲で関係式

速度 ~ui , 角速度 ω~ i で静止流体中を動く質量 mi の単体球状粒子に働く力は (a) 重力 f~Gi ≡ mi~g ,(b) 流体力
f~W i および (c) 粒子間接触力に分類され, , 球状粒子の運動方程式は)∑(
d~ui
f~c,ij + f~d,ij + f~W i + f~Gi , (1)=midt

CD (Re , Ω∗ ) = CD (Re , 0.5)
CL (Re , Ω∗ ) = CL (Re , 0.5)

Ii

d~
ωi
=
dt

(

~i
~τW i = −πρνd3i ω)~ ij × f~c,ij + f~d,ij + ~τW i ,R(8)で与えられる. これまでに,Re が大きい場合に球状粒子が流体から受けるトルクの一般形は知られていない.
従って, 以下では式 (8) を流体からのトルクとして採用
する.

(2)

j

で与えられる. ただし,f~c,ij および f~d,ij は粒子 i, j 間接
~ ij は粒子 i の中心
触力および粘性摩擦力を表す. また Rと粒子 j との接触点を結ぶベクトル, Ii は粒子 i の慣性/モーメントで, 粒径 di の球状粒子の場合, Ii ≡ mi d2i 10で表される.

2.1(7)を採用する. Re が十分小さい場合, 粒径 di , 角速度 ω~iの粒子が周囲流体から受けるトルクは, ストークス近似
により

j

∑

(6)

2.2

粒子間接触力本研究では, 砂礫粒子の粒径は 10−4 m から 0.1m の範囲を取るものとし, 電磁気力などの微小粉体に働くよ
うな力は考えないものとする. 粒径が di および dj の２つの砂礫粒子に働く接触力を接触面に対して鉛直な成分 f~n と水平な成分 f~t に分割すると, 鉛直成分 f~n球状粒子が周囲流体から受ける力およびトルク

速度 ~ui で静止流体中を動く粒径 di の球状粒子に働く流体力は, 流体の周辺粒子からの影響を考えないもの
とすると,

c,ij

c,ij

c,ij

(3)

は,2 粒子の中心を結ぶ単位ベクトルを ~nij ,2 粒子間距
離を rij として, ヘルツ則により( n )3/2nf~c,ij= Kij δij
~nij ,(9)で与えられる. ただし ρ は流体の密度, CD および CMはそれぞれ, 慣性係数を, また Vi と Si は 球の進行方//向に垂直な面への射影面積で, Si = πd2i 4, Vi = πd3i 6,と表される. また f~Bi = −Vi ρ~g は浮力である 20) . 角速¯ij ≡ (di + dj )/2 としてで与えられる. ただし d√/nδij ≡ d¯ij − rij , Kij = di dj (2 (di + dj ) D2 ) D ≡()/
3 1 − σ 2 (2E), E と σ はそれぞれ砂礫粒子のヤング率およびポアソン比である.

d~ui
1+ f~Bi ,f~W i = − CD ρSi |~ui | ~ui − CM ρVi
2
dt

2

t
t
接触力の接線成分 f~c,ij
は, 接線方向の変位 δij=)∫(~~vij · tij dt に依存する. ただし ~vij ≡ ~vj − ~vi ,~vw,ij ≡~j − w~ i で ~tij は相対速度 ~vij − ~vw,ij の方向w~j × R~i × Rを向き,~nij に垂直な単位ベクトル. 接触粒子間の滑りtを考慮した場合の接触力の接線方向成分 f~c,ijに関して
tは幾つかのモデルが提唱されている. 本研究では f~c,ij
t
は δij
を用いて,
t
n
~tij
f~c,ij= −µχ f~c,ij(/t方程式 (1),(2) の他に, これらの方程式をさらに時間微分した式を必要とする. また個々の粒子はこれらの物理量の他に, 粒子固有の時刻 tj および時間刻み ∆tj を持
つ. 従ってエルミート法は異なる時間スケールを持つ
複数の粒子が共存するような多体系の数値計算に適している. 以下では個々の時間刻みの比が２のべき乗で表されるような, ブロック時間刻みのエルミート法を用いる 14) . ブロック時間刻みエルミート法では, 各粒子の一ステップ後の時刻 tj + ∆tj が, 最小になるような粒子群 σ ∈ (σ1 , . . . , σM ) を選択することができる. この最小値 tσ + ∆tσ を tG として, 任意の粒子 j の時刻
tG での物理量の値は, 時刻 tj の周りのテイラー展開を
用いて

(10)
)3/2

で µは
とする6) . ただし χ ≡ 1 − 1 − δij δmaxt
摩擦係数,δmax は接線方向の変位 δijの最大値で, この
値を越えると接触粒子がお互いに滑り始める. 弾性理
t論から δijは()
2−σnδmax = µδij(11)
2 (1 − σ)で与えられる. 接触による粘性摩擦力の法線成分, 接線成分は
nnf~d,ij= ηij(~vij · ~nij ) ~nij
[(
)
]
t
t
~vij · ~tij ~tij − ~vw,ij
f~ = η
d,ij

ij

(12)
(13)

で与えられる. ただし mij を２粒子の換算質量, ~vw,ij ≡~j − w~ i としてw~j × R~i × R√
1/2nηij≡ α Kij mij d¯ij
√
1/2
(14)
η t ≡ β K m d¯ .
ij

ij

2

3

(tG − tj )(tG − tj ) ˙~aj +~aj + (tG − tj ) ~uj + ~rj
62
2(tG − tj ) ˙~uG,j =
~aj + (tG − tj ) ~aj + ~uj2
2(tG − tj ) ~˙~j +ω
Ωj + (tG − tj ) Ω~jω~ G,j =
2(15)
~rG,j =

で与えられる. 非選択粒子も含めた全粒子の時刻 tG での値 (15) を選択粒子に対する運動方程式 (1),(2) およ
びそれらを微分した式に代入することによって, 選択粒
子 σ ∈ (σ1 , . . . , σM ) の時刻 t = tG での物理量)(~˙ G,σ~ G,σ , ~a˙ G,σ , Ω(16)~aG,σ , Ωij ij

σ∈(σ1 ,...,σM )弾性理論に依ると, パラメータ α の β に対する比は縦方向弾性係数と横方向弾性係数の比に相当し, ポアソン√比を用いて, β = 2 (1 − σ) と記述される. この関係√
から,β は α を用いて β = α 1/2 (1 − σ) と記述される 24) . パラメータ α の値は跳ね返り係数との関係から決定される.

3.

が得られる.4 次のエルミート法では, 時刻 t = tG でこ
れらの値を取る補間多項式は以下のように与えられる.

~aσ (t) = ~a0,σ + (t − tσ ) ~a˙ 0,σ
2

3

2

3

(t − tσ ) ¨(t − tσ ) (3)~a0,σ +~a0,σ2
6
~˙ 0,σ~ σ (t) = Ω~ 0,σ + (t − tσ ) ΩΩ+多段格子エルミート法+

混合粒径砂礫粒子系の時間発展に用いる多段格子エルミート法は, 多体粒子系の時間発展を記述するエルミート法と, 混合粒径から成る粒子系の相互作用を記述する多段格子法に分けられる. 以下, これらの概要をそれぞれ述べる.

(t − tσ ) ~¨(t − tσ ) ~ (3)Ω0,σ +Ω0,σ ,
2
6(17)

ここでこれらの多項式中の係数は時刻 t = tG で, 式(17) が (16) の値を満たすように,()−6 (~a0,σ − ~aG,σ ) − ∆τσ 4~a˙ 0,σ + 2~a˙ G,σ¨0,σ =~a∆τσ2()
12 (~a0,σ − ~aG,σ ) + 6∆τσ ~a˙ 0,σ + ~a˙ G,σ(3)~a0,σ =
∆τ 3
(
) σ
(
)
~ 0,σ − Ω~ G,σ − ∆τσ 4Ω~˙ 0,σ + 2Ω
~˙ G,σ−6 Ω~¨ 0,σ =
Ω
∆τσ2
(
)
(
)
~ 0,σ − ~aG,σ + 6∆τσ Ω~˙ 0,σ + Ω~˙ G,σ12 Ω~ (3) =
Ω
0,σ
∆τσ3
(18)

3.1 エルミート法の概要
一般に, 予測子・修正子法が過去の複数の時点の情報
から現在の物理量を予測・修正するのに対し, エルミー
ト法は, 過去の一時点における高階微分の情報から現在の物理量を予測・修正する12),13) . したがって一般に4 次のエルミート法では, 各粒子 j は位置 ~rj , 速度 ~vj , 加~j ≡ ω
速度 ~aj ≡ ~v˙j , 角速度 ω~ j , 角加速度 Ω~˙j 以外に加速~˙j ）を持度および角加速度の時間微分量（それぞれ ~a˙j ,Ω
つ. またこれらの量の時間発展を記述するために, 運動
3

で与えられる. ただし ∆τσ ≡ tG − tσ . 多項式 (17) を時間積分し, 選択された粒子 σ ∈ (σ1 , . . . , σM ) の時刻 tGでの修正された物理量の値はる粒子の粒子番号()~jσ = j(σ,1) , j(σ,2) , . . . , j(σ,n) .

∆τσ4 ¨∆τσ5 (3)~a0,σ +
~a24
120 0,σ∆τσ3 ¨∆τσ4 (3)~uσ = ~uG,σ +~a0,σ +~a6
24 0,σ∆τσ3 ~¨∆τσ4 ~ (3)ω~σ = ω
~ G,σ +
Ω0,σ +Ω(19)
6
24 0,σで与えられる. 最後に, 個々の選択粒子 σ ∈ (σ1 , . . . , σM )の固有時間 tσ は tG に更新される. またこれらの粒子の固有時間刻み ∆tσ は自由に更新することが可能である. したがって, 他の粒子と接触している粒子に対しては,∆tσ をヘルツ則から導かれる以下の粒子衝突時間よりも小さく取る必要がある25) .
(
)1/5
m2σj
tcoll = 2.94
.
(20)
2 u
Kσj
σj

をもとめる. 各選択粒子の他粒子との接触判定に対しては, これらの近傍粒子に対してのみ行うことによって,接触判定の計算コストを大幅に軽減することが可能になる. 実際の沈降シミュレーションにおいては, 格子の段数 MS を上限値 Mcr に一致させることは必ずしも得策ではない. 計算領域内の粒子が非常に希薄な場合に多段格子の数を増やすと, 近傍格子点の数が近傍粒子の数よりも大きくなってしまい計算コストが大きくなってしまう. 実際,MS の最適値は, 粒度分布 f (d) および計算領域内における粒子の空間充填率 φ に依存する. 我々は, 各選択粒子の周りの平均格子点数および粒径が MS以下の平均粒子数を,f (d) および φ の汎関数として求め,MS の最適値をもとめるアルゴリズムを開発し, 真空中の 2 次元および 3 次元の砂粒子沈降シミュレーショ
ンにおいて、粒度分布が広がり小粒径粒子の割合が増加するほど、計算時間が短縮されることを示した9) .
以下の数値シミュレーションでは, この方法によってもとめた格子段数を用いて, 数値計算を実行するものとする.

~rσ = ~rG,σ +ただし mσj および uσj はそれぞれ 粒子 σ と粒子 j の換算質量および相対速度. 一方, 他粒子と接触していない粒子に対しては, 固有時間刻みの更新には注意が必要である. 特にある粒子が他の粒子に接近しているよ
うな場合は, あらかじめその粒子と接触するまでの時間を予想しておき, 接触予想時間よりも小さくなるように∆tσ を更新しなければならない. 我々はそのような非接触粒子の固有時間刻み ∆tσ の更新法を考案した9) . また, 以下では固有時間刻みは適当な値 tC を単位として,∆tσ = 2n tC のように離散化しておく. 以下の数値計算
では、tC = κc tcoll ,（κc = 0.04）とする。

3.2

4.

以下では多粒度分布に従う砂礫粒子の運動の記述に対して, 図-1 のような複数の格子サイズ (l1 , l2 , . . . , lMS )から成る多段格子を用いる. 以下の数値計算では, K 段
目格子の格子幅を(K = 1, 2, . . . , MS )

4.1

(21)

(22)

によって決定する. また粒径が lK−1 < d ≤ lK の範囲にある ~x の位置の粒子に対しては, ~x/lK の整数部分から成る K 段目格子の格子番号を割り振るものとする.
また定義から格子の段数 MS の上限 Mcr は, 与えられた粒度分布における最小粒径 dmin を用いて

lMcr +1 < dmin ≤ lMcr .

粒子流束の設定

以下の数値計算では, ランダムに粒子を計算領域の上方から一定初速度で落下させ, 計算領域内に堆積する粒子層の空隙率について考察する. 上方から落下する粒子に関しては, まず計算領域と同じ底面積 L2 から成る直方体内に, 重複しないようにランダムに粒子を生成し,それらの粒子に一定速度 ur を与えておく. その際生成された粒子群の空間充填率を φ0 とすると, 計算領域内に単位時間に降り注がれる粒子流の質量はのように与える. ただし L は沈降シミュレーションを行う計算領域の底面の一辺の長さ. また定数 K0 は dmaxを与えられた粒度分布における最大粒径として最大の格子の格子幅が dmax に一致するように, 条件

2−(K0 +1) L < dmax ≤ 2−K0 L,

沈降シミュレーションの解析条件および解析方法

以下の沈降シミュレーションでは底面積 L2 , 高さ zc
で, 側面に周期境界条件を用いた計算領域内で粒子を沈降させて堆積層を作り, その空隙率を評価する. 以下ではそのような沈降シミュレーションでの初期設定および空隙率の解析法について記述する.

多段格子法の概要

lK = 2−(K−1+K0 ) L(24)

F = ρ0 φ0 L2 ur(25)で与えられる. ただし ρ0 は砂礫粒子の質量密度である.
√
以下ではこの質量流束 F を自由落下時間 tf ≡ L/g
で規格化した無次元量 Q ≡ F tf を用いて, 生成堆積層の空隙率の質量流束依存性について調べる.

(23)

静止流体中での沈降粒子の速度は, 流体摩擦によって一定値に漸近する. 一般に粒子レイノルズ数 Re の関数
である水の抵抗係数 CD に対して, 粒子の最終速度はの条件から決定される. 多段格子法では, まず各選択粒子の近傍格子点をもとめ, 次にこれら近傍格子点に属す

4

図–1 多段格子法における粒子の配置の例 (MS = 3). 各粒子はその粒径サイズに合った段の格子に配属される. 沈降シミュ
レーションではこのような複数段の格子が計算領域全体に張り巡らされている.

f(d)

n(d)

1.0
0.8
0.6

60

S=0.0
S=0.5
S=1.0
S=2.0

S=0.0
S=0.5
S=1.0
S=2.0

50
40
30

0.4
20

0.2

10

d

d

log
3.0

2.5

2.0

log

L

1.5

3.0

(a) 粒径加積曲線 f (d)

2.5

2.0

1.5

L

(b) 個数密度分布 n(d),図–2 粒度分布 (29),Σ = 4 log 2 に対する (a) 粒径加積曲線 f (d) と (b) 個数密度分布 n(d).

以下のルビイ則で近似されることが知られている 26) .
√
√
2 36
36
ur
√
+
−
,
(26)
=
3
d
d∗
(η − 1) gdi∗

転の強さを与える無次元パラメータとして= di ω r.m.s / (2ur ) .
Ωr.m.s∗のデフォルト値としては 7.5 × 10−3を用いる. Ωr.m.s∗を用いる. また計算ボックスの幅 L は,L = 4dmax , 高さは zc = 2L とする. 他の物理量の値は表 1 に記載
する. これらの表の値を用いた場合, 最大粒径ボック
ス比 dmax /L = 0.15 および dmax = 10−2 [m] における単一粒度分布の質量流束 F/mj = 100 [個/s] はQ ≈ 20.34,φ0 ≈ 0.03 に相当する.

ここで η ≡ ρ0 /ρ は水に対する粒子の密度比で d∗ ≡/(η − 1) gd3i ν 2 . 以下では, 単一粒度分布ではその粒径からルビイ則 (26) で決まる最終速度を初期落下速度 urとして与える. また混合粒度分布では、一般に漸近速度は粒径ごとに異なるが, 一定流束を与えるために, ル
ビイ則 (26) で決まる最大粒径での最終速度を粒子の初期落下速度 ur として与えるものとする.

4.3

計算領域上方で一定速度 ur のままで鉛直落下する粒子が計算領域上面に到達するまでの時間は, 粒子の初期高さを z0j , 計算領域上面の高さを zc とすると,

tarr,j = (z0j − zc )/ur

以下では, 与えられた粒度分布における最大粒径 dmax
および最小粒径 dmin に対し, 多段格子における最大格
子サイズ l1 および最小格子サイズの下限値 lMcr をそれぞれ dmax および 2dmin に等しく取る. 従って条件
dmax = L/4 から, MS 段多段格子における各レベルの格子サイズ lK (21) は(27)で与えられる. 従って, 各粒子 j はグローバルな時刻 tG
が tG ≥ tarr,j を満たした瞬間に, 高さ zj = zc から, 速
度 ur で運動を開始するものとする.

4.2

粒度分布および多段格子サイズ

lK = 2−(K+1) L (K = 1, 2, . . . , MS ).
で与えられる.
f (d)∆d が沈降シミュレーションにおける物理パラメータ

以下の単一粒度分布の沈降シミュレーションでは, 粒子の初期回転角速度として, 各座標成分が標準偏差 ω r.m.sを取る正規分布に従ってランダムに与える. また初期回

f (d) = αS

(28)

以下では, [d, d + ∆d] の範囲の質量( d )S
L= αS dˆS(dmin ≤ d ≤ dmax )(29)

5

(a)Q = 6.78,t = 19.1tf(b) Q = 6.78,t = 57.3tf ,(c)Q = 67.8,t = 1.9tf(d) Q = 67.8,t = 7.6tf ,図–3 Talbot 型粒度分布 S = 1.0,Σ = 4 log 2 における沈降シミュレーションのスナップショット.

で表される粒度分布を考察する.dmin ,dmax と l1 , lMcrとの対応付けや f (d) の規格化条件から αS は

αS =

S+1
S+1
l1S+1 − lM
cr +1

.

る小粒径粒子の割合は, パラメータ S で特徴付けられる.
また粒径の広がりは, パラメータ Σ ≡ log (dmax /dmin )で特徴付けられる.Σ を固定して S を変えた場合の f (d)および n(d) は図-2 のようになる. 図-2 から明らかに,Σ
が一定の場合には S が小さいほど小粒径粒子の割合が大きいことが分かる. 以下では,S や Σ を変えて沈降シ
ミュレーションを行い, 堆積の空隙率のこれらのパラメータ依存性について考察する.

(30)と与えられる. またこの粒度分布に対する粒径加積曲線∫d
ˆ dˆ および規格化された個数密度分布
f (d) ≡ dmin f (d)d
/∫
l1
˜
f (d)
f (d)
( ) dd˜n (d) ≡(31)v (d)lM cr+1 v d˜

4.4

はそれぞれ

)
αS ( S+1d− lnS+1
S+1(S − 2) dS−3
n (d) = S−2
S−2
l1 − lM
cr+1
f (d) =

空隙率の評価以下の沈降シミュレーションでは, 計算領域の高さ Lまでの領域を幅 δz = 0.1L ごとに高さ z1 , z2 , z3 . . . z10に分割し, 高さ z < zj までの領域の空間充填率が 0.6を越えた段階で, 堆積層の高さが zj に達したと判断す
る. 数値計算は堆積層の高さが zj = 0.4L に達するまで行い,0.15L < z < 0.35L の領域での空隙率を評価す
る. また本研究ではモンテカルロ積分法を用い, この領

(32)
(33)

と表される. ただし v(d) は球状粒子 1 個の体積. したがって分布 (29) は指数が S+1 の Talbot 型粒度分布1),2)に相当する. また粒度分布 (29) では, 大粒径粒子に対す

6よび小角分離の割合が最大の状態を RLP と呼ぶ 4),5) .
φ-Z 面上での RCP および RLP の位置は, 統計力学的な手法によって得られることが示されている7) .

表–1 物理パラメータ値

変数名

記号

最大粒径

dmax

10−2 m粒子質量密度ρ02.65 × 103 kg/m

流体密度

ρ

1.0 × 103 kg/m

流体の動粘性係数

ν

1.0 × 10−6 m2 /s

ヤング率

E

1.0 GPa

ポアソン比

σ

0.17

単一粒度分布の場合に粒径の値および質量流束 Q を変えて, 静止流体中での沈降シミュレーションを行い,得られた堆積層の φ,Z の値を図-5 に重ねて示す. 解析方法の正当性を確認するため, 真空中での沈降シミュレー
ションを行い, Zhang らの結果6) と比較を行った. 図5 から、真空中の場合では、プロット点は Zhang らの得た曲線付近に分布することが確認できた. 一方静止流
体中でのシミュレーションでは、Q = 20.34 ではプロッ
ト点は RLP 上方に分布するのに対し（図- 5(a)）,Q を
203.4 まで上げると真空中での場合6) と異なり, プロッ
ト点の位置が図の RLP の位置に近づくことが分かった（図- 5(b)）. 実際, 図-6 より,Ωr.m.sの大きさに依らず,Q∗が大きくなるほど堆積層の空隙率は多きくなることが確認できる.

動径方向粘性係数

α

0.2

5.2

摩擦係数

µ

3

3

混合粒度分布の沈降シミュレーションと空隙率

単一粒度分布でのシミュレーションから, 空隙率の値は粒子の微量回転の大きさにほとんど依存しない. このことを踏まえ, 以下の混合粒度分布の沈降シミュレー
ションでは,Ωr.m.s= 0 として,Talbot 型粒度分布のパ∗ラメータ S および Σ と堆積層の空隙率との間の関係について図-7 に示す.

0.5

域内のランダムに置かれた点のうち, どの粒子とも交わらない点の割合として空隙率を評価する.

多段格子エルミート法を用いた沈降シミュレーショ
ンでは, 沈降粒子は計算領域の中を落下していき, 堆積層を形成していく（図-3）. このときの沈降粒子の割合
は質量流束の値によって決定される. 以下の沈降シミュ
レーションで質量流束の値はパラメータ Q で特徴付けられるが, 混合粒径の場合は分級が起こるため, 単一粒径の場合と空隙率の Q 依存性が異なることが予想される. 従って以下では堆積層の空隙率の値を単一粒径と混合粒径の場合とに分けて解析する.

図-7 から Q = 67.8(φ0 = 0.1) では, 堆積層の空隙率は S の値に依らず, 粒度の広がり Σ が大きくなるほど, 空隙率は小さくなることが分かる. しかしながらSulaimann らの対数正規型の粒度分布に対する粒子付着法による結果1),2) （図中の黒線）よりも空隙の下がりは小さくなる。これは実際の堆積過程で作られるブリッジが、粒子付着法では考慮せれないためであると考えられる。また空隙率の減少する度合いは,S が小さ
いほど, 即ち小粒径粒子の割合が大きいほど大きくなる
ことが分かる (図-7(a)). これは小粒径粒子が大粒径粒子の間の隙間を埋めるためで, 粒子付着法でも同様の結果が得られている 2),5) .

5.1 単一粒径分布の沈降シミュレーションと空隙率
一般に球状粒子の充填の特性は, 空隙率 λ（または空間充填率 φ = 1 − λ）の他に, 粒子１個が平均して何個の粒子に隣接しているかを表す数値である平均配位数
Z によって特徴付けられる. 一般にランダム充填における空間充填率と平均配位数の関係は, ブリッジの数およびブリッジを構成する粒子の角分離（ブリッジ角度）の大きさで決まる (図-4).Nolan らによれば、図-5 で示されるように単一粒度分布でもブリッジの割合が多くな
り角分離が小さくなるほど, 空隙率が上がり平均配位数が下がる4),5) . 同一粒度分布でブリッジの割合および小角分離の割合が最小の状態を RCP, ブリッジの割合お一方,Q = 6.78(φ0 = 0.01) では,S が小さいと,Σ を大きくしても空隙率は下がらない (図-7(d)). 実際, 質量流束を小さくすると, 堆積層が形成されるまでに時間が掛かり, 小粒径粒子が下方の空隙に沈降してしまい
上部の空隙が大きくなると考えられる. 堆積層の上層（0.25L < z < 0.35L）と下層（0.15L < z < 0.25L）における粒径 d が dmax /3.0 より小さな小粒径粒子の割合 Γの時間変化を図-8 に示す. 図-8 から Q = 67.8(φ0 = 0.1)では, 堆積層が形成されるまでに,Γ は, 上層と下層で同程度の値を取るのに対し,Q = 6.78(φ0 = 0.01) では, 時
刻が大きくなるにつれて, 下層での Γ が単調増加していることが分かる. このことは相対的に堆積層上層で空隙が増加することを意味する. したがって, 混合粒度分

5.

結果, 考察7

図–4 複数球によるブリッジの生成.(a) 角分離大のブリッジ,(b) 角分離小のブリッジ

Z

Z

φ

φ
(a) Q = 20.34,d ≤ 0.1m,(b) Q ≤ 203.4 ,d = 0.01m

図–5 単一粒度分布の沈降シミュレーションで得られるプロット点の φ-Z 面上での分布. （○）Ωr.m.s= 7.5 × 10−2 , （□）∗
7.5 × 10−3 , （△）7.5 × 10−4 . またアスタリスクは真空重力場中での沈降シミュレーションの最終状態.(a) 流束パラ
メータ Q = 20.34 , 粒径 d ≤ 0.1m,(b) 流束パラメータ Q ≤ 203.4 , 粒径 d = 0.01m の沈降シミュレーションの最終状
態. Q の値が大きい場合に, 最終状態は RLP に近づく. 図中の曲線は Zhang らによって得られた真空中の沈降シミュ
レーションによる結果6) .

布の沈降では,Q が小さい場合に, 単一粒度分布の場合には見られない分級の効果が空隙率 λ の値に影響を及ぼすことが分かる.

6.

考えられる.
複合粒径モデルでは, 分級の効果が現れるので空隙率の質量流束依存性は単一粒径モデルのものとは異なってくる. 本研究では Talbot 型の粒度分布に従う非回転粒子に対して, 小粒径の割合および粒度の広がりを表す
パラメータ（それぞれ S および Σ）および質量流束と空隙率の関係について考察した. 我々は質量流束が大きい場合, 小粒径の割合に依らず, 粒度の広がり Σ が大きくなるほど, 空隙率が小さくなるという結果を得た.
また空隙率の減少率は小粒径の割合が大きくなるほど（S が小さくなるほど）大きくなり,Σ ≈ 2.0 で空隙率は0.3 程度まで小さくなるという結果を得た. これらは粒子付着法を利用した空隙率の評価方法1),2),5) と同様の傾向を示している. 一方質量流束が小さい場合, 小粒径
粒子の割合が大きくても空隙率の下がりは見られなくなる. 我々は堆積層の上層と下層で, 小粒径粒子の割合の時間発展を比較し, 質量流束が小さい場合には, 下層
で小粒径粒子の割合が増加していることを確認した. こ
れまでも離散粒度分布での振動を伴う多体粒子計算で,分級が起こることが示されている15) . 本論文で, 我々は連続粒度分布における沈降シミュレーションでも分級が起こることを確認した. これまでの河床変動モデル
では空隙率の変化は粒度分布の変化のみに起因するとまとめ

本研究では粒度分布と間隙率の間の関係を求めるため, 多体粒子沈降シミュレーションによって得られる堆積層の空隙率評価を行った. 従来の DE 法を用いた多体粒子シミュレーションでは, 接触相互作用の評価に対する計算コストが大きくなるため,DE 法は連続多粒度分
布に対してはこれまでほとんど適用されていない. した
がって我々は多粒度分布における粒子間接触相互作用の評価を能率的に行うための計算スキームとして, 多段
格子エルミート法を開発した9) .
本研究ではこのような多段格子エルミート法を用い,
粒子沈降シミュレーションによって得られる堆積層の空隙率の質量流束および粒度分布依存性を考察した. 単
一粒径モデルでは, 粒子の回転や粒径の大きさに依らず, 質量流束が大きくなると, 空隙率が大きくなるという結果が得られた. これは, 同時に落下する粒子が多くなるほど, ブリッジのできる割合が大きくなり, 堆積層のランダム充填での状態が RLP に近づくためであると

8

λ

Z

(a) λ の Q 依存性 (Ωr.m.s= 7.5 × 10−4 )∗(b)Z の Q 依存性 (Ωr.m.s= 7.5 × 10−4 )∗λZ(c) λ の Q 依存性 ( Ωr.m.s= 7.5 × 10−2 )∗(d)Z の Q 依存性 (Ωr.m.s= 7.5 × 10−2 )∗図–6 単一粒度分布の沈降シミュレーションにおける空隙率 λ および平均配位数 Z の流束パラメータ Q 依存性. 図より Ωr.m.s∗の値に依らず質量流束が大きくなると空隙率が上がり平均配位数が下がる. 従って堆積層の状態は流束が大きくなるほ
ど RLP の状態に近づく.

考えられているが, 今回の結果から粒子濃度や運動状態の変化によって空隙率が変化することが考えられる.

nology 72 (1992) 149.
5) G. T. Nolan and P. E. Kavanagh, Computer simulation of random packings of spheres with log-normal
distributions.Powder Technology 76 (1993) 309.
6) Z. P. Zhang, L. F. Liu, Y. D. Yuan and A. B. Yu, A
simulation study of the eﬀects of dynamic variableson the packing of spheres. Powder Technology 116(2001) 23.
7) C. Song, P. Wang and H. A. Makse, A phase diagram
for jammed matter. Nature 453 (2008) 629.
8) L. E. Silbert, D. Ertas, G. S. Grest, T. C. Halsey
and D. Levine Geometry of frictionless and frictional
sphere packings., Phys. Rev. E 65 (2002) 031304.
9) Y. Sota and A. Watanabe,International Journal
for Numerical Methods in Engineering 投稿中 (2011).
10) M. P. Allen and D. J. Tildesley, Computer Simulationof Liquids ( Oxford Science Publications, 1989) .
11) H. Sigurgeirsson, A. Stuart and W. Wan, Algorithms
for particle-Field Simulations with Collisions. J. Comput. Phys. 172 (2001) 766.
12) J. Makino, Optimal Order and Time-step Criterion
for Aarseth-type N-Body Integrators. Astrophys J
369 (1991) 200.
13) J. Makino and S. J. Aarseth, On a Hermite Integrator
with Ahmad-Cohen Scheme for Gravitational ManyBody Problems. Publ. Astron. Soc. Japan 44 (1992)
141.
14) S. Konstantinidis and K. D. Kokkotas MYRIAD: A
new N-body code for simulations of Star Clusters. Astronomy and Astrophysics 2010; 522: A70.
15) 原田英治, 後藤仁志, 酒井哲郎：分級過程の三次元性に関する計算力学的アプローチ, 水工学論文集, 46 (2002)

本研究では沈降シミュレーションの第一段階として,砂粒子多体系の運動に対して以下のことを省略した.

1. 個々の粒子の運動が周囲流体におよぼす寄与
2. 流体によって伝達される粒子相互間に働く力
一般に高濃度の粒子系と流体の共存系に関して
は,Anderson らによる流体運動に対する局所平均の式およびそれに対する流体抵抗力の式が知られている19),27) .
これらの式で特に空隙率が小さい充填層に対しては, 式(3) は Ergun の式18) に置き換えられる. したがって、より厳密な空隙率の議論では粒子系の運動の流体への抗力を考える必要があり, 今後の課題と言える.
参考文献
1) M. Sulaiman, D. Tsutsumi, M. Fujita and K. Hayashi,Classification of Grain Size Distribution Curves of
Bed Material and the Porosity. Annuals of Disas.
Prev. Res. Inst., Kyoto Univ. 50B (2007) 615.
2) M. Sulaiman, 混合砂礫の空隙率と空隙率の変化を考慮
した河床変動モデルに関する研究 (京都大学学位論
文, 2008).
3) P. A. Cundall and O. D. L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29(1979) 47.
4) G. T. Nolan and P. E. Kavanagh, Computer simulation of random packing of hard spheres, Powder Tech-

9

λ

λ

0.45

0.45

0.40

0.40

0.35

0.35

0.30

0.30

0.25
0.0

0.5

1.0

1.5

2.0

2.5

Σ

0.25
0.0

Σ
0.5

a) Q = Q0

1.5

2.0

2.5

2.0

2.5

b) Q = 0.5Q0

λ

λ

0.45

0.45

0.40

0.40

0.35

0.35

0.30

0.30

0.25
0.0

1.0

0.5

1.0

1.5

2.0

2.5

Σ

0.25
0.0

0.5

(c) Q = 0.3Q0

1.0

1.5

Σ

(d) Q = 0.1Q0

図–7 Talbot 型粒度分布のパラメータ S および Σ と堆積層の空隙率 λ との間の関係. ただし Ωr.m.s= 0 かつ Q0 ≡ 67.8 で∗(a) Q = Q0 ,(b) 0.5Q0 , (c) 0.3Q0 ,(d) 0.1Q0 . それぞれの図中での曲線は S = 0.0(赤）,1.0(青）,2.0(緑）を表している. 図中の黒線は Sulaimann らの粒子付着法による結果1),2) 。

Γ

Γ

0.8

0.5

上層(0.25<z/L<0.35)
下層(0.15<z/L<0.25)

0.4

0.6

0.3

0.4
0.2

上層(0.25<z/L<0.35)

0.2

0.1

下層(0.15<z/L<0.25)
0.0
2

3

4

5

0.0
20

6

(a) Q = Q0

30

40

50

60

(b) Q = 0.1Q0

図–8 Talbot 型粒度分布沈降シミュレーションで得られる堆積層の上層（0.25L < z < 0.35L）および下層（0.15L < z < 0.25L）における小粒径粒子 (d < dmax /3.0) の割合. ただし S = 0.0 で Q0 ≡ 67.8 として (a) Q = Q0 ,(b) 0.1Q0 . それぞれの図中での曲線は上層が赤, 下層が青を表している.

16) 後藤仁志, 林 稔, 安藤 怜, 酒井哲郎： 固液混相流解
析のための DEM-MPS 法の構築, 水工学論文集 47(2003) .
17) 後藤仁志, 林 稔, 安藤 怜, 鷲見 崇, 酒井哲郎： 砂礫混合層を伴う混相流解析のための DEM-MPS 法マルチ
スケールリンクの開発, 海岸工学論文集 50 (2003) .
18) S. Ergun, Fluid flow through packed columns. Chemical Engineering Progress 48 (1952) 89.
19) T. B. Anderson and R. Jackson, A fluid mechanical description of fluidized beds.. Industrial and Engineering Chemistry Fundamentals 6 (1967) 527.
20) I. Kim, S. Elghobashi and W. A. Sirignano, On
the equation for spherical-particle motion: eﬀect ofReynolds and acceleration numbers. J. Fluid Mech.
367 (1998) 221.
21) R. Kurose and S. Komori, Drag and lift forces on arotating sphere in a linear shear flow. J. Fluid Mech.

384 (1999) 183.
22) C. Yih, FLUID MECHANICS (WEST RIVER
PRESS, 1979).
23) L. Schiller and Z. Naumann, Uber die grundlegenden berechungen bei der schwerkraftaufbereitung. Z. Ver Deutsch. Ing.77 (1933) 318.
24) 後藤仁志, 数値流砂水理学, 森北出版,2004 年
25) J.Duran 粉粒体の物理学（中西 秀, 奥村 剛訳）, 吉岡書店,2002 年.
26) W. W. Rubey, Setting velocities of gravels, sand andsilt particles. Journal of Science 25 (1933) 325.
27) 川口 寿裕, 田中 敏嗣, 辻 裕：有限要素法による流動層の数値シミュレーション, 日本機械学会論文集, 58(1992) 2119.
(2011 年 3 月 8 日 受付)

10