Homepage: ]https://staff.ul.ie/eugenebenilov/

# Nonisothermal evaporation

E. S. Benilov<EMAIL_ADDRESS>[ Department of Mathematics and
Statistics, University of Limerick, Limerick V94 T9PX, Ireland

###### Abstract

Evaporation of a liquid layer on a substrate is examined without the often-
used isothermality assumption – i.e., temperature variations are accounted
for. Qualitative estimates show that nonisothermality makes the evaporation
rate depend on the conditions the substrate is maintained at. If it is
thermally insulated, evaporative cooling dramatically slows evaporation down;
the evaporation rate tends to zero with time and cannot be determined by
measuring the external parameters only. If, however, the substrate is
maintained at a fixed temperature, the heat flux coming from below sustains
evaporation at a finite rate – deducible from the fluid’s characteristics,
relative humidity, and the layer’s depth (whose importance has not been
recognized before). The qualitative predictions are quantified using the
diffuse-interface model applied to a liquid evaporating into its own vapor.

##  I Introduction and preliminary estimates

Evaporation of liquids has been studied for over a century, since the
pioneering work of James Clerk Maxwell [Maxwell77, Maxwell90], and it is still
studied now – both theoretically (e.g., Refs.
[SobacTalbotHautRednikovColinet15, TalbotSobacRednikovColinetHaut16, Sazhin17,
WilliamsKarapetsasMamalisSefianeMatarValluri20, FinneranGarnerNadal21,
DallabarbaWangPicano21]) and experimentally (e.g., Refs.
[JakubczykKolwasDerkachovKolwasZientara12, HolystLitniewskiJakubczyk17,
SchweiglerBensaidSeifritzSelzerNestler17, PoosVarju20, VarjuPoos22]) – as
numerous issues have yet to be resolved.

Consider, for example, a flat liquid layer. It is generally believed that it
evaporates at a steady rate depending on the liquid’s parameters (temperature,
heat of vaporization, etc.) and the humidity of air. There are numerous
measurements of evaporation rates; a recent review of this work in application
to water can be found in Refs. [PoosVarju20, VarjuPoos22].

Fig. 1 shows a selection111Fig. 1 shows those empiric formulae listed in Table
1 of Ref. [VarjuPoos22] that involve only the relative humidity and
characteristics uniquely related to the temperature, such as the saturated
pressure and vaporization heat of water (calculated using Refs.
[WagnerPruss02, Hendersonsellers84], respectively). One of the formulae
includes also the wind speed which was set to zero. All the other empiric
results cited in Ref. [VarjuPoos22] involve further parameters (e.g., the
horizontal scale of the vessel), making a comparison with the low-parameter
formulae impossible. of measurements of the evaporation rate $E$, for water
evaporating into still air, as a function of the temperature $T$ within a
“room temperature” range. There is evident discord in these results,
suggesting that important factors vary from experiment to experiment.

The present paper identifies at least some of these factors. It is shown that
nonisothermal effects – e.g., the heat exchange between the liquid and
substrate (and side walls, if any) – can make $E$ depend on the distance
between the interface and substrate, and the material the latter is made of.

Figure 1: Evaporation rate vs. temperature, according to various empiric
formulae. Curves (1)–(6) correspond to Refs. [MoghimanJodat07,
MackayVanwesenbeeck14, SharpleyBoelter38, Shah12, BoelterGordonGriffin46,
TangPaukenJeterAbdelkhalik93] , respectively. In all cases, the relative
humidity is 50%.

To illustrate the importance of heat exchange with the boundaries, consider an
amount of liquid in a thermally insulated vessel – and let half of the liquid
evaporate. The temperature of the remaining half decreases due to evaporative
cooling – and the size of the decrease is easy to estimate. Assuming for
simplicity that vaporization heat $\Delta h$ and heat capacity $c_{P}^{(l)}$
of the liquid do not change significantly with $T$, one can approximate the
temperature decrease by

$\Delta T=\frac{\Delta h}{c_{P}^{(l)}}.$ (1)

For $\Delta h=2442\,\mathrm{J\,g}^{-1}$ and
$c_{P}^{(l)}=4.182\,\mathrm{J\,g}^{-1}\mathrm{K}^{-1}$ (which correspond to
water at $25^{\circ}\mathrm{C}$ [LindstromMallard97]), Eq. (1) yields a
somewhat unexpected result:

$\Delta T\approx 584\,\mathrm{K}.$

In reality, however, evaporation of liquid in an insulated vessel slows down
to a virtual standstill well before it half-evaporates. Since the dependence
of $E$ on $T$ at normal conditions is typically exponential, even a moderate
temperature decrease can reduce the evaporation rate by an order of magnitude.

Alternatively, let the vessel’s walls and bottom be kept at a fixed
temperature (such a setting has probably more applications). In this case, the
energy loss to vaporization is replenished by the incoming heat flux, which
can be readily calculated,

$-\kappa^{(l)}T^{\prime}=E\,\Delta h,$ (2)

where $\kappa^{(l)}$ is the liquid’s thermal conductivity and the temperature
gradient $T^{\prime}$ can be expressed through the temperature difference
$\Delta T$ between the interface and the nearest boundary, and the
corresponding distance $D$,

$T^{\prime}=-\frac{\Delta T}{D}.$ (3)

To determine $\Delta T$ for water at $25^{\circ}\mathrm{C}$, set
$\kappa^{(l)}=0.6065\,\mathrm{W\,m}^{-1}\mathrm{K}^{-1}$ [LindstromMallard97]
and

$E=0.025361\,\mathrm{g\,m}^{-2}\mathrm{s}^{-1},$ (4)

which is the average of empiric curves 2–6 in Fig. 1 (curve 1 cannot be used,
as $25^{\circ}\mathrm{C}$ is not in its range). With these values, Eqs.
(2)–(3) yield

$\Delta T=\frac{E\,\Delta
h}{\kappa^{(l)}}D\approx\frac{1\,\mathrm{K}}{1\,\mathrm{cm}}\times D.$ (5)

Evidently, this estimate is both qualitatively and quantitatively different
from that for insulated substrates.

The difference between insulated and fixed-$T$ vessels demonstrates the
importance of heat fluxes from boundaries and, generally, nonisothermal
effects. In this work, they are explored using the simplest setting:
evaporation of a liquid into _its own undersaturated vapor_. It is described
by a relatively simple model which does not include the diffusive mass flux
(pure fluids do not diffuse). Evaporation in this case occurs via advection
[Benilov22a], but heat conduction is similar to that in mixtures.

The described setting will be examined using the so-called diffuse interface
model (DIM). It was proposed in 1901 by Diederik Korteweg [Korteweg01] and has
been used since then in thousands of papers and for tens of applications (some
of this work is reviewed in Ref. [Benilov23a]). The DIM is particularly suited
to the problem at hand: it describes both liquid and vapor, as well as the
interfacial dynamics – as opposed to models built of ‘blocks’ describing one
item each. The use of the DIM is convenient but not crucial, however, as
nonisothermal effects can be introduced into any good model of evaporation.

In Sec. II of the present paper, the problem is formulated mathematically. In
Sec. III, evaporation is examined under the assumption of isothermality. This
case will be used as a yardstick for the full problem examined in IV. Other
effects potentially explaining the discord among the empiric curves in Fig. 1
are discussed in Sec. V. In Sec. VI, the results are summarized, plus it is
clarified there when the heat flux from air is weak and its effect on
evaporation, negligible (which is of interest in a broader context, not just
for the present work).

##  II Formulation

###  II.1 Thermodynamics

Thermodynamic properties of a fluid can be described by the dependence of its
specific (per unit mass) internal energy $e$ and specific entropy $s$ on the
density $\rho$ and temperature $T$. The functions $e(\rho,t)$ and $s(\rho,t)$
are supposed to be constrained by the Gibbs relation; if written in terms of
$e$ and $s$, it takes the form

$\frac{\partial e}{\partial T}=T\frac{\partial s}{\partial T}.$ (6)

The fluid’s equation of state, or the expression for the pressure $p$, is
given by

$p=\rho^{2}\left(\frac{\partial e}{\partial\rho}-T\frac{\partial
s}{\partial\rho}\right),$ (7)

and the chemical potential, or specific free energy, by

$G=e+\frac{p}{\rho}-Ts.$ (8)

It follows from (6)–(8) that

$\frac{\partial p}{\partial\rho}=\rho\dfrac{\partial G}{\partial\rho},$ (9)
$\dfrac{\partial p}{\partial T}=-\rho^{2}\frac{\partial s}{\partial\rho}.$
(10)

These two identities enable one to replace $p$ with $G$ or $s$, which happens
to be convenient in the problem at hand.

Define the heat capacity at constant volume,

$c_{V}=\frac{\partial e}{\partial T},$ (11)

the specific enthalpy,

$h=e+\frac{p}{\rho},$ (12)

and the heat capacity at constant pressure,

$c_{P}=\left(\frac{\partial h}{\partial T}\right)_{p=\operatorname{const}}.$

Expressing the derivative at constant pressure via the partial derivatives
with respect $\rho$ and $T$, and recalling Eq. (12), one obtains

$c_{P}=\frac{\partial e}{\partial T}-\left(\dfrac{\partial
e}{\partial\rho}-\frac{p}{\rho^{2}}\right)\dfrac{\partial p}{\partial
T}\left(\dfrac{\partial p}{\partial\rho}\right)^{-1}.$ (13)

###  II.2 Hydrodynamics

Consider a liquid layer on a horizontal solid substrate, and vapor above the
liquid. If the vapor in undersaturated, the liquid evaporates, giving rise to
a vertical flow. This setting is characterized by the velocity $w(z,t)$,
density $\rho(z,t)$, and temperature $T(z,t)$, where $z$ is the vertical
coordinate and $t$, the time.

####  II.2.1 Governing equations

It can be safely assumed that the Reynolds number associated with evaporation
is small, so that the contribution of inertia to the balance of momentum is
negligible. Thus, Stokes-flow (slow-flow) approximation can be employed.

Using the diffuse-interface model (DIM), one can write the Stokes-flow version
of the hydrodynamic equations in the form

$\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial
z}\underset{\text{mass flux}}{\underbrace{\left(\rho w\right)}}=0,$ (14)
$\underset{\text{pressure gradient}}{\underbrace{\frac{\partial p}{\partial
z}}}~{}=~{}\underset{\text{viscous
stress}}{\underbrace{\frac{\partial}{\partial z}\left(\eta\frac{\partial
w}{\partial z}\right)}}~{}~{}+\underset{\text{van der Waals
force}}{\underbrace{K\rho\frac{\partial^{3}\rho}{\partial z^{3}}}},$ (15)
$\frac{\partial\left(\rho e\right)}{\partial t}+\frac{\partial}{\partial
z}\underset{\text{heat flux}}{\underbrace{\left[w\left(\rho
e+p-\eta\frac{\partial w}{\partial z}\right)-\kappa\frac{\partial T}{\partial
z}\right]}}\\\ =\underset{\text{rate of work done by van der Waals
force}}{\underbrace{w\,K\rho\frac{\partial^{3}\rho}{\partial z^{3}}}},$ (16)

where the effective viscosity $\eta$ is related to the shear viscosity
$\mu_{s}$ and and bulk viscosity $\mu_{b}$ by

$\eta=\frac{4}{3}\mu_{s}+\mu_{b},$

$\kappa$ is the thermal conductivity and $K$, the so-called Korteweg
parameter. The expression for the van der Waals force comes from the DIM,
otherwise (14)–(16) are standard equations of compressible Stokes-flow
hydrodynamics. The three-dimensional, non-Stokes flow versions of Eqs.
(14)–(16) were derived in Ref. [Giovangigli20] from the Enskog–Vlasov kinetic
theory, and in Ref. [GalloMagalettiCasciola21] via nonequilibrium
thermodynamics. In this paper, a brief derivation of the DIM expression for
the van der Waals force is given in Appendix A.

Note that $\eta$ and $\kappa$ depend generally on $\rho$ and $T$, whereas $K$
is a constant. Its value is related to, and can be deduced from, the fluid’s
surface tension – for water, for example, $K\approx 1.9\times
10^{-17}\mathrm{m}^{7}\mathrm{s}^{-2}\mathrm{kg}^{-1}$ [Benilov23a].

####  II.2.2 Boundary conditions far above the interface

Assume that, far above the liquid–vapor interface, the viscous stress is zero,

$\frac{\partial w}{\partial z}\rightarrow 0\qquad\text{as}\qquad
z\rightarrow+\infty,$ (17)

and the vapor density and temperature tend to certain values,

$T\rightarrow T^{(v)}\qquad\text{as}\qquad z\rightarrow+\infty,$ (18)
$\rho\rightarrow\rho^{(v)}\qquad\text{as}\qquad z\rightarrow+\infty.$ (19)

For evaporation to occur, $\rho^{(v)}$ should be smaller than the saturated
vapor density $\rho^{(v.sat)}$ – which is determined, together with the
matching saturated liquid density $\rho^{(l.sat)}$, by the Maxwell
construction:

$p(\rho^{(l.sat)},T^{(v)})=p(\rho^{(v.sat)},T^{(v)}),$ (20)
$G(\rho^{(l.sat)},T^{(v)})=G(\rho^{(v.sat)},T^{(v)}),$ (21)

where it is implied that $\rho^{(l.sat)}\geq\rho^{(v.sat)}$.

Physically, the equalities of the pressure and chemical potential in the two
phases guarantee the mechanical and thermodynamic equilibria of the interface,
respectively. Mathematically, conditions (20)–(21) can be derived from the DIM
(to be elaborated later) or from any other good model describing a static flat
interface in an unbounded space.

Require also that the saturated vapor and liquid be thermodynamically stable,
which amounts to

$\left(\frac{\partial p}{\partial\rho}\right)_{\rho=\rho^{(v.sat)}}\geq
0,\qquad\left(\frac{\partial p}{\partial\rho}\right)_{\rho=\rho^{(l.sat)}}\geq
0,$ (22)

i.e., an increase in $\rho$ does not decrease the pressure. The Maxwell
construction (20)–(21) and requirements (22) uniquely define $\rho^{(v.sat)}$
and $\rho^{(l.sat)}$ as functions of $T$.

Before proceeding further, it is convenient to integrate the momentum equation
(15) and fix the constant of integration via boundary conditions (17)–(19),
which yields

$p=\eta\frac{\partial w}{\partial
z}+K\left[\rho\frac{\partial^{2}\rho}{\partial
z^{2}}-\frac{1}{2}\left(\frac{\partial\rho}{\partial
z}\right)^{2}\right]+p^{(v)},$ (23)

where $p^{(v)}=p(\rho^{(v)},T^{(v)})$.

####  II.2.3 Boundary conditions at the substrate

Let the substrate be located at $z=0$ and write the no-flow boundary condition
in the form

$w=0\qquad\text{at}\qquad z=0.$ (24)

Two different conditions for $T$ will be examined: one assuming that the
substrate is kept at a fixed temperature,

$T=T_{0}\qquad\text{at}\qquad z=0,$ (25)

and another corresponding to a thermally insulated substrate,

$\frac{\partial T}{\partial z}=0\qquad\text{at}\qquad z=0.$ (26)

Due to the presence of higher-order derivatives of $\rho$ in expression for
the van der Waals force, a separate boundary condition is required for the
density. Several different versions of such are used in the literature (e.g.,
Refs. [Seppecher96, PismenPomeau00, GalloMagalettiCasciola21]), describing
slightly different models of the fluid–substrate interaction. The specific
form of this boundary condition is important only if the thickness of the
liquid–vapor interface is comparable to the liquid layer’s depth – which is,
obviously, not the case for a _macroscopic_ layer considered in this work.

Thus, the simplest version of the boundary condition for $\rho$ will be used –
the one suggested in Ref. [PismenPomeau00],

$\rho=\rho_{0}\qquad\text{at}\qquad z=0,$ (27)

where $\rho_{0}$ characterizes the fluid–substrate interaction. This condition
reflects the balance of forces affecting fluid molecules adjacent to the
substrate, and it can be derived under the same assumptions as the DIM itself
[Benilov20a].

Note that none of the conclusions reported in this paper depends on the
specific value of $\rho_{0}$.

####  II.2.4 A flat interface in an unbounded fluid

In a sufficiently deep layer, the interface is not affected by the substrate.
Boundary condition (24) can be moved to minus-infinity,

$w\rightarrow 0\qquad\text{as}\qquad z\rightarrow-\infty,$ (28)

and conditions (25)–(27), replaced with

$T\rightarrow T^{(l)}\qquad\text{as}\qquad z\rightarrow-\infty,$ (29)
$\rho\rightarrow\rho^{(l)}\qquad\text{as}\qquad z\rightarrow-\infty.$ (30)

Note that Eq. (23) and boundary conditions (28)–(30) imply that the liquid’s
temperature $T^{(l)}$ and density $\rho^{(l)}$ are not entirely arbitrary, but
are related to the vapor parameters by

$p(\rho^{(l)},T^{(l)})=p(\rho^{(v)},T^{(v)}).$ (31)

Mathematically, this constraint is a result of the Stokes-flow approximation:
if the time derivative in the momentum equation were retained, (31) would not
hold. Physically, the constraint suggests that an adjustment of pressure
should occur (via fast acoustic waves) before the flow becomes truly slow.

###  II.3 Nondimensionalization

To nondimensionalize Eqs. (14)–(16), introduce characteristic scales of
pressure $P$, density $\varrho$, and viscosity $\bar{\eta}$. Using these
parameters, one can define a velocity scale and a spatial scale,

$V=\frac{Pl}{\bar{\eta}},\qquad l=\sqrt{\frac{K\varrho^{2}}{P}}.$

As seen later, this choice of $V$ and $l$ corresponds to an asymptotic regime
where the pressure gradient, viscous stress, and van der Waals force in Eq.
(15) are all of the same order. The interfacial scale $l$ will be referred to
as ‘microscopic’, and the depth of the whole liquid layer, ‘macroscopic’.

The following nondimensional variables will be used:

$z_{nd}=\frac{z}{l},\qquad t_{nd}=\frac{V}{l}t,$
$w_{nd}=\frac{w}{V},\qquad\rho_{nd}=\frac{\rho}{\varrho},\qquad
T_{nd}=\frac{\varrho RT}{P},$ $e_{nd}=\frac{\varrho e}{P},\qquad
s_{nd}=\frac{s}{R},$ $p_{nd}=\frac{p}{P},\qquad G_{nd}=\frac{\varrho
G}{P},\qquad h_{nd}=\frac{\varrho h}{P},$
$\left(c_{V}\right)_{nd}=\frac{c_{V}}{R},\qquad\left(c_{P}\right)_{nd}=\frac{c_{P}}{R},$

where $R$ is the specific gas constant of the fluid under consideration. It is
convenient to also nondimensionalize the viscosity and thermal diffusivity,

$\eta_{nd}=\frac{\eta}{\bar{\eta}},\qquad\kappa_{nd}=\frac{\kappa}{\bar{\kappa}},$

where $\bar{\eta}$ is a characteristic value of $\eta$, whereas $\bar{\kappa}$
is _not_ that of $\kappa$ – but is given by

$\bar{\kappa}=\frac{KR\varrho^{3}}{\bar{\eta}}.$

This choice of $\bar{\kappa}$ conveniently eliminates all nondimensional
parameters in the governing equations, but one should keep in mind that
$\kappa_{nd}$ can be large or small.

In terms of the new variables, Eqs. (14)–(16) take the form (the subscript nd
omitted)

$\frac{\partial\rho}{\partial t}+\frac{\partial\left(\rho w\right)}{\partial
z}=0,$ (32) $p=\eta\frac{\partial w}{\partial
z}+\rho\frac{\partial^{2}\rho}{\partial
z^{2}}-\frac{1}{2}\left(\frac{\partial\rho}{\partial z}\right)^{2}+p^{(v)},$
(33) $\frac{\partial\left(\rho e\right)}{\partial t}+\frac{\partial}{\partial
z}\left[w\left(\rho e+p-\eta\frac{\partial w}{\partial
z}\right)-\kappa\frac{\partial T}{\partial z}\right]\\\
=w\rho\frac{\partial^{3}\rho}{\partial z^{3}}.$ (34)

The nondimensional versions of the thermodynamics identities from Sec. II.1
look exactly as their dimensional counterparts, and so do the boundary
conditions from Secs. II.2.3–II.2.2, provided $\rho_{0}$ is nondimensionalized
by $\varrho$ and $T_{0}$ by $P/\varrho R$.

###  II.4 The van der Waals fluid

In what follows, general conclusions will be illustrated using the van der
Waals fluid, whose nondimensional internal energy and entropy are

$e=c_{V}T-\rho,\qquad s=c_{V}\ln T-\ln\frac{\rho}{1-\rho},$ (35)

where the heat capacity at constant volume $c_{V}$ is a given constant. Then,
Eqs. (7)–(8) yield the following expressions for the pressure and chemical
potential:

$p(\rho,T)=\frac{T\rho}{1-\rho}-\rho^{2},$ (36)
$G(\rho,T)=T\left(\ln\frac{\rho}{1-\rho}+\frac{1}{1-\rho}\right)-2\rho\\\
+c_{V}T\left(1-\ln T\right).$ (37)

To illustrate the properties of the van der Waals fluid, expressions (36)–(37)
and the Maxwell construction (20)–(21) were used to compute the saturated
densities $\rho^{(l.sat)}$ and $\rho^{(v.sat)}$. The results are shown in Fig.
2. Observe that, if $T>T_{cr}$, only one phase exists, so interfaces do not.

Figure 2: The saturated densities of liquid and vapor vs. the temperature, for
the van der Waals fluid. The nondimensional critical point in this case is
$T_{cr}=8/27$ (marked by the black dot and dotted line).

##  III Isothermal evaporation

The assumption of isothermality is used in many papers on evaporation (as
quantified by the 2 million results yielded by a Google search for
“isothermal” + “evaporation”). In terms of the present model, isothermality
corresponds to the limit $\kappa\rightarrow\infty$, in which case Eq. (34)
yields $\partial T/\partial z\rightarrow 0$. Assuming that $T$ is also
independent of $t$, one can treat the temperature in Eqs. (32)–(33) as a known
parameter. The isothermal reduction of the full DIM was first examined in Ref.
[PismenPomeau00].

###  III.1 Steady evaporation

Consider first a liquid–vapor interface in an unbounded space and assume that
it is _steadily_ receding due to evaporation. Its velocity is equal to
$-E/\rho^{(l)}$ where $E$ is the nondimensional evaporation rate and
$\rho^{(l)}$, the liquid’s nondimensional density.

Mathematically, the assumption of steadiness corresponds to the following
ansatz:

$\rho=\rho(\xi),\qquad w=w(\xi),$

where

$\xi=z+\frac{E}{\rho^{(l)}}t.$

In terms of $\xi$, the general boundary conditions (17), (19) and (28), (30)
become

$w\rightarrow 0\qquad\text{as}\qquad\xi\rightarrow-\infty,$ (38)
$\rho\rightarrow\rho^{(l)}\qquad\text{as}\qquad\xi\rightarrow-\infty,$ (39)
$\frac{\mathrm{d}w}{\mathrm{d}\xi}\rightarrow
0\qquad\text{as}\qquad\xi\rightarrow+\infty,$ (40)
$\rho\rightarrow\rho^{(v)}\qquad\text{as}\qquad\xi\rightarrow+\infty.$ (41)

When rewritten in terms of $\xi$, Eq. (32) and conditions (38)–(39) yield

$w=E\left(\frac{1}{\rho}-\frac{1}{\rho^{(l)}}\right),$ (42)

which automatically satisfies condition (40) as well. Then Eq. (33) becomes

$p(\rho,T)=-\frac{E\eta(\rho,T)}{\rho^{2}}\frac{\mathrm{d}\rho}{\mathrm{d}\xi}+\rho\frac{\mathrm{d}^{2}\rho}{\mathrm{d}\xi^{2}}-\frac{1}{2}\left(\frac{\mathrm{d}\rho}{\mathrm{d}\xi}\right)^{2}\\\
+p(\rho^{(v)},T).$ (43)

This equation and boundary conditions (39) and (41) are translationally
invariant – thus, to ensure the solution’s uniqueness, the following extra
requirement is imposed

$\rho=\frac{\rho^{(l)}+\rho^{(v)}}{2}\qquad\text{at}\qquad z=0.$ (44)

Eq. (43) and conditions (39), (41), and (44) constitute a boundary-value
problem for the function $\rho(\xi)$ and undetermined parameters $E$ and
$\rho^{(l)}$. The latter can be found straight away from constraint (31) with
$T^{(l)}=T^{(v)}=T$, which yields

$p(\rho^{(l)},T)=p(\rho^{(v)},T).$ (45)

This equation relates $\rho^{(l)}$ to the (known) density $\rho^{(v)}$ of the
vapor far above the interface. To make $\rho^{(l)}$ unique, one should require
that the liquid phase is stable,

$\left(\frac{\partial p}{\partial\rho}\right)_{\rho=\rho^{(l)}}\geq 0.$

It turns out that, if $T$ is sufficiently large and $\rho^{(v)}$ is
sufficiently small, Eq. (45) does not have any solutions. Such a situation is
illustrated in Fig. 3a: if the vapor pressure happens to be in the shaded
region, it cannot match _any_ value of the liquid pressure. Such a pattern
arises only if $T>1/4$ – otherwise the local minimum of $p(\rho)$ is negative
and, for any $\rho^{(v)}$, there exists a matching value of $\rho^{(l)}$.

Figure 3: Existence of physically meaningful solutions of Eq. (45) for the van
der Waals fluid. (a) The equation of state $p(\rho)$ for the particular case
of $T=0.27315$; the region where the vapor density does not have a match in
the liquid region is shaded. (b) The existence region (marked with a letter
“E”) of solutions of Eq. (45), in the $(T,\rho^{(v)})$ parameter plane. Curve
(1) separates the values of $\rho^{(v)}$ that have a match in the liquid
region from those that do not; curve 2 is the saturation curve (the vapor at
infinity should be undersaturated); curve 3 is $T=T_{cr}$.

Fig. 3(b) shows the region of the $(T,\rho^{(v)})$ plane where Eq. (45) admits
physically meaningful solutions. What happens if $(T,\rho^{(v)})$ are such
that no steady solution exists will be clarified in Sec. II.3.

Boundary-value problem (39), (41), (43)–(44) was solved numerically (using the
function bvp5c of MATLAB) for the van der Waals fluid (36). The viscosity was
assumed to be proportional to the density, with the proportionality
coefficient implied to be scaled out during the nondimensionalization,

$\eta=\rho.$ (46)

This is the simplest model qualitatively reflecting growth of a fluid’s
viscosity with density.

Fig. 4 shows typical profiles $\rho(\xi)$ for various values of $T$ [panels
(a)] and various values of relative humidity $\rho^{(v)}/\rho^{(v.sat)}$
[panel (b)]. The latter also illustrates that, unless $T$ is close to
$T_{cr}$, the liquid’s density $\rho^{(l)}$ is close to its saturated value
$\rho^{(l.sat)}$.

Figure 4: Solutions of boundary-value problem (39), (41), (43)–(44), for
various temperatures and relative humidities. Panel (a):
$\rho^{(v)}/\rho^{(v.sat)}=0.5$; (1) $T=0.1$, (2) $T=0.15$, (3) $T=0.2$. Panel
(b): $T=0.2$; $\rho^{(v)}/\rho^{(v.sat)}=0.1,~{}0.5,~{}0.9$; these curves are
not numbered, but can still be distinguished by comparing them to the dotted
straight line marking the density $\rho^{(v.sat)}$ of saturated vapor. The
other dotted line shows the density $\rho^{(l.sat)}$ of saturated liquid.

The dependence of $E$ on $T$ is illustrated in Fig. 5. Note that the curves
depicted could not be extended to $T=0$ due to computational difficulties (the
vapor density becomes too small). In the opposite limit, the curves are
truncated due to the nonexistence of solution of Eq. (45) as explained above.

Figure 5: Nondimensional evaporation rate $E$ vs. temperature $T$, as
described by boundary-value problem (39) and (41), (43)–(44). The curves are
marked with the corresponding values of the relative humidity,
$\rho^{(v)}/\rho^{(v.sat)}$. The dotted curves show the corresponding
asymptotic result, (54)–(56).

###  III.2 The limit of nearly-saturated vapor

If the vapor is close to saturation, evaporation must be slow, i.e.,

$\frac{\rho^{(v)}}{\rho^{(v.sat)}}\approx 1\qquad\Rightarrow\qquad E\ll 1.$

In this case, boundary-value problem (39), (41), (43)–(44) can be solved
asymptotically.

If the fluid temperature is far from its critical value, the problem involves
another small parameter, $\rho^{(v.sat)}/\rho^{(l.sat)}$, making the analysis
awkward. To remedy this, it is assumed that this parameter is order-one, then
the asymptotic results are shown numerically to work for
$\rho^{(v.sat)}/\rho^{(l.sat)}\ll 1$ as well.

Define the equilibrium solution $\rho^{(0)}(z)$ as that of Eq. (43) with
$E=0$, subject to boundary conditions (39), (41), and (44) with
$\rho^{(l)}=\rho^{(l.sat)}$ and $\rho^{(v)}=\rho^{(v.sat)}$, i.e.,

$\rho^{(0)}\frac{\mathrm{d}^{2}\rho^{(0)}}{\mathrm{d}\xi^{2}}-\frac{1}{2}\left(\frac{\mathrm{d}\rho^{(0)}}{\mathrm{d}\xi}\right)^{2}\\\
-p(\rho^{(0)},T)+p(\rho^{(v.sat)},T)=0,$ (47)
$\rho^{(0)}\rightarrow\rho^{(l.sat)}\qquad\text{as}\qquad\xi\rightarrow-\infty,$
(48)
$\rho^{(0)}\rightarrow\rho^{(v.sat)}\qquad\text{as}\qquad\xi\rightarrow+\infty,$
(49)
$\rho^{(0)}=\frac{\rho^{(l.sat)}+\rho^{(v.sat)}}{2}\qquad\text{at}\qquad\xi=0.$
(50)

It can be shown that, unless $\rho^{(l.sat)}$ and $\rho^{(v.sat)}$ satisfy the
Maxwell construction (20)–(21), the above boundary-value problem does not have
a solution (see Appendix B.1).

The solution $\rho$ of the full problem (39), (41), (43)–(44) is close, but
not equal, to the equilibrium solution $\rho^{(0)}$. The deviation of the
former from the latter is comparable to the evaporation rate $E$, so let

$\rho=\rho^{(0)}+\rho^{(1)}+\cdots,$

where $\rho^{(1)}=\mathcal{O}(E)$ and $\cdots=\mathcal{O}(E^{2})$. Assume also
that $E$ scales with the deviation of the relative humidity from unity,

$1-\frac{\rho^{(v)}}{\rho^{(v.sat)}}=\mathcal{O}(E),$

and that the liquid density is close to its saturated value by the same
margin,

$1-\frac{\rho^{(l)}}{\rho^{(l.sat)}}=\mathcal{O}(E).$

To zeroth order, Eq. (47) and boundary conditions (39), (41) are satisfied
identically. The first order yields, after straightforward algebra,

$\frac{\mathrm{d}}{\mathrm{d}\xi}\left(\frac{1}{\rho^{(0)}}\frac{\mathrm{d}\rho^{(1)}}{\mathrm{d}\xi}\right)+\frac{1}{\rho^{(0)2}}\left[\frac{\mathrm{d}^{2}\rho^{(0)}}{\mathrm{d}\xi^{2}}-\frac{\partial
p(\rho^{(0)},T)}{\partial\rho^{(0)}}\right]\rho^{(1)}=\frac{E\eta(\rho^{(0)},T)}{\rho^{(0)4}}\frac{\mathrm{d}\rho^{(0)}}{\mathrm{d}\xi}-\frac{1}{\rho^{(0)2}}\frac{\partial
p(\rho^{(v.sat)},T)}{\partial\rho^{(v.sat)}}\left(\rho^{(v)}-\rho^{(v.sat)}\right),$
(51)
$\rho^{(1)}\rightarrow\rho^{(l)}-\rho^{(l.sat)}\qquad\text{as}\qquad\xi\rightarrow-\infty,$
(52)
$\rho^{(1)}\rightarrow\rho^{(v)}-\rho^{(v.sat)}\qquad\text{as}\qquad\xi\rightarrow+\infty.$
(53)

The evaporation rate $E$ can be found without finding $\rho^{(1)}(\xi)$ (and
even considering $\rho^{(2)}(\xi)$ and the higher corrections). To do so,
observe that the operator on the left-hand side of Eq. (51) is self-adjoint,
and the homogeneous version of problem (51)–(53) is satisfied by
$\rho^{(1)}=\mathrm{d}\rho^{(0)}/\mathrm{d}\xi$ – hence, the full
(nonhomogeneous) version should be orthogonal to
$\mathrm{d}\rho^{(0)}/\mathrm{d}\xi$.

To derive the orthogonality condition, multiply (51) by
$\mathrm{d}\rho^{(0)}/\mathrm{d}\xi$ and integrate from $\xi=-\infty$ to
$\xi=+\infty$. Using boundary condition (48)–(49) for $\rho^{(0)}$ and
(52)–(53) for $\rho^{(1)}$, one obtains the desired expression for the
evaporation rate:

$E=\frac{H}{A},$ (54)

where

$\displaystyle H$
$\displaystyle=\left(1-\frac{\rho^{(v.sat)}}{\rho^{(l.sat)}}\right)\frac{\partial
p(\rho^{(v.sat)},T)}{\partial\rho^{(v.sat)}}\frac{\rho^{(v.sat)}-\rho^{(v)}}{\rho^{(v.sat)}},$
(55) $\displaystyle A$
$\displaystyle=\int_{-\infty}^{\infty}\frac{\eta(\rho^{(0)},T)}{\rho^{(0)4}}\left(\frac{\mathrm{d}\rho^{(0)}}{\mathrm{d}\xi}\right)^{2}\mathrm{d}\xi.$
(56)

Observe that the last factor in the expression for $H$ is the difference
between the relative humidity and unity.

Asymptotic result (54)–(56) is compared to the numeric solution of the exact
problem in Fig. 5. Interestingly, the agreement between the two solutions is
reasonably good even when the relative humidity is far from unity: the
relative error of the most part of the asymptotic curve for the 50% humidity
is less than $0.15$ and that for the 10% curve is less than $0.22$. These
errors are only exceeded near the terminal point, where the asymptotic
solution noticeably under-predicts the exact one.

The coefficient $A$ given by expression (56) always arises when the DIM is
used to examine evaporation [Benilov22a], condensation [Benilov22b], or
another setting where these processes play a role [Benilov20d]. It involves
the leading-order solution $\rho^{(0)}(z)$ which needs to be computed before
the integral in (56) can be evaluated. To avoid this extra computation, one
can express $A$ as a closed-form integral (see Appendix B.2)

$A=2^{1/2}\int_{\rho^{(v.sat)}}^{\rho^{(l.sat)}}\frac{\eta(\rho,T)}{\rho^{4}}\sqrt{\rho\left[G(\rho,T)-G(\rho^{(v.sat)},T)\right]-p(\rho,T)+p(\rho^{(v.sat)},T)}\mathrm{d}\rho.$
(57)

If the nondimensional temperature is low (which is the case for many common
fluids at “normal conditions” [Benilov20b]), expressions (57) can be
calculated asymptotically,

$A\approx
0.14219\frac{\eta(0,T)\,T^{1/2}}{\left(\rho^{(v.sat)}\right)^{5/2}}\qquad\text{if}\qquad
T\ll 1,$ (58)

where $\eta(0,T)$ is the small-density limit of the viscosity $\eta(\rho,T)$.
Since $\rho^{(v.sat)}\rightarrow 0$ as $T\rightarrow 0$, one can deduce from
the above expression and formula (54) that the evaporation rate vanishes as
$T\rightarrow 0$ (as it should do physically).

###  III.3 Unsteady evaporation

It remains to find out how the liquid evaporates if no steady solution exists
. In terms of Fig. 3b, this occurs when the pair $(T,\rho^{(v)})$ is outside
the existence region.

The steady and unsteady scenarios of evaporation were simulated using the set
of evolution equations (32)–(33) using the method of lines [Schiesser78].
Typical results are illustrated in Fig. 6.

Figure 6: Two scenarios of isothermal evaporation. In both cases, $T=0.26$ and
the near-substrate density $\rho_{0}$ equals, for simplicity, the value
prescribed by the initial condition. (a) $\rho^{(v)}/\rho^{(v.sat)}=0.4$ (a
steady solution exists and is used as the initial condition); (b)
$\rho^{(v)}/\rho^{(v.sat)}=0.3$ (steady solution does not exist, the initial
condition is given by (58)–(59)]. In both cases, the time $t$ of the snapshot
is related to the curve number $n$ by $t=22\,n$.

The parameters of Fig. 6(a) are such that a steady solution exists. Even
though it was obtained for an unbounded space, the presence of the substrate
does not change it until the interface is very close to the substrate (this
stage of the evolution is not shown in the figure). The run depicted in Fig.
6(a) originates from an initial condition obtained by solving the steady
problem (39), (41), (43)–(44) and shifting the solution to the right by a
distance of $\Delta z=65$ from the substrate.

For the parameters of Fig. 6(b), no steady solution exists. The following
initial condition was used:

$\rho=\rho^{(0)}(z-\Delta
z)\left[\left(\frac{\rho^{(v)}}{\rho^{(v.sat)}}+1\right)\right.\\\
+\left.\left(\frac{\rho^{(v)}}{\rho^{(v.sat)}}-1\right)\tanh\frac{z-z_{1}}{W}\right],$
(59)

where $\rho^{(0)}(z)$ is the equilibrium solution described by (47)–(50),
$\Delta z$ is the initial position of the interface, $z_{1}$ is the position
of a transitional region where the vapor density changes from $\rho^{(v.sat)}$
[as ‘prescribed’ by $\rho^{(0)}(z)$] to $\rho^{(v)}$ corresponding to the
chosen relative humidity, and $W$ is this region’s width . In Fig. 6(b), the
following parameter values were used:

$\Delta z=65,\qquad z_{1}=67,\qquad W=2.$ (60)

The difference between the two scenarios is evident: if no steady solution
exists, the whole layer ‘empties out’ – quickly and all at once. This occurs
because the pressure in the liquid cannot equilibrate with that in the vapor
(recall that condition (45) does not hold) – and the resulting pressure
gradient generates a strong evaporative flow. It can be conjectured that, in a
three-dimensional model – where bubbles may arise – the liquid in this regime
_boils_. If this conjecture is true, the impossibility of matching the
pressure in the vapor and liquid phases at a certain temperature provides a
criterion for a ‘near-vacuum boiling’ (which includes the ‘vacuum boiling’
proper as a limiting case).

Observe the near-substrate boundary layer in Fig. 6(b). It develops due to
boundary condition (27) forcing the near-substrate density to remain fixed.
Computations with various values of $\rho_{0}$ have shown that this boundary
layer has no impact on the global dynamics.

Note that evaporation of a pure fluid has been previously examined in Ref.
[BarbanteFrezzottiGibelli14] – and not only via the DIM, but also via
simulations molecular dynamics. The parameter range explored in this paper
happened to be inside the existence region of steady solutions, so the
unsteady regime was not observed.

##  IV Nonisothermal evaporation

If evaporation is steady, the density field near the interface is microscopic
– i.e., its dimensional spatial scale is comparable to $l$ introduced in Sec.
II.3. For the temperature field, however, no mechanism exists for keeping it
short-scale. As a result, evaporative cooling rapidly spreads out and
eventually become macroscopic.

Thus, the problem splits into two subproblems:

(A) an analysis of the _macroscopic_ temperature field, which yields the heat
fluxes coming into the interface from the vapor and liquid sides,

(B) an analysis of the _microscopic_ interfacial dynamics, which inter-relates
the heat fluxes and evaporation rate.

Subproblems (A) and (B) are examined in Secs. IV.1 and IV.2, respectively. The
former is solved three times: for an unbounded space (the simplest case), for
a semi-infinite space bounded by a fixed-temperature substrate (the case with
most applications), and for an insulated substrate.

Admittedly, the analyses presented below involve several assumptions, so the
results obtained will be verified via numerical simulations of the exact
governing equations in Sec. IV.3.

###  IV.1 Macroscopic solution

The evaporative cooling is unlikely to exceed, say, 10–20 degrees – otherwise,
as argued in the Introduction, the evaporation would effectively stop. Within
the normal-conditions range, 20 degrees is a small fraction of the absolute
temperature, so that $T(z,t)$ can be decomposed into a constant part $T_{0}$
and a small variation $\tilde{T}(z,t)$,

$T=T_{0}+\tilde{T}.$ (61)

For a fixed-$T$ substrate, $T_{0}$ is its temperature.

Consider liquid and vapor far from the interface, were the density is close to
either $\rho^{(l)}$ or $\rho^{(v)}$. The spatial scale of the temperature
field in these regions is macroscopic, which nondimensionally means “much
greater than unity”. Under such assumptions, the governing equations (32)–(34)
can be reduced to a heat conduction equation (see Appendix 101).

Let the position of the interface be $z=z_{i}(t)$. Its velocity can be
expressed in terms of the evaporation rate and the liquid’s density,

$\frac{\mathrm{d}z_{i}}{\mathrm{d}t}=-\frac{E}{\rho^{(l)}}.$

Changing to the co-moving reference frame
$\left(z,t\right)\rightarrow\left(\xi,t\right)$, where

$\xi=z-z_{i}=z+\int\frac{E(t)}{\rho^{(l)}}\mathrm{d}t,$ (62)

one can describe the temperature evolution by

$\displaystyle\rho^{(l)}c_{P}^{(l)}\left(\frac{\partial\tilde{T}}{\partial
t}+\frac{E}{\rho^{(l)}}\frac{\partial\tilde{T}}{\partial z}\right)$
$\displaystyle=\kappa^{(l)}\frac{\partial^{2}\tilde{T}}{\partial z^{2}}\hskip
11.38092pt\,\text{if}\hskip 11.38092ptz<0,$ (63)
$\displaystyle\rho^{(v)}c_{P}^{(v)}\left(\frac{\partial\tilde{T}}{\partial
t}+\frac{E}{\rho^{(l)}}\frac{\partial\tilde{T}}{\partial z}\right)$
$\displaystyle=\kappa^{(v)}\frac{\partial^{2}\tilde{T}}{\partial z^{2}}\hskip
11.38092pt\text{if}\hskip 11.38092ptz>0,$ (64)

where

$c_{P}^{(l)}=c_{P}(\rho^{(l)},T_{0}),\qquad
c_{P}^{(v)}=c_{P}(\rho^{(v)},T_{0}),$
$\kappa^{(l)}=\kappa(\rho^{(l)},T_{0}),\qquad\kappa^{(v)}=\kappa(\rho^{(v)},T_{0}),$

Eqs. (63)–(64) are to be solved with the following conditions:

$\left(\tilde{T}\right)_{z=+0}=\left(\tilde{T}\right)_{z=-0},$ (65)
$\kappa^{(v)}\left(\frac{\partial\tilde{T}}{\partial
z}\right)_{z=+0}-\kappa^{(l)}\left(\frac{\partial\tilde{T}}{\partial
z}\right)_{z=-0}=E\,\Delta h,$ (66)

where the right-hand side of (66) describes consumption of heat in the
interface due to evaporative cooling. The simplest initial condition will be
assumed,

$\tilde{T}=0\qquad\text{at}\qquad t=0,$ (67)

i.e., initially, the temperature field is uniform.

The problem formulated contains a hidden small parameter. To identify it,
denote the vertical spatial scale of the solution by $D$ (in case of a finite
layer, $D$ coincides with its depth) and observe that

Eq. (63) $\displaystyle\text{:}\qquad\frac{E\text{-involving
term}}{\text{right-hand side}}\sim ED\frac{c_{P}^{(l)}}{\kappa^{(l)}},$ (68)
Eq. (64) $\displaystyle\text{:}\qquad\frac{E\text{-involving
term}}{\text{right-hand side}}\sim
ED\frac{c_{P}^{(v)}}{\kappa^{(v)}}\frac{\rho^{(v)}}{\rho^{(l)}}.$ (69)

Thus, if

$D\ll\frac{\kappa^{(l)}}{Ec_{P}^{(l)}},$ (70)

the $E$-involving term in Eq. (63) can be neglected.

To understand how restrictive condition (70) is, one can estimate its right-
hand side for water at, say, $25^{\circ}\mathrm{C}$. Using Ref.
[LindstromMallard97] to obtain values for $\kappa^{(l)}$ and $c_{P}^{(l)}$ and
assuming estimate (4) for $E$, one obtains

$D\ll 5.7\,\mathrm{m}.$

Many physically important examples comply with this restriction, and those
that do not are discussed in Sec. IV.3.

Observe also the small factor $\rho^{(v)}/\rho^{(l)}$ in estimate (69) –
hence, its right-hand is likely to be smaller than that of (68). Estimating
the ratio of the two right-hand sides for water and air at
$25^{\circ}\mathrm{C}$, one obtains

$\frac{\left(c_{P}^{(v)}/\kappa^{(v)}\right)\left(\rho^{(v)}/\rho^{(l)}\right)}{\left(c_{P}^{(l)}/\kappa^{(l)}\right)}\approx
0.0066.$

Thus, if the $E$-involving term in Eq. (63) is negligible, its counterpart in
Eq. (64) is too.

Omitting these terms in Eqs. (63)–(64), one obtains

$\displaystyle\rho^{(l)}c_{P}^{(l)}\frac{\partial\tilde{T}}{\partial t}$
$\displaystyle=\kappa^{(l)}\frac{\partial^{2}\tilde{T}}{\partial
z^{2}}\qquad\,\text{if}\qquad z<0,$ (71)
$\displaystyle\rho^{(v)}c_{P}^{(v)}\frac{\partial\tilde{T}}{\partial t}$
$\displaystyle=\kappa^{(v)}\frac{\partial^{2}\tilde{T}}{\partial
z^{2}}\qquad\text{if}\qquad z>0.$ (72)

These equations do not explicitly involve $t$ – hence, can be solved via the
Laplace transformation (the presence of $E(t)$ in the matching condition (66)
does not pose a problem, as this is not a coefficient but a right-hand side).

####  IV.1.1 A flat interface in an unbounded space

This case corresponds to the following boundary conditions:

$\tilde{T}\rightarrow 0\qquad\text{as}\qquad z\rightarrow\pm\infty.$ (73)

Applying the Laplace transformation to Eqs. (71)–(72), then recalling initial
condition (67) and boundary conditions (65)–(66) and (73), one obtains

$\displaystyle\hat{T}$ $\displaystyle=-\frac{s^{-1/2}\hat{E}\,\Delta
h}{\sqrt{\rho^{(l)}c_{P}^{(l)}\kappa^{(l)}}+\sqrt{\rho^{(v)}c_{P}^{(v)}\kappa^{(v)}}}\exp\left(\sqrt{\frac{\rho^{(l)}c_{P}^{(l)}s}{\kappa^{(l)}}}z\right)\hskip
32.72049pt\text{if}\qquad z<0,$ $\displaystyle\hat{T}$
$\displaystyle=-\frac{s^{-1/2}\hat{E}\,\Delta
h}{\sqrt{\rho^{(l)}c_{P}^{(l)}\kappa^{(l)}}+\sqrt{\rho^{(v)}c_{P}^{(v)}\kappa^{(v)}}}\exp\left(-\sqrt{\frac{\rho^{(v)}c_{P}^{(v)}s}{\kappa^{(v)}}}z\right)\qquad\text{if}\qquad
z>0,$

where the variables with hats represent the transforms, e.g.,

$\hat{T}(s,t)=\int_{0}^{\infty}\tilde{T}(z,t)\operatorname{e}^{-st}\mathrm{d}t.$

The inverse transform of $\hat{T}(s,t)$ will be matched to the microscopic
solution obtained later. Three characteristics of the former will be needed
for the matching: the temperature at the interface and the heat fluxes toward
it. Calculating the inverse transforms of these characteristics only [which is
easier than calculating the whole $T(z,t)$], one obtains

$\left(\tilde{T}\right)_{z=0}=-\frac{\Delta
h}{\sqrt{\rho^{(l)}c_{P}^{(l)}\kappa^{(l)}}+\sqrt{\rho^{(v)}c_{P}^{(v)}\kappa^{(v)}}}\int_{0}^{t}\frac{E(t^{\prime})}{\sqrt{\pi\left(t-t^{\prime}\right)}}\mathrm{d}t^{\prime},$
(74) $\kappa^{(l)}\left(\frac{\partial\tilde{T}}{\partial
z}\right)_{z=-0}=-\frac{\sqrt{\rho^{(l)}c_{P}^{(l)}\kappa^{(l)}}E\,\Delta
h}{\sqrt{\rho^{(l)}c_{P}^{(l)}\kappa^{(l)}}+\sqrt{\rho^{(v)}c_{P}^{(v)}\kappa^{(v)}}},\qquad\kappa^{(v)}\left(\frac{\partial\tilde{T}}{\partial
z}\right)_{z=+0}=\frac{\sqrt{\rho^{(v)}c_{P}^{(v)}\kappa^{(v)}}E\,\Delta
h}{\sqrt{\rho^{(l)}c_{P}^{(l)}\kappa^{(l)}}+\sqrt{\rho^{(v)}c_{P}^{(v)}\kappa^{(v)}}}.$
(75)

Note that (75) [but not (74)] can be deduced from basic symmetries of the
heat-conduction equations (71)–(72) – which probably explains their relatively
simple form.

####  IV.1.2 Fixed-temperature substrate

In the coordinate system co-moving with the interface, a fixed-$T$ substrate
corresponds to the following boundary condition,

$\tilde{T}=0\qquad\text{at}\qquad z=-D,$ (76)

where $D(t)$ is the depth of the liquid layer. Since this condition is set at
a moving boundary, Eqs. (71)–(72) can no longer be solved exactly, but one can
still solve them asymptotically under the same assumption (70) which makes the
$E$-involving term negligible. It makes evaporation slower than the
temperature evolution – hence, one can ‘freeze’ $D$ – i.e., assume that it
depends on a slow-time variable different from $t$.

Now, apply the Laplace transformation to Eqs. (71)–(72) and boundary
conditions (65)–(66), (76) with ‘frozen’ $D$. After straightforward algebra,
one can show that the temperature field tends to

$\tilde{T}\rightarrow\left\\{\begin{tabular}[c]{ll}$-\dfrac{E\,\Delta
h}{\kappa^{(l)}}\left(z+D\right)\qquad$&$\text{if}\qquad-D<z<0,$\\\
$-\dfrac{E\,\Delta h}{\kappa^{(l)}}D$&$\text{if}\qquad
z>0$\end{tabular}\right.\qquad\text{as}\qquad t\rightarrow\infty.$ (77)

Solution (77) has a clear physical meaning: it describes a quasi-steady heat
flux from the substrate toward the interface, feeding evaporation. The vapor
above the interface has uniform temperature because the temperature variations
have spread out to $z\rightarrow+\infty$.

Note that the linear dependence of $T$ on $z$ in liquid has indeed been
observed experimentally [WardStanga01].

####  IV.1.3 Thermally insulated substrate

The boundary condition at the substrate in this case is

$\frac{\partial\tilde{T}}{\partial z}=0\qquad\text{at}\qquad z=-D.$

As before, the initial-boundary-value problem for $\tilde{T}(z,t)$ (with
$D(t)$ ‘frozen’), can be solved via the Laplace transformation – but the
resulting solution is much more cumbersome than those in the two previous
cases. One can still see that, in the liquid, the temperature field becomes
spatially uniform, with a small parabolic correction:

$\tilde{T}\approx\tilde{T}_{0}(t)\\\
+\frac{\mathrm{d}\tilde{T}_{0}}{\mathrm{d}t}\frac{\rho^{(l)}c_{P}^{(l)}}{\kappa^{(l)}}\frac{\left(D+z\right)^{2}}{2}\qquad\text{if}\qquad-D<z<0.$
(78)

Since the case of insulated substrate is cumbersome mathematically and trivial
physically, it will not be examined analytically in further detail (but will
be simulated numerically in Sec. IV.3).

###  IV.2 Microscopic solution

Strictly speaking, nonisothermal evaporation is not steady, but can be
regarded as _quasi_ steady – i.e., adjusted to the current temperature of the
interface and incoming heat fluxes, and slowly changing with them.

In this section, $E$ is calculated asymptotically for quasi-steady
nonisothermal evaporation under an additional assumption that the relative
humidity is close to 100%. Only the case of fixed-$T$ substrate will be
examined; those of insulated substrate and unbounded space will be briefly
discussed in the next subsection.

The assumption of quasi-steadiness amounts to rewriting the governing
equations (32)–(34) in the co-moving reference frame $\left(\xi,t\right)$
[$\xi$ is defined by (62)] and omitting the time derivatives, which yields

$\frac{E}{\rho^{(l)}}\frac{\partial\rho}{\partial\xi}+\frac{\partial\left(\rho
w\right)}{\partial\xi}=0,$ (79) $p=\eta\frac{\partial
w}{\partial\xi}+\rho\frac{\partial^{2}\rho}{\partial\xi^{2}}-\frac{1}{2}\left(\frac{\partial\rho}{\partial\xi}\right)^{2}+p^{(v)},$
(80) $\frac{E}{\rho^{(l)}}\frac{\partial\left(\rho
e\right)}{\partial\xi}+\frac{\partial}{\partial\xi}\left[w\left(\rho
e+p-\eta\frac{\partial w}{\partial\xi}\right)-\kappa\frac{\partial
T}{\partial\xi}\right]\\\ =w\rho\frac{\partial^{3}\rho}{\partial\xi^{3}}.$
(81)

These equations are to be solved with the ‘old’ boundary conditions for $w$,
(38) and (40), whereas those for $\rho$ and $T$ will be obtained later via
matching with the macroscopic solution.

As in the isothermal case, one can reduce Eq. (79) to expression (42) for $w$,
and substitute it into Eqs. (80)–(81) which become

$p(\rho,T)=-\frac{E\eta(\rho,T)}{\rho^{2}}\frac{\partial\rho}{\partial\xi}+\rho\frac{\partial^{2}\rho}{\partial\xi^{2}}-\frac{1}{2}\left(\frac{\partial\rho}{\partial\xi}\right)^{2}+p(\rho^{(v)},T_{0}),$
(82) $\kappa(\rho,T)\frac{\partial
T}{\partial\xi}=E\left[e(\rho,T)+\left(1-\frac{\rho}{\rho^{(l)}}\right)\frac{p(\rho,T)}{\rho}-\left(1-\frac{\rho}{\rho^{(l)}}\right)\frac{\partial^{2}\rho}{\partial\xi^{2}}-\frac{1}{2\rho^{(l)}}\left(\frac{\partial\rho}{\partial\xi}\right)^{2}-Q\right]+E^{2}\left(\frac{1}{\rho}-\frac{1}{\rho^{(l)}}\right)\frac{\eta}{\rho^{2}}\frac{\partial\rho}{\partial\xi},$
(83)

where Eq. (83) was integrated and $Q$ is the integration constant (to be fixed
later). As before, assume that the solution is close to equilibrium: the
leading-order temperature is uniform and equals $T_{0}$ [the same constant as
in the macroscopic representation (61)], and the leading-order density is
described by the solution $\rho^{(0)}$ of the equilibrium boundary-value
problem (47)–(50). These assumptions amount to

$\rho=\rho^{(0)}+\rho^{(1)}+\mathcal{\cdots},\qquad T=T_{0}+T^{(1)}+\cdots,$

where $\rho^{(1)}$ and $T^{(1)}$ are $\mathcal{O}(E)$. Upon substitution of
these expansions into Eq. (82), the zeroth order cancels out and the first-
order yields

$\frac{\partial}{\partial\xi}\left(\frac{1}{\rho^{(0)}}\frac{\partial^{2}\rho^{(1)}}{\partial\xi^{2}}\right)+\frac{1}{\rho^{(0)2}}\left[\frac{\partial^{2}\rho^{(0)}}{\partial\xi^{2}}-\frac{\partial
p(\rho^{(0)},T_{0})}{\partial\rho^{(0)}}\right]\rho^{(1)}\\\
=\frac{E\eta(\rho^{(0)},T_{0})}{\rho^{4}}\frac{\partial\rho^{(0)}}{\partial\xi}-\frac{1}{\rho^{(0)2}}\frac{\partial
p(\rho^{(v.sat)},T_{0})}{\partial\rho^{(v.sat)}}\left(\rho^{(v)}-\rho^{(v.sat)}\right)+\frac{1}{\rho^{(0)2}}\frac{\partial
p(\rho^{(0)},T_{0})}{\partial T_{0}}T^{(1)}.$ (84)

By comparison with its isothermal counterpart (51), Eq. (84) includes an extra
term (the last term on its right-hand side). Expanding, in turn, the
temperature equation (83), one obtains

$\kappa(\rho^{(0)},T_{0})\frac{\partial
T^{(1)}}{\partial\xi}=E\left[e(\rho^{(0)},T_{0})+\left(\frac{1}{\rho^{(0)}}-\frac{1}{\rho^{(l.sat)}}\right)p(\rho^{(0)},T_{0})-\left(1-\frac{\rho^{(0)}}{\rho^{(l.sat)}}\right)\frac{\partial^{2}\rho^{(0)}}{\partial\xi^{2}}-\frac{1}{2\rho^{(l.sat)}}\left(\frac{\partial\rho^{(0)}}{\partial\xi}\right)^{2}-Q\right].$

The derivatives of $\rho^{(0)}$ in this equation can be eliminated via
equalities (97) and (99), and after straightforward algebra involving the use
of definition (8) of the chemical potential, one obtains

$\kappa(\rho^{(0)},T_{0})\frac{\partial
T^{(1)}}{\partial\xi}=ET_{0}\left[s(\rho^{(0)},T_{0})-s(\rho^{(v.sat)},T_{0})\right]+e(\rho^{(v.sat)},T_{0})-Q.$
(85)

The solution of this equation and the outer (macroscopic) solution (77) match
only if the former satisfy the following boundary condition

$\displaystyle\frac{\partial T^{(1)}}{\partial\xi}$
$\displaystyle\rightarrow-\frac{E\,\Delta
h(T_{0})}{\kappa(\rho^{(l.sat)},T_{0})}\hskip
29.87547pt\text{as}\qquad\xi\rightarrow-\infty,$ (86) $\displaystyle T^{(1)}$
$\displaystyle\rightarrow-\dfrac{E\,\Delta
h(T_{0})}{\kappa(\rho^{(l.sat)},T_{0})}D\qquad\text{as}\qquad\xi\rightarrow+\infty.$
(87)

Condition (87) implies that the constant $Q$ in Eq. (85) is

$Q=e(\rho^{(v.sat)},T_{0}),$ (88)

and condition (86) yields the same result, but subject to

$\Delta
h(T_{0})=T_{0}\left[s(\rho^{(v.sat)},T_{0})-s(\rho^{(l.sat)},T_{0})\right],$
(89)

which is actually an identity (as shown in Appendix D).

Keeping in mind that asymptotics (86) predicts linear behavior of
$T^{(1)}(\xi)$ as $\xi\rightarrow-\infty$, one can show that Eq. (84) admits a
similar linear asymptotics for $\rho^{(1)}$, i.e.,

$\rho^{(1)}=\mathcal{O}(\xi)\qquad\text{as}\qquad\xi\rightarrow-\infty.$ (90)

In vapor, one should, as before, assume

$\rho^{(1)}\rightarrow\rho^{(v)}-\rho^{(v.sat)}\qquad\text{as}\qquad\xi\rightarrow+\infty.$
(91)

The evaporation rate $E$ can be found using the same procedure as that in the
isothermal case. Multiplying Eq. (84) by $\mathrm{d}\rho^{(0)}/\mathrm{d}\xi$,
one eliminates $\rho^{(1)}$ by integrating by parts, and recalling boundary
conditions (90)–(91) for $\rho^{(1)}$ and boundary-value problem (48)–(49) for
$\rho^{(0)}$. After straightforward algebra involving the use of thermodynamic
identity (10), one obtains

$E=\frac{H}{A+A^{\prime}},$ (92)

where $H$ and $A$ are given by the isothermal expressions (55)–(56),
respectively, and

$A^{\prime}=-\frac{1}{E}\int_{-\infty}^{\infty}T^{(1)}\frac{\mathrm{d}s(\rho^{(0)},T_{0})}{\mathrm{d}\xi}\mathrm{d}\xi$

can be shown to be independent of $E$.

To do so, recall that the function $T^{(1)}$ in $A^{\prime}$ is the solution
of boundary-value problem (85)–(88) which needs to be solved before
$A^{\prime}$ can be calculated. This can be by-passed, however, by rearranging
$A^{\prime}$ in the form

$A^{\prime}=-\frac{1}{E}\int_{-\infty}^{\infty}\left[T^{(1)}+\dfrac{E\,\Delta
h}{\kappa(\rho^{(l.sat)},T_{0})}D\right]\frac{\mathrm{d}\left[s(\rho^{(0)},T_{0})-s(\rho^{(l.sat)},T_{0})\right]}{\mathrm{d}\xi}\mathrm{d}\xi\\\
+\dfrac{\Delta
h(T_{0})}{\kappa(\rho^{(l.sat)},T_{0})}D\left[s(\rho^{(v.sat)},T_{0})-s(\rho^{(l.sat)},T_{0})\right],$

then integrating the first term by parts, recalling boundary conditions
(86)–(87), and rearranging the second term using identity (89). As a result,
the expression for $A^{\prime}$ involves only $\partial T^{(1)}/\partial\xi$,
which can be eliminated using Eqs. (85) and (88). Eventually, one obtains

$A^{\prime}=A_{1}^{\prime}+A_{2}^{\prime}D,$ (93)

where

$\displaystyle A_{1}^{\prime}$
$\displaystyle=T_{0}\int_{-\infty}^{\infty}\frac{\left[s(\rho^{(0)},T_{0})-s(\rho^{(l.sat)},T_{0})\right]\left[s(\rho^{(0)},T_{0})-s(\rho^{(v.sat)},T_{0})\right]}{\kappa(\rho^{(0)},T_{0})}\mathrm{d}\xi,$
$\displaystyle A_{2}^{\prime}$
$\displaystyle=\dfrac{T_{0}\left[s(\rho^{(v.sat)},T_{0})-s(\rho^{(l.sat)},T_{0})\right]^{2}}{\kappa(\rho^{(l.sat)},T_{0})}.$

Expressions (92)–(93) constitute the most general asymptotic result of the
present work. It can be used to find the dependence of the evaporation rate
$E$ and the liquid’s depth $D$ on the time variable $t$; to do so, one needs
to complement (92)–(93) with the ordinary differential equation

$\frac{\mathrm{d}D}{\mathrm{d}t}=-\frac{E}{\rho^{(l.sat)}}.$

With $E$ determined by expressions (92), the above equation can be readily
solved.

###  IV.3 Numerical simulations

The above predictions regarding the substrate’s effect on evaporation have
been tested via numerical integration of the exact governing equations
(32)–(34), for various initial conditions and parameter values. Note that, in
the numerics, the energy equation (34) was conveniently replaced with the
(mathematically equivalent) temperature equation (100).

Figs. 7–8 illustrate a typical evolution. It was computed for the viscosity
$\eta$ given by (46) and

$c_{V}=3,\qquad\kappa=\kappa_{0}\rho,$

where the former value corresponds to diluted water vapor, and the latter
reflects the general tendency of thermal conductivity to increase with
density. Comparing the nondimensionalization of this paper with that of Ref.
[Benilov23a], one can deduce that $\kappa_{0}$ is the reciprocal of the
“isothermality parameter” $\beta$ introduced in the latter: if, in a certain
region of the flow, $\kappa_{0}$ is large ($\beta$, small), the temperature
field in this region is almost uniform. According to the estimate of Ref.
[Benilov23a] for water under normal conditions, $\beta\approx 0.06$, which
approximately corresponds to

$\kappa_{0}=17.$

The simplest initial condition for the temperature was used,
$T=\operatorname{const}$, and the initial density field was represented by the
steady isothermal solution – i.e., that of (39), (41), (43)–(44) – for
$\rho^{(v)}/\rho^{(v.sat)}=0.9$, shifted to the right by a distance of $\Delta
z=20$ from the substrate. The boundary conditions for the vapor density and
temperature were moved from $z=+\infty$ to $z=400$, and the computation was
stopped well before the perturbations of $T$ and $\rho$ reached anywhere near
this point.

Figure 7: Numerical solutions for a fixed-temperature substrate (left-hand
panels) and insulated substrate (right-hand panels): (a) density field, (b)
temperature field, (c) heat flux. The initial conditions and parameters are
described in Sec. IV.3. The curves represent the snapshots taken at the
following times: (0) $t=0$, (1) $t=20$, (2) $t=200$, (3) $t=500$, (4)
$t=1400$. The density distribution evolves very slowly, to the extent that
curves (0)–(4) in the upper panels are indistinguishable. Note that this
figure shows only a fraction of the $z$ range of the underlying computation
which was carried for $z\in\left(0,400\right)$. Figure 8: The evaporation rate
vs. time, for the computations represented in Fig. 7: (i) insulated substrate,
(f) fixed-temperature substrate [the dotted curve shows the asymptotic
solution (92 )–(93)].

The following features of the evolution can be observed:

  * •

Fig. 7 shows that the evolution of the density field is much slower than the
temperature evolution (as assumed in the asymptotic analysis). This conclusion
applies to both kinds of substrates.

  * •

The two lower left-hand panels of Fig. 7 show that, for the fixed-$T$
substrate, the temperature field in the liquid settles into a quasi-steady
pattern with linear dependence of $T$ on $z$ and uniform heat flux [as
predicted by asymptotics (77)]. Fig. 8 shows that the evaporation rate is
changing very slowly in this case and is close to its asymptotic value
predicted by formulae (92)–(93).

  * •

The two right-hand lower panels of Fig. 7 show that, for an insulated
substrate, the liquid’s temperature becomes nearly-uniform [as predicted by
asymptotics (78)] and is rapidly decreasing, while $E$ tends to zero (see Fig.
8). Since $T$ cannot drop below the dew point (at which the initial density of
vapor coincides with the saturated density), it can only approach it from
above222Calculating the dew point for the parameters of Figs. 7–8, one obtains
$\left(T-T_{0}\right)\times 10^{3}\approx-4.294$. The convergence to this
value in the actual figures is extremely slow, making it difficult to extend
the simulation to the stage where it is almost reached..

  * •

This effectively means that, in the case of insulated substrate, the
evaporation rate depends on $t$ – hence, is not fully determined by a set of
external (measurable) parameters. It can only be measured directly.

  * •

For fixed-$T$ substrates, in turn, the evaporation rate does not tend to zero
and _is_ determined by external parameters (as evidenced by the good agreement
between the dotted and solid curves, both marked with (f), in Fig. 8).

###  IV.4 Can the problem be solved if condition (70) does not hold?

If the depth $D$ of the evaporating layer is large enough, the timescale of
the long-term temperature evolution is comparable to that of evaporation. This
invalidates the asymptotic approach used in Sec. (IV.1) for solving the
macroscopic problem. Another shortcoming of the asymptotic approach is that it
works only if the relative humidity is close to unity.

However, even though there is no timescale separation in this case, one can
still exploit the _spatial_ scale separation – i.e., the difference between
the microscopic scale of the interface and the macroscopic scale of the liquid
layer. This implies solving the macroscopic equations (63)–(64) numerically
and, at each time step, feeding the computed temperature of, and fluxes at,
the interface into the microscopic problem – then using it to compute $E$ from
the microscopic problem and feeding $E$ back into the macroscopic equations.

The approach outlined above will not be discussed in further detail. It is
worth implementing only after the present model is extended to evaporation
into air (as opposed to the liquid’s own vapor).

##  V Can other factors contribute to the discrepancies among experimental
results?

There are two factors, not included in the present model, which may
potentially contribute to the discord among the empiric formulae illustrated
in Fig. 1. These factors are convection in the liquid and diffusion of the
vapor in air.

  * •

In a layer on a fixed-$T$ substrate, evaporation can cause a large temperature
difference between the substrate and interface – which, in turn, can trigger
off convection. The size and number of the convection cells depend on the
depth $D$ of the tank with liquid, but also on its width $L$ – hence, two
experiments differing by the value of $L$ but identical otherwise may yield
different results. It is not a priory clear how large this difference is, but
it is telling that some of the empiric formulae listed in Refs. [PoosVarju20,
VarjuPoos22] do involve $L$.

  * •

It is well known that diffusion of vapor, or any other substance, in an
_infinite_ semispace does not have a steady regime: the one-dimensional
diffusion equation does not simply have such solutions. Physically, the vapor
accumulates near the interface in this case, making the evaporation rate tend
to zero with time – while the layer of nearly-saturated vapor slowly expands
toward infinity.
If, however, the vapor is collected at a _finite_ height above the interface,
a steady regime does exist, such that the vapor flux is spatially uniform and
vapor concentration is linear. This effectively means that the measured
evaporation rate should depend on the evaporation chamber’s height. Such a
dependence is indeed acknowledged in some papers (e.g., Ref.
[SparrowKratzSchuerger83]), but not mentioned in the others (including the
results presented in Fig. 1).

Both of the above factors can be incorporated into the diffuse-interface
model: to study convection, a three-dimensional version of the DIM should be
used, whereas diffusion of vapor in air can be examined via a multicomponent
version [Benilov23a]. The work on the latter problem is in progress.

Among other factors affecting evaporation, one might mention radiative
exchange of heat between liquid and air. It is often small, but can still be
important for the case of insulated vessels, where the heat fluxes from the
boundaries are zero.

##  VI Summary and concluding remarks

The following features of nonisothermal evaporation of liquid layers on a
substrate are reported in this paper:

  1. (a)

Qualitative estimates (1) and (5) show that nonisothermality affects
evaporation even under normal conditions (and probably more so in high-
temperature high-pressure industrial processes).

  2. (b)

Calculations and simulations show that, if the substrate is insulated, the
temperature decreases toward the dew point, while the evaporation rate $E$
tends to zero. This implies that $E$ cannot be deduced by measuring the
external parameters only.

  3. (c)

If the substrate is maintained at a fixed temperature, the heat flux coming
from the substrate supports evaporation at a finite rate. The heat flux in the
liquid is spatially uniform and the temperature profile is linear, which
agrees with measurements [WardStanga01]. Asymptotic formula (92) has been
obtained, relating $E$ to the fluid’s parameters, the layer’s depth, and
relative humidity (which was assumed to be close to unity).

Another conclusion has been drawn for the limit of isothermal evaporation:

  1. (d)

If the temperature is sufficiently high and the vapor density is sufficiently
low, the vapor pressure cannot be matched by that of the liquid (as
illustrated in Fig. 3a). In such cases, the low pressure strains liquid,
encouraging cavitation – and it can be conjectured that ‘near-vacuum boiling’
occurs.
If, however, the vapor pressure does have a match for a liquid state, the
pressure does equilibrate – but the chemical potentials of the two phases are
still different. In this case, evaporation occurs without boiling. The
existence of such regimes on the $(\rho^{(v)},T)$ plane is illustrated in Fig.
3b.

Conclusions (b)–(d) have been drawn using the _pure-fluid_ version of the
diffuse-interface model, where air is approximated by vapor of the same fluid.
These results should rather be viewed as a proof of concept, not an accurate
predictive tool. To obtain the latter, one needs to use the _multicomponent_
version of the DIM [Benilov23a].

There is one result, however, that can be extended to evaporation into air
with no further work – namely, expressions (75) for the heat fluxes toward an
interface in an unbounded space. Denoting these fluxes by $Q^{(air)}$ and
$Q^{(liquid)}$, one can deduce from (75) that

$\frac{Q^{(air)}}{Q^{(liquid)}}=-\sqrt{\frac{\left(\rho
c_{p}\kappa\right)^{(air)}}{\left(\rho c_{p}\kappa\right)^{(liquid)}}},$ (94)

where the minus reflects the opposite directions of the fluxes. Using Refs.
[TheEngineeringToolbox-AirDensity, TheEngineeringToolbox-AirSpecificHeat,
TheEngineeringToolbox-AirThermalConductivity] and [LindstromMallard97] to
estimate the numerator and denominator of the fraction on the right-hand side
of (94), one obtains for water and air at $25^{\circ}\mathrm{C}$

$\left|\frac{Q^{(air)}}{Q^{(liquid)}}\right|\approx 0.0035.$

Evidently, the heat flux coming from air should be accounted for only if it is
the _only_ heat flux (which it indeed is for liquids in an insulated vessel).

It is interesting to speculate how the results of this paper can be extended
to droplets.

Evaporation of a _sessile_ droplet is probably similar to that of a liquid
layer on a fixed-$T$ substrate (since the droplet is small, it cannot change
the substrate’s temperature). A _floating_ droplet, however, is likely to
behave differently. Since the only source of vaporization heat in this case is
air (whose density and thermal conductivity are small), evaporation should be
slow. Evaporative cooling for such droplets should be as important as for a
liquid layer on an insulated substrate.

###### Acknowledgements.

The author is grateful to Daniel Jakubczyk and Tibor Poós for elucidating
discussions of experimental results on evaporation.

##  Appendix A The van der Waals force as described by the diffuse-interface
model

The DIM is based on the following assumptions:

  1. (i)

The long-range attractive intermolecular force (van der Waals force) can be
modelled by a pair-wise isotropic potential $\Phi(r)$ where $r$ is the
distance between two molecules. The net force affecting a given molecule is
the algebraic sum of the forces exerted on it by the other molecules.

  2. (ii)

The spatial scale of the van der Waals force is much smaller than the
interfacial thickness.

Now, consider a three-dimensional fluid with molecular mass $m$, so that
$\rho/m$ is the number density. According to assumption (i), the density of
the collective force exerted by the molecules at a point $\mathbf{r}$, is

$\mathbf{F}(\mathbf{r},t)=\frac{\rho(\mathbf{r},t)}{m}\mathbf{\nabla}\int\frac{\rho(\mathbf{r}^{\prime},t)}{m}\Phi(|\mathbf{r}^{\prime}-\mathbf{r|})\,\mathrm{d}^{3}\mathbf{r}^{\prime}.$
(95)

Then, according to assumption (ii), let the spatial scale of
$\rho(\mathbf{r},t)$ be much larger than that of $\Phi(r)$, in which case
expression (95) can be simplified asymptotically. To do so, change in it
$\mathbf{r}^{\prime}\rightarrow\mathbf{r}^{\prime}+\mathbf{r}$ and then expand
$\rho_{j}(\mathbf{r}^{\prime}+\mathbf{r},t)$ about $\mathbf{r}^{\prime}$,
which yields

$\mathbf{F}_{i}(\mathbf{r},t)=\frac{\rho(\mathbf{r},t)}{m}\mathbf{\nabla}\int\left[\rho(\mathbf{r},t)+\mathbf{r}^{\prime}\cdot\mathbf{\nabla}\rho(\mathbf{r},t)\vphantom{\frac{1}{2}}\right.\\\
+\left.\frac{1}{2}\mathbf{r}^{\prime}\mathbf{r}^{\prime}:\mathbf{\nabla\nabla}\rho(\mathbf{r},t)+\cdots\right]\frac{\Phi(r^{\prime})}{m}\mathrm{d}^{3}\mathbf{r}^{\prime}.$

Since $\Phi(r^{\prime})$ is isotropic, the second integral in the above
expansion vanishes, and one obtains

$\mathbf{F}=\rho\left(C\mathbf{\nabla}\rho+K\mathbf{\nabla}\nabla^{2}\rho+\cdots\right),$
(96)

where

$C=\int\frac{\Phi(r^{\prime})}{m^{2}}\mathrm{d}^{3}\mathbf{r}^{\prime},\qquad
K=\int\frac{r^{\prime
2}}{2}\frac{\Phi(r^{\prime})}{m^{2}}\mathrm{d}^{3}\mathbf{r}^{\prime}.$

After the substitution of (96) into the hydrodynamic equations, the term
involving $C$ can be absorbed into the internal energy – i.e., eliminated by
changing

$e\rightarrow e+\frac{C}{2}\rho,\qquad p\rightarrow p+C\rho^{2}.$

This reflects the fact that the energy of molecular interactions is a kind of
internal energy.

Thus, without loss of generality, one can set in expression (96) $C=0$.
Omitting also the small terms hidden in “$\cdots$”, one obtains

$\mathbf{F}=K\rho\mathbf{\nabla}\nabla^{2}\rho.$

The one-dimensional reduction of this expression represents the van der Waals
force in Eqs. (15)–(16).

##  Appendix B Properties of boundary-value problem (47)–(49)

###  B.1 The Maxwell construction (20)– (21)

To verify Eq. (20), consider the limit $\xi\rightarrow-\infty$ in Eq. (47).
Taking into account boundary condition (48), one immediately obtains (20).

To verify Eq. (21), differentiate (47) and use thermodynamic identity (9) to
obtain

$\frac{\mathrm{d}^{3}\rho^{(0)}}{\mathrm{d}z^{3}}-\frac{\partial
G(\rho^{(0)},T)}{\partial\rho^{(0)}}\frac{\mathrm{d}\rho^{(0)}}{\mathrm{d}z}=0.$

Integrating this equality and fixing the constant of integration via boundary
condition (49), one obtains

$\frac{\mathrm{d}^{2}\rho^{(0)}}{\mathrm{d}\xi^{2}}-G(\rho^{(0)},T)=-G(\rho^{(v.sat)},T).$
(97)

Taking in the above equation the limit $\xi\rightarrow-\infty$ and recalling
boundary condition (48), one recovers Eq. (21).

###  B.2 Derivation of formulae (57)–(58)

First, observe that identity (9) implies that

$G=\frac{\partial\left(\rho G-p\right)}{\partial\rho}.$ (98)

Now, multiply Eq. (97) by $\mathrm{d}\rho^{(0)}/\mathrm{d}\xi$ and integrate.
Using identity (98) and boundary condition (49), one obtains, after
straightforward algebra,

$\frac{1}{2}\left(\frac{\mathrm{d}\rho^{(0)}}{\mathrm{d}\xi}\right)^{2}=\rho^{(0)}\left[G(\rho^{(0)},T)-G(\rho^{(v.sat)},T)\right]-p(\rho^{(0)},T)+p(\rho^{(v.sat)},T).$
(99)

Given that $\rho^{(0)}(\xi)$ is supposed to be a decreasing function, it
follows from (99) that (the superscript (0) omitted)

$\mathrm{d}\rho=-2^{1/2}\sqrt{\rho\left[G(\rho,T)-G(\rho^{(v.sat)},T)\right]-p(\rho,T)+p(\rho^{(v.sat)},T)}\,\mathrm{d}\xi.$

This result allows one to transform (56) into (57).

To obtain expression (58), assume that $T$ is small – hence, $\rho^{(v.sat)}$
is also small, and the main contribution to integral (57) comes from the
region near the lower limit of integration. Since the density is small there,
the general expressions for the chemical potential and pressure can be
replaced with their ideal-gas limits,

$G(\rho,T)\sim T\ln\rho,\qquad p(\rho,T)\sim T\rho,$

and the viscosity can be replaced with its small-density limit,

$\eta(\rho,T)\sim\eta(0,T).$

Since $\rho^{(l.sat)}\gg\rho^{(v.sat)}$, the upper limit of (57) can be
replaced with $\infty$, and one obtains

$A\sim\frac{2^{1/2}\eta(0,T)\,T^{1/2}}{\left(\rho^{(v.sat)}\right)^{5/2}}\int_{1}^{\infty}\frac{\sqrt{\hat{\rho}\ln\hat{\rho}-\hat{\rho}+1}}{\hat{\rho}^{4}}\mathrm{d}\hat{\rho},$

where $\hat{\rho}=\rho/\rho^{(v.sat)}$. The integral in this expression can be
evaluated numerically, yielding Eq. 58).

##  Appendix C The temperature equation

###  C.1 Reduction of the energy equation to the temperature equation

Replace $\partial e/\partial t$ in Eq. (34) with

$\frac{\partial e}{\partial t}=\frac{\partial
e}{\partial\rho}\frac{\partial\rho}{\partial t}+\frac{\partial e}{\partial
T}\frac{\partial T}{\partial t},$

and then eliminate $\partial\rho/\partial t$ using the density equation (32).
Recalling definition (11) of the heat capacity at constant volume, one obtains

$\rho c_{V}\left(\frac{\partial T}{\partial t}+w\frac{\partial T}{\partial
z}\right)+B\frac{\partial w}{\partial z}=\eta\left(\frac{\partial w}{\partial
z}\right)^{2}\\\ +\frac{\partial}{\partial z}\left(\kappa\frac{\partial
T}{\partial z}\right),$ (100)

where

$B=p-\rho^{2}\frac{\partial e}{\partial\rho}$ (101)

characterizes the production (consumption) of thermal energy due to mechanical
compression (expansion) of the fluid. The first term on the right-hand side of
(100) describes the production of heat by viscosity.

Eq. (100) is explicitly resolved with respect to $\partial T/\partial t$ and,
thus, is more convenient for computations than Eq. (34).

###  C.2 The heat conduction equation

Recall ansatz (61) where the temperature field was decomposed into a
background value $T_{0}$ and a small variation $\tilde{T}$. To use it in the
forthcoming formal derivation, introduce a small parameter $\epsilon$, and let

$T=T_{0}+\epsilon\tilde{T}.$

A similar ansatz is applied to the density field,

$\rho=\rho_{0}+\epsilon\tilde{\rho},$

the velocity will be scaled via

$w=\epsilon^{1/2}\tilde{w},$

and the free variables, via

$z=\epsilon^{-1/2}z_{new},\qquad t=\epsilon^{-1}t_{new}.$

Recalling that the original nodimensionalization in Sec. II.3 assumed the
spatial scale to be that of the interface, the stretched variable $z_{new}$
implies that now one considers a region _far above_ or _far below_ the
interface.

Rewriting Eqs. (32)–(33) and (100) in terms of the new variables, omitting the
subscript new, and keeping the leading-order only, one obtains

$\frac{\partial\tilde{\rho}}{\partial t}+\rho_{0}\frac{\partial w}{\partial
z}=0,$ (102) $\frac{\partial}{\partial z}\left[\frac{\partial
p(\rho_{0},T_{0})}{\partial\rho_{0}}\tilde{\rho}+\frac{\partial
p(\rho_{0},T_{0})}{\partial T_{0}}\tilde{T}\right]=0,$ (103)
$\rho_{0}c_{V}(\rho_{0},T_{0})\frac{\partial\tilde{T}}{\partial
t}+B(\rho_{0},T_{0})\frac{\partial w}{\partial
z}=\kappa(\rho_{0},T_{0})\frac{\partial^{2}\tilde{T}}{\partial z^{2}}.$ (104)

The parameter $\epsilon$ has now played its role of an ‘indicator’ of small
terms, as they have all been omitted. One can therefore set $\epsilon=1$, so
that the rescaled variables coincide with those used in the main body of the
paper.

Next, integration of Eq. (103) plus an assumption that the constant of
integration is zero (which amounts to the requirement that the pressure at
infinity does not vary in time) yield

$\tilde{\rho}=-\frac{\partial p}{\partial T}\left[\frac{\partial
p(\rho_{0},T_{0})}{\partial\rho_{0}}\right]^{-1}\tilde{T}.$

Using this equality and Eq. (102), one obtains

$\frac{\partial w}{\partial z}=\frac{1}{\rho_{0}}\frac{\partial
p(\rho_{0},T_{0})}{\partial T_{0}}\left[\frac{\partial
p(\rho_{0},T_{0})}{\partial\rho_{0}}\right]^{-1}\frac{\partial\tilde{T}}{\partial
t},$

so that Eq. (104) becomes (the subscript 0 omitted)

$\left[\rho c_{V}+\frac{B}{\rho}\frac{\partial p}{\partial
T}\left(\frac{\partial
p}{\partial\rho}\right)^{-1}\right]\frac{\partial\tilde{T}}{\partial
t}=\kappa\frac{\partial^{2}\tilde{T}}{\partial z^{2}}.$ (105)

Finally, definition (101) of $B$ and definition (13) of $c_{P}$ help one to
reduce (105) to the standard heat conduction equations, as required.

##  Appendix D Proof of identity (89)

By definition, the vaporization heat is

$\Delta h(T)=h(\rho^{(v.sat)},T)-h(\rho^{(l.sat)},T),$ (106)

where $\rho^{(v.sat)}$ and $\rho^{(l.sat)}$ are the densities of the saturated
vapor and liquid, respectively, and the enthalpy $h$ is given by (12).
Recalling definition (8) of the chemical potential (which implies that
$h=G+Ts$), one can transform (106) into

$\Delta h(T)=G(\rho^{(v.sat)},T)+Ts(\rho^{(v.sat)},T)\\\
-G(\rho^{(l.sat)},T)-Ts(\rho^{(l.sat)},T).$

Recalling equality (21) (the second part of the Maxwell construction), one
obtains (89), as required.