Nonisothermal evaporation

TL;DR

This paper investigates how temperature variations influence liquid evaporation rates on substrates, revealing that thermal conditions significantly alter evaporation dynamics and can be modeled using a diffuse-interface approach.

Contribution

It introduces a nonisothermal model for evaporation, highlighting the impact of substrate thermal conditions on evaporation rates, which was not previously well-understood.

Findings

Evaporation rate depends on substrate thermal conditions.

Insulation slows evaporation, potentially to zero.

Fixed temperature substrates sustain finite evaporation rates.

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.

Equations

$$\Delta T=\frac{\Delta h}{c_P^{(l)}}.$$

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

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

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

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

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

$$\frac{\partial e}{\partial T}=T\frac{\partial s}{\partial T}.$$

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

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

$$\frac{\partial p}{\partial \rho }=\rho \frac{\partial G}{\partial \rho },$$

$$\frac{\partial p}{\partial T}=-{\rho }^2\frac{\partial s}{\partial \rho }.$$

$$c_V=\frac{\partial e}{\partial T},$$

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

$$c_P={\left(\frac{\partial h}{\partial T}\right)}_{p=\mathrm{const}}.$$

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

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial z}\underset{\text{mass flux}}{\underbrace{\left(\rho w\right)}}=0,$$

$$\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}}},$$

$$\begin{gathered} \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{wK\rho \frac{{\partial }^3\rho }{\partial z^3}}}, \end{gathered}$$

$$\eta =\frac{4}{3}{\mu }_s+{\mu }_b,$$

$$\frac{\partial w}{\partial z}\to 0\text{as}z\to +\infty ,$$

$$T\to T^{(v)}\text{as}z\to +\infty ,$$

$$\rho \to {\rho }^{(v)}\text{as}z\to +\infty .$$

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

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

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

$$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)},$$

$$w=0\text{at}z=0.$$

$$T=T_0\text{at}z=0,$$

$$\frac{\partial T}{\partial z}=0\text{at}z=0.$$

$$\rho ={\rho }_0\text{at}z=0,$$

Full text

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.

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

[TABLE]

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:

[TABLE]

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,

[TABLE]

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$,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

and the chemical potential, or specific free energy, by

[TABLE]

It follows from (6)–(8) that

[TABLE]

[TABLE]

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,

[TABLE]

the specific enthalpy,

[TABLE]

and the heat capacity at constant pressure,

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

$\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,

[TABLE]

and the vapor density and temperature tend to certain values,

[TABLE]

[TABLE]

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:

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

and another corresponding to a thermally insulated substrate,

[TABLE]

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],

[TABLE]

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,

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

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,

[TABLE]

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:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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)

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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:

[TABLE]

[TABLE]

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.

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:

[TABLE]

where

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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)}$.

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,

[TABLE]

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)}$.

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.

III.2 The limit of nearly-saturated vapor

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

[TABLE]

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.,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

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:

[TABLE]

where

[TABLE]

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)

[TABLE]

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

[TABLE]

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.

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:

[TABLE]

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:

[TABLE]

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)$,

[TABLE]

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,

[TABLE]

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

[TABLE]

one can describe the temperature evolution by

[TABLE]

where

[TABLE]

[TABLE]

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

[TABLE]

[TABLE]

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,

[TABLE]

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

[TABLE]

Thus, if

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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 quasisteady – 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

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.,

[TABLE]

In vapor, one should, as before, assume

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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}$

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

A similar ansatz is applied to the density field,

[TABLE]

the velocity will be scaled via

[TABLE]

and the free variables, via

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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