3972
Langmuir 2005, 21, 3972-3980
Analysis of the Effects of Marangoni Stresses on the Microflow in an Evaporating Sessile Droplet Hua Hu* and Ronald G. Larson Department of Chemical Engineering, University of Michigan, Ann Arbor, Michigan 48109-2136 Received October 6, 2004. In Final Form: January 25, 2005 We study the effects of Marangoni stresses on the flow in an evaporating sessile droplet, by extending a lubrication analysis and a finite element solution of the flow field in a drying droplet, developed earlier.1 The temperature distribution within the droplet is obtained from a solution of Laplace’s equation, where quasi-steadiness and neglect of convection terms in the heat equation can be justified for small, slowly evaporating droplets. The evaporation flux and temperature profiles along the droplet surface are approximated by simple analytical forms and used as boundary conditions to obtain an axisymmetric analytical flow field from the lubrication theory for relatively flat droplets. A finite element algorithm is also developed to solve simultaneously the vapor concentration, and the thermal and flow fields in the droplet, which shows that the lubrication solution with the Marangoni stress is accurate for contact angles as high as 40°. From our analysis, we find that surfactant contamination, at a surface concentration as small as 300 molecules/µm2, can almost entirely suppress the Marangoni flow in the evaporating droplet.
1. Introduction Marangoni effects, manifested as “tears of wine”, were observed as early as the 1800s.2 In a wine glass, the evaporation of alcohol generates a surface tension gradient, which produces a traction on the wine surface causing the wine to climb up the side of glass where it forms a thin film. As the wine accumulates, a bulging rim of liquid forms along the top of the film, which eventually pinches into droplets which roll under their own weight, like tears, back into the wine. The Italian physicist Marangoni gave a detailed description of the movement of a liquid surface induced by a surface tension gradient, generated either by a composition or a temperature variation along the free surface. In 1901, Be´nard3 discovered the convection cells in a thin liquid film that were later named after him. Later, Block4 performed more careful experiments on a thin liquid film, and Pearson5 gave a detailed theoretical analysis of Be´nard’s observation, concluding that a surface tension gradient, i.e., the Marangoni stress, causes the convection patterns in a thin film. While much experimental and theoretical work has been performed to investigate various surface-tension-driven (i.e., Marangoni) flows in thin films or shallow pools,6-12 there are only a few papers reporting Marangoni flow in an evaporating sessile droplet. Zhang and Yang13 experimentally studied the natural convection in evaporating drops, where they observed a Marangoni flow. Davis and * Corresponding author. E-mail:
[email protected]. (1) Hu, H.; Larson, R. G. Langmuir 2005, 21, 3963. (2) Scriven, L. E.; Sternling, C. V. Nature 1960, 187, 187. (3) Be´nard, H. Ann. Chim. Phys., Ser. 7 1901, 23, 62. (4) Block, M. J. Nature 1956, 178, 650. (5) Pearson, J. R. A. J. Fluid Mech. 1958, 4, 489. (6) Scriven, L. E. Chem. Eng. Sci. 1960, 12, 98. (7) Nield, D. A. J. Fluid Mech. 1964 19, 341. (8) Fanton, X.; Cazabat, A. M. Langmuir 1998, 14, 2554-2561. (9) De Gennes, P. G. Eur. Phys. J. E 2001, 6, 421-424. (10) Stroock, A. D.; Ismagilov, R. F.; Stone, H. A.; Whitesides, G. M. Langmuir 2003, 19, 4358-4362. (11) Mancini, H.; Maza, D. Europhys. Lett. 2004, 66, 812-818. (12) Vogel, M. J.; Miraghaie, R.; Lopez, J. M.; Hirsa, A. H. Langmuir 2004, 20, 5651-5654. (13) Zhang, N. L.; Yang, W. J. Trans. ASME 1992, 104, 656-662.
co-workers14-16 studied both stationary and spreading droplets with consideration of evaporation, Marangoni stresses, and moving contact lines. They applied a lubrication analysis to obtain the evolution equation of the droplet free surface profile as a droplet spreads or resides on a heated substrate. The velocity fields in the droplet were thereafter obtained using the evolution equation for the droplet profile. Savino and co-workers17 numerically solved the axisymmetric steady-state Navier-Stokes equations taking into account the Marangoni stress at the liquid-air interface. They also performed experiments to map the velocity field in the droplets. Their theoretical results were not consistent with the experimental ones due to the errors in using the PIV technique to map the velocity field in a spherical cap droplet, which acts like a lens distorting the real flow field. Meanwhile, many researchers have taken advantage of Marangoni flow in drying droplets to affect the pattern of deposition from the droplet onto the underlying substrate.18-21 Wang and co-workers,18 for example, produced a patterned porous thin film on a substrate. Stebe and co-workers20,21 used a surfactant that during drying develops a concentration gradient along the droplet surface leading to a gradient in surface tension; by varying the initial surfactant concentration, they were able to vary the pattern of particle deposition on the substrate. Marangoni flow plays a key role in coating, thin film deposition, crystal growth, and production of photonic materials. Therefore, a thorough understanding of Marangoni flow in an evaporating droplet will be impor(14) Ehrhard P.; Davis, S. H. J. Fluid Mech. 1991, 229, 365. (15) Anderson D. M.; Davis, S. H. Phys. Fluids 1995, 7, 248. (16) Oron A.; Davis, S. H.; Bankoff, S. G. Rev. Mod. Phys. 1997, 69, 931. (17) Savino R.; Paterna, D.; Favaloro, N. J. Thermophys. Heat Transfer 2002, 16, 562-574. (18) Wang, H. T.; Wang, Zh. B.; Huang, L. M.; Mitra, A.; Yan Y. S. Langmuir 2001, 17, 2572-2574. (19) Maillard, M.; Motte, L.; Pileni, M. P. Adv. Mater. 2001, 13, 200204. (20) Nguyen, V. X.; Stebe, K. J. Phys. Rev. Lett. 2002, 22, 32823285. (21) Truskett, V.; Stebe, K. J. Langmuir 2003, 19, 8271-8279.
10.1021/la0475270 CCC: $30.25 © 2005 American Chemical Society Published on Web 03/26/2005
Microflow in an Evaporating Droplet
Langmuir, Vol. 21, No. 9, 2005 3973
tant both in fundamental research and in practical applications. Although Davis and co-workers and Savino and coworkers have developed theories for the Marangoni flow in an evaporating droplet, an analytical theory for the locally resolved axisymmetric flow in a slowly evaporating droplet with a pinned contact line is still lacking. Such a theory is needed, for example, to predict particle deposition from a drying droplet in the presence of the Marangoni flow when, as usual, the particles do not rapidly diffuse across the height of the droplet during flow and deposition. In what follows, in section 2, we develop a lubrication theory to describe the velocity field in an evaporating droplet in the presence of Marangoni stresses. To demonstrate the accuracy of the lubrication theory, we also develop a finite element method to solve simultaneously the thermal and flow fields in the evaporating droplet. In section 3, we present the results obtained from these methods and discuss the effect of surface-active contaminants on the Marangoni flow in the evaporating droplet. We finally summarize in section 4. 2. Theory 2.1. Expressions for the Velocity Field with a Thermal Marangoni Stress Boundary Condition. In our companion paper,1 we established the governing equations for an evaporating droplet without considering the thermal transfer due to latent heat of evaporation. Since in many experiments there is strong evidence of the thermal cooling affecting the flow pattern in the drying droplet, we here add the energy equation to the set of governing equations. Therefore, we have the Laplace equation for the vapor concentration distribution
∆c ) 0
(1)
The flow equations are
1 ∂(rur) ∂uz + )0 r ∂r ∂z µ
(( ) ) ( ( ) )
µ
(2)
∂2ur ∂P ∂ 1 ∂ (rur) + 2 ) ∂r r ∂r ∂r ∂z
(3)
∂2uz ∂P 1 ∂ ∂uz r + 2 ) r ∂r ∂r ∂z ∂z
(4)
where we have neglected inertial terms, since the Reynolds number Re ≡ Fu j rR/µ is small (0.003) for weak flow in the slowly evaporating droplet considered here. Here, we will also neglect the buoyancy-driven flow because of the small value of a dimensionless group B ≡ Fgh02C/7.1375β, which was introduced by Pearson5 to estimate the relative strength of the buoyancy-induced flow compared to that of Marangoni flow. We choose the water density F ) 1 g cm-3; water thermal expansion coefficient22 C ) 2.07 × 10-4 °C-1; the temperature coefficient of surface tension for water, β ) -0.1657 dyn cm-1 °C-1; g ) 980 cm s-2; and the droplet height h0 ) 0.04 cm. Thus, we obtain B ≈ 3 × 10-4, which shows that the buoyancy-induced flow is very weak compared with the Marangoni flow in the evaporating droplet. We also consider heat transfer in the droplet and the glass coverslip, on which the droplet rests. The energy equation is
(∂T∂t + u‚∇T) + k∆T ) 0
Fcp
(5)
where cp is the specific heat and k is the thermal conductivity. Here, we define an inverse Stanton number St-1 ) Fcpu j rR/k, which is a ratio of the convective to the thermal diffusive effects. In typical experiments with water droplets, such as those described in Hu and Larson,23 parameters in this group have the following approximate values: height-averaged radial velocity u j r ) 1 µm/s, contact line radius R ) 1 mm, k ) 1.4536 × 10-3 cal K-1 cm-1 s-1 (from Bird et al.24), Cp ) 1 cal g-1 K-1, and water density F ) 1 g cm-3, which gives St-1 about 0.02. This implies that the rate of the convective heat transfer is much smaller than the rate of the conductive heat transfer, and so we can neglect the convection term in the energy equation. In the finite element analysis presented shortly, we confirm that the heat convection term is negligible. We can also neglect the transient in the energy equation, which we show is valid in Appendix A by estimating the ratio of the relative rates of change of droplet height to that of temperature. Thus, the energy equation (5) simplifies to the Laplace equation
∆T ) 0
(6)
We solve eq 6 for a system containing both the droplet and a glass substrate. So the boundary conditions for eq 6 are as follows: 1. On the droplet surface S ) {h(r,t)|r e R}: Jh ) Hw(J‚n), where Hw is the latent heat of evaporation of water and Jh is the heat flux. 2. On the glass surface outside the droplet z ) 0: R < r < ∞: Jh ) 0. 3. At the lower boundary of the glass z ) -hg, 0 < r < ∞ and in the glass far away from the central axis -hg < z < 0, r f ∞: T ) constant where hg is the thickness of the coverslip. The constant-temperature boundary condition 3 on the underside of the glass coverslip is reasonable in experiments with a water-immersion objective (which we will describe in a future publication), since in that case the gap between the coverslip and the microscope objective is filled with water, which acts as a heat bath. Note that the boundary condition on the free surface of the droplet assumes no heat loss due to either conduction or natural convection in the air above the droplet. These thermal boundary conditions could readily be modified, if necessary to describe other experimental conditions. We now develop a semianalytical lubrication solution for the flow field produced by droplet evaporation under the additional approximation of a relatively flat droplet, h/R , 1. Our procedure is numerically to solve the vapor concentration and heat equations, eqs 1 and 6, to obtain the temperature profile along the liquid-air interface and to obtain the Marangoni stress, which becomes a boundary condition for the momentum equations (3) and (4), used in the lubrication solution as described below. We check the accuracy of this semianalytical solution using a finite element analysis to solve simultaneously the mass balance equation (1), the flow equations (2)-(4), and the heat equation (6). A Marangoni stress, which is a surface-tension gradient on a free liquid surface induced by a temperature or a surfactant-concentration gradient, must be balanced by (22) Perry’s Chemical Engineers’ Handbook, 7th ed.; McGraw-Hill: New York, 1997. (23) Hu, H.; Larson, R. G. J. Phys. Chem. B 2002, 106, 1334-1344. (24) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport phenomena; John Wiley & Sons: New York, 1960.
3974
Langmuir, Vol. 21, No. 9, 2005
Hu and Larson
a shear stress in the liquid, which drives a recirculating flow. From a force balance on the surface and assuming that the droplet is almost flat so that the tangential stress along the surface is approximately τrz|z)h, we obtain
τrz|z)h ) dσ/dr
(7)
where σ is the surface tension. Since the temperature on the surface of the evaporating droplet is nonuniformly depressed by evaporative cooling, the surface tension varies along the droplet free surface. In our finite element analysis (see below) we find that the temperature profile along the droplet surface is well fitted by
T/∆T0 ) ar˜ b + (1 - a)r˜ 2 + c
dσ/dT ) β
(9)
where β is a property of the liquid; for water, β ) -0.1657 dyn cm-1 °C-1. From the chain rule, the surface-tension gradient is then given by
dσ dT dσ ) R dr˜ dT R dr˜
(10)
Now, assuming a flat surface, we can write the boundary condition, eq 7, using the Newtonian expression for the shear stress of a viscous liquid, as
|
)z)h(r,t)
|
∂uz Ma (abr˜ b-1 + 2(1 - a)r˜ ) tf ∂r
z)h
(11)
where the thermal Marangoni number is defined as Ma ≡ -β∆T0tf/µR, which is a ratio of the Marangoni force to the viscous force, and µ is the solvent viscosity. We apply the new boundary condition (11) within the lubrication theory developed in our companion paper1 and obtain the flow equations with a Marangoni-stress boundary condition
u˜ r )
(
) ( )} ( )
3 1 1 z˜ z˜ 2 [(1 - r˜ 2) - (1 - r˜ 2)-λ(θ)] 2 - 2 + 8 1 - ˜t r˜ h ˜ h ˜
{
r˜ h02h ˜ 2
R
(J ˜ λ(θ)(1 - r˜ 2)-λ(θ)-1 + 1)
3 z˜ 2 z˜ 2h h ˜ ˜2
+
˜ Mah0h 3 z˜ 2 z˜ (abr˜ b-1 + 2(1 - a)r˜ ) (12) 2R 2h h ˜ ˜2
(
)
3 1 z˜ 2 z˜ 3 [1 + λ(θ)(1 - r˜ 2)-λ(θ)-1] + 4 1 - ˜t 3h ˜2 h ˜
(
)
z˜ 2 3 1 z˜ 3 [(1 - r˜ 2) - (1 - r˜ 2)-λ(θ)] - 3 h ˜ (0, ˜t) 2 2 1 - ˜t 2h ˜ 3h ˜ h02 z˜ 3 2 -λ(θ)-1 2 (J ˜ λ(θ)(1 r ˜ ) + 1) z ˜ + R2 h ˜ r˜ 2h02 z˜ 3 J ˜ λ(θ)(λ(θ) + 1)(1 - r˜ 2)-λ(θ)-2 z˜ 2 2 R h ˜ r˜ 2h02 z˜ 3 2 -λ(θ)-1 (J ˜ λ(θ)(1 r ˜ ) + 1) h ˜ (0, ˜t) R2 h ˜2 Mah0 Mah0 2 b-2 z˜ 3 (ab r˜ (abr˜ b + + 4(1 - a)) z˜ 2 + 4R 2R h ˜
{
(
)
( ) () }
( )
(8)
where a, b, and c are the fitting parameters representing the shape of the temperature profile, ∆T0 is the temperature difference between the edge and the top of the droplet, and r˜ ≡ r/R. For most liquids, the surface tension decreases with increasing temperature and does so linearly for small temperature variations, i.e.
∂ur ∂z
u˜ z )
()
2(1 - a)r˜ 2)
z˜3 h ˜ (0,t˜) (13) h ˜2
where we have defined the dimensionless variables: u˜ r ) ˜ ) h/h0. tf is urtf/R; u˜ z ) uztf/h0; ˜t ) t/tf; r˜ ) r/R; z˜ ) z/h0; h the drying time, R is the contact line radius, h0 is the ˜ ≡ -J(0,t)/Fh˙ (0,t) is the initial height of the droplet, and J dimensionless mass flux at the top of the droplet surface, where J(0,t) and h˙ (0,t) are given by eqs 11 and 23, respectively, in ref 1. When there is no Marangoni stress along the droplet surface, i.e., the Marangoni number Ma is zero, eqs 12 and 13 reduce to the solutions in our companion paper.1 When the velocity gradient ∂uz/∂r at the free boundary surface is neglected (by dropping the terms in braces, see the discussion in ref 1), eqs 12 and 13 reduce to the simplified “classical” expressions for the Marangoni boundary condition. Terms arising from ∂uz/∂r are normally dropped in a standard lubrication analysis but are included here because of the flow singularity at the droplet edge, as is discussed in our companion paper.1 2.2. General Expressions for the Velocity Field with Marangoni Stresses. In section 2.1, an approximate solution for the full velocity field with the Marangoni stress on the droplet free surface was derived, in which the Marangoni stress was generated by a temperature gradient. However in some situations, such as those examined experimentally by Nguyen and Stebe,20 the Marangoni stress may be generated by a surfactant concentration gradient or by both temperature and surfactant concentration gradients, and the resulting flow affects the solute deposition behavior. For an arbitrary Marangoni stress distribution, we define a function g(r,t), in the Marangoni boundary condition
∂ur ∂z
|
z)h(r,t)
)
|
dσ µ dr
z)h
|
∂uz ∂r
≡ g(r,t)
z)h
(14)
The function g(r,t) combines the surface tension gradient (produced by temperature, surfactant-concentration gradient, or both) and a boundary vertical velocity gradient (produced by evaporation) along the radial direction. Once these are determined, ∂ur/∂z|z)h(r,t) can be used as a boundary condition to solve the velocity field, giving
u˜ r )
( (
3 1 1 [(1 - r˜ 2) - (1 - r˜ 2)-λ(θ)] 8 1 - ˜t r˜ ˜ g(Rr˜ ,tf˜t)tfh0h 2R
) )
z˜ 2 z˜ -2 h ˜2 h ˜ 3 z˜ 2 z˜ (15) 2h h ˜ ˜2
Microflow in an Evaporating Droplet
Langmuir, Vol. 21, No. 9, 2005 3975
(
)
3 1 z˜ 2 z˜ 3 [1 + λ(θ)(1 - r˜ 2)-λ(θ)-1] + 2 4 1 - ˜t 3h ˜ h ˜ z˜ 2 3 1 z˜ 3 [(1 - r˜ 2) - (1 - r˜ 2)-λ(θ)] h ˜ (0, ˜t) + 2 1 - ˜t 2h ˜ 2 3h ˜3 g′(Rr˜ ,tf˜t)tfh0 2 z˜ 3 g(Rr˜ ,tf˜t)tfh0 2 z˜ 3 z˜ + z˜ 4r˜ R 4R h ˜ h ˜ r˜ g(Rr˜ ,tf˜t)tfh0 z˜ 3 h ˜ (0,t˜) (16) 2R h ˜2
u˜ z )
(
)
(
(
)
)
where g′(Rr˜ ,tf˜t) is the derivative of g(r,t) with respect to r˜ . For any specific case, if we know the surface tension distribution and boundary vertical velocity gradient along the free surface, then we can obtain g(r,t), which allows us to obtain u˜ r and u˜ z from the above equations. We must note, however, that while the temperature distribution along the droplet surface is dominated by heat diffusion and can therefore be obtained independently of the flow field, surfactant-concentration gradients are controlled by the flow. Hence the Marangoni stress boundary condition produced by surfactants must be obtained simultaneously with the velocity field. 2.3. Finite Element Model. A finite element model corresponding to the governing equations (1) to (6) is
Mc ) J
(17)
[C(u) + L]T ) F′
(18)
K(T)u ) F(T)
(19)
where the rate of mass transfer due to diffusion is represented by the matrix M, which is an assembly over ne all elements in a mesh, namely, M ) ∑e)1 me. The total rate of mass transfer at the boundary is represented by ne e j . The rate of heat transfer per unit J, where J ) ∑e)1 temperature due to convection is represented by the matrix C(u), which is an assembly over all elements in a mesh, ne namely, C(u) ) ∑e)1 c(u)e. The rate of heat transfer per unit temperature due to conduction is represented by L, ne e l . The total rate of heat transfer at the where L ) ∑e)1 ne f′e. The boundary is represented by F′, where F′ ) ∑e)1 temperature-dependent viscous diffusion term is reprene ke. The total sented by the matrix K(T), where K(T) ) ∑e)1 force acting on the boundary is represented by the term ne e f . The boundary conditions are F(T), where F(T) ) ∑e)1 described in section 2.1 and the finite element equations (18) to (20) are solved simultaneously. Details of our method of solving these equations can be found in ref 1. In general, if both temperature and surfactant concentration gradients are present along the droplet surface, eq 33 in ref 1 should be used to solve the finite element model (18) to (20). Thus, the Marangoni stress given by eq 33 in ref 1 contains both thermal and surfactantconcentration-gradient Marangoni stresses. The latter stresses can be obtained by incorporating a mass balance equation for the surfactant concentration in the droplet. However, this is beyond the scope of this work and will be considered in future work. 3. Results and Discussions 3.1. Temperature Field. To obtain the Marangoni stress arising from a surface-tension gradient due to evaporative cooling, we first obtain the evaporation flux profile on the droplet free surface from either the analytic formula of Hu and Larson23 or the finite element solution
Figure 1. Contour plots of the temperature fields (in °C) in the droplet and the glass coverslip obtained by a finite element analysis of the heat equation, eq 18, at two contact angles. The parameters used for solving the heat eq 18 are as follows: vapor latent heat, Hw ) 541 cal g-1; thermal conductivity of water, kw ) 1.4536 × 10-3 cal cm-1 s-1 K-1; and thermal conductivity of glass, kg ) 2.2976 × 10-3 cal cm-1 s-1 K-1. The dimension of the glass substrate is taken to be 1.3 mm in radius and 0.15 mm in thickness. The results are insensitive to increases in the radius of the substrate.
of eq 17. With this nonuniform evaporation-flux distribution, the latent heat flux can then be calculated and applied as a boundary condition to solve the finite element heat equation. In our finite element analysis, the parameters used for solving the heat eq 18 and the flow eq 19 are as follows: vapor latent heat of evaporation of water Hw ) 541 cal g-1 (from ref 25), thermal conductivity of water kw ) 1.4536 × 10-3 cal cm-1 s-1 K-1 (from Bird et al.24), thermal conductivity of glass kg ) 2.2976 × 10-3 cal cm-1 s-1 K-1 (from Bird et al.24), viscosity of water µ ) 0.01 P, and temperature coefficient of the surface tension of water β ) -0.1657 dyn cm-1 K-1 (β is obtained from ref 25). By applying the finite element analysis, the temperature distributions in a droplet with a contact-line radius of 1 mm and initial height 0.364 mm are obtained for 40° and 10° contact angles and plotted in Figure 1. These contour plots show that at the initial contact angle of 40° the temperature increases from the top to the bottom of the droplet, and from the center to the edge of the droplet, while at a contact angle of 10°, the temperature decreases from the center to the edge of the droplet. Thus, the thermal field in the droplet changes significantly with evaporation time, and its radial gradient even reverses direction. The reversal of temperature-gradient direction occurs because at early times, the longer conduction path from the bottom of the glass to the top of the droplet makes the temperature lower at the top of the droplet than elsewhere, while at long times the faster rate of evaporation at the droplet’s edge makes it cooler there. From the temperature fields, the droplet surface temperature profiles at different contact angles are extracted and plotted in Figure 2. When the contact angle decreases (25) CRC Handbook of Chemistry and Physics, 63rd ed.; CRC Press: Cleveland, OH, 1984.
3976
Langmuir, Vol. 21, No. 9, 2005
Hu and Larson
Figure 2. The temperature profiles along the droplet surface at contact angles of 40°, 35°, 30°, 25°, 20°, 15°, and 10°. The solid lines are the finite element results, and the dashed lines are the fitting results using eq 8.
below 14°, the radial temperature gradient along the surface exhibits a transition from positive to negative slope so that the coldest point along the droplet surface shifts from the center to the edge of the droplet. The surface temperature from the finite element analysis can be fitted accurately by eq 8; see Figure 2. The choice of the functional form given in eq 8 is based on two considerations. First, across the most of the droplet, the temperature distribution can be represented by a simple parabolic function. However, near the contact line the temperature varies rapidly and deviates from the parabolic curve. Therefore, we included a term of higher order in r in eq 8 to correct this deviation. The lack of a term that is first order in radius r is due to the symmetry condition that the first derivative of the temperature along the droplet surface must be zero at r ) 0. Parameters a and b in eq 8 represent the strength of the temperature gradient along the droplet surface. The constant c is just the temperature at top of the droplet, divided by ∆T0, the temperature difference between the edge and the top of the droplet. The effect of the temperature field and the resultant Marangoni stress on the flow will be seen in the following paragraphs. 3.2. Velocity Field. The time-dependent velocity fields coupled with the thermal fields obtained from the finite element analysis are plotted in Figure 3, in which at the initial contact angle of 40° the Marangoni number Ma ≡ -β∆T0tf/µR ) 841 is large and therefore produces a strong recirculation in the droplet. As the contact angle decreases, the temperature gradient on the droplet surface is attenuated so that the Marangoni number decreases. When the contact angle decreases below a critical value of 14°, the recirculation disappears and the velocity vectors all point toward the edge of the droplet, because at this contact angle, the shear stress at the surface changes sign from negative to positive. Thus, the sign of the Marangoni number Ma changes at a contact angle of 14°. For the Marangoni-stress boundary condition, we can also obtain the analytical results if we know the parameters a, b, and c in eqs 12 and 13. We therefore fit the finite element temperature distributions at different contact angles depicted by the solid lines in Figure 2 with
Figure 3. The time-dependent velocity fields calculated by the finite element method with the Marangoni-stress boundary condition. Plots a to d are for contact angles 40°, 30°, 20°, and 10°, respectively. Table 1. The Parameters a, b, c, and ∆T0 Obtained by Fitting the Computed Temperature Profile by the Phenomenological Expression, Equation 8
40° 35° 30° 25° 20° 15° 10°
a
b
c
∆T0
av rel error, %
max rel error, %
0.334064 0.286465 0.302891 0.448445 0.487827 0.807235 0.833545
8 10 12 16 20 26 8
1601.529 2087.617 2532.286 2618.861 4045.307 3681.104 -5342.43
0.015594 0.011965 0.009865 0.00954 0.006177 0.006788 -0.00468
0.017 0.009 0.015 0.021 0.012 0.013 0.019
0.11 0.04 0.08 0.25 0.43 0.22 0.44
the empirical expression eq 8, which gives the dashed lines, which lie almost on top of the solid lines. The fitting parameters a, b, and c in eq 8 are listed in Table 1 along with the average and largest relative errors between the fitted and the finite element results. The errors are tiny for all contact angles, confirming the accuracy of the fitting function (8). Hence, with this accurate analytical representation of the temperature profile along the free surface, we can compute the velocity field using the analytical solution equations (12) and (13). As we shall see, the analytical flow field obtained in this way is nearly same as the flow field from the finite element analysis, shown in Figure 3. The streamlines for a contact angle of 40° obtained from the finite element analysis are very similar to those from the analytical solution, as can be seen from the plots in Figure 4. Similar agreement is obtained at a contact angle of 10° (not shown). Furthermore, the ur(z) and uz(z) profiles computed by the analytical solution nearly overlap those computed by the finite element model for contact angles
Microflow in an Evaporating Droplet
Langmuir, Vol. 21, No. 9, 2005 3977
Figure 4. Streamline plots of the flow field from the finite element model for the Marangoni-stress boundary condition at contact angle of 40°, from (a) FEM solution and (b) analytical solution.
40° and 10°ssee Figures 5 and 6sand all contact angles between these two (not shown). The maximum average relative difference between the finite element and the analytical results for radii up to r˜ ) 0.9 is less than 3%. It is surprising that the finite element model and analytical solution agree with each other even better than in the absence of the Marangoni stress; see ref 1. This level of agreement is retained even when the classical lubrication analytic solution is used; i.e., when the terms in braces in eqs 12 and 13, due to the boundary term ∂uz/∂r|z)h, are dropped. Evidently, the Marangoni stress tends to swamp some of the errors in the classical lubrication approximation, leading to a more accurate solution than in the absence of Marangoni stresses. 3.3. Effects of Maragoni Stress on the Velocity Field. Comparing the velocity profiles of ur(z) and uz(z) in Figure 6 of this paper, with Marangoni stress, to that of Figure 6 in our companion paper1 which neglects the Marangoni force, we find that the differences in the radial velocity profiles ur(z) are particularly pronounced. When the Marangoni number is zero, the gradient ∂ur/∂z|z)h is close to zero at the free surface. When the Marangoni number is negative (as is the case for a contact angle of 10°), there is a positive velocity gradient ∂ur/∂z|z)h on the droplet free surface, and when the Marangoni number is positive (as is true for a contact angle of 40°), this gradient is negative. The axial gradient of the radial velocity profile ur(z) is dominated by the contribution of the Marangoni term to ∂ur/∂z|z)h, which is much higher than that of ∂uz/∂r|z)h in the boundary condition (11). The result in Figure 6a shows that a negative Marangoni number promotes motion of the fluid from the center toward the edge of the droplet surface, i.e., a clockwise (or reverse)
Figure 5. (a) Radial and (b) vertical velocities versus vertical position at different radial positions r ) 0.1, 0.2 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 0.9 mm, for a contact angle of 40° and a Marangoni-stress boundary condition. The solid lines are from the finite element method and the dashed lines from the analytic solution in eqs 12 and 13.
circulation. The transition from counterclockwise (positive Ma) to clockwise (negative Ma) recirculation, as the contact angle drops below 14°, can be seen in Figures 5 and 6. For large fixed r (near the droplet edge to the right of the stagnation point) the radial velocity changes from positive to negative in sign with increasing z for counterclockwise rotation (Figure 5a) while the opposite occurs for clockwise rotation (Figure 6a). In Figure 6a, the zone of negative ur(z) is tiny, indicating a very weak recirculation. 3.4. Surface-Active Contaminants. The theoretical results in Figure 3 predict a strong recirculation flow in the water droplet. However, in experiments with drying water droplets at most a weak recirculation flow is observed.26 Moreover, a negligible Marangoni flow in water droplets is consistent with the commonly observed “coffee ring” pattern produced by solute deposition from a drying droplet.26,27 However, our FEM and analytical flow fields predict a very strong thermally driven Marangoni flow in drying water droplets. It has been reported, however, that surface-active contaminants that collect on the free surface can almost entirely suppress Maranogni flow in an evaporating water droplet and that low concentrations of (26) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R.; Witten, T. A. Nature 1997, 389, 827-829. (27) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R. Phys. Rev. E 2000, 62, 756-765.
3978
Langmuir, Vol. 21, No. 9, 2005
Hu and Larson
and so resides only at the interface between water and air. Since we have a relatively flat droplet, we can obtain the advection-diffusion equation for the surfactant along the free surface as follows
∂ν 1 ∂ D ˜ ∂ ∂ν (u˜ r|z)hr˜ ν) ) r˜ r ˜ ∂r ˜ r ˜ ∂r˜ ∂r˜ ∂t˜
( )
(22)
where ν is the number of surfactant molecules per unit area; D ˜ is dimensionless diffusivity, which is D ˜ ) Dtf/R2, and u˜ r|z)h is the radial velocity along the droplet surface. If we assume that the transfer of surfactant along the droplet surface is a quasi-steady-state process, we obtain, from eq 22
∂ν u˜ r|z)hν ) ∂r˜ D ˜
(23)
For a dilute concentration of surface-active agent, the surface pressure induced by surfactants on the free surface is
σS ) -νkBT
(24)
Using the chain rule, we then have
dσS dν dσS )R dr˜ R dν dr˜
(25)
Substituting eqs 23 and 24 into eq 25, we obtain the surfactant contribution to the Marangoni stress Figure 6. The same as Figure 5 except for a contact angle of 10°.
these contaminants are almost unavoidable for water surfaces, which are easily contaminated.17,20 We will therefore use our theory to analyze how a small concentration of surface contaminant can affect the Marangoni flow in an evaporating droplet. We consider a simple case, in which a Marangoni stress is induced by a surface tension gradient generated by both temperature and surfactant concentration gradients on the free surface of a relatively flat water droplet. Thus, the total Marangoni stress on the free surface is dσ/R dr˜ , given by
dσT dσS dσ ) + R dr˜ R dr˜ R dr˜
(26)
Combining eqs 20, 21, and 26 with eq 14, we obtain the axial gradient of the radial velocity on the free surface, which is represented by the function g(r,t)
g(r,t) ) -
Ma (abr˜ b-1 + 2(1 - a)r˜ ) tf kBT u˜ r|z)hν ∂uz µR D ∂r ˜
|
z)h
(27)
(20)
where dσT/R dr˜ is the Marangoni stress induced by a temperature gradient and dσS/R dr˜ is the Marangoni stress induced by a surfactant concentration gradient. If we know the total Marangoni stress in eq 20, we can substitute it into general solutions (15) and (16) to obtain an analytical solution. From eqs 8-10, we can derive the thermal contribution to the Marangoni stress
dσT µMa )(abr˜ b-1 + 2(1 - a)r˜ ) R dr˜ tf
kBT u˜ r|z)hν dσS )R dr˜ R D ˜
(21)
We now derive the surfactant contribution to the Marangoni stress. First, we make the additional assumption that the surfactant is insoluble in the droplet fluid
Substituting eq 27 into the general eqs 15 and 16, we then derive an analytical solution for the flow field in the presence of both temperature and surfactant concentration gradients on the free surface for a flat droplet, namely
u˜ r )
(
)
z˜ 2 3 1 1 z˜ [(1 - r˜ 2) - (1 - r˜ 2)-λ(θ)] 2 - 2 + 8 1 - ˜t r˜ h ˜ h ˜ ˜ r˜ h02h z˜ 3 z˜ 2 (J ˜ λ(θ)(1 - r˜ 2)-λ(θ)-1 + 1) + 2 2h R h ˜ ˜2 Mah0h ˜ 3 z˜ 2 z˜ (abr˜ b-1 + 2(1 - a)r˜ ) + 2R 2h h ˜ ˜2 u˜ r|z)hν kBTtfh0h ˜ z˜ 3 z˜ 2 (28) 2 2h D ˜ 2µR h ˜ ˜2
{
(
(
(
)
)
)}
Microflow in an Evaporating Droplet
Langmuir, Vol. 21, No. 9, 2005 3979
(
) ) )
z˜ 3 3 1 z˜ 2 [1 + λ(θ)(1 - r˜ 2)-λ(θ)-1] + 2 4 1 - ˜t 3h ˜ h ˜ 3 1 z˜ 2 z˜ 3 [(1 - r˜ 2) - (1 - r˜ 2)-λ(θ)] - 3 h ˜ (0,t˜) 2 2 1 - ˜t 2h ˜ 3h ˜ h02 z˜ 3 2 -λ(θ)-1 2 (J ˜ λ(θ)(1 r ˜ ) + 1) z ˜ R2 h ˜ r˜ 2h02 z˜ 3 J ˜ λ(θ)(λ(θ) + 1)(1 - r˜ 2)-λ(θ)-2 z˜ 2 + 2 R h ˜ r˜ 2h02 z˜ 3 2 -λ(θ)-1 (J ˜ λ(θ)(1 r ˜ ) + 1) h ˜ (0,t˜) R2 h ˜2 Mah0 2 b-2 Mah0 z˜ 3 (ab r˜ (abr˜ b + + 4(1 - a)) z˜ 2 + 4R 2R h ˜ u˜ r|z)hν kBTtfh0 2 z˜ 3 ˜3 2 z 2(1 - a)r˜ ) 2 h ˜ (0,t˜) z˜ + h ˜ D ˜ 4r˜ µR2 h ˜ νkBTtfh0 2 z˜ 3 ∂ (u˜ r|z)h) z˜ ∂r˜ 4D ˜ µR2 h ˜ u˜ r|z)hν kBTr˜ tfh0 z˜ 3 h ˜ (0,t˜) (29) D ˜ 2µR2 h ˜2
u˜ z )
(
()
(
( ) () ( ) { ( ) ( ) () }
In eqs 28 and 29, the terms in braces are the contributions from the surfactants on the free surface. If we neglect the effect of surfactant contaminants, i.e., the terms in braces, then eqs 28 and 29 reduce to eqs 15 and 16. To obtain the term u˜ r|z)h in eqs 28 and 29, we set z ) h in eq 28, and solve for u˜ r|z)h to give the radial velocity on the free surface
[
u˜ r|z)h ) -
3 1 1 [(1 - r˜ 2) - (1 - r˜ 2)-λ(θ)] 8 1 - ˜t r˜
2 ˜ ˜ 1 Mah0h 1 r˜ h0 h (J ˜ λ(θ)(1 - r˜ 2)-λ(θ)-1 + 1) (abr˜ b-1 + 2 2 R 2 2R
]/(
2(1 - a)r˜ )
1+
)
˜ ν kBTtfh0h (30) 2 D ˜ 4µR
In above equation, we find that the numerator (in brackets) is the radial velocity on the free surface in the absence of the surfactant. The denominator (in parentheses) represents the suppressing effect of the surfactant contaminants on this radial free surface velocity. We can estimate from this denominator the concentration, ν, required to suppress the radial velocity by, say, a factor of 100 by setting
˜ ν kBTtfh0h ) 100 D ˜ 4µR2
(31)
Using typical values for the parameters for a water droplet (R ) 1 mm, tf ) 360 s, D ≈ 1 × 10-5 cm2 s-1, h0 ) 0.36 mm), we obtain
ν∼
100 × 4 × 0.01 × (1 × 10-5) (4 × 10-14)(0.036 × 108)
≈ 300 molecules/µm2
(about 5 × 10-22 mol/µm2). From this value, we can easily find that for a water droplet with a contact line radius of 1000 µm, the total amount of surfactant on the free surface needed to suppress radial flow along the surface by a factor of 100 is about 1 × 109 molecules. We noted earlier that
in the absence of either thermal or surfactant Marangoni stress the radial velocity on the free surface is nearly zero but becomes large when thermal Marangoni stresses are present. We now find, however, that a tiny concentration of contaminants can greatly suppress the radial surface velocity that is produced by the thermal Marangoni stress. Hence, the presence of a trace concentration of surface contaminant can largely cancel out the thermal Marangoni effect. In experiments with water droplets, it is very difficult to control the amount of surfactants below the low level that suppresses thermal Marangoni flow, even when the experiments are carried out in a clean room.28-31 Other volatile fluids, whose surfaces are not so easily contaminated, might be expected to show strong Marangoni flows due to evaporative cooling.17 Our theory, which includes Marangoni effects due to both temperature variations and insoluble surface active agents, could be used to predict not only the effects of unintentional surfactant contaminants, but also the effects due to deliberately introduced surfactants, such as those employed at high concentrations by Stebe and co-workers20 to control the deposition patterns of solutes from drying droplets. At high concentrations near surface saturation, the pressure-concentration isotherms depart greatly from the ideal form given by eq 24, and interesting deposition patterns can thereby be induced. In order apply our theory to analyze the particle deposition process in the presence of the surfactants in the droplet, an additional convectivediffusion equation should be introduced to calculate the surfactant concentration distribution. Once the surfactant concentration distribution is obtained along the droplet surface, the Marangoni force can be calculated and the velocity field is then obtained from the general equations (15) and (16). Since the convective flow and surfactant concentration distribution equations are in general coupled, an iterative scheme will be needed to find simultaneous solutions for both the flow and surfactant concentration fields. We expect that our lubrication theory will be of considerable help in explaining and controlling these deposition patterns. 4. Summary We have performed a thorough theoretical study of the effects of Marangoni stress on flow in an evaporating droplet. Marangoni stress due to thermal gradients along the free surface produced by latent heat of evaporation is introduced into the free surface boundary condition, allowing an analytical solution to be obtained for the flow field, using the lubrication approximation. The solution requires specification of the temperature gradient along the droplet surface, which can be obtained by a numerical solution of the thermal field. We developed a finite element model to solve the coupled thermal and velocity fields in the droplet, both to obtain the needed thermal field for the lubrication solution, and to confirm the validity and accuracy of the lubrication theory for the flow field. These solutions show that the heat of vaporization, the nonuniform path lengths for heat conduction, and the nonuniform evaporation rate lead to a nonuniform distribution of temperature along the air-liquid interface and hence a nonuniform surface tension, which drives a thermal Marangoni flow. The lubrication approximation (28) Ward, C. A.; Stanga, D. Phys. Rev. E 2001, 64, Art. No. 051509. (29) Barnes, G. T.; Hunter, D. S. J. Colloid Interface Sci. 1982, 88, 437-443. (30) Cammenga, H. K.; Schreiber, D.; Rudolph, B. E. J. Colloid Interface Sci. 1983, 92, 181-188. (31) Schreiber, D.; Cmmernga, H. K. Phys. Chem. Chem. Phys. 1981, 85, 909-914.
3980
Langmuir, Vol. 21, No. 9, 2005
Hu and Larson
provides an accurate analytical solution to the timedependent axisymmetric flow field, including the recirculation due to the Marangoni-stress boundary condition, when compared to a finite element solution. We predict that the nonuniformity in heat-conduction path lengths produces a positive Marangoni number at a large contact angle, leading to an inward radial flow along the droplet surface, while the nonuniformity in evaporation rate yields a negative Marangoni number at small contact angles, yielding an outward flow. The general structure of the lubrication solution admits inclusion not only of thermal Marangoni flows but also of flow induced by surfactant concentration gradients along the droplet surface. From our lubrication solution, we learn that a small concentration, as about 300 molecules/µm2, of insoluble surfactant on the free surface, significantly reduces the radial flow produced by the thermal Marangoni stresses. This explains why surfaces that are easily contaminated, such as water surfaces, typically do not show thermal Marangoni effects, while droplets with less easily contaminated surfaces do show pronounced thermal Marangoni flows. Appendix A Here we will show that the temperature field in the drying droplet is at quasi-steady state by estimating the ratio of the relative rates of change of droplet height to that of temperature
( ) 1 dh h0 dt 1 dT δT dt
(
)
Ah0FCp
(A2)
where A is the area of the substrate covered by the droplet and Hw is the heat of vaporization. The decrease in temperature continues until it is limited by conduction of heat from the substrate. Thus, the eventual fluid temperature drop δT can be estimated by a quasi-steadystate energy balance, yielding
dh kδT ∼ FHw h0 dt
(A3)
where the left side is the rate of which heat flows into the droplet from the substrate and the right side is the rate at which latent heat is lost by evaporation. Incorporating these two estimates in eqs A2 and A3 into eq A1 gives
1 dh CpFh0 dh/dt h0 dt ) 1 dT k δT dt
(A4)
Noting that dh/dt is approximately the vertical velocity, u j z, and that the ratio of the vertical to the radial velocity is of the order of the drop aspect ratio h0/R, we obtain
1 dh h0 dt h0 2 -1 = St 1 dT R δT dt
()
(A1)
where δT is the decrease in temperature of the liquid due to latent heat of vaporization resulting from evaporation. If the above ratio is small, then the droplet temperature equilibrates rapidly compared to the rate at which the droplet evaporates. The initial rate of temperature change can be estimated from a transient adiabatic energy balance, yielding
dT dh ∼ AFHw dt dt
(A5)
where St-1 ) FCpu˜ rR/k. Since h0/R < 1, and St-1 , 1, the droplet temperature field reaches steady state quickly compared to the rate of change in droplet height and the temperature field therefore remains at quasi-steady state. Acknowledgment. We thank NASA microgravity research division for supporting this study through Grant NAG3-2134 and NAG3-2708. LA0475270
