Fast Evaporation of Spreading Droplets of Colloidal Suspensions

ARTICLE pubs.acs.org/Langmuir

Fast Evaporation of Spreading Droplets of Colloidal Suspensions Kara L. Maki†,§ and Satish Kumar*,‡ †

Institute for Mathematics and Its Applications and ‡Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, Minnesota 55455, United States ABSTRACT: When a coffee droplet dries on a countertop, a dark ring of coffee solute is left behind, a phenomenon often referred to as the coffee-ring effect. A closely related yet lesswell-explored phenomenon is the formation of a layer of particles, or skin, at the surface of the droplet during drying. In this work, we explore the behavior of a mathematical model that can qualitatively describe both phenomena. We consider a thin axisymmetric droplet of a colloidal suspension on a horizontal substrate undergoing spreading and evaporation. In contrast to prior work, precursor films (rather than pinned contact lines) are present at the droplet edge, and evaporation is assumed to be limited by how quickly molecules can transfer out of the liquid phase (rather than by how quickly they can diffuse through the gas phase). The lubrication approximation is applied to simplify the mass and momentum conservation equations, and the colloidal particles are allowed to influence the droplet rheology through their effect on the viscosity. By describing the transport of the colloidal particles with the full convectiondiffusion equation, we are able to capture depthwise gradients in particle concentration and thus describe skin formation, a feature neglected in prior models of droplet evaporation. The highly coupled model equations are solved for a range of problem parameters using a finite-difference scheme based on a moving overset grid. The presence of evaporation and a large particle Peclet number leads to the accumulation of particles at the liquidair interface. Whereas capillarity creates a flow that drives particles to the droplet edge to produce a coffee ring, Marangoni flows can compete with this and promote skin formation. Increases in viscosity due to particle concentration slow down droplet dynamics and can lead to a reduction in the spreading rate.

’ INTRODUCTION When a coffee droplet dries on a countertop, a dark ring of coffee solute is left behind, a phenomenon often referred to as the coffee-ring effect. A closely related yet less-well-explored phenomenon is the formation of a layer of particles, or skin, at the surface of the droplet during drying. Figure 1 shows an image of skin formation from an experiment in which a droplet of 1-μmdiameter polystyrene spheres suspended in water evaporates on a microscope slide under ambient conditions. Similar experiments were performed by Deegan et al., who in some cases observed the growth of a cap of particles near the top of the droplet.1 This growth was found to happen at relatively short times and was ascribed to a flow that drives particles from the droplet edge to the top along the droplet surface. At later times, particles were found to move from the surface cap to the droplet edge, leading to a higher concentration of particles there. The coffee-ring effect has been studied extensively since the pioneering work of Deegan et al.2 because of its relevance in a number of materials- and biomedical-related applications.36 In most theoretical investigations of this phenomenon, it is assumed that the particle distribution is either uniform in the vertical direction (i.e., the particle distribution is independent of the vertical distance from the substrate)1,2,79 or that particle concentration gradients in the vertical direction are negligible (i.e., the particle distribution can be decomposed into a vertically averaged component and a small perturbation that includes the vertical dependence).10 However, describing skin formation requires accounting for vertical concentration gradients. The objective of the present work is to explore the behavior of a mathematical model that can do this. r 2011 American Chemical Society

Before discussing relevant prior work on droplet evaporation, we note two modeling studies that have examined skin formation in planar liquid films. Routh and Zimmerman considered particle distributions in a horizontal liquid layer evaporating at a constant rate.11 Bulk fluid flow was neglected. As evaporation occurs, the liquidair interface drops in height, leading to higher particle concentrations at the interface relative to the substrate. The sharpness of the particle concentration profile depends on the relative importance of the timescales of evaporation and particle diffusion. Yiantsios and Higgins built upon this work by investigating the effects of convection due to thermally induced Marangoni flows as well as the role of rheology.12 As a consequence of evaporation, an initially planar film develops surface undulations, and the particle concentration along the liquidair interface is governed by the interplay among evaporation, diffusion, and Marangoni effects. For droplets laden with colloidal particles, Deegan et al.1,2,7 proposed a mechanism for particle transport based on the observation that the droplet contact line is fixed during a majority of the drying time. An outward flow inside the droplet is produced when the liquid removed by evaporation from the edge of the droplet must be replenished to satisfy the fixedcontact-line constraint. This outward flow carries particles to the edge of the droplet. To explore this mechanism, Deegan et al.1,2,7 developed a model for fluid flow within the droplet which assumed that (i) the droplet shape is always a spherical cap with a fixed Received: June 3, 2011 Revised: August 11, 2011 Published: August 11, 2011 11347

dx.doi.org/10.1021/la202088s | Langmuir 2011, 27, 11347–11363

Langmuir

ARTICLE

Figure 2. Sketch of evaporating droplet on a flat, solid substrate.

Figure 1. Image of a drying droplet of 1-μm-diameter polystyrene spheres (Thermo Scientific) suspended in water evaporating on a microscope slide under ambient conditions. The scale bar represents 1 mm, and the volume fraction of particles is approximately 0.001. The light circular-like line is the edge of the droplet, and the dark cloud in the center is the skin of particles. The image is taken at a relatively early time; at later times (not shown), the skin breaks up and the particles collect at the droplet edge to produce a coffee ring.

contact line, (ii) the flow in the droplet is primarily in the radial direction, and (iii) evaporation is rate limited by the diffusion of liquid molecules into the surrounding air. The model indeed predicts that particles will collect at the contact line, but more quickly than is actually observed. The discrepancy was attributed in part to the collection of particles at the dropletair interface (i.e., skin formation), which the model does not describe. Numerous publications have built upon the original model proposed by Deegan et al. to investigate the spatial patterns left behind by evaporating colloidal droplets.8,9,1322 Although a comprehensive review of these is beyond the scope of this article, we highlight the particularly noteworthy observation that strong counter flows can be driven by thermal Marangoni effects. These flows can reverse the coffee-ring effect and produce a particle distribution with high concentrations in the center of the drop.9,14,15 The direction of the Marangoni flows is determined by the ratio of the substrate thermal conductivity to liquid thermal conductivity.23 Although prior work has contributed much to our fundamental understanding of the coffee-ring effect, it has not emphasized the phenomenon of skin formation. A common feature in prior work on drying droplets of colloidal suspensions is the assumption that the evaporation process is diffusion-limited. In general, the mass flux at the liquidair interface is determined by both the phase-transition process (i.e., the transfer of molecules across the interface) and the removal of vapor into the surrounding air by diffusion (assuming convection-free air).2427 When the rate-limiting step is the phase-transition process, another evaporation regime emerges where the diffusion of vapor away from the interface is rapid and the concentration of vapor in the air is uniform.24 We refer to this evaporation regime as transfer-rate-limited. Here, the evaporative flux can be determined by a classical Hertz Knudsen relation (derived from ideas in the kinetic theory of gases), and the dynamics of evaporation are governed by thermal diffusion in the liquid and solid. The transfer-rate-limited regime allows for decoupling of the liquid dynamics from the air dynamics because the latter can be ignored. This approach has been used to model spreading evaporating droplets of pure liquid on heated substrates.2830

As noted above, the focus of the present work is to examine the phenomenon of skin formation. However, rather than examining the regime of diffusion-limited (i.e., slow) evaporation (which appears to characterize many of the experiments discussed earlier), we consider the opposite limit of transfer-rate-limited (i.e., fast) evaporation. There are several reasons for taking this approach. First, most previous work on the coffee-ring effect has focused on diffusion-limited evaporation, but little seems to be known about what happens in the opposite limit. Second, the model problem that we will examine contains a precursor film at the droplet edge rather than pinned contact lines and so is relevant to experimental situations in which a thin adsorbed film (∼100 nm or less) is in equilibrium with the solid and vapor phases.29 One such practical situation includes a colloidal droplet drying on a heated substrate in a saturated vapor environment. As a consequence, the droplet will spread (due to capillary pressure gradients) as it evaporates, a situation that has also been overlooked in prior work (which emphasizes pinned contact lines). Third, as we will discuss later, the model problem that we consider can be treated consistently within the lubrication approximation, whereas the problem of diffusion-limited evaporation of a droplet with pinned contact lines and a spherical-cap shape (considered in previous work) cannot. We next present our mathematical model and the numerical methods applied to it, followed by results and discussion. Also different from prior work is that we consider the role of rheology by allowing the droplet viscosity to depend on the particle concentration. By describing the transport of the colloidal particles with the full convectiondiffusion equation, we are able to capture depthwise gradients in particle concentration and thus describe skin formation.

’ MATHEMATICAL MODEL Our model consists of an axisymmetric droplet of a Newtonian solvent containing colloidal particles, with the droplet in contact with ambient air and precursor films (Figure 2). The droplet, which has an initial radius R0 and height H0 in cylindrical coordinates (r0 , z0 ), rests on a solid substrate at z0 = 0. Its interface with the air is located at z0 = h0 (r0 , t0 ), and the associated velocity field is denoted by v0 = (v0 r, v0 z). Except for some of the material properties given in Table 1, the notation is such that primes indicate dimensional variables. To isolate the most basic effects of droplet rheology, we will model the colloidal particles as hard spheres. With this assumption, an increase in particle concentration will lead to an increase in viscosity, which we will describe with the KriegerDougherty relationship !2 ϕ ð1Þ μ ¼ μ0 1 ϕm where μ0 is the viscosity of the pure suspending liquid, ϕ is the particle volume fraction, and ϕm is the volume fraction at maximum packing31 (equal to 0.64 for hard spheres). Because the particle volume fraction is proportional to the particle 11348

dx.doi.org/10.1021/la202088s |Langmuir 2011, 27, 11347–11363

Langmuir

ARTICLE

concentration, the droplet equations of motion will be coupled to a particle transport equation through the viscosity relation (eq 1). Assuming that the droplet is small enough that gravitational effects can be neglected, the equations governing the flow inside the droplet are 0

0

∇ 3v ¼ 0 0 ∂v 0 0 0 F þ v v ∇ ¼ ∇ 0 3 T0 3 ∂t 0

ð2Þ

T0 ¼ p0 I þ μð∇0 v0 þ ∇0 v0T Þ

ð4Þ

ð3Þ

where k denotes the thermal conductivity and cp denotes the specific heat of the suspending liquid. On the solid substrate, z0 = 0, we assume a constant temperature, a no-slip condition, and impermeability such that 0

T 0 ¼ T 0 , vr0 ¼ 0, and vz0 ¼ 0

where T 0 is room temperature. At the interface, z0 = h0 (r0 , t0 ), we must impose boundary conditions to balance the mass, normal stress, tangential stress, and heat flux. We have adapted the one-sided evaporation model of Burelbach et al.,32 which decouples the dynamics of the ambient air from the liquid. The model is based upon the assumptions that the ambient air density, viscosity, and thermal conductivity are much smaller than those in the liquid. The mass balance is given by 0

J 0 ¼ Fðv0 vI Þ 3 n0

ð7Þ

where J0 is the evaporative mass flux through the liquidair interface, v0 I = (0, ∂h0 /∂t0 ) is the velocity of the interface, and ð ∂h0 =∂r 0 , 1Þ n0 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð∂h0 =∂r 0 Þ2 þ 1

ð8Þ

is the outward unit normal vector to the interface. The normal stress balance across the liquidair interface gives 0

pa n0 3 T0 3 n0 ¼ kσ Π0

ð9Þ

where p0 a represents the ambient air pressure, and the mean curvature k is given by k ¼ ∇ 3 ½ð1 þ j∇ h j Þ

0 0

∇h

ð10Þ

The normal stress balance includes a contribution from capillarity as well as a disjoining pressure, Π0 , modeling van der Waals interactions. In the tangential stress balance, t0 3 T0 3 n0 ¼ t0 3 ∇0 σ

ð12Þ

We assume that the surface tension of the droplet is linearly related to the interface temperature, T0 I = T0 (r0 , h0 (r0 , t0 ), t0 ), by 0

0

ð11Þ

ð13Þ

Parameter σ0 is the surface tension of the pure suspending liquid at room temperature T 0 0, and γ t dσ/dT 0 . For common liquids, γ > 0; therefore, the fluid on the interface flows from hot to cold regions. Because the fluid inside the droplet is viscous, it is dragged along with the surface flow, a phenomenon known as the thermocapillary effect. The heat-flux balance is given by kn0 3 ∇0 T0 ¼ Lm J 0

ð14Þ

where Lm is the latent heat of vaporization per unit mass. To close the fluid model, we must specify a constitutive equation for the evaporative mass flux J0 . Because we assume transfer-rate-limited evaporation, the evaporative mass flux is determined by the HertzKnudsen relation. We follow the work of Ajaev et al. and add a correction due to Laplace pressure so that 0

0

0

KJ 0 ¼ Rðp0 pa Þ þ ðT I T S Þ

ð6Þ

0

0 0 2 1=2

ð1, ∂h0 =∂r 0 Þ t0 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð∂h0 =∂r 0 Þ2 þ 1

σ ¼ σ 0 γðT I T 0 Þ

where t0 denotes time, F denotes density, and p0 denotes pressure. Because the droplet is evaporating, we must also describe the temperature T 0 in the droplet 0 ∂T 0 0 0 Fcp ∇ þ v T ð5Þ ¼ k∇02 T 0 3 ∂t 0

0

the unit tangent vector to the interface is

ð15Þ

where K and R are constants described in the next paragraph, and T 0 S is the saturation temperature.29,30,33 With this model, evaporation is driven by both deviations in the interfacial temperature as well as pressure. Moreover, this form of the mass flux allows for the formation of a precursor film where evaporation and van der Waals forces balance each other. The parameter K in eq 15 is derived from kinetic theory and measures the relative importance of kinetic effects. It is given by 0 3=2 pffiffiffiffiffiffiffiffiffiffiffi TS 2πR ð16Þ K¼ ωFv Lm where R is the universal gas constant divided by the molar mass of pure solvent vapor, ω is the accommodation coefficient describing the probability of phase change, and Fv is the density of the pure solvent vapor. We take ω = 7.5 10 1.34 Parameter R, explained in more detail below, is chosen such that a desired precursor film thickness is achieved. To simplify the governing equations, we assume the droplet is thin so that ε = H0 /R0 , 1, thereby allowing the application of the lubrication approximation.35 In prior studies on evaporating droplets,1,9,14,36 a further simplification is made because of the fact that the capillary number, given by Ca ¼

μ0 V 0 9:28 1010 ¼ σ0 R 0 ðmÞ

ð17Þ

where V 0 is the characteristic velocity, is small. Following Ajaev,29 the characteristic velocity scale is determined from the interfacial mass and energy balances so that V 0 = kT 0 S/ FLmR0 . Deegan and co-workers considered experiments with 5 10 4 m e R0 e 0.15 m leading to estimates for Ca of 1.8 10 6 or less.1 On the basis of this estimate, they argue that surface tension effects dominate over viscous effects and that the droplet has the shape of a spherical cap (assuming gravity is negligible). However, as discussed by Ajaev and co-workers,29,30 11349


Langmuir

ARTICLE

the effects of surface tension and viscosity become equally important for evaporating droplets if ε = H0 /R0 = O(Ca1/3). Thus, to explore this limit, we consider asymptotic expansions in powers of Ca1/3 and take ε = Ca1/3 so that Ca = Ca/ε3 = 1. We introduce the following scalings to nondimensionalize the governing equations and boundary conditions: 0

0

0

0

1=3 0

0

ð18Þ

Taking the limit ε f 0 leads to the following leading-order equations: 1 ∂ ∂vz ðrvr Þ þ r ∂r ∂z ∂p ∂ μ ∂vr 0¼ þ ∂r ∂z μ0 ∂z

0¼

ð19Þ ð20Þ

∂p ∂z

ð21Þ

∂2 T ∂z2

0¼

ð22Þ

with 0

0

T0 TS 0 Θ, vr ¼ 0, vz ¼ 0 Ca2=3 T S

T ¼

ð23Þ

at z = 0 and ∂h ∂h þ vz ¼ J ∂r ∂t 1 ∂ ∂h A̅ p ¼ pa r 3 r ∂r ∂r h

ð24Þ

vr

ð25Þ

μ ∂vr ∂Σ ¼ μ0 ∂z ∂r

ð26Þ

∂T ¼J ∂z

ð27Þ

where M is a modified Marangoni number,

1=3 0

r ¼ R r, z ¼ Ca R z, h ¼ Ca R h, 0 v0r ¼ V 0 vr , V z ¼ Ca1=3 V 0 vz , t 0 ¼ R 0 =V 0 t, μ0 V 0 μ0 V 0 0 p0 ¼ 2=3 p, Π ¼ Π, Ca R 0 Ca2=3 R 0 0 0 T 0 T S ¼ Ca2=3 T S T, J 0 ¼ FV 0 Ca1=3 J

0¼

Solving the temperature problem gives T(r, z, t) = J(r, t)z + Θ. By eq 13, we find that 0 ∂Σ CaγT ∂T ∂h ∂T ∂ ¼ 0 S þ ð31Þ ¼ M ðJhÞ ∂r ∂r ∂z ∂r V μ0 ∂r T γ M ¼ S σ0

first introduced by Gramlich et al.37 for thin-film flow over topographical features and subsequently used by Ajaev.29 The radial velocity component is found by integrating the radial component of the momentum conservation equation (eq 20), applying the tangential stress balance (eq 26) at the interface, and implementing the no-slip condition (eq 23) at the substrate: Z z Z z ∂p μ0 ∂ μ0 ðz hÞ dz þ M ðJhÞ dz vr ¼ ∂r 0 μ ∂r 0 μ ð33Þ Because ∂p/∂z = 0, the pressure is simply stated by the normal stress condition (eq 25) A̅ 1 ∂ ∂h p ¼ pa r ð34Þ 3 r ∂r ∂r h Substituting the interface temperature TI = Jh + Θ and eq 34 into the dimensionless evaporative mass flux (eq 30) gives A̅ 1 1 ∂ ∂h Θδ r þ 3 J ¼ ð35Þ K̅ þ h r ∂r ∂r h The evaporative mass flux varies along the dropletair interface and is largest near the contact-line region. The vertical velocity component can be found by integrating the conservation-of-mass equation (eq 19) to yield Z z 1 ∂ ðrvr Þ dz ð36Þ vz ¼ 0 r ∂r Finally, we arrive at the evolution equation for the interface by considering the kinematic condition (eq 24), ∂h 1 ∂ þ ðrQ Þ ¼ J ∂t r ∂r

29

at z = h(r, t). Note that following Ajaev we model the disjoining pressure with the expression Φ¼

A̅ h3

where

ð28Þ

Z

Q ¼

h

vr dz

ð37Þ

ð38Þ

0

The dimensionless parameters in the above equations are Σ¼

ð32Þ

Ca σ A and A̅ ¼ V 0 μ0 6πσ 0 R 02 Ca4=3 1=3

ð29Þ

where A is the Hamaker constant. The dimensionless constitutive relation for J is KJ ̅ ¼ δðp pa Þ þ T I 0

where K = KFV /Ca

1/3

T

ð30Þ 0

S

0

and δ = Rμ0V /Ca

4/3 0

0

R T S.

It is important to note that the lubrication approximation with the no-slip boundary condition breaks down near the contact line if the contact line is evaporating but remains stationary (i.e., no precursor film). Specifically, the kinematic condition (eq 24) is violated at the contact line. With ∂h/∂t = 0 (stationary contact line), vr = 0 (no slip), and vz = 0 (no penetration), the kinematic condition gives the contradiction 0 = J > 0. Even with the inclusion of a Navier slip law, the theory still breaks down because the substrate and interfacial boundary conditions must hold simultaneously at the contact line. For example, in the case 11350


Langmuir

ARTICLE

of a Navier slip law, vr j z ¼ 0 ¼ β

∂vr ∂Σ ∂h j ¼ β jr ¼ R ¼ βMJ jr ¼ R ∂r ∂r ∂z z ¼ 0

where β is a positive slip coefficient and R* is the location of the contact line. Therefore, at the contact line, ∂h/∂t = 0, ∂h/∂r < 0, J > 0, vr < 0, and vz = 0. Thus, the kinematic condition (eq 24) is violated. The above observations suggest that evaporating droplets with pinned contact lines cannot be rigorously described within the lubrication approximation. Although Deegan et al.1,2 do apply the lubrication approximation, they do not explicitly check whether their velocity field satisfies momentum conservation. Hu and Larson9,14 also apply the lubrication approximation and assume the form h(r, t) = h(0, t)(1 r2) for the droplet shape. Without van der Waals forces, ∂p/∂r = 0 by eq 25. Therefore, the first term in eq 33 for vr will be zero. If μ = μ0, then vr can at most be a linear function of z, but Hu and Larson took it to be quadratic in z. The success of these prior models in explaining the experimental data must thus be viewed as fortuitous, perhaps because the errors incurred are relatively small. Nevertheless, to develop a model fully consistent with the lubrication approximation, we consider here the case where the contact line is not pinned but is instead replaced by a thin precursor film. As noted in the Introduction, this model is also relevant to experimental situations where a thin adsorbed film exists at the droplet edge. The droplet-profile evolution equation (eq 37) requires an initial condition and four boundary conditions. The initial condition is a droplet of constant curvature on a precursor film: ( h0 ð1 r 2 Þ þ heq for r e 1, ð39Þ hðr, 0Þ ¼ otherwise heq The uniform equilibrium thickness, heq = (δA/Θ)1/3, is chosen such that the evaporative mass flux is zero in eq 35. Note that along the flat precursor film all radial derivatives are zero and therefore the expression in the square brackets in eq 35 is zero when h = heq. We chose the parameter δ such that a desired equilibrium thickness is achieved; for the current study, heq = 0.025. The length of the computational domain, denoted by the parameter L, is chosen to ensure that the contact-line region (i.e., the region where the droplet meets the precursor film) stays sufficiently inside the computational domain. Hence, near r = L, the substrate is covered by a precursor film of thickness heq. We use L = 2. The boundary conditions specify symmetry and smoothness and are expressed by ∂h ∂3 h ∂h ∂3 h ð0, tÞ ¼ 0, 3 ð0, tÞ ¼ 0, ðL, tÞ ¼ 0, 3 ðL, tÞ ¼ 0 ∂r ∂r ∂r ∂r

ð40Þ

Next, we turn our attention to the particle dynamics. We assume that the diameter of a colloidal particle is much smaller than the characteristic height of the droplet; therefore, a continuum model for the particles is used. In particular, the particle concentration c0 (r0 , z0 , t0 ), the mass of particles per unit volume, is governed by the convectiondiffusion equation ∂c0 1 ∂ 0 ∂c0 ∂ ∂c0 0 0 0 þ v 3∇ c ¼ 0 0 r D 0 þ 0 D 0 ð41Þ r ∂r ∂z ∂t 0 ∂r ∂z

where D is the particle diffusivity. If we scale eq 41 by the inverse of the particle density Fp, then we arrive at an expression for particle transport in terms of the particle volume fraction ϕ = c0 /Fp. The particle volume fraction, introduced in the viscosity equation (eq 1), ranges from 0 to the maximum packing fraction of ϕm. As in Yiantsios and Higgins,12 we assume that D is given by the generalized StokesEinstein equation: ! 1:85ϕ 6:55 d D ¼ D0 ð1 ϕÞ ð42Þ dϕ 0:64 ϕ The coefficient D0 is the Einstein diffusivity given by D0 = kBT*/6πμ0a, where kB is the Boltzmann constant, T* is a characteristic temperature taken to be T 0 S, and a is the particle radius. We consider 1-μm-diameter spheres and therefore a = 5 10 7m. The particle transport equation is supplemented with no-flux boundary conditions to reflect the assumption that particles do not leave the droplet: ∂c0 ∂c0 ∂c0 ¼ 0 at z0 ¼ 0, 0 ¼ 0 at r 0 ¼ 0, 0 ¼ 0 at r 0 ¼ R0 L ∂z0 ∂r ∂r

ð43Þ

There is also a condition on the diffusive flux at the liquidair interface: c0 ðv0 v I Þ 3 n0 ¼ Dð∇0 c0 Þ 3 n0 at z0 ¼ h0 0

ð44Þ

This condition plays an important role in skin formation because it leads to the collection of particles along the evaporating interface. Using the scalings introduced in eq 18, the leading-order particle transport equation reads ∂ϕ 1 ∂ ∂ 1 ∂ D ∂ϕ þ ðrvr ϕÞ þ ðvz ϕÞ ¼ ð45Þ ∂t r ∂r ∂z Pe ̅ ∂z D0 ∂z where Pe = Ca1/3Pe and Pe = V 0 Ca1/3R0 /D0. The leading-order boundary conditions are given by ∂ϕ ∂ϕ ∂ϕ ðr, 0, tÞ ¼ 0, ð0, z, tÞ ¼ 0, ðL, z, tÞ ¼ 0 ∂z ∂r ∂r

ð46Þ

1 D ∂ϕ ∂h ∂h ¼ ϕ vr þ vz at z ¼ h Pe ∂r ∂t ̅ D0 ∂z

ð47Þ

The initial particle concentration is assumed to be spatially uniform inside the droplet and zero in the precursor film. For the parameter values we explore, the precursor film is assumed to be thinner than the diameter of a particle. Thus, the initial function is given by 8 > for r e 0:7 < ϕ0 ð48Þ ϕðr, z, 0Þ ¼ bðrÞ for 0:7 < r < 1 > :0 otherwise where ϕ0 ∈ [10 3, 0.3] and b(r) is a fifth degree polynomial whose coefficients are chosen so that ϕ (r, z, 0) has continuous first and second derivatives. In summary, the evolution of the droplet shape is found by simultaneously solving eq 37 in the 1D stationary domain 0 < r < L and eq 45 in the 2D moving droplet-shape domain 0 < r < L, 11351


Langmuir

ARTICLE

Table 1. Material Parameters for Water at Room Temperaturea T 0 S (K) T

0

0

Table 2. Estimates of the Nondimensional Parameters Calculated Using the Material Parameter Values in Table 1

2.9799 102

Ca

9.28 l010[R0 (m)]1

2.9800 10

A

8.15 106 [R0 (m)]2/3

9.9700 10

M

6.95 101

2

(K)

F (kg/m ) 3

2

Fv (kg/m )

2.0000 10

K

1.21 101[R0 (m)]2/3

R(J/(kg K)) ̅ μ0 (Pa s)

2

4.6115 10

δ

6 76 l0°[R0 (m)]4/3

4

Θ

3. 53 l0°[R0 (m)]2/3

1

̅Pe

1.43 101[R0 (m)]2/3

3

2

9.0000 10

k (W/(m K))

6.0700 10

cp (J/(kg K)) Lm (J/kg)

4.1800 10 2.4450 106

K (K m2 s/kg)

7.5500 10°

σ0 (N/m)

7.2000 102

γ (N/(mK))

1.6800 104

A (J)

1.0000 1017 2

D0 (m /s)

3

4.8400 1013

a

Many of these values are from published sources,18,32 but values for A and T 0 S are larger than published data in order to achieve a larger precursor film thickness of O (103). This larger thickness allows for adequate resolution of the concentration field given the grid size in our finite-difference calculations.

0 < z < h. The two equations are coupled by the radial and vertical velocities, which both depend on droplet shape h and particle volume fraction ϕ. Estimates of the nondimensional parameters are given in Table 2.

’ NUMERICAL SCHEME A finite-difference-based method of lines is implemented to solve the coupled system of nonlinear partial differential equations given in eqs 37 and 45. In the method of lines, space and time are discretized independently, thereby replacing the original coupled partial differential equations with a large system of coupled ordinary differential equations. The droplet evolution equation (eq 37) is solved in MATLAB. The stationary domain, 0 < r < L, is discretized by an equally spaced grid, and the spatial derivatives are approximated by centered finite differences. To avoid complete coupling, we implement a semi-implicit backward Euler time-integration scheme. In particular, the integrals of the viscosity, where μ = μ(ϕ), in the expression (eq 38) for Q are treated explicitly. The particle transport equation (eq 45) is more challenging to approximate because its domain is moving. To do so, we implement the method of moving overset grids in Overture. Overture is an object-oriented environment for solving partial differential equations on overset time-dependent grids.3840 The framework, developed at Lawrence Livermore National Laboratory, includes support for grid generation, difference operators, boundary conditions, database access, and graphics. Moving overset grids have been recognized to be an attractive approach to handling problems with complex moving geometries. Specifically, they have been successfully used to solve flow problems with deforming boundaries.41 By definition, an overset grid (also known as composite overlapping grids or Chimera grids) is a collection of logically rectangular, curvilinear component grids covering a domain and overlapping where they meet.42 Each component grid is defined by a smooth mapping from parameter space (m, n) ∈ [0,1] to physical space (r, z), ðr, zÞ ¼ Gðm, n, tÞ

The component grids communicate through interpolation. Figure 3 displays the droplet-shaped grid built by overlapping grid generator Ogen in Overture.38 The overlapping grid is composed of a boundary-fitting curvilinear grid at the interface with two background Cartesian grids for the remaining area. On the overlapping grid, the spatial derivatives are approximated by second-order curvilinear finite differences. By curvilinear, we mean that approximations to the derivatives of ϕ(r, z, t) with respect to (r, z) are formed on the unit square (m, n) by an application of the chain rule. For example, ∂ϕ ∂m ∂ϕ~ ∂n ∂ϕ~ ¼ þ ∂r ∂r ∂m ∂r ∂n

ð50Þ

where ϕ~ðm, n, tÞ ¼ ϕðGðm, n, tÞ, tÞ

ð51Þ

and ∂m ∂n , ∂r ∂r are inverse vertex derivatives computed from the mapping of G. The curvilinear finite differences are obtained by discretizing the (m, n) derivatives in eq 50 by standard centered second-order finite differences.38 The integrals of the viscosity appearing in the expressions for vr and vz are approximated on a single Cartesian grid to avoid the complexities associated with the curvilinear component grids. To do so, μ is defined by evaluating an interpolant of μ built on the droplet-shape grid. If the point on the Cartesian grid is outside the droplet-shape grid, then the value is assigned by extrapolation. Next, the integrals of μ are calculated on the single Cartesian grid by applying the trapezoidal rule. Finally, the integral values on the droplet-shape grid are assigned by evaluating an interpolant built on the single Cartesian grid. To handle the moving domain, the evolution equation is solved in a frame moving with the grids. Specifically, the boundary-fitted grid moves according to the droplet dynamics while the background Cartesian grid remains fixed. When the moving grid support is used in Overture, at each time step a new overlapping composite grid is efficiently generated by Ogen.43 The time integration is done using a semi-implicit backward Euler method. To avoid nonlinearities, we treat the radial and vertical velocities as well as the diffusivity explicitly. For numerical stability, we retain the lower-order radial diffusion term in the particle transport equation44

∂ϕ 1 ∂ ∂ 1 ∂ D ∂ϕ ∂ D ∂ϕ þ ðrvr ϕÞ þ ðvz ϕÞ ¼ þ ε2 ∂t r ∂r ∂z Pe ∂r D0 ∂r ̅ ∂z D0 ∂z

ð52Þ

ð49Þ 11352


Langmuir

ARTICLE

Figure 3. Overlapping grid of the droplet-shaped domain in the nondimensional coordinate system. The three different component grids are shown in red (the boundary-fitting grid), blue, and green.

’ RESULTS AND DISCUSSION We begin by examining skin formation in an evaporating droplet with a constant viscosity, μ = μ0, and a constant diffusivity of particles, D = D0. We then systematically incorporate various additional effects to understand their influence on the particle distribution. Without the influence of rheology, the droplet and particle dynamics are independent. Thus, we start by describing typical droplet behavior without Marangoni stresses as shown in Figure 4.29 In what follows, dimensional variables are plotted in all figures to allow for easy comparison to prior work, and the droplet has an initial radius of R0 = 2 10 3 m. Unless otherwise stated, the parameters are taken to be the values given in Table 2. Initially, the droplet undergoes surface-tensiondriven spreading, as shown in Figure 4a. This is followed by a time period in which the droplet recedes, shown in Figure 4b,

where a nearly uniform curvature droplet (away from the contact-line region) slowly evaporates. The droplet is continuously losing mass by evaporation. Figure 4c,d shows the evaporative mass flux during the spreading and receding phases, respectively. The evaporative mass flux reaches a maximum near the contact-line region and then rapidly decays as the precursor film is approached because of the van der Waals forces. The addition of thermocapillarity, controlled by parameter M, hinders the surface-tension-driven spreading. Figure 4e shows the droplet evolution with M = 50. Marangoni stresses generate a surface flow from the contact-line region (where the temperature is higher) to the apex of the droplet (where the temperature is lower). This flow acts to slow down the droplet spreading and keeps the droplet upright. As a result of the change in the droplet shape, the evaporative mass flux is smaller (Figure 4f); therefore, it takes longer for the droplet to evaporate. Overall, the droplet dynamics are governed by a competition among capillary spreading, evaporation, and thermocapillarity. We now consider particle dynamics without any coupling to the droplet dynamics. A typical evolution of the particle distribution and velocity field, with M = 0, Pe = 10, and ϕ0 = 0.1, is shown in Figure 5. Snapshots of the particle concentration contours are displayed on the left, and the corresponding velocity fields are shown on the right. The arrows in the velocity field plots illustrate the direction of flow and are covered by contours of the velocity amplitude, v0 (in units of m/s). Therefore, there is no meaning attached to the length of the arrows. When looking at the figures, it is important to examine the color bar associated with each plot because the scale changes from figure to figure. As the droplet begins to evaporate, a front of particles collects at the dropletair interface (Figure 5ad). Because a peak in the evaporative mass flux occurs near the contact-line region (Figure 4c), the maximum particle concentration in the particle front occurs there. The particle concentration near the interface continues to grow as the droplet continues to loss mass. Simultaneously, the particles in the front both diffuse vertically toward the substrate and convect with the surface flow. In ) )

The simulation proceeds until the particle concentration reaches within 103 of the maximum volume fraction ϕm at some point in the domain. The Matlab code is validated by assuming that the droplet contains no particles (i.e., ϕ = 0 and μ = μ0). Craster et al. employed the same droplet evolution model (eq 37) to study the pinning, retraction, and terracing of evaporating droplets containing nanoparticles.10 In their work, they first considered the dynamics of the droplet in the absence of nanoparticles. We verified our Matlab code by reproducing the droplet profiles published in Figure 1 of their work.10 The Overture code is verified by the method of analytical functions. In the method of analytical functions, a forcing function is chosen such that the true or exact solution is known. It is important that the partial derivative of the exact solution with respect to r vanishes at r = 0 so that the particle transport equation can be evaluated at the boundary. With the true solution in hand, the numerical errors are measured and tracked. In addition, the 1D diffusion model developed by Routh and Zimmerman was used as a test, and the Overture code was able to reproduce Figures 2 and 3, published in ref 11.

11353


Langmuir

ARTICLE

Figure 4. Droplet dynamics in the absence of particles with parameter values given in Table 2 for R0 = 2 10 3m. (a, b) Evolution of droplet height in the spreading (by increments of 5.39 s) and receding phases (by increments of 26.95 s), respectively, when M = 0. The arrow points in the direction of increasing time, and the plots in b are for later times than those shown in a. (c, d) Evaporative mass flux in the spreading and receding phases, respectively, when M = 0. The plots are shown for times corresponding to those used in a and b. (e) Evolution of droplet height by increments of 53.90 s when M = 50. Note that the spreading phase is suppressed. (f) Evaporative mass flux corresponding to e. Here and elsewhere, the units of h0 and r0 are m, and those of J0 are kg/(m2 s).

general, the vertical diffusion causes the particle front to spread inside the droplet, and the behavior of the surface flow depends on the strength of Marangoni stresses. Therefore, the formation of a skin is determined by a competition among diffusion, convection, and evaporation. Inside the droplet, the particles convect toward the contactline region and begin to accumulate there, which is similar to the coffee-ring effect. The particles inside the droplet are carried to the contact-line region by surface-tension-driven spreading (i.e., flows driven by capillary pressure gradients). Note that this particle-transport mechanism is different from the one studied by Deegan and co-workers1,2,7 in their work on the coffee-ring effect. Deegan and co-workers consider the case of a sphericalcap droplet with pinned contact lines subject to diffusion-limited

evaporation. There, particle transport occurs because of flows generated by the evaporative mass flux (rather than by capillary pressure gradients). Because liquid is evaporating near the contact line, flow must occur toward the contact line in order to keep it pinned. For the problem of a spreading droplet that we consider, the maximum in the particle concentration eventually switches from the particle front near the contact-line region to the contact-line region itself (Figure 5ad), leading to a coffeering-like effect. When the particle concentration in the contactline region reaches the maximum volume fraction, the simulation is stopped. Effects of Peclet Number. The Peclet number measures the importance of convective transport to diffusive transport. Therefore, at low Peclet numbers, diffusion dominates and particle 11354


Langmuir

ARTICLE

Figure 5. Results from a computation with Pe = 10, M = 0, ϕ0 = 0.1, μ = μ0, and D = D0. (ad) Contours of the particle volume fraction at various times during droplet drying. (IIV) Corresponding velocity profiles. The computation is stopped when the particle volume fraction is within 10 3 of ϕm = 0.64. 11355


Langmuir

ARTICLE

Figure 6. Results from a computation with Pe = 50, M = 0, ϕ0 = 0.1, μ = μ0, and D = D0. (ad) Particle concentration contours inside the droplet at various times. (IIV) Comparison of particle concentrations at the interface obtained with different values of Pe. The center of the droplet is located at r0 = 0, and the contact-line position is the largest radial position plotted. The results presented are for Pe = 10 (O) and Pe = 50 (*). 11356


Langmuir

ARTICLE

Figure 7. Results from a computation with Pe = 50, M = 50, ϕ0 = 0.1, μ = μ0, and D = D0. (ad) Contours of the particle volume fraction at various times during droplet drying. (IIV) Corresponding velocity profiles. The speed of the liquid inside the precursor film is effectively zero. 11357


Langmuir motion is strongly driven by concentration gradients. The evaporating droplet initially forms particle concentration gradients near the dropletair interface. These vertical concentration gradients are then smoothed out by diffusion, causing the particle front to spread vertically toward the substrate. As the Peclet number increases, diffusion weakens and vertical particle concentration gradients become more pronounced. This is seen in Figure 6ad, where Pe = 50 and all other parameters are the same as in Figure 5. Consequently, the initial evaporationinduced particle front is more closely packed (Figure 6ad) and has a larger concentration at the interface (Figure 6IIV). In addition, the particle dynamics are more influenced by the fluid flow, which carries particles both inside the droplet and on the droplet surface toward the contact-line region. Therefore, the maximum volume fraction is reached in 45.22 s, whereas for Pe = 10 the maximum volume fraction is attained at 48.40 s. Note the inversion in the particle concentration profiles at the interface (Figure 6IIV). At early times, the concentration is nearly uniform along most of the interface before dropping sharply in the contact-line region. As time progresses, the concentration builds up in the contact-line region and drops over the rest of the interface. To summarize, weak diffusion promotes the formation of a skin. This observation is qualitatively consistent with the experimental results of Deegan and co-workers,1,2,7 who saw skin formation only in experiments with micrometer-sized particles.1 Because the diffusivity, D0, scales inversely with the particle radius, larger particles have smaller diffusivities. Once a skin is formed, its dynamics are then determined by the convection. Effects of Marangoni Number. When the intensity of the Marangoni stresses is increased, the fluid flow within the droplet develops a recirculation with the direction of the surface flow toward the apex of the droplet. The strong Marangoni flow prevents initial droplet spreading (Figure 4e). Figure 7 shows the particle concentration contours for Pe = 50 and M = 50. The initial local maximum of concentration for the interfacial particle front that forms near the contact-line region now persists throughout the simulation. In particular, the recirculation flow pulls the particles in the front toward the apex of the droplet and creates a localized region of high particle concentration. Competing with this are the convective flows inside the droplet, which collect particles in the contact-line region. Even though the surface flow becomes stronger as time proceeds, the convective flows inside the droplet eventually dominate and the local maximum concentration in the contactline region becomes larger than the local maximum on the interface. Soon after, the maximum volume fraction is reached in the contact-line region and the simulation is stopped. Figure 8 compares the final particle concentrations along the interface for Pe = 50 and different Marangoni numbers. In general, the convective forces pull either the local region of high particle concentration that forms near the contact-line region toward the apex of the droplet (M = 50) or the contact line (M = 0). For M = 0, the final concentration profile monotonically increases until the maximum volume fraction is reached in the contact-line region. For M = 50, there is a local maximum in the concentration profile away from the contact-line region. This local maximum would correspond to the formation of a skin. Thus, we can conclude that the skin forms with weak diffusion and persists with a strong Marangoni flow. Effects of Initial Particle Concentration. We explore the influence of the initial particle concentration on the particle

ARTICLE

Figure 8. Comparison of the final particle concentrations (end of the simulation) at the interface obtained with different values of M. The center of the droplet is located at r = 0, and the contact-line position is the largest radial position plotted. The results presented are for M = 0 (*) and M = 50 (O).

dynamics. Recall that when μ = μ0 and D = D0 the droplet dynamics are not influenced by the particle concentration. Therefore, a smaller initial particle concentration allows the numerical simulation to run longer. Figure 9 shows the particle dynamics with ϕ0 = 0.001. The evolution of the particle concentration in Figure 9a,b mimics the ϕ0 = 0.1 dynamics in Figure 7ad, except that the scales for the color bars are different. As time progresses, the surface flow grows in strength and pulls the particles toward the apex. Consequently, the local region of high particle concentration becomes the location of the maximum particle concentration and moves toward the apex. This localized region of high particle concentration eventually reaches the apex of the droplet. Thus, given enough time, the Marangoni flows can pull the skin that forms near the contact-line region to the apex of the droplet. Effects of Rheology. Rheology influences the evolution of the droplet as well as the particle dynamics. In general, as the particle concentration increases, the viscosity increases before diverging at the maximum volume fraction. This causes the fluid flow to slow down, especially near the interface. Figure 10 shows results from a simulation with the same parameters as in Figure 7 but with a concentration-dependent viscosity. A comparison of Figures 10IIV and 7IIV shows that the velocity along the interface is slower when the viscosity depends on the particle concentration. For example, in Figure 7II the maximum speed is 0.001 m/s, whereas in Figure 10II it is 0.0005 m/s. Therefore, the Marangoni effect is dampened and the weakened recirculation flow greatly influences the final particle concentration. In particular, the final profile (Figure 10d) resembles what is seen when M = 0 (Figure 8), where the profile monotonically increases until the maximum volume fraction is reached in the contact-line region. As a consequence, the skin does not persist throughout the simulation. In the case of a droplet with pinned contact lines undergoing diffusion-limited evaporation, Hu and Larson discussed the need to suppress Marangoni flows in order to produce a coffee ring in drying water droplets.14,15 They showed both experimentally and theoretically that surface contaminants could weaken the thermally induced Marangoni flows. In the case of a spreading droplet undergoing transfer-rate-limited evaporation considered here, 11358


Langmuir

ARTICLE

Figure 9. Results from a computation with Pe = 50, M = 50, ϕ0 = 0.001, μ = μ0, and D = D0. (ad) Contours of the particle volume fraction at various times during droplet drying. (IIV) Corresponding velocity profiles. The speed of the liquid inside the precursor film is effectively zero. 11359


Langmuir

ARTICLE

Figure 10. Results from a computation including the influence of rheology. (ad) Contours of the particle volume fraction at various times for Pe = 50, M = 50, ϕ0 = 0.1, μ = μ(ϕ), and D = D0. (IIV) Corresponding velocity profiles. The speed of the liquid inside the precursor film is effectively zero.

11360


Langmuir our results indicate that rheology is another mechanism that can dampen thermally induced Marangoni flows. The intensity of the rheological damping of the Marangoni flow depends on the initial particle concentration. In general, the rheological influence is delayed with smaller initial particle concentrations because it takes a longer period of time for the particle concentration to reach a level where the viscosity has

Figure 11. Comparison of the final particle concentrations (end of the simulation) at the interface obtained with and without the influence of rheology. The remaining nondimensional parameters are Pe = 50, ϕ0 = 0.001, M = 50, and D = D0. The results presented are for μ = μ0 (O) and μ = μ(ϕ) (*).

ARTICLE

a significant effect. In particular, Figure 11 compares the final particle concentration profiles for ϕ0 = 0.001 with both concentration-independent and concentration-dependent viscosities. Compared to Figure 10, where ϕ0 = 0.1, the weakening of the Marangoni flow is not significant enough to interfere with the development of the interfacial skin away from the contactline region. At later times, the diverging viscosity does slow down the movement of the skin, causing the peak of the interfacial particle concentration to be shifted toward the contact-line region and the concentration profile not to be as steep compared to the case of constant viscosity. The dynamics of the droplet, including the contact-line position, are also influenced by the rheology. Figure 12a,b shows the contact-line position as a function of time for Pe = 50, either M = 0 or M = 50, and different initial particle concentrations. We define the contact-line position to be the maximum in the interfacial curvature as done in Ajaev.29 With M = 0, the initial droplet dynamics are dominated by capillary spreading. As a result, when the particles begin to collect in the contact-line region, the droplet spreading is slowed down by the diverging viscosity. Because the number of particles arriving at the contactline region depends on the initial particle concentration, droplet spreading is slower with increasing initial concentration. In contrast, when M = 50, the effects of rheology are noticeable only when the initial particle concentration is increased to ϕ0 = 0.3. This behavior occurs because the smaller the initial particle concentration, the longer it takes for the concentration in the contact-line region to influence the viscosity significantly.

Figure 12. Comparison of the dynamics of the droplet with Pe = 50, either M = 0 or M = 50, different initial particle concentrations, and μ = μ(ϕ). Contact line position vs time for (a) M = 0 and (b) M = 50. The discontinuities in the contact line position vs time plots are numerical artifacts arising from the computational grid. Evolution of droplet height for (c) M = 0 and (d) M = 50. 11361


Langmuir By the time that happens, the contact-line dynamics are dominated by evaporation. Figure 12c,d shows a comparison of the droplet shape at 15 s. When M = 0, the capillary spreading is hindered with increasing initial particle concentration, causing the apex of the droplet to be larger and the contact line to be shifted to the left. In contrast, for M = 50, the apex of the droplet becomes smaller with increasing initial particle concentration. The Marangoni flow along the surface of the droplet is weakened by the growing viscosity; therefore, the droplet is not as upright. Effects of Diffusivity. The diffusivity D given by the generalized StokesEinstein equation (eq 42) is a nonmonotonic function. Initially, as the particle concentration increases, the diffusivity decreases because of hindrance effects. However, as the concentration continues to increase, the diffusivity begins to grow, eventually diverging at ϕm. In the exploration of the effects of the Peclet number, diffusion was found to play an important role in the formation of the interfacial particle front. As the effects of diffusion are dampened, a more compact particle front develops. With the nonconstant diffusivity, the formation of the interfacial particle front strongly depends on the particle concentration. For brevity, we simply summarize the key findings below. When 0 < ϕ < 0.3 or 0.57 < ϕ < ϕm, D/D0 > 1 and the diffusivity suppresses the skin formation. However, for 0.3 < ϕ < 0.57, D/D0 < 1 and the skin formation is accelerated. Because the diffusivity begins to increase rapidly when ϕ > 0.48, a sharp transition in diffusivity develops near the interface. Consequently, at later times, the numerical calculations cannot be trusted because the spatial approximations of the diffusivity derivatives are very poor.

’ CONCLUSIONS We have examined the problem of evaporation in a spreading droplet laden with colloidal particles. As the droplet evaporates, particles begin to accumulate at the dropletair interface, creating a particle front that is most concentrated in the contactline region because of the larger evaporation rates there. If the Peclet number is increased (i.e., diffusion is weakened), then the particle front becomes more compact and concentrated. In the early stage of this skin formation, the eventual location of interfacial particles is then determined by the Marangoni number. If the Marangoni number is sufficiently small, then surfacetension-driven spreading collects the particles along the interface in the contact-line region. If the Marangoni number is sufficiently large, then thermocapillarity generates a recirculation flow that competes with capillary spreading and drives particles toward the apex of the droplet. Therefore, skin formation is promoted by weak diffusion (large Peclet numbers) and a strong Marangoni flow (large Marangoni numbers). Lower initial particle concentrations also promote skin formation because they allow more time for Marangoni stresses to act before the maximum volume fraction is reached. The addition of rheological effects slows down the dynamics. The increases in viscosity due to the particle concentration are most noticeable along the interface, where the particle front develops, and in the contact-line region. As a result, the Marangoni flow is dampened and skin formation is hindered. Rheological effects also slow down the rate at which the droplet spreads, and this effect is stronger at smaller Marangoni numbers because of the weaker thermocapillary flows that tend to carry particles to the droplet apex. There is little effect on the receding

ARTICLE

rate because the droplet receding phase is governed by evaporation rather than viscosity. As noted in the Introduction, prior modeling work on the coffee-ring effect has not emphasized the phenomenon of skin formation, focused on diffusion-limited evaporation, and considered pinned contact lines. In contrast, the present work considers the opposite limit of transfer-rate-limited evaporation and replaces the contact line with a precursor film, thereby allowing the droplet to spread. The solution of the full convectiondiffusion equation allows us to describe vertical concentration gradients and thus model skin formation. The use of the precursor film makes the model most relevant to experimental situations in which a thin adsorbed film exists at the droplet edges. In addition, our model accounts for rheological effects, which have largely been overlooked in prior work. Although our work has considered transfer-rate-limited evaporation, some of the phenomena we observe are very similar to what is seen in the diffusion-limited evaporation regime. For example, we observe the formation of a coffee-ring and the development of a surface cap of particles at the droplet apex. Our model may thus provide guidance in the development of more rigorous and accurate models for droplet evaporation problems regardless of the evaporation regime. We note that accounting for diffusion-limited evaporation in our model would likely require solving for a concentration field in the gas phase simultaneously with the droplet and particle dynamics, a considerably more involved calculation. As discussed earlier, accounting for pinned contact lines cannot be done consistently within the lubrication approximation and so would require solutions of the full equations of motion or perhaps a matched asymptotic expansion, regardless of the evaporation regime. Finally, we note that there are a number of phenomena that are not accounted for in our model that may be relevant to various experimental situations. These include Marangoni effects due to the adsorption of particles at the dropletair interface, particle adsorption at the substrate, and the fluidsolid phase transition that occurs as the packing fraction approaches the maximum value. Nevertheless, as noted above, the present model accounts for several features not emphasized in prior work and may thus serve as a basis for the development of more advanced models that account for these complex but potentially important phenomena.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Present Addresses

§ School of Mathematical Sciences, Rochester Institute of Technology, Rochester, NY 14623. E-mail: [email protected].

’ ACKNOWLEDGMENT K.L.M. is grateful to Bill Henshaw from LLNL for many useful discussions on moving overset grids. We thank Wieslaw Suszynski of the Coating Process and Visualization Laboratory at the University of Minnesota for assistance in performing the experiment shown in Figure 1. K.L.M. acknowledges the support of this research by the Institute for Mathematics and Its Applications at the University of Minnesota, as well as support from an AWMNSF mentoring travel grant. S.K. acknowledges support from the Department of Energy under award no. DE-FG02-07ER46415. 11362


Langmuir

’ REFERENCES (1) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R.; Witten, T. A. Phys. Rev. E 2000, 62, 756–765. (2) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R.; Witten, T. A. Nature 1997, 389, 827–829. (3) Jing, J.; Reed, J.; Huang, J.; Hu, X.; Clarke, V.; Edington, J.; Housman, D.; Anantharaman, T. S.; Huff, E. J.; Mishra, B.; Porter, B.; Shenker, A.; Wolfson, E.; Hiort, C.; Aston, C.; Schwartz, D. C. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 8046–8051. (4) Derby, B. Annu. Rev. Mater. Res. 2010, 40, 395–414. (5) Brutin, D.; Sobac, B.; Loquet, B.; Sampol, J. J. Fluid Mech. 2011, 667, 85–95. (6) Wong, T.-S.; Chen, T.-H.; Shen, X.; Ho, C.-M. Anal. Chem. 2011, 83, 1871–1873. (7) Deegan, R. D. Phys. Rev. E 2000, 61, 475–485. (8) Fischer, B. J. Langmuir 2002, 18, 60–67. (9) Hu, H.; Larson, R. G. Langmuir 2005, 21, 3963–3971. (10) Craster, R. V.; Matar, O. K.; Sefiane, K. Langmuir 2009, 25, 3601–3609. (11) Routh, A. F.; Zimmerman, W. B. Chem. Eng. Sci. 2004, 59, 2961–2968. (12) Yiantsios, S. G.; Higgins, B. G. Phys. Fluids 2006, 18, 082103. (13) Popov, Y. Phys. Rev. E 2005, 71, 036313. (14) Hu, H.; Larson, R. G. Langmuir 2005, 21, 3972–3980. (15) Hu, H.; Larson, R. G. J. Phys. Chem. B 2006, 110, 7090–7094. (16) Yarin, A. L.; Szczech, J. B.; Megardis, C. M.; Zhang, J.; Gamota, D. R. J. Colloid Interface Sci. 2006, 294, 343–354. (17) Okuzono, T.; Kobayashi, M.; Doi, M. Phys. Rev. E 2009, 80, 021603. (18) Bhardwaj, R.; Xiaohua, F.; Attinger, D. New J. Phys. 2009, 11, 075020. (19) Bhardwaj, R.; Fang, X.; Somasundaran, P.; Attinger, D. Langmuir 2010, 26, 7833–7842. (20) Joshi, A. S.; Sun, Y. J. Displ. Technol. 2010, 6, 579–585. (21) Joshi, A. S.; Sun, Y. Phys. Rev. E 2010, 82, 041401. (22) Hu, H.; Larson, R. G. In Evaporative Self-Assembly of Ordered Complex Structures; Lin, Z., Ed.; World Scientific: Singapore, 2011. (23) Ristenpart, W. D.; Kim, P. G.; Domingues, C.; Wan, J.; Stone, H. A. Phys. Rev. Lett. 2007, 99, 234502. (24) Sultan, E.; Boudaoud, A.; Ben Amar, M. J. Fluid Mech. 2005, 543, 183–202. (25) Cazabat, A.; Guena, G. Soft Matter 2010, 6, 2591–2612. (26) Murisic, N.; Kondic, L. Phys. Rev. E 2008, 78, 065301. (27) Murisic, N.; Kondic, L. J. Fluid Mech. 2011, 679, 219–246. (28) Anderson, D. M.; Davis, S. H. Phys. Fluids 1995, 7, 248–265. (29) Ajaev, V. J. Fluid Mech. 2005, 528, 279–296. (30) Sodtke, C.; Ajaev, V. S.; Stephan, P. J. Fluid Mech. 2008, 610, 343–362. (31) Russel, W. B.; Saville, D. A.; Schowalter, W. R. Colloidal Dispersions; Cambridge University Press: Cambridge, U.K., 1989. (32) Burelbach, J. P.; Bankoff, S. G.; Davis, S. H. J. Fluid Mech. 1988, 195, 463–494. (33) Ajaev, V. S.; Homsy, G. M. J. Colloid Interface Sci. 2001, 240, 259–271. (34) Marek, R.; Straub, J. Int. J. Heat Mass Transfer 2001, 44, 39–53. (35) Oron, A.; Davis, S. H.; Bankoff, S. G. Rev. Mod. Phys. 1997, 69, 931–980. (36) Hu, H.; Larson, R. G. J. Phys. Chem. B 2002, 106, 1334–1344. (37) Gramlich, C. M.; Kalliadasis, S.; Homsy, G. M.; Messer, C. Phys. Fluids 2002, 14, 1841–1850. (38) Chesshire, G.; Henshaw, W. D. J. Comput. Phys. 1990, 90, 1–64. (39) Brown, D. L.; Henshaw, W. D.; Quinlan, D. J. In Scientific Computing in Object-Oriented Parallel Environments; Ishikawa, Y., Oldehoeft, R. R., Reynders, J. V. W., Tholburn, M., Eds.; Lecture Notes in Computer Science; Springer: Berlin, 1997; Vol. 1343, pp 177184. (40) Henshaw, W. D.; Chand, K. K. The Overture Project. http:// www.llnl.gov/casc/Overture.

ARTICLE

(41) Fast, P.; Shelley, M. J. J. Comput. Phys. 2004, 195, 117–142. (42) Meakin, R. L. In Handbook of Grid Generations; Thompson, J. F., Soni, B. K., Weatherill, N. P., Eds.; CRC Press: Boca Raton, FL, 1999; Chapter 11 (43) Henshaw, W. D.; Schwendeman, D. W. J. Comput. Phys. 2006, 216, 744–779. (44) Hundsdorfer, W.; Verwer, J. G. Numerical Solution of TimeDependent Advection-Diffusion-Reaction Equations; Springer-Verlag: Berlin, 2003.

11363


Fast Evaporation of Spreading Droplets of Colloidal Suspensions

Recommend Documents