Article pubs.acs.org/Langmuir
Electrostatic Suppression of the “Coffee Stain Effect” Alexander W. Wray,*,† Demetrios T. Papageorgiou,‡ Richard V. Craster,‡ Khellil Sefiane,§ and Omar K. Matar† †
Department of Chemical Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, U.K. Department of Mathematics, Imperial College London, South Kensington Campus, London SW7 2BZ, U.K. § Institute for Materials and Processes, School of Engineering, The University of Edinburgh, King’s Buildings, Mayfield Road, Edinburgh EH9 3JL, U.K. ‡
ABSTRACT: The dynamics of a slender, evaporating, particle-laden droplet under the effect of electric fields are examined. Lubrication theory is used to reduce the governing equations to a coupled system of evolution equations for the interfacial position and the local, depthaveraged particle concentration. The model incorporates the effects of capillarity, viscous stress, Marangoni stress, elecrostatically induced Maxwell stress, van der Waals forces, concentration-dependent rheology, and evaporation. Via a parametric numerical study, the one-dimensional model is shown to recover the expected inhomogeneous ring-like structures in appropriate parameter ranges due to a combination of enhanced evaporation close to the contact line, and resultant capillarity-induced flow. It is then demonstrated that this effect can be significantly suppressed via the use of carefully chosen electric fields. Finally, the three-dimensional behavior of the film and the particle concentration field is briefly examined. having a flux that decays exponentially near the contact line,15,16 or by relaxing the no-slip criterion there.17 Hu and Larson18 have shown that the introduction of Marangoni flow can in fact suppress the ring effect. However, they do note that the recirculation effect is significantly weaker in experiments than as predicted by theoretical calculations, at least for the case of water. This effect has since been extensively investigated.19,20 There have since been studies of consistent lubrication-theory models both in the absence21 and presence22 of thermal effects. These used a depth-averaged concentration, an assumption relaxed by Maki and Kumar,14 with resultant skin formation. There is a very large body of existing literature on this topic. An extensive review has been given recently by Larson,23 to which we refer the interested reader for additional details. As mentioned above, the suppression of the coffee stain effect holds significant potential interest in practical and industrial settings. Thus, recent, predominantly experimental, results on the use of electric fields to suppress the ring stain effect hold a lot of interest for potential modeling. These results indicate that the “coffee stain” effect can be suppressed using either AC or DC electric fields. In the case of AC potential, the alternation and cycles of voltage lead to the oscillations of the contact line which inhibit the deposition of the particles in the same position.24 Mampallil et al.25 performed additional experiments to demonstrate the same mechanism over a range of parameters. In the case of DC fields, on the other hand, there are no oscillations of the contact line. Instead, the
1. INTRODUCTION It is a commonly experienced phenomenon that the evaporation of a nanoparticle-laden droplet gives rise to a distinctly inhomogeneous residue, even when the original particle distribution was entirely uniform. This results in what has come to be known as the “coffee stain” or “coffee ring” effect. It is of course commonly seen in situations such as drying coffee drops, or indeed in the drying of watercolor paint.1 The effect is of significant interest in a wide array of practical contexts. While there are many contexts in which the ring-like inhomogeneity is a useful patterning technique,2−4 it is also found to be an undesirable effect. This is the case, for instance, in DNA microarrays,5,6 chemical recovery,7 matrixassisted laser desorption/ionization mass spectrometry,8 and nanofabrication.9 This phenomenon has been the subject of significant investigation,10,11 with Deegan et al. explaining that the most important factors are an increased flux near the (pinned) contact line and a resultant capillarity-induced restoring flow. This results in an almost-complete mass flux of the particles in the bulk toward the contact line. In particular, Deegan et al.1 gave an explicit ODE based on a lubrication approximation assuming a spherical cap. Similarly, Hu and Larson12 have also given modeling results and demonstrate good agreement with their corresponding experiments. Hu and Larson13 have used a lubrication-like analysis. However, as pointed out by Maki and Kumar,14 this and other similar analyses suffer two flaws. First, the posited ansätze do not themselves satisfy the lubrication form of the momentum equations. Second, the resultant velocity profiles are singular at the contact line, due to a singularity in the evaporative flux there. It is possible to artificially remove this singularity by © 2014 American Chemical Society
Received: February 28, 2014 Revised: May 2, 2014 Published: May 2, 2014 5849
dx.doi.org/10.1021/la500805d | Langmuir 2014, 30, 5849−5858
Langmuir
Article
Figure 1. Left: (a) Photograph of stain rings resulting from the evaporation of a drop containing 0.1% titanium nanoparticles in water on CYTOP. (b) Stick−slip behavior exhibited by the contact angle and base radius evolution (corresponding to individual rings in part a). Right: (a) Photograph of a uniform deposit from a similar drop with electrowetting on dielectrics (EWOD) at a DC potential of 10 V. (b) Corresponding receding contact line, from Orejon et al.26 The data for the contact angles and radii in part b are measured using a DSA100 instrument (from Krüss Optronic) monitoring the drops shown in part a from right to left of the panels of the deposits shown. The 500 μm scale shown in part a is valid for both horizontal and vertical axes.
agglomeration on fluid flow, thermal properties, and surface tension are ignored, but we retain the dependence of viscosity on local particle concentration. A second parallel electrode, assumed without loss of generality to have zero potential, is located above the drop. The potential of the lower electrode is allowed to vary as a function of space and time. The geometry of the situation is as given in Figure 2. We assume that the fluid is sufficiently viscous that the hydrodynamics may be described by the Stokes equations
contact line recedes monotonically. The uniform deposition of nanoparticles is believed to be the result of the combination of the contact line receding and electrophoretic force driving particles toward the substrate;26 this is shown in Figure 1. However, despite the evident interest and breadth of work in the literature, the existing studies of ring suppression have all been experimental, and we wish to redress the balance. Our aim, therefore, is to investigate situations where long-wave theory reproduces the coffee ring effect, investigating the mechanisms at work. We also wish to investigate how suppression of the coffee stain effect might be attempted using electric fields, via the use of leaky-dielectric theory,27 as has previously been used successfully in such situations.28 In section 2, we outline the problem that we shall be studying, before using a long-wave approximation to simplify its governing equations. In section 3, we discuss the reproduction of the ring effect in the current model, before studying its suppression and enhancement using electric fields. Our conclusions are given in section 4.
∇p = μ∇2 u ,
2. PROBLEM FORMULATION 2.1. Governing Equations and Nondimensionalization. We model a slender droplet of density ρ and viscosity μ deposited on an ultrathin film on a flat, horizontal electrode. The droplet is laden with nanoparticles of density cρ. We use the Cartesian coordinate system (x, y, z) with y oriented so as to be normal to the electrode, with the origin resting at the center of the droplet, and the interface residing at y = h. The velocities in the respective directions are u = (u, v, w). We assume that gravity is negligible. The effects of particle
∇·u = 0
(1)
Figure 2. Schematic of the system geometry. Above: a droplet resting on a substrate with a precursor film, with a parallel electrode above. Below: a possible wall potential ϕw. 5850
dx.doi.org/10.1021/la500805d | Langmuir 2014, 30, 5849−5858
Langmuir
Article
Table 1. Table of All Problem Parameters and Characteristic Values symbol ρl μ0 cρ Ts γ0 σ0 A ΔHv λ R ρv + σR ε0 εG εL 3 symbol
dimensional quantity
characteristic value
ε V
aspect ratio velocity
1.00 × 10−1 1.00 × 10−4 m s−1
(
dimensionless Hamaker constant
ε × 5.3 × 10−9
Ca͠ Pe
capillary number
ε−3 × 1.39 × 10−6
Péclet number
ε−1 × 2.07 × 105
superheat electric Weber number
1.00 × 10 εϕb2 × 8.85 × 10−3V−2
Marangoni number kinetic parameter pressure contribution to phase change
ε × 6.27 × 10 2.30 × 10−3 ε−2 × 1.00 × 10−13
Θ Eb Ma K Δ
n ·T·t = t·∇s γ
A 6πh3
λTs/(ρl ΔH v 3)
εA /(6π 32Vμ0 ) ε−3μ0V/γ0
ε−1V 3 /+ (Tw − Ts)/Ts εε0ϕb2 /(μ0 V 3) ε2σ0Ts/(Vμ0) (2πRTs)1/2ρlV/(ρvΔHv)
5
ht + uhx − v +
(2)
where n and t are the normal and tangential vectors to the interface, respectively, Π is the disjoining pressure, T = −pI + 4e + τ is the total stress tensor, and ∇s = ∇ − n(n·∇) is the surface gradient operator. Here, 4e is the Maxwell stress due to the presence of electric fields, τ = μ(∇u + ∇uT) is the usual viscous stress tensor, and γ is the surface tension. We allow the surface tension to vary as a linear function of the temperature22 (but neglect its dependence on particle concentration29) so that γ = γ0 − σ0(T − Ts), where γ0 is the surface tension at the saturation temperature Ts, and σ0 = −dγ/dT. The Maxwell stress is given by27,30 4e = εEE − (1/2)ε|E|2I, where ε is the local permittivity, E is the electric field, and I is the identity matrix. The viscosity μ is allowed to be a function of the local particle concentration c. In particular, following Craster et al.,22 we take μ = μ0(1 + c/cρ)2. We ignore particle interactions in the disjoining pressure, instead incorporating only van der Waals effects, giving
Π=
m grouping
ε−2Vμ0 /(ρl ΔH v 3)
In addition, at the interface the fluid satisfies the kinematic condition,
where p is the pressure. These are solved subject to no-slip and impermeability at the lower wall, and the stress conditions at the interface n ·T·n = Π − γ ∇s ·n ,
m−3 kg−1 s4 A2
0
2
unit kg m−3 Pa s J(kg K)−1 K N/m N(m K)−1 J J kg−1 W(m K)−1 J K−1 kg−1 kg m−3 m2 s−1
1.00 × 103 1.00 × 10−3 4.18 × 103 3.73 × 102 7.2 × 10−2 1.68 × 10−4 1.00 × 10−20 2.26 × 106 6.07 × 10−1 4.61 × 102 2 × 10−2 4.84 × 10−13 1.00 × 1012 8.85 × 10−12 1.00 × 10° 8.00 × 101 1.00 × 10−3 characteristic value
liquid density viscosity particle density saturation temperature surface tension surface tension temperature dependence Hamaker constant latent heat of vaporization thermal conductivity gas constant vapor density particle diffusion coefficient conductivity ratio permittivity of free space relative permittivity of gas phase relative permittivity of liquid phase height of droplet name
J =0 ρl
(4)
where ρl is the liquid density and J is the evaporative flux. The temperature distribution T is assumed to be conduction dominated, so that it satisfies Laplace’s equation ∇2 T = 0
(5)
This is solved subject to a thermal flux condition at the interface and temperature equality at the wall, given, respectively, by ΔH vJ = −λ n ·∇T ,
T |y = 0 = Tw
(6)
where ΔHv is the latent heat of vaporization, λ is the thermal conductivity, and Tw is the temperature of the wall, and we have neglected heat transfer effects (i.e., we are working in the limit of small Biot number). To close this system, the following nonequilibrium interfacial condition is applied:31,32 pve pv
(3)
where A is the Hamaker constant. This is in line with the modeling of Ajaev,31 with the incorporation of van der Waals forces allowing for the existence of a stable ultrathin precursor film, removing the singularity issue at the contact line.
−1=
⎞ ΔH v ⎛ Ti 1 (p − pv ) + ⎜ − 1⎟ RTs ⎝ Ts ρl RTs ⎠
(7)
where Ti is the interfacial temperature, R is the gas constant per unit volume, pv is the pressure of the vapor phase, and pve is the equilibrium vapor pressure. The saturation temperature Ts satisfies31,32 5851
dx.doi.org/10.1021/la500805d | Langmuir 2014, 30, 5849−5858
Langmuir
Article
⎛ ⎞2 ⎜ ⎟ J 2π ⎜ ⎟ Ts = ⎞ R ⎜⎜ ⎜⎛ pve ρv p − 1⎟ ⎟⎟ ⎠⎠ ⎝ ⎝ v
decoration on the variables and work in dimensionless variables hereafter. All the parameters in this problem and their characteristic values are given in Table 1. 2.2. Long-Wave Reduction. We look for a model in the long-wave regime, so we make the substitution x = x̃/ε, where x̃ is assumed to be of order unity and ε is an aspect ratio, given by the ratio of the height of the drop to its radius. We assume that the contact angle is sufficiently small for the long-wave reduction to be appropriate. In addition, we assume a capillarity driven flow, so we scale the capillary number Ca = μ0V/γ0 such that Ca = ε3Ca͠ , where Ca͠ is of order unity (and indeed we take Ca͠ = 1 throughout). Consideration of the continuity equation, kinematic equation, and normal stress suggest scaling v = εṽ, t = t/̃ ε, and p = 6(ε−1). Under these assumptions, the Stokes equations (eq 1) reduce to
(8)
where ρv is the vapor density. The ratios of the thermal conductivity, density, viscosity, etc., of the vapor to the liquid are all neglected in this one-sided model, and we do not consider explicitly the dynamics in the vapor. The particle concentration c is assumed to satisfy an advection−diffusion equation
Dc = +∇2 c Dt
(9)
where + is the diffusion coefficient of the nanoparticles. This is subject to no flux at the wall, cy = 0, and an interfacial diffusion flux due to concentration variation associated with evaporation: cJ +n ·∇c = c(u − us) ·n = ρl
u x ̃ + vỹ = 0,
where us is the velocity of the interface. During the course of evaporation, particles may be drawn into the precursor layer; in practice, the particle diameters may be at least comparable with the thickness of this ultrathin precursor film. Additional modeling of the dynamics in such a situation, especially with regard to impact on the contact line, could be of significant future interest; for the moment, such effects are neglected. The potentials in the liquid and gas regions, denoted by ϕL and ϕG, respectively, satisfy Laplace’s equation
ϕyyL,G = 0
−q = εGϕyG|y = h − εLϕyL|y = h
c t ̃ + uc x ̃ + vc̃ y =
c0t ̃ + uc0x ̃ + vc̃ 0y =
where ε0 is the permittivity of free space and εG and εL are the usual relative permittivities of the gas and liquid phases, respectively. The quantities are nondimensionalized using22,30
t=
3 ̂ t, V
T − Ts = T̂ , Ts μ = μ0 μ ̂
q=
ε0ϕb 3
c = cρc ̂,
3
q̂ ϕ
= ϕ bϕ ̂
L,G
,
1 (c0xx̃ ̃ + c1yy) Pe͠
This is then cross-sectionally averaged using upper interface, eq 10 gives, with J = εJ,̃
(21)
(1/h)∫ h0 dy.
c1y = hxc0x + c0J ̃Pe
p̂ ,
At the (22)
In order to discern the form of J, we note that, under the rescaling of the deviation temperature T = εT̃ , the temperature satisfies T̃ yy = 0 subject to (14)
L,G
(20)
where Pe = V 3/+ is the Péclet number. Following Jensen and Grotberg,33 we assume the Péclet number is small, setting Pe = εPe͠ , and assume rapid vertical diffusion, setting c = c0(x, t) + ε2c1(x, y, t) + .... Then, c is governed by
(13)
p=
1 2 (ε c xx̃ ̃ + cyy) εPe 22
where σR is the ratio of the electrical conductivity of the liquid to that of the gas. The charge q can be derived from the Gauss condition
u = V u,̂
(19)
The particle concentration c is governed by
(12)
Vμ0
(18)
The Gauss condition (eq 13) becomes
σR (n ·∇ϕL)|y = h = (n ·∇ϕG)|y = h
̂ (x , y , h) = 3(x ̂, y ̂, h),
(17)
σR ϕyL|y = h = ϕyG|y = h
subject to appropriate equipotentials at the upper and lower electrodes, respectively, ϕG|y=d = 0 and ϕL|y=0 = ϕw(x, t). This is complemented by the continuity of potential and current at the interface
−q = ε0(εG∇ϕG − εL∇ϕL)|y = h ·n
(16)
The continuity of potential at the interface (eq 12) remains unchanged; the continuity of current becomes
(11)
ϕG| y = h = ϕL| y = h ,
pỹ = 0
Before completing the exposition of the hydrodynamic problem, we consider the equations governing the other variables to discern their scalings. The scaled version of Laplace’s equation for the potentials (eq 11) becomes
(10)
∇2 ϕL,G = 0
px̃ ̃ = [μ(c)uy]y ,
Tỹ = 0 = (Tw − Ts)/Ts ≡ Θ,
J = ρl VJ ̂,
Tỹ |y = h = −J ̃
(23)
where Θ represents a dimensionless difference between the wall temperature and the equilibrium saturation temperature. Thus, finally, we complete the hydrodynamic problem. The velocities satisfy no-slip and impermeability at the wall, u = ṽ = 0. At the interface, the normal stress condition is given by
(15)
where ϕb is the maximum value of the potential of the wall electrode, μ0 is the viscosity of the liquid in the absence of particles, 3 is the initial height of the drop, and V = λTs/ ρlΔHv3 is a characteristic velocity.22 We now drop the hat
p ̃ = pv − 5852
2 2 h xx̃ ̃ E ( − b (εGϕyG − εLϕyL ) − 3 Ca͠ 2 h
(24)
dx.doi.org/10.1021/la500805d | Langmuir 2014, 30, 5849−5858
Langmuir
Article
where E b /ε = ε0ϕb 2 /(μ0 V 3) measures the relative strength of
p = pv − hxx − Π −
the electric field 30 and (/ε = A/(6π32 Vμ 0 ) is the dimensionless rescaling of the Hamaker constant,22 both scaled to ensure retention at leading order. The tangential stress condition is given by ̃ ̃ μuy|y = h = −E bq(ϕx ̃ + h xϕ ̃ y) − Ma(Tx ̃ + h x T ̃ y)
⎡ h2c ⎛ εGϕw ⎞ 1 3 ⎤ h cpx ⎥ ϕw ⎟ − (hc)t + ⎢ ⎜Ma(Jh)x + E b h − d x⎠ 3μ ⎣ 2μ ⎝ ⎦x
(25)
=
(27)
Θ + Δp K+h
⎡ h2c ⎛ ⎤ ⎞ εGϕw 1 3 ∇ϕw ⎟ − (hc)t + ∇·⎢ ⎜Ma∇(Jh) + E b h c ∇p ⎥ 3μ h−d ⎠ ⎣ 2μ ⎝ ⎦
(28)
=
where K = (2πRTs)1/2ρlV/(ρvΔHv) and Δ = ε−2Vμ0/(ρlHv3 ) are dimensionless constants that denote the relative significance of interfacial kinetic effects and of variation in the pressure on the phase-change interfacial temperature. Solving the leading order of Laplace’s eq 17 subject to eq 18 and equipotential at the electrodes gives ϕL = ϕw
σ R (y − d ) , σR (h − d ) − h
q = ϕw
εL − σR εG σ R (h − d ) − h
We take σR ≫ 1 so that ϕG = ϕw
y−d , h−d
q=
ϕw εG d−h
(38)
hc(x , 0) = 0.1 × h(x , 0)
(39)
0= (31)
⎛ Δ( ⎞1/3 Θ − Δ(/h∞3 ⎟ ⇒ h∞ = ⎜ ⎝ Θ ⎠ K + h∞
(40)
This has been solved using a custom numerical code that is first order semi-implicit in time and second order in space. For validation, the alternate systems of eqs 32 and 33 were solved using EPDCOL34 (with c(x, 0) = 0.1) and the results compared for accuracy. In addition, the grids in each case and the temporal error control parameter were varied to check for convergence. For the calculation of the 2D eqs 36 and 37, we use a semiimplicit scheme based on the pseudolinear scheme (pL1) of Witelski and Bowen,35 which is first order in time and second order in space. We use initial conditions which are the natural axisymmetric extension of eqs 38 and 39. 3.2. Reproduction of the Ring Effect. We first wish to validate whether or not the experimentally observed ring
Dropping the suffix of c0, the evolution equations are then ⎡ h2 ⎛ εGϕw ⎞ 1 3 ⎤ ϕw ⎟ − ht + ⎢ ⎜Ma(Jh)x + E b hp⎥ + J = 0 h − d x⎠ 3μ x ⎦ ⎣ 2μ ⎝ x (32)
⎡ h ⎛ εGϕw ⎞ 1 2 ⎤ ϕw ⎟ − ct + ⎢ ⎜Ma(Jh)x + E b h p ⎥cx x 3μ x ⎦ h−d ⎠ ⎣ 2μ ⎝ cJ 1 (hcx)x + = h hPe
2 3 ⎧ ⎪(1 − x ) + h ∞ |x | ≤ 1 h(x , 0) = ⎨ ⎪ |x | > 1 ⎩ h∞
This corresponds to a drop laden with particles released onto a thin film of uniform thickness h∞. The thickness of the precursor film h∞ is selected by considering eq 28 and arranging that, with a uniformly zero electric field, the flux J is zero, so that
(30)
We can then see that the tangential term τ := −Ebq(ϕx + hxϕy) = Eb[εGϕw/(h − d)]ϕwx. Next we solve eq 16 subject to no-slip and eq 25. This gives ⎞ ⎛ y2 μu = px ⎜ − hy⎟ + Ma(Jh)x y + τy ⎠ ⎝2
(37)
3. RESULTS 3.1. Numerical Procedure. We wish to solve eqs 32 and 35 subject to the initial conditions (29)
ϕL = ϕw ,
∇·(h∇c) Pe
We note that, for the purpose of investigating solute distribution, we find it more helpful not to consider c but rather its depth-integrated value ∫ h0 c dy = hc. The variation of this conserved quantity reflects the mass transport throughout the drop and shall be studied extensively in the next section.
σ R (h − d ) − (h − y ) , σR (h − d ) − h
ϕG = ϕw
(36)
=0
2.3. Solution of Long-Wave Equations. Equation 7 together with eq 23 can be solved to give J=
(35)
⎡ h2 ⎛ ⎞ εGϕw 1 3 ⎤ ∇ϕw ⎟ − ht + ∇·⎢ ⎜Ma∇(Jh) + E b h ∇p ⎥ + J 3μ h−d ⎠ ⎣ 2μ ⎝ ⎦
where ũ = (1/h)∫ h0 u dy. With these definitions, the kinematic condition is
T = Θ − yJ ,
(hcx)x Pe
These can naturally be extended in the expected way to three dimensions, giving
(26)
ht + (uh ̅ )x + J = 0
(34)
We choose, instead of 33, to rewrite this in terms of the conserved quantity hc:
where Ma/ε2 = σ0Ts/V is the Marangoni number.22 Putting all this together, and dropping the text decoration, gives 1 c0t + uc̅ 0x = (hc0x)x hPe
2 E bεG ⎛ ϕw ⎞ ⎜ ⎟ 2 ⎝d − h⎠
(33)
where J is given by eq 28 and p is given by 5853
dx.doi.org/10.1021/la500805d | Langmuir 2014, 30, 5849−5858
Langmuir
Article
Figure 3. Maps of the final distribution of the depth-integrated concentration hc as functions of x and Θ at the point of termination (h = h∞ everywhere). Thus, as seen in the inset images, which correspond to the respective dotted lines, a particular value of the vertical coordinate Θ gives a distribution hc as a function of x.
Figure 4. Top row: interfacial fluxes; left, interfacial flux due to the combination of superheat and van der Waals forces; right, interfacial flux due to capillarity-induced phase change due to the interfacial pressure jump. Bottom row: internal flow fields; left, the flow induced by capillarity; right, the flow that would be induced for Ma ≠ 0.
coordinate corresponds to a different value of Θ, while the horizontal axis corresponds to the spatial position x. We plot the value of hc corresponding to this (Θ, x) pair as a density map. It is seen that higher values of Θ give rise to greater maximal values of hc, while low values of Δ and K give rise to more pronounced ring formations, as seen in the right-hand column. For example, the inset of the image corresponding to K = 0.01, Δ = 0.001 shows a distinct transport of the solute away from
depositions can be reproduced by this model. Therefore, we investigate parametrically the effect of varying the superheat 0.05 ≤ Θ ≤ 1, and the relative significances of the interfacial pressure-induced phase change effects Δ ∈ {0.001, 0.1} and interfacial kinetic effects K ∈ {0.1, 0.01}. We therefore simulate the system up until the point where the film becomes flat (h = h∞ everywhere), and at this point consider the depth-integrated concentration hc. We demonstrate this value of hc at the point of termination in Figure 3. In each image, the vertical 5854
dx.doi.org/10.1021/la500805d | Langmuir 2014, 30, 5849−5858
Langmuir
Article
Figure 5. Maps corresponding to the final hc distribution in x−Θ planes with ϕw = ΦG (eq 41).
previously reported,18 Marangoni effects will in general oppose the ring formation. As we intend to suppress the effect using electric fields, we allow Ma = 0 from here on. We note that there is another phase of the evolution that these mechanisms do not really cover: the initial transient. In particular, the drop initially spreads/retracts slightly under the effect of surface tension to ensure an interface of uniform curvature, and to match to the precursor film. However, this phase exists in the purely hydrodynamic case and so is not of relevance here. 3.3. Parametric Study of Electric Fields. We wish to demonstrate the degree of control that the introduction of an electric field affords over the final distribution of hc. The option of an arbitrary choice of ϕw gives us an enormous degree of freedom: in order to simplify this, we restrict ourselves to particular families of solutions. In particular, as has been previously seen,28 fields with a Gaussian potential distribution allow for a significant degree of control. We therefore start by examining the spatially symmetric, two-parameter family of potentials
the center line to accumulate in two symmetric peaks. Thus, the right-hand column reproduces the expected ring effect to some extent. However, in contrast to experimental observations, we note that the peaks are at x ≈ 0.35, some distance from the initial contact line of x ≈ 0.9. This is due largely to the fact that, in our model, the contact line is not pinned, and indeed moves freely due to evaporation near the contact line. This discrepancy and potential resolutions are discussed further in section 3.5. In addition, the expected ring stain is absent from the left column. In order to elucidate on the mechanisms at work, we examine a drop in detail for Δ = 0.001, K = 0.01, Θ = 0.6, ( = 10−6, Ma = 0, and Pe = 103 at t = 0.26, which is typical of the drop evolution we have seen. First, we examine the fluxes in the top row of Figure 4. On the left, we can see the contribution of the evaporative flux due to the terms (Θ − ΔΠ)/(K + h), while, on the right, we see the contribution of the pressure to phase change −Δhxx/(K + h). As expected, we get a sharp increase in the magnitude of the fluxes at the edges of the drop. Via pressure-induced phase change, capillarity works to suppress evaporation. Thus, a smaller value of Δ results in flux being dominant close to the contact line, resulting in a more pronounced coffee stain effect (right-hand column of Figure 3). A larger value results in a more uniform evaporation profile; in combination with a receding contact line, this results in little evidence of a ring stain (left-hand column of Figure 3). We examine the flow field inside the drop in the bottom row of Figure 4. The most important constituent here is the flow outward due to capillarity. This is due to interfacial shape change due to evaporation, essentially as reported by Maki and Kumar.14 It is noted that, for the majority of the bulk, the contribution of van der Waals forces is weak, although in the precursor film it serves to prevent rupture. On the right, we plot the flow that would arise due to Marangoni forces. In this case, we have neglected them by setting Ma = 0, but the flow field that would be induced is easily computed for a given configuration of h and c, and this is what we plot here. However, the nature of the field seems to be generic: as
ϕw (x) = ΦG(x ; :G, 3 G) = exp( −:G(x − 3 G)2 ) + exp( −:G(x + 3 G)2 ) (41)
We take d = 2 throughout. We also have the freedom to select Eb as we choose. To allow the drop to settle from its artificial initial condition before it starts being strongly influenced by the field, we choose E b = ,(1 − exp( − 10t ))
(42)
We anticipate that, by having the Gaussians near the edge of the drop, the bulk of the fluid may be driven back toward the center. We therefore examine the parameter ranges 0.45 ≤ 3 G ≤ 0.75, :G ∈ {5, 15}, and , ∈ {150, 300}. We note that the first neglected electrostatic terms in eqs 17−19, 24, and 25 are of 6(ε 2). Thus, care should be taken, as a sufficiently large 5855
dx.doi.org/10.1021/la500805d | Langmuir 2014, 30, 5849−5858
Langmuir
Article
value of , potentially challenges the fact that these terms have been neglected. The results are given in Figure 5. Briefly, we note that with the inclusion of an electric field the interface may no longer be quite uniformly h = h∞ at the point where evaporation terminates, due to the balance of a nonuniform field with evaporative and van der Waals effects, though the difference is small. As is clearly visible, it is possible to significantly affect the final distribution of hc via the use of electric fields. We observe that the interface will deviate in such a way as to drive the bulk of the drop away from the peaks of the electric field. We note that this characteristic behavior is in fact principally due to the tangential component of the Maxwell stress induced at the interface; this tangential stress is caused by the variation of the wall potential ϕw, and is not exhibited when the wall potential is constant. Hence, as predicted, this “two-prong” field allows the drop to be forced toward the center of the domain. In order to quantify this, we consider the variance Var(hc) = ∫ x2hc dx. The parametric dependence of this on 3 G is given in Figure 6 for
the case :G = 5, , = 300. A decrease in Var(hc) is observed as 3 G is decreased from its initial value of 0.7. However, it is noted that after a certain point the value of 3 G lies within the body of the drop for a nontrivial portion of the period of evaporation. This results in some of the fluid being driven outward, hence increasing the variance once again. As a result, we get a minimum value for some intermediate value of 3 G . This gives rise to the lowest value of the variance for any variant of ϕG calculated, with Var(hc) ≈ 1.08 for 3 G ≈ 0.6. Clearly the shape of a pair of Gaussian curves is not optimal, and we would expect to have a greater degree of success without the outer “downslope” of the curves. We therefore investigate a new family, ϕw (x) = ΦT (x ; :T , 3T ) = 2 + tanh( −:T (x − 3T )2 ) + tanh( −:T (x + 3T )2 )
(43)
We examine the parameter ranges 0.4 ≤ 3T ≤ 0.75, :T ∈ {3, 10}, and , ∈ {150, 300} in Figure 7. It is seen that having the tanh curves less steeply sloped gives rise to a better suppression of the ring effect, while a stronger field and narrower placement of the curves almost augments this phenomenon. In particular, the best suppression is given by the case :T = 3, , = 300. The parametric dependence of Var(hc) on 3T is given in Figure 6. We note that the single best suppression is given with 3 = 0.4 , giving a variance of 0.898 as opposed to 1.727 in the electrically passive case: a reduction of 48%. 3.4. Enhancement of the Ring Effect Using Electric Fields. It is also of some practical interest to examine whether the ring effect can be enhanced via the use of electric fields. In particular, for the parameters
Figure 6. Parametric dependence of Var(hc) on 3 G with ϕw = ΦG for :G = 5, , = 300 (dotted line) and dependence on 3T with ϕw = ΦT for :T = 3, , = 300 (solid line). Plotted for comparison: Var(hc) in the electrically passive case (dashed line).
Figure 7. Maps corresponding to the final hc distribution in x−Θ planes with ϕw = ΦT (eq 43). 5856
dx.doi.org/10.1021/la500805d | Langmuir 2014, 30, 5849−5858
Langmuir
Article
Figure 8. Results of one-dimensional calculations for the parameters (eq 44) subject to the field (eq 45). Left, t = 0.5; right, t = 1.4. The drop is divided symmetrically in two by the field.
Figure 9. Results of the two-dimensional calculations for the parameters (eq 44) subject to the field (eq 46). Left, t = 0.5; middle, t = 0.9; right, t = 1.2. The drop, after initially dividing into a symmetric ring, subsequently collapses into a pair of droplets.
Pe = 200,
( = 10−6 ,
Θ = 0.2,
radically). This is not the mechanism observed by Eral et al.,24 where a moving contact line is indeed observed. Another alternative would have been to use a stick−slip model for the contact line, and this is an alternative that we would like to pursue in future. We also note that, where other studies1,10,13,18 imposed the condition that the drop interface assumed the shape of a spherical cap, we were able to relax this assumption. This is necessary because the multitude of flow effects, especially with the introduction of the electric fields, has rendered this assumption inappropriate. Thus, while in the electrically passive case capillarity typically acts to hold the drop in a surface of constant curvature (see Figure 4), this can be radically violated by electric fields (see Figure 9). Finally, we note that the suppression mechanism we have described in section 3.3 is somewhat different from that described in the experimental context of Orejon et al.26 Their mechanism is based to a large degree on electrophoretic effects (which we have not incorporated here) and a particular electric field that is somewhat different from the fields we have prescribed. We note that in fact the long-wave formulation we have used could not be used to model the field induced in the geometry of Orejon et al.26 (although it would be quite manageable to model that exact situation via the use of a conformal map to a simpler geometry, and solving the flow field exactly using complex methods).
Δ = 0.005, (44)
K = 0.2
we impose the fields ϕw = exp( −10x 2),
E b = 10 exp( −10t )
(45)
In our one-dimensional calculations, as seen in Figure 8, the electric field was observed to divide the drop exactly as expected. The process is entirely symmetric, and the two resultant drops evaporate individually. However, if we attempt the same procedure in two dimensions for the field, ϕw = exp( − 10(x 2 + y 2 ))
(46)
we find the results given in Figure 9: a capillary instability destabilizes the structure, and results in a collapse of the fluid ring, heading toward two individual droplets. This result is not entirely unsurprising but does highlight the necessity of performing the full three-dimensional calculations. Note that this does not exclude the potential possibility of using electric fields to selectively increase/decrease evaporation in certain regions of the drop. However, this would be difficult to manipulate without having a dramatic effect via the Maxwell stresses at the interface. 3.5. Comparison of Results to Previous Studies and Experiments. As expected, the mechanism that manifests here is closer to that described by Maki and Kumar14 than that of Deegan et al.10 This is for two main reasons: first, the fact that we have allowed for a precursor film instead of a contact line model and, second, Deegan et al.’s assumption of perfect interfacial sphericity. Using a precursor film model, it is very unlikely that a perfect ring formation will result (and indeed, the flux at the edge of the drop immediately forces it to start receding in practice). However, we note that it would have been inappropriate to use a fixed contact line model as used elsewhere,10,16 as this model almost inevitably forces complete mass transport to the contact line (unless the robust evaporation profile can be altered
4. CONCLUSIONS The behavior of a nanoparticle-laden droplet deposited onto a heated substrate has been studied, both in the presence and absence of electric fields. An asymptotic long-wave expansion has been used to derive a canonical nonlinear model that retains many physical effects (rheological dependence on local particle concentration, viscous effects, capillarity, Marangoni stresses, and Maxwell stresses). The leaky dielectric model has been reduced to the limit where the liquid phase is highly conducting to afford a 5857
dx.doi.org/10.1021/la500805d | Langmuir 2014, 30, 5849−5858
Langmuir
Article
(14) Maki, K. L.; Kumar, S. Fast evaporation of spreading droplets of colloidal suspensions. Langmuir 2011, 27, 11347−11363. (15) Masoud, H.; Felske, J. D. Analytical solution for Stokes flow inside an evaporating sessile drop: Spherical and cylindrical cap shapes. Phys. Fluids 2009, 21, 042102. (16) Fischer, B. J. Particle convection in an evaporating colloidal droplet. Langmuir 2002, 18, 60−67. (17) Petsi, A.; Burganos, V. Stokes flow inside an evaporating liquid line for any contact angle. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2008, 78, 036324. (18) Hu, H.; Larson, R. G. Marangoni effect reverses coffee-ring depositions. J. Phys. Chem. B 2006, 110, 7090−7094. (19) Girard, F.; Antoni, M.; Faure, S.; Steinchen, A. Evaporation and Marangoni driven convection in small heated water droplets. Langmuir 2006, 22, 11085−11091. (20) Dunn, G.; Wilson, S.; Duffy, B.; David, S.; Sefiane, K. The strong influence of substrate conductivity on droplet evaporation. J. Fluid Mech. 2009, 623, 329−351. (21) Matar, O. K.; Craster, R. V.; Sefiane, K. Dynamic spreading of droplets containing nanoparticles. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2007, 76, 056315. (22) Craster, R. V.; Matar, O. K.; Sefiane, K. Pinning, retraction, and terracing of evaporating droplets containing nanoparticles. Langmuir 2009, 25, 3601−3609. (23) Larson, R. G. Transport and deposition patterns in drying sessile droplets. AIChE J. 2014, 60, 1538−1571. (24) Eral, H. B.; Augustine, D. M.; Duits, M.; Mugele, F. Suppressing the coffee stain effect: how to control colloidal self-assembly in evaporating drops using electrowetting. Soft Matter 2011, 7, 4954− 4958. (25) Mampallil, D.; Eral, H. B.; van den Ende, D.; Mugele, F. Control of evaporating complex fluids through electrowetting. Soft Matter 2012, 8, 10614−10617. (26) Orejon, D.; Sefiane, K.; Shanahan, M. E. Evaporation of Nanofluid Droplets with Applied DC Potential. J. Colloid Interface Sci. 2013, 407, 29−38. (27) Saville, D. Electrohydrodynamics: the Taylor-Melcher leaky dielectric model. Annu. Rev. Fluid Mech. 1997, 29, 27−64. (28) Yeo, L.; Craster, R.; Matar, O. Drop manipulation and surgery using electric fields. J. Colloid Interface Sci. 2007, 306, 368−378. (29) Sefiane, K.; Skilling, J.; MacGillivray, J. Contact line motion and dynamic wetting of nanofluid solutions. Adv. Colloid Interface Sci. 2008, 138, 101−120. (30) Wray, A. W.; Papageorgiou, D. T.; Matar, O. K. Electrified coating flows on vertical fibres: enhancement or suppression of interfacial dynamics. J. Fluid Mech. 2013, 735, 427−456. (31) Ajaev, V. S. Spreading of thin volatile liquid droplets on uniformly heated surfaces. J. Fluid Mech. 2005, 528, 279−296. (32) Moosman, S.; Homsy, G. M. Evaporating menisci of wetting fluids. J. Colloid Interface Sci. 1980, 73, 212−223. (33) Jensen, O. E.; Grotberg, J. B. The spreading of heat or soluble surfactant along a thin liquid film. Phys. Fluids A 1993, 5, 58−68. (34) Keast, P.; Muir, P. H. Algorithm 688: EPDCOL: A More Efficient PDECOL Code. ACM Trans. Math. Software 1991, 17, 153− 166. (35) Witelski, T.; Bowen, M. ADI schemes for higher-order nonlinear diffusion equations. Appl. Numer. Math. 2003, 45, 331−351.
greater degree of analytical tractability. The resultant equations have allowed us to carry out a wide range of parametric studies in order to understand the observed phenomena. It has been demonstrated that the experimentally observed ring configurations can be recovered in our model. It was also observed that this effect could be manipulated via the use of carefully selected electric fields. However, it was noted that the one-dimensional calculations presented here do not always extend naturally to the expected axisymmetric analogue, so that the ability to derive, analyze, and simulate the governing equation in two dimensions is of great importance. We also point out that any effects ignored in our model (such as disjoining pressure due to particle concentration or imperfect conduction effects) can readily be incorporated. Finally, we reiterate that the simplicity of our model, together with its ability to reproduce experimental results, renders it very useful in potential industrial applications.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■ ■
ACKNOWLEDGMENTS Work partly supported by NSF grant DMS-0707339 to D.T.P. and an EPSRC Doctoral Training Award to A.W.W. REFERENCES
(1) 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: Stat., Nonlinear, Soft Matter Phys. 2000, 62, 756−765. (2) Huang, J.; Tao, A. R.; Connor, S.; He, R.; Yang, P. A general method for assembling single colloidal particle lines. Nano Lett. 2006, 6, 524−529. (3) Layani, M.; Gruchko, M.; Milo, O.; Balberg, I.; Azulay, D.; Magdassi, S. Transparent conductive coatings by printing coffee ring arrays obtained at room temperature. ACS Nano 2009, 3, 3537−3542. (4) Bodiguel, H.; Doumenc, F.; Guerrier, B. Stick- slip patterning at low capillary numbers for an evaporating colloidal suspension. Langmuir 2010, 26, 10758−10763. (5) Blossey, R.; Bosio, A. Contact line deposits on cDNA microarrays: a “twin-spot effect. Langmuir 2002, 18, 2952−2954. (6) Heim, T.; Preuss, S.; Gerstmayer, B.; Bosio, A.; Blossey, R. Deposition from a drop: morphologies of unspecifically bound DNA. J. Phys.: Condens. Matter 2005, 17, S703−S716. (7) Zhang, D.; Xie, Y.; Mrozek, M. F.; Ortiz, C.; Davisson, V. J.; BenAmotz, D. Raman detection of proteomic analytes. Anal. Chem. 2003, 75, 5703−5709. (8) Kim, Y.; Hurst, G. B.; Doktycz, M. J.; Buchanan, M. V. Improving spot homogeneity by using polymer substrates in matrix-assisted laser desorption/ionization mass spectrometry of oligonucleotides. Anal. Chem. 2001, 73, 2617−2624. (9) Norris, D. J.; Arlinghaus, E. G.; Meng, L.; Heiny, R.; Scriven, L. Opaline Photonic Crystals: How Does Self-Assembly Work? Adv. Mater. (Weinheim, Ger.) 2004, 16, 1393−1399. (10) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. R.; Witten, T. A. Capillary flow as the cause of ring stains from dried liquid drops. Nature 1997, 389, 827−829. (11) Adachi, E.; Dimitrov, A. S.; Nagayama, K. Stripe patterns formed on a glass surface during droplet evaporation. Langmuir 1995, 11, 1057−1060. (12) Hu, H.; Larson, R. G. Evaporation of a sessile droplet on a substrate. J. Phys. Chem. B 2002, 106, 1334−1344. (13) Hu, H.; Larson, R. G. Analysis of the microfluid flow in an evaporating sessile droplet. Langmuir 2005, 21, 3963−3971. 5858
dx.doi.org/10.1021/la500805d | Langmuir 2014, 30, 5849−5858
