粉体せん断流の理論的および数値的研究
齊藤 国靖2006 年 2 月 1 日

要旨
散逸粒子である粉体で構成される粉体流には通常の Newton 流体には見られない数多くの特徴がある。散逸的な斥力相互作用をする多体粒子系である
粉体流は、非ニュートン流体に分類され、非一様な分布をすることがほとん
どであり、境界や重力の影響を強く受ける。また、静的な領域と動的な領域が共存し、自由冷却や斜面流などといった特徴的な現象が数多くあり、未だに未解決な問題が多く、また、理論的な解釈が乏しいのが現状である。そのため、本研究では、粉体流の特徴のエッセンスを取り出す目的で、比較的単純な系であるせん断流を扱い、数値的な結果と理論的な解釈を報告する。

DEM（離散要素法）による分子シミュレーションにより、せん断粉体流はエネルギー的に一定な準安定状態の間にダイナミックな相構造を形成し、最終的に相分離した状態で定常化することが解った。定常状態については運動論を基に導出された連続方程式でその構造をよく記述することができ、高密度・低温な領域に対しても運動論が適用可能であるという結論と、粉体の回転自由度を法線方向の跳ね返り係数に繰り込むことができるという結果を報告する。また、通常の流体が見せる一様な解は粉体流においては不安定であるという線型安定性解析の結果も報告する。
ダイナミックな相構造が見られる緩和過程については、現在解析の途中であるが、これも運動論による連続方程式を数値的に解き、DEM による粒子シミュレーションと共通の結果を得ることができた。また、この数値解析によって、相分離の過程で重要となる項を摘出しようと考えている。

1

目次
第 1 章 はじめに

1.1

本研究の背景 . . . . . . . . . . . . . . . . . . . . . . . . . . .

4
4

1.2

本研究の動機と目的. . . . . . . . . . . . . . . . . . . . . . .

5. . . . . . . . . . . . . . . . . . . . . . .

77
8

第 2 章 粉体モデル

2.1
2.2

粉体粒子の相互作用

離散要素法（DEM） . . . . . . . . . . . . . . . . . . . . . . .

2.2.1
2.2.2

シミュレーションで用いる粉体モデル

2.2.3
2.2.4

２次元 . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

. . . . . . . . .
座標変換 . . . . . . . . . . . . . . . . . . . . . . . . . .

89

時間発展の差分スキーム . . . . . . . . . . . . . . . . . 13

第 3 章 跳ね返り係数

3.1
3.2

15
. . . . . . . . . . . . . . . . . . . . . . . 15
DEM における跳ね返り係数 . . . . . . . . . . . . . . . . . . . 17

3.3

理論値とシミュレーション値の比較 . . . . . . . . . . . . . . . 19

跳ね返り係数の定義

第 4 章 DEM による粉体せん断流の粒子シミュレーション

21

4.1
4.2

粉体せん断流のシミュレーション . . . . . . . . . . . . . . . . 21

4.3

シミュレーション結果 . . . . . . . . . . . . . . . . . . . . . . 23

せん断流モデル . . . . . . . . . . . . . . . . . . . . . . . . . . 21

4.3.1
4.3.2

粉体せん断流の様子 . . . . . . . . . . . . . . . . . . . 23
系の時間発展 . . . . . . . . . . . . . . . . . . . . . . . 25

第 5 章 回転自由度の繰り込み

5.1

35
連続方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

5.2
5.3

極性流体 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

5.4

回転自由度を繰り込んだ連続方程式 . . . . . . . . . . . . . . . 41

回転自由度の繰り込み . . . . . . . . . . . . . . . . . . . . . . 39

5.4.1
5.4.2

連続方程式 . . . . . . . . . . . . . . . . . . . . . . . . 41
スケーリング . . . . . . . . . . . . . . . . . . . . . . . 422

第 6 章 定常状態の解析

6.1
6.2
6.3
6.4
6.5

44

定常解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

DEM データと理論値の比較 . . . . . . . . . . . . . . . . . . . 47
定常状態での圧力について . . . . . . . . . . . . . . . . . . . . 48
弾性極限および強いせん断応力がかかる場合について . . . . . 48
高密度粉体ガス系の振る舞い . . . . . . . . . . . . . . . . . . . 49

第 7 章 Linear Stability Analysis

60

7.1

粉体せん断流の線型安定性解析 . . . . . . . . . . . . . . . . . 60

7.2
7.3

線型化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
安定性解析

7.3.1
7.3.2
7.3.3
7.4

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

無限系 . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

kx = 0 の場合 . . . . . . . . . . . . . . . . . . . . . . . 66
kx 6= 0 の場合 . . . . . . . . . . . . . . . . . . . . . . . 68

有限系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

第 8 章 緩和過程

71

8.1
8.2

緩和過程の数値解析の目的と方法 . . . . . . . . . . . . . . . . 71

8.3
8.4
8.5

境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

連続方程式の差分化. . . . . . . . . . . . . . . . . . . . . . . 72

完全な連続方程式の数値解析結果 . . . . . . . . . . . . . . . . 74
連続方程式の簡単化. . . . . . . . . . . . . . . . . . . . . . . 74

第 9 章 まとめ

9.1
9.2
9.3

81

議論―粉体流におけるバグノールド則 . . . . . . . . . . . . . . 81
結論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
展望と課題

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

3

第 1 章 はじめに

1.1

本研究の背景粉体およびその粒子集団である粉体流の性質は実に多様であり、土石流や砂嵐などの災害の例を持ち出すまでもなく、我々の身の回りでも粉体や粉体流に関する現象は数え切れないほど存在する。それらの多くはただサイエンスの対象としてではなく、薬学や工学、農学の分野でも粉体および粉体流が扱われている。粉体および粉体流は素粒子や原子核などの近代物理学の対象とは違い、古代からその存在を認識されてきた。しかし古くから認識されていたにも関わらず、物理学の対象として扱われ始めたのは比較的最近のことであり、そのメカニズムや統計的な性質は現在の統計物理学の手法をもってしても完全に明らかにされるものではない [1]。粉体および粉体流のメカニ
ズムの解明が困難なのは、そもそも構成要素である粉体についてあまりよく解っていないからである。粉体はまず有限の大きさと質量を持ち、形状も球体ではなく色々である。また、粉体同士の相互作用、つまり粉体同士の衝突を考えてみても、衝突の際に破損して質量が変わったり、剛体球のように瞬間的に衝突が終わるのではなく、接触時間が有限であり、その間に粉体粒子に作用する力が変化するなど、難しい問題を多々含んでいる。
これらの粉体および粉体流に関する困難をできるだけ単純化しながら、かつその性質を余すことなく記述しようとする試みから端を発し、現在世界中で研究が活発に行われており、議論の絶えることのない分野になっている。
まず、粉体の形ひとつでも、粉体のとり得る形状は無限にあるため話が進まない。従って、本論文では粉体は球状であるとして話を進める。また、粉体が衝突する際に破損したり変形したりする状況は考えないことにする。そうすると、粉体というのは球状を保ったまま運動する粒子として扱われる。しかし、粉体を単純化するにあたり、粉体が衝突する際に非弾性衝突をしてエネルギーが散逸されるという性質は欠かすことができない。これにより粉体
粒子は剛体粒子と一線を画すことができるからであり、この性質が粉体流の極めて特徴的な現象の大きな要因になってくるからである。従って、本研究でも、粉体とは非弾性衝突によりエネルギーを失う球形粒子であり、粉体流とは外部からエネルギーの供給があり、内部でエネルギーを散逸して定常化する非平衡系であると考える。
このような粉体および粉体流に対して粒子数密度が比較的に薄い場合には、

4

粉体ガス、すなわち粉体の集団を気体として扱ったときの理論がある程度確立されている。代表的なものは、非弾性ボルツマン（Boltzmann）方程式を用いて、チャップマン・エンスコッグ（Chapmann-Enskog）法による摂動論の枠組みで粉流体の流体方程式を導出するというものである。 [2] これは、粉体粒子の衝突を瞬間的なものと仮定し、非弾性衝突の効果を考慮した衝突項（collision term）を持つボルツマン方程式を、空間的に一様な解の回りで空間微分を展開パラメータとして展開する方法である。具体的にはブレイ（Brey）
らの報告に詳しい説明がある [3]。チャップマン・エンスコッグ法の枠組みでは、e → 1 の極限と、分子カオスの仮定が成立する程の低密度な粉体気体を前提としてはいるものの、粉体気体の一様冷却の性質を示したり、輸送係数の定量的予言を与えるなど、低密度な粉体流の性質を捉えることができる。
チャップマン・エンスコッグ法の枠組みでは、e → 1 や、分子カオスの仮定が成立する低密度な粉体気体であることが必要であったが、例えば自由冷却現象では、粉体粒子の非弾性衝突によって高密度なクラスター領域が生じる [4]。このように、高密度領域が形成されて空間的に非一様な状態になると、分子カオスの仮定が崩れ、粉体粒子間の相関が生じてくるため、非弾性
ボルツマン方程式の成立条件が怪しくなってくる。この粉体粒子間の相関を、
川崎とオッペンハイム（Oppenheim） [5] による多体相関の枠組みで捕らえた報告がノイエ（Noije）とエルンスト（Ernst）によってなされ [6–8]、粉体気体の長距離相関の存在を定量的に説明することができた。
一方、希薄系に対する粉体ガスのボルツマン方程式を有限濃度の系に拡張
した非弾性 Enskog 方程式があり、長距離相関が重要でない現象に対してはかなり高濃度までその適用が試みられ、半定量的に正しい結果を与えている。
また、粉体流の従う流体方程式の導出に関しては、チャップマン・エンスコッ
グ法によるもの [12] の他に、グラッド（Grad）展開によって粉体流の流体方程式を導出したジェンキンス（Jenkins）とリッチマン（Richmann）の論
文 [9]（第３章に詳述）などがあり、粉体流を解析的に理解しようという動きも盛んである [10, 11]。また、粉体流の数ある特徴的な現象、例えば重力下での粉体流の流れである斜面流や、本研究で扱ったような無重力状態での粉体流のせん断などを、このような流体方程式を解析的および数値的に扱うことで現象を理解しようという報告が数多くある [13, 16, 68]。

1.2

本研究の動機と目的先述したように、粉体流には通常の流体には見られない特徴的な現象が数多く存在する [17]。本研究の目的はこの粉体流の基本的な性質を数値的に扱い、理論的な理解を得ることである。特に粉体流は境界条件と外場に強く依存し、そのために静的な領域と動的な領域が共存する系も実現する。しかし、粉体流の基本的な性質を理解するためには、粉体流が強く依存してしまうよ5

うな外場、特に重力の影響を排除した方が明確であり、あまり特別な境界条
件を設けないことが必要である。また、粉体流の実験との対応をみることで、数値的な結果や、理論的な解釈の裏付けにもなるため、なるべく実験可能なモデルを想定するのがよい。このような理由から、まず無重力状態での２次
元粉体せん断流を扱う。また、実験での検証も考慮に入れて、有限なバルク
内でのせん断流を考えることにする。
粉体せん断流については多くの報告がある。ほとんどの場合が無重力状態
であり、せん断面が周期境界条件（Lees-Edwards 境界）か、本研究のように壁を設けてあるかでまず二分される。Lees-Edwards 境界ではキャンベル（Campbell） [18] やタン（Tan） [19] などが、後に述べる離散要素法（discreteelement method : DEM）による分子シミュレーションを行っている。また、
せん断面を壁にした分子シミュレーションについては、ポップケン（Popken）とクリアリィ（Cleary） [20] などが報告している。いずれも一様なせん断流は不安定であり、通常の流体におけるせん断流の様子とは大きく異なる結果を報告している。

6

第 2 章 粉体モデル

2.1

粉体粒子の相互作用

粉体ガスの構成要素は多数の粉体粒子であり、その相互作用を最も単純に考えると弾性体同士の斥力相互作用である。弾性体同士の衝突を考えるには、弾性体の接触理論が基礎になる。弾性体については古くから研究があり、最も有名なのがヘルツ（Hertz）による固体接触理論 [21] であり、様々な弾性論の教科書に記載されている [22, 23]。ヘルツの固体接触理論によると、弾性体同士が接触するときに働く力は、３次元球の接触時における２体問題によって導かれる。
まず、法線方向に働く力について概説する。両物体の圧縮力の合力を Fn とする。接触時の力のため両物体は変形し、その変形した物体の接触面を半径b の円とする。この半径 b および接触領域の δn は変形前の粒子半径 a、a0 を
使うと、
b
δn

aa0 1/3)
a + a0
1
1= Fn2/3 [D2 ( + 0 )]a a
= Fn1/3 (D

(2.1)
(2.2)

で与えられる。ここで、
0

3 1 − σ02
1 − σ02
(
+
)
(2.3)
4
E
E
であり、σ0 、E はポアソン比、ヤング率と呼ばれる物質定数である。従って、3/2
Fn ∝ δn という関係が得られる。
D=

次に、接線方向に働く力についてはミンドリン（Mindlin）により論じられている [24]。δt を弾性体が接触した瞬間からの接線方向の変位とすると、先にでてきた接触面の半径 b を用いて、

Ft= 4bδt ((2 − σ0 )(1 + σ0 ) (2 − σ00 )(1 + σ00 ) −1
+
)
E
E
1/2

で与えられる。従って、Ft ∝ δt δn(2.4)という関係が得られる。

さらに、粉体粒子の衝突の際には物質間の粘性による散逸力が働くことも知られており、接線方向に対しては、粉体粒子同士が接線方向に高速で衝突したときに、動摩擦係数 µ に対して接線方向の力が Ft > µFn に達して滑り始めて動摩擦力（クーロン摩擦）に切り替わることも考慮に入れなければな
らない。

7

2.2

離散要素法（DEM）

2.2.1

シミュレーションで用いる粉体モデル

先述したとおり、粉体粒子の接触の際には、弾性論より導かれる弾性力と、物質間の粘性による散逸力および動摩擦力も考慮に入れた摩擦力が働く。これらの粉体粒子の相互作用を再現するため、本シミュレーションでは、法線方向および接線方向に対して、弾性力は Hertz 力の代わりに線型バネによる斥力を用いる。これは、２次元系における Hertz 則は線型バネに対数補正が3/2付いたものに置き換えられるため、Fn ∝ δnという関係よりも線型バネの方がより現実に近いモデルになることと、計算の単純化と効率化を図るためである。ちなみに、粉体に斥力が働くことを本質とするため、斥力モデルの選択に結果の詳細は依らないという報告もある [25]。
また、物質間の粘性による散逸力には相対速度に比例する粘性抵抗で置き換え、接線方向に対しては Ft > µ|Fn | に達したときにクーロン摩擦に切り替
わるようになっている。さらに、粉体粒子は接触の際に変形するため、慣性

図 2.1: 粉体粒子の相互作用のモデル
モーメントが時間変化するはずであるが、シミュレーションでは簡単のため
接触中の粒子の慣性モーメントは一定と考える。また、衝突の際に粒子が破損し、質量が変わるといった状況もここでは簡単のため考慮に入れていない。
このような２体衝突のモデルを多体衝突に適用し、多粒子系を扱えるようにしたのが離散要素法（DEM）と呼ばれる方法である。DEM は岩盤力学の研究者であったクンドール（Cundall）によって提唱され、広く粉体粒子のシ
ミュレーションに用いられている方法である。従って、DEM とは各粒子の並進自由度と回転自由度のみを残し、それ以外の内部振動に伴う自由度の配分や気体への熱伝導などといった全ての自由度を散逸として系のパラメータに押し込めたモデルである。

8

以上により、i 番目と j 番目の粒子が衝突する際、i 番目の粒子に対して法ij
線方向および接線方向それぞれに働く力を Fij
n 、Ft とすると、

Fij
n
Fij
t

ij
= −kn uijn − ηn u̇n(ij−kt uijt − ηt u̇t=−µ|Fijn |tij(2.5)(Ft < µ|Fn |)(Ft > µ|Fn |)(2.6)となる。ここで、上付きのドットは時間微分を表し、tij は接触粒子同士の中心を結ぶ線分に対しての接線方向の単位ベクトルである。また、uij
n および

uij
t は、i 番目と j 番目の粒子の中心を結ぶ線分に対しての法線方向および接線方向の相対変位であり、

uij
n
uij
t

= xi − xj − (xi (0) − xj (0))σ= xi − xj − (xi (0) − xj (0)) + (∆φi + ∆φj )tij
2

(2.7)
(2.8)

である。ここで、xi は i 番目の粒子の位置であり、xi (0) は i 番目の粒子の接触直前の位置である。また、σ は粉体粒子直径、∆φi は i 番目の粒子の角度変位である。
従って、衝突によって i 番目の粒子が従う運動方程式は(Pij
mẍi = j (Fijn + Ft )P
I ω̇i = σ2 j Mij(2.9)

である。ここで、ωi は i 番目の粒子の角速度（つまり ωi = ∆φ̇i ）であり、Mij
は i 番目の粒子に働くモーメントMij= nij × Fijt(2.10)である。また、nij は接触粒子同士の中心を結ぶ線分に対しての法線方向の単位ベクトルである。

2.2.2

座標変換

DEM では、粒子の位置、速度および角速度はすべて直交座標（xyz 座標）で表され、接触したときに接触粒子同士の中心を結ぶ線分に対して法線方向と接線方向の相対変位と相対速度から力を算出し、それを再び直交座標に戻すというプロセスが必要である。ここでは接触している２粒子の直交座標での相対変位と相対速度を、接触粒子同士の中心を結ぶ線分に対して法線方向と接線方向の相対変位と相対速度に変換する座標変換の方法を述べる。

i 番目の粒子と j 番目の粒子が接触したとき、i 番目の粒子の中心を原点として相対変位と相対速度を球座標 (r, Θ, Φ) の各単位ベクトル（er ,eΘ ,eΦ ）の方向に変換すれば、法線方向が er 、２つの接線方向 t 、s がそれぞれ eΘ 、eΦの方向に対応するようになる。従って、直交座標での相対変位と相対速度を、

9

まず y 軸回りに Θ だけ回し、次に z 軸回りに Φ だけ回せば、法線方向および接線方向の相対変位と相対速度に変換できる。つまり、i 軸回りに角度 Θ
だけ回転させる回転行列を

Ri (Θ)(2.11)と書くと、球座標での単位ベクトルは、行列R(Φ, Θ)≡ Rz (Φ)Ry (Θ)
cos Θ cos Φ=  cos Θ sin Φ− sin Θ− sin Φcos Φ
0

sin Θ cos Φsin Θ sin Φ (2.12)

cos Θ

によって、

er

= R(Φ, Θ)ez

(2.13)eΘ= R(Φ, Θ)ex

(2.14)eΦ= R(Φ, Θ)ey

(2.15)と変換される。任意のベクトル A は
 
Ax
 
A = Ay 

Az= Ax ex + Ay ey + Az ez= Ar er + AΘ eΘ + AΦ eΦ
= Ar R(Φ, Θ)ez + AΘ R(Φ, Θ)ex + AΦ R(Φ, Θ)ey
 
AΘ
 
= R(Φ, Θ) AΦ (2.16)

Arと変換され、相対変位と相対速度の変換も決められる。直交行列の性質 Rt (Φ, Θ)R(Φ, Θ) =1 を用いれば、直交座標の各成分から球座標の各成分への変換則は
 
 
AΘ
Ax
 
 
t
=
R
(Φ,
Θ)
AΦ 
Ay 
Ar

(2.17)Azで表される。
次に接触点での粒子の回転を考える。i 番目の粒子に対して、角速度ベクト10ル ωi を用いて、接触点での速度 vi のうち、回転による寄与は

vi = ωi ×

σ
er
2

i
i
= (ωri er + ωΘeΘ + ωΦeΦ ) ×==σ ii(ω eΘ − ωΘ
eΦ )
2 Φ

 i
ω
σ  Φi −ωΘ 
20σer
2(2.18)

である。但し、eΘ × er = −eΦ 、eΦ × er = eΘ を用いた。

(1) 式から

i
ωΦ
i
ωΘ

=
=

−ωxi sin Φ + ωyi cos Φ(2.19)

ωxi cos Θ cos Φ + ωyi cos Θ sin Φ − ωzi sin Θなので、接触点での速度は次のように変換される。
 i 
 i
ω
vΘ
σ  Φi  i
−ωΘ 
vΦ  =
2i0vr
− sin Φcos Φσ=− cos Θ cos Φ − cos Θ sin Φ
2
0
0
 i
ωx
σ t i=L  ωy 
2ωzi(2.20)

 i
0
ωx
 i
sin Θ ωy 
0ωzi(2.21)

但し、Lt

− sin Φ= − cos Θ cos Φ0cos Φ− cos Θ sin Φ
0

0




sin Θ
0(2.22)である。
同様に、モーメントの逆変換は

Mijσer × Fijt2
σ=er × (Frij er + FΘij eΘ + FΦij eΦ )2σ(−FΦij eΘ + FΘij eΦ )=2σ((−FΘij sin Φ − FΦij cos Θ cos Φ)ex + (FΘij cos Φ − FΦij cos Θ sin Φ)ey=
2+FΦij sin Θez )

=

11

なので




Mxij
 ij 
My  =Mzij



− sin Φσ
 cos Φ20
 ij Fσ  Θij L FΦ 2
Frij=

− cos Θ cos Φ− cos Θ sin Φ
sin Θ  ij 0FΘ  ij 0 FΦ 
0Frij(2.23)

である。
以上の操作を、接触する２粒子 i、j に対し、粒子間距離q
Rij =(xj − xi )2 + (yj − yi )2 + (zj − zi )2(2.24)

および

αij

=

q
(xj − xi )2 + (yj − yi )2(2.25)を用いて、
zj −ziRijx −x
cos Φij = jαij icos Θij =αsin Θij = Rijijsin Φij =yj −yiαij(2.26)と定義される角度を用いて計算すればよい。なお、αij = 0 の時には

zj − zi ≥ 0 ならばΘ = 0

zj − zi ≤ 0 ならばΘ = π(2.27)

とする。

2.2.3

２次元

２次元の場合は、３次元の結果を xy 平面に射影すればよい（尚、xz 平面に射影すると、θ の変域が 0 ≤ θ ≤ π であるため、x ≤ 0 の領域を記述できないことに注意。）。３次元系における z 方向の相対変位および x 軸、y 軸回りの角度変位は 0 なので、

zj − zi = ∆φx = ∆φy = 0(2.28)となる。また、cos Θij = 0x −x
cos Φij = jαij isin Θij = 1
y −ysin Φij = jαij i12(2.29)

なので、Θ = 0 となり、回転行列はR(Φ,− sin Φ0π
) = 02−1となり、接触面への速度変換は
− sin Φt
L = 0cos Φsin Φ 
0

cos Φ
0

cos Φ
0

0

0



(2.30)

0
1

(2.31)

0となる。従って接触面での相対変位は  ij  ij 
ut
0 − sin Φ cos Φ∆x− sin Φ  ij  σ 
 ij 

+
=
∆y
0
cos
Φ
sin
Φ
u


 s

 0
2
ij
ij
un
−1
0
0
∆z
0


−∆z − r∆Φx sin Φ + r∆Φy cos Φ=  −∆x sin Φ + ∆y cos Φ + r∆Φz 

cos Φ
0
0



0
∆Φij
x


1 ∆Φij
y 
0

∆Φij
z

∆x cos Φ + ∆y sin Φ





0



= −∆x sin Φ + ∆y cos Φ + r∆Φz ∆x cos Φ + ∆y sin Φ(2.32)となる。相対速度についても同様に変換されるので、接触面での力は
Ã ! Ã
!Ã !
Fnij
cos Φ sin Φ
Fxij=(2.33)
Fsij− sin Φ cos Φ
Fyij

Ftij = −Fzij = 0と変換され、モーメントも
!
Ã
Mxij

Myij

=

Ã
!
−Ftij sin Φ
Ftij cos Φ

Mzij =(2.34)

=

Ã !
0
0

σ ij
F
2 s

(2.35)(2.36)と逆変換されて、よく知られた２次元系の結果と一致している。

2.2.4

時間発展の差分スキーム

粉体粒子の運動方程式

mv̇i

ij
= (Fijn + Ft )(2.37)

I ω̇iij(2.38)= M
13

を、時間刻み δt の２次の精度で差分する。時刻 δt たったときの i 番目の粒子の速度および角速度は、Adams-Bashforth の二段法を用いて

vi (t + δt) =ωi (t + δt) =δt(3Fi (t) − Fi (t − δt))2m
δt
ωi (t) + (3Mij (t) − Mij (t − δt))
2I
vi (t) +(2.39)(2.40)ijとなる。ただし、Fi ≡ Fij
n + Ft である。次に、位置はxi (t + δt)
φi (t + δt)δt(vi (t + δt) − vi (t − δt))2
δt
= φi (t) + (ωi (t + δt) − ωi (t − δt)) × nij
2= xi (t) +として時間発展する。これにより、本シミュレーションは、時間発展について δt2 の精度を保障する。

14

第 3 章 跳ね返り係数3.1

跳ね返り係数の定義

粉体粒子同士の衝突は非弾性衝突であり、跳ね返り係数は法線方向および接線方向で定義される。粉体粒子の跳ね返り係数、特に接線方向の跳ね返り
係数については様々な議論があり [26, 28]、その全容は明らかにされていな
い。ここでは法線方向および接線方向の跳ね返り係数の一般的な定義を述べ、DEM による跳ね返り係数の理論的な値と測定値について説明することにする。まず、図 3.1 のように、粉体粒子１と粉体粒子２の２体衝突を考える。衝突前後の位置、速度、角速度を
0

0

0

0

0

0

(x1 , v1 , ω1 ) → (x1 , v1 , ω1 )(x2 , v2 , ω2 ) → (x2 , v2 , ω2 )とする。但し、上付きの 0 は衝突後であることを示す。ここで、衝突時における接触面の法線方向の単位ベクトルとしてk=x2 − x1|x2 − x1 |(3.1)

を定義すると、衝突前後の接触面での相対速度は

G=g+s×k0

00

G =g +s ×k

図 3.1: ２体衝突

15

(3.2)
(3.3)

となる。ここで g と s はそれぞれ並進運動と回転運動の相対速度で、= v1 − v2σs =(ω1 + ω2 )
2

g

(3.4)
(3.5)

である (σ : 粒子直径)。 接線方向の単位ベクトルをj=(G × k) × k|G × k|(3.6)と定義すると、(s × k) · k = 0 などに注意して、法線方向と接線方向の跳ね
返り係数が以下のように定義できる。
0

G · k = −e(G · k)

(3.7)

0

G · j = −β(G · j)

(3.8)

法線方向の跳ね返り係数 e は物質により決る定数で、接線方向の跳ね返り係
数 β は衝突前の相対速度の接触面に対する入射角 γ によって変わり、現象論
的にはウォルトン（Walton）の式(
2−1 + µ(1 + e)(1 + mσ
4I ) cot γβ=β0(γ ≥ γ0 )(γ ≤ γ0 )(3.9)

で与えられる [31]。但し、m は粒子質量、I は粒子の慣性モーメント、µ は
クーロン摩擦係数である。ここで、γ は G と k のなす角で定義され、

G = G⊥ k + Gk j

(3.10)とした時のcot γ =

G⊥
Gk

(3.11)

として定義され、衝突前の相対速度の接触面に対する入射角を表す。つまり、衝突前の相対速度の接触面に対する入射角が臨界角 γ0 以下であれば接線方向の跳ね返り係数は β0 で一定であるが、γ0 を超えると接触面で滑りが起こり、接線方向の跳ね返り係数は cot γ に比例して大きくなる。つまり、臨界角 γ0は 接線方向に作用する力がクーロン摩擦

Fk = −µ|F⊥ |j(3.12)に切り替わる角度を表している。

(3.9) において、β0 および γ0 の値はウォルトンの現象論 [31] からは決めることができない。一方、理論的な β の取り扱いについてはモウ（Maw）等の理論 [34, 35] が知られており、國仲と早川 [27] もそれに従っている。しかし、

16

モウ等の理論は複雑なだけでなく、β を陽にきめることができないため、より簡略で物理的な理論が必要となっている。
また、(3.7) 式で定義した法線方向の跳ね返り係数 e は一般に衝突速度に強く依存する [36–41]。しかし、本研究で採用した線型バネのモデルでは法線方向の跳ね返り係数 e が定数となり、簡単化した取り扱いが可能になっている。

DEM における跳ね返り係数

3.2

DEM における粉体粒子が非弾性衝突を起こす場合、法線方向の跳ね返り
係数は、衝突する２つの粉体粒子間に働く力が線型バネによる弾性力と線型抵抗による粘性力であることから、減衰振動の方程式

M ün + ηn u̇n + kn un=0(3.13)

を解くことで求めることができる。但し、un は法線方向の相対変位であり、M は換算質量である。
また、接線方向の跳ね返り係数については、γ ≥ γ0 であれば (3.9) 式を使って β を計算でき、γ ≤ γ0 であれば衝突するときに滑りが無いと仮定し、法線方向と同様、減衰振動の方程式

M üt + ηt u̇t + kt ut= 0(3.14)

を解くことで β0 を求めることができる。但し、ut は接線方向の相対変位である。一般に減衰振動

M ü + k u̇ + ηu = 0の解は条件により３つに分けられる。

1. η

2

η
− 4M k < 0 の場合 a = 2M
、b =√(3.15)

4M k−η 2として2M

u(t) =

u(0) −atesin(bt)bu̇(t) =

u(0)e−at (cos(bt) −(3.16)asin(bt))b(3.17)

2. η 2 − 4M k = 0 の場合
u(t)= u(0)te−at

u̇(t)

=

(3.18)(1 − at)u(0)e−at(3.19)

3. η 2√− 4M k > 0 の場合
b0 =η 2 −4M kとして
2M

u(t) =

u(0) −atesinh(b0 t)bu̇(t) = u(0)e−at (cosh(b0 t) −17

(3.20)asinh(b0 t))
b0(3.21)

120%QWNQOD៺ᡂ

80

40
ធ⸅ജ
0
ᤨ㑆
図 3.2: DEM による分子シミュレーションで、２つの粉体粒子が接触している間に働く接線方向の力が時間的に変化する様子をプロットしたもの。入射角 cot γ が十分大きい場合であり、クーロン摩擦の判定条件により、接触中に滑ることはなく、摩擦力ではなく弾性力と粘性しか働かないことが解る。
従って、このときの接線方向の跳ね返り係数は β = β0 で定数であり、減衰振動の方程式を解くことで値を求めることができる。
２粒子が離れる時刻は、法線方向の相対変位 un (t) が 0 になる時で、この時
刻 tp での速度がちょうど衝突後の速度ということになる。従って、法線方向
および接線方向の跳ね返り係数は

e =
β

=

u̇n (tp )u̇n (0)u̇t (tp )−u̇t (0)−(3.22)(3.23)と表され、条件によって異なる形をとる。
例えば、２次元球の場合に、各パラメータの値を
kn =3.0 × 103 η =3.0n
kt =kn /4 η =η /2
t(3.24)nとすると、ηn2 − 4M kn < 0 となるので、1. より、tp は

tp =

π
bn

18

(3.25)

となり、従って、

e =
=

−

u̇n (tp )u̇n (0)πηnexp(− p)2mkn − ηn2(3.26)となる。一方、ηt2 − 4M kt = (ηn2 − 4M kn )/4 < 0 なのでu̇t (tp )u̇t (0)a0b0b0a0
exp(− t π)( 0t sin( t π) − cos( t π))
bn
bt
bn
bn

= −

β0

=

p
πAηt 2T /A − 1
πAηt1√= exp(− √
)[ p
sin(
)
ηn 2N − 1
ηn 2N − 1
2T /A − 1pπAηt 2T /A − 1√− cos()](3.27)
ηn 2N − 1となる。但し、

A=1+

mσ 2
4I

N=

πkn
4ηn2

T =

πkt
4ηt2

(3.28)である。実際に (3.24) 値を代入すると、

β0 ' 0.769235(3.29)である。

3.3

理論値とシミュレーション値の比較
I=

1
mσ 2
8

µ = 0.2(3.30)

として、(3.9) の理論値の傾きを求めるとµ(1 + e)(1 +

mσ 2
)
4I

=

7
µ(1 + e) ' 1.12442
2(3.31)である。以上の解析により求まった値とシミュレーション値を比較したのが
図 3.3 であり、理論とシミュレーションはよく一致している。また、ポリス
チレン球の衝突実験による測定結果 [32] と同じ特徴を捉えている。
このように、接線方向の跳ね返り係数は臨界角 γ0 を境に、衝突したときに粘着して滑らない状況から完全に滑っている状態に劇的に変化する。ここでの結果はウォルトンの式 (3.9) を比較的簡単な議論で初めて導いたものである。また、同様の議論が３次元系の衝突問題にも適用でき、モウ等の理論で用いられる複雑な議論を避けて、簡潔で物理的な理論から β を決定すること
ができる。

19

1
0.5
Ǫ 0
-0.5
-1
0

5

10
EQVǫ

15

20

図 3.3: ２次元球の場合の β の理論値（線）とシミュレーション値（データ点）の比較。kn = 3.0 × 103 , ηn = 3.0 , kt = kn /4 , ηt = ηn /2 。

20

第 4 章 DEM による粉体せん断流の粒子シミュレーション
4.1

粉体せん断流のシミュレーション

粉体流のシミュレーションには、DEM やイベントドリブンの剛体球などによる粒子シミュレーションや、流体方程式を差分法や有限要素法により数値的に解く方法などがある。分子シミュレーションによる方法は直接的である
ため、様々な状況を再現することができる。従って、連続体をベースにした数
値解析よりも実験的な側面が強い。一方、流体方程式を数値的に解く場合は境界条件や初期値の設定などに充分注意を払う必要があり、シミュレーショ
ン以前に系の条件を把握しておく必要がある。そのため、本研究では DEM
による粒子シミュレーションにより、粉体せん断流を再現するところから始める。
粉体せん断流の粒子シミュレーションに関する研究は盛んに行われており、キャンベル（Campbell） [18]、デイブ（Dave） [44]、タン（Tan） [19] など、数多い報告がされている。彼らはいずれもせん断をかける面を周期境界条件にした Lees-Edwards の境界条件を使っている。Lees-Edwards の境界条
件は、1972 年にリー（Lees）とエドワーズ（Edwards）により考案された方法で [45]、せん断流の分子シミュレーションの境界条件として広く用いられている。Lees-Edwards の境界条件は、いずれの境界も周期境界にとるため、幅が無限に広い系でのせん断流の一部をシミュレーションしていることになる。
一方、有限な幅におけるせん断流の分子シミュレーションもあり、ポップ
ケン（Popken）とクリアリィ（Cleary） [20] は摩擦が働く壁を動かすことで粉体ガスにせん断をかける分子シミュレーションを行い、有限要素法解析による結果との一致を報告している。
本研究では、実験との対応を容易にするという目的で、有限な幅での粉体せん断流を扱う。壁には粉体ガスを構成している粉体と同じ粉体を埋め込み、
平坦な壁よりもせん断がかかるようなモデル（バンピー境界）を適用した。
以下に具体的な設定と分子シミュレーションの結果を述べる。

4.2

せん断流モデル

無重力下で、２次元粉体ガスにせん断をかける。

21

図 4.1: せん断流モデル
粉体粒子の法線方向の跳ね返り係数は e = 0.85、接線方向の跳ね返り係数は初期状態
粉体粒子を容器内にランダムに配置し、Maxwell 分布に従ってそれぞれ初速度を与えた一様な粉体ガスを初期条件にした。

境界条件
せん断方向（x 軸方向）の境界は周期境界条件とし、粒子を埋め込んだ壁（∓y の位置）をそれぞれ x 方向に ±U/2 の速度で動かすことでせん断をかけた。

粒子を埋め込んだ壁は、静止した状態から、以下の運動方程式に従って時間発展する。

Mw

dudt= −γw (u − U/2) + F ex

(4.1)

但し、u は壁の速度、Mw は壁の質量であり、F ex は粉体粒子が壁に埋め込んだ粒子と衝突することにより壁に与えられる外力であり、緩和係数 γw (= 10)で壁の速度は終端速度 ±U/2 に到達する。実際には壁の質量密度が粉体粒子の 5 × 106 倍であるため、粉体粒子の衝突により与えられる外力の影響はほとんどなく、壁は連続的に終端速度に達する。

22

表 4.1: 系の基本的な値
粒子数

N = 5000

バルク幅

W = 180σ∆ = 180σν̄ = 0.12

容器高さ
平均面積占有率

表 4.2: DEM で用いる各パラメータの値

kn = 3.0 × 103[(U/2σ)2 ]ηn = 3.0
kt = kn /4[U/2σ][(U/2σ)2 ]ηt = ηn /2µ = 0.2[U/2σ]

スケーリング
長さを粉体粒子直径 σ で、時間スケールを 2σ/U でスケールして無次元化を行った。また、粉体粒子質量を m = 1 としてある。また、DEM で用いる
各パラメータの値は第２章の例で使用したものと同一であり、表 4.2 にまと
めておく。

4.3

シミュレーション結果

4.3.1

粉体せん断流の様子

図 4.2 はシミュレーションにおける粉体粒子を描いたもので、点は粉体粒子を表し、枠が容器に相当する。系は左上の番号順に時間発展し、時間間隔は一定である。
図 4.2 から読み取れる系の時間発展は以下の様にまとめられる。まず、初期の一様な状態において壁際の粒子が壁に埋め込まれた粒子と衝突し (0)、バ
ルク中央に向かって弾かれ、壁からやや離れた位置にある粒子と非弾性衝突を起こしてクラスターを形成する (1)。このクラスターは衝突の余力でそのままバルク中央へと移動していき (2,3)、やがてバルク中央で衝突、合体し(4)、１つの大きなクラスターを形成する (5)。

23













[

Z
図 4.2: シミュレーションの様子。y = ∓∆ の境界で、粒子を埋め込んだ壁がそれぞれ x 方向に速度 ±U/2 で動くことでせん断がかかっている。図中の番号は時間発展の順番を表す。

24

'PGTI[40002000   

0
0

100



200
VKOG

300

400

図 4.3: 系の全エネルギーの時間発展。早い段階にエネルギー的にほぼ一定な状態に落ち着き、その後ゆっくりと緩和しながら定常状態に移行していく。

4.3.2

系の時間発展エネルギー
系の全エネルギーはE= (運動エネルギー) + (回転エネルギー)
=

N
X
I
m
( vi2 + ωi2 )22n=1で定義され、全粉体粒子の並進運動のエネルギーと回転運動のエネルギーの和である。
図 4.3 は系の全エネルギーの時間発展であり、図中の番号は時間発展の順番を表し、図 4.2 のものに一致する。図 4.3 より、系は比較的早い段階にエネ
ルギー的にほぼ一定な状態に落ち着き、その後ゆっくりと緩和しながら定常状態に移行していくのが解る。今後、全エネルギーがほぼ一定となる時間領域（図 4.3 での 1,2,3,4）を準安定状態と呼び、そこでの流体変数の時間発展を緩和過程と呼ぶことにする。

25

[
ǻ
K ⇟⋡࠮࡞

ǻK

㕙Ⓧ #


 ⇟⋡࠮࡞
Z
図 4.4: セルの取り方。

面積占有率
面積占有率は、粉体粒子の存在領域の総面積 A に対する粉体粒子の面積が占める割合で定義され、

ν=

1 X πσ 2
A i 4

(4.2)

となる。ここで、和は面積 A の領域に存在する全ての粉体粒子についてとる。シミュレーションで測定する場合には、せん断方向への分布を平均してしまい、壁に垂直な方向への分布をとることにする。これは、図 4.4 の様にバルクの面積をせん断方向に Nb 分割したセルに対する粉体粒子の面積占有率を求めることで、壁に垂直な方向の面積占有率の分布をとるということである。つまり、せん断方向のバルク幅を ∆x 、壁に垂直な方向のバルク幅を∆y としたとき、∆y を Nb 個に等分割し、i 番目のセル（面積 A = δy ∆x 、但し δy = ∆y /Nb ）に対する粉体粒子の面積占有率を y = δy i での面積占有率として分布をとるということである。シミュレーションではバルクを Nb = 180
分割し、δy = σ にとっている。
図 4.5 より、系がエネルギー的に一定な準安定な状態の内に、壁からやや離れた位置にクラスターが形成され、その後クラスターがバルク中央へ移動していき、バルク中央で衝突・合体をする。最終的にバルク中央でクラスターを形成したまま定常状態に至る。このように、せん断によって粉体流中に相分離が生じる。

26

粉体温度
粉体温度は、１体分布関数 f (x, c, t) を用いて
Z
1
1 ∞
θ =
dx (c − u)2 f (x, c, t)n −∞ 2(4.3)

で定義される。但し、n は粒子数密度である。シミュレーションでは、面積占有率と同じ各セル内の粒子に対して

θ

=

1X1
(c − u)2
n i 2

(4.4)

を計算し、粉体温度の分布とした。
図 4.5 と図 4.6 を比較すると、粉体温度と面積占有率は、時間発展する間ほぼ一対一に対応し、クラスター領域では温度が極端に低くなる。つまり、面積占有率の極大値（クラスターのピーク）の位置は、粉体温度の極小値の位置に対応しながら時間発展していく。つまり、初期の一様な分布から出発し、壁よりやや離れた位置にクラスターが形成されると、その位置での粉体温度はやや低くなり、クラスターの移動に伴い、低温な領域も移動していく。また、クラスターがバルク中央で衝突・合体して定常状態に至ると、バルク中
央での粉体温度はほとんど 0 になってしまう。

速度場（せん断方向）
せん断方向の速度場は、面積占有率と同じセル内にある粒子の各速度のせ
ん断方向の成分を平均したものである。
図 4.7 より、平均速度が 0 の初期状態から、粒子を埋め込んだ壁に引きずられる形でせん断方向の速度場が時間発展する。次に、各壁からやや離れた位置にクラスターが形成されるが、クラスターに挟まれた領域（プラグ領域と呼ぶ）でせん断方向の速度場は殆ど 0 のままである。これは、壁に弾かれた粉体粒子がせん断方向の速度成分をもってバルク内に向かうも、クラスターに阻まれ、プラグ領域までせん断の影響が及ばないことによる。徐々にクラスターがバルク中央に向かって移動していき、プラグ領域が狭くなっていくに従い、平均速度 0 の領域も狭くなっていく。最終的にクラスターが衝突・合体すると、クラスター領域の速度場の傾きが希薄領域よりも低い状態になり、そのまま定常状態に至る。従って、相分離した状態が保たれると、クラスター領域の速度場の傾きが希薄領域よりも低い状態で一定となる。

速度場（壁に垂直方向）
壁に垂直な方向の速度場も、面積占有率と同じセル内にある粒子の各速度の壁に垂直な方向の成分を平均したものである。

27

図 4.8 より、平均速度が 0 の初期状態から始まり、壁際の粉体粒子が壁に弾かれるため、壁に垂直方向の速度場が ∓∆/2 のそれぞれの壁際で極大および極小となる。次に、壁よりもやや離れた位置にクラスターが形成されると、壁に垂直方向の速度場はクラスターの壁側の界面で極大となる。これは、な
おも壁に弾かれた粉体粒子がクラスターに向かって飛んできており、クラス
ターの壁側の界面でクラスターと衝突して速度を落とすためである。この衝突の余力で、クラスター全体はさらにバルク中央へと押されていく。最終的にクラスターがバルク中央で衝突・合体し、定常状態に至ると、クラスターを形成する粉体粒子はほぼ停留し、さらに壁に弾かれる粉体粒子がほとんどなくなるため、壁に垂直方向の速度場はほぼ 0 になる。

角速度場
図 4.9 より、角速度場の分布も、面積占有率と同じセル内にある粒子の各角速度を平均したものである。シミュレーションは２次元であるため、角速度は z 成分 ωz のみである。
初期のランダムな状態では平均 0 の角速度分布であるが、せん断の方向か
ら明らかなように、壁付近の粒子は壁に弾かれて負の角速度を持つようになる。また、バルク内の粒子は互いに接触することで次第に角速度を失っていき、バルク中央付近では角速度がほぼ 0 になる。壁に弾かれた粉体粒子の持つ角速度はバルク中央まで伝播することはなく、壁に弾かれることによる角速度の供給と、バルク内での粉体粒子同士の接触による角速度の消失とが均衡を保ったまま定常状態に至る。

回転温度
回転温度は

Tr

=

1
n

Z ∞

dx

−∞

I
(ωz − ω̄)2 f (x, v, t)
2m

(4.5)

で定義される。但し、シミュレーションは２次元であるため、角速度は z 成分 ωz のみであり、ω̄ は角速度場である。また、面積占有率と同じ各セル内の粒子に対してTr=

1X I(ωz − ω̄)2
n i 2m(4.6)

を計算し、回転温度の分布とした。
図 4.10 より、初期のランダムな状態では一様に分布しているが、早い段階から全体的に小さな値をとるようになり、クラスターの界面付近で比較的に値が大きくなるものの、定常状態に至るまで特に変化がない。

28

1

1





0.8

0.8

0.6
ǵ
0.4

0.6
ǵ
0.4

0.2

0.2

0
-90 -60 -30

0
[

30

60

90

0
-90 -60 -30

0
[

30

60

1

1




0.8

0.8

0.6
ǵ
0.4

0.6
ǵ
0.4

0.2

0.2

0
-90 -60 -30

0
[

30

60

90

1

0
-90 -60 -30

0
[

30

60

90

1



0.8

0.8

0.6
ǵ
0.4

0.6
ǵ
0.4

0.2

0.2

0
-90 -60 -30

90

0
[

30

60

90

0
-90 -60 -30

図 4.5: 面積占有率の時間発展。

29

0
[

30

60

90

12

12



10
8

8

ǰ6

ǰ6

4

4

2

2

0
-90 -60 -30



10

0
[

30

60

90

0
-90 -60 -30

0
[

30

60

12

12


10

8

ǰ6

ǰ6

4

4

2
0
-90 -60 -30

2
30

60



10

8

0
[

90

12

0
-90 -60 -30

0
[

30

60

12


10

8

ǰ6

ǰ6

4

4

2
0
-90 -60 -30

2
0-90 -60 -30

30

60

90



10

8

0
[

90

90

図 4.6: 粉体温度の時間発展。

30

0
[

30

60

90

0.8

0.8




0.4

0.4

ux 0

ux 0

-0.4

-0.4

-0.8
-90 -60 -30

0
y

30

60

90

0.8

-0.8
-90 -60 -30

0
y

30

60

0.8




0.4

0.4

ux 0

ux 0

-0.4

-0.4

-0.8
-90 -60 -30

0
y

30

60

0.8

90

-0.8
-90 -60 -30

0.4

ux 0

ux 0

-0.4

-0.4
0
y

30

60

0
y

30

60

0.8



0.4

-0.8
-90 -60 -30

90

90

-0.8
-90 -60 -30

図 4.7: せん断方向の速度場の時間発展。

31

90


0
y

30

60

90

0.08

0.08



0.04
W[

0.04
W[

0

0
-0.04

-0.04
-0.08
-90 -60 -30

0
[

30

60

-0.08
-90 -60 -30

90

0.08

0
[

30

60

0.08




0.04
W[

0.04
W[

0
-0.04

0
-0.04

-0.08
-90 -60 -30

0
[

30

60

-0.08
-90 -60 -30

90

0.08

0
[

30

60

90

0.08



0.04
W[

90

0.04
W[

0
-0.04

0
-0.04

-0.08
-90 -60 -30

0
[

30

60

90

-0.08
-90 -60 -30

図 4.8: 壁に垂直方向の速度場の時間発展。

32

0
[

30

60

90

0.4

0.4




0
ȁ\

0

-0.4

ȁ\

-0.4

-0.8

-0.8

-1.2

-1.2

-1.6
-90 -60 -30

0
[

30

60

-1.6
-90 -60 -30

90

0.4

0
[

30

60

0.4




0
ȁ\

0

-0.4

ȁ\

-0.4

-0.8

-0.8

-1.2

-1.2

-1.6
-90 -60 -30

0
[

30

60

-1.6
-90 -60 -30

90

0.4

0
[

30

60



0

0

-0.4

ȁ\

-0.4

-0.8

-0.8

-1.2

-1.2

-1.6
-90 -60 -30

90

0.4


ȁ\

90

0
[

30

60

90

-1.6
-90 -60 -30

図 4.9: 角速度場の時間発展。

33

0
[

30

60

90

1.2

1.2



0.8

0.8

6T

6T
0.4

0
-90 -60 -30

0.4

0
[

30

60

0
-90 -60 -30

90

1.2

0
[

30

60

90

1.2




0.8

0.8

6T

6T
0.4

0
-90 -60 -30

0.4

0
[

30

60

0
-90 -60 -30

90

1.2

0
[

30

60

90

1.2





0.8

0.8

6T

6T
0.4

0
-90 -60 -30

0.4

0
[

30

60

90

0
-90 -60 -30

図 4.10: 回転温度の時間発展。

34

0
[

30

60

90

第 5 章 回転自由度の繰り込み

5.1

連続方程式

粉体粒子を構成要素とする多体系を連続体として記述する方法としてはい
くつか報告がある [9–12] が、何れも運動論から出発し、チャップマン・エン
スコッグ法や、後に述べるグラッド（Grad）展開などの方法により分布関数から速度の各モーメントを計算することで連続方程式を導出するものである。
本研究の解析では、２次元系での解析がなされているジェンキンス・リッチ
マン（Jenkins and Richman [9]）による連続方程式を採用する。
まず、連続方程式の導出はボルツマン（Boltzmann）方程式∂f+ c · ∇f
∂t

= Ψc

(5.1)

のモーメントを計算することから始まる。ここで、f (x, c, ω, t) は１体分布関数
であり、x、c、ω はそれぞれ着目する粒子の位置、速度、角速度であり、回転自由度も含めた形であることに注意していただきたい。Ψc は衝突項（collisionterm）であり、２体分布関数 f (2) (x1 , c1 , ω1 ; x2 , c2 , ω2 ; t) を用いて、
ZZ
2Ψc = σdc1 dω1dk(c12 · k)
c12 ·k>0

×f (2) (x − σk, c1 , ω1 ; x, c2 , ω2 ; t)(5.2)

である。ここで、σ は粒子直径であり、k は着目している粒子 1 の中心から
衝突相手の粒子 2 の中心に向かう方向の単位ベクトルである。積分は粒子 1と粒子 2 の衝突が起こるために、これらの粒子の相対速度 c12 = c1 − c2 を用いて、c12 · k > 0 が満たされる範囲内で実行される。

次に、速度の i 次のモーメント ψ i をボルツマン方程式の両辺にかけて速度および角速度で積分することで、モーメントの平均の時間発展を決定する方R
程式を導く。１体分布関数による平均を <>≡ dcdωf で表すと、ボルツマ
ン方程式は∂< ψ i > +c · ∇ < ψ i > =
∂t

Ψc

(5.3)

となる。ここで、衝突項は、衝突による位相体積要素の収縮を無視すると
ZZ
0dk(c12 · k)(ψ2i − ψ2i )Ψc = σ 2 dc1 dc2 dω1 dω2
c12 ·k>0

×f (2) (x − σk, c1 , ω1 ; x, c2 , ω2 ; t)
35

(5.4)

である。モーメント ψ i に付けた上付きのダッシュ0 は衝突後のものであることを示し、下付きの番号は粒子番号である。ここで 1 と 2 の添字を入れ換え、

k を −k に書き換えることで、ZZ
Φc = σ 2 dc1 dc2 dω1 dω20c12 ·k>0

dk(c12 · k)(ψ1i − ψ1i )×f (2) (x, c1 , ω1 ; x + σk, c2 , ω2 ; t)(5.5)

としておく。
次に、２体分布関数

f (2) (x − σk, c1 , ω1 ; x, c2 , ω2 ; t)の引数 x−σk と x をひとくくりとして考え、テイラー（Taylor）展開すると、

f (2) (x − σk, c1 , ω1 ; x, c2 , ω2 ; t) = f (2) (x, c1 , ω1 ; x + σk, c2 , ω2 ; t)11
+ (−σk · ∇ + (σk · ∇)2 − (σk · ∇)3 + · · · )2!
3!
×f (2) (x, c1 , ω1 ; x + σk, c2 , ω2 ; t)
(5.6)
移項して、

f (2) (x, c1 , ω1 ; x + σk, c2 , ω2 ; t) = f (2) (x − σk, c1 , ω1 ; x, c2 , ω2 ; t)1
1
+ (σk · ∇ − (σk · ∇)2 + (σk · ∇)3 − · · · )2!
3!
×f (2) (x, c1 , ω1 ; x + σk, c2 , ω2 ; t)
(5.7)(5.7) 式を (5.5) 式に代入し、(5.4) 式と足すと
Z
Z
0
0
2Ψc =
dc1 dc2 dω1 dω2
dk(c12 · k)(ψ1i + ψ2i − ψ1i − ψ2i )f (2) (x − σk, c1 , ω1 ; x, c2 , ω2 ; t)c12 ·k>0ZZ01
1+d∇ · dc1 dc2 dω1 dω2
dk(c12 · k)k(ψ1i − ψ1i )(1 − σk · ∇ + (σk · ∇)2 − · · · )2!
3!
c12 ·k>0
×f (2) (x, c1 , ω1 ; x + σk, c2 , ω2 ; t)

(5.8)

よって

Ψc

= −∇ · χk + χc

(5.9)

但し、χc は collision part と呼ばれ、速度のモーメント ψ i に作用する演算子の形でiχc [ψ ] =

1
2

Z

×f

dc1 dc2 dω1 dω2

(2)

Z

0

c12 ·k>0(x − σk, c1 , ω1 ; x, c2 , ω2 ; t)360dk(c12 · k)(ψ1i + ψ2i − ψ1i − ψ2i )(5.10)

と書ける。また、χk は kinetic part と呼ばれ、速度のモーメント ψ i に作用する演算子の形で
Z
Z
0
σ
1
1
i
χk [ψ ] = −dc1 dc2 dω1 dω2
dk(c12 · k)k(ψ1i − ψ1i )(1 − σk · ∇ + (σk · ∇)2 − · · · )22!
3!
c12 ·k>0

×f (2) (x, c1 , ω1 ; x + σk, c2 , ω2 ; t)(5.11)と書ける。ここで、粒子直径が充分小さいことと、連続体では長波長極限をとることから、kinetic part は最初の項のみ残し、
ZZ0
σdk(c12 · k)k(ψ1i − ψ1i )χk [ψ i ] = −dc1 dc2 dω1 dω2
2
c12 ·k>0

×f (2) (x, c1 , ω1 ; x + σk, c2 , ω2 ; t)(5.12)となる。次に、２体分布関数に対しては、分子カオスの仮定（molecular chaos）
を用いて

f (2) (x − σk, c1 , ω1 ; x, c2 , ω2 ; t)

' g(x −σk)f (x − σk, c1 , ω1 , t)f (x, c2 , ω2 , t)2とする。但し、g(x) は２体相関関数（Radial distribution function）であり、ここでは位置ではなく、面積占有率の関数として

g(ν) =

16 − 7ν
16(ν − 1)2(5.13)としたものを使用する [42]。なお、面積占有率 ν とは、２次元の場合、粒子数密度 n を用いて ν = nπd2 /4 と定義され、系の面積に対して粒子が占める面積の割合を表すものである。以上が散逸のないときの Enskog 理論の考え方である。実際には衝突によって位相体積要素が収縮するので (5.4) 式が成立しない。このモデルは法線方向の跳ね返り係数 e が比較的 1 に近い場合のみ有効で、非弾性衝突の効果を衝突ルールにのみ押し込めた近似的なものであることに注意する必要がある。
次に、１体分布関数をマクスウェル分布

f (0) (c, ω, t)

=

n
Ci Ci
Ω2
exp(−
−
)
1/2
2T
2Θ
(2πT )(2πΘ)(5.14)のまわりで展開することを考える。ここで、C(= c − u) は速度のゆらぎ（uは平均速度）であり、Ω(= ω − ω̄) は角速度のゆらぎである（ω̄ は平均角速度）。また、T は粉体温度であり、T ≡ 21 < Ci Ci > と定義され、Θ は粉体の回転温度であり Θ ≡ 12 < Ωi Ωi > と定義される。

グラッド [46] の方法に従い、１体分布関数を以下のように展開し、

f (x, c, ω, t)=(1 − ai∂
ai j ∂ 2
ai jk
∂3
+
−
+ ···)∂ci
2! ∂ci ∂cj
3! ∂ci ∂cj ∂ck

×f (0) (c, ω, t)(5.15)
37

３次まで残したものを使って速度の各モーメントを計算すると、Dt = ∂t +u·∇
として回転自由度も含めた連続方程式

Dt ρ=−ρ∇ · u





 ρDt u=−∇ · P









nIDt ω=−∇ · L + χ(Iω)ρDt T =−P · ∇u − ∇ · q − χρDt Θ=−L · ∇ω − ∇ · w − ξ − ω · χc (Iω)

(5.16)となる。但し、ρ は質量密度であり、圧力テンソル P はP= ρ < CC > +χk [C](5.17)

であり、回転自由度に起因する圧力テンソル L は

nd2< CΩ > +χk [IΩ]16

(5.18)

1
1
ρ < C 2 C > +χk [ C 2 ]22(5.19)L =である。また、熱流 q は

q =であり、回転自由度に起因する熱流 w は

w

=

nσ 2
I
ρ < Ω2 C > +χk [ Ω2 ]16
2

(5.20)である。並進運動のエネルギー散逸率 χ はχ1= χc [ c2 ]2

(5.21)

であり、回転運動のエネルギー散逸率 ξ はξ= χc [

σ2 2
ω ]
16

(5.22)

である。

5.2

極性流体

(5.16) の方程式系は極性流体力学と呼ばれる。極性流体とは元来分子極性
等の内部自由度を連続体の方程式の一部として採用し、その振る舞いを記述する方法である。こうした考え方はむしろ弾性体の静力学において導入され [47, 48]、 粉体の静力学にも適用されている。コンピューターのディスプレー等に広く使われている液晶を例にとって説明しよう [49]。液晶は細長い
ほぼ剛体的な分子からなる系と考えられる。細長い剛体棒を高密度に充填す
ると, 棒が平行に揃った状態が安定になる。この様な液晶状態にある分子集団38が流れる場合は等方的な相や等方的な流体とは異なった流動特性を示す。連続体の式の上でも液晶流体では, 分子の配向ベクトルを含む異方的な構成方
程式が必要になる。また分子回転に伴うトルクが重要な役割を果たすのでス
トレステンソルの反対称部分が分子回転の自由度とカップルして残る。また
この系では角運動量保存の式が運動方程式と独立となり局所回転密度の時間変化は偶 (ぐう) 応力 (おうりょく)( couple stress) を通して発展する。極性流体あるいは極性弾性体はこのように極性をもった構成要素が何らかの意味で巨視的な運動方程式とカップルする場合に一般的に有効な枠組である。その結果、弾性体ではポアソン比 (Poisson s ratio) が負になり得るなどの顕著な影響があり [50, 51]、粉体が静止状態から破断する際にも局所的な回転を通して極性弾性体的な扱いが有効であると考えられる。一方で、構成要素を有限に留めて連続体近似を行う手法の概念的な分かりにくさや構成方程式の複雑化から来る解析の面倒さ、各具体例で考える際には一般論を特に必要としないこと等から極性連続体の考え方が浸透しているとは言い難い。
気体論を用いてスピン自由度を極性流体として特徴づける考え方は古くから Dahler 等によって試みられ [52–54] 非圧縮流体 [55] やその中のサスペン
ションのレオロジー [56] も議論されている。粉体流への適用は金谷によるもの [57] や Lun によるもの [58] が知られ、SPH 法による有効性の検証 [59] も試みられている。また境界ではスピンや偶応力テンソルが重要であることはジェンキンス（Jenkins）によって指摘されている [60]。最近、御手洗等 [61]が斜面流に極性流体を適用し、スピン効果を見事に定量的に説明することに成功した。しかし彼女等の解析を通して明らかになったことは、スピン自由度が流体の回転と独立に振る舞うのは境界近傍に限定され、また流体そのもの運動にはスピン自由度の影響は殆ど現れないというものであった。またバ
ビック（Babic）によれば [62], 球対称な粒子では衝突による偶応力テンソルは消えてしまうということであり、スピン自由度の影響はますます小さいことが明らかになっている。一方で、ごく最近にも [63] でも回転自由度は等分配則を充たさないことが示されたり、剪断流下での極性流体の線形安定性解析がなされて [64] おり,

その回転効果を定量的に議論することの必要性も増している。

5.3

回転自由度の繰り込み

しかしながら、回転自由度を含めた連続方程式 (5.16) は自由度の数が多く、実際には解析が困難である。また、既述のようにバルクでの影響は殆どない。
そこで、ジェンキンス（Jenkins）とザング（ Zhang）によって提唱された、
回転自由度を動摩擦係数 µ の展開という形で法線方向の跳ね返り係数に繰り込む理論 [65] を用いて、連続方程式を簡単化する。ここではその２次元版である [66] に従い、回転自由度を法線方向の跳ね返り係数に繰り込むことに

39

する。
彼らによると、並進運動のエネルギー散逸率 χ および回転運動のエネルギーΘ
散逸率 ξ は粉体温度と回転温度の比（一般に< 1）で展開することができ、
TX1Θχ = −χc [ c2 ] =an ( )(5.23)
2
T
n
および

ξ

=

−χc [

X
σ2 2
Θ
ω ]=
bn ( )
16Tn(5.24)とできる。さらに、Herbst [67] らは χ および ξ をマクスウェル分布で計算し、

1
χc [ c2 ] ∝2σ2 2χc [ ω ] '16

ν 2 g(ν)T 3/2 (−1 + e + χrot )(5.25)

0(5.26)であることを導いている。ここで、χrot は回転の自由度を含む項である。従っΘて、まず (5.24) 式で ξ = 0 と置いてについて解くと、TΘT= (1 + e)µ +(1 + e)2 (β0 − 2) 2
µ + O(µ3 )β0 + 1(5.27)となる。ここで、β0 は接線方向の跳ね返り係数で、粒子が衝突するときの入射角によって値を変えるため、この時点では T と Θ の関係は一意には決まらない。ところが、この関係式を並進運動のエネルギー散逸率 χ に代入すると、

1n2 g(ν) −1
1Tχc [ c2 ][]( )1/22322(2πT ) 2πΘ
16T Kπ ΘΘΘΘ= −1 + e + Q + M + ( )2 N + O(( )3 )TTT

(5.28)となる。ここで、

8π √
T πΘ
σ
Q = −µ + (1 + e)µ2 + O(µ3 )

K

M
N

=

= µ + O(µ3 )3= − µ + O(µ3 )2(5.29)(5.30)(5.31)(5.32)である。従って、上式の結果をマクスウェル分布で計算した結果と比べると、

χrot

= Q+

Θ
Θ
M + ( )2 NTT

(5.33)

であるから、有効跳ね返り係数を

eef f= e+Q+ΘΘΘ
M + ( )2 N + O(( )3 )
T
T
T
40

(5.34)

Θ
の２次まで繰
T
り込むことができる。また、動摩擦係数が小さい場合に、eef f を O(µ2 ) まで残すととすることで、回転の効果を法線方向の跳ね返り係数の中に

eef f' e − µ + 2(1 + e)µ2(5.35)となる。この時点で接線方向の跳ね返り係数が、衝突時における粒子の入射角によって変わることによる不定性は完全に省かれている。
本研究では、回転自由度を繰り込んだ跳ね返り係数 eef f を使うことで、ジェンキンス（Jenkins）とリッチマン（Richman） [9] により与えられた連続方
程式を簡単化して解析を行った。

5.4

回転自由度を繰り込んだ連続方程式

5.4.1

連続方程式

連続方程式の解析を簡単にするため、(5.16) 式のうち、角速度および回転温度の発展方程式は考えずに、並進自由度の発展方程式のみ考える。このとき、回転自由度の効果は法線方向の跳ね返り係数に繰り込むため、法線方向の跳ね返り係数を有効跳ね返り係数 eef f に置き換える。従って回転自由度を繰り込んだ連続方程式は Dt ρ=−ρ∇i ui
ρDt ui =−∇j PijρDt T =−Pij (∇i uj ) − ∇i qi − χ(5.36)となる。ここで、Pij は圧力テンソルであり、静水圧 p、粘性率 η 、ずり粘性率 µ を用いて、

Pij= [p − η(∇ · u)]δij − µD̂ij(5.37)

である。但し、δij はクロネッカーのデルタであり、D̂ij = (∂i uj + ∂j ui )/2 である。また、q は熱流であり、熱伝導率 κ を用いて、

q = −κ∇T − λ∇ρ(5.38)

である。但し、密度勾配の係数 λ は密度と粉体温度の関数である。また、χ
はエネルギー散逸率であり、具体的には

χ=r
1 − e2 2
π 1/2√ ρ g(ν)T 1/2 [4T − 3
T (∇ · u)]
2
4 2π

である。

41

(5.39)

5.4.2

スケーリング本研究では、長さを粒子直径 σ 、速度を U/2 でスケーリングした。従って、時間、位置、速度、粉体温度は次のように無次元化される。

t

=

2σ ∗
t
U

x = σx∗

u=

U ∗
u
2

T =

U2
θ
4d

但し、t∗ 、x∗ 、u∗ 、θ はそれぞれ無次元化された時間、位置、速度および粉体温度であり、d は次元である。ここから、上付きの ∗ によって無次元である

ことを表すことにする。次に、連続方程式の各項を無次元化する。まず、ラグランジュ微分 Dt は

Dt

=
=

U ∂
(
+ u∗ · ∇)
2σ ∂τ
U
Dt∗
2σ

と無次元化される。次に、静水圧 p、粘性率 η 、ずり粘性率 µ、熱伝導率 κ、
熱流の一部である密度勾配の係数 λ を無次元化すると、

ρp U 2
p(ν)θ
4
σρp U
η =
η(ν)θ1/2
2
σρp U
µ =
µ(ν)θ1/22
κ = σρp U κ(ν)θ1/2
p =(5.40)(5.41)(5.42)(5.43)

3

λ =

σρp U
λ(ν)θ3/2
8(5.44)となる。ジェンキンスとリッチマンによって求められた右辺の各関数は表 5.1にあげておく。これらの関数を用いると、圧力テンソル Pij は

Pijρp U 2 ∗Pij
4=(5.45)となり、無次元化された圧力テンソルPij∗∗= [p(ν)θ − η(ν)θ1/2 (∇∗ · u∗ )]δij − µ(ν)θ1/2 D̂ij

(5.46)

が得られる。熱流 qi は

q =

ρp U 3 ∗q
8(5.47)となり、無次元化された熱流q∗= −κ∗ (ν)θ1/2 ∇∗ θ − λ∗ (ν)θ3/2 ∇∗ ν

(5.48)

が得られる。エネルギー散逸率 χ はχ=ρp U 3 ∗χ
8σ
42(5.49)

表 5.1: 無次元化された各関数。

p(ν)

=

1
ν[1 + (1 + e)νg(ν)]
2

η(ν)

=

1
√ (1 + e)ν 2 g(ν)
2π

µ(ν)

=

r

κ(ν)

=

λ(ν)

=

√

π
1
(1 + e)(3e + 1)
(1 + e)(3e − 1)
1[g(ν)−1 +ν+(+ )(1 + e)ν 2 g(ν)]
2 7 − 3e4(7 − 3e)8(7 − 3e)π2π[

13(2e2 + e + 1)9(1 + e)(2e − 1)
1g(ν)−1 +ν+(+)(1 + e)ν 2 g(ν)](1 + e)(19 − 15e)8(19 − 15e)32(19 − 15e)2πrπ 3e(1 − e)
1 d(ν 2 g(ν))−[4g(ν)−1 + 3(1 + e)ν]2 16(19 − 15e)νdνとなり、無次元化されたエネルギー散逸率
∗

χ

=

r
π 1/2 ∗ ∗1 − e2 2
1/2√ ν g(ν)θ [4θ − 3θ (∇ · u )]2
4 2π

(5.50)

が得られる。
従って、これらを連続方程式に戻すと、ρp や U などといった次元を持つものは両辺で相殺するため、無次元化された連続方程式は以下のようになる。

Dt ν=−ν∇ · uνDt ui =−∇j Pij(5.51)
 1 νDt θ=−Pij (∇i uj ) − ∇ · q − χ2

但し、簡単のため無次元であることを表す ∗ は省いた。

43

第 6 章 定常状態の解析

6.1

定常解

本章では (5.51) 式で導いた流体方程式系を用いて、単純せん断流の定常解を求めてみる。定常状態では、系の対称性から各流体変数はせん断方向（x方向）への依存性がなく、壁に垂直な方向（y 方向）にのみ依存し、

ν = ν(y)ux = ux (y)uy = 0θ = θ(y)(6.1)と書ける。また、他の全ての関数も y のみの関数で書ける。

uy = 0 となるので、流体変数 X = (ν, ux , uy , θ)t は
Dt X = ∂t X = 0

(6.2)

∇·u=0

(6.3)

を満たし、さらになので、質量保存の式は自動的に満される。従って、定常状態での連続方程式で残るのは運動方程式とエネルギーの式のみで、
 0=∇y Pxy0=∇y Pyy
0=−Pyx ∇y ux − ∇y qy − χ(6.4)となる。ここで、(6.20) 式の第１、２式（運動方程式）から、定常状態では
圧力テンソル Pxy (= P yx), Pyy が空間的に一定であることが解る。

Pxy = τ(一定)

Py = p(一定)となる。従って、Pxy , Pyy の表式から
µ(ν) 1/2 dux


−
θ
=τ
2
dy

 1 ν[1 + (1 + e)νg(ν)]=p
2

(6.5)

(6.6)

という関係が得られる。これより、せん断方向の速度勾配と温度を ν で表すことができ、

dux
dy
θ

2τ
µ(ν)θ1/2p2p=ν[1 + (1 + e)νg(ν)]
α(ν)

= −

(6.7)

=

(6.8)

44

となる。但し、α(ν) ≡

2
ν[1 + (1 + e)νg(ν)](6.9)

である。従って、これらを (6.20) の第３式（エネルギー方程式）に代入すれば、ν のみの微分方程式が得られる。
ここで、(6.20) のエネルギー方程式の各項はPyx ∇y ux

qy

= −τ

2τ
µ(ν)θ1/2

τ 2 α(ν)1/2= −2 1/2
µ(ν)
p
κ(ν)
1
λ(ν)
= −p3/2 p∇y () − p3/2∇y να(ν)α(ν)3/2α(ν)=(6.10)dp3/2
κ(ν)[(1 + (1 + e) (ν 2 g))− λ(ν)]∇y ν(6.11)dνα(ν)2α(ν)3/2(1 − e)η(ν)θ3/2
η(ν)= p3/2 (1 − e)α(ν)3/2

χ =(6.12)となるので、エネルギー方程式は

d
d
1
κ(ν)
dν
[(1 + (1 + e) (ν 2 g))[− λ(ν)] ]dy α(ν)3/2
dν
2α(ν)
dy
=²

2α(ν)1/2
η(ν)− (1 − e)µ(ν)α(ν)3/2という ν に関する微分方程式となる。但し、
µ ¶2
τ
² =
p

(6.13)

である。

1
1
d
κ(ν)
[( + r (ν 2 g))− λ(ν)]dνα(ν)α(ν)3/2 2(6.14)

η(ν)2α(ν)1/2− (1 − e)µ(ν)α(ν)3/2(6.15)

F (ν)

=

G(ν)

= ²

として、ν に関する微分方程式をddν
[F (ν) ] =
dy
dy

G(ν)

d
と置く。いま、 dνH(ν) = F (ν) とすると、

d
H(ν)
dy

dH(ν) dνdν dydν= F (ν)dy=

45(6.16)

なので、H(ν) を使って、ν に関する微分方程式が H(ν) に関する２階の常微分方程式となる。

d dH(ν)
(
) =
dy
dy
この両辺に

G(ν)

dH(ν)
dy ∗ を掛けて、バルク中央（y = 0）から任意の位置 y まで積

分すると

dH(ν)
1 d dH(ν) 2() = G(ν)
2 dy
dy
dy
Z y
1 dH(ν) 2
dH(ν) 2dH(ν) 0
[(
) −(
) ] =
G(ν)
dy
2
dy
dy y=0dy 00となる。左辺において、

dν
|y=0
dy

=

0

(6.17)という条件を加えることで、(dνdH(ν) 2= (F (ν) |y=0 )2 = 0)
dy y=0dyとなる。また、右辺は積分変数を変換してZ H(ν(y))
Z ydH(ν) 0dy=
G(ν)dH(ν)G(ν)
dy 0
H(ν(0))0Z ν(y)=F (ν 0 )G(ν 0 )dν 0ν(0)dとなる。最後の式変形では、 dνH(ν) = F (ν) より dH(ν) = F (ν)dν となることを用いた。
以上により1 dH(ν) 2()
2 dy ∗
つまりとなるので、=Z ν(y)F (ν 0 )G(ν 0 )dν 0

ν(0)

s Z
ν(y)
dH(ν)
2
=
±
F (ν 0 )G(ν 0 )dν 0
dy ∗
ν(0)

(6.18)dH(ν)= F (ν) dνdy
dy により F (ν) に戻すと

F (ν)
±q R
dν = dyν(y)2 ν(0) F (ν 0 )G(ν 0 )dν 0となる。従って、これを再び 0 から y まで積分することにより、
Z ν(y)F (ν 0 )q R 0±dν 0 = yνν(0)2 ν(0) F (ν 00 )G(ν 00 )dν 00(6.19)(6.20)となり、ν を y の関数として求めることが出来る。また、(6.7) および (6.8)
式により、速度場と粉体温度も y の関数として求まる。

46

表 6.1: 定常状態において衝突している粉体粒子の内訳

6.2

２体衝突

97.6%

多体衝突

2.4%

DEM データと理論値の比較

上記の方法で定常解を数値的に求め、DEM による分子シミュレーションのデータと比較した。定常解を数値的に求める際に、ν(0) の値が不定値として残るが、ここでは算出した面積占有率 ν(y) をバルク全体で積分し、その値が

1
∆

Z ∆/2

ν(y)dy= ν̄

(6.21)−∆/2に十分近い値に収束するまで反復計算をしている。また、境界条件から決まる p
および τ の値はフィッティングパラメータとして使い、p = 0.10、τ = −0.0017とした。

粉体粒子の回転が無い場合
まず、粉体粒子の回転を排除した DEM による分子シミュレーションの定常状態を、回転自由度を除いた連続方程式系 (5.51) の定常解を上述の方法により求めたものと比較したのが図 6.1 、図 6.2、図 6.3 であり、それぞれ面積
占有率、粉体温度、せん断方向の速度場を比較したものである。これらはいずれも定量的な一致をみることができ、相分離した状態でのクラスター領域（高密度・低温）でも運動論により導かれた連続方程式が有効であることが解る。高密度・低温な領域でも運動論が適用できたことは、定常状態のクラス
ター領域（高密度・低温）内での粉体粒子は殆ど自由に飛んでおり（接触している粒子は全体の 1.014 パーセント）、また接触している粒子のうち、97パーセント以上が２体衝突であるという DEM による分子シミュレーションの実測値からも納得できる結果である。

粉体粒子の回転がある場合
次に、粉体粒子の回転を入れた DEM による分子シミュレーションの定常状
態を、回転自由度を繰り込んだ連続方程式系 (5.51) を上述の方法により求めたもの（つまり回転自由度を繰り込んだ跳ね返り係数 eef f = 0.798 を用いて計算した定常解）と比較したのが図 6.4 、図 6.5、図 6.6 である。なお、フィッティングパラメータ p と τ の値は図 6.1 で使用したものと同じで、p = 0.10、

τ = −0.0017 である。
47

表 6.2: 定常解のフィッティング値と DEM による実測値

p
τ

定常解

シミュレーション

0.10−1.7 × 10−3

2.5 × 10−3
−1.7 × 10−3

バルク中央における面積占有率のピークの値のずれが大きくなるものの、半定量的な一致をみている。また、粉体温度とせん断方向の速度場について
はいずれも定量的に一致する。これにより、回転する粉体粒子からなる粉体
せん断流の定常状態も、回転自由度を法線方向の跳ね返り係数に繰り込むことで半定量的に説明できたということがいえる。

6.3

定常状態での圧力について定常解を求める際に、圧力 p および τ はフィッティングパラメータとして扱ったが、その値は DEM による分子シミュレーションで実際に測定することもできる。ここでは [77–81] に従い、

Pαβ

=

n
1 X σ ij ij X i j(n F +Cα Cβ )
A
2 α β
i=1

(6.22)

i≤j

として圧力を測定した。但し、A は圧力を測定するセルの面積、nij は粒子 iの中心から粒子 j の中心に向かう方向の単位ベクトル、Fij は粒子 i が粒子 jに及ぼす合力、Ci は粒子 i の速度ゆらぎ Ci = ci − u である。

DEM による分子シミュレーションでも定常状態において p および τ は空間的に一様にはなることは確認していて、τ の平均値はフィッティングの値とほぼ同じになる。しかし、p の平均値はフィッティングの値と強度一桁異なる結果になってしまった（表 6.2）。この不一致の原因を明らかにすることは境界条件の議論も含めて今後の課題である。

6.4

弾性極限および強いせん断応力がかかる場合について弾性極限 (e → 1) では、被積分関数の分母のルートの中が負、すなわちG(ν) > 0 になってしまい、定常解が計算できない。また、² > 0.06 になる程強いせん断応力がかかる場合も、G(ν) > 0 になってしまい計算できない。
DEM においても、上述の極限においてはいずれもエネルギーが定常化せず、定常状態が得られない。従って、これらは非定常問題であり、新たな課題として残っている。

48

表 6.3: 高密度系の定常解のフィッティング値と DEM による実測値（回転が
無い場合）
定常解

p
τ

シミュレーション−3

9.3 × 10−4
−1.0 × 10−4

3.0 × 10−8.0 × 10−4

表 6.4: 高密度系の定常解のフィッティング値と DEM による実測値（回転が
ある場合）
定常解

p
τ

6.5

2.7 × 10シミュレーション−2

5.8 × 10−4

−9.0 × 10−4

−2.1 × 10−4

高密度粉体ガス系の振る舞い

平均面積占有率が ν̄ = 0.24 である高密度粉体ガス系のせん断流の面積占有率の時間発展の様子を図 6.7 に示す。ここではシステムサイズを保つために、
せん断方向のバルク幅を半分にすることで平均面積占有率を２倍にしている。
図 6.7 によると、２本のクラスターの間のプラグ領域に、新しく小さなクラスターが形成されているのが解る。これは２本のクラスターがバルク中央に移動する際に、プラグ領域の粉体粒子が２本のクラスターに弾かれて、中央
で粉体粒子同士が非弾性衝突を起こして滞留するためである。高密度粉体ガス系では、この様な２次的な現象も観測できる。また、粉体温度の時間発展の様子（図 6.8）と比較すると、クラスターが形成される領域で粉体温度が比較的に低くなることが解る。
図 6.11、図 6.12、図 6.13 および図 6.14、図 6.15、図 6.16 は粉体粒子の回転が無い場合および粉体粒子の回転がある場合の定常状態を、式 (6.21) がν̄ = 0.24 に収束するように計算した定常解と重ねてプロットしたものである。
高密度の場合でも式 (6.20) から計算される定常解により半定量的な一致を得ることができる。
なお、フィッティングパラメータ p、τ の値と DEM による粒子シミュレー
ションから実測した値は表 6.3、6.4 にまとめておく。

49

1
0.8
0.6
ǵ
0.4
0.2
0
-90

-60

-30

0
[

30

60

90

図 6.1: 回転を排除した DEM による分子シミュレーションの定常状態における面積占有率 ν （ドット）と回転自由度を除いた定常解（実線）の比較。定
量的に一致している。

12

8
ǰ
4

0
-90

-60

-30

0
[

30

60

90

図 6.2: 回転を排除した DEM による分子シミュレーションの定常状態における粉体温度 θ （ドット）と回転自由度を除いた定常解（実線）の比較。境界付近のゆらぎは大きいが、定量的に一致している。

50

0.8

0.4
WZ

0

-0.4
-0.8
-90

-60

-30

0
[

30

60

90

図 6.3: 回転を排除した DEM による分子シミュレーションの定常状態におけるせん断方向の速度場 ux （ドット）と回転自由度を除いた定常解（実線）の比較。定量的に一致している。

1
0.8
0.6
ǵ
0.4
0.2
0
-90

-60

-30

0
[

30

60

90

図 6.4: 粉体粒子の回転を入れた DEM による分子シミュレーションの定常状態における面積占有率 ν （ドット）と回転自由度を繰り込んだ定常解（実線）の比較。バルク中央のピークのずれが大きい。低密度な領域では定量的に一
致する。

51

16
12
ǰ

8
4
0
-90

-60

-30

0
[

30

60

90

図 6.5: 粉体粒子の回転を入れた DEM による分子シミュレーションの定常状態における粉体温度 θ （ドット）と回転自由度を繰り込んだ定常解（実線）の比較。境界付近のゆらぎは大きいが、定量的に一致している。

0.8

0.4
WZ

0

-0.4

-0.8
-90

-60

-30

0
[

30

60

90

図 6.6: 粉体粒子の回転を入れた DEM による分子シミュレーションの定常状態におけるせん断方向の速度場 ux （ドット）と回転自由度を繰り込んだ定常解（実線）の比較。定量的に一致している。

52

1

1


0.8



0.8

0.6

0.6

ǵ

ǵ
0.4

0.4

0.2

0.2

0
-90 -60 -30

0
[

30

60

0
-90 -60 -30

90

1

0
[

30

60

90

1


0.8



0.8
0.6

0.6
ǵ

ǵ
0.4

0.4

0.2

0.2

0
-90 -60 -30

0
[

30

60

0
-90 -60 -30

90

0
[

30

60

90

1

1


0.8



0.8

0.6

0.6
ǵ

ǵ
0.4

0.4

0.2

0.2

0
-90 -60 -30

0
[

30

60

0
-90 -60 -30

90

0
[

30

60

90

図 6.7: 高密度粉体ガス系における面積占有率 ν の時間発展。

53

5

5



4

4

3

3

ǰ

ǰ
2

2

1

1

0
-90 -60 -30

0
[

30

60

0
-90 -60 -30

90

5

0
[

30

60

90

5


4



4

3

3

ǰ

ǰ
2

2

1

1

0
-90 -60 -30

0
[

30

60

0
-90 -60 -30

90

5

0
[

30

60

90

5


4


4

3

3

ǰ

ǰ
2

2

1

1

0
-90 -60 -30

0
[

30

60

0
-90 -60 -30

90

0
[

30

60

図 6.8: 高密度粉体ガス系における粉体温度 θ の時間発展。

54

90

0.8

0.8




0.4

0.4

WZ 0

WZ 0

-0.4

-0.4

-0.8
-90 -60 -30

0
[

30

60

-0.8
-90 -60 -30

90

0.8

0
[

30

60

0.8



0.4

0.4

WZ 0

WZ 0

-0.4

-0.4

-0.8
-90 -60 -30

0
[

30

60

-0.8
-90 -60 -30

90

0.8

0
[

30

60

90

0.8




0.4

0.4

WZ 0

WZ 0

-0.4

-0.4

-0.8
-90 -60 -30

90

0
[

30

60

90

-0.8
-90 -60 -30

0
[

30

60

90

図 6.9: 高密度粉体ガス系におけるせん断方向の速度場 ux の時間発展。

55

0.08

0.08




0.04

0.04
W[

W[

0

-0.04

-0.04
-0.08
-90 -60 -30

0

0
[

30

60

-0.08
-90 -60 -30

90

0.08

0
[

30

60

0.08




0.04
W[

0.04
W[

0

-0.04
-0.08
-90 -60 -30

0

-0.04
0
[

30

60

-0.08
-90 -60 -30

90

0.08

0
[

30

60



0.04

0.04
W[

0

-0.04
-0.08
-90 -60 -30

90

0.08


W[

90

0

-0.04
0
[

30

60

-0.08
-90 -60 -30

90

0
[

30

60

90

図 6.10: 高密度粉体ガス系における壁に垂直な方向の速度場 uy の時間発展。

56

1
0.8
0.6
ǵ
0.4
0.2
0
-90 -60 -30

0
[

30

60

90

図 6.11: 高密度粉体ガス系。回転を排除した DEM による分子シミュレーショ
ンの定常状態における面積占有率 ν （ドット）と回転自由度を除いた定常解（実線）の比較。

4
3
ǰ2
1
0
-90 -60 -30

0
[

30

60

90

図 6.12: 高密度粉体ガス系。回転を排除した DEM による分子シミュレーショ
ンの定常状態における粉体温度 θ （ドット）と回転自由度を除いた定常解（実線）の比較。境界付近のゆらぎは大きいが、定量的に一致している。

57

0.4
0.2
WZ 0
-0.2
-0.4
-90 -60 -30

0
[

30

60

90

図 6.13: 高密度粉体ガス系。回転を排除した DEM による分子シミュレーショ
ンの定常状態におけるせん断方向の速度場 ux （ドット）と回転自由度を除いた定常解（実線）の比較。定量的に一致している。

1
0.8
0.6
ǵ
0.4
0.2
0
-90 -60 -30

0
[

30

60

90

図 6.14: 高密度粉体ガス系。粉体粒子の回転を入れた DEM による分子シミュ
レーションの定常状態における面積占有率 ν （ドット）と回転自由度を繰り込んだ定常解（実線）の比較。バルク中央のピークのずれが大きい。低密度な領域では定量的に一致する。

58

4
3
ǰ2
1
0
-90 -60 -30

0
[

30

60

90

図 6.15: 高密度粉体ガス系。粉体粒子の回転を入れた DEM による分子シミュ
レーションの定常状態における粉体温度 θ （ドット）と回転自由度を繰り込んだ定常解（実線）の比較。境界付近のゆらぎは大きいが、定量的に一致し
ている。

0.4
0.2
WZ 0
-0.2
-0.4
-90 -60 -30

0
[

30

60

90

図 6.16: 高密度粉体ガス系。粉体粒子の回転を入れた DEM による分子シミュ
レーションの定常状態におけるせん断方向の速度場 ux （ドット）と回転自由度を繰り込んだ定常解（実線）の比較。定量的に一致している。

59

第 7 章 Linear Stability
Analysis7.1

粉体せん断流の線型安定性解析

第３章のはじめに述べたように、粉体せん断流の解析には、連続体をベースにした数値解析などがある。しかし、粉体流の従う連続方程式を解析的に解く手段はなく、連続体をベースにした数値解析も実験的な役割の方が大きいのが現状である。そこで、粉体せん断流の場合に、通常の流体のせん断流に見られるような単調な解の不安定性を議論することで、粉体せん断流の性質を洗い出そうとする試みが数多く報告されている。本研究で扱ったような２次元粉体ガスのせん断流に関する線型安定性解析は、アラム（Alam）とノット（Nott） [68–70] をはじめ、バビック（Babic）やサーベッジ（Savage）な
ど [71–74]、多くの報告がなされている。多くの場合、バルク幅が無限である
として線型安定性解析をしているが、アラムらの報告 [68] では有限なバルク幅での線型安定性解析を行っており、バルク幅が引き起こす粉体せん断流の定常解の分岐について詳しい説明がある。さらに、最近は回転自由度のある
流体系の安定性解析も詳しく調べられている。
また、粉体ガスではなく、重力の影響も考慮した粉体流の線型安定性解析を行ったものに、斜面流の線型安定性解析の報告 [13, 15] などもある。
ここでは本シミュレーションで用いたパラメータを使用し、Lees-Edwards
境界における粉体せん断流に関する線型安定性解析をすることで、相分離した解の安定性を議論する。尚、本章の結果はほぼ文献 [70, 71] と一致する。

7.2

線型化

連続方程式には一様解ν = ν0 (一定)



du

x0


=
a (一定)


dy


uy0 =
0


θ
=
θ
(
一定)


s 0




2(1 − e)η0 1/2


θ0
 a =
µ0

60

(7.1)

も自明な解として存在するが、実際には相分離した解が得られる。ここでは一様解の不安定性を議論するために、線型安定性解析を行う。
まず、各流体変数を列ベクトルν(x, y, t)ux (x, y, t)

X(x, y, t) = uy (x, y, t)
θ(x, y, t)(7.2)

で表したとき、一様解 X0 (y) と微小な摂動 X̂(x, y, t) に分解すれば、

X(x, y, t) =X0 (y) + X̂(x, y, t)(7.3)

であるから、連続方程式に代入して X̂ に関して線型な項のみを残す。まず、Lagrange 微分 Dt = ∂t + u · ∇ は以下の様に線型化される。

(∂t + u · ∇)ν =(∂t + a∂x )ν̂ ν(∂ + u · ∇)u =(∂ + a∂ )û + aû
t
x
t
x x
y

ν(∂t + u · ∇)uy =ν0 (∂t + a∂x )ûy ν(∂ + u · ∇)θ =ν (∂ + a∂ )θ̂
t

0

t

(7.4)

x

ここで添え字の 0 は一様解の場合であることを表し、以下同様の表現を採用する。次に、各関数を線型化する。まず２体相関関数はg(ν)= g0 + gν ν̂

(7.5)と線型化される。但し、gν は g(ν) の１次の摂動 ν̂ の係数であり、一様解の値で決まる定数である。以下同様に、静水圧は

p = p0 + pν ν̂ + pθ θ̂

(7.6)

粘性率のうち面積占有率による部分は

η(ν)= η0 + ην ν̂

(7.7)

ずり粘性率のうち面積占有率による部分はµ(ν) =µ0 + µν ν̂

(7.8)と線型化される。但し、各関数の摂動の係数は表 7.1 にまとめておく。次に、
圧力テンソルを線型化する。圧力テンソルは、一様解の場合のもの Pij0 および摂動部分 Pˆij に分けることができ、

Pij= Pij0 + P̂ij

(7.9)

となる。ここで、

Pij0 =

Ã
p0
τ0

τ0
p0

!

0
Dij
=

61

Ã

0

a/2

!
a/2
0

(7.10)

表 7.1: 各関数の１次の摂動の係数gν =7ν0 − 25
16(ν0 − 1)3

pν =

1
1
θ0 [1 + (1 + e)ν0 g0 ] + ν0 θ0 (1 + e)(ν0 gν + g0 )
2
2

pθ =

1
[1 + (1 + e)ν0 g0 ]ν02

1+eην = √ (ν02 gν + 2ν0 g0 )
2π
µν =

r

π 16(ν0 − 1)(25 − 7ν0 ) (1 + e)(3e + 1)
(1 + e)(3e − 1)
1[++(
+ )(1 + e)(ν02 gν + 2ν0 g0 )]2
2 (7 − 3e)(7ν0 − 16)
4(7 − 3e)8(7 − 3e)πであり、摂動部分に関しては

P̂ij

D̂ij

1/2

1/2

0
= [p̂ − η0 θ0 ∇ · û]δij − µ̂θ0 Dij
1
1/2
−1/2
0
− µ0 θ0− µ0 θ0 Dˆij
θ̂Dij
2
Ã
!
1
∂x ûx
(∂
û
+
∂
û
)
y
x
x
y
2
=
1
∂y ûy2 (∂y ûx + ∂x ûy )(7.11)(7.12)である。従って、圧力テンソルの各成分の具体形は
1/2

1/2

1/2

1/2

P̂xx

= pν ν̂ + pθ θ̂ − (η0 + µ0 )θ0 ∂x ûx − η0 θ0 ∂y ûyP̂yy= pν ν̂ + pθ θ̂ − (η0 + µ0 )θ0 ∂y ûy − η0 θ0 ∂x ûx
1
1
1
1
1/2
−1/2
1/2
1/2
= − aµν θ0 ν̂ − aµ0 θ0θ̂ − µ0 θ0 ∂y ûx − µ0 θ0 ∂x ûy2422= P̂xyP̂xyP̂yxである。これにより、圧力テンソルの空間微分は以下の様に線型化される。

−∇j Pij

00
00= −Pxx∂x ûx − Pxy∂x ûy − (Pyx∂y ûx + aPˆyx ) − Pyy∂y ûy= (−p0 ∂x − τ0 ∂y )ûx + (−p0 ∂y − τ0 ∂x )ûy − aPˆyx
1 2
1
1/2
−1/2
=
a µν θ0 ν̂ + a2 µ0 θ0θ̂24+(−p0 ∂x − 2τ0 ∂y )ûx + (−p0 ∂y − 2τ0 ∂x )ûy62

(7.13)

次に、熱流を線型化する。熱流も一様解の場合のもの q0 と摂動部分 q̂ に分けることができ、

q = q0 + q̂

(7.14)となる。ここで、

q0

= 0

q̂ =

(7.15)

1/2
3/2−κ0 θ0 ∇θ̂ − λ0 θ0 ∇ν̂(7.16)である。従って、熱流の勾配は以下の様に線型化される。
3/2

1/2

−∇ · q = λ0 θ0 [∂x2 + ∂y2 ]ν̂ + κ0 θ0 [∂x2 + ∂y2 ]θ̂(7.17)

次に、散逸項を線型化する。散逸項も一様解の場合のもの χ0 と摂動部分 χ̂
に分けることができ、

χ = χ0 + χ̂

(7.18)となる。ここで、

χ0

=

χ̂ =

3/2

(1 − e)η0 (ν0 )θ0

3
3/2
1/2
(1 − e)ην θ0 ν̂ + (1 − e)η0 θ0 θ̂
2
r
3 π
−
(1 − e)η0 θ0 ∇ · û
4 2

(7.19)(7.20)である。
以上により、線型化された連続方程式は以下の様になる。

質量保存則

∂t ν̂

= −a∂x ν̂ − ν0 (∂x ûx + ∂y ûy )(7.21)

運動方程式∂t ûx1
1/2= ν0−1 ( aµν θ0 ∂y − pν ∂x )ν̂2
1−1/2+ν0−1 ( aµ0 θ0∂y − pθ ∂x )θ̂
4
1
1/2
1/2
+(ν0−1 [(η0 + µ0 )θ0 ∂x2 + µ0 θ0 ∂y2 ] − a∂x )ûx2
µ0 1/2)θ ∂x ∂y − a)ûy+(ν0−1 (η0 +2 063(7.22)

∂t ûy1
1/2= ν0−1 ( aµν θ0 ∂x − pν ∂y )ν̂2
1−1/2+ν0−1 ( aµ0 θ0∂x − pθ ∂y )θ̂4µ0 1/2+ν0−1 (η0 +)θ ∂x ∂y ûx
2 0
1
1/2
1/2
+(ν0−1 [(η0 + µ0 )θ0 ∂y2 + µ0 θ0 ∂x2 ] − a∂x )ûy2

(7.23)

エネルギー保存則

∂t θ̂

3/2

1/2

3/2

= ν0−1 [2λ0 θ0 (∂x2 + ∂y2 ) + a2 µν θ0 + (1 − e)ην θ0 ]ν̂
3
1
−1/2
1/2
1/2
+ (1 − e)η0 θ0 ] + a∂x )θ̂+(ν0−1 [2κ0 θ0 (∂x2 + ∂y2 ) + a2 µ0 θ0
2
2
r
3 π
+ν0−1 [(−
(1 − e)η0 θ0 − 2p0 )∂x − 4τ0 ∂y ]ûx
4 2
r
3 π
+ν0−1 [(−
(1 − e)η0 θ0 − 2p0 )∂y − 4τ0 ∂x ]ûy(7.24)
4 2

以上をまとめて行列の形に書けば

∂t X̂
ここで、


−a∂x

al∂y − (pν /ν0 )∂x
L = al∂ − (p /ν )∂
ν
0 y
 x
2
2
β1 (∂x + ∂y ) + α1

= LX̂−ν0 ∂x2q∂x + m∂y2 − a∂xq 0 ∂x ∂yw∂x − (4τ0 /ν0 )∂y(7.25)−ν0 ∂y0
q ∂x ∂y − aq∂y2 + m∂x2 − a∂xw∂y − (4τ0 /ν0 )∂x

0an∂y − (pθ /ν0 )∂x

an∂x − (pθ /ν0 )∂yβ2 (∂x2 + ∂y2 ) + α2 + a∂xであり、l, m, n, q, q 0 , α1 , α2 , β1 , β2 は表 7.2 の通り、一様解の値のみに依る関
数である。

7.3

安定性解析

7.3.1

無限系

流体変数の時間部分はラプラス変換、空間部分はフーリエ変換し、X(t, x, y)= X(0, 0, 0)eσt eikx eiky

(7.26)

として線型化された連続方程式を書き直せば

σ X̂ = LX̂
64(7.27)








表 7.2: 行列 L の成分における空間微分の係数

l=

1 −1
1/2
ν µν θ 0
2 0

m=
1/2

1 −1
1/2
ν µ0 θ02 0
1/2

q = ν0−1 (η0 + µ0 )θ0

q 0 = ν0−1 (η0 + µ0 /2)θ01−1/2
n = ν0−1 µ0 θ043
w = ν0−1 (−43/2

α1 = 2a2 l + (1 − e)ν0−1 ην θ03/2

β1 = 2ν0−1 λ0 θ0r

π(1 − e)η0 θ0 − 2p0 )231/2
α2 = a2 n + (1 − e)ν0−1 η0 θ02
1/2

β2 = 2ν0−1 κ0 θ0となり、行列 L は以下のようになる。
−iakx−iν0 kx ialky − i(pν /ν0 )kx −qkx2 − mky2 − iakx
L = 
 ialk − i(p /ν )k−q 0 kx ky
x
ν
0 y

−β1 (kx2 + ky2 ) + α1 iwkx − i(4τ0 /ν0 )ky−iν0 ky0−q 0 kx ky − a−qky2 − mkx2 − iakx
iwky − i(4τ0 /ν0 )kx


ianky − i(pθ /ν0 )kx

iankx − i(pθ /ν0 )ky22−β2 (kx + ky ) + α2 + iakx

ここで、x 軸方向にせん断がかかる系においては、時間とともに x 軸方向に
引き伸ばされた空間構造をとる [75]。従って、時刻 t = 0 での波数を kx0 , ky0
として、波数を

k = (kx0 , ky0 − kx0 t)(7.28)とする。つまり波数は時間の線型関数になる。これにより行列 L は

L = A0 + A 1 t + A2 t 2(7.29)と t の関数になる。ここで、行列 A0 は行列 L において kx → kx0 および

ky → ky0 としたものに等しく、行列 A1 と A2 は

0
0
iν0

1/2
−1
0
 −il
ν
µ
θ
k
q
kx0
0 0
y0
0
A1 = kx0 iν −1 p
0
q kx0
2qky0
 0 ν
−1
2β1 ky04iν0 τ0
−iw

65


0

−in 

iν0−1 pθ 

2β2 ky0





である。

7.3.2

0


2  0
A2 = kx0
 0

−β1

0
−m
0
0

0
0
−q
0


0

0 

0 

−β2

(7.30)

kx = 0 の場合

kx = 0 であるため、A1 = A2 = 0 となり、摂動はσ X̂ = A0 X̂

(7.31)となる。つまり、kx = 0 の場合は、せん断がかかっても時間的に変化しない
恒常的な不安定性をみることになる。但し、
00−iν0 ky0

2

−a
ilk
−mk
y0
y0
A0 = 
−i(p /ν )k
2
0
−qky0
ν
0 y0

2
+ α1 −4i(τ0 /ν0 )ky0
iwky0
−β1 ky0

0
inky0





 (7.32)−i(pθ /ν0 )ky0 2+ α2−β2 ky0である。一様解を不安定化させるモードを取り出すため、式 (7.32) が自明でない解を持つ条件

det[σI − A0 ] = 0(7.33)

を解く。但し、I は４行４列の単位行列である。(7.33) 式は σ の４次方程式

σ 4 + ω3 σ 3 + ω 2 σ 2 + ω1 σ + ω0=0(7.34)

であり、各係数はそれぞれ

ω3
ω2

ω1

ω0

2
= −(m + q + β2 )ky0
− α2=

(7.35)

4(mq + qβ2 + β2 m)ky0+(mα2 + qα2 + 4nτ0 /ν0 + wpθ /ν0 − pν )ky0

(7.36)

2+(pν α2 − (α1 + 4aτ0 /ν02 )pθ )ky0

(7.37)

4+[(4τ0 l/ν0 + mα1 )pθ − (4nτ0 /ν0 + mα2 )pν ]ky0

(7.38)

6
4= −mqβ2 ky0+ [(m + β2 )pν + (nv − mα2 )q − (β1 + mw/ν0 )pθ ]ky06
= −m(pν β2 − pθ β1 )ky0である。本シミュレーションで使用した値に合わせ、平均面積占有率は ν̄ =0.12、跳ね返り係数は e = 0.798 として式 (7.34) 数値的に解を求めたのが図

66

4G ȁ
0.2

4G ȁ
0.002

0
0

-0.2
-0.002-0.4-0.004-0.6-0.8-0.006-1-0.008

-1.2

-0.01

-1.4

-0.012
2
0

0

0.5

1
M[

1.5

0.05

0.1
M[

0.15

図 7.1: 決定方程式を数値的に求めた解の実部。右は小さな波数領域での拡大図。小さな波数領域では２つの独立なモードと１組の縮退した負のモードが得られ、波数が少し大きくなると、独立だったモードが縮退して負になる。
小さな波数領域で独立であったモードの内の一方はほとんど 0 であり、もう一方の独立な解が支配的なモードであることが解る。また、さらに大きな波数領域では縮退していたモードの縮退が解けるものの、いずれも負の値をとるため、一様解の不安定性には寄与しない。

Im(ȁ)
0.04
0.02
0
-0.02
-0.04
0

0.5

1
ky

1.5

2

図 7.2: 決定方程式を数値的に求めた解の虚部。４つのモードがそれぞれ対称に現れる。図中のギャップのある部分では σ が実数解になっている。

67

0. 2

2

WPUVCDNG

1.6

UVCDNG

1.2

M[

0.8

0.4




0
0

0.2

0.4

0.6

G

0.8

1

図 7.3: 平均面積占有率 ν̄ = 0.12 の場合の跳ね返り係数 e に対する不安定な波数領域。本シミュレーションで使用した跳ね返り係数 e = 0.798 ではすべての波数領域で一様解が不安定であることが解る。e → 1 の極限で、一様解の安定領域が拡がってくる。

7.1 である。このとき、小さな波数領域で正になる不安定モードが現れるのが解る。
次に、平均面積占有率を ν̄ = 0.12 に固定した場合に、跳ね返り係数 e の値に対する不安定モードの実部の値を波数を縦軸にしてマッピングしたものが
図 7.3 である。本シミュレーションで使用した跳ね返り係数 e = 0.798 では
すべての波数領域で一様解が不安定である。また、弾性極限 e → 1 で一様解の安定領域が拡がってくるという特徴が見られる。

次に、跳ね返り係数を e = 0.798 に固定した場合に、平均面積占有率 ν の値に対する不安定モードの実部の値を波数を縦軸にしてマッピングしたものが図 7.4 である。本シミュレーションで使用した平均面積占有率 ν̄ = 0.12 ではすべての波数領域で一様解が不安定である。また、希薄極限 ν̄ → 0 で一様解の安定領域が拡がってくるという特徴が見られる。

7.3.3

kx 6= 0 の場合

e = 0.798、ν = 0.12 の場合の不安定領域の時間発展を図 7.6 に示す。この場合、不安定領域はせん断方向の波数 kx の方向に引き伸ばされていく。

68

1.6

UVCDNG

0.00051.2

0.001

M[ 0.80.0015

00.00002WPUVCDNG

0.4

0.00001

0
0

0.3

ǵ

0.6

0.9

図 7.4: 跳ね返り係数 e = 0.798 の場合の平均面積占有率 ν̄ に対する不安定な波数領域。本シミュレーションの平均面積占有率 ν̄ = 0.12 ではすべての波数
領域で一様解が不安定であることが解る。ν̄ → 0 の極限で、一様解の安定領域が拡がってくる。

0.9


0.8





0.7
0.6


ǵ 0.5



UVCDNGWPUVCDNG

0.4

0.3



0.2




0.1
0
0

0.2

0.4

G

0.6

0.8

1

図 7.5: 不安定モードが現れる小さな波数 ky = 0.05 に対して、跳ね返り係数

e と平均面積占有率 ν̄ を軸にして不安定領域をマッピングしたもの。弾性・希薄極限で一様解の安定領域が拡がってくる。

69

0.2

0.2
V

V

0.1

0.1

M[ 0

M[ 0

-0.1

-0.1

-0.2
-0.2

-0.1

0
MZ

0.1

-0.2
-0.2

0.2

0.2

-0.1

0
MZ

0.1

0.2

-0.1

0
MZ

0.1

0.2

0.2
V

V

0.1

0.1

M[ 0

M[ 0

-0.1

-0.1

-0.2
-0.2

-0.1

0
MZ

0.1

-0.2
-0.2

0.2

図 7.6: e = 0.798、ν0 = 0.12 での不安定領域の時間発展。図中の t の値は時間を表す。島状の不安定領域がせん断によって引き伸ばされていく様子が
解る。

7.4

有限系

有限系では波数は

kn

=

2nπ
∆

(n = 0, 1, 2, · · · )(7.39)という離散的な値をとる。このとき、∆ が小さければ離散性の効果が顕著に現れ、波数 k → 0 の極限に存在しなかった線型安定な領域が存在するようになる（因みに kn = 0 は中立安定である）。従って、有限周期系にせん断をかけたときには一様せん断流が実現しうる。有限系での議論は、アラム（Alam）とノット（Nott）によって詳しい論文がある [68]。彼らの論文によると、バルク幅 ∆ の逆数をオーダーパラメータにして定常解が分岐することが報告されており、∆ が小さければ一様解が安定になるということも示されている。

70

第 8 章 緩和過程

8.1

緩和過程の数値解析の目的と方法

DEM による分子シミュレーションの結果、粉体せん断流は比較的に早い段階でエネルギー的にほぼ一定な準安定状態に達した後、壁からやや離れた
位置に２つのクラスターを形成し、バルク中央に向かって移動していき、バ
ルク中央での衝突・合体の後、相分離した状態で定常に至る。この間の各流
体変数の緩和の様子は第３章に示した通りである。第４章では定常状態の解析を行い、DEM による分子シミュレーションの結果を再現することができた。これは第５章で述べたように、定常状態においては連続方程式が面積占有率の微分方程式に落とすことができ、粉体温度と速度場を面積占有率の分布から一意に決めることができるためである。
しかし、粉体せん断流が定常状態に至るまでの緩和過程を解析するには、定常解を求めた時のように方程式を簡単化することが難しい。そこで、まずは連続方程式を数値的に解き、その結果が DEM による分子シミュレーショ
ンの結果とよく一致することを確認しておく。次に、数値的に解いた連続方程式の項のうち、緩和過程で見られたクラスターの形成など、粉体せん断流の特徴的な現象に大きく寄与する項を残していくことで、相分離に至るまでの緩和過程の様子を記述できる最も簡単化された連続方程式を求める。そして、このように半実験的に求められた連続方程式の解析解を求めるというのが本章の目的である。
流体の連続方程式の数値解法については数多くの参考書があり [83–87]、通常の流体に関しての解法は確立している。粉体流の場合も、連続方程式の形は違うものの、偏微分方程式を数値的に解く方法は通常の流体に対してとられている方法と特に変わる点はない。本シミュレーションでは、空間微分については空間刻みの２次の精度を保ち、時間発展については４次のルンゲ・クッタ法を用いている。また、境界条件については、粉体せん断流の場合について [88–92] に報告があり、壁と衝突する粉体粒子の運動量保存と壁でのエネルギー収支の関係を定式化してある。本研究で用いたバンピー境界の場合については、ヒューイ（Hui）らによって詳しく説明されている [93]。

71

8.2

連続方程式の差分化

第３章で与えられた、回転自由度を繰り込んだ連続方程式を差分化する。
また、各流体変数 u をせん断方向（x 方向）に平均化し、壁に垂直な方向（y方向）への変化のみ考える。空間微分は空間刻みを h として中央差分することで、打ち切り誤差を O(h2 ) の程度に揃えている。なお、本シミュレーショ
ンでは h = 1 にとってある。

du
dy
d2 u
dy 2

ui,j+1 − ui,j−12hui,j+1 − 2ui,j + ui,j−1
h2

=
=

(8.1)
但し、u の添え字の i は時間ステップを表し、j は空間ステップを表している。
連続方程式は２階の空間微分まで含んでいるので、空間微分を含む項をまとめて f (u, ∇y u, ∇2y u) とすると、∂t u= f (u, ∇y u, ∇2y u)(8.2)と書ける。従って、上述の差分化により

ui+1,j= ui,j + δtf (ui,j , ui,j−1 , ui,j+1 )(8.3)とすることで各流体変数の時間ステップを１つ進ませることができる。時間
発展には４次の Runge-Kutta 法を採用しており、時間刻み δt は 0.01 とした。

8.3

境界条件

境界条件は、ジョンソン（Johnson）とジャクソン（Jackson）によるもの [89]を簡単化して用いる。まず、境界面での粉体流の速度と壁の速度の差をスリッ
プ速度 usl = u − uw として定義する。ここで、uw は壁の速度である。スリッ
プ速度 usl の方向の単位ベクトルを t とし、バルク内に向かう方向へ壁と垂直に立てた単位ベクトルを n とすると、壁で衝突する粉体粒子と壁の運動量保存則を考え、−n · P · t =πφΩ(ν, θ)|usl |4が境界で成立する。但し、P は境界での圧力テンソルであり、

(8.4)
π
は２次元粉
4

πd2
に起因する。また、φ は壁の荒さを表すパラメー
4タであり、Ω(ν, θ) は壁と粉体粒子の衝突頻度で、境界での面積占有率と粉体
温度できまり、

体粒子の質量 m = ρpΩ(ν, θ) ∝72

νg(ν)θ1/2(8.5)

である [94]。ここでは比例定数を (8.4) 式におけるパラメター φ に押し込めて Ω(ν, θ) = νg(ν)θ1/2 とする。次に、境界におけるエネルギー収支を考え、

n · q = −usl · P · n − Γ(ν, θ)(8.6)

が境界で成立する。ここで、q は境界での熱流であり、Γ(ν, θ) は粉体粒子と壁との衝突によるエネルギー散逸率である。壁でのエネルギー散逸率は、壁の硬さを表すパラメータ Φ と衝突頻度と境界での粉体温度の積に比例し、ππ
Γ(ν, θ) =ΦΩ(ν, θ)θ = Φνg(ν)θ3/2
(8.7)
4
4
である。なお、本シミュレーションで使用したパラメータの値は φ = 0.20、

Φ = 0.24 である。これらの境界条件を
fb (u, ∇y u)

=

0

(8.8)

と書いたとき、差分化により

fb (ui,j , ui,j−1 , ui,j+1 ) =

0

(8.9)

と変形できる。これを両壁際での境界値 ui,0 、ui,N について解くことにより、

ui,0= fb1 (ui,1 , ui,2 )(8.10)ui,N= fb1 (ui,N −1 , ui,N −2 )(8.11)

として境界値を決定する。但し、hN = ∆ であり、j = 0 は y = −∆/2 での値を表し、j = N は y = ∆/2 での値を表している。つまり、各流体変数を離散化し、面積占有率は νi,j 、速度場は ux,i,j および uy,i,j 、粉体温度は θi,j で表すと、(8.4) 式は圧力テンソルの具体的な表式から

ux,i,0=ux,i,N=−ux,i,2 + 4φhνi,1 g(νi,1 )/µ(νi,1 )1 + 4φhνi,1 g(νi,1 )µ(νi,1 )
ux,i,N −2 + 4φhνi,N −1 g(νi,N −1 )/µ(νi,N −1 )1 + 4φhνi,N −1 g(νi,N −1 )µ(νi,N −1 )(8.12)(8.13)となる。ここで面積占有率は境界で連続であると考え、

νi,0= 2νi,1 − νi,2

(8.14)

νi,N= 2νi,N −1 − νi,N −2(8.15)とする。次に (8.6) 式は熱流の具体的な表式から

θi,0

θi,N

θi,1 λ(νi,1 )(νi,2 − νi,0 )κ(νi,1 )[φ(1 + ux,i,0 )2 − φθi,1 ]νi,1 g(νi,1 )+2hκ(νi,1 )θi,N −1 λ(νi,N −1 )(νi,N −2 − νi,N )=κ(νi,N −1 )[φ(1 − ux,i,N )2 − φθi,N −1 ]νi,N −1 g(νi,N −1 )+2hκ(νi,N −1 )=73

(8.16)(8.17)

となり、境界での粉体温度の値が決まる。また、壁を通して粒子の移動は無
いため、境界では uy = 0 であるから

uy,i,0= 0(8.18)

uy,i,N= 0(8.19)となる。

8.4

完全な連続方程式の数値解析結果

まず、回転自由度を繰り込んだ連続方程式 (5.51) を y にのみ依存するとして１次元モデルに縮約し、それを完全な形で差分化した場合の数値計算結果を示す。DEM による粒子シミュレーションにおいて、系の全エネルギーがほぼ一定になる直前のデータを多項式でフィッティングしたものを初期条件とする（図 8.1）。

DEM による粒子シミュレーションの面積占有率の時間発展が図 8.2 であり、回転自由度を繰り込んだ連続方程式を完全な形で数値的に解いて求めた
面積占有率の時間発展が図 8.3 である。いずれも等間隔な時間ステップでプロットしてあり、連続方程式を数値的に解いた結果が DEM による粒子シミュ
レーションの結果とよく一致することが解る。
また、回転自由度を繰り込んだ連続方程式を完全な形で数値的に解いて求めた粉体温度の時間発展は図 8.4、せん断方向の速度場の時間発展は図 8.5、壁に垂直な方向の速度場の時間発展は図 8.6 であり、いずれも図中の番号は同じ時間ステップを表す。これらを DEM による粒子シミュレーションの結果である第４章の図 4.6、図 4.7 と比較すると、いずれもよく一致することが解る。但し、壁に垂直な方向の速度場については、第４章の図 4.8 と比べると絶対値が一桁異なっている。

8.5

連続方程式の簡単化

回転自由度を繰り込んだ連続方程式のうち、η(ν) は微小な量として 0 にする。また、エネルギー散逸率のうち速度場の内積に依る項はほとんど無視できる。さらに、運動方程式（速度場の発展方程式）において移流項は無視す

74

0.4

0.16

0.2

0.12

WZ [

ǵ [

0

0.08

-0.2
0.04
-80

-40

0
[

40

-0.4-80

80

-40

0
[

40

80

-40

0
[

40

80

50

0.04

40

W[ [

ǰ [

0

30
20

-0.04
-80

-40

0
[

40

80

10
-80

図 8.1: 初期条件として採用した多項式（実線）。ドットは DEM による粒子シ
ミュレーションにおいて、系の全エネルギーがほぼ一定になる直前のデータ。
るとすれば、連続方程式は主要な項のみ残した形に簡約でき、

∂t ν
∂t ux
∂t uy
∂t θ

= −∂y (νuy )µ(ν) 1/2
θ ∂y ux ]= ∂y [
2= ∂y [µ(ν)θ1/2 ∂y uy − p(ν)θ](8.20)(8.21)(8.22)= −uy ∂y θ − ν −1 µ(ν)θ1/2 (∂y ux )2+2ν −1 µ(ν)θ1/2 (∂y uy )2 − 2ν −1 ∂y qy√− 2π(1 − e2 )νg(ν)θ3/2(8.23)となる。これを数値的に解いた結果は図 8.7、図 8.8、図 8.9、図 8.10 であり、
回転自由度を繰り込んだ連続方程式を完全な形で解いた場合とほとんど同じであることが解る。これにより、粉体せん断流の発展方程式は簡単化された形 (8.20)〜(8.23) で記述できることが解った。

75

1
0.8



0.6
ǵ
0.4







0.2
0
-90 -60 -30

0
[

30

60

90

図 8.2: DEM による分子シミュレーションの面積占有率の時間発展。



0.8

0.6

ǵ
0.4




0.2

0
-90

-60

0
[

-30

30

60

90

図 8.3: 連続方程式を完全な形で数値的に解いた場合の面積占有率の時間発
展。DEM による分子シミュレーションの緩和過程をよく再現している。

76

12
ǰ8
4



0
-90
-60 -30

0
[

30

 
60

90

図 8.4: 連続方程式を完全な形で数値的に解いた場合の粉体温度の時間発展。

DEM による分子シミュレーションの緩和過程をよく再現している。

0.8
0.4




WZ



0
-0.4
-0.8
-90 -60 -30

0
[

30

60

90

図 8.5: 連続方程式を完全な形で数値的に解いた場合のせん断方向の速度場の時間発展。DEM による分子シミュレーションの緩和過程をよく再現して
いる。

77

0.8
0.4
W[ 0




-0.4


-0.8
-90 -60 -30

0
[

30

60

90

図 8.6: 連続方程式を完全な形で数値的に解いた場合の壁に垂直な方向の速度
場の時間発展。DEM による分子シミュレーションの緩和過程をよく再現し
ている。



0.8
0.6
ǵ



0.4





0.2
0
-90

-60

-30

0
[

30

60

90

図 8.7: 簡単化した連続方程式を数値的に解いた場合の面積占有率の時間発展。連続方程式を完全な形で数値的に解いた場合とほぼ同じである。

78

12
8
ǰ
4
0
-90 -60 -30

0
[





30

60

90

図 8.8: 簡単化した連続方程式を数値的に解いた場合の粉体温度の時間発展。
連続方程式を完全な形で数値的に解いた場合とほぼ同じである。

0.8
0.4




WZ 0



-0.4
-0.8
-90 -60 -30

0
[

30

60

90

図 8.9: 簡単化した連続方程式を数値的に解いた場合のせん断方向の速度場の時間発展。連続方程式を完全な形で数値的に解いた場合とほぼ同じである。

79

0.8
0.4


W[ 0


-0.4



-0.8
-90 -60 -30

0
[

30

60

90

図 8.10: 簡単化した連続方程式を数値的に解いた場合の壁に垂直な方向の速度場の時間発展。連続方程式を完全な形で数値的に解いた場合とほぼ同じである。

80

第 9 章 まとめ

9.1

議論―粉体流におけるバグノールド則

傾きのある板の上に粉体をのせると、粉体は重力の影響を強く受け、互いに接触しながらほぼ層状に流れ落ちる。これを斜面流といい、粉体流の特徴的な現象の一つとして研究が盛んである [95–97]。１９５４年、バグノールド（Bagnold）は斜面近傍の流れを無視して考えた場合に、粉体流の接線応力 τが速度勾配 γ̇ の２乗に比例するというバグノールド則（Bagnold scaling）

τ ∝ γ̇ 2(9.1)

が成立することを提唱した [95]。これは通常の流体の接線応力が速度勾配に比例するのに対し、粉体流の大きな特徴の一つと考えられている。その後、ポリクエン（Pouliquen）によってバグノールド則は実験によって確認され [96]、シルベルト（Silbert）らによる DEM の粒子シミュレーションでも実証されている [97]。また、カルマン（Kurman）はバグノールド則の成立にはバー
ネットオーダー（Barnett order）まで必要であることを指摘し、サントス（Santos）、ガルゾ（Garzo）、ダフティ（Dufty）は粉体ガスの場合、一様状態とせん断下では輸送係数に違いが現れるものの、定常状態においてバグノー
ルド則が成立することを示すなど、バグノールド則をめぐっては様々な議論がある。
昨年、御手洗（Mitarai）と中西（Nakanishi）により、斜面流におけるバ
グノールド則をニュートン流体の構成方程式から導くという興味深い報告がなされた [14]。彼等によると、まずエネルギー方程式−∇y q + τ γ̇ − χ = 0(9.2)

において ∇y q ' 0 という仮定をおくことで、τ γ̇ = χ(9.3)という関係式を得る。これをニュートン流体の枠内で圧力、熱流、エネルギー81散逸率を流体変数で展開した方程式系

f1 (ν)σ −d T

p =

(9.4)

1/2 1−d

σ

T

1/2

τ

=

f2 (ν)mqy=−f3 (ν)m−1/2 σ 1−d T 1/2 ∇y T

(9.6)

f4 (ν)m(9.7)−1/2 −1−d

χ =

σ

T

γ̇

3/2

(9.5)

の (9.7) 式に代入し、(9.5) 式を用いることで、粉体温度に対してT =

f2 (ν)
mσ 2 γ̇ 2
f3 (ν)(9.8)という関係を得る。但し、f1 (ν), f2 (ν), f3 (ν), f3 (ν) はいずれも面積占有率 νの関数であり、T は粉体温度である。これを再び (9.5) に代入することでτ=

f2 (ν)3/2
mσ 2−d γ̇ 2
f3 (ν)1/2(9.9)

を得る。彼等は斜面流の場合を想定しており、面積占有率 ν は一様とするため、上式はバグノールド則 τ ∝ γ̇ 2 である。

一方、本研究の第６章で得られた粉体せん断流の定常解においても、定常状態で面積占有率 ν が一様であれば、エネルギー方程式 (6.16) よりs(1 − e2 )g(ν) 1/2γ̇ =νθ(9.10)µ(ν)という関係が得られ、これにより

τµ(ν)θ1/2γ̇2µ(ν)3/2= −γ̇ 22(1 − e2 )1/2 νg(ν)1/2= −(9.11)となり、バグノールド則 τ ∝ γ̇ 2 が導き出せるのである。ところが、既に第７章でも述べた様に、粉体せん断流の場合には面積占有率が一様な解は不安定であり、せん断方向の速度勾配 γ̇ は一定ではない。従って、粉体せん断流では (9.11) 式が成立せず、バグノールド則が破れるのである。

9.2

結論

本研究における DEM による粉体せん断流の分子シミュレーションの結果、
せん断面同士の距離が有限である場合に、系のエネルギーがほぼ一定となる
ような準安定状態を経て相分離した定常状態へと緩和する様子が解った。このことは、本論文には記述しなかった Lees-Edwards 境界による無限系の場
合でも同様な結果を得ている。

82

また、運動論から導出された連続方程式を定常状態の場合に解き、DEM による分子シミュレーションとの定量的な一致をみただけでなく、回転自由度
を法線方向の跳ね返り係数に繰り込むことが実際に良い近似になっていることを示すことができた。これにより、粉体せん断流の相分離した定常状態を、
高密度・低温な領域も含めて連続体として説明することができた。
さらに、有限系および無限系の粉体せん断流においては通常の流体のせん断に見られるような一様解が不安定であることを確認し、有限系の場合にはバルク幅がある程度小さければ一様解が安定化することも解った。これにより、有限系では法線方向の跳ね返り係数をパラメータにして一様解が相分離した解に遷移していくという体系的な捉え方ができるようになった。
また、回転自由度を法線方向の跳ね返り係数に繰り込んだ連続方程式を数値的に解くことにより、粉体せん断流の時間発展を定量的に捉えることに成功し、クラスターが形成されて相分離した定常状態に緩和するのに必要な発展方程式を絞り込むという段階に差し掛かっている。

9.3

展望と課題

連続方程式の定常解を DEM による分子シミュレーションから得られた定常状態とフィッティングさせるときに用いた静水圧の値 p が、DEM による分子シミュレーションでの実測値と大きく異なってしまったことは、境界条件の議論も含めて今後の課題である。
また、DEM による分子シミュレーションで定常状態に至らないような非定常状態における解析も今後の研究対象である。また、緩和過程の解析では、絞り込んだ発展方程式はまだ解析的に扱える範囲ではなく、今後その解析手法も含めて議論していかなければならない。

83謝辞
本研究を進めるにあたり、指導教員である早川尚男先生には理論全般および数値解析全般に関して手厚いご指導をいただき、ここに深く感謝申し上げます。また、早川先生と武末先生にセミナー等でご指導いただいた非平衡系の物理学に関する基礎的な知識が本研究の支えになりました。本当にありがとうございました。DEM の内容については当研究室の占部さん、河原田さんに適切なご指摘を数多くいただき大変感謝しております。さらに、DEM の基本的な内容に関して、同志社大学の下坂先生による解説を参考にさせていただきました。大変ありがとうございました。また、光藤さんには数値計算に関する部分や、違った視点から的確な助言を下さり、心より感謝しております。そして、同じ学年で研究を続けてきた西野くんは、研究だけでなく様々な面で励みになり、心より感謝しています。また、富田先生と武末先生と佐野先生をはじめ、人間環境学研究科のみなさんや、当研究室の OB である國仲さんと金さんには本研究の内容に関して適切な指示を頂き、大変感謝しております。また、御手洗さんと波多野さんには粉体流に関する理論的な議論をしていただき、大変感謝しております。最後に、今回宇宙フォーラム事業の一環として本研究に対して多くの援助をして下さったことに、関係者の皆様方に深くお礼申し上げます。

84

関連図書
[1] 早川 尚男 『散逸粒子系の力学』岩波講座 物理の世界 物理と数理
４ （岩波書店, 2003）[2] N. V. Brilliantov and T. Pöschel, Kinetic Theory of Granular Gases,
Oxford Univ. Press, (2004)[3] J. J. Brey, J. W. Dufty, C. S. Kim and A. Santos, Phys. Rev. E 58,
4638, (1998)[4] Goldhirsch, I., and G. Zanetti, Phys. Rev. Lett. 70, 1619, (1993)[5] Kawasaki, K. and I. Oppenheim, Phys. Rev. 139, A1763, (1965)[6] van Noije, T. P. C. and M. H. Ernst, Phys. Rev. Lett. 79, 411, (1997)[7] van Noije, T. P. C. and M. H. Ernst and R. Brito, Physica A 251, 266,
(1998)[8] van Noije, T. P. C. and M. H. Ernst, Phys. Rev. E 61, 1765, (2002)[9] J. T. Jenkins and M. W. Richman, Phys. Fluids 28, 3485, (1985)[10] C. K. K. Lun, S. B. Savage, D. J. Jeffrey and N. Chepurniy, J. Fluid
Mech. 140, 223, (1984)[11] C. K. K. Lun, J. Fluid Mech. 233, 539, (1991)[12] V. Garzo and J. W. Dufty, Phys. Rev. E 59, 5895, (1998)[13] N. Mitarai and H. Nakanishi, J. Fluid Mech. 507, 309, (2004)[14] N. Mitarai and H. Nakanishi, Phys. Rev. Lett. 94, 128001, (2005)[15] Y. Forterre and O. Pouliquen, JFM 467, 361, (2002)[16] Goldhirsch, I., Ann. Rev. Fluids Mech, 35, 267, (2003)[17] H. M. Jaeger, S. R. Nagel and R. P. Behringer, Rev. Mod. Phys. 68,
1259, (1996)[18] C. S. Campbell, J. Fluid Mech. 348, 85, (1997)85

[19] Tan, M.-L. and Goldhirsch, I. Phys. Fluids, 9, 856, (1997)[20] L. Popken and P. W. Cleary, Journal of Computational Physics, 155,
1, (1999)[21] H. Hertz, J. Reine Angew. Math. 92, 156, (1881)[22] L. D. Landau and E. M. Lifshitz, Theory of Elasticity (3rd Ed.) (Pergamon, New York, 1960)[23] A. E. H. Love, A Treatise on the Mathematical Theory of Elasticity(4th Ed.) (Dover, New York, 1944)[24] R. D. Mindlin, J. Appl. Mech. Trans. ASME 16, 259, (1949)[25] K. M. Aoki and T. Akiyama, Phys. Rev. Lett 77, 4166, (1996)[26] H. Kuninaka and H. Hayakawa, Phys. Rev. Lett 93, 154301, (2004)[27] H. Kuninaka and H. Hayakawa, JPSJ 72, 1655, (2003)[28] H. Hayakawa and H. Kuninaka, Phase Transition 77, 889, (2004)[29] M. Y. Louge and M. E. Adams, Phys. Rev. E 65, 021303, (2002)[30] W. Gddsmith, Impact, (Dover, 1960, 2001)[31] O. R. Walton and R. L. Braun, J. Rheol. 30, 949, (1986)[32] L. Labous, A. D. Rosato and R. N. Dava, Phys. Rev. E 56, 5717,(1997)[33] Stronge, W.J. Impact Mechanics. Cambridge University Press, Cambridge.
[34] N. Maw, J. R. Barber and J. N. Fawcett, Wear[35] N. Maw, J. R. Barber and J. N. Fawcett, Wear[36] F. G. Bridge, A. Hatzes and D. N. Lin, Nature 309, 333, (1984)[37] G. Kuwabara and K. Kohno, Jpn J. Appl. Phys, 2, 1230, (1987)[38] N. Brillianton, F. Spahn, J. M. Hertzsch and T. Pöschel, Phys. Rev. E,
53, 5382, (1996)[39] W. A. Morgado and I. Oppenheim, Phys. Rev. E, 55, 1940, (1997)[40] R. Ramirez, T Pöschel, N. Brillianton and T. Schwagers, Phys. Rev.
E, 60, 4465, (1999)86

[41] K. L. Johnson, Contact Mechanics, (Cambridge University Press, 1985)[42] L. Verlet, D. Levesque, Mol. Phys. 46, 969, (1982)[43] P. Résibois, M. De Leener, Classical Kinetic Theory of Fluids, (Wiley
1977)[44] R. N. Dave, A. D. Rosato and K. Bhaswan, Mechanics Reserch Communications, 22, 335, (1995)[45] A. W. Lees and S. F. Edwards, J. Phys. C : Solid State Phys. 5, 1921,(1972)[46] H. Grad, Comm. Pure Appl. Math. 2, 331, (1949)[47] E.Cosserat and F. Cosserat, Theorie des Crops Deformables (A. Herrmann, Paris 1960)[48] A.C Eringen, J. Math. Mech. 16, 1, (1966)[49] P. G. de Gennes and J. Prost, The Physics of Liquid Crystals: Second
Edition (Oxford Univ. Press, Oxford, 1993)[50] R. S. Lakes, Science, 235, 1038, (1987)[51] Sergey V. Dmitriev*, Takeshi Shigenari and Kohji Abe, J.Phys. Soc.
Jpn. 70, 1431, (2001)[52] D.W. Condiff and J. S. Dahler, Phys. Fluids, 7, 842, (1964)[53] B.J. McCoy, S.I. Sandler and J. S. Dahler, J. Chem. Phys., 45, 3485,(1966)[54] J. S. Dahler and M. Theodosopulu, Adv. Chem. Phys., 31, 155, (1975)[55] G. Lukaszewicz, Micropolar Fluid:(Brikhaeuser, Boston, 1999)

Theory and Applications[56] H. Hayakawa, Phys. Rev. E, 61, 5477, (2000)[57] K. Kanatani, Trans. Jpn. Soc. Mech. Eng. B, 45, 507, (1979)[58] C.K. K. Lun, J. Fluid Mech., 2333, 4539, (1991)[59] J. Kano, A.Shimosaka and J. Hidaka, J. Soc. Powder Technol. Jpn.,
33, 95, (1996)[60] J. T. Jenkins, J. Appl. Mech., 59, 120, (1992)

87

[61] A. Mitarai, H. Hayakawa, and H. Nakanishi, Phys. Rev. Lett, 88,174301, (2002)[62] , Int. J. Eng. Sci., 35, 523, (1997)[63] I. Goldhirsch, S. H. Noskowicz and O.Bar-Lev, Phys. Rev. Lett, 95,
068002, (2005)[64] B. Gayen and M. Alam, preprint[65] J. T. Jenkins and C. Zhang, Phys. Fluids 14, 1228, (2002)[66] D. K. Yoon and J. T. Jenkins, Phys. Fluids 17, 083301, (2005)[67] O. Herbst, M. Huthmann, and A. Zippelius, Granular Matter 2, 211,
(2000)[68] P. R. Nott, M. Alam, K. Agrawal, R. Jackson and S. Sundaresan, J.
Fluid Mech. 397, 203, (1999)[69] M. Alam and P. R. Nott, J. Fluid Mech. 377, 99, (1998)[70] M. Alam and P. R. Nott, J. Fluid Mech. 343, 267, (1997)[71] M. Babic, J. Fluid Mech. 254, 127, (1993)[72] S. A. Kinnas and N. E. Fine, J. Fluid Mech. 254, 151, (1993)[73] C. Wang, R. Jackson and S. Sundaresan, J. Fluid Mech. 308, 31, (1996)[74] S. B. Savage, J. Fluid Mech. 241, 109, (1992)[75] O. M. Phillips, Ann. Rev. Fluid Mech. 1, 245, (1969)[76] 日本流体力学会 編集 水島 二郎・藤村 薫 著 『流れの安定性』流体力学シリーズ５（朝倉書店, 2003）[77] F. da Cruz, Ph.D. thesis, Ecole Nationale des Ponts et Chaussees, 2004,
http://pastel.paristech.org/archive/00000946.
[78] M. Latzel, S. Lunding and H. J. Hermann, Granular Matter 2, 123,
(2000)[79] M. Prochnow, Ph.D. thesis, Ecole Nationale des Ponts et Chaussees,
2002, http://pastel.paristech.org/archive/00000321.
[80] O. R. Walton and R. L. Braun, J. Rheol. 30, 949, (1986)[81] J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic
Press, London, 1986).
88[82] Frederic da Cruz, Sacha Emam, Michael Prochnow, Jean-Noel Roux,and Francois Chevoir, Phys. Rev. E 72, 021309, (2005)[83] 山崎 郭滋 『偏微分方程式の数値解法入門』（森北出版株式会社, 1995）[84] 登坂 宣好 大西 和栄 『偏微分方程式の数値シミュレーション』第５刷（東京大学出版会, 1995）[85] 数値流体力学編集委員会編『圧縮性流体解析』（東京大学出版会,

1995）[86] 河村 哲也 『流体解析 I 』第６刷（朝倉出版, 2002）[87] スハス V パタンカー著 水谷 幸夫・香月 正司 共訳 『コンピュータによる熱移動と流れの数値解析』[88] J. T. Jenkins and M. W. Richman, J. Fluid Mech. 171, 53, (1986)[89] P. C. Johnson and R. Jackson, J. Fluid Mech. 176, 67, (1987)[90] M. Y. Louge, Phys. Fluids 6, 2253, (1994)[91] M. W. Richman, Acta Mechanica 75, 227, (1988)[92] J. T. Jenkins and M. Y. Louge, Phys. Fluids 9, 2835, (1997)[93] K. Hui, P. K. Haff, J. E. Ungar and R. Jackson, J. Fluids Mech., 145,
223, (1984)[94] Y. Forterre and O. Pouliquen, J. Fluid Mech. 467, 361, (2002)[95] R. A. Bagnold, Proc. Roy. Soc. London A, 225, 49, (1954)[96] O. Pouliquen, Phys. Fluids, 11, 542, (1999)[97] L. E. Silbert, D. Ertas, G. S. Grest, T. C. Halsey, D. Levine and S. J.
Plimton, Phys. Rev. E, 64, 053102, (2001)89