2.1. Second-gradient theory: Korteweg's stress and capillary fluid equations

In the second-gradient theory as introduced by van der Waals (Reference van der Waals1894), free energy per unit volume is expressed as

(2.1)\begin{equation} \mathcal{A}_{vdW} = \mathcal{A} + \tfrac{1}{2}\kappa {\lvert \boldsymbol{\nabla}\rho\lvert}^2, \end{equation}

where $\mathcal {A}$ is the bulk free energy per unit volume, $\rho$ is the density and $\kappa$ is the capillary coefficient. The second term represents the interface energy while the bulk free energy is solely a function of the local density and temperature (Cahn & Hilliard Reference Cahn and Hilliard1958; Giovangigli Reference Giovangigli2020). The equilibrium state of the corresponding system is obtained by minimizing the free energy in a given volume under the constraint of constant total mass, leading to the stress tensor (Anderson et al. Reference Anderson, McFadden and Wheeler1998),

(2.2)\begin{equation} \boldsymbol{\mathsf{T}}_{K} = \boldsymbol{\nabla}\otimes\frac{\partial\mathcal{L}}{\partial(\boldsymbol{\nabla}\rho)}- \mathcal{L}\boldsymbol{\mathsf{I}}, \end{equation}

where $\boldsymbol{\mathsf{I}}$ is unit tensor and $\mathcal {L}$ is the Lagrange function,

(2.3)\begin{equation} \mathcal{L} = \mathcal{A} + \tfrac{1}{2}\kappa {\lvert \boldsymbol{\nabla}\rho\lvert}^2 - \lambda\rho, \end{equation}

and $\lambda$ is the Lagrange multiplier for the mass constraint, or chemical potential,

(2.4)\begin{equation} \lambda=\frac{\partial\mathcal{A}}{\partial\rho}-\kappa\boldsymbol{\nabla}^2\rho, \end{equation}

where $\boldsymbol {\nabla }^2$ is the Laplace operator. This in turn leads to the following Korteweg's stress tensor (Korteweg Reference Korteweg1901):

(2.5)\begin{equation} \boldsymbol{\mathsf{T}}_{K} = (P-\kappa\rho\boldsymbol{\nabla}^2\rho - \tfrac{1}{2}\kappa {\lvert \boldsymbol{\nabla}\rho\lvert}^2 )\boldsymbol{\mathsf{I}} + \kappa \boldsymbol{\nabla}\rho\otimes\boldsymbol{\nabla}\rho, \end{equation}

where

(2.6)\begin{equation} P=\rho \frac{\partial\mathcal{A}}{\partial\rho}-\mathcal{A} \end{equation}

is the thermodynamic pressure, or equation of state. From the local balance equations for mass and momentum, one obtains the macroscopic governing laws for an isothermal capillary fluid:

(2.7)$$\begin{gather} \partial_t\rho + \boldsymbol{\nabla}\boldsymbol{\cdot}\rho \boldsymbol{u} = 0, \end{gather}$$

(2.8)$$\begin{gather}\partial_t\rho\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{\cdot}\rho \boldsymbol{u}\otimes\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}= 0, \end{gather}$$

where $\boldsymbol {u}$ is the fluid velocity and the stress tensor $\boldsymbol{\mathsf{T}}$ is

(2.9)\begin{equation} \boldsymbol{\mathsf{T}} = \boldsymbol{\mathsf{T}}_K + \boldsymbol{\mathsf{T}}_{NS}. \end{equation}

The Navier–Stokes viscous stress tensor reads

(2.10)\begin{equation} \boldsymbol{\mathsf{T}}_{NS}={-}\mu\boldsymbol{\mathsf{S}}-\eta (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u})\boldsymbol{\mathsf{I}}, \end{equation}

where $\boldsymbol{\mathsf{S}}$ is the trace-free rate-of-strain tensor,

(2.11)\begin{equation} \boldsymbol{\mathsf{S}}=\boldsymbol{\nabla}\boldsymbol{u} + {\boldsymbol{\nabla}\boldsymbol{u}}^{{{\dagger}}} -\tfrac{2}{3}(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u})\boldsymbol{\mathsf{I}}, \end{equation}

and $\mu$ and $\eta$ are the dynamic and the bulk viscosity, respectively.

The momentum balance equation in (2.8) can be recast in the following form:

(2.12)\begin{equation} \partial_t\rho\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{\cdot}\rho \boldsymbol{u}\otimes\boldsymbol{u} + \boldsymbol{F}_{K} +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{NS}=0, \end{equation}

where Korteweg's force $\boldsymbol {F}_{K}$ is the divergence of the Korteweg pressure tensor,

(2.13)\begin{equation} \boldsymbol{F}_{K}=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{K}. \end{equation}

The latter can be written in the following form:

(2.14)\begin{equation} \boldsymbol{F}_{K}=\boldsymbol{\nabla}P_0+\boldsymbol{\nabla}(P-P_0) -\kappa\rho\boldsymbol{\nabla}(\boldsymbol{\nabla}^2\rho), \end{equation}

where we have introduced a reference pressure $P_0$. Navier–Stokes momentum equations with Korteweg's force in (2.14) shall be a target for reconstruction by a suitable kinetic model.

2.2. Kinetic model for non-ideal fluid

To introduce a kinetic model for non-ideal fluid, we begin with the first Bogolioubov–Born–Green–Kirkwood–Yvon (BBGKY) equation,

(2.15)\begin{equation} \partial_t\, f + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f = \mathcal{J}=\iint \boldsymbol{\nabla} V(\lvert \boldsymbol{r}-\boldsymbol{r}_1\rvert) \boldsymbol{\cdot}\frac{\partial}{\partial \boldsymbol{v}}f_{2}(\boldsymbol{r},\boldsymbol{v}, \boldsymbol{r}_1,\boldsymbol{v}_1,t) \,{\rm d}\boldsymbol{v}_1\,{\rm d}\boldsymbol{r}_1, \end{equation}

where $f(\boldsymbol {r},\boldsymbol {v},t)$ and $f_{2}(\boldsymbol {r},\boldsymbol {v},\boldsymbol {r}_1,\boldsymbol {v}_1,t)$ are the one- and the two-particle distribution functions, respectively, $\boldsymbol {r}$, $\boldsymbol {r}_1$ and $\boldsymbol {v}$, $\boldsymbol {v}_1$ are particles’ position and velocity, while $V$ is a potential of a pair interaction. The local equilibrium state is defined by the Maxwellian $f^{eq}$ at constant temperature $T$, parametrized by the local values of density $\rho$ and flow velocity $\boldsymbol {u}$,

(2.16)\begin{equation} f^{eq}=\frac{\rho}{(2{\rm \pi} RT)^{3/2}}\exp\left[-\frac{(\boldsymbol{v}-\boldsymbol{u})^2}{2RT}\right], \end{equation}

where $R$ is the gas constant. Furthermore, let us introduce a projector $\mathcal {K}$ onto local equilibrium at constant temperature (Gorban & Karlin Reference Gorban and Karlin2005),

(2.17)\begin{equation} \mathcal{K}\mathcal{J}= \left(\dfrac{\partial f^{eq}}{\partial \rho}-\frac{1}{\rho}\boldsymbol{u}\boldsymbol{\cdot}\frac{\partial f^{eq}}{\partial \boldsymbol{u}}\right)\int \mathcal{J} \,{\rm d}\boldsymbol{v} +\frac{1}{\rho}\dfrac{\partial f^{eq}}{\partial \boldsymbol{u}}\boldsymbol{\cdot} \int \boldsymbol{v}\mathcal{J} \,{\rm d}\boldsymbol{v}. \end{equation}

The projector property, $\mathcal {K}^2=\mathcal {K}$, can be verified by a direct computation. With the projector in (2.17), the interaction term in (2.15) is split into two parts by writing an identity,

(2.18)\begin{equation} \mathcal{J} = (1-\mathcal{K})\mathcal{J} + \mathcal{K} \mathcal{J}. \end{equation}

The first term,

(2.19)\begin{equation} \mathcal{J}_{loc}=(1-\mathcal{K})\mathcal{J}, \end{equation}

satisfies the local conservation of both mass and momentum,

(2.20)\begin{equation} \mathcal{K}\mathcal{J}_{loc}=0. \end{equation}

It is conventional to model the locally conserving part of the interaction with a single relaxation time Bhatnagar–Gross–Krook (BGK) approximation,

(2.21)\begin{equation} \mathcal{J}_{loc} \to \mathcal{J}_{BGK} ={-}\frac{1}{\tau}(1-\mathcal{K})f ={-}\frac{1}{\tau}(f- f^{eq}), \end{equation}

where the relaxation time $\tau$ is a free parameter. The second term in the identity of (2.18),

(2.22)\begin{equation} \mathcal{J}_{nloc}=\mathcal{K}\mathcal{J}, \end{equation}

satisfies the local mass but not the local momentum conservation. Indeed, after integration by part in the velocity $\boldsymbol {v}$ and neglecting boundary integrals, we arrive at

(2.23)\begin{equation} \mathcal{J}_{nloc} ={-}\frac{1}{\rho}\frac{\partial f^{eq}}{\partial \boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{F}_{nloc}, \end{equation}

where the force $\boldsymbol {F}_{nloc}$ reads

(2.24)\begin{equation} \boldsymbol{F}_{nloc}=\iiint \boldsymbol{\nabla}V(\lvert \boldsymbol{r}-\boldsymbol{r}_1\rvert)f_{2}(\boldsymbol{r},\boldsymbol{v}, \boldsymbol{r}_1,\boldsymbol{v}_1,t) \,{\rm d}\boldsymbol{v}_1\,{\rm d}\boldsymbol{r}_1\,{\rm d}\boldsymbol{v}. \end{equation}

Collecting the BGK approximation together with the non-local contribution, a generic kinetic model may be written as

(2.25)\begin{equation} \partial_t\, f + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f ={-}\frac{1}{\tau}(f - f^{eq}) - \frac{1}{\rho}\frac{\partial f^{eq}}{\partial \boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{F}_{nloc}. \end{equation}

Evaluation of the force in (2.24) requires us to specify the particles interaction. It is customary to invoke the Enskog–Vlasov model (Enskog Reference Enskog1921; Vlasov Reference Vlasov1961), where both hard-sphere collisions and a weak long-range attraction potential contribute to a non-local momentum transfer. For the hard-sphere Enskog part, a de-localization of the collision is responsible for a non-vanishing contribution of momentum transfer through the distance between the centres of the spheres upon their impact while the Vlasov approximation contributes non-locally to the momentum transfer from a distributed mean-field force. Evaluation of both the Enskog and Vlasov contributions to the force in (2.24) proceeds along familiar lines (Chapman & Cowling Reference Chapman and Cowling1939; He et al. Reference He, Shan and Doolen1998; He & Doolen Reference He and Doolen2002; Martys Reference Martys2006) and is reported in Appendices A and B for completeness,

(2.26)\begin{equation} \boldsymbol{F}_{EV} = \boldsymbol{F}_{E}+\boldsymbol{F}_{V}. \end{equation}

The first term is the lowest-order contribution of the collisional momentum transfer from the Enskog hard-sphere model,

(2.27)\begin{equation} \boldsymbol{F}_{E} = \boldsymbol{\nabla}b\rho^2\chi RT+ O(\boldsymbol{\nabla}^3\rho). \end{equation}

Here $b=4v_{HS}$, with $v_{HS}=\mathcal {V}_{HS}/m$ the specific volume of a hard sphere of diameter $d$ and mass $m$, while $\mathcal {V}_{HS}={\rm \pi} d^3/6$ is the volume of the sphere. Moreover, $\chi$ is the equilibrium two-particle correlation function, evaluated at the local density reduced by the specific volume of the hard sphere, $\chi =\chi (b\rho (\boldsymbol {r},t))$. To the lowest order, $\chi =1+(5/8)b\rho +O((b\rho )^2)$, cf. Chapman & Cowling (Reference Chapman and Cowling1939). The second term in (2.26) is the contribution of a long-range attraction potential $V$ in the mean field Vlasov approximation. To third order in the gradient of density,

(2.28)\begin{equation} \boldsymbol{F}_{V} ={-}\boldsymbol{\nabla}a\rho^2 - \kappa \rho\boldsymbol{\nabla}\boldsymbol{\nabla}^2\rho +O(\boldsymbol{\nabla}^5\rho), \end{equation}

where parameters $a$ and $\kappa$ are

(2.29)$$\begin{gather} a ={-}2{\rm \pi}\int_{d}^{\infty} r^2V(r) \,{\rm d}r, \end{gather}$$

(2.30)$$\begin{gather}\kappa={-}\frac{2{\rm \pi}}{3}\int_{d}^{\infty} r^4V(r) \,{\rm d}r. \end{gather}$$

Thus, with the approximations specified, the non-local force in (2.26) becomes

(2.31)\begin{equation} \boldsymbol{F}_{EV} = \boldsymbol{\nabla}(P_{EV}-P_0) - \kappa \rho\boldsymbol{\nabla}\boldsymbol{\nabla}^2\rho, \end{equation}

where the reference ideal gas pressure $P_0$ is provided by the local Maxwellian in (2.16),

(2.32)\begin{equation} P_0=\frac{1}{3}\int \lvert\boldsymbol{v}-\boldsymbol{u}\rvert^2 f^{eq}\,{\rm d}\boldsymbol{v}=\rho R T, \end{equation}

while the equation of state of the Enskog–Vlasov gas is of van der Waals type,

(2.33)\begin{equation} P_{EV}=\rho RT(1 + b\rho\chi) -a\rho^2. \end{equation}

This allows us to extend the Enskog–Vlasov kinetic model, and a phenomenological equation of state $P$ in (2.6) can be used instead of $P_{EV}$ in (2.33). Moreover, the reference pressure $P_0$ can be made selective by rescaling the local equilibrium,

(2.34)\begin{equation} f^{eq}=\frac{\rho}{(2{\rm \pi} P_0/\rho)^{D/2}}\exp\left[-\frac{(\boldsymbol{v}-\boldsymbol{u})^2}{2P_0/\rho}\right], \end{equation}

where $D$ is the space dimension. While the Enskog–Vlasov partition above corresponds to selecting $P_0 = \rho R T$, an alternative is provided by Reyhanian, Dorschner & Karlin (Reference Reyhanian, Dorschner and Karlin2020), where $P_0 = P$ is chosen. Another construction was proposed (Martys Reference Martys1999, Reference Martys2001) where all contributions to the pressure tensor are introduced with the local attractor, rescaling temperature and introducing a velocity shift to match the targeted system of moments. Using the rescaled equilibria in (2.34), a family of kinetic models parametrized by the reference pressure reads

(2.35)\begin{equation} \partial_t\, f + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f ={-}\frac{1}{\tau}(f - f^{eq}) - \frac{1}{\rho}\frac{\partial f^{eq}}{\partial \boldsymbol{u}}\boldsymbol{\cdot}[\boldsymbol{\nabla}(P-P_0)-\kappa \rho\boldsymbol{\nabla}\boldsymbol{\nabla}^2\rho]. \end{equation}

The kinetic equation in (2.35) shall be considered as a semi-phenomenological model of a non-ideal fluid, with the relaxation time $\tau$, the capillarity coefficient $\kappa$, pressure $P$ and reference pressure $P_0$ as phenomenological input parameters, while the Enskog–Vlasov realization will serve as a representative example for estimates of various flow regimes.

The analysis of the kinetic model in (2.35) under the conventional scaling of a small deviation from a uniform state (Chapman & Cowling Reference Chapman and Cowling1939),

(2.36a,b)\begin{equation} \boldsymbol{\nabla}\to \epsilon\boldsymbol{\nabla},\quad \partial_t\to\epsilon\partial_t, \end{equation}

is detailed in Appendix C. To second order in space derivatives, the resulting momentum balance equation reads

(2.37)\begin{equation} \partial_t\rho\boldsymbol{u} + \epsilon\boldsymbol{\nabla}\boldsymbol{\cdot}\rho \boldsymbol{u}\otimes\boldsymbol{u} + \epsilon\boldsymbol{\nabla}P +\epsilon\boldsymbol{\nabla}\boldsymbol{\cdot}\epsilon\boldsymbol{\mathsf{T}}_{NS}+O(\epsilon^3) =0, \end{equation}

where the dynamic viscosity $\mu$ and the bulk viscosity $\eta$ in the Navier–Stokes stress tensor in (2.10) are defined by the reference pressure ($D=3$),

(2.38)$$\begin{gather} \mu=\tau P_0, \end{gather}$$

(2.39)$$\begin{gather}\eta=\left(\frac{5}{3}-\frac{\partial\ln P_0}{\partial\ln\rho}\right)\tau P_0. \end{gather}$$

Thus, the momentum balance equation in (2.37) is form-invariant with respect to the choice of reference pressure, provided $P_0$ satisfies a sub-isentropic condition,

(2.40)\begin{equation} {P_0}\le C\rho^{5/3}, \end{equation}

for some $C>0$. With (2.40), the bulk viscosity in (2.39) is positive and vanishes when the reference pressure follows an isentropic process for ideal monatomic gas, $P_0=C\rho ^{5/3}$. For example, any polytropic process, $P_0=A\rho ^n$, $1\le n\le 5/3$ satisfies the sub-isentropic condition and results in $\eta =(5/3-n)\tau P_0$. A special case of the isothermal process $n=1$ returns $\eta =(2/3)\tau P_0$, and the viscous stress tensor becomes

(2.41)\begin{equation} \boldsymbol{\mathsf{T}}_{NS}={-}\tau P_0(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{\dagger}). \end{equation}

However, when compared to the two-phase momentum equation in (2.12), the macroscopic limit recovers only the non-ideal gas component thereof while missing Korteweg's capillarity contribution. Indeed, the third-order term, $\sim \epsilon ^3\rho \boldsymbol {\nabla }\boldsymbol {\nabla }^2\rho$ in (2.28) and (2.35), does not contribute to the momentum equation in (2.37) under the scaling of (2.36a,b). This is consistent with the well-known results from kinetic theory (Chapman & Cowling Reference Chapman and Cowling1939) and is not surprising. The scaling of (2.36a,b) is essentially based on the Knudsen number, which overrides the relative contribution of the capillarity term by two orders, cf. Appendix C. Thus, under the weak non-uniformity assumption in (2.36a,b), the capillarity terms are seen as higher-order Burnett-level contributions, and cannot appear in the main (first and second) orders in the momentum balance equation of (2.37). In fact, the condition in (2.36a,b) rules out situations at an interface between phases where gradients of density become large over a relatively short distance. Therefore, for the kinetic model in (2.35) to recover in-full the momentum balance in (2.12), a different scaling needs to be applied.

2.3. Scaling by velocity increment and lattice Boltzmann equation

2.3.1. Time step and forcing

A rescaling of the kinetic model in this section shall be maintained by introducing a time step $\delta t$. As a preliminary consideration, we evaluate the contribution of the force term over the time step. To that end, as noted by Kupershtokh, Medvedev & Karpov (Reference Kupershtokh, Medvedev and Karpov2009), for a generic force $\boldsymbol {F}$, we can write the action of the force on the distribution function as a full derivative in a frame that moves with the local fluid velocity,

(2.42)\begin{equation} \frac{1}{\rho}\frac{\partial f^{eq}}{\partial \boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{F} = \frac{{\rm d} f^{eq}}{{\rm d}t}. \end{equation}

Introducing the velocity increment,

(2.43)\begin{equation} \delta\boldsymbol{u}= \frac{\boldsymbol{F}}{\rho}\delta t, \end{equation}

and integrating in time, leads to an approximation,

(2.44)\begin{equation} \mathcal{F}=\int_{t}^{t+\delta t} \frac{1}{\rho}\frac{\partial f^{eq}}{\partial \boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{F} \,{\rm d}t \approx f^{eq}(\boldsymbol{u}+\delta\boldsymbol{u}) - f^{eq}(\boldsymbol{u}). \end{equation}

This so-called exact difference method (EDM) becomes accurate for a gravity force, $\boldsymbol {F}/\rho ={\rm const.}$, otherwise it often provides a reliable estimate for the force term and is widely used. In what follows, the scaling to be applied assumes smallness of the velocity increment in (2.43) rather than smoothness of the spatial distribution of the force. Since the velocity increment is based on a time step, it is natural to proceed with a lattice Boltzmann realization of the kinetic equation.

2.3.2. Standard lattice and product form

The lattice Boltzmann model shall be realized with the standard discrete velocity set $D3Q27$, where $D=3$ stands for three dimensions and $Q=27$ is the number of discrete velocities,

(2.45a,b)\begin{equation} \boldsymbol{c}_i=(c_{ix},c_{iy},c_{iz}),\quad c_{i\alpha}\in\{{-}1,0,1\}. \end{equation}

We first define a triplet of functions in two variables, $\xi _{\alpha }$ and $\zeta _{\alpha \alpha }$,

(2.46)$$\begin{gather} \varPsi_{0}(\xi_{\alpha},\zeta_{\alpha\alpha}) = 1 - \zeta_{\alpha\alpha}, \end{gather}$$

(2.47)$$\begin{gather}\varPsi_{1}(\xi_{\alpha},\zeta_{\alpha\alpha}) = \frac{\xi_{\alpha} + \zeta_{\alpha\alpha}}{2}, \end{gather}$$

(2.48)$$\begin{gather}\varPsi_{{-}1}(\xi_{\alpha},\zeta_{\alpha\alpha}) = \frac{-\xi_{\alpha} + \zeta_{\alpha\alpha}}{2}, \end{gather}$$

and consider a product form associated with the discrete velocities $\boldsymbol {c}_i$ in (2.45a,b),

(2.49)\begin{equation} \varPsi_i= \varPsi_{c_{ix}}(\xi_x,\zeta_{xx}) \varPsi_{c_{iy}}(\xi_y,\zeta_{yy}) \varPsi_{c_{iz}}(\xi_z,\zeta_{zz}). \end{equation}

All pertinent populations below shall be determined by specifying the parameters $\xi _\alpha$ and $\zeta _{\alpha \alpha }$ in the product form in (2.49). A two-dimensional version of the model on the $D2Q9$ lattice is obtained by omitting the $z$-component in all formulae below. Finally, we use the notation $\delta \boldsymbol {r}_{i}=\boldsymbol {c}_i\delta t$ for the lattice links, and denote the grid spacing in any direction $\alpha =x,y,z$ as $\delta r=\vert c_{i\alpha }\vert \delta t$, $c_{i\alpha }\ne 0$. For the $D3Q27$ discrete velocity set in (2.45a,b), the lattice spacing is the same in all Cartesian directions.

2.3.3. The lattice Boltzmann equation

The local density $\rho$ and flow velocity $\boldsymbol {u}$ are defined using the populations $f_i$,

(2.50)$$\begin{gather} \rho(\boldsymbol{r},t)=\sum_{i=0}^{Q-1}f_i(\boldsymbol{r},t), \end{gather}$$

(2.51)$$\begin{gather}\rho\boldsymbol{u}(\boldsymbol{r},t)=\sum_{i=0}^{Q-1}\boldsymbol{c}_i\, f_i(\boldsymbol{r},t). \end{gather}$$

With a reference pressure $P_0$, and by setting the parameters,

(2.52)$$\begin{gather} \xi_{\alpha}=u_{\alpha}, \end{gather}$$

(2.53)$$\begin{gather}\zeta_{\alpha\alpha}=\frac{P_0}{\rho}+u_{\alpha}^2, \end{gather}$$

the local equilibrium populations are represented with the product-form in (2.49),

(2.54)\begin{equation} f_i^{eq}= \rho\prod_{\alpha=x,y,z}\varPsi_{c_{i\alpha}}\left(u_\alpha,\frac{P_0}{\rho}+u_{\alpha}^2\right). \end{equation}

The LBGK equation, supplemented with a forcing term, is written as

(2.55)\begin{equation} f_i(\boldsymbol{r}+\boldsymbol{c}_i \delta t, t+\delta t) - f_i(\boldsymbol{r}, t) = \omega(f_i^{eq} - f_i) + (f_i^*-f_i^{eq}), \end{equation}

where $\omega$ is a dimensionless relaxation parameter, the equilibrium populations are provided by (2.54), while the extended equilibrium populations $f_i^*$ are defined by the product-form in (2.49) with the following assignment for the parameters:

(2.56)$$\begin{gather} \xi_{\alpha}^{*} = u_{\alpha}+\frac{F_{\alpha}\delta t}{\rho}, \end{gather}$$

(2.57)$$\begin{gather}\zeta_{\alpha\alpha}^{*} = \frac{P_0}{\rho} +u_{\alpha}^2 + \frac{\varPhi_{\alpha\alpha}}{\rho}, \end{gather}$$

where $\varPhi _{\alpha \alpha }/\rho$ is a correction term for the diagonals of the non-equilibrium momentum flux tensor,

(2.58)\begin{equation} \varPhi_{\alpha\alpha} = \left(1-\frac{\omega}{2}\right) \delta t\partial_{\alpha}\left(\rho u_{\alpha} \left(u_{\alpha}^2 + \frac{3P_{0}}{\rho}-3\varsigma^2\right)\right), \end{equation}

where $\varsigma = \delta r/\sqrt {3}\delta t$ is the so-called lattice speed of sound. The form of the correction term is a result of the multi-scale analysis presented in Appendix D. Thus, the extended equilibrium is explicitly written as

(2.59)\begin{equation} f_i^*=\rho\prod_{\alpha=x,y,z}\varPsi_{c_{i\alpha}}\left(u_\alpha+\frac{F_{\alpha}\delta t}{\rho},\frac{P_0}{\rho}+u_{\alpha}^2+\frac{\varPhi_{\alpha\alpha}}{\rho}\right). \end{equation}

Comments are in order.

• • If the correction term of (2.58) is omitted in (2.57), the population in (2.59) becomes the equilibrium with the increment $\delta \boldsymbol {u}$ in (2.43) due to force action added to the flow velocity $\boldsymbol {u}$. This corresponds to the EDM forcing maintained by the second bracket in the right-hand side of the LBGK equation in (2.55).

• • The correction term in (2.56) is added following a proposal by Saadat et al. (Reference Saadat, Hosseini, Dorschner and Karlin2021). Its purpose is to compensate for the bias of the $D3Q27$ velocities in (2.45a,b), $c_{i\alpha }^3=c_{i\alpha }$, and to restore the Galilean invariance of the normal components of the viscous stress tensor.

The LBGK model in (2.55) is generic with respect to the choice of reference pressure $P_0$ and the force $\boldsymbol {F}$. We now proceed with the special case of Korteweg's force to establish a representation thereof matched to the lattice Boltzmann system.

2.3.4. Pseudo-potential and capillarity

Following the representation in (2.14), Korteweg's force includes two distinct parts, the term supplying the non-ideal gas equation of state and the capillarity term responsible for the surface tension. Introducing a pseudo-potential $\psi$,

(2.60)\begin{equation} \psi=\begin{cases} \sqrt{ P-P_0} & \text{if } P>P_0,\\ \sqrt{P_0-P} & \text{if } P \leq P_0, \end{cases} \end{equation}

Korteweg's force is written as

(2.61)\begin{equation} \boldsymbol{F} =\begin{cases} 2\psi\boldsymbol{\nabla}\psi - \kappa\rho\boldsymbol{\nabla}\boldsymbol{\nabla}^2\rho & \text{if } P>P_0,\\ -2\psi\boldsymbol{\nabla}\psi - \kappa\rho\boldsymbol{\nabla}\boldsymbol{\nabla}^2\rho & \text{if } P \leq P_0. \end{cases} \end{equation}

In the lattice Boltzmann setting, the pseudo-potential part is represented as a linear combination of the first- and second-neighbours contributions,

(2.62)\begin{equation} \delta t \psi(\boldsymbol{r})\boldsymbol{\nabla}\psi(\boldsymbol{r})= \psi(\boldsymbol{r})\sum_{i=0}^{Q-1} \frac{w_i}{\varsigma^2} \boldsymbol{c}_i [\mathcal{G}_1\psi(\boldsymbol{r}+\boldsymbol{c}_i \delta t) + \mathcal{G}_2\psi(\boldsymbol{r}+2\boldsymbol{c}_i \delta t)] + {O}([\delta r\boldsymbol{\nabla}]^5), \end{equation}

where the weights $w_i$ are defined by the product-form in (2.49) at $\xi _{\alpha }=0$, $\zeta _{\alpha \alpha }=\varsigma ^2$,

(2.63)\begin{equation} w_i=\prod_{\alpha=x,y,z}\varPsi_{c_{i\alpha}}(0,\varsigma^2), \end{equation}

and where the parameters $\mathcal {G}_1$ and $\mathcal {G}_2$ satisfy the conditions,

(2.64)$$\begin{gather} \mathcal{G}_1+2\mathcal{G}_2=1, \end{gather}$$

(2.65)$$\begin{gather}\mathcal{G}_1+8\mathcal{G}_2=0. \end{gather}$$

The condition in (2.64) maintains the equation of state, while the condition in (2.65) eliminates the third-order error. Non-compliance with the first and/or second condition would introduce respectively deviations of order ${O}([\delta r\boldsymbol {\nabla }])$ and/or ${O}([\delta r\boldsymbol {\nabla }]^3)$ from the target Korteweg stress tensor.

The capillarity part of Korteweg's force in (2.61) is represented in a similar way but using the density instead of the pseudo-potential,

(2.66)\begin{align} \delta t \tilde{\kappa}\rho(\boldsymbol{r})\boldsymbol{\nabla}\boldsymbol{\nabla}^2\rho(\boldsymbol{r}) &= \rho(\boldsymbol{r})\sum_{i=0}^{Q-1} \frac{w_i}{\varsigma^2} \boldsymbol{c}_i [\mathcal{G}_3\rho(\boldsymbol{r}+\boldsymbol{c}_i \delta t) + \mathcal{G}_4 \rho(\boldsymbol{r}+2\boldsymbol{c}_i \delta t)]\nonumber\\ &\quad + {O}([\delta r\boldsymbol{\nabla}]^5), \end{align}

where $\tilde {\kappa }=\kappa \delta r^2$ and the parameters $\mathcal {G}_3$ and $\mathcal {G}_4$ satisfy the conditions,

(2.67)$$\begin{gather} \mathcal{G}_3+2\mathcal{G}_4=0, \end{gather}$$

(2.68)$$\begin{gather}\mathcal{G}_3+8\mathcal{G}_4=6\kappa. \end{gather}$$

The condition in (2.67) cancels the first-order derivative, while the condition in (2.68) maintains the capillarity contribution. Combining both contributions of the pseudo-potential in (2.62) and the capillarity in (2.66), we obtain the lattice Boltzmann form of Korteweg's force in (2.61) as

(2.69)\begin{align} \delta t\boldsymbol{F} &={\pm} 2\psi(\boldsymbol{r})\sum_{i=0}^{Q-1} \frac{w_i}{\varsigma^2} \boldsymbol{c}_i \left[\frac{4}{3}\psi(\boldsymbol{r}+\boldsymbol{c}_i\delta t) - \frac{1}{6} \psi(\boldsymbol{r}+2\boldsymbol{c}_i\delta t)\right]\nonumber\\ &\quad +\tilde{\kappa}\rho(\boldsymbol{r})\sum_{i=0}^{Q-1} \frac{w_i}{\varsigma^2} \boldsymbol{c}_i [2\rho(\boldsymbol{r}+\boldsymbol{c}_i\delta t) - \rho(\boldsymbol{r}+2\boldsymbol{c}_i\delta t)] + {O}([\delta r\boldsymbol{\nabla}]^5), \end{align}

where the sign convention follows (2.61).

In the context of the lattice Boltzmann method, various pseudo-potential representations have long been in use, and comments are to make a distinction to the present formulation.

• • By setting $\mathcal {G}_2=0$ in (2.62) and ignoring Korteweg's capillarity term in (2.66) by choosing $\mathcal {G}_3=\mathcal {G}_4=0$, one recovers a model first proposed by Shan & Chen (Reference Shan and Chen1993) for special equations of state (see (G1), (G2), (G3) in Appendix G) and extended to a general equation of state by Yuan & Schaefer (Reference Yuan and Schaefer2006). Unlike the above condition in (2.65) which eliminates the third-order error in the pseudo-potential part of Korteweg's force in (2.14), the model of Shan & Chen (Reference Shan and Chen1993) requires the third-order Korteweg-like term to be retained so that it mimics surface tension effects. Consequently, the force in the models of Shan & Chen (Reference Shan and Chen1993), Yuan & Schaefer (Reference Yuan and Schaefer2006) becomes

  (2.70)\begin{equation} \boldsymbol{F}_{SC}=2\psi\boldsymbol{\nabla}\psi+\tfrac{1}{3}\psi\boldsymbol{\nabla}\boldsymbol{\nabla}^2\psi. \end{equation}
  The force in (2.70) neither conforms with the van der Waals second-gradient theory and Korteweg's stress of § 2.1 (unless $\psi =A\rho$ which, however, does not lead to phase separation), nor can it be derived from the BBGKY equation with a central-force particles interaction of § 2.2. Another relatively minor issue is the fixed parameter $1/3$, which is mimicking the capillarity coefficient and is the result of the third-order error retained in the expansion.

• • By imposing the condition in (2.64) while still discarding Korteweg's capillarity contribution in (2.66), one arrives at a dual-range force model of Sbragaglia et al. (Reference Sbragaglia, Benzi, Biferale, Succi, Sugiyama and Toschi2007), Shan (Reference Shan2008),

  (2.71)\begin{equation} \boldsymbol{F}_{DR}=2\psi\boldsymbol{\nabla}\psi+(\tfrac{4}{3}-\mathcal{G}_1) \psi\boldsymbol{\nabla}\boldsymbol{\nabla}^2\psi. \end{equation}
  The model in (2.71) is an improvement on the first-neighbour model in (2.70), in that it allows for a variable capillarity-like coefficient. At the same time, it does not resolve the inconsistency with Korteweg's stress tensor.

• • Various other modifications of the original models of Shan & Chen (Reference Shan and Chen1993) were proposed to improve on the main inconsistency by introducing and tuning ad hoc numerical coefficients tailored to a selected equation of state (Kupershtokh et al. Reference Kupershtokh, Medvedev and Karpov2009; Li, Luo & Li Reference Li, Luo and Li2013; Huang, Yin & Killough Reference Huang, Yin and Killough2019; Luo, Fei & Wang Reference Luo, Fei and Wang2021). However, to the best of our knowledge, none of these proposals came to recognize that the problem lies in the fact that it is impossible to represent both parts of Korteweg's force, the non-ideal equation of state and the capillarity term, while using a pseudo-potential alone. This fact follows both from the phenomenological thermodynamics reviewed in § 2.1, as well as from a more microscopic approach of § 2.2. The pseudo-potential in both cases represents only the equation of state, while the capillarity term requires the density field to be used, as featured by (2.66).

• • Pseudo-potential is a convenient form of representing the equation of state contribution to Korteweg's force, tailored to the lattice Boltzmann setting. If other numerical methods are used to evaluate the force in (2.61), such as higher-order finite difference, the model is usually renamed to a free energy approach.

With the generic LBGK model in (2.55) and Korteweg's force in (2.69) both specified, we proceed to the analysis of the hydrodynamic limit under a suitable scaling.

2.3.5. Hydrodynamic limit under small velocity increment scaling

Chapman–Enskog analysis of the LBGK equation in (2.55) was performed under the following scaling. With the characteristic values of the flow velocity $\mathcal {U}$, the flow scale $\mathcal {L}$, the density $\rho$, the force $\mathcal {F}$ and the velocity increment $\delta u$, the following conditions apply:

(2.72)$$\begin{gather} \frac{\delta u}{\mathcal{U}}\sim \frac{\delta t \mathcal{F}}{\rho \mathcal{U}}\sim\varepsilon, \end{gather}$$

(2.73)$$\begin{gather}\frac{\delta r}{\mathcal{L}}\sim\varepsilon. \end{gather}$$

The first scaling condition in (2.72) refers to a smallness of velocity increment, that is, to the smallness of the force action over time $\delta t$. The second scaling condition in (2.73) is a resolution requirement. Both conditions are assumed to hold simultaneously. Details of the analysis are provided in Appendix D while the result is summarized below.

Let us introduce a transformed velocity $\boldsymbol {U}$ by shifting the local velocity $\boldsymbol {u}$ by half of the velocity increment $\delta \boldsymbol {u}$ in (2.43),

(2.74)\begin{equation} \boldsymbol{U}=\boldsymbol{u}+\varepsilon\frac{\delta t\boldsymbol{F}}{2\rho}. \end{equation}

Then the following mass and momentum balance equations to second order in $\varepsilon$ are recovered when the force in (2.69) is used under the scaling in (2.72) and (2.73),

(2.75)$$\begin{gather} \partial_t \rho + \varepsilon\boldsymbol{\nabla}\boldsymbol{\cdot}\rho\boldsymbol{U} + {O}(\varepsilon^3) = 0, \end{gather}$$

(2.76)$$\begin{gather}\partial_t \rho \boldsymbol{U} + \varepsilon \boldsymbol{\nabla}\boldsymbol{\cdot} \rho \boldsymbol{U}\otimes\boldsymbol{U} + \varepsilon\boldsymbol{F}_{K} + \varepsilon\boldsymbol{\nabla}\boldsymbol{\cdot}\varepsilon\boldsymbol{\mathsf{T}}_{NS} + {O}(\varepsilon^3) = 0, \end{gather}$$

where the dynamic and the bulk viscosity in the Navier–Stokes stress tensor in (2.10) are related to the relaxation parameter $\omega$ and the reference pressure $P_0$ as follows:

(2.77)$$\begin{gather} \mu= \left(\frac{1}{\omega} - \frac{1}{2}\right)\delta t P_0, \end{gather}$$

(2.78)$$\begin{gather}\eta= \left(\frac{5}{3} - \frac{\partial\ln P_0}{\partial\ln\rho}\right) \left(\frac{1}{\omega} - \frac{1}{2}\right)\delta t P_0. \end{gather}$$

Unlike the previous result in (2.37), the momentum balance in (2.76) includes not only the non-ideal gas pressure but also the capillarity term, and is thus consistent with Korteweg's force in the momentum balance. It should be pointed out that the scaling in (2.72) refers to the smallness of the increment of the flow velocity rather than to the smallness of either the time step or of the force. Thus, rescaling the kinetic model in (2.35) based on the smallness of the flow velocity increments results in both the non-ideal gas equation of state and the capillarity revealed at the Euler level $O(\varepsilon )$ of the momentum balance in (2.76). This is in contrast to the conventional scaling, which is tied to the non-uniformity and surface tension that would appear only at a Burnett level $O(\epsilon ^3)$.

2.3.6. Code structure

The structure of the algorithm is given in figure 1. The algorithm consists of four main building blocks: (a) evaluation of macroscopic variables $\rho$, $\boldsymbol {u}$ and $P_0$ via (2.50); (b) evaluation of non-local macroscopic contributions $\varPhi _{\alpha \alpha }$ and $F_{\alpha }$ via (2.58), (2.69) and (2.60); (c) evaluation of equilibrium and extended equilibrium functions via (2.54) and (2.59) and (d) discrete space/time evolution following (2.55).

Figure 1. Flow chart for a typical code based on the proposed two-phase lattice Boltzmann model.

In the remainder of this paper, and without loss of generality, we set the reference pressure $P_0=\varsigma ^2\rho$ to minimize the correction term in (2.58). This choice of reference pressure follows the third-order quadrature leading to the family of lattices used here and is quite common in both ideal and non-ideal fluid lattice Boltzmann simulations (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017). Other possibilities of choosing the reference pressure are illustrated in Appendix H. In the next section, the model is scrutinized by a set of numerical tests probing various aspects of thermo- and hydrodynamic consistency.

3.1. Liquid–vapour coexistence

We begin with the validation of liquid–vapour coexistence. Two-dimensional flat interface simulations were performed on a grid $800\times 10$, filled with the vapour phase of a fluid with a specified equation of state and periodic boundary conditions. A column of the liquid phase over 400 grid-points was placed at the centre of the domain. Simulations were ran until the steady state was reached. The steady state was monitored via an $L_0$ norm convergence criterion based on the liquid density $\rho _l$ at the centre of the drop and the vapour density $\rho _v$ at a location outside the drop. A theoretical prediction for the coexistence density ratio $\rho _l/\rho _v$ is readily obtained via the equilibrium condition leading to the Maxwell equal-area construction,

(3.1)\begin{equation} \int_{\rho_v}^{\rho_l} \frac{P_{sat}-P}{\rho^2}\,{\rm d}\rho=0, \end{equation}

where $P_{sat}(T)$ is the saturation pressure at which the liquid and vapour phases coexist at a given temperature $T$ below the critical point.

Initially, four generally adopted equations of state (EoS) were considered: the van der Waals EoS (van der Waals Reference van der Waals1873),

(3.2)\begin{equation} P = \frac{\rho R T}{1-b\rho} - a \rho^2, \end{equation}

where parameters $a$ and $b$ are related to critical temperature $T_c$ and pressure $P_c$ as

(3.3a,b)\begin{equation} a=\frac{27}{64}\frac{R^2T_c^2}{P_c},\quad b=\frac{1}{8}\frac{RT_c}{P_c}; \end{equation}

the Peng–Robinson EoS (Peng & Robinson Reference Peng and Robinson1976),

(3.4)\begin{equation} P = \frac{\rho R T}{1-b\rho} - \frac{a \alpha(T) \rho^2}{1+2\rho b - b^2 \rho^2}, \end{equation}

with

(3.5)\begin{equation} \alpha(T) = [1 + (0.37464 + 1.54226\omega' - 0.26992 \omega'^2) (1-\sqrt{T/T_c})]^2, \end{equation}

where $\omega '$ the acentric factor ($\omega '=0.344$ for water), and

(3.6a,b)\begin{equation} a=0.45724\frac{R^2T_c^2}{P_c},\quad b=0.0778\frac{R T_c}{P_c}; \end{equation}

the Riedlich–Kwong–Soave EoS (Redlich & Kwong Reference Redlich and Kwong1949; Soave Reference Soave1972),

(3.7)\begin{equation} P = \frac{\rho R T}{1-b\rho} - \frac{a \alpha(T) \rho^2}{1+\rho b}, \end{equation}

with

(3.8)\begin{equation} \alpha(T) = [1 + (0.480 + 1.574\omega' - 0.176 \omega'^2) (1-\sqrt{T/T_c})]^2, \end{equation}

and

(3.9a,b)\begin{equation} a=0.42748\frac{R^2T_c^2}{P_c},\quad b=0.08664\frac{R T_c}{P_c}; \end{equation}

and the Carnahan–Starling EoS (Carnahan & Starling Reference Carnahan and Starling1969),

(3.10)\begin{equation} P = \rho R T \frac{1+b\rho/4 + {(b\rho/4)}^2 - {(b\rho/4)}^3}{{(1-b\rho/4)}^3} - a\rho^2, \end{equation}

with

(3.11a,b)\begin{equation} a=0.4963\frac{R^2T_c^2}{P_c},\quad b=0.18727\frac{R T_c}{P_c}. \end{equation}

Figure 2 demonstrates that the stationary density ratios $\rho _l/\rho _v$ obtained from the simulation are in excellent agreement with the theoretical coexisting liquid–vapour density ratios that are defined by Maxwell's equal-area rule in (3.1), for all four equations of state and for ratios as high as at least $\rho _l/\rho _v\sim 10^{11}$. It is noted that high coexisting density ratios were obtained without any tuning parameters, universally for all equations of state considered. Therefore, and without loss of generality, in the remainder of this article, we only consider the van der Waals equation of state in (3.2).

Figure 2. Liquid–vapour coexistence for various equations of state. Grey lines, Maxwell's equal-area construction in (3.1). Red symbol, simulation. (a) van der Waals (3.2) ($a=0.000159$, $b=0.0952$); (b) Peng–Robinson (3.4) ($a=0.000159$, $b=0.0952$); (c) Carnahan–Starling (3.10) ($a=0.000868$, $b=4$); (d) Riedlich–Kwong–Soave (3.7) ($a=0.000159$, $b=0.0952$). For all simulations, $\tilde {\kappa }=0.02$.

3.2. The principle of corresponding states and thermodynamic convergence

A discussion on the principle of corresponding states and the necessity of adherence to it in the context of the van der Waals theory of a capillary fluid in the simulation of realistic systems at large density ratios is in order. According to Guggenheim (Reference Guggenheim1945), the principle of corresponding states is ‘the most useful byproduct of van der Waals’ equation of state’. The principle maintains that all properties that depend on inter-molecular forces are related to the critical properties of the substance in a universal way, regardless of the molecular compound of interest. While for real fluids, as shown by Guggenheim (Reference Guggenheim1945), the principle is an approximation, it is an exact property within the van der Waals capillary fluids theory. This is, most notably, shown by the absence of substance-dependent parameters in the non-dimensional form of the van der Waals equation of state and mean-field scaling laws. For the equation of state, the principle of corresponding states implies that the reduced pressure $P_r=P/P_c$ is a universal function of the reduced temperature $T_r=T/T_c$ and of the reduced density $\rho _r=\rho /\rho _c$,

(3.12)\begin{equation} \frac{P}{P_c} = f\left(\frac{T}{T_c}, \frac{\rho}{\rho_c}\right).\end{equation}

The universality of the reduced pressure in (3.12) can be used to write Maxwell's equal-area rule in reduced form,

(3.13)\begin{equation} \int_{\rho_{r,v}}^{\rho_{r,l}} \frac{P_{r,sat}(T_r) - P_r}{\rho_r^2} \,{\rm d}\rho_r=0. \end{equation}

The coexistence density ratio $\rho _l/\rho _v$ at a given reduced temperature $T_r$ is therefore also universal. Deviation from these universal values can be rooted either in numerical artefacts or inconsistent thermodynamics of the model. For a consistent discretization in the sense of Lax's equivalence theorem (Lax & Richtmyer Reference Lax and Richtmyer1956), the former should be vanishing in the limit of a well-resolved simulation. The substance-agnostic form of the non-dimensional equations, i.e. the principle of corresponding states, allows us to probe convergence by systematically thickening the interface.

To probe the consistency with the principle of corresponding states in our model, we performed simulations for different values of the weak attraction parameter $a$ of the van der Waals equation of state in (3.2). Reduced coexistence densities are shown in figure 3 for four different values of $a$. It is observed that, down from the critical point to $T_r\approx 0.4$, coexistence densities essentially overlap for all the four values of $a$, in agreement with the principle of corresponding states. Deviations from the principle of corresponding states are most pronounced on the vapour side, as observed in figure 3. This can be explained as follows. The local value of the scaling parameter $\varepsilon$ in (2.72) is estimated as $\varepsilon _{l,v}\sim \delta t \mathcal {F}/\rho _{l,v}\mathcal {U}$ in the liquid and in the vapour, respectively. Hence, their relative magnitude scales as the inverse of the density ratio, $\varepsilon _{v}/\varepsilon _l\sim \rho _l/\rho _v$. Thus, even if the scaling condition in (2.72) is satisfied for the liquid, $\varepsilon _l\ll 1$, it is still prone to violation for the vapour, if the density ratio $\rho _l/\rho _v$ becomes sufficiently large. Conversely, if the scaling condition is satisfied for vapour, $\epsilon _v\ll 1$, it is also satisfied for liquid. Furthermore, due to the scaling relation between interface width and temperature detailed in § 3.6, the interface gets thinner at lower temperatures. This in turn means that for a given $\delta r$, fewer grid-points resolve the interface with decreasing temperature. Whenever the parameters ${\delta r}/{W}$ and ${\delta t F}/{\rho \mathcal {U}}$ increase, contributions of higher orders in $\varepsilon$ become significant and lead to deviations from the analytical predictions. The characteristic interface width scales as $W \propto 1/\sqrt {a}$ (Jamet et al. Reference Jamet, Lebaigue, Coutris and Delhaye2001). As such, lower values of $a$ lead to larger $W$, in parallel with smaller force increments over $\delta t$, which restores dominance of order $\varepsilon$ terms over Burnett and super-Burnett level contributions, and therefore the corresponding states principle. This last point is demonstrated by the convergence of the coexistence density ratio to the analytical predictions with decreasing $a$. For well-resolved simulations, as shown in figure 2, the model correctly recovers the coexistence densities and thus complies with the principle of corresponding states. Below, we refer to the convergence of the scheme to the principle of corresponding states as the thermodynamic convergence.

Figure 3. Convergence to the principle of corresponding states. Coexistence densities as obtained from (grey lines) Maxwell's construction and (markers) simulations with different choices of $a$: (green $\times$) $a=0.0102$; (magenta $+$) $a=0.0051$; (blue triangles) $a=0.0026$; (red squares) $a=0.0013$ and (black circles) $a=0.00064$.

To further clarify the utility of the principle of corresponding states in probing convergence, we have also conducted simulations similar to those reported in figure 3, with the first-neighbour pseudo-potential model (Shan & Chen Reference Shan and Chen1993), (2.70), supplemented with the van der Waals equation of state (Yuan & Schaefer Reference Yuan and Schaefer2006). The results are shown in figure 4. Two points are worth noting here. First, a notion of convergence is present as co-existence densities do converge to a limit. However, that limit does not correspond to the capillary fluid model. Deviations become pronounced at lower temperatures where density ratios are larger. These differences are explained by the thermodynamic consistency issue where the mechanical stability condition governing the classical pseudo-potential model is recovered as (Benzi et al. Reference Benzi, Biferale, Sbragaglia, Succi and Toschi2006; Shan Reference Shan2008)

(3.14)\begin{equation} \int_{\rho_{r,v}}^{\rho_{r,l}} (P_{r,sat}(T_r) - P_r ) \frac{1}{\psi_r^2}\frac{\partial \psi_r}{\partial \rho_r} \,{\rm d}\rho_r=0, \end{equation}

which is different from (3.13). Here $\psi _r=\psi /\psi _c$, where $\psi _c$ is the pseudo-potential at critical state. This in turn naturally leads to observations of figure 4 where the discrete solver is convergent, but to a continuum-level behaviour dictated by (3.14) rather than the capillary fluid of (3.13). It is worth nothing that a number of numerical modifications have been proposed to improve on this issue. A detailed overview of such changes is out of the scope of the present work. Interested readers are referred to Chen et al. (Reference Chen, Kang, Mu, He and Tao2014) for a review.

Figure 4. Coexistence densities as obtained from (grey lines) Maxwell's construction and (markers) simulations with different choices of $a$: (green $\times$) $a=0.0102$; (magenta $+$) $a=0.0051$; (blue triangles) $a=0.0026$; (red squares) $a=0.0013$ and (black circles) $a=0.00064$. Simulations are conducted using the first-neighbour pseudo-potential model in (2.70).

3.3. Temperature dependence of the surface tension near the critical point

Surface tension at the liquid–vapour interface decreases with increasing temperature and vanishes at the critical point (Guggenheim Reference Guggenheim1945). For the van der Waals equation of state, the surface tension coefficient $\sigma$ follows a scaling law as $T_r\to 1$ (van der Waals Reference van der Waals1894; Blokhuis & Kuipers Reference Blokhuis and Kuipers2006),

(3.15)\begin{equation} \sigma = \frac{16a}{27b^2}\sqrt{\frac{\kappa}{a}} {(1-T_r)}^{3/2}. \end{equation}

To probe the consistency of the proposed lattice Boltzmann model, the temperature dependence of the surface tension was numerically measured in two different configurations, the flat interface and the circular drop, in a temperature interval $T_r\in [0.85,\, 1]$.

In the first configuration, the surface tension coefficient was evaluated using its definition for the flat interface (Kirkwood & Buff Reference Kirkwood and Buff1949),

(3.16)\begin{equation} \sigma = \int_{-\infty}^{+\infty} (P_{xx} - P_{yy})\,{\rm d} x, \end{equation}

where the interface is normal to the $x$-axis in a two-dimensional simulation setup. The normal $P_{xx}$ and the tangential $P_{yy}$ components of the discrete pressure tensor were computed using a formalism developed in Appendix E, following a proposal by Shan (Reference Shan2008).

In the second configuration, simulations of circular liquid drops surrounded with vapour at the centre of a square domain were conducted. At each temperature, four different initial drop radii were considered, $R_0\in \{45,55,65,75\}$, chosen in such a way that the interface width is sufficiently small compared to the initial drop radius. In the simulation, we used $W/R_0 \le 0.1$, where $W=(\rho _l-\rho _v)/\max \lvert \boldsymbol {\nabla }\rho \lvert$ is the interface width, to minimize the curvature dependence of surface tension. The corresponding surface tension coefficient was evaluated using the Laplace law ($D=2$) in a form,

(3.17)\begin{equation} \Delta P = \frac{(D-1)\sigma}{R_e}, \end{equation}

where the radius $R_e$ corresponds to the equimolar dividing surface (Gibbs Reference Gibbs1874).

A brief reminder of Gibbs’ theory of dividing surfaces is in order. The total mass in both the diffuse and sharp interface pictures can be written as

(3.18)\begin{equation} \int_{V} \rho \,{\rm d}V = \rho_l V_l + \rho_v V_v + \varGamma, \end{equation}

where $\rho _l V_l$ and $\rho _v V_v$ are the masses in the bulk liquid and vapour phases in the sharp interface picture, while $\varGamma$ is the excess mass on a dividing surface $\varSigma$, or mass adsorbance (Gibbs Reference Gibbs1874). By requiring that no mass be stored on the dividing surface, we get the definition of the equimolar surface,

(3.19)\begin{equation} \varGamma = 0. \end{equation}

The family of dividing surfaces in the case of drop or bubble are concentric spheres ($D=3$) or concentric circles ($D=2$) parametrized by their radius $R$. In particular, for a two-dimensional drop, the mass adsorbance can be written as a function of the radius of the dividing circle,

(3.20)\begin{equation} \varGamma(R)= \int_0^{2{\rm \pi}}\int_0^{\infty}(\rho(r) - \rho_v)r \,{\rm d}r\,{\rm d}\varphi - \int_0^{2{\rm \pi}}\int_0^{R}(\rho_l - \rho_v)r \,{\rm d}r\,{\rm d}\varphi, \end{equation}

while the zero-adsorption condition in (3.19), $\varGamma (R_e)=0$, implies the equimolar radius $R_e$,

(3.21)\begin{equation} R_e = \sqrt{\frac{\displaystyle\int_0^{\infty}(\rho(r) - \rho_v)r \,{\rm d}r}{(\rho_l - \rho_v)}}. \end{equation}

The drop configuration along with the scaling of the pressure difference across the interface with drop radius are shown in figure 5.

Figure 5. (a) Circular $D=2$ drop configurations; (b) pressure difference scaling with drop radius for $T_r=0.99$, 0.98, 0.97 and 0.96. The pressure difference is defined as $\Delta P=P_{in}-P_{out}$. The slope of the straight line is the surface tension coefficient.

The results as obtained from both the flat interface and the drop configurations are shown in figure 6. It is clearly observed that the surface tensions as obtained from the proposed formulation (using either one of the considered configurations) agree very well with (3.15), provided that $W \ll R_{e}$. Discussion of the curvature dependence of surface tension shall be continued in § 3.5.

Figure 6. Temperature dependence of the surface tension coefficient near the critical point. Dashed grey line, theory using (3.15); red circles, simulation results using Laplace's law in (3.17); blue squares, surface tension coefficient computed by (3.16).

3.4. Control of surface tension

The present formulation allows us to select the surface tension in the simulation via the capillarity parameter $\kappa$, independently from the density ratio and temperature. To demonstrate this feature, flat interface simulations were performed at three reduced temperatures, $T_r=\{0.99, 0.98, 0.97\}$, for different values of $\tilde {\kappa }$. The surface tension was evaluated using (3.16). Results in figure 7 demonstrate that the surface tension can be effectively tuned using $\tilde {\kappa }$ and that changes in surface tension do not affect the equilibrium density ratios. This is further detailed in table 1 where the values of surface tension and deviations of the vapour and liquid densities from theoretical values are given as a function of $\tilde {\kappa }$ for $T_r=0.97$. It is clearly observed that while the surface tension covers two orders of magnitude, the deviation from the predicted density of vapour is at most $0.041$ percent. Furthermore, in the limit of vanishing $\kappa$, the surface tension also vanishes, as expected from the theory, (3.15).

Figure 7. (a) Surface tension as a function of capillarity parameter $\tilde {\kappa }$; (b) liquid/vapour densities as a function of $\tilde {\kappa }$. Results of simulation shown at three different reduced temperatures: blue square, $T_r=0.99$; red circles, $T_r=0.98$; black triangles, $T_r=0.97$. Dashed grey lines, theoretical coexistence densities from the Maxwell construction.

Table 1. Effect of the choice of $\tilde {\kappa }$ on surface tension and deviations in equilibrium vapour and liquid phases densities for $T_r=0.97$.

3.5. Effect of curvature on surface tension: Tolman length

The drop simulation in § 3.3 made use of the Laplace law, relying on the equimolar dividing surface of radius $R_e$ in (3.21). Further discussion on the non-uniqueness of the choice of the dividing surface and the curvature dependence of surface tension is in order. Following Gibbs (Reference Gibbs1874), the free energy of a drop or bubble separated from the surrounding vapour or liquid by a dividing circle ($D=2$) or sphere ($D=3$) of length or area $\varSigma$ is $A=U-TS+\sigma \varSigma$, where $U$ and $S$ are the internal energy and entropy of bulk phases while the last term is the adsorbance of free energy. The equilibrium condition requires vanishing of the variation $\delta A$; for the isothermal case, we have

(3.22)\begin{equation} \delta A={-}P_{l,v}\delta V_{l,v}-P_{v,l}\delta V_{v,l} +\varSigma\delta \sigma+\sigma\delta \varSigma=0, \end{equation}

where $P_{l,v}$ and $P_{v,l}$ are the pressures inside and outside the liquid drop or vapour bubble. Using $\delta V_{l,v}=-\delta V_{v,l}=2(D-1){\rm \pi} R^{D-1} \delta R$ and $\delta \varSigma =2(D-1)^2{\rm \pi} R^{D-2} \delta R$ leads to a generalized Laplace law,

(3.23)\begin{equation} \Delta P = \frac{(D-1)\sigma(R)}{R} + \frac{{\rm d} \sigma(R)}{{\rm d}R}. \end{equation}

The derivative of surface tension ${\rm d}\sigma /{\rm d}R$ is termed a notional derivative by some authors (Blokhuis & Bedeaux Reference Blokhuis and Bedeaux1992) to stress that it refers to arbitrariness of the dividing surface. Apart from the equimolar surface in (3.21), the surface of tension is another possible choice to lift the ambiguity of the dividing surface. The notional derivative vanishes at the surface of tension,

(3.24)\begin{equation} \left.\frac{{\rm d}\sigma}{{\rm d}R}\right|_{R=R_s}=0, \end{equation}

thereby reducing the generalized Laplace law in (3.23) to a standard form,

(3.25)\begin{equation} \Delta P = \frac{(D-1)\sigma(R_s)}{R_s}. \end{equation}

Integrating (3.23) from $R_s$ to $R$, and eliminating $\Delta P$ using (3.25), one obtains an analytic expression for the notional surface tension $\sigma (R)$ relative to its minimum $\sigma _s$ at the surface of tension $R_s$,

(3.26)\begin{equation} \frac{\sigma(R)}{\sigma_s} = \frac{1}{D}{\left(\frac{R_s}{R}\right)}^{D-1} + \frac{D-1}{D}\left(\frac{R}{R_s}\right). \end{equation}

Equation (3.26) provides for a simple way of identifying the surface of tension and the corresponding surface tension. Two-dimensional simulations at $T_r=0.98$ were conducted with different initial drop and bubble sizes, $R_0\in \{30,40,50,60,70,80,100,120,140\}$. For each dividing surface of radius $R$, the corresponding surface tension was evaluated as (Blokhuis & Bedeaux Reference Blokhuis and Bedeaux1992)

(3.27)\begin{equation} \sigma(R) = \int_{0}^{\infty} {\left(\frac{r}{R}\right)}^{D-1}[P_{{\perp}}(r,R) - P_{{\parallel}}(r)]\,{\rm d}r, \end{equation}

where $P_{\parallel }(r)$ is the tangential component of the pressure tensor, computed via the discrete pressure tensor detailed in Appendix E, and

(3.28)\begin{equation} P_{{\perp}}(r,R) = P_{in} - (P_{out} - P_{in}) H(r-R) \end{equation}

is the normal pressure component in the sharp interface system, where $H$ is the Heaviside step function. For each drop or bubble, surface tension $\sigma (R)$ in (3.27) was probed at seventeen equidistant dividing circles between $R_{min}=5\delta r$ and $R_{max}=165\delta r$. The discrete values $\sigma (R)$ obtained in these simulations were fitted with (3.26), with $\sigma _s$ and $R_s$ as the free fitting parameters. Figure 8 shows that the data for all drops and bubbles collapsed on a single master curve are in excellent agreement with the theory in (3.26). Thus, the proposed model correctly identifies the surface of tension for the van der Waals fluid.

Figure 8. Variation of surface tension $\sigma$ with the radius of the dividing circle $R$ at $T_r=0.98$. Grey dashed line, theory in (3.26); symbol, simulations with van der Waals equation of state; blue circles, drop simulations; red diamonds, bubble simulations. Shading from dark to light, drops and bubbles of increasing sizes are shown.

While the notional derivative ${\rm d}\sigma (R)/{\rm d}R$ vanishes at the surface of tension, the same does not hold for the surface tension at the surface of tension, ${\rm d}\sigma _s/{\rm d}R_s\neq 0$. In other words, surface tension $\sigma _s$ depends on the curvature of the surface of tension. In a seminal paper, Tolman (Reference Tolman1949) characterized the curvature dependence of surface tension by the Tolman length. For sufficiently large $R_s$, the leading-order curvature dependence of the surface tension may be written as (Tolman Reference Tolman1949)

(3.29)\begin{equation} \sigma(R_s) \approx \sigma_0\left(1 \mp \frac{(D-1)\delta_T}{R_s}\right), \end{equation}

where $\delta _T$ is the Tolman length and $\sigma _0$ is the flat interface surface tension coefficient. Here, the negative (positive) sign corresponds to drops (bubbles).

We first compare the leading-order Tolman model in (3.29) with the simulation. The values of $R_s$ obtained in the previously described drops and bubbles simulations of figure 8 are plotted against the pressure difference for different drop and bubble sizes in figure 9. It is clear that for smaller drops and bubbles, the pressure difference deviates from the Laplace law with constant $\sigma _s=\sigma _0$, indicating a curvature-dependent surface tension. Fitting the data points with (3.29), the Tolman length can be extracted from the simulation, here $\delta _T=9\delta r$ for both drops and bubbles, at the reduced temperature $T_r=0.98$.

Figure 9. Pressure difference scaling with surface of tension radius $R_s$ for liquid drops and vapour bubbles. Symbol, simulation data for van der Waals equation of state with $a=0.18$ and $b=0.095$l; black line, Laplace law with $\sigma =\sigma _0$; blue lines, best fit with (3.29) used to compute Tolman length $\delta _T$; red lines, best fit with the second-order Helfrich expansion in (3.30).

While the leading-order Tolman correction in (3.29) improves the agreement with data at moderate $R_s$, deviation persists for smaller drops and bubbles at $\delta r/R_s> 0.03$. Higher-order terms in the inverse powers of $R_s$, important for droplets or bubbles of small radius, are neglected by (3.29) and are addressed by Helfrich (Reference Helfrich1978)

(3.30)\begin{equation} \sigma = \sigma_0 \mp \sigma_0\frac{(D-1)\delta_T}{R_s} + \frac{k{(D-1)}^2}{2R_s^2} + \frac{\bar{k}(D-2)}{R_s^2} + \cdots, \end{equation}

where $k$ and $\bar {k}$ are the bending and Gaussian rigidities; note that the latter vanishes for $D=2$. Taking the second-order term in (3.30) into account, the best fit in figure 9 results in bending rigidity $k=1.049\times 10^{5}\sigma _0\delta r^2$.

Finally we consider the limit of a flat interface where the Tolman length can be derived independently from the above considerations. In this case, the location of the surface of tension $X_s$ can be found as the normalized first-order moment of the normal stress difference (Rao & Berne Reference Rao and Berne1979),

(3.31)\begin{equation} X_s = \frac{\displaystyle\int_{-\infty}^{\infty} x(P_{xx} - P_{yy})\, {\rm d} x}{\displaystyle\int_{-\infty}^{\infty} (P_{xx} - P_{yy}) \,{\rm d} x}. \end{equation}

With the dividing surface as a vertical straight line at $X$ in two dimensions, the mass adsorbance $\varGamma (X)$ is defined as (see example in figure 10)

(3.32)\begin{equation} \varGamma(X)= \int_{-\infty}^{\infty} [\rho(x) - \rho_v - (\rho_l-\rho_v) H(x-X)]\,{\rm d} x. \end{equation}

Similar to the case of cylindrical symmetry considered in § 3.3, the equimolar surface is found by annihilating the mass adsorbance, $\varGamma (X_e)=0$.

Figure 10. Example of mass adsorbance for a flat interface. Black continuous line, density profile at $T_r=0.98$; red dashed line, sharp interface profile with the dividing surface at $X/\delta r=70$. Grey area represents the mass adsorbance $\varGamma (X)$ in (3.32).

The Tolman length in the limit of the flat interface is the distance between the surface of tension and the equimolar surface (Blokhuis & Bedeaux Reference Blokhuis and Bedeaux1992; Blokhuis & Kuipers Reference Blokhuis and Kuipers2006),

(3.33)\begin{equation} \delta_T=X_e - X_s. \end{equation}

An example from the simulation is presented in figure 11. As shown by Blokhuis & Bedeaux (Reference Blokhuis and Bedeaux1992), the Tolman length scales with the reduced temperature as

(3.34)\begin{equation} \delta_T\propto{(1-T_r)}^{{-}1}. \end{equation}

To validate the scaling of (3.34) in our model, flat interface simulations were conducted over a range of temperatures in the vicinity of the critical state. The Tolman length in (3.33) for various temperatures is shown in figure 11. The results obtained from simulations are in excellent agreement with the theoretical scaling by (3.34). A similar study using the Shan–Chen equations of states and the pseudo-potential model was presented in (Lulli et al. Reference Lulli, Biferale, Falcucci, Sbragaglia and Shan2021). Furthermore, the flat interface simulations lead to a $\delta _T=9.2\delta r$ at $T_r=0.98$, in agreement with the value obtained from drop and bubble simulations, $\delta _T=9\delta r$ at the same temperature. The relatively small discrepancy can be attributed to higher-order terms in the curvature, neglected in the fitting process.

Figure 11. (a) Temperature dependence of Tolman length $\delta _T$. Results from the flat interface simulation for van der Waals fluid are shown with blue squares while the grey dashed line represents theoretical scaling by (3.34). (b) Surface of tension, equimolar surface and Tolman length for a flat interface. Continuous black line, density profile at $T_r=0.98$; red dashed line, sharp approximation with the equimolar surface as the dividing surface; blue dotted line, sharp approximation with the surface of tension as the dividing surface. Distance between the surface of tension and the equimolar surface is the Tolman length. In all simulations, $a=0.18$ and $b=0.095$.

3.6. Interface width: temperature scaling and control

In the present work, we use a definition for the interface width bearing numerical information as to how well the stiff gradients are resolved on a given mesh, making it directly related to the velocity increment per time step:

(3.35)\begin{equation} W = \frac{\rho_l-\rho_v}{\max\lvert\boldsymbol{\nabla}\rho\lvert}, \end{equation}

where $\rho _l$ and $\rho _v$ are the densities of saturated liquid and vapour, respectively. It can readily be observed that in the limit of a sharp interface, i.e. resolved with $\delta r$, $\delta r/W\to 1$. As previously noted for the co-existence densities and surface tension, the model recovers thermodynamical properties/scalings of the second-gradient fluid in the limit $\varepsilon \to 0$, akin to $\delta r/W\to 0$. As such, in the limit of $\delta r/W \to 1$, one expects numerical effects to dominate and observe deviations from thermodynamics of the van der Waals fluid.

Surface tension vanishes as the temperature approaches the critical, cf. § 3.3, (3.15) and figure 6, while the interface diverges as $T\to T_c$. As noted by Widom (Reference Widom1965), the van der Waals theory predicts the temperature scaling of the interface width as

(3.36)\begin{equation} W(T_r)\propto (1-T_r)^{{-}1/2}. \end{equation}

To validate the consistency of the proposed model, simulations of flat interfaces were carried out in a range of reduced temperatures $T_r$ near the critical point, and corresponding interface widths $W(T_r)$ in (3.35) were measured. Given the effect of the choice of $a$ on interface thickness, simulations were also carried out with different $a$. The results obtained from simulations near the critical point are shown in figure 12. As noted by many authors in the literature (Jamet et al. Reference Jamet, Lebaigue, Coutris and Delhaye2001), the parameter $a$ can be used to control the interface thickness, as done in the present work, at a given density ratio, leaving the ratios and Maxwell construction unaffected. In agreement with the equivalent states theory, the right-hand side plot in figure 12 points to the universality of the scaling of the interface width near the critical point regardless of the choice of $a$. In addition, it is interesting to note that for a fixed grid-size $\delta r$, as $(1-T_r)\to 1$, $\delta r/W\to 1$ (equivalent to the scaling parameter $\varepsilon$ introduced in the multi-scale analysis), indicating deviation from the thermodynamically converged state. This is illustrated by the deviation of the numerical interface thickness, starting at $T_r\approx 0.98$ from the theoretical predictions. Lowering the value of $a$, i.e. rescaling the interface by a factor $1/\sqrt {a}$ and therefore lowering $\varepsilon$, it is observed that interface is again well resolved and the scaling in (3.36) is restored. In agreement with previous sections, in the limit of $\varepsilon \to 0$, the model is shown to be thermodynamically converged and recovers the properties of the second-gradient fluid.

Figure 12. (a) Interface width as a function of temperature. Blue square, simulation with $a=0.184$; red circle, simulation with $a=0.02$; grey dashed line, theoretical scaling by (3.36). (b) Effect of the choice of $a$ on the interface width. Simulation for $T_r= 0.98$ (diamond), 0.99 (square), 0.995 (circle).

4.1. Shear stress: layered Poiseuille flow

The setup consists of a rectangular domain filled with the liquid phase at the bottom and the vapour phase on top. The flow is driven by a body force. Top and bottom are subject to no-slip boundary conditions while the inlet and outlet are fixed by periodicity. The momentum balance at the steady state reduces to a well-posed system of ordinary differential equations,

(4.1a)$$\begin{gather} \partial_y \mu_l \partial_y u + \rho_l g = 0, \quad \forall y: 0\leq y \leq h_l, \end{gather}$$

(4.1b)$$\begin{gather}\partial_y \mu_v \partial_y u + \rho_v g = 0, \quad \forall y: h_l\leq y \leq H, \end{gather}$$

closed by the following boundary conditions:

(4.2a)$$\begin{gather} u\lvert_{y=0}\, = 0, \end{gather}$$

(4.2b)$$\begin{gather}u\lvert_{y=H}\, = 0, \end{gather}$$

(4.2c)$$\begin{gather}u\lvert_{y=h_l^{-}}\, = u\lvert_{y=h_l^{+}}, \end{gather}$$

(4.2d)$$\begin{gather}\mu_l\partial_y u\lvert_{y=h_l^{-}}\, = \mu_v \partial_y u\lvert_{y=h_l^{+}}. \end{gather}$$

Here $H$ is the height of the channel, $h_l$ the height of the section filled with the liquid phase, $\mu _l$ and $\mu _v$ are the dynamic viscosities of the liquid and vapour, respectively, and $g$ is the acceleration due to the body force. An analytical solution is given in Appendix F. This configuration is of particular interest as it probes the ability of a two-phase model to capture jump conditions at the interface which happens only if deviatoric components of the viscous stress tensor in (2.11) are correctly recovered.

Two sets of parameters were considered: (a) $\mu _l/\mu _v=10.1$, $\rho _l/\rho _v=10.1$ ($T_r=0.77$) and (b) $\mu _l/\mu _v=11.3$, $\rho _l/\rho _v=1030$ ($T_r=0.36$). Simulations were conducted on a grid of size $5\times 200$, with a constant acceleration $g=10^{-8}\delta r/\delta t^2$. Furthermore, the same simulations were performed with the conventional second-order equilibrium (Succi Reference Succi2002) instead of the product form in (2.54). The steady-state results are compared to the analytical solution in figure 13. In contrast to the conventional second-order equilibrium, the use of the product form of (2.54) in our model recovers the correct jump conditions and therefore results in continuous velocity profiles, also for a large density ratio, in excellent agreement with analytical solution. A look at the maximum non-dimensional velocities in the domains, respectively $5\times 10^{-4}$ and $0.04$ for density ratios 10 and 1000, further shows that in contrast to single-phase flows where deviations become only important for large non-dimensional velocities, here, due to the presence of large density gradients at the interfaces, third-order contributions to the equilibrium must be retained even for small velocities.

Figure 13. Steady-state velocity profiles for the layered Poiseuille flow. Configuration (a) is where $T_r=0.77$, $\rho _l/\rho _v=10.1$, $\mu _l/\mu _v=10.1$. Configuration (b) is where $T_r=0.36$, $\rho _l/\rho _v=1030$ and $\mu _l/\mu _v=11.3$. Grey plain line, analytical solution; black dashed line, LBM with product-form equilibrium; red dashed line, LBM with conventional second-order equilibrium.

4.2. Normal stress: dissipation of acoustic waves

To assess Galilean invariance of the diagonal components of the viscous stress tensor and the effect of the correction term, the dissipation of acoustic waves was measured numerically. Acoustic waves were initialized as a small perturbation of density with an amplitude $\delta \rho$ and subject to uniform background velocity $U_0$ representing a moving reference frame,

(4.3)\begin{equation} \rho(x,0) = \rho_0 + \delta \rho \sin(2{\rm \pi}/\lambda), \end{equation}

where $\rho _0$ is the density of saturated liquid or vapour at a given temperature and $\lambda$ is the wavelength of the perturbation. The maximum velocity in the domain was monitored over time and fitted to an exponential function of the form $\max (u) - U_0 = \exp (-\varLambda t)$, where the coefficient $\varLambda$ is tied to the dissipation rate of the normal modes as

(4.4)\begin{equation} \eta = \frac{\varLambda}{{(2{\rm \pi}/\lambda)}^2}. \end{equation}

Simulations were performed on a periodic domain of size $L=128\delta r\times \delta r$ with $\lambda =64\delta r$ at the temperature $T_r=0.36$, corresponding to a density ratio $\rho _l/\rho _v\approx 10^3$. The perturbation amplitude was set to $\delta \rho = 10^{-4}\rho _0$ while the reference frame velocity varied in the range $U_0\delta t/\delta r\in [0, 0.3]$. The dissipation rate from (4.4) with and without the correction term are shown in figure 14. It is observed that the correction term for the third-order diagonal moments of the equilibrium populations in (2.58) restores the Galilean invariance (independence from the reference frame velocity) of the dissipation rate of normal modes.

Figure 14. Dissipation rate of acoustic mode at $T_r=0.36$. Circle, saturated liquid; triangle, saturated vapour. Blue, simulation of the model with the correction term; red, simulation of the model without the correction term.

4.3. Viscosity-capillarity coupling: oscillating drop

Rayleigh's classical study of a freely oscillating drop of non-viscous liquid under small deformations laid ground for the understanding of capillary waves (Rayleigh Reference Rayleigh1879). Assuming purely axisymmetric oscillation modes, Rayleigh's oscillation frequency of the $n$th mode, $n=2,3,\dots$, for large density ratios, $\rho _l/\rho _v\gg 1$, is given by

(4.5)\begin{equation} f_n = \frac{1}{2{\rm \pi}}\sqrt{\frac{\sigma n (n-1)(n+2)}{\rho_l R_0^3}}. \end{equation}

Two-dimensional simulations were performed for a drop of initial radius $R_0=40\delta r$, for a relatively low kinematic viscosity, $\nu _l\delta t/\delta r^2=0.01$. Modes of order $n=2,3,4,5,6$ were initiated with a monochromatic perturbation, $R(\theta,0)=R_0(1 + a_n \cos (n\theta ) - \frac {1}{4}na_n^2)$, with the amplitude $a_n=0.1$. The reduced temperature of the van der Waals fluid was set to $T_r=0.36$, corresponding to a density ratio $\rho _l/\rho _v\approx 10^3$. The initial and mid-cycle shapes of Rayleigh's modes are shown in figure 15. Simulations were performed over eight oscillation cycles and corresponding oscillation periods were identified. Figure 15 demonstrates an excellent comparison of the oscillation frequencies measured in the simulation with the corresponding theoretical values from (4.5), even for higher-order modes. While Rayleigh's analysis was restricted to non-viscous liquids, a number of analytical solutions for viscous oscillating drops have been obtained over the years. In particular, Aalilija, Gandin & Hachem (Reference Aalilija, Gandin and Hachem2020) derived a time-dependent solution of the following form:

(4.6)\begin{equation} R(\theta,t) = R_0 ( 1 + \epsilon_n \cos(n\theta) - \tfrac{1}{4}n\epsilon_n^2), \end{equation}

where $\epsilon _n = a_n \exp (-\lambda _n t) \cos (2{\rm \pi} f_n t)$, while $\lambda _n$ is the damping rate of the $n$th mode,

(4.7)\begin{equation} \lambda_n = 2n(n-1){\nu_l}{R_0^{{-}2}}. \end{equation}

Figure 15. (a) Drop oscillation period for modes $n=2,3,4,5,6$. Red circle, simulation; blue square, Rayleigh's modes from (4.5). (b) Corresponding shapes of the drop at $t=1/2\delta t f_n$ and $t=1/\delta t f_n$.

To validate the damping rate of capillary waves, the previously described two-dimensional drop was simulated in a range of kinematic viscosities, $\nu _l\delta t/\delta r^2\in [0.01, 0.1]$, subject to a monochromatic perturbation with the lowest mode $n=2$ and the initial perturbation amplitude was set to $a_2=0.2$. The envelope of oscillations was fitted using the exponential function $\exp (-\varLambda t)$, see figure 16. The effective damping coefficient $\varLambda$ was used along with (4.7) to evaluate the apparent kinematic viscosity. As shown in figure 16, the dissipation rates obtained by simulation agree well with the analytical solution from (4.7).

Figure 16. (a) Time evolution of the drop radius in directions $\theta =0$ and $\theta ={\rm \pi} /2$ for kinematic viscosity $\nu \delta t/\delta r^2=0.02$. Grey line, analytical solution from (4.6); symbol: simulation; dashed black line, the $\exp (-\varLambda t)$ fit used to compute the dissipation rate in (4.7). (b) Apparent viscosity as computed via (4.7). Grey line, analytical solution from (4.7); symbol, simulation. All results correspond to the oscillation mode $n=2$.

4.4. Isothermal speed of sound

Finally, we validate the speed of sound at different temperatures in both the liquid and vapour phases. For the van der Waals fluid, the isothermal sound speed is

(4.8)\begin{equation} c_s = \sqrt{\frac{\partial P}{\partial \rho}\bigg|_{T}} = \sqrt{\frac{RT}{{(b\rho-1)}^2} - 2a\rho}. \end{equation}

The sound speed was measured by monitoring the position of a pressure front over time in a quasi-one-dimensional simulation at different temperatures. The obtained results are compared to theory using (4.8) in figure 17. It is observed that the simulations accurately capture the speed of sound in both liquid and vapour phases, for the entire range of density ratios on the coexistence diagram in figure 2, up to at least $\rho _l/\rho _v\sim 10^{11}$. Results for other equations of state are compared in Appendix G.

Figure 17. Isothermal sound speed for van der Waals fluid at various temperatures. Line, theory using (4.8); symbol, simulation. Upper branch, saturated liquid; lower branch, saturated vapour.

4.5. Interaction with solid boundaries: equilibrium contact angle

A comprehensive review of various implementations of the contact angle on solid boundaries in the lattice Boltzmann setting can be found in Li et al. (Reference Li, Luo, Kang and Chen2014). Here we follow a proposal by Benzi et al. (Reference Benzi, Biferale, Sbragaglia, Succi and Toschi2006) where a virtual density, and therefore a corresponding virtual pseudo-potential, is attributed to the solid nodes. The calculations of the force using (2.69) are then carried at all fluid nodes. The no-slip condition is imposed via the modified bounce-back scheme of Bouzidi, Firdaouss & Lallemand (Reference Bouzidi, Firdaouss and Lallemand2001) for curved boundaries. The virtual density attributed to the solid nodes is the free parameters controlling the contact angle. The ability of this approach to correctly reproduce the dynamics of the contact line in a high-density ratio regime will be validated below.

4.5.1. Young–Laplace equation

As a first validation of the solid/fluid interaction, the case of a two-dimensional channel of height $H$ and length $L$ is considered. Initially a rectangular column of liquid of length $L/3$ and height $H$ surrounded by vapour was placed at the centre of the channel. Simulations were performed until the system reaches steady state, and the angle between the liquid/vapour interface and the solid wall at the triple-point was measured. The static contact angle measured at the triple-point should verify the Young–Laplace law for this configuration:

(4.9)\begin{equation} \Delta P = \frac{2\sigma \sin \theta}{H}, \end{equation}

where $\theta$ is the equilibrium contact angle and $\sigma$ the liquid/vapour surface tension.

Here, to validate the static angle, the reduced temperature of van der Waals fluid was set to $T_r=0.36$, which corresponds to a density ratio $\rho _l/\rho _v=10^3$. The height of the channel is set to $H=75\delta r$ and the length to $L=300\delta r$. The pressure difference $\Delta P$ is computed as the difference in pressure between two monitoring points located respectively inside and outside the liquid column. The convergence criterion used to assess steady state is based on the time variations of pressure at these two monitoring points. The contact angles were measured using the low-bond axisymmetric drop shape analysis module (Stalder et al. Reference Stalder, Melchior, Müller, Sage, Blu and Unser2010) in ImageJ. The results obtained for the entire contact angle range are shown in figure 18. The measured contact angles show excellent agreement with the Young–Laplace equation. Furthermore, the results show that for $\rho _w\to \rho _v$, the contact angle tends to $\theta \to 180^{\circ }$, while as $\rho _w\to \rho _l$, the contact angles tend to $\theta \to 0^{\circ }$ therefore covering the entire range of contact angles.

Figure 18. Static contact angles as (blue squares) obtained from the Young–Laplace equation and (red circles) measured directly from the simulations.

4.5.2. Wettability-driven drop motion in a channel

A difference in contact angle on opposite sides of a liquid drop placed on a surface with a non-uniform wettability may initiate a motion of the drop. As a further validation of the solid/liquid interaction of the proposed model, the case of a liquid column placed in a channel and subjected to a wettability step function is considered. The configuration, illustrated in figure 19, consists of a channel of length $L$ and width/height of $H$. Initially a liquid column of length $L_l$ is placed at the channel centre. Once the liquid column has reached equilibrium, the wettability on the right-half of the channel is changed, which results in unequal contact angles on the different sides of the liquid column. An approximate analytical solution for the centre-of-mass position $X$ and centroid velocity $U$ was found by Esmaili, Moosavi & Mazloomi (Reference Esmaili, Moosavi and Mazloomi2012),

(4.10)$$\begin{gather} X(t) = \tau_{\infty} U_{\infty} ( {\rm e}^{-{t}/{\tau_{\infty}}} + (t/{\tau_{\infty}}) - 1 ), \end{gather}$$

(4.11)$$\begin{gather}U(t) = U_{\infty} ( 1 - {\rm e}^{-{t}/{\tau_{\infty}}}), \end{gather}$$

where the saturated centroid velocity $U_{\infty }$ and transition time $\tau _{\infty }$ are computed respectively as

(4.12)$$\begin{gather} U_{\infty} = \frac{\sigma H (\cos \theta^{-} - \cos \theta^{+} )}{6[ \rho_l \nu_l L_l + \rho_v \nu_v (L - L_l)]}, \end{gather}$$

(4.13)$$\begin{gather}\tau_{\infty} = \frac{H^2[ \rho_l L_l + \rho_v (L - L_l)]}{12[ \rho_l \nu_l L_l + \rho_v \nu_v (L - L_l)]}. \end{gather}$$

Here, $H=70\delta r$ while $L=1260\delta r$, $L_l=420\delta r$, $\nu _l=\nu _v=0.075\delta r^2/\delta t$ and the density ratio is set to $10^3$ as in the previous cases. Initially, the contact angles on both sides of the water column were set to $\theta ^{-}=59.3^{\circ }$. Once a steady state was reached, the contact angle on the right-hand side of the column was changed to $\theta ^{+}=63.4^{\circ }$, which resulted in a net force acting on the liquid. The liquid column centroid velocity and centre-of-mass position are compared to the analytical solutions in figure 20. Both quantities closely match analytical predictions.

Figure 19. Schematics of the wettability-driven liquid column. Arrow indicates the direction of motion of the liquid column when the contact angle $\theta _{+}$ (right) exceeds the contact angle $\theta _{-}$ (left).

Figure 20. (a) Centroid velocity and (b) centre-of-mass position as obtained from (grey line) analytical solution and (red markers) LB simulations for the wettability-driven motion case.

5.1. Thermodynamic convergence and resolution requirements

Two-phase lattice Boltzmann models are routinely used to simulate the impact of liquid drops on solids. While experiments are typically conducted with water and other liquids in air, the applicability of the two-phase LBM is justified whenever the dynamic effects are concerned. However, the majority of LB models do not reach the experimentally relevant density ratios of, say, the water/air density. The proposed LBM was demonstrated to be valid and thermodynamically well-posed at high density ratios. However, thermodynamic convergence, especially at higher density ratios comes at a cost in terms of interface resolution requirements. For instance, at $T_r=0.36$ resulting in a density ratio of $10^3$, to guarantee deviation below one percent from theory in the vapour phase, one must have $W>8\delta r$. In large-scale hydrodynamic-dominated cases targeted in the present section, full thermodynamic convergence is not necessary. The only manifestations of deviation from the thermodynamically converged state that must be kept under control are spurious currents. Here, for density ratios $10^3$, the interface thickness, in (3.35), is fixed at $W=5\delta r$ resulting in a seven percent deviation of the vapour phase density from theory and spurious currents below $10^{-3}\delta r/\delta t$. The Tolman length for this choice of parameter is $\delta _T/\delta r=2.4$ which, for the resolutions considered below ($R_0=75\delta r$), results in $\delta _T/R_0=0.032$ therefore minimizing the impact of curvature on effective surface tension. Below, we consider benchmark simulations of dynamic effects at realistic density ratios of increasing complexity to probe the accuracy and numerical stability of our model in the dynamic setting.

5.2. Contact time on flat non-wetting surfaces

Extensive studies of drop impact dynamics have shown that the contact time on so-called super-hydrophobic surfaces is independent of the Weber number, $We=\rho _l D_0 U_0^2/\sigma$, and only scales with the inertio-capillary time, $\tau _i=\sqrt {\rho _l D_0^3/8\sigma }$, i.e. for a given drop initial diameter $D_0$, the contact time is not affected by the impact velocity (Gauthier et al. Reference Gauthier, Symon, Clanet and Quéré2015). The different stages of the impact process are illustrated in figure 21 through water drops impacting two surfaces with different contact angles. Simulations carried out under the same conditions are shown and already point to very good agreement between simulations and experiments. To quantify the ability of the proposed model to capture the dynamics of drop impact and to investigate the Weber independence of the contact time on non-wetting surfaces, simulations were performed for a wide range of Weber numbers $\textit{We}\ \in [1,40]$, resulting in Reynolds numbers in the range of $\textit{Re}\ \in [400, 1400]$, with $\textit{Re}=U_0 D_0/\nu _l$. Simulations have been performed in boxes of size $4D_0\times 4D_0\times 4D_0$ with $D_0=150\delta r$. The temperature was set to $T_r=0.36$ resulting in a density ratio of $10^3$. The obtained data show very good agreement with experimental contact times measured by Gauthier et al. (Reference Gauthier, Symon, Clanet and Quéré2015) and confirm the Weber independence of the contact times. Both numerical and experimental data are shown in figure 22.

Figure 21. Drop impacting flat solid surface with different contact angles (first and second row) $\theta =180$ and (third and fourth row) $\theta =90$. The Weber and Reynolds number are respectively 3.5 and 750. Experimental data from Vadillo et al. (Reference Vadillo, Soucemarianadin, Delattre and Roux2009) are shown in the first and third rows. Density ratios in both experiments and simulations are $10^3$ and kinematic viscosity ratios are set to 15.

Figure 22. Drop contact times on flat non-wetting surface (contact angle $\theta =165$ for different Weber numbers as obtained from simulations and experiments. Simulations results are shown with red circular markers while experimental data reported by Gauthier et al. (Reference Gauthier, Symon, Clanet and Quéré2015) are illustrated with blue square markers. The dashed grey line represents the average contact time as obtained from simulations, ${\bar {t_c}}/\tau _i=2.4$.

5.3. Reducing the contact time: pancake bouncing

To further reduce the drop contact times, a number of different strategies have been devised during the past decades. Recently, Liu et al. (Reference Liu, Moevius, Xu, Qian, Yeomans and Wang2014) proposed to use macroscopic structures, in the form of tapered posts, to reduce the contact time. It has been shown that above a certain threshold Weber number, these structures can decrease the contact time by approximately $75$ %. This mechanism is also known as pancake bouncing, due to the pancake-like shape of the drop at take off. A detailed numerical study of pancake bouncing using LBM has been presented in (Mazloomi Moqaddam, Chikatamarla & Karlin Reference Mazloomi Moqaddam, Chikatamarla and Karlin2017) for a limited range of densities. Here, to further demonstrate the versatility of our model, we consider a realistic density ratio of $10^3$. The geometry consists of a flat substrate populated by tapered posts of base and tip radii $R_B$ and $R_b$ and height $h$ with a centre-to-centre distance $w$ in a simple square arrangement. The simulations were run following the geometrical configurations considered in experiments (Liu et al. Reference Liu, Moevius, Xu, Qian, Yeomans and Wang2014). Curved wall boundaries were implemented via the modified bounce-back scheme introduced by Bouzidi et al. (Reference Bouzidi, Firdaouss and Lallemand2001). The geometry is illustrated in figure 23. As for the flat substrate simulations, the domain size was set to $4D_0\times 4D_0\times 4D_0$ with $D_0=150\delta r$. Snapshots from the different stages of impact for two different Weber numbers, one below the pancake bouncing threshold and one above, from both simulations and experiments are shown in figure 24. The results show excellent agreement between simulations and experiments. The contact times for different Weber numbers, as obtained from simulations, are compared to the experimental results from Liu et al. (Reference Liu, Moevius, Xu, Qian, Yeomans and Wang2014) in figure 25. It is shown that the numerical simulations not only accurately capture the contact time reduction due to pancake bouncing and the shape of the drop at take off, but they also correctly predict the onset of pancake bouncing. The shape of the drop at take off is assessed via the pancake quality parameter $Q$ defined as the ratio of the drop radius at take off to its radius at maximum spreading (Liu et al. Reference Liu, Moevius, Xu, Qian, Yeomans and Wang2014).

Figure 23. Illustration of the geometry of tapered posts. Configuration follows the experiment setup by Liu et al. (Reference Liu, Moevius, Xu, Qian, Yeomans and Wang2014).

Figure 24. Drop impacting tapered posts at different Weber numbers (first and second rows) $We=28.2$ with pancake bouncing and (third and fourth rows) $We=14.2$. The first and third rows are experiments from Liu et al. (Reference Liu, Moevius, Xu, Qian, Yeomans and Wang2014) while the second and fourth rows are from simulations.

Figure 25. (a) Drop contact times and (b) pancake quality at rebound on tapered posts for different Weber numbers as obtained from simulations and experiments. Simulations results are shown with red circular markers while experimental data reported by Liu et al. (Reference Liu, Moevius, Xu, Qian, Yeomans and Wang2014) are illustrated with blue square markers.

5.4. Extreme density ratios: inertia-dominated coalescence of mercury drops

The sudden and pronounced topological changes involved in the coalescence of drops, especially the formation and subsequent evolution of the neck between them, have been of interest and subject to study for many years. The dynamics of the dimensionless neck radius $r_{neck}/R_0$ (with $R_0$ the drop's initial radius and $r_{neck}$ the neck radius) have been shown to belong to one of two regimes (Eggers, Lister & Stone Reference Eggers, Lister and Stone1999) characterized by the Ohnesorge number, $Oh=\mu /\sqrt {2\rho _l\sigma R_0}$: highly viscous Stokes or inertial. The neck radius evolution over time has been shown to scale either as $r_{neck}/R_0\propto t/\tau _i$ (in the viscous regime) or $r_{neck}/R_0\propto \sqrt {t/\tau _i}$ (in the inertial regime). Here, to better illustrate the ability of the proposed scheme to model flows with extreme density ratios, we will focus on the first regime, more specifically, the coalescence of mercury droplets, with Ohnesorge numbers of the order of $10^{-4}$ (for 1 g drops). We consider a case following the experimental configuration of Menchaca-Rocha et al. (Reference Menchaca-Rocha, Martínez-Dávalos, Núñez, Popinet and Zaleski2001) with two 1 g mercury droplets. To match the proper density ratio of a mercury/air system ($\rho _{Hg}=13\,600\ \text {kg}\ \text {m}^{-3}$), the temperature is set to $T_r=0.267$, resulting in a density ratio of $\rho _l/\rho _v\sim 11\,500$. The two drops, resolved with 150 points on their radii are placed in a rectangular domain of size $400\times 400\times 800$ with a centre-to-centre distance of 304 points. The drops are connected initially via a neck of 4-points radius. The evolution of the drops’ shape over time are compared to experiments from Menchaca-Rocha et al. (Reference Menchaca-Rocha, Martínez-Dávalos, Núñez, Popinet and Zaleski2001) in figure 26. To further illustrate the agreement of numerical simulations with experiments, the evolution of the neck radii $r_{neck}$ over time are shown in figure 27. Both qualitative comparison of the drops’ shape and quantitative comparison of the neck radius show that the presented solver is able to correctly model multi-phase flow dynamics with extreme density ratio. To the best of the authors’ knowledge, this is the first lattice Boltzmann simulation at such high density ratios, regardless of the formulation, i.e. pseudo-potential, free energy, phase-field etc.

Figure 26. Sequential images from different stages of the mercury drops coalescence and sub-sequential capillary waves propagation. First and third-row images are from experiments of Menchaca-Rocha et al. (Reference Menchaca-Rocha, Martínez-Dávalos, Núñez, Popinet and Zaleski2001), while second and fourth rows are from LB simulations. Snapshots are taken at $\Delta t=3.5$ ms intervals.

Figure 27. Time evolution of the neck radii from both (grey line) simulations and (red circular markers) experiments as reported by Menchaca-Rocha et al. (Reference Menchaca-Rocha, Martínez-Dávalos, Núñez, Popinet and Zaleski2001). The dashed blue line represents the $\sqrt {t/\tau _i}$ scaling with $\tau _i=\sqrt {\rho R_0^3/\sigma }$.