Competitive evaporation of multiple sessile droplets

An asymptotic model is derived for the competitive diffusion-limited evaporation of multiple thin sessile droplets under the assumption that the droplets are well separated. Exact solutions of the model are obtained for a pair of and for a polygonal array of identical droplets, and the model is found to perform well even outside its formal range of validity, up to and including the limit of touching droplets. The shielding effect of droplets on each other is demonstrated, and the model is used to investigate the effect of this shielding on droplet evolutions and lifetimes, as well as on the coffee-ring effect. The theoretical predictions of the model are found to be in good agreement with recent experimental results for seven relatively closely spaced droplets, suggesting that the model could be a useful tool for studying a wide range of other droplet configurations.


Introduction
Motivated in part by its numerous applications in industry (such as, for example, in printing, nanofabrication, spray cooling and biochemical assays), the evaporation of sessile droplets has been extensively studied experimentally, numerically and analytically (Routh 2013;Larson 2014;Brutin & Starov 2018;Giorgiutti-Dauphiné & Pauchard 2018).In particular, in recent years considerable attention has been paid to aspects of evaporation such as the coffee-ring effect (Deegan et al. 1997) and droplet lifetimes (Stauber et al. 2014(Stauber et al. , 2015)).
The vast majority of the previous work has concerned a single droplet, but in practice most droplets do not occur in isolation, and so interactions between droplets are of great practical and scientific interest.In the present contribution we investigate the competitive diffusion-limited evaporation of multiple thin sessile droplets.In this case the evaporative flux is determined by solving a challenging three-dimensional mixed boundary value problem for Laplace's equation.The problem is also of interest as it has analogues in a wide range of other physical situations, including the dissolution of immersed nanobubbles and nanodroplets (Lohse & Zhang 2015; 884 A45-2 A. W. Wray, B. R. Duffy and S. K. Wilson Qian, Arends & Zhang 2019), the electrical contact resistance between elastic bodies (Argatov 2011), the electric field around electrodes at constant potential (Jackson 1999, § 3.13), the thermal field around plates held at constant temperature (Lee & Chien 1994) and flow through a porous membrane (Fabrikant 1985).
As a result, there is significant interest in configurations in which this mixed boundary value problem can be solved exactly.Unfortunately, only a few exact solutions are known, notably those for a flat disk (Weber 1873) and a spherical cap (Lebedev 1965).For studies of more complicated geometries, such as droplets with elliptical or triangular contact lines, authors have had to resort to numerical methods (Sáenz et al. 2015(Sáenz et al. , 2017)).
Previous studies of the evaporation of multiple sessile droplets have used a variety of experimental, numerical and analytical approaches (Lacasta et al. 1998;Schäfle et al. 1999;Kokalj et al. 2010;Sokuler et al. 2010;Carrier et al. 2016;Castanet et al. 2016;Shaikeea et al. 2016;Shaikeea & Basu 2016a,b;Hatte et al. 2019;Khilifi et al. 2019;Schofield et al. 2020).The critical difference in comparison with what is observed with a single droplet is the occurrence of 'shielding', in which the presence of other droplets causes a droplet to evaporate more slowly than it would in isolation.Physically this is because each droplet feels an increased vapour concentration in the surrounding atmosphere due to the presence of the other droplets, reducing its rate of evaporation.Also relevant here is the extensive body of recent work on the dissolution of immersed nanobubbles and nanodroplets (Lohse & Zhang 2015;Qian et al. 2019), in particular, the work on the interactions between multiple nanobubbles or nanodroplets (Weijs, Seddon & Lohse 2012;Weijs & Lohse 2013;Dollet & Lohse 2016;Laghezza et al. 2016;Michelin, Guérin & Lauga 2018;Zhu et al. 2018;Xie & Harting 2019) which also exhibit an analogous shielding effect due to an increased concentration of nanobubble gas or nanodroplet liquid in the immersing fluid.
Two of the relatively few analytical studies relevant to the present problem are those by Dollet & Lohse (2016) and Fabrikant (1985).In their work on the dissolution of surface nanobubbles, Dollet & Lohse (2016), adopting an approach similar to that used by Argatov (2011) in a different physical context, considered the dissolution of two well-separated nanobubbles which are treated as dissolving independently at leading order, with the shielding effect entering as a spatially uniform first-order correction to the local volume fluxes from the nanobubbles.On the other hand, in his work on potential flow through a porous membrane, Fabrikant (1985) used a series of integral transformations to reduce the problem to a system of integral equations from which he obtained approximations to the integral volume fluxes through the pores.It is this latter approach that we shall follow in the present work.While it is restricted to thin droplets, this approach allows us to obtain explicit expressions for the spatially varying local fluxes from any given configuration of droplets.In § 2 we define precisely the mathematical problem that we are solving.In § 3 we put the approximate approach of Fabrikant on an asymptotic basis, before extending the procedure to obtain explicit expressions for the local fluxes from the droplets.In § 4 we calculate the evolutions, and hence the lifetimes, of multiple droplets in two particular configurations.In § 5 we explain qualitatively what effect the presence of another droplet has on the coffee-ring effect, while in § 6 we compare the theoretical predictions of the present asymptotic model with recent experimental results.Finally, we summarise our conclusions in § 7.
Competitive evaporation of multiple sessile droplets 884 A45-3 are thin, so that to a leading-order approximation they appear to be flat (Dunn et al. 2008;Schofield et al. 2018), and that the evaporation is in the diffusion-limited regime (Sultan, Boudaoud & Ben Amar 2005), so that the local mass flux from each of the droplets, J k for k = 1, 2, . . ., N, may be calculated from the vapour concentration c (Popov 2005).Far from the droplets c approaches the constant far-field value c ∞ , while on the surfaces of the droplets, denoted by S k for k = 1, 2, . . ., N, it takes the constant saturation value c sat .(Note that Dunn et al. (2008), Dunn et al. (2009) and Schofield et al. (2018) treated the case of an isolated droplet when c sat depends on temperature.)We non-dimensionalise the system using a characteristic radius of the droplets, a ref , according to where D denotes the diffusion coefficient of vapour in the atmosphere and hats denote dimensionless quantities.Immediately dropping the hats, c satisfies Laplace's equation (2.4) Introducing local polar coordinates (ρ, φ) centred at (x k , y k ), the flux from S k is denoted by J k (ρ, φ) and is given by Note that in order to determine the evolution of the droplets it is only the J k that we require: for our present purpose the solution for c is of secondary importance.
Using an ingenious series of integral transformations, Fabrikant (1985) showed that for the problem (2.2)-(2.4),J k (denoted by −σ k in his notation) satisfies where J 0 = J 0 (ρ) given by is the flux from the kth droplet in isolation.
Note that, up to this point, Fabrikant's results are exact.

Explicit asymptotic expressions for the fluxes J k
Our aim is to derive explicit expressions for the fluxes J k given by (2.6) in the asymptotic limit of well-separated droplets.Our strategy is as follows : Fabrikant (1985) demonstrated a way of using (2.9) to obtain approximations to the integral fluxes F k .We put this approach on an asymptotic basis, before using (2.6) to derive explicit expressions for the fluxes J k .These expressions are, in general, spatially varying, and will be shown to be in excellent agreement with the exact solution of the problem for a pair of identical droplets.
Let r k,n = (x n − x k ) 2 + (y n − y k ) 2 denote the distance between the centres of the kth and nth droplets (which must, by definition, satisfy r k,n a k + a n ), and let ψ k,n defined by tan ψ k,n = (y n − y k )/(x n − x k ) denote the angle that the line joining the centres makes with the θ = 0 direction, as shown in figure 1.We assume that the radius of each droplet a k is small relative to the distance between the centre of that droplet and any other, i.e.
In the integrands of (2.6) and (2.9) we write ρ = r k,n + R n , where |R n | a n r k,n , and φ = ψ k,n + Ψ n .Note that 1 a n /r k,n sin Ψ n ∼ Ψ n , so that both R n /r k,n and Ψ n are small across the range of integration.
3.1.Equations for the integral fluxes F k At leading order in , equation (2.9) becomes This is exactly the asymptotically rigorous derivation of Fabrikant's equation ( 15) with his 'central estimate' (denoted by b k,n = r k,n in his notation).Equation (3.2) constitutes a set of N linear algebraic equations for the integral fluxes F k which may, in principle, be solved exactly for any given configuration of droplets.
3.2.Explicit expressions for the fluxes J k (ρ, φ) At leading order in , equation (2.6) becomes With the integral fluxes F k obtained from (3.2), equation (3.3) gives the fluxes J k explicitly.We now analyse two particular configurations, namely for a pair of identical droplets in § 3.3, and for a polygonal array of N identical droplets in § 3.4.

3.3.
Solution for a pair of identical droplets Consider a pair of identical droplets with the same radii a 1 = a 2 = a and centres a distance b ( 2a) apart located at (−b/2, 0) and (b/2, 0).From (3.2) the two droplets have the same integral flux Hence (3.3) gives the flux from the droplet centred at (−b/2, 0), denoted by J 1 (ρ, φ), to be The corresponding expression for the flux from the droplet centred at (+b/2, 0), denoted by J 2 (ρ, φ), is obtained by replacing φ with φ + π in (3.5).Note that the term ρ 2 + b 2 − 2ρb cos φ appearing in (3.5) is the square of the distance from (ρ, φ) to the centre of the other droplet, which is always strictly positive.The second term in the square brackets in (3.5) represents the shielding effect of the other droplet and is inversely proportional to the distance between the centres of the droplets in the limit of well-separated droplets, i.e. it is O(1/b) in the limit b → ∞.However, note that this term is always finite, and thus the fluxes J 1 and J 2 have the same square-root singularity as J 0 (i.e. the same singularity as an isolated droplet) at the contact lines of the droplets.

Solution for a polygonal array of N identical droplets
Consider N identical droplets with the same radii a k = a and centres (x k , y k ) for k = 1, 2, . . ., N located at the vertices of a regular N-gon.The circumradius (i.e. the distance between the centre of the polygon and the centres of the droplets) is denoted by b/2, and so the distance between the centres of the kth and nth droplets is (3.9) Numbering the droplets anticlockwise with the Nth one at (−b/2, 0), we find and (3.3) gives the flux from the Nth droplet, J N = J N (ρ, φ), to be Once again, the second term in the square brackets represents the shielding effect of all the other droplets on the Nth droplet, and the flux J N has the same square-root singularity as J 0 at the contact line of the Nth droplet.

Model validation
As a check on the validity of the present asymptotic model, we consider the pair of identical droplets with radius a a distance b apart described in § 3.3 when a = 1.
First we compare the leading-order asymptotic approximations for the integral flux F given by (3.4) and the spatially varying local fluxes J 1 given by (3.5) and J 2 with the corresponding exact solutions, which we denote by F K , J 1K and J 2K , respectively.We computed the exact solutions using the formulation due to Kobayashi (1939) (see also Sneddon (1966, pp. 259-264); note a small but significant typographical error in this reference: the first Bessel function appearing in equation (8.4.9) should be of order h ± m, not h + m), increasing the number of modes in the truncated series expansion until convergence to at least five decimal places was achieved, and thereby recovering the impressively accurate numerical values of F K obtained by Kobayashi (1939, table 1) himself 80 years ago.Figure 2(a) shows the exact solution F K (solid line) and the asymptotic approximation F given by (3.4) (dashed line) as functions of b.In particular, figure 2(a) shows that, as a consequence of the shielding effect, F K increases monotonically with b from the value F K 3.011 when b = 2 to F K → 4 − in the limit b → ∞.For all values of b the asymptotic solution F is in excellent agreement with the exact solution F K , indeed so good that the two are almost indistinguishable in figure 2 and has a magnitude of approximately 0.367 %, in agreement with the corresponding results of Fabrikant (1985, table 2).As figure 2 shows, despite the asymptotic model being formally valid only for well-separated droplets, the agreement with the exact solution is excellent up to and including the limit of touching droplets.
Secondly, for the particular case b = 2.2, corresponding to two relatively closely spaced droplets, figure 3 shows a comparison of the contours of the exact solutions for the fluxes J 1K and J 2K (figure 3a) with the corresponding contours of the asymptotic solutions for J 1 and J 2 (figure 3b), and comparisons of J 1K and J 2K (solid lines) with J 1 and J 2 (dashed lines) along the cross-sections x = ±1.1 (−1 y 1) and y = 0 (−2.1 x −0.1 and 0.1 x 2.1) through the centres of the droplets (figure 3c).As figure 3 shows, the asymptotic model accurately captures the shielding effect: evaporation is depressed on the side of the droplet that is closer to the other droplet, resulting in asymmetric profiles of J 1 and J 2 along the line y = 0 (but symmetric profiles along the lines x = ±b/2 = ±1.1).Again, despite the asymptotic model being formally valid only for well-separated droplets, the agreement with the exact solution in this case of two relatively closely spaced droplets is excellent.In particular, in this case the prediction for the minimum value of the flux is in error by only 1.65 %.
Comparison of (a) the contours of the exact solutions for the fluxes J 1K and J 2K with (b) the corresponding contours of the asymptotic solutions for J 1 and J 2 , and (c) comparisons of J 1K and J 2K (solid lines) with J 1 and J 2 (dashed lines) along the crosssections x = ±1.1 (−1 y 1) and y = 0 (−2.1 x −0.1 and 0.1 x 2.1) through the centres of the droplets for a pair of identical droplets when a = 1 and b = 2.2.
Finally, we note that for further validation we also developed a finite-differencebased code to solve the original governing equations (2.2)-(2.4)numerically.However, despite using O(10 5 ) points we found that convergence to the desired number of decimal places could not be attained, due to the singularities at the contact lines.Indeed, in the present case of two identical droplets we found that the asymptotic model is generally accurate to more decimal places than the solution obtained via our finite-difference-based code.

Evolutions and lifetimes of droplets
Encouraged by the excellent agreement between the predictions of the asymptotic model and exact solution for a pair of identical droplets described in § 3.5, we now use the solutions of (3.2) for the integral fluxes F k for k = 1, . . ., N to determine the evolutions, and hence the lifetimes, of any configuration of droplets.
Denoting the (small) contact angle of the kth droplet by θ k , its free surface profile by z = h k , and its volume by V k , we scale θ k , h k , V k and time t according to where θ ref is a characteristic contact angle, ρ fluid is the constant density of the fluid and the hats are again immediately dropped.For droplets small enough that surface-tension effects dominate gravitational effects, we have where we now allow both their radii a k = a k (t) and their contact angles θ k = θ k (t), and hence their free surface profiles h k = h k (ρ, t) and their volumes V k = V k (t), to be functions of time.For simplicity, we assume that the distances between the centres of the droplets remain fixed, but other regimes are, in principle, treatable within the same framework.The rate of decrease of the volume of the kth droplet is given by the integral flux F k , i.e. −dV k /dt = F k , and hence 4.1.Evolution and lifetime of a pair of identical droplets For the pair of identical droplets each with radius a = a(t), contact angle θ = θ (t) and volume V = V(t) a distance b apart analysed in § 3.3, substituting the solution for F given by (3.4) into (4.3)gives We consider three modes of evaporation (Picknett & Bexon 1977;Stauber et al. 2014Stauber et al. , 2015)), namely the extreme constant radius (CR) mode for which a ≡ ā and θ = θ (t) subject to the initial condition θ(0) = θ, the extreme constant angle (CA) mode for which θ ≡ θ and a = a(t) subject to the initial condition a(0) = ā and the intermediate stick-slide (SS) mode for which a ≡ ā and θ = θ (t) until θ reaches its critical (receding) value θ = θ (0 θ θ ) at the de-pinning time t = t , and thereafter θ ≡ θ and a = a(t) subject to the initial condition θ (0) = θ, where ā and θ are constants (both of which may, without loss of generality, be set equal to unity, but are retained in what follows for clarity).Note that the SS mode reduces to the CR mode in the limiting case θ = 0 and to the CA mode in the limiting case θ = θ.which readily yields the exact solutions for the evolutions of θ and V, namely In particular, the lifetime of the droplets in the CR mode, denoted by t CR , is determined by setting V = 0 in (4.6) to obtain where t CR∞ = πā 2 θ/16 is the lifetime of the same droplets evaporating in isolation, which is recovered when the distance between the droplets becomes large, i.e.
which yields the implicit solutions for the evolution of a and hence of V = πa 3 θ /4, namely where the function I = I(a) is defined by and hence is given by In particular, the lifetime of the droplets in the CA mode, denoted by t CA , is given by i.e.
where t CA∞ = 3πā 2 θ/32 is the lifetime of the same droplets evaporating in isolation, i.e. t CA → t CA + ∞ as b → ∞.Note that in the opposite limit of initially touching droplets, b → 2ā + , then t CA → (2π In the SS mode, the droplets initially evaporate in a CR phase satisfying (4.5) according to (4.6) until θ = θ , which occurs at t = t , where t is given by Thereafter, the droplets evaporate in a CA phase satisfying (4.8) with θ replaced by θ according to where the function I = I(a) is again given by (4.11).In particular, the lifetime of the droplets in the SS mode, denoted by t SS , is given by i.e.
where t SS∞ = πā 2 (2 θ + θ * )/32 is the lifetime of the same droplets evaporating in isolation, i.e. t SS → t SS + ∞ as b → ∞.Note that in the opposite limit of initially touching droplets, b → 2ā + , then Figure 4(a) shows t CR , t CA and t SS given by (4.7), (4.13) and (4.17), respectively, plotted as functions of b ( 2) when ā = 1 and θ = 1.In particular, figure 4(a) shows that, due to the shielding effect, t CR , t CA and t SS are all monotonically decreasing functions of b for fixed θ (i.e. the lifetime of a pair of droplets is always longer than that of the same droplets in isolation), but monotonically increasing functions of θ for fixed b (i.e. the lifetime of the droplets in the SS mode always lies between that of the same droplets in the extreme modes, t CR t SS t CA ).Furthermore, as figure 4(a) also shows, while the relative impact of the shielding effect on the lifetime of the droplets decreases slightly with θ , even touching droplets evaporating in the CR mode (i.e.t CR → πā 2 θ /12 − in the limit b → 2ā + ) still have a shorter lifetime than isolated droplets evaporating in the CA mode (i.e.t CA → 3πā 2 θ /32 + in the limit b → ∞).

Evolution and lifetime of a polygonal array of N identical droplets
Following the same approach as in § 4.1, we can obtain the corresponding expressions for the evolution and lifetimes of the polygonal array of N identical droplets analysed in § 3. and hence given by in the appropriate places.In particular, the lifetime of the droplets in the SS mode (from which the lifetime of the droplets in the extreme modes can readily be deduced) Competitive evaporation of multiple sessile droplets 884 A45-13 is given by (4.22) Figure 4(b) shows t CR , t CA and t SS for a polygonal array of N = 5 droplets plotted as functions of b ( 2/ sin(π/5) 3.403) when ā = 1 and θ = 1.In particular, figure 4(b), which is typical of the corresponding results for all values of N 3, shows that the behaviour of a polygonal array of N identical droplets is qualitatively similar to that of a pair of identical droplets (i.e. to that in the special case N = 2 shown in figure 4a).However, as figure 4(b) illustrates, touching droplets evaporating in the CR mode having a shorter lifetime than isolated droplets evaporating in the CA mode is peculiar to the special case N = 2.

Radially integrated evaporative flux R k
We now examine how the radially integrated evaporative flux from the kth droplet, denoted by R k = R k (φ) and given by varies with the local polar angle φ.As discussed below, this quantity is of particular interest because of its connection to the coffee-ring effect.
For the pair of identical droplets analysed in § 3.3, substituting the solution for J 1 given by (3.5) into (5.1)gives the radially integrated flux in the droplet centred at (−b/2, 0), denoted by R 1 (φ), to be where k = a/b, R 0 = 2a/π is the radially integrated flux from an isolated droplet, and F is given by (3.4).The corresponding expression for the radially integrated flux in the other droplet centred at (b/2, 0), denoted by R 2 (φ), is obtained by replacing φ with φ + π in (5.2).In particular, in the limit k = a/b → 0 we have In order to assess the accuracy of the expressions (5.2) and (5.3), figure 5 shows a comparison of them with the exact solution for the radially integrated flux R 1 computed using the formulation due to Kobayashi (1939) as described in § 3.5 for three different values of b.In particular, figure 5 shows that the asymptotic solution (5.2) is quite accurate even for relatively large values of a/b, but, as expected, the expansion of the asymptotic solution (5.3) is a useful approximation only for sufficiently small values of a/b.Sáenz et al. (2017) demonstrated that, for a droplet containing nanoparticles evaporating in the CR mode, the distribution of the final residue is strongly related to the radially integrated fluid flux within the droplet, which is in turn closely related to the radially integrated evaporative flux R k .In particular, the radial directions with the greatest values of R k have the greatest fluid flux within the droplet, giving the greatest concentration of residue at the contact line in that direction.Of course, because of other effects, such as particle diffusion and azimuthal flow, R k does not give a perfect characterisation of the final residue distribution.However, Sáenz et al. (2017) observed that even in the highly asymmetric case of a triangular droplet the flow is predominantly radial, and in many situations R k will indeed provide a good indication of what final residue will occur.The solutions for R 1 for a pair of identical droplets shown in figure 5 reveal that, as expected, the azimuthal variation of R 1 is most pronounced when the droplets are closest together (i.e. when the shielding effect is most significant), and so we would anticipate that this situation would give rise to the most non-homogeneous coffee rings, with least residue accumulating where the contact lines are closest together and most residue accumulating where they are furthest apart.

FIGURE 2 .
FIGURE 2. Plots of (a) the exact solution for the integral flux F K (solid line) and the (almost indistinguishable) asymptotic approximation F given by (3.4) (dashed line), and (b) the associated relative error e rel defined by (3.12) as functions of b for a pair of identical droplets a distance b ( 2) apart when a = 1. https://doi.org/10.1017/jfm.2019.

FIGURE 5 .
FIGURE 5. Comparison between the exact solution (solid lines), the asymptotic solution given by (5.2) (dashed lines) and the expansion of the asymptotic solution given by (5.3) (dotted lines, indistinguishable from the solid line for b = 10) for the radially integrated evaporative flux R 1 from the droplet centred at (−b/2, 0) for the pair of identical droplets analysed in § 3.3 plotted as functions of the scaled local polar angle φ/π when a = 1 for b = 2, 6 and 10.The direction of the arrow corresponds to increasing values of b.The horizontal line R 1 = R 0 = 2a/π 0.637 corresponds to the limiting value of R 1 in the limit b → ∞.

6.
Comparison with the experimental results ofKhilifi et al. (2019) The ultimate test for the value of the present asymptotic model is, of course, how well its theoretical predictions compare with experimental results.To this end, we compare the theoretical predictions of the present model with the recent experimental results ofKhilifi et al. (2019) for seven relatively closely spaced droplets.Khilifi et al. (2019) studied the evaporation of sessile droplets of water at an initial temperature of 20 • C on a sapphire substrate into an enclosed chamber with 50 % humidity at an ambient temperature of 21 • C. Initial experiments using an isolated droplet revealed that the droplets evaporate in an SS mode with an initial contact angle https://doi.org/10.1017/jfm.2019.
Fabrikant then related the integral of the flux over S k , denoted by F k (also denoted by F k in his notation) and given by