Molecular Dynamics Simulations on Evaporation of Droplets with Dissolved Salts

Molecular dynamics simulations are used to study the evaporation of water droplets containing either dissolved LiCl, NaCl or KCl salt in a gaseous surrounding (nitrogen) with a constant high temperature of 600 K. The initial droplet has 298 K temperature and contains 1,120 water molecules, 0, 40, 80 or 120 salt molecules. The effects of the salt type and concentration on the evaporation rate are examined. Three stages with different evaporation rates are observed for all cases. In the initial stage of evaporation, the droplet evaporates slowly due to low droplet temperature and high evaporation latent heat for water, and pure water and aqueous solution have almost the same evaporation rates. In the second stage, evaporation rate is increased significantly, and evaporation is somewhat slower for the aqueous salt-containing droplet than the pure water droplet due to the attracted ion-water interaction and hydration effect. The Li+-water has the strongest interaction and hydration effect, so LiCl aqueous droplets evaporate the slowest, then NaCl and KCl. Higher salt concentration also enhances the ion-water interaction and hydration effect, and hence corresponds to a slower evaporation. In the last stage of evaporation, only a small amount of water molecules are left in the droplet, leading to a significant increase in ion-water interactions, so that the evaporation becomes slower compared to that in the second stage.


Introduction
The physics of droplet evaporation in an infinite space has attracted much interest due to the crucial role played in energy engineering, in chemical engineering as well as in environmental processes [1], and has been under investigation for many years by different methods: hydrodynamics [2], kinetic theory [3], and molecular simulations [4][5][6][7][8][9][10][11].Molecular dynamics attempts to simulate the real behavior of Nature by identifying each atom and following their motion in time through the basic laws of classical mechanics [7].The system behavior and temporal evolution of its thermodynamic and transport properties can be obtained by statistically averaging the results of all molecular motions.Molecular dynamics simulating an evaporation process has no need of some assumptions made by CFD (computational fluid dynamics), so this method was adopted to study the droplet evaporation [4][5][6][7][8][9][10][11] and the evaporation of flat thin liquid films on solid surfaces [12][13][14].These studies focused on evaporation of droplets consisting of one component under various conditions, and compared the evaporation rates simulated by molecular dynamics and predicted by classical kinetic theory, such as the D 2 law [5,6,9].The D 2 law was derived based on the droplet evaporation in an infinite space [15], and predicts that the derivative of the square of the droplet diameter with respect to time is constant, or dD 2 /dt = -K, where K is the evaporation constant.However, Semenov et al. [16] presented that for the evaporation of sessile drops on hydrophobic substrates, the evaporation rate is proportional to the radius of the three phase line instead of being proportional to the area of the surface of the droplet.Semenov et al. [17] also investigated the effect of the influence of kinetic effects on evaporation of pinned sessile water droplets of submicrometer size placed on a heat conductive substrate.Their computer simulation model took into account the following phenomena: influence of curvature of the droplet's surface on the saturated vapor pressure above the surface (Kelvin's equation), the effect of latent heat of vaporization, thermal Marangoni convection, and Stefan flow inside an air domain above the droplet.
The evaporation of droplets of dissolved salts (aqueous droplets) has extensive applications in many industrial processes such as crystallization [2,18], electrospraying [19], electrospinning [20], and atmospheric science [21][22][23][24][25]. Starov and Churaev [2] investigated the crystallization process of aqueous solutions in a thin capillary using hydrodynamics, and they presented that the evaporation flux differs significantly from that predicted by the classical solution with a one-component liquid.Recently, some studies have used molecular dynamics simulations to investigate the evaporation of a small water cluster dissolving a single kind of charged ion [19,[26][27][28].Caleman and Spoel [26] simulated the evaporation from a water cluster (N = 216 and 512) containing either Cl .Daub and Cann [27] studied evaporation and condensation of a small cluster (N = 10, 20, 30 and 40) of water or methanol containing one single Ca 2+ , Na + or Cl − ion in either a vacuum or under argon gas.Daub and Cann [19] also studied the evaporation of a water cluster (N = 10, 15 or 20) containing one Na + ion or one Ca 2+ ion under the action of an electric field.Daub and Cann's results demonstrated that the interaction between ions and water molecules affects the evaporation of the cluster [19,27].Kӧhler [21] investigated the process of formation of liquid cloud drops based on equilibrium thermodynamics, where water vapor condensed with existence of nucleus (solutes), and proposed the well-known Kӧhler theory expressed as: where p w is the water vapor pressure outside the droplet, p° is the corresponding saturation vapor pressure over a flat surface, σ w is the droplet surface tension, ρ w is the density of pure water, n s is the moles of solute, M w is the molecular weight of water, and D is the droplet diameter.According to Kӧ hler's theory, the droplet diameter, water surface tension, and molar concentration of the solute significantly affect the water vapor pressure and hence the droplet evaporation rate.Generally, an aqueous solution is electrically neutral and it includes the same number of cations and anions for monovalent salts, the interaction between cations and anions may influence the evaporation properties of the aqueous droplet, however, few molecular dynamics simulations have been carried out considering this process.The salt crystallization process from an evaporating NaCl aqueous solution has been studied by Mucha and Jungwirth [18] using molecular dynamics simulations, but they did not focus on evaporation rates.This work uses molecular dynamics simulations to investigate the evaporation of water droplets containing either dissolved LiCl, NaCl or KCl salt.The droplet is surrounded and heated by a nitrogen gas atmosphere at a constant temperature.The effects of the salt concentration and salt type on the evaporation rate are examined.By analyzing the spatial position and the interaction energy of ions and water molecules, the differences between evaporation rates for various cases are explained in detail.

Interatomic Potential and Initial Configuration
For molecular dynamics simulations, selecting a proper intermolecular potential function and constructing a correct initial configuration of system are important to correctly describe the physical process concerned.This work simulates the evaporation of water droplets with or without dissolved salts.Three salts, LiCl, NaCl or KCl, at various concentrations are added to the water droplet.The droplet is surrounded and heated by nitrogen gas at a constant temperature.Since the system includes nitrogen molecules, water molecules, Li + , Na + , K + , and Cl − ions, the long-range Coulombic force between ions must be considered.Thus, a combined potential model of Lennard-Jones 12-6 potential and Coulombic potential is adopted here, which can be expressed as [29,30]: ( where the subscripts i and j denote ith and jth particles (atoms or ions), q is the charge of particle, r is the distance between particles, σ and ε are the minimum energy and the zero energy separation distance.The water molecules are characterized by the SPC/E model [30].The values of the potential parameters q, σ and ε for the same particles are summarized in Table 1 [29,30].The following mixing rules are adopted to describe the potential parameters between different particles, or: (3) The truncated distances for short-range and long-range forces are taken as 10 Å and the PPPM summation technique [31] is used to modify the long-range Coulomb interaction.The equations of motion are integrated using the Velocity-Verlet algorithm [5].The method of constraints [29] is applied for nitrogen molecules and water molecules to maintain their bond lengths and angles.
The initial configuration of the system is shown in Figure 1, where the nitrogen molecules, water molecules, and ions are distinguished with different colors.The droplet is placed in the center of a cubic box with a side length of 12 nm, and nitrogen molecules surround the droplet.The number of water molecules in the droplet is 1,120, and the number of nitrogen molecules is 600, corresponding to a 16.47 kg•m −3 gas density.The number of salt molecules is assumed to be 0, 40, 80, and 120, respectively, to analyze the effect of salt concentration on droplet evaporation.It is noted that the maximum salt mole concentration is 9.7%, which is less than its saturation concentration, and hence salt crystallization cannot occur.The droplet radius is fixed to 2 nm for all cases.This value corresponds to a density of 1 g•cm −3 as the droplet is composed of pure water.

Preparation of Initial Equilibrium State
Before the onset of evaporation a well-defined system has to be prepared.Initial velocities of particles in both the gaseous phase (nitrogen) and the droplet are generated by assuming a Maxwell-Boltzmann distribution based on the initial temperature of 298 K. Periodic boundary condition is applied to the three coordinates of the box.A time step length δt = 1 fs is used for all cases.The system with the initial configuration is simulated in an NVT ensemble and it reaches an equilibrium state after 100,000 time steps.For each time step, the velocities of the gaseous phase and the droplet are separately rescaled to maintain a constant temperature of 298 K.It is noted that a very small amount of water molecules (less than 10 water molecules for the case in Figure 2a) escape from the droplet and occur in the surrounding vapor when the initial equilibrium state is reached, as shown Figure 2a.
To analyze the temporal evolution of the evaporation rate, one must define whether a water molecule belongs to the droplet or to the vapor.The method originally proposed by Shigeo et al. [32] is adopted here, which is based on counting the number of neighboring water molecules around each water molecule.Neighbor molecules are determined as the molecules within a distance of 4.34 Å from the molecule of interest.A molecule is considered to be in the droplet if its neighbor molecule number n neighbor ≥ 9, in the vapor phase if n neighbor ≤ 1, or in interface region if 2 ≤ n neighbor ≤ 8.The interface is ignored in the present work, so n neighbor ≥ 4 is used as a threshold value to determine the droplets.

Droplet Evaporation
The temperature of the gaseous phase is abruptly increased to 600 K to trigger the droplet evaporation and this temperature is kept by the velocity-rescaling method [4] for nitrogen molecules for every time step throughout the evaporation process for all cases.Total evaporation time is set as 1,600 ps.The instantaneous positions and velocities of particles are recorded every 1 ps and all quantities of interest are calculated statistically for 10 recorded results to reduce statistical fluctuations.Figure 2 presents snapshots for the water droplets with 120 dissolved NaCl molecules and the number Entropy 2013, 15 1237 densities of water molecules, Na + and Cl − ions at different evaporation instants.It can be seen that the droplets deviate from the initial spherical shape and their volume gradually decreases with time, however, Na + and Cl − ions cannot escape from the droplet and finally crystallize as the droplet evaporates completely.

Effect of Salt Concentration
To analyze the effect of the salt concentration on the droplet evaporation, 0, 40, 80, and 120 LiCl molecules are added into 1,120 water molecules, respectively, to prepare aqueous droplets.The initial temperature of the droplet is 298 K, the droplet is heated by 600 K nitrogen gas and evaporates.The temporal evolution of the water molecule number in the droplet is shown in Figure 3.
Figure 3 shows that the evaporation process can be divided into three stages.At the beginning of evaporation the evaporation rates for pure water and three aqueous solutions are low, since only a small amount of heat is transferred to the droplet and water has high latent heat of evaporation, and no visible difference is observed for four droplets which means that Li + and Cl − have not yet affected the evaporation.Later, the evaporation rates increase because more heat is transferred to the droplet and the difference between four droplets occurs, the water droplet evaporates faster than three aqueous droplets and vanishes about at t = 1,150 ps.The aqueous droplet with high LiCl concentration has a lower evaporation rate than that with low LiCl concentration.In the last stage of evaporation (at about t > 1,200 ps), the evaporation rates for aqueous droplets decrease compared to that in the second stage.Figure 4 shows that hydration number of Li + is 3.80 at t = 0 ps, 3.71 at t = 200 ps, 3.45 at t = 400 ps, as well as 3.34 at t = 600 ps, with only 9.2% decrease from t = 0 ps to t = 600 ps; however, hydration number of Cl − is 8.70 at t = 0 ps, 8.68 at t = 200 ps, 8.67 at t = 400 ps, and 8.65 at t = 600 ps.Therefore, only 52 and four water molecules escape the confinement of Li + and Cl − , respectively.At the same period, about 270 water molecules escape from the droplet due to evaporation (Figure 3).The results above demonstrates that the free water molecules with a weak interaction with ions made the biggest contribution to evaporation rate at the beginning of evaporation, and hence no visible difference is observed for pure water and aqueous solution with various LiCl concentrations.The hydration numbers of Li + and Cl − at t = 600 ps for aqueous droplets with 40, 80, and 120 LiCl molecules are listed in Table 2.Although high LiCl concentration leads to a small hydration number, the hydration effect is enhanced because the total number of water molecules bounded by Li + and Cl − is increased.The addition of Li + and Cl − into the water droplet also affects the interaction between water molecules.Table 3 lists the coordination number of water molecular at t = 0 ps, which is defined as the average number of water molecules in a sphere with 0.35 nm radius around a water molecule.The value of 0.35 nm chosen here is based on the fact that it is a standard length to determine the formation of hydrogen bonds between water molecules [26].Table 3 shows that the coordination number of water molecular is reduced for high LiCl concentration, thus, the interaction between water molecules becomes less with increased LiCl concentration.The average interaction energies between water molecules and ions (Li + and Cl

−
) for various LiCl concentrations are calculated by Equation ( 2) and shown in Figure 5.The negative value means that water molecules are attracted by ions.The interaction energy is stronger for high LiCl concentration at t < 1,200 ps.Based on results above, the low evaporation rate of the droplet with high LiCl concentration can be attributed to stronger hydration effect and stronger attractive force to water imposed by Li + and Cl − as compared to that with low LiCl concentration.Figure 5 also shows that the interaction energy is significantly elevated at t > 1,200 ps, because less and less water molecules are left in the droplet, hence, the evaporation becomes slower at t > 1,200 ps.The vapor pressure is low at the beginning stage of evaporation, and it gradually increases as water molecules escape from the droplet.High vapor pressure means larger evaporation resistance, which is another important factor for the slower evaporation rate in the last stage of evaporation than that in the second stage.
The aqueous droplet with 40 dissolved LiCl molecules has the highest evaporation rate at t < 1,200 ps (Figure 3), thus, less water molecules are left in the droplet compared to the droplets with 80 and 120 LiCl, so that each water molecule in the droplet is surrounded by more ions at t > 1,200 ps, which leads to a stronger ion-water interaction for the droplet with 40 LiCl.Therefore, the crossover of curves of ion-water interaction for various LiCl concentrations is observed at t = 1,200 s.

Effect of Salt Category
Due to the difference of interactions between water molecules and various ions, the evaporation rates of droplets with dissolved different salts may be different.Three common salts LiCl, NaCl and KCl are used to analyze this effect.Figure 6 shows temporal evolution of the number of water molecules in the droplets with 120 LiCl, NaCl or KCl molecules.Again, three stages are observed during evaporation for all the three aqueous droplets.The evaporation rates of aqueous droplets are lower than that of pure water droplet, and the slowest is LiCl aqueous droplet, then NaCl and KCl.  Figure 7 shows the distribution of water molecules in a sphere with 0.45 nm radius around Li + , Na + , and K + at evaporation instants of 0, 300, and 1,000 ps.Oxygen atoms are closer to cations than hydrogen atoms due to attracted Coulombic interaction.The number of water molecules around Li + , Na + , and K + differs significantly, more water molecules occur around Li + , then Na + and K + .As the droplet evaporates, the water molecules around cations are gradually reduced.a2) and (a3) KCl at 0, 300, and 1,000 ps; (b1), (b2) and (b3) NaCl at 0, 300, and 1,000 ps; (c1), (c2) and (c3) LiCl at 0, 300, and 1,000 ps.(White balls: H, red balls: O, purple balls: Cl − blue ball: K + , Na + or Li + ).
The radial distribution functions g(r) Cation-Cl − , g(r) Cation-O , g(r) Cl − -O for LiCl, NaCl, and KCl aqueous droplets at various evaporation instants of 0, 600, 1,300, and 1,600 ps are show in Figure 8, where the subscript "Cation" denotes K + , Na + , or Li + , respectively.The positions of first peak of g(r) Cation-Cl − are 0.24 nm, 0.27 nm and 0.33 nm for LiCl, NaCl, and KCl aqueous droplets (Figure 8a), and the positions are almost unchanged throughout the evaporation process.However, the peak values of g(r) Cation-Cl − are elevated with the time, which means that more and more cations and chloride ions aggregate together.Eventually, a crystal will form when all water molecules in the droplet evaporate completely.The peak values of g(r) Cation-O (Figure 8b) and g(r) Cl − -O (Figure 8c) for LiCl aqueous droplet are the largest throughout the evaporation process, then for NaCl and the smallest for KCl.Thus, the strongest hydration effect occurs for LiCl aqueous droplet according to Equation (4).
Comparison of Figure 8b,c indicates that the difference of g(r) Cation-O for three aqueous droplets is more significant than that of g(r) Cl − -O .Therefore, only the average interaction energy between water molecules and cations (K + , Na + , or Li + ) is calculated by Equation ( 2) and is plotted in Figure 9.The attractive force between water molecules and Li + is the strongest, while the weakest is for K + .The results confirm again that the strong hydration effect and attractive force are responsible for the slow evaporation.The results can also be connected to the Hoffmeister series effect [33] in term of structure breakers or structure enhancer cations.Hofmeister series is a classification of ions in order of their ability to salt out.The order of cations is usually given as: K + > Na + > Li + in Hofmeister series, therefore, the present results are in good agreement with the Hofmeister series effects.In Nature, many aqueous solutions include two or more solutes, so it is necessary to discuss the evaporation of water droplet simultaneously dissolved various kinds of salts.The evaporation of KCl + LiCl aqueous droplet is simulated and the results are shown in Figure 10, where KCl20+LiCl60 means that droplet dissolves simultaneously 20 KCl molecules and 60 LiCl molecules.With the same salt concentration, the evaporation rates of KCl20 + LiCl60 and KCl40 + LiCl40 aqueous droplets are between the ones of KCl and LiCl aqueous droplets, and faster evaporation occurs at KCl40 + LiCl40 since K + has weaker hydration effect and smaller attractive force towards water molecules than Li + .

Conclusions
The evaporations of pure water droplets as well as NaCl, KCl and LiCl aqueous droplets are studied by molecular dynamics simulations.The droplets are placed in a gaseous nitrogen surrounding and heated by the surrounding with a constant high temperature of 600 K.The evaporation of aqueous droplet can be divided into three stages with different evaporation rates.The rate is slow at the beginning of evaporation, because only a small amount of heat is transferred to the droplet and water has high latent heat of evaporation.The rate is increased in the second stage as more and more heat is transferred to the droplet, however, the rate is again decreased at the last stage of evaporation due to the much stronger ion-water interaction.
The addition of salts into water droplet results in a slower evaporation rate compared to pure water droplets, which can be attributed to the strong hydration effect and strong attractive force on the water imposed by cations and anions.The evaporation rates for various aqueous droplets are LiCl < NaCl < KCl, and the evaporation becomes slower for high salt concentrations, due to the stronger hydration effect and attractive force occur in LiCl aqueous droplets and at high salt concentration.
The interaction potential model of particles is important for MD simulation results, however, the present study focuses on the effect of the salts concentration and category on the evaporation of the aqueous droplets, the conclusions are expected to be still applicable when a different interaction potential model is adopted.Furthermore, the previous MD studies [5,6,9] showed that because there was no "bulk" liquid for nano-scale droplets, its evaporation rate deviated from the classical D 2 law, which was regarded as a good description of evaporation of micro-and milli-scale droplets.Therefore, it is worth studying further that whether the present results may extend to micro-and milli-scale droplets.
, and their results showed a somewhat slower evaporation rate for clusters with Cl − and Na + than those with H 2 PO 4 − and NH 4 +

Figure 1 .
Figure 1.Initial configuration of system: green balls are N, white balls are H, red balls are O, blue balls are positive ions, and purple balls are chloride ions.

Figure 3 .
Figure 3. Temporal evolution of the water molecule number in the droplet with various salt concentrations.

Figure 5 .
Figure 5. Average interaction energies between water molecules and ions (Li + and Cl − ) for various LiCl concentrations.

Figure 6 .
Figure 6.Temporal evolution of the water molecule number in the droplets with various salts.

Figure 9 .Figure 10 .
Figure 9. Average interaction energies between water molecules and ions (Li + and Cl − ) for various salts.

Table 1 .
Values of potential parameters.

Table 2 .
hydration number of Li + and Cl − for different cases at t=600 ps.

Table 3 .
Coordination number of water molecule at t=0 ps.