2.1. Basic conservation equations

Consider a one-component two-phase system occupied in a domain $\varOmega \in \mathbb{R}^d$ , $d=2,3$ . We have three basic physical quantities to describe the dynamics of the systems, which are the number density $n(\boldsymbol{x},t)$ , the internal energy density $e(\boldsymbol{x},t)$ and the velocity $\boldsymbol{u}(\boldsymbol{x},t)$ . They are functions of space $\boldsymbol{x} \in \varOmega$ and time $t$ . When there is a two-phase interface, the gradient contribution of $n$ to the energy should be considered (Onuki Reference Onuki2007). This leads to a modified internal energy density function

(2.1) \begin{equation} \hat {e}(\boldsymbol{x},t) = e(\boldsymbol{x},t) + \frac {K}{2} |\boldsymbol{\nabla }n|^2, \end{equation}

where $K$ is a positive parameter. By further considering the kinetic energy, the total energy density is given by

(2.2) \begin{equation} \hat {e}_{\textit{tot}}(\boldsymbol{x},t) =\hat {e}(\boldsymbol{x},t) + \frac {m_0 n }{2} \boldsymbol{u}^2, \end{equation}

where $m_0$ is the mass of one molecule. Thus $m_0 n$ defines a density of the fluids, denoted as $\rho (x,t)$ for later use.

The general conservation equations for the number density $n$ , momentum density $m_0 n \boldsymbol{u}$ and total energy density $e_{\textit{tot}}$ are given by

(2.3) \begin{align} \frac {\partial n}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(n \boldsymbol{u}) &= 0, \end{align}

(2.4) \begin{align} m_0\left ( \frac {\partial n \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(n \boldsymbol{uu})\right ) &= \boldsymbol{\nabla }\boldsymbol{\cdot }\hat { \boldsymbol{\sigma }}, \end{align}

(2.5) \begin{align} \frac {\partial e_{\textit{tot}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(e_{\textit{tot}} \boldsymbol{u}) &= \boldsymbol{\nabla }\boldsymbol{\cdot }(\hat { \boldsymbol{\sigma }} \boldsymbol{u}) - \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{j}_{\!e}, \end{align}

where $\boldsymbol{j}_{\!e}$ is the heat flux, and $\hat { \boldsymbol{\sigma }}$ is the stress tensor. These three equations correspond to the laws of mass, momentum and energy conservation in classical mechanics. The system is not closed, as $\boldsymbol{j}_{\!e}$ and $\hat {\boldsymbol{\sigma }}$ are not yet known.

In addition, we need to consider proper boundary conditions. For simplicity, we assume that $\varGamma _S=\partial \varOmega$ is the solid–fluid interface. We also suppose that the inner energy densities on the solid boundary are functions of $n$ given by $e_b(n)$ as in Onuki (Reference Onuki2007). Then the energy conservation equation on the boundary is given by

(2.6) \begin{align} \frac {\partial e_b}{\partial t} + \boldsymbol{\nabla} _{\varGamma } \boldsymbol{\cdot }(e_b \boldsymbol{u}) = \boldsymbol{j}_{\!e} \boldsymbol{\cdot }\boldsymbol{\nu } - \boldsymbol{j}_{\!e, s} \boldsymbol{\cdot }\boldsymbol{\nu }, \end{align}

where $\boldsymbol{j}_{\!e, s}\boldsymbol{\cdot }\boldsymbol{\nu }$ is the thermal flux from the boundary to solid phase, and $\boldsymbol{j}_{\!e}\boldsymbol{\cdot }\boldsymbol{\nu }$ is the thermal flux from the liquid phase to the boundary. Here, $\boldsymbol{\nu }$ is the unit outward normal vector on $\varGamma _S$ . The equation implies that the difference between two fluxes leads to the change of surface energy density $e_b$ . More explanations on the equation can be found in Liu et al. (Reference Liu, Qian and Xu2012). Notice that $e_b$ is a function of $n$ . The above equation is reduced to

(2.7) \begin{equation} e_b^{\prime}(n)\, \frac {\partial n}{\partial t} + \boldsymbol{\nabla} _{\varGamma } \boldsymbol{\cdot }(e_b \boldsymbol{u}) = \boldsymbol{j}_{\!e} \boldsymbol{\cdot }\boldsymbol{\nu } - \boldsymbol{j}_{\!e, s} \boldsymbol{\cdot }\boldsymbol{\nu }. \end{equation}

In addition, we also assume that $\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nu } = 0$ on $\varGamma _S$ , i.e. the boundary is impenetrable. However, the tangential velocity $\boldsymbol u_{\tau }= \boldsymbol{P}\boldsymbol{u}=\boldsymbol{u}-(\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nu })\boldsymbol{\nu }$ may not be zero due to the slip effect. Here, we define the tangential projection operator $\boldsymbol{P} := \boldsymbol{I} - \boldsymbol{\nu } \boldsymbol{\nu }^{\rm T}$ .

To close the system, the specific expressions for the stress tensor and heat flux, along with the necessary boundary conditions, must be determined. These will be derived in the following subsections using an extended Onsager variational principle.

2.2. The entropy functional

Before introducing the Onsager principle for non-isothermal systems, we first specify the entropy density $s(n,e)$ in the system, which is a function of the number density $n$ and the internal energy density $e$ . Standard thermodynamic theory defines the temperature $T(n,e)$ , i.e.

(2.8) \begin{align} \frac {1}{T}=\frac {\partial s(n,e)}{\partial e}. \end{align}

For example, for the van der Waals fluids (see Onuki Reference Onuki2007), we have

(2.9) \begin{align} s (n,e) & = - k_B n \!\left (\! d \ln \hbar + \frac {d}{2} \ln \frac {d \pi k_B n }{ m_0} - \frac {d}{2} \ln (e + \epsilon a^d n^2 ) + \ln \frac {n}{1 - a^d n} \!\right ) + \frac {k_B n (d + 2)}{2},\nonumber\\ T (n,e) & = \frac {2 (e + \epsilon a^d n^2 )}{ d k_B n}, \end{align}

where $k_B$ is the Boltzmann constant, $\hbar$ is the Planck constant, $a$ is the radius of a molecule, and $\epsilon$ is the magnitude of attracting potential. Notice that in the literature, it is common to take the number density $n$ and the temperature $T$ as primary unknowns, with the internal energy $e$ and entropy density $s$ being functions of $n$ and $T$ . A straightforward calculation then yields the formula used in our work. Other examples include the widely used Peng–Robinson equation of state. The expressions for the entropy $s$ and internal energy $e$ can be found in Smejkal & Miky (Reference Smejkal and Miky2017) and Kou & Sun (Reference Kou and Sun2018b ).

Similar to the case for internal energy density, we must also incorporate a gradient contribution into the entropy density when a liquid–vapour interface is present. The modified entropy density is then given by (cf. Onuki Reference Onuki2007)

(2.10) \begin{equation} \hat {s}(n,e) = s (n, e) - \frac {C}{2} |\boldsymbol{\nabla }n|^2, \end{equation}

where $C$ is a positive parameter. In addition, the entropy density on the solid boundary $\varGamma _S$ is also a function of $n$ , given by $s_b(n)$ . Then the total entropy is calculated by

(2.11) \begin{align} \mathcal{S} & := \int _{\varOmega } \hat {s}(n(\boldsymbol{x},t),e(\boldsymbol{x},t))\, \textrm{d}x + \int _{\varGamma _S} s_b (n(\boldsymbol{x},t))\, \textrm{d}A, \end{align}

which is a functional with respect to $n$ and $e$ .

2.3. The extended Onsager principle for non-isothermal systems

To use the Onsager principle for non-isothermal systems, we define a Rayleigh functional,

(2.12) \begin{align} \mathcal{R}( \boldsymbol{u}, \boldsymbol{j}_{\!e}, \boldsymbol{u}_\tau |_{\varGamma _S} ,\dot {n}|_{\varGamma _S} ):=\varPhi -\dot {\mathcal{S}}_{\textit{prod}}, \end{align}

where $\varPhi (\boldsymbol{u}, \boldsymbol{j}_{\!e},\boldsymbol{u}_\tau |_{\varGamma _S}, \dot {n}|_{\varGamma _S})$ is the entropy production function, and $\dot {\mathcal{S}}_{\textit{prod}}$ is the time changing rate of the entropy produced in the system, which is the total entropy minus that from the environment (noticing that the system is not closed due to heat transfer on the solid boundary). The entropy production function of the system is given by

(2.13) \begin{align} \varPhi := &\int _{\varOmega } \left (\frac {\eta }{T} \left\| \frac {\boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u}^{\rm T}}{2} - \frac {1}{d} (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}) \boldsymbol{I} \right\|_F^2 + \frac {\xi }{2 T} (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u})^2 + \frac {1}{2 \lambda } \left|\frac {\boldsymbol{j}_{\!e}}{T}\right|^2\right ) \textrm{d} x \nonumber \\ &{}+ \int _{\varGamma _S} \Big (\frac {\alpha }{2 T} | \boldsymbol{u}_{\tau } |^2 + \frac {\beta }{2 T} \dot {n}^2\Big )\, \textrm{d} A, \end{align}

where $\boldsymbol{I}$ is the unit matrix, $\|{\cdot }\|_F$ denotes the Frobenius norm, $\dot {n}$ is the material derivative of $n$ , $d$ is the spacial dimension, $\lambda$ is thermal conductivity, $\alpha$ is the slip coefficient, $\beta$ is the phenomenological relaxation parameter, and $\eta$ and $\xi$ are the shear and bulk viscosities, respectively. The entropy producing rate is given by

(2.14) \begin{equation} \dot {\mathcal{S}}_{{\textit{prod}}} := \dot {\mathcal{S}} + \int _{\varGamma _S} \frac {\boldsymbol{j}_{\!e, s} \boldsymbol{\cdot }\boldsymbol{\nu }}{T}\, \textrm{d} A-\int _{\varGamma _S}\frac {1}{T} \boldsymbol{P}\hat { \boldsymbol{\sigma }}\boldsymbol{\nu } \boldsymbol{\cdot }\boldsymbol{u}_{\tau }\, \textrm{d} A, \end{equation}

which is the rate of increase of the total entropy in the bulk, plus outgoing entropy flux through the domain boundary, and minus the entropy production rate due to the work done at the boundary. Here, $\boldsymbol{\nu }$ is the outer normal direction on $\varGamma _S$ . For adiabatic and no-slip boundary conditions, the last two terms vanish.

We propose an extended version of the Onsager principle as follows. The constitutive equations and the boundary conditions can be obtained by minimising the Rayleigh functional, i.e.

(2.15) \begin{align} &\min _{\boldsymbol{u}, \boldsymbol{j}_{\!e},\boldsymbol{u}_\tau |_{\varGamma },\dot {n}|_{\varGamma }} \mathcal{R}\quad \text{s.t.}\quad \textrm{(2.3){-}(2.7)}. \end{align}

In the following, we derive in detail the expressions for $\hat { \boldsymbol{\sigma }}$ and $\boldsymbol{j}_{\!e}$ , along with the boundary conditions characterising the wetting and thermal effects.

2.4. Derivations of the constitutive relations and boundary conditions

We first use (2.3)–(2.7) to express $\dot {\mathcal{S}}$ as a functional of $\boldsymbol{u}$ , $\boldsymbol{j}_{\!e}$ , $\boldsymbol{u}_\tau |_{\varGamma _S}$ and $\dot {n}|_{\varGamma _S}$ . Direct calculations lead to

(2.16) \begin{align} \dot {\mathcal{S}} & = {\int _{\varOmega }} \left(\frac {\partial \hat {s}}{\partial t} +\boldsymbol{\nabla }\boldsymbol{\cdot }(\hat {s}\boldsymbol{u} ) \right)\textrm{d} x + \int _{\varGamma _s}\left( \frac {\partial s_b}{\partial t}+\boldsymbol{\nabla} _{\varGamma }\boldsymbol{\cdot }(s_b\boldsymbol{u}_{\tau } ) \right)\textrm{d} A= {\int _{\varOmega }} \frac {\partial \hat {s}}{\partial t}\, \textrm{d} x + \int _{\varGamma _s} \frac {\partial s_b}{\partial t}\, \textrm{d} A\nonumber \\ &= \int _{\varOmega } \left (\frac {\partial s}{\partial e} \frac {\partial e}{\partial t} + \frac {\partial s}{\partial n} \frac {\partial n}{\partial t} - \frac {1}{2} \frac {\partial }{\partial t} (C\, |\boldsymbol{\nabla }n|^2)\right ) \textrm{d} x+ \int _{\varGamma _s} s_b^{\prime}(n)\, \frac {\partial n}{\partial t}\, \textrm{d} A\nonumber \\ & = \int _{\varOmega } \left (\frac {\partial s}{\partial e} \frac {\partial \hat {e}}{\partial t} - \frac {1}{2} \frac {\partial s}{\partial e} \frac {\partial }{\partial t} (K\, |\boldsymbol{\nabla }n|^2)+ \frac {\partial s}{\partial n} \frac {\partial n}{\partial t} - \frac {1}{2} \frac {\partial }{\partial t} (C\, |\boldsymbol{\nabla }n|^2)\right ) \textrm{d} x+ \int _{\varGamma _s} s_b^{\prime}(n)\, \frac {\partial n}{\partial t}\, \textrm{d} A\nonumber \\ & = \int _{\varOmega } \left (\frac {1}{T} \frac {\partial \hat {e}}{\partial t} + \frac {\partial s}{\partial n} \frac {\partial n}{\partial t} - \frac {M(T)}{T}\, \boldsymbol{\nabla }n \boldsymbol{\cdot }\boldsymbol{\nabla }\frac {\partial n}{\partial t}\right ) \textrm{d} x+ \int _{\varGamma _s} s_b^{\prime}(n)\, \frac {\partial n}{\partial t}\, \textrm{d} A. \end{align}

Here, we have used the definition of the temperature, i.e. $ {1}/{T} = {\partial s}/{\partial e}$ , and $M(T):={K+CT}$ .

Consider (2.5). The time derivative of total energy density $e_{\textit{tot}}$ can be evaluated as

(2.17) \begin{align} \frac {\partial e_{\textit{tot}}}{\partial t} = \frac {\partial \hat {e}}{\partial t} + \frac {m_0\, | \boldsymbol{u} |^2}{2} \frac {\partial n}{\partial t} + m_0 n \boldsymbol{u} \boldsymbol{\cdot }\frac {\partial \boldsymbol{u}}{\partial t} = \frac {\partial \hat {e}}{\partial t} - \frac {m_0\, | \boldsymbol{u} |^2}{2} \frac {\partial n}{\partial t} + m_0 \boldsymbol{u} \boldsymbol{\cdot }\frac {\partial (n \boldsymbol{u})}{\partial t} . \end{align}

Direct calculations lead to

(2.18) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }(e_{\textit{tot}} \boldsymbol{u}) &= \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\left(\hat {e} + \frac {m_0 n}{2}\, | \boldsymbol{u} |^2\right) \boldsymbol{u}\right)\nonumber \\ &=\boldsymbol{\nabla }\boldsymbol{\cdot }(\hat {e}\boldsymbol{u})+\frac {m_0 n\, | \boldsymbol{u} |^2}{2} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} + \frac {m_0\, | \boldsymbol{u} |^2}{2} \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }n + m_0 n \boldsymbol{u} \boldsymbol{\cdot }\big ((\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}\big )\nonumber \\ &=\boldsymbol{\nabla }\boldsymbol{\cdot }(\hat {e}\boldsymbol{u})- \frac {m_0\, | \boldsymbol{u} |^2}{2} \boldsymbol{\nabla }\boldsymbol{\cdot }(n \boldsymbol{u}) + m_0 \boldsymbol{u} \boldsymbol{\cdot }\big (\boldsymbol{\nabla }\boldsymbol{\cdot }(n \boldsymbol{uu})\big ). \end{align}

Substituting the above two equations into (2.5) yields

(2.19) \begin{align} &\frac {\partial \hat {e}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\hat {e} \boldsymbol{u}) - \frac {m_0\, | \boldsymbol{u} |^2}{2} \left(\frac {\partial n}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(n \boldsymbol{u})\right) + m_0\boldsymbol{u} \boldsymbol{\cdot }\left(\frac {\partial (n \boldsymbol{u})}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(n \boldsymbol{uu})\right)\nonumber\\&\quad = \boldsymbol{\nabla }\boldsymbol{\cdot }(\hat { \boldsymbol{\sigma }} \boldsymbol{u}) - \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{j}_{\!e}. \end{align}

When combined with (2.3) and (2.4), the above equation reduces to

(2.20) \begin{equation} \frac {\partial \hat {e}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\hat {e} \boldsymbol{u}) = \hat { \boldsymbol{\sigma }} : \boldsymbol{\nabla }\boldsymbol{u} - \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{j}_{\!e}, \end{equation}

where we have used the relation $ \boldsymbol{\nabla }\boldsymbol{\cdot }(\hat { \boldsymbol{\sigma }} \boldsymbol{u})=(\boldsymbol{\nabla }\boldsymbol{\cdot }\hat { \boldsymbol{\sigma }}) \boldsymbol{\cdot }\boldsymbol{u}+ \hat { \boldsymbol{\sigma }} : \boldsymbol{\nabla }\boldsymbol{u}$ .

By substituting (2.20) into (2.17) and integrating by parts, we have

(2.21) \begin{align} \dot {\mathcal{S}} = &\int _{\varOmega } \left(\left(\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\right) + \frac {\partial s}{\partial n}\right) \frac {\partial n}{\partial t} + \frac {1}{T} \bigl(- \boldsymbol{\nabla }\boldsymbol{\cdot }(\hat {e} \boldsymbol{u}) - \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{j}_{\!e} + \hat { \boldsymbol{\sigma }} : \boldsymbol{\nabla }\boldsymbol{u}\bigr)\right) \textrm{d} x \nonumber \\ &{}+ \int _{\varGamma _S}\left( - \frac {M}{T} \frac {\partial n}{\partial \boldsymbol{\nu }} \frac {\partial n}{\partial t} + s_b^{\prime}(n)\, \frac {\partial n}{\partial t} \right)\textrm{d} A. \end{align}

We introduce a notation on the chemical potential, i.e.

(2.22) \begin{align} \hat {\mu } := - T\left( \frac {\partial s}{\partial n} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\right)\right)\!. \end{align}

Then we substitute (2.3) into (2.21) and integrate by parts. This gives

(2.23) \begin{align} \dot {\mathcal{S}} =& \int _{\varOmega } \left (\frac {\hat {\mu }}{T} \boldsymbol{\nabla }\boldsymbol{\cdot }(n \boldsymbol{u}) + \frac {1}{T} \bigl(- \boldsymbol{\nabla }\boldsymbol{\cdot }(\hat {e} \boldsymbol{u}) - \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{j}_{\!e} + \hat { \boldsymbol{\sigma }} : \boldsymbol{\nabla }\boldsymbol{u}\bigr) \right ) \textrm{d} x \nonumber \\ &{}+ \int _{\varGamma _S}\left( - \frac {M}{T} \frac {\partial n}{\partial \boldsymbol{\nu }} + s_b^{\prime}(n) \right)\frac {\partial n}{\partial t}\, \textrm{d} A\nonumber \\ = &\int _{\varOmega } \left (- n \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\frac {\hat {\mu }}{T} + {\hat {e}} \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\frac {1}{T} + \boldsymbol{j}_{\!e} \boldsymbol{\cdot }\boldsymbol{\nabla }\frac {1}{T} - \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \frac {\hat { \boldsymbol{\sigma }}}{T} \right ) \boldsymbol{\cdot }\boldsymbol{u}\right ) \textrm{d} x\nonumber \\ & {}+ \int _{\varGamma _S} \left ( \frac {1}{T} \boldsymbol{u}_{\tau } \boldsymbol{\cdot }(\boldsymbol{P} \hat { \boldsymbol{\sigma }} \boldsymbol{\nu }) - \frac {\boldsymbol{j}_{\!e} \boldsymbol{\cdot }\boldsymbol{\nu }}{T} +\left( - \frac {M}{T} \frac {\partial n}{\partial \boldsymbol{\nu }} + s_b^{\prime}(n) \right)\frac {\partial n}{\partial t} \right ) \textrm{d} A. \end{align}

In the derivation, we have used the relation $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nu }=0$ and the notation for the tangential projection operator $\boldsymbol{P} = \boldsymbol{I} - \boldsymbol{\nu } \boldsymbol{\nu }$ .

We then compute the entropy production rate $\dot {\mathcal{S}}_{ {\textit{prod}}}$ in the whole system:

(2.24) \begin{align} \dot {\mathcal{S}}_{ {\textit{prod}}} ={} & \dot {\mathcal{S}} + \int _{\varGamma _S} \frac {\boldsymbol{j}_{\!e, s} \boldsymbol{\cdot }\boldsymbol{\nu }}{T}\, \textrm{d} A-\int _{\varGamma _S}\frac {1}{T} \boldsymbol{P}\hat { \boldsymbol{\sigma }}\boldsymbol{\nu } \boldsymbol{\cdot }\boldsymbol{u}_{\tau }\, \textrm{d} A \nonumber \\ ={}& \int _{\varOmega } \left (\boldsymbol{u} \boldsymbol{\cdot }\left(- n\, \boldsymbol{\nabla }\frac {\hat {\mu }}{T} + {\hat {e}}\, \boldsymbol{\nabla }\frac {1}{T} - \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \frac {\hat { \boldsymbol{\sigma }}}{T} \right )\right) + \boldsymbol{j}_{\!e} \boldsymbol{\cdot }\boldsymbol{\nabla }\frac {1}{T} \right )\textrm{d} x \nonumber \\ & {}+ \int _{\varGamma _S} \left ( - \frac {(\boldsymbol{j}_{\!e} - \boldsymbol{j}_{\!e, s}) \boldsymbol{\cdot }\boldsymbol{\nu }}{T} +\left( - \frac {M}{T} \frac {\partial n}{\partial \boldsymbol{\nu }} + s_b^{\prime}(n) \right)\frac {\partial n}{\partial t} \right ) \textrm{d} A=:I_1+I_2. \end{align}

The first term $I_1$ on the right-hand side of the equation is already a linear functional with respect to $\boldsymbol{u}$ and $\boldsymbol{j}_{\!e}$ . We further compute the second term $I_2$ . By substituting (2.7) into $I_2$ and integrating by parts, we obtain

(2.25) \begin{align} I_2 =& \int _{\varGamma _S} \left ( - \frac {1}{T} \left ( e_b^{\prime}(n)\, \frac {\partial n}{\partial t} + \boldsymbol{\nabla} _{\varGamma } \boldsymbol{\cdot }(e_b \boldsymbol{u}_{\tau }) \right )+\left( - \frac {M}{T} \frac {\partial n}{\partial \boldsymbol{\nu }} + s_b^{\prime}(n) \right)\frac {\partial n}{\partial t} \right ) \textrm{d} A \nonumber \\ = & \int _{\varGamma _S} \left ( - \frac {1}{T}\, \boldsymbol{\nabla} _{\varGamma } \boldsymbol{\cdot }(e_b \boldsymbol{u}_{\tau }) -\frac {1}{T} \left (M \frac {\partial n}{\partial \boldsymbol{\nu }} + e_b^{\prime}(n) -T\, s_b^{\prime}(n) \right )\frac {\partial n}{\partial t}\right ) \textrm{d} A \nonumber \\ =& \int _{\varGamma _S} \left ( - \frac {1}{T^2}\, \boldsymbol{\nabla} _{\varGamma } T \boldsymbol{\cdot }(e_b \boldsymbol{u}_{\tau }) - \frac {L (n)}{T} \dot {n} + \frac {L (n)}{T} \boldsymbol{u}_{\tau } \boldsymbol{\cdot }\boldsymbol{\nabla} _{\varGamma } n \right ) \textrm{d} A . \end{align}

Here, in the last equation, we have used $\dot {n}= {\partial n}/{\partial t}+\boldsymbol{u}_{\tau } \boldsymbol{\cdot }\boldsymbol{\nabla} _{\varGamma } n$ on $\varGamma _S$ , and $L(n):=M ({\partial n}/{\partial \boldsymbol{\nu }}) + e_b^{\prime}(n) -T\, s_b^{\prime}(n)$ . Notice that the Helmholtz free energy on the surface can be defined as $f_b = e_b - T s_b$ . Then we have $L (n) = M ({\partial n}/{\partial \boldsymbol{\nu }}) + {\partial f_b (n,T)}/{\partial n}|_T$ . Substituting the expression for $I_2$ in (2.24) gives

(2.26) \begin{align} \dot {\mathcal{S}}_{{\textit{prod}}} = &\int _{\varOmega } \left (\boldsymbol{u} \boldsymbol{\cdot }\left(- n\, \boldsymbol{\nabla }\frac {\hat {\mu }}{T} + {\hat {e}}\, \boldsymbol{\nabla }\frac {1}{T} - \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \frac {\hat { \boldsymbol{\sigma }}}{T} \right )\right) + \boldsymbol{j}_{\!e} \boldsymbol{\cdot }\boldsymbol{\nabla }\frac {1}{T} \right )\textrm{d} x \nonumber \\ &{}+ \int _{\varGamma _S} \left. \left ( - \frac {e_b}{T^2}\, \boldsymbol{\nabla} _{\varGamma } T + \frac {L (n)}{T}\, \boldsymbol{\nabla} _{\varGamma } n \right ) - \frac {L (n)}{T} \dot {n} \right ) \textrm{d} A. \end{align}

Now we see that $\dot {\mathcal{S}}_{{\textit{prod}}}$ is reduced to a linear functional with respect to $\boldsymbol{u}$ , $\boldsymbol{j}_{\!e}$ , $\boldsymbol{u}_\tau |_{\varGamma _S}$ and $\dot {n}|_{\varGamma _S}$ .

Based on the above results, the constraint optimisation problem (2.16) is reduced to an unconstrained one, i.e.

(2.27) \begin{equation} \min _{\boldsymbol{u}, \boldsymbol{j}_{\!e}, \boldsymbol{u}_\tau |_{\varGamma _S},\dot {n}|_{\varGamma _S}} \mathcal{R}=\varPhi -\dot {S}_{\textit{prod}}. \end{equation}

From (2.14), we see that $\varPhi$ is a quadratic functional with respect to $\boldsymbol{u} $ , $\boldsymbol{j}_{\!e} $ , $\boldsymbol{u}_\tau |_{\varGamma _S}$ and $\dot {n}|_{\varGamma _S}$ . It is easy to know that the Rayleigh functional $\mathcal{R}(\boldsymbol{u}, \boldsymbol{j}_{\!e}, \dot {n}|_{\varGamma _S},\boldsymbol{u}_\tau |_{\varGamma _S})$ has a unique minimiser, which is derived by setting the first derivatives of $\mathcal{R}$ to be zero. Direct calculations give

(2.28) \begin{equation} \delta \varPhi = \int _{\varOmega }-\boldsymbol{\nabla }\boldsymbol{\cdot }\frac { \boldsymbol{\sigma }}{T}\boldsymbol{\cdot }\delta \boldsymbol{u}+\frac {\boldsymbol{j}_{\!e}\boldsymbol{\cdot }\delta \boldsymbol{j}_{\!e}}{\lambda T^2}\,\textrm{d}x +\int _{\varGamma _S}\frac {1}{T}(\alpha \boldsymbol{u}_\tau +\boldsymbol{P} \boldsymbol{\sigma }\boldsymbol{\nu })\boldsymbol{\cdot }\delta \boldsymbol{u}_\tau +\frac {\beta }{T}\dot {n}\boldsymbol{\cdot }\delta \dot {n}\, \textrm{d}A, \end{equation}

where $ \boldsymbol{\sigma }$ is the dissipative stress tensor defined as

(2.29) \begin{equation} \boldsymbol{\sigma } := \left ( \eta \left(\boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u}^{\rm T}\right) + \left(\xi - \frac {2 \eta }{d}\right) \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{uI} \right )\!. \end{equation}

In the end, we derive the following constitutive equations and boundary conditions:

(2.30) \begin{align} \dfrac {\delta \mathcal{R}}{\delta \boldsymbol{j}_{\!e}} = 0 \quad \Rightarrow & \quad \boldsymbol{j}_{\!e} = - \lambda \boldsymbol{\nabla }T \quad \textrm {in}\ \varOmega , \end{align}

(2.31) \begin{align} \dfrac {\delta \mathcal{R}}{\delta \boldsymbol{u}} = 0 \quad \Rightarrow & \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {1}{T} \hat { \boldsymbol{\sigma }}\right) = \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {1}{T} \boldsymbol{\sigma }\right) - \left(n\, \boldsymbol{\nabla }\frac {\hat {\mu }}{T} - \hat {e}\, \boldsymbol{\nabla }\frac {1}{T}\right) \quad \textrm {in}\ \varOmega , \end{align}

(2.32) \begin{align} \dfrac {\delta \mathcal{R}}{\delta \boldsymbol{u}_{\tau }|_{\varGamma _S}} = 0 \quad \Rightarrow & \quad \beta \dot {n} + L (n) = 0 \quad \textrm {on}\ \varGamma _S, \end{align}

(2.33) \begin{align} \dfrac {\delta \mathcal{R}}{\delta \dot {n}|_{\varGamma _S}} = 0 \quad \Rightarrow & \quad \alpha \boldsymbol{u}_{\tau } + \boldsymbol{P} { \boldsymbol{\sigma }} \boldsymbol{\nu } - \frac {e_b}{T}\, \boldsymbol{\nabla} _{\varGamma } T - L (n)\, \boldsymbol{\nabla} _{\varGamma } n = 0 \quad \textrm {on}\ \varGamma _S. \end{align}

To further simplify (2.31), we use the expression for pressure $ \hat {p} := n \hat {\mu } - \hat {e} + T \hat {s}$ (Onuki Reference Onuki2007). In the following, we will show that $n\, \boldsymbol{\nabla }({\hat {\mu }}/{T}) - \hat {e}\, \boldsymbol{\nabla }({1}/{T}) $ can be rewritten as a divergence form. Notice that $\hat {\mu }=- T (({\partial s}/{\partial n}) + \boldsymbol{\nabla }\boldsymbol{\cdot }(({M}/{T})\, \boldsymbol{\nabla }n) )=\mu -T\,\boldsymbol{\nabla }\boldsymbol{\cdot }(({M}/{T})\, \boldsymbol{\nabla }n)$ , where we denote $\mu =-T({\partial s}/{\partial n}) $ for simplicity in expressions. Direct calculations give

(2.34) \begin{align} \begin{split} \hat {p} & = n (\mu - T \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n)\right) - e - \frac {K}{2}\, |\boldsymbol{\nabla }n|^2 + T \left(s - \frac {C}{2}\, |\boldsymbol{\nabla }n|^2\right)\\ & = n \mu - e + Ts - nT\, \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\right) - \frac {M}{2}\, |\boldsymbol{\nabla }n|^2\\ & = p - nT\, \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\right) - \frac {M}{2}\, |\boldsymbol{\nabla }n|^2. \end{split} \end{align}

Here, we set $p=n \mu - e + Ts$ . By the definition of $T$ and $\mu$ , we easily have $\boldsymbol{\nabla }s= -({\mu }/{T})\, \boldsymbol{\nabla }n+({1}/{T})\, \boldsymbol{\nabla }e$ . This further leads to

(2.35) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {p}{T} \boldsymbol{I}\right)=\boldsymbol{\nabla }\frac {p}{T}= - e\, \boldsymbol{\nabla }\frac {1}{T} + n\, \boldsymbol{\nabla }\frac {\mu }{T} . \end{align}

Direct calculations yield

(2.36) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {1}{T} \hat {p} \boldsymbol{I}\right) & = \boldsymbol{\nabla }\frac {\hat {p}}{T}= \boldsymbol{\nabla }\left(\frac {p}{T} - n\, \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\right) - \frac {M}{2 T}\, |\boldsymbol{\nabla }n|^2 \right)\nonumber\\ & = \boldsymbol{\nabla }\frac {p}{T} - n\, \boldsymbol{\nabla }\left ( \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\right)\right ) - \left (\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\right)\right ) \boldsymbol{\nabla }n \nonumber\\ &\quad - \frac {M}{2 T}\, \boldsymbol{\nabla }|\boldsymbol{\nabla }n|^2 - |\boldsymbol{\nabla }n|^2\, \boldsymbol{\nabla }\frac {M}{2 T}\nonumber\\ & = \boldsymbol{\nabla }\frac {p}{T} - n\, \boldsymbol{\nabla }\left ( \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\right)\right ) - \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\, \boldsymbol{\nabla }n\right) - \frac {| \boldsymbol{\nabla }n|^2}{2}\, \boldsymbol{\nabla }\left(\frac {K}{T}+C\right)\nonumber\\ & = n\, \boldsymbol{\nabla }\frac {\mu }{T} - e\, \boldsymbol{\nabla }\frac {1}{T} - n\, \boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\right) - \frac {K}{2}\, | \boldsymbol{\nabla }n|^2\, \boldsymbol{\nabla }\frac {1}{T} - \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\, \boldsymbol{\nabla }n\right)\nonumber\\ & = n\, \boldsymbol{\nabla }\frac {\hat {\mu }}{T} - \hat {e}\, \boldsymbol{\nabla }\frac {1}{T} - \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\, \boldsymbol{\nabla }n\right) . \end{align}

Thus we have

(2.37) \begin{equation} n\, \boldsymbol{\nabla }\frac {\hat {\mu }}{T} - \hat {e}\, \boldsymbol{\nabla }\frac {1}{T}=\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {1}{T} \boldsymbol{\varPi }\right)\!, \end{equation}

where

(2.38) \begin{align} \boldsymbol{\varPi }: = \hat {p} \boldsymbol{I} + M\, \boldsymbol{\nabla }n\, \boldsymbol{\nabla }n. \end{align}

The tensor $\boldsymbol{\varPi }$ is the reversible stress tensor containing the gradient contributions (Onuki Reference Onuki2007). By substituting (2.37) into (2.31), we have

(2.39) \begin{equation} \hat { \boldsymbol{\sigma }} = \boldsymbol{\sigma } - \boldsymbol{\varPi }, \end{equation}

where $ \boldsymbol{\sigma }$ is the standard viscous stress defined in (2.29), and $\boldsymbol{\varPi }$ is given in (2.38).

2.5. The two-phase flow equations

Previously, we have derived the constitutive equations (2.30) and (2.39) together with the boundary conditions (2.32) and (2.33) by using the extended Onsager principle. We substitute them into the momentum conservation equation (2.4) and the energy conservation equation (2.5). By combining all the conservation equations and the boundary conditions, we derive a complete system for the one-component two-phase flows:

(2.40) \begin{align} \smash{\raise-65pt\hbox{$\bigg\{$}}&\frac {\partial n}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(n \boldsymbol{u})= 0 \quad & \hbox{in } \varOmega, \end{align}

(2.41) \begin{align} &m_0 \left (\frac {\partial n \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(n \boldsymbol{uu})\right ) = \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{\sigma } - \boldsymbol{\varPi })\quad & \hbox{in } \varOmega, \end{align}

(2.42) \begin{align} &\frac {\partial e_{\textit{tot}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(e_{\textit{tot}} \boldsymbol{u}) = \boldsymbol{\nabla }\boldsymbol{\cdot }(( \boldsymbol{\sigma } - \boldsymbol{\varPi }) \boldsymbol{u}) +\boldsymbol{\nabla }\boldsymbol{\cdot }(\lambda\, \boldsymbol{\nabla }T)\quad & \hbox{in } \varOmega, \end{align}

(2.43) \begin{align} &\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nu }=0,\quad \alpha \boldsymbol{u}_{\tau } + \boldsymbol{P} \boldsymbol{\sigma } \boldsymbol{\nu } + \frac {e_b}{T}\, \boldsymbol{\nabla} _{\varGamma } T - L (n)\, \boldsymbol{\nabla} _{\varGamma } n = 0 \quad & \textrm {on}\ \varGamma _S, \end{align}

(2.44) \begin{align} &\beta \left(\frac {\partial n}{\partial t} + \boldsymbol{u}_\tau \boldsymbol{\cdot }\boldsymbol{\nabla} _{\varGamma } n\right) + L (n) = 0 \quad & \textrm {on} \ \varGamma _S, \end{align}

(2.45) \begin{align} &e_b^{\prime}(n)\, \frac {\partial n}{\partial t} + \boldsymbol{\nabla} _{\varGamma } \boldsymbol{\cdot }(e_b \boldsymbol{u}_{\tau }) = -\lambda\, \boldsymbol{\nabla }T \boldsymbol{\cdot }\boldsymbol{\nu } - \boldsymbol{j}_{\!e, s} \boldsymbol{\cdot }\boldsymbol{\nu } \quad & \textrm {on}\ \varGamma _S, \end{align}

where

(2.46) \begin{align} &e_{\textit{tot}} = e + \frac {K}{2}\, |\boldsymbol{\nabla }n|^2+\frac {m_0 n }{2}\,|\boldsymbol{u}|^2, \quad \boldsymbol{\sigma } = \left ( \eta \left(\boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u}^{\rm T}\right) + \left(\xi - \frac {2 \eta }{d}\right) \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{uI} \right )\!,\quad \nonumber \\ & \boldsymbol{\varPi } = \hat {p} \boldsymbol{I} + M\, \boldsymbol{\nabla }n\, \boldsymbol{\nabla }n,\quad \frac {1}{T} = \frac {\partial s(n,e)}{\partial e},\quad \hat {p} = n \hat {\mu } - \hat {e} + T \hat {s}, \nonumber \\ & \hat {\mu } = - T\left ( \frac {\partial s}{\partial n} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {M}{T}\, \boldsymbol{\nabla }n\right)\right )\!, \quad L (n) = M \frac {\partial n}{\partial \boldsymbol{\nu }} + \frac {\partial f_s (n,T)}{\partial n}. \end{align}

Let $n(x,t)$ , $\boldsymbol{u}$ and $e(x,t)$ be the solutions of the system (2.40)–(2.45) for given initial values. If we substitute (2.30)–(2.33) into (2.26), integration by parts leads to

(2.47) \begin{align} \frac {\textrm{d} \mathcal{S}_{\textit{prod}}}{\textrm{d} t}&=2\varPhi (\boldsymbol{u}, \boldsymbol{j}_{\!e}, \boldsymbol{u}_\tau |_{\varGamma _S},\dot {n}|_{\varGamma _S})\geqslant 0. \end{align}

The term $\varPhi$ is non-negative by definition.

Instead of (2.42), we can also derive the equation for $\hat {e}=e + ({K}/{2})\, |\boldsymbol{\nabla }n|^2$ , i.e.

(2.48) \begin{equation} \frac {\partial \hat {e}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\hat {e} \boldsymbol{u}) = ( \boldsymbol{\sigma } - \boldsymbol{\varPi }) : \boldsymbol{\nabla }\boldsymbol{u} +\boldsymbol{\nabla }\boldsymbol{\cdot }(\lambda\, \boldsymbol{\nabla }T). \end{equation}

Equation (2.42) can also be replaced by an equation for entropy,

(2.49) \begin{equation} \frac {\partial \hat {s}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\hat {s} \boldsymbol{u}) = - \boldsymbol{\nabla }\boldsymbol{\cdot }\widehat {\boldsymbol{J}} + \frac {1}{T} \boldsymbol{\sigma } : \boldsymbol{\nabla }\boldsymbol{u} + \frac {\boldsymbol{j}_{\!e}}{T^2} \boldsymbol{\cdot }\boldsymbol{\nabla }T, \end{equation}

where $\widehat {\boldsymbol{J}} = [M\, \boldsymbol{\nabla }n (\partial n / \partial t + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }n) - \boldsymbol{j}_{\!e}] / T$ .

When the solid boundary is adiabatic, i.e. $ \boldsymbol{j}_{\!e, s}=0$ , (2.45) actually gives a Neumann-type boundary condition for the temperature $T$ . In many situations, we fix the temperature of the boundary, then (2.45) will be replaced by a Dirichlet boundary condition $T (x, t) |_{\varGamma _S} = T_s (x, t)$ . When $\alpha$ goes to infinity, the boundary condition (2.43) will reduce to the standard no-slip boundary condition $\boldsymbol{u}=0$ on $\varGamma _S$ . When $\beta$ goes to zero, condition (2.44) will reduce to $L(n)=0$ , i.e. $M({\partial n}/{\partial \boldsymbol{\nu }}) =- {\partial f_s}/{\partial n}$ . This corresponds to the Young’s equation for the equilibrium contact angle. In applications, the slip coefficient $\alpha$ can be directly obtained from molecular dynamics simulations (Qian et al. Reference Qian, Wang and Sheng2003; Lācis et al. Reference Lācis, Johansson, Fullana, Hess, Amberg, Bagheri and Zaleski2020, Reference Lācis, Pellegrino, Sundin, Amberg, Zaleski, Hess and Bagheri2022). The phenomenological relaxation parameter $\beta$ can be chosen by optimally matching the flow field with that in molecular dynamics simulations (see more details in Qian et al. Reference Qian, Wang and Sheng2003).

2.6. Reformulation for the system

To facilitate the derivation of the lattice Boltzmann scheme, we reformulate the system (2.40)–(2.45) as follows. First, we will use $\rho =m_0 n$ as the fundamental unknown instead of $n$ . Consequently, the mass equation (2.40) is reformulated to

(2.50) \begin{align} \frac {\partial \rho }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u})= 0. \end{align}

Notice the relation $\boldsymbol{\nabla }\boldsymbol{\cdot }(({1}/{T}) (\hat {p} \boldsymbol{I} + M\, \boldsymbol{\nabla }n\, \boldsymbol{\nabla }n)) = n\, \boldsymbol{\nabla } ({\hat {\mu }}/{T}) - \hat {e}\, \boldsymbol{\nabla }({1}/{T})$ . By the definition $\hat {p} = n\hat {\mu } - \hat {e} + T \hat {s}$ , direct calculations lead to the relation

(2.51) \begin{align} \boldsymbol{\nabla }\hat {p} + \boldsymbol{\nabla }\boldsymbol{\cdot }M (\boldsymbol{\nabla }n\, \boldsymbol{\nabla }n)=n\, \boldsymbol{\nabla }\hat {\mu } + \hat {s}\, \boldsymbol{\nabla }T +\frac {1}{T}M(\boldsymbol{\nabla }n\boldsymbol{\cdot }\boldsymbol{\nabla }T)\,\boldsymbol{\nabla }T. \end{align}

Using this equation, the momentum (2.4) can be regrouped as

(2.52) \begin{equation} \frac {\partial \rho \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{uu}) = - n\, \boldsymbol{\nabla }\hat {\mu } - \hat {s}\, \boldsymbol{\nabla }T -\frac {1}{T}M(\boldsymbol{\nabla }n\boldsymbol{\cdot }\boldsymbol{\nabla }T)\,\boldsymbol{\nabla }T + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }. \end{equation}

In addition, the total energy transport (2.42) can be rewritten as

(2.53) \begin{align} &\frac {\partial e_{\textit{tot}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(e_{\textit{tot}} \boldsymbol{u}) = - \boldsymbol{\nabla }\boldsymbol{\cdot }((\hat {p} \boldsymbol{I} + M\,\boldsymbol{\nabla }n\, \boldsymbol{\nabla }n) \boldsymbol{\cdot }\boldsymbol{u}) +\boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{\sigma } \boldsymbol{\cdot }\boldsymbol{u}) + \boldsymbol{\nabla }\boldsymbol{\cdot }(\lambda\, \boldsymbol{\nabla }T)\nonumber \\ &\quad = -(\hat {p} \boldsymbol{I} + M\,\boldsymbol{\nabla }n\, \boldsymbol{\nabla }n) : \boldsymbol{\nabla }\boldsymbol{u} - \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\cdot }(\hat {p} \boldsymbol{I} + M\,\boldsymbol{\nabla }n\, \boldsymbol{\nabla }n) + \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{\sigma } \boldsymbol{\cdot }\boldsymbol{u}) + \boldsymbol{\nabla }\boldsymbol{\cdot }(\lambda\, \boldsymbol{\nabla }T)\nonumber \\ &\quad = -(\hat {p} \boldsymbol{I} + M\,\boldsymbol{\nabla }n\, \boldsymbol{\nabla }n) : \boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{\sigma } \boldsymbol{\cdot }\boldsymbol{u}) + \boldsymbol{\nabla }\boldsymbol{\cdot }(\lambda\, \boldsymbol{\nabla }T)\nonumber \\ &\qquad {}- \boldsymbol{u} \boldsymbol{\cdot }(n\, \boldsymbol{\nabla }\hat {\mu } + \hat {s}\, \boldsymbol{\nabla }T+\frac {1}{T}M(\boldsymbol{\nabla }n\boldsymbol{\cdot }\boldsymbol{\nabla }T)\,\boldsymbol{\nabla }T). \end{align}

The boundary conditions can also be updated accordingly.

2.7. Discussions

To conclude this section, we remark on the connections between the variational approach and classical methods in the literature – e.g. Guo & Lin (Reference Guo and Lin2015), Magaletti et al. (Reference Magaletti, Gallo, Marino and Casciola2016), Kou & Sun (Reference Kou and Sun2018a ) and Gallo et al. (Reference Gallo, Magaletti, Cocco and Casciola2020) for non-isothermal systems, and Ren & E (Reference Ren and E2011) and Magaletti et al. (Reference Magaletti, Gallo, Perez, Carrillo and Kalliadasis2022) for isothermal cases. The classical approach begins with the conservation equations and assumes that the system is in a linear response regime. Constitutive relations are then derived by analysing the entropy production rate and invoking the Clausius–Duhem inequality. In principle, these approaches can yield identical systems of equations, provided that all reversible and irreversible terms are properly specified. While the classical approach is well suited for isothermal systems or those with simple thermomechanical coupling, it requires the careful application of various thermodynamic relations for more complex scenarios (Onuki Reference Onuki2007). In contrast, the extended Onsager principle provides a unified framework for handling systems with diverse effects – such as wetting, phase transition and cross-diffusion in multi-component fluids. This approach starts from the definition of the Rayleigh functional, from which the constitutive and thermodynamic relations are derived directly through variational calculus. The resulting system automatically satisfies the Onsager reciprocal relations and the second law of thermodynamics.

We now distinguish the present model in (2.40)–(2.45) from other models in the literature. A key difference lies in the formulation of the heat flux. In our model, the heat flux is proportional to the temperature gradient, with the effect of the density gradient implicitly handled by the coupled system. In contrast, the models proposed by Magaletti et al. (Reference Magaletti, Gallo, Marino and Casciola2016) and Gallo et al. (Reference Gallo, Magaletti, Cocco and Casciola2020) explicitly incorporate a mass gradient term in the heat flux. Furthermore, our model differs from the quasi-incompressible models in Guo & Lin (Reference Guo and Lin2015) and Sun et al. (Reference Sun, Wu, Jiang, Wise and Guo2024). Those works employ a phase-field function within a Cahn–Hilliard framework to describe phase separation. Our model, however, is developed within a fully compressible framework and is capable of describing both phase separation and phase transitions in a single-component, two-phase flow system.

3.1. The well-balanced discrete velocity kinetic model for velocity distribution function

Building upon the established kinetic model for decoupled hydrodynamic and energy equations in Guo et al. (Reference Guo, Zheng, Shi and Zhao2007), and inspired by the well-balanced LBM for two-phase flow (Guo Reference Guo2021), we develop a discrete velocity kinetic model for non-isothermal two-phase flow through appropriate modifications to the source distribution functions, as detailed below.

The evolution equations of the DVM for velocity distribution function read

(3.1) \begin{equation} \partial _t f_i + \boldsymbol{c}_i \boldsymbol{\cdot }\boldsymbol{\nabla }f_i = -\frac {1}{\tau _f}\left(f_i-f_i^{\textit{eq}}\right) + F_i, \end{equation}

where $f_i$ is the velocity distribution function after the discretisation of the molecule velocity space, and $f_i^{\textit{eq}}$ is the modified local equilibrium distribution function, defined as

(3.2) \begin{equation} f_i^{\textit{eq}} = {\begin{cases} {\rho - (1 - {\omega _0})\rho _0/(R T_0) + {\omega _0}{s_0}(\boldsymbol{u})}, &i = 0,\\ {{\omega _i}\rho _0/(R T_0) + {s_i}(\boldsymbol{u}),}&i \ne 0. \end{cases}} \end{equation}

Here, $s_i(\boldsymbol{u})=\omega _i (({\boldsymbol{c}_i \boldsymbol{\cdot }(\rho \boldsymbol{u})}/{R T_0}) + ({(\rho \boldsymbol{uu}):(\boldsymbol{c}_i \boldsymbol{c}_i -c_s^2 \boldsymbol{I})}/{2 (R T_0)^2}) )$ , $R$ is the gas constant, $T_0$ is the reference temperature, $\rho _0$ is a numerical constant, and the weights $\omega _i$ can be found in the Appendix of Guo (Reference Guo2021). In our simulation, the value of $\rho _0$ is set to 0 for simplicity. Also, $F_i$ is the force term defined as

(3.3) \begin{equation} F_i = \omega _i \left (\frac {\boldsymbol{c}_i \boldsymbol{\cdot }\boldsymbol{F}_{\textit{non}}}{R T_0} + \frac {\boldsymbol{u} (\boldsymbol{F}_{\textit{non}} + c_s^2 \,\boldsymbol{\nabla }\rho ):(\boldsymbol{c}_i \boldsymbol{c}_i - R T_0 \boldsymbol{I})}{(R T_0)^2} + \frac {1}{2} \left(\frac {|\boldsymbol{c}_i|^2}{R T_0}-d\right)(\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\rho )\right )\!, \end{equation}

where $\boldsymbol{F}_{\textit{non}} = - n\, \boldsymbol{\nabla }\hat {\mu } - \hat {s}\, \boldsymbol{\nabla }T - ({1}/{T})M(\boldsymbol{\nabla }n\boldsymbol{\cdot }\boldsymbol{\nabla }T)\,\boldsymbol{\nabla }T$ .

By using the isotropic constraints, the following moments can be obtained:

(3.4) \begin{equation} \left \{\!\! \begin{array}{l} \sum _i f_{i}^{\textit{eq}} = \sum _i f_{i} = \rho ,\\[5pt] \sum _i \boldsymbol{c}_i f_{i}^{\textit{eq}} = \sum _i \boldsymbol{c}_i f_{i} = \rho \boldsymbol{u},\\[5pt] \sum _i \boldsymbol{c}_i \boldsymbol{c}_i f_{i}^{\textit{eq}} = \rho _0 \boldsymbol{I} + \rho \boldsymbol{uu},\\[5pt] \sum _i \boldsymbol{c}_i \boldsymbol{c}_i \boldsymbol{c}_i f_{i}^{\textit{eq}} = R T_0 \rho [\boldsymbol{uI}],\\[5pt] \sum _i F_{i}(\boldsymbol{x},t) = 0,\\[5pt] \sum _i \boldsymbol{c}_i F_{i} (\boldsymbol{x},t) = \boldsymbol{F}_{\textit{non}},\\[5pt] \sum _i \boldsymbol{c}_i \boldsymbol{c}_i F_{i}= [\boldsymbol{uF}_{\textit{non}}] + R T_0 [\boldsymbol{u}\, \boldsymbol{\nabla }\rho ] + R T_0 \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\rho \boldsymbol{I}, \end{array} \right . \end{equation}

where $[\boldsymbol{uI}]_{\alpha \beta \gamma } = u_\alpha \delta _{\beta \gamma } + u_\beta \delta _{\alpha \gamma } + u_\gamma \delta _{\alpha \beta }$ , $[\boldsymbol{uF}]_{\alpha \beta } = u_{\alpha } F_{\beta } + u_{\beta } F_{\alpha }$ and $[\boldsymbol{u}\, \boldsymbol{\nabla }\rho ]_{\alpha \beta } = u_\alpha\, \partial _\beta \rho + u_\beta\, \partial _\alpha \rho$ .

3.2. Chapman–Enskog analysis of the well-balanced discrete velocity kinetic model

We will derive the equations satisfied by the moments (3.4) by Chapman–Enskog analysis. We introduce the following Chapman–Enskog expansions for the evolution (3.1):

(3.5) \begin{align} \partial _{t} &= \varepsilon\, \partial _{t_1} + \varepsilon ^2\, \partial _{t_2} + \cdots ,\quad \boldsymbol{\nabla }= \varepsilon \boldsymbol{\nabla} _1,\quad F_i = \varepsilon F_i^{(1)}, \end{align}

(3.6) \begin{align} f_i &= f_i^{0} + \varepsilon f_i^{(1)} + \varepsilon ^2 f_i^{(2)} + \cdots . \end{align}

Here, $\varepsilon$ is a small parameter proportional to the Knudsen number, $t_1=\varepsilon t$ , $t_2=\varepsilon ^2 t$ and $x_1=\varepsilon x$ are corresponding macroscopic scales, and $\boldsymbol{\nabla} _1 = (\partial _{x_1}, \partial _{y_1}, \partial _{z_1})$ .

Substituting the expansions into (3.1) yields the following equations at different scales:

(3.7) \begin{align}&\,\,\, \text{for}\ O(\varepsilon ),\quad \partial _{t_1} f_{i}^{(0)} + \boldsymbol{c}_i \boldsymbol{\cdot }\boldsymbol{\nabla} _1 f_{i}^{(0)} = - \frac {1}{ \tau _f} f_{i}^{(1)} + F_{i}^{(1)}, \end{align}

(3.8) \begin{align}& \text{for}\ O(\varepsilon ^2),\quad \partial _{t_2} f_{i}^{(0)} + \partial _{t_1} f_{i}^{(1)} + \boldsymbol{c}_i \boldsymbol{\cdot }\boldsymbol{\nabla} _1 f_{i}^{(1)} = - \frac {1}{\tau _f} f_{i}^{(2)}. \end{align}

Before recovering the macroscopic equation, the following equations are derived from the first two equations of the moment system (3.4):

(3.9) \begin{align} \sum _i f_i^{(0)}=\rho ,\quad \sum _i c_i f_i^{(0)}=\rho \boldsymbol{u}. \end{align}

By taking the zeroth- and first-order velocity moments of (3.7), we obtain the following hydrodynamic equations at $O(\varepsilon )$ :

(3.10) \begin{align} \frac {\partial \rho }{\partial t_1} + \boldsymbol{\nabla} _1 \boldsymbol{\cdot }(\rho \boldsymbol{u})&=0, \end{align}

(3.11) \begin{align} \frac {\partial (\rho \boldsymbol{u})}{\partial t_1} + \boldsymbol{\nabla} _1 \boldsymbol{\cdot }(\rho _0 \boldsymbol{I} + \rho \boldsymbol{u}\boldsymbol{u})&= \boldsymbol{F}_{\textit{non}}^{(1)}. \end{align}

Similarly, by taking the zeroth- and first-order velocity moments of (3.8), we get the hydrodynamic equations at $O(\varepsilon ^2)$ :

(3.12) \begin{align} \frac {\partial \rho }{\partial t_2} = 0, \end{align}

(3.13) \begin{align} \frac {\partial (\rho \boldsymbol{u})}{\partial t_2} + \boldsymbol{\nabla} _1 \boldsymbol{\cdot }\sum _i \boldsymbol{c}_i \boldsymbol{c}_i f_i^{(1)}= 0. \end{align}

From (3.7), we have

(3.14) \begin{align} -\frac {1}{\tau _f} \sum _i \boldsymbol{c}_i \boldsymbol{c}_i f_i^{(1)} & = \partial _{t_1} (\rho _0 \boldsymbol{I} + \rho \boldsymbol{uu}) + \boldsymbol{\nabla} _1 \boldsymbol{\cdot }(R T_0\rho [\boldsymbol{uI}]) - \left[\boldsymbol{uF}_{\textit{non}}^{(1)}\right] \nonumber \\ &\quad {}- R T_0 [\boldsymbol{u}\, \boldsymbol{\nabla} _1 \rho ] - R T_0 \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla} _1 \rho \boldsymbol{I}\nonumber \\ &=R T_0 \rho (\boldsymbol{\nabla} _1 \boldsymbol{u} + \boldsymbol{\nabla} _1 \boldsymbol{u}^{\rm T} + \boldsymbol{\nabla} _1 \boldsymbol{\cdot }(\boldsymbol{uI})) + O\left(\boldsymbol{u}^3\right), \end{align}

where we have used the expressions

(3.15) \begin{align} \partial _{t_1}(\rho u_\alpha u_\beta ) &= \rho u_\alpha\, \partial _{t_1} u_\beta + u_\beta\, \partial _{t_1} (\rho u_\alpha )\nonumber\\ &= u_\alpha F_{non \beta }^{(1)} - u_\alpha u_\beta\, \partial _{t_1} \rho - u_\alpha\, \boldsymbol{\nabla} _{1\beta } ( \rho u_\beta u_\gamma ) + u_\beta F_{non \alpha }^{(1)} - u_\beta \boldsymbol{\nabla} _{1\alpha } (\rho u_\alpha u_\gamma )\nonumber\\ &= \left[\boldsymbol{uF}_{\textit{non}}^{(1)}\right] + u_\alpha u_\beta \boldsymbol{\nabla} _{1\gamma } \rho u_\gamma - u_\alpha \boldsymbol{\nabla} _{1\beta } (\rho u_\beta u_\gamma )- u_\beta\, \boldsymbol{\nabla} _{1\alpha } (\rho u_\alpha u_\gamma )\nonumber\\ &= \left[\boldsymbol{uF}_{\textit{non}}^{(1)}\right] - \boldsymbol{\nabla} _1 \boldsymbol{\cdot }(\rho \boldsymbol{uuu}) \end{align}

and

(3.16) \begin{align} \begin{split} \boldsymbol{\nabla} _1 \boldsymbol{\cdot }(R T_0 \rho [\boldsymbol{uI}]) & = \boldsymbol{\nabla} _{1\alpha } (R T_0 \rho (u_\alpha \delta _{\beta \gamma } + u_\beta \delta _{\alpha \gamma } + u_\gamma \delta _{\alpha \beta }))\\ &= \boldsymbol{\nabla} _{1\alpha }(R T_0 \rho u_\alpha \delta _{\beta \gamma }) + \boldsymbol{\nabla} _{1\gamma } (R T_0 \rho u_\beta ) + \boldsymbol{\nabla} _{1\beta } (R T_0 \rho u_\gamma ). \end{split} \end{align}

Substituting (3.14) into (3.13) and enforcing $\eta = R T_0\tau _f \rho$ , we obtain

(3.17) \begin{align} \text{for}\ O(\varepsilon ^2 ),\quad \frac {\partial (\rho \boldsymbol{u})}{\partial t_2} = \boldsymbol{\nabla} _1 \boldsymbol{\cdot }\eta ( \boldsymbol{\nabla} _1 \boldsymbol{u} + \boldsymbol{\nabla} _1 \boldsymbol{u}^{\rm T} + \boldsymbol{\nabla} _1 \boldsymbol{\cdot }(\boldsymbol{uI})) +O(\boldsymbol{u}^3) . \end{align}

Combining the equations at both orders $O(\varepsilon )$ and $O(\varepsilon ^2)$ , the following hydrodynamic equations are recovered:

(3.18) \begin{align} \frac {\partial \rho }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u}) &= 0, \end{align}

(3.19) \begin{align} \frac {\partial \rho \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{uu}) &= \boldsymbol{\nabla }\boldsymbol{\cdot }\eta ( \boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u}^{\rm T} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{uI})) + \boldsymbol{F}_{\textit{non}}. \end{align}

Here, the term $O(\boldsymbol{u}^3)$ is not Galilean invariant and is therefore omitted for flows with small Mach numbers.

3.3. The discrete velocity kinetic model for the total energy distribution function

The evolution equations of the DVM for total energy distribution function $h_i$ read

(3.20) \begin{equation} \partial _t h_i + \boldsymbol{c}_i \boldsymbol{\cdot }\boldsymbol{\nabla }h_i = -\frac {1}{\tau _h}(h_i - h_i^{\textit{eq}}) + \frac {Z_i}{ \tau _{h}}(f_i-f_i^{\textit{eq}}) + q_i. \end{equation}

Here, $Z_i=\boldsymbol{c}_i\boldsymbol{\cdot }\boldsymbol{u} - |\boldsymbol{u}|^2/2$ , the equilibrium distribution function $h_i^{\textit{eq}}$ is defined as

(3.21) \begin{equation} h_i^{\textit{eq}} = e_{\textit{tot}} f_i^{\textit{eq}} + \omega _i \frac {T}{2}\left(\frac {c_i^2}{(R T_0)^2} - D\right)\!, \end{equation}

and $q_i$ is the source term distribution function defined as

(3.22) \begin{equation} q_i=\omega _i e_{\textit{tot}} \frac {\boldsymbol{c}_i\boldsymbol{\cdot }\boldsymbol{F}_{\textit{non}}}{(R T_0)^2} + f_i \boldsymbol{c}_i \boldsymbol{\cdot }\hat {\boldsymbol{F}}, \end{equation}

where $\hat {\boldsymbol{F}} = \tilde {\boldsymbol{F}} + \boldsymbol{F}_{\textit{non}}/\rho$ , with $\tilde {\boldsymbol{F}} = (((-\hat {p} - M n_x^2)u_x - M n_x n_y v_x)/\rho u,(-M n_x n_y u_y + v_y(-\hat {p} - M n_y^2))/\rho v)^{\rm T}$ .

Similarly, $h_i$ and $q_i$ satisfy the following moments conditions:

(3.23) \begin{equation} \left \{\!\! \begin{array}{l} \sum _i h_{i}^{\textit{eq}} = \sum _i h_{i} = e_{\textit{tot}},\\[5pt] \sum _i \boldsymbol{c}_i h_{i}^{\textit{eq}} = e_{\textit{tot}} \boldsymbol{u},\\[5pt] \sum _i \boldsymbol{c}_i \boldsymbol{c}_i h_{i}^{\textit{eq}} = (R T_0)^2 T \boldsymbol{I} + e_{\textit{tot}}\boldsymbol{uu},\\[5pt] \sum _i q_{i} = \rho \boldsymbol{u} \boldsymbol{\cdot }\hat {\boldsymbol{F}},\\[5pt] \sum _i \boldsymbol{c}_i q_{i}= e_{\textit{tot}} \boldsymbol{F}_{\textit{non}} + \sum _i \boldsymbol{c}_i \boldsymbol{c}_i f_i \boldsymbol{\cdot }\hat {\boldsymbol{F}},\\[5pt] \sum _i \boldsymbol{c}_i \boldsymbol{c}_i q_{i}= 0. \end{array} \right . \end{equation}

3.4. Chapman–Enskog analysis of the discrete velocity kinetic model for the total energy distribution function

As in § 3.2, applying the Chapman–Enskog expansion to the evolution (3.20) yields the following equations at different scales:

(3.24) \begin{align} \text{for}\ O(\varepsilon ),\quad \partial _{t_1} h_{i}^{(0)} + \boldsymbol{c}_i \boldsymbol{\cdot }\boldsymbol{\nabla} _1 h_{i}^{(0)} &= - \frac {1}{ \tau _h} h_{i}^{(1)} + \frac {Z_k}{\tau _{h}} f_i^{(1)} + q_{i}^{(1)}, \end{align}

(3.25) \begin{align} \text{for}\ O(\varepsilon ^2),\quad \partial _{t_2} h_{i}^{(0)} + \partial _{t_1} h_{i}^{(1)} + \boldsymbol{c}_i \boldsymbol{\cdot }\boldsymbol{\nabla} _1 h_{i}^{(1)} &= - \frac {1}{\tau _h} h_{i}^{(2)} + \frac {Z_k}{\tau _{h}} f_i^{(2)} + q_i^{(2)}. \end{align}

Similarly, the following equations can be derived from (3.23):

(3.26) \begin{align} \sum _i h_{i}^{(m)} = 0,\quad m \geqslant 1. \end{align}

By taking the zeroth-order moments of (3.24) and (3.25), we obtain the following macroscopic equations at orders $O(\varepsilon )$ and $O(\varepsilon ^2)$ :

(3.27) \begin{align} \text{for}\ O(\varepsilon ),\quad \frac {\partial e_{\textit{tot}}}{\partial t_1} + \boldsymbol{\nabla} _1 \boldsymbol{\cdot }(e_{\textit{tot}} \boldsymbol{u})&= \rho \boldsymbol{u} \boldsymbol{\cdot }\hat{\boldsymbol{F}}^{(1)}, \end{align}

(3.28) \begin{align} \text{for}\ O(\varepsilon ^2 ),\quad \frac {\partial e_{\textit{tot}}}{\partial t_2} + \boldsymbol{\nabla} _1 \boldsymbol{\cdot }\sum _i \boldsymbol{c}_i h_i^{(1)}&= 0. \end{align}

From (3.24), we obtain

(3.29) \begin{align} -\frac {1}{\tau _h} \sum _i \boldsymbol{c}_i h_i^{(1)} & = \partial _{t_1} \sum _i \boldsymbol{c}_i h_i^{\textit{eq}} + \boldsymbol{\nabla} _1 \boldsymbol{\cdot }\sum _i \boldsymbol{c}_i \boldsymbol{c}_i h_i^{\textit{eq}} - \frac {1}{\tau _{h}} \sum _i \boldsymbol{c}_i \boldsymbol{c}_i f_i^{(1)}\boldsymbol{\cdot }\boldsymbol{u} - \sum _i \boldsymbol{c}_i q_i^{(1)}\nonumber \\ &=\partial _{t_1} (e_{\textit{tot}} \boldsymbol{u}) + \boldsymbol{\nabla} _1 \boldsymbol{\cdot }(e_{\textit{tot}} \boldsymbol{uu} + (R T_0)^2 T \boldsymbol{I}) + \frac {1}{\tau _h} \boldsymbol{\sigma }^{(1)} \boldsymbol{\cdot }\boldsymbol{u} \nonumber \\ &\quad {}- \rho \boldsymbol{uu} \boldsymbol{\cdot }\hat {\boldsymbol{F}}^{(1)} - e_{\textit{tot}} \boldsymbol{F}_{\textit{non}}^{(1)}. \end{align}

Notice that

(3.30) \begin{align} \partial _{t_1}(e_{\textit{tot}} \boldsymbol{u}) & = \boldsymbol{u}\, \partial _{t_1} e_{\textit{tot}} + e_{\textit{tot}}\, \partial _{t_1} \boldsymbol{u}\nonumber\\ &=\rho \boldsymbol{uu} \boldsymbol{\cdot }\hat{\boldsymbol {F}}^{(1)} - \boldsymbol{u}\,\boldsymbol{\nabla} _1 \boldsymbol{\cdot }(e_{\textit{tot}} \boldsymbol{u}) + e_{\textit{tot}}\, \partial _{t_1} \boldsymbol{u}\nonumber\\ &=\rho \boldsymbol{uu} \boldsymbol{\cdot }\hat{\boldsymbol {F}}^{(1)} - \boldsymbol{u}\,\boldsymbol{\nabla} _1 \boldsymbol{\cdot }(e_{\textit{tot}} \boldsymbol{u}) + e_{\textit{tot}} \boldsymbol{F}_{\textit{non}}^{(1)} - e_{\textit{tot}} \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla} _1 \boldsymbol{u}\nonumber\\ &=\rho \boldsymbol{uu} \boldsymbol{\cdot }\hat{\boldsymbol {F}}^{(1)} + e_{\textit{tot}} \boldsymbol{F}_{\textit{non}}^{(1)} - \boldsymbol{\nabla} _1 \boldsymbol{\cdot }(e_{\textit{tot}} \boldsymbol{uu}). \end{align}

Substituting this equation into (3.29) yields

(3.31) \begin{equation} -\frac {1}{\tau _h} \sum _i \boldsymbol{c}_i h_i^{(1)} = \boldsymbol{\nabla} _1 \boldsymbol{\cdot }(R T_0)^2 T \boldsymbol{I} + \frac {1}{\tau _h} \boldsymbol{\sigma }^{(1)} \boldsymbol{\cdot }\boldsymbol{u}. \end{equation}

By combining the equations at both orders $O(\varepsilon )$ and $O(\varepsilon ^2)$ , and imposing $\lambda = (R T_0)^2\tau _h$ , the following total energy transport equation is recovered:

(3.32) \begin{equation} \frac {\partial e_{\textit{tot}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(e_{\textit{tot}} \boldsymbol{u}) = \boldsymbol{\nabla }\boldsymbol{\cdot }\lambda\, \boldsymbol{\nabla }T + \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{\sigma } \boldsymbol{\cdot }\boldsymbol{u}) + \rho \boldsymbol{u} \boldsymbol{\cdot }\hat{\boldsymbol {F}}. \end{equation}

3.5. Lattice Boltzmann scheme

Building on the discrete velocity kinetic models established in the previous subsections, we develop a lattice Boltzmann scheme for the non-isothermal two-phase fluid system through characteristic discretisation and trapezoidal quadrature, leading to the following implicit scheme:

(3.33) \begin{align} f_i(\boldsymbol{x} + \boldsymbol{c}_i\, \delta t, t+\delta t) - f_i(\boldsymbol{x},t) ={} &\frac {\delta t}{2} [\varOmega _f(x+\boldsymbol{c}_i\, \delta t, t+\delta t) + F_i(\boldsymbol{x} + \boldsymbol{c}_i\, \delta t, t+ \delta t)] \nonumber \\ &{}+ \frac {\delta t}{2}[\varOmega _f(x, t) + F_i(\boldsymbol{x} , t)], \end{align}

(3.34) \begin{align} h_i(\boldsymbol{x} + \boldsymbol{c}_i\, \delta t, t+\delta t) - h_i(\boldsymbol{x},t) ={} & \frac {\delta t}{2} [\varOmega _h(x+\boldsymbol{c}_i\, \delta t, t+\delta t) + S_i(\boldsymbol{x} + \boldsymbol{c}_i\, \delta t, t+ \delta t)] \nonumber \\ &{}+\frac {\delta t}{2}[\varOmega _h(x, t) + S_i(\boldsymbol{x} , t)], \end{align}

where $\delta t$ is the time step, $\varOmega _f = (f_i^{\textit{eq}} - f_i)/\tau _f$ , $\varOmega _h = (h_i^{\textit{eq}} - h_i)/\tau _h$ and $S_i=-Z_i \tau _f \varOmega _f/\tau _h + q_i$ .

The implicit nature of the scheme is removed by defining the distribution functions

(3.35) \begin{align} \hat {f}_i &= f_i - \frac {\delta t}{2}(\varOmega _f + F_i), \end{align}

(3.36) \begin{align} \hat {h}_i &= h_i - \frac {\delta t}{2}(\varOmega _h + S_i). \end{align}

By substituting (3.35) and (3.36) into the discrete evolution equations (3.33) and (3.34), we obtain

(3.37) \begin{align}& \hat {f}_i(\boldsymbol{x} + \boldsymbol{c}_i\, \delta t, t+\delta t) - \hat {f}_i(\boldsymbol{x},t) = - \omega _f[\hat {f}_i(\boldsymbol{x},t) - f_i^{\textit{eq}}(\boldsymbol{x},t)] + \delta t\left(1-\frac {\omega _f}{2}\right) F_i, \end{align}

(3.38) \begin{align} &\hat {h}_i(\boldsymbol{x} + \boldsymbol{c}_i\, \delta t, t+\delta t) - \hat {h}_i(\boldsymbol{x},t) \nonumber\\ &\quad= -\omega _h\!\left(\hat {h}_i(\boldsymbol{x},t) - h_i^{\textit{eq}}(\boldsymbol{x},t)\!\right) + \delta t\left(\!1-\frac {\omega _h}{2}\!\right)\!\left(q_i + \frac {Z_i}{\tau _h}\left(f_i - f_i^{\textit{eq}}\right)\!\right)\nonumber\\ &\quad = -\omega _h\!\left(\hat {h}_i(\boldsymbol{x},t) - h_i^{\textit{eq}}(\boldsymbol{x},t)\!\right) + \delta t\!\left(\!1-\frac {\omega _h}{2}\!\right)\!q_i + \omega _h Z_i\!\left(\!1+\frac {\delta t}{2\tau _f}\!\right)^{-1}\!\!\left(\!\hat {f}_i - f_i^{\textit{eq}} + \frac {\delta t}{2} F_i\!\!\right) \nonumber\\ &\quad = -\omega _h\left(\hat {h}_i(\boldsymbol{x},t) - h_i^{\textit{eq}}(\boldsymbol{x},t)\!\right) + \delta t\!\left(\!1-\frac {\omega _h}{2}\!\right)q_i + \omega _h\!\left(\!1-\frac {\omega _f}{2}\!\right)Z_i \!\left(\!\hat {f}_i - f_i^{\textit{eq}} + \frac {\delta t}{2} F_i\!\right)\!, \end{align}

where $\omega _f=2\,\delta t/(2\tau _f + \delta t)$ and $\omega _h=2\,\delta t/(2\tau _h + \delta t)$ . The fluid density, velocity and total energy can be computed from the new distribution function as

(3.39) \begin{equation} \rho =\sum _i \hat {f}_i,\quad \rho \boldsymbol{u} = \sum _i \boldsymbol{c}_i \hat {f}_i + \frac {\delta t}{2} \boldsymbol{F}_{\textit{non}} \end{equation}

and

(3.40) \begin{equation} e_{\textit{tot}} = \sum _i \hat {h}_i + \frac {\delta t}{2} \rho \boldsymbol{u} \boldsymbol{\cdot }\hat{\boldsymbol {F}}. \end{equation}

In a standard lattice Boltzmann scheme, the computational domain is discretised into a uniform lattice. For example, in the D2Q9 model, the lattice spacing is set to $\delta x = c\, \delta t$ , with $c=\sqrt {3R T_0}$ . Consequently, the lattice size is closely dependent on the discrete velocity set, meaning that the spatial discretisation is coupled with the velocity discretisation.

3.5.1 Boundary condition treatment

Treating boundary conditions in the LBM requires a method to transform conditions on fluid variables into equivalent conditions for the distribution functions. This study employs dynamic boundary conditions that incorporate thermal and wetting effects. Following a boundary treatment method similar to that in Huang & Wang (Reference Huang and Wang2018), the boundary values for number density $n$ , velocity $\boldsymbol{u}$ and temperature $T$ are obtained from the dynamic boundary conditions in (2.43), (2.44) and (2.45). For a mesh point $(i,0)$ on the bottom wall, the number density at $m$ time step $n^m(i,0)$ is calculated using a finite difference scheme

(3.41) \begin{align} n^m(i,0) ={} & n^{m-1}(i,0) - \frac {\delta t}{\beta }\, u_x^{m-1}(i,0)\left(\frac {3n^{m-1}(i,0) - 4n^{m-1}(i,1) + n^{m-1}(i,2)}{2}\right) \nonumber \\ & {}- \frac {\delta t}{\beta }\, L\left(n^{m-1}(i,0)\right)\!, \end{align}

where $\delta t$ is the time step. The velocity $\textbf{u}^m$ and temperature $T^m$ are updated similarly. With the updated macroscopic variables on the boundary, unknown distribution functions at boundary nodes are calculated by the non-equilibrium extrapolation scheme (Guo, Zheng & Shi Reference Guo, Zheng and Shi2002).

In numerical simulations, it is essential to impose appropriate outflow boundary conditions at open boundaries. This study employs a constant-pressure boundary condition to model the outflow. The unknown density distribution functions are determined by the method described in Huang, Wang & Yang (Reference Huang, Wang and Yang2024).

4.1. Validation of the scaling law

Figure 1.

The density contours of droplet evaporation at different times: $t=0.1$ , $t=0.35$ , $t=0.75$ .

Figure 2.

Time evolution of the square of the dimensionless diameter for the present model.

This subsection presents a numerical study of the $D^2$ law for droplet evaporation and the Cox–Voinov law. We begin by examining droplet evaporation in an open environment. The process follows the $D^2$ law, which states that the square of the droplet diameter $D$ decreases linearly with time, i.e.

(4.1) \begin{equation} \left(\frac {D}{D_0}\right)^2 = 1-\kappa t, \end{equation}

where $D_0$ is the initial droplet diameter, and $\kappa$ is the evaporation constant. In our simulations, the computational domain is a square with grid resolution $200\times 200$ . A droplet with initial diameter $D_0 = 60\, \delta x$ is placed at the centre of the domain. The droplet’s temperature is set to the saturation temperature $T_{s} = 0.86 T_c$ , where $T_c$ stands for the critical temperature. The surrounding vapour is maintained at a higher temperature $T_g=1.0 T_c$ , while, subject to this temperature gradient at the liquid–vapour interface, the droplet undergoes evaporation.

For the temperature field, a Dirichlet condition $T=T_g$ is applied, while periodic boundary conditions are used for the velocity field. Figure 1 shows density contours at different dimensionless times, illustrating the droplet’s gradual evaporation driven by the interfacial temperature gradient. Quantitatively, figure 2 plots the temporal evolution of the squared droplet diameter, showing excellent agreement with the $D^2$ law.

Figure 3.

Numerical validation of the Cox–Voinov law.

Next, we consider a droplet deposited on a smooth substrate with initial contact angle $\theta _0 \gt \theta _{\textit{eq}}$ , which spreads and relaxes towards an equilibrium shape determined by $\theta _{\textit{eq}}$ . The spreading dynamics can be quantitatively analysed using the Cox–Voinov law, which relates the apparent contact angle $\theta _a$ to the velocity $U$ of the moving contact line via

(4.2) \begin{align} \theta _a^3 - \theta _{\textit{eq}}^3 \propto Ca, \end{align}

where the capillary number is defined as $Ca = \mu U /\gamma$ . Here, $\mu$ is the dynamic viscosity, and $\gamma$ is the surface tension. We set $\mu = 0.167$ $\text{N}\ \text{s}\ \text{m}^{-2}$ in our simulations.

The simulations are conducted on a $256\times 256 \times 256$ lattice. Initially, a droplet with contact angle $\theta _0 = 2/3\pi$ is placed at the centre of the substrate, and the equilibrium contact angle is taken as $\theta _{\textit{eq}} = \pi /3$ . To extract the dynamic contact angle from the simulation data, we assume that the droplet is closely approximated by a spherical-cap shape after an initial transition period. Under this assumption and by enforcing volume conservation, the contact angle can be determined numerically. The volume $V$ of a spherical cap is given by

(4.3) \begin{equation} V=\frac {\pi }{3} R^3(2+\cos \theta _a)(1-\cos \theta _a)^2. \end{equation}

Since our numerical scheme conserves mass, the droplet volume $V$ remains constant. By measuring the droplet height $h_d(t)$ and contact line diameter $2\,r(t)$ , we compute the time-dependent radius of the droplet as

(4.4) \begin{equation} R(t) = \frac {r(t)^2 + h_d(t)^2}{2\,h_d(t)}, \end{equation}

then solve (4.3) for the contact angle $\theta (t)$ . Figure 3 shows a plot of $\theta ^3(t) - \theta _{\textit{eq}}^3$ versus $Ca$ from a numerical simulation of a spreading droplet. The linear scaling observed confirms excellent agreement with the Cox–Voinov law.

Figure 4.

Dependence of the apparent contact angle on wall temperature.

In addition, we also investigated the effect of wall temperature on the apparent contact angle $\theta$ . The numerical results are illustrated in figure 4. It can be observed that the apparent contact angle gradually decreases with increasing wall temperature, which is consistent with theoretical predictions (Adamson Reference Adamson1973; Sullivan Reference Sullivan1981; Diaz et al. Reference Diaz, Savage and Cerro2017).

4.2. Leidenfrost droplet impact on a flat plate

Droplet impact on a heated plate above the Leidenfrost temperature is a phenomenon that has attracted considerable research interest. In this subsection, the developed model and LBM are validated through simulations of a droplet impacting a heated plate above the Leidenfrost temperature. The initial droplet radius is set as $ R_0=50$ . The liquid and vapour saturation temperature, denoted as $ T_s$ , is maintained at $ 0.8 T_c$ , leading to the corresponding coexistence densities $ \rho _l=7.2$ and $ \rho _g=0.197$ . Under these conditions, the measured surface tension $ \gamma$ is determined to be 0.11235.

To facilitate a more precise comparison with experimental data, several dimensionless numbers are introduced to align with realistic conditions. The Prandtl number for the liquid phase is defined as

(4.5) \begin{align} {\textit{Pr}}_l = \frac {\mu _l c_{\nu l}}{\lambda _l}, \end{align}

while the Prandtl number for the gas phase is given by

(4.6) \begin{align} {\textit{Pr}}_g = \frac {\mu _g c_{\nu g}}{\lambda _g}, \end{align}

where $\mu _l$ , $\mu _g$ , $c_{\nu l}$ , $c_{\nu g}$ , $\lambda _l$ and $\lambda _g$ are the dynamic viscosities, the specific heat capacities, and the thermal conductivities of the liquid phase and the gas phase, respectively. Additionally, the Jakob number is incorporated as

(4.7) \begin{align} \textit{Ja} = \frac {c_{\nu l} (T_w - T_s)}{h_{fg}}, \end{align}

where $T_w$ is the wall temperature, and $h_{fg}$ is the latent heat, which is calculated by $h_{fg} = h_g-h_l$ , with $h_g$ and $h_l$ the enthalpy values of the gas and liquid phases, respectively. The Weber number and Bond number are also considered to represent the experimental conditions, given by

(4.8) \begin{align} \textit{We} = \frac {D_0 \rho _l U^2}{ \gamma }, \quad Bo = \frac {\rho _l g D_0^2}{ \gamma }, \end{align}

where $U$ is the droplet initial velocity. In our simulations, the selected parameters are $ c_{\nu l} = 7.4$ , $ c_{\nu g} = 3.4$ , $ \lambda _l = 0.25$ and $ \lambda _g = 0.1$ , which result in Prandtl numbers $ {\textit{Pr}}_g=0.94$ and $ {\textit{Pr}}_l=1.49$ , respectively.

For the boundary conditions, the top of the domain is defined as an outflow boundary, while the bottom solid wall is treated as a no-slip boundary maintained at a constant temperature $ T_w$ . Periodic conditions are applied to the lateral boundaries. Given that the contact temperature between the droplet and the heated plate rapidly equilibrates to $ T_w$ upon impact, the temperature of the adjacent fluid layer is also set to $ T_w$ .

The proposed model is validated against the experimental study in Biance et al. (Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006), where a water droplet with initial radius $ R_0 = 1$ mm impacts a flat plate heated to $ 300\,^\circ\rm C$ at initial velocity $0.8\ \text{m}\ \text{s}^{-1}$ . The corresponding dimensionless parameters in this case are set to be ${\textit{Ja}} = 0.50$ , ${\textit{We}} = 21.6$ , $Bo = 0.16$ .

Figure 5.

The experiment snapshots and simulation results for droplet impact on a heated plate.

A comparison of the simulation results with experimental snapshots is presented in figure 5. The numerical results show excellent agreement with experimental observations, confirming that the developed lattice Boltzmann model reliably predicts Leidenfrost droplet impingement dynamics.

For a more quantitative validation, we further examine the experimental study by Wachters et al. (Reference Wachters, Smulders, Vermeulen and Kleiweg1966), where the dimensionless numbers are specified as $ We = 16$ and $ Ja = 0.61$ . The evolution of the droplet’s gravitational centre height, denoted as $ \epsilon _{\textit{mass}}$ , is computed and compared with the corresponding experimental data. The gravitational centre height of the droplet is determined using the expression

(4.9) \begin{equation} \epsilon _{\textit{mass}} = \frac {\sum _{\rho \gt \rho _{\textit{avg}}}\rho (\boldsymbol{x})\,h(\boldsymbol{x})}{\sum _{\rho \gt \rho _{\textit{avg}}}\rho (\boldsymbol{x})}, \end{equation}

where $ h(\boldsymbol{x})$ denotes the height of each computational cell, and $ \rho (\boldsymbol{x})$ represents the local density within the cell.

Figure 6.

Comparison of the evolution of the height of the dimensionless droplet gravity centre at ${\textit{We}}=16$ , ${\textit{Ja}}=0.61$ .

The numerical results obtained from the present model, along with the experimental data reported by Wachters et al. (Reference Wachters, Smulders, Vermeulen and Kleiweg1966), are depicted in figure 6. As illustrated in the figure, the simulation results exhibit strong agreement with both the experimental data and the numerical results in Wang et al. (Reference Wang, Fei, Lei, Wang and Luo2022a ), further validating the accuracy and reliability of the proposed model in capturing the Leidenfrost droplet dynamics.

4.3. Droplet evaporation on a heated surface

Figure 7.

Comparison of the experiment snapshots and simulation results of droplet evaporation on a heated surface: three-dimensional simulation.

Figure 8.

Profiles of the normalised radius $R^*$ with the dimensionless time $T^*$ .

The evaporation of droplets on heated surfaces represents a fundamental heat transfer process that is prevalent across numerous engineering applications. In this subsection, we aim to investigate the phenomenon of droplet evaporation on a heated surface, and validate the numerical results through comparison with the experimental data reported in Dash & Garimella (Reference Dash and Garimella2014). Within the numerical framework, the Jakob number is employed to characterise the ratio of sensible heat to latent heat during the phase transition between liquid and vapour.

In the present numerical simulation, a three-dimensional computational domain with grid resolution $ 200 \times 200 \times 200$ is employed. Initially, a droplet with diameter 70 lattice units is positioned at the centre of the heated wall, with its temperature specified as $ T_s = 0.86 T_c$ . Periodic boundary conditions are imposed in the $ x$ - and $ y$ -directions, whereas a constant-pressure boundary condition is applied at the top boundary. Furthermore, the temperatures of the bottom and top boundaries are maintained at $ T_w$ and $ T_s$ , respectively, throughout the simulations. Under the condition $ c_{\nu l} = 5.0$ , the wall temperature is set to $ T_w = 0.8931 T_c$ .

Figure 7 shows the evolution of the droplet profile at different times, as predicted by the proposed model, alongside corresponding experimental results from Dash & Garimella (Reference Dash and Garimella2014) for validation. The droplet gradually evaporates due to the vapour pressure difference between the liquid–vapour interface and the outlet. Furthermore, comparison with the experimental data confirms that the present model accurately captures the morphological evolution of the evaporating droplet.

Figure 8 shows the temporal evolution of the droplet radius from the present three-dimensional numerical model, with experimental data and numerical results from Huang et al. (Reference Huang, Wang and Yang2024) also shown for comparison. In this figure, the droplet radius and time are normalised using a reference case with droplet diameter 20 lattice units. The results from our three-dimensional simulations are in excellent agreement with both the experimental measurements and the prior numerical data.

4.4. Bubble nucleation and departure

Figure 9.

Distributions of the density contours in nucleate boiling and the comparison with experimental.

Figure 10.

Variation of bubble departure diameter with gravity acceleration.

To further evaluate the performance of the proposed model and numerical algorithm, we conduct simulations for the liquid–vapour phase change process, specifically examining bubble nucleation and departure in nucleate boiling. The computational domain is discretised with mesh resolution $100 \times 100 \times 150$ . Initially, the cavity is entirely filled with saturated liquid, and the system’s temperature is maintained at $T_s = 0.86 T_c$ . A localised heating region, represented by a square area of 4 lattice units in length at the bottom centre, is maintained at fixed temperature $T_w = 0.96 T_c$ , while the remaining walls are held at the saturated temperature $T_s = 0.86 T_c$ .

The equilibrium contact angle is prescribed as $100^{\circ }$ . A no-slip boundary condition is enforced at the solid walls, whereas periodic boundary conditions are applied along the $x$ - and $y$ -directions. The fluid properties, including the specific heat at constant volume, remain consistent with those used in the previous subsection, and the thermal conductivity is set to $2/3$ . Furthermore, a buoyant force, defined as $\boldsymbol{F}_b = (\rho - \rho _{ave})\boldsymbol{g}$ , is incorporated into the velocity field, where $\rho _{ave}$ denotes the average density over the entire system, and $\boldsymbol{g}$ represents the gravitational acceleration.

Figure 9 depicts the evolution of bubble morphology as predicted by the present thermal lattice Boltzmann model, alongside experimental observations reported in Mukherjee & Dhir (Reference Mukherjee and Dhir2004). The numerical results effectively capture the key stages of bubble nucleation, growth and departure, demonstrating strong agreement with experimental findings in terms of bubble dynamics and growth patterns.

To further quantify the numerical predictions, the power-law relationship between the bubble departure diameter and gravitational acceleration is examined. Figure 10 presents the variation of the bubble departure diameter $D_d$ as a function of gravity acceleration, where the solid line represents the theoretical relationship $D_d=0.166g^{-0.5}$ . The comparison reveals that the departure diameters predicted by different numerical models are in close agreement, and the observed dependency of $D_d$ on $g^{-0.5}$ aligns well with the theoretical results proposed in Dhir (Reference Dhir2006). These findings further validate the effectiveness of the present model and the lattice Boltzmann scheme in simulating liquid–vapour phase change phenomena.

Finally, it is worth noting that the phase transition in our deterministic model is induced by the spinodal instability of the free energy. In fact, the nucleation threshold for superheating given by the deterministic models is larger than observed in experiments. The genuine nucleation of a vapour bubble in a liquid is induced by thermal stochastic fluctuations (Lohse & Prosperetti Reference Lohse and Prosperetti2016; Menzl et al. Reference Menzl, Gonzalez, Geiger, Caupin, Abascal, Valeriani and Dellago2016). Recently, nucleation in liquid–vapour phase transitions has been studied using a fluctuation hydrodynamic theory (Gallo et al. Reference Gallo, Magaletti, Cocco and Casciola2020, Reference Gallo, Magaletti and Casciola2021, Reference Gallo, Magaletti, Georgoulas, Marengo, De Coninck and Casciola2023). In these studies, the stochastic model demonstrates excellent agreement with experimental data regarding the onset temperature, and provides a promising approach for accurately predicting nucleation phenomena.