3.1 Evolution of the droplet radius

Obtaining an evolution equation for the droplet radius, $r(t)$ , is carried out in the same fashion as in related studies involving contact line dynamics without evaporation (see e.g. Hocking Reference Hocking1983; Savva & Kalliadasis Reference Savva and Kalliadasis2009; Vellingiri, Savva & Kalliadasis Reference Vellingiri, Savva and Kalliadasis2011; Savva & Kalliadasis Reference Savva and Kalliadasis2012, Reference Savva and Kalliadasis2013), whereby we treat the bulk of the liquid drop separately from the region in the vicinity of the contact line. In the central region of the droplet, the outer region, slip and evaporation-induced free-surface deformations are neglected, but they become important near the contact line. Gravitational effects are important at the onset in the outer region only, since they diminish as the droplet becomes smaller due to evaporation.

3.1.1 Outer region

The analysis of the outer region is nearly identical to the analysis undertaken by Hocking (Reference Hocking1983), the main difference being the time variation of the droplet volume. Here we reiterate the main results, albeit with a slightly different notation, referring the interested reader to the original work of Hocking for further details. The rate of spreading ${\dot{r}}$ is assumed to be small, so that we may employ the quasi-static approximation, introducing the expansion

(3.1) $$\begin{eqnarray}h_{\mathit{out}}\sim h_{0}(x,r,v)+{\dot{r}}h_{1}(x,r,v)+\cdots ,\end{eqnarray}$$

i.e. time enters the problem through $v$ , $r$ and ${\dot{r}}$ . Anticipating the scalings for ${\dot{r}}$ and $\dot{v}$ determined later on, we find that, as $\unicode[STIX]{x1D706}\rightarrow 0$ , ${\dot{r}}=O(1/|\!\ln \unicode[STIX]{x1D706}|)$ (as in the non-volatile case) and $\dot{v}=O(\unicode[STIX]{x1D706}|\!\ln \unicode[STIX]{x1D706}|)$ (as will follow from the analysis of (2.20)). Hence, although ${\dot{r}}$ is small, it is nonetheless generally much greater than $\unicode[STIX]{x1D706}$ and $\dot{v}$ and is accordingly treated as such. Thus, we neglect the evaporation and slip effects in the outer region, so that $h_{\mathit{out}}$ satisfies at $O(\unicode[STIX]{x1D706}^{0})$

(3.2) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}h_{\mathit{out}}+\frac{1}{x}\unicode[STIX]{x2202}_{x}\left\{h_{\mathit{out}}^{3}x\unicode[STIX]{x2202}_{x}\left(\unicode[STIX]{x2202}_{x}^{2}h_{\mathit{out}}+\frac{1}{x}\unicode[STIX]{x2202}_{x}h_{\mathit{out}}-Bh_{\mathit{out}}\right)\right\}=0.\end{eqnarray}$$

Collecting powers of ${\dot{r}}$ in (3.2) gives

(3.3) $$\begin{eqnarray}\unicode[STIX]{x2202}_{x}\left(\unicode[STIX]{x2202}_{x}^{2}h_{0}+\frac{1}{x}\unicode[STIX]{x2202}_{x}h_{0}-Bh_{0}\right)=0\end{eqnarray}$$

to $O({\dot{r}}^{0})$ . This is a third-order differential equation, and is solved subject to (2.18a ), (2.18d ) and (2.19). The solution to (3.3) is given in terms of modified Bessel functions as

(3.4) $$\begin{eqnarray}h_{0}=\frac{\unicode[STIX]{x1D703}}{\sqrt{B}}\frac{I_{0}(c)-I_{0}(x\sqrt{B})}{I_{1}(c)},\end{eqnarray}$$

where $c=r\sqrt{B}$ and $\unicode[STIX]{x1D703}$ is the apparent contact angle, defined as

(3.5) $$\begin{eqnarray}\unicode[STIX]{x1D703}=\frac{2v\sqrt{B}I_{1}(c)}{r^{2}I_{2}(c)}.\end{eqnarray}$$

Note that $\unicode[STIX]{x1D703}$ corresponds to a rescaled angle, which needs to be multiplied by $\unicode[STIX]{x1D717}_{s}$ to obtain the true one. When $B=0$ , equations (3.4) and (3.5) reduce to

(3.6a,b ) $$\begin{eqnarray}h_{0}=\frac{\unicode[STIX]{x1D703}r}{2}\left(1-\frac{x^{2}}{r^{2}}\right)\quad \text{and}\quad \unicode[STIX]{x1D703}=\frac{8v}{r^{3}},\end{eqnarray}$$

respectively. From the $O({\dot{r}})$ terms in (3.2) we find that $h_{1}$ must satisfy

(3.7) $$\begin{eqnarray}\unicode[STIX]{x2202}_{r}h_{0}+\frac{1}{x}\unicode[STIX]{x2202}_{x}\left\{h_{0}^{3}x\unicode[STIX]{x2202}_{x}\left(\unicode[STIX]{x2202}_{x}^{2}h_{1}+\frac{1}{x}\unicode[STIX]{x2202}_{x}h_{1}-Bh_{1}\right)\right\}=0,\end{eqnarray}$$

obtained by taking $\unicode[STIX]{x2202}_{t}h_{0}={\dot{r}}\unicode[STIX]{x2202}_{r}h_{0}$ . Multiplying both sides by $x$ and integrating yields

(3.8) $$\begin{eqnarray}\unicode[STIX]{x2202}_{x}\left(\unicode[STIX]{x2202}_{x}^{2}h_{1}+\frac{1}{x}\unicode[STIX]{x2202}_{x}h_{1}-Bh_{1}\right)=\frac{r^{3}I_{1}(c)I_{2}(c)[xI_{1}(c)-rI_{1}(x\sqrt{B})]}{4v^{2}[I_{0}(c)-I_{0}(x\sqrt{B})]^{3}}.\end{eqnarray}$$

This is also a third-order differential equation solved subject to homogeneous boundary conditions, namely $h_{1}=0$ at $x=r(t)$ , $\unicode[STIX]{x2202}_{x}h_{1}=0$ at $x=0$ and $\int _{0}^{r}xh_{1}\,\text{d}x=0$ . Here we are merely interested in obtaining the behaviour of the slope of $h_{1}$ as the contact line is approached. Following Hocking (Reference Hocking1983), the slope has the asymptotic form

(3.9) $$\begin{eqnarray}-\unicode[STIX]{x2202}_{x}h_{1}\sim \frac{1}{\unicode[STIX]{x1D703}^{2}}\ln \left[\frac{\text{e}^{2}}{\unicode[STIX]{x1D6FF}}\left(1-\frac{x}{r}\right)\right]\end{eqnarray}$$

as $x\rightarrow r$ , where $\unicode[STIX]{x1D6FF}$ is determined from

(3.10) $$\begin{eqnarray}\ln \unicode[STIX]{x1D6FF}=\int _{0}^{1}x\left\{\frac{cI_{1}^{3}(c)[xI_{1}(c)-I_{1}(xc)]^{2}}{I_{2}^{2}(c)[I_{0}(c)-I_{0}(xc)]^{3}}-\frac{1}{1-x}\right\}\text{d}x,\end{eqnarray}$$

which changes with $r$ and needs to be evaluated numerically. In (3.9), $O(x-r)$ terms are neglected and, just as in many instances throughout the manuscript, the constant terms are absorbed into the logarithmically diverging term. It is important to emphasise that although the above analysis concerns sessile droplets (with $B>0$ ), the extension to pendant droplets (with $B<0$ ) is trivial and one needs to replace the modified Bessel functions with Bessel functions, take $c=r\sqrt{-B}$ throughout and change the sign of the first term in braces in (3.10), keeping also in mind that sufficiently large pendant droplets may possibly detach from the substrate. Thus, since our model does not account for the detachment dynamics of pendant droplets, we restrict our treatment to pendant droplets for which $c$ is always less than the limiting value of $c$ , the first positive root of the Bessel function $J_{1}$ . When $c$ attains this value, the denominator of (3.4) vanishes when modified appropriately for pendant droplets. Figure 2 shows a plot of $\ln \unicode[STIX]{x1D6FF}$ as a function of $r\sqrt{|B|}$ for both pendant and sessile droplets. As expected, we see that as the droplet evaporates so that $r\rightarrow 0$ , the effect of gravity on the value of $\unicode[STIX]{x1D6FF}$ diminishes and $\unicode[STIX]{x1D6FF}\rightarrow 1/2$ , the zero Bond number limit. For pendant droplets, there is a vertical asymptote near 3.83, the first positive root of  $J_{1}$ .

Figure 2. Dependence of $\ln \unicode[STIX]{x1D6FF}$ on $r\sqrt{|B|}$ for a sessile ( $B>0$ ; solid curve) and a pendant ( $B<0$ ; dashed curve) droplet. As the droplet evaporates, $r\rightarrow 0$ , the curves approach asymptotically the value of $\ln \unicode[STIX]{x1D6FF}=-\text{ln}\,2$ , corresponding to the $B=0$ case.

Hence, the asymptotics of the slope of the free surface as $x\rightarrow r$ can be readily obtained by combining (3.4) and (3.9):

(3.11) $$\begin{eqnarray}-\unicode[STIX]{x2202}_{x}h\sim \unicode[STIX]{x1D703}+\frac{{\dot{r}}}{\unicode[STIX]{x1D703}^{2}}\ln \left[\frac{\text{e}^{2}}{\unicode[STIX]{x1D6FF}}\left(1-\frac{x}{r}\right)\right].\end{eqnarray}$$

This behaviour is to be matched, within an appropriate overlap region, with the corresponding behaviour of the inner region.

3.1.2 Inner region

Whilst considering the solution in the outer region, we did not use the contact angle condition, equation (2.18b ), which can only be applied in the inner region, where slip and evaporation effects become important. From the inner region near the contact line, we anticipate a behaviour that allows the matching with the outer-region dynamics as we move away from the contact line.

To investigate the inner-region dynamics, introduce the change of variables

(3.12a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D702}=\frac{r-x}{\unicode[STIX]{x1D706}}\quad \text{and}\quad \unicode[STIX]{x1D719}=\frac{h}{\unicode[STIX]{x1D706}},\end{eqnarray}$$

which allows us to write (2.14) as

(3.13) $$\begin{eqnarray}{\dot{r}}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D719}+\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}[\unicode[STIX]{x1D719}^{2}(\unicode[STIX]{x1D719}+1)\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}^{3}\unicode[STIX]{x1D719}]=-\frac{E}{\unicode[STIX]{x1D719}+K}\quad \text{for}~\unicode[STIX]{x1D702}>0,\end{eqnarray}$$

where we have retained $O(\unicode[STIX]{x1D706}^{0})$ terms only, neglecting $O(\unicode[STIX]{x1D706}^{2}B)$ terms for moderate gravitational effects. The boundary conditions to be applied are $\unicode[STIX]{x1D719}=0$ and $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D719}=1$ at $\unicode[STIX]{x1D702}=0$ and, in order to match with (3.11), $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D719}$ must be no more than logarithmically large as $\unicode[STIX]{x1D702}\rightarrow \infty$ . Note that (3.13) arises in the two-dimensional geometry as well (Savva et al. Reference Savva, Rednikov, Colinet, Karayiannis, König and Balabani2014), which is not surprising given that the contact line appears to be a straight line in its immediate vicinity when $r\gg \unicode[STIX]{x1D706}$ .

Generally, the evaporative flux terms are comparable to the capillary terms and need to be retained in the leading-order equations (see also appendix A, where the case $E\ll 1$ is explored). To proceed, just like the outer region, consider a small- ${\dot{r}}$ expansion in (3.13) of the form

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D719}\sim \unicode[STIX]{x1D719}_{0}+{\dot{r}}\unicode[STIX]{x1D719}_{1}+\cdots .\end{eqnarray}$$

At leading order in ${\dot{r}}$ , (3.13) gives

(3.15) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}[\unicode[STIX]{x1D719}_{0}^{2}(\unicode[STIX]{x1D719}_{0}+1)\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}^{3}\unicode[STIX]{x1D719}_{0}]=-\frac{E}{\unicode[STIX]{x1D719}_{0}+K},\end{eqnarray}$$

which is solved for $\unicode[STIX]{x1D702}>0$ subject to

(3.16a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{0}=0\quad \text{and}\quad \unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D719}_{0}=1\end{eqnarray}$$

at $\unicode[STIX]{x1D702}=0$ and requiring that $\unicode[STIX]{x1D719}_{0}$ behaves linearly at infinity, namely

(3.16c ) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{0}\sim \unicode[STIX]{x1D703}_{m}\unicode[STIX]{x1D702}\quad \text{as}~\unicode[STIX]{x1D702}\rightarrow \infty ,\end{eqnarray}$$

where $\unicode[STIX]{x1D703}_{m}$ is the (rescaled with $\unicode[STIX]{x1D717}_{s}$ ) macroscopic Young’s angle modified by evaporation. As it turns out, $\unicode[STIX]{x1D703}_{m}$ is a degree of freedom arising from the far-field expansion of the problem and needs to be determined for given values of $E$ and $K$ . Thus, solving (3.15) with (3.16) numerically allows us to extract the value of $\unicode[STIX]{x1D703}_{m}$ (see appendix B.1 for details).

Figure 3. (a) Plots of $\unicode[STIX]{x1D703}_{m}$ as a function of $E$ for various values of $K$ . (b) Plots of $\unicode[STIX]{x1D703}_{m}$ as a function of $K$ for various values of $E$ . The solid curves correspond to the values of $\unicode[STIX]{x1D703}_{m}$ obtained by solving (3.15)–(3.16); the dashed curves correspond to the asymptotic result for $E\ll 1$ , equation (A 5) with (A 8); the black dotted curves show the asymptotic behaviour $\unicode[STIX]{x1D703}_{m}\sim \unicode[STIX]{x1D701}(K)E^{1/4}$ as $E\rightarrow \infty$ ; the grey dotted curves show the result of Hocking’s asymptotics, equation (3.17). The asymptotic results are only shown for the cases where their applicability can reasonably be expected.

Figure 3 summarises some representative computations performed for $\unicode[STIX]{x1D703}_{m}$ and its asymptotics. Figure 3(a) shows the dependence of $\unicode[STIX]{x1D703}_{m}$ with $E$ for different values of $K$ . For small values of $E$ , the curves obtained from the full problem (3.15)–(3.16) are tangential to the asymptotic result deduced in appendix A for weakly modified  $\unicode[STIX]{x1D703}_{m}$ , see (A 5) with (A 8), and is shown as dashed lines. Figure 3(b) shows plots of $\unicode[STIX]{x1D703}_{m}$ as $K$ is varied between 0.1 and 25 for different values of $E$ . As $K$ increases, $\unicode[STIX]{x1D703}_{m}$ gradually approaches Young’s angle. For the smaller values of $E$ we have good agreement between the full numerics (solid curves) and the asymptotic result (A 5) with (A 8) (dashed curves), which improves as $K$ is further increased. When $E=10$ and $E=100$ , the small- $E$ asymptotic result, equation (A 5) with (A 8), over-predicts $\unicode[STIX]{x1D703}_{m}$ and for this reason these curves are discarded. This is to be expected given that these values of $E$ cannot possibly conform with the conditions stated in appendix A, which are necessary for the asymptotic result to hold. From these plots it is clear that $\unicode[STIX]{x1D703}_{m}$ increases as $E$ increases and as $K$ decreases. Thus evaporation enhances $\unicode[STIX]{x1D703}_{m}$ but its effect is diminished if the kinetic resistance effects become too large. These results are consistent, at least qualitatively, with those of Rednikov et al. (Reference Rednikov, Rossomme and Colinet2009), who used a disjoining pressure model for the contact line dynamics, the main difference being the presence of a weak singularity of $\unicode[STIX]{x1D703}_{m}$ in the present model, manifested as $K\rightarrow 0$ .

As $E\rightarrow \infty$ , we find that $\unicode[STIX]{x1D703}_{m}\gg 1$ and $\unicode[STIX]{x1D703}_{m}$ tends asymptotically to $\unicode[STIX]{x1D701}(K)E^{1/4}$ , where $\unicode[STIX]{x1D701}(K)$ is a function of $K$ to be determined numerically (see appendix B.2). In this limit, we find that the original, non-scaled, evaporation-modified Young’s angle, $\unicode[STIX]{x1D717}_{m}=\unicode[STIX]{x1D703}_{m}\unicode[STIX]{x1D717}_{s}$ , will no longer depend on the actual value $\unicode[STIX]{x1D717}_{s}$ , because $E$ scales with $\unicode[STIX]{x1D717}_{s}^{-4}$ . The results corresponding to the strong-evaporation result, $\unicode[STIX]{x1D703}_{m}=\unicode[STIX]{x1D701}(K)E^{1/4}$ , are represented in figure 3 as black dotted curves and appear to yield a very good estimate even for $E=O(1)$ , provided that $K$ is sufficiently small. Noteworthy is also that similar $E^{1/4}$ power-law behaviours in evaporation-modified angles have also been reported in different evaporation settings, e.g. by Morris (Reference Morris2001) in a setting which is equivalent to taking $K\gg 1$ (see also Todorova et al. Reference Todorova, Thiele and Pismen2012), but not in an explicit context of weak or strong evaporation, and also by Rednikov et al. (Reference Rednikov, Rossomme and Colinet2009) in the weak evaporation limit for a different model. Interestingly, the recent analysis of Saxton et al. (Reference Saxton, Whiteley, Vella and Oliver2016) uncovered a remarkably similar, $E^{2/7}$ power-law behaviour in the diffusion-limited evaporation scenario with $E\gg 1$ .

It is important to note that evaporation-modified contact angles have also been considered in a similar setting by Hocking (Reference Hocking1995) and it is thus of interest to put his findings in the present context. More specifically, motivated by the work of Anderson & Davis (Reference Anderson and Davis1995), Hocking focused on steady two-dimensional droplet shapes in which the fluid lost due to evaporation was replaced by a flux of fluid through the substrate so that there was no contact line motion and no volume variations. He used the same ingredients for evaporation as here, but retained the evaporative term in his outer-region analysis, which is neglected here. Another key difference is the assumed order of magnitude of the kinetic resistance to evaporation. Here, this non-equilibrium effect is treated to be small and only appreciable in the vicinity of the contact line, just as slip is. In contrast, Hocking (Reference Hocking1995) formally treated the kinetic resistance to be significant at the macro-scale, although in the end he aimed at the small- $K$ limit (in his terms). Hence, unsurprisingly, Hocking needed to obtain the evaporation-modified angles from the coupling of the micro- and macro-scales, unlike here, where $\unicode[STIX]{x1D703}_{m}$ is solely determined from the inner-region asymptotics. Clearly, Hocking’s asymptotic treatment and ours correspond to two distinct parameter regimes, neither of which is a particular case of the other. Nevertheless there can exist a region in the parameter space where the two asymptotic cases overlap, when the kinetic resistance is small compared to the baseline case of Hocking (Reference Hocking1995), but large in terms of the present study, i.e. when $K\gg 1$ and $\unicode[STIX]{x1D706}K\ll 1$ . The overlap region in question is described by Hocking’s (Reference Hocking1995), equation (14), where we note a typo in the denominator of the logarithmic term, with the factor 2 to be replaced by 3. When rendered in our notation and scalings and with the typo rectified, equation (14) in Hocking (Reference Hocking1995) can be written as

(3.17) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{m}^{4}=1+\frac{4E}{K}\ln \frac{\text{e}K}{\unicode[STIX]{x1D703}_{m}},\end{eqnarray}$$

which is a transcendental equation for $\unicode[STIX]{x1D703}_{m}$ . In accordance with Hocking’s derivation, equation (3.17) is valid, in our terms, for $K\gg 1$ and $E/K\ll 1$ and generally for $(E/K)\ln K=O(1)$ , allowing for $O(1)$ deviations from Young’s angle, scaled to unity here. In a sense, equation (3.17) can be viewed as the counterpart of the work of Morris (Reference Morris2001) carried out in the limit $K\gg 1$ . Although, as we mentioned, equation (3.17) still holds when $\unicode[STIX]{x1D703}_{m}-1=O(1)$ , it simplifies to $\unicode[STIX]{x1D703}_{m}=1+(E/K)\ln (\text{e}K)$ in the particular case $\unicode[STIX]{x1D703}_{m}-1\ll 1$ . This two-term asymptotic result matches perfectly with our result for weakly evaporation-modified angles, equation (A 5), when $\unicode[STIX]{x1D6FC}$ , given by (A 8), is expanded for $K\gg 1$ . To demonstrate this overlap, Hocking’s asymptotic result, equation (3.17), is represented in figure 3 by grey dotted curves and is seen to agree well within the domain of its validity with the full computation result (solid curves).

Looking at the $O({\dot{r}})$ terms of (3.13) using (3.14) gives a linear problem to determine  $\unicode[STIX]{x1D719}_{1}$ ,

(3.18) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D719}_{0}+\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}[\unicode[STIX]{x1D719}_{0}^{2}(\unicode[STIX]{x1D719}_{0}+1)\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}^{3}\unicode[STIX]{x1D719}_{1}+\unicode[STIX]{x1D719}_{0}(3\unicode[STIX]{x1D719}_{0}+2)\unicode[STIX]{x1D719}_{1}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}^{3}\unicode[STIX]{x1D719}_{0}]=\frac{E\unicode[STIX]{x1D719}_{1}}{(\unicode[STIX]{x1D719}_{0}+K)^{2}},\end{eqnarray}$$

which is treated with homogeneous conditions at $\unicode[STIX]{x1D702}=0$ ,

(3.19a ) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{1}=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D719}_{1}=0\end{eqnarray}$$

and by requiring that $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D719}_{1}$ is no more than logarithmically large as $\unicode[STIX]{x1D702}\rightarrow \infty$ , where the last condition and the asymptotics of (3.18) and $\unicode[STIX]{x1D719}_{0}$ dictate that

(3.19b ) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D719}_{1}\sim \frac{1}{\unicode[STIX]{x1D703}_{m}^{2}}\ln \unicode[STIX]{x1D6FD}\unicode[STIX]{x1D702}\quad \text{as}~\unicode[STIX]{x1D702}\rightarrow \infty .\end{eqnarray}$$

Here, $\unicode[STIX]{x1D6FD}$ is a degree of freedom which needs to be determined as part of the solution for given values of $E$ and $K$ (see appendix B.3 for a brief discussion on the implementation of the numerical scheme for  $\unicode[STIX]{x1D719}_{1}$ ).

Combining (3.16c ) and (3.19b ) and going back to the original variables, we find

(3.20) $$\begin{eqnarray}-\unicode[STIX]{x2202}_{x}h\sim \unicode[STIX]{x1D703}_{m}+\frac{{\dot{r}}}{\unicode[STIX]{x1D703}_{m}^{2}}\ln \left(\unicode[STIX]{x1D6FD}\frac{r-x}{\unicode[STIX]{x1D706}}\right)\quad \text{as}~(r-x)/\unicode[STIX]{x1D706}\rightarrow \infty ,\end{eqnarray}$$

which is to be matched with the asymptotics of the outer solution, (3.11). It is clear from (3.20) that $\unicode[STIX]{x1D706}/\unicode[STIX]{x1D6FD}$ is a measure of the size of the inner region and it is thus of interest to explore how $\unicode[STIX]{x1D6FD}$ varies as a function of $E$ and $K$ . Figure 4 shows plots of $\ln \unicode[STIX]{x1D6FD}$ as $E$ and $K$ are varied and, unlike the corresponding plots of $\unicode[STIX]{x1D703}_{m}$ , the dependence on these parameters is non-monotonic. In figure 4(a) $\ln \unicode[STIX]{x1D6FD}$ exhibits an initial decrease followed by an increase as $E$ becomes large, but this behaviour is delayed for large $K$ . A similar behaviour is observed in figure 4(b) as $K$ is varied while $E$ is kept constant, but we now have that $\ln \unicode[STIX]{x1D6FD}\rightarrow 1$ as $K\rightarrow \infty$ , which corresponds to the non-volatile case (Hocking Reference Hocking1983). The approach to this asymptotic behaviour becomes more gradual for larger values of  $E$ .

Figure 4. (a) Plots of $\ln \unicode[STIX]{x1D6FD}$ as a function of $E$ for various values of  $K$ . (b) Plots of $\ln \unicode[STIX]{x1D6FD}$ as a function of $K$ for various values of  $E$ .

At this point, we note that the differences in the contact line model employed, be it a precursor or a slip model, will only manifest themselves in the inner-region dynamics, when the droplet is sufficiently large. Thus, had we used, for example, a disjoining pressure model in our formulation we would have anticipated the same behaviour as in (3.20) with the slip length being replaced by the associated micro-scale of the disjoining pressure model, say, the thickness of the precursor film (Savva & Kalliadasis Reference Savva and Kalliadasis2011) or a length scale obtained by the Hamaker constant, surface tension and the contact angle (Colinet & Rednikov Reference Colinet and Rednikov2011). Similar conclusions can also be drawn with other popular contact line models (see, e.g. Sibley et al. Reference Sibley, Savva and Kalliadasis2012, Reference Sibley, Nold, Savva and Kalliadasis2013, Reference Sibley, Nold, Savva and Kalliadasis2015b ).

3.1.3 Matching

One readily observes that the outer solution, equation (3.11), cannot match with the inner one, equation (3.20), since the $x$ -dependent logarithmic terms have different coefficients. This issue can be resolved with an intermediate layer sandwiched between the inner and outer regions as shown in figure 1 (Hocking Reference Hocking1983). As noted in previous studies (see, e.g. Savva & Kalliadasis Reference Savva and Kalliadasis2009; Sibley et al. Reference Sibley, Nold and Kalliadasis2015a ), this intermediate region ultimately justifies why matching can be carried out simply by considering the cubes of the outer and inner slopes, equations (3.11) and (3.20), respectively. Matching inner and outer solutions can be accomplished without an intermediate region by an alternative path originally proposed by Lacey (Reference Lacey1982). In this case, matching is performed not only for the first few terms of the asymptotic expansion, but effectively for the infinite series of the inner and outer solutions (for a discussion placed in the context of contact line motion with mass transfer, see the work of Oliver et al. Reference Oliver, Whiteley, Saxton, Vella, Zubkov and King2015). In a more recent development, Sibley et al. (Reference Sibley, Nold and Kalliadasis2015a ) discuss this alternative method of matching inner and outer solutions for contact line motion in a more general setting, explaining how it can be performed in a straightforward manner for more complicated problems for which simply considering the cubes of the slopes argument does not work. As expected, however, both approaches yield the same results.

In the present problem, considering the cubes of (3.11) and (3.20) does allow us to cancel the $x$ -dependent logarithmic terms and by matching the terms that are constant in $x$ we obtain the following evolution equation for the droplet radius

(3.21) $$\begin{eqnarray}{\dot{r}}=\frac{\unicode[STIX]{x1D703}^{3}-\unicode[STIX]{x1D703}_{m}^{3}}{3\ln \displaystyle \frac{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}r}{\unicode[STIX]{x1D706}\text{e}^{2}}},\end{eqnarray}$$

which is valid for ${\dot{r}}=O(1/|\!\ln \unicode[STIX]{x1D706}|)$ and is determined with $O(1/|\text{ln}\,\unicode[STIX]{x1D706}|^{3})$ error as $\unicode[STIX]{x1D706}\rightarrow 0$ . Although a detailed consideration of an intermediate region is omitted here, as done, for example in other works studying contact line motion in different settings (see, e.g. Eggers Reference Eggers2005b ; Savva & Kalliadasis Reference Savva and Kalliadasis2014; Xu & Jensen Reference Xu and Jensen2016), the matching procedure can be rigorously justified by following nearly identical arguments as the ones presented in the above-mentioned works.

Equation (3.21) tells us that the contact line will move an order-unity distance within $t=O(|\!\ln \unicode[STIX]{x1D706}|)$ as $\unicode[STIX]{x1D706}\rightarrow 0$ , but we will also argue heuristically later on (see § 4.2) that it can be extended to the evaporation stage as well. Before continuing with our analysis, we remark that a similar analysis undertaken by Anderson & Davis (Reference Anderson and Davis1995) with a one-sided evaporation model assumed as a starting point a functional relation between ${\dot{r}}$ and $\unicode[STIX]{x1D703}$ ; here (3.21) is the outcome of the matched asymptotic analysis we have undertaken. Moreover, recall that the $\unicode[STIX]{x1D703}_{m}^{4}$ law obtained through matching by Hocking (Reference Hocking1995), equation (3.17), is for evaporation-modified angles and not for motion-modified ones, which results to a $\unicode[STIX]{x1D703}^{3}$ term, as seen in (3.21). Lastly, we note that an analogous result was obtained by Saxton et al. (Reference Saxton, Whiteley, Vella and Oliver2016) in the diffusion-limited evaporation case, but without the $O(1/\ln ^{2}\unicode[STIX]{x1D706})$ terms as retained here.

On the other hand, an asymptotically consistent expansion for ${\dot{r}}$ in which $O(1/\ln ^{2}\unicode[STIX]{x1D706})$ terms are retained requires, in principle, the inclusion of an $O(1/|\!\ln \unicode[STIX]{x1D706}|)$ correction to $r(0)$ to be used as an initial condition to (3.21), as compared to the value of $r(0)$ for the full problem. This correction would arise from a (mostly numerical) treatment of the capillary action stage which occurs for $t=O(1)$ (see also Oliver et al. Reference Oliver, Whiteley, Saxton, Vella, Zubkov and King2015; Ren et al. Reference Ren, Trinh and Weinan2015; Saxton et al. Reference Saxton, Whiteley, Vella and Oliver2016). Yet, although our simulations repeatedly show that the capillary action phase and the initial droplet profile are unimportant for the subsequent dynamics, one can take the $r(0)$ correction to be a priori negligibly small by considering initial drop shapes which are already sufficiently close to the quasi-static profiles, equation (3.4).

3.2 Evolution of the droplet volume

Equation (3.21) depends on the droplet volume, $v(t)$ , through the expression for $\unicode[STIX]{x1D703}$ , equation (3.5). Hence, to close the system, we will utilise (2.20) to obtain an evolution equation for $\dot{v}$ in terms of $r(t)$ , $v(t)$ and the system parameters, ${\mathcal{E}}$ , ${\mathcal{K}}$ , $\unicode[STIX]{x1D706}$ and  $B$ .

As a first approximation, we may assume that the fine details coming from the contact line region are negligible. Hence, we may take

(3.22) $$\begin{eqnarray}\dot{v}=-\int _{0}^{r}\frac{{\mathcal{E}}x}{h_{0}+{\mathcal{K}}}\,\text{d}x,\end{eqnarray}$$

where $h_{0}$ is the leading-order outer solution, equation (3.4). This expression is not amenable to further analysis and one needs to compute this integral at each time, unless we consider the $B=0$ case, so that we may use (3.6a ) for $h_{0}$ in (3.22) to obtain

(3.23) $$\begin{eqnarray}\dot{v}=-\frac{{\mathcal{E}}r}{\unicode[STIX]{x1D703}}\ln \left(\frac{\unicode[STIX]{x1D703}r}{2{\mathcal{K}}}+1\right).\end{eqnarray}$$

From (3.23) we can easily deduce that $\dot{v}=O(\unicode[STIX]{x1D706}|\!\ln \unicode[STIX]{x1D706}|)$ as $\unicode[STIX]{x1D706}\rightarrow 0$ for moderate values of $E$ and $K$ . However, equation (3.23) is not expected to work well for all the parameter regimes of interest. The reason is because the influence of the inner region near the contact line is not always accounted for in sufficient detail. To properly do this, we can pursue further our asymptotic approach that relies on the assumption that the droplet size is much larger than ${\mathcal{K}}$ and $\unicode[STIX]{x1D706}$ . Ultimately, however, as the droplet volume diminishes, this formula will fail to describe the long-term behaviour of the volume, just as (3.21) will fail in this limit.

Unlike the matching carried out to obtain ${\dot{r}}$ , to estimate $\dot{v}$ we need to consider all three regions, inner, intermediate and outer. This is particularly true if (2.20) is to be evaluated at the spreading stage, although, as alluded to earlier, $v(t)$ is not altered appreciably during this stage. Evaporation effects manifest themselves during the third stage, where $\unicode[STIX]{x1D703}\approx \unicode[STIX]{x1D703}_{m}$ . As the transition to the evaporation stage occurs, the extent of the intermediate region shrinks and can be omitted altogether in the last two stages. However, in what follows we retain the presence of the intermediate region contributions for the sake of completeness.

Figure 5. A schematic diagram showing the three regions in the vicinity of the contact line at $x=r$ (not drawn to scale), distinguishing between the apparent, macroscopic angle, $\unicode[STIX]{x1D703}$ , the microscopic, static angle, $\unicode[STIX]{x1D703}_{s}=1$ (normalised to unity in our chosen non-dimensionalisation) and the evaporation-modified Young’s angle,  $\unicode[STIX]{x1D703}_{m}$ .

Based on the above, we expect that all scales contribute to the net evaporative flux at leading order. To proceed, we split $\dot{v}$ into three parts, each consisting of an integral carried out in each region, namely

(3.24) $$\begin{eqnarray}\dot{v}=q_{1}+q_{2}+q_{3},\end{eqnarray}$$

where

(3.25a-c ) $$\begin{eqnarray}q_{1}=-\int _{\tilde{r}}^{r}\frac{x\unicode[STIX]{x1D706}E}{h+\unicode[STIX]{x1D706}K}\,\text{d}x,\quad q_{2}=-\int _{r_{\ast }}^{\tilde{r}}\frac{x\unicode[STIX]{x1D706}E}{h+\unicode[STIX]{x1D706}K}\,\text{d}x\quad \text{and}\quad q_{3}=-\int _{0}^{r_{\ast }}\frac{x\unicode[STIX]{x1D706}E}{h+\unicode[STIX]{x1D706}K}\,\text{d}x,\end{eqnarray}$$

with $r_{\ast }$ and $\tilde{r}$ being the radii where the intermediate region matches within appropriate overlap regions with the solution in the outer and inner regions, respectively (see figure 5), i.e. such that $\tilde{\unicode[STIX]{x1D702}}=(r-\tilde{r})/\unicode[STIX]{x1D706}\gg 1$ and $r-r_{\ast }\ll 1$ , but, at the same time, we also have ${\dot{r}}\ln \tilde{\unicode[STIX]{x1D702}}\ll 1$ and ${\dot{r}}\ln (r-r_{\ast })\ll 1$ . Based on these requirements for the asymptotics of $\tilde{r}$ and  $r_{\ast }$ , each integral in (3.25) is estimated by making use of the smallness of ${\dot{r}}$ and the appropriate asymptotic limits, retaining only the leading algebraic-order terms in  $\unicode[STIX]{x1D706}$ .

The integral over the inner region is estimated using

(3.26) $$\begin{eqnarray}q_{1}=-\int _{\tilde{r}}^{r}\frac{x\unicode[STIX]{x1D706}E}{h_{\mathit{in}}+\unicode[STIX]{x1D706}K}\,\text{d}x\sim -\int _{0}^{\tilde{\unicode[STIX]{x1D702}}}\frac{r{\mathcal{E}}}{\unicode[STIX]{x1D719}_{0}+K}\,\text{d}\unicode[STIX]{x1D702},\end{eqnarray}$$

as $\unicode[STIX]{x1D706}\rightarrow 0$ , where we took $h_{\mathit{in}}=\unicode[STIX]{x1D706}\unicode[STIX]{x1D719}_{0}(\unicode[STIX]{x1D702})$ , with $\unicode[STIX]{x1D702}=(r-x)/\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D719}_{0}$ determined from (3.15). Using the asymptotics of $\unicode[STIX]{x1D719}_{0}$ as $\unicode[STIX]{x1D702}\rightarrow \infty$ in the left-hand side of (3.15) and integrating, we find that (3.26) behaves like

(3.27) $$\begin{eqnarray}q_{1}\sim -\frac{r{\mathcal{E}}}{\unicode[STIX]{x1D703}_{m}}\ln (\tilde{\unicode[STIX]{x1D702}}\tilde{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D703}_{m}\text{e}^{-3/2})\quad \text{as}~\tilde{\unicode[STIX]{x1D702}}\rightarrow \infty ,\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D6FD}}$ is a parameter appearing in the next term of the asymptotics of $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D719}_{0}$ as a function of $\unicode[STIX]{x1D719}_{0}$ , namely

(3.28) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D719}_{0}\sim \unicode[STIX]{x1D703}_{m}-\frac{E}{2\unicode[STIX]{x1D703}_{m}^{3}\unicode[STIX]{x1D719}_{0}}\ln (\unicode[STIX]{x1D719}_{0}\tilde{\unicode[STIX]{x1D6FD}})\quad \text{as}~\unicode[STIX]{x1D719}_{0}\rightarrow \infty\end{eqnarray}$$

and is determined using similar techniques as those discussed in appendix B. Had we merely used $h_{\mathit{in}}=\unicode[STIX]{x1D706}\unicode[STIX]{x1D703}_{m}\unicode[STIX]{x1D702}$ in (3.26), which is the leading-order behaviour as we approach the intermediate region, we would have obtained

(3.29) $$\begin{eqnarray}q_{1}\sim -\frac{r{\mathcal{E}}}{\unicode[STIX]{x1D703}_{m}}\ln \frac{\unicode[STIX]{x1D703}_{m}\tilde{\unicode[STIX]{x1D702}}}{K}\quad \text{as}~\tilde{\unicode[STIX]{x1D702}}\rightarrow \infty ,\end{eqnarray}$$

which is equivalent to the inner-region contributions of (3.22) with $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{m}$ .

Figure 6. (a) Solid curves: plots of $\ln \tilde{\unicode[STIX]{x1D6FD}}$ as a function of $E$ for various values of $K$ and $0.001\leqslant E\leqslant 30$ ; dashed curves: plots of $\ln (\text{e}^{3/2}/K)$ corresponding to each of the values of $K$ shown as solid curves, to compare (3.27) and (3.29). (b) Plots of $\ln \tilde{\unicode[STIX]{x1D6FD}}$ as a function of $K$ for various values of  $E$ . The curves from bottom to top correspond to the values $E=0.1$ , 1, 10 and 100, respectively.

In figure 6 we show a few representative plots of $\tilde{\unicode[STIX]{x1D6FD}}$ for various values of $E$ and  $K$ . In figure 6(a) we provide plots for $10^{-3}\leqslant E\leqslant 30$ , to contrast the values of $\tilde{\unicode[STIX]{x1D6FD}}$ with $\text{e}^{3/2}/K$ as a means to compare (3.27) and (3.29). It is readily seen that (3.29) gives a rough estimate for $q_{1}$ which works better for smaller values of $E$ and larger values of $K$ , i.e. for weaker evaporation effects. This is not surprising given that in this limit $\unicode[STIX]{x1D719}_{0}\approx \unicode[STIX]{x1D702}$ with $\unicode[STIX]{x1D703}_{m}\approx 1$ (cf. appendix A) is no different from the simplified form used to obtain (3.29). In particular, it is interesting to note that the calculations presented in figure 6(a) suggest that $\tilde{\unicode[STIX]{x1D6FD}}\rightarrow \text{e}^{3/2}/K$ as $E\rightarrow 0^{+}$ .

To estimate $q_{2}$ , perform a change of variable to the stretched coordinate $\unicode[STIX]{x1D702}$ and consider to leading order as $\unicode[STIX]{x1D706}\rightarrow 0$

(3.30) $$\begin{eqnarray}q_{2}\sim -\int _{\tilde{\unicode[STIX]{x1D702}}}^{\unicode[STIX]{x1D702}_{\ast }}\frac{\unicode[STIX]{x1D706}^{2}rE}{h_{\mathit{int}}}\,\text{d}\unicode[STIX]{x1D702},\end{eqnarray}$$

where $\unicode[STIX]{x1D702}_{\ast }=(r-r_{\ast })/\unicode[STIX]{x1D706}$ and $h_{\mathit{int}}$ is the intermediate solution

(3.31) $$\begin{eqnarray}h_{\mathit{int}}=\unicode[STIX]{x1D706}\unicode[STIX]{x1D702}\left(\unicode[STIX]{x1D703}_{m}^{3}+3{\dot{r}}\ln \frac{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D702}}{\text{e}}\right)^{1/3}\end{eqnarray}$$

as expressed in terms of the inner variable $\unicode[STIX]{x1D702}$ . The derivation of (3.31) is nearly identical to the intermediate solution of Hocking (Reference Hocking1983) for non-volatile droplets, arguing, as in Saxton et al. (Reference Saxton, Whiteley, Vella and Oliver2016), that the evaporation term is negligible within the intermediate region and matching the solution with appropriate conditions at the far field of the inner region. In (3.30) we have neglected the kinetic resistance effects which is an acceptable approximation for ${\mathcal{K}}\ll 1$ . Hence, we obtain

(3.32) $$\begin{eqnarray}q_{2}\sim -\frac{r{\mathcal{E}}}{2{\dot{r}}}\left[\left(\unicode[STIX]{x1D703}_{m}^{3}+3{\dot{r}}\ln \frac{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D702}_{\ast }}{\text{e}}\right)^{2/3}-\left(\unicode[STIX]{x1D703}_{m}^{3}+3{\dot{r}}\ln \frac{\unicode[STIX]{x1D6FD}\tilde{\unicode[STIX]{x1D702}}}{\text{e}}\right)^{2/3}\right].\end{eqnarray}$$

Next, we need to find the asymptotic behaviour of $q_{2}$ as we have both $\tilde{\unicode[STIX]{x1D702}}\rightarrow \infty$ and $r_{\ast }\rightarrow r$ . Using (3.21) we modify the first term in the square brackets of (3.32) so that we can eventually take the limit $r_{\ast }\rightarrow r$ . Hence we obtain

(3.33) $$\begin{eqnarray}q_{2}\sim -\frac{r{\mathcal{E}}}{2{\dot{r}}}\left\{\left[\unicode[STIX]{x1D703}^{3}+3{\dot{r}}\ln \left(\frac{\text{e}}{\unicode[STIX]{x1D6FF}}\frac{r-r_{\ast }}{r}\right)\right]^{2/3}-\left(\unicode[STIX]{x1D703}_{m}^{3}+3{\dot{r}}\ln \frac{\unicode[STIX]{x1D6FD}\tilde{\unicode[STIX]{x1D702}}}{\text{e}}\right)^{2/3}\right\},\end{eqnarray}$$

noting that the limits $\tilde{\unicode[STIX]{x1D702}}\rightarrow \infty$ and $r_{\ast }\rightarrow r$ must be taken with the understanding that the $O({\dot{r}})$ terms of the asymptotic expansions are still of higher order compared to the $O(1)$ terms. Hence, expanding into a series in ${\dot{r}}$ gives as $\unicode[STIX]{x1D706}\rightarrow 0$

(3.34) $$\begin{eqnarray}q_{2}\sim -r{\mathcal{E}}\left\{\frac{3(\unicode[STIX]{x1D703}+\unicode[STIX]{x1D703}_{m})}{2(\unicode[STIX]{x1D703}^{2}+\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}_{m}+\unicode[STIX]{x1D703}_{m}^{2})}\ln \frac{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}r}{\unicode[STIX]{x1D706}\text{e}^{2}}+\frac{1}{\unicode[STIX]{x1D703}}\ln \left[\frac{\text{e}}{\unicode[STIX]{x1D6FF}}\left(1-\frac{r_{\ast }}{r}\right)\right]-\frac{1}{\unicode[STIX]{x1D703}_{m}}\ln \frac{\unicode[STIX]{x1D6FD}\tilde{\unicode[STIX]{x1D702}}}{\text{e}}\right\}.\end{eqnarray}$$

Lastly, to find the leading-order expression for the integral within the outer region, $q_{3}$ , we neglect, just as in the intermediate region, kinetic resistance effects and consider

(3.35) $$\begin{eqnarray}q_{3}\sim -\int _{0}^{r_{\ast }}\frac{x{\mathcal{E}}}{h_{0}}\,\text{d}x\end{eqnarray}$$

as $\unicode[STIX]{x1D706}\rightarrow 0$ , where $h_{0}$ is the leading-order outer solution given by (3.4). Clearly the integrand has a logarithmic singularity as $r_{\ast }\rightarrow r$ , since $h_{0}$ vanishes as $r_{\ast }\rightarrow r$ . Hence, we can remove the singularity as $r_{\ast }\rightarrow r$ outside the integrand by writing

(3.36) $$\begin{eqnarray}q_{3}\sim -\frac{{\mathcal{E}}}{\unicode[STIX]{x1D703}}\left\{\int _{0}^{r}\left[\frac{x\sqrt{B}I_{1}(c)}{I_{0}(c)-I_{0}(x\sqrt{B})}-\frac{r}{r-x}\right]\text{d}x-r\ln \left(1-\frac{r_{\ast }}{r}\right)\right\}\end{eqnarray}$$

as $r_{\ast }\rightarrow r$ . In the limit $B\rightarrow 0$ , equation (3.36) simplifies to

(3.37) $$\begin{eqnarray}q_{3}\sim \frac{r{\mathcal{E}}}{\unicode[STIX]{x1D703}}\ln \left[2\left(1-\frac{r_{\ast }}{r}\right)\right].\end{eqnarray}$$

Adding the contributions from (3.27), (3.34) and (3.36) cancels the singular terms involving $\ln \tilde{\unicode[STIX]{x1D702}}$ and $\ln (r-r_{\ast })$ and we obtain an expression for $\dot{v}$ determined up to and including $O(\unicode[STIX]{x1D706})$ terms,

(3.38) $$\begin{eqnarray}\dot{v}=-r{\mathcal{E}}\left\{\frac{1}{\unicode[STIX]{x1D703}_{m}}\ln \frac{\tilde{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D703}_{m}}{\unicode[STIX]{x1D6FD}\text{e}^{1/2}}+\frac{1}{\unicode[STIX]{x1D703}}\ln \frac{\text{e}}{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}+\frac{3(\unicode[STIX]{x1D703}+\unicode[STIX]{x1D703}_{m})}{2(\unicode[STIX]{x1D703}^{2}+\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}_{m}+\unicode[STIX]{x1D703}_{m}^{2})}\ln \frac{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}r}{\unicode[STIX]{x1D706}\text{e}^{2}}\right\},\end{eqnarray}$$

where

(3.39) $$\begin{eqnarray}\ln \unicode[STIX]{x1D6FE}=\int _{0}^{1}\left[\frac{1}{1-x}-\frac{xcI_{1}(c)}{I_{0}(c)-I_{0}(xc)}\right]\text{d}x,\end{eqnarray}$$

noting that $\unicode[STIX]{x1D6FE}\rightarrow 2$ as $r\sqrt{B}\rightarrow 0$ (see figure 7 for plots of $\unicode[STIX]{x1D6FE}(c)$ for pendant and sessile droplets). To compute $\unicode[STIX]{x1D6FE}$ for pendant droplets we take $c=r\sqrt{-B}$ , replace the modified Bessel functions by Bessel functions and change the sign of the second term in (3.39). In the stages dominated by evaporation, where we can take $\unicode[STIX]{x1D703}\approx \unicode[STIX]{x1D703}_{m}$ , equation (3.38) simplifies further to

(3.40) $$\begin{eqnarray}\dot{v}=-\frac{{\mathcal{E}}r}{\unicode[STIX]{x1D703}}\ln \frac{\tilde{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D703}r}{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D706}\text{e}^{3/2}}.\end{eqnarray}$$

The arguments leading to the derivation of (3.38) reveal that it is able to capture the evolution of $v(t)$ both in the spreading and the evaporation stages as $\unicode[STIX]{x1D706}\rightarrow 0$ and for $O(1)$ values for $E$ , $K$ and $B$ . As it turns out, using the simpler expression (3.40) to simulate the full dynamics does not typically yield appreciably different results from using the more complete expression (3.38). This is mainly due to the fact that spreading takes place on a shorter time scale compared to the long evaporation stage during which $\unicode[STIX]{x1D703}\approx \unicode[STIX]{x1D703}_{m}$ .

Figure 7. Dependence of $\ln \unicode[STIX]{x1D6FE}$ on $r\sqrt{|B|}$ for a sessile ( $B>0$ ; solid curve) and a pendant ( $B<0$ ; dashed curve) droplet. As the droplet evaporates, $r\rightarrow 0$ , the curves approach asymptotically the value of $\ln \unicode[STIX]{x1D6FE}=\ln 2$ , corresponding to the $B=0$ case.

As already anticipated from our earlier discussion, equation (3.38) and, perhaps more transparently, equation (3.40) reveal that, in the distinguished limit under consideration, $\dot{v}=O(\unicode[STIX]{x1D706}|\!\ln \unicode[STIX]{x1D706}|)$ as $\unicode[STIX]{x1D706}\rightarrow 0$ . Comparing this scaling with the scaling ${\dot{r}}=O(1/|\!\ln \unicode[STIX]{x1D706}|)$ as $\unicode[STIX]{x1D706}\rightarrow 0$ , which was deduced from (3.21), confirms a posteriori the assumptions underpinning the analysis of the second stage, namely that $\unicode[STIX]{x1D706}\ll {\dot{r}}\ll 1$ and $\dot{v}\ll {\dot{r}}$ . However it is clear that our analysis will break down when the argument of the logarithm in (3.40) tends to unity with $r\rightarrow r_{c}$ , where

(3.41) $$\begin{eqnarray}r_{c}=\frac{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D706}\text{e}^{3/2}}{\tilde{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D703}_{m}}=O({\mathcal{K}}),\end{eqnarray}$$

i.e. when the size of the droplet becomes comparable to the kinetic length scale. When this happens, equation (3.40) incorrectly predicts that $\dot{v}\rightarrow 0$ with $v\neq 0$ . This indicates that ${\mathcal{K}}$ can no longer be assumed to be smaller than the droplet size which invalidates the assumptions put forth for our analysis to hold.

The derivation of an evolution equation for $v(t)$ concludes the asymptotic analysis we have undertaken. Although we concede that the arguments we presented above are centred around the spreading time scale, we did not deem necessary a separate treatment of the evaporation stage to produce a composite expansion valid across both stages, as done, for example in the analysis of Saxton et al. (Reference Saxton, Whiteley, Vella and Oliver2016). Later on, we will offer some heuristic arguments as to why doing so is acceptable (see § 4.2), but, more importantly, throughout the following section we will provide rather compelling numerical evidence from over one hundred simulations of the governing partial differential equation, which are found to be in excellent agreement with the predictions of the equations (3.21) and (3.38) obtained from matching.

4.1 The four-stage dynamics

Figure 8 shows a typical calculation which is performed for parameters that are close to the values for water reported in § 2. We see that the solutions of the full (solid curves) and reduced problems (dashed curves) are visually indistinguishable and the solutions are on top of each other; there are some barely noticeable differences at early times (see, e.g. figure 8 b), which are solely attributed to our choice of  $\unicode[STIX]{x1D716}$ . Indeed, if we repeat the calculation with a much smaller value of  $\unicode[STIX]{x1D716}$ , the dynamics is indistinguishable throughout. We should also note that using the simpler expression for  $\dot{v}$ , equation (3.23), instead of the more rigorously obtained (3.38) the agreement between the full and reduced problems degrades and the solutions no longer exhibit such excellent agreement, particularly during the late stages of the dynamics (see also § 4.3).

From figure 8 we see that the contact line advances initially up to $t=O(10)=O(1/|\!\ln \unicode[STIX]{x1D706}|)$ and then gradually recedes due to evaporation, with the apparent contact angle, defined by (3.6b ), being close to the evaporation-modified Young’s angle (see figure 8 b), as expected from our analysis. Looking at the various plots in figure 8, we can identify the four stages for the evaporation dynamics of the contact line, anticipated earlier in the discussion, at the beginning of § 3. The first stage, is typically very brief and corresponds to the relaxation of the droplet shape in the bulk to the leading-order outer solution (in this calculation the first stage lasts approximately up to $t=O(0.1)$ ). In the second stage, the dynamics is driven primarily by (3.21) with $v(t)\approx v(0)$ . This is evidenced in figure 8(d), where we show for comparison the evolution of the droplet radius in the non-volatile case. Here we see good agreement with the non-volatile case up to $t=O(1)$ , but the agreement degrades when $\unicode[STIX]{x1D703}^{3}$ becomes of the same order of magnitude as $\unicode[STIX]{x1D703}_{m}^{3}$ , which, interestingly, occurs long before the volume changes become appreciable. Due to the brief duration of the first two stages, the volume of the droplet hardly changes (see, e.g. the droplet profiles ‘ $a$ ’–‘ $c$ ’ in figure 8 a and the plots of $r$ and $v$ against a logarithmic time scale in figure 8 d).

The distinguishing characteristic of the third stage is that the apparent contact angle remains close to $\unicode[STIX]{x1D703}_{m}$ . During this stage, evaporation effects dominate the dynamics and most of the liquid of the droplet evaporates into the saturated atmosphere (see, e.g. the droplet profiles ‘ $d$ ’–‘ $h$ ’ in figure 8 a and the plots of $r$ and $v$ against a linear time scale in figure 8 c). At the fourth stage, the droplet becomes too small and all assumptions upon which our analysis is based no longer hold. Hence, this stage is not captured by the reduced model; the computations of the full problem suggest that during this final stage $\unicode[STIX]{x1D703}$ decreases abruptly below $\unicode[STIX]{x1D703}_{m}$ (see figure 8 b). These final moments of the droplet’s lifetime are admittedly rather difficult to resolve with high accuracy due to the computation being prone to round-off errors and the aforementioned difficulties associated with the stiffness of the full problem. Nevertheless, it is clear that this stage will be rather short in duration and is thus deemed of lesser importance in our discussion.

Figure 9. Evaporating droplet with $E=9800$ , $K=14$ , $B=0$ , $\unicode[STIX]{x1D706}=2\times 10^{-5}$ , $\unicode[STIX]{x1D716}=10^{-2}$ and $v(0)=r(0)=1$ . (a) Evolution of the droplet radius and volume. (b) Evolution of the apparent contact angle, equation (3.6b ). In (a,b) the solutions to both the full (solid curves) and the reduced problems (dashed curves) are virtually indistinguishable.

In a second example, we test the limits of applicability of our analysis by considering the case of an FC-72-like droplet discussed in § 2. More specifically, if we consider the case when $\unicode[STIX]{x0394}\unicode[STIX]{x1D6E9}=10$  K, $\unicode[STIX]{x1D717}_{s}=4^{\circ }$ , $d=2.2$  mm, while all other parameters are kept unchanged, we have approximately the following parameters, $E=9800$ ( ${\mathcal{E}}=0.20$ ), $K=14$ and $\unicode[STIX]{x1D706}=2\times 10^{-5}$ . This is a rather extreme scenario given that the value of ${\mathcal{E}}$ is comparatively large and perhaps beyond the regime of applicability of our asymptotic theory. Yet, as figure 9 shows, our theory still applies and accurately predicts the dynamic behaviour of the full problem. Even in this case we have nearly indistinguishable dynamics throughout. However, here we see that the droplet recedes from the outset due to the evaporation-modified Young’s angle being larger than the initial apparent contact angle (as also remarked by Todorova et al. Reference Todorova, Thiele and Pismen2012). This calculation is also an example which illustrates that the distinction and duration of the previously discussed stages are strongly dependent on the parameters of the problem. The estimated evaporation time in this case is $t_{\ast }\approx 8.8$ . Based on the parameters of the problem, we can conclude that a droplet of volume $4.7~\unicode[STIX]{x03BC}\text{l}$ would evaporate completely within 10 s, a figure which is in the right ballpark based on the discussion in the work of Raj et al. (Reference Raj, Kunkelmann, Stephan, Plawsky and Kim2012).

4.2 The evaporation stage

From the results presented above, we saw that (3.21) and (3.38) are able to adequately describe the dynamics throughout. We also saw that the droplet evaporates completely in finite time after a comparatively long third stage, which is the stage during which the apparent contact angle is roughly equal to the evaporation-modified Young’s angle. Thus, if we simply take $\unicode[STIX]{x1D703}\approx \unicode[STIX]{x1D703}_{m}$ , equation (3.21) gives ${\dot{r}}\approx 0$ , i.e. the contact line velocity is very small and the contact line motion is slaved to the slow process of evaporation, so that ${\dot{r}}\approx 2\dot{v}/(3\unicode[STIX]{x1D703}_{m}^{1/3}v^{2/3})$ (when $B=0$ ). This means that at the onset of the third stage, where $v=O(1)$ , both ${\dot{r}}$ and $\dot{v}$ are of the same order and hence ${\dot{r}}$ ceases to be the greatest smallness parameter, and, self-consistently, one can neglect its active role when evaporative effects dominate. However, as $v\rightarrow 0^{+}$ , it follows that, just as in the second stage, we have that ${\dot{r}}\gg \dot{v}$ . This means that (3.21), which arises from the quasi-static analysis of the second stage may become relevant again. Thus we can reasonably expect that both the reduced problem, equations (3.21) and (3.40), and the modified problem in which (3.40) is solved together with $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{m}$ (see also Saxton et al. Reference Saxton, Whiteley, Vella and Oliver2016) will yield comparable results. Undertaking a more rigorous approach to justify this argument is likely to be rather unwieldy, so instead we resorted to making comparisons of the solutions of the reduced and modified problems with the full problem when $B=0$ with $r(0)=2/\unicode[STIX]{x1D703}_{m}^{1/3}$ and $v(0)=1$ for various choices of $E$ and $K$ . The particular choice for $r(0)$ was made to ensure that the droplet does not undergo a spreading stage. For all cases considered, the reduced problem was visually indistinguishable from the full problem, whereas with the modified problem there were small differences which became more noticeable as $r\rightarrow 0^{+}$ . This provided us with sufficient numerical evidence that the validity of (3.21) can be extended into the evaporation stage as well, as it is able to capture the transition of $\unicode[STIX]{x1D703}\rightarrow \unicode[STIX]{x1D703}_{m}$ and dynamics of the evaporation stage as $v\rightarrow 0^{+}$ , provided that $r\gg \unicode[STIX]{x1D706}$ .

On the other hand, by replacing (3.21) with $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{m}$ in the reduced problem, we can decouple the system of equations for $r(t)$ and $v(t)$ so that we can utilise just the equation for $\dot{v}$ . Since the droplet loses an order unity of its volume during the slow evaporation stage, this ultimately allows us to obtain an estimate of the time it takes a droplet of initial volume $v(0)=v_{0}$ to evaporate completely when $B=0$ . The evaporation time, $t_{\ast }$ , can be computed from (3.40) and using $r=2v^{1/3}\unicode[STIX]{x1D703}_{m}^{-1/3}$ , namely

(4.1) $$\begin{eqnarray}t_{\ast }=\frac{\unicode[STIX]{x1D703}_{m}^{4/3}}{2{\mathcal{E}}}\unicode[STIX]{x2A0D}_{0}^{v_{0}}\frac{\text{d}v}{v^{1/3}\ln \displaystyle \frac{\tilde{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D703}_{m}^{2/3}v^{1/3}}{\unicode[STIX]{x1D706}\text{e}^{3/2}}}.\end{eqnarray}$$

It is important to note that since (3.40) does not allow the volume to vanish completely, the integral for estimating $t_{\ast }$ only exists in the Cauchy principal value sense and as such is given above. Similar considerations can be followed to obtain $t_{\ast }$ when $B\neq 0$ , but the integral expression is cast with the radius being the variable of integration. Unlike the case with $B\neq 0$ , equation (4.1) can be evaluated analytically to obtain

(4.2) $$\begin{eqnarray}t_{\ast }=\frac{3\unicode[STIX]{x1D706}^{2}\text{e}^{3}}{2{\mathcal{E}}\tilde{\unicode[STIX]{x1D6FD}}^{2}}\text{Ei}\left(\ln \frac{v_{0}^{2/3}\tilde{\unicode[STIX]{x1D6FD}}^{2}\unicode[STIX]{x1D703}_{m}^{4/3}}{\unicode[STIX]{x1D706}^{2}\text{e}^{3}}\right),\end{eqnarray}$$

where $\text{Ei}(x)$ is the exponential integral, defined as

(4.3) $$\begin{eqnarray}\text{Ei}(x)=-\int _{-x}^{\infty }\frac{\text{e}^{-y}}{y}\,\text{d}y\end{eqnarray}$$

for $x>0$ . Using the large-argument expansion $\text{Ei}(x)=\text{e}^{x}/(x-1)+O(\text{e}^{x}x^{-3})$ , we can write (4.2) as

(4.4) $$\begin{eqnarray}t_{\ast }=\frac{3}{4{\mathcal{E}}}\frac{v_{0}^{2/3}\unicode[STIX]{x1D703}_{m}^{4/3}}{\ln {\displaystyle \frac{\tilde{\unicode[STIX]{x1D6FD}}v_{0}^{1/3}\unicode[STIX]{x1D703}_{m}^{2/3}}{\unicode[STIX]{x1D706}\text{e}^{2}}}}.\end{eqnarray}$$

Although ${\mathcal{K}}$ does not appear explicitly in (4.4), it is clear from our earlier discussion that ${\mathcal{K}}$ together with ${\mathcal{E}}$ and $\unicode[STIX]{x1D706}$ influence the values of both $\unicode[STIX]{x1D703}_{m}$ and $\tilde{\unicode[STIX]{x1D6FD}}$ which appear in (4.4) above.

Figure 10. (a) Evaporation times for droplets with $K=14$ , $B=0$ , $\unicode[STIX]{x1D706}=2\times 10^{-5}$ and $v(0)=10^{n/3}$ , $r(0)=2\unicode[STIX]{x1D703}_{m}^{-1/3}[v(0)]^{1/3}$ for $n=0$ , $\pm 1$ , $\pm 2$ , $\pm 3$ , when $E=4$ (black circles) and $E=24$ (grey circles), as predicted by the solution to the full problem. The solid curves correspond to the predictions of (4.4) for each set of parameters. (b,c) Plots of $v$ and $r$ as functions of $t_{\ast }-t$ for each solution to the full problem from (a), respectively. For each set of parameters (grey curves for $E=24$ and black curves for $E=4$ ), the data ultimately collapse onto the same dotted curve, which corresponds to the theoretical prediction of (4.5); the dashed curves in (b) and (c) correspond to the classical power laws $v\sim (t_{\ast }-t)^{3/2}$ and $r\sim (t_{\ast }-t)^{1/2}$ , respectively.

Figure 10 shows the dependence of $t_{\ast }$ on the initial droplet volume for two different sets of parameters. The values of $t_{\ast }$ obtained from solutions to the full problem are in excellent agreement with the values predicted with (4.4). Unlike the typical scaling $t_{\ast }\sim v_{0}^{2/3}$ , which can be deduced by simple arguments, e.g. by assuming the integral evaporation flux to be proportional to the droplet radius, in this set of parameters we have found better agreement with the scaling $t_{\ast }\sim v_{0}^{0.63}$ for both sets of calculations (note that other exponents can be obtained for different values of the system parameters). This result highlights the importance of the presence of the logarithmic terms in (4.4) in modulating the power law. An important consequence of the calculation presented above is that the evolution of $v(t)$ and $r(t)$ during the evaporation stage can be approximately predicted using the implicit formulae

(4.5a,b ) $$\begin{eqnarray}t_{\ast }-t=\frac{3}{4{\mathcal{E}}}\frac{v^{2/3}\unicode[STIX]{x1D703}_{m}^{4/3}}{\ln {\displaystyle \frac{\tilde{\unicode[STIX]{x1D6FD}}v^{1/3}\unicode[STIX]{x1D703}_{m}^{2/3}}{\unicode[STIX]{x1D706}\text{e}^{2}}}}\quad \text{and}\quad t_{\ast }-t=\frac{3}{16{\mathcal{E}}}\frac{r^{2}\unicode[STIX]{x1D703}_{m}^{2}}{\ln {\displaystyle \frac{\tilde{\unicode[STIX]{x1D6FD}}r\unicode[STIX]{x1D703}_{m}}{2\unicode[STIX]{x1D706}\text{e}^{2}}}}.\end{eqnarray}$$

The above predictions are readily verified in figures 10(b) and 10(c), where we plot $v$ and $r$ , respectively, as functions of $t_{\ast }-t$ for all 14 solutions of the full problem (solid curves), which were performed in preparation of figure 10(a). Looking at figure 10(b), all data from the solutions to the full problem collapse onto two (dotted) curves, corresponding to the theoretical result, (4.5a ), for the two different sets of parameters under consideration. The same happens when we plot $r$ as a function of $t_{\ast }-t$ (see figure 10 c), apart from the brief initial transients in the first two stages of the dynamics, which appear as nearly vertical lines in the plot. For comparison, we included in figure 10(b,c) the simple power-law scalings, $v\sim (t_{\ast }-t)^{3/2}$ and $r\sim \sqrt{t_{\ast }-t}$ , respectively, to demonstrate that they do not always predict the behaviour of the system accurately. In fact better fit to the data with the radius is provided by power laws with an exponent of approximately 0.53, i.e. slightly larger than $1/2$ . Interestingly, this is in qualitative agreement with what has been reported in the literature for diffusion-limited evaporation of water droplets (Cazabat & Guéna Reference Cazabat and Guéna2010). However, as will be discussed later, our model predicts such apparent exponents to be always larger than $1/2$ , which does not match experiments with other liquids. This is not surprising given the different evaporation regime (limited mostly by heat transfer) considered here.

Figure 11. (a) Influence of $K$ on the evaporation time, for droplets with $\unicode[STIX]{x1D706}=10^{-4}$ , $B=0$ , $v(0)=r(0)=1$ for two different values of $E$ , namely $E=1$ (black) and $E=8$ (grey). The solid curves correspond to the theoretical prediction, equation (4.4); the circles represent the results obtained by solving the full problem when $K=10^{-2n/3}$ , $n=0$ , $\pm 1$ , $\pm 2$ , $\pm 3$ . (b) Influence of $E$ on the evaporation time, for droplets with $\unicode[STIX]{x1D706}=10^{-4}$ , $B=0$ , $K=v(0)=r(0)=1$ . The solid curve corresponds to (4.4); the circles to solutions of the full problem for $E=10^{n/3}$ , $n=0,1,2,\ldots ,9$ .

In figure 11(a) we show typical plots demonstrating the dependence of $t_{\ast }$ on $K$ . We find that $t_{\ast }$ increases with $K$ provided that $K$ is large enough. This dependence, however, is not monotonic and there exists a moderate value of $K$ for which $t_{\ast }$ attains a minimum. This can be attributed to the different mechanisms that operate at disparate values of $K$ . The increase of $t_{\ast }$ for large values of $K$ is to be expected, since the evaporative flux diminishes as $K$ increases. For small $K$ , decreasing $K$ towards $0$ increases $\unicode[STIX]{x1D703}_{m}$ (and the apparent contact angle as a result), which in turn decreases the evaporative flux, given that the droplet becomes thicker with a smaller evaporation surface and a smaller contact area with the substrate. Nevertheless, these effects are less dramatic for larger values of $E$ , as the comparison of the two sets of calculations in figure 11(a) reveals. The dependence of $t_{\ast }$ on $E$ is more straightforward to describe, since, as expected, we always have a monotonic decrease of the evaporation time with $E$ (see figure 11 b). We should note that although $E$ appears explicitly at the denominator of $t_{\ast }$ , see (4.4), $\unicode[STIX]{x1D703}_{m}$ and $\tilde{\unicode[STIX]{x1D6FD}}$ also depend on $E$ , giving rise to non-trivial power-law dependencies (best fit to the data in figure 11 b was achieved with $t_{\ast }\sim E^{-0.69}$ ).

Figure 12. Comparison of the asymptotic predictions with the solution to the full problem for the evolution of $r(t)$ (solid curves) with $\unicode[STIX]{x1D706}=2\times 10^{-5}$ , $E=40$ , $B=0$ , $v(0)=r(0)=1$ and for $K=14$ , 140, 1400 and 14 000. The dashed curves show solutions to (3.21) with (3.38); the dotted curves show the solutions to (3.21) with (3.23). As $K$ increases, the agreement of the dashed curves with the full solution degrades; the opposite happens with the dotted curves.

4.3 The large- $K$ limit

The large- $K$ limit can be of relevance if $f_{a}$ is rather small for the reasons outlined in § 2. For example, some experimental studies with water reported values of $O(10^{-3})$ for $f_{a}$ (see Marek & Straub Reference Marek and Straub2001, and the references therein). In this limit, we expect that the evolution of $v(t)$ will not be well described by (3.38) because the analysis ceases to work for relatively large radii, see (3.41).

At the same time, one can also infer from figure 6(a) and the discussion about $q_{1}$ that simply approximating the inner region as a wedge of slope $\unicode[STIX]{x1D703}_{m}$ suffices to capture well the inner-region contributions to  $\dot{v}$ . This, in turn, means that utilising the simpler expression for $\dot{v}$ , equation (3.22), or, when $B=0$ , the exact expression (3.23), will provide a better approximation for $\dot{v}$ , noting that it does not suffer from the cutoff radius limitation of (3.38). By following similar arguments as in the derivation of (4.4) we obtain an estimate for $t_{\ast }$ that pertains in this limit, namely

(4.6) $$\begin{eqnarray}t_{\ast }=\frac{3v_{0}^{2/3}\unicode[STIX]{x1D703}_{m}^{4/3}}{4{\mathcal{E}}\ln \displaystyle \frac{v_{0}^{1/3}\unicode[STIX]{x1D703}_{m}^{2/3}}{{\mathcal{K}}\text{e}^{1/2}}},\end{eqnarray}$$

where ${\mathcal{K}}$ appears explicitly, unlike (4.4) where the dependence on $K$ is implicit.

Sample computations with large values of $K$ are shown in figure 12 for water-like parameters (see also figure 8). The evolution of $v(t)$ is not shown here because the differences are more subtle compared to those seen in the evolution of $r(t)$ . Here we compare the solutions to the full problem (solid curves) and solutions to the reduced problems using (3.21) either with (3.38) (dashed curves) or with (3.23) (dotted curves) for various values of $K$ . We demonstrate that indeed using (3.38) is not appropriate for the larger values of $K$ because $r_{c}$ becomes unacceptably large. In contrast, equation (3.23) appears to perform well throughout, albeit with less favourable agreement for moderately large values of  $K$ , particularly as $t\rightarrow t_{\ast }$ (see, e.g. the calculation when $K=14$ ). Note also that the non-monotonic variation of $t_{\ast }$ with $K$ (the simulation with $K=14$ requires a longer time for the droplet to evaporate compared to $K=140$ and $K=1400$ ) is consistent with what has been observed in figure 11(a).

Figure 13. Comparison of pendant (grey curves) and sessile (black curves) droplet dynamics when $|B|=3$ ; the rest of the parameters are as in figure 8. (a) Snapshots of pendant (left) and sessile (right) droplets at different times. Top to bottom profiles correspond to $t=10$ , 500, 1000, 1500 and 1800, respectively; arrows indicate the direction of gravity. (b,c) Evolution of the droplet radius and volume, respectively. Solid curves show the solutions to the full problem; the solutions to the reduced problem (dashed curves) are indistinguishable from those to the full problem. For comparison, the solution when $B=0$ is also shown (dotted curves).

4.4 The influence of gravity; pendant droplets

As previously mentioned, the effect of gravity diminishes as the droplet shrinks in size due to evaporation. Although the qualitative features of the dynamics remains unaltered, a number of important observations can be made in relation to the effect of gravity on the evaporation time. Figure 13 summarises the results of two calculations using the same parameters and initial condition as in figure 8, but with non-zero Bond number. More specifically, we study the case when $B=3$ (sessile droplet) and $B=-3$ (pendant droplet) and the results are contrasted with the calculation of figure 8, where $B=0$ . First, we note the excellent agreement of our theory with the solutions to the full problems, confirming that our analysis is applicable when $B\neq 0$ as well. Second, we see from these plots that the sessile droplet evaporates faster than the pendant one, with the droplet in the zero-gravity case completely evaporating at some instant in time between the two. The reason for the enhanced evaporation for the sessile droplet can be readily deduced in figure 13(a), where we see that gravity flattens a sessile droplet, which, apart from increasing the contact area, also boosts the evaporation rate density, see (2.10). This ultimately increases the evaporation flux, particularly at early times, when gravity manifests itself more strongly. Exactly the opposite happens when considering pendant droplets due to the fattening of the droplet near its axis of symmetry (see left-hand side of figure 13 a).

If the evolution for the volume shown in figure 13(c) is plotted as a function of $t_{\ast }-t$ (see figure 14 a) we see that for gravitational effects to be negligible, namely for the curves to approach the dashed curve corresponding to the zero gravity case, a rather significant amount of time must elapse until the volume of the droplet becomes sufficiently small. Indeed, a similar observation can be made if we plot the evaporation time as a function of the initial volume when $B=\pm 3$ (see figure 14 b). For smaller volumes, the estimate given by (4.4), the dotted curve, is reasonable. However, for larger sessile droplets, there is a significant deviation of the evaporation time from the value predicted by (4.4). This deviation can be estimated by numerical means by following the same principles as in § 4.2 which were invoked in estimating $t_{\ast }$ when $B=0$ – see the dashed curves in figure 14(b). We should note here that for pendant droplets only four data points were collected by solving the full problem for volumes between 0.1 and 1 due to the previously mentioned critical value of $c\approx 3.83$ , beyond which pendant droplets cannot be suspended from the substrate and droplet detachment/breakup may occur.

Figure 14. (a) Plot of the evolution of the volume shown in figure 13(c) as a function of $t_{\ast }-t$ . (b) Plot of the evaporation time, $t_{\ast }$ , as a function of the initial volume when $B=3$ (black circles) and $B=-3$ (grey circles). The rest of the parameters are as in figure 13. The circles correspond to data extracted by solving the appropriate full problems. The dotted curve corresponds to the theoretical estimate, equation (4.4); the dashed curves in (b) are numerically computed evaporation times, using techniques which were discussed in § 4.2.

4.5 The influence of slip

From the preceding discussion, we saw that $\unicode[STIX]{x1D706}$ can influence the parameters of the inner-region asymptotics appearing in (3.21) and (3.38), namely $\unicode[STIX]{x1D703}_{m}$ , $\unicode[STIX]{x1D6FD}$ and $\tilde{\unicode[STIX]{x1D6FD}}$ if we keep the other parameters of the system fixed (i.e. for fixed ${\mathcal{E}}$ and ${\mathcal{K}}$ ). For example, in figure 15 we present the results of two computations where we keep all system parameters the same apart from the values of $\unicode[STIX]{x1D706}$ . We see that when $\unicode[STIX]{x1D706}=10^{-3}$ , the droplet spreads more and, as a consequence of the flattening of the droplet, total evaporation occurs faster, nearly twice as fast compared to the case when $\unicode[STIX]{x1D706}=10^{-5}$ . Figure 16(a) shows that $\unicode[STIX]{x1D703}_{m}$ decreases as $\unicode[STIX]{x1D706}$ increases for fixed ${\mathcal{E}}$ and ${\mathcal{K}}$ , thus explaining why in the computations of figure 15 the droplet with the larger $\unicode[STIX]{x1D706}$ spreads to a larger radius before entering the evaporation stage. Lastly, figure 16(b) shows solutions of the full problem (circles) compared with the predictions of (4.4) (solid curves). In all computations considered, we see that $t_{\ast }$ decreases nearly linearly as $\ln \unicode[STIX]{x1D706}$ increases with a slope which depends on the other parameters of the system, although for comparatively larger values of ${\mathcal{E}}$ the dependence of $t_{\ast }$ on $\unicode[STIX]{x1D706}$ becomes weaker.

Figure 15. (a) and (b) depict the evolution of the droplet radius and volume, respectively for two disparate values of the slip length, $\unicode[STIX]{x1D706}$ , $10^{-5}$ (grey curves) and $10^{-3}$ (black curves). The rest of the parameters are ${\mathcal{E}}=10^{-3}$ , ${\mathcal{K}}=10^{-4}$ , $B=0$ with $v(0)=r(0)=1$ . The solutions to the full problems are shown as solid curves and are indistinguishable from those of the reduced problems (plotted as dashed curves).

Figure 16. (a) Dependence of $\unicode[STIX]{x1D703}_{m}$ on $\unicode[STIX]{x1D706}$ and for three sets of pairs of values $({\mathcal{E}},{\mathcal{K}})$ : $(10^{-3},10^{-4})$ – curves marked by black circles; $(10^{-3},10^{-3})$ – curves marked by white circles; $(10^{-2},10^{-3})$ – curves marked by grey circles. (b) Dependence of $t_{\ast }$ on $\unicode[STIX]{x1D706}$ for the same set of parameters as in (a), $B=0$ and $v(0)=r(0)=1$ . The circles correspond to values of $t_{\ast }$ obtained from the full problem; the solid curves correspond to the estimate given by (4.4).

4.6 The influence of substrate wettability

The wettability of a substrate is characterised by the value of Young’s equilibrium contact angle, $\unicode[STIX]{x1D717}_{s}$ . It is thus of interest to explore its effect on the dynamics and, more specifically on the evaporation time. One may naturally expect that a lower value for $\unicode[STIX]{x1D717}_{s}$ would cause the droplet to spread more, thus increasing its contact area with the wall, so that, just as gravity, shorter evaporation times would be observed.

To demonstrate this, consider a calculation of a fluid with parameters $E=0.125$ , $K=14$ , $\unicode[STIX]{x1D706}=2\times 10^{-5}$ and a Young’s angle of $\unicode[STIX]{x1D717}_{s}=20^{\circ }$ , subject to the initial condition $r(0)=1.5$ and $v(0)=1$ . This set of parameters, denoted as set A, corresponds to parameter values which are close to those of FC-72. This calculation is repeated with all physical parameters kept the same, apart from $\unicode[STIX]{x1D717}_{s}$ which is taken to be $\unicode[STIX]{x1D717}_{s}=10^{\circ }$ (set B). To maintain the same scaling for the radius and keep the same dimensional initial volume as in set A, $E$ must be increased by 16 times, whereas $v(0)$ and $\unicode[STIX]{x1D706}$ must undergo a twofold increase, whilst the rest are kept the same (hence, $E=2$ , $v(0)=2$ and $\unicode[STIX]{x1D706}=4\times 10^{-5}$ for set B). This is because $E$ scales with $\unicode[STIX]{x1D717}_{s}^{-4}$ , whereas $\unicode[STIX]{x1D706}$ and $V_{0}$ scale with $\unicode[STIX]{x1D717}_{s}^{-1}$ if the characteristic length scale, $d$ , is to be kept constant. Hence, decreasing $\unicode[STIX]{x1D717}_{s}$ further to $5^{\circ }$ (set C), we need to take $E=32$ , $\unicode[STIX]{x1D706}=8\times 10^{-5}$ and $v(0)=4$ , with the remaining parameters kept the same as in set A. The results of these calculations are shown in figure 17(a). Since time scales differently across the three set of parameters (left-hand side of figure 17 a), a proper comparison of the dynamics in the three cases would be made by matching the time scales of each set. Hence, on the right-hand side of figure 17(a), we match the time scales of the sets A and C to that of set B and we see that the actual evaporation time is approximately 76 % higher in set A and 12 % lower in set C compared to the evaporation time of set B.

Figure 17. Influence of $\unicode[STIX]{x1D717}_{s}$ (set A: light grey curves $\unicode[STIX]{x1D717}_{s}=20^{\circ }$ ; set B: grey curves $\unicode[STIX]{x1D717}_{s}=10^{\circ }$ ; set C: black curves $\unicode[STIX]{x1D717}_{s}=5^{\circ }$ ). (a)  $E=0.125\times 16^{i}$ , $K=14$ , $\unicode[STIX]{x1D706}=2^{i+1}\times 10^{-5}$ , $r(0)=1.5$ , $v(0)=2^{i}$ where $i=0,1,2$ for $\unicode[STIX]{x1D717}_{s}=20^{\circ }$ , $10^{\circ }$ , $5^{\circ }$ , respectively, which corresponds to keeping the initial geometry and all physical parameters the same for all sets apart from $\unicode[STIX]{x1D717}_{s}$ . (b) As in (a) but with $E=1.5\times 16^{i}$ . (c) As in (a) but with $E=18\times 16^{i}$ . Left plots: $r$ as a function of the dimensionless time, $t$ ; right plots: $r$ as a function of the dimensionless time for set B, $t^{\prime }$ . In all plots solid and dashed curves correspond to the solutions to the full and reduced problems, respectively.

By increasing the values of $E$ in all three sets of parameters by 12 times and producing their corresponding plots in figure 17(b), we find that the curves for sets B and C nearly collapse on top of each other when time is rescaled, with the dimensional evaporation time in set B differing from those of sets A and C by about 17 % and 1 %, respectively (right-hand side of figure 17 b). A further increase in the values of $E$ in all three sets of parameters of figure 17(a) by 144 times (i.e. 12 times higher compared to the calculations of figure 17 b) yields the plots shown in figure 17(c). We readily observe that all curves nearly collapse on top of each other when time is rescaled appropriately as in panels (a) and (b) (see right-hand side of figure 17 c). In this case, we find that the dimensional evaporation time in set B differs from those of A and C by approximately 2 % and 0.1 %, respectively.

From these calculations, we see that for small values of $E$ there is a clear dependence of the dynamics on the value of $\unicode[STIX]{x1D717}_{s}$ , which aligns with the expectation that lower values for $\unicode[STIX]{x1D717}_{s}$ yield shorter evaporation times. However, we have also seen that for sufficiently large values of $E$ the dependence of the evaporation time on $\unicode[STIX]{x1D717}_{s}$ diminishes. This can be rationalised by using (4.4) to cast the dimensional evaporation time in terms of the parameters of the problem as

(4.7) $$\begin{eqnarray}\frac{3\unicode[STIX]{x1D70C}LV_{0}^{2/3}\unicode[STIX]{x1D717}_{m}^{4/3}}{4k\unicode[STIX]{x0394}\unicode[STIX]{x1D6E9}(2\unicode[STIX]{x03C0})^{2/3}}\left[\ln {\displaystyle \frac{\tilde{\unicode[STIX]{x1D6FD}}V_{0}^{1/3}\unicode[STIX]{x1D717}_{m}^{2/3}}{3b\text{e}^{2}(2\unicode[STIX]{x03C0})^{1/3}}}\right]^{-1},\end{eqnarray}$$

where $V_{0}$ is the initial dimensional droplet volume. First, we note that $\unicode[STIX]{x1D717}_{s}$ does not appear explicitly in the expression above. The presence of $\unicode[STIX]{x1D717}_{s}$ is implicit, however, due to its influence on $\unicode[STIX]{x1D717}_{m}$ . At the same time, we have seen in § 3.1.2 that this dependence diminishes in the strong evaporation limit ( $E\gg 1$ ), which allows us to conclude that under such conditions the evaporation time becomes practically independent of the wetting properties of the substrate, which effectively corresponds to the completely wetting regime. Indeed, this is a direct consequence of the fact that the inner region, where the influence of $\unicode[STIX]{x1D717}_{s}$ features more prominently, influences the dynamics through $\unicode[STIX]{x1D717}_{m}$ , whose value is dominated by the evaporation-induced contribution when $E\gg 1$ . Hence as the apparent contact angle becomes close to $\unicode[STIX]{x1D717}_{m}$ during the comparatively longer evaporation stage, the influence of Young’s angle, $\unicode[STIX]{x1D717}_{s}$ , becomes negligible.

4.7 Apparent power laws

The section concludes with some remarks on the apparent power-law behaviours characterising the evolution of the radius and volume of the droplet close to complete evaporation. In a typical experiment, it is observed that the time remaining until extinction is approximately proportional to the droplet radius raised to some power (see, e.g. Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000; Shahidzadeh-Bonn et al. Reference Shahidzadeh-Bonn, Rafaï, Azouni and Bonn2006; Cazabat & Guéna Reference Cazabat and Guéna2010, and the references therein). The theoretical prediction of these exponents, which sometimes deviate from the classical $r^{2}$ -law (i.e. when the time remaining until extinction is proportional to the radius squared), was identified as one of the ‘hot topics’ discussed in the review of Bonn et al. (Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009). As previously alluded to, the logarithmic dependencies which arose from our analysis also predict non-trivial (apparent) power laws. For example, if we consider $\text{d}(\ln t_{\ast })/\text{d}(\ln v_{0})$ using (4.4), we deduce that $t_{\ast }\sim v_{0}^{n}$ , where $n$ is given by

(4.8) $$\begin{eqnarray}n=\frac{2}{3}-\frac{1}{3}\left(\ln \frac{\tilde{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D703}_{m}^{2/3}v_{0}}{\unicode[STIX]{x1D706}\text{e}^{1/2}}\right)^{-1}\approx \frac{2}{3}-\frac{1}{3}\left(\ln \frac{\tilde{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D703}_{m}^{2/3}}{\unicode[STIX]{x1D706}\text{e}^{1/2}}\right)^{-1}.\end{eqnarray}$$

Clearly the ‘local’ exponent $n$ depends on $v_{0}$ itself, but we may neglect $v_{0}$ appearing inside the logarithm by invoking a key assumption of our analysis, namely that the droplet is macroscopically large, or, equivalently, $r\gg r_{c}$ . In all cases considered, this approximation works rather well. For example, the differences in the exponents obtained from fitting the data of figure 10(a) and the exponents computed from (4.8) were within  $O(10^{-4})$ .

Figure 18. Plots of $v$ as a function of $t_{\ast }-t$ for different sets of parameters. The solutions to the full problems and the corresponding reduced problems are shown by solid and dashed curves, respectively. In all plots, the grey curves correspond to the data obtained with ${\mathcal{E}}=8\times 10^{-5}$ , ${\mathcal{K}}=28\times 10^{-5}$ , $\unicode[STIX]{x1D706}=2\times 10^{-5}$ , $B=0$ and $r(0)=v(0)=1$ , also plotted in figure 8. The black curves depict calculations where only one parameter is varied compared to the grey curves: (a)  ${\mathcal{K}}=28\times 10^{-4}$ – the classical power laws are also shown for visual comparison; (b)  ${\mathcal{E}}=8\times 10^{-3}$ ; (c)  $\unicode[STIX]{x1D706}=2\times 10^{-4}$ ; (d)  $B=10$ . The dotted curves in (a–c) are based on power laws of the form $v\sim (t_{\ast }-t)^{a}$ where $a=1/n$ . The inset in (d) shows the volume evolution at early times.

From (4.8), we also see that for realistic values of the various parameters we always have $n<2/3$ (recall that $n=2/3$ corresponds to the classical $r^{2}$ -law). Nevertheless, having exponents $n>2/3$ is likely to become more feasible with the inclusion of additional effects which are neglected in the simple model of our study. Although a detailed parametric investigation of these apparent power laws is beyond the scope of this work, for the sake of completeness, we show a few representative calculations in figure 18 to reveal some qualitative trends about the dependence of the exponent $a$ in $v\sim (t_{\ast }-t)^{a}$ on the system parameters using the calculation in figure 8 as a reference system. Looking at figure 18(a), we find that larger values of ${\mathcal{K}}$ increase $a$ slightly, but they manifest themselves more strongly as the droplet shrinks, thus altering the power law as $t\rightarrow t_{\ast }$ . Indeed, one would expect a power law of the form $v\sim (t_{\ast }-t)^{3}$ in the kinetically limited regime, where the global evaporation rate is proportional to the surface area of the droplet. This behaviour is observed here because using the larger value of ${\mathcal{K}}$ moves the onset of this regime to larger droplet sizes. From figure 18(b) we deduce that increasing ${\mathcal{E}}$ decreases $a$ , but the change is generally negligibly small. Moreover, as we have seen in the preceding section, increasing $\unicode[STIX]{x1D706}$ speeds up the evaporation process and as a result we see an increase in the prefactor (figure 18 c). In figures 18(a–c) we also plotted the apparent power law predicted by using $a=1/n$ where $n$ is given by the approximation (4.8). We readily see that the theoretical estimate predicts the overall trend quite well, provided that the volume of the droplet does not become too small. In figure 18(d) we explore the influence of gravity. Although its inclusion changes the exponent only at the onset (see the inset of figure 18 d), this change can yield appreciable differences in $t_{\ast }$ (for example, for the two simulations shown in figure 18 d, the evaporative dynamics is approximately 1.5 times slower when $B=0$ compared to $B=10$ ).

A final observation that can be made by looking at the curves of figure 18 is the clear breakdown of our analytical theory (the dashed curves) when the volume of the droplet becomes small, but this only occurs when typically more than 99 % of the fluid of the droplet evaporates. This departure of the theory from the numerics is strongly dependent on the parameters of the problem (compare, e.g. the simulations shown in figure 18 d), but ultimately, as expected from (3.40), the theory incorrectly predicts that the droplet stops evaporating when $v$ reaches a (small) value in finite time.