Evaporating Droplets - Langmuir (ACS Publications)

Evaporating Droplets - Langmuir (ACS Publications)pubs.acs.org/doi/abs/10.1021/la0204646CachedSimilarby M Cachile - ‎2...

1 downloads 46 Views 191KB Size


Langmuir 2002, 18, 8070-8078

Evaporating Droplets M Cachile,† O. Be´nichou, C. Poulard, and A. M. Cazabat* Laboratoire de Physique de la Matie` re Condense´ e, Colle` ge de France, 11 Place Marcelin Berthelot, 75231 Paris Cedex 05, France Received May 18, 2002. In Final Form: July 22, 2002 The dynamics of evaporation of drops of simple, completely wetting liquids is investigated experimentally and theoretically in the situation where no pinning of the contact line is occurring. When the contact angle is low enough, it is shown that a convenient approximation for the equation governing the dynamics of the drop is to keep only the evaporation term. The simplified equation has been studied mathematically and numerically. The agreement with experimental data is extremely good, allowing us to address more complex cases involving gradient and hydrodynamic effects.

I. Introduction Experimental and theoretical studies of liquids evaporating on heated or isothermal substrates have been recently performed by many authors.1-3,5,6,9-11,15-17,20-22,24-27 Anderson and Davis1 have shown that the contact angle of a drop put on a heated substrate, and continuously fed to keep a constant volume, is increased, in agreement with analysis by Hocking16 and Morris.20 A similar result is obtained by Garoff and co-workers with dip-coated films of evaporating liquids.24 Deegan9-11 and Parisse21 have more specifically studied drops of water or colloidal dispersions evaporating on isothermal substrates, when the contact line is pinned and evaporation proceeds through a decrease of the contact angle. Mixed situations, * Corresponding author. E-mail: [email protected] † Present address: Grupo de Medios Porosos, Facultad de Ingenierıa - UBA, Paseo Colo´n 850 (1063), Buenos Aires, Argentina. (1) Anderson, D. M.; Davis, S. H. Phys. Fluids 1995, 7, 248. (2) Betterton, M.; Brenner, M.; Stone, H. 2002. Preprint, in preparation. (3) Birdi, K. S.; Vu, D. T.; Winter, A. J. Phys. Chem. 1989, 93, 3702. (4) Blake, T. D. AIChE Spring meeting, New Orleans, LA, 1988; Paper I.a. (5) Blossey, R.; Bosio, A. Langmuir 2002, 18, 2952. (6) Bourge`s-Monnier, C.; Shanahan, M. E. R. Langmuir 1995, 11. 2820. (7) Cachile, M.; Be´nichou, O.; Cazabat, A. M. 2002. Langmuir 2002, in press. (8) Crank, J. The mathematics of Diffusion, 2nd ed.; Oxford science publications: Oxford, England, 1975. (9) Deegan, D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R.; Witten, T. A. Nature 1997, 389 827. (10) Deegan, D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R.; Witten, T. A. Phys. Rev. E 2000, 62, 757. (11) Deegan, R. D. Phys. Rev. E 2000, 61, 475. (12) Derjaguin, B. V.; Churaev, N. V. Wetting films; Nauka: Moscow, 1984. (13) Derjaguin, B. V.; Churaev, N. V.; Muller, V. M. Surface forces; Consultant bureau: New York and London, 1987. (14) Derjaguin, B. V.; Churaev, N. V.; Muller, V. M. Adv. Colloid Interface Sci. 1988, 28, 197. (15) Ehrard, P.; Davis, S. H. J. Fluid Mech. 1991, 229, 365. (16) Hocking, L. M. Phys. Fluids 1995, 7, 2950. (17) Hu, H.; Larson, R. G. J. Phys. Chem. 2002, 106. 1334. (18) Jackson, J. D. Classical Electrodynamics, 2nd ed.; Wiley: New York, 1975. (19) Lyklema, J. Fundamentals of Interface and Colloid Science; Academic Press: London, 1991; Vol. I. (20) Morris, S. J. S. J. Fluid. Mech. 2001, 432, 1. (21) Parisse, F.; Allain, C. Langmuir 1996, 13, 3598. (22) Picknett, R. G.; Bexon, R. J. Colloid Int. Sci. 1977, 61, 336. (23) Pomeau, Y. C. R. Acad. Sci. Paris 2000, 328 (IIb), 411. (24) Qu, D.; Rame´, E.; Garoff, S. Phys. Fluids 2002, 14, 1154.

with stick-slip motion of the receding contact line, were reported by Shanahan,6,26 Birdi,3 and Hu.17 The complexity of the problem is noticeably due to the interplay between evaporation, hydrodynamic flow, contact line pinning, disjoining pressure for thin films, gravity for thick ones, and Marangoni flows if temperature or concentration gradients develop during the process.2,25 Therefore, it is interesting to define simple situations, where only a limited number of these effects come into play. For example, a drop of pure, completely wetting liquid with no specific interactions with the substrate is expected to recede without pinning, and with a very small contact angle. No significant temperature gradient can develop in such thin, flat drops, and gravity does not play at late times, when the drops are small enough. Therefore, pinning, gravity, and Marangoni flow will play no role. That situation is met with drops of alkanes, spreading on oxidized silicon wafers and mica. As a detailed report of experimental setup and data is presented elsewhere,7 we shall give only a brief summary of the results. The aim of the present paper is to examine the consequences of the mathematical structure of the evaporation term, and identify which effects are accounted for by considering only that term. When this is done, it is easier to recognize the other contributions. II. Experiments and Basic Equations The drops of alkanes are put on flat, smooth, isothermal, horizontal substrates. The drops first spread, then stop for some time, while the contact angle decreases, then recede and finally disappear completely. Let t0 be the vanishing time. The radius R(t) of each drop is measured as a function of time t, and plotted vs the time (t0 - t) to vanishing. The data for different drops collapse on a single curve at late times, and a log-log plot is a straight line, the slope y of which is slightly less than 0.5; see Figure 1. The contact angle (cf. Table 1) can be measured if equal thickness interference fringes are visible at the edge of the drop.28 On silicon wafers, this is the case for octane and nonane over the whole range of times. For heptane, the fringes are visible only just before drop vanishing, for hexane, no fringe is seen at all. On mica, the optical contrast is less and the fringes are visible only in the last second of the drop’s life. (25) Rafaı¨, S.; Bonn, D. Preprint, in preparation, 2002. (26) Shanahan, M. E. R. Langmuir 1995, 11, 1041. (27) Wang, G. AIP Conference proceeding 197, Drops and bubbles; American Institute of Physics: Melville, NY,1988. (28) To avoid disturbances, the contact angle is usually measured with the same magnification as the drop itself. Therefore, only small angles ( 0


To complete the specifications imposed on the Φ function, we note that eq 16 implies that the evaporation rate of the whole drop may be written as

∫01uΦ(u) du

dV ) -2πJ0R(t) dt


which is thus proportional to R(t). For this evaporation rate to be finite, we have to specify that Φ is such that

∫01uΦ(u) du < ∞


From now on, we search asymptotic solutions in power law form

R(t) ∼ C(t0 - t)y, with y > 0


∂h (R(t), t) ∼ D(t0 - t)x ∂r


B. Asymptotic Behavior in the Limit t f t0. 1. Study of the Radius R(t). Taking the limit t f t0 in eq 20 leads to

h0(0) ) J0

dt′ ∫0t R(t′) 0


We note that this last equation implies that, when t f t0, R(t) goes to 0 slower than (t0 - t); that is, y < 1. Subtracting eqs 31 and 20 gives

The method proposed here for determining the exponent y consists of evaluating and then identifying the leading asymptotic contributions of each of the two members of eq 32.

Evaporating Droplets

Langmuir, Vol. 18, No. 21, 2002 8075

2. Study of the Slope (∂h/∂r)(R(t),t). We analyze eq 21 when t f t0. In this limit, we have

The lhs of eq 32 can be written as

h0(0) - h0(R(t)) ) J0 1 1 h (0) - h0(0) + h′′0(0)R2(t) + . . . J0 0 2


)] (33)


h′0(R(t)) ∼ C(t0 - t)yh′′0(0) and




1 2 C h′′0(0)(t0 - t)2y (a) ∼ 2J0



dt′ y 0 - t′)


( )




∂h (36)

The asymptotic analysis of term (c) is presented in Appendix A. It is shown that y > 1/3 and that

∞ -u(y-1) e (1 0

1 (c) ∼ (t0 - t)1-y C


- Φ(e

)) du (37)

Using eqs 34, 36, and 37 in eq 32, and knowing that 2y > 1 - y (because y > 1/3), we obtain the implicit integral equation for y

1 ) 1-y

∫0∞e-u(y-1)(Φ(e-uy) - 1) du


nγnRn(t)∫0 ∑ n)2




(t0 - t)1-2y




∑ nγn(n + 2)y - 1



Knowing that y > 1/3, we obtain at leading order from eq 21

that is

1 (t - t)1-y (b) ∼ C(1 - y) 0



The first term of the rhs of eq 32 behaves like

1 t dt′ ∼ ∫t ∫tt R(t′) C (t



(R(t), t) ∼ -





(t0 - t)1-2y ∑ 2 n)2(n + 2)y - 1



On the other hand ∞


) ∫0 e-u(2y-1)Φ′(e-uy) du ∑ n)2(n + 2)y - 1 ∞


∫01Φ′(v)v1-(1/y) dv

1 y


so that we finally get

J0 ∂h (R(t), t) ∼ - 2( ∂r yC

∫01Φ′(v)v1-(1/y) dv) (t0 - t)1-2y


which gives

y A(y) ≡ 1-y

∫0 (Φ(v) - 1)v 1


Finally, if the integral

dv ) 0


Since A(y) is a continuous increasing function of y with

lim A(y) ) -∞




lim A(y) ) +∞ yf1


we are thus sure that there exists a unique solution y of the implicit equation A(y) ) 0 and that

1 0, then y <

∫0tΦ′(v)v1-(1/y) dv


exists, we obtain the general relation

x ) 1 - 2y


C. Specific Choices of Evaporation Function Φ. Comparison with Numerical Results. We now apply these general results to specific choices of the evaporation function Φ, and we compare them with a numerical analysis of eq 16. We begin with a short presentation of the integration scheme used. 1. Numerical Integration. Since there are no hydrodynamical or capillary effects, the time evolution of a drop only due to the evaporation as given by eq 16 can be easily integrated numerically. The simplest scheme, with a uniform stationary grid, has been tested to work well enough to determine the exponent y of R vs (t0 - t). The number of nodes (N) of the grid has been chosen to be 105 after testing that for 105 e N e 5 × 105 the difference in the values of R was less than 0.01%. The time step was set to vary as dt λn, with λ e 1 and n the number of iterations. This diminution in the time step dt is useful to avoid loosing information at the latest stages of evolution, where an acceleration of evaporation occurs. Anyway, we have found no dependence on λ except at the very end of the evaporation. Now, we will apply the results found in the general case, to special choices of the evaporation function Φ. The


Langmuir, Vol. 18, No. 21, 2002

Cachile et al.

first case corresponds to the original function introduced in refs 9-11 by Deegan, eq 5, mentioned in the Introduction. 2. Evaporation Function Proposed by Deegan.9-11 It is very easy to check that, in this case, A(1/2) ) 0, where A(y) has been defined by eq 39. In this case, the radius behaves like R(t) ∼ C(t0 - t)1/2. As for the slope at the edge of the drop, eq 47 leads to

lim tft0

2J0 ∂h (R(t), t) ) - 2 ∂r C


∫01(1v- v2)3/2 dv ) -∞


which means that the contact angle tends to the limiting value (π/2). Numerical integrations of this function have been performed, using spherical caps as initial condition for the profile. R0 ≡ R(t ) 0) was kept the same in all the initial conditions. When an initial profile evolves, after some transient time, the profile is no longer a spherical cap but a semiellipsoid. This fact implies that the contact angle is = π/2. Once the angle attains π/2, it remains constant, and the profiles remain elliptic. For long timessusually when the angle has already attained π/2sthe radius vs (t0 - t) starts to follow a power law, with an exponent 0.5, and keeps the same behavior until the drop disappears. More precisely, the simulation is stopped when R e 10-4R0, but the analysis is usually performed up to R e 10-2R0. With the notation of eq 49, the numerical results can be written as follows: x ) 0 (the angle remains constant) and y ) 0.5 ( 0.001, in excellent agreement with the analytical predictions. From a physical point of view, this last result does not really make sense, since the expression of the evaporation term Φ1 is valid only for low contact angles (cf. eqs 3 and 4). We consider now a Φ function that enables to remove the divergence of the evaporation term and which leads to much more reasonable results for the asymptotical behavior of the contact angle. 3. Nondivergent Evaporation Function. We have chosen the function:

Φ2(x) )


1 - e-m 1-x -m x1 - x2 1 - e 1



where m is a positive parameter. Φ2(0) ) 1 and Φ2(1) ) m/(1 - e-m)). Physically relevant values for m must be large enough to play only close to the drop’s edge. This is related to the fact that J2(x, t) ) (J0/R(t)) Φ2(x) is not a self-consistent boundary condition for the electrostatic problem. We give here some values of the exponent y, obtained by solving the implicit integral equation eq 39, and compare them with numerical results. On one hand, the integral given by eq 48 is finite

∫0 Φ′2(v)v1-(1/y) dv < ∞ 1


which implies that the exponent x is given by

x ) 1 - 2y


On the other hand, as was done for the case of Φ1, numerical integrations for Φ2 were performed. Spherical

Figure 6. Profile evolution from the numerical integration with m ) 3. The initial profile is a spherical cap with contact angle θ0 ) π/10. In the inset a log-log plot of R vs (t0 - t) is shown. Table 4. Comparison of the Exponent y Obtained by Asymptotical Analysis (Equation 39) and Numerical Integration m y from eq 39 y from numerical study x from eqs 53 and 39 2 3 5 ∞

0.415 0.438 0.464 0.5

0.415 ( 0.005 0.435 ( 0.005 0.465 ( 0.005 0.500 ( 0.001

0.17 0.124 0.072 0

caps were used again as initial conditions, keeping for every integration the same value of R0. The numerical integrations were performed for several values of m, ranging from m ) 2 to m ) 5. In every case, at long times, the radius follows a power law of the form (t0 - t)y where the value of y depends on m. The comparison between the analytical and numerical results for the exponent y are shown in Table 4. The profiles in this case are neither elliptic nor spherical, and the contact angle tends to 0 when (t0 - t) f 0. In Figure 6, we show the evolution of a profile with m ) 3 and θ0 ) π/10. In the inset a log-log plot of R vs (t0 - t) is shown, where a power law dependence is evident. In this case, the value of y is 0.435 ( 0.005. D. Results Concerning the Second Evaporation Term. Let us assume that the evaporation term is of the form shown in eq 15. Following exactly the same method as exposed in section III.B and doing the same kind of hypothesis on the Φ function, we obtain, when t f t0, R(t) ∝ (t0 - t)y and (∂h/∂r) (R(t), t) ∝ (t0 - t)x, where y is implicitly the solution of the integral equation


φ(v) - 1

∫01 v1+(1/y)






This relation stands for the equivalent of eq 49 in the case where the evaporation rate is proportional to the square of the radius of the drop and not to the radius of the drop. If an explicit function of the form (51) is used, m has to be small.

Evaporating Droplets

For example:

m ) 0.5 leads to y ) 0.55 x ) 0.45 m ) 0.25 leads to y ) 0.53 x ) 0.47 Note that a direct calculation for J ) cst (m f 0), assuming a spherical cap shape, gives x ) y ) 0.5 IV. A Summary of Results The comparison between the experimental results reported in Table 1 and the theoretical ones reported in Table 4 shows that our model works extremely well when no film is present and when the contact angle is very small. As the mathematical analysis is strictly exact, this allows us to conclude about the consequences of the shape chosen for the evaporation term. In particular, the relation between the exponents for radius and angle, 2y + x ) 1, which results empirically from the proportionality between evaporation rate and radius for a spherical cap, is found to have a wider validity; i.e., it holds as soon as eq 14 is valid. The mathematical analysis shows in which cases y is larger or smaller than 0.5. In particular, if one just smoothes at the edge the function J(r,R) proposed by Deegan,9-11 y is always less than 0.5, whatever the exact shape of the smoothing factor. Significant changes in J(r,R) over the whole range r < R are needed to obtain y > 0.5, and there is no physics behind this point. There is certainly some physics in the relation between y and m in the case of the smoothed function used eq 51: the exponent y is larger when the drop’s edge is steeper. However, the measured values are too close and no trend emerges. This is not surprising, as the study was done on linear alkanes, precisely because they had no specific interactions with the substrate. Complementary experiments with more complex liquids are clearly needed at this point. But the most interesting result is the fact that the contact angles are predicted to decrease in a model which does not include hydrodynamical contributions and, in a way which is controlled by the structure of J(r,R) at the very edge of the drop. All these conclusions become questionable when the contact angle is larger. The assumption of an isothermal receding motion is weakened. Worse, even the sign of a thermal gradient is not obvious. If we assume that the bottom of the drop is in thermal equilibrium with the substrate and that the heat diffusion is very fast, a stationary field of temperatures is achieved, with the center of the drop colder than the edges, which are at the temperature of the substrate. Then the surface tension gradient is directed inward, and the drop is supposed to recede faster. However, as the evaporation rate is very large at the edges, they might well be colder, in which case the surface tension gradient would be directed outward. This is the case for the main part of the receding motion for hexane, because we clearly see the edge of the drop feeding the film. In fact, evaporation is so fast with hexane that the initial outwardly directed surface tension gradient probably never fades out. In the presence of a film, the only general conclusion is that the measured value of y can be larger than 0.5 even with a decreasing contact angle (x > 0). To say more, a careful study of both drop and film, with specific attention paid at the edge of the drop, is needed. HMDZ on silicon wafers is a simple case, where the contact angle is small, the drop’s edge is quite visible, and in the part of the receding motion where y is measured, the drop does not

Langmuir, Vol. 18, No. 21, 2002 8077

feed the film. Actually, the predicted relation between the exponents, y + x ) 1, is well obeyed. However, the general situation is more complex, and this is obvious with hexane. V. General Conclusion As a general conclusion, our simple analysis, where only evaporation is considered, accounts remarkably well for the behavior of evaporating drops of alkanes and HMDZ on oxidized silicon wafers. The discrepancies observed on the contact angle close to drop vanishing are well understood. The association of rigorous mathematical techniques with numerical procedures to treat exactly a model system, itself obtained from a real situation through a preliminary careful physical analysis, appears to be extremely powerful. Now we are in the situation to go further. First, if the power laws describing the dynamics are found, it is not the case for the prefactors. More precisely, the value of the contact angle when the drops start to recede has to be determined by a more complete analysis.2 The problem is certainly complex, because the substrate seems to play an important role.7 In the framework of our simple model, a better description of the evaporation rate at the drop’s edge is needed. This rate is proportional to Pv(equ), therefore can be calculated from the disjoining pressure Π(h), which contains all the specific liquidsubstrate interaction. Then, a convenient approximation of J(r,t) with the structure of eq 14 could be chosen. Afterward, the resulting values for the contact angle could be checked against experiment. Second, some data are not yet explained. For cyclohexane or water on mica for example, y is larger than 0.5 (y = 0.57 for cyclohexane and 0.68 for water11), although the drops recede with a relatively large contact angle and certainly without film. Thermal gradients could possibly come into play. However, one must be very careful with the interpretation of large y values: if some significant contact line pinning takes place at the microscopic scale, the receding motion slows down. A larger value of t0 results. As the parameter for the analysis is t0 - t, one sees immediately that the measured value of y is larger in that case. In other words, a large value of y is not necessarily due to a faster receding dynamics. A slowing down due to pinning would do exactly the same. A systematic study of the more complex cases is therefore needed. Acknowledgment. The study was initiated by D. Bonn, who brought the problem to our attention. We gratefully acknowledge helpful discussions with M. Brenner, H. Stone, and M. Betterton, who gave access to unpublished work. We thank the reviewers for judicious comments. Appendix A: Asymptotical Analysis of the Term (c) of Equation 32 In this appendix, we give details on the asymptotic evaluation of the term (c) of eq 32 and show that y > 1/3. All along this appendix we denote the lhs (rhs, respectively) the left-hand side (right-hand side, respectively) of eq 32. First, we start by using the series expansion of Φ in eq 22 to get ∞

γnRn(t)∫0 ∑ n)2

(c) ) -








Langmuir, Vol. 18, No. 21, 2002

Cachile et al.

We begin with assuming successively that y < 1/3 and y ) 1/3 and then we show that these conditions lead to inconsistent results. (a) If we assume that y < 1/3, then


(c) ∼ -γ2


dt′ C2 (t0 - t)2y R3(t′)



(b) If we assume that y ) 1/3, then



dt′ (t0 - t)2y 3 R (t′)



γ2 ln(t0 - t) (t0 - t)2/3 C


rhs ∼ -

γ2 ln(t0 - t) (t0 - t)2/3 C


h′′0(0) 2 C (t0 - t)2/3 2J0


(A2) Thus

As a consequence,

rhs ∼ -γ2C2

(c) ∼ -




lhs ∼ lhs ∼ -

1 2 C h′′0(0)(t0 - t)2y 2J0


which implies that


h′′0(0) ) 2J0γ2


dt′ R3(t′)



(c) ∼ -


γn(t0 - t)1-y

Cn)2(n + 1)y - 1

∞ 1 ) - (t0 - t)1-y γn C n)2

∞ ∑ ∫0 e-u((n+1)y-1) du

Knowing that

h′′0(0) < 0, γ2 > 0, J0 > 0,


Since γ2 * 0, these behaviors are not compatible, which implies that y * 1/3. (c) Last, we assume that y > 1/3. In this case

>0 ∫0t Rdt′ 3 (t′)

the condition (A5) is never satisfied. Thus, y g 1/3.


∫0∞e-u(y-1)(1 - φ(e-uy)) du

1 ) (t0 - t)1-y C


which is the equation (eq 37) given in the main text. LA0204646