2.2.1. Gradient dynamics form

Having set the stage for formulating mesoscopic models as gradient dynamics on an underlying energy functional, we next introduce such a model for an evaporating sessile liquid drop with profile $h(\boldsymbol {r},t)$ in a gap of height $d$ (see figure 1). We employ two fields, on the one hand, the amount $\psi _1(\boldsymbol {r},t)$ of the substance in liquid state in the drop per substrate area and, on the other hand, the amount $\psi _2(\boldsymbol {r},t)$ of the substance in vapour form in the gas phase also per substrate area. The field $\psi _1(\boldsymbol {r},t)$ is proportional to the thickness of the liquid film, i.e.

(2.4)\begin{equation} \psi_1(\boldsymbol{r},t)= \rho_{liq} h(\boldsymbol{r},t),\end{equation}

where $\rho _{liq}$ is the constant liquid density (to be specified later). All employed densities are number densities, i.e. are given in units of particles per volume. The field $\psi _2(\boldsymbol {r},t)$ is proportional to the height-averaged vapour density $\rho (\boldsymbol {r},t)$, namely,

(2.5)\begin{equation} \psi_2(\boldsymbol{r},t)= \rho(\boldsymbol{r},t) [d-h(\boldsymbol{r},t)]. \end{equation}

The gas phase can either consist of pure vapour or of vapour in an inert gas (called ‘air’ in the following). The latter has height-averaged density $\rho _{air}(\boldsymbol {r},t)$ and the resulting total gas density is $\rho _{tot}=\rho +\rho _{air}$. The literature distinguishes the following three main cases.

1. (i) Evaporation into pure vapour treated isothermally. In this case no diffusion occurs, all dynamics in the gas phase is due to pressure equilibration via convective motion. However, this is a very fast process as it occurs with the speed of sound (Maxwell Reference Maxwell1890). On the time scale of evaporation one then assumes uniform vapour, i.e. gas pressure. This implies that $\rho (\boldsymbol {r},t)$ is uniform.

2. (ii) Evaporation into air treated isothermally. In this case vapour diffusion is very important. As in case (i), the total pressure equilibrates fast, i.e. $\rho _{tot}$ is constant and uniform, e.g. $p=k_B T \rho _{tot}$ for an ideal gas, where $k_B$ is the Boltzmann constant and $T$ the temperature. In consequence,

   (2.6)\begin{equation} \rho_{air}(\boldsymbol{r},t)=\rho_{tot}-\rho(\boldsymbol{r},t) = \rho_{tot}-\frac{\psi_2(\boldsymbol{r},t)}{d-h(\boldsymbol{r},t)}. \end{equation}

3. (iii) As case (i) or case (ii), but not treated isothermally. Then, heat diffusion becomes a limiting factor as the heat used up as latent heat during the liquid–vapour phase transition needs to be transported to the interface. It also becomes important that normally a jump in the heat flux is considered at the interface that ultimately controls the evaporation flux. Here, we will not discuss this case.

Next, we write the model for the coupled dynamics of the local amounts of liquid and vapour in the gradient dynamics form (2.3). In particular, the $\psi _i$ follow the mixed conserved and non-conserved dynamics

(2.7) \begin{equation} \left.\begin{gathered} \partial_t \psi_1 = \boldsymbol{\nabla}\, \boldsymbol{\cdot}\, \left({\mathsf{Q}}_{11}\boldsymbol{\nabla}\frac{\delta \mathcal{F}}{\delta \psi_1} +{\mathsf{Q}}_{12}\boldsymbol{\nabla}\frac{\delta \mathcal{F}}{\delta \psi_2}\right)- {\mathsf{M}}_{11}\frac{\delta \mathcal{F}}{\delta \psi_1} - {\mathsf{M}}_{12} \frac{\delta \mathcal{F}}{\delta \psi_2},\\ \partial_t \psi_2 = \boldsymbol{\nabla}\, \boldsymbol{\cdot}\, \left({\mathsf{Q}}_{21}\boldsymbol{\nabla}\frac{\delta \mathcal{F}}{\delta \psi_1}+ {\mathsf{Q}}_{22}\boldsymbol{\nabla}\frac{\delta \mathcal{F}}{\delta \psi_2}\right) -{\mathsf{M}}_{21}\frac{\delta \mathcal{F}}{\delta \psi_1}-{\mathsf{M}}_{22}\frac{\delta \mathcal{F}}{\delta \psi_2}. \end{gathered}\right\}\end{equation}

Note, however, that $h$ is the relevant field for many of the physical effects and that often it may be more convenient to write the governing equations in terms of $h$ (2.4) and $\psi =\psi _2$.

First, we focus on case (ii) before discussing the amendments needed for case (i). As on the considered time scales, the total gas pressure is uniform, there is no pressure gradient that would drive convective transport in the gas phase. In consequence, there is no dynamic coupling between the liquid and gas layer, i.e. ${\mathsf{Q}}_{12}={\mathsf{Q}}_{2 1}=0$. Transport in the gas phase is then limited to vapour diffusion within the air. As a result, the conserved mobilities are

(2.8)\begin{equation} \boldsymbol{\mathsf{Q}}=\left(\begin{array}{cc} \dfrac{1}{\rho_{liq}}\,\dfrac{\psi_1^3}{3 \eta} & 0 \\[8pt] 0 & \tilde{D} \psi_2 \end{array}\right),\end{equation}

where $\eta$ is the dynamic viscosity of the liquid and $\tilde {D}$ is a diffusive mobility constant. We obtain ${\mathsf{Q}}_{11}$ by considering the standard long-wave evolution equation (Oron et al. Reference Oron, Davis and Bankoff1997; Craster & Matar Reference Craster and Matar2009). We multiply its gradient dynamics formulation (Thiele Reference Thiele2010) by the liquid density $\rho _{liq}$ and replace $h$ by $\psi _1/\rho _{liq}$ to obtain

(2.9)\begin{equation} \partial_t \psi_1 = \rho_{liq}\partial_t h = \rho_{liq} \boldsymbol{\nabla}\, \boldsymbol{\cdot}\, \left(\frac{h^3}{3 \eta}\boldsymbol{\nabla} \frac{\delta \mathcal{F}}{\delta h}\right) = \boldsymbol{\nabla}\, \boldsymbol{\cdot}\, \left({\mathsf{Q}}_{11} \boldsymbol{\nabla}\frac{\delta \mathcal{F}}{\delta \psi_1}\right). \end{equation}

The form of the diffusive mobility ${\mathsf{Q}}_{22}$ is for dilute solutions discussed by Thiele (Reference Thiele2011) and Xu, Thiele & Qian (Reference Xu, Thiele and Qian2015). The presently employed form $\tilde {D} \psi _2=\tilde {D} (d-h)\rho$ is the direct equivalent for the present case of vapour diffusion in the gap between the liquid and upper wall. It can be related to the usual diffusion constant $D$ of the vapour particles in air by $D = k_B T \tilde {D}$.

In the simplest case, the non-conserved mobilities are

(2.10)\begin{equation} \boldsymbol{\mathsf{M}} = \tilde{M}\left(\begin{array}{cc} 1 & -1 \\ -1 & 1 \end{array} \right),\end{equation}

where $\tilde {M}$ is an evaporation rate constant, which can be estimated, e.g. from the Hertz–Knudsen equation (Knudsen Reference Knudsen1915; Librovich et al. Reference Librovich, Nowakowski, Nicolleau and Michelitsch2017), which linearly incorporates an accommodation (or ‘sticking’) coefficient of the gas molecules. Note that more intricate models for phase change may be implemented via the non-conserved terms, e.g. by employing a variable evaporation rate $\tilde {M}(\psi _1, \psi _2)$. Our choice in (2.10) represents the simplest case for an evaporation that is purely driven by the differences in chemical potentials. Moreover, the particular matrix (2.10) ensures that the total particle number $\int (\psi _1+\psi _2)\,\mathrm{d}\kern0.7pt x\,\mathrm {d} y$ is conserved, i.e. $\psi _1+\psi _2$ satisfies the continuity equation $\partial _t (\psi _1+\psi _2)= \boldsymbol {\nabla }\boldsymbol {{\cdot }}(\,\boldsymbol {j}_{\psi _1}+\boldsymbol {j}_{\psi _2})$.

The free energy in long-wave approximation is

(2.11)\begin{equation} \mathcal{F} \!=\! \int_\varOmega \left[\underbrace{\gamma\left(1+\frac{|\boldsymbol{\nabla} h|^2}{2}\right)}_{\mathrm{surface}} +\! \underbrace{g(h)}_{\text{wetting}} +\!\underbrace{h f_{liq}(\rho_{liq})}_{\text{liquid bulk}} +\underbrace{(d\!-\!h) f_{vap}(\rho)}_{\text{vapour bulk}} +\underbrace{(d\!-\!h) f_{air}(\rho_{air})}_{\text{air bulk}} \right] \mathrm{d}\kern0.7pt x\,\mathrm{d} y, \end{equation}

where the $f_i$'s are bulk liquid, vapour and air energies per volume, $\gamma$ is the liquid–gas interface tension and $g(h)$ is the wetting energy per area. The functional is accompanied by the constraint of particle number conservation across the two phases. The condition is

(2.12)\begin{equation} 0 = \int_\varOmega \left[\psi_1 + \psi_2 - \bar{n}\right]\mathrm{d}\,\boldsymbol{x} = \int_\varOmega \left[h \rho_{liq} + (d-h) \rho -\bar{n}\right] \mathrm{d}\kern0.7pt x\,\mathrm{d} y, \end{equation}

where $\bar {n}$ is the mean (liquid and vapour) particle number per substrate area. Here, particle flux through the boundaries of the domain $\varOmega$ is excluded, however, it can be easily incorporated. Also, gravity is neglected as we consider small droplets, but may be added in the form of potential energy.

The presented long-wave formulation of the mesoscopic model is best suited for shallow droplets. However, as briefly discussed in § 3 of Thiele (Reference Thiele2018), it is an advantage of the gradient dynamics formulation that one may separately discuss improvements of the energy functional and the mobilities. There, it is argued that a good strategy for model improvements consists in making the energy as exact as possible and keeping the mobilities as simple as possible. Here, in particular, we may use a full-curvature trick similar to Gauglitz & Radke (Reference Gauglitz and Radke1988), Snoeijer et al. (Reference Snoeijer, Andreotti, Delon and Fermigier2007) and von Borries Lopes, Thiele & Hazel (Reference von Borries Lopes, Thiele and Hazel2018) that allows us to study the evaporation of droplets with larger contact angles. In practice, this is achieved by writing the surface energy term in (2.11) with the full metric factor $\gamma \sqrt {1+(\boldsymbol {\nabla } h )^2}$. Here, we also improve the non-conserved mobilities by accordingly replacing the evaporation rate per surface area $\tilde {M}$ in (2.10) by $\tilde {M}\sqrt {1+(\boldsymbol {\nabla } h )^2}$. Below we will refer to this amendment as the full-curvature formulation of the mesoscopic model. For brevity, in the following we refer to the standard long-wave formulation of the mesoscopic model as the ‘long-wave model’. Although experience shows that the described energy-focused strategy can even correct qualitatively incorrect results of a long-wave description (von Borries Lopes et al. Reference von Borries Lopes, Thiele and Hazel2018), and is also applied in descriptions like the Cahn–Hilliard model (Cahn Reference Cahn1965) for the decomposition of a binary mixture (Thiele Reference Thiele2018), it should be noted that the approach may be considered as not being ‘rational’ in an asymptotic sense: part of the neglected terms are of the same order in the smallness parameter of the long-wave approximation as the included improvements of the energy.

Furthermore, we remark that (2.11) is in a somewhat mixed form as it considers at the same time equations of state encoded in the $f$'s for liquid and vapour that may result in phase change, but also uses the film height $h$ as it is the most important quantity for mesoscopic hydrodynamics. As stated above, we treat $\rho _{liq}$ as constant and will consider the vapour as an ideal gas. This is the result of a two-step procedure, explained in the next section (that may be skipped by the reader who is mainly interested in the resulting system of equations).

2.2.2. From real to ideal vapour

An equation of state for a real gas, e.g. a van der Waals gas, predicts coexisting liquid and vapour densities at equilibrium. In an out-of-equilibrium hydrodynamic description, liquid and vapour densities may vary in space and time. However, there is no unique way to express the three fields $\rho _{liq}(\boldsymbol {r},t), \rho (\boldsymbol {r},t)$ and $h(\boldsymbol {r},t)$ in terms of the two fields $\psi _1 (\boldsymbol {r},t)$ and $\psi _2 (\boldsymbol {r},t)$ that our gradient dynamics (2.7) is based on. Therefore, we introduce the following two-step procedure of simplifications that allows us to treat $\rho _{liq}$ as a constant.

(i) We assume a thick homogeneous sharp-interface flat film of thickness $h$, i.e. the energy functional (2.11) becomes

(2.13)\begin{equation} \mathcal{F}_0[\rho_{liq}, \rho_{vap}] = \int_\varOmega \left[ h f_{liq}(\rho_{liq}) +(d-h) f_{vap}(\rho_{vap}) \right] \mathrm{d}\kern0.7pt x\,\mathrm{d} y, \end{equation}

where we also dropped the constant $\gamma$. To calculate $\rho _{liq}$ and $\rho _{vap}$ at coexistence (as functions of temperature $T$), we minimise

(2.14)\begin{align} \mathcal{F}_1[\rho_{liq}, \rho_{vap},h,\xi] &= \int_\varOmega \left[ h f_{liq}(\rho_{liq}) +\xi f_{vap}(\rho_{vap}) \right] \mathrm{d}\,x\,\mathrm{d} y\nonumber\\ &\quad - \tilde \mu \int_\varOmega \left[h \rho_{liq} + \xi \rho -\bar{n}\right]\mathrm{d}\,x\,\mathrm{d} y + \tilde{p} \int_\varOmega \left[h + \xi - d\right]\mathrm{d}\kern0.7pt x\,\mathrm{d} y \end{align}

with respect to $\rho _{liq}, \rho _{vap},h$ and $\xi =d-h$ that we treat as an independent field. Here, $\tilde {\mu }$ and $\tilde {p}$ are Lagrange multipliers for mass and volume conservation. We obtain

(2.15)\begin{equation} \left.\begin{gathered} \tilde{\mu} = f'_{liq}(\rho_{liq}) = f'_{vap}(\rho_{vap}),\\ \tilde{p} = f_{liq}(\rho_{liq}) - \tilde{\mu}\rho_{liq} = f_{vap}(\rho_{vap}) - \tilde{\mu}\rho_{vap}, \end{gathered}\right\} \end{equation}

i.e. the standard Maxwell construction for phase coexistence. These we use to obtain the coexisting $\rho _{liq}$ (and, therefore, $f_{liq}$) and $\rho _{vap}$ analytically or (most likely) numerically. (To do so, we either employ a single function $f=f_{vap}=f_{liq}$ that allows for a liquid–gas phase transition (e.g. the van der Waals free energy) or combine a purely entropic $f_{vap}(\rho ) = k_B T \rho [\log (\varLambda ^3 \rho )-1]$ with an $f_{liq}(\rho )$ that allows for a phase transition.) For the considered isothermal system, this fixes $\rho _{liq}$, $\rho _{vap}$ and $f_{liq}$ at coexistence values.

(ii) Next we approximate the equation of state in the vapour and the liquid phase. For the liquid phase, we neglect any compressibility, i.e. we fix the density in the liquid layer to $\rho _{liq}$. In other words, the ‘liquid branch’ of the equation of state $p(\rho )$ is replaced by a vertical line at $\rho =\rho _{liq}$. The ‘gas branch’ of the equation of state is either directly replaced by the ideal gas law given above or, alternatively, is expanded up to linear order about $\rho =\rho _{vap}$. This results in a shifted and scaled ideal gas law. The latter approach has the advantage that coexistence pressure and concentration are exactly as in the equation of state one started with.

In this way the relation $h=\psi _1/\rho _{liq}$ introduced above becomes meaningful and the relation of the variations with respect to $h$ and to $\psi _1$ alluded to earlier is justified. Also using $\rho =\psi _2/(d-h)$, the free energy functional (2.11) is written as

(2.16)\begin{align} \mathcal{F}[\psi_1,\psi_2] &= \int_\varOmega \left[\frac{\gamma}{2 \rho_{liq}^2} {(\boldsymbol{\nabla} \psi_1)}^2 + g\left(\frac{\psi_1}{\rho_{liq}}\right) +\frac{f_{liq}}{\rho_{liq}} \psi_1 +\left(d-\frac{\psi_1}{\rho_{liq}}\right) f_{vap}\left(\frac{\psi_2}{d-\psi_1/\rho_{liq}}\right) \right.\nonumber\\ &\quad \left. +\left(d-\frac{\psi_1}{\rho_{liq}}\right) f_{air}\left(\rho_{tot}-\frac{\psi_2}{d-\psi_1/\rho_{liq}}\right) \right] \mathrm{d}\,x\,\mathrm{d} y. \end{align}

Here, it only depends on $\psi _1$ and $\psi _2$. Note that $f_{liq}$ is now a constant given, for instance, by the Maxwell construction as explained above. Alternatively, one may directly deduce the value of $f_{liq}$ from the saturation pressure of the liquid by considering the equilibrium state of a liquid film in a saturated atmosphere as done here in § 2.2.4, resulting in (2.29).

2.2.3. Evolution equations

Here, we bring all the information discussed in the previous sections together and present the resulting long-wave model.

The energy functional (2.16) in terms of $\psi _1$ and $\psi _2$ is now minimised together with the particle number constraint (2.12) (Lagrange multiplier $\mu$) with respect to variations in the two fields. This gives

(2.17)\begin{align} \frac{\delta \mathcal{F}}{\delta \psi_1} &= \frac{1}{\rho_{liq}} \left[ -\gamma\Delta h + g'\left(h\right) +f_{liq} -f_{vap}(\rho)+ \rho f_{vap}'(\rho)\right.\nonumber\\ &\quad\left.\vphantom{-\gamma\Delta h + g'\left(h\right) +f_{liq} -f_{vap}(\rho)+ \rho f_{vap}'(\rho)} - f_{air}(\rho_{tot}-\rho) - \rho f_{air}'(\rho_{tot}-\rho)\right] - \mu, \end{align}

where the brackets contain Laplace pressure, Derjaguin pressure, liquid energy and vapour pressure. Note that the somewhat unusual final contribution to the vapour pressure is a direct consequence of the assumption of constant $p_{tot}$. Here and in the following, we use $h$ and $\rho$ as abbreviations where appropriate.

The variation with respect to $\psi _2$ gives

(2.18)\begin{equation} \frac{\delta \mathcal{F}}{\delta \psi_2} =f_{vap}'(\rho) - f_{air}' (\rho_{tot}-\rho) - \mu, \end{equation}

i.e. a difference in chemical potentials. If we now specify vapour and air to be ideal gases, $f_{vap}=k_B T \rho [\log (\varLambda ^3 \rho )-1]$ and $f_{air}=k_B T \rho _{air}[\log (\varLambda ^3 \rho _{air})-1]$ with the mean free path length $\varLambda$, we obtain

(2.19)\begin{equation} \frac{\delta \mathcal{F}}{\delta \psi_1} = \frac{1}{\rho_{liq}} \left\{-\gamma\Delta h+ g'(h) +f_{liq} + k_B T\rho_{tot} - k_B T\rho_{tot} \log[\varLambda^3 (\rho_{tot}-\rho)]\right\} - \mu. \end{equation}

The gradient of this pressure drives the dynamics in the liquid film. The constant liquid energy density $f_{liq}$ and the constant ideal pressure $k_BT\rho _{tot}$ do not contribute to the conserved dynamics, but the final term in the curly parenthesis does. It is a direct consequence of treating case (ii), i.e. of imposing a uniform gas pressure, that is, constant $\rho _{tot}$. The other variation is then

(2.20)\begin{equation} \frac{\delta \mathcal{F}}{\delta \psi_2} =k_B T \{\log(\varLambda^3 \rho) - \log[\varLambda^3 (\rho_{tot}-\rho)]\} - \mu. \end{equation}

Also here the second term in the curly parenthesis is a consequence of the imposed uniform pressure in case (ii).

Introducing the obtained expressions into the general two-field gradient dynamics (2.7) with (2.8) and (2.10) gives

(2.21) \begin{gather} \partial_t \psi_1 = \boldsymbol{\nabla}\, \boldsymbol{\cdot}\, \left(\frac{\psi_1^3}{3 \rho_{liq}^2\eta}\boldsymbol{\nabla} \left\{ -\gamma \Delta h + g'\left(h\right)- k_B T\rho_{tot} \log[\varLambda^3(\rho_{tot}-\rho)] \right\}\right) - \tilde{M}E,\nonumber\\ \partial_t \psi_2 = \boldsymbol{\nabla}\, \boldsymbol{\cdot}\, \left(\tilde{D} k_B T \frac{d-h}{1-\rho/\rho_{tot}}\boldsymbol{\nabla}\rho\right) + \tilde{M} E,\nonumber\\ \mathrm{where}\ E=\frac{1}{\rho_{liq}} \left[-\gamma \Delta h + g'(h)+f_{liq}\right]\nonumber\\ + k_B T \left\{ \frac{\rho_{tot}}{\rho_{liq}} +\left(1- \frac{\rho_{tot}}{\rho_{liq}}\right) \log[\varLambda^3(\rho_{tot}-\rho)] - \log(\varLambda^3 \rho) \right\}. \end{gather}

This is the final result for case (ii).

Then, the saturation vapour density $\rho _{sat}$ above a flat thick film is obtained by setting the transfer term to zero ($E=0$) and dropping capillarity and wettability influences, i.e.

(2.22)\begin{equation} E = \frac{f_{liq}}{\rho_{liq}} + k_B T \left\{ \frac{\rho_{tot}}{\rho_{liq}} +\left(1- \frac{\rho_{tot}}{\rho_{liq}}\right) \log[\varLambda^3(\rho_{tot}-\rho_{sat})] - \log(\varLambda^3 \rho_{sat}) \right\}=0. \end{equation}

The corresponding film height $H$ is determined by the conservation of mass ${\bar {n} = \psi _1+\psi _2 = H \rho _{liq} + (d-H)\rho _{sat}}$.

Note that the somewhat unexpected terms related to $\rho _{tot}$ in the two kinetic equations (2.21) turn out to be well behaved: the additional contributions in the equation for $\psi _2$ provide a factor $1/[1-\rho /\rho _{tot}]$ to the generalized diffusion constant. Normally, $\rho /\rho _{tot}\ll 1$ and might be neglected. If, however, $\rho /\rho _{tot}\to 1$, we approach the limit of pure vapour, where diffusion is not the proper transport process any more: the factor in question diverges, i.e. $\rho$ becomes instantaneously uniform. The additional logarithmic term in the equation for $\psi _1$ corresponds to a flux proportional to $(\boldsymbol {\nabla }\rho )/(1-\rho /\rho _{tot})$. For $\rho /\rho _{tot} \ll 1$, this is the gradient in partial vapour pressure that drives some flow in the adjacent liquid layer, one can see it as an ‘osmotic coupling’. It is normally very small as compared with the other pressure gradients. The local evaporation rate $\tilde {M}E$ contains additional terms proportional to the ratio of total gas density and liquid density $\rho _{tot}/\rho _{liq}$ that we expect to be small. For dry air and water, the ratio is about $10^{-3}$, and humid air has an even lower density than dry air.

For $\rho /\rho _{tot}\ll 1$ and $\rho _{tot}/\rho _{liq}\ll 1$, (2.21) written in terms of $h$ and $\rho$ reduce to

(2.23)\begin{equation} \left.\begin{gathered} \partial_t h ={-}\boldsymbol{\nabla}\, \boldsymbol{\cdot}\, \left\{\frac{h^3}{3 \eta}\boldsymbol{\nabla} \left[ \gamma \Delta h - g'\left(h\right) \right]\right\}- j_{evap},\\ \partial_t[(d-h)\rho] = \boldsymbol{\nabla}\, \boldsymbol{\cdot}\, \left[D (d-h) \boldsymbol{\nabla}\rho\right] + \rho_{liq} j_{evap},\\ \mathrm{where}\ j_{evap}/M ={-}\gamma \Delta h + g'(h)+f_{liq} - \rho_{liq} k_B T \log\left(\frac{\rho}{\rho_{tot}-\rho}\right). \end{gathered}\right\}\end{equation}

Here, we have introduced the diffusion constant $D = \tilde {D} k_B T$ and the evaporation rate constant of $M = \tilde {M} / \rho _{liq}^2$, thus converting the local particle evaporation rate $\tilde {M} E$ into a volume rate $j_{evap} = \tilde {M} E / \rho _{liq}$. Equation (2.23) allows one to study the crossover between transition-limited drop evaporation dynamics (small $\rho$) and diffusion-limited evaporation dynamics ($\rho$ at drop close to saturation $\rho _{sat}$).

Before discussing limiting cases and applying our model to study evaporating drops, we briefly compare it to existing models. Equation (2.21) or the simplified (2.23) represents a seamless reduced model for the evaporation of a sessile droplet in a gap. It captures both limiting cases individually considered by Murisic & Kondic (Reference Murisic and Kondic2011). It is of lower complexity than Sultan et al. (Reference Sultan, Boudaoud and Ben Amar2005) as in the description of liquid and vapour the vertical dimension has been eliminated. Its gradient dynamics form makes it transparent and allows for versatile adaptation to many situations, as further discussed in the conclusion.

In case (ii) considered up to here, convective motion in the gas layer is neglected assuming a uniform total gas pressure $p_{tot}$, i.e. a uniform total gas density ${\rho _{tot}=\rho _{air}+\rho }$ that is an externally controlled parameter of the system. However, this directly couples $\rho _{air}$ to the vapour density $\rho$ and implies the discussed small additional contributions to the pressure in the liquid, to the evaporation/condensation rate and vapour diffusion. In other words, these terms are direct consequences of the gradient dynamics structure. Although, normally, one may neglect them due to their smallness, one needs to keep in mind that ultimately this breaks the thermodynamic consistency and, therefore, may result in unphysical behaviour.

If we consider case (i), evaporation into pure vapour, diffusion is excluded right from the beginning. As in case (ii), we assume that convective motion is very much faster than evaporation, what directly implies that the vapour density $\rho$ is uniform in the entire gas layer. It is set by the conditions at the lateral boundaries of the gap. In other words, the vapour pressure is then a given constant and the governing equation in gradient dynamics form is (2.2). The corresponding energy is (2.11) without the air energy, without the constant $\gamma$ and with constant $\rho =\rho _{vap}$, i.e.

(2.24)\begin{equation} \mathcal{F} = \int_\varOmega \left[ \frac{\gamma}{2} {(\boldsymbol{\nabla} h)}^2 + g(h) +h f_{liq}(\rho_{liq})+(d-h) f_{vap}(\rho_{vap}) \right] \mathrm{d}\,x\,\mathrm{d} y.\end{equation}

Minimisation with respect to film thickness gives

(2.25)\begin{equation} \frac{\delta \mathcal{F}}{\delta h} ={-}\gamma\Delta h+ g'(h) + f_{liq} - f_{vap}, \end{equation}

i.e. (2.2) becomes

(2.26)\begin{equation} \partial_t h={-}\boldsymbol{\nabla}\, \boldsymbol{\cdot}\, \left\{Q\boldsymbol{\nabla} [\gamma\Delta h - g'(h)] \right\} + M\left[\gamma\Delta h - g'(h) - f_{liq} + f_{vap} + p_{vap}\right]. \end{equation}

Grouping all constants ( $f_{liq}, f_{vap}, p_{vap}$) in the evaporation term into a single one, the equation corresponds to the long-wave description in the transition-limited case (Pismen Reference Pismen2004; Thiele Reference Thiele2010) where evaporation/condensation is controlled by the difference between the sum of Laplace and Derjaguin pressure in the liquid and the constant and uniform imposed vapour pressure.

2.2.4. Specific functions and parameters

One may now proceed by non-dimensionalising, thereby using sensible values of $\rho _{liq}$ and $f_{liq}$ as parameters. However, as they are not independent, this may result in artefacts. The better approach is to employ a specific equation of state/free energy that allows for liquid–gas phase transition and calculate these values for a particular temperature $T$.

An example is the van der Waals equation of state, see, e.g. Landau & Lifshitz (Reference Landau and Lifshitz1980, §§ 76 and 84) giving the pressure

(2.27)\begin{equation} p_{vdW}=\frac{k_B T\rho}{1-b\rho}-a\rho^2, \end{equation}

where $a$ and $b$ are constants (attraction strength and effective excluded volume, respectively). Approximated, in the gas phase with $a,b\to 0$ it becomes the ideal gas law $p_{id}=k_B T\rho$. The corresponding Helmholtz free energy per volume is

(2.28)\begin{equation} f_{vdW}={-}k_B T \rho\left[\log\left(\frac{1-\rho b}{\varLambda^3\rho}\right) + 1 \right]-a\rho^2, \end{equation}

becoming $f_{id}=k_B T \rho [\log (\varLambda ^3\rho ) - 1 ]$ for $a,b\to 0$. Note that $p=-f+\rho f'= \rho ^2\partial_\rho (f(\rho )/\rho )$.

Furthermore, the saturation vapour density $\rho _{sat}$ follows from the equilibrium condition, e.g. in the simple case of a thick flat liquid film (no Laplace and Derjaguin pressure) setting the evaporation rate to zero in (2.23),

(2.29)\begin{equation} j_{evap}/M = f_{liq} - \rho_{liq} k_B T \log\left(\frac{\rho_{sat}}{\rho_{tot}-\rho_{sat}}\right) = 0.\end{equation}

The saturation vapour pressure $p_{sat}$ resulting from the equation of state may then be used to express the humidity as $\phi = p / p_{sat}$.

If one assumes a specific wetting potential $g(h)$ that allows for partial wetting, a second spatially homogeneous equilibrium state is found typically at very small height $h_a$ that corresponds to a microscopic adsorption layer

(2.30)\begin{equation} g'(h_a) = \rho_{liq}k_B T \log\left(\frac{\rho}{\rho_{tot}-\rho}\right) - f_{liq}.\end{equation}

The height $h_a$ depends on the vapour density $\rho$ and will be slightly shifted from the height $h_p$ where the Derjaguin pressure vanishes ($g'(h_p) = 0$), i.e. the adsorption layer height of the saturated case $h_p = h_a(\rho _{sat})$. In fact, this allows for a set of spatially homogeneous equilibrium states, where the vapour density is an arbitrary constant $0 < \rho \leqslant \rho _{sat}$ controlled externally, e.g. by the humidity of the laboratory, while the substrate is macroscopically dry (the ‘moist case’ of de Gennes Reference de Gennes1985) with $h = h_a$. Note that due to the saturation conditions (2.29) and (2.30), we find that $g'(h_a) \leqslant 0$, i.e. the adsorption layer height is shifted to lower values $h_a \leqslant h_p$. When, on the other hand, an initial $\rho (\boldsymbol {r}, t) \geqslant \rho _{sat}$ is given, no equilibrium is guaranteed to exist (depending on the choice of the wetting potential). The evaporation rate $E$ will then become negative leading to a transfer of material from the vapour phase to the liquid film, i.e. condensation.

In the following, we use the simple wetting potential for partially wetting liquids

(2.31)\begin{equation} g(h) = \frac{H_A}{2h^2} \left(\frac{2 h_p^3}{5 h^3} - 1\right),\end{equation}

with the Hamaker constant $H_A \approx 5/3\gamma h_p^2 \theta _e^2$ that relates to the macroscopic equilibrium contact angle $\theta _e$ by a mesoscopic Young's law: $\cos \theta _e = 1 + g(h_p)/\gamma$ (Churaev Reference Churaev1995).

For the simulation of the long-wave model in § 3, we use the following set of realistic parameters (unless stated otherwise): initial drop height $h_0 = h(r=0, t=0)={1}\ \textrm {mm}$, gap height $d={3}\ \textrm {mm}$, precursor film thickness $h_p = h_0/1000$, diffusion constant $D = \tilde {D} k_B T = 2.8\times 10^{-5}\ \textrm {m}^2\ \textrm {s}^{-1}$ (Lee & Wilke Reference Lee and Wilke1954; Marrero & Mason Reference Marrero and Mason1972), temperature $T={22}\,^{\circ }\textrm {C}$, contact angle $\theta _e = {44.38}^{\circ }$ and a domain radius of $L={25}\ \textrm {mm}$. The liquid properties are chosen according to water at room temperature: viscosity $\eta = 8.9\times 10^{-4}\ \textrm {Pa}\ \textrm {s}$, surface tension $\gamma = 72\times 10^{-3}\ \textrm {N}\ \textrm {m}^{-1}$, particle density according to mass density and molar mass $\rho _{liq} = ({997}\ \textrm {kg}\ \textrm {m}^{-3})/({18}\ \textrm {kg}\ \textrm {mol}^{-1})\,N_A$, and the saturation vapour pressure $p_{sat}={2643}\ \textrm {Pa}$. The parameter values are extracted from Lide (Reference Lide2004). The constant of phase change $M$ is employed as a free parameter that allows us to move between the different limiting cases.

2.2.5. Radial symmetry and numerical implementation

For an efficient numerical simulation of the dynamics, we consider a radial symmetry, i.e. for the long-wave model, $h(\boldsymbol {x}, t) = h(r, t)$ and $\psi (\boldsymbol {x}, t) = \psi (r, t)$. This allows us to reduce (2.23) to a spatially one-dimensional model in polar coordinates,

(2.32) \begin{gather} \partial_t h = \left(\partial_r + \frac{1}{r}\right) \left\{\frac{h^3}{3 \eta}\partial_r \left[-\gamma \left( \partial_r^2 h + \frac{1}{r} \partial_r h \right) + g'\left(h\right) \right]\right\} - j_{evap},\nonumber\\ \partial_t[(d-h)\rho] = \left(\partial_r + \frac{1}{r}\right) \left[D (d-h) \partial_r \rho\right] + \rho_{liq} j_{evap},\nonumber \\ \mathrm{where}\ j_{evap}/M ={-}\gamma \left( \partial_r^2 h + \frac{1}{r} \partial_r h \right) + g'(h)+f_{liq}- \rho_{liq}k_B T \log\left(\frac{\rho}{\rho_{tot}-\rho}\right).\end{gather}

The dynamic equations are then solved using the finite-element method implemented in the C$++$ library oomph-lib (Heil & Hazel Reference Heil and Hazel2006), which employs a second-order backward differentiation formula for temporal integration. In particular, the library provides both spatial and temporal adaptivity, which is essential due to the multi-scale character of the dynamics, e.g. when spatially resolving the fields in the bulk and near the contact line. In the simulations, the domain is typically discretised with ${\approx }1000$ grid points with spacings that range from $2^{-7}$ to $2^{-18}$ times the domain size.

In these calculations, we employ an adsorption layer of height $h_p$ that is by a factor of $10^3$ smaller than the initial drop height $h(r=0, t=0)$. Larger ratios are possible but result in a strongly increasing numerical effort. As the effective adsorption layer height is coupled to the gas phase (see, e.g. (2.30)), the overall pressure balance can cause liquid transport through the layer during equilibration. The resulting fluxes are negligible if the adsorption layer is very thin compared with the droplet size. The effect can be further suppressed by modulating the evaporation rate coefficient $M(h)$ with a (smooth) step function, effectively disabling evaporation from the adsorption layer. This is particularly important on the slow time scale of droplet evaporation.

To ensure smoothness of all fields at the drop centre ($r=0$), we employ natural (homogeneous) Neumann boundary conditions

(2.33a–c)\begin{equation} \left.\partial_r h\right|_{r=0} = 0, \quad\left.\partial_r \frac{\delta \mathcal{F}}{\delta \psi_1}\right|_{r=0} = 0 \quad\text{and}\quad \left.\partial_r [(d-h)\rho]\right|_{r=0} = 0.\end{equation}

These conditions also ensure zero liquid and vapour flux through the computational domain boundary at $r=0$. At the outer boundary far from the drop the gap between the plates is open, i.e. air and vapour can freely be exchanged with the surrounding laboratory. Thus, we assume a constant lab humidity (or vapour concentration) $\rho _{lab}$ corresponding to a Dirichlet boundary conditions at $r=L$. Together with natural Neumann boundary conditions for the film profile $h$ and the chemical potential $\delta \mathcal {F} / \delta \psi _1$, this gives the remaining boundary conditions

(2.34a–c)\begin{equation} \left.\partial_r h\right|_{r=L} = 0, \quad\partial_r \left.\frac{\delta \mathcal{F}}{\delta \psi_1}\right|_{r=L} = 0 \quad\text{and}\quad \rho|_{r=L} = \rho_{lab}.\end{equation}

In particular, these conditions allow for a non-zero vapour flux through the outer boundary into the laboratory environment. The external vapour concentration is here chosen such that it corresponds to a low relative humidity of, i.e., $\rho _{lab}/\rho _{sat} = 0.1$.

Results obtained with the developed long-wave model are presented in § 3. In the subsequent § 4 they are compared with Stokes-equation results and corresponding experiments.