Lifetime of a Nanodroplet: Kinetic Effects and Regime Transitions

A transition from a d to a d law is observed in molecular dynamics (MD) simulations when the diameter (d) of an evaporating droplet reduces to the order of the vapor’s mean free path; this cannot be explained by classical theory. This Letter shows that the d law can be predicted within the Navier-Stokes-Fourier (NSF) paradigm if a temperature-jump boundary condition derived from kinetic theory is utilized. The results from this model agree with those from MD in terms of the total lifetime, droplet radius, and temperature, while the classical d law underpredicts the lifetime of the droplet by a factor of 2. Theories beyond NSF are also employed in order to investigate vapor rarefaction effects within the Knudsen layer adjacent to the interface.

A sound knowledge of the evaporation of nanodroplets is of great importance to many applications, such as combustion and spray drying [1], and the design of next-generation evaporative cooling nanodevices [2]. Understanding the mechanisms underlying the evaporation in nanodevices requires observation at spatial and temporal resolutions that challenge current experimental techniques. While computational atomistic descriptions, such as molecular dynamics (MD) simulations, are capable of modeling evaporation at nanoscales, it is unrealistic to use them to study multiscale problems-processes spanning a wide range of time and length scales-because of memory and computational time limitations. This motivates the consideration of continuum models that go beyond the classical Navier-Stokes-Fourier (NSF) description.
From a fundamental viewpoint, as well as for benchmarking, the process of evaporation of a spherical droplet is important. The continuum description, offered by the NSF equations, agrees with a molecular description when the droplet radius (a ¼ d=2) is much larger than the mean free path (λ) in the vapor [1,3,4], i.e., for sufficiently small Knudsen number (Kn ¼ λ=a).
For an isothermal drop below the critical temperature, the NSF equations without temperature-jump boundary conditions predict the time rate of change of the square of the droplet radius (or diameter) to be constant, also known as the d-squared (d 2 ) law of evaporation [5]. Interestingly, it has been observed in MD simulations [3,6,7] that when the Knudsen number Kn ¼ λ=a ≳ 1, the droplet radius evolves linearly in time. Notably, as will be shown in this Letter, the transition from d 2 to d law cannot be explained through the NSF equations with classical boundary conditions, i.e., with the Hertz-Knudsen-Schrage (HKS) relation [8] and an assumption that the temperature is continuous across the liquid-vapor interface.
In the literature, the transition from the d 2 to a d law is usually discussed in the MD framework [3,7] and explained theoretically by introducing ad hoc corrections due to rarefaction effects, offering very little insight into the macroscopic processes that dictate this changeover.
Rarefaction manifests itself through a temperature jump and kinetic boundary (Knudsen) layer, which have been observed experimentally [9][10][11] and predicted theoretically [12,13]. Notably, even at Kn ≈ 10 −3 -where the d-squared scaling is still seen-the NSF equations with classical boundary conditions are unable to give a good quantitative prediction of the total evaporation time of micrometer size droplets [9].
In this Letter, we show that the crossover from the d 2 to a d law is caused by a prominent temperature jump at the liquid-vapor interface. Using MD as a benchmark, we show that for nanosized droplets the NSF equations along with the temperature-jump boundary conditions, give good predictions for the evaporation process when the initial Knudsen number is below ≲0.5. In specific limits, analytic progress allows us to predict the transition from the d 2 to d law, with an explicit expression given between droplet radius and time.
Typically, there is a difference in velocity distribution function between molecules ejected from the condensed phase and those coming from the vapor, which causes a strong nonequilibrium in a thin layer adjacent to the interface, the Knudsen layer. The result is an actual jump in temperature at the liquid-vapor interface alongside a variation in temperature across the Knudsen layer, predicted accurately by the linearized Boltzmann equation (LBE) [14]. As the NSF equations cannot capture the Knudsen layer, the temperature jump boundary conditions include both the actual and apparent contributions to the jump [15,16]. In contrast, the regularized 26 moment (R26) equations [13] are able to approximate the Knudsen layer (and hence the temperature profile near the interface) so that the boundary condition does not include the apparent component.
Modeling the liquid drop.-Consider a spherical incompressible liquid droplet immersed in its own vapor. The droplet and the vapor flow surrounding the droplet retain spherical symmetry during the lifetime of the droplet, so that all the equations are in a spherical coordinate system with r being the radial distance from the center of the droplet. Let aðtÞ be the radius of the droplet at time t, and the temperature and pressure of the vapor at a distance far from the droplet be T ∞ and p ∞ , respectively. Throughout this Letter, properties relating to the liquid have the subscript l, and those relating to the vapor have the subscript v.
Owing to the conservation of mass, the velocity inside the droplet is zero, i.e., the radial velocity u l ¼ 0. The temperature T l ðt; rÞ inside the droplet is given by the energy balance equation where ρ l and c l are the density and specific heat of the liquid, respectively. The heat flux q l in the liquid is given by Fourier's law, i.e., q l ¼ −κ l ∂T l =∂r, where κ l is the thermal conductivity.
Modeling the surrounding vapor.-We consider a monoatomic ideal gas, assume the interface separating the liquid and vapor is infinitely thin, and the gradients in field variables (such as temperature) across the interface are modeled by setting the jump boundary conditions at the interface, see, for example, Refs. [8,13,16,17]. We study slow evaporation, which means small changes in the pressure and temperature from their equilibrium state p ∞ and T eq , where T eq is the saturation temperature at pressure p ∞ . Accordingly, only terms that are linear in deviations from the reference equilibrium state (p ∞ , T eq ) are considered.
The process in the vapor can be further simplified due to the high density ratio of the liquid phase to vapor phase. This ensures that the timescales for heat and mass diffusion within the vapor are very small compared to the liquid phase, which creates a quasisteady process in the vapor. The linearized conservation laws for mass and energy in the vapor are then At the interface (r ¼ a), mass and energy conservation give respectively, where H 0 is the specific heat of evaporation and j is the mass flux through interface. In equilibrium, both the mass flux j and the heat flux q v are zero, and the chemical potential and temperature are continuous across the liquid-vapor interface. These conditions lead to the Clausius-Clapeyron-Kelvin relation: Here, γ is the surface tension coefficient and Þ=RT eq is the saturation pressure for a planar surface with R being the specific gas constant.
In what follows, we consider the NSF equations with (i) classical boundary conditions and (ii) temperature-jump boundary conditions. Theories beyond NSF, in particular, the R26 equations [13], and the linearized Boltzmann equation (LBE) [14] are employed to investigate Knudsen layer effects and the results are compared with the MD simulations.
Model 1: NSF without jump.-When assuming equality of liquid and vapor temperatures at the interface, the HKS relation can be applied alongside NSF constitutive relations to give where ϑ is the condensation-evaporation coefficient, Model 2: NSF with temperature-jump.-When including temperature jumps, linear irreversible thermodynamics [16,17] gives that where T v ¼ T ∞ þ ðq v =κ v Þða 2 =rÞ is the temperature in vapor. In Eq. (5), the elements of the Onsager resistivities matrix r αβ are taken from Ref. [15] as These coefficients were obtained from kinetic theory by assuming that the evaporation or condensation coefficient does not depend on the impact energy and that all noncondensing vapor particles are being thermalized. Furthermore, for all numerical computations, we assume ϑ ¼ 1 (an assumption supported by the MD simulations [18]).
The results for ϑ ≠ 1, which are similar, are discussed in the Supplemental Material [19].
Comparison between models and MD results.-Equations (1) and (3) The property values, such as H 0 , T eq , and ρ l , are all evaluated for a Lennard-Jones gas as in Ref. [21], in order to compare the results with the MD results given in Ref. [5]. The full list of flow parameters are given in the Supplemental Material [19].
The MD simulation data (circle in Fig. 1) predict a short initial growth of the droplet, due to condensation of the surrounding hot vapor on the cold liquid drop. As shown by Fig. SM7 in the Supplemental Material [19], relaxing the assumption of quasistatic flow in the vapor and accounting for the finite size of the domain, this growth can be captured in our model. However, as this initial transient is essentially an artifact of the initial MD configuration, we do not focus on this behavior in the present work. From Fig. 1, the NSF equations with classical boundary conditions (dashed black line) predict a linear evolution for the radius squared, but this macroscopic model does not give a good prediction of the rate of evaporation (the slope of the curve) or the evaporation time, about half of MD predictions. On the other hand, the NSF equations with jump boundary conditions (solid blue lines) give a good prediction for the rate of evaporation and total lifetime of the droplet. Notably, in the MD simulations a departure from the d 2 law is observed when the Knudsen number  Fig. SM8 in the Supplemental Material [19]), in the inset of Fig. 1 the slope jda 2 =dtj vs a is compared for macroscopic theories and MD simulation [22].
Revisiting d 2 law for nanodroplets.-An explicit analytical relation between radius of the droplet and time can be obtained from the NSF equations with no-jump boundary conditions (4) after some simplifying assumptions. The MD studies show that the evaporation process of the droplet consists of essentially two stages [5,22]. During the first stage, heat flows from the hot surrounding vapor to the droplet, causing the droplet temperature to rise from the initial temperature T 0 l to T Ã l (≃T eq , see Fig. SM2 in the Supplemental Material [19]), and at the same time, the radius increases due to the condensation of hot vapor at the colder droplet surface. After this relatively negligible initial growth, the MD simulations in Ref. [5] show an almost isothermal behavior for the droplet. During this stage, heat supplied by the hot vapor q v contributes only to the evaporation, the mass flux is given by j ¼ −q v =H 0 . Hence, for model 1, Eqs. (3) and (4) yield From this analytic result (denoted by × in Fig. 1), which agrees well with numerical results, one can see clearly that a 2 evolves linearly in time. This is the well-known d 2 law of evaporation [4]. The total evaporation time t F given by the d 2 law (6) is t F ¼ ρ l H 0 a 2 0 =2κ v ðT ∞ − T eq Þ, which is smaller than predicted by the MD simulations [5], as shown in Fig. 1.
One can obtain an analytic solution analogous to Eq. (6) by taking the temperature jump into account, if one neglects the Kelvin's correction, which does not play a significant role on evaporation, see Fig. SM4 of the Supplemental Material [19]. From the jump boundary conditions (5), along with j ¼ −q v =H 0 , one can obtain a solution for q v , j, and T Ã l . Neglecting the Kelvin's correction, these equations become linear and yield KnÞ is the effective heat conductivity in the vapor. Substituting j in Eq. (3) and integrating, we obtain where the coefficient The temperature jump Utilizing Eq. (9), for a millimeter sized drop the temperature jump is seen to be negligible (δT ≈ 5 × 10 −4 K), for micrometer sized drops it becomes moderate (δT ≈ 0.5 K), while for a ten nanometer drop the jump is large (δT ≈ 40 K). Clearly from Eq. (7), when the droplet radius is significantly bigger than λ, the radius squared evolves linearly with time. However, when the droplet radius is much smaller than λ, a evolves linearly in time, a behavior that cannot be described without the second term on the left-hand side of Eq. (7), which comes from the jump boundary condition. In Fig. 1, the results from our numerical simulations are compared to Eq. (7); there is a close agreement. The total evaporation time t F obtained from Eq. (7) is given by Clearly, t F depends on α 0 and Kn 0 , and is larger than that predicted by the classical d-squared law (6). Figure SM3 in the Supplemental Material [19] presents t F as a function of Kn 0 for different values of the evaporation coefficient ϑ.
Again, a good agreement between numerical simulations and Eq. (10) is observed. For 0.5 ≤ ϑ ≤ 1, values typically obtained from MD simulations [18], α 0 and hence t F change within less than 10%. Temperature profiles.-The temperature vs the inverse distance (x ¼ a=r) is shown in Fig. 2 at t ¼ 20.92 and 32.63 ns. The interface is at x ¼ 1, with vapor on the left and liquid on the right. The MD results from Ref. [5] are denoted by circles when a ¼ 10.819 nm (at t ¼ 20.92 ns) and 9.327 nm (at t ¼ 32.63 ns), respectively. It is difficult to differentiate the temperature profiles in MD at these times, because of thermal fluctuations, so the same symbols are used to represent the temperature at all times. The results of the NSF with classical boundary conditions are depicted with black dashed lines. Note that this model predicts a ¼ 8.257 and 4.057 nm with average liquid temperature T l ¼ 105.3 and 107.6 K, respectively. Therefore, as well as substantial differences in drop radii, the model also fails to capture the temperature profile, by missing all MD data points in the vapor and overestimating the temperature in the liquid. While results from the NSF equations with jump boundary conditions (solid red lines) are consistent with MD predictions for the liquid temperature, the agreement between the NSF and MD for the vapor temperature is not satisfactory, particularly close to the interface.
Kinetic effects.-The NSF equations cannot produce correct temperature profiles in the vapor because of (a) their inability to capture the Knudsen layer and (b) the finite thickness of the interface. However, the Knudsen layer is significantly larger than the interfacial width, which is about 8.6 Å [21], so the assumption of a sharp interface is justified.
The LBE provides accurate solutions for the entire range of Knudsen numbers. However, the coefficients tabulated in Ref. [14] do not provide any information about the flow fields. To recover this information we apply the R26 equations, which provide a significant improvement over the NSF equations, accurately approximating LBE results up to Kn ≲ 1 [13] (See Fig. SM1 in Supplemental Material [19] for a comparison of the models). The analytic expressions for the R26 equations derived in Ref. [13] allow one to decompose the temperature profiles in the vapor into a classical Fourier contribution (dashed lines in Fig. 3) for t ¼ 20.92 ns, plus a contribution from the Knudsen layer. Notably, the R26 theory provides a better qualitative and quantitative agreement with the MD data compared to the NSF results.
Finally, we consider a case where kinetic effects dominate from the outset, with Kn 0 ¼ 1.33, as studied by Refs. [3,7] using MD simulations. In Fig. 4, we plot (a=a 0 ) vs time with initial droplet radius a 0 ¼ 2.3 nm and initial liquid temperature T 0 l ¼ 93 K. The far-field temperature T ∞ and pressure p ∞ are 500 K and 0.4 MPa, respectively. Indeed, the NSF with jump (solid line) and the PHYSICAL REVIEW LETTERS 123, 154501 (2019) LBE (dot-dashed lines) confirm the expected linear decay of the radius with time, which is also seen in the MD results (circles), until t ¼ 1.2 nm, at which point the droplet radius becomes comparable to the interfacial width [21] as well as the cutoff radius (2.5σ LJ ) used in the MD simulations.
The main advantage of continuum based models is that they can readily be incorporated in mutiscale multiphysics softwares used for engineering design optimization. In summary, it has been shown that to predict scaling transitions occurring when a ≈ λ the NSF must be supplemented with a temperature jump boundary condition, but to capture finer flow characteristics one must go beyond NSF. Motivated by these findings, future work could be to study mixtures, where concentration jumps will play a vital role, and to consider processes close to the critical point, where the Enskog-Vlasov formulation [23] allows one to capture nonideal effects in the vapor and gives a unified description of liquid and vapor.  [3] (circles). Symbols × and þ concur with the d 2 law (6) and Eq. (7), respectively.