大阪大集中講義
Subject ：散逸粒子系の統計力学
Subtitle： 散逸粒子を例にした非平衡統計物理学入門

Speaker ： 早川尚男 (Dept. Phys., Kyoto Univ.)
Date： 2004年6月16日-6月18日

Location ： 大阪大学

粉体の新しい本 （早川尚男）• 粉体とは何だろうか• 非弾性衝突• 気体論的アプローチ• 遅い流れ• 静力学• 散逸粒子系の物理の私的展望目次（１）• I．粉体物理とは– I-A. 粉体の様々な現象– I-B. 数値計算法• II．非平衡の現状：粉体の果たす役割– II-A. 非平衡の本質： 不可逆性について– II-B. Langevin方程式とFokker-Planck方程式

目次（２）• III. Boltzmann方程式とその応用– III-A. Boltzmann方程式と流体方程式– III-B. Burnett orderのBoltzmann方程式の解– III-C. 非弾性気体の統計力学– III-D. 実験で観測できる非弾性気体と非平衡統計• 非弾性衝突について（巨視的散逸性についての考察）– IV-A. 準静的理論– IV-B. 弾性円盤の衝突シミュレーション– IV-C. はねかえり係数が１を越える現象• まとめ

Ｉ．Introduction - What is granular physics?• Granular Materials： Inelastic Repulsive
Interactions, Excluding Volume Effects• Rich Examples in Daily Life: Important inIndustries• Interdisciplinary ： Soil Mech., Civil Eng., Chem.
Eng., Earth Sci., etc.
• Strong Fluctuation : Not easy to control• A Good Subject for Theoretical Physics• Important for Education of Physics

I-A.Physics of inelastic collisions ofmacroscopic materials
You do not know the details of inelastic collisions of macroscopic materials.
For example, the coefficient of restitution becomes easily larger than 1.

M.Y.Louge and M.E.Adams,
Phys. Rev. E, 65, 021303 (2002).

Shape of sand piles•••

Static shape is from the frictional
force. (Cf: No shape of fluid).
Force distribution is not uniform.
The shape and the force
distribution depend on history ofmaking sand piles.
Pressure distribution obeys an
exponential function.

Force propagates along force
chains• Strong fluctuations• Concentration ofstresses in small
number of particles• collapse of a silo by
stress concentration!

Force chains (2)• Force chains can be observed in an arrayof photo-elastic disks (below).

Frictional flow in seeds of mustard (photo by
Chicago group)•••

coexistence of solid-likeregion and fluid-likeregionSeeds in solid-likeregion has very slow
creep motion.
Fluid-like region obeys
Bagnold’s law (shear
stress is proportional to
the square of velocity
gradient).

Granular Gases (What can we obtain if molecules ingases are dissipative?)（１） The system is a model ofintersteller dusts.
（２）Theoretical analysis is possible,

Goldhrisch andZannetti (1993 PRL)

粉体振動層のパターン形成
An oscillon observed in the experiment in Texas group

Typical Patterns

Multiple oscillons

Surface patterns in vibrating beds

Typical Patterns (2)Multi-rolls convections

Segregation

Radial segregation
Axial segregation

粉体振動層の準２次元的挙動• 左図：鵜川、佐野（東京農工大）による実験• 弾性板の振動と類似• 半定量的に有効性は確認されている。
• 斥力系で弾性的挙動
は興味深いUgawa and Sano

実験で観測されたもの（２）I-B.Standard method of simulation for granularparticles (DEM)• DEM(=discrete or distinct element method) is astandard technique to simulate granular assemblies.
• This model is a kind of that for molecular dynamics
simulation based on a soft sphere model.
• The model can describe both static friction anddynamical situations which is superior to high speed
simulators such as inelastic hard sphere models.

Outline of DEM
•

•

•
•

DEM includes both contact
repulsive force and dynamical
friction force.
DEM includes the tangential
contact force as well as thenormal contact force.
DEM usually contains therotational motion of particles.
For the simplicity of contact
judgment we usually adopt soft
spheres but we can introduce
more realistic models.

Equations of motion for DEM

DEMの問題点• ただでさえ多いパラメーターの中で決まらないパラメーターが多い。散逸に関するパラメーターは全く未知。
• 形状に依存する。また（後述のように）表面状態によって接触力は全く異なる。
• 大概の粉体は丸くないので接触（摩擦）力は大きく、転がりにくい。球状のモデルでは転がりすぎ。
• 遅い。スケールアップできない。
• シミュレーションできたとしても何も理解が進まない。

動的シミュレーション• 粉体粒子シミュレーション（ＤＥＭ）• 粒径分布 [0.8,1.0]の一様分布• 層高 8, 横幅 90, 奥行き2.25,粒子数 1342
• ばね定数 2 x 10^4, 横ばね定数 1 x 10^4 , 摩
擦係数 0.4• 壁： 摩擦係数が倍、硬さも粒子の倍• 散逸パラメータ：1.0 (対臨界減衰), 時間刻み 2 x
10^{-4}• 角振動数 2.5 振幅 1.6 時刻 1x10^3まで• 粒径 1 mmとすると周波数40 Hz, 振幅 2.0 mm,
加速度振幅 10• 但し最大粒径d, 重力加速度gで無次元化

３次元DEMシミュレーションによるundulationの再現

動的シミュレーション（Ｍｏｖｉｅ）• 粉体粒子シミュ
レーション• 再現される。

Ripple Formation•

振動数が低く容器幅が大きい場合にはrippleになる。（容器の横壁に到達できない）

粉体物理の現状と問題点• 各論としては重要で面白い問題が多い。応用上も重要。
• しかし教養のない物理学者には難しい。
• 結果に普遍性がないために、その結果にあまり興味を持たれない。
• 工学畑に（それも多くの学科に分散して）専門家集団が存在するので後発の物理学者にできることは多くない。
• 些細なモデルの違いで結果が大きく変わる。
– 実験>>シミュレーション>>>>理論• 理論物理の対象としては難しすぎる。
– 理想的かつ簡単な問題を扱うべき。

II. 非平衡物理の現状• Boltzmann

線形応答理論（枠組みの完成）– 入れ物を計算する。
• 臨界現象輸送係数の異常 (Mode Coupling, 動的RG)• Dense Gas の輸送係数の有限性– 新しい方向を探る• Pattern formation (1980-1992)• Chaos (1977-1990)• 複雑系 (1990-2000)– 各論 （粉体も含めて）、五十年一日の如く生物物理

II-A 非平衡の本質： 熱力学の不可逆性

熱力学的自然観と力学的自然観

Ludwig Boltzmann

Boltzmann方程式の導出

H定理（１）

Boltzmann方程式

Ｈ定理（２）平衡分布

問題点

計算機実験によるLoschmidtの反論の実証

Boltzmannの再反論• 実際には統計力学はこの問題を完全に解決していない。
• 平衡統計力学はこの問題を回避してアンサンブルの理論を構成したが、非平衡統計力学は不可逆性の問題に正面から取り組む必要がある。

力学と統計力学は整合するか。

エントロピーを増やすためには粗視化が必要

粗視化とエントロピーについての直感的考察

粗視化についてのイメージ

II-B: Ｌａｎｇｅｖｉｎ方程式とFokker-Planck方程式

Fokker-Planck 方程式の導出

Fokker-Planck方程式とFDT

安定性

加速度効果ランジュバン方程式の導出

一般化Langevin方程式と粗視化

この式は可逆粗視化が必要①初期条件に最も確からしい分布を取る②メモリー効果を切断し、マルコフ近似を取る

III. Boltzmann方程式とその応用• 気体論が非平衡の基本– 粒子の接触がないので最も簡単かつ汎用性がある• Boltzmann方程式の応用を調べることは重要– III-A. Boltzmann方程式と流体力学– III-B. Burnett order solution and its application(金賢得、早川)– III-C. 非弾性ガスの統計力学– III-D. 現実的な粉体ガスとその統計力学（河原田、早川）• 新しいLangevin方程式の応用

III-A. Boltzmann方程式と流体力学

流体力学変数と流体方程式

Pressure tensor and heat current : C-E method

局所平衡流体とBoltzmann方程式

輸送係数と線形応答理論
線形化Boltzmann方程式

輸送係数の続き

但し詳細は省略するが、ハードコア系の計算の結果、熱伝導率と粘性率は以下の通り：

尚、熱伝導率と粘性率の定義は以下で与えられる：III-B.Burnett order solution of Boltzmann
equation (with H-D. Kim)
• The distribution function for hard spheres is expanded as• To solve the Boltzmann equation Burnett assumed

Results (1)Results (2)Comparision between the Burnett solution with MD

Results (3)• This suggests that there are no universal properties as
expected by SST or extended thermodynamics.
• Thermodynamic quantities exhibit qualitative differencesamong hard spheres, Maxwell molecules(=BGK model) andinformation theory.

Thermodynamic quantities• Pressure tensor• Kinetic temperature• Shanon’s entropy

Application to chemical reaction• From our result we can calculate chemical reaction rate
under the heat conduction. Note that the Burnett order isessential for the calculation.
• In Present’s model, the rate for chemical reaction is• Here the cross section is

The result of calculation of reaction rateSummary and Discussion for III-B• We have obtained a systematic perturbative solution of
Boltzmann equation at Burnett order.
• The solution is useful for the test of extended
thermodynamic argument or the calculation of chemical
reaction.
• The effect of boundaries is crucial in this problem. (inprogress)• The correlation effects may not be negligible. We should
check it.

III-C. 非弾性衝突をする粒子系の統計力学• 統計力学の基本は気体論• 非弾性衝突をする気
体はどういう統計則に従うか• 統計的定常状態は保てないので外力で補う必要（空気中のサス
ペンション等）非弾性ガスの特徴• 一様状態が不安定
構造形成• 長距離相関がある• エネルギー散逸によって運動状態が維持できない• 容易に実験できない
Goldhrisch and Zannetti (1993)

非弾性衝突するガス

したがってshear modeはつねに不安定 (河原（九大理））

Haff’s law

Chapman-Enskog method for inelastic gases

続きEvaluation of transport coefficients

輸送係数の計算結果（はねかえり係数依存性）

The results are quantitatively accurate.
The result is inconsistent with theprediction by Green-Kubo formula!
Green-Kubo formula can be used only
elastic case.

The linear stability analysis for uniform state• As suggested the shear mode is always unstable.
• Other three modes are degenerated but one of three is unstable,

不安定

Beyond Haff’s law（河原（九大理）の３次元シミュレーションより）

Many Body Effects in Granular Gases• Correlation effects are always important in granular gases,
because uniform states are unstable.
• Intrinsic long-range correlation plays important roles.
• The velocity correlation function is defined by

BBGKY hierarchy

Long range correlation

Driven Granular Gases by Heat Bath
• 熱浴で励起された非弾性ガス
サスペンションの簡単なモデルであり定常化する• 高速流動層として応用上も重要（非弾性衝突＋流体）

Model and framework• We start from the inelastic Boltzmann equation

where J[f,f] is

g=v-v_1 is the relative velocity andwhere e is the restitution constant.

The Fokker-Planck operator L_{FP} is defined throughwhere the first and the second terns represent friction andthe activation.
Hydrodynamic variables are (V=v-u)

The loss of energy in each collision isThus, the macroscopic loss of the energy appears as

The balance equations for hydrodynamic variables arewhere P_{ij} and q are the pressure tensor and heat flux,
respectively.

Homogeneous States• The basic solution is that for homogeneous state. We
assume

The scaling form of the inelastic Boltzmann equation is given by

where

are the function ofHere we note the distribution function andandcan be obtained as the functions of time and e.

The steady value ofis plotted as a function of e.

Transport coefficients
The transport coefficients such as viscosity and heat
conductivity are introduced as

where
These are obtainable. For example, the viscosity is given by

Is the viscosity of e=1.

Results (1)Results (2)一様状態の線形安定性(in progress)• 弱い熱浴: usual instabilities• 中位の熱浴– shear modeは低波長側で常に不安定(e<1)– 残りの温度、圧縮、密度モードは熱浴の導入で安定化し得る。
渦の発生のみ• 強い熱浴– 安定化• 定量的評価はまだやっていない。
• 相関効果を論じる必然性がある。(in progress)

Application of granular gas kinetic theory for flow
on a slope or in a pipe
• Collisional granular
flow is similar to
granular gases.
• In the vicinity of theboundary, rotation ofparticles is excited,• This flow can becharacterized by
microploar fluid
mechanics.

Validity of micropolar gas-kinetic theory• Figure shows comparison
between MD simulation andanalytic expression of
micropolar fluid mechanics
for intrinsic spin field.
• Agreement is almost perfect.
• Unfortunately, the effect ofspins to translational velocity
is not large.
• See, N.Mitarai, HH and H.
Nakanishi, PRL 88, 174301(2002).

ＩＩＩ－Ｄ．Experimental Accesible GranularGases and its Statistical Mechanics(with A. Kawarada)• Motivation• Purpose• Simulation setup• Simulation Result and Analysis• Discussion• Conclusion

Motivation• Granular gas is
Clustering of granular gases
G.Zanetti,Phys. Rev.Lett.
easiest model of granules. I.Goldhrirsch and70,1619(1993)• There are many researches by theoreticaland numerical methods.
• Most of systems cannot be regarded as
idealistic granular gases.

Purpose• We propose an experimental accessible
system of granular gases.
• We find characteristic behavior of velocity
distribution function (VDF) in such systems.
• We make a theory for the system.
• Parts of the contents in this section has been
accepted for publication in J. Phys. Soc.Jpn.
(cond-mat/0404633)

Setup of our Simulation
width

depth

ShakingheightSide view

Top view

Fixed spherical particles (scatters):3500, 10500, 35000
Free spherical particles:1000, 3000, 10000

Container size :height=1.8d, width=depth=110d , 198.1d , 330d[d:diameter of particle]
Area fraction :0.06

Volume fraction 0.02

Amplitude:1.2d, Frequency:

0 .5 g / d,Acceleration amplitude:11.5

Quasi-two dimension

Simulation Result andAnalysis

Energy
4300

2800

E

E

4300

2800

1300

1300

0

0
0

25
t/T

0

50

10

20

30

40

t/T
(a)

(b)

(c)

(d)

The system reaches a statistical “steady state”
after about 25 cycles of the oscillation.

Granular temperature
0.08

1 N
2
m
(
V
(
t
)
−
V
(
t
))
θd =
∑ i
0
dN i =10.06

θθ1 0.04
3

1 N
V0 (t ) =Vi ( t )∑N i =1

0.02

0

0

10

20

30

40

50

t/T

The velocity fluctuation reaches the statistical“steady state” .

Density
0.24
0.16
0.08
0
33
11
density

x

22
0.15

11
33
0.1

22
y
0.05

Density: at 45 period

The system does not have any definite clustersVortex and Velocity filed
0.3
0.15
0
-0.15
-0.3

30
25
20
15

33

10
5

11

0
0

5

10

15

20

25

30

Velocity filed: at 46 periodvortex
0.2

x

11

22
0.1
0

33

22
y
-0.1
-0.2

Vortex: at 46 period

The system does not have any systematic flows

Velocity Distribution Function (VDF)
various area fraction
ψ=6%
ψ=30%
ψ=50%

1

ψ=7%

~ 0.1
f(c)

0.01

0.001
-3

-2

-1

0

−2

c
~

1

2

f (v, t) = nv0 (t) f (v / v0 ), v / v0 = c

3Velocity Distribution Function(VDF) Dependence on space
Dependence on time
(A)
(B)

1000

(A)
(B)
(C)
(D)

10

100

n

n

1

10

0.1
1
-2

-1.5

-1

-0.5

0

0.5

1

v

Phase:（A）π、（B）3π/2

1.5

2

-2

-1

0

1

2

v

（A) corner、(B)（C）along wall、(D) center。

A little dependence on time and space exists
at high energy tail.

Flatness (FN)
VDF in the cooling process is nearly Gaussian.
FN ≡

c4
c

2

VDF can be approximated by an exponentialfunction.

2

VDF without Coulomb’s friction
undulation
cooling state
no tangential
no friction

1

0.1~ 0.01

f(c)

0.001

0.0001

-4

-3

-2

-1

0

1

2

3

4

c

The essence to produce the large flatness
in VDF is Coulomb’s friction.

Analysis VDF• Langevin equationduu= −γ +ηdtu• Fluctuation-dissipation relationη α (t) = 0, η α (t) η β (t' ) = 2 D δ αβ δ(t − t' )• Fokker-Planck (Kramers equation)γ ∂ u∂ P( x , u , t)
∂2
∂t

={

⋅

m ∂u u+D• Solution (Steady state) ∴ P (u ) = 2
st

∂u

2

} P( x , u , t)m
exp  −
3T


3m 
u
T 

The langevin equation with Coulomb’s friction law
obey the exponential VDF.

Diffusion coefficientmean square displasement

50

f(x)
date

40

< ( r ( t ) − r ( 0 )) 2 >≈ Dt a
a = 1 .094 ± 0 .002403( t ≥ 100 )

30
20
10
0
10

20

30

40

50

t/T

The diffusion coefficient is evaluated from the
2
(
)
r
(t)
−
r
(t
)
= 4 D(t − t 0 )0relation

Fluctuation-Dissipation Relation
0.035

Fitting function

0.03~~ 0.452
D = 0 .0735 T− 0 .007~~
D = 0 .0744 T − 0 .0054

0.025

~
D 0.02

~~
D(T)
(a)
(b)

0.015

D =γ T

0.01

0.005
0.15

0.2

0.25

0.3

~
T

0.35

0.4

0.45

3m

0.5

Particle number: (a)3000 ,(b)1000

The friction γ is independent of T.

⇒ γ = 0 .129 g

Free-Cooling state

VDF in a free-cooling state(× d / g )
6.5

50
250
600

1

6

5.5

0.1

~
f(c)

5

FN
4.5

0.01

4
0.001
-2

3.5

-1.5

-1

-0.5

0

c

0.5

1

1.5

2

0

100

200

300

400

500

t

VDF in cooling process is independent of time

600

(

d/g

)

Energy decay
10000

T(0)
E(t) =+cα(1 + t/t 0 )

data
f(x)
1000

T(0) = 255.0 ± 0.6295
t 0 = 267.0 ± 2.399

E100

α = 1.267 ± 0.004988

10

1
0.05

0.5

5

50

t

c = − 0.4764 ± 0.0307

500

(

d/g )

Energy decay is slower than Haff’s law .

Coulomb friction and Langevin equation

Coulomb’s fluid

VDF

VDF (2)Discussion• 日本宇宙フォーラムでは宇宙環境実験の対象として粉体を重視している。
–９月に微小重力物理の代表例として粉体物理の紹介をする。
–実験の申請も佐野さん（東大）と一緒に行っている。
• ヒアリングが来週ある。

• 宇宙物理との絡み：小惑星物理や惑星のリング形成における非弾性ガス物理の役割

IV．Physics of Inelastic collisions by elastic
materials（Kuninaka & HH）
Inelastic collisions are characterized by the− v
coefficient of restitution defined by e =

r

/ v i

Plastic deformation
Quasi-static region

1 − e ∝ v1/5

非弾性衝突の物理• 巨視的物体の非弾性衝突は散逸や不可逆性の起源を考えるのに適した現象（かもしれない）。
• 粉体物理の最も基礎的なプロセス• 非弾性衝突自体あまりよく研究されているとは言い難い– 正面衝突でのはねかえり係数の衝突速度依存性– ななめ衝突でのはねかえり係数– e>1となる異常な現象

ＩＶ－Ａ．準静的な衝突理論• Hertzの接触理論– 圧縮量hと弾性力F_{el}の関
係

R,R’: 半径： Y,Y_0: Young 率,
σ_0, ポアソン比One point lesson for Hertzian contact theory

無散逸の衝突論と接触時間

準静的な衝突理論（２）• 摩擦力• 運動方程式• 無次元（）

準静的な理論

IV-B 弾性円盤の
弾性円盤の数値シミュレーション
数値シミュレーション
これまでの研究： 一体問題（衝突物体のみで壁の内部自由度なし）
法線方向のはねかえり係数の変化
Hisao Hayakawa and Hiroto Kuninaka"Simulation and theory of the impact of two-dimensional elastic disks"
Chem. Eng. Sci. 57 (2002) 239-252:正面衝突時の正面衝突時の衝突速度と衝突速度と法線方向のはねかえり
法線方向のはねかえり係数のはねかえり係数の係数の解析

等温完全弾性体モデルでは不十分

An example of simulation of inelastic collision
by a completely elastic disk and a wall.

壁と円盤の円盤の衝突モデル
衝突モデル（モデル（無散逸）
無散逸）

2R

質点数：円盤は1600で、一様乱数で初期配置
質点間相互作用：1 2 1
V ( x) = k1 x + k 2 x 4 (k 2 / k1 ≅ 10 −3 )2
4
（Xは自然長からの変位）

ポアソン比： ν ≅ 7 .5 × 10

8R
（円盤の並進エネルギー）＝ （円盤の並進エネルギー）＋（内部振動）
衝突前

衝突後

−2

モデルの説明• ランダム格子モデル：有限サイズで格子配向による異方性を避ける• 壁は固定境界（横と縦）Numerical Modeling of interaction between a wall and a disk

表面のroughness

roughnessなし
分散 δ = 2 .15 × 10 −2 R で正規乱数 dri , dy j を発生させる
平均は０
円盤：

( R, θ i ) ⇒ ( R + dri , θ i )壁：( x j , y j ) ⇒ ( x j , y j + dy j )表面粒子の相互作用のモデル図

Comparision with Hertzian contact theory

直線は2次元Hertz理論。両者が一致するには散逸が必要

Comparision with quasi-static theory

d_1が１になっていない！

接線方向のはねかえり
接線方向のはねかえり係数のはねかえり係数

衝突角度

衝突前

vt'
ψ2
β ≡− =−vt
ψ1−1 ≤ β ≤ 1

衝突後vtψ 1 = − = tan γvnvc ' = v 'i − R ω ' × nv 't
ψ2 = −vn

For β: See Kuninaka and HH, JPSJ 72, 1655 (2003).

Waltonのの衝突理論

O.R.Walton and R.L.Braun:
J. of Rheol., 30 (1986) 949.

slip
stick

µ 0 : 摩擦係数

− β 0ψ 1sticke : はねかえり係数mR 2
ψ 1 − µ 0 (1 + e)(1 +
)
I

slip

mR 2
1+
I

=3 （2D）=7/2 (3D)

cot γ

cot γ VS. β

Outline of derivation of Walton’s theory
•
•
•
•
•
•
•

衝撃
トルクの変化
接点の速度
法線方向の速度変化
接線方向の速度変化
滑っているときの摩擦法則
これらにcotγ＝（n・J)/|n×J|を• 用いるとJ=-m(1+e)(n.v_c)+μm(1+β） cotγ[v_c-(n.v_c) n]を得る。そこからWalton則が得られる。

Comparision between Walton’s theory and an experimentψ 1 は Vn = const.

にして変化させる・実線はWaltonの理論式
・Waltonの理論との比較
から、
stick

slip

µ0 ≅ 0.18
β 0 ≅ 0.56
以下の実験結果に近い値
L.Labous et al.:Phys. Rev. E, 56 (1997) 5717.

・表面のroughnessが重要
である。
slip

stick

Mawらによる衝突理論との比較
N.Maw,J.R.Barber,and J.N.Fawcett(1976),
W.J.Stronge:ImpactMechanics(2000)

zx

ただしz方向にはtκ

stiffness= π / 2tc
β x = 3.02, β z = 1

m:換算質量k r :慣性半径

tf

tc

0

κ / e2= (π / 2ηtc ) β x / β zη = (2 −ν ) / 2(1 −ν )ν：ポアソン比

歪の時間発展をx,z方向について解く

・の範囲でfrequencyがΩ/eになることに注意して、境界条件 vz (tc ) = 0 vz (t f ) = −evz (0) で解くと、
v z + u z (t ) = 0 とし、Z方向の成分

時刻 t1 , t 2 , t3の導入

γ大

X方向の成分

時刻 t1 , t2 , t3 の求め方

３つの衝突角度領域の境界での極限をとれば、上の式から、

tan γ

よってこの衝突角度領域で常にstickした状態になるResults

Maw理論
理論との理論との比較
との比較

eは 0 ≤ ψ ≤ 0.67での平均値を用いた。
1

μはWalton理論との比較から決めた（ µ = µ0 ）

IV-C.Anomalous behavior of e>1• はねかえり係数が１を越える現象• 斜め衝突では別に不思議ではない– 普遍的な現象（例：水切り）• ミクロな散逸とマクロな散逸の関係が問題になる。

Comparison of 2D simulation and 2D experiment
Our result

Calsamiglia et al: J.Appl.Mech. 66 (1999) 146.

1200 particles, v_i=0.1c, with gravity

1.3
1.25
1.2
1.15
e

1.1
1.05
1
0.95
0.9
0.85

0

0.5

1

1.5

2

2.5

3

3.5

tan γ

tan 4.5 くらい

k1壁 = 10 × k1円盤

壁が硬いと表面突起に散乱されe>1となりうる4

4.5

☆e>1ののシミュレーション（
シミュレーション（２）−2

k1 = 10 × k1
壁

円盤．．．壁の方を軟らかくする。

(cf. 実験に対応)
衝突速度 V = 0 . 1c に固定して、 ψ 1 を変化させる。

ただし 実験より遙かに高速衝突
壁の底面は固定境界
重力はなし

シミュレーション結果
円盤と壁の密度比を1にして計算を行った。

ρ d ： diskの密度

ρ w ：壁の密度シミュレーション実験

Impact with intermediate γ (tanγ=2)•

The disk isbounded on asoft wall. It is justlike an impact ofa stone on water
surface.

Impact with large γ ( tan γ=10)• Disk slips on
a wall. No
reflection for
normaldirection.

解析●重要だと思われる要因：
壁の局所的な変形と、それに伴うはねかえり係数の補正●弾性理論に基づいて壁の局所変形を記述し、それによるはねかえり係数の補正を考える。

壁の局所変形D. A. Hills et al.: Mechanics of Elastic Contacts、ガーリン『弾性接触論』

f(x)で表される円盤下半分が半平面に押し付けられた
ときの接触面両端の座標→x=a,bc: 頂点のx座標
diskは2次関数で近似( x − c) 2
f ( x) =
+ δh
2R

c
1
1 − 2ν
−1= 1 − tanlπµ (2 − 2ν )

ν : ポアソン比µ : 摩擦係数

tan α ≡f (b ) − f ( a )lと定義する。αは摩擦によって変形した壁の有効傾きαによるeの補正は、となる。ただしまた座標変換からcはシミュレーションのデータから、c～0.56R
さらに

に注意して

を用いている。ここでは掘り起こしがないときの摩擦係数。

Ψ1 vs. μで計算可能Ψ1 vs. μのメカニズム• 突起と円盤が衝突した際にエネルギーロスがある。
• エネルギーロス/元のエネルギーは突起に衝突する毎に一定->何回衝突するかでロスが分かる。
• ロスから横方向の衝撃が分かり、そこから摩擦係数が• 分かる。

Phenomenology (1)Phenomenology (2)Ψ1 vs. e（跳ね返り係数）

現象論で結果を再現メカニズムの本質は把握できた

Discussion• Static theory of elasticityが有効か？を現在チェック中• 内部散逸の効果： 何故バネビーズモデルでは散逸をとらえきれないのか。当然という気もするが、物足りない。
• 初速の効果。実験に比べて速い衝突である。そのためにはねかえり係数のピークが出る。遅い衝突ではピークがlarge γにずれる。定量的な見積もりはまだ。
• ３次元効果。シミュレーションはできなくはないが大変。
• 当初の目的に比べてOutputが物足りない。

非弾性衝突のまとめ
ランダム格子による衝突モデルで、２次元の斜め衝突
シミュレーションを行った。
ランダム格子による衝突モデルでも、衝突角度と、円盤と壁の硬さの比によっては、はねかえり係数は１を超えることがある。
→ただし実験の様に線形には増加しない→衝突角度が大きいと壁が円盤を表面でキャッチできない
弾性論から壁の局所変形による、衝突角度補正を考えた→はねかえり係数の補正式に代入→１を越えるeの定性的振る舞いは再現３日目の概要• 非弾性衝突の項目は思いの外早く紹介が終わった。
• 非弾性気体も同様• より粉体らしい現象の紹介も必要• 遅い流れと静力学の問題を紹介する。

遅い流れ• 斜面上の流れでは粒子が互いに接触しながらゆっくりと流れる。
• 粒子は流れ方向に拘束されて準１次元的な流れとなっている。
• そこではNavier-Stokes流体とは大きく異なる性質を持つ。

Bagnold’s scaling

Bagnold’s scaling (2)接線応力は単位面積あたりに働く力だからτは速度差に反比例するので

Scaling in surface flow profile

(Angle of repose)

The critical line between no flow region and flow region

Pouliquenの実験の解釈

解釈（２）• 破断しやすい方向は斜面上では斜面に沿ってであり、応力鎖が分裂したときに斜面に沿った方向に分裂すると破壊限界を越えるので斜面が崩落する。したがって主軸が曲がる。
• 摩擦係数の低下は応力主軸が曲がることで理解できる。
• と考えている。

VI. 静力学• 応力鎖の統計力学は重要。粉体らしさが最も顕著に現れる。
• 横壁との摩擦も重要。(Janssen’s law)

応力分布は非自明• 砂山の頂点の直下が圧力の極大とは限らない。
• 履歴依存FPA model and dip

粉体は塑性体か弾性体か？
• 塑性体：力はまっすぐ伝わる。破壊寸前（伝統的な見方。
境界条件をつけられない）：双曲型• 弾性体：安定性あり。（境界条件に依存）：楕円型
x

z
f
θ

f

f1

f2

1

f2

1

θ1

f

θ2

(b)

(a)

応力鎖の散乱ルール

数値モデルから再構成した応力鎖塑性論から弾性論へそしてまた弾性論へ• 応力鎖のBoltzmann方程式（分裂のみ）

l：応力鎖の平均自由行程 f: 応力鎖の強度, n:方向ベクトル

Boltzmann equation for the stress fieldChapman-Enskog method for stress chains• We can obtain a static equation for stress balance.

弾性論の回復
弾性論的ストレステンソルと弾性率を決定可能。Poisson
比は熱力学的安定性を満たさない。
深さhでの垂直応力の応答。
実線は弾性論による計算（弾性率に依存しない）。

塑性論から弾性論、そして更に弾性論
•

•
•
•

左図： 粉体の一点からの刺激がどのように伝播するか（ダブルピークからシングルピーク
へ）
ダブルピーク： 塑性論的（直線的）応答
シングルピーク：弾性論的（拡散的）応答
更に応力鎖の融合を考慮すると深い極限で塑性論的応答が回復する。

Fluctuations in granular flows••

Study of Granular Flow（Solid-Gas Mixtures） １９９１－
Experiments by Matsushita’s group ： At the sametime
– They are interested in power spectra.

– Exp. at initial stage a=1.5? , Peng and Herrmann a=4/3, My
prediction: a=4/3 (1996)– The experiments by Matsushita’s group confirms a=4/3.

Experimental Setup of Granular Flowin a Pipe
HopperLaser Sensor

Laser

Glass Pipe

Out of Air

AD converter

Flow Meter

Computer

Flask

Cock

• From laser sensor, theyobserve the time evolutionof density fluctuation.
• The effect of air is crucial
to produce power laws. So
they control the air flow.
• The paper is published in
PRL (Moriyama et al.
1998).

Results and interpretationx = 10 cm
0

10

P(f)

10-2

x = 40 cm

10-4

10-6

10-10
10-1

~ f - 4/ 3

x = 80 cm

10-8

100

101

102

f (Hz)

103

104• They haveconfirmed that theexponent is a=4/3.
• I reproduce theresult from thesimulation of anonlinear Langevinequation.

Analysis Based on Langevin
Equation
10-2

P(f)

10-4

~ f - 4/ 3

10-6

10-8
100

101

102

f

103

104

• From the relaxation theory of local
fluctuation we can derive a=4/3 forweakly stable region.
• But, theory is unsatisfied.

Application of ASEP to granular flows• ASEP (asymmetric simple exclusion process)is an exactly solvable model but I feel that themodel contains some reality. At least, themodel is regarded as a simple one for traffic
flows and granular flows.



~
1

2

~
3

  

~ - Æ~  LSimulation of ASEP
6

t

t

6

j

(

-

A) low density phase= 0:2;= 1:0

j-C ) maximal current phase= 1:0;
= 1:0

(

β

6

C

t

t

6

• Depending on αandβthere are four phases.
• The exact phase diagram(see below).

A
1/2
j

-

j

coexistence line

: ;

=02

:

=02

-

B ) high density phase= 1:0;
= 0:2

B

(

0

1/2

α

Domain wall theory of ASEP on thecoexistence line
• S. Takesue, HH and T. Mitsudo (PRE68015103R) demonstrate that power spectrum
on the coexistence line can be explained by
the domain wall theory.
• On the line, the density is saturated near theexit, but it is almost zero near the entrance.
So the random walk of the domain wall isessential.

DW theory predicts
simulation
a=3/2 (TASEP)
for the power
3 1/2[2α(1−α)(1−2α)]
spectrum!
ω−3/2

103

2L

I(ω)

10

2

101
1
−1

10

10−5

10−4

10−3 ω

• The upper figure
shows comparison of
DW theory andsimulation of ASEP
for power spectrum.
• The lower figure
shows thecomparison of DW
simulation and−2−1
expression.
1
10analytic
10

From a toy model to a realistic model• ASEP is exactly solvable but it is too simple.
Analysis of DW theory is just like exercise.
• If we introduce multi-lanes for flows, themodel may become realistic. If the rule oflane change is not complicated, we may
keep most of good properties of ASEP.

Multi-lane ASEP by T. Mitsudo and HH(inprogress)• The rule of lane-change: if a vehicle exists infront of the moving vehicle, it prefers changing
the lane.
• The external field can be added.
• The result is similar to granular flow on a slope.
• The velocity obeys an exponential function of
depth.

結論• 非弾性衝突の常識には嘘が多い。研究は黎明期にある。はねかえり係数は１を越える。
• 気体力学は基本的なことが分かっていない.
• 非弾性気体の統計力学は応用上も重要で概念的にも重要• 粉体流の問題ではNewton流体とBagnold流体の遷移がよく分かっていない。実際の流れは混合的である。

結論（２）• 粉体は静力学においても双曲的、楕円的な性
質を併せ持つ。
• その静力学の性質が動力学でも重要。摩擦だ
けでなくundulationをはじめとする様々なパターン形成で本質的• ASEPなどの新しい非平衡統計力学の手法を使って揺らぎの解析が可能。

まとめ• 粉体は非平衡物理のスタンダードかつオーソ
ドックスな対象である• 気体物理の自然な延長に非弾性ガスの統計力学がある• 非弾性衝突の物理はまだ発展途上にある。