粉体気体系の
粉体気体系の分子シミュレーション
分子シミュレーション
ー乱流化と乱流化と非平衡定常状態ー
非平衡定常状態ー

名古屋工業大学 礒部雅晴
isobe@nitech.ac.jp
http://stat.fm.nitech.ac.jp/~isobe

アウトライン1.
イベントドリブン
1. イベントドリブン型分子動力学法
イベントドリブン型分子動力学法
型分子動力学法
イベントドリブン型分子動力学法（（Event-Driven Molecular））
Event-Driven
Molecular Dynamics
Dynamics）
2.
次元粉体気体の次元粉体気体乱流化
2. 22次元粉体気体
次元粉体気体の
の乱流化と
乱流化と
と統計則
統計則
次元粉体気体の乱流化と(Granular(Granular Turbulence
Turbulence inin TwoTwo Dimensions)
Dimensions)3.
重力下粉体気体系
重心運動
についての
しい
3. 重力下粉体気体系の
重力下粉体気体系のの重心運動についての重心運動についてのについての新新しい現
しい現
現
重力下粉体気体系の重心運動についての新についての新しい現象論と象論
揺動散逸関係象論とと揺動散逸関係の
揺動散逸関係のの破破れれ象論と揺動散逸関係の(Granluar(Granluar Gas
Gas under
under Gravity)
Gravity)
4.
今後
4. 今後の
今後の
の展望
展望
今後の

１．イベントドリブン型分子動力学法
イベントドリブン型分子動力学法(Event-Driven Molecular Dynamics)

イベントドリブン型分子動力学法
イベントドリブン
イベントドリブン型分子動力学法
型分子動力学法 (Event(EventDriven MD)
MD)
イベントドリブン型分子動力学法(Event--Driven•• 分子動力学法
分子動力学法ニュートンの
ニュートン
運動方程式
計算機
って
直接解
ニュートンのの運動方程式を運動方程式を
を計算機を
計算機を
を使
使って直接解って直接解直接解くことによ
くことによニュートンの運動方程式を計算機をって直接解くことによ
直接解くことによりり多体系の
多体系
動力学
研究
する
重要
多体系の
の動力学を
動力学を
を研究する
研究する
する重要
重要な
な方法論
方法論
多体系の
動力学を
研究する重要
する重要な
重要な

••
••

イベントドリブン型分子動力学法
イベントドリブン
イベントドリブン型分子動力学法
型分子動力学法
イベントドリブン型分子動力学法(Alder&Wainwright(Alder&Wainwright1957)
1957)
タイムステップドリブン型分子動力学法
タイムステップドリブン
タイムステップドリブン型分子動力学法
型分子動力学法
タイムステップドリブン型分子動力学法(Rahman(Rahman1964,
1964,Verlet
Verlet1967)
1967)•• イベント
イベントとはとはイベントとはとは何か??
イベントとは何
とは何
何か

粒子間衝突
境界
弾性壁
熱壁
など
との
衝突
粒子間衝突,
境界（
弾性壁、
熱壁など
など）
との衝突
衝突、
サブセル
粒子間衝突, 境界（
境界（（弾性壁、
弾性壁、、熱壁など
熱壁など）など））との衝突との衝突、衝突、、サブセル
通過など
通過
通過など
など…
通過など
など…

•• イベントスケジューリング
イベントスケジューリング

効率化の
効率化
効率化の
のエッセンス
エッセンス:
効率化のエッセンス:
粒子間衝突の
粒子間衝突
時間
効率
よく
粒子間衝突の
の中
中で
で最
最もも短
短い
い時間を
時間を
を効率よく
効率よく
よく探
探すこと
すこと
粒子間衝突の
時間を
効率よく探
よく探

礒部雅晴,
礒部雅晴 剛体球系分子動力学シミュレーション
剛体球系分子動力学シミュレーション －方法論の
方法論の歴史的発展と歴史的発展と大規模シミュレーショ
大規模シミュレーショ
ンへの応用
への応用－
応用－, 日本物理学会誌 Vol.10 pp.748-753 (2007).

Historical
Historical Review
Review ofof Event
EventEvent--Driven
Driven MD
MD
Alder and Wainwright, J.Chem.Phys. 31, 459 (1959). Pioneer Paper of MD
Studies in Molecular Dynamics. I. General Method

Generalized algorithm of Event-Driven MD including Square-Well potential.

イベントドリブンＭＤ
イベントドリブン
ＭＤ
単純
イベントドリブンＭＤ
ＭＤの
の最
最も
も単純な
単純な
なアルゴリズム
アルゴリズムイベントドリブンＭＤのＭＤの
単純な
剛体球ポテンシャル
剛体球ポテンシャル

衝突則

粒子位置と
粒子位置と速度からすべての速度からすべての粒子からすべての粒子ペア
粒子ペア (i , j)について衝突について衝突する
計算する.
衝突する時間する時間を時間を計算するする O(N2)

最短の最短の時間で時間で衝突する
衝突する粒子する粒子ペア粒子ペアをペアを探す tmin. O(N)

すべての粒子すべての粒子を粒子を上で求めた最短時間
めた最短時間だけ最短時間だけ進だけ進める tmin. O(N)
衝突則を衝突則を使って、って、衝突後の衝突後の粒子の
粒子の速度
を計算する
計算する.
する O(1)

衝突に衝突に関与した
関与した粒子した粒子ペア
粒子ペアにペアに関して、して、次の衝突の衝突の時間を再計算する
再計算する.
する O(N)

イベントドリブンＭＤ
イベントドリブン
ＭＤ
アルゴリズム発展
イベントドリブンＭＤ
ＭＤののアルゴリズム発展
アルゴリズム発展
発展のの歴史
歴史
イベントドリブンＭＤのＭＤのアルゴリズム発展の
発展の
D.C.Rapaport, J. Comp. Phys. 34, 184 (1980).

最小衝突時間の最小衝突時間の計算量を
計算量を劇的に劇的に減らした:らした
O(N) →O(log N)

The Art of Molecular Dynamics Simulation - 2nd edn
D. C. Rapaport (Cambridge University Press, 2004)

２分探索木を
分探索木を使用：
使用：イベントドリブン分子動力学法
イベントドリブン分子動力学法の分子動力学法の本質的な
本質的なブレークスルー

イベントドリブンＭＤ
イベントドリブン
ＭＤ
アルゴリズム発展
イベントドリブンＭＤ
ＭＤののアルゴリズム発展
アルゴリズム発展
発展のの歴史
歴史
イベントドリブンＭＤのＭＤのアルゴリズム発展の発展の
しかし、
しかし、Rapaport (1980) のアルゴリズムはアルゴリズムはとても複雑
とても複雑すぎて
複雑すぎて理解
すぎて理解しづらく
理解しづらく、
しづらく、プロ
グラミングコード作成
年代にの論
グラミングコード作成には
作成には大変には大変な
大変な困難を
困難を要する。
する。90年代
年代に入り、Rapaportの文を基礎に基礎にデータ構造
データ構造を
構造を単純化させ
単純化させ簡略化
させ簡略化させるいくつかの簡略化させるいくつかの改良させるいくつかの改良が改良が行われた。
われた。

B.D.Lubachevsky, J.Comp.Phys. 94,255 (1991).

Heap Algorithm

M. Marin, D. Risso and P. Cordero, J.Comp.Phys. 109, 306 (1993)
M. Marin and P. Cordero, Comp.Phys.Commun. 92, 214 (1995).

More efficient and simple
algorithm by using Complete
Binary Tree.

イベントドリブンＭＤ
イベントドリブン
ＭＤ
アルゴリズム発展
イベントドリブンＭＤ
ＭＤののアルゴリズム発展
アルゴリズム発展
発展のの歴史
歴史
イベントドリブンＭＤのＭＤのアルゴリズム発展の
発展の
M.Isobe, Int. J. Mod. Phys.C 10, 1281 (1999).
•Extended Exclusive Grid Particle
Method•Dynamical Upper Time Cut-off•Cell-Crossing Eventをを考慮しない考慮しない。
しない。
非平衡系に
非平衡系に適用可能な
適用可能な単純で単純で極めて効率的な
効率的な アルゴリズムがアルゴリズムが完成した
完成した。
した。

S. Muller and S. Luding,
J.Comput.Phys. 193, 306 (2004)

並列化への
並列化への試への試み

A. Donev, S. Torquato and F.H.Stillinger,
J. Comput. Phys. 202, 737 (2005)

液晶系への
液晶系への応用
への応用

Development
Development ofof Efficiency
Efficiency ofof Event
EventDriven MDMD
Event--Driven

50 years ago… Alder & Wainwright could perform only 2000
collisions per one (CPU) hour with 100 particles.

Now, we can perform a few billion (~109)collisions per one
hour on the personal computer with one million particles.

EventくのEventDriven MD
MDのの多
多くの応用例
くの応用例
応用例
Event--Driven
MDの
くの応用例
•• イベントドリブン
イベントドリブン型分子動力学
は型分子動力学
シミュレーション
イベントドリブン型分子動力学
型分子動力学シミュレーション
シミュレーションは
は、、
イベントドリブン型分子動力学シミュレーション型分子動力学シミュレーションはシミュレーション
非平衡統計物理学のひとつの非平衡統計物理学のひとつの重要
非平衡統計物理学のひとつののひとつの重要
重要な
なツール
ツール
非平衡統計物理学のひとつの重要のひとつの重要な
重要な

••
••
••
••
••
••

２２次元融解アルダー
次元融解(アルダー転移(アルダー
アルダー転移転移)
アルダー転移転移)２２成分ガラス
成分
ガラス
成分ガラス
ガラスの
のシミュレーション
シミュレーション
成分ガラスの
ガラスの
流体力学の
の
諸問題
流体力学
レイリーベナール
対流
カルマン
流体力学の
(レイリーベナール
レイリーベナール対流
対流，
カルマン渦
渦, ,レイリー
レイリー
流体力学の諸問題(レイリーベナール対流
レイリーベナール対流，
対流，，カルマン渦
カルマン渦
テイラー不安定性
テイラー
テイラー不安定性
不安定性,
など..)
テイラー不安定性
不安定性, など
など..)
化学反応系
化学反応系
粉体気体
粉体
シア
粉体気体(GranularGas),粉体シア
粉体シア
シア流
流, ,粉体振動層
粉体振動層
粉体気体(Granular Gas),
粉体シア流
シア流
微小重力下での
微小重力下
での
コロイド
微小重力下でのでのコロイド
コロイドの
の結晶化
結晶化
微小重力下でのコロイドでのコロイドの
コロイドの
衝撃波
衝撃波
低次元熱伝導問題と低次元熱伝導問題
熱伝導率
低次元熱伝導問題と低次元熱伝導問題とと熱伝導率の
熱伝導率の
熱伝導率の
の発散
発散
ロングタイムテールとロングタイムテール
輸送現象
輸送係数
発散
ロングタイムテールと
輸送現象（
輸送係数の
の発散）
発散）
ロングタイムテールとと輸送現象（
輸送現象（（輸送係数の
輸送係数の
発散））

1.1.
2.2.
3.3.
4.4.
5.5.

微小重力場中の微小重力場中粉体系
分岐
相転移
微小重力場中の
の粉体系の
粉体系の
の分岐（
分岐（
相転移）
(PRE2001他
他))
微小重力場中の
粉体系の分岐（（相転移）相転移））(PRE2001
２２次元粉体気体の次元粉体気体次元粉体気体のの乱流化
乱流化(PRE2003(PRE2003他
他))
次元粉体気体の２２次元ロングタイムテール
次元
ロングタイムテール次元ロングタイムテール
ロングタイムテール問題
問題(PRE(PRE2008
2008他
他))
次元ロングタイムテール問題
ロングタイムテール問題
シアストレス自己相関関数
シアストレス
自己相関関数
シアストレス自己相関関数
自己相関関数とモラセステール(Mol.Sim.2009(Mol.Sim.2009他
他))
シアストレス自己相関関数と自己相関関数ととモラセステール
粉体系の粉体系
非平衡定常状態
しい
粉体系のの非平衡定常状態と非平衡定常状態と新しい現象論
しい現象論
現象論(JPSJ2008(JPSJ2008他
他））
粉体系の非平衡定常状態とと新
しい現象論

••
••
••

２．２次元粉体気体系の次元粉体気体系の乱流化と乱流化と統計則(Granular Turbulence in Two Dimensions)

名古屋工業大学 礒部雅晴
isobe@nitech.ac.jp
http://stat.fm.nitech.ac.jp/~isobe粉体気体
粉体気体(Granular(Granular Gas)
Gas) 再訪
再訪
粉体気体(Granular

1993年の Goldhirsh&Zanetti(PRL.70,1619)の記念碑的論文以来、記念碑的論文以来、粉体自身の
粉体自身の運動の時間スケール
時間スケールに比べ外場（
外場（重力、
重力、振動板など
振動板など）など）で特徴づけられる
特徴づけられる運動
づけられる運動の運動の時間スケー
時間スケー
ルがはるかに小さい系（すなわち微小重力、
微小重力、高励起）に対し、非平衡統計物理学の伝統的手法を拡張した系統的研究が展開され、一大分野として確立されている。
この分野は特に、粉体気体(Granular Gas)と呼ばれ、何冊かの本や系統的レビューも出版されている。



T. Pöschel & S. Luding (Eds.), Granular Gases (Lecture Notein Physics, 564), Springer (2001).; T. Pöschel & N. Brilliantov(Eds.), Granular Gas Dynamics (Lecture Notes in Physics,624), Springer (2003).; N. Brilliantov & T. Pöschel, Kinetic
Theory of Granular Gases, Oxford University Press (2004).; T.
Pöschel & T. Schwager, Computational Granular Dynamics,
Springer (2004).;I. Goldhirsch, Raipd Granular Flows,
Annu.Rev.Fluid.Mech.35, 267 (2003).



粉体気体(Granular Gas)を2次元乱流(Two-Dimensional
Turbulence)の観点から研究を遂行する。

極限系においては、系統的な研究はほとんどない。

自由冷却粉体気体とは
自由冷却粉体気体とは
自由冷却粉体気体とは
とは何
何か
か？
？
自由冷却粉体気体とは何
とは何
外場がない
初期状態が熱平衡状態
非弾性剛体球系

システムサイズ（粒子数）

熱力学極限

密度（占有率）

非圧縮極限

反発係数

準弾性極限

粒子間の非弾性衝突によりいくつかの時空
時空スケール
時空スケールがスケールが出現
時間発展

空間パターン

３つの階層による主な解析手段

一様冷却状態(Homogeneous Cooling State)
シアリング不安定性 (Shearing Instability)

マクロ（流体力学的現象論）
マクロ（流体力学的現象論）

クラスター不安定性(Clustering Instability)

メソ（気体運動論）
メソ（気体運動論）

最終アトラクター状態(Final Attractor State)
非弾性コラプス(Inelastic Collapse)

ミクロ（分子動力学）
ミクロ（分子動力学）

密度場
密度場

速度場
速度場

Calculation
Calculation(2010.6)(2010.6)

N=4.19N=4.19million
millionνν=0.60=0.60
r=0.996368
r=0.996368
C/N=6956
C/N=6956

温度場
温度場

圧縮場
圧縮場

渦度場
渦度場

時間発展
時間発展

HCS

Shearing

Clustering

Final Attractor
2
N=2048
N=20482νν=0.60=0.60
r=0.996r=0.996

密度場
密度場

温度場、渦度場のパワースペクトルがエネルギースペ
クトル、渦度スペク
トル

速度場
速度場

温度場
温度場

E(k)

渦度場
渦度場

Z(k)

Initial State

Onset of 3rd Stage

Final Steady State２２次元乱流現象を
次元乱流現象
こすひとつの
単純次元乱流現象をを起起こすひとつの単純
こすひとつの単純
単純な
なモデル
モデル
次元乱流現象をこすひとつの単純な
単純なIsobe, PRE 68, 040301 (2003).
2003年に行われたEvent-Driven MDを用いた大規模計算では、 自由冷却粉体気体の２次元非弾性剛体球（ＩＨＳ）
次元非弾性剛体球（ＩＨＳ）系（ＩＨＳ）系に決定論的熱浴（
決定論的熱浴（速度スケーリング
速度スケーリング）
スケーリング）を課したモデル系において、２次元乱流に特有ないくつもの現象が発見され、２次元流体乱流で
次元流体乱流で予測されたスケーリング
されたスケーリング指数
スケーリング指数そのものが指数そのものが（そのものが（初めて）めて）得られた。
２次元非弾性剛体球（
次元非弾性剛体球（剛体円板）
剛体円板）モデルはとても
モデルはとても単純はとても単純な単純なモデルである
デルである→
である→
初期条件 & 境界条件
熱平衡状態（密度は一様、速度はマックスウェル分布）
２次元周期境界条件、無重力状態

発達した
発達した乱流した乱流を
乱流を研究するため研究するため、するため、メソ、
メソ、マクロな
マクロな理論解析が
理論解析が
不可能な
不可能な領域に領域に着目し着目し、初期から最終状態に至るまでの長時間かつ広範なパラメーター空間における乱流化とエネ
空間構造
ルギー＆渦度スペクトル、散逸スケールなどの空間構造の形成や形成や相関関数の相関関数の時間発展、
時間発展、スケーリングのスケーリングの性質に関性質
して系統的に調べることを目標。

M.Isobe, Phys. Rev. E 68, 040301 (2003).

自由冷却粉体気体の
自由冷却粉体気体
相図
反発係数
密度空間
自由冷却粉体気体の
の相図（
相図（
反発係数－
－密度空間）
密度空間）
自由冷却粉体気体の相図（（反発係数－
反発係数－密度空間））
Point Particles

密度（占有率）Close-packed Density
Alder (KTHNY) Transition

熱平衡
(r=1)
HCS
マクロ（流体力学的現象論）
マクロ（流体力学的現象論）
運動論－シア
運動論－シア領域境界
領域境界

シア－コラプス
シア－コラプス
領域境界領域境界

メソ（気体運動論）
メソ（気体運動論）
HCS
Shearing
Clustering
Final Attractor

ミクロ（分子動力学）
ミクロ（分子動力学）
ＵｎｋｏｗｎHCS
Shearing
Clustering
Inelastic Collapse
システムサイズ依存

反発係数
ｒ→０

系統的な理論解析が行われている範囲
は主に低密度、弾性領域のごくわずかである。

低密度系(ν<0.55?)

自由冷却粉体気体の
自由冷却粉体気体
時間発展
自由冷却粉体気体のの時間発展の
時間発展の
のシナリオ
シナリオ
自由冷却粉体気体の
時間発展の

時間発展

平衡状態
(r=1)
運動論－シア
運動論－シア
領域境界
領域境界

HCS
マクロ（流体力学的現象論）
マクロ（流体力学的現象論）
メソ（気体運動論）
メソ（気体運動論）
HCS
Shearing

ミクロ（分子動力学）
ミクロ（分子動力学）
Clustering

シア－コラプス
シア－コラプス
領域境界
領域境界Final Attractor
ＵｎｋｏｗｎHCS

Shearing

Clustering

Inelastic Collapse

システムサイズ依存

反発係数
ｒ→０

系統的な理論解析が行われている範囲
は時間発展の初期段階のみである。

高密度系（排除体積効果）はよくわからない。

高密度系での高密度系
での最終高密度系でのでの最終最終アトラクター
アトラクター
高密度系での最終での最終アトラクター
最終アトラクター
McNamara
McNamara&&Young,
Young,PRE
PRE53,
53,50895089(1996).
(1996).
高密系にてモードの存在がEDMDシミュレーションの結果から示唆
高密系にて22Vortex
Vortexモードの存在がEDMDシミュレーションの結果から示唆
された。
された。
Soto,
Soto,Mareschal
Mareschal&&Mansour,
Mansour,PRE
PRE62,
62,38363836(2000).
(2000).
振幅方程式を用いた現象論：シアバンドと2
Vortexの間のモードの違いなど、
振幅方程式を用いた現象論：シアバンドと2 Vortexの間のモードの違いなど、研究された。
研究された。
Hill
Hill&&Mazenko,
Mazenko,PRE
PRE63,
63,031303
031303(2001);(2001);PRE
PRE67,061302
67,061302(2003).
(2003).
高密系でも扱えるように修正された現象論的流体方程式が提案され、DNSに高密系でも扱えるように修正された現象論的流体方程式が提案され、DNSにより最終アトラクターが状態が研究された。
より最終アトラクターが状態が研究された。

2
N=512
N=5122νν=0.75=0.75

Temperature
TemperatureField
Field(↓)EDMDで新しく見つかった1(↓)EDMDで新しく見つかった1Vortex
Vortexアトラクター状態
アトラクター状態

11Vortex
Vortexr=0.9968
r=0.9968

Shear
ShearBandBand
r=0.9958
r=0.9958

22Vortex
Vortexr=0.9953
r=0.9953

Bose
BoseBose--Condensed
Condensed State(?)
State(?) in
in Granular
Granular Gas
Gas




1 Vortex Mode: Only particles near the edge of the void-crystal boundaries are
moving actively, therefore, it looks like a ring shape in the temperature field. The
prediction of final attractors seems to be difficult since different final attractors(1Vortex, Shear, 2 Vortex) appear even if the restitution coefficient is almost
the same. Of course, these results are out of expectation of the previous
phenomenological theory.
Bose-Einstein condensation: In the dense system, other physical mechanisms
might be relevant, such as the condensed state in terms of 2D turbulence. Indeed, inthe study of 2D turbulence, when the minimal correlation scale (such as
Kolmogorov or dissipation scale) is larger than the system size, the one or two
vortex state appears and a type of Bose-Einstein condensation was pointed out by
Kraichnan (Tabeling, Phys.Rep. 362 (2002)).
Some
SomeExamples
ExamplesofofBose-Condensed
Bose-CondensedState
StateininFluid
FluidTurbulence
Turbulence

Experiment
Experiment

DNS
DNS

Energy
Energy Spectrum
Spectrum on
on 2D
2D Turbulence
Turbulence
TheoryKraichnan-Batchelor Theory

Energy Spectrum E(k) [m3s-2]
Enstrophy Transfer Rateη[s-3]Wave number k [m-1]

••
••

••

••
••

R.H.Kraichnan,
R.H.Kraichnan,Phys.Fluids
Phys.Fluids10,
10,14171417(1967).
(1967).

The
Theinverse
inverseenergy
energycascadecascadeand
andforward
forwardenstrophy
enstrophycascadecascade(2D(2Dforced
forcedturbulence)
turbulence)

G.Batchelor,
G.Batchelor,Phys.Fluids
Phys.FluidsSuppl.
Suppl.IIII12,
12,233
233(1969).
(1969).

The
Theforward
forwardenstrophy
enstrophycascadecascade(2D(2Dfreely
freelydecaying
decayingturbulence)
turbulence)From
thedimensional
analysis(likeaKolmogorov
scaling
From the dimensional analysis (like a Kolmogorov scalingin
in3D
3Dturbulence),
turbulence), ––33power
powerlaw
lawin
inenergy
energyspectrum
spectrumisisderived
derivedunder
underthetheassumptions
assumptionsthat
thatthetheenergyenergyspectrum
spectrumonly
onlydepends
dependsononenstrophy
enstrophytransfer
transferandandwavewavenumber.
number.

エネルギースペクトルのエネルギースペクトル指数
スケーリング
エネルギースペクトルのの指数と指数ととスケーリング則
スケーリング則
則
エネルギースペクトルの指数とスケーリング則••2nd
2nd Stage~3に移行する際、エネStage~3rdrdStage
Stageに移行する際、エネ
ルギースペクトルの低波数側は、
2次元乱ルギースペクトルの低波数側は、2次元乱流におけるKraichnan-Leith流におけるKraichnan-LeithBatchelor(KLB)
理論の予測と同じ指数-3をBatchelor(KLB)理論の予測と同じ指数-3を
とる。
とる。
••K
Kdから最小散逸スケールが見積もれる。
dから最小散逸スケールが見積もれる。
ldl/ /σσ~36.9=0.60)~36.9(ν(ν=0.60)d••散逸領域（2次元）には、ほぼ散逸領域（2次元）には、ほぼ3粒子が関与している。
N=4
ν/ (l(ld/σ)2/π∼102/π∼103粒子が関与している。
N=4ν/
/σ)
d
流体力学の議論をするには、最小散逸ス
流体力学の議論をするには、最小散逸ス
ケールよりも大きなシステムサイズが必要。
ケールよりも大きなシステムサイズが必要。

・・最小散逸スケールはσ/(1-r2)20.5(σ:粒子直径)
最小散逸スケールはσ/(1-r)0.5でスケーリングできる。
でスケーリングできる。(σ:粒子直径)・（スケールされた）最小散逸スケールは、粒子数、反発係数の依存・（スケールされた）最小散逸スケールは、粒子数、反発係数の依存性はほとんどなく、粒子占有率に依存。
ν=0.60では、ほぼ2の値をとる。
性はほとんどなく、粒子占有率に依存。ν=0.60では、ほぼ2の値をとる。
(ν(ν=0.25では、=0.25では、 KKd ~0.7
、 ν=0.70ではKd ~3).粒子占有率でのスケーリン
d ~0.7、 ν=0.70ではKd ~3).粒子占有率でのスケーリン
グは現在、確認中。
グは現在、確認中。
・・慣性領域
慣性領域は
慣性領域
反発係数
準弾性極限
レイノルズ
慣性領域は
は、、反発係数が反発係数がが準弾性極限で準弾性極限でで大大。。→→高高レイノルズ数
レイノルズ数
数。。
慣性領域は
反発係数が
準弾性極限でレイノルズ数自己組織化する
自己組織化
する
乱流
パターン
最小散逸
自己組織化する
する乱流
乱流パターンパターンとと最小散逸スケール
最小散逸スケール
スケール
自己組織化する乱流
する乱流パターン
乱流パターンとパターンと最小散逸スケール
N=4.19N=4.19million
millionνν=0.60=0.60
r=0.996368
r=0.996368
C/N=6956
C/N=6956

温度場
温度場

McNamara
McNamara&&Young
Young(1996)(1996)
Peng
Peng&&Ohta
Ohta(1997)(1997)
Radjai&Radjai&Roux(2003)
Roux(2003)
過去のこれらの研究は、数過去のこれらの研究は、数千粒子の規模で研究され千粒子の規模で研究されている。しかし、そのサイズている。しかし、そのサイズはほとんど空間相関の最はほとんど空間相関の最小散逸スケールと同じであ小散逸スケールと同じであ
る（！）
る（！）

N~103
最小散逸スケール

エネルギー緩和
エネルギー
緩和
クラスター
による
逸脱
エネルギー緩和
緩和（
Haff則
則と
とクラスター化
クラスター化化による逸脱による逸脱
逸脱）エネルギー緩和（
緩和（（Haff
Haff則クラスター化
による逸脱）
逸脱））

-1
•3
•3rd
rdStageは、E(t)~tStageは、E(t)~t-1が理論予測されていが理論予測されているが、高密度では逸脱する。さらに、tるが、高密度では逸脱する。さらに、ts~800
s~800付近から、べき緩和からも逸脱する。
付近から、べき緩和からも逸脱する。

べき緩和
緩和から
べき
緩和
から
の逸脱

・・シミュレーションは、全運動エネルギー（温シミュレーションは、全運動エネルギー（温度）を一定の拘束条件で、時間をスケール度）を一定の拘束条件で、時間をスケールしている。tとtStageではHaffしている。tとtsの関係は、1~2sの関係は、1~2 StageではHaff
則の結果と一致し、クラスタ化後(3stage~)則の結果と一致し、クラスタ化後(3stage~)2は、t~t
は、t~ts s2の関係になる。
の関係になる。

衝突率の
衝突率
時間発展
クラスター
衝突率の
の時間発展、
時間発展、
クラスター衝突領域
衝突領域
衝突率の時間発展、、クラスター衝突領域
クラスター衝突領域
圧縮場
圧縮場

N=4.19N=4.19million
millionνν=0.60=0.60
r=0.996368
r=0.996368
C/N=6956
C/N=6956

相対拡散リチャードソン
相対拡散ととリチャードソン則
リチャードソン則
則
相対拡散とリチャードソン則

TheoryJ.P.L.C.
5-32.
J.P.L.C.Salazar&L.R.Collins,
Salazar&L.R.Collins,Annu.Rev.Fluid
Annu.Rev.FluidMech.
Mech.41,
41,(2009)(2009)40
405-32.

Richardson Law→ R(t)~t3 (L.F.Richardson, Proc.R.Soc.London Sect.A 110, 709 (1926).)

1.
2.
3.

Dissipation Subrange:
Inertial Subrange:
Diffusion Subrange:

Inertial Subrange (2D)

•第第3ステージ
ステージは
ステージは、r(ts) >ld となったと同時となったと同時に同時に始まる。
まる。
ldがDissipation SubrangeとIntertial Subrangeを分ける(Kolmogorovスケールのような)よい指標になっていることがわかる。
•システムサイズLが十分に大きいため、第
第4ス
ス
テージ（
テージ（クラスター衝突
クラスター衝突）衝突）は、粒子の粒子の相対拡散r(t相対拡散 s)
がLに
に到達する
到達する前
する前で移行すること移行することがわかる。つますることり、第4ステージは周期境界条件の影響ではなく、熱力学極限で存在すると考えられる。

テイラーマイクロスケールとテイラーマイクロスケール
レイノルズ
テイラーマイクロスケールとテイラーマイクロスケールととレイノルズ数
レイノルズ数
レイノルズ数
数• テイラーマイクロスケール: λ流体乱流では、テイラーマイクロスケールがよく使用されており、 2次元非圧縮性等方性乱流では、
λ=(8ν<uλ=(8ν< 2>/<ε>)1/2となる。
となる。

エンストロフィー散
エンストロフィーエンストロフィー散散
エンストロフィー散
逸率の
の
逸率
逸率の
逸率
の時間発展
時間発展

• レイノルズ
レイノルズ数
数: R λ=<u2>1/2 λ/νテイラーマイクロスケールを使ってレイノルズ数を導出すると、実験室規模では、O(102~103)、地球物理規模ではO(104)で乱流となることが知られている。
粉体気体のエネルギー散逸率<ε>は、時間と共に減少しシミュレーションから見積もれる。 動粘性係数νはr=1のEnskog理論から求まる。 ts~0ではR λ∼12.4であ∼12.4
るが、エンストロフィー散逸率が最大となるクラスタリン
グ領域直前ts~252で R λ∼445.4 に増大し、乱流化する。

粉体気体系の粉体気体系の準弾性極限における準弾性極限におけるレイノルズにおけるレイノルズ数
レイノルズ数は初期平衡状態から
初期平衡状態から時間
から時間と時間と共に増え、非圧縮一様等方性の圧縮一様等方性の条件が条件が適用可能と適用可能と思われる第われる第2ステージのステージの終わりにおいては、
わりにおいては、発達し発達した乱流の乱流の値まで増加まで増加する増加する。
する。粉体気体系の粉体気体系の乱流化を乱流化を示す定量的な定量的な根拠が根拠が得られた。
られた。

自由冷却粉体気体の
自由冷却粉体気体
相図
反発係数
密度空間
自由冷却粉体気体の
の相図（
相図（
反発係数－
－密度空間）
密度空間）
自由冷却粉体気体の相図（（反発係数－
反発係数－密度空間））

密度（占有率）Close-packed Density

Point Particles

Alder (KTHNY) Transition

熱平衡
(r=1)
HCS
運動論－シア
運動論－シア領域境界
領域境界

シア－コラプス
シア－コラプス
領域境界
領域境界

HCS
Shearing
Clustering
エンストロフィーカスケード
Final Attractor

Bose-Einstein凝縮

HCS
Shearing
Clustering
Inelastic Collapse
システムサイズ依存
ミクロ（分子動力学）
ミクロ（分子動力学）

反発係数
ｒ→０

システムサイズ依存？

EDMDによる粉体気体系
による粉体気体系は粉体気体系は乱流の乱流の統計則並びに統計則並びに乱流化
びに乱流化の乱流化のミクロな起源を起源を探究する
探究する重要する重要な
重要な方法論を
方法論を提供していると提供していると思していると思われる。
われる。

３．重力下粉体気体系の
重力下粉体気体系の重心運動についての重心運動についての新しい現象論
しい現象論と現象論と揺動散逸関係の
揺動散逸関係の破れ(Granular Gas under Gravity)

Miyakonojo National College of Technology
Jun’ichi Wakou
Nagoya Institute of Technology Masaharu IsobeDrivenDriven Granular
Granular Gas
Gas under
under Gravity
Gravity

Introduction•• Energy
Energy Balance
Balance•• SupplySupply(External(External Driving):
Driving):VibratingVibratingWall,
Wall,Thermal
ThermalWall
Wall…
…

•• Dissipation
Dissipation(Inelastic(InelasticCollisions):
Collisions):••

DrivenDrivenGranular
GranularGas
Gasis
isone
oneof
ofprototype
prototypeto
toinvestigate
investigate
thetheNonequilibrium
NonequilibriumSteady
SteadyState
State (NESS)(NESS)

Vibrating Wall

ThermalWall

Collective
Collective Behavior
Behavior in
in NESSNESS

Introduction

Various
Various Interesting
Interesting Phenomena
Phenomena
••
••
••
••

2
Convection
Convection(Faraday(Faraday1831,
1831,Evesque,
Evesque,Rajchenbach
Rajchenbach89)
89)b
2Waves(Fauve
Surface
SurfaceWaves(Fauveetetal.
al.89,
89,Umbanhowar
Umbanhowaret
etal.
al.96)
96)
Size
SizeSegregation
Segregation(Duran(Duranet
etal.
al.93)
93)
Transition
Transitionfrom
fromaacondensed
condensedstate
statetotoaafluidized
fluidizedstate
state(Luding(Ludingetetal.
al.94,
94,Goldshtein
Goldshteinet
etal.
al.95).
95).

dZ
M= −Mg + F
dt

Macroscopic
Macroscopic Characterization
Characterization••••

ItItisisimportant
importantto
tostudy
studyhowhowmacroscopic
macroscopicquantities
quantitiesin
inthe
the
systemsystem(Center(Centerof
ofmass,mass,Granular
Granulartemperature,temperature,etc..)
etc..)
dependdependononthe
thesystem
systemparametersparameters(N,(N,r,r,A,
A,ω...
ω...)ω...)Scaling
Scalingrelationshiprelationshipfor
forthethemacroscopic
macroscopicquantities
quantities(Luding(Ludingetetal.
al.94,
94,Warr
Warretetal.
al.95).
95).

Purpose
Purpose ofof thisthis Study
Study (Our
(Our Goals)
Goals)

Introduction

 We
Wepropose
proposeaanovel
novelapproach
approachbased
basedononaaLangevinLangevinequation
equationformotion2
forfluctuating
fluctuating
motionof
ofthe
thecenter
centerofofmassmassofofdriven
drivendZgranulargas
gravity.
granularundergravity.
M gasunder
= −Mg
+F

dt

2

b

 In
Inthisthisframework,
framework,thetheanalyticalanalyticalsolution
solutionof
ofthe
theLangevinLangevinequation
equationisisused
usedtotoderive
deriveanalytic
analyticexpression
expressionfor
forseveral
severalmacroscopic
macroscopicquantities
quantitiesandandthethepower
powerspectrum
spectrumfor
forthethecentercenterofofmass(COM).
mass(COM).
 In
Inorderordertototesttestourourtheory,
theory,weweperformed
performedextensive
extensiveeventeventdriven
drivenmolecularmoleculardynamics
dynamicssimulations
simulationsfor
foroneone-andandtwotwodimensional
dimensionalsystems
systemsforforaawide
widerange
rangeofofparametersparametersand
andcomparecomparewith
withthethetheoreticaltheoreticalpredictions.
predictions.

MotionMotion ofof thethe Center
Center ofof Mass
Mass

Theory•• Equation
Equation ofof motion
motion ofof Center
Center ofof Mass
Mass (COM)
(COM)
2
2 d Z
dM
Z 2 = −Mg + Fb
M 2 dt= −Mg + Fb

dt ofparticles
M
M::total
totalmassmassofparticles
Z(t)Z(t)::heightheight ofofthe
theCOMCOM
FFbb(t)(t)::force
forceexerted
exertedby
bythe
thebottom
bottom plate
plate•• Important
Importanttime
timescales:
scales:τ：τ：oscillation
oscillationperiod
periodof
ofthethebottom
bottomplate,plate,
ττrelrel::macroscopic
macroscopicrelaxation
relaxationtime.
time.
Fluidizedττ≪τ
Fluidizedstate
state≪τrelrelg

Fb

J. Wakou, A. Ochiai and M. Isobe, J. Phys. Soc. Jpn. 77, 034402 (2008)

Theoretical
Theoretical Approach
Approach(Course(Course(Course--Graining
Grainingfrom
fromMicro
MicroMicro--totoMacroscopic)
Macroscopic)

Introduction•• Molecular
Molecular Dynamics
Dynamics

Microscopic(GallasLuding et al.94) etc..
(Gallasetetal.92,
al.92,
2 Luding et al.94) etc..

dZ
M 2 = −Mg + Fb•• Kinetic
ThoeryKineticdt
Thoery (Jenkins(JenkinsSavage
Savage83)
83)
Velocity
VelocityDistribution
Distribution(Bernu(Bernuetetal.
al.94)
94)
Scaling
ScalingRelations
Relations(Kumaran(Kumaran98)
98)•• Hydordynamics
Hydordynamics (Jenkins,(Jenkins,Richman
Richman85)
85)

Density,Density,Flow,
Flow,Temperature
TemperatureProfiles
Profiles(Brey(Breyetetal.01)
al.01)
Scaling
ScalingRelations
Relations(Lee(Lee95)95)•• LangevinLangevin equation
equation (Wakou,(Wakou,Ochiai,
Ochiai,Isobe
Isobe08)08)(which(whichaims
aimstotodescribe
describemacroscopic
macroscopicproperties)
properties)

CenterCenterofofMass,Mass,Global
GlobalTemperature
Temperature

MacroscopicAnalogy
Analogy with
with thethe Brownian
Brownian Motion
MotionAssumption:
Assumption:
FFbb consists
consists ofof thethestatic pressurepressure &&I static
1.
1. Random
Random force
force
2.
2. Frictional
Frictional force
force
3.
3. Elastic
Elastic force
force(D(QTEGCEVKPIQP
VJGDQVVQORNCVG

Analogy

TheoryLangevinLangevin Equation
Equation

Theory•• LangevinLangevin Equation
Equation ofof thethe COM
COM
2

d Z
dZ
2
= − Ω ( Z − Z ) + fs − µ
+R
2
dt
dt
1. Elastic ForceFrictional Force
2. Elastic Force
Random ForceElastic
ElasticForce
Force duedue to:to:1.
1. Expansion
Expansion andand contraction
contraction ofof thethe granular
granular gas
gas
2.
2. Excitation
Excitation ofof aa soundsound wave
waveˆ / τrel , µ = µˆ / τrel•• Macroscopic
Macroscopic Parameters:
Parameters: Ω = Ωττrelrel≡≡c/g
c/g(c:(c:thermal
thermalvelocity)
velocity)Properties
Properties ofof Force
Force

Theory•• Elastic
Elastic Force
Force duedue toto excitation
excitation ofof soundsound wave
wave
Hydrodynamic-wave theory
Hydrodynamicsoundsound-wave
theoryififvvbottom<< c, thenbottom << c, thengfs (t) ≅ σˆ A0ωcosωtc(if>> c, f might become nonlinear (shock-like) )(ifvvbottom
bottom >> c, fss might become nonlinear (shock-like) )•• Random
Random Force
Force

Gaussian
Gaussianwhite
whitenoise
noise
Fluctuation-Dissipation Relation
Fluctuation-Dissipation
RelationR( t )R( t' ) = 2µk BT / M × δ( t − t' )(T(Tisisgranular
granulartemperature)temperature)Solution
Solution ofof LangevinLangevin equation
equation

Theoryt

Z(t) − Z = A 0ζ sin(ωt + θ) + ∫ G(t − t′)R(t′)dt′−∞① Power injected by the bottom plate

(
σˆ
A ω)
Pb = Fb v 0 = −M ωZ v 0 = Mg
1+ O ω
ˆ −2
2
c
τrel


ˆ ≡ ωτ rel ≅>> 1 
ω
τ


2

(

( ))

② Power spectrum for the height of the COM
ICM(ω
ˆ ′)

2

2µˆ
π  Agζ 
c
′
′
ˆ +ω
ˆ ) + δ(ωˆ −ωˆ )) + 2= N 2  (δ(ω
3
ˆ −ω
Ng
2  c (Ω
ˆ ′2 )2 + (µˆ ωˆ ′)25Delta functional peaks at thefrequency of the bottom plate

Continuous spectrum

EventEventDriven MD
MD :1D:1D System
SystemEvent--Driven
\

㨯㨯㨯

\0

••
••

\
I



\
\
\
\

••
••
••
••

Model

N
N inelastic
inelastic particle
particle under
under gravity
gravity
Bottom
Bottom plate
plate oscillate
oscillate sinusoidally
sinusoidally
ZZ00(t)=Aωωtt(t)=A00sin
sinω

2
A
ω
Acceleration:
Acceleration:
Γ= 0g
Inelastic
Inelastic collisions
collisions between
between particles
particles
withwith thethe restitution
restitution coefficient
coefficient rr
Elastic
Elastic collisions
collisions between
between particles
particles
andand thethe bottom
bottom plate
plate
NoNo rotationalrotational degrees
degrees ofof freedom
freedom

S. Luding et al., Phys. Rev. E 49, 1634 (1994) and many papers.
M. Isobe and A. Ochiai, Mol. Sim. 33, 147 (2007).
J. Wakou, A. Ochiai and M. Isobe, J. Phys. Soc. Jpn. 77, 034402 (2008).

Macroscopic
Macroscopic Characterization
Characterization

Introduction

1D-Driven MD
1D Event
Event-Driven
MD (Luding(Luding et
et al.
al. 94)
94)
••

Order
-r)N
Order2Parameter:
Parameter:X=(1
X=(1-r)N

••

dZ
2
ω
A
0
M
=
−
Mg
+
F
Acceleration:Γ=bAcceleration:
2
g
dt

••

X>X
X>XCC ::Condensed
CondensedState
State

••

X<X
ΓΓsmall)
X<XCC::((Γsmall)Condensed
CondensedState
State
((Γ
ΓΓlarge)large)Fluidized
FluidizedState
State

••

Averaged
Averagedheightheightof
ofCOM
COM
2
V
in
inFluidized
FluidizedState
State
Z∝ 0

V0：Maximum Velocityof the vibrating plate

X

Kinetic
KineticTheory(2D,3D)
Theory(2D,3D) (Warr(Warretet al.
al.95)(Kumaran
95)(Kumaran98)
98)••

Theoretical
TheoreticalDerivation
Derivationfor
forthetheHeight
HeightofofCOM
COMPowerPowerSpectrum
Spectrumfor
forthetheheight
heightofof thetheCOM
COM

2µˆ
c5
ICM(ω
ˆ ′)
= 2
3
ˆ −ω
Ng continuous (Ωˆ ′2 )2 + (µˆ ωˆ ′)2ˆ : curve-fit parametersµˆ, Ωspectrum

Theoretica l prediction sˆ = 1 .5 .
with µˆ = 2.0 and Ω

The continuous part of the power spectra
at various system parameters collapse ona single master curve, which agree wellwith the theoretical prediction.

ResultPower spectrum ICM(ωω’)
without scaling

Behavior at high frequency
is out of our theory

Summary
Summaryof
of 1D
1DSystem
System with
withVibrating
VibratingWall
Wall

Discussion•• Comparison-DrivenComparisonwith
withthetheextensive
extensivemicroscopic
microscopicEvent
Event-Driven
MD
MDsimulation
simulationfor
foraawide
widerange
rangeparametersparameters•• ττ<<: fluidization<<ττrelrel : fluidization•• AA00ωω<<<<cc::hydrodynamichydrodynamicsound-wave
sound-wavetheory
theoryisisvalid
valid
••
••
••
••

Amplitude
Amplitudeof
ofOscillation
Oscillationof
ofthe
theCOMCOM
Power
PowerInjected
Injectedby
bythe
theBottom
BottomPower
PowerSpectrum
Spectrumfor
forthetheheight
heightof
ofthe
theCOMCOM
Kinetic
KineticEnergy
Energyof
ofthe
theCOM
COM(Violation(ViolationofofEquipartition)
Equipartition)•• Theoretical-driven
Theoreticalprediction
predictionisisgood
goodagreementagreementwith
withthetheevent
event-driven
MD
MDsimulation
simulationby
bythe
thefollowing
followingthree
threeparameters:
parameters:

ˆ = 1 .5σˆ = 1.5, µˆ = 2.0, ΩEventEventDriven MD
MD :2D:2D System
SystemEvent--Driven
N=200, r=0.996
2500

fraction•• PackingNparticles
N inelastic
inelastic
particles with
with diameter
diameter dd•• Periodic
Periodic boundary
boundary condition
condition inin thethehorizontal
horizontal direction
direction•• Bottom
Bottom plate
plate isis aa thermal
thermal wall
wall withwithtemperaturetemperature TT00
k B T0•• Driving
Driving Intensity:
Intensity: Λ =
mgd

2000

1500

1000

500

0

0

10 0

Model

0.0005•• Inelastic
Inelastic collisions
collisions between
between particles
particles
withwith thethe restitution
restitution coefficient
coefficient rr•• No
No rotationalrotational degrees
degrees ofof freedom
freedom. .

M. Isobe and H. Nakanishi, J. Phys. Soc. Jpn. 68, 2882 (1999).
M. Isobe, Phys. Rev. E 64, 031304 (2001).
J.Wakou and M.Isobe, AIP Conf. Proc. 1145,171 (2009).;1217,135 (2010).

LangevinLangevin Equation
Equation

Theory•• LangevinLangevin Equation
Equation ofof thethe COM
COM
2

d Z
dZ
2
= − Ω ( Z − Z ) + fs − µ
+R
2
dt
dt
1. Elastic ForceFrictional Force
2. Elastic Force
Random ForceElastic
ElasticForce
Force duedue to:to:1.
1. Expansion
Expansion andand contraction
contraction ofof thethe granular
granular gas
gas
2.
2. Excitation
Excitation ofof aa soundsound wave
waveˆ / τrel , µ = µˆ / τrel•• Macroscopic
Macroscopic Parameters:
Parameters: Ω = Ωττrelrel≡≡c/g
c/g(c:(c:thermal
thermalvelocity)
velocity)PowerPowerSpectrum
Spectrumof
ofthe
theheight
heightofof thetheCOMCOM
Theory:Theory:

µˆ
c5
ICM(ω
ˆ ′)
=
ˆ2 −ω
Ng3 (Ω
ˆ ′2 )2 + (µˆ ωˆ ′)2ResultContinuous spectrum

Nearly Elastic Case
Theoretical predictions
ˆ = 1.5.
with µˆ = 2.0 and Ω

In the case of nearly elastic particles, the power spectrum (with same
parameters as 1D case) agrees well with our theoretical prediction.

PowerPowerSpectrum
Spectrumof
ofthe
theheight
heightofof thetheCOMCOM
Theory:Theory:

µˆ
c5
ICM(ω
ˆ ′)
=
ˆ2 −ω
Ng3 (Ω
ˆ ′2 )2 + (µˆ ωˆ ′)2ResultContinuous spectrum

Largely Inelastic Case
Theoretical predictions
ˆ = 1.5.
with µˆ = 2.0 and ΩWhat is the causeof the deviation?

We find that systematic deviation from the theoretical predictions as
inelasticity increases , which have high intensity.

CauseCauseofofDeviation
Deviationat
at Large
LargeInelasticity
Inelasticity

Discussion

Failure of the Law of Equipartition

KCM / ( kBT / 2 )100

Λ=10001000
500
200
equipartition

100.340 * ( kBT / 2 )KCM: Kinetic energy of the COM
-1.07

1 N 1
k B T global ≡m vi2∑
N i =1 2

1

0.1
0.01

0.1
kBT / 2

1

KCM ≡

k T
1
MV 2 ≠ B global
2
2

We found the failure of the law of equipatition at low temperature.
This is not surprising because the system is far from equilibrium.
No theoretical predictions on the behavior of KCM / (kBT/2).

The
TheViolation
Violationof
ofthe
theLaw
Lawof
ofEquipartition
Equipartition

Discussion

Development of the density inversion(Simulation: Isobe&Nakanishi 99; Theory: Brey et al. 01, Bromberg et al. 03)

π1 / 2
2 1/ 2
Nh: number of layers of particles at rest
λ≡Nh 1 − r2
Density inversion occurs at λ > λc()

λ C = 1.06569

Violation
Violationof
of Fluctuation
FluctuationFluctuation--Dissipation
DissipationTheorem
Theorem

Discussion

Fluctuation-Dissipation Theorem (of the first kind)

R(t)R(t' ) = I δ(t − t' )Langevin equationk B Tglobal
M
K CM =
I=
4µ2⇒

I = 2 µ kT global M

Equipartition

Is the Langevin equation still useful?
4µ
We put I =K CM and determine KCM from simulation dataM

重心運動についての
重心運動
についての
重心運動についてのについての現象論
現象論
重心運動についての現象論についての現象論
J.Wakou, A.Ochiai, M.Isobe
J.Phys.Soc.Jpn.77(08)034402

重心高さZに対するランジュバン方程式
d2Z
Ω： 固有振動の角振動数
2
=
−
Ω
(
Z
−
Z
)
−
µ
V
+
r
μ： 摩擦係数
dt 2
摩擦力
摩擦力 ランダム力
ランダム力

線形弾性力

r ( t )r ( t' ) = I δ( t − t' ) ただし、 I = 2µk B TCM M

重心高さの
パワースペクトル理論：
理論
 c  3µˆ
ICM (ωˆ )   2MkBTCM  =
2
2
ˆ2 −ω
ˆ 2 + (µˆ ωˆ) g  Ω()

ˆ = 1.5
フィッティングパラメタ µˆ = 0.6, Ω

揺動散逸関係においては
揺動散逸関係においては、においては、TCMが適切な
適切な「温度」
温度」変数

２２次元系の次元系
結果
次元系のの結果のまとめ
結果のまとめのまとめ
次元系の結果のまとめ•• 重心運動
重心運動の
と粉粒体の
粉粒体
について
重心運動の
の温度
温度T
の温度
温度T
について、
エネル
重心運動の
温度TCM
粉粒体の
温度Tggについて、について、、エネル
CMと粉粒体のギー等分配則≠T が観測された
ギー
等分配則
観測
された
ギー等分配則
等分配則の
の破
破れ
れTTCM
された。
ギー等分配則の等分配則の観測された。
された。。
CM≠Tggが観測された•• エネルギー
エネルギー等分配則
等分配則
エネルギー等分配則
等分配則の
の破
れは
エネルギー等分配則の
等分配則の
破れは
密度反転現象が密度反転現象密度反転現象がが進進むにつれてむにつれて
密度反転現象が
顕著になる
顕著
になる
顕著になる
になる。
顕著になる。
になる。。
•• 密度反転
密度反転が
きても
重心
密度反転が
が起
きても重心
重心の
の速
速
密度反転が
起きても重心きても重心の
重心の
度分布は
度分布
ガウス
分布
度分布は
はガウス分布
ガウス分布
分布。
度分布は
ガウス分布。
分布。。
•• 重心運動
重心運動をするランジュバン
種
記述
する
ランジュバン
方程式
えた
。（
重心運動を
を記述する
記述する
するランジュバン
ランジュバン方程式
方程式を
を考
考えた。（
えた。（
。（第
第22種
種））
重心運動を
記述
するランジュバン方程式
ランジュバン方程式を
方程式を
えた。（第
。（第
揺動散逸関係においてを用いるとシミュレーション
揺動散逸関係いるとシミュレーション
結果
揺動散逸関係においてにおいてT
シミュレーション結果
結果とと揺動散逸関係においてにおいてTCMいるとシミュレーション結果
シミュレーション結果と結果とCMを用いるとシミュレーション
ある程度良
ある
程度良
一致
ある程度良
ある程度良い
程度良い
い一致。
一致。
一致。。
程度良い

Concluding
Concluding Remarks
Remarks

Conclusion

 AALangevinLangevinequation
equationfor
forthethemotion
motionofofCOM
COMwaswasformulated
formulatedto
todescribe
describe
thethemacroscopic
macroscopicproperties
propertiesof
ofthe
thefluidized
fluidizedstate
stateofofgranular
granulargas
gasunderundergravity.
gravity. Comparison
Comparisonwith
withthetheresultresultof
ofextensive
extensivesimulations:
simulations:Vibrating Wall

Thermal Wall

1D System

excitation of the soundwave (*)

Under
considerationGood
agreement

2D System

Under consideration

No density inversion

Under consideration

density inversion (**)

(nearly elastic)

2D System(Large inelasticity)

Deviation

 TheTheoriginoriginofofViolation
Violationof
ofthe
theLaw
Lawof
ofEquipartition
Equipartitionisis
different
differentbetween
between(*)(*)andand(**)(**)
 Future
FutureWorks:
Works:Extended
Extendedtheoreticaltheoreticalframework
frameworkfor
forthethepresence
presenceofofdensity
densityinversion
inversion(modified(modifiedFluctuation-dissipation
Fluctuation-dissipationrelation)
relation)

今後の
今後
今後の
の展望
展望
今後の

AlderとWainwrightによって拓
によって拓かれたイベントドリブン型分子
イベントドリブン型分子動力学法は、現在、現在、数百万粒子系の数百万粒子系の分子動力学シミュレーショ
分子動力学シミュレーショ
ンがパソコンレベルでパソコンレベルで可能になった可能になった。
になった。
 分子動力学法は
分子動力学法は、ミクロなミクロな粒子レベル
粒子レベルとレベルとマクロなマクロな現象論（
現象論（流体力学）力学）をつなげる 非平衡統計物理学の非平衡統計物理学の構築（
構築（粉体系への
粉体系への応用への応用など）など）に不可欠な不可欠な重要な重要なツールとしてますますツールとしてますます重要としてますます重要になっていく
重要になっていくと思われる。
われる。


御静聴、御静聴、ありがとうございました。
ありがとうございました。