Article pubs.acs.org/Langmuir
Sagging of Evaporating Droplets of Colloidal Suspensions on Inclined Substrates Leonardo Espín and Satish Kumar* Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, Minnesota 55455, United States ABSTRACT: A droplet of a colloidal suspension placed on an inclined substrate may sag under the action of gravity. Solvent evaporation raises the concentration of the colloidal particles, and the resulting viscosity changes may influence the sag of the droplet. To investigate this phenomenon, we have developed a mathematical model for perfectly wetting droplets based on lubrication theory and the rapid-vertical-diffusion approximation. Precursor films are assumed to be present, the colloidal particles are taken to be hard spheres, and particle and liquid dynamics are coupled through a concentration-dependent viscosity and diffusivity. Evaporation is assumed to be limited by how rapidly solvent molecules can transfer from the liquid to the vapor phase. The resulting one-dimensional system of nonlinear partial differential equations describing the evolution of the droplet height and particle concentration is solved numerically for a range of initial particle concentrations and substrate temperatures. The solutions reveal that the interaction between evaporation and non-Newtonian suspension rheology gives rise to several distinct regimes of droplet shapes and particle concentration distributions. The results provide insight into how evaporation and suspension rheology can be tuned to minimize sagging as well as the well-known coffee-ring effect, an outcome which is important for industrial coating processes.
1. INTRODUCTION Industrial coating processes often have at their core the deposition of a thin film of a colloidal suspension onto a solid substrate. Film deposition is then followed by solvent evaporation, leaving behind a particulate coating.1−4 Whereas deposition and evaporation of colloidal suspensions on horizontal substrates have been extensively studied, the situation where the substrate is inclined has received less attention. Inclined substrates are particularly relevant to spray coating operations, where droplets impact a substrate and then spread and flow under the action of surface tension and gravity. Gravity-driven flow causes droplets to sag such that they are thicker downstream. These thickness variations can be locked in by evaporation, which raises the viscosity of the coating and slows surface-tension-driven leveling. Such thickness variations are generally undesirable, and understanding how they develop in colloidal suspensions is a particularly complex problem since particle concentrations, and thus liquid viscosity, can vary in space and time. Thickness variations can occur in films and droplets of colloidal suspensions even in the absence of evaporation due to rheological effects. A prominent example of this is when regions of the droplet having lower particle concentration flow faster than regions having higher concentration.5,6 For droplets of colloidal suspensions on horizontal substrates, solvent evaporation can produce flows that cause particles to accumulate at the liquid−vapor interface (“skin formation”)7 or in the contact-line region.8−10 The latter phenomenon is the well-known “coffee-ring” effect. © 2014 American Chemical Society
The purpose of the present study is to develop a mathematical model to study the sagging of evaporating droplets of colloidal suspensions on inclined substrates. The model is an extension of an earlier one developed by us to study the forced spreading of films and droplets of colloidal suspensions in the absence of evaporation.6 Lubrication theory and the rapid-vertical-diffusion approximation are used to derive a coupled system of one-dimensional partial differential equations describing the evolution of the droplet height and particle concentration. The colloidal particles are taken to be hard spheres, and particle and liquid dynamics are coupled through a concentration-dependent viscosity and diffusivity. We focus on the case of perfectly wetting droplets where precursor films are present at the droplet edge. Evaporation is assumed to be limited by how rapidly solvent molecules can transfer from the liquid to the vapor phase, allowing us to use a so-called “one-sided” model.11−18 (The implications of this assumption are discussed in section 2.2.) The evolution equations are solved numerically, and the results allow us to develop a map relating the suspension volume fraction and evaporation rate to coating quality (characterized by the amount of sag and particle concentration distribution). The governing equations and scalings for our model are described in section 2, followed by a discussion in section 3 of the numerical method we use. Various aspects of the results are Received: August 13, 2014 Revised: September 16, 2014 Published: September 17, 2014 11966
dx.doi.org/10.1021/la503229z | Langmuir 2014, 30, 11966−11974
Langmuir
Article −2 ⎛ Φ(x′, z′, t ′) ⎞ ⎟ μ′ = ⎜1 − ⎝ ⎠ 0.64
discussed in sections 4, 5, and 6. Finally, concluding remarks are offered in section 7.
2. GOVERNING EQUATIONS AND SCALINGS 2.1. Hydrodynamics. We consider the two-dimensional flow of a droplet spreading on a uniformly heated substrate at temperature T0 that forms an angle 0 < α ≤ π/2 with respect to the horizontal (Figure 1). The droplet consists of a Newtonian
which has been used in several prior studies. Equation 5 indicates that the suspension viscosity rapidly increases with the particle concentration, and it diverges when the volume fraction approaches the maximum value for hard spheres, Φ = 0.64. In the present work we focus on the case of perfectly wetting droplets and assume that thin precursor films are present at the droplet edge. In particle-free systems, precursor films have been observed and have thicknesses that range from 10 nm to 1 mm.28 The use of a precursor film is also relevant for situations where a droplet flows over a prewetted plane.29 An additional benefit of using precursor films is that they naturally alleviate the stress singularity that arises at moving contact lines. Finally, we note that by incorporating a suitable disjoining pressure, partially wetting systems can be modeled.30,31 2.2. Evaporation. 2.2.1. General Remarks. Before proceeding, it is important to make some general remarks about the modeling of evaporation. Solvent evaporation involves transfer of solvent molecules from the liquid droplet into the vapor phase and transport of those molecules away from the liquid−vapor interface. In general, modeling this process involves solving equations that describe solvent transport in both the liquid and vapor phases. In addition, the location of the liquid−vapor interface is typically unknown, leading to a highly nonlinear free-boundary problem that is a formidable computational challenge. To simplify matters, prior work has often focused on two limiting cases. In one case, evaporation is assumed to be limited by how quickly solvent molecules can diffuse in the vapor phase. The droplet shape is assumed to be a fixed spherical cap, and a diffusion equation is solved in the vapor phase to determine the solvent concentration there. The resulting evaporation flux at the droplet−vapor interface can then be used as a boundary condition to determine the fluid flow and solute transport inside the droplet. This approach has been widely used to study the evaporation into stagnant air of liquid droplets resting on horizontal substrates.10,32,33 However, the assumption of a spherical-cap interface may no longer hold if the droplet is on an inclined surface. In the other limiting case, evaporation is assumed to be limited by how quickly solvent molecules can transfer from the liquid phase to the vapor phase. The advantage of this approach is that only the solvent concentration in the liquid phase needs to be accounted for, leading to a so-called “one-sided” model of evaporation.11−18 When coupled with the lubrication approximation, free-boundary problems can readily be solved. However, these one-sided models are typically based on the assumption that the vapor phase is pure (i.e., solvent molecules only, no air),11,34,35 which may not be the case in many industrial coating processes. It should also be noted that, for single-component systems, evaporation corresponds to boiling, which is usually undesirable in coating-related applications because the resulting vapor bubbles may produce defects. More detailed discussion of the assumptions underlying the one-sided evaporation model can be found in refs 16 and 36. Although the assumptions of the one-sided evaporation model may appear to be rather restrictive, recent work has shown that it can be adapted to describe the evaporation of water droplets into an ambient air atmosphere in the absence of contact-line pinning. Murisic and Kondic16 combined the
Figure 1. Schematic of a droplet spreading on a heated inclined substrate.
solvent in which hard, spherical, colloidal particles are dispersed. The suspension has constant density ρ and a particle-concentration-dependent viscosity μ(Φ), where Φ(x, z, t) is the particle volume fraction, x is the coordinate parallel to the substrate, z is the coordinate normal to the substrate, and t is time. The droplet is assumed to be long and thin such that its characteristic thickness h0 and length l satisfy ε = h0/l ≪ 1, allowing the lubrication approximation to be invoked. The characteristic draining speed of the flow is u* = [h02ρg sin(α)]/ μ0, where g is the gravitational acceleration and μ0 is the solvent viscosity.6,19−21 The time scale is defined as t* = l/u*, and the pressure scale, p* = lρg sin(α), can be obtained from the xcomponent of the momentum balance equations. With these scales we define the following dimensionless variables (denoted with primes): x = lx′
z = h 0 z′
p = p*p′
u = u*u′
l t= t′ u*
v = εu*v′
μ = μ0 μ′
At leading order in ε, the equations governing momentum and mass conservation reduce to 0 = −p′x ′ + (μ′u′z ′)z ′ + 1
(1)
0 = p′ z ′
(2)
u′x ′ + v′z ′ = 0
(3)
and these are subject to no-slip and no-penetration conditions at the substrate: u′(x′, z′ = 0, t ′) = v′(x′, z′ = 0, t ′) = 0
(5) 6,7,24−27
(4)
The boundary conditions at the liquid−vapor interface will be discussed in the following section. We assume that the suspension viscosity is described by the well-known Krieger−Dougherty equation22,23 11967
dx.doi.org/10.1021/la503229z | Langmuir 2014, 30, 11966−11974
Langmuir
Article
where β = κTSu0 /[3h0 3ρ2 g sin(α)] and Θ = (T0 − TS)/εTS is the dimensionless substrate temperature or superheat. With the defined scalings for T and J we can render eqs 6 and 7 dimensionless:
lubrication approximation and a one-sided evaporation model to describe the spreading of a water droplet on a silicon wafer. To determine the saturation temperature, which is a key parameter in the one-sided model, they assumed that the saturation pressure was equal to the partial pressure of water and then used the Clasius−Clapeyron equation (which relates vapor pressure and temperature for a single-component system) to calculate the saturation temperature. Thus, the saturation temperature was determined by assuming that the entire system is effectively a single-component system held at the partial pressure of water. This assumption might be expected to be more accurate the higher the partial pressure is relative to the system pressure (which is usually atmospheric pressure). Predictions from the model of Murisic and Kondic16 of the water droplet size as a function of time agreed well with their experimental observations. However, when applied to the case of a 2-propanol droplet, the agreement between theory and experiment was much less satisfactory. Interestingly, the ratio of the partial pressure (taken to be the saturation pressure) to atmospheric pressure is ∼10−2 for water, whereas it is ∼10−3 for 2-propanol. This is consistent with the observation made above that better agreement might be expected for larger values of this ratio. Alternatively, it has been pointed out that if the water surface is contaminated, the contaminated layer can inhibit evaporation, making the rate-limiting step the transfer of solvent molecules to the vapor phase.35 Given the simplicity and versatility of the one-sided model, we adopt it in the present work. Its use here is also appropriate because we neglect contact-line pinning (which was also absent in the work of Murisic and Kondic16) and consider heated substrates.35 Below, we present the key equations of the onesided model. 2.2.2. Model Details. To decouple the dynamics of the vapor phase from those of the liquid, the one-sided model assumes that the density, viscosity, and thermal conductivity of the vapor are much smaller than those of the liquid. Under these conditions, it can be shown that all of the heat conducted to the liquid−vapor interface, z = h(x, t), is used to vaporize the liquid,11 leading to the following interfacial energy balance: 3 = −κTz
at
z = h(x , t )
J ′ = − βT ′ z ′
at
z = h(x , t )
(11)
σ = σ0 − γ(T − TS)
(12)
where σ0 is the surface tension of the solvent at the saturation temperature TS. The lubrication form of the tangential-stress balance equation at the interface is then μuz = σx + σzhx, which has the following nondimensional form: μ′u′z ′ = −Ma(T ′x ′ + T ′z ′h′x ′)
at
z = h(x , t )
(13)
where the Marangoni number Ma = γTS/[l ρg sin(α)]. In the contact-line region, van der Waals forces suppress evaporation, resulting in the formation of a thin precursor film.13 van der Waals interactions are modeled with a disjoining pressure term Π = (/h3, where ( is the Hamaker constant.13 The disjoining pressure is incorporated through the interfacial normal-stress balance, which has the lubrication form −(p − pv) = σ0hxx + Π, where pv is the pressure in the vapor phase (assumed to be pure). The scaling for the disjoining pressure term, Π = lμ0u*Π′/ h02, is obtained by balancing with the pressure term of the normal-stress balance. The resulting dimensionless disjoining pressure Π′ = (′/h′3, where the dimensionless Hamaker constant satisfies (′ = (/[h0 3lρg sin(α)]. Thus, the dimensionless normal-stress balance is 2
− (p′ − p′ v ) =
1 h′x ′ x ′ + Π′ Ca ̂
(14)
where the parameter Ca ̂ = [l 3ρg sin(α)]/σ0h0 is a rescaled capillary number. According to this definition, the small-slope ̂ ≪ ε−1. restriction of lubrication theory6 requires that Ca 1/3 The evaporation model is closed by specifying a constitutive equation for the mass flux J. As in prior works, we use a a modified Hertz−Knudsen relation: 2πRTS ρv
(7)
J=
3 1 (p − pv ) + (T − TS) ρ TS
at
z = h(x , t ) (15)
Under the lubrication approximation, the energy conservation equation in the liquid is Tzz = 0
at
(10)
where we have defined h′ = h/h0. The surface tension of the solvent is taken to be a linear function of the temperature:7,11,13,18
(6)
z = h(x , t )
z′ = h′(x′, t ′)
−h′t ′ − u′h′x ′ + v′ = J ′
where J = J(x, t) is the mass flux at the interface, 3 is the latent heat of vaporization per unit mass, and κ is the thermal conductivity of the liquid. We obtain the scaling for the evaporative flux, J = ερu*J′, by balancing the terms in the kinematic boundary condition at the liquid−vapor interface: ρ( −ht − uhx + v) = J
at
where R is the gas constant per unit mass and ρv is the vapor density. Equation 15 reflects the assumption that departures from thermodynamic equilibrium (p ≠ pv or T ≠ TS) drive evaporation in the one-sided model.11,16 Equation 15 has the following dimensionless form:
(8)
As a reference temperature, we choose the saturation temperature, TS, at the given vapor pressure.7,11,13,18 We then define a dimensionless temperature T′ via the equation T = TS + εTST′. Integration of eq 8 subject to eq 6 and the condition T(x, z = 0, t) = T0 yields the dimensionless droplet temperature
KJ ′ = δ(p′ − p′v ) + T ′
where
at
z = h(x , t )
K = [h0 2ρ2 g sin(α)(2πRTS)1/2 ]/ρv 3μ0
(16)
and
2
T′ = Θ −
J′ z′ β
δ = [l g sin(α)]/h0 3 . 2.3. Particle Concentration. Transport of the colloidal particles is governed by a convection-diffusion equation:
(9) 11968
dx.doi.org/10.1021/la503229z | Langmuir 2014, 30, 11966−11974
Langmuir
Article
∂c + u ⃗ ·∇c = ∇·(D∇c) ∂t
⎡ h2 ⎛ ⎤ ⎞ h Ma h (Jh)x ⎥ϕx ϕt + μ−1(ϕ)⎢ ⎜1 + xxx + Πx⎟ + ⎠ β 2 ⎣3⎝ ⎦ Ca ̂ Jϕ 1 (D(ϕ)hϕx)x + = h h(Pe) (26)
(17)
where c = c(x, z, t) is the particle concentration. The diffusivity of the colloidal particles is described by a generalized Stokes− Einstein equation:6,24,37 6.55
D = D0(1 − Φ)
d ⎛ 1.85Φ ⎞ ⎜ ⎟ dΦ ⎝ 0.64 − Φ ⎠
where the dimensionless evaporative flux is given by
(18)
J=
where D0 = kBTS/6πμ0a, with kB being Boltzmann’s constant and a the radius of the colloidal particles. Equation 17 is complemented by no-flux conditions at the substrate and liquid−vapor interfaces:
cz|z = 0 = 0 and at
z = h(x , t )
(20)
Table 1. Order-of-Magnitude Estimates of Physical Parameters
Equation 17 is coupled to the momentum conservation equations through the particle volume fraction, defined through c(x, z, t) = ρpΦ(x, z, t), where ρp is the density of the particles. The dimensionless lubrication form of the convection-diffusion equation is ε 2Pe(Φt ′ + u′Φx ′ + v′Φz ′) = ε 2(D′Φx ′)x ′ + (D′Φz ′)z ′
(21)
where Pe = u*l/D0 is the Peclet number. To further simplify eq 21, we consider the rapid-verticaldiffusion limit,6,18,25,27 in which diffusion in the z-direction is assumed to be dominant so that the strongest concentration gradients are in the x-direction. We thus decompose the particle volume fraction as its z-average plus a small correction: Φ(x′, z′, t ′) = ϕ(x′, t ′) + ε 2Peϕ1(x′, z′, t ′)
(22)
where ε2Pe ≪ 1 is assumed. Using eq 22 with eqs 19−21 allows elimination of zdependence in the particle transport via vertical averaging in that direction. At leading order in ε2Pe the equation for the zaveraged particle volume fraction ϕ is 1 h(ϕt + u ̅ ϕx) = (Dhϕx)x + Jϕ Pe
1 h
∫0
(23)
h
f (x , z , t ) d z
(24)
and the prime notation has been dropped. 2.4. Evolution Equations and Parameter Values. It follows from eq 22 that, at leading order in ε2Pe, the equations of fluid motion couple with the particle concentration through a z-independent viscosity term, μ(ϕ). This allows us to explicitly integrate eqs 1−4 with respect to z and obtain a pair of coupled evolution equations for the droplet height and zaveraged particle volume fraction:
⎤⎫ Ma h2 (Jh)x ⎥⎬ − J β 2 ⎦⎭
typical value 10−5 10−2 103 10−3 0.0728 10−8 297 10−2 461.15 0.607 4.18 × 103 106 10−4 10−20−10−17
3. SOLUTION METHOD Equations 25 and 26 are solved numerically on the domain 0 ≤ x ≤ L subject to boundary conditions
h(0, t ) = b
⎡ h3 ⎞ ∂ ⎧ h3 ⎛ h ht = − ⎨μ−1(ϕ)⎢ + ⎜ xxx + Πx⎟ ̂ ⎠ ∂x ⎩ 3 ⎝ Ca ⎣3 +
parameter droplet thickness, h0 (m) droplet length, l (m) density, ρ (kg/m3) viscosity, μ0 [kg/(m s)] surface tension, σ (kg/s2) particle diameter, a (m) saturation temperature, TS (K) vapor density, ρv (kg/m3) gas constant, R [J/(kg K)] thermal conductivity, κ [W/(m K)] specific heat, cp [J/(kg K)] latent heat, 3 (J/kg) γ [N/(m K)] Hamaker constant, ( (J)
have a diffusivity D0 ≈ 10−11 m2/s at the saturation temperature given in Table 1. Using this and other values in Table 1, we obtain the estimates Ca ̂ ≈ 104, Pe ≈ 105, Ma ≈ 0.05, β = O(1), δ ≈ 10−5, ( ≈ 10−4 , and K ≈ 10−3. In addition, b = 0.005 and corresponds to a value of 50 nm when the superheat Θ = O(10−2) (b = 0.001 if Θ = O(1)). It follows that, for the values of the parameters shown in Table 1, the small-slope restriction ̂ ≪ ε−1 and the fast-diffusion restriction ε2Pe ≪ 1 are Ca 1/3 satisfied. We note that different choices of the saturation temperature could be made, but we do not expect this to qualitatively affect the results we present later.
where z-averaged quantities are defined as f ̅ (x , t ) =
(27)
We note that when the droplet height becomes thin enough, reaching the value h ≡ b = ((δ /Θ)1/3, the evaporative flux (eq 27) reduces to zero, allowing for the formation of a thin precursor film.13 Table 1 provides order-of-magnitude estimates of various physical parameters.7,11,25 Particles with diameter ∼10−8 m
(19)
−D(n·̂ ∇)c = 0
⎛ hxx ⎞⎤ β ⎡ + Π⎟⎥ ⎢Θ − δ ⎜ ⎝ Ca ̂ ⎠⎦ βK + h ⎣
h(L , t ) = b
hx(0, t ) = 0 hx(L , t ) = 0
(28)
ϕx(0, t ) = 0
(29)
ϕx(L , t ) = 0
The initial condition is given by a fourth-order polynomial which satisfies the conditions in eq 28 imposed when L = 1 and the extra condition of having unit area between x = 0 and x = 1. As described in section 2, the quantity b = ((δ /Θ)1/3 represents the thickness of the precursor film. The initial particle concentration is ϕ(x, 0) = ϕI for 0 ≤ x ≤ 1 and ϕ(x, 0) = 0 for x > 1.
(25) 11969
dx.doi.org/10.1021/la503229z | Langmuir 2014, 30, 11966−11974
Langmuir
Article
The numerical solutions are obtained with a fully implicit finitedifference scheme. All spatial derivatives are approximated with second-order centered differences, and time integration is performed with the DDASPK iterative solver.38 To validate our numerical solution, we compare it against the self-similar approximation h(x , t ) =
K differ from the estimates in section 2, they were chosen for simplicity and in such a way that the products βδ, βK, and Ma/β are consistent with those obtained from the estimated values. (Note that the evaporation parameter β only enters the governing eqs 25−27 through these products; in particular, Ma/β = 0.0068 with the values estimated in section 2.) Computations are carried out for different values of the dimensionless substrate temperature Θ and initial particle concentration ϕI. If the substrate inclination angle α is changed, the values of the relevant dimensionless parameters (Ca ,̂ Pe, β, δ, K, and ( ) are adjusted according to the definitions given in section 2. As shown in ref 6, the width of the precursor film b affects the speed and the long-term shape evolution of the interface of a spreading droplet of a colloidal suspension. To facilitate comparison of solutions having different parameter values, we fix the width of the precursor film to the value b = 0.01 and adjust the dimensionless Hamaker constant according to the value of the substrate temperature Θ. Finally, as noted above, we set Ma = 0 for all computations, with the exception of those of section 6, where we briefly investigate how Marangoni effects influence the results.
(x /t )1/2 1−
ϕI 0.64
(30)
Equation 30 can be obtained for a droplet with uniform particle concentration in the whole computational domain 0 ≤ x ≤ L when evaporation is absent, Θ = δ = Ma = 0. Equation 30 is the self-similar solution of Huppert,39 with an extra term that accounts for the increase in viscosity due to the presence of particles. Equation 30 is obtained by ignoring the effects of surface tension and, as a consequence, is valid only away from the contact-line region.39 As shown in Figure 2, the excellent agreement between the selfsimilar approximation (eq 30) and the numerical solution indicates
4. QUANTIFYING SAG In Figure 3a, we show how changes in the initial particle concentration modify the liquid drainage, or sag. To do this, we fix the substrate temperature at Θ = 0.5 and let the droplet evolve until a final state is reached for a given initial volume fraction ϕI. If enough liquid has been lost through evaporation, the particle volume fraction will reach its maximum possible value, ϕ(x, tM) = 0.64, for some value of x at time tM = tM(ϕI, Θ), and the computation has to be stopped. In practice, the computations are stopped when the volume fraction is equal to or larger than 98.44% of ϕ(x, tM). We note that, despite the appearance of the droplets shown in Figure 3, the actual aspect ratio is h0/l = 10−3 according to the estimated values shown in Table 1. In Figure 3b, we show the concentration (volume fraction) profiles for ϕI = 0.05 and ϕI = 0.5. In the former case, the particles accumulate at the rear of the droplet. In the latter case, concentration spikes appear at both the front and back of the droplet. These concentration profiles will be discussed in more detail in section 5. As shown in Figure 3a, for lower initial particle concentrations the droplets spread rapidly due to the relatively low viscosity. The leading front can reach distances several times the initial diameter of the droplet. Also, for these lowconcentration droplets, the volume reduction due to evapo-
Figure 2. Numerical solution and self-similar approximation (eq 30) at time t = 1084 for the case Ca ̂ = 1000, Pe = 104, b = 0.01, and ϕI = 0.1 when no evaporation is present, Θ = δ = Ma = 0. that the numerically computed droplet is thinning and spreading at the correct speed. Since eq 30 is obtained by ignoring the effects of surface tension, it is valid only away from the contact-line region as is evident from Figure 2. To facilitate the presentation and interpretation of the results, it is useful to consider a substrate inclination angle of α = 90° and the following values of the dimensionless parameters: Ca ̂ = 103, Pe = 104, β = 10−3, δ = 1, Ma = 0, and K = 1. Although the values for β, δ, Ma, and
Figure 3. (a) Droplet profiles for different values of the initial volume fraction ϕI and fixed substrate temperature Θ = 0.5. Profiles are shown at time t = tM(ϕI, Θ) with a substrate inclination angle of α = 10°. (b) Concentration profiles for ϕI = 0.05 (dashed line) and ϕI = 0.5 (solid line). 11970
dx.doi.org/10.1021/la503229z | Langmuir 2014, 30, 11966−11974
Langmuir
Article
concentration plane. The reference value sc has been chosen to correspond to the amount of sag displayed by the ϕI = 0.3 droplet, shown with a dotted line in Figure 3. Thus, values of S larger than sc correspond to the regions below the curves of Figure 4 (for a given substrate inclination angle), and values smaller than sc are in the region above. For reference, the values of sag for the droplets shown in Figure 3 are S(0.1, 0.5) = 8.115 and S(0.05, 0.5) = 16.84 (which are above the continuous line) and S(0.5, 0.5) = 0.054 (which is below the continuous line). Thus, the curves in Figure 4 indicate that if fluid sagging is to be avoided, the substrate temperature increase has to be considerably larger for the cases of droplets with low initial particle concentration, whereas only moderate increases are required to halt sagging for more concentrated suspensions. Also, as the inclination angle of the substrate, and as a consequence the draining speed, is increased, the temperatures required to avoid sagging must increase as well. However, the increase is not uniform; much larger increases are required for droplets with lower initial particle concentrations.
rative liquid loss near x = 0 can be so severe that evaporative displacement of the back contact line occurs, together with the formation of “dry patches”, where the droplet has thinned until it reaches the width of the precursor film. In contrast, the more concentrated suspension with ϕI = 0.3 only slightly deforms before the maximum volume fraction is reached. Thus, the amount of sag is the result of a competition between gravitydriven liquid drainage and evaporation, which counters drainage by increasing the viscosity of the suspension through loss of liquid. To quantify sag, we define the following quantity: S(ϕI , Θ) =
1 L
∫0.5 h(x , 0) dx
tM
L
∫0 ∫0.5 h(x , t ) dx dt
(31)
Equation 31 is a normalized measure of droplet volume past x = 0.5 over the droplet lifetime (0 ≤ t ≤ tM). This definition is motivated by analogous ones that have appeared in prior work.40−42 In Figure 4 we show curves corresponding to S(ϕI = 0.3, Θ = 0.5) = 2.114 s c in the superheat−initial particle
5. PARTICLE DISTRIBUTIONS We now discuss in more detail the distribution of particles inside the evaporating droplets and how this depends on the problem parameters. Figure 5 shows the droplet height and particle concentration in the presence and absence of evaporation. Here, the initial concentration and superheat have relatively small values (ϕI = 0.1, Θ = 0.1). Gravity-driven drainage causes the formation of a bulky region at the front of the droplet with a much thinner region at the back. This, together with the fact that the evaporative flux (eq 27) is largest near the contact lines (which follows from the observation that the flux decreases with height and increases with interfacial curvature7,18), implies that the largest loss of liquid occurs near x = 0. As a consequence, the largest concentrations occur in the tail region of the droplet. For a discussion of the behavior of the velocity field in the case of spreading on horizontal substrates, see ref 7. This large loss of liquid also leads to a displacement of the contact line at the rear of the droplet, seen in Figure 5a, or in Figure 3 for the curves ϕI = 0.05 and ϕI = 0.1. (Note that evaporation is suppressed once the liquid height reaches the size of the precursor film b.) As can be observed in Figure 5b, this displacement is accompanied by a gradual increase in viscosity due to the increase in concentration, which eventually
Figure 4. For a given inclination angle the continuous or dotted line separates the regions of no sag (for superheats above the corresponding curve) and sag (for superheats below the corresponding curve). The curves indicate the values of substrate superheat Θ and initial particle concentration ϕI where the observed sag (eq 31) satisfies S = sc.
Figure 5. (a, left) Droplet height h(x, t) and (b, right) particle concentration ϕ(x, t) at time tM = 20.62 for the cases when the substrate superheat is Θ = 0.1 (red dashed curves) and Θ = δ = 0 (black curves). The solutions correspond to the case ϕI = 0.1 and α = 10°. 11971
dx.doi.org/10.1021/la503229z | Langmuir 2014, 30, 11966−11974
Langmuir
Article
7b). These large concentrations raise the droplet viscosity and cause this part of the droplet to flow more slowly than the rest of the droplet. The concentration spikes at both the front and rear contact lines seen in Figure 7b when Θ = 0.2 are similar to the wellknown coffee-ring effect that has been observed for droplets evaporating on horizontal substrates.7−9,43 Figure 3b shows another example of a coffee-ring-like pattern in a different region of parameter space (ϕI = 0.5, Θ = 0.5). Evaporation raises the particle concentration near both contact lines, halting droplet spreading. In the rest of the droplet, the concentration is nearly uniform. These concentration profiles form when evaporation is faster than liquid drainage. Droplet drying proceeds inward from the edges since evaporation is fastest at the contact lines. The above results can be summarized by the map shown in Figure 8, which identifies the different types of concentration
halts the displacement and causes the back contact line to become pinned when the volume fraction approaches 0.64. For larger initial concentrations and the same superheat, the large concentration spike at the rear of the droplet is considerably diminished, as shown in Figure 6. In this regime, the higher
Figure 6. Droplet height h(x, t) (black curve) and particle concentration ϕ(x, t) (red curve) for the case Θ = 0.1, ϕI = 0.5, and α = 10°. The dashed line shows the initial droplet height h(x, t = 0).
initial concentration causes the droplet to spread more slowly and the maximum packing fraction to be reached more quickly. As a result, the concentration distribution is more uniform and the droplet sags considerably less. As alluded to above, the concentration increases in the rear of the droplet first and then progressively through the rest of the droplet. This behavior is similar to that of “lateral” consolidation fronts that have been observed experimentally in films of colloidal suspensions drying on horizontal substrates.4 We now consider what happens when the superheat is increased to a slightly larger value. Figure 7 shows calculations for ϕI = 0.55 and two different values of Θ. The Θ = 0.1 case is similar to that shown in Figure 6. However, when Θ is increased to 0.2, the droplet sags less (Figure 7a) and spikes in the concentration profiles are observed near both the front and back contact lines (Figure 7b). The kinks in the height profile near the rear of the droplet for Θ = 0.1 (Figure 7a) are a consequence of the large particle concentrations there (Figure
Figure 8. Map showing different types of concentration distributions in the ϕI−Θ plane. The dotted line is redrawn from Figure 4 for the case α = 10°. Below the solid curve, a concentration spike is observed in the back only. Below the dashed curve, the concentration distribution is relatively uniform.
distributions that are observed in the ϕI−Θ plane. For reference, we include the curve from Figure 4 for α = 10° (shown here with a dotted line); below this line the droplet sags, whereas above this line no sag occurs according to the definition given earlier.
Figure 7. (a) Droplet height h(x, t = tM) and (b) particle concentration ϕ(x, t = tM) for the case α = 10° and ϕI = 0.55 when the superheat Θ = 0.1 (continuous curve) and Θ = 0.2 (dashed curve). 11972
dx.doi.org/10.1021/la503229z | Langmuir 2014, 30, 11966−11974
Langmuir
Article
7. CONCLUSIONS Motivated by the need to improve fundamental understanding of the physical mechanisms that influence the behavior of thinfilm coatings of colloidal suspensions, we have investigated how the interaction between evaporation and suspension rheology modifies the behavior of droplets sliding down an inclined substrate. Our results make clear that there is a trade-off between sagging and the coffee-ring effect. Faster evaporation rates lead to less sagging but also less uniform particle concentration distributions. Larger initial particle concentrations and smaller evaporation rates tend to minimize sagging and yield a relatively uniform concentration distribution. Whereas much prior work on evaporating droplets of colloidal suspensions has focused on horizontal substrates, the novelty of the present contribution is the consideration of inclined substrates, a situation of considerable industrial importance (e.g., spray coating of discrete objects). Although our model makes several simplifying assumptions, it provides a platform for exploring a number of additional effects. These include partial wetting,12,14 particle adsorption and desorption at the liquid−vapor interface and in the contact-line region,44 different evaporation modes (e.g., diffusion-limited evaporation), more complex rheology,45,46 the effects of structural disjoining pressure,18,26,47−50 particle concentration gradients in the vertical direction,7 and three-dimensional dynamics. Our study (and extensions of it) can also be used to guide the design and interpretation of experiments on evaporating droplets of colloidal suspensions, which display many fascinating yet unexplained phenomena.10
Also shown in this plot are lines that indicate where (i) a spike in the concentration is observed only at the back of the droplet (cf. Figure 5) and (ii) the concentration distribution is relatively uniform (cf. Figure 6). These lines were drawn by setting threshold values for a concentration spike such that the lines nearly fall on the sag−no sag curve (dotted line). In this way, sagging can be correlated with different types of concentration distributions. Figure 8 makes clear that there is a trade-off between sagging and the coffee-ring effect. Larger values of the superheat lead to less sagging, but also tend to lead to coffee-ring-like concentration distributions (cf. Figure 7b). More uniform concentration distributions can be achieved by lowering the superheat, and sagging can be reduced by raising the initial particle concentration.
6. MARANGONI EFFECTS We now briefly discuss the influence of Marangoni flows on droplet behavior. Variations in surface tension due to temperature gradients drive interfacial flows from warmer regions (lower surface tension) toward colder regions (higher surface tension). As noted in prior studies involving horizontal substrates,7,18 these flows tend to elongate the droplet in the vertical direction, slow spreading, and reduce the evaporative flux. As a consequence, the time required to reach the maximum packing fraction increases. Our results for inclined substrates are consistent with these observations (Figure 9). Increasing the value of Ma from 0 to
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported through the Industrial Partnership for Research in Interfacial and Materials Engineering at the University of Minnesota. S.K. acknowledges support from the Department of Energy under Award No. DE-FG0207ER46415. We thank Professor Lorraine Francis for helpful discussions.
Figure 9. Droplet height plots for the case ϕI = 0.05, Θ = 10, and α = 90° for several values of the Marangoni number. The black curves correspond to the same computation at two different times: the black dashed curve shows the computation at the same time as the red, continuous curve, and the dashed−dotted curve shows the computation at time tM(ϕI = 0.05, Θ = 10, Ma = 0.005) = 2.065. Note that tM(ϕI = 0.05, Θ = 10, Ma = 0) = 1.7313.
■
REFERENCES
(1) Cardinal, C. M.; Jung, Y. D.; Ahn, K. H.; Francis, L. F. Drying regime maps forparticulate coatings. AIChE J. 2010, 56, 2769−2780. (2) Kim, S.-H.; Lee, S. Y.; Yi, G.-R.; Pine, D. J.; Yang, S.-M. Microwave-assisted self-organization of colloidal particles in confining aqueous droplets. J. Am. Chem. Soc. 2006, 128, 10897−10904. (3) Aliseda, A.; Berchielli, A.; Doshi, P.; Lasheras, J. Chemical Engineering in the Pharmaceutical Industry: R&D to Manufacturing; Wiley: Hoboken, NJ, 2010; pp 781−799. (4) Routh, A. F. Drying of thin colloidal films. Rep. Prog. Phys. 2013, 76, No. 046603. (5) Buchanan, M.; Molenaar, D.; Villiers, S. D.; Evans, R. Pattern formation in draining thin film suspensions. Langmuir 2007, 23, 3732−3736. (6) Espı ́n, L.; Kumar, S. Forced spreading of films and droplets of colloidal suspensions. J. Fluid Mech. 2014, 742, 495−519. (7) Maki, K. L.; Kumar, S. Fast evaporation of spreading droplets of colloidal suspensions. Langmuir 2011, 27, 11347−11363.
0.005 tends to increase the height of the droplet during the earlier stages of evaporation. However, the reduced evaporative flux means that gravity has more time to act on the droplet, and this leads to an increase in sag. We note that the value of Ma/β used is these calculations is 5, which is much larger than the value estimated in section 3. Thus, we would not expect Marangoni flows to have a significant influence on the results presented in the prior section in the parameter range we have examined. 11973
dx.doi.org/10.1021/la503229z | Langmuir 2014, 30, 11966−11974
Langmuir
Article
(8) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G. Capillary flow as the cause of ring stains from dried liquid drops. Nature 1997, 389, 827−829. (9) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R.; Witten, T. A. Contact line deposits in an evaporating drop. Phys. Rev. E 2000, 62, 756−765. (10) Larson, R. Transport and deposition patterns in drying sessile droplets. AIChE J. 2014, 60, 1538−1571. (11) Burelbach, J. P.; Bankoff, S. G.; Davis, S. H. Nonlinear stability of evaporating/condensing liquid films. J. Fluid Mech. 1988, 195, 463− 494. (12) Ajaev, V. S. Viscous flow of a volatile liquid on an inclined heated surface. J. Colloid Interface Sci. 2004, 280, 165−173. (13) Ajaev, V. S. Spreading of thin volatile liquid droplets on uniformly heated surfaces. J. Fluid Mech. 2005, 528, 279−296. (14) Ajaev, V. S.; Gambaryan-Roisman, T.; Stephan, P. Static and dynamic contact angles of evaporating liquids on heated surfaces. J. Colloid Interface Sci. 2010, 342, 550−558. (15) Sodtke, C.; Ajaev, V.; Stephan, P. Dynamics of volatile liquid droplets on heated surfaces: theory versus experiment. J. Fluid Mech. 2008, 610, 343−362. (16) Murisic, N.; Kondic, L. On evaporation of sessile drops with moving contact lines. J. Fluid Mech. 2011, 679, 219−246. (17) Oron, A.; Davis, S. H.; Bankoff, S. G. Long-scale evolution of thin liquid films. Rev. Mod. Phys. 1997, 69, 931−980. (18) Craster, R. V.; Matar, O. K.; Sefiane, K. Pinning, retraction, and terracing of evaporating droplets containing nanoparticles. Langmuir 2009, 25, 3601−3609. (19) Goodwin, R.; Homsy, G. M. Viscous flow down a slope in the vicinity of a contact line. Phys. Fluids 1991, 515−528. (20) Moriarty, J. A.; Schwartz, L. W.; Tuck, E. O. Unsteady spreading of thin liquid films with small surface tension. Phys. Fluids 1991, 3, 733−742. (21) Eres, M. H.; Schwartz, L. W.; Roy, R. V. Fingering phenomena for driven coating films. Phys. Fluids 2000, 12, 1278−1295. (22) Brady, J. F. The rheological behavior of concentrated colloidal dispersions. J. Chem. Phys. 1993, 99, 567−581. (23) Krieger, I. M.; Dougherty, T. J. A mechanism for nonNewtonian flow in suspensions of rigid spheres. J. Rheol. 1959, 3, 137− 152. (24) Yiantsios, S. G.; Higgins, B. G. Marangoni flows during drying of colloidal films. Phys. Fluids 2006, 18, No. 082103. (25) Warner, M. R. E.; Craster, R. V.; Matar, O. K. Surface patterning via evaporation of ultrathin films containing nanoparticles. J. Colloid Interface Sci. 2003, 267, 92−110. (26) Matar, O.; Craster, R.; Sefiane, K. Dynamic spreading of droplets containing nanoparticles. Phys. Rev. E 2007, 76, No. 056315. (27) Tsai, B.; Carvalho, M. S.; Kumar, S. Leveling of thin films of colloidal suspensions. J. Colloid Interface Sci. 2010, 343, 306−313. (28) Popescu, M. N.; Oshanin, G.; Dietrich, S.; Cazabat, A.-M. Precursor films in wetting phenomena. J. Phys.: Condens. Matter 2012, 24, No. 243102. (29) Ye, Y.; Chang, H.-C. A spectral theory for fingering on a prewetted plane. Phys. Fluids 1999, 11, 2494−2515. (30) Schwartz, L. W.; Eley, R. R. Simulation of droplet motion on low-energy and heterogeneous surfaces. J. Colloid Interface Sci. 1998, 202, 173−188. (31) Schwartz, L. W. Hysteretic effects in droplet motions on heterogeneous substrates: direct numerical simulation. Langmuir 1998, 14, 3440−3453. (32) Erbil, H. Evaporation of pure liquid sessile and spherical suspended drops: a review. Adv. Colloid Interface Sci. 2012, 170, 67− 86. (33) Stauber, J.; Wilson, S.; Duffy, B.; Sefiane, K. On the lifetimes of evaporating droplets. J. Fluid Mech. 2014, 744, R2. (34) Haut, B.; Colinet, P. Surface-tension-driven instabilities of a pure liquid layer evaporating into an inert gas. J. Colloid Interface Sci. 2005, 285, 296−305.
(35) Cazabat, A.-M.; Guéna, G. Evaporation of macroscopic sessile droplets. Soft Matter 2010, 6, 2591−2612. (36) Sultan, E.; Boudaoud, A.; Ben Amar, M. Evaporation of a thin film: diffusion of the vapour and Marangoni instabilities. J. Fluid Mech. 2005, 543, 183−202. (37) Russel, W. B.; Saville, D. A.; Schowalter, W. R. Colloidal Dispersions; Cambridge University Press: Cambridge, U.K., 1995. (38) Brown, P. N.; Hindmarsh, A. C.; Petzold, L. R. Using Krylov methods in the solution of large-scale differential-algebraic systems. SIAM J. Sci. Comput. 1994, 15, 1467−1488. (39) Huppert, H. E. Flow and instability of a viscous current down a slope. Nature 1982, 300, 427−429. (40) Overdiep, W. S. The effect of a reduced solvent content of solvent-borne solution paints on film formation. Prog. Org. Coat. 1986, 14, 1−21. (41) Overdiep, W. S. The levelling of paints. Prog. Org. Coat. 1986, 14, 159−175. (42) Croll, S. G.; Kisha, L. W. Observations of sagging in architectural paints. Prog. Org. Coat. 1992, 20, 27−52. (43) Wray, A.; Papageorgiou, D.; Craster, R.; Sefiane, K.; Matar, O. Electrostatic suppression of the “coffee stain effect”. Langmuir 2014, 30, 5849−5858. (44) Fraštia, L.; Archer, A. J.; Thiele, U. Modelling the formation of structured deposits at receding contact lines of evaporating solutions and suspensions. Soft Matter 2012, 8, 11363−11386. (45) Jeong, H. J.; Hwang, W. R.; Kim, C.; Kim, S. J. Numerical simulations of capillary spreading of a particle laden-droplet on a solid surface. J. Mater. Process. Technol. 2010, 210, 297−305. (46) Hewitt, D. R.; Balmforth, N. J. Thixotropic gravity currents. J. Fluid Mech. 2013, 727, 56−82. (47) Trokhymchuk, A.; Henderson, D.; Nikolov, A.; Wasan, D. T. A simple calculation of structural and depletion forces for fluids/ suspensions confined in a film. Langmuir 2001, 17, 4940−4947. (48) Kondiparty, K.; Nikolov, A. D.; Wasan, D.; Liu, K. Dynamic spreading of nanofluids on solids. Part I: experimental. Langmuir 2012, 28, 14618−14623. (49) Liu, K.; Kondiparty, K.; Nikolov, A. D.; Wasan, D. Dynamic spreading of nanofluids on solids. Part II: modeling. Langmuir 2012, 28, 16274−16284. (50) Wasan, D.; Nikolov, A.; Kondiparty, K. The wetting and spreading of nanofluids on solids: role of the structural disjoining pressure. Curr. Opin. Colloid Interface Sci. 2011, 16, 344−349.
11974
dx.doi.org/10.1021/la503229z | Langmuir 2014, 30, 11966−11974
