2.1. Hydrodynamics

The liquid is Newtonian with density $\rho$, viscosity $\mu$, and surface tension $\sigma '$. The characteristic vertical and lateral length scales of the droplet are given by $H$ and $L$, respectively. Their ratio, $\epsilon =H/L\ll 1$, is assumed small so that lubrication theory may be leveraged. The liquid motion is subject to no-slip and no-penetration conditions at the solid substrate, and stress balances at the liquid–gas interface.

We use the scale $p^*=\sigma _0\epsilon ^2/H$ for the pressure $p'$, $v_r^*=\epsilon ^3\sigma _0/\mu$ for the $r$ velocity $v_r'$, $v_z^*=\epsilon v_r^*$ for the $z$ velocity $v_z'$, and $t^*=H/v_r^*\epsilon$ for the time $t'$. These scales are derived from balancing terms in the mass and momentum balances (Larsson & Kumar Reference Larsson and Kumar2021). The quantity $\sigma _0$ is the surface tension at the saturation temperature $T_s$ and is used to scale $\sigma '$. Finally, for the evaporative flux $J'$, we use a general scale $J^*$ (defined in § 2.2) and non-dimensionalize the temperature $T'$ as $T=(T'-T_{s})/(T_b-T_s)$, where $T_b>T_s$ is the temperature of the heated substrate. For the remainder of this work, we use dimensionless variables (indicated without a prime superscript):

(2.1)\begin{equation} \left.\begin{gathered} r'=\epsilon^{-1}Hr,\quad z'=Hz,\quad h'=Hh,\quad T'=(T_b-T_s)T+T_{s},\\ v_r' = v_r^*v_r,\quad v_z'=\epsilon v_r^*v_z,\quad p'=p^*p,\quad t'=t^*t,\quad J'=J^*J. \end{gathered}\right\} \end{equation}

Surface tension $\sigma$ is assumed to vary linearly with the dimensionless interface temperature $T_I=T\vert _{z=h}$:

(2.2)\begin{equation} \sigma=1-\epsilon^2\,Ma\,T_I, \end{equation}

where the Marangoni number is $Ma=(\partial \sigma '/\partial T' )(T_b-T_s)/\sigma _0\epsilon ^2$.

Under the lubrication approximation, the Navier–Stokes equations become

(2.3)$$\begin{gather} \frac{\partial^2 v_r}{{\partial z}^2}=\frac{\partial p}{\partial r}, \end{gather}$$

(2.4)$$\begin{gather}\frac{\partial p}{\partial z}=0, \end{gather}$$

(2.5)$$\begin{gather}\frac{1}{r}\,\frac{\partial}{\partial r}(rv_r)+\frac{\partial v_z}{\partial z}=0, \end{gather}$$

with the boundary conditions $v_r=v_z=0$ at $z=0$, and

(2.6)$$\begin{gather} p=-\frac{1}{r}\,\frac{\partial}{\partial r}\left( r\,\frac{\partial h}{\partial r} \right)-\varPi, \end{gather}$$

(2.7)$$\begin{gather}\frac{\partial v_r}{\partial z}=-Ma\,\frac{\partial T_I}{\partial r}, \end{gather}$$

at $z=h$, where $\varPi$ is the disjoining pressure.

As stated, the no-slip condition imposed at the substrate does not allow for contact-line motion. There are several approaches to ameliorating this, and in this work we use a precursor-film approach because the resulting equations are much easier to solve (Savva & Kalliadasis Reference Savva and Kalliadasis2011; Kumar & Charitatos Reference Kumar and Charitatos2022). It is assumed that the substrate is covered with a thin precursor film and the liquid pressure is modified with a disjoining pressure to account for van der Waals interactions. As in previous work (Schwartz Reference Schwartz1998; Murisic & Kondic Reference Murisic and Kondic2011; Espín & Kumar Reference Espín and Kumar2017; Pham & Kumar Reference Pham and Kumar2017, Reference Pham and Kumar2019; Charitatos & Kumar Reference Charitatos and Kumar2020), we use a two-term disjoining pressure that includes both attractive and repulsive interactions,

(2.8)\begin{equation} \varPi(h)=A_1\left( \left( \frac{A_2}{h} \right)^3-\left( \frac{A_2}{h} \right)^2 \right), \end{equation}

where $A_1$ is a dimensionless Hamaker constant, and $A_2$ controls the thickness at which the disjoining pressure vanishes. The constants $A_1$ and $A_2$ also control the scaled equilibrium contact angle $\theta _{eq}$ through the relation (Schwartz Reference Schwartz1998)

(2.9)\begin{equation} \theta_{eq}\approx\sqrt{A_1A_2}. \end{equation}

This is related to the lab-frame contact angle $\theta _{eq}'$ by $\theta _{eq}'\approx \epsilon \theta _{eq}$. The precursor-film thickness $h_p$ is determined by the condition $J(h_p)=0$, which will be discussed in §§ 2.3.1 and 2.3.2. We have found that using a one-term attractive model for $\varPi$ (corresponding to a perfectly wetting liquid) results in smaller apparent contact angles, but does not change qualitatively the findings presented in this work (Ajaev Reference Ajaev2005; Eggers & Pismen Reference Eggers and Pismen2010; Maki & Kumar Reference Maki and Kumar2011; Karpitschka et al. Reference Karpitschka, Liebig and Riegler2017).

Solving (2.3) and (2.4) subject to boundary conditions (2.6) and (2.7) yields the radial velocity

(2.10)\begin{equation} v_r=\left( \frac{1}{2}\,z^2-hz \right)\frac{\partial p}{\partial r}-z\,Ma\,\frac{\partial T_I}{\partial r}. \end{equation}

Finally, mass conservation at the interface $z=h$ coupled with the continuity equation (2.5) gives the kinematic condition governing the film height:

(2.11)\begin{equation} \frac{\partial h}{\partial t}=-\frac{1}{r}\,\frac{\partial}{\partial r}\int_0^hrv_r\,\textrm{d}{z}-EJ=\frac{1}{3r}\,\frac{\partial}{\partial r} \left( rh^3\,\frac{\partial p}{\partial r}\right)+\frac{Ma}{2r}\,\frac{\partial}{\partial r}\left( rh^2\,\frac{\partial T_I}{\partial r}\right)-EJ, \end{equation}

where the pressure $p$ is given by (2.6), the evaporative number is $E=J^*/\epsilon \rho v_r^*$ (the ratio of evaporative and convective fluxes), and specific expressions for the evaporative flux $J$ will be given in § 2.3. This equation is subject to the symmetric boundary condition $h(r)=h(-r)$ as well as matching onto the precursor film $h\to h_p$ as $r\to \infty$.

2.2. Energy transport

The solid substrate is held at a constant dimensionless temperature $T=1$ on its bottom face at $z=-b$, where $b=b'/H$ is the thickness of the solid substrate. The solid and liquid have thermal conductivities $k_s$ and $k_\ell$, respectively. Because the droplet is thin, the temperature is governed by

(2.12)\begin{equation} \frac{\partial^2 T}{{\partial z}^2}=0 \end{equation}

in both phases, and subject to the boundary conditions

(2.13a,b)$$\begin{gather} \left.T\right\vert_{z=-b}=1,\quad \left.T\right\vert_{z=0^-}=\left.T\right\vert_{z=0^+}, \end{gather}$$

(2.14a,b)$$\begin{gather}\left.\frac{\partial T}{\partial z}\right\vert_{z=0^-}=\frac{\kappa}{b}\left.\frac{\partial T}{\partial z}\right\vert_{z=0^+},\quad -\left.\frac{\partial T}{\partial z}\right\vert_{z=h}=J, \end{gather}$$

where $\kappa =bk_\ell /k_s$ is a scaled thermal conductivity ratio, and we use the evaporative flux scale $J^*=\Delta T'\,k_\ell /L_m H$, where $L_m$ is the latent heat of vaporization.

Solving these equations yields the temperature field

(2.15)\begin{equation} T = \begin{cases} 1-\kappa J\left( 1+\dfrac{z}{b} \right), & -b \leq z < 0, \\ 1-J\left( z+\kappa \right), & 0 \leq z \leq h. \end{cases} \end{equation}

Of chief concern is the temperature at the liquid–gas interface,

(2.16)\begin{equation} T_I=\left.T\right\vert_{z=h}=1-J(h+\kappa), \end{equation}

which depends on the evaporative flux $J$, the film height $h$, and the conductivity ratio $\kappa$.

2.3. Evaporation

We now present expressions for the evaporative flux $J$ that appears in (2.11) and (2.16) for both one-sided and diffusion-limited evaporation.

2.3.1. One-sided evaporation

The one-sided evaporation model gives a constitutive relation for the evaporative flux $J$, which is (in non-dimensional form) (Burelbach et al. Reference Burelbach, Bankoff and Davis1988; Ajaev Reference Ajaev2005)

(2.17)\begin{equation} KJ=\delta p+T_I, \end{equation}

where

(2.18a,b)\begin{equation} K=\frac{k_\ell\sqrt{2{\rm \pi} R_gT_s^3}}{\rho_vL_m^2\alpha H},\quad \delta=\frac{p^*T_s}{\rho L_m\,\Delta T}. \end{equation}

Here, $R_g$ is the specific gas constant, $\rho _v$ is the density of the vapour, and $\alpha$ is the accommodation coefficient. The parameter $K$ gives a measure of kinetic effects and the volatility of the liquid, and can be interpreted as the inverse of an effective mass-transfer coefficient (Ajaev Reference Ajaev2005; Parrish & Kumar Reference Parrish and Kumar2020). The parameter $\delta$ measures the effect of changes in local pressure on the phase-change temperature (Ajaev Reference Ajaev2005). Substituting (2.16) into (2.17) yields the one-sided evaporative flux

(2.19)\begin{equation} J=\frac{1+\delta p}{K+h+\kappa}, \end{equation}

which is an explicit, local function of the film height $h$.

The precursor-film thickness $h_p$ is determined by the condition $J=0$ for a flat film. With the pressure given by (2.6), (2.19) reduces to finding the roots of a cubic polynomial:

(2.20)\begin{equation} 0=1+\delta p=1-\delta\,\varPi(h_p) \quad\Rightarrow\quad 0=h_p^3+\delta A_1A_2^2h_p-\delta A_1A_2^3. \end{equation}

The discriminant of this polynomial is negative, thus the single real root is

(2.21)\begin{equation} h_p=2A_2\sqrt{\frac{\delta A_1}{3}}\sinh\left[ \frac{1}{3}\ \text{asinh} \left(\frac{3}{2}\sqrt{\frac{3}{\delta A_1}} \right)\right], \end{equation}

where $\sinh$ is the hyperbolic sine function, and $\textrm {asinh}$ is its inverse. This relation is used to determine the precursor-film thickness in the initial condition to the numerical method described in § 3 when using one-sided evaporation.

2.3.2. Diffusion-limited evaporation

Under diffusion-limited evaporation, the evaporative flux $J$ is determined by the diffusive flux of vapour at the liquid–gas interface $z=h$. It has been shown that in many scenarios, transient effects as well as convection are negligible in the gas phase (Larson Reference Larson2014), so the transport problem reduces to pseudo-steady-state diffusion. The droplet interface is assumed to be in thermodynamic equilibrium, allowing us to relate the thermodynamic states of the liquid and gas phases.

Figure 2 depicts the problem for the dimensionless gas-phase concentration $c_g$ above the droplet. We use the vertical coordinate $z^*=z'/L$ since there is no apparent vertical scale (unlike inside the droplet), and $c_g$ is defined on the semi-infinite domain $z^*\in (\epsilon h, \infty )$ and $r\in (0,\infty )$. Because the droplet is thin and $\epsilon h\ll 1$, we can neglect the droplet height and instead use the constant domain $z^*\in (0,\infty )$ (shown formally in Appendix A). It is assumed that the vapour is at an ambient concentration $c_\infty$ far from the droplet, so we non-dimensionalize $c_g'$ as $c_g=(c_g'-c_\infty ')/c_{s}$, where $c_s$ is the concentration of saturated vapour. We then have $c_g\to 0$ as $r\to \infty$ or $z^*\to \infty$, as shown in figure 2.

Figure 2. Schematic of the vapour concentration $c_g$ problem in the semi-infinite domain $r\in (0,\infty )$ and $z^*\in (0,\infty )$. The thin liquid droplet (and precursor film) is in thermodynamic equilibrium with the vapour phase, so we use the Kelvin equation to describe $c_g$ at $z^*=0$.

The governing equation for $c_g$ is

(2.22)\begin{equation} 0=\nabla^2c_g=\frac{\partial^2 c_g}{{\partial z^*}^2}+\frac{1}{r}\,\frac{\partial}{\partial r}\left( r\,\frac{\partial c_g}{\partial r} \right). \end{equation}

For boundary conditions in $r$, we have the symmetric condition $c_g(r)=c_g(-r)$, which implies that $\partial c_g/\partial r$ and ${\partial ^3 c_g}/{\partial r^3}$ vanish at $r=0$. We also enforce matching to the ambient concentration $c_g\to 0$ as $r\to \infty$. In $z^*$, we match to the ambient concentration as $z^*\to \infty$ and then enforce thermodynamic equilibrium with the liquid droplet at $z^*=0$. For this, we use the Kelvin equation, which accounts for the effect of interface curvature and disjoining pressure on thermodynamic equilibrium (Eggers & Pismen Reference Eggers and Pismen2010):

(2.23)\begin{equation} \left.c_g\right\vert_{z^*=0}={\rm e}^{Zp}-c_\infty, \end{equation}

where $c_\infty = c_\infty '/c_s<1$, $Z=Mp^*/\rho R_g^* T$ is an effective compressibility factor, $M$ is the molecular weight of the liquid, and $R_g^*$ is the universal gas constant. If $Z=0$, then we obtain $c_g\vert _{z^*=0}=1-c_\infty$, or in dimensional terms, $c_g'\vert _{z^*=0}=c_s$, which is consistent with the classical diffusion-limited evaporation model. However, note that the term $Zp$ in (2.23) is never negligible in the precursor film due to large disjoining pressures (see (2.6) and (2.8)) that prevent the film from evaporating away. Later, we will adjust the values of $c_\infty$ and $Z$ to match the droplet lifetime with that from one-sided evaporation.

To obtain the evaporative flux $J$, we employ Fick's law so that

(2.24)\begin{equation} J=-\frac{J_D^*}{J^*}\left.\frac{\partial c_g}{\partial z^*}\right\vert_{z^*=0}, \end{equation}

where $J_D^*=c_s D\epsilon /H$ is the evaporative-flux scale arising from diffusion (opposed to $J^*$ from energetic considerations), and $D$ is the diffusivity of the vapour. For the purpose of this work, we assume the ratio $J_D^*/J^*$ to be near unity and neglect it, but in practice it may be absorbed into the evaporative number $E$ without loss of generality. A noteworthy result is the analytical solution for $J$ developed by Eggers & Pismen (Reference Eggers and Pismen2010); considering a problem analogous to that shown in figure 2, they showed that (in the current notation)

(2.25)\begin{equation} J(r)=-\int_0^\infty\mathcal{K}(r,r')\,\frac{\partial c_g\vert_{z^*=0}}{\partial r'}\,\textrm{d}{r'}, \end{equation}

where the kernel $\mathcal {K}$ involves the complete elliptic integrals of the first and second kind. As will be discussed in § 3, we have found numerical solution of (2.22) and (2.24) to be more stable in our nonlinear simulations, but the analytical results encapsulated in (2.25) offer valuable insight into the behaviour of $c_g$ and $J$, and will be leveraged in § 3. Furthermore, (2.25) demonstrates that $J$ is a global function of $h$ – distinct from the local function of $h$ obtained under one-sided evaporation (2.19).

The precursor-film thickness $h_p$ is determined by $J=0$, which is equivalent to $(\partial c_g/\partial z^*)\vert _{z^*=0}=0$ by (2.24). For a flat film at steady state (e.g. a precursor film), (2.22) implies that $c_g$ depends linearly on $z^*$. Since $c_g$ vanishes far from the film, and since the first derivative of $c_g$ vanishes at the film surface, $c_g=0$ everywhere for a flat film at steady state. Equation (2.23) then yields the cubic equation

(2.26)\begin{equation} 0=h_p^3+\frac{ZA_1A_2^2}{-\ln c_\infty}\,h_p-\frac{ZA_1A_2^3}{-\ln c_\infty}, \end{equation}

which is equivalent to (2.20) with $\delta$ replaced by $Z/(-\ln c_\infty )$ (note that $c_\infty <1$). Thus through the same method employed to obtain (2.21), we obtain

(2.27)\begin{equation} h_p=2A_2\sqrt{\frac{ZA_1}{-3\ln c_\infty}}\sinh\left[ \frac{1}{3}\ \text{asinh} \left( \frac{3}{2}\sqrt{-\frac{3\ln c_\infty}{ZA_1}} \right) \right]. \end{equation}

This relation is used to determine the precursor-film thickness in the initial condition to the numerical method described in § 3 when using diffusion-limited evaporation.

2.4. Parameter values

The effects of the key dimensionless groups shown in table 1 will be investigated in §§ 5 and 6. The remaining parameters – $A_1$ and $A_2$ for the disjoining pressure, $\delta$ and $K$ for one-sided evaporation, and $Z$ and $c_\infty$ for diffusion-limited evaporation – are fixed for the remainder of this work. The constants $A_1=10^2$ and $A_2=10^{-3}$ are chosen close to those used in previous work (Pham & Kumar Reference Pham and Kumar2017; Charitatos et al. Reference Charitatos, Pham and Kumar2021) and give a scaled equilibrium contact angle $\theta _{eq}\approx 18^\circ$ (see (2.9)). Following Murisic & Kondic (Reference Murisic and Kondic2011), we choose $K=10$ and $\delta =10^{-3}$ (determined by fitting predictions to experimental measurements of evaporating water droplets), giving $h_p\approx 3.9\times 10^{-4}$. As will be discussed in § 4, we choose $Z=7.69\times 10^{-5}$ and $c_\infty =9.26\times 10^{-1}$ to match the total evaporation time as well as $h_p$ under diffusion-limited evaporation. Note that at this small value of $Z$, the vapour concentration is very nearly constant over the droplet surface as discussed in § 2.3.2.

Table 1. Important dimensionless parameters and typical values.

3.1. One-sided evaporative flux

Under one-sided evaporation, $J$ given by (2.19) is an explicit, local function of the liquid thickness $h$ and is thus easily computed alongside the finite-difference method used to solve (2.11). There are no additional numerics required.

3.2. Diffusion-limited evaporative flux

Under diffusion-limited evaporation, $J$ is obtained from the solution of (2.22) and is thus an implicit, global function of $h$ (see (2.25)). When approached numerically, commonly (2.22) is solved with a finite-element method (Diddens Reference Diddens2017; Diddens et al. Reference Diddens, Kuerten, van der Geld and Wijshoff2017a; Pahlavan et al. Reference Pahlavan, Lisong, Bain and Stone2021; Thayyil Raju et al. Reference Thayyil Raju, Diddens, Li, Marin, van der Linden, Zhang and Lohse2022), but many analytical simplifications have been developed. Eggers & Pismen (Reference Eggers and Pismen2010) derived (2.25), which removes the necessity of solving a partial differential equation (PDE) at each time step. However, the kernel $\mathcal {K}$ is ill-conditioned and singular at $r'=r$, so computation of the integral is numerically complex and prone to error that can induce instability in the time-stepping algorithm. An alternative is to assume that the droplet is a spherical cap, in which case the evaporative flux is given by

(3.2)\begin{equation} J= \frac{2}{\rm \pi}\,\frac{1-c_\infty}{\sqrt{R^2-r^2}}, \end{equation}

which removes any numerical complications when computing the evaporative flux (note that a similar approximation can also used for non-thin droplets) (Wilson & D'Ambrosio Reference Wilson and D'Ambrosio2022). Equation (3.2) is simple and useful for gaining physical insight, so we will reference it when discussing our results. However, we seek to develop understanding for more general droplet shapes, so to compute $J$ in our nonlinear simulations, we solve (2.22) using the novel hybrid spectral finite-difference method detailed below.

Since the value of $h$ is required in the equilibrium boundary condition for $c_g$ (2.23), it is convenient to choose the same discretization for $r$ for both (2.11) and (2.22) (second-order centred finite differences). This necessitates a finite $r$ domain despite the boundary conditions for (2.22) as $r\to \infty$. It has been shown that, away from the contact line, $h$ decays to the precursor-film thickness as $1/r$ (Eggers & Pismen Reference Eggers and Pismen2010), so an $r$ domain that is significantly larger than the droplet radius is required to resolve the evaporative flux properly. We choose $r\in (0,100)$ because scaling relations (Eggers & Pismen Reference Eggers and Pismen2010), as well as numerical results from our simulations, show that this is adequate to resolve the evaporative flux within numerical tolerances.

We must also discretize (2.22) in the $z^*$ domain. It is possible to similarly choose a finite distance (Diddens Reference Diddens2017; Loussert et al. Reference Loussert, Doumenc, Salmon, Nikolayev and Guerrier2017) and solve in a truncated $z^*$ domain by finite elements or finite differences. However, unlike the $r$ domain, the proper distance in the $z^*$ domain is unclear. Furthermore, we are not restrained by boundary data as in the $r$ domain, so it is possible to use a different discretization that is capable of resolving semi-infinite domains. Thus we use a spectral method in $z^*$ with Laguerre basis functions on the domain $z^*\in (0,\infty )$. The Laguerre functions (Boyd Reference Boyd2001) $L_n(z^*)$ are a modification of the Laguerre polynomials (Abramowitz & Stegun Reference Abramowitz and Stegun1964) $\ell _n(z^*)$ to a semi-infinite domain:

(3.3)\begin{equation} L_n(z^*)=\exp({-z^*/2})\,\ell_n(z^*). \end{equation}

We choose these basis functions over other candidates because they decay exponentially, which matches the exponential decay of $c_g$ shown by Eggers & Pismen (Reference Eggers and Pismen2010).

This hybrid spectral finite-difference method represents the gas-phase concentration as the expansion

(3.4)\begin{equation} c_g(r,z^*,t)\approx\sum_{n=0}^{M-1}a_n(r,t)\,L_n(z^*), \end{equation}

where the coefficients $a_n$ depend on $r$ and $t$. Inserting this into (2.22) and applying orthogonality of the Laguerre functions gives the PDE

(3.5)\begin{equation} a_n''+\frac{1}{r}\,\frac{\partial}{\partial r}\left( r\,\frac{\partial a_n}{\partial r} \right)=0, \end{equation}

for each $n$, where $a_n''(a_1,\ldots,a_{n})$ gives the coefficients of the second $z^*$ derivative of $c_g$. This PDE can be discretized further in $r$ with the same finite-difference method, and solved for $a_n$ alongside (2.11). The evaporative flux is then obtained from the relation

(3.6)\begin{equation} J(r,t)=-\sum_{n=0}^{M-1}a_n'(r,t), \end{equation}

which follows from the identity $L_n(0)=1$. Here, $a_n'$ are the coefficients of the first $z^*$ derivative of $c_g$.

The end result of this numerical method is the reduction of (2.22) to an $NM\times NM$ banded system of equations with bandwidth $4M$, where $M$ is the number of Laguerre functions, and $N$ is the number of finite-difference cells, that is solved efficiently by standard banded-system solvers. This requires computational time and memory similar to a full finite-difference/element method, but is capable of resolving the entire semi-infinite domain with spectral accuracy in $z^*$. To the best of our knowledge, this is the first numerical method to employ Laguerre functions in the context of evaporation. See the supplementary material available at https://doi.org/10.1017/jfm.2023.873 for a detailed discussion of implementing this method.

For accuracy, we must have high resolution near the droplet contact line. However, this level of resolution is not necessary far from the droplet. Thus, near the droplet, we use cells with constant width, but those far from the droplet become wider. Formally, the cells at a radius less than a critical radius, $r< r_c$, have constant width, and each cell beyond $r_c$ is a factor $c$ wider than the previous one. For a given number of nodes $N$, the number of nodes inside the critical radius $N_i$ is given by

(3.7)\begin{equation} N_i=-\alpha+\frac{1}{\ln c}\,{W}(\alpha c^{N+\alpha}\ln c), \end{equation}

where $\alpha =cr_c/(c-1)(d-r_c)$, $d$ is the size of the $r$ domain, and $W$ is the principal branch of the Lambert-$W$ function (Corless et al. Reference Corless, Gonnet, Hare, Jeffrey and Knuth1996). In this work, we use $r_c=2$, $c=1.05$, $d=100$ and $N=4000$ nodes, giving $N_i\approx 3807$. This approach covers the entire domain $r\in (0,100)$ with only minor cell stretching, while placing over $95\,\%$ of the cells near or inside the droplet ($r<2$).

This hybrid spectral finite-difference method also allows computation of the vapour field $c_g$ above the droplet, whereas direct expressions for $J$ such as (2.25) and (3.2) provide no information about $c_g$. Figure 3 shows $c_g$ above a droplet with radius $R\approx 1.4$ (located at $z^*=0$) where we present $z^*$ on a logarithmic scale to show that the vapour concentration $c_g$ decays exponentially. Vapour is concentrated near the droplet at $z^*=0$ and approaches the scaled ambient value $c_g=0$ as $z^*\to \infty$ or $r\to \infty$. Note that establishing this vapour field is critical to the diffusion-limited evaporation model; processes such as convection can interfere and cause the diffusion-limited model to become inaccurate (Shahidzadeh-Bonn et al. Reference Shahidzadeh-Bonn, Rafaï, Azouni and Bonn2006). As will be discussed in § 5, one-sided evaporation may offer a more accurate description in such a scenario.

Figure 3. (a) Example of the vapour concentration $c_g$ above the droplet (located at $z^*=0$). (b) Plots of the vapour concentration $c_g$ at fixed $r$. Straight dashed lines are provided for reference to show that $c_g$ decays exponentially ($\ln c_g\sim z^*$) for sufficiently large $z^*$.

5.1. Comparison of droplet radius

We now compare the two evaporation models to understand the differences in the predictions of the droplet radius. Turning first to diffusion-limited evaporation, the droplet radius shown in figure 5(a) (black lines) is well-approximated by (4.5) because the droplet shape is well-approximated by a spherical cap. Note that this results in (4.4) for the total evaporative flux $J_{tot}\sim R$ which scales with the droplet circumference. This result holds regardless of the parameters $c_\infty$ and $Z$ as long as the droplet shape is well-approximated by a spherical cap. For other situations, such as those where gravity is important, this approximation may not hold. Note that the nonlinear simulations are not restricted to the case of spherical-cap droplets.

Now focusing on one-sided evaporation and assuming a spherical-cap shape (4.2), we may rewrite $J$ given by (2.19) in terms of $\varGamma$ (neglecting the pressure contribution because it is unimportant in the bulk):

(5.1)\begin{equation} J=\frac{\left(K+\kappa\right)^{-1}}{1+\varGamma\left( 1-\left( \dfrac{r}{R} \right)^2 \right)}. \end{equation}

In the limit that $\varGamma =R\tan \theta _a/2(K+\kappa )\ll 1$, this reduces to the constant evaporative flux $J\approx ( K+\kappa )^{-1}$. The total evaporative flux then becomes

(5.2)\begin{equation} J_{tot}=\int_0^R rJ\,\textrm{d}{r}\approx\frac{R^2}{2(K+\kappa)}. \end{equation}

The condition $\varGamma \ll 1$ states that the droplet shape does not influence the evaporative flux, hence giving a nearly constant $J$ throughout the bulk. In contrast to the simple analytical relations that one obtains for diffusion-limited evaporation (4.4), (5.2) is valid only in parameter regimes where $\varGamma \ll 1$.

With $R\sim 1$, $\theta _a\approx 23^\circ$ and $K+\kappa =10$, we have $\varGamma \ll 1$ for results presented in this work – as evidenced by the nearly constant one-sided evaporative flux in figure 4(b). Substituting (5.2) into (4.1) results in a linear scaling for the droplet radius under one-sided evaporation:

(5.3)\begin{equation} R(t)=\frac{4E}{3\tan\theta_a(K+\kappa)}\,(t_0^{(o)}-t). \end{equation}

Note that this assumes a constant apparent contact angle $\theta _a$, and despite $\theta _a$ varying as shown in figure 5(b), well-approximates the radius evolution shown in figure 5(a) (red lines). Equivalently, (5.3) can be obtained by simplifying (4.8) in the limit $\varGamma \ll 1$. We note that a similar linear relationship can be derived for non-thin droplets undergoing constant evaporation (Masoud & Felske Reference Masoud and Felske2009).

In experimental studies of drying droplets, it has been observed that $J_{tot}$ is proportional to the surface area (${\sim }R^2$) in some systems (such as forced convection over the droplet) (Shahidzadeh-Bonn et al. Reference Shahidzadeh-Bonn, Rafaï, Azouni and Bonn2006), and diffusion-limited evaporation fails to capture the evolution because it predicts $J_{tot}$ proportional to the circumference of the droplet ($J_{tot}\sim R$ as shown in (4.4)). In such a scenario, one-sided evaporation may be able to capture the droplet evolution since it predicts $J_{tot}\sim R^2$ when $\varGamma \ll 1$; by treating the value of $K$ in relations (5.2) and (5.3) as a fitting parameter, one can adjust the evaporation rate to fit experimental data.

5.2. Comparison of droplet contact angle

We now discuss trends in the apparent contact angle $\theta _a$ and elucidate fundamental differences in the structure of the two evaporation models near the droplet contact line. Figure 5(b) shows that after the initial spreading phase, $\theta _a$ remains nearly constant under diffusion-limited evaporation (black line) but decreases steadily under one-sided evaporation (red line). To understand this difference, we must look more closely at the dynamics near the droplet contact line that give rise to $\theta _a$. Figure 6 shows a scaled-up view of the droplet contact line, where capillary effects tend to spread the droplet while evaporation tends to contract it. Eggers & Pismen (Reference Eggers and Pismen2010) showed that under diffusion-limited evaporation, the size of the contact-line region $w$ scales as $w\sim R^{1/3}$, so $\xi$ depicted in figure 6 is a scaled distance away from the contact line such that $\xi \sim 1$ in the contact-line region:

(5.4)\begin{equation} \xi\sim\frac{R-r}{R^{1/3}}. \end{equation}

Figure 6. Schematic of droplet contact-line region. Evaporation competes with capillary flow to establish the apparent contact angle $\theta _a$. The coordinate $\xi$ is a scaled distance from the droplet contact line given by (5.4). The bulk pressure is given by (5.5).

The driving force for spreading arises from the difference in capillary pressure in the bulk and disjoining pressure in the precursor film. For a spherical-cap droplet, substituting (4.2) into (2.6) (and neglecting disjoining pressure in the bulk) gives

(5.5)\begin{equation} p\approx-\frac{1}{r}\,\frac{\partial}{\partial r}\left( r\,\frac{\partial h}{\partial r} \right)=\frac{2\tan\theta_a}{R(t)}. \end{equation}

While the pressure in the precursor film $p=-\varPi (h_p)$ remains constant and negative, the bulk pressure scales as $p\sim R^{-1}$, which grows as the droplet shrinks. Thus the driving force for spreading becomes larger as the droplet evaporates, which tends to decrease $\theta _a$ over time.

Opposing this capillary spreading is evaporation, where a non-uniform evaporative flux that is focused at the contact line will tend to increase $\theta _a$. Under diffusion-limited evaporation (assuming a spherical cap for simplicity), we can use (3.2) and (5.4) to show that near the contact line

(5.6)\begin{equation} J\sim\frac{1}{\sqrt{R^2-r^2}}\sim\frac{1}{\sqrt{R(R-r)}}\sim\frac{1}{R^{2/3}\sqrt{\xi}}. \end{equation}

First, we see that $J$ grows like the classical square-root singularity (Eggers & Pismen Reference Eggers and Pismen2010) $\xi ^{-1/2}$ as we move towards the contact line ($\xi \to 0$), obtaining its maximum (singular) value at $r=R$. Second, note that as the droplet shrinks and $R$ decreases, $J$ increases and becomes more non-uniform (the peak shown in figure 4(b) becomes steeper). Thus diffusion-limited evaporation tends to increase $\theta _a$ as the droplet shrinks, balancing the increasing capillary pressure and giving the nearly constant $\theta _a$ shown in figure 5(b) (black line).

For one-sided evaporation, note that near the contact line we have $h\approx h_p$, and the pressure is dominated by the disjoining pressure, so we have from (2.19) that

(5.7)\begin{equation} J=\frac{1+\delta p}{K+h+\kappa}\approx\frac{1}{K+h_p+\kappa}\left( 1-\delta\,\varPi(h) \right). \end{equation}

Thus, as we move towards the contact line and precursor film where $1-\delta \,\varPi (h)\to 0$ (see (2.20)), the evaporative flux $J$ decreases (shown in figure 4b). This behaviour is opposite to that for diffusion-limited evaporation, resulting in the steady decrease in $\theta _a$ shown in figure 5(b) (red line). This may explain the decreasing contact angle observed experimentally by Murisic & Kondic (Reference Murisic and Kondic2011) for evaporating water droplets well-described by one-sided evaporation.

These effects are amplified as the evaporative number $E$ is increased. Figure 7 shows the apparent contact angle $\theta _a$ over time for both one-sided and diffusion-limited evaporation as $E$ is increased. With slow evaporation, the droplet will assume its equilibrium shape, and the contact angle will remain nearly constant at $\theta _{eq}$ under both evaporation models (though the radius evolution remains qualitatively different, as shown in figure 5a). This is what we observe for $E=10^{-3}$ in figure 7 (black lines). Now, as $E$ is increased under diffusion-limited evaporation, $\theta _a$ increases but remains nearly constant throughout most of the droplet lifetime. As $E$ is increased, the effect of evaporation increasing $\theta _a$ is intensified (see figure 6), leading to a larger $\theta _a$. This in turn increases the bulk capillary pressure (5.5), maintaining the balance between capillary spreading and non-uniform evaporation. However, one-sided evaporation tends to decrease $\theta _a$. So while increasing $E$ results in a net increase in $\theta _a$ due to faster overall evaporation, there is a steeper decrease in $\theta _a$ over time due to increasing capillary pressure.

Figure 7. The apparent contact angle $\theta _a$ over time for (a) diffusion-limited and (b) one-sided evaporation as $E$ is increased by an order of magnitude for each line. The black line is for $E=10^{-3}$. Note that $\theta _a$ approaches the equilibrium contact angle $\theta _{eq}\approx 18^\circ$ (see § 2.4) as $E$ is decreased because the droplet approaches its equilibrium shape.