Langmuir 2005, 21, 3963-3971
3963
Analysis of the Microfluid Flow 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 The axisymmetric time-dependent flow field in an evaporating sessile droplet whose contact line is pinned is studied numerically and using an analytical lubrication theory with a zero-shear-stress boundary condition on the free surface of the droplet at low capillary and Reynolds numbers. A finite element algorithm is developed to solve simultaneously the vapor concentration and flow field in the droplet under conditions of slow evaporation. The finite element solution confirms the accuracy of the lubrication solution, especially when terms of higher order in the droplet flatness ratio (the ratio of droplet height to radius, h/R) are included in the lubrication theory to account more accurately for the singular flow near the contact line.
1. Introduction Upon drying, a droplet of liquid typically leaves a ring of solute on the substrate on which it rested. This phenomenon, which is caused by an evaporation-driven flow,1,2 has been used as the basis of a high-throughput DNA mapping technique.3 Thus, when a droplet of DNA solution is deposited onto a chemically treated substrate, the flow inside the droplet induced by evaporation carries DNA molecules toward the droplet’s edge. As the droplet dries out, DNA molecules, which are negatively charged, become bound to the substrate (if the substrate is positively charged) and are stretched out along the radial direction. The stretched DNA molecules can then be subjected to a restriction enzyme digestion, and the small gaps left by enzymatic digestion serve as an optical map of the restriction sites. Similarly, Blossey and Bosio4 reported that DNA microarrays, made by drying multiple droplets onto a glass surface, form ring patterns of DNA on the substrate, which negatively affects the accuracy of the subsequent fluorescence imaging of DNA/RNA hybridization. It was also recently reported5-8 that colloidal crystals can be made by deposition, layer by layer, of small particles initially suspended in a drying droplet. Modeling of these important deposition processes requires a solution (preferably analytical) for the velocity field in a drying droplet that is both accurate and simple enough to use in analysis of mass transport and deposition within the droplet. In general, droplet evaporation occurs in two stages, as shown in Figure 1. In the first stage, the contact line is pinned and the contact angle decreases during drying. When the contact angle reaches a critical value, which for water droplets on glass is about 2-4° as reported in our * Corresponding author. E-mail:
[email protected]. (1) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R.; Witten, T. A. Nature 1997, 389, 827-829. (2) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R. Phys. Rev. E 2000, 62, 756-765. (3) Jing, J. P.; Reed, J.; Huang, J.; Hu, X.; Clarke, V.; Edington, J.; Housman, D.; Anantharaman, T. S.; Huff, E. J.; Mishra, B.; Porter, B.; Shenkeer, A.; Wolfson, E.; Hiort, C.; Kantor, R.; Aston, C.; Schwartz, D. C. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 8046-8051. (4) Blossey, R.; Bosio, A. Langmuir 2002, 18, 2952-2954. (5) Velikov, K. P.; Christova, C. G.; Dullens, R. P. A.; van Blaaderen, A. V. Science 2002, 296, 106-109. (6) Cong, H. L.; Cao, W. X. Langmuir 2003, 19, 8177-8181. (7) Fan, F. Q.; Stebe, K. J. Langmuir 2004, 20, 3062-3067. (8) Cong, H. L.; Cao, W. X. Langmuir 2004, 20, 8049-8053.
Figure 1. A schematic of the droplet-evaporation process.
previous paper,9 the contact line starts to recede; this is the second stage. During the first stage, which typically occupies most of the drying time, the evaporation flux at the edge of the droplet is much higher than that at the center, leading to more solvent loss at the edge of the droplet than at the center. To keep the contact line pinned, the solvent must flow from the droplet center toward the edge to compensate for the solvent loss. Consequently, a flow is generated in the evaporating droplet whose characteristics depend on the rate of evaporation and the evaporation-flux distribution over the droplet surface. This is the mechanism by which nonuniform solute deposition occurs from an evaporating droplet, leading, for example, to the commonly observed “coffee ring” phenomenon first explained by Deegan and co-workers.1,2 Deegan and co-workers1,2 used a particle-tracking method to measure the height-averaged radial velocity in an evaporating droplet. They derived a height-averaged radial velocity and showed it to be largely consistent with their experimental results. They also measured and successfully predicted the time-dependent radial distribution of solute deposition, although their predictions can be questioned since they assumed rapid vertical diffusion of solute particles, which cannot be justified for the sizes of droplets and solute particles they used. Fischer10 applied a lubrication theory to investigate particle deposition patterns on the substrate under various evaporation conditions. Although he predicted the particle deposition patterns, the evaporation flux distribution on the droplet surface is not correctly described in his model, leading to (9) Hu, H.; Larson, R. G. J. Phys. Chem. B 2002, 106, 1334-1344. (10) Fischer, B. J. Langmuir 2002, 18, 60-67.
10.1021/la047528s CCC: $30.25 © 2005 American Chemical Society Published on Web 03/26/2005
3964
Langmuir, Vol. 21, No. 9, 2005
an inaccurate prediction of particle deposition from the drying droplet. Davis and co-workers in a series of papers11-13 studied both stationary and spreading droplets with consideration of evaporation and contact line speed. 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. Despite the successful theories of Deegan et al. and Davis et al., a theory for the complete locally resolved axisymmetric flow in a slowly evaporating droplet with a pinned contact line is still lacking. A full axisymmetric analytical solution for the flow field in the evaporating droplet is needed, for example, to predict DNA stretching and deposition in the high-throughput mapping application mentioned above, and is required for a proper understanding of particle deposition when, as usual, the particles do not rapidly diffuse across the height of the droplet during flow and deposition. Because particle or polymer deposition from a drying droplet can most readily be simulated by Brownian dynamics simulations, which require tiny time steps, it is important to develop accurate analytical expressions for the time-dependent velocity field. It would be cumbersome and computationally very expensive if one needed to obtain these velocities from a numerical solution to this time-dependent flow problem at each time step needed for the Brownian dynamics calculations. In what follows, in section 2, we develop a lubrication theory to describe the velocity field in the evaporating droplet neglecting the surface-tension gradient, i.e., the Marangoni stress, on the droplet free surface. To demonstrate the accuracy of the lubrication theory, we also develop a finite element method to solve simultaneously the vapor concentration and flow fields in the evaporating droplet. In section 3, we present and discuss the results obtained from these methods. In a later publication, we will describe the effects of Marangoni stress on this flow. 2. Theory 2.1. Governing Equations. In general, the shape of a sessile drop on a substrate is controlled by a normalstress balance on the free surface, which is influenced by gravity, fluid flow within the droplet, surface tension, and “vapor recoil”12 due to acceleration of liquid upon evaporation. When fluid flow is very slow, the capillary number Ca ) µu j r/σ is small and the surface tension dominants the normal stress balance, so that the shape is that of a static droplet. For small droplets, gravity is also negligible, and this shape is then simply a spherical cap, achieved at small values of Bond number Bo ) FgRh0/σ, which accounts for the balance of surface tension and gravitational force. Here F is the fluid density, g is the gravitational constant, R is the droplet radius, h0 is the initial height of the droplet, σ is the surface tension of the liquid, µ is the liquid viscosity, and u j r is the average radial velocity. For small water droplets with contact-line radii of 0.8-1.0 mm and heights of 0.3-0.4 mm and slow flows (around 1 µm/s), the Bond number is in the range 0.03-0.04 and the capillary number is around 10-8, so that the spherical-cap approximation is a good one. Consider then in Figure 2 a small droplet on a glass surface, whose shape is that of a spherical cap. The local (11) Ehrhard, P.; Davis, S. H. J. Fluid Mech. 1991, 229, 365-388. (12) Anderson, D. M. and Davis, S. H. Phys. Fluids 1995, 7, 248265. (13) Oron, A.; Davis, S. H.; Bankoff, SG. Rev. Mod. Phys. 1997, 69, 931-980.
Hu and Larson
Figure 2. A sessile spherical-cap droplet on the substrate in a cylindrical coordinate system.
radial coordinate and the local height are r and h(r,t), respectively. S ) {h(r,t) ) (R2/sin2 θ - r2)1/2 - R/tan θ|r e R} defines the surface of the droplet. At the interface between the liquid and the vapor, the vapor concentration c is assumed to be equal to the saturation value cv. Far above the droplet, the vapor concentration approaches an ambient value Hcv (where H is the humidity of the ambient air). For stagnant air (no air currents) the difference in water vapor concentration between the droplet surface and the far-field drives the evaporation of water into the air, according to the diffusion equation. In our previous work,9 we reported that because of the low evaporation rate, the dimensionless group R2/Dtf ) cv(1 - H)/F is small (about 10-5), where tf is the time for the droplet to evaporate and D is the vapor diffusivity. Hence, the vapor concentration field above the droplet can be considered to be at a quasi-steady state. We therefore can neglect the transient term in the diffusion equation and obtain a simple Laplace equation for the vapor concentration distribution
∆c ) 0
(1)
The boundary conditions for eq 1 have been described in our previous work.9 In the liquid, we need to consider mass and momentum balances. In a cylindrical coordinate system (Figure 2), these yield the continuity and Stokes equations
1 ∂(rur) ∂uz + )0 r ∂r ∂z
µ
(2)
(( ) ) ( ( ) )
µ
∂2ur ∂ 1 ∂ ∂P (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. The boundary conditions for eqs 2, 3, and 4 are as follows. The total shear stress at the free surface is given by
τrz|z)h(r,t) ) µ
( | ∂ur ∂z
+
z)h(r,t)
| )
∂uz ∂r
z)h(r,t)
(5)
For the zero-shear-stress boundary condition τrz|z)h(r,t) ) 0, we obtain
∂ur ∂z
|
z)h(r,t)
)-
|
∂uz ∂r
z)h(r,t)
(6)
Microfluid Flow
Langmuir, Vol. 21, No. 9, 2005 3965
Here we include the term ∂uz/∂r|z)h(r,t) that is dropped in a standard lubrication analysis. Although it is of higher order in the ratio of the droplet height to radius, this term is relatively large, even for rather flat droplets, evidently because of the singularity in evaporation flux J(r,t). This flux influences uz through the kinematic condition
uz|z)h(r,t) )
∂h J(r,t) + ∂t F
(7)
where ∂h/∂t is the downward movement of the free surface due to evaporation and J(r,t)/F is the solvent velocity induced by evaporation flux. From the relationship between the height and the volume of the droplet, using the parabolic approximation to the droplet free surface profile, we obtain
2m ˘ (t) ∂h ) h˙ (0,t)(1 - r˜ 2) ) (1 - r˜ 2) 2 ∂t FπR
(8)
where h˙ (0,t) is the rate of motion of the surface at the top of the droplet, m ˘ (t) is the rate of mass loss from the droplet due to evaporation, and r˜ ≡ r/R. Since the contact line is pinned during droplet evaporation, the nonslip boundary condition is applied on the substrate. 2.2. Analytical Solutions for Full Velocity Field. The rate of droplet flattening and the velocity field within the droplet depend on the evaporation rate. The evaporation flux at the surface is J(r,t) ) D(n‚∇c)|z)h(r), where n is a unit vector normal to the droplet free surface. In our previous study,9 we showed that for contact angles θ between 0° and 90° the evaporation flux along the droplet surface can be well approximated by the simple form given by Deegan and co-workers1
J(r,t) ) J0(θ)(1 - r˜ 2)-λ(θ)
(10)
where θ has units of radians. To a good approximation, J0(θ) is expressed by9
J0(θ) )
Dcv(1 - H) (0.27θ2 + 1.30)(0.6381 R 0.2239(θ - π/4)2) (11)
From eq 9 and a mass balance equation, Deegan and co-workers1 derived the dimensionless height-averaged radial velocity u j˜ r
u j˜ r ≡
ur ) Az2 + Bz + C
(13)
where the constants A, B, C are determined from the boundary conditions. (2) The shear stress on the droplet surface is approximated by τrz|z)h(r,t). (3) We take
|
∂uz ∂r
z)h(r,t)
≈
∂ (u ∂r z
|
z)h(r,t))
All three of these approximations should improve as the droplet becomes flatter. The second and third approximations become exact when the free surface is horizontal so that the z axis is perpendicular to the free surface. Our theory has improved on that of Deegan et al. by considering a shear stress at the free surface, which will allow us to model the Marangoni flow that occurs when the temperature is not uniform along the droplet surface. Most importantly, we have derived the full axisymmetric velocity field in the evaporating droplet, whereas Deegan et al. only presented a height-averaged radial velocity, and no axial velocity field. The full flow field will be needed to predict particle or DNA deposition during droplet evaporation. Combining eqs 7, 8, and 9, we have
uz|z)h(r,t) ) h˙ (0,t)(1 - r˜ 2) +
J(0,t) (1 - r˜ 2)-λ(θ) F
(14)
Differentiating both sides of eq 14 with respect to r gives
(9)
where λ(θ) is a parameter reflecting the uniformity of evaporation. When the contact angle is 90°, we obtain λ ) 0, implying that the evaporation flux is uniform along the surface of the droplet, while for small contact angles, λ is close to 0.5, yielding a singularity in evaporation flux at the edge of the droplet. We also found J0(θ) and λ(θ) as functions of the contact angle θ. λ(θ) has a linear dependence on θ, namely
1 θ λ(θ) ) 2 π
(1) Applying the lubrication theory to the flow equations (3) and (4), we obtain that the radial velocity profile ur(z) at each value of r is parabolic in z
u j rtf 1 1 1 ) [(1 - r˜ 2)-λ(θ) - (1 - r˜ 2)] (12) R 4 1 - ˜t r˜
From eq 12, the full velocity field can be derived using the lubrication approximation, valid for relatively flat droplets. This analytical solution is derived using the following approximations:
|
∂uz ∂r
z)h(r,t)
≈
-
∂(uz|z)h) ) ∂r
1 h˙ (0,t)2r˜ (1 + J ˜ (0,t)λ(θ)(1 - r˜ 2)-λ(θ)-1) (15) R
where J ˜ (0,t) ≡ -J(0,t)/Fh˙ (0,t) and J(0,t) ) J0(θ). Hence, substituting eq 15 into eq 6, we have
∂ur ∂z
|
) z)h(r,t)
1 h˙ (0,t)2r˜ (1 + J ˜ (0,t)λ(θ)(1 - r˜ 2)-λ(θ)-1) R (16)
Using the no-slip condition on the substrate (where z ) 0), we easily obtain C ) 0, and then eq 13 becomes
ur ) Az2 + Bz
(17)
From the above equation, we can calculate the derivative of ur with respect to z at z ) h(r,t), which is
∂ur ∂z
|
z)h(r,t)
) 2Ah(r,t) + B
(18)
From eqs 16 and 18 we can solve for B, which is
B ) -2Ah(r,t) +
2h˙ (0,t)r˜ (1 + J ˜ (0,t)λ(θ)(1 - r˜ 2)-λ(θ)-1) R (19)
Substituting (19) into (17) and integrating (17) from 0 to
3966
Langmuir, Vol. 21, No. 9, 2005
Hu and Larson
h with respect to z, we can derive an expression for u jr
u jr )
(
2h˙ (0,t)r˜ 1 Ah2(r,t) + -2Ah(r,t) + (1 + 3 R J ˜ (0,t)λ(θ)(1 - r˜ 2)-λ(θ)-1)
)
h(r,t) (20) 2
Hence, solving for A gives
A)
(
h˙ (0,t)r˜ -3 (1 + J ˜ (0,t)λ(θ)(1 u jr 2 R 2h (r,t)
)
r˜ 2)-λ(θ)-1)h(r,t) (21) For simplicity, we use h to represent h(r,t) and J ˜ to represent J ˜ (0,t˜) in the following. Inserting both A and B into eq 17, gives
(
)
h˙ (0,t˜)r˜ tfho 3 z˜ 2 z˜ j˜ r 3 (1 + J ˜ λ(θ)(1 u˜ r ) u 2 2 h ˜ h ˜ R2 3 z˜ 2 r˜ 2)-λ(θ)-1) z˜ (22) 2h ˜
(
)
where we have defined the dimensionless variables: u˜ r ) urtf/R; u˜ z ) uztf/h0; ˜t ) t/tf; r˜ ) r/R; z˜ ) z/h0; h ˜ ) h/h0. tf is the drying time, R is the contact line radius, and h0 is the initial height of the droplet and we have used h(r,t) ) h(0,t)(1 - r˜ 2). In our previous paper,9 we concluded that the evaporation rate and the rate of decrease of the droplet height remain nearly constant during evaporation. Thus, we can approximate the rate of the decrease of the droplet height by a simple form
h˙ (0,t˜) ≈ -
ho tf
(23)
Substituting eqs 12 and 23 into eq 22, the radial velocity in the evaporating droplet is derived
u˜ r )
(
) )}
z˜ 2 3 1 1 z˜ [(1 - r˜ 2) - (1 - r˜ 2)-λ(θ)] 2 - 2 + 8 1 - ˜t r˜ h ˜ h ˜ r˜ h02h ˜ 3 z˜ 2 z˜ 2 -λ(θ)-1 (J ˜ λ(θ)(1 r ˜ ) + 1) (24) 2h R2 h ˜ ˜2
{
(
With the known radial velocity u˜ r shown in eq 24, we are now able to derive the vertical velocity u˜ z through the continuity equation, giving
u˜ z )
(
)
3 1 z˜ 2 z˜ 3 + [1 + λ(θ)(1 - r˜ 2)-λ(θ)-1] 2 4 1 - ˜t 3h ˜ h ˜
(
)
3 1 z˜ 2 z˜ 3 [(1 - r˜ 2) - (1 - r˜ 2)-λ(θ)] h ˜ (0,t˜) 2 1 - ˜t 2h ˜ 2 3h ˜3 h02 z˜ 3 (J ˜ λ(θ)(1 - r˜ 2)-λ(θ)-1 + 1) z˜ 2 + 2 R h ˜ r˜ h02 z˜ 3 2 -λ(θ)-2 2 J ˜ λ(θ)(λ(θ) + 1)(1 r ˜ ) z ˜ R2 h ˜ r˜ h02 z˜ 3 (J ˜ λ(θ)(1 - r˜ 2)-λ(θ)-1 + 1) 2 h ˜ (0,t˜) (25) 2 R h ˜
{
(
)
( ) () }
The equations derived for the ordinary lubrication analysis (without the second boundary term ∂uz/∂r|z)h(r,t)) can be obtained from eqs 24 and 25 by dropping the terms in braces. Note that these terms in braces are higher order in droplet aspect ratio (by a factor of h02/R2) than the terms not in braces and so would be dropped in a standard lubrication analysis. However, the singularity at the droplet edge, where r˜ approaches unity, is higher order in the terms in braces, and so these terms must eventually become large close to the contact line, no matter how flat the droplet is, and so we obtain a more accurate solution when we retain the terms in braces. 2.3. Finite Element Model. 2.3.1. System of Linear Equations. In the last section, we developed a lubrication theory to obtain approximate analytical solutions for the radial and the vertical velocities in the evaporating droplet. To check the validity and accuracy of the approximate analytical solutions, we here also describe a finite element analysis to solve the flow field in the evaporating droplet. However, we are dealing with a moving-boundary, freesurface problem, which cannot be solved by commercial finite element software. Therefore, to obtain a numerical solution for eqs 1, 2, 3, and 4, we developed our own finite element code. By using the conventional variational method and the Galerkin approximation, the finite element model for the thermal and flow fields is then expressed as
Mc ) J
(26)
Ku ) F
(27)
where the rate of mass transfer due to diffusion is represented by the matrix M, which is an assembly over ne me. The total all elements in a mesh, namely, M ) ∑e)1 rate of mass transfer at the boundary is represented by ne e j . The viscous diffusion term is J, where J ) ∑e)1 ne ke. The total represented by the matrix K, where K ) ∑e)1 force acting on the boundary is represented by the term ne e f . (In the absence of Marangoni F, where F ) ∑e)1 stresses, this force is set to zero.) Here ne is the total number of elements in the finite element mesh. The superscript e means that these terms are calculated at the element level. The detailed descriptions and formulas are given in Appendix A. The derivation of eqs 26 and 27 can be found in Hughes.14 Equations 26 and 27 are, respectively, the finite element equations for the vapor concentration field given by eq 1 and for the flow field given by eqs 2-4. 2.3.2. Boundary Conditions. For the vapor concentration eq 26, the boundary conditions are described in section 2.1, and in ref 9. For the flow eq 27, the no-slip boundary condition holds along the bottom of the droplet at the liquid-glass interface. On the upper free surface of the droplet, however, the boundary conditions are complicated and involve a moving boundary. For an accurate computation of the flow field, in our finite element model we do not assume that the droplet is flat even though near the end of the drying process the contact angle is very small. Nevertheless, the shape of the droplet can still be regarded as a spherical cap during evaporation because the capillary number is miniscule (Ca ∼ 10-8). Thus, the boundary conditions are imposed along both tangential and normal directions on the droplet’s spherical free surface. In other words, along the tangential direction (14) Hughes, T. J. R. The finite element method: linear static and dynamic finite element analysis; Prentice Hall: Englewood Cliffs, NJ, 1987.
Microfluid Flow
Langmuir, Vol. 21, No. 9, 2005 3967
on the free surface a zero shear stress is imposed, and along the normal direction a kinematic boundary condition is used. First, let us consider the normal direction. The kinematic boundary condition with a phase change is
[ ∂h∂t ] + (n‚J) F
(n‚u) ) [nr, nz]‚ 0,
(28)
r sin θ r sin θ , 1R R
[
2 1/2
)) ]
( (
(29)
The first term [nr, nz]‚[0, ∂h/∂t] on the right-hand side of eq 28 is the projection of the droplet-surface-boundary motion onto the normal direction. For a spherical cap, the surface motion ∂h/∂t can be exactly calculated by the following equation
∂h ) ∂t
[
((
h4(0,t) - R4
) )
h2(0,t) + R2 4h (0,t) 2h(0,t) 3
2
1/2
-r
+
2
]
R2 + h2(0,t) 2h2(0,t)
h˙ (0,t) (30)
where h˙ (0,t) is the rate of motion of the surface at the top of the droplet, which can be calculated by the finite element model, eq 26. Alternatively, from eq 8, h˙ (0,t) can also be approximated by
h˙ (0,t) )
2m ˘ (t) FπR2
(31)
Second, let us consider the tangential direction. There are two kinds of shear-stress boundary conditions, namely, a stress-free boundary condition and a Marangoni-stress boundary condition induced by the nonuniform distribution of temperature along the droplet surface. For the shear-stress-free boundary condition, the following equation must be satisfied on the boundary
(t‚τ) ) 0
(32)
where t is a unit vector along the tangential direction given by t ) [-nz nr]. For the Marangoni stress boundary condition, a shear stress in the liquid must be generated to balance the surface tension gradient on the droplet surface. Thus, we have
(t‚τ) )
QT‚K‚Q‚u ) QT‚F
(34)
where the rotation matrix QT is
where (n‚J) is the evaporation flux along the direction normal to the free surface and is given by eq 12, where n is a unit vector along the normal direction, and nr and nz are the r and z components of the unit normal vector n, respectively. It is not difficult to calculate n, which is
[nr, nz] )
Therefore, in the FEM analysis, for the nodes on the droplet surface, a coordinate transformation of the flow equation (27) must be performed so that both tangential and normal boundary conditions in eqs 28, 32, and 33 can be precisely applied
dσ dr
(33)
In the above equation (33), the surface tension gradient dσ/dr along the droplet surface is generated by a variation of either temperature or composition. The results for the Marangoni-stress boundary condition are discussed in our companion paper in this issue.15 Because the droplet surface is a spherical cap, the surface is not parallel to either of the global coordinates. (15) Hu, H.; Larson, R. G. Langmuir 2005, 21, 3972.
[
·
· · ‚‚‚ Q) ‚‚‚ ‚‚‚
‚‚‚ nzi -nri
‚‚‚ nri nzi
‚‚‚
‚‚‚
‚‚‚ ‚‚‚ ‚‚‚ · · ·
]
(35)
for ith node on the droplet surface, where nri and nz1 are the r and z components of the normal unit vector at the droplet free surface for node i. 3. Results and Discussion The FEM model equations (26) and (27) and the analytical solutions (24) and (25) are used to compute the velocity field in a water droplet, whose contact line radius and initial height are 1.0 and 0.364 mm, respectively. The initial contact angle for this droplet is 40°, and we assume that it takes 360 s to dry out. (These conditions were chosen to match typical experiments described in Hu and Larson.9) To compute the time-dependent velocity fields, the finite element meshes are generated for different contact angles θ of 40°, 30°, 20°, and 10°, which approximately correspond, respectively, to reduced times of ˜t ≡ t/tf ) 0, 0.25, 0.5, and 0.75 through a linear relationship between θ and ˜t. Although we can also use eqs 7 and 8 to calculate the relationship between contact angle and evaporation time, as has been discussed in our previous work,9 the relationship between θ and ˜t is accurately approximated by a linear equation. The finite element mesh is used to provide points at which the analytical flow field is evaluated. With a zero-shear-stress boundary condition on the free surface of the droplet, the velocity-vector field as a function of time or contact angle is computed by the finite element model and plotted for different contact angles in Figure 3. Near the end of drying, the radial velocity is much higher than the vertical velocity almost everywhere in the droplet, while these two components are more nearly equal in magnitude at the beginning of drying, which is expected since the ratio of ur/uz scales with droplet aspect ratio R/h. At different contact angles or different times, we observe similar flow patterns in that near the droplet’s center the vertical velocities are higher than the radial velocities while around the edge of the droplet the radial velocities dominate. The streamlines obtained by the finite element method at the two typical contact angles of 40° and 10° are plotted in Figure 4; they have a roughly hyperbolic shape. The analytical solutions for the velocities expressed by eqs 24 and 25 are qualitatively identical to the finite element results plotted in Figures 3 and 4. A quantitative comparison is given below. We plot the radial and vertical velocities ur and uz against the vertical position z at different radial positions r in Figure 5a-b for a contact angle of 40° and in Figure 6 for a contact angle of 10°. These figures show that the approximate analytical results from eqs 24 and 25 are reasonably consistent with the finite element results, especially for the smaller contact angles. The differences in the radial and vertical velocities between two methods
3968
Langmuir, Vol. 21, No. 9, 2005
Figure 3. (a-d) The time-dependent velocity fields calculated by the finite element method with the shear-stress-free boundary condition. Plots a to d are for contact angles 40°, 30°, 20°, and 10°, respectively. The droplet has a contact line radius of 1.0 mm, initial height of 0.364 mm, and drying time of 360 s. Unless indicated otherwise, the sizes and the drying times of the droplet for each of the following figures are the same.
increase with the radial position. The functions of ur(z) and uz(z) obtained by the approximate analytical solution are almost identical to those obtained by the finite element model near the center of the droplet, where we find that the average relative differences in the radial velocities ur between these methods are about 0.1% and 0.05% for contact angles of 40° and 10°, respectively, while the average relative differences in vertical velocities uz are 0.5% for both contact angles. However, the curves of ur(z) and uz(z) obtained from the approximate analytical solution deviate more from those obtained by the finite element model near the edge of the droplet, where we find that at r ) 0.9 mm the average relative differences in the radial velocities ur are about 8% and 10% for contact angles 40° and 10°, respectively, while the average relative differences in vertical velocity uz are about 2% and 4%, respectively. These discrepancies in the radial and vertical velocities near the edge of the droplet arise presumably in part because the droplet surface is more tilted near the edge than at the center of the droplet, thus lessening the accuracy of the lubrication theory. Nevertheless, it is remarkable how accurate the lubrication solution is for a contact angle as large as 40°. In fact, near the edge of the droplet, the lubrication solution is actually more accurate for a 40° contact angle than for a 10° angle, probably because the singularity in the flow at the edge of the droplet is stronger for a flatter droplet, and this degrades both the lubrication and numerical solutions, so that it is not clear whether the lubrication or the numerical solution is more accurate at the droplet edge.
Hu and Larson
Figure 4. Streamline plots of the flow field from the finite element model at a contact angle of (a) 40° and (b) 10°.
If we apply a conventional lubrication solution to the drying droplet, in which the terms in braces in eqs 24 and 25 are dropped, we obtain the results that are compared to the finite element results in Figures 7 and 8 for contact angles of 40° and 10°. These two figures show that while the standard lubrication solutions are similar to the finite element results near the center of the droplet, as one moves away from center, the difference between the two solutions increases. While the analytical solution containing the extra free-surface terms (the terms in braces) in eqs 24 and 25 predicts a nonzero gradient in ur(z) at the free surface, the standard (or “classical”) lubrication solution in Figures 7 and 8 predicts a zero gradient. The validity of the classical lubrication theory thus requires not only a flat surface but also a uniform motion of that surface so that the term ∂uz/∂r in the shear-stress boundary condition (6) is zero. However, from the expression for uz in eq 7, we can see that both terms ∂h/∂t and J(r,t)/F on the right-hand side are functions of the radial position r. When r is small, both ∂h/∂t and J(r,t)/F are nearly uniform so that uz is almost uniform, and the gradient ∂uz/∂r is then nearly zero near the center of the droplet for any contact angle. In addition, for any contact angle the surface of the droplet is almost horizontal near the center. Since near the center of the droplet both requirements (i.e., small ∂uz/∂r and a flat surface) for validity of the lubrication theory are satisfied, the results obtained by both sets of analytical solutions are almost the same as the results obtained by finite element analysis (see Figures 5-8). However, as the radial position increases, from eq 7, we can see that (1) a large gradient ∂uz/∂r develops along the droplet surface and (2) the surface tips away from the horizontal. Thus, for large contact angles, both requirements for the classical lubrication theory fail near the
Microfluid Flow
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°. The solid lines are from the finite element method and the dashed lines from the analytic solution (27) and (28).
edge of the droplet, and then the predicted results from the two sets of analytical solutions are not equal. Even for a small contact angle (see Figures 6a and 8a), for which the droplet is very flat, there remains a difference between the two sets of the analytical solutions near the droplet edge, because although the surface is now nearly horizontal everywhere, the gradient ∂uz/∂r remains large and in fact is singular at the edge of the droplet. By including the higher order ∂uz/∂r term in the boundary condition, Figures 5 and 6 show that we obtain surprisingly accurate results even for droplets that are not flat. A simplified version of analytical solutions (24) and (25) has been successfully used to simulate the DNA deposition and stretching process in the evaporating droplet 16 by incorporating the flow equations directly into Brownian dynamics simulations of the DNA stretching. This has led to good agreement of the predicted stretch of DNA molecules deposited onto a surface during droplet drying with the corresponding experimental measurements. From the flow field solution in the absence of the Marangoni stress, we can now understand the mechanisms of the DNA stretching, coffee ring deposition, and other related phenomena in the evaporating droplet. (16) Chopra, M.; Li, L.; Hu, H.; Burns, M.; Larson, R. G. J. Rheol. 2003, 47, 1111-1132.
Langmuir, Vol. 21, No. 9, 2005 3969
Figure 6. The same as Figure 5 except for a contact angle of 10°.
However, due to thermal cooling or solute composition variation along the droplet free surface, a surface tension gradient, i.e., the Marangoni stress, is expected to be generated along the droplet free surface. This Marangoni stress causes additional convection currents in the droplet and can affect the particle deposition patterns.17 In a companion paper,15 we will account for the effects of Marangoni stress on the flow in an evaporating droplet. 4. Summary A lubrication analytical solution was obtained for the complete axisymmetric, time-dependent, flow field in the absence of the Marangoni stress. We also developed a finite element model for the velocity field in the droplet, which confirms the validity and accuracy of the lubrication theory under conditions of small capillary number where the interface remains that of a spherical cap. Because of the presence of a singularity at the edge of the droplet, the inclusion of terms of higher order in the droplet aspect ratio improves the accuracy of the lubrication solution. Using the flow field derived here, Brownian dynamics simulations of solute deposition or DNA stretching and deposition can be carried out conveniently. Appendix A. The Finite Element Method In the finite element model, eqs 26 and 27, the matrixes M, J, K, and F are summations over all the elements in (17) Nguyen, V. X.; Stebe, K. J. Phys. Rev. Lett. 2002, 88, 164501.
3970
Langmuir, Vol. 21, No. 9, 2005
Hu and Larson
Figure 7. (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°. The solid lines are from the finite element model while the dashed lines are from the “classical” analytical solution (27) and (28), obtained when the terms in braces are dropped.
a mesh. Their expressions at the element level are described in the following. In eq 26, the mass diffusion is described by
me ) [mab] ) [
∫ΩDNa,iNb,ir dΩ]
(A1)
e
where m is the rate of mass-diffusion matrix at the element level and D is the diffusivity of vapor. a and b are the indices of the element in the matrix, a, b ) 1, 2, 3, 4 for a quadralateral element. The subscript i ) 1,2 represents r and z, respectively, and the Einstein convention is used to sum over repeated subscripts. Ω is the integration domain in a element. The boundary conditions on the mass transfer equation are given by
je ) [ja] ) [
∫Γ
h
nijiNar dΓh]
(A2)
where je is the total rate of mass transfer over the boundary Γh at the element level, and ji is the mass flux along the ith direction. Na is a bilinear node function for a quadralateral element. The detailed expression for Na can be found in Hughes.14 In eq 27, the viscous diffusion term is described by
ke ) [kpq] ) [
∫Ω BaT‚(D + D′)‚Bbr dΩ]
(A3)
Figure 8. The same as Figure 7 except at a contact angle of 10°.
where p and q are indices of the components in the viscousdiffusion matrix ke at the element level, and p, q ) nsd(a - 1) + i, nsd ) 2, a ) 1, 2, 3, 4, and i ) 1, 2. BaT, D, and D′ are
BaT )
[
N Na,1 0 Na,1 a r Na,2 0 Na,2 0
[ ] [ ]
]
(A4)
2 0 D)µ 0 0
0 2 0 0
0 0 1 0
0 0 0 2
(A5)
1 1 D′ ) λ 0 1
1 1 0 1
0 0 0 0
1 1 0 1
(A6)
In eq A6, the value of the penalty coefficient λ must be large to satisfy closely the constraint of the continuity equation. A normal criterion for choosing the value of λ is that the ratio λ/µ ≈ 107-1013. Here, we choose λ/µ to be 1010. Equation A5 calculates the viscous behavior of the liquid, while eq A6 is the penalty term. The force at the free surface is described by
Microfluid Flow
fe ) [fp] ) [
Langmuir, Vol. 21, No. 9, 2005 3971
∫Γ
h
niσijNar dΓh]
(A7)
where fe is the total force over the boundary at the element level, and σij is the stress at the boundary of element. The kinematic boundary condition shown in eq 28 will be included in FEM eqs 26 and 27 when the elements are assembled together, see Hughes14 for details. Appendix B. Finite Element Mesh and Convergence From the analysis of the dimensionless groups in section 2.1, we know that mass and momentum transfer can be regarded as quasi-steady-state processes. Therefore, as long as we can solve the vapor concentration and velocity fields at any contact angle, the time evolution of the vapor concentration and the velocity can be determined through a simple linear relationship between contact angle and evaporation time, which has been reported in our companion paper.15 However, we still must refine the finite element mesh differently at each contact angle to obtain the converged results. The convergence of the vapor concentration field can be found in our previous paper.9 We performed a numerical test to solve the velocity field at a contact angle of 40° using successively refined finite element meshes with 2500, 4100, 8900, 13400, and 19000
elements. From this numerical test we find that when there are more than 10000 elements, the relative error defined in eqs B1 below is less than our convergence criterion of 0.01. For smaller contact angles, more elements, as many as 30000, are needed for accurate solutions. We therefore use finite element meshes that have between 15000 and 30000 elements, depending on the contact angle, and these meshes are generated by a commercial software package, Ansys 6.0 (Ansys, Inc). The relative errors in the velocity and temperature fields are defined as
euerr
)
|
n+1 n u j r|r)0.9 -u j r|r)0.9 n+1 u j r|r)0.9
|
< 0.01
(B1)
n where u j r|r)0.9 is the height-averaged radial velocity at r˜ ) 0.9 (or r ) 0.9 mm for a droplet with radius R ) 1 mm) for the nth finite element mesh refinement. Since the velocity field is axisymmetric, we can save some computational time by solving the problem over only half of the spherical-cap domain.
Acknowledgment. We thank NASA microgravity research division for supporting this study through Grant NAG3-2134 and NAG3-2134. LA047528S
