数理解析研究所講究録 1305 巻 2003 年 68-72

68
非弾性剛体球系の一様冷却状態のダイナミクス
九州大学理学研究院中西秀 (Hiizu Nakanishi)

Department of Physics, Kyushu University
序:

粉粒休系の単純なモデルとして、非弾性粒子系がしばしば用いられる。その中でももっとも簡単な、無重力中で外部からの駆動のない非弾性剛体系の動力学を調べる。
弾性剛体球系はボルツマン気体と呼ばれ、 19 世紀の統計力学の揺藍期には気体のモデルとして研究された。それが、マックスウェル・ボルツマン分布など平衡の統計力学の確立に重要な役割を果たしたことは良く知られている。 ここでの問題意識は、ボルツマン気体のエネルギー保存則が少し破れた時にどのようなことが起こるか、 ということで、粉粒体系の振る舞いとの関連以外に、 統計力学的興味も大きい。
非弾性剛体球系 :ここで調べる非弾性ボルツマン気体系は、 同じ直径および質量の剛体球系からなり、粒子間衝突は法線方向の反発係数 のみによって特徴づけられるとする。
即ち、回転の自由度は考えない。衝突する 2 粒子の衝突直前および直後の速度を、 それぞとすると、衝突則は、および ,れ$r$$\vec{v}_{1},$$\mathfrak{h}$$v_{1}^{\neg}$$v_{2}^{\neg}$$v_{1}^{\prec}$$=\vec{v}_{1}-1\star r(\vec{n}\cdot\vec{v}_{12})\vec{n}$

(1)

$\{$

$v_{2}^{\neg}$$=\vec{v}_{2}+\mathrm{T}(\underline{1}+\underline{r}\vec{n}\cdot\overline{v}_{12})\overline{n}$となる。 ここで、 n\rightar ow は衝突の瞬間の 2 粒子の中心を結ぶ方向を向いた単位ベクトルで、である。
系の初期状態としては弾性系の平衡状態にとり、そこから時間発展させる。 しばらくの間、系は一様なまま非弾性衝突によって徐々にエネルギーを失って行く (一様冷却状態)$\vec{v}_{12}\equiv\vec{v}_{1}-\vec{v}_{2}$[1]。 しかし、非弾性衝突によって相対速度を失う為に、速度場の非一様性が生じ渦構造が見えるようになる (ずれ不安定性) [2]。 さらに時間がたつと、密度場の非一様性が発達し (クラスター不安定性) [3]、密度の濃い領域の粒子クラスターが互いに衝突したり分裂一定反発係数の剛休球モデルでは、いずれ、密度の高い領域で、有限時間内に無限衝突が生じる、いわゆる非弾性崩壊が起こり $[5, 6]$ このモデルの範囲内では、したりする$[4]_{\text{。}}$、それ以上時間発展を追いかけることができなくなる。
一様冷却状態での温度の冷却は Haff 則と呼ばれ、その時間依存性は、以下のように簡単に導かれる。 まず、 粉体温度 を、 通常の分子運動の温度に習って、一様冷却状態 :$T$$\frac{d}{2}T=\frac{1}{2}m<v>\triangleleft$と定義する。 ここで、 は空間の次元、 $m$ は粒子の質量で、 $<\ldots>$ は統計平均を表す。反発係数 で非弾性衝突するたびに法線方向の相対速度が 倍に減少するので、粒子の運動
$d$

$r$

$r$

エネルギーは、一回衝突あたり割合として$2 \gamma_{0}=\frac{1-r^{2}}{d}$69
図 1: エネルギーの依存性。反発係数
実線はr=0.9、粒子数密度$\hat{\mathrm{o}\check{\backslash \wedge \mathrm{h}\check{\mathrm{I}\mathrm{r}}\triangleright}}$$\tau$$n=0.25_{\text{。}}$

通常の剛体系のエネルギー減衰を示し、破線は速度交換した系のもの (本文参照)。細線は Haff 則 (4) を示すが、破線と重ならないように少しずらせてある。
$\tau$だけ減少する。 その結果、衝突頻度を度は、$\omega$とすると、系が一様な場合には、平均として温(2)$\frac{dT(t)}{dt}=-2\gamma_{0^{\{}}vT(t)$に従って減少するであろう。一方、衝突頻度比例するので、 (2) 式は、$\omega$

は粒子の平均速度に比例し、それは$\sqrt{T}$に$\frac{dT(t)}{dt}=-CT(t)^{3/2}$となる。 ここで、 $C$ は時間によらない比例係数。 この微分方程式の解は直ちに求まり、(3)$T(t)= \frac{T_{0}}{(1+t/t_{0})^{2}}$を得る。時間が経つとともに粒子速度がゆつくりになって行くので、時間を実際の時間ではなく個々の粒子の平均衝突回数 ではかると都合の良いことが多い。 と$t$$\tau$$\tau$$t$の関係は、$\omega dt=d\tau$で与えられるが は時間の関数なので、 両者は比例関係ではない。
で時間をはかると、 Haff 則 (3) は単純な指数関数$\omega$$\tau$(4)$T(\tau)=T_{0}\exp[-2\gamma_{0}\tau]$で表される。
実際のシミュレーションを見てみると、 確かに最初は指数関数的に減衰し、 (4) に良くあっているが、 クラスター不安定性が生じるとともに減衰は緩やかになる (図 y。

ボルツマン・エンスコッグ方程式 : 一様冷却状態での一粒子分布関数相関を無視すると、 ボルツマン. エンスコッグ方程式、$\rho(\vec{r},\vec{v}, t)$

は、 2 体$\frac{\partial\rho(\vec{r}_{1},\tilde{v}_{1},t)}{\partial t}+\vec{v}_{1}\cdot\vec{\nabla}_{1}\rho+\frac{\vec{f_{1}}}{m}\cdot\frac{\partial\rho(\vec{r}_{1},\vec{v}_{1},t)}{\partial\vec{v}_{1}}$$= \chi\sigma^{d-1}\int d\vec{v}_{2}\int_{(\tilde{v}_{12}\cdot\hat{\sigma})>0}d\hat{\sigma}(\vec{v}_{12}\cdot\hat{\sigma})[\rho(\vec{r}_{1}, v_{1}’)\prec\rho(\tilde{r}_{1}, v_{2}’)\prec-\rho(\vec{r}_{1},\vec{v}_{1})\rho(\vec{r}_{1},\vec{v}_{2})]$

(5)

70

Sonine 係数 $a\ell(2\leq\ell\leq 5)$ の

図 2:

時間依存性。反発係数は r=0.9(上),095(中), 098(下) で、密度は全て $n=$0.25。データは 100 回の平均値を表す。

で与えられる。 ここで、 は外力、 は接触距離での 2 体相関、 は粒子直径である。まは衝突後に ,た、 は衝突時の 2 粒子の相対位置の方向を示す単位ベクトルで、$\vec{f_{1}}$$\sigma$$\chi$$\hat{\sigma}$$\ovalbox{\tt\small REJECT}$$\vec{v}_{1}$$\overline{v}_{1)}^{\mathrm{v}\prime}v_{2}^{\tau\prime}$になる 2 粒子の衝突前速度である。

弾性系の場合には、 この方程式の一様定常解としてマックスウェル・ボルツマン分布が得られる。
非弾性系の場合には、 もちろん定常解はないが、平均粒子速度 $v_{0}(t)\equiv\sqrt{<v^{2}>/d}$ でに対する分布関数 ( , t)、スケールした速度ど$\vec{c}$$\equiv\overline{v}/v_{0}(t)$$\tilde{\rho}$(6)$\rho(\vec{v},t)\equiv\frac{n}{v_{0}(t)^{d}}\overline{\rho}(\vec{C}_{)}t)$は、定常になるかも知れない。 ここで、 は粒子密度で、分布の位置依存性は無視した。
のみにあり、 の 依存性がないとして、 (6) 式をボルツ
速度分布の時間依存性は
またそれは、マン. エンスコッグ方程式 (5) に代入すると定常解を得ることができ $[7, 8]$速度の大きいところで指数関数的な振る舞いをすることが予想されて [9]、実際、 シミュレーションでもそのような傾向が現れている [11]。
$n$$v_{0}(t)$$\tilde{\rho}$

$t$、

71

図 3: a2 の最小値の反発係数依存性。
粒子数密度は n=0.25。ボルツマン. エ蕾$[mathring]_{\mathrm{Q}}^{1}$$r$

ンスコッグ方程式の最低近似の結果 [8]を破線、 6 次の近似の結果 [11] を実線
は速度交換をした系の値。
で示す。
$(\cross)$

ー様冷却状態での速度分布関数 :速度分布の時間発展を調べるには、分布関数 を直交$\overline{\rho}$関数系に展開してその係数を調べるのが便利である。分布がガウス関数からそんなに外れていない時には、 Sonine 多項式を用いて、(7)$\tilde{\rho}(\vec{c}, t)=\frac{1}{\prime\pi}e^{-c^{2}}\sum_{\ell=0}^{\infty}a_{\ell}(t)S_{\ell}(c^{2})$のように展開する。 ここで、
$S_{0}(x)=1$ ,$S_{\ell}(x)$は Sonine 多項式で、 低次のものは、$S_{1}(x)=-x+ \frac{1}{d}$,$S_{2}(x)= \frac{1}{2}x^{2}-\frac{1}{2}(d+2)x$ 十 $\frac{1}{8}d(d+2)$,が分布の形を与えるが、 全確率が 1 になるこ. とから a0=1、 速度ので与えられる。
スケールから $a_{1}=0$ となるので、 2 次より高次の係数が実際の分布の形の情報を与える。の時間発展は、展開 (7) をボルツマン. エンスコッグ方程式 (5) に代入することによって得られ、得られた方程式が数値的に解かれている $[10, 11]$ 。その結果、 マツクスウエル$a\ell(t)$$a\ell$

分布、 即ち $a_{\ell}=0$ (\ell >0)、 を初期状態をとっても、一粒子あたり数回の衝突という非常に短い時間のうちに、ゼロでない定常値に落ち着くことが示されている。即ち、非弾性系
でも、一様冷却状態ではスケールされた分布関数は時間によらない定常な形に極めて短時間で落ち着き、 それはマックスウエル分布から外れている、 ということがボルツマン. エンスコッグ方程式から得られている。
しかしながら、実際に 2 次元系でシミュレーションをしてみる
シミュレーション結果 :は $\tau\leq 10$ である値に達するが、 そこで定常にはと、 図 2 のような結果を得る。即ち、ならずに、 それからクラスター不安定性が現れる時刻まで、徐々にゼロにむかつて絶対値の最大が減少して行く。 この傾向は、反発係数 を変えてみても同様である。初期の値は、 が 1 に近い時には、ボルツマン. エンスコッグ方程式で得られる値に近い (図 3)。
$a_{\ell}$$r$$|a\ell|$$r$

この$a_{\ell}$が一様冷却状態でも定常値を示さないという性質は、 ボルツマン. エンスコッグ方程式で予想される振る舞いとは異なり、そこで無視されたものの効果、即ち速度相関によるもののはずである。 このことを確かめる為に、剛体球系のダイナミクスを人為的に変更して、一定時間間隔で粒子の速度をランダムに交換して、計算機実験を行なった。即ち、速度分布は保存したまま、速度相関を取り除いたのである。 その場合の a2 の時間発展を、 $r=0.9$ でいくつかの密度での振る舞いとともに、 図 4 に示す。 どの密度でも、 も72

図 4: Sonine 係数 a2 の時間依存性。反発係数は $r=0.9$ で、密度は $n=0.25$ ,0.1111, 00625。速度交換を行なったものを破線で示す。データは 100 回の平均値を表す。
$\tau$とのダイナミクスで計算したものは初期緩和の後に徐々に絶対値が減少しているが、速度相関を乱したものは定常値に落ち着き、 その値もボルツマン. エンスコッグ方程式による予想値に近い (図 3)。
結論 :ボルツマン. エンスコッグ方程式に基づいて、非弾性系であっても一様冷却状態

では、平均速度でスケールした一体分布関数は定常形を持ち、その形はマックスウェル分布とは異なるものになると信じられていた。本研究では、 2 次元系でのシミュレーションを行なうことにより、実際には分布関数は定常形を持たないことを示した。 このことは、非弾性形においては、一様状態であっても速度相関が発達するので、 2 体相関を無視したボルツマン. エンスコッグ方程式の有効性は、速度相関のない弾性平衡系の場合に比べてかなり限られていることを示唆する。 より詳しい結果は、 [12] に発表されている。

参考文献. Haff, J. Fluid Mech. 134 (1983) 401.
. Ernst, Phys. Rev.
[2] T.P.C. van Noije and[1]$\mathrm{P}.\mathrm{K}$$\mathrm{M}.\mathrm{H}$$\mathrm{E},$$61$(2000) 1765.

[3] I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70 (1993) 1619.

. Herrmann, Chaos, 9(1999) 673.
(1992) 496.
. Young, Phys. Fluids. [5] S. McNamara and[4] S. Luding and$\mathrm{H}.\mathrm{J}$$\mathrm{W}.\mathrm{R}$[6] S. McNamara and$\mathrm{A},$$4$. Young, Phys. Rev. $\mathrm{E}50(1994)$ R28.

$\mathrm{W}.\mathrm{R}$[7] A. Goldshtein and M. Shapiro, J. Fluid Mech. 282 (1995) 75.

[8] T.P.C. van Noije and[9][10]$\mathrm{S}.\mathrm{E}$$\mathrm{M}.\mathrm{H}$. Ernst, Granular Matter $1(1998)$

$57$

.

. Esipov and T. P\"oschel, J. Stat. Phys. 86 (1997) 1385.

$\mathrm{N}.\mathrm{V}$. Brilliantov and T. P\"oschel, Phys. Rev.

[11] M. Huthmann,$\mathrm{J}.\mathrm{A}$$\mathrm{E}61$(2000) 2809.

G. Orza, and R. Brito, Granular Matter 2(2000) 189.

[12] H. Nakanishi, Phys. Rev.

$\mathrm{E}$, to be published (2002).