Evaporation of liquid droplets of nano- and micro-meter size as a function of molecular mass and intermolecular interactions: experiments and molecular dynamics simulations

Transport of heat to the surface of a liquid is a limiting step in the evaporation of liquids into an inert gas. Molecular dynamics (MD) simulations of a two component Lennard-Jones (LJ) fluid revealed two modes of energy transport from a vapour to an interface of an evaporating droplet of liquid. Heat is transported according to the equation of temperature diffusion, far from the droplet of radius R. The heat flux, in this region, is proportional to temperature gradient and heat conductivity in the vapour. However at some distance from the interface, Al, (where l is the mean free path in the gas), the temperature has a discontinuity and heat is transported ballistically i.e. by direct individual collisions of gas molecules with the interface. This ballistic transport reduces the heat flux (and consequently the mass flux) by the factor R/(R + Al) in comparison to the flux obtained from temperature diffusion. Thus it slows down the evaporation of droplets of sizes R B Al and smaller (practically for sizes from 10 nm down to 1 nm). We analyzed parameter A as a function of interactions between molecules and their masses. The rescaled parameter, A(kBTb/e11) , is a linear function of the ratio of the molecular mass of the liquid molecules to the molecular mass of the gas molecules, m1/m2 (for a series of chemically similar compounds). Here e11 is the interaction parameter between molecules in the liquid (proportional to the enthalpy of evaporation) and Tb is the temperature of the gas in the bulk. We tested the predictions of MD simulations in experiments performed on droplets of ethylene glycol, diethylene glycol, triethylene glycol and tetraethylene glycol. They were suspended in an electrodynamic trap and evaporated into dry nitrogen gas. A changes from B1 (for ethylene glycol) to approximately 10 (for tetraethylene glycol) and has the same dependence on molecular parameters as obtained for the LJ fluid in MD simulations. The value of x = A(kBTb/e11) 1/2 is of the order of 1 (for water x = 1.8, glycerol x = 1, ethylene glycol x = 0.4, tetraethylene glycol x = 2.1 evaporating into dry nitrogen at room temperature and for Lennard-Jones fluids x = 2 for m1/m2 = 1 and low temperature).


Introduction
An evaporation process occurs at the interface between liquid and gas in the region of nanometer size. This region is very hard to access in direct thermodynamic experiments and therefore profiles of thermodynamic parameters governing evaporation in the vicinity of the interface are not easily available from experiments. 1,2 Computer simulations, on the other hand, give full information about evaporation in nano-sized systems, but unfortunately cannot be easily scaled to macroscopic systems studied in experiments. 1 In this paper we combine molecular dynamics (MD) simulations in nano-meter size liquid-gas systems and experiments performed on micro-meter size droplets evaporating into inert gas to study evaporation. We obtain a simple formula for evolution of liquid droplet radius R(t) at time t, during evaporation as a function of thermodynamic and material parameters characterizing the system.
Liquids evaporate into inert gas at the expense of internal energy of the gas. The heat transfer from the gas to the liquid is the limiting step in the process of evaporation. During evaporation the bulk temperature of the gas, T b , is higher than the temperature of the evaporating liquid, T liq . Far from the interface heat is transported towards the interface via temperature diffusion. Close to the interface temperature has a discontinuity (at the scale of mean free path in the gas phase). This effect was first measured by Ward and coworkers [3][4][5][6] in water close to the triple point, where the mean free path was of the order of micrometers. This effect was later observed in molecular dynamics simulations 7 of Lennard-Jones molecules in small nanoscopic systems. Due to the temperature discontinuity heat is transported in a ballistic manner in this region. The ballistic transport modifies the standard formula for evaporation where only temperature diffusion is taken into account. In our previous paper 8 we have established in experiments and MD simulations the following formula for liquid droplet radius R(t) as a function of time t during evaporation into inert gas: where l is the mean free path in the gas phase, Dh = h(T b ) À h(T liq ) is the total change in enthalpy upon evaporation, k v is the heat conductivity in the gas phase at T liq and r liq is the liquid density at T liq . All parameters in this equation apart from A can be simulated for model systems (see the ESI †) or taken from the NIST database for any given gas and liquid. Eqn (1) follows from the conservation of energy: the left hand side of eqn (1) gives the enthalpy consumed by the liquid droplet per unit time and surface area and the right hand side gives the heat flux from the vapour to the liquid. In this equation the ballistic mechanism of energy transfer reduces the heat flux from the vapour to the liquid near the interface by a factor R/(R + Al) in comparison to heat flux obtained purely from the temperature diffusion mechanism. Parameter A is a fit parameter. In the present paper we establish the dependence of A on thermodynamic and molecular parameters of the system (temperature of the gas phase, interaction between liquid molecules, molecular mass of liquid molecules and their molecular size). These parameters are important for energy transfer in molecular collisions between gas molecules and the liquid interface. In particular energy transfer from vapour to liquid is strongly reduced in molecular collisions when the molecular mass of gas molecules is much smaller than the molecular mass of liquid molecules. A similar role is being played by molecular interactions in the liquid. The temperature of gas (proportional to the kinetic energy of gas molecules) acts in the opposite manner. We performed computer simulations of a two-component Lennard-Jones system. Both components differed in mass and size. We performed MD simulations at a few temperatures in the gas phase. Additionally we re-analysed our previous experimental data 8 for very small liquid droplets (i.e. comparable to Al) to obtain reliable values of A. For R c Al eqn (1) is not sensitive to the ballistic transport of energy, and therefore the experimental data should be fit separately for large and small R values. We get DT = T b À T liq from fits to eqn (1) at large R. Next we fix DT and fit A parameter at small R.
In the computer simulations component 2 played the role of inert gas and component 1 formed a liquid droplet. During evaporation component 1 also entered into the gas phase, but component 2 did not dissolve in the liquid phase. The temperature of the bulk gas was set between the critical temperatures of the components. The gas component was well above its critical point and the liquid phase was below its critical point. The liquid droplet maintained a constant temperature during evaporation, with T liq being smaller than T b . After a short initial transient state, the temperature profile T(r,t) was quasi-stationary i.e. it quickly adapted to the changes of the droplet radius R(t) at time t. Thus we had T(r,t) = T(r/R(t)) during the evaporation. A schematic drawing of our system is shown in Fig. 1.
The system studied using computer simulations mimicked the typical two-component systems studied by us experimentally i.e. evaporation of ethylene glycol, diethylene glycol, triethylene glycol and tetra-ethylene glycol into a dry nitrogen atmosphere. We combined our MD simulations and experimental studies of these four compounds to find A as a function of molecular parameters characterizing the molecules in the system.

Computer simulations
Molecular Dynamics (MD) simulations were performed using the Verlet ''leapfrog'' method 9,10 to solve the Newton equations of motion that describe the time evolution of a two-component system composed of low volatile (the first component indexed by 1) and high volatile (the second component indexed by 2) molecules. The molecules of the same kind (ii = 11 or 22) interact via the short range Lennard-Jones (LJ) based potential given in the following form: and the repulsive interaction between components 1 and 2 was given by and 0 for d/s 12 4 2 1/6 . Here i = 1, 2 is the component index and d is the inter-molecular distance. The A 0 , B 0 , C 0 , and C 0 0 constants were adjusted to make the potential differentiable in the whole space. The particles of component 2 (forming an inert gas phase) had always the same mass m 2 = 1, the size parameter s 22 = 1 and the energy parameter e 22 = 1/4. The temperature was always given in 4e 22  Simulations of droplet evaporation were performed in the three dimensional sphere of radius R b . The list of all simulation runs is given in Table S1 of the ESI. † A schematic picture of the system is shown in Fig. 1. The distance from the centre of the sphere is denoted by r. Initially, the sphere was filled with the gas of component 2 at density r 2 = r 0 2 and the liquid droplet that is composed of molecules of component 1 was placed with the centre of the mass at r = 0. Two systems were considered: the first of R b = 340.6, r 0 2 = 0.015, s 11 = 1.0 with the total number of particles N = 2.91 Â 10 6 and the second of R b = 557.3, r 0 2 = 0.0065, s 11 = 2.0, N = 4.91 Â 10 6 . The number of molecules of component 1 was always N 1 = 2.15 Â 10 5 . The particles of component 2 were reflected elastically and specularly from the boundary of the sphere at r = R b . The way the particles are reflected has no influence on the evaporation process. (R b À R(t))/l is many times larger than 1. As a consequence, the particle distribution function close to the liquid surface does not depend on the way the conditions at R b are imposed. The boundary temperature T b was imposed by appropriate scaling of velocities of molecules of component 2 for r 4 R b À ks 22 where k = 3 to 5 depending on r 0 2 . The droplet volume was negligible if compared to the total volume of the system. Far from the droplet surface the density, r 2 , was constant and equal to the initial value, r 0 2 . To make simulation conditions closer to real ones the molecules of component 1 could go out as well as come back into the simulation sphere. For r 4 R b the molecules of component 1 moved independently in stochastic gas (see Fig. 1) characterized by the temperature, T b , and the self-diffusion constant, 9 D 2 . The stochastic medium was simulated introducing the Langevin force into the Verlet scheme (here with no 11 interaction). The random force was taken from the Euler-Maruyama approximation for the stochastic term. 11 The approximation was used to generate random force to move molecules of component 1 at distances r 4 R b . Also for r 4 R b gas 2 was replaced by the Brownian gas. For r o R b no random force was used and the time evolution was always performed using a classical MD method. Molecules of component 1 could cross (in both sides) the R b border and spread all over the space as in real systems for which movement of particles 1 is not restricted (i.e. like for R b -N).
The molecules of component 1 formed the droplet of radius R of the number liquid density, r liq , and the liquid temperature, T liq . The droplet radius, R, was defined via the total number density i.e. r(r = R) = r liq /2. Based on additional simulations we estimated the critical temperature for the interactions given by eqn (2) to T C E 0.65e ii (ESI †). The pressure value as a function of r was calculated during simulations using the standard formula (ref. 9). In the simulations we considered only a quasi-stationary stage during which the droplet evaporated at constant temperature and density. The pressure maintained its value in the initial non-stationary regime (ref. The liquid was under thermodynamic conditions well below its critical point while the gas of component 2 (e 22 = 1/4) formed an inert gas well above its corresponding T C . Such simulations correspond to the experimental case of evaporation of liquids (e.g. water) into the inert gas e.g. air. The values of k v in eqn (1) were determined from additional MD simulations (ESI †) applying the method of Muller-Plathe 12 and Dh was calculated from the simulation data. As in our previous paper, 8 l was estimated from The values of the self-diffusion coefficient of molecules of component 2, D 2 , were determined from the Einstein formula 9 by performing additional MD simulations for the gas of particles of component 2 at T 2 = T b and r 2 = r 0 2 (ESI †). The value of A was determined from the fit to solution of eqn (1) given by eqn (5).

Kinetic parameter A from computer simulations
The evaporation process observed in simulations proceeds via an initial strongly non-equilibrium process characterized by the fast decrease of R, rapid decrease of T liq and increase of r liq .
Thus initially the enthalpy needed for evaporation comes from the decrease of temperature of the droplet. This stage is very short. In the second quasi-stationary stage the temperature in the gas phase develops a profile adjusting to the change in the radius of the droplet R(t), T(r,t) = T(r/R(t)). At this stage of evaporation the temperature of the droplet, T liq , was constant. During this process we observed (Fig. 2) that molecules of component 1 leaving the droplet-gas interface had markedly lower temperature than the temperature of the gas.
In order to determine A we solved eqn (1): where B = k V (T b À T liq )/(Dhr liq ). The dependence of A on m 1 is a consequence of a decrease in efficiency of energy  Fig. 2 in which the temperature discontinuity at the gasliquid interface is much higher (i.e. the energy transfer is slower) for m 1 /m 2 = 8 than for m 1 /m 2 = 1. The temperature jump for the gas temperature profile at the interface is given by This discontinuity increases with A. Thus by comparing the upper and lower panels in Fig. 2 we find that A is a function of molecular mass. In Fig. 3 we plot A as a function of molecular mass for different molecular sizes. We removed any temperature dependence of A in Fig. 3, since we observed that A is inversely proportional to the square root of gas bulk temperature as shown in Fig. 4. The final approximate formula for A as a function of molecular parameters and temperature is as follows (Fig. 3 and 4): where a 1 and a 0 (two dimensionless parameters of the order of 1) weakly depend on their molecular size (range of interactions and structure of the liquid). The main mode of energy transfer between gas and liquid is via direct collisions of gas molecules with the surface of a liquid and exchange of kinetic energy. The interaction parameter between gas and liquid e 12 does not affect this transfer appreciably in rarified gas.   We have verified this physical observation in computer simulations. In our simulations gas 2 is always a rarefied gas well above its critical temperature. We set (ref. 8) e 12 = e 22 and e 22 = ge 11 where the interaction parameter between liquid molecules e 11 = 1. We changed g from 0.05 to 0.0005 without any observable changes in A. We have also verified that A does not depend on the strength of interactions between molecules of component 2 which forms a gas phase. Finally the dependence of A on molecular interaction between molecules of component 1 (forming a liquid) and their masses were tested in experiments performed on a series of similar chemical compounds: ethylene glycol, diethylene glycol, triethylene glycol and tetra-ethylene glycol evaporating into dry nitrogen gas.

Experimental determination of parameter A
We have used previously obtained data 8 on diethylene glycol and triethylene glycol and performed additional experiments on ethylene glycol and tetraethylene glycol.
In the case of real molecules, the dependence of evaporation coefficient A on molecular parameters is not as straightforward as in computer simulations. In MD simulations we could independently change the masses, interaction parameters and sizes of spherical molecules and in this manner we find A as a function of each of these parameters separately. In real molecules their masses, sizes and interaction parameters are inter-dependent. Eqn (7) shows which of the molecular parameters most strongly influence the value of A. Here we further verify eqn (7) in experiments.
Liquids of all these compounds were introduced into the droplet injector with a sterile syringe. Next the droplets were injected into the electrodynamic trap. They were charged during injection and thus could be trapped by appropriate combination of AC and DC electrodynamic fields. The initial droplet radius was several mm. In the electro-dynamic trap single droplets levitated in a small volume set by the trapping geometry and parameters of the system. The evolution of the droplet radius R(t) in time was obtained via the recorded angle-resolved light scattering and analysis based on the Mie scattering theory. Details of the experimental set-up are given in our previous publication. 8 In Fig. 5 we show a radius of a droplet of triethylene glycol (very small droplet) and of ethylene glycol (a slightly large droplet) as a function of time. The data were fit to solution of eqn (1) (see eqn (5)). There were two free parameters in the fit: the kinetic parameter, A, and the temperature difference between the vapour temperature far from the droplet, T b , and the droplet temperature, T liq . The second of the fitting parameters was important for all sizes of studied droplets. The parameter A was important only for small droplets i.e. at times when the droplet size decreased to a value comparable to Al (where the mean free path in nitrogen was l = 66 nm). Therefore in order to find both T b À T liq and A reliably, a droplet, which during evaporation reached a size below B2 mm, was selected. The presence of impurities influences the results very severely 1,12,13 and therefore we used pure liquids. The ballistic mode of evaporation is utterly negligible for droplets larger than B6 mm. Thus, fitting the model to the radius evolution above 6 mm was used to find the T b À T liq parameter only (see the lower panel of Fig. 5). Extending the fit to further evaporation enables finding both T b À T liq and A reliably. T b À T liq should not change in the extended fit and serves as The radius R(t) as a function of time for triethylene glycol droplets evaporating into the dry nitrogen atmosphere. The following parameters characterized the system: Dh = 80 kJ mole À1 , r liq = 1121 kg m À3 , k v = 0.0258 J m À1 s À1 K À1 , l = 66 nm. The blue line is a fit with T b À T liq as the only fitting parameter and A = 0. The red line is a fit with suitably chosen A and T b À T liq (see eqn (5)). Lower panel: The radius R(t) as a function of time for ethylene glycol droplets evaporating into the dry nitrogen atmosphere. The following parameters characterized the system: Dh = 62 kJ mole À1 , r liq = 1109 kg m À3 , k v = 0.0258 J m À1 s À1 K À1 , l = 66 nm. The red line is a fit with suitably chosen A and T b À T liq in eqn (5). The blue line shows the region of droplet sizes affected only by the temperature difference and not by the kinetic parameter, A.

Soft Matter Paper
Open Access Article. Published on 26 July 2017. Downloaded on 3/1/2020 10:48:28 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

View Article Online
a check against the influence of impurities. In short, after the fit of T b À T liq we fixed its value and for longer times (and thus smaller droplets) fitted only A.
In Fig. 6 we show the dependence of A on molecular mass, temperature and intermolecular interactions of liquid molecules. The interaction between molecules is proportional to the enthalpy of evaporation e 11 B Dh (Dh is given in the ESI † for all studied compounds). The dependence of A(k B T/e 11 ) is nearly linear as a function of m 1 /m 2 in very good agreement with MD results presented in eqn (7). The value of A for TTEG (tetraethylene glycol) is loaded with relatively high uncertainty since a TTEG droplet is more prone to contamination (e.g. with water) due to its hour-long evaporation. The dependence of A upon temperature, as predicted by MD calculations, is hardly accessible in our experiments. For the 273-298 K temperature range studied experimentally the square root of temperature changes is B4% while the accuracy of A is at least several %.

Conclusions
Kinetic effects and ballistic transport of energy is being currently studied in computer simulations. [16][17][18][19][20][21][22][23][24][25] The main difficulty in the quantitative prediction of the evaporation flux is the appropriate estimate of the energy flux from the gas to the surface of the liquid in the close vicinity of the interface. In the ballistic region close to the interface gas molecules individually collide with molecules in the liquid, but the collisions involve all nearest neighbours of a given molecule in the liquid and become multiple particle processes that cannot be simply reduced to binary collision. The finite range of intermolecular-interactions couples many liquid molecules during collisions with a gas molecule. Therefore we cannot describe this effect as a pure binary collision with well-known energy transfer. The phenomenological formula for A (eqn (7)) contains information about such processes. The practical message from our study is that slowly evaporating liquids consisting of high molecular mass are affected by the ballistic mode of energy transfer even for large, micro-meter size droplets. In the nanometric size range of droplets, evaporation of all liquids is very strongly influenced by parameter A. Eqn (1) together with eqn (7), the main result of this work, should apply to evaporation of all liquids from LJ liquids to water and for all sizes of droplets. The details of molecular structure affect the evaporation process. In our equations interactions between liquid molecules are mainly included in the enthalpy of evaporation, Dh, (eqn (1)) and in parameter A (eqn (7)). The origin of two parameters (a 0 and a 1 ) in eqn (7) is not known. Most probably they depend on the range of interactions and structure of the liquid and tell us how well molecules in the liquid are coupled during collisions with gas molecules. Therefore eqn (7) i.e. dependence of A on molecular mass, can be used for a series of compounds having similar chemical structures and the corresponding structure of the liquid as we have shown for a series of Lennard-Jones liquids (Fig. 3) and a series of glycols (Fig. 6).
The temperature discontinuity, DT, across an interface crossed by a thermal flux, j, is related to the Kapitza resistance, R k = DT/j. The temperature discontinuity at the liquid dropletgas interface (eqn (6)) (first observed by Ward and coworkers [2][3][4][5][6] and explained by Bond and Struchtrup 26 ) allows us to quantify the Kapitza resistance using energy flux from eqn (1). We find the Kapitza resistance for our system: R k = Al/k v . The typical values of R k are 10 À5 to 10 À6 m 2 K W À1 ; for water 1.3 Â 10 À5 , glycerol 10 À5 , TEG 1.3 Â 10 À6 , TTEG 2.5 Â 10 À6 m 2 K W À1 and similarly for the Lennard-Jones liquid.
The ballistic mode of energy transfer discussed in this paper is important for various phenomena, when characteristic sizes in the system become comparable to the mean free path in the gas. For example this mode of energy transfer arises when gas bubbles are formed in liquid by e.g. laser heating 27 or for the flow of gases through nanoholes of sizes comparable to the mean free path. 28 In this paper we have not discussed the issue of nucleation and growth just after the nucleation. All studied systems (LJ or glycols) were under thermodynamic conditions far from the region where nucleation plays a role. 29 Fig. 6 Parameter A as a function of molecular mass, m 1 , divided by the mass of nitrogen, m 2 . The interaction energy e 11 = 2Dh/(sN A ), where Dh is the enthalpy of evaporation 14 (per mole; see also the ESI †) at 273 K, s is the number of nearest neighbours in the liquid (we set s = 4 from computer simulations of ethylene glycol (EG) 15 ) and N A is the Avogadro number. The dashed line is the linear fit to eqn (7). From the fit to eqn (7) we find a 1 = 0.3 and a 0 = À0.2, to be compared with a 1 = 0.42 to 0.65 and a 0 = 1.40 obtained in simulations of LJ fluid. The Lennard-Jones liquid of spherical molecules studied in MD simulations (Fig. 3) is very much different from the linear chain molecules studied experimentally. Nonetheless eqn (7) works very well for both systems giving very similar values of a 1 (compare Fig. 6 with Fig. 3). However a 1 and a 0 are not universal and vary from system to system depending probably on subtle differences in the range of intermolecular interactions and the structure of the liquid. In general we do not have a physical interpretation of these parameters. For example for water we found A(k B T b /e 11 ) 1/2 = 1.8(m 1 /m 2 = 0.64) and for glycerol A(k B T b /e 11 ) 1/2 = 1(m 1 /m 2 = 3.3). More simulations and experiments on different compounds are needed to elucidate the origin of a 1 and a 0 . Nonetheless, since both a 1 and a 0 are of the order of 1 eqn (7) seems to grasp most of the physical ingredients of ballistic transport of energy in the system and its influence on the evaporation flux.