Thermal Science & Engineering Vol.15 No.2 (2007)

分子動力学法による微小液滴の∗固体壁面への衝突シミュレーション
松本 充弘†

中澤 伸之†

Molecular Dynamics Simulation of Droplet Impingement on Solid Wall

Mitsuhiro MATSUMOTO†

and Shinji NAKAZAWA†

Abstract
Impingement process of a small droplet on solid wall was investigated with molecular dynamics simulation
technique. A liquid droplet consisting of about 14,000 Lennard-Jones particles was thrown with a given speed onto asmooth wall, and the change of its size and shape was analyzed. After the collision, the droplet spreads on a“hydrophilic,” or strongly interacting, wall, but bounces on a “hydrophobic” wall. The spreading behavior issuppressed by surrounding gas.
Key Words: Droplet impingement, Nano droplet, Nanofluidics, Contact angle, Molecular dynamics
simulation

記 号
T:
V:
r:
x, y, z ::
α:
ε:
φ:
ρ:
σ
添字
LJ:
sc:
w:
w-sc:
１ンタ, スピンコーティングなど様々な分野に利用さ設定温度
液滴の初期速度
粒子間距離
粒子の座標
粒子-壁相互作用パラメタ
エネルギーパラメタ
相互作用ポテンシャル
粒子数密度
粒子サイズパラメタ

れ[1]，最近ではリソグラフィを使わずに半導体回路[K][m/s][m][m]を製造する技術への応用[2]も注目されている．壁面への液滴衝突過程について，これまでにさまざまな可視化実験（例えば[3]）や数値流体計算（例えば[4]）が行われている．最近，衝突後の液滴変形ダイナミ

[J]
[J]
[m-3]
[m]

クスが液滴周囲の雰囲気（蒸気以外のガスの存在）に強く影響されることが Xu らにより報告された[5]．
それによれば，液滴が衝突後に飛沫（splash）を生成するのは，固体壁上を前進する液膜の先端がガスにLennard-Jones 相互作用
ソフトコア相互作用
液滴粒子-固体壁相互作用
ガス粒子-固体壁相互作用

よって遮られるためであると結論づけている．このように，衝突液滴の挙動は多くの外的要因に支配されるため，その制御のために検討すべきことは多い．
現在はピコリットル（直径 10 μm 程度）の液体微粒化技術が確立しているが，将来的にはサブミクロ

緒 言ンからナノメートル領域の微小液滴の生成・制御技液滴をうまく制御して固体壁に衝突させる技術術が必要とされることは間違いないであろう．このは，ミスト冷却，微粒子製造，インクジェットプリ

スケールにおいては，液滴の挙動において表面張力†

受理日: 2006 年 3 月 12 日, 第 43 回日本伝熱シンポジウムより受付，担当エディター: 小原 拓
京都大学大学院工学研究科 機械理工学専攻 (〒606-8501 京都市左京区吉田本町)- 55 -© 2007 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.15 No.2 (2007)が支配的となる．本研究ではこのようなナノメート⎛σ ⎞φsc ( r ) = 4ε ⎜ ⎟ルスケールの微小液滴について，その衝突挙動に影(3)⎝r⎠

響を与える外的要因を検討する出発点として，分子動力学シミュレーションを行った．
２12を仮定し，ガスが液滴の変形を抑制する効果のみに注目して解析を行った．なお，ガス－固体壁間も同計算方法

じ soft core ポテンシャルに由来するポテンシャル
2.1 計算モデル
計算を簡単にするため，液滴を構成する粒子間にφ w− sc ( z ) =は，次の Lennard-Jones（LJ）12-6 相互作用⎣⎢⎝ r ⎠9(4)を仮定した．このため，ガス粒子は固体表面に全く⎡⎛ σ ⎞⎛σ ⎞ ⎤φLJ ( r ) = 4ε ⎢⎜ ⎟ − ⎜ ⎟ ⎥124π⎛σ ⎞εσ 3 ρ ⎜ ⎟90⎝z⎠6(1)⎝ r ⎠ ⎦⎥

吸着しない．実際の系においては，吸着したガス分子が衝突後の液滴の挙動にある程度の影響を与えるを仮定した．
この LJ 粒子からなる球状液滴を作成し，ある設定温度 T で熱平衡状態をつくった後に，一定ことが考えられるが，本研究では平滑な固体壁を仮定しているので，吸着するガス粒子モデルを用いたとしても，その影響を調べることはできない．の並進速度 V を与えて固体壁に衝突させる．本研究

以下では，長さを σ ，エネルギーを ε ，質量を粒では，計算時間を節約するために，固体壁として，LJ 粒子が一様に最密充填密度 ρ で分布している平子質量 m ，時間を τ = σ m ε により無次元化した滑面を仮定した．このとき，固体壁の法線方向を z 軸単位系で諸量を表すが，Lennard-Jones 粒子系の代表とし， z ≤ 0 の半無限領域に広がる固体壁を考える．

としてよく用いられるアルゴン原子の諸量

距離 z ( ≥ 0 ) にある LJ 粒子との間には次の相互作

σ = 0.34 nm
ε = 1 .67 × 10 −21 J
m = 6.63 × 10 −26 kg
τ = 2.14 ps

用がはたらくことがわかる：
0+∞+∞−∞−∞−∞φw ( z ) = αρ ∫ dz ' ∫ dx ' ∫ dy ' φ LJ

( x' + y ' +( z − z' ) )22⎡ 1 ⎛ σ ⎞9 1 ⎛ σ ⎞3 ⎤= 4παεσ 3 ρ ⎢ ⎜ ⎟ − ⎜ ⎟ ⎥⎢⎣ 90 ⎝ z ⎠ 12 ⎝ z ⎠ ⎥⎦2( 2)を用いて換算した値も付記する．

ここで， ε と σ はそれぞれ LJ 相互作用のエネルギ

2.2 計算条件

ーパラメタとサイズパラメタである．また，液滴と各粒子の運動方程式を leap frog 法により数値積分

固体壁の相互作用の強さ（親水性，疎水性）を調節する無次元パラメタ α を新しく導入した．

することで，粒子の軌跡を追跡した．時間刻みは0.005 (アルゴン換算で約 10 fs) とした．約 14,000

さらに，ガス雰囲気の影響を検討するため，液滴個の LJ 粒子を球状に配置して液滴を作成し，サイズを構成する粒子との相互作用が弱く，また設定温度が 100σ × 100σ × 100σ の計算セル中央部に置いて温において「気体」としてふるまう粒子を液滴周囲に度 T で平衡化した．平衡化の間に若干の粒子が蒸発導入する．この目的のためには，気液臨界温度が十するが，結果として得られた液滴直径は約 30σ (約分に低い（すなわち液滴粒子に比べて小さい ε を持10 nm)であった．この静止液滴を構成する全粒子につ）LJ 粒子を利用するのが 1 つの方法である．しか同一の初期並進速度 V を与え，計算セルの一端し，これは新たなパラメタ（ガス粒子の ε ）を導入することになる．また，液滴へのガスの吸着や溶解( z = 0 ) に設けた固体壁に垂直に衝突させた．ただが起きることになるが，実際系（例えば空気中の水し，この計算中には温度制御は行わない．また，x, y滴）を考えると，吸着や溶解が液滴の衝突挙動に与方向には周期境界条件を仮定した．制御パラメタのえる影響は副次的なものと思われる．このため，本すべての組み合わせ ( 2 × 8 × 3 × 2 = 96 通り) につい研究では，ガス粒子間ならびにガス－液滴粒子間の相て各々200,000 ステップ (約 2 ns)の計算を行い，1000互作用として反発力のみをもつ soft core ポテンシャル

ステップごとに蓄積した粒子データを解析した．

- 56 -© 2007 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.15 No.2 (2007)

Fig. 1 Density profile (contour map) of the droplet at rest on solid wall; (left) T=0.65, (right) T=0.75.

・平衡化温度 T ：LJ 流体の三重点温度近傍である

実験[5]などに比べるとかなり高密度であるが，ガス0.65 (アルゴンでは 78 K 相当)と，それより少し高温粒子の平均自由行程が液滴サイズ以下の場合にのみの 0.75 (90 K 相当)の 2 種類とした．

ガスの影響が顕著にあらわれると考えたためである．

・初期並進速度 V ：ここでは，0.2, 0.4, 0.6, 0.8, 1.0,

この条件下でのガス粒子の平均自由行程は，約 3.5 σ

1.25, 1.5, 2.0 の 8 通り (32～320 m/s 相当)を選んだ．

(液滴直径の 1/8 程度) で，液滴との相互作用は十分・固体壁との相互作用パラメタ α ：1.0, 0.5, 0.1 の３にあると思われる．

種類とした．式(2)からわかるように， α = 1.0 は，固体壁が液滴と同じ“材質”で構成されていること

３

シミュレーション結果と考察に対応するので，ここでは“親水性 hydrophilic”と呼ぶことにする．一方，α = 0.1 は，固体壁と液滴の3.1 固体壁の“濡れ性”相互作用が非常に小さいことに対応し，“疎水性まず，固体壁の相互作用パラメタ α の影響を調べhydrophobic”と称することにする．その中間 α = 0.5

る基礎データとして，固体壁上に静置した液滴の数は“semi-hydrophobic”と表記する．

密度分布を Fig. 1 に示す．ただし，液滴の軸対称性・ガス雰囲気：上述のようにして構成した LJ 粒子だを仮定して解析を行った．

けの系と，液滴部分を除く計算セル全体に液滴と同マクロスコピックな濡れ性は“接触角”によってじ初期温度の 64,000 個のガス (soft core) 粒子を均記述することが一般的であるが，分子スケール／ナ一に分布させた系について比較する．これは 20 atm

ノスケールでは，気液界面の揺らぎや固液界面での程度のガス圧に対応する．これは，報告されている液体の構造化のために接触角を一意的かつ精確に求- 57 -© 2007 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.15 No.2 (2007)

Table 1 Evaluated contact angle of the droplet on
solid wall.

義した．Fig. 3 に液滴サイズと接触面積の時間変化を，典型的な条件について示した．これらの解析結

αT=0.65

T=0.750.1150°±10°

150°±10°0.5

100°±5°90°±5°

は大きく，飛散（蒸発）による液滴サイズの減少は1.0

40°±5°

20°±5°

著しい．また，当然ながら液滴は壁面上に大きく広がる．特に，前述のように V = 1.5, 2.0 の場合に，真果から，以下のことがわかる．
3.3.1

衝突速度依存性

どの条件においても，衝突速度 V が大きいほど衝突時に液滴が受ける衝撃

空中では，多くの条件において液滴は衝突後にいくつかのクラスターに分裂する様子が見られた．

めることは難しい[6]．ここでは仮に，数密度が 0.7

3.3.2 固体壁の性質による違い Fig. 3 (a), (b), (c)の等密度面と固体壁とのなす角を接触角とみなすこを比較すると明らかであるが，親水壁（a）上では液とにして，Table 1 のように接触角を読み取った．こ滴は壁に沿って液膜状に大きく広がるのに対し，疎の結果から， α = 0.5 の固体壁は予想通り中立的な水性が強い（c）場合には，広がりかけた液膜は表面性質を持つと言えるだろう．

張力によって再び収縮した後，リバウンドする．これは，Fig. 2 のスナップショットからもうかがえる．

3.2 衝突後の液滴の観察

慣性力に比べて表面張力が支配的である微小液滴の液滴の形状変化のスナップショットを，典型的な3 つの条件について，Fig. 2 に示した．多くの場合，特徴がよくあらわれている． α = 0.5 (b)は，やはり100,000 ステップ (約 1 ns 相当) 程度でほぼ定常とな中間的な挙動を示す．

る．見やすくするために，各粒子の持つ運動エネル

3.3.3

ギーに応じて色分けして示している．本研究では，ことにより，衝突後の液滴挙動は，定性的には初期粒子は並進自由度のみをもつ質点であり，また粒子の液滴温度には強くは依存しないことがわかる．た間距離のみのポテンシャル関数 (1) を仮定していだし，スナップショットなどを詳細に検討すると，液滴温度依存性

Fig. 3(a)と(d)を比較するる．このため，液滴から脱離した粒子はポテンシャ

高温では液滴の表面張力が小さい（ T = 0.75 と 0.65

ルエネルギー相当分だけエネルギーを失い，運動エでの Lennard-Jones 液体の表面張力の比は 0.8 程度

ネルギーが低下していることが色の変化からわかる．

[7]）ために，液滴がやや分裂しやすいことがわかる．

これは，熱平衡状態での蒸気（液滴の粒子と蒸気のこのように，温度の影響が小さいのは，本研究の計粒子は温度が等しい）とは異なる点である．

算条件においては，衝突により並進運動エネルギー一連のスナップショットから，ガスがない場合にから転換される熱エネルギーが大きいためであろう．

は，衝突の衝撃により，液滴から多くの LJ 粒子が飛たとえば，初期並進速度が V = 1.0 の場合に，仮にそ散（あるいは蒸発）していることがわかる．この図のすべての並進運動が衝突によりランダムな運動にには示していないが，液滴の衝突速度が大きく（ V = 1.5, 2.0 ），かつ濡れにくい（ α = 0.1, 0.5 ）場なったと仮定すると，温度は 0.33 程度上昇すること

合には，数個から 100 個程度の粒子で構成される数

この計算では並進自由度しか持たない LJ 粒子系を多くのクラスターが液滴から直接に脱離する様子も扱っているが，多原子分子のように回転・振動など見られた．他方，ガス雰囲気中では，このような飛の他の自由度をもつ場合には，相対的に液滴温度の散や液滴の変形は抑制される．

影響は大きくなると予想される．になり，初期温度の違い（0.1）を圧倒する．ただし，3.3.4 雰囲気による違い 親水壁の場合 [Fig. 3(a)3.3 データ解析と(e)] には，真空中（蒸気のみ）での液滴衝突に比液滴の変形挙動を，(1) 液滴サイズ（構成粒子数），べて，ガスの存在下では液滴サイズの減少が小さく，(2) 液滴形状（高さと投影面積），および (3) 接触面液滴の広がりも小さい．これは，スナップショット積，に注目して定量的に解析した．ここでは，あるFig. 2 からも明らかである．ガスが液滴の形状変化，閾値（ − 3ε とした）より低いポテンシャルエネルギ

特に壁に沿って広がる液膜のフロントラインの進行ーをもつ LJ 粒子を“液滴を構成している粒子”と定を妨げているためであり，実験で報告されている挙- 58 -© 2007 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.15 No.2 (2007)###Fig. 2 Sequential snapshots of the droplet collision process; (left) on hydrophilic wall in vacuum, (middle) on
hydrophobic wall in vacuum, (right) on hydrophilic wall in gas, where gas particles are not shown. Theimpinging speed is V=1.0 and the initial temperature is T=0.65. The color of each particle represents its
kinetic energy.

- 59 -© 2007 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.15 No.2 (2007)

Fig. 3 Examples of droplet dynamics after the impingement; change of the droplet size (number of constituent LJparticles) and the contact area during the droplet collision.

動[5]に類似している．ただし，本計算の液滴は極め3.3.5て小さく（直径 10 nm 程度），表面張力が支配的で数密度分布から液滴の形状変化を追跡した Fig. 4あるため，実験（液滴直径は 3 mm 程度）で報告さにその一部を示す．接触が始まった後，液滴の中心れているような飛沫生成は見られない．

部が固体壁に近づくにつれて液滴周辺部を押し出す

液滴形状の時間変化

軸対称性を仮定した

他方，疎水壁の場合 [Fig. 3(c)と(f)] には，ガスがことにより液膜が広がっていく．前述のように，濡液滴の蒸発を抑える効果は見られるが，液滴の広がれ性が低く衝突速度が大きい場合 [Fig. 4 (d)] には，り（接触面積）には大きな差はない．これは，もと液膜先端がちぎれてクラスター化する．しかし，１

もと疎水壁との衝突による液滴の広がりが小さいた

粒子ずつの飛散（蒸発）は液膜の場所によらず至るめであろう．

ところで起きていることが，粒子の詳細な軌道解析- 60 -© 2007 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.15 No.2 (2007)(a) α = 1.0, V = 1.0(b) α = 1.0, V = 2.0(c) α = 0.1, V = 1.0(d) α = 0.1, V = 2.0

Fig. 4 Change of the droplet shape (density contour map) in vacuum at T = 0.65 .
によりわかった．

本研究の一部は，京都大学 21 世紀 COE プログラム「動的機能機械システムの数理モデルと設計論―

４

おわりに複雑系の科学による機械工学の新たな展開―」の援ナノメートルスケールの微小液滴が固体表面に衝助を受けて実施された．

突する際の動的挙動を分子動力学法により調べ，衝突速度・液滴温度・固体壁の性質・雰囲気の影響を参考文献検討した．計算の都合により，並進自由度のみをもつ Lennard-Jones 相互作用粒子を使用し，一様かつ滑らかな固体壁を仮定するなど，非常にラフなモデル系を用いたため，現実系との直接的な比較はまだ困難であるが，衝突後の微小液滴の変形ダイナミクスが固体壁の性質（濡れ性）や雰囲気（ガス圧）に大きく依存するという知見は，微小液滴の利用技術を開発する上で重要だと思われる．今後，構造を持った壁面への衝突解析や衝突角度への依存性解明など，多方面への展開が可能である．
衝突液滴の変形は連続体計算[4]により解析することが一般に行われているが，蒸発や分裂を伴う微小液滴の複雑な挙動の解析には分子動力学法が有用であろう．

- 61 -[1] Nasr, G.G., Yule, A.J. and Bendig, L., Industrial
Sprays and Atomization: Design, Analysis, andApplications (2002) Springe-Verlag.
[2] 下田達也, “インクジェットで Si-TFT を形成”,
日経マイクロデバイス, July-253 (2006) 55-61.
[3] Mundo, C., Sommerfeld, M. and Tropea, C., Int. J.
Multiphase Flow, 21 (1995) 151-173.
[4] Kamnis, S. and Gu, S., “Numerical modelling ofdroplet impingement”, J. Phys. D: Appl. Phys., 38(2005) 3664-3673.
[5] Xu, L., Zhang, W., and Nagel, S.R., “Dropsplashing on a dry smooth surface”, Phys. Rev. Lett.,
94 (2005) 184505.
[6] Maruyama, S. et al., “Liquid droplet in contact witha solid surface”, Microscale Thermophys. Eng., 2(1998) 49-62.
[7] Toxvaerd, S., “Perturbation theory for nonuniform
fluids: surface tension”, J. Chem. Phys., 53 (1971)
3116-3120.
© 2007 The Heat Transfer Society of Japan