数理解析研究所講究録 1413 巻 2005 年 176-185

176

Lennard-Jones 粒子系における気液界面構造
小串典子 (Fumiko Ogushi)
湯川論 (Satoshi Yukawa)伊藤伸泰 (Nobuyasu Ito)伊藤伸泰研究室

東京大学大学院工学系研究科物理工学専攻

Department of Applied Physics, School of Engineering,
The Univercity of Tokyo概要
3 次元 Lennard-JOn6 粒子系における気液界面構造について, 非平衡分子動力学法を用いて研究する。
系は $L_{X}\cross L_{Y}\mathrm{x}Lz$ の立方体てあり, この中に Lennard-Jones ポテンシャノレを持つ粒子が詰められてい熱浴をつけ, 部分的に温度を制御するる。 方向の両端から一定の幅の領域に異なる温度の方向については弾性壁とし, y-, z- 方向には周期境界とすことて系に温度勾配を与える。 境界条件は,る。 粒子の密度及び熱浴の温度を調整することて, 系は単相の状態も混相の状態も実現することが可能てある。
気泡の臨界核は数ミクロン以下てあり, このような大きさの系における熱伝導は明らか ない。 そこてます, 粒子数数百程度の系においてどのような熱伝導が生じているかを調べる。 臨界点から 200% 程度
高温・高密度の超臨界流状態における熱伝導率の系のサイズ に対する依存性を調べたところ, 剛体粒子系などと同様に $1/\sqrt{L}$ のサイズ依存性を示し, 粒子の大きさに対して数百程度の長さでは巨視的な系における熱伝導率から 1 割程度減少することが分かった。
次に, 気液共存状態を実現し界面における密度変化の様子を調べた。気相, 液相単相における密度に対して 1/3 乗て比例する長さを考え, この長さを単位として界面ての密度変化の式を考えたところ, これ$\mathrm{N}\mathrm{o}\mathrm{s}\mathrm{e}’- \mathrm{H}\mathrm{o}\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}$

$x$

$x$

$\grave{\text{て}}$$L$とシミュレーションの結果が良い一致を見せた。

1

はじめに沸騰現象は身近な現象であり,さらに熱機関などにおいては非常に重要な現象である。
沸騰現象を説明する理論に古典核生成論がある [1]。 気泡核はその半径が一定の大きさ r。
以下では成長することが出来ず潰れ, r。以上の半径を持つ場合は成長する。 この半径 r。
を臨界核半径という。 古典核生或論では臨界核半径を気泡のある場合とない場合の自由エネルギー差から求める。 3 次元の場合, 気泡核が生成された場合の自由エネルギー $f(r)$は表面とバルクの自由エネルギーから次のように計算される。

(1)$f(r)=4 \pi r^{2}f_{S}-\frac{4}{3}\pi r^{3}f_{v}$は単位体積あたりの自由エネルギーは単位面積あたりの自由エネルギー,であることである。 臨界核では $df/dr=0$ であるのでこれから臨界核半径はが分かる。 古典核生成論は, 微視的には界面の構造を考慮しておらず不十分であり, 実験値とすれることが知られている。 コンピューターシミュレーションや実験, また様々な理論的なアプローチ [2] も行われているが, 未だ核生成については十分な理解は得られていここで,$f_{S}$$f_{v}$$r_{\mathrm{c}}=2f_{8}/f_{v}$ない。 Oxtoby ら [3] は, 自由エネルギーを求める際に界面の構造を考慮し古典核生成論を補正することで, 古典核生成論よりも実際の値に近い結果を得ている。 ここで, 界面における密度変化は密度汎関数法 [4, 5, 6] から求められる。
このように, 沸騰現象を理解するためには臨界核近傍での気泡核の振る舞いを理解する必要があり, そのためには気液界面の構造を理解することが重要である。
沸騰現象は, 非平衡定常状態の熱伝導が重要な役割を果たしている現象の一つである。
巨視的な系における熱伝導現象は Fourier の法則に従っており, 熱伝導率は熱流と温度勾177
Heat BathHeat BathLocal box

図 1: 系の形状配の比で表されるが,微視的なダイナミクスに基づく理論的理解は未だ十分得られておらす, 臨界核近傍及びこれ以下の大きさの系では, どのような熱伝導が生じているか明らかてないため,ますこれについて調べる必要がある。
ある種のカオスモデルでは Fourier 則が実現することが知られており [7, 8, 9], 剛体粒子系 $[10, 11]$ や Fer 面-Pasta-Ulam 格子模型 $[8, 12]$ を用いた研究から, 熱伝導率が系の-Uram 格
サイズ及び次元に依存することも知られている。 島田ら [8] は, Fer 面子模型及ひ剛体球系を用いて, 熱伝導率が系のサイズ 及び次元 に依存しており, 3 次$\beta$$- \mathrm{P}\mathrm{a}\mathrm{s}\mathrm{t}\cdot \mathrm{a}$$\beta$$d$$L$の依存性を持つことを示した。 この依存性は久保公式 $[13, 14]$ から求めら
元系ではにはれる熱伝導率の依存性と一致している。 久保公式から, 自己相関関数の依存性をの依存性を持ち 3 次元では従って減少する。 従って, 熱伝導率は$1/\sqrt{L}$$t^{-^{\underline{d}}}\sim$$\langle j(0)j(t)\rangle$’$\sqrt{L}$$L^{-\frac{d}{2}+1}$持つ。 ここで次元 は 2 以上である。
村上らは [11], 非平衡分子動力学法を用いて剛体粒子系における熱伝導現象を調べ, 3次元では $1/\sqrt{L}$ のサイズ依存性を持ち, 固相・流体相の共存する系では, 各相の温度勾配$d$は異なる傾きを示した。
同様に気相・液相の共存する系について調べるためには, 剛体粒子系では気相と液相の区別をつけることは出来ないので, 引力相互作用を持つ粒子モデルを用いる。一般に, 気相や液相を扱う場合 Lennard-Jones ポテンシャルが用 1) られる。 Lennard-Jones 粒子を用
次元系において臨界点近くでのいた気泡核に関する研究も多くなされており $[15, 16]_{0}$$2$熱伝導率の挙動について調べたものもある [17]。
臨界核半径の大きさは数ミクロン以下であり, このような犬きさの系では上で述べたように熱伝導そのものも未だ明らかでない。 さらに, Lennard-Jones 相互作用は剛体粒子や非線形格子に比べ相互作用範囲が長く: 局所的な構造はさらに複雑である。 従って,とは異なる依存性を持つ可能性
Lennard-Jones 粒子系特有の依存性を持ち, 単純な$1/\sqrt{L}$もある。 そこで,まず粒子数数百粒子程度の系での熱伝導率につ 1) て剛体粒子系などと同様の依存性を持つかについて調べ, 次に気液共存系を実現し, その界面における密度変化について調べる。

178
2model系は立方体 $L_{X}\cross L_{Y}\cross Lz$ であり, その中に Lennard-Jones ポテンシャルを持つ粒子を配置する。 方向に一定の幅で系を分割し, それぞれのセルの中で局所的な温度・密度を計算する。 系の形状を Fig. 1 に示す
粒子 と の間における Lennard-Jones(12-6) ポテンシャルは次のように定義される。
$x$

$i$

$j$

(2)

$\phi(r_{ij})=4\epsilon\{(\frac{\sigma}{r_{ij}})^{12}-(\frac{\sigma}{r_{ij}})^{6}\}$は粒子と粒子 間の距離である。

は粒子の直径,はポテンシャルの強さを決める。
及ひ各粒子の質量はそれぞれ 1 とし, 計算時間の短縮のため, cut-Off を 3 とした。境界条件は,方向は, 系の両端に Lennard-Jones ポテンシャルの斥力部分を持つ弾性壁をつ {九及び 方向には周期境界を課した。 系に温度勾配を与えるため,方向の両端から一定幅の領域に異なる温度の Nos\’e-Hoover 熱浴 [18] をつけた。 熱浴の温度は, 一方を高温, もう一方を低温とする。 熱浴における運動方程式を以下に示す。

$r_{ij}$$\sigma,$

$i$

$j$

$\sigma$$\epsilon$$\epsilon$

$x$

$y$

$x$

$z$

$T_{H}$$T_{L}$$m \frac{d^{2}r_{i}}{dt^{2}}=F_{i}+\xi m\frac{dr_{i}}{dt}$(3)$\frac{d\xi}{dt}=\frac{2(K_{i}-K_{\mathrm{H}\mathrm{B}})}{\tau^{2}}$(4)

そしては熱浴内にいる粒子の番号である。
はそれぞれ熱浴変数, 時定数, 粒子 の運動エネルギーそして熱浴 (of ) 内にいる粒子のもつ運動エネルギーである。
は以下全てのシミュレーションにおいて 1 とした。 系の更新には Verlet 法を用い, 熱浴変数 の計算には Euler 差分を用いた。 これら二つの差分法は一定の時間刻みを用$i$$\xi,$$i$$T_{H}$$\tau,$$K_{i}$$K_{\mathrm{H}\mathrm{B}}$$T_{L}$$\tau$$10^{-4}$$\xi$

及びいて計算した。温度及び密度はそれぞれを用いて無次元化した値を用いる。
系は Verlet 法を用いて更新し, 熱浴変数の計算には Euler 差分を用いる。 どちらの差分法にも一定の時間刻みを用いる。
$\sigma^{-2}$$\epsilon/k_{B}$$10^{-4}$

331

熱伝導率
シミュレーション

熱浴の温度をそれぞれ $T_{H}=4.0,$ $T_{L}=2.8$ と , 系の平均粒子密度 は 0446, 0562, 及
び 0652 とする。 3 次元 Lenn と -Jones 粒子系における臨界点の値 [20] はおよそ温度 132,密度 0302 と知られている。パラメーターは臨界点よりおよそ 200% 高温・高密度の領域$\llcorner$$\rho$$\mathrm{d}$にある。 この温度域では, 固相には 13 以上の密度が必要であるので, このとき系は超臨は 4 に固定し, 平界流状態にある。熱伝導率を求めるシミュレーションでは,及び均密度一定のまま,を変化させる。初期条件では, 粒子は立方格子に配置する。 また,物理量は 5 本の独立したサンプルから求め, 誤差は標準偏差を用いている。
$L_{Y}$$L_{X}$$L_{Z}$

179

図 2:

32

$x$

方向における温度変化 : 温度分布は温度勾配の方向に線形になっている

局所温度

局所平衡は各粒子の熱化のスケールと巨視的な質量流れのスケールが異なっていることに基づいている。 しかし, このことは微視的な系においては明らかではない。村上ら [11]は非平衡定常状態での剛体球系において局所平衡が成立するかどうかについて調べ, 局所的な速度分布が近似的に Maxwe 垣分布になっていることを確認した。 このことから, 局
所平衡は成立していると考えられ, 局所的な温度を平衡状態における温度同様に定義することが出来る。 Lennard-Jones 粒子系においても同様に定義することが出来ると考え, 各セルにおける局所温度を次のように定義する。各セルにおける局所温度を次のように定義する。

(5)$T(x)=m \frac{\langle\vec{v_{i}}^{2}\rangle_{\mathrm{c}\mathrm{e}11}}{3k_{\mathrm{B}}}$

ここで$k_{\mathrm{B}},\vec{v_{i}}$はそれぞれ Boltzmann 定数及び粒子 の速度である。
$i$$\langle\ldots\rangle_{\mathrm{c}\mathrm{e}11}$は局所セルにおける長時間平均である。 Fig. 2 に局所温度を示す。 温度分布は線形になっており, 温度勾配は一定になっていることが分かる。

33

有限サイズにおける熱伝導率$\kappa(L_{X})$

次に熱浴変数が系にした全仕事から熱伝導率を求める。時刻 $t=0$ から までの間に, 熱を次のように定義する。
浴変数が系に働いた全仕事 $E_{H}(t)$ 及び$t$$E_{L}(t)$$E_{H}(t)= \int_{0}^{t}dW_{H}$(6)$E_{L}(t)= \int_{0}^{t}dW_{L}$(7).

はそれぞれ高温側, 低温側の仕事である。 $dW_{H}$ および $dW_{L}$ は単位時間当たりに熱浴変数が系にした仕事である。 $dW_{H}$ 及び $dW_{L}$ はそれぞれ高温側, 低温側の熱浴$E_{H}(t),$$E_{L}(t)$変数による仕事であり以下のように定義する。
$dW_{H}=- \sum_{i}\xi_{H}P_{i}dr_{i}$

(8)

180

図 3: 熱浴変数の系にする全仕事 : 熱浴側の仕事はそれぞれ正・負の一定の傾きを持ち, 系には定常な熱が流れている(9)$dW_{L}=- \sum_{i}\xi_{L}P_{i}dr_{i}$は単位時間に粒子 が動いは熱浴内にいる粒子の番号てある。 は粒子 の運動量,た距離である。
系に定常な熱が流れている場合, 高温側の熱浴変数のした仕事 $dW_{H}$ は系に一定の熱エネルギーを与え, $dW_{H}$ は正の値を持つ。低温側の熱浴変数のした仕事 $dW_{L}$ は系から一定の熱エネルギーを吸い取り,は負の値を持つ。 このとき, $E_{H}(t)$ 及び $E_{L}(t)$ を時間ににあたり, それぞれ正の値, 負の値を持つことが分かる。定常
対する変化率は状態では熱の収支は釣り合っており, $E(t)$ の傾きの絶対値は高温側 低温側で同じ値を$i$$P_{i}$$i$$i$$dr_{i}$$d\nu V_{L}$$d7\mathrm{T}^{\gamma_{H}},dW_{L}$$\circ$とる。
及ひ $E_{L}(t)$ の振る舞いを Fig. 3 に示す。 時刻 $t=0$ の時, 系は初期条件により与えられた粒子配置と温度を持っている。 この初期状態は, 定常状態における粒子配置, 温度$E_{H}(t,)$

分布とは異なっており, 定常状態において系の持つエネルギーと比べ余分なエネルギーを持っている。従ってシミュレーションが始まると, 系は先ず余分なエネルギーを捨て (あるいは吸収し) 初期状態から定常状態に近づいていく。 Fig. 3 を見ると, $E(t)$ は $E(0)=0$から高温側・低温側ともに急激に減少し, 負になっている。 これは, 初期状態において系が定常状態と比べ余分なエネルギーを持っており, シミュレーションが始まると, 急激な吸熱を生じながら余分なエネルギーを捨て, 定常状態に近づいていると考えられる。 この初期状態からの緩和の後, $E_{H}(t)$ は正の傾きを持ち, $E_{L}(t)$ は負の傾きを持つようになっている。 この時の傾きの値は, その絶対値がほぼ同じ値となっており, これは, 系には定常な熱が流れていることを示している。
次に $E_{H}(t),$ $E_{L}(t)$ から熱伝導率を求める。 熱伝導率の比て表される。

$\kappa(L_{X})$は熱流$J$と温度勾配 $dT/dx$(10)$\kappa(L_{X})=J/\frac{dT}{dx}$定常状態での値から, 熱伝導率を求めるため異なるスタート時刻$t_{0}$を用いて$\kappa(L_{X}, t, t_{0})$を定義する。
$\kappa(L_{X,\backslash }tt_{0})’=(\frac{dT}{dx})^{-1}\cdot\frac{E(t)-E(t_{0})}{t-t_{0}}$(11)181) : 系が定常状態に達した後のデータから求められた熱から求められた熱伝導率 \kappa (L い ,伝導率は似たような振る舞いを示し, いつれも一定値に近づく
図 4: 異なる$t,$$t_{0}$$t_{0}$

図 5: 熱伝導率のサイズ依存性 : 熱伝導率は $1/\sqrt{L}$ の依存性を示し, サイズが数百程度では巨視的な系における熱伝導率から 1 割程度減少するまで間に系に流れる熱流に当たる。 $E(t,)$ 及び $E(t_{0})$ の値を示す。定常状態に達には $E_{H}(t)$ と $E_{L}(t)$ の絶対値の平均を用いる。 Fig. 4 に$E(t)-E(t_{0})/t-t_{0}$ は時刻$t_{0}$

から$t$$\kappa(L_{X}, t, t_{0})$

した後では熱流は一定であり, 定常状態のデータから求めた$\kappa(L_{X}, t, t_{0})$は$t_{0}$によらずほほ同様の振る舞いを示しながら, 一定値に近づいていく。従って $E(t)$ から求めたに依存した振る舞いを示さなくなった時点で, 系は定常状態に達していると考えらがのから, 有限系における熱伝導率 \kappa (L え) はれ, 定常状態における$\kappa.(L_{X\backslash }, t, t_{0})$$t_{0}$$\kappa(L_{J}\backslash ’, t, t_{0})$$\kappa(L_{X}, t, t_{0})$$tarrow\infty$34での値として求められる。

結果との振る舞いをについて
以 [の解析から求められた熱伝導率仮定しこの 3 パラメータを求める。 その結果, /3 の値は 05 と考えられ, この結果は剛体粒子系における熱伝導率のサイズ依存性 $\kappa+a/\sqrt{L_{X}}[8,11]$ と一致している。 $l\mathrm{J}=0.5$ と$\kappa+a/L_{X}^{\beta}$$\kappa(L_{X})$$\kappa(L_{X})$及び を用いて fit した結果を Fig. 5 に示す パラメータは一定値である。 パラメータ と は同じ
は巨視的な系における熱伝導率 に相当し,し,$\kappa(L_{\lambda’})$を 2 パラメータ$\kappa$オーダーである。従って,少する。

$L_{X}$$\kappa$$a$$\kappa$$\kappa$$a$が数百程度の系では熱伝導率$\kappa(L_{X})$は$\kappa$$a$から 10% 程度減

182

図 6: 巨視的な系における熱伝導率

Fig. 6 に4$\kappa$を示す。 熱伝導率$\kappa$$\kappa$:$\kappa$は密度に比例する

は密度に比例している。

気液界面

4.1

シミュレーション

次に, 熱浴の温度・平均密度をかえ, 気液共存系を実現する。 系の形状及び境界条件,差分法については, 全て熱伝導率を求めた場合と同様である。 平均密度 0446, 熱浴の温度 $T_{H}=1.32,$ $T_{L}=0.8$ として, シミュレーションを行った。 Fig. 7 に局所温度を示す。温度勾配は $x=20$ 付近を境に大きく異なっており, この境の両側では温度勾配は一定になっていることが分かる。 高温側が気相, 低温側が液相であり, 気液共存系を実現している。

4.2

界面における密度変化

次に, 密度変化について調べる。均質核生成を考え, 密度の気泡がの液相に密度生じたとする。 この時, 密度汎関数法から密度を $\tanh$ の形を持つ位置の関数として求めることが出来る。 これは市 ffusc.-interface または Chan-Hilliard モデル [4] として知られて$D_{g}$$D_{l}$いる。 この界面を用いて核が生成された後の成長の様子を調べる方法を phase-field 法という。 この手法では, 厚みを持った界面を記述でき, さらに成長過程の複雑な形状についても良く再現することが出来る。 しかし, 密度を求めるためにはパラメータとして界面の厚みが残っている。
は気を考える。
及ひ
今, 各相における密度の 1/3 乗に比例する特徴的な長さは液相側の界面の長さである。
相側の,$X_{g}$$X_{l}$$X_{g}$$X_{l}$$D_{g}=a_{g}X_{g}^{-3}$(12)$D_{l}=a_{l}X_{l}^{-3}$(13)は比例係数である。 この長さを用いて, 気相側と液相側の密度の関数をそれぞれ次のように書く。 は気相側と液相側の界面の密度がつながる位置である。
ここで$a_{g}$

及ひ$a_{l}$$x_{0}$$D_{g}(x)=D_{g}+ \beta_{g}tanh(\frac{x-x_{0}}{X_{g}})$

(14)

183

図 7: 粒子数 448,

$L_{\lambda}\cdot=60,$$T_{H}=1.32,$ $T_{L}=0.8$(15)$D_{l}(x)=D_{l}+ \beta_{l}tanh(\frac{x-x_{0}}{X_{l}})$及びは係数である。仮に $a_{g}=a_{l}$ とするととるので及びは各相での密度と を用いて表される。
次に, 密度は界面中の位置において $D_{g}(x)=D_{l}(x)$ かっ
ここで$\mathcal{B}_{g}$$\beta_{l}$$X_{g}$の和は界面の厚み$L$

であ$L$$X_{l}$$X_{g}$$X_{l}$$x_{0}$$\underline{dD}_{d\mathrm{L}^{x)\frac{dD_{l}(x}{dx}}}d.x=$とすると, 密度はそれぞれ次のように書かれる。

ここで$X_{g},$$X_{l}$$D_{\mathit{9}}(x)=D_{\mathit{9}}+ \frac{D_{l}-D_{g}}{X_{g}+X_{l}}X_{g}(1+tanh(\frac{x-x_{0}}{X_{\mathit{9}}}))$(16)$D_{l}(x)=D_{l}- \frac{D_{l}-D_{g}}{X_{g}+X_{l}}X_{l}(1-tanh(\frac{x-x_{0}}{X_{l}}))$(17)

は各相における密度と界面の厚み$L$で表されるので, この関数の形は$L$のみで決まる。

4.3

結果

式 (16),(17) は界面における密度変化を, 密度の 1/3 乗に比例するとして求めた気相側と液相側の界面の厚みて界面の厚み と位置$L$ので,及ひを用いて表すものである。 $D_{g}(x),$ $D_{L}(x)$ はパラメータとしを含んでいるが,は関数を水平方向に移動させるだけである$X_{g},$$x_{0}$$X_{l}$$x_{0}$

この関数の表す曲線は界面の厚みのみをパラメータとして決まる。 Fig. 8 にをシミュレーションから得られたデータに重ねて表示した。

$D_{g}(x)$$D_{l}(x)$シミュレーションから得られたデータと, $D_{g}(x),$ $D_{l}(x)$ は非常に良く一致してぃることが分かる。平均密度, 両端の温度, 系の長さ, 太さを変え同様のシミュレーションを行ったが, いづれの場合においても $D_{g}(x),$ $D_{l}(x)$ はシミュレーション結果ととても良い一致を示した。

5

まとめ
臨界点から 200% 程度高温・高密度の超臨界流状態において熱伝導率の系のサイズに対する依存性を調べた 0 定常状態からはなれた初期状態からシミュレーションを開始する

184

図 8: 粒子数 446, $L_{x}\mathrm{x}L_{Y}\mathrm{x}Lz=60\mathrm{x}4,$ $T_{H}=1.32,$ $T_{L}=0.8$と, 系ははじめに持っていた余分なエネルギーを捨てながら緩和し, 定常状態に至る。定常状態では, 系には定常な熱流が流れており局所温度は温度勾配の方向に対し線形になっていた。熱伝導率は系のサイズに依存し, 系が短くなるに従い小さくなる。 この振る舞い
は $\kappa+a/\sqrt{L_{X}}$ の形に書) 九 これは以前の研究 $[8, 11]$ 及び long-time ta 垣と一致していた。と が同じオーダーであることから,パラメータはがおよそ数百で巨視的な熱伝導率 の値から 10% 程度減少する。 また, 今回のシミュレーション領域において》巨視的な系における熱伝導率 は系の平均密度に比例していた。
次に, 同様の系のにおいて気液共存系を実現し, 界面での密度変化を調べた。 気液界面における密度変化を記述する長さとして, 密度の 1/3 に比例する長さを考え, これを用いて密度汎関数法を基に界面での密度 $D_{g}(x)$ 及び $D_{l}(x)$ を求めた。 これとシミュレーションによる結果を比較したところ, 系の温度・密度また系の大さによらず非常に良く一致した。共存域における温度に対して界面の厚さを調べると, 温度が臨界点に近づくにっれて, 界面の厚みは大きくなっている。 Chain-Hillar モデルなどから求められる界面と異なり, 気相側と液相側においてそれぞれの関数をっくることで, $D_{g}(x),$ $D_{l}(x)$ は非対称界$\kappa$$\kappa(L_{X})$$a$$L_{X}$$\kappa$$\kappa$$X_{g},$$X_{l}$

面を表すことが出来る。
しかし, この関数もパラメータに界面の厚みが残っており, 密度を界面の厚みをパラメータとして含まない形で表すことは今後の課題である。共存域の温度が臨界点に近づいてぃく場合の界面の厚みの変化の様子を調べることも今後の課題である。 さらに, 界面における温度変化も非常に興味深く重要な課題である。

参考文献[1] M. Volmer and H. Z. Flood: Z. Phys. Chem. A190273 (1934).
R. Becker and W. D\"oring: Ann. Phys. 24719 (1935).
J. Frenkel: Kinetic Theory of Liquids (Theodre Steinkopff Dresden, 1939)[2] F. F. Abraham: Homogeneous Nucleation Theory (Academic, New York, 1979)
A. Laaksonen, V. Talanqucr and D. W. Oxtoby: Annu. Rev. Phys. Chem. 46489(1995).

185
J.M.Howe: Philos. Mag. A74761(1996)[3] D. W. Oxtoby: J. Phys. Condensed Matter 47627 (1992).
[4] J. W. Chan and J. E. Hilliard: J. Chem. Phys. 28258 (1958). J. Chem. Phys. 31688 (1959).
[5] D. W. Oxtoby and R. Evans, J. Chem. Phys. 897521 (1988).
X. C. Zeng and D. W. Oxtoby: J. Chem. Phys. 944472 (1991).
F. F. Abraham: J. Chem. Phys. 60246 (1974). Phys. Rep. 5393 (1979)[6] K. Laasonen, S. Wonczak, R. Strey and A. Laaksonen: J. Chem. Phys. 1139741(2000).
, T. Pusztai and P. F. James: J. Chem. Phys. 1176157 (2002).
L.
$\mathrm{G}\mathrm{r}\acute{\mathrm{a}}\mathrm{n}\mathrm{a}^{j}\mathrm{s}\mathrm{y}$[7] J. L. Lebowitz and H. Spohn: J. Stat. Phys. 19633(1978).
G. Casati, J. Ford, F. Vivaldi and W. M. Visscher: Phys. Rev. Lett. 521861(1984).
B. Hu, B. Li and H. Zhao: Phys. Rev. E. 572992(1998).
D. Alonso, R. Artuso, G. Casati and I. Guarneri: Phys. Rev. Lett. 821859(1999).
[8] T. Shimada, T. Murakami, S. Yukawa, K. Saito and N. Ito: J. Phys. Soc. Jpn. 623150 (2000).

[9] B. Li, G. Casati and J. Wang: Phys. Rev. E. 67021204 (2003).
[10] B. J. Alder and T. E. Wanwright: Phys. Rev. A 118 (1970).
[11] T. Murakami, T. Shimada, S. Yukawa and N. Ito: J. Phys. Soc. Jpn. 721049(2003).

[12] J. R. Dorfman and E. G. D. Cohen: Phys. Rev. A 6776(1972). Phys. Rev. A 12292 (1975).

[13] R. Kubo: J. Phys. Soc. Jpn. 12570(1957).
[14] R. Kubo, M. Yokota, and S. Nakajima: J. Phys. Soc. Jpn. 121203(1957).
[15] T. Kinjo and M. Matsumoto: Fluid Phase Equil. 144 (1998) 343.

[16] T. Kinjo and M. Matsumoto: Fluid Phase Equil. 14 (1999) 138[17] H. Okumura and N. Ito: Phys. Rev. E67045301(2003).
[18] T. Hamanaka, R. Yamamoto and A. Onuki: in printing in Phys. Rev. E (2004).

[19] W. G. Hoover: Phys. Rev. A 311695(1985).
S. Nose’:J. Chem. Phys. 81511(1984).

[20] A. Lotfi, J. Vravec and J. Fischer: Mol. Phys. 761319(1992). R. Agrawal and D. A.
Kofke: Mol. Phys. 8543(1995).