2.1. Formulation

Consider a mixture of two fluids with partial densities $\rho _{i}$, where $i=1,2$. The first fluid will be referred to as either ‘liquid’ or ‘vapour’ (depending on the phase it is in), and the second fluid will be referred to as ‘air’.

Under the assumption of spherical symmetry, $\rho _{i}$ generally vary with the radial coordinate $r$ and time $t$. In the classical model, however, one assumes that inside the drop

(2.1)\begin{equation} \rho_{1}=\rho_{1}^{(l)},\quad\rho_{2}=0\quad\text{if}\ r\in(0,r_{d}) , \end{equation}

where $\rho _{1}^{(l)}$ is a known constant, and the drop's radius $r_{d}$ depends on $t$.

Outside the drop, both densities do vary. The air density is much larger than that of the vapour, but their variations are of the same order (as, physically, the pressure should be spatially uniform). As a result, the variations of $\rho _{2}$ are relatively weak, so one assumes

(2.2)\begin{equation} \rho_{2}=\rho_{2}^{(a)}\quad\text{if}\ r\in( r_{d},\infty), \end{equation}

where the dry-air density $\rho _{2}^{(a)}$ is a known parameter. Note that, physically, $\rho _{1}^{(l)}$ and $\rho _{2}^{(a)}$ are inter-related by the requirement that the pressure difference between the inside and outside of the drop equals the capillary correction due to the curvature of the interface.

The classical model of evaporation also involves the liquid's thermal conductivity $\kappa ^{(l)}$, air's thermal conductivity $\kappa ^{(a)}$, saturated vapour density $\rho _{1}^{(v.sat)}$, vaporisation heat $\varDelta$, and vapour-in-the-air diffusivity $\mathcal {D}$. The dependence of these parameters on the temperature $T$ is supposed to be known.

The classical model is usually combined with the quasi-steady approximation, based on an assumption that the diffusivity and thermal conductivity are sufficiently large (which they are indeed for water under normal conditions). Then the term $\partial \rho _{1}/\partial t$ in the diffusion equation and the term $\partial T/\partial t$ in the heat-conduction equation can both be neglected. The simplest version of the classical model also neglects the Stefan flow (which carries the vapour away from the surface of the liquid), the Soret effects (the mass flux due to the temperature gradient), and the Dufour effect (the heat flux due to the density gradient).

Then, for the spherically symmetric case, one obtains

(2.3)\begin{gather} \frac{\partial}{\partial r}\left( r^{2}\,\frac{\partial T}{\partial r}\right) =0\quad\text{if}\ r\in( 0,r_{d}) , \end{gather}

(2.4)\begin{gather}\frac{\partial}{\partial r}\left( r^{2}\,\frac{\partial\rho_{1}}{\partial r}\right) =0,\quad \frac{\partial}{\partial r}\left( r^{2}\,\frac{\partial T}{\partial r}\right) =0\quad\text{if}\ r\in( r_{d} ,\infty) . \end{gather}

At the drop's centre, the usual boundary condition applies:

(2.5)\begin{equation} \frac{\partial T}{\partial r}=0\quad\text{at}\ r=0.\end{equation}

To formulate the boundary conditions at the drop's surface, introduce

(2.6a,b)\begin{equation} T^{({\pm})}=\lim_{r\rightarrow r_{d}\pm0}T,\quad\left( \frac{\partial T}{\partial r}\right)^{({\pm})}=\lim_{r\rightarrow r_{d}\pm0} \frac{\partial T}{\partial r}, \end{equation}

where $^{(+)}$ and $^{(-)}$ mean ‘just outside’ and ‘just inside’ the drop, respectively. Then the continuity of the temperature and heat flux imply

(2.7)\begin{gather} T^{(-)}=T^{(+)}\quad\text{at}\ r=r_{d}, \end{gather}

(2.8)\begin{gather}\kappa^{(l)}\left( \frac{\partial T}{\partial r}\right)^{(-)}-\kappa ^{(a)}\left( \frac{\partial T}{\partial r}\right)^{(+)}=\dot{r}_{d}\rho _{1}^{(l)}\varDelta\quad\text{at}\ r=r_{d}, \end{gather}

where

(2.9)\begin{equation} \dot{r}_{d}=\frac{\mathrm{d}r_{d}}{\mathrm{d}t}. \end{equation}

As discussed in the Introduction, the classical model assumes the Maxwell boundary condition for the vapour density at the drop's surface, i.e.

(2.10)\begin{equation} \rho_{1}=\rho_{1}^{(v.sat)}\quad\text{at}\ r=r_{d},\end{equation}

and the following conditions should be imposed on the vapour–air mixture far away from the drop:

(2.11)\begin{equation} \rho_{1}\rightarrow\rho_{1}^{(a)},\quad T\rightarrow T^{(a)}\quad \text{as}\ r\rightarrow\infty.\end{equation}

To ensure that the liquid is evaporating – not condensating or being in equilibrium – let the vapour at infinity be undersaturated:

(2.12)\begin{equation} \rho_{1}^{(a)}<\rho_{1}^{(a.sat)(a)}, \end{equation}

where $\rho _{1}^{(a.sat)(a)}=\rho _{1}^{(a.sat)}(T^{(a)})$. Finally, the velocity of the drop's (receding) surface is related to the evaporative flux by

(2.13)\begin{equation} \rho_{1}^{(l)}\dot{r}_{d}=\mathcal{D}\,\frac{\partial\rho}{\partial r} \quad\text{at}\ r=r_{d}.\end{equation}

Boundary-value problem (2.3)–(2.13) governs $\rho _{1}(r,t)$, $T(r,t)$ and $r_{d}(t)$.

2.2. The dependence of $\rho _{1}^{(v.sat)}$ on $T$

One's experience with raindrops sitting on ones's skin suggests that evaporative cooling does not exceed several degrees. This seems to enable one to neglect the dependence of the external parameters on the temperature – say, ‘freeze’ them at $T=T^{(a)}$. Such an approximation is indeed applicable to all parameters but one – namely, the saturated vapour density $\rho _{1}^{(v.sat)}$.

This claim is illustrated for water in figure 1. One can see that the relative change of $\rho _{1}^{(v.sat)}$ across the room temperature range exceeds those of the remaining parameters.

Figure 1.

The temperature dependencies of the saturated vapour density $\rho _{1}^{(v.sat)}$, vapour-in-the-air diffusivity $\mathcal {D}$, thermal conductivity $\kappa ^{(a)}$ of the air, and vaporisation heat $\varDelta$ (all non-dimensionalised on their reference values): (1) $\rho _{1}^{(v.sat)}(T)/\rho _{1}^{(v.sat)}(25\,^{\circ }\mathrm {C})$, (2) $\mathcal {D}(T)/\mathcal {D}(25\,^{\circ }\mathrm {C})$, (3) $\kappa ^{(a)}(T)/\kappa ^{(a)}(25\,^{\circ }\mathrm {C})$, and (4) $\varDelta (T)/\varDelta (25\,^{\circ }\mathrm {C})$. Curves (1)–(4) correspond to the empirical formulae of Wagner & Pruss (Reference Wagner and Pruss1993), Engineering ToolBox (2018), White (Reference White2005) and Lindstrom & Mallard (Reference Lindstrom and Mallard1997), respectively.

In this paper, the dependence $\rho _{1}^{(v.sat)}(T)$ is approximated by Antoine's equation (Antoine Reference Antoine1888). If written in terms of the density – as opposed to the pressure, as in the original formulation – it takes the form

(2.14)\begin{equation} \frac{\rho_{1}^{(v.sat)}}{\rho_{A}}=\frac{T_{A}}{T}\exp\left( -\frac{T_{A} }{T}\right) ,\end{equation}

where, for water,

(2.15a,b)\begin{equation} \rho_{A}=6.6326\times10^{4}\ \mathrm{kg\ m}^{{-}3},\quad T_{A}=5292\ \mathrm{K}.\end{equation}

These values have been obtained by fitting (2.14) to the high-accuracy IAPWS empirical formula (Wagner & Pruss Reference Wagner and Pruss1993) for the room temperature range. The relative deviation of the two formulae is less than 0.23 %, which is much better than the accuracy of Antoine's equation with $\rho _{A}$ and $T_{A}$ taken from a thermodynamics handbook.

Not only does Antoine's equation provide an accurate approximation, but it also has physical meaning. As shown by Benilov (Reference Benilov2020), (2.14) describes the generic van der Waals fluid in the low-temperature limit – i.e. when $T$ is much closer to the freezing point than to the critical point. Observe also that $T_{A}$ is large, which is what makes the dependence of $\rho _{1}^{(v.sat)}$ on $T$ strong.

Equation (2.14) can be rearranged conveniently using a reference temperature (as done by Sobac et al. Reference Sobac, Talbot, Haut, Rednikov and Colinet2015). Choosing that to be $T^{(a)}$, one obtains

(2.16)\begin{equation} \frac{\rho_{1}^{(v.sat)}}{\rho_{1}^{(v.sat)(a)}}=\frac{T^{(a)}}{T}\exp\left( -\frac{T_{A}}{T}+\frac{T_{A}}{T^{(a)}}\right) ,\end{equation}

where $\rho _{1}^{(v.sat)(a)}$ is the saturated vapour density at infinity. One can further assume $( T^{(a)}-T) /T^{(a)}\ll 1$ and reduce (2.16) to

(2.17)\begin{equation} \rho_{1}^{(v.sat)}=\rho_{1}^{(v.sat)(a)}\exp\left[ \frac{T_{A}}{T^{(a)2} }\,( T-T^{(a)}) \right] .\end{equation}

Observe that even if $T$ is close to $T^{(a)}$, the expression in the square brackets can still be order-one (because $T_{A}\gg T^{(a)}$).

2.3. Solution of boundary-value problem (2.3)–(2.13)

Equations (2.3)–(2.4) can be solved readily as

(2.18)\begin{gather} T=\textrm{const}_{1}+\frac{\textrm{const}_{2}}{r} \quad\text{if}\ r\in( 0,r_{d}) , \end{gather}

(2.19)\begin{gather}\rho_{1}=\textrm{const}_{3}+ \frac{\textrm{const}_{4}}{r},\quad T=\textrm{const}_{5}+ \frac{\textrm{const}_{6}}{r} \quad\text{if}\ r\in( r_{d},\infty), \end{gather}

where the constants of integration can be fixed via boundary conditions (2.5)–(2.11):

(2.20) \begin{gather} \left.\begin{gathered} \textrm{const}_{1}=T^{(a)}+r_{d}\dot{r}_{d}\,\frac{\rho_{1}^{(l)}\varDelta }{\kappa^{(a)}},\quad\textrm{const}_{2}=0,\quad\textrm{const} _{3}=\rho_{1}^{(a)},\\ \textrm{const}_{4}=r_{d}\left[ \rho_{1}^{(v.sat)}\left( T^{(a)} +r_{d}\dot{r}_{d}\,\frac{\rho_{1}^{(l)}\varDelta}{\kappa^{(a)}}\right) -\rho _{1}^{(a)}\right], \quad\textrm{const}_{5}=T^{(a)},\\ \textrm{const}_{6}=r_{d}^{2}\dot{r}_{d}\,\frac{\rho_{1}^{(l)}\varDelta }{\kappa^{(a)}}. \end{gathered}\right\} \end{gather}

Recalling (2.16) for $\rho _{1}^{(v.sat)}$ and (2.13) for $\dot {r}_{d}$, one obtains

(2.21)\begin{equation} \underset{\text{diffusion}}{\underbrace{\frac{\rho_{1}^{(l)}}{\mathcal{D} \rho_{1}^{(v.sat)(a)}}\,\dot{r}_{d}r_{d}}}+\exp \underset{\text{non-isothermality}}{\underbrace{\dfrac{T_{A}\rho_{1}^{(l)}\varDelta}{T_{{}}^{(a)2}\kappa^{(a)}}\,\dot{r}_{d}r_{d}}}=H,\end{equation}

where

(2.22)\begin{equation} H=\frac{\rho_{1}^{(a)}}{\rho_{1}^{(v.sat)(a)}} \end{equation}

is the relative humidity of air far away from the drop.

Ordinary differential equation (2.21) determines $r_{d}(t)$; the terms on its left-hand side are labelled according to their physical meanings.

Equation (2.21) yields the following expression for a drop's time of evaporation:

(2.23)\begin{equation} t_{e}=r_{d}^{2}(0)\,\frac{x\rho_{1}^{(l)}}{2\mathcal{D}\rho_{1}^{(v.sat)(a)} },\end{equation}

where $r_{d}(0)$ is the drop's initial radius, and $x$ is the solution of the non-dimensional algebraic equation

(2.24)\begin{equation} {-x}+\exp( -\alpha x) =H,\end{equation}

with

(2.25)\begin{equation} \alpha=\frac{T_{A}\mathcal{D}\rho_{1}^{(v.sat)(a)}\varDelta}{T^{(a)2}\kappa^{(a)}}.\end{equation}

Formula (2.23) is often referred to as the ‘$d^{2}$ law’, where $d$ is the drop's diameter.

It is instructive to estimate $\alpha$ for, say, the midpoint of the room temperature range, $T^{(a)}=25\,^{\circ }\mathrm {C}$, and pressure $1$ atm. Letting the liquid be water, one can use the NIST online calculator (Lindstrom & Mallard Reference Lindstrom and Mallard1997) to obtain

(2.26)\begin{equation} \rho_{1}^{(l)}=997.05\ \mathrm{kg\ m}^{{-}3},\quad\rho_{1}^{(v.sat)(a)} =0.023075\ \mathrm{kg\ m}^{{-}3},\quad \varDelta=2441.7\ \mathrm{kJ\ kg}^{{-}1}, \end{equation}

whereas for air, Engineering ToolBox (2003, 2009, 2018) yield

(2.27) \begin{equation} \rho_{2}^{(a)}=1.184\ \mathrm{kg\,m}^{{-}3},\quad\kappa^{(a)} =0.02624\ \mathrm{W}\ \mathrm{m}^{{-}1}\ \mathrm{K}^{{-}1},\quad\mathcal{D} =2.49\times10^{{-}5}\ \mathrm{m}^{2}\ \mathrm{s}^{{-}1}. \end{equation}

Substituting values (2.26)–(2.27) into (2.25), and recalling the value in (2.15) for $T_{A}$, one obtains

(2.28)\begin{equation} \alpha\approx3.1828.\end{equation}

3.1. The governing equations

Let $z$ be the vertical Cartesian coordinate, and let $w$ be the vertical velocity of the mixture as a whole. Introduce also the pressure $p$, and keep in mind that it is not an independent unknown: since the DIM assumes the fluid to be compressible, $p$ is a known function of $\rho _{1}$, $\rho _{2}$ and $T$ (the equation of state). The shear and bulk viscosities – $\mu _{s}$ and $\mu _{b}$, respectively – are also known functions of the same arguments.

Introduce the partial chemical potentials $G_{1}$ and $G_{2}$. Their thermodynamic definition is given in Appendix A.1 – but technically, one can simply treat them as given functions of $\rho _{1}$, $\rho _{2}$ and $T$, such that

(3.1)\begin{equation} \frac{\partial p}{\partial\rho_{j}}=\sum_{i}\rho_{i}\,\frac{\partial G_{i} }{\partial\rho_{j}}.\end{equation}

In the low-density limit, $p$ and $G_{i}$ tend to their ideal gas approximations:

(3.2)\begin{equation} p\sim\sum_{i}R_{i}\rho_{i},\quad G_{i}\sim R_{i}T\ln\rho_{i}\quad \text{as}\ \rho_{i}\rightarrow0,\end{equation}

where $R_{i}$ are the specific gas constants. For water and air, they are

(3.3a,b)\begin{equation} R_{1}=461.53\ \mathrm{m}^{2}\ \mathrm{s}^{{-}2}\ \mathrm{K}^{{-}1},\quad R_{2}=289.41\ \mathrm{m}^{2}\ \mathrm{s}^{{-}2}\ \mathrm{K}^{{-}1}, \end{equation}

with $R_{2}$ calculated as the $70/30$ weighted average of the gas constants of nitrogen and oxygen.

The DIM comprises the standard hydrodynamic equations involving an undetermined external force $F_{i}$, which will be identified later with the van der Waals force. Under the slow-flow (Stokes) approximation, a compressible isothermal binary fluid is governed by

(3.4a,b)\begin{gather} \frac{\partial( \rho_{1}+\rho_{2}) }{\partial t}+\frac {\partial[ ( \rho_{1}+\rho_{2}) w] }{\partial z}=0,\quad \frac{\partial\rho_{1}}{\partial t}+\frac{\partial( \rho_{1}w+J) }{\partial z}=0, \end{gather}

(3.5)\begin{gather}\frac{\partial p}{\partial z}=\frac{\partial}{\partial z}\left( \eta\,\frac{\partial w}{\partial z}\right) +\rho_{1}F_{1}+\rho_{2}F_{2}, \end{gather}

where the effective viscosity is

(3.6)\begin{equation} \eta=\frac{4\mu_{s}}{3}+\mu_{b}, \end{equation}

the diffusive flux is (e.g. Giovangigli Reference Giovangigli2021)

(3.7)\begin{equation} J=D\left( F_{1}-\frac{\partial G_{1}}{\partial z}-F_{2}+\frac{\partial G_{2} }{\partial z}\right) , \end{equation}

and $D$ is the diffusion coefficient.

To understand the relationship of $D$ with the diffusivity $\mathcal {D}$ from the classical model, assume that the fluid is diluted, so that the van der Waals force is negligible ($F_{i}=0$) and the chemical potentials are given by (3.2). As a result, (3.7) reduces to

(3.8)\begin{equation} J\sim D\left( -\frac{R_{1}T}{\rho_{1}}\,\frac{\partial\rho_{1}}{\partial z}+\frac{R_{2}T}{\rho_{2}}\,\frac{\partial\rho_{2}}{\partial z}\right) \quad\text{as}\ \rho_{i}\rightarrow0. \end{equation}

Let the vapour density be much smaller than that of the air ($\rho _{1}\ll \rho _{2}$), but let their variations be comparable ($\partial \rho _{1}/\partial z\sim \partial \rho _{2}/\partial z$) – so that the above expression reduces to

(3.9a,b)\begin{equation} J\sim{-}D\,\frac{R_{1}T}{\rho_{1}}\,\frac{\partial\rho_{1}}{\partial z} \quad\text{as } \rho_{i}\rightarrow0, \frac{\rho_{1}}{\rho_{2} }\rightarrow0. \end{equation}

A comparison of this expression with that for the classical diffusive flux yields

(3.10a,b)\begin{equation} D\sim\frac{\rho_{1}}{R_{1}T}\,\mathcal{D}\quad\text{as}\ \rho _{i}\rightarrow0, \frac{\rho_{1}}{\rho_{2}}\rightarrow0.\end{equation}

For a dense mixture of fluids with comparable densities, $D$ cannot be deduced from $\mathcal {D}$, but its dependence on $( \rho _{1},\rho _{2},T)$ is supposed to be known from measurements.

For a multicomponent flow, the (one-dimensional) DIM approximation of the van der Waals force is (Benilov Reference Benilov2023a)

(3.11)\begin{equation} F_{i}=\sum_{j}K_{ij}\,\frac{\partial^{3}\rho_{j}}{\partial z^{3}}, \end{equation}

where the Korteweg matrix $K_{ij}$ characterises the interaction of molecules of components $i$ and $j$ (Newton's third law implies $K_{ij}=K_{ji}$). For the water–air combination, the Korteweg matrix is (see Appendix B)

(3.12)\begin{equation} \begin{bmatrix} K_{11} & K_{12}\\ K_{21} & K_{22} \end{bmatrix}=\begin{bmatrix} 1.87810 & 0.84978\\ 0.84978 & 1.36880 \end{bmatrix} \times10^{{-}17}\ \mathrm{m}^{7}\ \mathrm{s}^{{-}2}\ \mathrm{kg}^{{-}1}. \end{equation}

Next, the pressure term in the momentum equation (3.5) can be rewritten conveniently in terms of the chemical potentials. Using identity (3.1), one obtains

(3.13)\begin{equation} \rho_{1}\,\frac{\partial G_{1}}{\partial z}+\rho_{2}\,\frac{\partial G_{2} }{\partial z}=\frac{\partial}{\partial z}\left( \eta\,\frac{\partial w}{\partial z}\right) +\rho_{1}F_{1}+\rho_{2}F_{2}. \end{equation}

This equation and (3.7) can be viewed as a set of algebraic equations for $F_{1}-\partial G_{1}/\partial z$ and $F_{2}-\partial G_{2}/\partial z$. Solving these equations and recalling (3.11) for $F_{i}$, one obtains

(3.14)\begin{gather} \frac{\partial}{\partial z}\left( \sum_{j}K_{1j}\,\frac{\partial^{2}\rho_{j} }{\partial z^{2}}-G_{1}\right) ={-}\frac{1}{\rho_{1}+\rho_{2}}\left[ \frac{\partial}{\partial z}\left( \eta\,\frac{\partial w}{\partial z}\right) -\frac{\rho_{2}J}{D}\right] , \end{gather}

(3.15)\begin{gather}\frac{\partial}{\partial z}\left( \sum_{j}K_{2j}\,\frac{\partial^{2}\rho_{j} }{\partial z^{2}}-G_{2}\right) ={-}\frac{1}{\rho_{1}+\rho_{2}}\left[ \frac{\partial}{\partial z}\left( \eta\,\frac{\partial w}{\partial z}\right) +\frac{\rho_{1}J}{D}\right] . \end{gather}

Equations (3.4) and (3.14)–(3.15) determine the unknowns $\rho _{1}$, $\rho _{2}$, $w$ and $J$.

3.2. Non-dimensionalisation

In what follows, the vertical coordinate $z$ will be non-dimensionalised on the interfacial thickness $\bar {z}$, which implies that distances should be measured from the current height $Z(t)$ of the interface. Mathematically, this corresponds to the change of variables $( z,t) \rightarrow ( z_{new},t_{new})$, where

(3.16)\begin{equation} z_{new}=z-Z(t),\quad t_{new}=t. \end{equation}

Rewriting (3.4) in terms of the new variables, introducing $\dot {Z}=\mathrm {d}Z/\mathrm {d}t$, and omitting the subscript $_{new}$, one obtains

(3.17a)\begin{gather} \frac{\partial( \rho_{1}+\rho_{2}) }{\partial t}-\dot{Z}\,\frac{\partial\rho_{2}}{\partial z}+\frac{\partial[ ( \rho_{1}+\rho_{2}) w] }{\partial z}=0, \end{gather}

(3.17b)\begin{gather}\frac{\partial\rho_{1} }{\partial t}-\dot{Z}\,\frac{\partial\rho_{1}}{\partial z}+\frac{\partial( \rho_{1}w+J) }{\partial z}=0, \end{gather}

whereas (3.14)–(3.15) remain the same as before.

There are three density scales in the problem – the liquid density, the air density and the saturated vapour density. Since they differ from each other by orders of magnitude, no single scaling of $\rho _{i}$ would fit all the asymptotic regions. Still, a ‘generic’ non-dimensionalisation is needed, so let the density scale be that of air, $\bar {\rho }=\rho _{2}^{(a)}$.

Introducing the air viscosity scale $\bar {\eta }$, define the characteristic pressure and velocity by

(3.18a,b)\begin{equation} \bar{p}=\frac{\bar{K}\bar{\rho}^{2}}{\bar{z}^{2}},\quad \bar{w}=\frac{ \bar{p}\bar{z}}{\bar{\eta}},\end{equation}

respectively. These choices for $\bar {p}$ and $\bar {w}$ correspond to an asymptotic regime where the pressure gradient, viscous stress and van der Waals force in the momentum equation (3.5) are of the same order.

Let $\bar {t}$ be the evaporation time scale, and introduce

(3.19)\begin{equation} \varepsilon=\sqrt{\frac{\bar{z}}{\bar{t}\bar{w}}}. \end{equation}

This parameter is responsible for the quasi-steady approximation – hence it is expected to be small. Another important parameter is the advection/diffusion ratio, defined by

(3.20)\begin{equation} \delta=\frac{\bar{w}\bar{\rho}^{2}\bar{z}}{\bar{D}\bar{p}}, \end{equation}

where $\bar {D}$ is the characteristic value of the diffusion coefficient $D$.

The following non-dimensional variables will be used:

(3.21a–j) $$\begin{gather} \left.\begin{gathered} t_{nd}=\frac{t}{\bar{t}},\quad z_{nd}=\frac{z}{\bar{z}},\\ ( \rho_{i})_{nd}=\frac{\rho_{i}}{\bar{\rho}},\quad w_{nd} =\frac{w}{\varepsilon\bar{w}},\quad J_{nd}=\frac{J}{\varepsilon\bar{\rho} \bar{w}},\quad Z_{nd}=\varepsilon Z,\\ \eta_{nd}=\frac{\eta}{\bar{\eta}},\quad D_{nd}=\frac{D}{\bar{D}} ,\quad ( K_{ij})_{nd}=\frac{K_{ij}}{\bar{K}},\quad (G_{i})_{nd}=\frac{\bar{\rho}}{\bar{p}}G_{i}. \end{gathered}\right\} \end{gather}$$

Rewriting (3.14)–(3.17) in terms of the non-dimensional variables, and omitting the subscript $_{nd}$, one obtains

(3.22)\begin{gather} \frac{\partial}{\partial z}\left( \sum_{j}K_{1j}\,\frac{\partial^{2}\rho_{j} }{\partial z^{2}}-G_{1}\right) ={-}\frac{\varepsilon}{\rho_{1}+\rho_{2}}\left[ \frac{\partial}{\partial z}\left( \eta\,\frac{\partial w}{\partial z}\right) -\delta\,\frac{\rho_{2}J}{D}\right] , \end{gather}

(3.23)\begin{gather}\frac{\partial}{\partial z}\left( \sum_{j}K_{2j}\,\frac{\partial^{2}\rho_{j} }{\partial z^{2}}-G_{2}\right) ={-}\frac{\varepsilon}{\rho_{1}+\rho_{2}}\left[ \frac{\partial}{\partial z}\left( \eta\,\frac{\partial w}{\partial z}\right) +\delta\,\frac{\rho_{1}J}{D}\right] , \end{gather}

(3.24)\begin{gather}\varepsilon\,\frac{\partial( \rho_{1}+\rho_{2}) }{\partial t} +\frac{\partial[ ( \rho_{1}+\rho_{2}) ( w-\dot {Z}) ] }{\partial z}=0,\quad\varepsilon\,\frac{\partial\rho_{1} }{\partial t}+\frac{\partial[ \rho_{1}( w-\dot{Z}) +J] }{\partial z}=0. \end{gather}

Evidently, the non-dimensional equations (3.22)–(3.24) can be transformed back to their dimensional versions, (3.14)–(3.17), by setting $\varepsilon =\delta =1$. This shortcut will be used for re-dimensionalising the results obtained.

In what follows, one will need the non-dimensional versions of the low-density expressions (3.2) for $p$ and $G_{i}$, and (3.10) for $D$. To obtain these, introduce

(3.25)\begin{equation} T_{nd}=\frac{\bar{R}\bar{\rho}}{\bar{p}}\,T,\quad\mathcal{D}_{nd}=\frac {\bar{\rho}^{2}}{\bar{D}\bar{p}}\,\mathcal{D},\end{equation}

where $\bar {R}$ is the characteristic gas constant. Rewriting (3.2) and (3.10) in terms of $T_{nd}$ and $\mathcal {D}_{nd}$, and omitting the subscript $_{nd}$, one can verify that the dimensional and non-dimensional versions of these expressions look exactly the same.

3.3. Non-dimensional parameters

Let the characteristic value of the Korteweg parameter be that of dry air, $\bar {K}=K_{22}$, given by (3.12). The characteristic values of the rest of the parameters are given by (2.26)–(2.27) and (3.3).

For common fluids, the room temperature range is much closer to the freezing point than the critical point (e.g. for water, the former is $0\,^{\circ }\mathrm {C}$ and the latter is $374\,^{\circ }\mathrm {C}$). In such cases, the interfacial thickness $\bar {z}$ is comparable to the liquid's molecular size – i.e. several ångströms – as confirmed by both experiments (e.g. Verde et al. Reference Verde, Bolhuis and Campen2012; Pezzotti et al. Reference Pezzotti, Galimberti and Gaigeot2017) and molecular simulations (e.g. Liu et al. Reference Liu, Harder and Berne2005; Pezzotti et al. Reference Pezzotti, Galimberti and Gaigeot2017; Dodia et al. Reference Dodia, Ohto, Imoto and Nagata2019; and references therein).

To obtain an estimate of $\bar {z}$ from within the DIM, one can use the particular case of an equilibrium interface. It is examined briefly in Appendix A.2 for the parameters specific to the water–air combination listed in Appendix B. The resulting profiles of the water and air densities are shown in figure 3, suggesting that the interfacial thickness is

(3.26)\begin{equation} \bar{z}\sim5 \unicode{x00C5}.\end{equation}

This value exceeds the water's molecular size by a factor of 2–3 (depending on how it is defined).

Figure 3.

A schematic illustrating the asymptotic structure of the solution. The curves $\rho _{1}(z)$ and $\rho _{2}(z)$ describe a flat equilibrium water/air interface at $T=25\,^{\circ }\mathrm {C}$, computed using the DIM (see Appendix A.2). The scale of the horizontal axes in both panels is logarithmic. The height of the non-shaded region is used as an estimate for the interfacial thickness $\bar {z}$.

According to (3.19), $\varepsilon$ depends on the time scale $\bar {t}$ – which, in turn, depends on the evaporation rate – hence $\varepsilon$ can be estimated only in application to a specific setting. Assuming that to be a drop of radius $r_{d}$, one can use the continuity of the advective mass flux to relate $\dot {r}_{d}$ to $\bar {w}$, we have

(3.27)\begin{equation} \dot{r}_{d}\rho_{1}^{(l)}\sim\bar{w}\rho_{2}^{(a)}. \end{equation}

Expressing $\bar {w}$ from this equality, substituting it into definition (3.19) of $\varepsilon$, and approximating $\dot {r}_{d}$ by $r_{d} /\bar {t}$, one obtains

(3.28)\begin{equation} \varepsilon=\sqrt{\frac{\rho_{2}^{(a)}\bar{z}}{\rho_{1}^{(l)}r_{d}}}. \end{equation}

Luckily, this expression no longer involves $\bar {t}$ (which is unknown without further calculations). Then for $r_{d}$ between 1 $\mathrm {\mu }$m and 1 nm, the above expression shows that $\varepsilon$ is between $0.771\times 10^{-3}$ and $2.44\times 10^{-2}$ – i.e. it is small.

To estimate $\delta$, estimate the diffusion coefficient by $\bar {D}=\mathcal {D}\bar {\rho }^{2}/\bar {p}$, where $\mathcal {D}$ is the standard vapour-in-the-air diffusivity. Then (3.18) and (3.20) yield

(3.29)\begin{equation} \delta=\frac{\bar{K}\bar{\rho}^{2}}{\mathcal{D}\bar{\eta}}.\end{equation}

Estimating the viscosity of air at $25\,^{\circ }\mathrm {C}$ using the results of Shang et al. (Reference Shang, Wu, Wang, Yang, Ye, Hu, Tao and He2019), one obtains

(3.30)\begin{equation} \bar{\eta}=4.20\times10^{{-}5}\ \mathrm{kg}\ \mathrm{m}^{{-}1}\ \mathrm{s}^{{-}1},\end{equation}

for which (3.29) yields $\delta \approx 1.83\times 10^{-8}$. Thus one can assume that

(3.31)\begin{equation} \delta\ll\varepsilon\ll1. \end{equation}

In addition to $\varepsilon$ and $\delta$, the problem involves two ‘hidden’ small parameters:

(3.32)\begin{equation} \epsilon^{(a/l)}=\frac{\rho_{2}^{(a)}}{\rho_{1}^{(l)}}\approx1.19\times 10^{{-}3},\quad\epsilon^{(v/a)}=\frac{\rho_{1}^{(v.sat)}}{\rho_{2}^{(a)} }\approx1.95\times10^{{-}2}. \end{equation}

The presence of four independent small parameters makes the analysis awkward: they arise in various combinations, and each time, one has to check whether these are large, small or order-one. In particular, the following estimates will be needed:

(3.33)\begin{gather} \frac{\varepsilon^{2}\epsilon^{(a/l)}\epsilon^{(v/a)}}{\delta}\ll \frac{\varepsilon^{2}\epsilon^{(a/l)}}{\delta}\approx3.84\times10^{{-}2} \ll1, \end{gather}

(3.34)\begin{gather}\frac{\delta}{\epsilon^{(v/a)}}\approx0.938\times10^{{-}6}\ll1. \end{gather}

Thus all of the above parameter combinations will be assumed small.

3.4. The boundary/matching conditions

When applied to the interface (inner region), (3.22)–(3.24) need boundary conditions at large $z$. These are to be obtained via matching the inner solution to those in the liquid and air (outer regions). The former is trivial: since the liquid is homogeneous and at rest, the interfacial solution should satisfy

(3.35)\begin{gather} \rho_{i}\rightarrow\rho_{i}^{(l)}\quad\text{as}\ z\rightarrow -\infty, \end{gather}

(3.36)\begin{gather}w\rightarrow0,\quad J\rightarrow0\quad\text{as}\ z\rightarrow -\infty. \end{gather}

Unlike the classical model where $\rho _{i}^{(l)}$ are external parameters, the DIM determines them as functions of the pressure and humidity of air. The DIM also describes dissolution of air in liquid – hence $\rho _{2}^{(l)}\neq 0$. This quantity is very small and thus unimportant both dynamically and thermodynamically.

The matching between the interface and air is less trivial: it can be carried out only after examining the large-distance behaviours of the inner solution and choosing the one that matches the outer solution. The properties of the latter are easy to guess: it should be such that $\rho _{1}\ll \rho _{2}$ (small vapour concentration in the air), and the variations of $\rho _{1}$ and $\rho _{2}$ should be such that the total pressure is spatially uniform (e.g. Ricci & Rocca Reference Ricci and Rocca1984; Mills & Coimbra Reference Mills and Coimbra2016).

To verify that (3.22)–(3.24) admit such an asymptotic solution, let

(3.37)\begin{align} \rho_{i}=\rho_{i}^{(+)}+\varepsilon\delta\rho_{i}^{\prime(+)}z+{O} ( \varepsilon^{2}\delta,\varepsilon\delta^{2}) ,\quad w=w^{(+)}+{O}(\varepsilon),\quad J=J^{(+)}+{O}(\varepsilon ),\end{align}

where $\rho _{i}^{(+)}$, $\rho _{i}^{\prime (+)}$, etc. depend on $t$ (but not on $z$) and may involve the hidden small parameters $\epsilon ^{(a/l)}$ and $\epsilon ^{(v/a)}$. It turns out that only $\rho _{i}^{(+)}$ is independent in these expansions, with the other coefficients being expressible through $\rho _{i}^{(+)}$.

To find how $w^{(+)}$ and $J^{(+)}$ depend on $\rho _{i}^{(+)}$, observe that to leading order, (3.24) can be integrated. Using boundary conditions (3.35)–(3.36) to fix the constants of integration, one obtains

(3.38)\begin{equation} w={-}\dot{Z}\left( \frac{\rho_{1}^{(l)}+\rho_{2}^{(l)}}{\rho_{1}+\rho_{2} }-1\right) +{O}(\varepsilon),\quad J={-}\dot{Z}\,\frac{\rho_{1}^{(l)}\rho_{2}-\rho_{1}\rho_{2}^{(l)}}{\rho_{1}+\rho_{2}}+{O} (\varepsilon),\end{equation}

which entails

(3.39)\begin{equation} w^{(+)}={-}\dot{Z}\left( \frac{\rho_{1}^{(l)}+\rho_{2}^{(l)}}{\rho_{1}^{(+)}+\rho_{2}^{(+)}}-1\right) ,\quad J^{(+)}={-}\dot{Z}\,\frac{\rho_{1}^{(l)}\rho_{2}^{(+)}-\rho_{1}^{(+)}\rho_{2}^{(l)}}{\rho_{1}^{(+)}+\rho_{2}^{(+)}}. \end{equation}

To calculate $\rho _{i}^{\prime (+)}$, one should substitute the expansions of $w$ and $J$ into (3.22)–(3.23). Since air is a diluted fluid, the effective viscosity does not depend on the density (e.g. Ferziger & Kaper Reference Ferziger and Kaper1972), i.e.

(3.40)\begin{equation} \eta\sim\eta_{0}\quad\text{as}\ \rho_{i}\rightarrow0,\end{equation}

where $\eta _{0}$ depends only on $T$. This property implies that when expressions (3.38) are substituted into (3.22)–(3.23), the linear term in the expansion for $w$ cancels out. The quadratic term does not, but its contribution is much smaller than that of $J^{(0)}$ (due to estimates (3.33)), and to leading order, one obtains the following algebraic equations for $\rho _{i}^{\prime (+)}$:

(3.41)\begin{gather} \sum_{i}\left( \frac{\partial G_{1}}{\partial\rho_{i}}\right)^{(+)}\rho _{i}^{\prime(+)} ={-}\frac{\rho_{2}^{(+)}J^{(+)}}{D^{(+)}( \rho_{1}^{(+)}+\rho_{2}^{(+)}) }, \end{gather}

(3.42)\begin{gather}\sum_{i}\left( \frac{\partial G_{2}}{\partial\rho_{i}}\right)^{(+)}\rho _{i}^{\prime(+)} =\frac{\rho_{1}^{(+)}J^{(+)}}{D^{(+)}(\rho_{1}^{(+)}+\rho_{2}^{(+)}) }, \end{gather}

where $D^{(+)}=D(\rho _{1}^{(+)},\rho _{2}^{(+)},T)$ and

(3.43)\begin{equation} \left( \frac{\partial G_{j}}{\partial\rho_{i}}\right)^{(+)}=\frac{\partial G_{j}(\rho_{1}^{(+)},\rho_{2}^{(+)},T)}{\partial\rho_{i}^{(+)}}. \end{equation}

Given the smallness of the air-to-liquid and vapour-to-air density ratios (3.32), (3.41)–(3.42) can be simplified. Using the diluted-fluid expressions (3.2) and (3.10) for $G_{i}$ and $D$, respectively, and keeping the leading order only, one obtains

(3.44)\begin{gather} \rho_{1}^{\prime(+)} =\frac{\rho_{1}^{(l)}\dot{Z}}{\mathcal{D} }, \end{gather}

(3.45)\begin{gather}\rho_{2}^{\prime(+)} ={-}\rho_{1}^{\prime(+)}\,\frac{R_{1}}{R_{2}}. \end{gather}

Equality (3.44) is the desired matching condition; it can also be viewed as the relationship of the speed $\dot {Z}$ of the interface to the vapour flux $\mathcal {D}\rho ^{\prime (+)}$, whereas (3.45) is the standard condition of constant pressure for diffusion in diluted fluids. As mentioned before, both conditions are obvious physically, but it is still important to verify that they are intrinsic to the DIM formalism.

3.5. The interfacial asymptotic region

Three asymptotic zones exist within the inner (interfacial) region,

(3.46)\begin{equation} \left.\begin{array}{ll@{}} \text{Zone 1\ (mostly air):} & \rho_{1}\ll\rho_{2}\ll\rho^{(l)},\\ \text{Zone 2\ (air--vapour mixture):} & \rho_{1}\sim \rho_{2}\ll\rho^{(l)},\\ \text{Zone 3\ (mostly liquid):} & \rho_{2}\ll\rho_{1}\sim \rho^{(l)}. \end{array}\right\} \end{equation}

Note that Zone 1 is not part of the air region: in the latter, the air density is nearly uniform, $\rho _{2}\approx \rho _{2}^{(a)}$, whereas in the former, variations of $\rho _{2}$ are comparable to $\rho _{2}^{(a)}$.

A schematic illustrating the asymptotic structure of the problem can be found in figure 3 (which is computed for the equilibrium interface, i.e. that where the relative humidity of air is $H=1$). Its most counter-intuitive features are the local maximum and minimum of the air density $\rho _{2}(z)$ in Zone 2. Note that they both disappear in the limit $K_{12}\rightarrow 0$ (see figure 11 of Benilov Reference Benilov2023a). One can thus assume that the maximum of $\rho _{2}$ is created by the van der Waals force exerted by the bulk of the liquid by attracting a certain amount of air towards the interface. The existence of the minimum of $\rho _{2}$ is harder to understand, but it is dynamically unimportant anyway. The air density is extremely small there – even smaller than that of the air dissolved in the liquid – and they both have a minuscule effect on the dynamics of the rest of the fluid.

It turns out that only Zone 1 contributes significantly to the evaporation rate, whereas the contributions of Zones 2–3 are negligible. To understand why, recall that the fluid density in the latter two zones exceeds that of the vapour component of air by an order of magnitude. As a result, Zones 2–3 are close to equilibrium (the fact that the ‘outside’ air is undersaturated cannot perturb them too much). Thus Zone 1 is the only region out of equilibrium – hence it is the only place where irreversible processes, such as evaporation, can occur.

The calculations associated with Zones 2–3 can be avoided via a certain workaround developed by Benilov (Reference Benilov2022b, Reference Benilov2023b) for evaporation of pure fluids. The workaround involves the following two steps.

1. (i) Equations (3.22)–(3.24) can be rearranged in such a way that the velocity $\dot {Z}$ of the interface is expressed in terms of the fluid's known parameters (relative humidity, etc.) and a certain integral of the solution.

2. (ii) The contribution of Zone 1 to this integral exceeds the contributions of Zones 2–3 by an order of magnitude. Thus to evaluate this integral asymptotically, it is sufficient to find the solution only in Zone 1.

3.5.1. Calculation of $\dot {Z}$

Recall that matching conditions (3.44)–(3.45) emerged from the governing equations at the order of ${O}(\varepsilon \delta )$. Finding $\dot {Z}$ does not require going that far, as the terms in (3.22)–(3.23) that involve $\dot {Z}$ are ${O} (\varepsilon )$. The terms ${O}(\delta )$ in these equations can also be omitted (as suggested by the estimates in § 3.3).

Substitute expressions (3.38) for $w$ and $J$ into (3.22)–(3.23), and, keeping the terms up to ${O} (\varepsilon )$, obtain

(3.47)\begin{gather} \frac{\partial}{\partial z}\left( \sum_{j}K_{1j}\,\frac{\partial^{2}\rho_{j} }{\partial z^{2}}-G_{1}\right) =\frac{\varepsilon\dot{Z}}{\rho_{1}+\rho_{2} }\,\frac{\partial}{\partial z}\left[ \eta\,\frac{\partial}{\partial z}\left( \frac{\rho_{1}^{(l)}+\rho_{2}^{(l)}}{\rho_{1}+\rho_{2}}\right) \right] , \end{gather}

(3.48)\begin{gather}\frac{\partial}{\partial z}\left( \sum_{j}K_{2j}\,\frac{\partial^{2}\rho_{j} }{\partial z^{2}}-G_{2}\right) =\frac{\varepsilon\dot{Z}}{\rho_{1}+\rho_{2} }\,\frac{\partial}{\partial z}\left[ \eta\,\frac{\partial}{\partial z}\left( \frac{\rho_{1}^{(l)}+\rho_{2}^{(l)}}{\rho_{1}+\rho_{2}}\right) \right] . \end{gather}

Keeping the same accuracy in boundary condition (3.37), one obtains

(3.49)\begin{equation} \rho_{i}\rightarrow\rho_{i}^{(+)}\quad\text{as}\ z\rightarrow +\infty.\end{equation}

Now, integrate (3.47) and (3.48) with respect to $z$ from $-\infty$ to $+\infty$. Recalling boundary conditions (3.35) and (3.49), one obtains

(3.50)\begin{equation} G_{i}^{(l)}-G_{i}^{(+)}={-}\varepsilon\dot{Z}A,\end{equation}

where

(3.51) \begin{equation} A=\int_{-\infty}^{+\infty}\frac{1}{\rho_{1}+\rho_{2}}\,\frac{\partial}{\partial z}\left[ \frac{( \rho_{1}^{(l)}+\rho_{2}^{(l)}) \eta}{\left( \rho_{1}+\rho_{2}\right)^{2}}\,\frac{\partial( \rho_{1}+\rho_{2}) }{\partial z}\right] \mathrm{d}z.\end{equation}

Next, consider the combination

(3.52)\begin{equation} \int_{-\infty}^{+\infty}[ \rho_{1}\times \text{(3.47)} +\rho_{2}\times \text{(3.48)} ] \,\mathrm{d}z. \end{equation}

Integrating by parts the terms involving the third derivatives, take into account boundary conditions (3.35) and (3.49), and recall identity (3.1) to obtain

(3.53)\begin{equation} p^{(l)}-p^{(+)}=0.\end{equation}

Thus the pressure in the liquid coincides with that in the air (as expected physically).

Equalities (3.50) and (3.53) relate the (unknown) liquid densities $\rho _{i}^{(l)}$ and the velocity $\dot {Z}$ of the interface to the air parameters (marked with $^{(+)}$). In addition, $\dot {Z}$ is also related to the air parameters by equality (3.44). If two of these three equalities are used to eliminate $\rho _{i}^{(l)}$ and $\dot {Z}$, then the remaining one will inter-relate the air parameters and thus turn into the generalised Maxwell boundary condition (GMBC).

To eliminate $\rho _{i}^{(l)}$ from equalities (3.50) and (3.53), recall that equations of state of liquids at normal conditions are typically such that even a minuscule change of the density gives rise to a huge change of the pressure. This means effectively that $\rho _{i}^{(l)}$ is close to its saturated value, $\rho _{i}^{(l)}\approx \rho _{i}^{(l.sat)}$ – so that (3.50) becomes

(3.54)\begin{equation} G_{i}^{(l.sat)}-G_{i}^{(+)}={-}\varepsilon\dot{Z}A. \end{equation}

In the state of saturation, the chemical potentials of the liquid and vapour are equal (see Appendix A.2) – hence one can change in the above equality $G_{1}^{(l.sat)}\rightarrow G_{1}^{(v.sat)}$. Then, using for both $G_{1}^{(v.sat)}$ and $G_{1}^{(+)}$ the diluted-fluid expression (3.2), one obtains

(3.55)\begin{equation} \rho_{1}^{(l.sat)}R_{1}T\ln\frac{\rho_{1}^{(+)}}{\rho_{1}^{(v.sat)} }=\varepsilon\mathcal{D}\rho_{1}^{\prime(+)}A.\end{equation}

This is, essentially, the GMBC.

With $A$ calculated (see the next subsubsection), equality (3.55) inter-relates the near-interface vapour density $\rho _{1}^{(+)}$ and its gradient $\rho _{1}^{\prime (+)}$. If $\varepsilon =0$, then (3.55) yields $\rho _{1}^{(+)}=\rho _{1}^{(v.sat)}$, which is the classical Maxwell boundary condition. If, however, $\varepsilon$ is small but not zero, then $\rho _{1}^{(+)}$ can differ significantly from $\rho _{1}^{(v.sat)}$ – because $A$ can actually be large, so that $\varepsilon A\sim 1$. This does not violate the employed asymptotic approach: it requires only $\varepsilon \ll 1$, which makes the time derivatives in (3.24) small and thus justifies the quasi-steady approximation.

3.5.2. Zone 1 of the interfacial asymptotic region

The coefficient $A$ arises in all problems where the DIM is applied to evaporation or condensation – but, so far, it has been calculated only for pure fluids and the case where the relative humidity is close to unity (Benilov Reference Benilov2020, Reference Benilov2022a,Reference Benilovb, Reference Benilov2023b). The calculation of $A$ for a binary mixture and arbitrary humidity is more complex technically, but the underlying idea is the same.

It is based on the observation that the main contribution to (3.51) comes from the zone where $\rho _{1}+\rho _{2}$ is small, but its derivatives are relatively large – i.e. near the air region, but not in it. This means effectively that (3.51) can be calculated asymptotically by considering Zone 1 only (see its definition at the beginning of § 3.5).

Before presenting the scaling of Zone 1, it is instructive to apply the ‘obvious’ approximations – i.e. take advantage of the diluted-fluid expressions (3.2), (3.10) and (3.40), plus neglect $\rho _{1}$ and $\rho _{2}^{(l)}$ if they appear next to $\rho _{2}$ and $\rho _{1}^{(l)}$, respectively. As a result, (3.47)–(3.48) and (3.51) become

(3.56)\begin{gather} \frac{\partial}{\partial z}\left( K_{12}\,\frac{\partial^{2}\rho_{2}}{\partial z^{2}}-R_{1}T\ln\rho_{1}\right) =\varepsilon\rho_{1}^{(l)}\dot{Z}\,\frac {\eta_{0}}{\rho_{2}}\,\frac{\partial^{2}}{\partial z^{2}}\left( \frac{1} {\rho_{2}}\right) , \end{gather}

(3.57)\begin{gather}\frac{\partial}{\partial z}\left( K_{22}\,\frac{\partial^{2}\rho_{2}}{\partial z^{2}}-R_{2}T\ln\rho_{2}\right) =\varepsilon\rho_{1}^{(l)}\dot{Z}\,\frac {\eta_{0}}{\rho_{2}}\,\frac{\partial^{2}}{\partial z^{2}}\left( \frac{1} {\rho_{2}}\right) , \end{gather}

(3.58)\begin{gather}A=\eta_{0}\rho_{1}^{(l)}\int_{-\infty}^{+\infty}\frac{1}{\rho_{2}}\,\frac{\partial}{\partial z}\left( \frac{1}{\rho_{2}^{2}}\,\frac{\partial \rho_{2}}{\partial z}\right) \mathrm{d}z. \end{gather}

The most general asymptotic limit is when the three terms in (3.57) are of the same order, which implies the rescaling

(3.59)\begin{gather} z=z_{0}+\left( \frac{K_{22}\rho_{2}^{(+)}}{R_{2}T}\right)^{1/2} z_{new},\quad \rho_{i}=\rho_{i}^{(+)}\left( \rho_{i}\right)_{new} ,\quad \dot{Z}={-}\frac{E_{0}}{\varepsilon\rho_{1}^{(l)}}\,\xi, \end{gather}

(3.60)\begin{gather}A=\frac{\rho_{1}^{(l)}R_{2}T}{E_{0}}\,A_{new}, \end{gather}

where $z_{0}$ is the ‘location’ of Zone 1, $E_{0}$ is given by (1.3) of the Introduction, and the minus in the expression for $\dot {Z}$ is introduced to make $\xi$ positive (recall that evaporation corresponds to $\dot {Z}<0$). Since the problem under consideration is translationally invariant, the actual value of $z_{0}$ does not feature in, and is not needed for, any leading-order results. One does not need (3.56) either, because now $A$ depends only on $\rho _{2}$.

Substitution of (3.59)–(3.60) into (3.57), boundary condition (3.49), and (3.58) yield (with the subscript $_{new}$ omitted)

(3.61)\begin{gather} \frac{\partial}{\partial z}\left( \frac{\partial^{2}\rho_{2}}{\partial z^{2} }-\ln\rho_{2}\right) ={-}\frac{\xi}{\rho_{2}}\,\frac{\partial^{2}}{\partial z^{2}}\left( \frac{1}{\rho_{2}}\right) , \end{gather}

(3.62)\begin{gather}\rho_{2}\rightarrow1\quad\text{as}\ z\rightarrow+\infty, \end{gather}

(3.63)\begin{gather}A=\int_{-\infty}^{+\infty}\frac{1}{\rho_{2}}\,\frac{\partial}{\partial z}\left( \frac{1}{\rho_{2}^{2}}\,\frac{\partial\rho_{2}}{\partial z}\right) \mathrm{d}z. \end{gather}

As shown in Appendix C, boundary condition (3.62) is sufficient to uniquely fix the solution of (3.61), making it unnecessary to examine the other asymptotic zones of the interfacial region.

It can be derived readily from (3.61) that

(3.64)\begin{equation} \rho_{2}\sim z^{2}\ln({-}z) \quad\text{as}\ z\rightarrow -\infty. \end{equation}

The growth of $\rho _{2}(z)$ as $z\rightarrow -\infty$ in Zone 1 agrees with the fact that this function has a local maximum in Zone 2 (see figure 3). Asymptotics (3.64) also guarantees that (3.63) converges at the lower limit.

The function $A(\xi )$ was computed numerically using the algorithm described in § C.2; the result is plotted in figure 4. For practical use, however, it is more convenient to use the approximation formula (1.4) (for more details, see § C.3).

Figure 4.

The function $A(\xi )$. Note that in the right-hand half of the figure, the scale of $\xi$ is logarithmic.

3.5.3. Summarising the GMBC for flat interfaces

Recall that the dimensional variables can be recovered by setting $\varepsilon =\delta =1$. Then equalities (3.44) and (3.59) yield a dimensional relationship of the near-interface air density $\rho _{1}^{(+)}$ and its gradient $\rho _{1}^{\prime (+)}$,

(3.65)\begin{equation} \xi=\frac{-\mathcal{D}\rho_{1}^{\prime(+)}}{E_{0}},\end{equation}

and (3.55) yields

(3.66)\begin{equation} \rho_{1}^{(+)}=\rho_{1}^{(v.sat)}\exp\left[ -\frac{R_{2}}{R_{1}}\,\xi\, A(\xi)\right] .\end{equation}

Equalities (3.65) and (3.66) inter-relate the near-surface vapour density $\rho _{1}^{(+)}$ and its normal derivative $\rho _{1}^{\prime (+)}$, hence they constitute the generalised Maxwell boundary condition (for flat interfaces), as required.

3.5.4. Curved interfaces

Interfacial curvature influences evaporation via the so-called Kelvin effect, i.e. by changing the effective saturated vapour pressure near the liquid's (curved) surface. This effect can be incorporated readily into the analysis under the assumption that the curvature radius is much larger than the interfacial thickness. This has been done, and the result has turned out to be physically obvious: the Kelvin effect adds a new term to the pressure condition (3.53),

(3.67)\begin{equation} p^{(l)}-p^{(+)}=\varepsilon C\sigma, \end{equation}

where $C$ is the interfacial curvature, and $\sigma$ is the surface tension. (The DIM expression for $\sigma$ is given by (B8) of Appendix B.)

The extra term in the pressure equation entails an extra term in the GMBC: instead of (3.66), it becomes (1.1) of the Introduction, whereas definition (3.65) of the relative evaporation rate $\xi$ coincides with (1.2) if one sets $\rho _{1}^{\prime (+)}=( \boldsymbol {n} \boldsymbol {\cdot }\boldsymbol {\nabla }\rho _{1})^{(+)}$.