分子動力学計算による気泡生成過程の研究
鈴木将 1 ，渡辺宙志 2
1

名古屋大学 COE 計算科学フロンティア，2 名古屋大学情報科学研究科
suzuki@phys.cs.is.nagoya-u.ac.jp

概要: 沸騰現象の非経験的な取り扱いを目指し、分子動力学法による全粒子計算を行った。全粒子計算を行うことで界面の形成と相転移は自発的に行われ、相間の摩擦や熱伝導率といった物性値も外から与えずに数値計算を行うことができる。我々は臨界沸騰近傍において系を膨張、収縮させることにより、気泡の多重臨界核の発生とその成長、及び消滅ダイナミクスを調べた。これらの結果を化学ポテンシャルに基づいた議論により評価し、気泡核生成時の界面張力の影響、及び臨界核半径について議論した。粒子を用いた直接計算による混相流研究の展望についても述べる。

1

はじめに

2

核生成とその成長ダイナミクス

沸騰現象は発電所などの熱機関において中心的な役割

沸騰をミクロに見れば、それまで液相にあった系の一を果たすが、強い非平衡現象であること、相転移を含む部に突然気泡が生じ、その気泡が成長することでおきる。

移動境界問題であることから扱いがきわめて難しい。特ある大きさ未満の気泡は界面張力により消滅し、それよに軽水炉などで問題となる気液二相流は、熱伝導や流動り大きい場合は成長する。この半径は臨界核半径と呼ばなどの輸送現象に加え沸騰という相転移も複雑にかられ、古典核生成論と呼ばれる理論により取り扱われる。

みあっており、スケールの異なる現象の同居が研究の妨

まずマクロな液相に半径 r の微小気泡が静的に存在するげとなってきた。このようなマルチスケールな問題につとき、気泡内部および外部の圧力 pg , pl は界面張力を γ

いては、現象をそのスケールにより分類し、計算するスとして力学的平衡条件からケールよりも小さいスケールについては粗視化された構

pg − pl成式を用いることで研究が行われてきた。気液二相流の2γr=(1)場合は、ミクロな現象である相転移と、マクロな現象でを満たす。また気液共存状態近傍の気相側及び液相側のある流動を同時に扱うため、相間の熱輸送、質量輸送、 化学ポテンシャルをそれぞれ µ (p ) および µ (p ) とす
ggll

摩擦を表現する構成式を用意することで計算が行われてると熱力学平衡条件より次式が成り立たなくてはならいる [1]。しかし、気液界面における熱やエネルギー輸ない。

送、質量変化の構成式は局所熱平衡及び線形非平衡を仮µg (pg ) =定して導出されたものがほとんどであり、沸騰現象といµl (pl ).

(2)

う強い非平衡系においてその妥当性に疑問が残る。
そこで、流体を構成する全ての粒子の運動を分子動力学法により計算することで沸騰現象を再現する。粒子法
においては界面は自発的に形成され、かつ気液界面における摩擦や熱輸送といった実験では調べづらいパラメータを入力として与える必要も無い。一般にナビエ・ストー
クス方程式などの支配方程式を差分化して時間発展を追う格子法に比べ、分子動力学法などの粒子法は計算量が大きくなるが、近年の計算能力の増大により、気液界面や気泡生成を直接全粒子計算により研究することが可能となった [2, 3, 4]。
本稿では分子動力学法による沸騰現象の全粒子計算のこれまでの成果と、その可能性について報告する。

ただし、化学ポテンシャルは気液共存圧力 p0 においてµg (p0 ) = µl (p0 ) となるべき量である。上記の力学的および熱力学的条件の釣り合いから臨界核半径 rc が決まる。密度揺らぎなど、なんらかの原因で rc を上回る半径の気泡が自発形成された場合、この半径は不安定である
ために気泡は暴走的に成長し、マクロな沸騰現象が起こると解釈される。しかし、古典核生成論が記述するのはあくまで液相中に気泡が静的に存在する場合であり、気泡の生成の瞬間、およびその後どのように成長していくかといったダイナミクスはまったく分かっていない。そこで、我々は分子動力学法により気泡生成および消滅を再現し、そのダイナミクスについて検討した。粒子間の相互作用としては Lennard-Jones ポテンシャルを採用した。粒子数は 64 万個とし、温度は Gaussian Thermostat法で制御した。
まず、系が単一の液相となるように温度及び密度を制御しておき、その後系を膨張させることで気泡を生成させた。系の密度が臨界密度を下回ると系は不安定化し、その後に同時多発的に気泡核が生成される (図 1)。この後、単一の気泡のみが成長し、他の気泡は収縮して消滅する Ostwald Ripening のような過程が観察された。これは、成長する気泡が局所的な密度と圧力を介してお互いに相互作用し、一番大きな気泡のみが選択的に成長するためと考えられる。
また、臨界核半径程度の非常に小さいスケールにおいて (1) 式及び (2) 式のいわばマクロな表現をそのまま適用することの妥当性を検討するため、気-液共存状態の温度・密度において単一の球状気泡を制御的に生成した上で、安定化した気泡半径を評価した。条件式 (1)、(2)より決まる pg , pl を用いて、与えられた系の平均密度 ρに対する気泡半径 r は

N=4 3πr ρg (pg ) +3()N4− πr3 ρl (pl )ρ3(3)と書くことができる。ただし、N は全粒子数、ρg および ρl はそれぞれ気相と液相の密度である。ここから評価される半径は広い領域でシミュレーションと一致したものの、特に半径の小さい領域では有意なずれが存在した。これは表面張力 γ のサイズ依存性や気泡形状の球からのずれを示唆していると考えられ、より詳細な計算を進めている。

3まとめと今後の展望図 1: 分子動力学法により自発形成された気泡核。上図は粒子の位置を表示したもの。局所的に密度が薄い気相状態と密度が高い液相状態が共存していることがわかる。下図は、局所的に気相密度にある領域を仮想粒子法により表示したもの。
クラスタリング判定により、異なる気泡を色分けして表示している。大小の気泡がいくつか生成されていることがわかる。

Lennard-Jones 粒子を用いた全粒子計算により、気泡 の発展により、マルチスケールな問題を直接計算する全の生成及び消滅過程を再現した。沸騰初期における気泡 粒子計算の果たす役割がますます大きくなるものと期待核の同時出現及び選択成長を確認し、また、制御された する。
本研究は名古屋大学 COE 計算科学フロンティア、及気泡の収縮過程を気相及び液相の化学ポテンシャルを用いて解析し、気泡成長のダイナミクスを明らかにした。 び科学研究費補助金 若手研究 (B) 課題番号 19740235 の本研究におおいて扱われた粒子数は 106 個程度であり、 助成を受けて行われた。
自発形成する気泡は数個程度であった。これにより沸騰の初期段階における気泡生成は再現できたが、その後の参考文献

多数の気泡が競合的に成長するマクロな沸騰現象を再現[1] 秋山守、有冨正憲、監修 「新しい気液二相流数値解析」、
コロナ社、2005 年.

するには、長さスケールで 10 倍、粒子数にして 103 倍程度の系のシミュレーションが必要となる。さらに、気泡が流れとカップルした気泡流の計算を行うためには、さらにもう一桁上、すなわち粒子数にして 1012 個の計算

[2] F. Ogushi, S.Yukawa, and N.Ito, ”Asymmetric Structure of Gas-Liquid Interface” J. Phys. Soc. Jpn. vol. 75
No.7 (2006) 073001.

が必要となる。これは以前では不可能な計算規模であっ [3] H. Okumura and N. Ito, ”Nonequilibrium molecular dyたが、近年のスーパーコンピュータ、たとえば京速計算機にとっては十分に到達可能な規模となっている。
これまでの流体ダイナミクスに関する計算手法は、いわばマクロからのアプローチとしてナビエ・ストークス方程式などの流体力学方程式を差分化して解く手法が主流であった。これに対し、粒子による計算はミクロからの構成論的なアプローチといえよう。計算能力の今後

namics simulations of a bubble”, Phys. Rev. E vol. 67(2003) 045301(R).

[4] S. Tsuda, S. Takagi, and Y. Matsumoto, ”A study onthe growth of cavitation bubble nuclei using large-scale
molecular dynamics simulations”, Fluid Dynamics Research, vol. 40 (2008) 606.