粉体の運動論モデルと流体力学的性質
第 28 回 Grain Formation Workshop, 平成 22 年度銀河のダスト研究会

矢野良輔、鈴木宏二郎
東京大学大学院先端エネルギー工学専攻

Sept. 3, 2010

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

1 / 10内容• 粉体（ダスト）を”流体”として考える上で、運動論的手法は重要となる．

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

2 / 10内容• 粉体（ダスト）を”流体”として考える上で、運動論的手法は重要となる．
• 通常、運動論的手法が適応できる領域は、希薄な粉体に限定される．濃い粉体の運動論は、未着手．

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

2 / 10内容• 粉体（ダスト）を”流体”として考える上で、運動論的手法は重要となる．
• 通常、運動論的手法が適応できる領域は、希薄な粉体に限定される．濃い粉体の運動論は、未着手．
• 希薄な粉体の運動は、非弾性ボルツマン方程式で記述される．

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

2 / 10内容• 粉体（ダスト）を”流体”として考える上で、運動論的手法は重要となる．
• 通常、運動論的手法が適応できる領域は、希薄な粉体に限定される．濃い粉体の運動論は、未着手．
• 希薄な粉体の運動は、非弾性ボルツマン方程式で記述される．
• 粉体の運動論モデルは、Brey ら (PRE, 1998) により幾つか提案されているが、流れ場ある下で、DSMC 法と比較していないために実用性は不明．弾性極限で、BGK モデルに漸近するので不十分．

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

2 / 10内容• 粉体（ダスト）を”流体”として考える上で、運動論的手法は重要となる．
• 通常、運動論的手法が適応できる領域は、希薄な粉体に限定される．濃い粉体の運動論は、未着手．
• 希薄な粉体の運動は、非弾性ボルツマン方程式で記述される．
• 粉体の運動論モデルは、Brey ら (PRE, 1998) により幾つか提案されているが、流れ場ある下で、DSMC 法と比較していないために実用性は不明．弾性極限で、BGK モデルに漸近するので不十分．
• Brey らのモデルとは別に、14 モーメント方程式から、
BGK(Bhatnagar-Gross-Krook) 型の運動論モデルを提案し Brey らのモデルと比較する．

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

2 / 10はじめに粉体は、非弾性衝突をすることからエントロピー開放系の物理として、昨今、多くの研究がなされている．工学的応用から理学的関心に至るまでその適応範囲は広い．
粉体はその大きさが、流れ場のスケールに対して無視できない場合も多く、そのような場合は、体積粘性などの効果が出現し、通常の流体の取り扱いとは少々異なる．ここでは、体積粘性が余り問題とならないようなできる限り薄い粉体を考える．即ち、Enskog 方程式でなく、Boltzmann 方程式を適応できる範囲の話に限定する．粉体同士の相互作用は、衝突のみとし、粉体の表面は極めて滑らかゆえに、粉体の回転運動は無視する．また、反発係数は、相対速度に寄らず一定である．
運動論方程式とは、メゾスケールの粒子運動を記述した方程式で、ここでは、Boltzmann方程式の簡易方程式を指す事とする．
運動論方程式を考える意義は、数学的性質が明確な方程式を使う事で現象を説明しやすくする事やモンテカルロ法に比べて数値揺らぎが無い事から、非定常な問題を解くのに適している．特に、揺らぎのような物を扱う際は、数値揺らぎと熱的揺らぎを区別せねばならず、モンテカルロ法による方法は、ある意味破綻する．

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

3 / 10粉体の運動論方程式
粉体は、非弾性衝突するので、ボルツマン方程式の衝突項は、非弾性衝突の効果を入れて書き換えなくてはならない．( は、反発係数))
Z (∂f (v)∂f (v)
1+ vi=f (v0 )f (v01 ) − f (v)f (v1 ) gσ sin χdχddv31 ,
∂t
2
∂xi
空間一様な系では、温度のみの関数となり、衝突項は次に書かれる．(ζ は冷却率)!
df
1 ∂vi fdT
= ζ. →= −ζT
2 ∂vi
dtdt
Brey らは、この性質から次のような運動論モデルを提案している．

∂f (v)∂t+ vi∂f (v)∂xi

 1 ∂Ci f= νel f (0) (T) − f + ζ2 ∂vi|{z} |{z }elastic

inelasticf (0) は、Maxwell-Boltzmann 分布ここで、Brey らのモデルは、ボルツマン方程式から直接導出したモデルではない事に注意されたい．非平衡条件下ではどうか？

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

4 / 10非弾性衝突ボルツマン方程式は、次に書きなおすことが出来る．
Z∂f (v)∂f (v)+ vi= ψ (G − f ) , ∧ψ =f (v1 )gσ sin χdχddv1 ,∂t∂xi
G を冷却率、輸送係数が合うように、14 モーメントで展開する．
!
1+G = f13
= A (0) + Ai (1) Hi(1) + Aij(2) Hij(2) + Ai (3) Hi(3) + aNE H(4) fM ,8(n)13 個の未定パラメタ A は、非弾性ボルツマン方程式から導かれる Grad の 13 モーメンiト方程式を満たすように決定される．
Mass conservation:

ZV3Momentum conservation:

ZV3

 +

ψ f13− f dv = 0, +
Ci ψ f13− f dv = 0,

Equation of pressure deviator:

Z

V3

Equation of heat flux:

 +Ci Cj ψ f13− f dv = −νpij + δij nνa(2) T,

Z
V3Ci C 22 +ψ f13− f dv = ν ($ − 1) qi ,矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

5 / 10総括：• Brey らのモデルは、非弾性衝突部分の取り扱いに於いて、非平衡系での正しい取り扱いが出来ているか不明である．少なくとも 13 モーメントの挙動は、合っていない．一方で、Yano らのモデルでは、13 モーメントの挙動は、ボルツマン方程式によるものと一致する
ものの、それ以上のモーメントの挙動に関しては、冷却の効果が考慮できていない．一方
で、Brey らのモデルは、13 モーメント以上のモーメントにもある程度冷却の効果が反映されていると期待できる．また、分布関数の正値性が常に保証される．一方で、Yano らのモデルでは、正値性が保証されない．
→ 数値計算で検証してみる．

Shock wave/shock wave

Heating wall
Tw=5Tп

Stagnation Streamline (SSL)
Figure: 円柱周りに形成される衝撃波の実験（左図）Y. Amarouchene et al., Phys. Fluids, 18, 031707(2006)、円柱周りに計算される衝撃波の解析（右図）

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 20106 / 10

反発係数 () による流れ場の変化〜DSMC 法の結果

• 反発係数は、 = 1.0, 0.975, 0.95, 0.925 と変化させる．
• Kn=0.1 として、一様流のマッハ数は、M = 8 とする．
• 逆べき分子として、F = −φ/rl として、l = 5, 10 を考える．κ = 5 は、マクスウェル分子と呼ばれる．

2
l= 5 1

Density
2-X/R

341051021011006Temperature 54Solid: ε=0.925
Dashed: ε=0.95
3Dash-Dot: ε=0.975
Dash-Dot-Dot:ε=1 2Density

1.5-X/R2

l= 1012.50320151050NORMALIZED TEMPERATURE

5Temperature 3157Velocity

NORMALIZED VELOCITY

1020

NORMALIZED DENSITY

157Velocity
6Solid: ε=0.925
ε=0.95
Dashed:5Dash-Dot: ε=0.975
Dash-Dot-Dot:ε=1 4

NORMALIZED TEMPERATURE

20NORMALIZED VELOCITYNORMALIZED DENSITY

Density, velocity and temperature profiles along SSL Density, velocity and temperature profiles along SSL• 外部からエネルギーを与え続けるので、定常状態が実現できる．
• 衝撃波背後の温度は、反発係数を小さくすると減少する．
• 衝撃波離脱距離は、反発係数を小さくすると減少する．
• 壁面で加熱しているので、κ = 10 ∧  = 0.95 では、壁面付近で再加熱領域が出来るが、
κ = 10 ∧  = 0.925 では、壁面付近で、温度は 0.12 に漸近する．

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

7 / 10 = 0.95, 0.975 の淀み流線上の密度、速度、温度分布〜運動論方程式の結果

12345017Velocity

6Solid: Present
5Dashed: BGK
Symbols: DSMC 415103l= 10
ε=1.0Density

2Temperature

51

1.52-X/R

2.53103.5

NORMALIZED VELOCITYNORMALIZED DENSITY/TEMPERATURE-X/R
Density, velocity and temperature profiles along SSL
2023408642-X/R
Density, velocity and temperature along SSL20718Velocity

166Solid: Present
5Dahed: Brey et al.
Symbols: DSMC 41412103l=10ε=0.975

86Density

42Temperature

1211.522.5

30-X/R

145Solid: PresentDahed: Brey et al.
Symbols: DSMC 4l= 53ε=0.95
Temperature 2Density
1121086462
11.5

22.5

3

3.5 4

NORMALIZED VELOCITY

210Velocity

018161412108642

NORMALIZED TEMPERATURE

141216-X/R
Density, velocity and temperature along SSL
1816141210864240Velocity

35306Solid: Present
Dahed: Brey et al. 4
Symbols: DSMC
l= 10ε=0.95
2TemperatureDensity25201510511.522.5

3NORMALIZED VELOCITY

5614718181614121086420

NORMALIZED TEMPERATURE

l= 5
ε=1.0
2
Temperature

816

NORMALIZED DENSITY

310NORMALIZED VELOCITY

Density

1218

NORMALIZED TEMPERATURE

10Velocity 6Solid: Present
Dahed: Brey et al.
Symbols: DSMC 4
l= 5ε=0.975
Temperature 2Density

14NORMALIZED VELOCITY

Solid: Present
5Dashed: BGKSymbols: DSMC 41516NORMALIZED TEMPERATURENORMALIZED DENSITY

6Density, velocity and temperature along SSLDensity, velocity and temperature along SSL
NORMALIZED DENSITY

Velocity

NORMALIZED VELOCITY720

NORMALIZED DENSITYNORMALIZED DENSITY/TEMPERATURE

Density, velocity and temperature profiles along SSL-X/R

Maxwell 分子 (l = 5) では、提案したモデルと DSMC 法による結果は良く合う．一方で、
l = 10 では、 = 0.95 で壁面や衝撃波背後の温度が合わない．

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

8 / 10淀み流線上の分布関数10-1l= 5ε=1.010-310-4-510-61010-110-2-3-3-1010-60cx/(2T∞ /m ) 0.5

l= 5ε=1.010-7-80100.5cx/(2T∞ /m)l= 5ε=0.975

1010-710

Distribution functions at -X/R=2.05
Solid: Present-1
10 Dashed: Brey et al.
Symbols: DSMC

10-6

l= 10
ε=1.0-10

10-70cx/(2T∞ /m)0.5

10-810Distribution functions at -X/R=2.05
10-1 Solid: Present

10-1

10-2 Symbols: DSMC

10-2-10-50cx/(2T∞ /m)0.5

1010-410-410-5510

Distribution functions at -X/R=1.75
l= 510-1 ε=0.95-210-31010-3-410-5

l= 5ε=0.95

10-1000.5cx/(2T∞ /m )Solid: Present
Dashed: Brey et al.
Symbols: DSMC

10-4
10-510-10-50cx/(2T∞ /m)

0.551010

Solid: Present
Dashed: Brey et al.
Symbols: DSMC

10-6
00.5cx/(2T∞ /m )10-1000.5cx/(2T∞ /m )10

Distribution functions at -X/R=1.75

Distribution functions at -X/R=2.05
Solid: Present
10-1
Dashed: Brey et al.
-210 Symbols: DSMC
1010-410-510-610-710-8-20.5

l= 10ε=0.975

10-5

l= 10ε=0.975-1010-110-2-3100cx/(2T∞ /m)-31010-6l= 10
ε=1.0-10

Distribution functions at -X/R=1.75-3Solid: Present
Dashed: Brey et al.
Symbols: DSMC

10-4
10-5

Dashed: Brey et al.

10-4

l= 5ε=0.975

10-6

Solid: Present
Dashed: BGK
Symbols: DSMC

10-5



10-610Solid: Present
Dashed: Brey et al.
Symbols: DSMC

1010-4

Distribution functions at -X/R=1.75

10-210-510 Dashed: BGK
Symbols: DSMC
-3-210-5-10

10-1

10-410-1 Solid: Present

10 Dashed: BGK
Symbols: DSMC
-310

10

Distribution functions at -X/R=2.05

Distribution functions at -X/R=2.49

10-1 Solid: Present-210-2
cx/(2T∞ /m)10

Distribution functions at -X/R=2.89

Distribution functions at -X/R=2.49

Distribution functions at -X/R=2.89
-1

10-2
10-3
10-4
10-510-610-7
10 Solid: Present-8
10 Dashed: BGK
Symbols: DSMC
-910-100l= 10ε=0.95

10-310-4
10-510-6l= 10ε=0.95-1000.5cx/(2T∞ /m )矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

10Solid: Present10-7 Dashed: Brey et al.
-8 Symbols: DSMC
10-100cx/(2T∞ /m)

0.510Sept. 3, 2010

9 / 10議論♠ BGK 型の運動論モデルを提案し、弱非弾性領域においてその性質を解析した．

♠ 非平衡領域の冷却率への効果が、速度依存型の衝突頻度を持つ逆べき分子の計算から確認された．
♠ 分布関数の正値性を満たさないために、反発係数を小さくすると分布関数の負の速度テイルに分布関数が負になる領域が確認された．
♠ 非弾性ボルツマン方程式における生成項の効果は、反発係数を下げると、負の速度テイルが弾性衝突の場合に比べて、衝撃波隆起位置で、過剰分布となる事が確認された．
♠ 今後、シェア流などに適応して、非弾性衝突による流れの不安定性の解析をする予定である．

矢野良輔、鈴木宏二郎 (東京大学大学院先端エネルギー工学専攻) 粉体の運動論モデルと流体力学的性質

Sept. 3, 2010

10 / 10