7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Molecular Dynamics Study of Vaporliquid Equilibrium Condition for Nanodroplet
H. Yaguchi *9
T. Yano F, M. Watanabe and S. Fujikawa *
Division of Mechanical and Space Engineering, Hokkaido University, Sapporo 0608628, Japan
t Department of Mechanical Engineering, Osaka University, Suita 5650871, Japan
Shyaguchi~,mechme.eng. hokudai.ac~jp
Keywords: molecular dynamics, thermodynamics, equilibrium condition, argon, vapor, droplet, chemical potential
Abstract
Molecular dynamics (MD) simulations of equilibrium system of single argon nanodroplet and its surrounding argon
vapor are carried out to address a fundamental issue whether the thermodynamic description is applicable to the
nanoscale inhomogeneous system. The profiles of both density and pressure in the transition layer are computed.
Numerical results show that the chemical potentials of liquid and vapor phases are not equal when a droplet is so
small that the number of molecules consisting the transition layer may be comparable to that in the droplet. The phase
equilibrium condition in thermodynamics fails in such a case, and therefore, all thermodynamics relations based on
this condition, such as Kelvin equation connecting vapor pressure to droplet radius, do not hold.
Introduction
/Inter~face
Liquid film
Interface
Droplet
A singlecomponent vaporliquid twophase system
consisting of a droplet and its surrounding vapor in an
equilibrium state is a fundamental system in various
fields of science and technology, such as ultrafine par
ticle formation (Adachi et al. 2003), cloud formation
(Chuang 2006), uptake of chemical species into aerosols
(Morita & Garrett 2008), and semiconductor cleaning
(Watanabe et al. 2009). According to thermodynamics,
the equilibrium condition of the vaporliquid system is
described as
Figure 1: (a) Nanodroplet of Rs=3.91 nm and (b) pla
nar liquid film in vaporliquid equilibrium state at 85 K.
radius Rs such that Eq. (1) holds.
Method of Molecular Dynamics Simulations
Figures 1 shows typical snapshots of the vaporliquid
equilibrium systems at 85 K for (a) a nanodroplet and
(b) a planar liquid film. The periodic boundary condi
tions are imposed for both cases in all three directions of
the simulation cells. The simulations of all nanodroplets
are performed in cubic simulation cells with dimensions
L x Lx L:; L for each simulation is given in Table 1.
where <, and <, are chemical potentials of bulk vapor
and liquid, respectively, and the temperature should be
uniform in the whole system. Equation (1) is of fun
damental importance in thermodynamics of phase equi
librium; in fact, all thermodynamics relations in the
phase equilibrium, such as ClausiusClapeyron relation,
Kelvin equation, and so on, can be derived from Eq. (1).
However, thermodynamics is a macroscopic or contin
uum theory for equilibrium systems. The continuum
theory implicitly assumes that the system concerned
comprises sufficiently large number of molecules. Our
aim of this study is to investigate the validity of thermo
dynamics for an argon nanodroplet in equilibrium with
argon vapor by molecular dynamics (MD) simulations.
More precisely, we obtain the lower bound of droplet
No. T IN L Rs e pe p, Pe p, T
p,RT
(K) (nm) (kg/m3) (MPa) (N/m)
1 84.8 8000 21.0 3.91 20.6 1412 6.21 6.6 0.106 0.0127 0.96
2 84.9 4000 16.5 2.99 18.9 1416 6.77 8.4 0.115 0.0124 0.96
3 85.0 2000 13.0 2.27 17.0 1420 7.55 10.6 0.126 0.0119 0.95
4 85.1 1000 10.5 1.70 14.9 1427 8.58 12.8 0.144 0.0107 0.95
5 85.0 730 10.0 1.43 13.9 1429 9.19 14.8 0.154 0.0104 0.95
6 85.3 515 9.5 1.16 12.0 1428 10.72 16.2 0.177 0.0093 0.93
7 85.0 450 11.0 0.84 10.3 1430 12.46 17.7 0.205 0.0073 0.93
8 84.9 515 12.0 0.74 9.5 1421 13.42 18.8 0.219 0.0069 0.92
9 84.9 500 12.0 0.63 8.8 1414 14.55 16.7 0.234 0.0052 0.91
10 84.0 515 12.0 0.51 8.0 1407 16.07 15.9 0.254 0.0040 0.90
11 85.4 700 13.0 0.41 6.8 1371 18.83 12.3 0.298 0.0025 0.90
12 84.6 515 11.5 0.27 6.6 1379 19.33 11.9 0.303 0.0016 0.89
T Zs e pe p p, TopR
(K) (nm) (kg/m3) (MPa) (N/m)
84.9 3.44 26.9 1394 (1410) 4.76 (4.59) 0.083 (0.079) 0.0136 (0.0131) 0.98
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: The parameters and numerical results of vaporliquid equilibrium states of an argon nanodroplet and its
vapor.
Table 2: The numerical results of vaporliquid equilibrium state of an argon planar liquid film and its vapor. The
values in parentheses are experimental ones at 85 K (Tegeleret al. 1999; Miller et al. 1976).
The simulation of planar liquid film is in rectangular cell
with L, x L, x L,; L, L,=5 nm and L,=30 nm. In
preparatory computations, we have confirmed that sta
tistical properties such as densities, mass fluxes, etc. are
not affected by the size of the simulation cells, the liner
dimension of which is several multiples of droplet ra
dius.
The number of all molecules N, the total energy E in
the system, and the volume V of simulation cell are kept
constant during the present MD simulation (NVE sim
ulation). The following 126 LennardJones potential is
used as intermolecular potential;
( 2
where rij is a distance between the centers of a molecule
i and a molecule j, E and a are potential parameters. As
is well known, different values of e and a give prop
erties of different monoatomic materials. In this study,
we choose argon to obtain physical perspectives from
the results. They are, respectively, given aS E/k=119.8 K
(k is Boltzmann constant) and a=0.3405 nm for argon
(Michels et al. 1949).
Equations of molecular motion are integrated by the
leapfrog method with the time step of 5 fs. First,
preparatory computations are performed with the tem
perature control by using the velocity scaling method
for more than 100 ns to realize vaporliquid equilibrium
at 85 K, and then, without the temperature control for
200 ns to eliminate influence of velocity scaling. The
criterion for the realization of the equilibrium states is
that the temperature and droplet size become steady and
the velocity distribution function becomes Maxwellian
everywhere in the quiescent system. After the care
ful confirmation of the equilibrium state, these config
urations in the system are adapted as the initial states:
thereafter the main computations are performed during
201.6 ns without the temperature control.
The instantaneous shapes of smaller nanodroplets are
far from a spheres, because the shapes are fluctuating
with the time. Nevertheless, it is confirmed that time
averaged spatialdistributions of macroscopic quantities
such as density, mass fluxes, velocity, and etc., are spher
ically symmetric. These quantities are averaged with re
spect to I and as well as the time t and are expressed as
functions of the radial distance r alone in the spherical
coordinate system (r, 8, ~) in which the origin coincides
Isnn. . . . . . . .
15000
500 o
( a) oo
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
loC o
( a) o
_L1
S100
,x50
r (nm)
r (nm)
Figure 2: Profiles of (a) density and (b) pressure in the
nanodropletvapor system at 85 K: Rs = 2.99 nm.
with the gravity center of nanodroplet.
In the simulation of the argon planar liquid film and
its vapor system, the computational and data processing
methods are almost the same as those in the preceding
study by Ishiyama et al. (21li14. The total number
of molecules is N=4000 and the temperature is set at
85 K. The preparatory computations to realize the equi
librium state are performed as well as in the simula
tions of nanodroplets. Cartesian coordinates (x, y, z),
the origin of which coincides with the gravity center
of the liquid film, are adopted. The physical quanti
ties are averaged with respect to x: and y as well as the
time t. There are two interfaces in the cell as shown in
Fig. 1(b); hence physical quantity profiles from liquid to
vapor are expressed as functions of the absolute value of
zcoordinate.
In this study, we calculate macroscopic physical quan
tities, such as the pressure and surface tension, from
the molecular motion on the basis of mechanical argu
ment. The detailed information of definition and cal
culation method for these macroscopic physical quanti
ties is available on the textbooks (Rowlinson & Widom
1982; Allen & Tildesley 1987) and our previous paper
(Yaguchi et al. 2010).
Figure 3: Profiles of (a) density and (b) pressure in the
nanodropletvapor system at 85 K: Rs = 1.43 nm.
Results and Discussion
Figures 2(a) and 3(a) are the density profiles in the r
direction for nanodroplets, and Fig. 4(a) is the den
sity profiles as functions of the absolute value of z
coordinate for planar liquid film. The high density re
gions on the lefthand side of profiles are the bulk liquid
phase whose density is pe, whilst the low density regions
are the bulk vapor phase whose density is p,. The region
with abrupt changes in density and pressure is the tran
sition layer.
Figures 2(b) and 3(b) show the profiles of the pres
sures pr, po, and pe, which are acting on the surfaces
normal to r, 8, and directions for nanodroplets. The
pressure ps and pp are identical to the pressure pr in
the bulk vapor and liquid phase, sufficiently outside the
transition layer. Thus, we define the pressure in the bulk
liquid as the liquid pressure pe and the pressure in the
bulk vapor as the vapor pressure p,, respectively. The
liquid pressure pe is defined as the pressure at r = 0, cal
culated by using the cubic spline interpolation with zero
gradient at r = 0, shown as solid curves in Fig. 2(b) and
3(b). As shown in Fig. 2(b) and 3(b), the liquid pressure
pe is larger than the vapor pressure p, due to the surface
tension ys.
141111. . . . . . . .
1 A o p,
20 
(b)
 z (nm)
Figure 4: Profiles of (a) density and (b) pressure in the
liquid filmvapor system at 85 K.
We adopt the radius of surface of tension as the
droplet radius Rs. The surface of tension is defined as
the surface on which the surface tension acts. The ra
dius Rs is determined simultaneously with the surface
tension ys from the pressure profiles in the transition
layer. The detailed calculation method and results of ys
and Rs are discussed in our previous paper (Yaguchi et
al. 2010).
Figure 4(b) shows the pressures p, r and pe,which
are acting on the surfaces normal to .r, y, and z direc
tions for planar liquid film. The pressure p, and r are
identical to the pressure p, in the bulk vapor and liquid
phase. The saturated vapor pressure p, is defined as the
pressure in the bulk vapor phase. Except for the transi
tion layer, the pressure p, is distributed with uniformity
in bulk vapor and liquid because contribution of surface
tension y, vanish in the case of planar surface. The sur
face of tension Zs is also determined with y, for liquid
film as well as nanodroplet.
Table 1 shows these numerical results of vaporliquid
equilibrium states of an argon nanodroplet and its va
por. The qualitative tendency agrees with the results of
preceding studies (Thompson et al. 1984: Nijmeijer et
al. 1992: Vrabec et al 2006). For reference, the mean
free path I = m/( Jyro pe,) of vapor molecules and the
compressibility factor, p /(pe,RT), are shown, where at
is mass of one molecule and R is gas constant. Table 2
also shows result of planar liquid film. The values of
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
MD in Table 2 are in good agreement with experimental
ones shown in parentheses (Tegeleret al. 1999; Miller
et al. 1976).
Now, we discuss the chemical potential in bulk vapor
and liquid phase to verify the thermodynamics equilib
rium condition, Eq. (1). Under a constant temperature,
i.e., clT 0, the differential of chemical potential is
10. o
(>o 9
S100
.r 50
t
where script / is v for vapor phase or I for liquid phase.
Integration of Eq. (3) with respect to Rs from oc to Rs
results in
1 (R,)
dliPi
where p' is the chemical potential both in vapor and
liquid for planar surface. Making use of Eq. (4), we can
evaluate an excess chemical potential of vapor p p'
and that of liquid p, p' at a given temperature T from
results ofMD simulations.
Figure 5(a) shows the relation between the specific
volume 1/pe, and pe, in the vapor phase, where the num
bers shown near the symbols correspond to these of the
data in Table 1. The solid curve in Fig. 5(a) is fitting
function, 1/pe, = /pe, + B, where 24 is a fitting pa
rameter and B is determined so that the curve passes
through the point of planar surface at pe,=0.083 MPa.
The integration of the righthand side of Eq. (4) is cal
culated by using this fitting function to obtain pe, p'
Figure 5(b) is 1/p, versus pe in the liquid phase. The
specific volume 1/p, is found to hardly depend on the
pressure pe: hence we introduce approximate function,
1/p, 1/, :~ where , is the specific volume in the
case of planar surface at pe=0.083 MPa. Substituting this
approximate function into the righthand side of Eq. (4),
we can evaluate p, p'.
Figure 6 shows the comparison betiveen pe, p' and
in. l' The vapor chemical potential pe, at 85 K mono
tonically increases with the decrease in nanodroplet ra
dius Rs, while the liquid chemical potential p, has a
maximum around radius of about 0.8 nm. The two
chemical potentials agree for droplets with radii larger
than 1.5 nm, and they disagree for nanodroplets with
radii smaller than 1.5 nm. That is, at 85 K, the phase
equilibrium condition in thermodynamics, Eq. (1), fails
for Rs less than 1.5 nm, and therefore, all thermody
namics relations based on Eq. (1), such as Kelvin equa
tion, do not hold. The number of molecules consisting
the transition laver may be comparable to that in such
a small nanodroplet. We here emphasize that, even for
nanodroplets with radius smaller than 1.5 nm, the equi
librium state in the sense that the velocity distribution
0 5 10 15
Liquid pressure pe (MPa)
i( u . .
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.0012
0.0010
0.0008
0.0006
0.0004
0.0002
Planar surface (Liquid film)
1 2 3 45 7 8
(b)
Vapor pressure p, (MPa)
Figure 5: Relations between specific volume and pressure in (a) vapor phase and (b) liquid phase.
function of molecules in small volume is Maxwellian,
is achieved everywhere in the system, and the mechani
cal equilibrium, as Laplace equation, also holds.
85 K g yv po
hO~e~
P
B
25000
20000
15000
10000
5000
Conclusions
Molecular dynamics (MD) simulations of equilibrium
system of single argon nanodroplets and its surround
ing argon vapor have been performed to investigate
the validity of the equilibrium condition in thermody
namics, Eq. (1), in the nanoscale inhomogeneous sys
tem. The chemical potentials of bulk vapor and liquid
phase disagree for nanodroplets with radius smaller than
1.5 nm, where the number of molecules consisting the
transition layer may be comparable to that in the nan
odroplet. Therefore, all thermodynamics relations based
on Eq. (1), such as the Kelvin equation, do not hold.
Acknowledgements
This work is supported by Japan Society for the Pro
motion of Science, GrantinAid for Scientific Research
(A)(No. 21246031) and (B)(No. 19360077).
References
Adachi, M., Tsukui, S., and Okuyama, K., Nanoparticle
Formation Mechanism in CVD Reactor with Ionization
of Source Vapor, J. Nanoparticle Res., Vol. 5, pp. 3137,
2003.
Chuang, P. Y., Sensitivity of Cloud Condensation Nuclei
Activation Processes to Kinetic Parameters, J. Geophys
ical Res., Vol. 111, DO9201, 2006.
0 1 2 3
Droplet radius Rs (nm)
Figure 6: Verification of chemical potential equilibrium
between vapor and liquid.
Morita, A. and Garrett, B. C., Molecular Theory of Mass
Transfer Kinetics and Dynamics at GasWater Interface,
Fluid Dyn. Res., Vol. 411\'8), pp. 459473 (2008).
Watanabe, M., Sanada, T., Hayashida, A., and Isago,
Y., Cleaning Technique Using HighSpeed SteamWater
Mixed Spray, Solid State Phenomena, Vol. 145146,
pp. 4346, 2009.
Michels, A., Wijker, H., and Wijker, Hk., Isotherms of
Argon between 00 C and 1500 C and Pressuures up to
2900 Atmospheres, Physica, Vol. 15(7), pp. 627633,
1949.
Ishiyama, T., Yano, T., and Fujikawa, S., Molecular Dy
namics Study of Kinetic Boundary Condition at an In
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
terface between Argon Vapor and its Condensed Phase,
Phys. Fluids, Vol. 16(8), pp. 28992906, 2004.
Rowlinson, J. S. and Widom, B., Molecular Theory of
Capillarity, Clarendon Press, Oxford, 1982.
Allen, M. P. and Tildesley, D. J., Computer Simulation
of Liquids, Clarendon Press, Oxford, 1987.
Yaguchi, H., Yano, T., and Fujikawa, S., Molecular Dy
namics Study of VaporLiquid Equilibrium State of an
Argon Nanodroplet and Its Vapor, J. Fluid Sci. Tech.,
Vol. 5(2), pp. 180191, 2010.
Thompson, S. M., Gubbins, K. E., Walton, J. P. R. B.,
Chantry, R. A. R., and Rowlinson, J. S., A Molecu
lar Dynamics Study of Liquid Drops, J. Chem. Phys.,
Vol. 81(1), pp. 530542, 1984.
Nijmeijer, M. J. P., Bruin, C., van Woerkom, A. B.,
Bakker, A. F., and van Leeuwen, J. M. J., Molecular
Dynamics of the Surface Tension of a Drop, J. Chem.
Phys., Vol. 96(1), pp. 565576, 1992.
Vrabec, J., Kedia, G. K., Fuchs, G., and Hasse, H., Com
prehensive Study of the VapourLiquid Coexistence
of the Truncated and Shifted LennardJones Fluid In
cluding Planar and Spherical Interface Properties, Mol.
Phys., Vol. 104(9), pp. 15091527, 2006.
Tegeler, Ch., Span, R., and Wagnera, W., A New Equa
tion of State for Argon Covering the Fluid Region
for Temperatures From the Melting Line to 700 K at
Pressures up to 1000 MPa, J. Phys. Chem. Ref. Data,
Vol. 28(3), pp. 779850, 1999.
Miller, Jr., J. W., and Yaws, C. L., Correlation Constants
for Liquids; Surface Tensions of Liquids, Chem. Eng.,
Vol. 83(23), pp. 127129, 1976.
