Direct Molecular Dynamics Simulation of Nucleation during

Subscriber access provided by - Access paid by the | UCSB Libraries

Statistical Mechanics

Direct molecular dynamics simulation of nucleation during supersonic expansion of gas to vacuum Martin Klíma, and Jiri Kolafa J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00066 • Publication Date (Web): 22 Mar 2018 Downloaded from http://pubs.acs.org on March 23, 2018

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 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 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.

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

Journal of Chemical Theory and Computation

Direct molecular dynamics simulation of nucleation during supersonic expansion of gas to vacuum Martin Klíma and Jiˇrí Kolafa∗ Department of Physical Chemistry, University of Chemistry and Technology, Prague, Technická 5, 166 28 Praha 6, Czech Republic E-mail: [email protected]

Abstract We develop a methodology for direct molecular-level simulation of adiabatic expansion of gas through a small orifice to vacuum. The gas gains supersonic speeds, cools down, and nucleates. The proposed approach combines equations of frictionless fluid dynamics with molecular dynamics simulation in an expanding periodic box. There are two key components of the proposed algorithm: (i) a time-reversible integrator tailored to an expanding system, and (ii) an iterative procedure employed to satisfy the condition of steady flow. For a conical nozzle (opening angle 60◦ ) the simulations with argon and water vapor predict cluster sizes in agreement with the experiment. Clusters of irregular shapes observed in the experiment [J. Lengyel et al. Phys. Rev. Lett. 2014, 112, 113401] are not reproduced. The role of friction, turbulence, and sonic boom originating at sharp nozzle edge is discussed.

1 INTRODUCTION Cold clusters of gases or water can be obtained in the laboratory by adiabatic expansion. In this experiment, vapor expands through a nozzle to vacuum, reaches a supersonic speed, adia∗ To

whom correspondence should be addressed

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

batically cools, nucleates to clusters of different sizes, and eventually freezes. 1–4 Such clusters may serve as a proxy to atmospheric aerosol particles 5 and other open questions in the physical chemistry of ice. 6 The process of pickup of various molecules on clusters in the beam, accompanied by molecular dynamics of the pickup process, 7 can provide cluster cross sections. These cross sections are in a certain range of conditions larger than for a ball. 8 The shape of the clusters is therefore far from regular, supposedly as a result of secondary coagulation of cold spherical clusters. Here we try to address this unsolved problem by simulation techniques. Homogeneous nucleation from gas has been intensively studied by atomistic molecular simulations for years. Most of these works have focused on the nucleation kinetics, testing the nucleation theories, 9,10 characterization of particles, 11 or detailed investigation of the nucleation mechanism. Also lattice models may answer fundamental questions of nucleation mechanisms. 12 Argon has been a model molecule in large-scale studies 13 or investigation of finite-size effects. 14 Large-scale studies continued by studying SPC/E water. 15 Typically, constant-volume simulation cells and adiabatic conditions 16 or various cooling protocols 11,17 have been applied. Cooling is often mediated by a carrier gas 18–21 instead of a simpler cooling by a stochastic-type thermostat. If temperature is low enough, nucleation from gas to liquid is followed by crystallization. It has been now well established that freezing of freshly nucleated droplets 22,23 occurs near surface. 24 This observation was supported also by molecular dynamics simulations. 25 At low temperatures the clusters may have fractal structure composed of frustrated grains. 24 For water, we can speak about a frozen nucleus with ice Ih structure from sizes of about 275 molecules. 26 Methodology more relevant to our task has been employed in the adiabatic expansion studies of the nucleation of naphthalene 27 and pharmaceuticals 28 from supercritical carbon dioxide. A cluster-heat-capacity model for study of homogeneous condensation in supersonic water vapor expansions was developed by Borner et al. 29 It combines a mathematical description of collisions, heat balance of growing and evaporating clusters (based on Monte Carlo result of heat capacities of small clusters) and is supported by molecular dynamics of cluster evaporation.

2


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


The temperature and velocity of ideal gas expanding adiabatically through a nozzle can be obtained from the equations of fluid dynamics. 30 The calculations could be done numerically for real gas if we knew its equation of state. However, nucleation is an irreversible process; therefore, the equation of state is ill-defined and/or time-dependent. Using molecular dynamics for obtaining the missing information seems to be natural. However, molecular dynamics can simulate only a small (periodic) box. The key methodology problem can be formulated as: The simulated expanding box does not know anything about the steady flow condition. More precisely, let us start from the ideal gas solution for the time dependence of density and simulate an expanding box. Two quantities are determined as functions of time: pressure and energy. From integrated pressure–volume work during the expansion (which is related to the total energy) we can calculate the velocity of the expanding box. The flux along the trajectory is then obtained from the fluid dynamics equations. However, the flux is not stationary for a nonideal gas. It is necessary to modify the time dependence of density and repeat calculations until a sufficiently steady flow is obtained.

2 THEORY 2.1

Nozzle

Figure 1: Nozzle geometry.

In the experiment, pressurized gas flows through a nozzle to vacuum. We assume slip wall conditions and a steady flow. The position of a tiny sample of gas can be given by one spatial

3



Page 4 of 27

coordinate which is a function of time, z(t). The geometry of a nozzle, see Figure 1, can be then described by a varying cross section, A (z), which is decreasing for z < 0 and increasing for z > 0; the minimum at the nozzle throat is A (0). The limits are A (±∞) = ∞. To be compatible with the motivation experiment, 7 we use the nozzle drilled in a 2 mm thick metal plate with the inlet diameter 2r = 0.1 mm. The nozzle shape is conical with opening angle α. Since the orifice is much smaller than the plate thickness, it can be assumed that the nozzle length is infinite. We estimate the cross section in the outlet area simply by a disk area; this is accurate enough for small opening angles α. The inlet geometry is approximated by a disk of radius r and a quarter of a torus. Thus

A (z) =

   π(r2 − πrz + 2z2 )

for z ≤ 0,

  π r + z tan(α/2)2 for z ≥ 0.

2.2 Supersonic expansion of ideal gas The calculations for ideal gas are a straightforward modification of textbook results. 30 Let gas flow steadily without friction from z = −∞ through the nozzle to z = +∞. The thermodynamic state of the gas along the z-axis is described by temperature, T (z), density, ρ(z), and pressure, p(z). The flow velocity is v(z). The initial values are T (−∞) = T0 , T (−∞) = p0 , ρ(−∞) = ρ0 , and v(−∞) = 0. These four functions obey four equations: (i) the ideal gas equation of state, (ii) the polytropic equations, (iii) the continuity equation

Jw = vρA ,

(1)

where Jw is the total mass flux, which is assumed constant, and (iv) the continuity equation for energy (Bernoulli equation for a stationary non-curl flow of a compressible fluid),

d( 21 v2 ) + dp/ρ = 0.

(2)

Calculation details can be found in the Supporting Information. As the result, we are able to obtain the density, temperature, and other quantities as a func-

4


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


tion of time. The gas reaches the speed of sound (M = 1) at the nozzle throat, z = 0. In addition, all four functions have a singularity in the first derivative at the throat because function A (z) for a conical nozzle has a jump in derivative at the minimum z = 0. This singularity disappears for nozzles with smooth walls (as Laval nozzles, see the Supporting Information for an example). The calculated density–time dependence serves as the first iteration of the expansion protocol in molecular dynamics.

2.3

Real gas

2.3.1 Simulation along given V (t) time dependence It is not possible to simulate directly a small (periodic) box carved out of flowing real gas because we do not know the flux in advance. We thus will conduct a simulation in periodic boundary conditions in a volume which is a function of time, V (t). Then we will calculate the flow from the simulation data and if the flow is not constant, V (t) will be adjusted. Since our task is to perform an adiabatic simulation, energy conservation is crucial; simple box rescaling is not sufficient. The proposed algorithm uses the ideas of the extended Lagrangian barostat by Andersen. 31 Let us consider a system of N point particles with masses mi , positions ~ri , and velocities ~vi = ~r˙i (the dot denotes time derivative) interacting by potential U(rN ), rN = {~r1 , . . . ,~rN }. Note that quantities with a subscript correspond to individual molecules whereas quantities without a subscript apply to the whole box. The sample is in a cubic periodic box of side L, which is a known function of time. The Newton’s equations of motion are modified by ~fi ∂U ~v˙i = − vf~vi , ~fi = − , mi ∂~ri where vf =

d ln L 1 d lnV = dt 3 dt

is calculated from given V (t) dependence (in the original algorithm 31 it is a dynamic variable). The right-hand side contains velocities which are not known in advance if the Verlet integrator (or its clones, as leap-frog) is used. There are three solutions of this problem available:

5



1. A time-reversible method based on the Trotter decomposition of the Liouville evolution operator. 32 2. Iterations until the final velocities are self-consistent with the Verlet integrator. 33 3. Predicting velocities using the TRVP (Time-Reversible Velocity Predictor). 34 We use the third method because it is simpler to code than option 1 and of the same accuracy as option 2. One step of the proposed algorithm for general models with constraints (fixed bond lengths) treated by SHAKE is N

~fi (t) := − ∂U(r (t)) , ∂~r i 1 L(t + h/2) vf (t) := ln , h L(t − h/2) ~vi (t) := TRVPk (~vi (t − h/2),~vi (t − 3h/2), . . . ), ~fi (t) − vf (t)~vi (t), ~ai (t) := mi ~vi (t + h/2) := ~vi (t − h/2) + h~ai (t) (leap-frog), ~ri (t + h) := ~ri (t) + h~vi (t + h/2) (leap-frog), λ (t + h/2) := [V (t + h)/V (t)]1/3 , ~ri (t + h) := SHAKE(rN ; bonds/λ (t + h/2) , ~vi (t + h/2) := ~ri (t + h) −~ri (t), Ekin (t) := Pcfg (t) :=

1 mi ~vi (t − h/2)2 +~vi (t + h/2) , ∑ 4 i " # 1 2Ekin (t) + ∑~ri (t) · ~fi (t) + · · · , 3V (t) i

~ri (t + h) := λ~ri (t + h) (final rescaling to V (t + h)).

Several comments may be useful: 1. The TRVP algorithm predicts the velocity from several previous values (the history). It depends on an integer parameter, k, which determines its time reversibility, O(h2k+3 ). We use k = 2 which guarantees sufficient energy conservation.

6


Page 6 of 27

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


2. The SHAKE algorithm is iterative. Precision of the bond lengths at least 10−7 Å is needed. 3. Rescaling by λ is split into two parts: (i) rescaling ri (t + h) from the leap-frog step by λ (t + h/2) which means that |~ri −~r j | inside SHAKE is rescaled by λ (which is implemented by the target bond lengths scaled by 1/λ ), and (ii) final rescaling of all coordinates by λ which will recover the correct bond lengths. Note that simple rescaling the whole configuration before SHAKE by λ would add terms proportional to v f~ri to the velocities used in the evaluation of the kinetic temperature (calculated from the kinetic energy of molecules). 4. The configurational pressure Pcfg (t) is written schematically. The real code contains sums over pair interactions with the nearest-image distances, constraint forces (calculated during SHAKE), etc. 35 The kinetic energy of the flow, Eflow = 21 mv2 , can be obtained from eq (2) rewritten to the form E˙flow = −V p˙

(3)

by time-reversible second-order numerical integration Eflow (t + h) = Eflow (t) −

V (t) +V (t + h) 2

(4)

× [p(t + h) − p(t)] + O(h ). 2

The sum of this energy and the enthalpy, H(t), should be constant

Etot = H + Eflow = const.

(5)

To achieve this result with high precision, the simulation must not use any corrections (LennardJones cutoff correction nor pressure correction to account for conserved momentum). The code is implemented in package MACSIMUS. 36 A testing run with the expansion of 100 Ar atoms from 120 to 60 K and the ideal-gas-based result for ρ(t) gave the maximum deviation in the Etot (t) curve 0.26 J/mol while the enthalpy

7



20 1000(H+Eflow) + const

H energy / (kJ mol−1)

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 8 of 27

H+Eflow 10

0

−0.2

Eflow 0

0.2

0.4

t / µs

Figure 2: Energy conservation in a test simulation of expansion of 300 SPC/E molecules (T0 = 409 K). H denotes the enthalpy and Eflow is the kinetic energy of the flow. The data are blocked by 1 ns. drops by 1.20 kJ/mol; the relative energy conservation error is 2.2·10−5 . The time step was 4 fs and one cluster formed during the expansion. For water, the relative energy conservation error is higher, 1.6 · 10−3 (see Figure 2), but still acceptable. The relative accuracy of the SHAKE iterative scheme was set to 10−7 because 10−6 had led (in shorter runs) to an increase of the total energy by about 3 · 10−3 relatively. In addition, a change of the algorithm to a “naive” version with an ad hoc rescaling of the box instead of the version given above leads to energy conservation error of several tens %. Since the integration (4) is tedious and we have verified that the total energy, eq (5), is well conserved, we will use simply Eflow = Etot − H for the box kinetic energy. From the box kinetic energy, one can calculate the box velocity as a function of time, p v = 2Eflow /m, and by integration the position of the flowing box in the nozzle, Z t

z(t) = z0 +

v(t 0 )dt 0 ,

(6)

0

where z0 is the initial position, z0 < 0. The flux is given by the continuity equation (1). At the end of the day, flux Jw should not depend on time.

8


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


2.3.2 Initial condition The simulation is started with temperature and density corresponding (using ideal gas formulas) to a time 0.2–0.4 µs before the expected flyby through the nozzle. The temperature at this point is already by a few kelvin lower that the initial temperature T0 . The flux velocity is not zero; it will be determined after the simulation.

2.3.3 Adjusting the flux There are two parameters in the above formulas, z0 and Etot . The direct integration of eq (6) with simulated v(t) leads in general to reaching the throat (z = 0) at a wrong time. One possibility to fix this problem would be a separate integration backward and forward from the time of flyby through the throat. For technical reasons (graphical adjusting the total energy) we consider z0 as a free parameter which is set so that the v(z) curve has an inflection point at z = 0. If the value of Etot is taken directly from simulation, quantity Eflow fluctuates as a consequence of a small sample size and may become negative; then, the velocity is undefined. This problem arises in the intake area where the gas velocity is low. We have tried five possibilities (which may be combined) to overcome this problem: 1. Smooth the Eflow time series. The data are already blocked (by 50 ps) which alleviates the problem, but for water not sufficiently. Suppressing fluctuations by smoothing and blocking best mimics the thermodynamic limit N → ∞. Since square root is a concave function, more smoothing leads to higher velocities. However, too much smoothing hides details on the v(t) and Jw curves, especially near the throat. 2. Define v = 0 for Eflow < 0. This approach is to some extent similar to smoothing, albeit now with a bias towards smaller values of v. p 3. Define v = − 2Eflow /m for Eflow < 0. This “arbitrary” options serves for testing purposes and enhances the bias towards smaller velocities. 4. Use such Etot so that Eflow < 0 never happens.

9



5. The total energy is considered as an adjustable parameter. If its value is slightly higher than the simulation-based constant energy, the probability of negative Eflow decreases; it corresponds to nonzero inlet velocity far from the nozzle throat.

2.3.4 Modifying V (t) to guarantee constant flux Jw If the flux calculated by eqs (6) and (1) is not constant enough (which happens for water), the V (t) (or ρ(t)) dependence controlling the simulation should be modified. We rescale the time t by a suitable increasing function s(t),

V (t) := V (s(t)).

The simplest possibility, s(t) = t for t < tnucl and s(t) = at for t > tnucl (where a > 0) led only to a moderate improvement of the flux. After several tests, we have arrived to the following algorithm. First, we smooth the flux for t > tnucl by fitting it to Jw (t) = a(t − tnucl )2 + b(t − tnucl ) + c,

where a, b, c,tnucl are parameters. For t < tnucl , s(t) = 1 (ideal gas) is sufficient. The scaling function for t > tnucl is s(t) = 1 + S

Jw (t) , maxt>0 Jw (t)

where S is a free parameter which is adjusted graphically. Several iterations are required. Instead of the above ad hoc procedure, one might use a time-dependent equation of state of the expanding gas and feed it back to the equations of fluid dynamics. One possibility could be an ideal gas composed of clusters (with certain size distribution) with the heat of condensation taken into account. 37

2.4

Simulation details

Argon is modeled by the traditional (although not very precise) potential. 38 The cutoff is set to 3.52σ = 12 Å. Both the potential and the force are cut and smoothed; 36 the potential is not

10


Page 10 of 27

Page 11 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


shifted. The timestep is 4 fs. Water is modeled by the popular simple SPC/E 39 potential. The cutoff is set to 16 Å. This includes both the Lennard-Jones cutoff and the site-site-based smoothed-and-shifted truncated Coulomb potential with the onset of smoothing (in the units of cutoff) at c = 0.7 40 (denoted as α in Ref. 36). Note that the standard Ewald summation is inefficient for very large boxes; in addition, maintaining the Ewald summation efficient over a range of box sizes leads to problems with energy conservation. The timestep of 1.42 fs is used for water. The connectivity of argon atoms to determine clusters is defined by the Ar–Ar distance up to the first minimum on the radial distribution function, 5.2 Å. For water, the intermolecular H–O radial distribution function minimum at 2.45 Å (corresponding to the hydrogen bond) is used. The averaged cluster sizes are calculated as the geometric mean of all observed clusters except very small ones. Clusters smaller than about 10 molecules cannot be observed in the experiment. We therefore chose 11 as the smallest cluster in the calculation of the averages and mark cases where this choice is inappropriate because the population around 11 is significant. The cross section of a cluster is calculated by sampling a projection to a randomly oriented plane with grid 1 Å×1 Å. Each atom is replaced by its Lennard-Jones diameter combined with the diameter of the smallest possible helium probe (1.25 Å). Averaging over all directions is performed using 16 directions with weights corresponding to Gaussian quadrature on sphere. 41 The simulations ran on a linux PC cluster. A single simulation for water required about two months on a single CPU. Unfortunately, parallelization based on regular domain decomposition, as implemented in MACSIMUS, is inefficient for highly heterogeneous systems with clusters.

2.5

Sonic boom cone

The inlet edge of the conical nozzle is the main obstacle to the flow. We expect that the abrupt change in the flow or even a shock wave is a source of sound waves that expand down the supersonic stream in the form of sonic boom. The shape of this wave is conical with intersection at a point on the nozzle axis.

11



Page 12 of 27

To calculate the position of the sonic boom wave, we adopt the same assumptions as in Section 2.2. Huygens’ wavefronts spread from the nozzle edge and are dragged faster than the speed of sound; the envelope corresponds to the sonic boom wave. The wavefronts are not spherical and we will calculate them numerically. Let us define a position in the nozzle by vector ~r = (z, y), where zˆ is the nozzle axis and yˆ is perpendicular (xˆ need not be considered because of rotational symmetry). The inlet edge is described by point (0, r0 ), the nozzle wall is given by y = r0 + z tan(α/2) for z > 0. As a simplification, we assume that the flux velocity is ~v = v 1,

y r0 + z tan α2

.

(7)

This approximation is good for slippery smooth walls and small α. A ray of sound emanates from the nozzle edge, is dragged by the flow, and refracts because of dispersion (caused by temperature gradient). Let the ray at point ~r and time t be described by the direction ~n, |~n| = 1. The initial condition is~r(0) = (0, r0 ) (nozzle edge). The algorithm (first order Euler integration in time, h is the time step) is:

r(t + h) := r(t) + h[c(t)~n(t) +~v(t)] (propagation), c(t + h) (refraction), c(t) ~n(t + h) (normalization), |~n(t + h)|

nz (t + h) := nz (t) ~n(t + h) :=

where v(t) means v(~r(t)) (flux velocity at position~r(t) as calculated in Section 2.2), and similarly for temperature T (t). Note that in our simplified model v and T do not depend on the y-coordinate. The speed of sound is given by c(t) = (γRT (t)/M)1/2 . The above integration is performed for a number of initial directions ~n = (cos θ , sin θ ), θ ∈ [π + α/2, 2π + α/2]. The results are drawn in the form of wavefronts.

12


Page 13 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


3 RESULTS AND DISCUSSION 3.1

Flux adjustment

We used a nozzle with angle 60◦ in all simulations. Our first task is to verify or adjust the expansion protocol, V (t), so that the calculated flux is constant.

3.1.1 Argon The argon simulations started from the temperature and pressure of real experiment; the initial density was calculated using the ideal gas equation of state. All initial temperatures are above the critical temperature, which is 151 K for real argon and 158 K for the used argon force field. 42 Under these conditions the back-calculated flux is constant within 3% so that no rerun with a modified V (t) dependence is needed.

Table 1: Experimental and simulated cluster sizes for a nozzle with α = 60◦ . T0 and ρ0 are the respective initial temperatures and density, Nmax is the largest cluster detected at simulation end for given condition, hNi and hSi denote the averaged cluster sizes and cross sections for N > 10, respectively. Psat and ρsat are the saturated pressure and density. The uncertainties of cluster sizes are 20–25% in the experiment, smaller in the argon simulations, and difficult to estimate for water.

Ar

H2 O

T0 ρ0 (K) (kg/m3 ) 173 2.78 193 7.48 203 11.85 382 0.79 409 1.75 424 2.43

Psat 43 (bar)

1.38 3.21 4.87

experiment ρsat 43 hNi 4 (kg/m3 ) 25 100 192 0.80 45 1.76 134 2.60 234

hSi 8 Psat 44,45 (Å2 ) (bar) † 260 556 790 245 0.61 690 1.56 1100 2.49

†

∗

simulation Nmax hNi 81 289 543 131 207 210

20∗ 92 234 38∗ 70∗ 202

hSi (Å2 ) 164 413 715 194 282 531

Extrapolated. The geometric mean strongly depends on the minimum cluster included, see Figures 5 and 6.

13



3.1.2 Water In real experiment, the inlet water vapor is saturated. We thus calculated the initial density from the ideal gas equation of state using the experimental saturated vapor pressure. Since the saturated vapor pressure of the model SPC/E water is lower than the experimental one (see Table 1), the source vapor is moderately supersaturated. (Homogeneous nucleation needs supersaturation of about 4 in real experiment and even more in short simulations.) In addition, the real nozzle is preheated to a temperature by about 5 K higher than T0 to avoid condensation of frost. The calculated flux of water vapor fluctuates much more than of argon with the same number of molecules. The effect can be best explained using eq (3). Pressure of gas can be in statistical thermodynamics expressed by several equivalent ways. 35,36 If we describe the system by centers of mass of water molecules and their orientations, the instantaneous (kinetic and configurational) pressure is given by

p=

∑i mi v2i ∂U − , 3V ∂V

where the sum is over all molecules, the partial derivative is interpreted using scaled coordinates (of the centers of mass with fixed molecular geometry and orientations, as in the virtual volume change method 46 ). The second term is small in the gas phase because it corresponds to intermolecular interactions. The first term contains only the translational degrees of freedom; therefore, it fluctuates because of energy transfer between rotational and translational degrees of freedom. Hence the kinetic energy of flow integrated by eq (3) fluctuates, too, especially in small systems. No such fluctuation is possible for argon. Our data are blocked by 50 ps. In order to determine the V (t) function, we have run many simulations with 100 molecules. Without further smoothing, the requirement of nonzero Eflux (point 4 of Section 2.3.3) leads to a high total energy and the averaged inlet velocities of 200 m/s and a large initial flux. With smoothing (over a rectangular window 1.5 ns wide), the inlet flux velocity drops to 100 m/s (see the dotted lines in Figure 3 left). Replacing the undefined velocity by zero (point 2) or a negative value (point 3) leads to fluxes indistinguishable within

14


Page 14 of 27

Page 15 of 27

4

Jw/(mg s−1)

1 v/(km s−1)

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


ideal gas adjusted

0

0

1

3

2

1

2

ideal gas adjusted

0

t/µs

1

2 t/µs

Figure 3: Velocity v(t) (left) and flux Jw (right) as functions of time for a simulation box with 100 water molecules (T0 = 409 K) flying through a nozzle. Blue lines show the results obtained using the ideal gas expansion protocol, red lines are the results with modified V (t). Solid lines are results with point 2 of the list of Section 2.3.3 applied, dotted lines with point 4. the scale of graphs. Similarly, using unsmoothed (only blocked) data gives the same results (the adjustment is more tedious, though); therefore, we do not show them in Figure 3. The flux in the inlet part still fluctuates (before final smoothing to draw the figure even more), but the more important part with supersonic speeds after adjustment is constant to 8%. In addition, the flux in the supersonic region is only marginally affected by flux details in the inlet part. Even though the inlet velocities of about 100 m/s may seem large in comparison with final 1.3 km/s, one should compare the corresponding kinetic energies. The dependence of density on time adjusted using 100 particles was used for final runs using 1000 molecules. The final simulations were performed with N = 1000 where the problems of negative kinetic energy are much less severe; however, we cannot afford several microsecond simulations just to adjust the flux. Nevertheless, the flux is constant with a lower precision of 12%, to be compared to 35% with the original ideal gas expansion protocol. Figure 4 shows the final time development of temperature during expansion. The temperature is higher than the ideal gas prediction because of the condensation heat released. In addition, the rotational degrees of freedom cool down slower than the translational ones. (The used rigid model of water does not describe vibrations.) This is because the system contains a lot of free molecules; in clusters, both motions equilibrate within a few picoseconds.

15



400 MD: rotations MD: average MD: translations ideal gas model

T/K

300

200

100

0

0

1

2

3

t/µs

Figure 4: Time development of temperature during expansion of water vapor (T0 = 382 K, p0 = 0.79 bar) and the equipartition of translational and rotational degrees of freedom.

500

80 T = 203 K T = 193 K T = 173 K

400

T = 203 K T = 193 K T = 173 K

60

300

n(N)

Nmax

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 16 of 27

40

200 20

100

0

0 0

1

2

1

2 log10(N)

t/µs

Figure 5: Left: Evolution of the size of the largest cluster for argon. The run with the largest final cluster from five simulations with 10,000 atoms was selected. Sampled by 10 ns, the enlarged inset by 1 ns. Right: Histogram of the logarithm of cluster sizes for argon, based on the last 100 ns of five trajectories; n(N) is the averaged number of clusters in a histogram bin.

16


Page 17 of 27

3.2

Cluster analysis

3.2.1 Argon For every initial condition, five independent simulations at least 1.8 µs long were performed with 10,000 atoms. The evolution of the largest cluster obtained is shown in Figure 5 (left). The distribution of clusters calculated from merged final 100 ns chunks of all simulations is shown in Figure 5 (right). It is approximately log-normal in agreement with the experiment. Also the predicted averaged cluster sizes, Table 1, match well the experiment. Large clusters still keep slowly growing after 2 µs of nucleation, but not so much to violate the observed agreement with the experiment. The typical mechanism of cluster growth is accretion. Less often two clusters collide and merge, then the released surface energy heats the cluster and a few atoms evaporate, which is documented in the inset in Figure 5.

3.2.2 Water 10 200

T = 424 K T = 409 K T = 382 K n(N)

Nmax

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


100

5

T = 424 K T = 409 K T = 382 K 0

0 0

1

2

1

2 log10(N)

t/µs

Figure 6: Left: Evolution of the size of the largest cluster for water. The run with the largest final cluster from five simulations with 1000 molecules was selected. Sampled by 10 ns. Right: Histogram of the logarithm of cluster sizes for water, based on the last 100 ns of five trajectories; n(N) is the averaged number of clusters in a histogram bin.

The water simulations were at least 2.6 µs long. The predicted cluster size for the hottest system is in a good agreement with the experiment. For two colder systems the averaged cluster sizes strongly depend on the threshold size; e.g., for T0 = 409 K and clusters with

17



N < 16 omitted we get hNi = 129 instead of 70. We can conclude that the simulation results are in a weak agreement with the experiment. There are two main sources of inaccuracy: 1. Except for the coldest system, there is only one large cluster in the box, see Figure 6. The cluster statistics is not reliable and larger clusters cannot grow at all. As seen from the histogram, the colder systems contains more smaller clusters; however, the distribution is far from log-normal. A larger number of molecules and longer simulations are needed to obtain more relevant results. 2. The initial system is supersaturated. Thus, larger clusters will likely form. The less common water model SPC with better description of the vapor saturation curve would have been more suitable for the simulation. However, preliminary results with the TIP4P/2005 water model 47 are not too far from the SPC/E results. On the other hand, decreasing the density to be on the model saturation curve would decrease the nucleation rate and cluster sizes and would be much less realistic than using the supersaturated initial gas. Both biases are of opposite sign and may have fortuitously cancelled each other.

3.2.3 Clusters of irregular shape?

Figure 7: Typical final cluster of 195 molecules obtained during simulation in a box with 1000 water molecules and T0 = 424 K. All obtained water clusters are almost spherical in disagreement with our motivation experiment, 8 see Figure 7. Although we are in the range of sizes where the observed anomalous increase of the cross section is small, there is no evidence within our model for any process that could lead to irregular cluster sizes. We have inspected several collisions of simulated

18


Page 18 of 27

Page 19 of 27

water clusters. They always resulted in a spherical shape in a few picoseconds. To be more specific, for water at the initial temperature of 409 K after a 2.5 µs expansion the averaged translational temperature of the system reaches 140–150 K and the rotational temperature is 170–180 K. However, the largest clusters are even warmer, 230–240 K, because evaporative cooling takes time. We also let collide two pairs of such final clusters; the merged clusters were spherical and the temperature increased to about 250 K. Even though SPC/E water has much lower freezing temperature than real water, even real water stays in supercooled liquid state at this condition. Clusters will continue to cool down by evaporation and may eventually stick together in the glassy or frozen state. But because of a low density of the jet, this process has a low probability without any disturbance of the flow.

200

y/µm

50 y/µm

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


100

0

0 0

50

0

z/µm

100

200

300

z/µm

Figure 8: Huygens’ wavefronts and the sonic boom in nozzles with the opening angles α = 30◦ and 60◦ for water (T0 = 409 K, saturated pressure) and the ideal gas approximation. The top half of the nozzle longitudinal section is shown. There are several phenomena that may violate the assumptions of our model. Let us assume first that the flow is not turbulent. Friction at the boundary layer near walls causes heating; in addition, the nozzle is kept at a higher temperature. In this case, it would be feasible and more accurate to use a two-dimensional flow model taking into account the center–wall velocity gradient. On the other hand, heat conduction is slow compared to supersonic speeds and the outer parts of the beam are skimmed. We therefore believe that our model is accurate enough for the central part of the stream provided that the flow is not turbulent. However, the sharp edge of the conical nozzle may cause turbulence, recirculation near the throat, or flow separation. Let us discuss consequences of sound waves generated at the edge.

19



They spread downstream in the form of a sonic boom. The calculated shape of the boom as the envelope of Huygens’ wavefronts is shown in Figure 8. The boom cone intersects at nozzle center in time t = 0.185 µs for α = 30◦ and 0.33 µs for α = 60◦ (in the latter case the accuracy is lower because of approximation (7)). The boom wave passing through the vapor causes an abrupt pressure and temperature increase followed by a decrease. At the conical intersection we therefore expect a large local disturbance in the density, pressure and temperature. In order to estimate the affected volume (size of the focus), we use Stokes’ formula (Ref. 48, page 302) for the exponential attenuation time of sound wave intensity,

τ = 3ρλ 2 /8π 2 µ,

where µ ≈ 0.012 mPa s is the dynamic viscosity of steam and λ the sound wavelength. For too short wavelengths (comparable to the mean free path in water vapor which is approximately 0.1 µm) the sound is too attenuated, for longer wavelengths the intersection region is too wide. Therefore, the order-of-magnitude estimate of the affected region is given by the wavelength for which τ equals the time of sound propagation t (the wave amplitude is attenuated to 1/e); numerically λ ≈ 10 µm (or slightly more because there is already some percentage of clusters present). This volume is small enough for the sonic boom to affect the flying clusters and gas. On the other hand, at time of intersection the clusters are still relatively small and hot so that they relax to a spherical shape if they collide. We therefore conclude that the presumed accretion of clusters is not caused by interacting sonic boom waves. The sonic boom waves may reflect and meet again farther in the nozzle, but they will be likely too weak. In order to assess the role of turbulence further, we calculated the Reynolds number for the initial steam temperature 409 K. (For the calculation, we used the ideal gas approximation, see the Supporting Information, eq (9), giving the speed at the throat 464 m/s and temperature 351 K.) The result is Re = 4300. The turbulence onset is Re > 2300 for an incompressible fluid in a long cylindrical pipe. Although any comparison can be only qualitative because supersonic flow in a cone differs from a cylinder, the sharp edge of our nozzle further enhances turbulence. We therefore conclude that vortices and density fluctuations may form and spread downstream.

20


Page 20 of 27

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


Collisions of already frozen clusters farther in the nozzle might be thus possible, which might lead to larger clusters of irregular shapes.

4 CONCLUSIONS The simulated cluster sizes for both argon and water match well the experiment within available precision (which is low for water). We therefore believe that our approach is soundly based and can describe nucleation in an adiabatic supersonic expansion, including all molecular-level effects like accretion of particles, collision of clusters, and evaporative cooling. On the other hand, the observed anomalous water or ice clusters of irregular shapes 8 were not found. The sonic boom cone expanding from the nozzle inlet edge intersects at a point where it might increase the aggregation rate, but the clusters here are still in a liquid state. Our model currently does not take into account friction nor turbulence. Whereas friction is probably negligible, the anomalous clusters might be be caused by turbulence spreading downstream from the sharp inlet edge. The simulated systems for water are too small for interestingly large clusters to develop. More sophisticated hardware (as GPUs) is needed to allow for repeated calculations with thousands water molecules for several microseconds in a narrow (α = 30◦ ) nozzle.

Acknowledgement This work would have been impossible without motivation and fruitful discussions with Michal Fárník and Juraj Fedor from the Heyrovský Institute of Physical Chemistry, Prague. Financial support by the Czech Science Foundation, Grant No. 15-12386S, is acknowledged.

Supporting Information Available SI.pdf: Derivation of all needed equations describing adiabatic flow of ideal gas through a nozzle and comparison of a conical nozzle with a parabolic one. free of charge via the Internet at http://pubs.acs.org/.

21


This material is available


References (1) Christen, W.; Rademann, K.; Even, U. Supersonic beams at high particle densities: Model description beyond the ideal gas approximation. J. Phys. Chem. A 2010, 114, 11189–11201. (2) Johnston, R. L. Atomic and molecular clusters; Taylor & Francis: London, 2002. (3) Buck, U.; Krohne, R. Cluster size determination from diffractive He atom scattering. J. Chem. Phys. 1996, 105, 5408–5415. (4) Bobbert, C.; Schutte, S.; Steinbach, C.; Buck, U. Fragmentation and reliable size distributions of large ammonia and water clusters. Eur. Phys. J. D 2002, 19, 183–192. (5) Kolb, C. E.; Worsnop, D. R. Chemistry and composition of atmospheric aerosol particles. Ann. Rev. Phys. Chem. 2012, 63, 471–491. (6) Bartels-Rausch, T. Ten things we need to know about ice and snow. Nature 2013, 494, 27–29. (7) Lengyel, J.; Koˇcišek, J.; Poterya, V.; Pysanenko, A.; Svrˇcková, P.; Fárník, M.; Zaouris, D. K.; Fedor, J. Uptake of atmospheric molecules by ice nanoparticles: Pickup cross sections. J. Chem. Phys. 2012, 137, 034304. (8) Lengyel, J.; Pysanenko, A.; Poterya, V.; Slavíˇcek, P.; Fárník, M.; Koˇcišek, J.; Fedor, J. Irregular shapes of water clusters generated in supersonic expansions. Phys. Rev. Lett. 2014, 112, 113401. (9) Karthika, S.; Radhakrishnan, T. K.; Kalaichelvi, P. A review of classical and nonclassical nucleation theories. Cryst. Growth Des. 2016, 16, 6663–6681. (10) Sosso, G. C.; Chen, J.; Cox, S. J.; Fitzner, M.; Pedevilla, P.; Zen, A.; Michaelides, A. Crystal nucleation in liquids: Open questions and future challenges in molecular dynamics simulations. Chem. Rev. 2016, 116, 7078–7116.

22


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


(11) Lümmen, N.; Kraska, T. Investigation of the formation of iron nanoparticles from the gas phase by molecular dynamics simulation. Nanotech. 2004, 15, 525–533. (12) Binder, K.; Virnau, P. Overview: Understanding nucleation phenomena from simulations of lattice gas models. J. Chem. Phys. 2016, 145, 211701. (13) Diemand, J.; Angelil, R.; Tanaka, K. K.; Tanaka, H. Large scale molecular dynamics simulations of homogeneous nucleation. J. Chem. Phys. 2013, 139, 074309. (14) Wedekind, J.; Reguera, D.; Strey, R. Finite-size effects in simulations of nucleation. J. Chem. Phys. 2006, 125, 214505. (15) Angelil, R.; Diemand, J.; Tanaka, K. K.; Tanaka, H. Homogeneous SPC/E water nucleation in large molecular dynamics simulations. J. Chem. Phys. 2015, 143. (16) Kraska, T. Molecular-dynamics simulation of argon nucleation from supersaturated vapor in the NVE ensemble. J. Chem. Phys. 2006, 124, 054507. (17) Römer, F.; Fischer, B.; Kraska, T. Investigation of the nucleation and growth of methanol clusters from supersaturated vapor by molecular dynamics simulations. Soft Mat. 2012, 10, 130–152. (18) Lümmen, N.; Kraska, T. Homogeneous nucleation of iron from supersaturated vapor investigated by molecular dynamics simulation. Aerosol Sci. 2005, 36, 1409–1426. (19) Lümmen, N.; Kraska, T. Molecular dynamics investigation of homogeneous nucleation and cluster growth of platinum clusters from supersaturated vapour. Nanotech. 2005, 16, 2870–2877. (20) Römer, F.; Kraska, T. Homogeneous nucleation and growth in supersaturated zinc vapor investigated by molecular dynamics simulation. J. Chem. Phys. 2007, 127, 234509. (21) Lümmen, N.; Kraska, T. Homogeneous nucleation and growth in iron-platinum vapour investigated by molecular dynamics simulation. Eur. Phys. J. D 2007, 41, 247–260.

23



(22) Manka, A.; Pathak, H.; Tanimura, S.; Woelk, J.; Strey, R.; Wyslouzil, B. E. Freezing water in no-man’s land. Phys. Chem. Chem. Phys. 2012, 14, 4505–4516. (23) Bhabhe, A.; Pathak, H.; Wyslouzil, B. E. Freezing of heavy water (D2 O) nanodroplets. J. Phys. Chem. A 2013, 117, 5472–5482. (24) Modak, V. P.; Amaya, A. J.; Wyslouzil, B. E. Freezing of supercooled n-decane nanodroplets: from surface driven to frustrated crystallization. Phys. Chem. Chem. Phys. 2017, 19, 30181–30194. (25) Vrbka, L.; Jungwirth, P. Homogeneous freezing of water starts in the subsurface. J. Phys. Chem. B 2006, 110, 18126–18129. (26) Pradzynski, C. C.; Forck, R. M.; Zeuch, T.; Slavíˇcek, P.; Buck, U. A fully size-resolved perspective on the crystallization of water clusters. Science 2012, 337, 1529–1532. (27) Römer, F.; Kraska, T. Molecular dynamics simulation of naphthalene particle formation by rapid expansion of a supercritical solution. J. Phys. Chem. C 2009, 113, 19028–19038. (28) Römer, F.; Kraska, T. Molecular dynamics simulation of the formation of pharmaceutical particles by rapid expansion of a supercritical solution. J. Supercrit. Fluids 2010, 55, 769–777. (29) Borner, A.; Li, Z.; Levin, D. A. Development of a molecular-dynamics-based clusterheat-capacity model for study of homogeneous condensation in supersonic water-vapor expansions. J. Chem. Phys. 2013, 138, 064302. (30) Landau, L. D.; Lifshitz, E. M. Course of theoretical physics. 6. Fluid mechanics, 2nd ed.; Pergamon Press: Oxford, 1987. (31) Andersen, H. C. Molecular dynamics simulations at constant pressure and/or temperature. J. Chem. Phys. 1980, 72. (32) Martyna, G.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit reversible integrators for extended systems dynamics. Mol. Phys. 1996, 87, 1117–1157.

24


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


(33) Toxvaerd, S. Algorithms for canonical molecular dynamics simulations. Mol. Phys. 1991, 72, 159–168. (34) Kolafa, J.; Lísal, M. Time-reversible velocity predictors for Verlet integration with velocity-dependent right-hand side. J. Chem. Theory Comput. 2011, 7, 3596–3607. (35) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Clarendon Press: Oxford, 1986. (36) Kolafa, J. MACSIMUS. http://www.vscht.cz/fch/software/macsimus, 2018; accessed March 17. (37) Wyslouzil, B.; Wilemski, G.; Beals, M.; Frish, M. Effect of carrier gas pressure on condensation in a supersonic nozzle. Phys. Fluids 1994, 6, 2845–2854. (38) Dumitrescu, L. R.; Smeulders, D. M. J.; Dam, J. A. M.; Gaastra-Nedea, S. V. Homogeneous nucleation of water in argon. Nucleation rate computation from molecular simulations of TIP4P and TIP4P/2005 water model. J. Chem. Phys. 2017, 146, 084309. (39) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269. (40) Kolafa, J.; Mouˇcka, F.; Nezbeda, I. Handling electrostatic interactions in molecular simulations: A systematic study. Collect. Czech. Chem. Commun. 2008, 73, 481–506. (41) Stroud, A. H. Approximate calculation of multiple integrals; Prentice-Hall: Englewood Cliffs, NJ, USA, 1971. (42) Thol, M.; Rutkai, G.; Koester, A.; Lustig, R.; Span, R.; Vrabec, J. Equation of state for the Lennard-Jones fluid. J. Phys. Chem. Ref. Data 2016, 45, 023101. (43) Wagner, W.; Pruss, A. The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use. J. Phys. Chem. Ref. Data. 2002, 31, 387–535.

25



(44) SAT-TMMC:

Liquid-vapor

coexistence

properties

–

Page 26 of 27

SPC/E

water

(LRC).

https://www.nist.gov/mml/csd/informatics/sat-tmmc-liquid-vapor-coexistenceproperties-spce-water-lrc, 2018; accessed January 5. (45) Fugel, M.; Weiss, V. C. A corresponding-states analysis of the liquid-vapor equilibrium properties of common water models. J. Chem. Phys. 2017, 146, 064505. (46) Harismiadis, V. I.; Vorholz, J.; Panagiotopoulos, A. Z. Efficient pressure estimation in molecular simulations without evaluating the virial. J. Chem. Phys. 1996, 105, 8469– 8470. (47) Abascal, J. L. F.; Vega, C. A general purpose model for the condensed phases of water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (48) Stokes, G. G., On the theories of the internal friction in fluids in motion, and of the equilibrium and motion of elastic solids. Trans. Cambridge Philos. Soc. 1845, 8, 287– 342.

26


Page 27 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


Graphical Abstract

27


Direct Molecular Dynamics Simulation of Nucleation during

Recommend Documents