Example of a Fluid-Phase Change Examined with MD Simulation

Subscriber access provided by JAMES COOK UNIVERSITY LIBRARY

Article

An Example of Fluid Phase Change Examined with MD Simulation: Evaporative Cooling of a Nano-scale Droplet Takashi Ao, and Mitsuhiro Matsumoto Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.7b02059 • Publication Date (Web): 22 Aug 2017 Downloaded from http://pubs.acs.org on August 26, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Langmuir is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

An Example of Fluid Phase Change Examined with MD Simulation: Evaporative Cooling of a Nano-scale Droplet Takashi Ao and Mitsuhiro Matsumoto∗ Department of Mechanical Engineering and Science, Kyoto University, Kyoto 615-8540, Japan E-mail: [email protected] Phone: +81 (75) 383 3655 Abstract We carried out a series of molecular dynamics simulations in order to examine evaporative cooling of a nano-scale droplet of Lennard-Jones liquid. After thermally equilibrating a droplet at a temperature Tini /Tt ' 1.2 (Tt : triple point temperature), we started the evaporation into vacuum by removing vaporized particles, with monitoring the change of droplet size and the temperature inside. As free evaporation proceeds, the droplet reaches a deep supercooled liquid state of T /Tt ' 0.7. The temperature was found to be uniform in spite of the fast evaporative cooling on the surface. The time evolution of the evaporating droplet properties was satisfactorily explained with a simple one-dimensional phase change model. After a sufficiently long run, the supercooled droplet was crystallized into a polycrystalline fcc structure. The crystallization is a stochastic nucleation process. The time and the temperature of inception were evaluated over 42 samples, which indicate the existence of a stability limit. (Submitted to Langmuir: Manuscript revised on August 11, 2017)

1 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction Molecular simulations have been widely used in science and engineering. In chemistry and chemical engineering, for example, they are utilized for many purposes, such as to determine phase behavior of materials, to investigate detailed static and dynamic properties in condensed phases, and to design new molecules with novel functions. Having focused on various types of fluid phase change in our group, such as evaporation, condensation, vapor nucleation, and cavitation, we owe a lot to Professor Gubbins, among the pioneering researchers in this field; his monograph 1 and reviews 2,3 on the role of molecular simulations have been good guides for us. Being intended to be a tribute to him, this paper describes an example of how molecular simulations work in examining fluid phase change of nano-scale droplets; note that a pioneering work on properties of droplets was also reported by his group. 4 Evaporation from droplets is a key issue of fluid phase change in many fields, such as inkjet printing, mist cooling, spray dry technology, combustion of liquid fuels, and meteorological problems. It is a dynamic process, where heat and mass transfers are closely coupled. A number of experimental and theoretical studies on evaporation dynamics of droplets under various conditions have been reported (see Ref. 5 and references therein). More recently molecular dynamics (MD) simulations have been utilized to obtain detailed information on droplet evaporation; since the pioneering work of Long et al., 6 evaporation of a single component droplet of simple fluids 7–11 and water 12,13 as well as a droplet of aqueous solutions 14 were examined. The target of these researches ranges widely, such as evaluating the evaporation rate, 6,14 comparing the droplet size change with macroscopic models, 8,11 investigating the shape change, 9 and clarifying the evaporation mechanism at the molecular scale. 13 Their simulation setup also varies, from direct temperature control 6,7 to introduction of ambient gas or vapor 9–11,14 through which the temperature and the pressure are controlled. In this paper we focus on evaporative cooling of single-component droplets to elucidate its fundamental behaviors. As evaporation from droplets proceeds, the surface temperature falls 2 ACS Paragon Plus Environment

Page 2 of 27

Page 3 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

due to the latent heat loss, sometimes reaching supercooled states. The evaporation speed depends on various control parameters, such as the droplet size and temperature and the atmosphere, as well as thermodynamic properties of the fluid such as saturated vapor density, heat of vaporization, and surface tension. Since the pioneering work of Ranz and Marshall, 15 a number of kinetic models have been developed for numerical analysis of evaporation, 16,17 but they require empirical or adjustable parameters of fluid properties which are sometimes difficult to obtain for supercooled states. 18 MD simulations of course provide a useful tool to examine this phenomenon in some details, although the size and the time of the system is essentially limited to microscales. Here we report our recent results of MD simulations for a nano-scale droplet of Lennard-Jones (LJ) liquid to investigate the dynamic process of evaporative cooling. During a preliminary simulation, 19 we have observed a spontaneous crystallization from the deeply supercooled liquid state. Ice formation through evaporative cooling of water droplets has been investigated experimentally 20,21 and theoretically. 22–24 For example, it was reported 21 that mm-size water droplets are cooled through evaporation under reduced pressure of 70 Pa, reaching a supercooled state of about −15 ◦ C, and finally become crystallized in several seconds. Crystallization dynamics is of course an attractive target for molecular simulations. A number of researches have been reported on freezing of bulk LJ liquid under canonical ensemble (i.e., temperature-controlled) conditions. 25–29 Also crystallization in LJ clusters have been studied and size dependent transitions are reported via Monte Carlo simulations 30 or canonical ensemble MD simulations. 31 In this paper, we want to investigate the freezing dynamics during evaporative cooling in a way as natural as possible; no artificial control of temperature or heat flux is preferable. Thus the simulations were carried out without any temperature control, i.e., under a “quasi” microcanonical ensemble condition; details are given in the next section.


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 27

Simulation Method and Conditions We performed a series of MD simulations to investigate dynamic behavior of nano-scale droplets. Briefly described in this section are the model and the procedure of the simulations.

Model We adopted a model of single-component monatomic fluid with a truncated Lennard-Jones (LJ) 12–6 potential as

φ(r) =

 [ ]    4 (σ/r)12 − (σ/r)6

(r ≤ rc )

   0

(r > rc )

(1)

between particles of distance r apart, with the cut-off distance rc = 4.0 σ. A conventional metric system is used in the following description; all quantities are expressed in reduced units with particle diameter σ, energy parameter , and particle mass m. The unit time is √

then τ = σ m/. The temperature and the pressure are expressed in the unit of /kB and /σ 3 , respectively, where kB is the Boltzmann constant. The MD simulations were executed using a leapfrog algorithm with time step ∆t = 0.001 τ .

Procedure First we prepare a liquid droplet of radius Rini . Using a velocity scaling technique, the droplet is equilibrated at temperature Tini in a sufficiently large cubic simulation cell with periodic boundary conditions for all three directions. During the equilibration run, some particles spontaneously evaporate from the droplet surface, the droplet size being reduced. After sufficient time (typically 105 steps or 100 τ ) the droplet reaches a thermal equilibrium with the surrounding saturated vapor. Then we start the main run, i.e., “evaporation into vacuum”; for that purpose we simply remove the particles in the vapor region instantly. We define a removal radius Rrm from the 4 ACS Paragon Plus Environment

Page 5 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

droplet center; once a particle passes through the boundary, we remove it from the system. In this simulation, we use Rrm = 70 σ. At low temperatures where saturated vapor density is also very low, we expect that the evaporation dynamics is not much affected by the position of the removal boundary Rrm , since the evaporated molecules move almost ballistically into the removal region. During this main run, no temperature control is imposed; thus the mean temperature of the droplet gradually decreases as evaporated particles carry away a part of energy, which is the latent heat of vaporization. This simulation has no exact correspondence to any conventional statistical ensembles because we remove particles from the system regardless of their chemical potential. It would be possible, however, to suppose that the simulation is done under a microcanonical ensemble if we continue to trace the evaporated particles in infinite space in a similar way to Ref. 12, instead of the removal. We store the thermodynamic properties of the system (e.g., number of particles, potential energy, temperature) and particle data (e.g., positions and velocities) with an appropriate interval (typically 103 steps) for later analyses. The simulation is continued up to 107 steps (104 τ ).


Langmuir

Results and Discussion We have carried out several combinations 19 of droplet size Rini and initial temperature Tini , among which the case of Rini = 15 σ and Tini = 0.8 is mainly discussed in this paper. Note that the triple point temperature Tt of LJ fluid is known to be 0.68 ∼ 0.69 in the reduced

Number of Particles

20000

Temperature [ε/kB]

unit, 32,33 below which liquid can exist only as a supercooled state.

0.8

Potential Energy [ε]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 27

-3.0

15000 10000 5000 0

2000

4000 6000 Time [τ]

8000

10000

0

2000

4000 6000 Time [τ]

8000

10000

0

2000

4000

8000

10000

0.7 0.6 0.5

-4.0 -5.0 -6.0 -7.0 6000

Time [τ]

Figure 1: Change of droplet properties during the evaporation.

Overall change As basic properties of the droplet, the number of constituent particles N , the temperature T evaluated from the mean kinetic energy of particles, and the potential energy per particle u are plotted in Fig. 1. Initial decrease of N and T is very rapid; within 40 τ , N reduced from the initial value ∼ 20, 000 to less than 15, 000 and T from 0.8 to below 0.7. After the droplet reaches a supercooled state T ≤ 0.67 at ∼ 200 τ , the evaporation continues and the 6 ACS Paragon Plus Environment

Fig. 1

Page 7 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Figure 2: Sectional snapshots of the evaporating droplet. The time is shown on each figure in reduced unit τ . temperature keeps decreasing. After t ' 4, 000 τ , we observed a sudden increase in T and a drop in u, which indicates a first-order phase change, most probably solidification; latent heat released during the phase change raises the droplet temperature. This kind of temperature recovery during solidification in bulk systems is well known in many fields, such as metallurgy and energy storage with phase change materials. Also in experiments of evaporating droplets partially insulated from their environments, similar temperature recovery is reported. 21,34 As LJ crystal has fairly a large saturated vapor pressure below the triple point temperature, 35,36 the evaporation continues after the solidification and the droplet temperature starts to decrease again. Examples of snapshot (sectional view) are shown in Fig. 2. As the evaporation proceeds, the droplet shrinks in size partly due to the reduction of the number of constituent particles and partly by the density increase at lower temperature. Rapid solidification occurs between 4, 000 τ and 4, 200 τ , and the whole droplet seems to be crystallized already at 4, 500 τ .


Fig. 2

1.2 1 0.8 0.6 0.4 0.2 0

0.7

-3

Number Density [σ ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Temperature [ε/kB ]

Langmuir

Page 8 of 27

2000 τ 4000 τ 5000 τ 7000 τ

0

5

10 r [σ]

15 2000 τ 4000 τ 5000 τ 7000 τ

0.6 0.5 0.4 0

5

10

15

r [σ]

Figure 3: Change of droplet profiles; (top) number density and (bottom) local temperature. Next we investigate structures of the droplet. Shown in Fig. 3 are the change of density profile and temperature profile as functions of radial distance r from the center of the droplet; for these analyses we assume spherical symmetry for the droplet. Although we expected some inhomogeneity in the droplet, such as lower temperature near the droplet surface due to the latent heat loss, the droplet was found to be surprisingly uniform; both the density and the temperature are essentially constant inside the droplet. We suppose that the droplet size is so small (R = 10 ∼ 12 σ) that heat is rapidly conducted, leading to the homogeneous local temperature. As solidification proceeds between 4, 000 τ and 5, 000 τ , the density and the temperature inside the droplet increase almost uniformly. The temperature reduces again since evaporation continues.

One-dimensional cooling model We try to understand the change of the droplet size and the mean temperature based on a simple one-dimensional thermodynamic model. By assuming a homogeneous spherical


Fig. 3

Page 9 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

droplet, the change rate of its radius R is described by the mass conservation law as [

(

d 4π 3 ρL · R dt 3

)]

= −4πR2 J

(2)

where ρL (T ) is the liquid density of the droplet, and J(T ) is the flux of evaporation into vacuum; both are functions of temperature T . Similarly, the energy conservation law gives an expression for temperature change as [

(

)

]

d 4π 3 ρL · R cp T = −4πR2 JL dt 3

(3)

where L(T ) is the latent heat of vaporization per unit mass, cp (T ) the isobaric heat capacity. We numerically solve eqs. (2) and (3) simultaneously for R(t) and T (t) with appropriate expressions for ρL (T ), J(T ), cp (T ), and L(T ). For that purpose, we employ an empirical expression of Lotfi et al. 37 for the density of saturated LJ liquid as

ρL (T ) = ρc + 0.477 (Tc − T )1/3 + 0.2124 (Tc − T ) − 0.01151 (Tc − T )3/2

(4)

with the critical point temperature Tc = 1.310 and density ρc = 0.314 as the reference state point. Equation (4) has been proposed for normal liquid in the temperature range 0.70 ≤ T ≤ 1.30, but we extrapolate it to supercooled states T < 0.70. Since we treat low temperature and low vapor pressure conditions, the evaporation flux J can be evaluated with assumptions of ideal gas for vapor and the evaporation coefficient 38,39 being unity. In the reduced units, it is expressed as √

J(T ) =

T ρV 2π

(5)

For the saturated vapor density ρV we also employ an empirical expression in Ref. 37 as ρV (T ; R∞ ) = ρc − 0.477 (Tc − T )1/3 + 0.053333 (Tc − T ) + 0.1261 (Tc − T )3/2 9 ACS Paragon Plus Environment

(6)

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 27

for bulk (i.e., R → ∞) saturated vapor. For the case of microscale droplets, however, the vapor pressure is known to be strongly curvature dependent; the Kelvin equation 40 is utilized to take account of the curvature effect, as (

2γ ρV (T ; R) = ρV (T ; R∞ ) exp ρL R0 T R

)

(7)

where γ is the surface tension, R0 the gas constant. Note that eq. (7) may fail for a nanoscale droplet. 41 We now need the temperature-dependent surface tension γ(T ); a linear function obtained by least square fitting of reported data for LJ liquid of normal state 42

γ(T ) = 2.470 − 2.046 T

(8)

is extrapolated to supercooled states. For the heat capacity, we assume a temperature-independent value cp = 4.0 for the low temperature range. 43 Finally an empirical expression for the heat of vaporization (or the enthalpy difference) is also found in Ref. 37 as

L(T ) = 6.4925 (Tc − T )1/3 + 2.9644 (Tc − T )2/3 − 1.8502 (Tc − T )3/2

(9)

Putting them all together in eqs. (2) and (3), we can trace the time evolution of R and T with appropriate initial conditions. Results for Rini = 13 and Tini = 0.8 are shown in Fig. 4 in comparison with our simulation data of similar initial conditions. The overall agreement before starting the solidification (i.e., t < 4000 τ ) is satisfactory, considering that some of the empirical expressions may become questionable in the range of such highly supercooled states. It is possible, in principle, to construct a similar model for the cooling process after solidification (e.g., t > 4500 τ in Fig. 1), since thermodynamic properties of the LJ system along the solid-vapor coexistence line have been reported. 35 However, no reliable data are available


Fig. 4

0.8

MD Model

0.7 0.6 0.5 0

1000

2000 3000 Time [τ]

15.0 Radius [σ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Temperature [ε/kB]

Page 11 of 27

4000

MD Model

14.0 13.0 0

1000

2000

3000

4000

Time [τ]

Figure 4: Comparison between the MD simulation data and a simple evaporative cooling model. for surface properties of LJ crystal. Thus we stop further investigation here, just pointing out that the similarity of temperature change (1000 τ ≤ t ≤ 3000 τ for the supercooled liquid droplet and 4500 τ ≤ t ≤ 6500 τ for the crystalline particle in Fig. 1, for example) may indicate some quantitative resemblance of their thermodynamic properties.

Solid structure As seen in the property change (Fig. 1) and the snapshots (Fig. 2), solidification starts between 4, 000 τ and 4, 200 τ , and the whole droplet seems to be crystallized at 4, 500 τ . To check the change of local ordering, a radial distribution function (RDF) at various steps is calculated for particles inside the droplet, i.e., particles with its distance from the droplet center being less than 8 σ, well inside the droplet surface of R ∼ 12 σ. Results are shown in Fig. 5, compared with a perfect fcc crystal at lower temperature T = 0.20. The shape indicates a liquid structure for early stages t ≤ 4, 000 τ , but the second peak suddenly splits to show a typical fcc-type RDF; structural change during the slow annealing for 4, 500 τ ≤ 11 ACS Paragon Plus Environment

Fig. 5

Langmuir

7

t = 500τ (T=0.60) 2000τ (T=0.52) 4000τ (T=0.48) 4500τ (T=0.57) 8000τ (T=0.47) Bulk FCC (T=0.20)

6 5 RDF

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 27

4 3 2 1 0 0

1

2

3

4

r [σ]

Figure 5: Time evolution of the radial distribution function; time average is taken for 10,000 steps (10 τ ). t ≤ 8, 000 τ is not clearly revealed in the RDF analysis.

Crystal nucleation The solidification is a typical nucleation process, which essentially has a stochastic character. 44–46 Starting the simulation from the macroscopically same state with microscopically different particle configurations, we expect to have different inception times for solidification. Examples are shown in Fig. 6, where results of ten different MD runs are plotted; the

Fig. 6

inception time is widely distributed. For more quantitative statistical analysis, the “inception time” ti is defined here as follows (Fig. 7):

Fig. 7

1. We define the mid-step of nucleation as the time t0 when the potential energy per particle becomes less than a threshold u0 . Looking at data shown in Fig. 6, we adopted u0 = −6.15 . 2. In all cases analyzed here, the droplet at t0 − 200 τ is in a liquid state. 3. We define the inception time ti as the time when the two tangential lines of the potential 12 ACS Paragon Plus Environment

-6.0 -6.1 -6.2 -6.3 -6.4 -6.5 -6.6 4000

5000

6000

7000

Time [τ]

Figure 6: Examples of potential energy change for ten MD runs. -6.00 -6.05 P.E. [ε]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Potential energy per particle [ε]

Page 13 of 27

-6.10 -6.15 -6.20 -6.25 “Inception Time” -6.30 3600

3800

4000

4200

4400

Time [τ]

Figure 7: Definition of inception time for solidification. energy curves at t0 and t0 − 200 τ are crossed. 4. As the potential energy fluctuates with time, data in the range of ±10 τ are used for least square fitting to determine each tangential line. With these steps, the inception time of crystallization ti is uniquely determined for each MD run, although the parameters u0 and the period 200 τ are rather arbitrarily chosen. We continued the MD simulation of a droplet in a canonical ensemble to prepare different initial configurations, with which we conducted 42 independent simulations of evaporation into vacuum. The obtained distributions of ti and the inception temperature Ti at ti are 13 ACS Paragon Plus Environment

Number of Cases

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Number of Cases

Langmuir

12 10 8 6 4 2 0 3000

4000

5000

Page 14 of 27

6000

7000

Inception Time t i [τ]

12 10 8 6 4 2 0

0.46 0.48 0.5 0.52 0.54 0.56 0.58 Inception Temperature Ti [ε/kB ]

Figure 8: Distribution of incubation time and temperature for droplet solidification. shown in Fig. 8. Although the number of samples (i.e., 42) is still insufficient for good statistics, the stochastic character is apparent. The distribution of ti is comparatively symmetric, but that of Ti is much skewed; small number of droplet samples become crystallized at higher temperatures, but chances are that nucleation starts when the temperature reduces below T ' 0.48 (T /Tt ∼ 0.7). It has been reported 47,48 that a stability limit, where the free energy barrier disappears, exists for supercooled liquid clusters. We expect that a similar stability limit exists for the LJ droplets. The average values over 42 samples are

Inception time hti i = 4950 [τ ] Potential energy at ti hui i = −6.13 [] Temperature at ti hTi i = 0.487 [/kB ]

(10) (11) (12)

Note that polycrystalline fcc structures are observed for all 42 cases, and no glassy states exist at the final step (1, 000 τ ), which confirms that simple fluids are very easy to be crys-


Fig. 8

Page 15 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

tallized. 49 As the crystallization proceeds, the droplet heats up due to the latent heat release. Modelling of the temperature increase is an interesting target, but investigation of the dynamics, i.e., how fast the temperature recovers, requires more detailed data, such as the relation between atomic motions and the development of local crystalline order. It is simple, in contrast, to evaluate the magnitude of temperature increase ∆Ts during the solidification. As shown in Fig. 1, the solidification proceeds fast enough, so we can assume that it is a (quasi) adiabatic process. Thus the potential energy reduction of ∆u = 0.12 ∼ 0.15 [] per particle, shown in Fig. 6 as examples, leads to the temperature increase of ∆Ts ' (2/3)∆u = 0.08 ∼ 0.1, which explains the behavior in Fig. 1. Thermodynamic argument gives another possible way to estimate ∆Ts . When most of the released latent heat of solidification is consumed to raise the droplet temperature, the energetics gives

∫

Ti +∆Ts Ti

cp dT ' Ls

(13)

where cp (T ) is the isobaric heat capacity, Ls is the heat of crystallization. Although cp should be an appropriately averaged value of supercooled liquid state and crystalline one, let us just use cp ' 4.0 of the supercooled liquid. 43 We also adopt the value Ls ' 1.0 of liquid close to the triple point. 36 This argument leads to ∆Ts = cp /Ls ' 0.25, which is apparently a highly overestimated value. We suppose that loose structure of the crystallized droplets brings this discrepancy; while Ref. 35 reported the energy per particle u = −6.14 for liquid phase and −7.26 for solid under coexisting conditions close to the triple point, we have −6.10 ∼ −6.05 for liquid droplet and −6.30 ∼ −6.20 for solid just after the crystallization. This indicates that the droplet with crystalline order is still in an energetically excited state, and gradually lose the extra energy in a long (typically 1000 τ ) period.


Page 16 of 27

-6.0 -6.2 -6.4 -6.6 2000 3000 4000 5000 6000 7000 8000 Time [τ] 3000 2000

Nc

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Potential Energy [ε]

Langmuir

1000 0 2000 3000 4000 5000 6000 7000 8000 Time [τ]

Figure 9: Number increase of “crystalline particles”.

Where crystallization starts It is interesting to know where the solidification is initiated; for experiments under normal situations, solidification almost always starts from walls and boundaries, as typically seen in Ref. 50. Molecular simulations provide an ideal setup to examine where and how the phase change starts; it was already reported with MD simulations 51 that droplet surface enhances the crystallization of small gold clusters. So how are the situations in our case? There are several indices for inspecting the local crystalline order, such as potential energy and bond orientations. 52,53 Since we know that an fcc crystalline structure appears in the droplet (Fig. 5), we here adopt a simple way of evaluating the number of nearest neighbors, just to look at where crystalline nuclei are born; a particle is judged to have a crystalline order if it has exactly twelve neighbors within the range of distance rth . The threshold value is chosen to be rth = 1.25 σ, which corresponds to the first minimum position of RDF for an fcc crystal. The number Nc of thus judged “crystalline particles” is plotted in Fig. 9 for the MD run shown in Fig. 1. As Nc increases, the potential energy per particle rapidly reduces, indicating 16 ACS Paragon Plus Environment

Fig. 9

Page 17 of 27

Figure 10: Sequential snapshots; “crystalline atoms” are shown with larger sphere. -3

Density of Cristalline Atoms [σ ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

0.8

t-t i = -100 τ 0τ 50 τ 100 τ 500 τ

0.6 0.4 0.2 0 0

5

10

15

r [σ]

Figure 11: Radial distribution of crystalline atoms. Average is taken over five MD runs. that this criterion for crystalline particles works well. Shown in Fig. 10 are sequential snapshots which indicate the positions of thus defined crystalline atoms. Crystalline atoms scatter in a liquid droplet (t ≤ 4, 000 τ ); a cluster of

Fig. 10

crystalline atoms suddenly appears at t ' 4, 100 τ , which is supposed to be a nucleus for crystallization. For randomly selected five MD runs, the radial distribution of crystalline atoms is investigated. The results are shown in Fig. 11, where the density of crystalline atoms is plotted at several relative times t − ti , where ti is the previously defined inception time. The distribution clearly indicates that crystal nucleation starts more often in inner parts of droplets than surface region.


Fig. 11

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 27

Polycrystalline structure The crystal growth is so rapid due to the deep supercool that we very often observed a polycrystalline structure in the final state; a typical example is seen in Fig. 2, where we see at least three crystalline grains, although only a single nucleus appears in another MD run as shown in Fig. 10. In some cases, several nuclei are generated almost simultaneously, as seen in Fig. 12.

Fig. 12

Figure 12: Example of polycrystalline structure, observed during an MD run different from Fig. 10.

Size dependence So far we have investigated droplets of the same size, the initial radius Rini of which is about 15 σ. Droplets with different size may show different behaviors. As an example, behavior of a smaller droplet with Rini = 10 σ is compared with that of Rini = 15 σ in Fig. 13. We have done only a single MD run for the smaller droplet and no statistical average is obtained, but the inception time is certainly shorter than that of larger droplets. There should be two possible explanations. (1) The ratio of surface area to volume is larger for smaller droplets, leading to more rapid cooling for smaller ones. (2) Internal pressure is higher for smaller droplets, which enhances the solidification of LJ fluid. 25 According to the concept of Laplace pressure, the internal pressure Pint of a spherical droplet with radius R


Fig. 13

Temperature [ε/kB]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

0.8

Internal Pressure [ε/σ 3 ]

Page 19 of 27

1.2 1.0 0.8 0.6 0.4 0.2 0.0

Rini =10σ Rini =15σ

0.7 0.6 0.5 0.4 0

2000

4000 Time [τ]

6000

8000

Rini =10σ Rini =15σ

0

2000

4000

6000

8000

Time [τ]

Figure 13: Comparison between a smaller droplet with initial radius Rini = 10 σ and a larger one with 15 σ. is larger than the ambient pressure Pamb due to the surface tension γ as

Pint = Pamb + 2γ/R

(14)

In our situation of evaporation into vacuum, Pamb ' 0 and thus we expect that Pint is inversely proportional to radius R, although γ has a slight dependence on R. These conjectures are supported by Fig. 13 showing the change of the temperature and Pint , where Pint is evaluated with the virial expression for particles within r = 8 σ from the droplet center for the larger droplet of Rini = 15 σ and r = 5 σ for the smaller of Rini = 10 σ; time average for 10 τ (10, 000 steps) is taken to suppress the fluctuations in Pint . The smaller droplet in liquid state has Pint ∼ 0.3 (its typical radius is R ∼ 8 σ), while the larger one with R ∼ 12 σ has Pint ∼ 0.2. Thus, for nano-scale droplets, the Laplace pressure contributes much to the droplet properties; it is very probable that the large Pini accelerates the crystallization through the shift of solid-liquid phase boundary. An opposite size effect is


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 27

Figure 14: Sectional view at the final step of a smaller droplet with Rini = 10 σ, indicating very loosened surface structure. expected for nano droplets of water, because the melting temperature is lowered upon higher pressure. This may explain the reason that crystallization did not occur in the evaporative cooling simulation of water droplets. 12 It is interesting to note that Pint of the smaller droplet shows a sudden increase (upto ∼ 1.2) after solidification, while little change of Pint is observed between pre- and postsolidification for the larger droplet. This implies that smaller solid clusters have much larger surface tension, which can be brought by the large surface curvature with highly loosened or strained structure, shown in Fig. 14 for example.


Fig. 14

Page 21 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Summary To demonstrate how molecular simulations work as a tool to investigate phase change dynamics, we showed an evaporative cooling process of a nano-scale droplet examined with a molecular dynamics simulation technique. The droplet consists of about 20, 000 particles of monatomic fluid with the Lennard-Jones interaction. During evaporation into vacuum, the temperature gradually decreases and the droplet reaches a deep supercooled state of T /Tt ∼ 0.7. The droplet is so small that the inside temperature is uniform. The change of temperature and droplet size is quantitatively explained with a simple one-dimensional model, based on empirical fluid properties and reasonable extrapolations to supercooled states. After some time in the supercooled state, the droplet becomes crystallized into a polycrystalline fcc structure. The solidification process is a typical homogeneous nucleation of a stochastic character, and the time and the temperature of inception were investigated with 42 samples. The crystal nuclei are generated more often at the central part of the droplet than in the surface region.

Acknowledgement We are grateful to Institute for Information Management and Communication, Kyoto University, for providing their supercomputer facilities.

Supporting Information Available Following movie files can be downloaded. • Movie 1 (File:EvaporativeCooling.m1v): Animation of droplet sectional view • Movie 2 (File:nucleation.m1v): Animation of crystalline particles in the droplet


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1) Gubbins, K. E.; Quirke, N. Molecular Simulation Industria; Routledge, 1997. (2) Gubbins, K. E. The role of computer simulation in studying fluid phase equilibria. Mol. Sim., 1989, 2, 4–6. (3) Gubbins, K. E. Applications of Molecular Simulaton. Fluid Phase Equil., 1993, 83, 1–14. (4) Thompson, S. M.; Gubbins, K. E.; Walton, J. P. R. B.; Chantry, R. A. R.; Rowlinson, J. S. A molecular dynamics study of liquid drops. J. Chem. Phys., 1984, 81, 530–542. (5) Sirignano, W. Fluid Dynamics and Transport of Droplets and Sprays; Cambridge Univ. Press, 1999. (6) Long, L. N.; Micci, M. M.; Wong, B. C. Molecular dynamics simulations of droplet evaporation. Computer Phys. Comm., 1996, 96, 167–172. (7) Bhansali, A. P.; Bayazitoglu, Y.; Maruyama, S. Molecular dynamics simulation of an evaporating sodium droplet. Int. J. Therm. Sci., 1999, 38, 66–74. (8) Walther, J. H.; Koumoutsakos, P. Molecular Dynamics Simulation of Nanodroplet Evaporation. J. Heat Transf., 2001, 123, 741–748. (9) Consolini, L.; Aggarwal, S. K.; Murad, S. A molecular dynamics simulation of droplet evaporation. Int. J. Heat Mass Transf., 2003, 46, 3179–3188. (10) Sumardiono, S.; Fischer, J. Molecular simulations of droplet evaporation processes: Adiabatic pressure jump evaporation. Int. J. Heat Mass Transf., 2006, 49, 1148–1161. (11) Landry, E. S.; Mikkilineni, S.; Paharia, M.; McGaughey, A. J. H. Droplet evaporation: A molecular dynamics investigation. J. Applied Phys., 102, 2007, 124301. 22 ACS Paragon Plus Environment

Page 22 of 27

Page 23 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

(12) Caleman, C.; van der Spoel, D. Temperature and structural changes of water clusters in vacuum due to evaporation. J. Chem. Phys., 2006, 125, 154508. (13) Mason, P. E. Molecular Dynamics Study on the Microscopic Details of the Evaporation of Water. J. Phys. Chem. A, 2011, 115, 6054–6058. (14) Wang, B.; Wang, X.; Chen, M.; Xu, J. Molecular Dynamics Simulations on Evaporation of Droplets with Dissolved Salts. Entropy, 2013, 15, 1232–1246. (15) Ranz, W.; Marshall, W. Evaporation from Drops. 1. Chem. Eng. Progress, 1952, 48, 141–146. (16) Vynnycky, M.; Mitchell, S. L. Evaporative cooling and the Mpemba effect. Heat Mass Transf., 2010, 46, 881–890. (17) Pati, S.; Chakraborty, S.; Som, S. Influenced of ambient vapor concentration on droplet evaporation in a perspective of comparison between diffusion controlled model and kinetic model. Int. J. Heat Mass Transf., 2011, 54, 5480–4584. (18) Debenedetti, P. G. Metastable Liquids; Princeton Univ. Press, 1996. (19) Matsumoto, M.; Tatsumi, J.; Hosoda, M. Evaporation dynamics of microdroplets. Proc. 15th International Heat Transfer Conference, 2014, IHTC15–9429. (20) Shin, H. T.; Lee, Y. P.; Jurng, J. Spherical-shaped ice particle production by spraying water in a vacuum chamber. Applied. Therm. Eng., 2000, 20, 439–454. (21) Satoh, I.; Fushinobu, K.; Hashimoto, Y. Freezing of a water droplet due to evaporation —heat transfer dominating the evaporation-freezing phenomena and the effect of boiling on freezing characteristics. I. J. Refrigeration, 2002, 25, 226–234. (22) Zhang, Z.; Gao, J.; S., Z. Droplet Vacuum Freezing Process Based on the Diffusioncontrolled Evaporation and Phase Transition Mechanism. Sci. Rep. 2016, 6, 35324. 23 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(23) Wang, C.; Xu, R.; Song, Y.; Jiang, P. Study on water droplet flash evaporation in vacuum spray cooling. Int. J. Heat Mass Transf., 2017, 112, 279–288. (24) Mirabedin, S. M.; Farhadi, F. Numerical investigation of solidification of single droplets with and without evaporation mechanism. Int. J. Refrigeration, 2017, 73, 219–225. (25) Rein ten Wolde, P.; Ruiz-Montero, M. J.; Frenkel, D. Numerical calculation of the rate of crystal nucleation in a Lennard-Jones system at moderate undercooling. J. Chem. Phys., 1996, 104, 9932–9947. (26) Moroni, D.; Rein ten Wolde, P.; Bolhuis, P. G. Interplay between Structure and Size in a Critical Crystal Nucleus. Phys. Rev. Lett., 2005, 94, 235703. (27) Trudu, F.; Donadio, D.; Parrinello, M. Freezing of a Lennard-Jones Fluid: From Nucleation to Spinodal Regime. Phys. Rev. Lett., 2006, 97, 105701. (28) Wang, H.; Gould, H.; Klein, W. Homogeneous and heterogeneous nucleation of Lennard-Jones liquids. Phys. Rev. E, 2007, 76, 031604. (29) Peng, L.-J.; Morris, J. R.; Lo, Y. C. Temperature-dependent mechanisms of homogeneous crystal nucleation in quenched Lennard-Jones liquids: Molecular dynamics simulations. Phys. Rev. B, 2008, 78, 012201. (30) Polak, W. Formation of regular polyicosahedral and defected crystalline structures in growing Lennard-Jones clusters. J. Crys. Growth, 2014, 401, 44–50. (31) Malek, S. M. A.; Morrow, G. P.; Saika-Voivod, I. Crystallization of Lennard-Jones nanodroplets: From near melting to deeply supercooled. J. Chem. Phys., 2015, 142, 124506. (32) Hansen, J.-P.; Verlet, L. Phase Transitions of the Lennard-Jones System. Phys. Rev., 1969, 184, 151–161.


Page 24 of 27

Page 25 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

(33) Levesque, D.; Verlet, L. Computer “Experiments” on Classical Fluids. IV. Transport Properties and Time-Correlation Functions of the Lennard-Jones Liquid near Its Triple Point. Phys. Rev. A, 1973, 7, 1690–1700. (34) Cogné, C.; Nguyen, P. U.; Lanoisellé, J. L.; Van Hecke, E.; Clausse, D. Modeling heat and mass transfer during vacuum freezing of puree droplet. Int. J. Referigeration, 2013, 36, 1319–1326. (35) Barroso, M. A.; Ferreira, A. L. Solid-fluid coexistence of the Lennard-Jones system from absolute free energy calculations. J. Chem. Phys., 2002, 116, 7145–7150. (36) Agrawal, R.; Kofke, D. A. Thermodynamic and Structural Properties of Model Systems at Solid-Fluid Coexistence 2. Melting and Sublimation of the Lennard-Jones System. Mol. Phys. 1995, 85, 43–59. (37) Lotfi, A.; Vrabec, J.; Fischer, J. Vapour liquid equilibria of the Lennard-Jones fluid from the N pT plus test particle method. Mol. Phys. 1992, 76, 1319–1333. (38) Yasuoka, K.; Matsumoto, M.; Kataoka, Y. Evaporation and condensation at a liquid surface. I. Argon. J. Chem. Phys., 1994, 101, 7904–7911. (39) Ishiyama, T.; Yano, T.; Fujikawa, S. Moleculear dynamics study of kinetic boundary condition at an interface between argon vapor and its condensed phase. Phys. Fluid., 2004, 16, 2899–2906. (40) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Clarendon Press, 1982. (41) Yaguchi, H.; Yano, T.; Fujikawa, S. Molecular dynamics study of vapor-liquid equilibrium state of an argon nanodroplet and its vapor. J. Fluid Sci. Tech., 2010, 5, 180–191. (42) Salomons, E.; Mareschal, M. Surface tension, adsorption and surface entropy of liquidvapour systems by atomistic simulation. J. Phys.: Cond. Matt., 1991, 3, 3645–3661. 25 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(43) Boda, D.; Lukás, T.; Liszi, J.; Szalai, I. The isochoric-, isobaric- and saturation-heat capacities of the Lennard-Jones fluid from equations of state and Monte Carlo simulations. Fluid Phase Equil., 1996, 119, 1–16. (44) Abraham, F. Homogeneous Nucleation Theory: The Pretransition Theory of Vapor Condensation; Academic Press, 1974. (45) Kashchiev, D. Nucleation; Butterworth, 2000. (46) Kelton, K.; Greer, A. Nucleation in Condensed Matter ; Pergamon, 2010; Vol. 15. (47) Mendez-Villuendas, E.; Saika-Voivod, I.; Bowles, R. K. A limit of stabilitiy in supercooled liquid clusters. J. Chem. Phys., 2007, 127, 154703. (48) Bartell, L. S.; Wu, D. T. Do supercooled liquids freeze by spinodal decomposition? J. Chem. Phys., 2007, 127, 174507. (49) Gotze, W.; Sjogren, L. Relaxation processes in supercooled liquids. Rep. Prog. Phys., 1992, 55, 241–376. (50) Jung, S.; Tiwari, M. K.; Doan, N. V.; Poulilkakos, D. Mechanism of supercooled droplet freezing on surfaces. Nature Commun., 2012, 3:615, 1–8. (51) Mendez-Villuendas, E.; Bowles, R. K. Surface Nucleation in the Freezing of Gold Nanoparticles. Phys. Rev. Lett., 2007, 98, 185503. (52) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Icosahedral Bond Orientational Order in Supercooled Liquids. Phys. Rev. Lett., 1981, 47, 1297–1930. (53) Errington, J. R.; Debenedetti, P. G.; Torquato, S. Quantification of order in the Lennard-Jones system. J. Chem. Phys., 2003, 118, 2256–2263.


Page 26 of 27

Page 27 of 27

Graphical TOC Entry Potential Energy Temperature

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

0.8

Evaporative Cooling Droplet

0.7 0.6 0.5 -3.0 -4.0 -5.0 -6.0 -7.0 0

2000

4000 Time

6000

8000


Example of a Fluid-Phase Change Examined with MD Simulation

Recommend Documents