Analytical Model for Diffusive Evaporation of ... - ACS Publications

May 14, 2018 - controlled droplet evaporation is conveniently solved in the toroidal coordinate system by applying the method of separation of variabl...
0 downloads 0 Views 3MB Size
Article Cite This: Langmuir XXXX, XXX, XXX−XXX

pubs.acs.org/Langmuir

Analytical Model for Diffusive Evaporation of Sessile Droplets Coupled with Interfacial Cooling Effect Tuan A. H. Nguyen,* Simon R. Biggs, and Anh V. Nguyen* School of Chemical Engineering, The University of Queensland, Brisbane, Queensland 4072, Australia S Supporting Information *

ABSTRACT: Current analytical models for sessile droplet evaporation do not consider the nonuniform temperature field within the droplet and can overpredict the evaporation by 20%. This deviation can be attributed to a significant temperature drop due to the release of the latent heat of evaporation along the air−liquid interface. We report, for the first time, an analytical solution of the sessile droplet evaporation coupled with this interfacial cooling effect. The two-way coupling model of the quasi-steady thermal diffusion within the droplet and the quasi-steady diffusioncontrolled droplet evaporation is conveniently solved in the toroidal coordinate system by applying the method of separation of variables. Our new analytical model for the coupled vapor concentration and temperature fields is in the closed form and is applicable for a full range of spherical-cap shape droplets of different contact angles and types of fluids. Our analytical results are uniquely quantified by a dimensionless evaporative cooling number Eo whose magnitude is determined only by the thermophysical properties of the liquid and the atmosphere. Accordingly, the larger the magnitude of Eo, the more significant the effect of the evaporative cooling, which results in stronger suppression on the evaporation rate. The classical isothermal model is recovered if the temperature gradient along the air−liquid interface is negligible (Eo = 0). For substrates with very high thermal conductivities (isothermal substrates), our analytical model predicts a reversal of temperature gradient along the droplet-free surface at a contact angle of 119°. Our findings pose interesting challenges but also guidance for experimental investigations.



INTRODUCTION Diffusive evaporation of a droplet is driven by the spatial gradient of vapor between the air−liquid interface and ambient air.1−4 For a steady-state diffusion-controlled process, the solution of the Laplace equation is widely used to predict the interfacial evaporation dynamics. Typically, the Laplace equation is solved either analytically or numerically together with the assumption that the atmosphere just above the air− liquid interface is saturated with vapor and the vapor saturation concentration is constant along the interface. Mathematically, there exists an analogy between these diffusive vapor concentration fields and the electrostatic potential fields around a spherical-cap conductor at constant potential because they are both governed by the Laplace equation. The analytical solution for the electrostatic potential of a charged conductor was first derived by Lebedev in 1965.5 Deegan et al.1,6 later used this solution to explain the nonuniform evaporative flux at the air− liquid interface of sessile droplets with acute contact angles. This classical analytical model does not take into account the effect of interfacial heat transfer, and temperature is considered uniform across the whole droplet. Thus, the saturated vapor © XXXX American Chemical Society

concentration along the air−liquid interface also does not change. This model is, therefore, called the “isothermal model”. Accurate measurements of droplet temperature, however, show a significant change in temperature inside the droplet and along the air−liquid interface.2,7,8 This nonuniform temperature field could be attributed to not only the nonuniform evaporative flux at the interface but also to the heat transfer process within the droplet. In particular, the higher the evaporative flux, the greater the latent heat of evaporation that is released, lowering the air−liquid interfacial temperature locally. Furthermore, the local temperature also varies because heat conduction length from the substrate to the air−liquid interface changes radially in accordance with the droplet height from the center to the edge. The combination of these heat and mass transfer processes leads to a nonuniform saturated vapor concentration along the air−liquid interface, different from the assumption of the “isothermal model”. Neglecting the cooling Received: November 8, 2017 Revised: May 8, 2018 Published: May 14, 2018 A

DOI: 10.1021/acs.langmuir.7b03862 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

equilibration of the vapor concentration field. Using similar arguments for the vapor phase, the criterion for attaining a quasi-steady state of heat transfer within a droplet is the ratio of heat equilibrium time in a droplet and the total evaporation time (droplet lifetime) theat/tf ≤ 0.1 (detailed discussion can be found in ref 15). Generally, the quasi-steady assumption is satisfied for a droplet if the averaged droplet temperature is smaller than ∼75, 38, 36, and 29 °C for water, isopropanol, ethanol, and methanol, respectively.16 The relative importance of convection heat transfer can be justified by the Péclet number (Pe), which is defined as the ratio of the rate of advection to thermal diffusion.15,17 It is reported that heat convection only becomes appreciable once Pe exceeds approximately 10 or so;15 thus, the assumption of neglecting thermal (Marangoni) convection inside a water droplet is valid for low Péclet number, say Pe < 10. Once this happens, the temperature field throughout the liquid droplet is determined by pure heat conduction. With above assumptions, mass conservation for vapor evaporation by diffusion is described by second Fick’s law. The scale analysis indicates that the evaporation of sessile droplets can be described by the quasi-steady state with neglecting the transient term in Fick’s law,1,4,5,18 giving ∇2C = 0. Heat transfer in a drying droplet can also be considered as quasi-steady state2 and is dominated by the heat conduction over the convection term.19 Therefore, the temperature field within a slowly evaporating droplet can be described by the Laplace equation, ∇2T = 0. These Laplace equations for C and T can be conveniently solved for the vapor concentration and temperature fields by applying the method of separation of variables in the toroidal coordinate system (α, β, φ)4,5 (Figure 1). Accordingly, the solution for vapor concentration is

effect and decoupling the thermal transport mechanism in the isothermal model may lead to considerable discrepancy in predicting the droplet evaporation dynamics and associated phenomena. For instance, Dash and Garimella9 reported that the evaporation rate of a 160° droplet is overpredicted by ∼20% when using the isothermal model. This discrepancy was attributed to a large temperature drop along the air−liquid interface due to the cooling effect and was confirmed by the numerical simulations.8,10,11 Chandramohan et al.7 recently reported a large temperature gradient due to the evaporative cooling effect along the air−liquid interface of a sessile droplet on a copper substrate. An improved model for the droplet evaporation, therefore, needs to allow the change in saturation vapor concentration along the air−liquid interface as a function of temperature. This “nonisothermal model” has also to be coupled with the equation of heat transfer inside the droplet and solved in an iterative manner. This coupled problem of diffusive vapor and heat conduction has been solved numerically by Dunn et al.,12 Xu and Ma.8 Pan et al.10,11 also carried out a numerical modeling to study the competing effects of external natural convection and the evaporative cooling. In another attempt to explain the discrepancy between experimental and theoretical results, Gleason and Putnam13 used the analytical solution of the isothermal model with replacing the vapor concentration along the air−liquid interface by a function of temperature. Hu and Larson2 used a finite-element method to compute the temperature field inside evaporating droplets whose boundary condition at the air−liquid interface was the analytical expression for the evaporative flux from the isothermal model. This one-way coupling model also shows that a nonuniform evaporation leads to a nonuniform distribution of temperature along the air−liquid interface. Despite these efforts, the complete theory for the effects of the interface cooling on the evaporation of liquid droplets is still lacking.



THEORETICAL ANALYSIS In this paper, we report, for the first time, closed-form, analytical solutions of vapor diffusion from a sessile droplet into the surrounding air by taking into account the interfacial cooling. The analytical solutions can be used for a full range of small spherical-cap droplets of different contact angles and types of fluids, resting on a flat substrate heated at constant temperature Tsub. This assumption of an isothermal surface is applicable for substrates made of materials having high thermal conductivity, such as copper (398 W/m K), gold (315 W/m K), aluminum (247 W/m K),14 and so forth. For simplicity, we ignore the effects of varying the substrate temperature, natural vapor convection, and convective heat transfer within the droplet. That being said, the present model is a quasi-steady diffusion-controlled droplet evaporation model. As discussed by Larson,15 the validity of the quasi-steady-state assumption for vapor-phase mass transfer is governed by the ratio of the characteristic time for the vapor field to adjust to the change in droplet shape to the characteristic time of droplet shrinkage (which is just the drying time), tvap/tf ≤ 5 × 10−3. This ratio can be approximated as five times the ratio of mass density of vapor and mass density of droplet liquid (∼5ρvap/ρL) which is quite small; thus, the vapor concentration above the surface of the droplet reaches a quasi-steady state well before the completion of drying. Although the vapor phase quasi-steady-state assumption is normally very well justified, thermal equilibration throughout the liquid droplet is much less rapid than

Figure 1. Schematic of a sessile droplet (spherical cap) evaporating from a flat substrate in rotationally symmetric cylindrical (r, z)and toroidal (α, β)coordinates. These coordinates are linked through complex mapping z + ir = iR coth(α + iβ/2), where R is the droplet base radius and i = − 1 . The substrate is heated at constant temperature Tsub ≥ T∞, which triggers the thermal conduction within the internal domain of π − θ ≤ β ≤ π, where θ is the contact angle of the liquid phase. Vapor diffuses from the droplet-free surface through the ambient air confined within the external domain of 2π ≤ β ≤ 3π − θ.

bounded by the physical domain: 0 ≤ α < ∞ and 2π ≤ β ≤ 3π − θ and satisfies the boundary condition C(0, 2π) = C∞ at infinity of ambient air at temperature T∞ (T∞ ≤ Tsub), (∂C/ ∂α)α=0 at the axis of symmetry and the zero flux of vapor diffusion (∂C/∂β)β=2π = 0 at the solid−vapour interface. On the other hand, the solution for temperature is bounded by 0 ≤ α < B

DOI: 10.1021/acs.langmuir.7b03862 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

Table 1. Magnitude of the Dimensionless Number (Eo) Characterizing the Interfacial Cooling Effect for Different Liquids at T = 295 K and p = 99.8 kPaa liquid

D [m2 s−1]

L [m2 s−2]

k [kg m s−3 K−1]

b [kg m−3 K−1]

Eo [−]

−5

2.45 × 10 1.20 × 106 5.49 × 105

0.604 0.203 0.161

1.11 × 10−3 9.47 × 10−3 2.84 × 10−2

0.11 0.84 1.03

2.44 × 10 1.50 × 10−5 1.06 × 10−5

water methanol acetone

6

a

D is the coefficient of vapor diffusion in air, L is the liquid latent heat of vaporization, k is the liquid thermal conductivity, and b is the thermal gradient of vapor saturation concentration at T.

∞ and π − θ ≤ β ≤ π and satisfies the boundary conditions of T(α, π) = Tsub at the droplet base and (∂T/∂α)a=0 = 0 at the axis of symmetry, giving C ̃ (α , β ) =

∫0

1 dF 3 dθ

{(τF coth θτ −

2 cosh α − 2 cos β

EC(τ )Piτ − 1/2(cosh α) cosh[(2π − β)τ ] dτ T (α , β) − Tsub = Tsub − T∞

(1)

ET =

1 dF 3 dθ

o

1 dF 3 dθ

)}

1 dF 3 dθ

1 dF 3 dθ



E T(τ )Piτ − 1/2(cosh α) sinh[(π − β)τ ] dτ

o

1 dF 3 dθ

)} (6)

(2)

where F(θ, τ) = sinh θτ/(sinh πτ sin θ) and dF/dθ = (τ cosh θτ sin θ − sinh θτ cos θ)/(sinh πτ sin2 θ). Eo = bLD/k is an evaporative cooling number and comes from the relationship between the fluxes at the droplet surface q(α) = Lj(α), which can be represented in the nondimensional form as ⎯→ ⎯ ⎯→ ⎯ Eo(∇C̃ · iβ )β = 3π − θ = (∇T̃ · iβ )β = π − θ . With the given mapping C̃ = 1 + T̃ along the droplet surface, the left-hand side of the above equation can be considered as a “quasi-temperature” (C̃ ) field across the droplet surface with a “thermal” conductivity Eo, whereas the right-hand side is an interfacial diffusive heat transfer with a thermal conductivity of 1. Hence, the evaporative cooling number Eo could be interpreted as the ratio between the interfacial diffusive and conductive heat transfer rates across the droplet surface, and in a broad sense, Eo could possibly be assigned as the Nusselt number. Generally, Eo characterizes the strength of the evaporative cooling effect on droplet evaporation. Its magnitude is determined by the ratio of the interfacial diffusive to conductive heat transfer rates and depends on the diffusion coefficient of liquid vapor in the atmosphere, the temperature dependence of the vapor saturation concentration, and the thermal properties of the fluid12,20−22 (Table 1). The larger the magnitude of Eo, the more significant the effect of evaporative cooling. When the effect of evaporative cooling can be neglected (i.e., Eo = 0), the problem of vapor diffusion in the atmosphere is decoupled from the problem of the heat transfer in the droplet and the present model reduces to the classical “isothermal model”23,24

⎯→ ⎯

⎯→ ⎯

Eo cosh θτ cosh πτ sinh θτ

{τF tanh[(θ − π)τ] − } {(τF coth θτ − ) − E (τF tanh[(θ − π)τ] −

2 cosh α − 2 cos β

j(α) = (D∇C· iβ )β = 3π − θ and q(α) = (k∇T · iβ )β = π − θ , yielding

q(α) =

{τF coth θτ − } ) − E (τF tanh[(θ − π )τ ] −

(5)

where Ce = Csat(Te) is the saturated vapor concentration of the liquid at temperature Te of the droplet edge (α → ∞ and Te = Tsub), τ is the integration dummy, and Piτ−1/2(cosh α) is the toroidal function (i.e., the first-kind Legendre function of the complex half-integral degree and the argument of the hyperbolic cosine function). These solutions are independent of φ because of the rotational symmetry. EC(τ) and ET(τ) in eqs 1 and 2, respectively, are functions of the integration dummy (independent of the toroidal coordinates α and β), which can be determined from the boundary conditions at the air−liquid interface as follows. The net of the heat and mass fluxes at the droplet surface depends on the toroidal coordinate and are determined by

j(α) =

cosh θτ cosh πτ cosh[(θ − π )τ ]



T ̃ (α , β ) =

∫0

C(α , β) − C∞ = Ce − C∞

EC =

∞ Ce − C∞ 2 (cosh α + cos θ)3/2 × EcPiτ − 1/2(cosh α) 0 R /D ⎡ cosh[(θ − π )τ ] sin θ ⎤ − τ sinh[(θ − π )τ ]⎥ dτ ⎢ ⎣ 2(cosh α + cos θ) ⎦ (3)



∞ Tsub − T∞ 2 (cosh α + cos θ)3/2 E TPiτ − 1/2(cosh α) 0 R /k ⎡ sinh θτ sin θ ⎤ − τ cosh θτ ⎥ dτ ⎢ ⎣ 2(cosh α + cos θ) ⎦ (4)



where k is the liquid thermal conductivity and D is the vapor diffusion coefficient. The local energy balance on the air−liquid interface establishes a relationship between the evaporative mass flux and the heat flux by q(α) = Lj(α), where L is the liquid latent heat of vaporization. Additionally, the atmosphere just above the air−liquid interface of the droplet is considered saturated with vapor, whose temperature dependence can be linearized locally as follows:8 C(α, 3π − θ) = Cs(T) = a + bT, where b = dCsat/dT = (Ce − C∞)/(Tsub − T∞) and a = Ce − bTsub. Thus, C̃ (α, β) = 1 + T̃ (α, β) along the droplet surface. Given these conditions, EC(τ) and ET(τ) can now be determined as follows (with detailed derivations in the Supporting Information)

C̃ iso(α , β) =

∫0

Piτ − 1/2(cosh α)

jiso (α) =

∫0

2 cosh α − 2 cos β



cosh θτ cosh[(2π − β)τ ] dτ cosh πτ cosh[(π − θ)τ ]

Cs − C∞ ⎧ sin θ ⎨ + R /D ⎩ 2



Piτ − 1/2(cosh α)

(7)

2 (cosh α + cos θ)3/2

τ cosh θτ tanh[(π − θ)τ ] ⎫ dτ ⎬ ⎭ cosh πτ (8)

C

DOI: 10.1021/acs.langmuir.7b03862 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir and T̃ iso = 0. In these equations, the subscript “iso” stands for the isothermal condition at the air−liquid interface.

the vapor concentration and temperature fields. For example, Figure 3 shows the theoretical predictions for the concentration contours of the vapor field above the air−liquid interface of a droplet with a contact angle of 60° using eqs 1 and 5. When Eo = 0, the current model reduces to the classical model of isothermal; thus, the saturated vapor concentration just above the free surface of the droplet is constant. Closed contour lines are concentrated near the free interface of the droplet and gradually diluted outward in a uniform pattern. However, when the interfacial cooling effect is coupled (i.e., Eo > 0), the vapor concentration near the free interface gradually decreases from the droplet apex toward the contact line. Thus, the isoconcentration lines near the free surface are not in closed contours. The temperature field inside a drying droplet is presented in Figure 4 in terms of isotherms calculated from eq 2. A cooler zone near the interface at the droplet apex is predicted in comparison with the zone near the droplet edge. This cold zone is corresponded to the low saturated vapor concentration zone observed in Figure 3. This can be explained by the heat conduction distance between the heated substrate and the air− liquid interface which becomes greater from the droplet edge toward the center (with acute contact angles). In other words, a relatively tall droplet (large contact angle) has a large effective thermal resistance between the substrate and the air−liquid interface. The higher the value of Eo used in eq 2, the stronger the cooling effect. Temperature Profile along the Air−Liquid Interface Due to the Cooling Effect. To visualize the cooling effect in more detail, the temperature profile along the air−liquid interface is plotted in Figure 5 for droplets of different contact angles and evaporative cooling numbers. The interfacial temperature profile is found to be dependent on the contact angle of the droplet and on the evaporative cooling number Eo. When the contact angle is smaller than 90° (Figure 5A), the interfacial temperature decreases monotonically from the droplet edge to the droplet center. The greatest interfacial temperature drop observed at the droplet center is due to the longest heat conduction distance between the heated substrate and the air−liquid interface there. Therefore, the larger the contact angle, the greater the temperature drop can be. Together with the droplet geometry, fluid type also has a significant influence on the interfacial temperature profile: a



RESULTS AND DISCUSSION Model Validation. Figure 2 shows the comparison between our analytical results and the published numerical results8 for

Figure 2. Comparison between the analytical solution obtained from this study and the numerical results presented in ref 8. Evaporative flux along the air−liquid interface is normalized by dividing by D(Ce − C∞)/R); the normalized interfacial temperature T̃ (α, π − θ) is given by eqs 2 and 6. The contact angle between the droplet and the substrate is 20°.

evaporative mass flux and temperature profile along the droplet surface. The comparison shows good agreement for different levels of cooling effect, that is, Eo = 0.5, 2.0, and 10. Additionally, the classical isothermal model is recovered from our new model if the temperature gradient along the air−liquid interface is negligible, Eo = 0 (eqs 7 and 8, red lines in Figure 2). These numerical results were already validated against available experimental evaporation data of sessile droplets resting on a highly thermal conductive aluminum substrate. Thus, these agreements verify the validity of our analytical model for the sessile droplet evaporation coupled with the interfacial cooling effect. It is worth noting that the advancement of our analytical coupling model is its applicability for a full range of droplets with different contact angles and types of liquids. Spatial Distribution of Vapor and Heat Fields. The close-form analytical solutions in eqs 1 and 2 now can be used to predict and investigate the effect of evaporative cooling on

Figure 3. Contour plots of vapor distribution above an evaporating droplet, θ = 60°, with (Eo = 0.5, 2 and 5) and without (Eo = 0) cooling effect. The color bar represents the normalized vapor concentration (C − C∞)/(Ce − C∞). D

DOI: 10.1021/acs.langmuir.7b03862 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

Figure 4. Contour plots of the temperature field within an evaporating droplet, θ = 60°, with (Eo = 0.5, 2 and 5) and without (Eo = 0) cooling effect. The color bar represents the normalized temperature (T − Tsub)/(Tsub − T∞).

Figure 5. Temperature distribution along the air−liquid interface of evaporating sessile droplets of different contact angles and evaporative cooling effect. Temperatures are calculated by eqs 2 and 6 and normalized as (T − Tsub)/(Tsub − T∞). The inserted graph in (B) is T̃ (α, π − θ) plotted at Eo = 1.0 with contact angles near the critical value (θ ≈ 119°).

numerical calculations8,10,26 have confirmed the interfacial cooling effect by which the interfacial temperature gradually cools down toward the droplet apex, the inversion of the interfacial cooling (iii) around θ ≈ 119° has not been reported before for isothermal substrates. It is probably because of the nonuniform vapor field that could lead to local readsorption of the vapor phase to the droplet interface (negative local evaporative flux). With the sake of the full analytical solutions (eqs 2 and 6), our prediction for a possible reversal of temperature gradient at ∼119° will pose important challenges for experimentalists to probe and confirm its validity. For the droplet evaporation from nonisothermal substrates with a finite thickness, different studies came up with different critical contact angles, for example, 9°,29 14°,2 and 31°,19 at which the reversal of temperature gradient occurs. To clarify these critical values, we need to develop a coupling model for nonisothermal substrate in the first instance. We will address this interesting development in a forthcoming publication. It is also worth noting that the numerical integration of eqs 1 and 2 is highly oscillated for large contact angles, that is, θ ≈ 160°; thus, a dedicated study is still needed to further explore this phenomenon for droplets drying from highly hydrophobic (superhydrophobic) surfaces. Additionally, for those droplets with large contact angles, the convective vortex motion inside the droplets could be substantial and should not be completely ignored. Because the heights of the droplet (height-to-radius aspect ratio) grow in time and the fluid velocities are greater (due to the thermal Marangoni effect) than those in the isothermal cases, the numerator of the Pe number (as defined and discussed in the previous sections) will increase. Thus, the advection due to the thermal Marangoni effect could play a

higher evaporative cooling number Eo causes a greater temperature drop. For droplets with large contact angles (on hydrophobic surfaces, Figure 5B), on the other hand, the evaporative cooling effect leads to more complex profiles of the interfacial temperature. Especially, near the droplet waist and the contact line region, there is a counteracting/combining effect of the local geometric confinement for mass transfer and the short distance for heat conduction from the substrate. Similar to the case of droplets with acute contact angles, the fluid near the contact line region experiences high local temperature because of the short heat conduction path from the substrate and thus supports high local saturated vapor concentration at the air− liquid interface. However, because of the geometry of obtuse contact angle droplets, the confined space between the waist and the substrate could lead to a nonuniform gradient of vapor concentration (i.e., nonclosed isoconcentration contours) near the droplet surface (similar to Figure 3 with Eo ≠ 0). The combination of these two effects leads to interesting predictions as shown in Figure 5B: (i) there is a narrow range of contact angle θ ≈ 119° (red lines in the inset of Figures 5B and S2) in which temperatures near the contact line and the droplet apex are almost similar; (ii) for θ < 119°, the trend of interfacial temperature profile is similar to that of acute droplets (the contact line region has the highest local temperature); and (iii) for θ ≥ 120°, the droplet surface near the contact line is colder than that at the droplet apex. For substrates with very high thermal conductivities (close to the isothermal condition used in our study), previous studies have shown that there is no inversion of the temperature gradient along the air−liquid interface in the range of 0 < θ < π/2. It agrees with our findings presented in Figure 5A. Although experimental results7,16,25 and E

DOI: 10.1021/acs.langmuir.7b03862 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

Figure 6. (A,B) Variation of normalized evaporative fluxes Jiso (dashed lines, eq 8) and J (eq 9) along the droplet surface on both hydrophilic (θ = 30°, 50°, 70°, 80°, and 90°) and hydrophobic (θ = 100°, 105°, 110°, 115°, and 120°) substrates. (C,D) Differences in the diffusive mass fluxes (ΔJ ≡ Jiso − J, eq 11) showing evaporation suppression due to the interfacial cooling effect (Eo = 0 and 1.0).

Figure 7. Variation of nondimensional evaporative flux (normalized by dividing by D(Ce − C∞)/R) along the droplet surface corresponding to different evaporative cooling numbers Eo = 0.0 (classical model), 0.1, 0.3, 0.5, 1.0, 2.0, and 5.0. (Left) θ = 60° and (right) θ = 120°.

dominant role in the heat transfer within an evaporating

J3 =



droplet, especially for organic liquid droplets.17 Evaporative Flux with the Effect of Interfacial Cool-



ing. Using the Mehler−Fock integral transform and after some

0

Ce − C∞ sin θ Cs̃ (α) + R /D 2

{

∫0

2 (cosh α + cos θ)3/2



ECPiτ − 1/2(cosh α)τ sinh[(π − θ)τ ] dτ

}

⎧ ⎪ ⎨ ⎪ ⎩

(

Eo τF tanh[(θ − π )τ ] −

1 dF 3 dθ

)

{(τF coth θτ − ) − E (τF tanh[(θ − π)τ] − 1 dF 3 dθ

o

⎫ τ cosh θτ tanh[(π − θ)τ ] ⎪ ⎬ dτ × Piτ − 1/2(cosh α) cosh πτ ⎪ ⎭

algebra, eq 3 can be displayed in a familiar form as follows j(α ) =

2 (cosh α + cos θ)3/2

1 dF 3 dθ

)}

(10)

Therefore, the evaporative cooling effect modifies the evaporative flux as

(9)

ΔJ = Jiso − J = (1 − Cs̃ )

where C̃ s = C̃ (α, 3π − θ) is the saturated vapor along the

sin θ − J3 2

(11)

Because C̃ (α, β) = 1 + T̃ (α,β) along the droplet surface, differences in the evaporative fluxes can also be represented as a function of the local air−liquid interface temperature: ΔJ = −0.5T̃ (α, π − θ) sin θ − J3. In the case of acute contact angles, both models show that the evaporative mass flux increases from the droplet surface center (r = 0) to the droplet edge (r = R), where the flux goes to infinity (Figure 6A). As previously commented,27 the singularity of the evaporative flux is due to the difference between molecules diffusing away from and back to the liquid phase per unit time. At the edge of the droplet with an acute contact angle, the probability of molecules diffusing back to the

droplet-free surface. The classical (isothermal) model, eq 8, predicts that the normalized evaporative flux (by dividing by D(Ce − C∞)/R) is a sum of two terms Jiso = J1 + J2, whereas the analytical model developed in this paper, on the other hand, predicts it as a sum of three terms in the evaporative flux equation, that is, J = C̃ sJ1 + J2 + J3. In these equations, J1 and J2 are the first and the second term in the bracket of eq 8, respectively, and J3 is a term associated with Eo extracted from eq 9 F

DOI: 10.1021/acs.langmuir.7b03862 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

concentration on the temperature, our analytical coupling model gives good agreement with previous numerical results on evaporation of sessile droplets resting on a highly thermal conductive substrate. For fast drying droplets on heated substrates, however, more realistic expressions for the variation of vapor concentration on temperature (higher order of polynomials) should be used instead. Furthermore, when the substrate is not isothermal, that is, poorly conductive substrates, a three-domain coupling solution needs to be developed based on the theoretical framework presented in this work.

droplet interface decreases and consequently the evaporative flux by diffusion increases toward the droplet edge. Additionally, the present model described by eq 3 generally gives lower predictions for the local surface flux and is dependent on the intensity of the cooling effect Eo (Figure 7). As can be seen from the Figure 6C, the difference between the isothermal mass flux and the nonisothermal one is uniform across the droplet surface except near the contact line area. ΔJ quickly decreases toward zero when r/R > 0.8 and becomes negative as r/R ≈ 1.0 (eq 9). The discrepancy ΔJ becomes even more apparent at higher contact angle and higher evaporative cooling number. On hydrophobic substrates, the isothermal model of eq 8 shows that the mass flux starts decreasing near the droplet waist and asymptotically approaches zero at the contact line (Figure 6B). That is because the local geometric confinement retards the vapor diffusion to ambient. However, the isothermal model fails to predict the trend of local evaporative flux on heated substrates. Predictions from our present model show that the fluid near the contact line region experiences the highest local temperature (θ < 119°, Figure 5) and thus supports the highest local saturated vapor concentration at the air−liquid interface. Consequently, the evaporative mass flux is the highest near the contact line as presented in Figure 6B. It is also worth noting that the geometric confinement effect mentioned above still plays a role in the current model. Evidently, for a small range of contact angle, that is, 105 < θ < 120°, our model predicts a special profile of the evaporative mass flux: J starts decreasing near the droplet waist (r/R > 1.0) and then rapidly increases in the proximity of the droplet edge (r → R). For a droplet of higher contact angles, on the other hand, both models predict that the flux asymptotically approaches zero at the droplet edge (r = R) as shown in Figures 6B and 7 for θ = 120°. Equation 11 is also used to quantify the reduction in the evaporation flux due to the interfacial cooling effect for droplets on hydrophobic surfaces. This reduction is similar to the case of acute contact angles and can be visualized from Figure 6D (θ ≤ 120°). Evidently, these significant changes in the local flux are predicted as a combination of the third (new) term associated with Eo as described by eq 10 and the distribution of the interfacial temperature T̃ (α, π − θ). However, on the hydrophobic surfaces, the latter term seems to play a dominant role28 (Figure S3).



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.langmuir.7b03862. Coordinate system; governing equations and solutions for mass and heat transfer; heat and mass fluxes at the droplet surface; and coupling conditions at the droplet surface (PDF)



AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (T.A.H.N.). *E-mail: [email protected] (A.V.N.). ORCID

Tuan A. H. Nguyen: 0000-0002-3024-2081 Anh V. Nguyen: 0000-0001-6703-2291 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We would like to thank Professor Suresh Garimella and reviewers for their helpful discussions on this research.



REFERENCES

(1) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R.; Witten, T. A. Capillary flow as the cause of ring stains from dried liquid drops. Nature 1997, 389, 827−829. (2) Hu, H.; Larson, R. G. Analysis of the Effects of Marangoni Stresses on the Microflow in an Evaporating Sessile Droplet. Langmuir 2005, 21, 3972−3980. (3) Nguyen, T. A. H.; Nguyen, A. V. Transient Volume of Evaporating Sessile Droplets: 2/3, 1/1, or Another Power Law? Langmuir 2014, 30, 6544−6547. (4) Nguyen, T. A. H.; Nguyen, A. V.; Hampton, M. A.; Xu, Z. P.; Huang, L.; Rudolph, V. Theoretical and experimental analysis of droplet evaporation on solid surfaces. Chem. Eng. Sci. 2012, 69, 522− 529. (5) Lebedev, N. N. Special functions and their applications, Revised English ed.; Prentice-Hall: Englewood Cliffs, 1965; Chapters 7 and 8. (6) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R.; Witten, T. A. Contact line deposits in an evaporating drop. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2000, 62, 756−765. (7) Chandramohan, A.; Weibel, J. A.; Garimella, S. V. Spatiotemporal infrared measurement of interface temperatures during water droplet evaporation on a nonwetting substrate. Appl. Phys. Lett. 2017, 110, 041605. (8) Xu, X.; Ma, L. Analysis of the effects of evaporative cooling on the evaporation of liquid droplets using a combined field approach. Sci. Rep. 2015, 5, 8614.



CONCLUSIONS To summarize, the process of diffusion-driven evaporation of sessile droplets on a solid planar surface has been reconsidered and modeled by taking into account the interfacial cooling effect. This nonisothermal vapor diffusion model is solved with the physically consistent boundary condition for saturated vapor density at the droplet-free surface and in coupling with the thermal diffusion within the liquid phase of the droplet. Our new analytical model is in closed form and is applicable for a full range of spherical-cap shape droplets of different contact angles and types of fluids. Accordingly, the intensity of the evaporative cooling effect is quantified by the dimensionless number, Eo, whose magnitude is the combination of effects of the diffusion coefficient of liquid vapor in the atmosphere, the temperature gradient of the vapor saturation concentration, and the thermal properties of the fluid. The larger the magnitude of Eo, the more significant the effect of the evaporative cooling effect. The classical isothermal model is recovered from the new predictions if the interface cooling is negligible, Eo = 0. Despite the use of a simple linear dependence of the saturated vapor G

DOI: 10.1021/acs.langmuir.7b03862 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir (9) Dash, S.; Garimella, S. V. Droplet Evaporation Dynamics on a Superhydrophobic Surface with Negligible Hysteresis. Langmuir 2013, 29, 10785−10795. (10) Pan, Z.; Dash, S.; Weibel, J. A.; Garimella, S. V. Assessment of Water Droplet Evaporation Mechanisms on Hydrophobic and Superhydrophobic Substrates. Langmuir 2013, 29, 15831−15841. (11) Pan, Z.; Weibel, J. A.; Garimella, S. V. Influence of Surface Wettability on Transport Mechanisms Governing Water Droplet Evaporation. Langmuir 2014, 30, 9726−9730. (12) Dunn, G. J.; Wilson, S. K.; Duffy, B. R.; David, S.; Sefiane, K. A mathematical model for the evaporation of a thin sessile liquid droplet: Comparison between experiment and theory. Colloids Surf., A 2008, 323, 50−55. (13) Gleason, K.; Putnam, S. A. Microdroplet Evaporation with a Forced Pinned Contact Line. Langmuir 2014, 30, 10548−10555. (14) Chung, D. D. L. Materials for thermal conduction. Appl. Therm. Eng. 2001, 21, 1593−1605. (15) Larson, R. G. Transport and deposition patterns in drying sessile droplets. AIChE J. 2014, 60, 1538−1571. (16) Kumar, M.; Bhardwaj, R. A combined computational and experimental investigation on evaporation of a sessile water droplet on a heated hydrophilic substrate. Int. J. Heat Mass Transfer 2018, 122, 1223−1238. (17) Chandramohan, A.; Dash, S.; Weibel, J. A.; Chen, X.; Garimella, S. V. Marangoni Convection in Evaporating Organic Liquid Droplets on a Nonwetting Substrate. Langmuir 2016, 32, 4729−4735. (18) Popov, Y. O. Evaporative deposition patterns: Spatial dimensions of the deposit. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2005, 71, 036313. (19) Ristenpart, W. D.; Kim, P. G.; Domingues, C.; Wan, J.; Stone, H. A. Influence of Substrate Conductivity on Circulation Reversal in Evaporating Drops. Phys. Rev. Lett. 2007, 99, 234502. (20) Perry, R. H.; Green, D. W.; Maloney, J. O. Perry’s Chemical Engineers’ Handbook, 7th ed; McGRAW-Hill, 1997. (21) Poling, B. E.; Prausnitz, J. M.; O’Connell, J. P. The Properties of Gases and Liquids, 5th ed.; McGRAW-Hill, 2001. (22) Raznjevic, K. Handbook of Thermodynamic Tables; Begell House, 1995. (23) Nguyen, T. A. H.; Nguyen, A. V. On the Lifetime of Evaporating Sessile Droplets. Langmuir 2012, 28, 1924−1930. (24) Nguyen, T. A. H.; Nguyen, A. V. Increased Evaporation Kinetics of Sessile Droplets by Using Nanoparticles. Langmuir 2012, 28, 16725−16728. (25) Patil, N. D.; Bange, P. G.; Bhardwaj, R.; Sharma, A. Effects of Substrate Heating and Wettability on Evaporation Dynamics and Deposition Patterns for a Sessile Water Droplet Containing Colloidal Particles. Langmuir 2016, 32, 11958−11972. (26) Dash, S.; Garimella, S. V. Droplet evaporation on heated hydrophobic and superhydrophobic surfaces. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2014, 89, 042402. (27) Poulard, C.; Guéna, G.; Cazabat, A. M. Diffusion-driven evaporation of sessile drops. J. Phys.: Condens. Matter 2005, 17, S4213. (28) Tran, H. V.; Nguyen, T. A. H.; Biggs, S. R.; Nguyen, A. V. On the predictions for diffusion-driven evaporation of sessile droplets with interface cooling. Chem. Eng. Sci. 2018, 177, 417−421. (29) Xu, X.; Luo, J.; Guo, D. Criterion for Reversal of Thermal Marangoni Flow in Drying Drops. Langmuir 2010, 26, 1918−1922.

H

DOI: 10.1021/acs.langmuir.7b03862 Langmuir XXXX, XXX, XXX−XXX