2.1. Thermodynamics

When studying hydrodynamics of a compressible fluid, one has to deal with its thermodynamic properties. In this section they are described briefly and in a self-contained form (for the benefit of readers specializing is incompressible hydrodynamics).

Consider an $N$-component compressible non-ideal fluid, characterised by the temperature $T$ and partial mass densities $\rho _{i}$, where $i=1\ldots N$. The fluid's thermodynamic properties are fully described by two functions: the internal energy $e(\rho _{1}\ldots \rho _{N},T)$ and entropy $s(\rho _{1}\ldots \rho _{N},T)$, both specific, i.e. per unit mass. The dependence of the fluid pressure $p$ on $(\rho _{1}\ldots \rho _{N},T)$, or the equation of state, is defined by

(2.1)\begin{equation} p=\rho\sum_{i}\rho_{i}\left( \frac{\partial e}{\partial\rho_{i}} -T\frac{\partial s}{\partial\rho_{i}}\right) , \end{equation}

(e.g. Giovangigli & Matuszewski Reference Giovangigli and Matuszewski2013), where

(2.2)\begin{equation} \rho=\sum_{i}\rho_{i} \end{equation}

is the total density (here and hereinafter, the summation is implied to be from $1$ to $N$ unless stated otherwise). The partial chemical potentials, in turn, are given by

(2.3)\begin{equation} G_{i}=\frac{\partial( \rho e) }{\partial\rho_{i}}-T\frac {\partial( \rho s) }{\partial\rho_{i}}. \end{equation}

Note that $e(\rho _{1}\ldots \rho _{N},T)$ and $s(\rho _{1}\ldots \rho _{N},T)$ are not fully arbitrary, but should satisfy the fundamental thermodynamic relation, or the Gibbs relation, which can be written in the form

(2.4)\begin{equation} \frac{\partial e}{\partial T}=T\frac{\partial s}{\partial T}. \end{equation}

The equivalence of this equality and the standard form of the Gibbs relation are shown in Appendix A.

Using (2.1)–(2.4), one can derive the identities

(2.5)$$\begin{gather} \frac{\partial p}{\partial T}=\sum_{i}\rho_{i}\frac{\partial G_{i}}{\partial T}+\rho s, \end{gather}$$

(2.6)$$\begin{gather}\frac{\partial p}{\partial\rho_{j}}=\sum_{i}\rho_{i}\frac{\partial G_{i} }{\partial\rho_{j}}, \end{gather}$$

then rewrite (2.6) in the form

(2.7)\begin{equation} G_{i}=\frac{\partial}{\partial\rho_{i}}\left( \sum_{j}\rho_{j}G_{j}-p\right) . \end{equation}

In what follows, one also needs the specific heat capacity at constant volume,

(2.8)\begin{equation} c=\frac{\partial e}{\partial T}, \end{equation}

where the traditional subscript $_{V}$ is omitted. In this paper $c$ is assumed to be positive, which is indeed the case for neutral fluids (Lynden-Bell Reference Lynden-Bell1999).

Define also

(2.9)\begin{equation} a_{i}={-}\frac{\partial e}{\partial\rho_{i}}, \end{equation}

which will be referred to as the generalized van der Waals parameter (of the $i$th species), and

(2.10)\begin{equation} B={-}p-\rho\sum_{i}\rho_{i}a_{i}. \end{equation}

Here $B$ is not one of the standard thermodynamic functions, but is convenient to use when thermodynamics is coupled to hydrodynamics; as seen later, $B$ characterises the release (consumption) of heat due to the fluid's mechanical compression (expansion). Using equation of state (2.1) and definition (2.9) of the van der Waals parameter, one can rearrange expression (2.10) in the form

(2.11)\begin{equation} B=\rho T\sum_{i}\rho_{i}\frac{\partial s}{\partial\rho_{i}}. \end{equation}

2.2. Examples: the Enskog–Vlasov and van der Waals fluids

In the low-density limit, the specific internal energy and entropy of any fluid should tend to those of an ideal gas, i.e.

(2.12a,b)\begin{equation} e\sim\frac{T}{\rho}\sum_{i}c_{i}\rho_{i},\quad s\sim\frac{\ln T}{\rho} \sum_{i}c_{i}\rho_{i}-\frac{1}{\rho}\sum_{i}R_{i}\rho_{i}\ln\rho_{i} \quad \text{as}\quad \rho_{i}\rightarrow0, \end{equation}

where $c_{i}$ is the specific partial heat capacity and $R_{i}=k_{B}/m_{i}$ is the specific gas constant, with $m_{i}$ the molecular mass of the $i$th species and

(2.13)\begin{equation} k_{B}=1.380649\times10^{{-}23}\ \mathrm{m}^{2}\ \mathrm{s}^{{-}2} \mathrm{kg\ K}^{{-}1} \end{equation}

the Boltzmann constant. Note also that one can replace $\ln T\rightarrow \ln T/\bar {T}$ and $\ln \rho \rightarrow \ln \rho /\bar {\rho }$, where $\bar {T}$ and $\bar {\rho }$ are suitable dimensional scales. This would make the arguments of the logarithms non-dimensional (with none of physically measurable parameters depending on $\bar {T}$ and $\bar {\rho }$).

A simple description of non-ideal fluids is delivered by the Enskog–Vlasov (EV) model, according to which

(2.14)$$\begin{gather} e=\frac{T}{\rho}\sum_{i}c_{i}\rho_{i}-\frac{1}{\rho}\sum_{ij}a_{ij}\rho _{i}\rho_{j}, \end{gather}$$

(2.15)$$\begin{gather}s=\frac{\ln T}{\rho}\sum_{i}c_{i}\rho_{i}-\frac{1}{\rho}\sum_{i}R_{i}\rho _{i}\ln\rho_{i}-\varTheta(\rho_{1}\ldots \rho_{N}), \end{gather}$$

where $\varTheta (\rho _{1}\ldots \rho _{N})$ is an arbitrary analytic function and $a_{ij}$ characterises the long-range attraction between molecules of the $i$th and $j$th species (with the implication that $a_{ij}=a_{ji}$). For a pure fluid, $a_{11}$ coincides with the generalized van der Waals parameter defined by (2.9) – hence, the notation.

One can readily verify that expressions (2.14) and (2.15) satisfy the Gibbs relation (2.4). They can be viewed as two-term expansions in small $T$, under an extra assumption that the heat capacity $c$ and the generalized van der Waals parameter $a_{i}$ are independent of $T$ and are zeroth-degree homogeneous functions of $\rho _{i}$ (for a pure fluid, $c$ and $a_{1}$ do not depend on $\rho _{1}$ at all).

The EV fluid model originates from the EV kinetic theory (de Sobrino Reference de Sobrino1967; Grmela Reference Grmela1971) and, as such, is naturally suited for use with the DIM that can be viewed as a hydrodynamic approximation of the EV kinetic equation (Giovangigli Reference Giovangigli2020, Reference Giovangigli2021). Equations (2.14) and (2.15) work very well for inert gases (Benilov & Benilov Reference Benilov and Benilov2019), they will be shown to work reasonably well for water, nitrogen and oxygen (see § 8 of this paper).

The EV involves too many parameters to be used as an illustration of general results – so, in such cases, the simpler van der Waals model will be employed. It is a particular case of (2.14) and (2.15) with

(2.16)\begin{equation} \varTheta={-}\frac{1}{\rho}\sum_{i}R_{i}\rho_{i}\ln\left( 1- {\sum_{j}} b_{j}\rho_{j}\right) , \end{equation}

where $b_{i}$ is the reciprocal of the maximum density of the $i$th species. Physically, $1/b_{i}$ can be interpreted as density of the closest packing.

Given expressions (2.14)–(2.16), definition (2.1) of the pressure and definition (2.3) of the chemical potential yield

(2.17)$$\begin{gather} p=\frac{T}{1- {\displaystyle\sum_{j}} b_{j}\rho_{j}}\sum_{i}R_{i}\rho_{i}-\sum_{ij}a_{ij}\rho_{i}\rho_{j}, \end{gather}$$

(2.18)$$\begin{gather}G_{i}=TR_{i}\ln\frac{\rho_{i}}{1- {\displaystyle\sum_{j}} b_{j}\rho_{j}}+\frac{Tb_{i} {\displaystyle\sum_{j}} R_{j}\rho_{j}}{1- {\displaystyle\sum_{j}} b_{j}\rho_{j}}-2\sum_{j}a_{ij}\rho_{j}+T( R_{i}+c_{i}-c_{i}\ln T) . \end{gather}$$

For a pure fluid ($N=1$), (2.17) reduces to the classical van der Waals equation of state (van der Waals Reference van der Waals1893).

2.3. Hydrodynamics

Consider a fluid flow characterised by the species densities $\rho _{i}(\boldsymbol {r},t)$, mass-averaged velocity $\boldsymbol {v}(\boldsymbol {r},t)$ and temperature $T(\boldsymbol {r},t)$, where $\boldsymbol {r}=[ x,y,z]$ is the position vector and $t$ the time. Let the species be affected by forces $\boldsymbol {F}_{i}$ (which will be later identified with the van der Waals forces), and the fluid as a whole, affected by viscosity. The shear viscosity $\mu _{s}$ and bulk viscosity $\mu _{b}$ depend generally on $\rho _{i}$ and $T$.

Let the flow be governed by following equations:

(2.19)$$\begin{gather} \frac{\partial\rho_{i}}{\partial t}+\boldsymbol{\nabla} \boldsymbol{{\cdot}}( \rho _{i}\boldsymbol{v}+\boldsymbol{J}_{i}) =0, \end{gather}$$

(2.20)$$\begin{gather}\frac{\partial( \rho\boldsymbol{v}) }{\partial t}+\boldsymbol{\nabla} \boldsymbol{{\cdot}}( \rho\boldsymbol{vv}) =\boldsymbol{\nabla} \boldsymbol{{\cdot}}\boldsymbol{\varPi }+\sum_{i}\rho_{i}\boldsymbol{F}_{i}-\boldsymbol{\nabla} p, \end{gather}$$

(2.21)$$\begin{gather} \frac{\partial}{\partial t}\left( \rho e+\frac{1}{2}\rho\vert \boldsymbol{v}\vert^{2}\right) +\boldsymbol{\nabla} \boldsymbol{{\cdot}}\left[ \left( \rho e+\frac{1}{2}\rho\vert \boldsymbol{v}\vert^{2}+p\right) \boldsymbol{v}-\boldsymbol{\varPi}\boldsymbol{{\cdot}}\boldsymbol{v}+\boldsymbol{Q}\right] \nonumber\\ =\sum_{i}\boldsymbol{F}_{i}\boldsymbol{{\cdot}}(\rho_{i}\boldsymbol{v}+ \boldsymbol{J}_{i}). \end{gather}$$

Here, the viscous stress tensor is

(2.22)\begin{equation} \boldsymbol{\varPi}=\mu_{s}[ \boldsymbol{\boldsymbol{\nabla} } \boldsymbol{v}+( \boldsymbol{\boldsymbol{\nabla} }\boldsymbol{v})^{\text{tr}}-\tfrac{2}{3}\boldsymbol{I}( \boldsymbol{\boldsymbol{\nabla} } \boldsymbol{{\cdot}}\boldsymbol{v}) ] +\mu_{b}\boldsymbol{I}( \boldsymbol{\boldsymbol{\nabla} }\boldsymbol{{\cdot}}\boldsymbol{v}) , \end{equation}

where the dotless product of two vectors (e.g. $\boldsymbol {\boldsymbol {\nabla }} \boldsymbol {v}$) produces a second-order tensor and the superscript $^{\text {tr}}$ denotes transposition. The diffusion fluxes $\boldsymbol {J}_{i}$ and the heat flux $\boldsymbol {Q}$ are related to the forces $\boldsymbol {F}_{i}$, temperature $T$ and chemical potentials $G_{i}$ by

(2.23)$$\begin{gather} \boldsymbol{J}_{i}=\sum_{j}D_{ij}\left[ \boldsymbol{F}_{j}-T\boldsymbol{\nabla} \left( \frac{G_{j}}{T}\right) \right] -\frac{\zeta_{i}}{T}\boldsymbol{\nabla} T, \end{gather}$$

(2.24)$$\begin{gather}\boldsymbol{Q}=\sum_{j}\zeta_{j}\left[ \boldsymbol{F}_{j}-T\boldsymbol{\nabla} \left( \frac{G_{j}}{T}\right) \right] -\frac{\kappa}{T}\boldsymbol{\nabla} T, \end{gather}$$

where $D_{ij}$, $\zeta _{i}$ and $\kappa$ are the transport coefficients (all three depend generally on $\rho _{i}$ and $T$).

To understand the physical meaning of expressions (2.23) and (2.24), we rewrite them in the form

(2.25)$$\begin{gather} \boldsymbol{J}_{i}=\sum_{j}D_{ij}\boldsymbol{F}_{j}- \underset{\text{diffusion} }{\underbrace{\sum_{j}D_{ij}^{\prime}\boldsymbol{\nabla} \rho_{j}}} -\underset{\text{Soret effect}}{\underbrace{\zeta_{i}^{\prime} \rho\boldsymbol{\nabla} T}}, \end{gather}$$

(2.26)$$\begin{gather}\boldsymbol{Q}=\sum_{j}\zeta_{j}\boldsymbol{F}_{j}- \underset{\text{Dufour effect}}{\underbrace{\sum_{j}\zeta_{j}\sum_{k}\frac{\partial G_{j}} {\partial\rho_{k}}\boldsymbol{\nabla} \rho_{k}}}-\underset{\text{heat conduction}}{\underbrace{\kappa^{\prime}\boldsymbol{\nabla} T}}, \end{gather}$$

where

(2.27a,b)$$\begin{gather} D_{ik}^{\prime}=\sum_{j}D_{ij}\frac{\partial G_{j}}{\partial\rho_{k}} ,\quad \zeta_{i}^{\prime}=\frac{\zeta_{i}}{T\rho}+\sum_{j}\frac{D_{ij}}{\rho} \left( \frac{\partial G_{j}}{\partial T}-\frac{G_{j}}{T}\right) , \end{gather}$$

(2.28)$$\begin{gather} \kappa^{\prime}=\frac{\kappa}{T}+\sum_{j}\zeta_{j}\left( \frac{\partial G_{j}}{\partial T}-\frac{G_{j}}{T}\right) \end{gather}$$

are the standard diffusivities, thermodiffusivities and thermal conductivity, respectively. The second term in expression (2.25) corresponds to the classical Fick law (the fluxes depend linearly on the density gradients), and the last term in expression (2.26) characterises heat conduction described by the Fourier law. The last term in (2.25) describes the Soret effect ($\boldsymbol {\nabla } T$ gives rise to diffusion) and the second term in (2.26), the Dufour effect ($\boldsymbol {\nabla } \rho _{j}$ gives rise to heat conduction).

The same four effects – diffusion, heat conduction, the Soret and Dufour effects – are described, obviously, by the original expressions (2.23) and (2.24), albeit in a form where the terms cannot be matched to a single effect each.

One might think that representing the fluxes in terms of $\boldsymbol {\nabla } \rho _{j}$ would be more natural than using $\boldsymbol {\nabla } ( G_{j}/T)$. Observe, however, that the coefficient of $( \boldsymbol {\nabla } T) /T$ in expression (2.23) coincides with the coefficient of $[ \boldsymbol {F}_{j}-T\boldsymbol {\nabla } ( G_{j}/T) ]$ in (2.24). This symmetry reflects the so-called Onsager reciprocal relations (Ferziger & Kaper Reference Ferziger and Kaper1972), which also imply that

(2.29)\begin{equation} D_{ij}=D_{ji}, \end{equation}

i.e. the diffusion of an $i$th species in a $j$th species occurs the same way as that of the $j$th species in the $i$th species.

It should also be assumed that the extended transport matrix,

(2.30)\begin{equation} D_{ij}^{(ext)}= \begin{bmatrix} & & & \zeta_{1}\\ & D_{ij} & & \vdots\\ & & & \zeta_{N}\\ \zeta_{1} & \ldots & \zeta_{N} & \kappa \end{bmatrix}, \end{equation}

is positive semidefinite ($D_{ij}^{(ext)}\succeq 0$), i.e.

(2.31)\begin{equation} \sum_{i=1}^{N+1}\sum_{j=1}^{N+1}d_{i}\,D_{ji}^{(ext)}d_{j}\geq0 \end{equation}

for all $( N+1)$-dimensional arrays $d_{i}$. As seen later, this property is essential for the entropy principle to hold.

Furthermore, the transport coefficients should be such that

(2.32a,b)\begin{equation} \sum_{i}D_{ij}=0,\quad \sum_{i}\zeta_{i}=0. \end{equation}

As a result, the density (2.19) and expressions (2.23) for the diffusion fluxes imply that

(2.33)\begin{equation} \frac{\partial\rho}{\partial t}+\boldsymbol{\nabla} \boldsymbol{{\cdot}}( \rho \boldsymbol{v}) =0, \end{equation}

where $\rho$ is the total density given by (2.2). Observe that, for a pure fluid, restrictions (2.32a,b) can be satisfied only if $D_{11}=0$ and $\zeta _{1}=0$, which means that pure fluids neither diffuse nor thermodiffuse.

Equations (2.19)–(2.24) have been first derived from the thermodynamics of irreversible processes (Meixner Reference Meixner1941). They were also derived from statistical mechanics (Bearman & Kirkwood Reference Bearman and Kirkwood1958; Mori Reference Mori1958) and non-equilibrium statistical thermodynamics (Keizer Reference Keizer1987); see (Giovangigli Reference Giovangigli1999) for more references. A derivation of the small-density version of (2.19)–(2.24) from the Boltzmann kinetic equation can be found in any textbook on kinetic theory (e.g. Ferziger & Kaper Reference Ferziger and Kaper1972).

In all these cases, the derived expressions for the transport coefficients automatically satisfy the Onsager relations and the rest of the properties listed above.

Note also that a reduction of the above equations for a binary fluid with no Soret and Dufour effects ($N=2$, $\zeta _{i}=0$) was used by Liu, Amberg & Do-Quang (Reference Liu, Amberg and Do-Quang2016) to show that such a model is able to describe the phase equilibrium for a real binary mixture of $\mathrm {CO}_{2}$ and ethanol.

2.4. Alternative forms of the momentum and energy equations

Since the transport fluxes (2.23) and (2.24) are expressed in terms of $G_{i}$ and $T$ (not $\rho _{i}$ and $T$), it is convenient to do the same for the pressure gradient in the momentum equation (2.20). Recalling identities (2.6)–(2.5) which imply that

(2.34)\begin{equation} \boldsymbol{\nabla} p=\sum_{i}\rho_{i}\boldsymbol{\nabla} G_{i}+\rho s\boldsymbol{\nabla} T, \end{equation}

and using (2.33) to simply the left-hand side of the momentum equation (2.20), one reduces it to

(2.35)\begin{equation} \rho\left[ \frac{\partial\boldsymbol{v}}{\partial t}+( \boldsymbol{v} \boldsymbol{{\cdot}}\boldsymbol{\nabla} ) \boldsymbol{v}\right] =\boldsymbol{\nabla} \boldsymbol{{\cdot}}\boldsymbol{\varPi}+\sum_{i}\rho_{i} ( \boldsymbol{F}_{i}-\boldsymbol{\nabla} G_{i}) -\rho s\boldsymbol{\nabla} T. \end{equation}

The energy equation, in turn, can be conveniently rewritten in terms of the temperature (which is a measurable quantity, unlike the internal energy $e$). Replacing, thus, in (2.21),

(2.36)\begin{equation} \frac{\partial e}{\partial t}=\sum_{i}\frac{\partial e}{\partial\rho_{i}} \frac{\partial\rho_{i}}{\partial t}+\frac{\partial e}{\partial T} \frac{\partial T}{\partial t}, \end{equation}

one should use the density equation to eliminate $\partial \rho _{i}/\partial t$. Then using (2.19) and (2.35), and recalling identities (2.8), (2.9) and (2.11), one obtains

(2.37)\begin{equation} \rho c\left( \frac{\partial T}{\partial t}+\boldsymbol{v}\boldsymbol{{\cdot}}\boldsymbol{\nabla} T\right) +\boldsymbol{\nabla} \boldsymbol{{\cdot}}\boldsymbol{Q}=\boldsymbol{\varPi}:( \boldsymbol{\nabla} \boldsymbol{v}) +B\boldsymbol{\nabla} \boldsymbol{{\cdot}}\boldsymbol{v}+\sum_{i}( \boldsymbol{F}_{i}-\rho a_{i}\boldsymbol{\nabla} ) \boldsymbol{{\cdot}}\boldsymbol{J}_{i},\end{equation}

where the symbol ‘$:$’ denotes the double scalar product of two tensors.

The first term on the right-hand side of (2.37) describes heat production by viscosity and the second term that by fluid compression or expansion.

2.5. The van der Waals force

Assume that a molecule of a $j$th species exerts on a molecule of an $i$th species an isotropic force with a potential $\varPhi _{ij}(r)$, where $r=( x^{2}+y^{2}+z^{2})^{1/2}$. Assuming for simplicity that the fluid is unbounded, one can write the mass-averaged force affecting the $i$th species in the hydrodynamic equations (2.20) and (2.21) in the form

(2.38)\begin{equation} \boldsymbol{F}_{i}(\boldsymbol{r},t)=\boldsymbol{\nabla} \sum_{j}\int\frac{\rho _{j}(\boldsymbol{r}^{\prime},t)}{m_{i}m_{j}} \varPhi_{ij}(|\boldsymbol{r}^{\prime}- \boldsymbol{r|})\,\mathrm{d}^{3}\boldsymbol{r}^{\prime}, \end{equation}

where $m_{i}$ is the molecular mass and the integration is implied to be over the whole space. To guarantee the convergence of the integral in (2.38) and those arising later, the potential $\varPhi _{ij}(r)$ is assumed to decay exponentially as $r\rightarrow \infty$.

Next, let the spatial scale of $\rho (\boldsymbol {r},t)$ be much larger than that of $\varPhi _{ij}(r)$, in which case expression (2.38) can be simplified asymptotically. To do so, change in it $\boldsymbol {r}^{\prime }\rightarrow \boldsymbol {r}^{\prime }+\boldsymbol {r}$ and then expand $\rho _{j}(\boldsymbol {r}^{\prime }+\boldsymbol {r},t)$ about $\boldsymbol {r}^{\prime }$, which yields

(2.39)\begin{align} \boldsymbol{F}_{i}(\boldsymbol{r},t)=\boldsymbol{\nabla} \sum_{j}\int\left[ \rho _{j}(\boldsymbol{r},t)+\boldsymbol{r}^{\prime}\boldsymbol{{\cdot}}\boldsymbol{\nabla} \rho_{j} (\boldsymbol{r},t)+\frac{1}{2}\boldsymbol{r}^{\prime}\boldsymbol{r}^{\prime} :\boldsymbol{\nabla}\boldsymbol{\nabla} \rho_{j}(\boldsymbol{r},t)+\cdots \right] \frac{\varPhi _{ij}(r^{\prime})}{m_{i}m_{j}}\,\mathrm{d}^{3}\boldsymbol{r}^{\prime}. \end{align}

Given the isotropy of $\varPhi _{ij}(r^{\prime })$, the second integral in the above expansion vanishes, and one obtains

(2.40)\begin{equation} \boldsymbol{F}_{i}=\sum_{j}W_{ij}\boldsymbol{\nabla} \rho_{j}+\sum_{j}K_{ij} \boldsymbol{\nabla} \nabla^{2}\rho_{j}+\cdots , \end{equation}

where

(2.41a,b)\begin{equation} W_{ij}=\int\frac{\varPhi_{ij}(r^{\prime})}{m_{i}m_{j}}\,\mathrm{d}^{3} \boldsymbol{r}^{\prime},\quad K_{ij}=\int\frac{r^{\prime2}}{2}\frac{\varPhi _{ij}(r^{\prime})}{m_{i}m_{j}}\,\mathrm{d}^{3}\boldsymbol{r}^{\prime}. \end{equation}

Since Newton's third law implies that $\varPhi _{ij}=\varPhi _{ji}$, the matrices $W_{ij}$ and $K_{ij}$ are symmetric.

Once expansion (2.40) is substituted into the hydrodynamic equations (2.20) and (2.21), its first term can be absorbed into the internal energy, i.e. eliminated by the change

(2.42a,b)\begin{equation} e\rightarrow e+\frac{1}{2\rho}\sum_{ij}W_{ij}\rho_{i}\rho_{j},\quad G_{i}\rightarrow G_{i}+\sum_{j}W_{ij}\rho_{j}. \end{equation}

This does not come as a surprise, as the energy associated with potential interactions of molecules can be viewed as a kind of internal energy; in fact, the second term of expansion (2.40) could also be (and sometimes is) absorbed into $e$. This is not done in this paper, however, as it would make $e$ a functional (instead of a function), with the implication that all the thermodynamic definitions and identities in § 2.1 would need to be rewritten in terms of functional derivatives.

Thus, without loss of generality, one can set in expression (2.40), $W_{ij}=0$. Omitting also the small terms hidden in ‘$\ldots$’, one obtains

(2.43)\begin{equation} \boldsymbol{F}_{i}=\sum_{j}K_{ij}\boldsymbol{\nabla} \nabla^{2}\rho_{j}, \end{equation}

which is the multicomponent extension of the standard DIM formula for the van der Waals force (e.g. Pismen & Pomeau Reference Pismen and Pomeau2000). The matrix $K_{ij}$ is the extension of the so-called Korteweg parameter for pure fluids, and it will be referred to as the Korteweg matrix. It should be positive definite, $K_{ij}\succ 0$, as the van der Waals force should be attractive, not repulsive.

Since the original representation (2.38) was a phenomenological model and, thus, the pairwise potential $\varPhi _{ij}$ cannot be measured, the Korteweg matrix should be viewed as a set of adjustable parameters. As seen later, they can be deduced from the measurements of the equation of state and surface tension. One should keep in mind, however, that the Korteweg matrix should not depend on the temperature. Such a dependence would be physically unjustified, as the intermolecular attraction (characterised by $K_{ij}$) should not depend on the molecules’ velocities (characterised by $T$).

2.6. Boundary conditions at a solid wall

Let the fluid occupy a domain $\mathcal {D}$, bounded by a smooth surface $\partial \mathcal {D}$. For simplicity, the so-called Navier slip is disallowed in this work, so that the boundary condition for the velocity is

(2.44)\begin{equation} \boldsymbol{v}=\boldsymbol{0}\ \text{at}\ \boldsymbol{r}\in\partial\mathcal{D}. \end{equation}

To ensure mass conservation, one should require

(2.45)\begin{equation} \boldsymbol{n}\boldsymbol{{\cdot}} \boldsymbol{J}_{i}=0\ \text{at}\ \boldsymbol{r}\in\partial\mathcal{D}, \end{equation}

where $\boldsymbol {n}$ is the outward-pointing unit normal to $\partial \mathcal {D}$.

The boundary condition for the temperature, in turn, depends on the problem at hand. Since this paper is concerned, inter alia, with the entropy principle and energy conservation, it will be assumed that no heat escapes through the boundary,

(2.46)\begin{equation} \boldsymbol{n}\boldsymbol{{\cdot}} \boldsymbol{Q}=0\ \text{at}\ \boldsymbol{r}\in\partial\mathcal{D}. \end{equation}

Boundary conditions (2.44)–(2.46) would be sufficient for the standard compressible multicomponent hydrodynamics, but the DIM requires an extra condition (due to the presence of higher-order derivatives of the density field in expression (2.43)).

The most common version of such a condition – prescribing a linear combination of the boundary value of the density and its normal gradient – ascends to the paper by Cahn & Hilliard (Reference Cahn and Hilliard1958). In application to pure fluids, the Neumann reduction of this condition was proposed by Seppecher (Reference Seppecher1996) and the Dirichlet reduction by Pismen & Pomeau (Reference Pismen and Pomeau2000). As shown by Benilov (Reference Benilov2020b), the latter follows from the assumptions under which the whole DIM is derived (pairwise intermolecular interactions, slowly varying density field), and so it is used in the present paper. Thus, we require that

(2.47)\begin{equation} \rho_{i}=\rho_{0,i}\ \text{at}\ \boldsymbol{r}\in\partial\mathcal{D},\end{equation}

where the constant $\rho _{0,i}$ is specific to the fluid/solid combination under consideration. The general version of the boundary condition for $\rho _{i}$ is discussed briefly in Appendix B.

To clarify the physical meaning of condition (2.47), consider the van der Waals forces acting on a fluid molecule located infinitesimally close to the wall: the solid attracts it towards the wall, while the other fluid molecules pull it away. The former force is fixed, whereas the latter grows with the near-wall density – so that the balance is achieved when the density assumes a certain value – namely, the parameter $\rho _{0,i}$ in condition (2.47). This argument suggests that a smaller value of $\rho _{0,i}$ corresponds to a hydrophobic wall (characterised by a large contact angle) and larger $\rho _{0,i}$ to a hydrophilic one.

According to its physical meaning, $\rho _{0,i}$ does not depend on the temperature. As shown later, its value can be deduced from a measurement of the contact angle.

3.1. Conservation laws and the H-theorem

It can be verified that the governing equations and boundary conditions introduced above conserve the mass of each species

(3.1)\begin{equation} M_{i}=\int_{\mathcal{D}}\rho_{i}\,\mathrm{d}^{3}\boldsymbol{r}, \end{equation}

and the total energy

(3.2)\begin{equation} E=\int_{\mathcal{D}+ }\left[ \rho e+\frac{1}{2}\rho\vert \boldsymbol{v} \vert^{2}+\frac{1}{2}\sum_{ij}K_{ij}( \boldsymbol{\nabla} \rho _{i}) \boldsymbol{{\cdot}}( \boldsymbol{\nabla} \rho_{j}) \right]\, \mathrm{d}^{3}\boldsymbol{r}.\end{equation}

The three terms in expression (3.2) represent the internal, kinetic and van der Waals energies.

The governing equations and boundary conditions also satisfy an H-theorem, reflecting the fact that the net entropy of a fluid in a thermally insulated container cannot decrease. To prove this, consider the following combination of the governing equations:

(3.3) \begin{equation} (2.37)+T\rho\sum_{i}\frac{\partial s}{\partial\rho_{i}} \times(2.19)+Ts\times (2.26). \end{equation}

After straightforward algebra involving the use of the thermodynamic identities presented in § 2.1 and expressions (2.23) and (2.24) for the fluxes, one obtains

(3.4)\begin{align} &\frac{\partial( \rho s)}{\partial t}+\boldsymbol{\nabla} \boldsymbol{{\cdot}}\left( \rho s\boldsymbol{v}-\sum_{i}\frac{G_{i}}{T}\boldsymbol{J}_{i}+ \frac{\boldsymbol{Q}}{T}\right) =\frac{\boldsymbol{\varPi}: ( \boldsymbol{\nabla} \boldsymbol{v}) }{T}\nonumber\\ &\quad +\left\{ \sum_{ij}D_{ij}\left\vert \boldsymbol{F}_{j}-T\boldsymbol{\nabla} \left( \frac{G_{j}}{T}\right) \right\vert^{2}+2\sum_{i}\zeta_{i}\left[ \boldsymbol{F}_{i}-T\boldsymbol{\nabla} \left( \frac{G_{i}}{T}\right) \right] \boldsymbol{{\cdot}}\left( -\frac{\boldsymbol{\nabla} T}{T}\right) \phantom{\left\vert -\frac{\boldsymbol{\nabla} T}{T}\right\vert^{2}}\right. \nonumber\\ &\quad +\left.\vphantom{\sum_{ij}\left\vert \left( \frac{G_{j}}{T}\right) \right\vert^{2}}\kappa T\left\vert -\frac{\boldsymbol{\nabla} T}{T} \right\vert^{2}\right\} . \end{align}

The first term on the right-hand side of this equation is non-negative due to the (easily verifiable) identity

(3.5)\begin{equation} \boldsymbol{\varPi}:( \boldsymbol{\nabla} \boldsymbol{v}) =\tfrac{1}{2}\mu _{s}\vert ( \boldsymbol{\nabla} \boldsymbol{v}) +( \boldsymbol{\nabla} \boldsymbol{v})^{\text{tr}}-\tfrac{2}{3}\boldsymbol{I} \boldsymbol{\nabla}\boldsymbol{{\cdot}}\boldsymbol{v}\vert^{2}+\mu_{b}( \boldsymbol{\nabla} \boldsymbol{{\cdot}}\boldsymbol{v})^{2}, \end{equation}

whereas the expression in curly brackets is non-negative because the extended transport matrix is positive semidefinite (see § 2.3). Thus, integrating (3.4) over $\mathcal {D}$, and using boundary conditions (2.44)–(2.46), one obtains

(3.6)\begin{equation} \frac{\mathrm{d}S}{\mathrm{d}t}\geq0, \end{equation}

where

(3.7)\begin{equation} S=\int_{\mathcal{D}}\rho s\,\mathrm{d}^{3}\boldsymbol{r}. \end{equation}

Inequality (3.6) is the desired H-theorem.

It follows from (3.4) and (3.5) that the exact equality in (3.6) can only hold if the velocity field is spatially uniform; together with the no-flow boundary conditions, this requirement amounts to $\boldsymbol {v}=\boldsymbol {0}$ (i.e. the fluid is static).

3.2. Stability via the entropy principle

The most common way to examine the stability of a steady solution of a set of equations consist in linearising these equations with respect to a small perturbation, assuming the harmonic dependence of the perturbation on $t$, and solving the resulting eigenvalue problem. In the problem at hand, however, it is much simpler to examine stability using the entropy principle.

If, at a certain steady state, the total entropy $S$ has a local maximum constrained by the conditions of fixed energy $E$ and mass $M_{i}$, this state is stable. The inverse is also true: if $S$ does not have a maximum, the corresponding state is unstable, because a perturbed solution with a higher entropy cannot evolve ‘back’. Neutrally stable oscillations are also prohibited by the entropy principle – hence, the system can only evolve further away, towards a steady state where the entropy does have a maximum.

Let a fluid be enclosed in a container (which can be later assumed to be infinitely large, if need be) and seek a maximum of $S$, constrained by the conditions of fixed $M_{i}$ and $E$. This problem amounts to finding the stationary points of the functional

(3.8)\begin{equation} H[\rho_{1}\ldots \rho_{N},T,\boldsymbol{v}]=S+\sum_{i}\eta_{i}M_{i}+\lambda E, \end{equation}

where $\lambda$ and $\mu _{i}$ are the Lagrange multipliers, and $S$, $E$ and $M_{i}$ are given by (3.7), (3.2) and (3.1), respectively. Varying $H$ with respect to $\boldsymbol {v}$ and equating the variation to zero, one obtains

(3.9)\begin{equation} \boldsymbol{v}=\boldsymbol{0}, \end{equation}

i.e. a steady state must be (unsurprisingly) static. Next, varying $H$ with respect to $T$, one obtains

(3.10)\begin{equation} \frac{\partial s}{\partial T}+\lambda\frac{\partial e}{\partial T}=0. \end{equation}

Comparison of this equality with the Gibbs relation (2.4) yields

(3.11)\begin{equation} \lambda={-}\frac{1}{T}. \end{equation}

Since $\lambda$ is a constant, (3.11) implies that $T$ is spatially uniform, i.e. a steady state must be isothermal.

Finally, varying $H$ with respect to $\rho _{i}$, recalling expressions (3.9)–(3.11) for $\boldsymbol {v}$ and $\lambda$, and keeping in mind definition (2.3) of the chemical potential $G_{i}$, one obtains

(3.12)\begin{equation} \frac{1}{T}\left( G_{i}-\sum_{j}K_{ij}\nabla^{2}\rho_{j}\right) +\eta_{i}=0. \end{equation}

This equation describes all steady states of the governing equations, and it will be extensively used in the rest of this paper. The temperature $T$ in (3.12) should be treated as a known parameter, whereas the Lagrange multiplier $\eta _{i}$ is to be deduced from the boundary conditions. The latter will be illustrated for the case of an infinite domain, plus an assumption that the fluid at infinity be spatially uniform and characterised by a coordinate-independent chemical potential $G_{i}=G_{\infty,i}$. In this case, (3.12) yields

(3.13)\begin{equation} \eta_{i}={-}\frac{G_{\infty,i}}{T}, \end{equation}

as required.

Note also that (3.12) can also be recovered by adapting the governing equations for the state of equilibrium. To do so, one should set $\boldsymbol {v}=\boldsymbol {0}$, $T=\text {const.}$ and $\partial \rho _{i}/\partial t=0$ in (2.19) and (2.20) and, recalling expression (2.23) for the diffusion fluxes, obtain (3.12), as required.

To examine a solution of (3.12) for stability, one needs to calculate the second variation of $H$. Omitting the algebra (involving the use of identities (2.3) and (2.4), definition (2.8) of the heat capacity $c$ and expression (3.11) for $\lambda$), one obtains

(3.14)\begin{align} \delta^{2}H=\frac{1}{T}\int_{\mathcal{D}}\left\{ \sum_{ij}\left[ -\frac{\partial G_{i}}{\partial\rho_{j}}( \delta\rho_{i}) ( \delta\rho_{j}) -K_{ij}( \boldsymbol{\nabla} \delta\rho_{i}) \boldsymbol{{\cdot}}( \boldsymbol{\nabla} \delta\rho_{j}) \right] -\frac{\rho c} {T}( \delta T)^{2}-\rho\vert \delta\boldsymbol{v}\vert^{2}\right\}\, \mathrm{d}^{3}\boldsymbol{r}. \end{align}

Evidently, perturbations of the temperature and velocity are negative and can only lower the total entropy – hence, the type of stationary point (maximum versus saddle) is fully determined by the variation of the density. Thus, setting $\delta T=0$ and $\delta \boldsymbol {v}=\boldsymbol {0}$, one obtains

(3.15)\begin{equation} \delta^{2}H=\frac{1}{T}\int_{\mathcal{D}}\sum_{ij}\left[ -\frac{\partial G_{i}}{\partial\rho_{j}}( \delta\rho_{i}) ( \delta\rho _{j}) -K_{ij}( \boldsymbol{\nabla} \delta\rho_{i}) \boldsymbol{{\cdot}}( \boldsymbol{\nabla} \delta\rho_{j}) \right] \,\mathrm{d}^{3}\boldsymbol{r}. \end{equation}

Expression (3.15) is the main tool used in this paper for studying the stability properties of steady states described by (3.12).

5.1. Spatially uniform states

Consider a uniform fluid where there is no flow and all species are in vapour phase. If the temperature drops, one of the species may become overcooled, giving rise to rapid condensation. A similar instability may occur when all or some of the species are in liquid phase and the temperature increases, giving rise to rapid evaporation.

To determine exactly which states are unstable (thus, do not occur in real world), one could perform the usual linear analysis. For a pure fluid ($N=1$), this is a straightforward task yielding the following stability criterion:

(5.1)\begin{equation} \frac{\partial p}{\partial\rho}>0. \end{equation}

Thus, instability occurs if an increase in density lowers the pressure, so that the flow generated by the pressure gradient brings even more fluid to this region.

For an arbitrary $N$, however, the analysis of linearised equations is extremely cumbersome. Instead, it will be examined via the entropy principle, i.e. using expression (3.15). Taking the limit $\mathcal {D}\rightarrow \mathbb {R} ^{3}$ (unbounded fluid), assuming that $\rho _{i}$ is spatially uniform and recalling that $K_{ij}\succ 0$ and $c>0$, one can deduce from (3.15) that $H$ has a maximum if

(5.2)\begin{equation} \frac{\partial G_{i}}{\partial\rho_{j}}\succ0. \end{equation}

This is the standard stability criterion for a spatially uniform state of a multicomponent fluid (Glansdorff & Prigogine Reference Glansdorff and Prigogine1971).

The following four comments are in order.

1. (i) For a physically meaningful $G_{i}(\rho _{1}\ldots \rho _{N},T)$, condition (5.2) holds for a sufficiently rarefied vapour or a sufficiently dense liquid.

2. (ii) The states with marginal stability are sometimes referred to as ‘spinodal points’ and stable vapour as ‘subspinodal vapour’.

3. (iii) Interestingly, the viscosity and transport coefficients do not appear in criterion (5.2). The corresponding effects can only slow the instability down, but not eliminate it.

4. (iv) To reconcile criterion (5.2) with its pure-fluid counterpart (5.1), note that, for $N=1$, $\partial G/\partial \rho$ and $\partial p/\partial \rho$ have the same sign (as implied by identity (2.6)).

5.2. Two-phase saturated states

Consider an interface separating liquid and vapour of the same pure fluid. If they are in equilibrium, their temperatures are equal, and the rest of the parameters satisfy the so-called Maxwell construction (Maxwell Reference Maxwell1875),

(5.3a,b)\begin{equation} G(\rho^{(l)},T)=G(\rho^{(v)},T),\quad p(\rho^{(l)},T)=p(\rho^{(v)},T), \end{equation}

where the superscripts ${(l)}$ and ${(v)}$ mark the parameters of the liquid and vapour, respectively. The former and latter equalities in (5.3a,b) guarantee thermodynamic and mechanical equilibria of the interface, respectively.

Assume also that both phases are stable,

(5.4a,b)\begin{equation} \left( \frac{\partial G}{\partial\rho}\right) _{\rho=\rho^{(l)}} >0,\quad \left( \frac{\partial G}{\partial\rho}\right) _{\rho=\rho^{(v)} }>0, \end{equation}

and the density of the liquid exceeds that of the vapour,

(5.5)\begin{equation} \rho^{(l)}>\rho^{(v)}. \end{equation}

Subject to these conditions, the Maxwell construction (5.3a,b) uniquely determines the saturated densities $\rho ^{(l)}$ and $\rho ^{(v)}$ as functions of $T$.

In what follows, the multicomponent version of the Maxwell construction will be shown to follow from the DIM's entropy principle.

Consider an insulated container with fluid subdivided between two states, liquid and vapour. If the liquid/vapour and fluid/wall interfaces are sufficiently thin, the corresponding full entropy, mass and energy can be approximately written in the form

(5.6)$$\begin{gather} S=V^{(l)}\rho^{(l)} s(\rho_{1}^{(l)}\ldots \rho_{N}^{(l)},T^{(l)})+V^{(v)} \rho^{(v)} s(\rho_{1}^{(v)}\ldots \rho_{N}^{(v)},T^{(v)}), \end{gather}$$

(5.7)$$\begin{gather}M_{i}=V^{(l)}\rho_{i}^{(l)}+V^{(v)}\rho_{i}^{(v)}, \end{gather}$$

(5.8)$$\begin{gather}E=V^{(l)}\rho^{(l)} e(\rho_{1}^{(l)}\ldots \rho_{N}^{(l)},T^{(l)})+V^{(v)} \rho^{(v)} e(\rho_{1}^{(v)}\ldots \rho_{N}^{(v)},T^{(v)}), \end{gather}$$

where $V^{(l)}$ and $V^{(v)}$ are the volumes of the liquid and vapour phases, respectively. We also introduce the full volume of the container,

(5.9)\begin{equation} V=V^{(l)}+V^{(v)}. \end{equation}

The Maxwell construction can be derived by maximizing $S$ subject to the constraints of fixed $M_{i}$, $E$ and $V$ (which are now functions, as opposed to being functionals in the previous subsection). Straightforward calculations show that the maximum of entropy is achieved if $T^{(l)}=T^{(v)}$ (isothermality) and

(5.10)$$\begin{gather} G_{i}(\rho_{1}^{(l)}\ldots \rho_{N}^{(l)},T)=G_{i}(\rho_{1}^{(v)}\ldots \rho_{N} ^{(v)},T), \end{gather}$$

(5.11)$$\begin{gather}p(\rho_{1}^{(l)}\ldots \rho_{N}^{(l)},T)=p(\rho_{1}^{(v)}\ldots \rho_{N}^{(v)},T). \end{gather}$$

The following four comments are in order.

1. (i) Since the solution describing coexistence of two phases satisfies the entropy principle, it is automatically stable.

2. (ii) The Maxwell construction can sometimes yield a solution with negative $V^{(l)}$ or $V^{(v)}$. In such cases, the two-phase equilibrium is irrelevant, and the fluid evolves towards the one-phase state with the same masses of the species and total energy.

3. (iii) The multicomponent Maxwell construction (5.10) and (5.11) comprises $N+1$ equations for $2N$ unknowns – hence, does not uniquely fix the liquid and vapour densities.

   To close set (5.10) and (5.11), one should assume the masses $M_{i}$ and the total energy $E$ to be known and view equalities (5.7)–(5.9) as additional equations. They bring the total number of equations to $2N+3$ (making the ‘expanded’ Maxwell construction appear overdetermined), but in this formulation $V^{(l)}$, $V^{(v)}$ and $T$ should also be viewed as unknowns.

   Physically, if a certain amount of fluid, with a certain amount of energy, is placed in a box, the entropy principle uniquely determines the final temperature and the proportion in which the box is subdivided between the liquid and vapour phases.

4. (iv) The closure of the Maxwell construction described in the previous bullet is, obviously, inapplicable to containers of infinite volume. To understand how conditions (5.10) and (5.11) should be closed in this case, consider the interface between the Earth's atmosphere and ocean. For this setting, one should prescribe the (atmospheric) pressure and composition of dry air above the ocean's surface. With these parameters given, (5.10) and (5.11) yield the saturated moisture content of the air and the saturated amounts of nitrogen, oxygen, etc. dissolved in the oceanic water.

To illustrate the use of the Maxwell construction, consider a van der Waals fluid, whose pressure and chemical potential are described by expressions (2.17) and (2.18). Assume for simplicity that the fluid is pure ($N=1$) and monatomic ($c_{i}=3R_{i}/2$), and let the scales $\varrho$, $P$ and $R$ used for non-dimensionalisation be such that

(5.12a–c)\begin{equation} a_{11}=1,\quad b_{1}=1,\quad R_{1}=1. \end{equation}

If $T<8/27$ (which is a subcritical temperature of the van der Waals fluid), two states exist representing the liquid and vapour phases, with some (but not all) of the states in between being unstable (see figure 1). If the temperature is supercritical, only one phase exists and interfaces do not.

Figure 1. The pressure (a) and chemical potential (b) of a pure van der Waals fluid as functions of the density, for three values of the temperature: (1) $T=0.26$ (subcritical), (2) $T=8/27$ (critical), (3) $T=0.33$ (supercritical). The critical point on curve (2) is marked with a filled circle; the saturation and spinodal points on curve (1) are shown by filled and empty squares, respectively. The labels ‘undersat(urated) vapour’, ‘oversat(urated)v(apour)’, etc. apply only to curve (1).

The simplest model describing a water/air interface is that with $N=2$. The Maxwell construction in this case should be complemented with one extra condition, setting the pressure above the interface equal to its atmospheric value,

(5.13)\begin{equation} p(\rho_{1}^{(v)},\rho_{2}^{(v)},T)=p_{A}. \end{equation}

The difference between a pure and a multicomponent fluid is illustrated in figure 2 for parameters (5.12a–c) and

(5.14a–e)\begin{equation} a_{22}=0.2,\quad a_{12}=0,\quad b_{2}=1,\quad R_{2}=0.6,\quad p_{A}=0.03. \end{equation}

These values reflect a compromise between simplicity, illustrative purposes (the curves with $N=1$ and $N=2$ should be visibly different) and an attempt to loosely match the parameters of water and air (loosely, because the van der Waals model includes few adjustable constants). In particular, parameters (5.14a–e) make the critical temperature of the second species noticeably smaller than that of the first species. As a result, a temperature range exists where the former is definitely vapour, whereas the latter can be either vapour or liquid (as is indeed the case with water under normal conditions).

Figure 2. The saturated densities (a) and pressure (b) versus the temperature, for a pure (solid line) and two-component (dotted line) van der Waals fluid, with the parameters (5.12a–c)–(5.14a–e). The empty circles mark the boiling point, the filled circle marks the critical point. The saturated densities of the second species (of the two-component fluid) are not shown.

The most important difference between a multicomponent fluid at fixed pressure and a pure fluid is that the former exhibits a boiling point. It occurs when the saturated pressure of the liquid-phase species matches the applied (atmospheric) pressure: as a result, the second species is completely replaced by the vapour of the first species, and so the fluid becomes pure. If it is heated beyond the boiling point, the pressure can no longer be kept fixed, but ought to increase (to match the saturated pressure of the first species, which grows with $T$).

One should keep in mind that the small-$T$ part of figure 2 is physically meaningless due to freezing (which is not described by the DIM).

6.1. Flat liquid/vapour interfaces

Consider an unbounded fluid involving liquid and vapour phases in equilibrium, separated by a static flat interface. Its spatial structure is described by (3.12) derived from the entropy principle. To adapt it specifically for a flat interface, let $\rho _{i}$ depend on a single coordinate – say, $z$ – so that the interface is horizontal. Let the fluid below the interface be liquid and above vapour,

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

(6.2)$$\begin{gather}\rho_{i}\rightarrow\rho_{i}^{(v)}\quad \text{as}\ z\rightarrow+\infty. \end{gather}$$

Taking in (3.12) the limit $z\rightarrow +\infty$, one can determine the constant $\eta _{i}$ and rewrite (3.12) in the form

(6.3)\begin{equation} \sum_{j}K_{ij}\frac{\mathrm{d}^{2}\rho_{j}}{\mathrm{d}z^{2}}=G_{i}-G_{i} (\rho_{1}^{(v)}\ldots \rho_{N}^{(v)},T). \end{equation}

The boundary-value problem (6.1)–(6.3) is invariant with respect to the change $z\rightarrow z+\textrm {const.}$ – hence, its solution is not unique. To make it such, an extra condition is needed, say,

(6.4)\begin{equation} \rho_{1}=\tfrac{1}{2}( \rho_{1}^{(l)}+\rho_{1}^{(v)}) \quad \text{at}\ z=0. \end{equation}

As shown before, the parameters $\rho _{i}^{(l)}$ and $\rho _{i}^{(v)}$ are not arbitrary but have to satisfy the Maxwell construction, i.e. conditions (5.10) and (5.11). Thus, they must be intrinsic to the boundary-value problem (6.1)–(6.4) and can be derived directly from it, as a condition for the existence of solution.

Indeed, consider (6.3) in the limit $z\rightarrow -\infty$ and recall (6.1), which immediately yields the first half of the Maxwell construction, (5.10). To derive the second half, consider

(6.5)\begin{equation} \sum_{i}\int_{-\infty}^{\infty}(6.3) \times \frac{\mathrm{d}\rho_{i}}{\mathrm{d}z}\,\mathrm{d}z. \end{equation}

The integral in this equation can be evaluated via identity (2.7), and the constant of integration can be fixed via boundary condition (6.2). After straightforward algebra, one obtains

(6.6)\begin{equation} \frac{1}{2}\sum_{ij}K_{ij}\frac{\mathrm{d}\rho_{i}}{\mathrm{d}z} \frac{\mathrm{d}\rho_{j}}{\mathrm{d}z}=\sum_{i}\rho_{i}[ G_{i} -G_{i}(\rho_{1}^{(v)}\ldots \rho_{N}^{(v)},T)] -p+p(\rho_{1}^{(v)} \ldots \rho_{N}^{(v)},T). \end{equation}

Taking in this equation the limit $z\rightarrow -\infty$ and recalling boundary condition (6.1), one recovers the second part of the Maxwell construction, (5.11), as required.

The boundary-value problem (6.1)–(6.4) was solved numerically for a two-component van der Waals fluid, using the MATLAB function BVP4c based on the three-stage Lobatto IIIa formula (Kierzenka & Shampine Reference Kierzenka and Shampine2001). Examples of numerical solutions with parameters (5.12a–c)–(5.14a–e) and

(6.7a,b)\begin{equation} K_{11}=K_{22}=1,\quad K_{12}=0 \end{equation}

are shown in figure 3. Observe that an increase in temperature makes the interface thicker, and it becomes infinitely thick when the critical temperature is reached.

Figure 3. Comparison between liquid/vapour interfaces in pure (solid line) and two-component (dotted line) van der Waals fluids with parameters (5.12a–c)–(5.14a–e) and (6.7a,b). Curves (1)–(3) correspond to $T=0.17,0.23,0.28132$. The last of these values is the (approximate) boiling point, so that the solid and dotted curves virtually coincide and $\rho _{2}$ is virtually zero.

As shown in Appendix C.1, all solutions describing liquid/vapour interfaces are stable as long as $\mathrm {d}\rho _{i} /\mathrm {d}z$ never vanishes ($\rho _{i}(z)$ is strictly monotonic) for all $i$.

6.2. One-dimensional drops and bubbles

Apart from liquid/vapour interfaces, (3.12) admits spatially localised solutions, describing drops or bubbles floating in vapour or liquid, respectively. In this subsection the simplest – one-dimensional (1-D) – solutions of this kind are discussed, describing a flat layer of increased or decreased density. They both correspond to the following boundary condition:

(6.8)\begin{equation} \rho_{i}\rightarrow\rho_{\infty,i}\ \text{as}\ z\rightarrow\pm\infty. \end{equation}

Here $\rho _{\infty,i}$ is the density of the vapour (liquid). Unlike the interface solutions examined above, the fluid outside the drop (bubble) does not have to be saturated vapour (liquid); it can actually be any stable or metastable state.

The 1-D reduction of the steady state (3.12) and expression (3.13) for $\eta _{i}$ can be written in the form

(6.9)\begin{equation} \sum_{j}K_{ij}\frac{\mathrm{d}^{2}\rho_{j}}{\mathrm{d}z^{2}}=G_{i} -G_{\infty,i}, \end{equation}

where $G_{\infty,i}=G_{i}(\rho _{\infty,1}\ldots \rho _{\infty,N},T)$ is the chemical potential at infinity.

Figure 4 illustrates typical solutions of the boundary-value problem (6.8) and (6.9) computed for a pure fluid. The following features can be observed in panel (a) depicting some 1-D drop solutions.

1. (i) Drop solutions exist only if $\rho _{\infty,1}$ corresponds to oversaturated vapour.

2. (ii) As $\rho _{\infty,1}\rightarrow \rho _{1}^{(v)}+0$, the drop becomes increasingly thick.

3. (iii) Once $\rho _{\infty,1}$ passes $\rho _{1}^{(v)}$, drop solutions cease to exist. Physically, this is because liquid drops surrounded by undersaturated vapour evaporate.

Figure 4. Solutions of boundary-value problem (6.8) and (6.9) for a pure van der Waals fluid with $T=0.2$. (a) Drop solutions with $\rho _{\infty,1}=0.1, 0.05, 0.03587, 0.0358485$ (curves (1)–(4), respectively). (b) Bubble solutions with $\rho _{\infty,1} =0.63, 0.68, 0.7286, 0.728671$ (curves (1)–(4), respectively). The labels ‘u(ndersaturated) v(apour)’, ‘o(versaturated) v(apour)’, etc. characterise $\rho _{\infty,1}$.

Similar tendencies have been observed for 1-D bubble solutions, illustrated in figure 4(b).

Most importantly, all drops (bubbles) in oversaturated vapour (undersaturated liquid) are likely to be unstable. For pure fluids ($N=1$), the instability can be proven rigorously; as shown in Appendix C.2, the entropy has a saddle (not maximum) on these solutions.

To understand how slightly perturbed drops and bubbles evolve, they were simulated numerically using the simplest asymptotic version of the DIM, the one based on (4.18), (4.19) and (4.9). It is adapted for $N=1$ and a single spatial coordinate in Appendix D, which also outlines the numerical method used.

Various initial conditions have been simulated for the van der Waals pure fluid, with only two patterns of evolution observed. If the initial condition includes a steady-drop solution $\rho _{1}^{(sd)}(z)$ plus some extra fluid, spontaneous condensation is typically triggered off, giving rise to two shock waves propagating away from the drop's centre. This behaviour is illustrated in figure 5(a) for

(6.10a,b)\begin{equation} T=0.2,\quad \rho_{\infty,1}=0.05, \end{equation}

and the initial condition

(6.11)\begin{equation} \rho_{1}=\rho_{1}^{(sd)}(1.01\,z),\quad w=0\ \text{at} t=0,\end{equation}

where $w$ is the vertical velocity (the other two velocity components are zero due to the 1-D nature of the flow).

Figure 5. Evolution of a perturbed drop, described by (D3) and (D4) and boundary conditions (D5a,b) and (D6a,b), for a pure van der Waals fluid and parameters (6.10a,b). (a) Initial condition (6.11) (the steady solution plus extra fluid); the time $t$ of a snapshot and the corresponding curve number $n$ are inter-related via $t=5000\,n$. (b) Initial condition (6.12) (the steady solution minus extra fluid); $t=500\,n$.

If, however, the initial condition includes less fluid than $\rho _{1}^{(nm)}(z)$, the drop would typically evaporate. This pattern is illustrated in figure 5(b) for parameters (6.10a,b) and the following initial condition:

(6.12)\begin{equation} \rho_{1}=\rho_{1}^{(sd)}(0.99\,z),\quad w=0\quad \text{at}\ t=0. \end{equation}

The two patterns can be understood physically on the basis of the fact that the steady-drop solutions exist only for drops surrounded by oversaturated metastable vapour, which is stable with respect to infinitesimally small perturbations, but may be unstable with respect to finite ones. One can thus assume that the drop solution provides a lower bound for the mass of perturbations capable of triggering off instability of the surrounding vapour.

It is also interesting to see how a drop would evolve if it is surrounded by saturated (not oversaturated) vapour.

The fact that the thickness of a steady drop becomes infinite as $\rho _{\infty,1}\rightarrow \rho _{1}^{(v)}$ suggests that a finite-size liquid layer, surrounded by saturated vapour, evaporates. The mechanism of this evaporation, however, is not immediately clear. Two-dimensional and three-dimensional drops, for example, evaporate because the curvature of their boundary gives rise to a liquid-to-vapour mass flux (Benilov Reference Benilov2020d, Reference Benilov2021, Reference Benilov2022b), but the boundaries of 1-D drops are flat. One can only assume that they evaporate due to a long-range interaction of the drop's upper and lower interfaces.

To explore this hypothesis, 1-D drops floating in saturated vapour were simulated numerically (using the same model and numerical method as before). It turned out that these drops do evaporate, and their evolution can be subdivided into the following three distinct stages.

1. (i) The boundaries of the drop rapidly assume the profile of a steady liquid/vapour interface, of the kind examined in the previous subsection.

2. (ii) The drop begins to get thinner, but the density of the drop's core remains close to $\rho _{1}^{(l)}$.

3. (iii) Once the thickness of the drop becomes comparable to the interfacial thickness, the density at the drop's centre begins to rapidly decrease, and the drop disappears.

The main characteristic of such behaviour is the evaporation time $t_{e}$, which can be defined as the interval over which the density at the drop's centre falls by a factor of, say, $10$, i.e.

(6.13)\begin{equation} \frac{\rho_{1}(0,t_{e})}{\rho_{1}(0,0)}=0.1. \end{equation}

Here $t_{e}$ was computed as a function of the drop's initial size $W$ for the following initial condition:

(6.14)\begin{equation} \rho_{1}=\rho_{1}^{(v)}+\frac{\rho_{1}^{(l)}-\rho_{1}^{(v)}}{1+\left( \dfrac{2z}{W}\right)^{6}},\quad w=0\quad \text{at}\ t=0. \end{equation}

It turned out that, for $W$ changing from $0.2$ to $10$, the evaporation time $t_{e}$ grows from $48.5$ to $10^{6}$, i.e. exponentially. This agrees with the hypothesis that the drop evaporation is caused by long-range interaction of the drop's boundaries and the fact that the density in liquid/vapour interfaces tends to $\rho _{1}^{(v)}$ and $\rho _{1}^{(l)}$ (as $z\rightarrow \pm \infty$) exponentially quickly.

6.3. Solid/fluid interfaces

Let the fluid be bounded below by a flat horizontal wall (substrate) located at $z=0$, so that the boundary condition (2.47) reduces to

(6.15)\begin{equation} \rho_{i}=\rho_{0,i}\quad \text{at}\ z=0. \end{equation}

Far above the substrate, the fluid is homogeneous,

(6.16)\begin{equation} \rho_{i}\rightarrow\rho_{\infty,i}\quad \text{as}\ z\rightarrow+\infty, \end{equation}

where $\rho _{\infty,i}$ corresponds to a stable or metastable state. Finally, (6.9) used in the previous subsection for 1-D drops and bubbles applies to the present case as well.

The boundary-value problem (6.9), (6.15) and (6.16) was solved numerically for a van der Waals pure fluid with

(6.17a,b)\begin{equation} T=0.2,\quad \rho_{0,1}=0.3. \end{equation}

Its typical solutions are illustrated in figure 6 ($\rho _{0,i}$ varies, $\rho _{\infty,1}$ is fixed) and in figure 7 ($\rho _{0,i}$ is fixed, $\rho _{\infty,1}$ varies). The following features can be observed.

1. (i) For a given $\rho _{\infty,1}$, there exists a certain pool of $\rho _{0,i}$ that can be ‘connected’ to this $\rho _{\infty,1}$. If, for example, $\rho _{\infty,1}$ equals the saturated vapour density (the case illustrated in figure 6a), then $\rho _{0,1}$ must be smaller than the saturated liquid density (and vice versa: if $\rho _{\infty,1}=\rho _{1}^{(l)}$, then $\rho _{0,i}>\rho _{1}^{(v)}$; see figure 6b). When $\rho _{\infty,1}$ approaches the pool's boundary, a clearly visible liquid/vapour interface emerges in the solution and gradually moves away from the substrate.

2. (ii) It follows from above that, for some pairs $( \rho _{0,1} , \rho _{\infty,1})$, no solution exists. As illustrated in figure 7, some such pairs involve unstable values of $\rho _{\infty,1}$ (and, thus, are unimportant), but there are also ones with a metastable $\rho _{\infty,1}$.

   

   Figure 6. Examples of solid/fluid interfaces for a pure van der Waals fluid with parameters (6.17a,b), for a given $\rho _{\infty,1}$: (a) $\rho _{\infty,1}=\rho _{1}^{(v)}$, $\rho _{0,1}$ varies from $0$ to $\rho _{1}^{(l)}$; (b) $\rho _{\infty,1} =\rho _{1}^{(l)}$, $\rho _{0,1}$ varies from $1$ to $\rho _{1}^{(v)}$. The arrows show the direction towards the boundary of the pool of allowable values of $\rho _{0,1}$.

   

   Figure 7. Examples of solid/fluid interfaces for a pure van der Waals fluid with parameters (6.17a,b). The fluid at infinity is: (1) undersaturated vapour, (2) oversaturated vapour with two possible solutions (the unstable one is shown as a thin line), (3) oversaturated vapour with a single solution, (4) undersaturated liquid with a single solution, etc. Observe that there are no solutions between curves (3) and (4).

   To clarify what happens in such cases, the evolution was simulated numerically (using the model and method used in the previous two subsections). It turned out that, if the boundary-value problem (6.9), (6.15) and (6.16) does not have a solution, the substrate triggers off a spontaneous phase change. This result can be readily interpreted physically: if a sufficiently hydrophilic substrate (with a sufficiently small $\rho _{0,1}$) touches oversaturated vapour, it triggers off spontaneous condensation, whereas a sufficiently hydrophobic substrate touching undersaturated liquid triggers off evaporation. In either case, a shock wave of phase change propagates away from the substrate, and the solid/fluid interface cannot be stationary.

3. (iii) For some pairs $( \rho _{0,1}, \rho _{\infty,1})$, there are two different solutions: a monotonic and non-monotonic one (see figure 7). It is unclear which of the two occurs in reality.

   It turned out that only monotonic solutions of the boundary-value problem (6.9), (6.15) and (6.16) maximise the entropy, whereas non-monotonic solutions do not (the former claim is proved in the general case, but the latter, only for $N=1$; see Appendix C.3). That is, non-monotonic solutions cannot be ruled out with certainty for $N\geq 2$ – yet the mere fact that the general stability proof that works for monotonic solutions cannot be extended to non-monotonic ones seems to resolve the dilemma in favour of the former.

4. (iv) Another physically important conclusion follows from the fact that the boundary-value problem (6.9), (6.15) and (6.16) has no more than one stable solution, and, thus, does not admit solutions describing a liquid layer of a finite thickness on a substrate, with saturated vapour above it. All such layers, regardless of their thickness, evaporate, and the Kelvin effect cannot be responsible for this effect (because the liquid/vapour interface is flat). The evaporation in this case can only be caused by long-range interaction between the interface and the substrate.

5. (v) It is interesting to compare the interfacial profiles shown in figures 6 and 7 to those computed by Evans, Stewart & Wilding (Reference Evans, Stewart and Wilding2017) via the density functional theory and Monte Carlo method for a Lennard-Jones fluid bounded by a single wall or contained between two parallel walls. The single-wall profiles of Evans et al. (Reference Evans, Stewart and Wilding2017) (see their figure 3) are qualitatively similar to those computed in this paper, but their two-wall profiles are riddled with short-scale oscillations (see figures 6, 14 and 20 of Evans et al. Reference Evans, Stewart and Wilding2017). It is not clear whether the oscillations are caused by ‘interference’ of the walls: on the one hand, the distance between the walls exceeds the spatial scale of wall–molecule interaction by a factor of $30$ (hence, the ‘interference’ should be weak), but, on the other hand, no oscillations occur when this distance is infinite.

   Whatever the nature of the oscillations is, one should not expect the DIM to describe them due to the omission of the higher-order derivatives of $\rho _{i}$ when obtaining expression (2.43) for the van der Waals force.

8.1. The van der Waals parameter of pure water, nitrogen and oxygen

For a pure fluid, the EV expression (2.14) for the internal energy yields (the subscript ${1}$ is omitted)

(8.1)\begin{equation} e=cT-a\rho, \end{equation}

which suggests that the van der Waals parameter $a$ can be determined as the linear fit of the empiric dependence of $cT-e$ on $\rho$. The heat capacity $c$ in this expression should be the same as that in the EV kinetic theory, i.e. $3R$ for water and $5R/2$ for nitrogen and oxygen. For simplicity, the fitting was carried out using only the data on the critical isobar $p=p_{cr}$ (as done by Benilov & Benilov Reference Benilov and Benilov2019), but the resulting straight line fits the isobars $p=p_{cr}/2$ and $p=2p_{cr}$ reasonably well too (see the top panels of figure 9). The parameters of the critical points for the fluids under consideration, as well as the other parameters needed here and hereinafter, can be found in tables 1 and 2.

Figure 9. Thermodynamic properties of $\mathrm {H}_{2}\mathrm {O}$ (left-hand panels) and $\mathrm {N}_{2}$ (right-hand panels). The non-connected symbols show the empiric data from Lindstrom & Mallard (Reference Lindstrom and Mallard1997) presented in isobaric form, for three values of the pressure $p$ relative to its critical value, $p_{cr}$ (see the legend). The gap between $\rho ^{(v)}$ and $\rho ^{(l)}$ in the empiric data for $p=p_{cr}/2$ reflects the impossibility (difficulty) of measurements in an unstable (metastable) fluid. The dotted lines show the parameters calculated via the EV fluid model.

Table 1. The parameters of the triple and critical points (subscripts $_{tr}$ and $_{cr}$, respectively) of $\mathrm {H}_{2}\mathrm {O}$, $\mathrm {N}_{2}$ and $\mathrm {O}_{2}$ (Lindstrom & Mallard Reference Lindstrom and Mallard1997).

Table 2. The parameters of $\mathrm {H}_{2}\mathrm {O}$, $\mathrm {N}_{2}$, $\mathrm {O}_{2}$ and air: $m$ is the molecular mass, $R$ is the specific gas constant, $a$ is the van der Waals parameter, $K$ is the Korteweg constant. The parameters of air are calculated as the 79-21 weighted averages of the corresponding parameters of nitrogen and oxygen, respectively.

The resulting values of the van der Waals parameter $a$ for $\mathrm {H} _{2}\mathrm {O}$, $\mathrm {N}_{2}$ and $\mathrm {O}_{2}$ are presented in table 2. Interestingly, $a$ of water exceeds significantly those of nitrogen and oxygen. This is likely to be caused by the difference in the molecular structure of these fluids: the water molecule has a non-zero dipolar moment, whereas nitrogen and oxygen molecules are symmetric – hence, do not. Since the van der Waals force is of an electric nature, one can assume that dipolar molecules interact stronger than neutral ones.

8.2. The equation of state of pure $\mathrm {H}_{2}\mathrm {O}$, $\mathrm {N}_{2}$, and $\mathrm {O}_{2}$

Pure fluids will be described by the EV model (2.14) and (2.15) with

(8.2)\begin{equation} \varTheta=R\left[{-}q^{(0)}\ln\left( 1-0.99\frac{\rho}{\rho_{tp}}\right) +\sum_{n=1}^{4}q^{(n)}\left( \frac{\rho}{\rho_{tp}}\right)^{n}\right] , \end{equation}

where $q^{(0)}\ldots q^{(4)}$ are undetermined coefficients and $\rho _{tp}$ is the fluid's density at the triple point ($\rho _{tp}$ is just a convenient density scale; the fact that, at the triple point, the liquid and vapour are in equilibrium with the solid phase is irrelevant, as solids are not described by the DIM). Note that the first term in expression (8.2) sets the maximum density at $\rho _{tp}/0.99$.

The coefficients $q^{(n)}$ were determined for $\mathrm {H}_{2}\mathrm {O}$, $\mathrm {N}_{2}$ and $\mathrm {O}_{2}$ by ensuring that the expressions for $p(\rho,T)$ and $G(\rho,T)$ corresponding to (8.2) yield the ‘correct’ – i.e. empiric – values for the critical density, temperature and pressure, as well as the liquid and vapour densities at the triple-point temperature $T=T_{tp}$ (five equations for the five unknown coefficients). All the necessary empiric data can be found in tables 1 and 2, and the computed values of $q^{(n)}$ are listed in table 3. Such an approach to calibrating the EV fluid model is a slight modification of that of Benilov (Reference Benilov2020b), which is, in turn, a modification of the approach of Benilov & Benilov (Reference Benilov and Benilov2018).

Table 3. The coefficients of the equation of state (8.2) for $\mathrm {H}_{2}\mathrm {O}$, $\mathrm {N}_{2}$ and $\mathrm {O}_{2}$.

To illustrate how well the EV model describes real fluids, the isobaric (with $p$ fixed) dependence of the temperature on the density is plotted for $\mathrm {H}_{2}\mathrm {O}$ and $\mathrm {N}_{2}$ in the lower panels of figure 9. The results for $\mathrm {O}_{2}$ (not presented) are similar to those for $\mathrm {N}_{2}$.

Note that the DIM has been coupled with realistic equations of state before (e.g. Caupin Reference Caupin2005; Magaletti, Gallo & Casciola Reference Magaletti, Gallo and Casciola2021); for pure water, this was typically done using the IAPWS-95 equation (Wagner & Pruß Reference Wagner and Pruß2002). The EV model used here is undoubtedly less accurate than the IAPWS-95 model, but it allows one to describe consistently all of the fluids involved ($\mathrm {H} _{2}\mathrm {O}$, $\mathrm {N}_{2}$ and $\mathrm {O}_{2}$).

8.3. The Korteweg parameter of pure $\mathrm {H}_{2}\mathrm {O}$, $\mathrm {N}_{2}$, and $\mathrm {O}_{2}$

Benilov (Reference Benilov2020b) proposed to deduce $K$ from the requirement that the DIM predicts the correct value of the surface tension $\sigma ^{(l/v)}$ of the liquid/vapour interface at the triple point. The same was done in the present work: $\sigma ^{(l/v)}$ was calculated via the DIM formula (7.12) with $G$ and $p$ of the EV fluid, and the value of $K$ was chosen for which (7.12) yields the same result as the empiric formula of Somayajulu (Reference Somayajulu1988). The resulting $K$ for $\mathrm {H}_{2}\mathrm {O}$, $\mathrm {N}_{2}$ and $\mathrm {O}_{2}$ can be found in table 2.

It should be emphasised that the values of $K$ computed via the above approach depend on the chosen fluid model. This explains the difference between the result for $K$ in this paper and the one computed by Benilov (Reference Benilov2020b): the former is based on the EV fluid model and the latter on the van der Waals model. The resulting 30 % difference in $\sigma ^{(l/v)}$ reflects the fact that the latter model is much less accurate.

As a test of the whole approach based on the DIM and EV models, the theoretical values for the saturation characteristics and $\sigma ^{(l/v)}$ have been compared with their empiric counterparts for the whole temperature range where liquid/vapour interfaces exist, i.e. between the triple and critical points. The results are shown in figure 10: evidently, the theoretical predictions are reasonably accurate.

Figure 10. Comparison of the results obtained via the DIM for an EV fluid (solid line) with the corresponding empiric data (dashed line): (a) the densities of saturated vapour and liquid, (b) saturated pressure, (c) surface tension of liquid/vapour interface. The empiric data in panels (a) and (b) are from Lindstrom & Mallard (Reference Lindstrom and Mallard1997), and those in panel (c) are from Somayajulu (Reference Somayajulu1988).

8.4. Parameters of water/air interaction

Generally, the parameters of a fluid should be deduced from measurements of its equation of state, surface tension, etc., but these are rarely known for multicomponent fluids. The water/air mixture at normal conditions is an exception in this respect: the density of air is small in this case, and its equation of state is close to that of an ideal gas. In addition, air will be treated as a pure fluid with parameters equal to the weighted averages of those of nitrogen and oxygen (see table 2).

Thus, let the function $\varTheta$ (the non-ideal part of the entropy of an EV fluid) be independent of the air density $\rho _{2}$, so that expressions (2.17) and (2.18) with $N=2$ yield

(8.3)\begin{equation} p=T\left[ R_{1}\rho_{1}+R_{2}\rho_{2}+( \rho_{1}+\rho_{2}) \frac{\mathrm{d}\varTheta(\rho_{1})}{\mathrm{d}\rho_{1}}\right] -a_{11}\rho _{1}^{2}-2a_{12}\rho_{1}\rho_{2}, \end{equation}

(8.4)\begin{align} G_{1}&=T\left[ R_{1}\ln\rho_{1}+\varTheta(\rho_{1})+( \rho_{1}+\rho _{2})\frac{\mathrm{d}\varTheta(\rho_{1})}{\mathrm{d}\rho_{1}}\right]\nonumber\\ &\quad -2(a_{11}\rho_{1}+a_{12}\rho_{2}) +T[ R_{1}+c_{1}( 1-\ln T) ] , \end{align}

(8.5)\begin{equation} G_{2}=T[ R_{2}\ln\rho_{2}+\varTheta(\rho_{1})] -2a_{12}\rho _{1}+T[ R_{2}+c_{2}( 1-\ln T) ] . \end{equation}

Here, the function $\varTheta (\rho _{1})$ is the one given by (8.2), with $\rho$ changed to $\rho _{1}$, and with its coefficients corresponding to water.

Before using expressions (8.3)–(8.5), one needs to fix the non-diagonal term $a_{12}$ of the matrix $a_{ij}$, responsible for the interaction of water and air molecules. It can be deduced from $\rho _{2} ^{(l)}$ (the amount of air dissolved in water): one should choose such $a_{12}$ that the Maxwell construction based on (8.3)–(8.5) predicts the correct value of $\rho _{2}^{(l)}$. Since $a_{12}$ is supposed to not depend on $T$, such a calculation should be done for a single temperature and the atmospheric pressure. For $T=25\,^{\circ }\mathrm {C}$ and $p=1\ \mathrm {atm}$, for example, $\rho _{2}^{(l)}=0.0227\ \mathrm {kg\ m}^{-3}$ (The Engineering Toolbox 2004). It can then be deduced that the Maxwell construction based on (8.3)–(8.5) yields the correct value of $\rho _{2}^{(l)}$ if

(8.6)\begin{equation} a_{12}=208\ \mathrm{m}^{5}\ \mathrm{s}^{{-}2}\ \mathrm{kg}^{{-}1}, \end{equation}

i.e. approximately equal to $a_{22}$ and ten times smaller than $a_{11}$.

The accuracy of expressions (8.3)–(8.5) can be illustrated by the corresponding value of the boiling point ($T\approx 109\,^{\circ }\mathrm {C}$, as opposed to the exact value of $T=100\,^{\circ }\mathrm {C}$). Furthermore, at ‘room temperature’ (say, $T=25\,^{\circ }\mathrm {C}$), the error should be even smaller, because the room temperature is close to the triple point of water ($T\approx 0\,^{\circ }\mathrm {C}$) where (8.3)–(8.5) were calibrated.

It still remains to determine the non-diagonal term $K_{12}$ of the Korteweg matrix $K_{ij}$.

Since $K_{ij}$ is supposed to be positive definite, $K_{12}$ should satisfy

(8.7)\begin{equation} -1<\frac{K_{12}}{( K_{11}K_{22})^{1/2}}<1. \end{equation}

One would think that $K_{12}$ could be determined by comparing the surface tension of liquid-water/air interface to that of liquid-water/vapour-water interface. It turns out, however, that the difference between these is so small that the existing experimental techniques cannot detect it (see Vargaftik, Volkov & Voljak (Reference Vargaftik, Volkov and Voljak1983) and references therein). This can be due to smallness of $K_{12}$, but more likely, because the density of air is small. Either way, the determination of $K_{12}$ via measurements of surface tension is impossible.

To at least confirm that air cannot affect significantly the surface tension of water, the boundary-value problem (6.1)–(6.4) was solved for the parameters of water and air at $T=25\,^{\circ }\mathrm {C}$, and $K_{12}$ varying through range (8.7). It turned out that the resulting variation of the surface tension is only $7\,\%$.

Even though the dependence of the surface tension (which is a global characteristic) on $K_{12}$ is weak, this parameter does influence the local structure of the interface.

Figure 11 depicts the density profiles $\rho _{1}(z)$ and $\rho _{2}(z)$ for the water/air interface at room temperature (i.e. the above-mentioned solution of the boundary-value problem (6.1)–(6.4)). One can see that, as $K_{12}$ approaches the right endpoint of range (8.7), the profile of the water density is getting steeper, and a layer of high air density is developing inside the interface. For curve 5 of figure 11(b), which is half-way through range (8.7), the maximum of $\rho _{2}(z)$ exceeds the atmospheric air density by a factor of approximately $20$. Since the amount of air drawn into the interface grows with $K_{12}$, one concludes that the air is ‘pulled in’ by the van der Waals attraction exerted by the bulk of the liquid.

Figure 11. The spatial structure of water/air interfaces for $T=25\,^{\circ }\mathrm {C}$ and $K_{12}=( K_{11}K_{22})^{1/2}\times n/10$, where $n$ is the curve number: (a) density of water; (b) density of air. The horizontal scales of panels (a) and (b) differ by a factor of 5; the latter and the horizontal scale of the shaded inset in panel (b) differ by a factor of 1500. The inset shows that the concentration of the gas dissolved in liquid is non-zero.

This suggests a possibility of deducing $K_{12}$ from empiric data on evaporation of drops (which depends on both global and local characteristics of the interface Benilov Reference Benilov2022b). Such an approach, however, requires an extension of the evaporation model of Benilov (Reference Benilov2022b) to multicomponent fluids, which is currently in progress.

Observe that the interfaces depicted in figure 11 are such that $\rho _{2}(z)$ is not monotonic – hence, the sufficient stability criterion from § 6.1 and Appendix C.1 does not apply. This does not, however, mean that this interface is unstable. It is, in fact, difficult to imagine that a microscopic non-monotonicity of the density of air dissolved in water could destabilise the whole interface, but this issue is still of interest theoretically and, thus, deserves further investigation.

8.5. The viscosity and thermal conductivity of water/air mixture

There is a large body of work on the viscosity of multicomponent fluids, with even the simplest theoretical models, e.g. that of Enskog–Chapman for dilute gases (e.g. Ferziger & Kaper Reference Ferziger and Kaper1972, yielding a fairly complicated dependence of $\mu _{s}$ and $\mu _{b}$ on $\rho _{i}$. Phenomenological results, on the other hand, tend to include many ad hoc parameters specific to the fluid under consideration (e.g. Davidson (Reference Davidson1993) and the references therein). Overall, the simplest option seems to be the expression for the shear viscosity proposed on a phenomenological basis by Hind, McLaughlin & Ubbelohde (Reference Hind, McLaughlin and Ubbelohde1960) and justified, under certain approximations, via statistical mechanics by Bearman & Jones (Reference Bearman and Jones1960),

(8.8)\begin{equation} \mu_{s}=\mu_{s,1}f_{1}^{2}+( \mu_{s,1}+\mu_{s,2}) f_{1}f_{2} +\mu_{s,2}f_{2}^{2}, \end{equation}

where $\mu _{s,i}$ is the shear viscosity of the $i$th species and

(8.9)\begin{equation} f_{i}=\frac{\rho_{i}/m_{i}}{\rho_{1}/m_{1}+\rho_{2}/m_{2}} \end{equation}

is its mole fraction. Expression (8.8) does not include any fluid-specific parameters, but is capable of predicting the shear viscosity with a reasonable accuracy. If, for example, air is treated as a mixture of $\mathrm {N}_{2}$ and $\mathrm {O}_{2}$, the error of (8.8) for the range $T=0-100\,^{\circ }\mathrm {C}$ is less than 1 %. A similar ‘mixture rule’ can be assumed for the bulk viscosity and thermal conductivity,

(8.10)$$\begin{gather} \mu_{b}=\mu_{b,1}f_{1}^{2}+( \mu_{b,1}+\mu_{b,2}) f_{1}f_{2} +\mu_{b,2}f_{2}^{2}, \end{gather}$$

(8.11)$$\begin{gather}\kappa=\kappa_{1}f_{1}^{2}+( \kappa_{1}+\kappa_{2}) f_{1} f_{2}+\kappa_{2}f_{2}^{2}. \end{gather}$$

Generally, various mechanical properties of a multicomponent fluid are often described by the same mixture rule, in which case it is referred to as ‘generalized’.

It still remains to fix the dependence of $\mu _{s,i}$, $\mu _{b,i}$ and $\kappa _{i}$ on $( \rho _{1},\rho _{2},T)$. In application to air, which can be treated as a dilute gas, one can assume these parameter depend only on $T$ (as predicted by the kinetic theory of dilute gases), i.e.

(8.12a–c)\begin{equation} \mu_{s,2}(\rho_{2},T)=\mu_{s,2}(0,T),\quad \mu_{b,2}(\rho_{2},T)=\mu _{b,2}(0,T),\quad \kappa_{2}(\rho_{2},T)=\kappa_{2}(0,T), \end{equation}

where $\mu _{s,2}(0,T)$, $\mu _{b,2}(0,T)$ and $\kappa _{2}(0,T)$ are the small-density limiting values of the corresponding parameters.

For water, whose liquid phase cannot be treated as a dilute gas, such an approximation is inapplicable. Aiming again for simplicity, one can approximate both viscosities by a parabola passing through two reference points, the zero-density limit and the saturated liquid state,

(8.13)$$\begin{gather} \mu_{s,1}(\rho_{1},T)=\mu_{s,1}(0,T)+[ \mu_{s,1}(\rho_{1}^{(l)} ,T)-\mu_{s,1}(0,T)] \left[ \frac{\rho_{1}}{\rho_{1}^{(l)}(T)}\right] ^{2}, \end{gather}$$

(8.14)$$\begin{gather}\mu_{b,1}(\rho_{1},T)=\mu_{b,1}(0,T)+[ \mu_{b,1}(\rho_{1}^{(l)} ,T)-\mu_{b,1}(0,T)] \left[ \frac{\rho_{1}}{\rho_{1}^{(l)}(T)}\right] ^{2}. \end{gather}$$

A similar approximation, but by a linear function, will be adopted for the thermal conductivity,

(8.15)\begin{equation} \kappa_{1}(\rho_{1},T)=\kappa_{1}(0,T)+[ \kappa_{1}(\rho_{1} ^{(l)},T)-\kappa_{1}(0,T)] \frac{\rho_{1}}{\rho_{1}^{(l)}(T)}.\end{equation}

The accuracy of the above approximations for $\mu _{s,1}$ and $\kappa _{1}$ are illustrated in figure 12.

Figure 12. The parameters of pure water at $T=327\,^{\circ }\mathrm {C}$: (a) shear viscosity, (b) thermal conductivity. The empiric data (Lindstrom & Mallard Reference Lindstrom and Mallard1997) are shown as a solid line and approximations (8.13) and (8.15) as a dashed line. The dots mark the reference points. The gap in the empiric data reflects the impossibility (difficulty) of measurements in an unstable (metastable) fluid.

To use (8.8)–(8.15), one needs to know how the viscosities and thermal conductivities of water and air depend on $T$. For $\mu _{s}$ and $\kappa$, such data are widely available (e.g. Lindstrom & Mallard Reference Lindstrom and Mallard1997; White Reference White2005). Measurements of $\mu _{b}$, on the other hand, are scarce, but can still be found in Holmes, Parker & Povey (Reference Holmes, Parker and Povey2011), Cramer (Reference Cramer2012) and Shang et al. (Reference Shang, Wu, Wang, Yang, Ye, Hu, Tao and He2019) for liquid water, vapour water and air, respectively.

The empiric formulae proposed in these papers will not be discussed in detail. The characteristic values they yield for normal conditions are listed in table 4.

Table 4. The empiric viscosity and thermal conductivity of liquid water, water vapour and air at $25\,^{\circ }\mathrm {C}$, and the corresponding references.

8.6. The transport coefficients

Recalling restrictions (2.32a,b), one can deduce that the extended transport matrix for a two-component fluid is

(8.16)\begin{equation} D_{ij}^{(ext)}= \begin{bmatrix} D_{11} & -D_{11} & \zeta_{1}\\ -D_{11} & D_{11} & -\zeta_{1}\\ \zeta_{1} & -\zeta_{1} & \kappa \end{bmatrix}. \end{equation}

Evidently, it involves only three independent coefficients, one of which (the thermal conductivity $\kappa$) has already been discussed in § 8.5. The other two, $D_{11}$ and $\zeta _{1}$, are discussed below.

It can be safely assumed that diffusion is of importance only where the water density is comparable to that of air. In the region where the former is high, the latter is miniscule, and so is its influence on the global dynamics. This effectively means that the diffusivities can be represented using the Enskog–Chapman method; even though it does not apply to liquid water, the error due to using it anyway is negligible.

According to the leading-order Enskog–Chapman formula (e.g. Ferziger & Kaper Reference Ferziger and Kaper1972),

(8.17)\begin{equation} D_{11}=\frac{\rho_{1}\rho_{2}( \rho_{1}m_{2}+\rho_{2}m_{1}) }{( \rho_{1}+\rho_{2})^{2}k_{B}T}\mathfrak{D}(T), \end{equation}

where $\mathfrak {D}(T)$ does not depend on $\rho _{1}$ and $\rho _{2}$. Its dependence on $T$ can be found in The Engineering Toolbox (2018), from which the following characteristic value can be deduced,

(8.18)\begin{equation} \mathfrak{D}(25\,^{{\circ}}\mathrm{C})=2.49\times10^{{-}5}\ \mathrm{m}^{2} \ \mathrm{s}^{{-}1}. \end{equation}

The Enskog–Chapman expression for the thermodiffusivity is fairly bulky; as a result, one often uses the simpler formula of de Groot & Mazur (Reference de Groot and Mazur1962),

(8.19)\begin{equation} \zeta_{1}=( \rho_{1}+\rho_{2}) \frac{( \rho_{1} /m_{1}) ( \rho_{2}/m_{2}) }{( \rho_{1}/m_{1} +\rho_{2}/m_{2})^{2}}\mathfrak{U}(T). \end{equation}

Unfortunately, no data on $\mathfrak {U}(T)$ for water/air mixture are available in the literature; the author of this paper was able to find only an estimate for a single temperature value (Lidon, Perrot & Stroock Reference Lidon, Perrot and Stroock2021),

(8.20)\begin{equation} \mathfrak{U}(21\,^{{\circ}}\mathrm{C})={-}4.98\times10^{{-}6}\ \mathrm{m}^{2} \mathrm{s}^{{-}1}. \end{equation}

It is shown later, however, that thermodiffusion in water/air interfaces is negligible, so the lack of data on $\mathfrak {U}(T)$ is unimportant.

In the next subsection characteristic values of the coefficients $D_{11}$ and $\zeta _{1}$ will be needed. These will be estimated for the particular case $\rho _{1}=\rho _{2}=1.17\ \mathrm {kg\ m}^{-3}$, i.e. when the water density matches that of air at normal conditions. Substituting these values and the molecular masses of water and air from table 2 into (8.17)–(8.20), one obtains

(8.21a,b)\begin{equation} D_{11}\approx1.37\times10^{{-}10}\ \mathrm{m}^{3}\ \mathrm{s\ kg},\quad \zeta _{1}\approx{-}2.76\times10^{{-}6}\ \mathrm{m}^{{-}1}\ \mathrm{s}^{{-}1}\ \mathrm{kg}. \end{equation}

Even though these values apply to different temperatures ($T=25\,^{\circ }\mathrm {C}$ and $T=21\,^{\circ }\mathrm {C}$), they will be used simultaneously in the same qualitative estimate (under the implied assumption that the four-degree difference does not matter).

8.7. The non-dimensional parameters

To understand which effects are important for water/air interfaces under normal conditions, one should estimate the non-dimensional parameters (4.13a–e). They will be calculated using the characteristics of water from tables 2–4 and estimates (8.21a,b). The characteristic pressure scale will be assumed to be $P=a\varrho ^{2}$, where $\varrho =997\ \mathrm {kg\ m}^{-3}$ is the density of liquid water at $25\,^{\circ }\mathrm {C}$ and $a$ equals the van der Waals parameter of water from table 2. The following expression for the viscosity scale will be used:

(8.22)\begin{equation} \mu=\tfrac{4}{3}\mu_{s}+\mu_{b} \end{equation}

that arises naturally in problems involving liquid films (Benilov Reference Benilov2020c, Reference Benilov2022a) and drops (Benilov Reference Benilov2022b). In the estimates for this paper, the viscosities $\mu _{s}$ and $\mu _{b}$ are those for water at $25\,^{\circ }\mathrm {C}$.

The following values of parameters (4.13a–e) have been obtained:

(8.23a–e)\begin{equation} \alpha\approx0.00139,\quad \tau\approx0.0653,\quad \beta\approx0.0593,\quad \nu\approx0.0321,\quad \delta\approx0.0176. \end{equation}

The smallness of the microscopic Reynolds number $\alpha$ suggests that inertia plays little role in the dynamics of water/air interfaces (hence, one can use the Stokes approximation). The smallness of the non-dimensional temperature $\tau$ does not have physical implications, but enables one to use asymptotic tools when calculating, say, the profile of the equilibrium interface (Benilov Reference Benilov2020c). The smallness of the Nusselt number $\nu$ implies that the Soret and Dufour effects are negligible, and so the lack of empiric data on them for water and air is unimportant. The smallness of $\beta$, in turn, implies that the flow is almost isothermal, whereas the small value of $\delta$ suggests that diffusion dominates advection.

Similar estimates have also been carried out for the parameters of air at $p=1\ \mathrm {atm}$ and $T=25\,^{\circ }\mathrm {C}$. It turned out that $\alpha$, $\beta$, $\nu$ and $\delta$ are even smaller than those for water, but the non-dimensional temperature $\tau$ for air is order one. The latter is not surprising, as the room temperature is much higher than the freezing temperatures of nitrogen and oxygen, but is close to that for water.

All these observations should be helpful when using the DIM to examine asymptotically the problems of evaporation of drops and moving contact lines.

One should keep in mind, however, that parameters (4.13a–e) are sensitive to the choice of the viscosity scale $\mu$. Benilov (Reference Benilov2020a), for example, chose $\mu$ equal to the half-sum of the shear viscosities of liquid water and vapour water; as a result, $\alpha$ and $\beta$ were noticeably larger than those calculated above. More generally, one should choose the viscosity scale, as well as the other characteristic scales, to essentially reflect the physics of the setting at hand.

Note also that estimates (8.23a–e) apply to water under normal conditions, but in industrial applications the governing non-dimensional parameters can be very different. In steam turbines, for example, the temperature can be as high as $540\,^{\circ }\mathrm {C}$, and the pressure can exceed $230\ \mathrm {atm}$ (Vasserman & Shutenko Reference Vasserman and Shutenko2017).