3.1. Kinetic model

In this section, we introduce a kinetic model tailored to recover the target hydrodynamic equations of a compressible non-ideal fluid, (2.6), (2.7) and (2.8), in the hydrodynamic limit. To this end, we follow the so-called double distribution function approach and consider two velocity distribution functions, $f(\boldsymbol{v}, \boldsymbol{x}, t)$ and $g(\boldsymbol{v},\boldsymbol{x}, t)$ , where $\boldsymbol{v}$ is the velocity of a particle. The idea of a double distribution function kinetic model was first proposed by Rykov (Reference Rykov1975) for polyatomic molecules, where the second distribution function represents the rotational–vibrational contribution to the internal energy. This approach was later adopted in the LBM (He, Chen & Doolen Reference He, Chen and Doolen1998; Guo et al. Reference Guo, Zheng, Shi and Zhao2007; Li et al. Reference Li, He, Wang and Tao2007; Karlin, Sichau & Chikatamarla Reference Karlin, Sichau and Chikatamarla2013), where the second distribution function has been used to represent different forms of energy. Here, we propose a model where the bulk energy $\rho E$ (2.1) is represented by the second distribution function, although other choices are possible. This specific choice leads to kinetic equations without complicated source terms and can be efficiently tackled by classical first-neighbour discrete-velocity lattices. We refer interested readers to Hosseini, Bhadauria & Karlin (Reference Hosseini, Bhadauria and Karlin2024) and Strässle et al. (Reference Strässle, Hosseini and Karlin2025) for an in-depth discussion.

Thus, the $f$ -distribution function defines the fluid density and momentum, while the $g$ -distribution defines the bulk energy (2.1)

(3.1) \begin{align} \int \{m,m\boldsymbol{v}\} f {\rm d}\boldsymbol{v}&= \{\rho , \rho \boldsymbol{u}\}, \end{align}

(3.2) \begin{align} \int mg {\rm d}\boldsymbol{v}&= \rho E, \end{align}

where $m$ is the mass of the particle. The kinetic model is defined by the coupled kinetic equations

(3.3) \begin{align} \partial _t f + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla } \kern-1pt f &= \mathcal{F}, \end{align}

(3.4) \begin{align} \partial _t g + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla } \kern-1pt g &= \mathcal{G}, \end{align}

where the collision terms on the right-hand side are sought in the following forms:

(3.5) \begin{align} \mathcal{F} &= \frac {1}{\tau }\left (f^{\textit{eq}} - f\right ) + \frac {1}{\lambda }\left (f_\lambda ^\star -f^{\textit{eq}}\right )\!, \end{align}

(3.6) \begin{align} \mathcal{G} &= \frac {1}{\tau }\left (g^{\textit{eq}} - g\right ) + \frac {1}{\lambda }\left (g_\lambda ^\star -g^{\textit{eq}}\right )\!. \end{align}

Here, the pair $\{f^{\textit{eq}}, g^{\textit{eq}}\}$ represents a local-equilibrium attractor, while $\{f_\lambda ^\star , g_\lambda ^\star \}$ represents an intermediate quasi-equilibrium attractor. Furthermore, $\tau$ and $\lambda$ are the corresponding relaxation times, and we proceed to define the equilibrium and the quasi-equilibrium distribution functions.

To describe the equilibrium distribution functions, we introduce a reference temperature parameter $\theta$ , with dimensions of $RT$ . This parameter fixes the width of the Maxwellian equilibria independent of the thermodynamic internal energy. The equilibrium associated with mass and momentum is then written as

(3.7) \begin{equation} f^{\textit{eq}}(\rho ,\boldsymbol{u},\theta ) = \frac {\rho }{m{(2\pi \theta )}^{D/2}} \exp \left [-\frac {|\boldsymbol{v}-\boldsymbol{u}|^2}{2\theta }\right ], \end{equation}

where $\rho$ and the $D$ components of $\boldsymbol{u}$ are fixed by the local conservation constraints (3.1), while $\theta$ is specified independently. Thus, $f^{\textit{eq}}$ forms a $D+2$ -parameter family.

The equilibrium distribution associated with the energy population is defined as

(3.8) \begin{equation} g^{\textit{eq}}(\rho ,\boldsymbol{u},T,\theta ) = \left ( \frac {\boldsymbol{v}^2}{2} +e(\rho ,T) -\frac {D\theta }{2} \right ) f^{\textit{eq}}(\rho ,\boldsymbol{u},\theta ), \end{equation}

so that the thermodynamic temperature $T$ enters through the specific internal energy $e(\rho ,T)$ , whereas $\theta$ enters through the Maxwellian kernel. The equilibrium distributions (3.7) and (3.8) satisfy the conservation laws of mass, momentum and bulk energy

(3.9) \begin{align} \int \{m,m\boldsymbol{v}\} f^{\textit{eq}}\, {\rm d}\boldsymbol{v}&= \{\rho , \rho \boldsymbol{u}\}, \end{align}

(3.10) \begin{align} \int mg^{\textit{eq}}\, {\rm d}\boldsymbol{v}&= \rho E. \end{align}

Other relevant higher-order moments of the equilibrium distributions are presented in table 5 in Appendix A.

The quasi-equilibrium state $\{f_\lambda ^\star , g_\lambda ^\star \}$ is constructed using shifted values of the flow velocity, reference temperature and thermodynamic temperature

(3.11) \begin{align} \boldsymbol{u}^\star _\lambda &= \boldsymbol{u}+\lambda \frac {\boldsymbol{F}}{\rho }, \end{align}

(3.12) \begin{align} \theta ^\star _\lambda &= \theta + \lambda \alpha \theta \left (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}\right )\!, \end{align}

(3.13) \begin{align} T_\lambda ^\star &= T - \lambda ^2\frac {\boldsymbol{F}\boldsymbol{\cdot }\boldsymbol{F}}{2\rho ^{2} c_v}, \end{align}

where $\boldsymbol{F}$ is the Korteweg force (2.11), and $\alpha$ is a non-dimensional parameter to be specified below. The quasi-equilibrium distribution $f_\lambda ^\star$ therefore is a shifted Maxwellian

(3.14) \begin{equation} f_\lambda ^\star (\rho ,\boldsymbol{u}_\lambda ^\star ,\theta _\lambda ^\star ) = f^{\textit{eq}}(\rho ,\boldsymbol{u}_\lambda ^\star ,\theta _\lambda ^\star ). \end{equation}

The corresponding quasi-equilibrium for the energy population is defined as

(3.15) \begin{equation} g^\star _\lambda (\rho ,\boldsymbol{u}_\lambda ^\star , T_\lambda ^\star ,\theta _\lambda ^\star ) = \left ( \frac {\boldsymbol{v}^2}{2} + e(\rho , T_\lambda ^\star ) - \frac {D\theta _\lambda ^\star }{2} + \frac {\boldsymbol{q}^{{c}}_\lambda \boldsymbol{\cdot }(\boldsymbol{v}-\boldsymbol{u})}{\rho \theta } \right ) f^{\textit{eq}}\left(\rho ,\boldsymbol{u}_\lambda ^\star ,\theta _\lambda ^\star \right)\!. \end{equation}

The vector $\boldsymbol{q}^{{c}}_\lambda$ appearing in the last term is defined by

(3.16) \begin{equation} \boldsymbol{q}^{{c}}_\lambda = \lambda P\left (\boldsymbol{\nabla }h - \frac {k}{\mu }\boldsymbol{\nabla }T\right )\!, \end{equation}

where $h$ is the specific enthalpy:

(3.17) \begin{equation} h = e + \frac {P}{\rho }. \end{equation}

Moments of the quasi-equilibrium distributions $f_\lambda ^\star$ (3.14) and $g^\star _\lambda$ (3.15) are presented in tables 5 and 6 in Appendix A.

The construction above leaves two scalar parameters to be specified: the reference temperature $\theta$ and the non-dimensional coefficient $\alpha$ in (3.12). To specify these free parameters, we perform a Chapman–Enskog analysis of the hydrodynamic limit of the kinetic equations. Details are given in Appendix A. The analysis shows that, to recover the target hydrodynamic set (2.6), (2.7) and (2.8), the reference temperature must be set as the thermodynamic flow work

(3.18) \begin{equation} \theta =\frac {P}{\rho }. \end{equation}

Recovery of the Navier–Stokes viscous stress tensor (2.9) requires the relaxation time $\tau$ be related to the shear viscosity $\mu$ as

(3.19) \begin{equation} \tau =\frac {\mu }{P}, \end{equation}

and the parameter $\alpha$ to have the form

(3.20) \begin{equation} \alpha =\left (\frac {D+2}{D}-\frac {\rho c_s^2}{P}-\frac {\eta }{\mu }\right )\!. \end{equation}

Note that the relaxation time $\lambda$ does not affect the hydrodynamic limit, thus remaining a free parameter which can be specified conveniently in the subsequent discretisation. We close this section with several explanatory comments.

1. (i) In the ideal-gas limit, $P\to \rho RT$ , the reference temperature (3.18) becomes proportional to the thermodynamic temperature, $\theta \to RT$ . For a generic non-ideal compressible fluid, the formulation of the reference temperature (3.18), as a function of $P$ and $\rho$ , rather than thermodynamic temperature $T$ , was first proposed by Reyhanian, Dorschner & Karlin (Reference Reyhanian, Dorschner and Karlin2020) in the context of the particles-on-demand method (Dorschner et al. Reference Dorschner, Bösch and Karlin2018).

2. (ii) A special case of kinetic models (3.3) and (3.4) was recently derived from the Boltzmann–Enskog–Vlasov kinetic theory by Karlin & Hosseini (Reference Karlin and Hosseini2026) using a projection operator technique. Here, we introduced a more general kinetic model, combining the aforementioned approach of Reyhanian et al. (Reference Reyhanian, Dorschner and Karlin2020) with the quasi-equilibrium representation of the non-local Korteweg force.

3. (iii) In the absence of the term (3.16), the non-equilibrium heat flux recovered in the hydrodynamic limit is proportional to the enthalpy gradient, $\boldsymbol{q}\propto \boldsymbol{\nabla }h$ , as first observed by Reyhanian et al. (Reference Reyhanian, Dorschner and Karlin2020). In the ideal-gas limit this reduces to Fourier’s law (2.12), since $\boldsymbol{\nabla }h\propto \boldsymbol{\nabla }T$ . For a general non-ideal fluid, however, $\boldsymbol{\nabla }h$ also contains density-gradient contributions, leading to a spurious non-Fourier component of the heat flux. The correction flux (3.16) compensates for this non-Fourier contribution, thereby recovering (2.12).

4. (iv) Setting $\alpha =0$ in (3.12) recovers a fixed bulk viscosity $\eta =\mu ((({D+2})/{D})-( {\rho c_s^2}/{P}) )$ . Shifting the reference temperature with the parameter $\alpha$ (3.20) renders the bulk viscosity an independent, tuneable and positive-definite parameter.

5. (v) From the Chapman–Enskog analysis, the zeroth moment of $g_\lambda ^\star$

   (3.21) \begin{equation} \int m g_\lambda ^\star {\rm d}\boldsymbol{v} = \frac {1}{2}\rho \boldsymbol{u}^2 + \lambda \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{F} + \frac {\lambda ^2}{2\rho }\boldsymbol{F}\boldsymbol{\cdot }\boldsymbol{F} + \rho e(\rho ,T_\lambda ^\star ), \end{equation}
   must differ from the corresponding moment of $g^{\textit{eq}}$ only by the term $\lambda \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{F}$ . The additional kinetic energy introduced by the shifted velocity $\boldsymbol{u}_\lambda ^\star$ must therefore be compensated by a commensurate reduction in the internal energy. Introducing $\delta T = T_\lambda ^\star - T$ , we can expand
   (3.22) \begin{equation} e(\rho ,T_\lambda ^\star ) = e(\rho ,T) + \left .\frac {\partial e}{\partial T}\right |_{\rho } \delta T + \mathcal{O}\left (\delta T^2\right )\!. \end{equation}
   To recover the correct moment up to order $\lambda ^2$ , one therefore requires
   (3.23) \begin{equation} \left .\frac {\partial e}{\partial T}\right |_{\rho } \delta T = -\frac {\lambda ^2}{2\rho ^2}\boldsymbol{F}\boldsymbol{\cdot }\boldsymbol{F}. \end{equation}
   This yields the expression (3.13) for the shifted temperature. Since $.\partial e/\partial T|_{\rho }=c_v\gt 0$ , the shifted temperature $T_\lambda ^\star$ is systematically lower than $T$ .

6. (vi) The role of the quasi-equilibrium state can be clarified by considering the limit $\lambda \to 0$ , in which relaxation toward the shifted state becomes equivalent to introducing explicit kinetic source terms

   (3.24) \begin{align} \lim _{\lambda \to 0}\frac {1}{\lambda }(f_\lambda ^\star -f^{\textit{eq}})&= \frac {1}{P}\boldsymbol{F}\boldsymbol{\cdot }(\boldsymbol{v}-\boldsymbol{u})f^{{ eq}}+\alpha (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u})\left (\frac {{\rho (\boldsymbol{v}-\boldsymbol{u})}^2}{2P}-\frac {D}{2}\right )f^{\textit{eq}}, \end{align}
   (3.25) \begin{align} \lim _{\lambda \to 0}\frac {1}{\lambda }(g_\lambda ^\star -g^{\textit{eq}})&= \frac {1}{P}\boldsymbol{F}\boldsymbol{\cdot }(\boldsymbol{v}-\boldsymbol{u})\left (\frac {\boldsymbol{v}^2}{2} + e - \frac {DP}{2\rho }\right )f^{\textit{eq}}\nonumber\\ &\quad + \alpha (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u})\left (\left (\frac {\boldsymbol{v}^2}{2} + e - \frac {DP}{2\rho }\right )\left (\frac {{\rho (\boldsymbol{v}-\boldsymbol{u})}^2}{2P} - \frac {D}{2}\right )- \frac {DP}{2\rho }\right )f^{\textit{eq}}\nonumber\\ &\quad +{\left (\boldsymbol{\nabla }h - \frac {k}{\mu }\boldsymbol{\nabla }T\right )\boldsymbol{\cdot }(\boldsymbol{v}-\boldsymbol{u})}f^{\textit{eq}}. \end{align}

   The first term in (3.24) provides the Korteweg force contribution to the momentum balance, while its counterpart, the first term in (3.25), generates the corresponding work term in the bulk-energy balance (2.8). Furthermore, the second term in (3.24) and its counterpart in (3.25) define the bulk viscosity contributions to the momentum and bulk-energy-balance equations. Finally, the last term in (3.25) is the heat-flux correction required to recover Fourier’s law, as discussed above.

7. (vii) The double-relaxation form of the kinetic (3.3), (3.4), (3.5) and (3.6) has the advantage of admitting the standard lattice Boltzmann space–time discretisation. Integration along characteristics combined with trapezoidal quadrature of the collision terms (He et al. Reference He, Chen and Doolen1998; Ansumali et al. Reference Ansumali, Arcidiacono, Chikatamarla, Prasianakis, Gorban and Karlin2007) yields an explicit, second-order accurate scheme (see details in Appendix B)

   (3.26) \begin{align} f(\boldsymbol{x}+\boldsymbol{v}\delta t,t+\delta t)=&{f}+ 2\beta \left (f^{\textit{eq}} - {f}\right )+ {\frac {\delta t}{\lambda }}\left (1-\beta \right ) \left (f_\lambda ^\star - f^{\textit{eq}}\right )\!, \end{align}
   (3.27) \begin{align} g(\boldsymbol{x}+\boldsymbol{v}\delta t,t+\delta t)=&{g}+ 2\beta \left (g^{\textit{eq}} - {g}\right )+ {\frac {\delta t}{\lambda }}\left (1-\beta \right ) \left (g_\lambda ^\star - g^{\textit{eq}}\right )\!. \end{align}

   Here, all quantities on the right-hand side are evaluated at $(\boldsymbol{x}, t)$ , $\delta t$ is the time step and $\beta \in [0,1]$ , is the transformed relaxation parameter

   (3.28) \begin{equation} \beta =\frac {\delta t}{2\tau +\delta t}. \end{equation}
   As a consequence of the trapezoidal time integration, the transformed populations are related to the hydrodynamic fields $\rho$ , $\boldsymbol{u}$ and $E$ as follows:
   (3.29) \begin{align} \int m {f} {\rm d}\boldsymbol{v} & = \rho , \end{align}
   (3.30) \begin{align} \int m \boldsymbol{v} {f} {\rm d}\boldsymbol{v} &= \rho \boldsymbol{u}-\frac {\delta t}{2}\boldsymbol{F}, \end{align}
   (3.31) \begin{align} \int m {g} {\rm d}\boldsymbol{v} &= \rho E - \frac {\delta t}{2} \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{F}. \end{align}

This concludes the presentation and analysis of the kinetic model for compressible non-ideal fluids. The lattice Boltzmann realisation follows.

3.2. Lattice Boltzmann realisation

In the lattice Boltzmann realisation of the kinetic model introduced above, we set $\lambda =\delta t$ . We consider the standard $D3Q27$ discrete-velocity set $\boldsymbol{v}_i=c\boldsymbol{c}_i$ in $D=3$ dimensions, with $Q=27$ velocities

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

The $D3Q27$ lattice (3.32) is characterised by the lattice speed of sound

(3.33) \begin{equation} \varsigma =\frac {1}{\sqrt {3}}c. \end{equation}

In the following, we use lattice units by setting $c=1$ . The discrete-velocity equations for the populations $f_i(\boldsymbol{x},t)$ and $g_i(\boldsymbol{x},t)$ , $i=1,\ldots , Q$ , follow from (3.26) and (3.27) as

(3.34) \begin{align} &f_i(\boldsymbol{x}+\boldsymbol{v}_i\delta t, t+\delta t)= f_i+ 2\beta \left (f_i^{\textit{eq}} - f_i\right ) + \left (1-\beta \right ) \left (f_i^{\star } - f_i^{\textit{eq}}\right )\!, \end{align}

(3.35) \begin{align} &g_i(\boldsymbol{x}+\boldsymbol{v}_i\delta t, t+\delta t)= g_i+ 2\beta \left (g_i^{\textit{eq}} - g_i\right ) + \left (1-\beta \right ) \left (g_i^{\star } - g_i^{\textit{eq}}\right )\!. \end{align}

It remains, then, to define the equilibrium populations $\{f_i^{\textit{eq}}, g_i^{\textit{eq}}\}$ and the shifted-equilibrium populations $\{f_i^\star , g_i^\star \}$ . To this end, we follow the product-form formalism (Karlin & Asinari Reference Karlin and Asinari2010) and introduce functions in two variables, $\xi _{\alpha }$ and $\zeta _{\alpha \alpha }$

(3.36) \begin{equation} \varPsi _{{i\alpha }}(\xi _{\alpha },\zeta _{\alpha \alpha })=1-c_{i\alpha }^2+ \frac {1}{2}\left [\left(3c_{i\alpha }^2-2\right)\zeta _{\alpha \alpha }+c_{i\alpha }\xi _{\alpha }\right ],\ i=1,\ldots ,Q,\ \alpha =x,y,z. \end{equation}

The equilibrium populations $f_i^{\textit{eq}}$ are defined by setting the parameters in the functions (3.36) as follows:

(3.37) \begin{align} &\xi _{\alpha }^{\textit{eq}}=u_{\alpha }, \end{align}

(3.38) \begin{align} &\zeta _{\alpha \alpha }^{\textit{eq}}={\theta } +u_{\alpha }^2, \end{align}

where $\theta$ is given by the thermodynamic flow work (3.18). With the definitions (3.37) and (3.38) in the functions (3.36), the local-equilibrium populations are written in product form

(3.39) \begin{equation} f_i^{\textit{eq}}= \rho \prod _{\alpha }\varPsi _{{i\alpha }}\left (u_\alpha ,{\theta }+u_{\alpha }^2\right )\!. \end{equation}

With $\lambda =\delta t$ , the shifted flow velocity (3.11) and shifted reference temperature (3.12) are

(3.40) \begin{align} &\boldsymbol{u}^\star = \boldsymbol{u}+\delta t\left (\frac {\boldsymbol{F}}{\rho }\right )\!, \end{align}

(3.41) \begin{align} & \theta ^\star = \theta + \delta t\alpha \theta \left (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}\right )\!. \end{align}

For the shifted-equilibrium populations $f_i^\star$ , the parameters $\xi _\alpha$ and $\zeta _{\alpha \alpha }$ in the functions (3.36) are set as follows:

(3.42) \begin{align} &\xi _{\alpha }^{\star } = {u_{\alpha }^\star }, \end{align}

(3.43) \begin{align} &\zeta _{\alpha \alpha }^{\star } = {\theta ^\star } +{\left ({u_{\alpha }^\star }\right )^2} + {\delta t}\varPhi _{\alpha \alpha }. \end{align}

Compared with the continuous-velocity kinetic model of § 3.1, the discrete-velocity formulation requires the additional correction term

(3.44) \begin{equation} \varPhi _{\alpha \alpha } = -\frac {1}{\rho }\partial _{\alpha }\left (\rho u_{\alpha }^3 + {3\rho u_\alpha (\theta -\varsigma ^2)} \right )\!. \end{equation}

This is the standard correction necessary to restore Galilean invariance in the hydrodynamic limit on first-neighbour discrete-velocity lattices (Prasianakis & Karlin Reference Prasianakis and Karlin2007; Li et al. Reference Li, Luo, He, Gao and Tao2012; Hosseini, Darabiha & Thévenin Reference Hosseini, Darabiha and Thévenin2020). Combining (3.42) and (3.43) with the same product-form construction utilising (3.36), the shifted-equilibrium populations may be written as

(3.45) \begin{equation} f_i^\star =\rho \prod _{\alpha }\varPsi _{{i\alpha }}\left ({u_{\alpha }^\star }, {\theta ^\star } { + \left ({u_{\alpha }^\star }\right )^2}+{\delta t}\varPhi _{\alpha \alpha }\right )\!, \end{equation}

which completes the definition of the forcing term in (3.34).

For the $g$ -populations, we follow the generating-function representation introduced in Karlin et al. (Reference Karlin, Sichau and Chikatamarla2013) and Saadat et al. (Reference Saadat, Hosseini, Dorschner and Karlin2021b ). The generating function is the bulk energy per unit mass

(3.46) \begin{equation} E(\rho ,\boldsymbol{u},T) = e(\rho ,T)+\frac {{u}^2}{2}, \end{equation}

consistent with (2.1). The corresponding equilibrium populations are constructed by repeated application of the operators

(3.47) \begin{equation} { \mathcal{O}_\alpha (\theta )E = \theta \dfrac {\partial E}{\partial u_\alpha }+ {u}_{\alpha } E,} \end{equation}

whose dependence on the reference temperature is indicated explicitly. The discrete equilibrium $g_i^{\textit{eq}}$ is defined by setting $\xi _\alpha =\mathcal{O}_{\alpha }$ and $\zeta _{\alpha \alpha }=\mathcal{O}^2_{\alpha }$ in the functions (3.36) and interpreting the product-form as an operator acting on the generating function (3.46)

(3.48) \begin{equation} g_i^{\textit{eq}}\left (\rho ,\boldsymbol{u},T,{\theta }\right )= \rho \prod _{\alpha }\varPsi _{{i\alpha }}\left (\mathcal{O}_\alpha (\theta ),[\mathcal{O}_\alpha (\theta )]^2\right ) E(\rho ,\boldsymbol{u},T). \end{equation}

Shifted-equilibrium populations $g_i^\star$ are defined using the equilibrium product-form (3.48) evaluated at shifted values (3.40) and (3.41), and adding a correction

(3.49) \begin{equation} g_i^\star = g_i^{\textit{eq}}\left (\rho ,\boldsymbol{u}^\star , T^\star , {\theta ^\star } \right ) + \begin{cases} \dfrac {1}{2}\boldsymbol{c}_i\boldsymbol{\cdot }\boldsymbol{q}^{{c}}, & c_i^2=1,\\ 0, & \text{otherwise}, \end{cases} \end{equation}

where the non-equilibrium energy flux $\boldsymbol{q}^c$ and shifted temperature $T^\star$ are given by (3.16) and (3.13), respectively, with $\lambda =\delta t$

(3.50) \begin{align} &T^\star = T - \delta t^2\left (\frac {\boldsymbol{F}\boldsymbol{\cdot }\boldsymbol{F}}{2\rho ^{{2}} c_v}\right )\!, \end{align}

(3.51) \begin{align} &{\boldsymbol{q}^{{c}}=\delta t P\left (\boldsymbol{\nabla }h - \frac {k}{\mu }\boldsymbol{\nabla }T\right )}. \end{align}

Finally, the locally conserved fields entering the equilibrium and shifted-equilibrium populations, namely the density $\rho$ , the momentum $\rho \boldsymbol{u}$ and the bulk energy $\rho E$ , are defined from the zeroth- and first-order moments of the populations, as in (3.29), (3.30) and (3.31), with sums over discrete velocities replacing velocity-space integrals

(3.52) \begin{align} &\rho =\sum _{i=1}^Q f_i, \end{align}

(3.53) \begin{align} &\rho \boldsymbol{u}=\sum _{i=1}^Q \boldsymbol{c}_i f_i + \frac {\delta t}{2} \boldsymbol{F}, \end{align}

(3.54) \begin{align} &\rho E =\sum _{i=1}^Q g_i + \frac {\delta t}{2}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{F}. \end{align}

Once the bulk energy $E$ is computed from (3.54), the thermodynamic temperature follows from the caloric equation obtained by integrating (2.5). For a van der Waals fluid, as the equation of state (2.31) is linear in $T$ , $c_v$ exhibits no density dependence. Integrating (2.5) therefore gives the explicit relation

(3.55) \begin{equation} T = \frac {E-\boldsymbol{u}^2/2+a\rho }{c_v}. \end{equation}

This explicit inversion is specific to equations of state with this simple caloric structure. For more general equations of state, such as the Peng–Robinson equation, the temperature must instead be recovered at each node by solving the corresponding nonlinear caloric relation.

Evaluating $\boldsymbol{F}$ , $\varPhi _{\alpha \alpha }$ and $\boldsymbol{q}^{{c}}$ requires spatial derivatives of $\rho$ , $T$ , $h$ and $\boldsymbol{u}$ . In the numerical applications below, all derivatives are computed using standard central second-order accurate finite differences, except for the first-order derivative in (3.44). This term is evaluated using an upwind-biased approximation, for instance along the $x$ -axis

(3.56) \begin{equation} \varPhi _{xx} = \frac {1+\textrm {sgn}(u_x)}{2}\frac {\varLambda _x\lvert _x - \varLambda _x\lvert _{x-\delta x} }{\delta x} + \frac {1-\textrm { sgn}(u_x)}{2}\frac {\varLambda _x\lvert _{x+\delta x} - \varLambda _x\lvert _{x} }{\delta x} ,\end{equation}

where $\textrm {sgn}(u_x)$ is the sign of $u_x$ and

(3.57) \begin{equation} \varLambda _x = \rho u_x^3 + {3\rho u_x(\theta -\varsigma ^2)}. \end{equation}

This upwind-biased approximation, while maintaining the formal order of accuracy of the solver, has been shown to improve stability in higher Mach number simulations (Hosseini et al. Reference Hosseini, Darabiha and Thévenin2020; Renard et al. Reference Renard, Feng, Boussuge and Sagaut2021; Saadat et al. Reference Saadat, Dorschner and Karlin2021a ).

Figure 2.

Overall structure of the proposed algorithm for the simulation of compressible non-ideal flows.

Figure 2. Long description

The flowchart begins with the initialization of variables. It then proceeds to compute density and force. Following this, it calculates velocity, energy, temperature, pressure, a factor, and a term. The next step involves computing intermediate distributions. Subsequently, it calculates equilibrium distributions. The process continues with the computation of a factor, followed by collision and streaming steps. Finally, the time is incremented, and the loop returns to the computation of density and force.

The overall structure of the proposed algorithm is shown in figure 2. A detailed multi-scale analysis of the hydrodynamic limit is presented in Appendix C, demonstrating that the present lattice Boltzmann model recovers the hydrodynamic set (2.6), (2.7) and (2.8). Numerical applications and validation of the model are presented in the next section.

4.1. Consistency: dispersion and dissipation of hydrodynamic modes

As a first step, we probe the dispersion and dissipation properties of the hydrodynamic shear, normal and entropic modes in the limit of a resolved flow. These benchmarks will consider: (a) speed of sound, (b) shear wave dissipation, (c) shear stress, viscous heating and entropic-mode dissipation and (d) the normal-mode dissipation rate. In all cases, and without loss of generality, we consider a van der Waals fluid fitted to the critical properties of nitrogen $\textrm {N}_{2}$ , which are listed in table 1. Below, the subscript $r$ refers to reduced variables normalised by their value at the critical point.

Table 1.

Critical properties of nitrogen $\textrm {N}_{2}$ . The critical density used here comes from fitting critical temperature and pressure from Jacobsen, Stewart & Jahangiri (Reference Jacobsen, Stewart and Jahangiri1986) to the van der Waals equation of state.

Table 1. Long description

The table presents critical properties of nitrogen, detailing the specific gas constant, critical pressure, critical density, and critical temperature. It contains one row and five columns. The columns are labeled as Substance, R over c v, P c in pascals, rho c in kilograms per cubic meter, and T c in kelvin. The row lists nitrogen with values 0.4, 3.4 times 10 to the power of 6 pascals, 241.96 kilograms per cubic meter, and 126.2 kelvin respectively.

The set-up for investigating the speed of sound consists of a one-dimensional domain of size $L_x = 0.1\,\rm m$ discretised with the spacing $\delta x=10 \,\unicode{x03BC}{\rm m}$ . The initial conditions are

(4.1) \begin{align} P(x) &= \begin{cases} P_0, & x \leqslant L_x/2,\\ P_0 + \delta P, & x \gt L_x/2, \end{cases}, \end{align}

(4.2) \begin{align} u_x(x) &= u_{x0}, \end{align}

(4.3) \begin{align} T(x) &= T_0. \end{align}

Once the initial conditions are set, the system is left to evolve, resulting in two oppositely moving pressure fronts propagating at a speed that becomes constant after a short initial transition time. For sufficiently weak perturbations, this corresponds to the speed of sound in the system. A series of such cases with different initial conditions were run on both liquid and vapour branches of the saturation curve for $T_0\in [88.34,126.2]\,\textrm {K}$ corresponding to $T_0/T_c\in [0.7,1]$ and compared with the analytical speed of sound (2.17) for the van der Waals equation of state (2.31)

(4.4) \begin{equation} c_s = \sqrt {\frac {R T(1+R/c_v)}{{\left (1-b\rho \right )}^2} - 2 a\rho }. \end{equation}

Note that for a saturated vapour/liquid, fixing $T_0$ sets $\rho _0$ and $P_0$ via Maxwell’s equal-area construction and the equation of state. The results are shown in figure 3 and are in excellent agreement.

Figure 3.

Speed of sound for nitrogen $\textrm {N}_{2}$ on the saturated liquid and vapour branches. Line: analytical solution from (4.4); markers: simulations.

Figure 3. Long description

A line graph showing the speed of sound for nitrogen on the saturated liquid and vapour branches. The x axis is labeled as T r and ranges from 0.7 to 1.0. The y axis is labeled as c s in meters per second and ranges from 200 to 500. The graph features a single data line with red square markers representing simulation data points. The line starts at approximately 500 meters per second at a reduced temperature of 0.7 and decreases steadily to around 200 meters per second at a reduced temperature of 1.0. The analytical solution from equation 4.4 is represented by a solid black line that closely follows the simulation data points. All values are approximated.

The next test is a measurement of the effective shear viscosity. For this, we set up a pseudo-one-dimensional, periodic domain of size $L_x=0.1\,\rm m$ , discretised with $N_x\times N_y=100\times 1$ grid points. The initial conditions are set as follows:

(4.5) \begin{align} u_y(x) &= {\textit{Ma}}_y c_s(\rho _0, T_0) + \delta {\textit{Ma}}_y c_s(\rho _0, T_0) \sin {\left (\frac {2\pi x}{L_x}\right )}, \end{align}

(4.6) \begin{align} u_x(x) &= 0, \end{align}

(4.7) \begin{align} \rho (x) &= \rho _0, \end{align}

(4.8) \begin{align} T(x) &= T_0. \end{align}

The maximum amplitude of the perturbation, $u_y-{\textit{Ma}}_y c_s(\rho _0, T_0)$ , in the domain is monitored throughout the simulation, and its evolution over time is fitted with an exponential decay function

(4.9) \begin{equation} u_{y}^{{max}}(t) \propto \exp {\left (-{\frac {4\pi ^2}{L_x^2}} \frac {\mu }{\rho } t\right )}. \end{equation}

The viscosity measured from simulations is compared with that predicted from the multi-scale analysis. The results, shown in figure 4, demonstrate excellent agreement and Galilean invariance of the measured viscosity.

Figure 4.

Kinematic viscosity as measured from shear wave decay simulations at different Mach numbers. Solid black line: analytical viscosity; square markers: viscosity measured from simulations.

Figure 4. Long description

A line graph showing kinematic viscosity measured from shear wave decay simulations at different Mach numbers. The x axis represents Mach numbers ranging from 0 to 1. The y axis represents kinematic viscosity in units of meters squared per second, ranging from 2 to 6 times 10 to the power of negative 4. The solid black line represents analytical viscosity, while the square markers represent viscosity measured from simulations. The data points are evenly distributed along the x axis, indicating measurements at regular intervals of Mach numbers. All values are approximated.

Next, we probe the thermal conductivity by monitoring the dissipation rate of temperature perturbations. The configuration consists of a pseudo-one-dimensional, periodic domain of size $L_x\times L_y=0.1\,\rm m\times 0.001\,\rm m$ . Initial conditions are set as follows:

(4.10) \begin{align} \rho (x) &= \rho _0 + \delta \rho \sin {\left (\frac {2\pi x}{L_x}\right )}, \end{align}

(4.11) \begin{align} P(x) &= P_0(\rho _0, T_0), \end{align}

(4.12) \begin{align} T(x) &= T(\rho (x), P_0), \end{align}

(4.13) \begin{align} u_x(x) &= {\textit{Ma}}_x c_s(\rho _0, T_0), \end{align}

(4.14) \begin{align} u_y(x) &= 0, \end{align}

where $T(\rho (x),P_0)$ is obtained from the equation of state. We then monitor the maximum of the temperature in the domain $T^{ {max}}$ and extract the thermal conductivity by fitting the data to

(4.15) \begin{equation} \lvert T^{{max}}-T_0\rvert \propto \exp \left (-\frac {4\pi ^2}{L_x^2}\frac {k}{\rho _0 c_p(\rho _0,T_0)}t\right )\!. \end{equation}

The measured thermal conductivities are then compared with the imposed value. Results are shown in figure 5 and are found to be in excellent agreement with the imposed value.

Figure 5.

Thermal conductivity as measured from simulations at different Mach numbers. Plain black line: analytical thermal conductivity; square markers: thermal conductivity measured from simulations.

To validate both viscous heating and the dissipation rate of the entropic modes, we next consider the two-dimensional thermal Couette flow. The case consists of a pseudo-one-dimensional domain of size $L_x$ with walls at $x=0$ and $x=L_x$ and periodic boundary conditions along the $y$ -direction. The flow is subject to the following boundary conditions:

(4.16) \begin{align} \{u_x, u_y, T\}(x=0) &= \{0, 0, T_w\}, \end{align}

(4.17) \begin{align} \{u_x, u_y, T\}(x=L_x) &= \{0, U_w, T_w\}. \end{align}

The analytical steady-state solution to this configuration can readily be derived as

(4.18) \begin{align} u_y(x) &= U_w \frac {x}{L_x}, \end{align}

(4.19) \begin{align} T(x) &= T_w + \frac {\mu U_w^2}{2k}\frac {x}{L_x}\left (1-\frac {x}{L_x}\right )\!. \end{align}

To validate the model, we consider a domain of size $L_x=1\,\textrm {mm}$ discretised with 100 grid points. Simulations were performed for ${ Pr}\in \{0.6, 1.2, 4.9\}$ and ${\textit{Ma}}\in \{0.8, 1.2, 1.6\}$ . The results are displayed in figure 6 and show excellent agreement with the analytical solutions.

Figure 6.

Left panel: temperature and density distribution across the channel for the thermal Couette flow at different Prandtl numbers. Triangle, square and circular markers are analytical results for ${Pr}\in \{0.6, 1.2, 4.9\}$ respectively. Solid and dashed lines are temperature and density profiles from simulations. Here, ${\textit{Ma}}=0.8$ for all cases. Right panel: temperature and density distribution for different Mach numbers. Triangle, square and circular markers are analytical results for ${\textit{Ma}}\in \{0.8, 1.2, 1.6\}$ respectively. Solid and dashed lines are temperature and density profiles from simulations. Here, ${Pr}=1.2$ for all cases.

Figure 6. Long description

Two graphs illustrate temperature and density distribution across a channel for thermal Couette flow. The left graph shows distribution at different Prandtl numbers, with triangle, square, and circular markers representing analytical results. Solid and dashed lines indicate temperature and density profiles from simulations. The right graph displays distribution at different Mach numbers, with similar markers and lines. Both graphs use dual y-axes: the left y-axis for density and the right y-axis for temperature. The x-axis represents the normalized channel width. All values are approximated.

Finally, we examine the dissipation rate of normal modes, i.e. acoustics. To do this, we set up a pseudo-one-dimensional domain of size $L_x$ with periodic boundary conditions in both the $x$ - and $y$ -directions. Defining a uniform background state, $(\rho _0, T_0, P_0)$ , we add a small perturbation to it at $t=0$

(4.20) \begin{align} u_x(x) &= {\textit{Ma}}_{x} c_s(\rho _0, T_0), \end{align}

(4.21) \begin{align} u_y(x) &= 0, \end{align}

(4.22) \begin{align} P(x) &= P_0 + \delta P \sin \left (\frac {2\pi x}{L_x}\right )\!. \end{align}

The density and temperature fields can be computed using isentropic relations for the van der Waals fluid (Kouremenos, Antonopoulos & Kakatsios Reference Kouremenos, Antonopoulos and Kakatsios1988; Nederstigt & Pecnik Reference Nederstigt and Pecnik2023)

(4.23) \begin{align} \frac {P_0}{\rho _0^{\gamma _{P\rho }^0}} &= \frac {P}{\rho ^{\gamma _{P\rho }}}, \end{align}

(4.24) \begin{align} \frac {T_0}{\rho _0^{\gamma _{T\rho }^0-1}} &= \frac {T}{\rho ^{\gamma _{T\rho }-1}}, \end{align}

where

(4.25) \begin{align} \gamma _{P\rho } &= \frac {\rho }{P}\frac {c_p}{c_v} \left (\frac {\partial P}{\partial \rho }\right )_T, \end{align}

(4.26) \begin{align} \gamma _{T\rho } &= \frac {1}{\rho c_v}\left (\frac {\partial P}{\partial T}\right )_{\kern-2pt \rho} +1. \end{align}

Here $c_p$ is the specific heat at constant pressure

(4.27) \begin{equation} c_p = c_v + T{\left (\frac {\partial P}{\partial T}\right )}_v {\left (\frac {\partial v}{\partial T}\right )}_P, \end{equation}

which for the van der Waals equation of state leads to

(4.28) \begin{equation} c_p = c_v + \frac {R^2T}{RT-2a\rho {(1-b\rho )}^2}. \end{equation}

In the ideal-gas limit $P\to \rho R T$ , both exponents reduce to $(\gamma _{P\rho },\gamma _{T\rho })\to {c_p}/{c_v}$ . We leave the system to evolve over time and monitor the acoustic energy (Landau & Lifshitz Reference Landau and Lifshitz1987)

(4.29) \begin{equation} E_{{acoustic}} = \frac {1}{2}\int \left [ \rho _0{|\boldsymbol{u}-\boldsymbol{u}_0|}^2 + \frac {{\rho '}^2 c_s^2(\rho _0, P_0)}{\rho _0} \right ] {\rm d}x, \end{equation}

where $\rho '=\rho -\rho _0$ . It can readily be shown that the decay for a propagating plane wave is proportional to

(4.30) \begin{equation} E_{ {acoustic}} \propto \exp {\left (-{\frac {4\pi ^2}{L_x^2}} \sigma t\right )}, \end{equation}

where (Landau & Lifshitz Reference Landau and Lifshitz1987)

(4.31) \begin{equation} \sigma = \frac {2(D-1)}{D} \frac {\mu }{\rho }+\frac {\eta }{\rho } + \frac {k}{\rho c_p}\left (\frac {c_p}{c_v}-1\right )\!. \end{equation}

Simulations were conducted for different initial velocities and effective dissipation rates measured. Results are shown in figure 7. The results show very good agreement with the analytical predictions.

Figure 7.

Normal dissipation rate $\sigma$ as measured from normal wave decay simulations at different Mach numbers. Plain black line: analytical dissipation rate, square markers: dissipation rate measured from simulations.

4.2. Multi-phase regime

4.2.1. Liquid–vapour co-existence

Turning to the two-phase regime, we first probe the liquid–vapour co-existence densities as a validation of thermodynamic and mechanical consistency. Simulations are conducted in a pseudo-one-dimensional domain of size $L_x=0.4\,\textrm {mm}$ with periodic boundary conditions. The domain is filled with saturated vapour, with a column of saturated liquid in the centre. Simulations are evolved until the density field converges. Simulations were run for $T_r\in [0.3,0.99]$ . As shown in figure 8, numerical results agree closely with the reference values obtained using Maxwell’s equal-area construction.

Figure 8.

Liquid–vapour co-existence densities for nitrogen $\textrm {N}_{2}$ . Line: Maxwell’s equal-area rule; symbols: simulation.

Figure 8. Long description

A line graph showing liquid-vapour co-existence densities for nitrogen. The x-axis represents density in kilograms per cubic meter, ranging from 0.4 to 600. The y-axis represents temperature in Kelvin, ranging from 40 to 140. The line graph follows Maxwells equal-area rule, and the symbols represent simulation data points. The graph shows an increasing trend in temperature with increasing density up to a peak, followed by a decreasing trend. All values are approximated.

4.2.2. Interface consistency and convergence

To illustrate the consistency of the proposed solver, we repeat the nitrogen liquid–vapour co-existence density test described above, assessing the convergence of the interface with increasing resolution. All physical parameters are unchanged: the nitrogen properties are those listed in table 1, the reduced temperature is $T_r=0.9$ with corresponding co-existence densities $(\rho ^l,\rho ^v)=(399.8,102.71)\,\textrm { kg}\,\textrm {m}^{-3}$ , and the capillarity coefficient is set to $\kappa =10^{-10}\,\textrm {m}^7\,\textrm {kg}^{-1}\,\textrm {s}^{-2}$ . The computational domain has a length of $L_x=0.5\,\textrm {mm}$ . Simulations are conducted at the various resolutions detailed in table 2. The results obtained from simulations are compared with data from a high-resolution iterative finite-difference solver for

(4.32) \begin{equation} \partial _x P = \kappa \rho \partial _x^3 \rho , \end{equation}

and are shown in figure 9. The results demonstrate both excellent agreement with the reference solution and convergence under grid refinement. To further quantify solver convergence, we compare the simulated vapour density with the reference value obtained from Maxwell’s equal-area construction. The results, shown in figure 10, exhibit second-order convergence.

Table 2.

Grid properties for the nitrogen $\textrm {N}_{2}$ liquid–vapour interface simulations.

Figure 9.

Liquid–vapour interface for nitrogen $\textrm {N}_{2}$ at $T_r=0.9$ . Black lines are converged results from implicit finite-difference solver and red markers from LBM simulations. Top left panel: $\delta x=5\,\unicode{x03BC}{\rm m}$ and top right panel: $\delta x=1\,\unicode{x03BC} \textrm {m}$ . Bottom left panel: $\delta x=0.5\,\unicode{x03BC} \textrm {m}$ . Bottom right panel: $\delta x=0.1\,\unicode{x03BC} \textrm {m}$ .

Figure 9. Long description

The image contains four line graphs labeled (a), (b), (c), and (d), each showing the liquid-vapour interface for nitrogen at different resolutions. The x-axis represents the position in micrometers, ranging from 0 to 10. The y-axis represents the reduced density, ranging from 0.5 to 1.5. Black lines indicate converged results from an implicit finite-difference solver, while red markers represent data from Lattice Boltzmann Method (LBM) simulations. Each graph shows a decreasing trend in reduced density as the position increases. All values are approximated.

Figure 10.

Convergence of the vapour-phase density for nitrogen $\textrm {N}_{2}$ at $T_r=0.9$ . Markers are results from simulations while the dashed line indicates second-order convergence.

Figure 10. Long description

A scatter plot showing the convergence of the vapour-phase density for nitrogen at a specific condition. The plot features several red square markers representing results from simulations and a black dashed line indicating second-order convergence. The x-axis represents the spatial step size in meters, ranging from 10 to the power of -6 to 10 to the power of 0. The y-axis represents the normalized density difference, ranging from 10 to the power of 6 to 10 to the power of 0. The data points show a trend where the normalized density difference decreases as the spatial step size decreases, following the second-order convergence line. All values are approximated.

4.3. Compressible configurations

4.3.1. One-dimensional non-ideal shock tubes

Until the early 1980s, shock-tube experiments were limited to gases exhibiting classical wave behaviour. Borisov et al. (Reference Borisov, Borisov, Kutateladze and Nakoryakov1983) first reported a shock-tube experiment aimed at investigating non-classical wave phenomena in a dense gas, i.e. near the thermodynamic critical point. The nonlinear dynamics of gases is characterised, to a large extent, by the fundamental derivative of gas dynamics (Thompson Reference Thompson1971)

(4.33) \begin{equation} \varGamma = 1 + \frac {\rho }{c_s}\left (\frac {\partial c_s}{\partial \rho }\right )_s. \end{equation}

For simple waves, $\varGamma$ represents the rate of change of the convected sound speed with respect to density. When $\varGamma \gt 0$ , the flow exhibits positive nonlinearity, i.e. disturbances steepen forward to form compression shocks. In contrast, when $\varGamma \lt 0$ , negative nonlinearity occurs and disturbances steepen backward, leading to expansion shocks. In regions of negative nonlinearity, gases display distinct non-classical phenomena.

Following Argrow (Reference Argrow1996), we illustrate non-classical wave fields using three non-ideal shock-tube cases, each with flow regions wholly or partly in the regime of negative nonlinearity. In all three cases, simulations consist of a one-dimensional domain of size $L_x=1\textrm {m}$ , initially divided into left and right halves. The initial conditions set for each half are listed in table 3. In all cases, the grid size is set to $\delta x = 0.001\,\rm m$ . Other simulation-specific parameters are listed in table 4. Note that the Courant–Friedrichs–Lewy (CFL) is here defined as the convective CFL

Table 3.

List of initial conditions for shock-tube cases.

Table 3. Long description

The table presents initial conditions for three shock-tube cases, each with flow regions in the regime of negative nonlinearity. It contains three rows and four columns. The columns are labeled Case, R/cv, (Pr, )left, and (Pr, )right. Row 1: Case I, R/cv 0.0125, (Pr, )left 1.09, 0.879, (Pr, )right 0.885, 0.562. Row 2: Case II, R/cv 0.329, (Pr, )left 1.6077, 1.01, (Pr, )right 0.8957, 0.594. Row 3: Case III, R/cv 0.0125, (Pr, )left 3.00, 1.818, (Pr, )right 0.575, 0.275.

Table 4.

Numerical parameters for the shock-tube cases.

(4.34) \begin{equation} \textrm {CFL} = \frac {\|\boldsymbol{u}\|}{\delta x/\delta t}, \end{equation}

Figure 11.

Reduced density and pressure fields for shock tube I at time $t=0.45 L_x\sqrt {\rho _c/P_c}$ . Solid lines are reference data from Guardone & Vigevano (Reference Guardone and Vigevano2002) and markers denote simulation results.

Figure 11. Long description

The image contains two line graphs. The left graph plots reduced density ( ) against position (x) in meters. The right graph plots reduced pressure ( ) against position (x) in meters. Both graphs feature solid lines representing reference data from Guardone & Vigevano (2002) and markers denoting simulation results. The x-axis ranges from 0 to 1 meters in both graphs. The y-axis of the left graph ranges from 0.5 to 0.9, while the y-axis of the right graph ranges from 0.85 to 1.10. The graphs show distinct step-like changes in both density and pressure fields.

where $\|\boldsymbol{u}\|=\sqrt {\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}}$ . In all cases, the initial discontinuity is located at half the length of the domain. The results for case I are shown and compared with the reference data in figure 11. This configuration is a typical example of non-classical gas dynamics: first, a compression front is observed moving from the high-pressure region to low-pressure region of the domain. In classical gas dynamics, one expects a sharp compression front, which is not observed here. A second front, moving in the opposing direction, from the low-pressure region to the high-pressure region, known as a rarefaction front, is also observed. Contrary to classical gas-dynamic expectations, this front is sharp. Such non-classical wave fronts have recently been observed experimentally by Colonna et al. (Reference Colonna, Medipati and Mercier2026).

Figure 12.

(a) Fundamental derivative $\varGamma$ distribution for shock tube I at $t=0.45 L_x\sqrt {\rho _c/P_c}$ . Solid lines are reference data from Guardone & Vigevano (Reference Guardone and Vigevano2002) and markers are from simulations. (b) Fundamental derivative $\varGamma$ iso-contours in $P_r$ – $\rho _r$ plane. Red square symbols represent the state of the shock tube shown on the left panel. The solid black line is the co-existence curve. Dashed black line represents $\Gamma=1$ .

Figure 12. Long description

The image consists of two line graphs. The left graph (a) depicts the fundamental derivative distribution for shock tube I at a specific condition. The x-axis represents the spatial coordinate in meters, ranging from 0 to 1, and the y-axis represents the fundamental derivative, ranging from 1.5 to 0. The solid lines represent reference data from Guardone & Vigevano (2002), and the markers represent data from simulations. The right graph (b) shows fundamental derivative iso-contours in a plane. The x-axis represents the density ratio, ranging from 0.4 to 1.0, and the y-axis represents the pressure ratio, ranging from 0.6 to 1.2. Red square symbols represent the state of the shock tube shown in the left graph in this space. The solid black line is the co-existence curve. The color bar on the right indicates the fundamental derivative values, ranging from 2 to 6. All values are approximated.

These observations can be explained by looking into the fundamental derivative of gas dynamics (4.33), which, in the case of a van der Waals fluid, can be computed explicitly as

(4.35) \begin{equation} \varGamma (P,\rho ) = \frac {\left (R/c_v+1\right )\left (R/c_v+2\right )\frac {P+a\rho ^2}{{(1/\rho -b)}^2} - 6a\rho ^4}{2\left (R/c_v+1\right )\frac {P+a\rho ^2}{{(1/\rho -b)}}-4a\rho ^4}. \end{equation}

It is observed in figure 12 that all points in the domain at $t=0.45 L_x\sqrt {\rho _c/P_c}$ fall within the $\varGamma \lt 0$ region. Physically, $\varGamma$ measures how the sound speed changes with compression or expansion along an isentrope. It directly controls the nonlinearity of acoustic waves. Consider first the compression front; along the front, both pressure and density decrease. When the $\varGamma \lt 0$ , the characteristic speeds decrease with increasing pressure (or density). Consequently, the local propagation speed is smaller in the higher-pressure region behind the front than in the lower-pressure region ahead of it. As a result, characteristics diverge across the compression region, causing the initially sharp front to spread. The compression wave therefore evolves into a compression fan rather than steepening into a discontinuity. Next consider the rarefaction wave; along this direction, pressure and density also decrease. For $\varGamma \lt 0$ , the characteristic speeds increase as pressure (or density) decreases. This leads to a convergence of characteristics, causing the rarefaction wave to steepen and the front to sharpen. In figure 11 one can clearly observe a rarefaction shock moving from right to left, i.e. low to high pressure. Additionally, one also observes a compression fan propagating into the low-pressure region.

In the second configuration, the fronts appear to follow the classical gas-dynamic behaviour. As shown in figure 13, in agreement with reference data, both pressure and density fields show a compression front moving towards the low-pressure side and a rarefaction wave moving in the opposite direction. To further confirm the classical characteristics of this shock-tube case, we plot the distribution of $\varGamma$ at $t=0.2 L_x\sqrt {\rho _c/P_c}$ in figure 14.

Figure 13.

Reduced density and pressure fields for shock tube II at time $t=0.2 L_x\sqrt {\rho _c/P_c}$ . Solid lines are reference data from Guardone & Vigevano (Reference Guardone and Vigevano2002) and markers are from simulations.

Figure 13. Long description

Two line graphs showing reduced density and pressure fields for shock tube II at time. The left graph plots reduced density against position in meters, with the y-axis ranging from 0.6 to 1.0. The right graph plots reduced pressure against position in meters, with the y-axis ranging from 0.8 to 1.6. Solid lines represent reference data from Guardone & Vigevano (2002), and markers represent simulation data. Both graphs show a step-like decrease in values as position increases. All values are approximated.

Figure 14.

(a) Fundamental derivative $\varGamma$ distribution for shock tube II at $t=0.2 L_x\sqrt {\rho _c/P_c}$ . Solid lines are reference data from Guardone & Vigevano (Reference Guardone and Vigevano2002) and markers are from simulations. (b) Fundamental derivative $\varGamma$ iso-contours in $P_r$ – $\rho _r$ plane. red square symbols represent the state of the shock tube shown on the left in the $P_r$ – $\rho _r$ space. The solid black line is the co-existence curve.

Figure 14. Long description

The image contains two line graphs. The first graph (a) shows the fundamental derivative distribution for shock tube II. The x-axis represents the position in meters, ranging from 0 to 1. The y-axis represents the fundamental derivative, ranging from 1.4 to 2.2. The graph features two lines: a red line with markers and a black line without markers. The red line with markers represents simulation data, while the solid black line represents reference data from Guardone & Vigevano (2002). The second graph (b) shows fundamental derivative iso-contours in the plane. The x-axis represents the reduced density, ranging from 0 to 2. The y-axis represents the reduced pressure, ranging from 0 to 2. The graph features iso-contours in various colors, with red square symbols representing the state of the shock tube shown on the left in the space. The solid black line represents the co-existence curve. All values are approximated.

Figure 15.

Reduced density and pressure fields for shock tube III at time $t=0.15 L_x\sqrt {\rho _c/P_c}$ . Solid lines are reference data from Guardone & Vigevano (Reference Guardone and Vigevano2002) and markers are from simulations.

Figure 15. Long description

The image contains two line graphs. The left graph shows the reduced density field for shock tube III at a specific time. The x-axis is labeled 'x (meters)' and ranges from 0 to 1. The y-axis is labeled ' ' and ranges from 0 to 2. The solid line represents reference data from Guardone & Vigevano (2002), and the markers represent simulation data. The right graph shows the reduced pressure field for shock tube III at the same time. The x-axis is labeled 'x (meters)' and ranges from 0 to 1. The y-axis is labeled 'P ' and ranges from 0 to 3. The solid line represents reference data from Guardone & Vigevano (2002), and the markers represent simulation data. All values are approximated.

In agreement with the reference data, $\varGamma$ remains positive throughout the domain. As in the previous shock-tube configuration, the thermodynamic states sampled by the flow lie in the supercritical region and close to the saturation line. However, because of the considerably lower specific heat capacity compared with case I, no negative- $\varGamma$ region is encountered in the present case. This indicates that non-classical BZT effects are strongly influenced by the specific heat capacity and are therefore more likely to arise in fluids composed of larger, more complex molecules with higher specific heat capacities.

Finally, in the third shock-tube experiment both initial states remain in the classical region. Figure 15 shows the density and pressure profiles for this shock tube at $t=0.15 L_x\sqrt {\rho _c/P_c}$ . As in the previous configurations, a rarefaction front propagates leftward towards the dense region. Its structure, however, is non-classical: the front initially resembles a rarefaction fan, but sharpens near $x=0.5\,\textrm {m}$ . This behaviour can be explained using the fundamental derivative shown in figure 16. For $x\lesssim 0.5\,\textrm {m}$ , $\varGamma \gt 1$ and the expansion behaves classically. Near $x=0.5\,\textrm {m}$ , $\varGamma$ becomes negative; consequently the rarefaction wave steepens into a shock. A second sign change of $\varGamma$ occurs on the compression side, where the flow crosses back from the non-classical regime, associated with a compression fan, into the classical regime, where the compression front steepens into a shock.

Figure 16.

(a) Fundamental derivative $\varGamma$ distribution for shock tube III at $t=0.15 L_x\sqrt {\rho _c/P_c}$ . Solid lines are reference data from Guardone & Vigevano (Reference Guardone and Vigevano2002) and markers are from simulations. (b) Fundamental derivative $\varGamma$ iso-contours in $P_r$ – $\rho _r$ plane. Red square symbols represent the state of the shock tube shown in the left in the $P_r$ – $\rho _r$ space. The solid black line is the co-existence curve.

Figure 16. Long description

The image contains two line graphs. The first graph (a) shows the fundamental derivative distribution for shock tube III. The x-axis represents the distance in meters, ranging from 0 to 1. The y-axis represents the fundamental derivative, ranging from -2 to 4. The solid lines are reference data from Guardone & Vigevano (2002), and the markers represent data from simulations. The second graph (b) displays fundamental derivative iso-contours in a plane. The x-axis represents the density ratio, ranging from 0 to 2. The y-axis represents the pressure ratio, ranging from 0 to 3. Red square symbols represent the state of the shock tube shown in the left graph in the space. The solid black line is the co-existence curve. The color bar on the right indicates the fundamental derivative values, ranging from -2 to 8. All values are approximated.

All three configurations show excellent agreement with reference data and demonstrate that the model properly captures the behaviour of a non-classical compressible gas.

4.4. Shock–liquid-column interaction

In our final application, we showcase the suitability of the model for compressible regimes, by considering the case of a circular liquid column interacting with a planar shock wave. The case consists of a two-dimensional domain of size $L_x\times L_y$ , here resolved with $800\times 800$ grid points, divided into subdomains via a shock positioned at $x_s$ . On the right-hand side of the shock front, the pre-shock state $(\rho _1,T_1,u_{x,1})$ is set to that of a saturated vapour at $T_r=0.9$ , following Reyhanian et al. (Reference Reyhanian, Dorschner and Karlin2020). The post-shock state to the left of the shock front, $(\rho _2,T_2,u_{x,2})$ , is derived using the Rankine–Hugoniot conditions. Furthermore, a saturated liquid column of radius $R$ at $T_r=0.9$ , resolved with $65$ grid points, is placed at $(x_c,y_c)$ in the pre-shock domain. Times are non-dimensionalised by the characteristic scale $t_0 = ( {2R_0 c_s^v}/{u_s c_s^l}){\sqrt { {\rho ^l}/{\rho ^v}}}$ where $u_s$ is the shock speed, $\rho ^{l,v}$ are the liquid–vapour densities and $c_s^{l,v}$ are the saturated liquid–vapour speeds of sound. The shock speed is defined via the shock Mach number ${\textit{Ma}}_{\textrm {s}}$ as $u_s = {\textit{Ma}}_s c_s^{v}$ .

To further stabilise simulations, especially near sharp fronts, we use a nonlinear numerical-viscosity scheme as devised in Cook & Cabot (Reference Cook and Cabot2004) and Fiorina & Lele (Reference Fiorina and Lele2005). This amounts to adding a numerical contribution to the transport coefficients as follows:

(4.36) \begin{equation} \{\mu ,\eta ,k\}^{ \textit{eff}} = \{\mu ,\eta ,k\} + \{\mu ,\eta ,k\}^{ {num}}. \end{equation}

The numerical contributions are defined as

(4.37) \begin{equation} \{\mu ,\eta \}^{\textit{num}} = \rho C_{\{\mu ,\eta \}} \delta x^{r+1}\overline {\lvert \boldsymbol{\nabla }^{r-1}S\lvert }, \end{equation}

where

(4.38) \begin{equation} S=\frac {1}{2}\sqrt {\left (\boldsymbol{\nabla }\boldsymbol{u}+\boldsymbol{\nabla }\boldsymbol{u}^\dagger \right ):\left (\boldsymbol{\nabla }\boldsymbol{u}+\boldsymbol{\nabla }\boldsymbol{u}^\dagger \right )}, \end{equation}

Figure 17.

Schlieren images of shock–liquid-column interaction case at, from top to bottom, $t/t_0=0$ , $t/t_0=0.3$ and $t/t_0=0.7$ . Schlieren images are generated as $\phi = \exp (-a( {\|\boldsymbol{\nabla }\rho \|}/{\textrm {max}(\|\boldsymbol{\nabla }\rho \|)}))$ with $a=100$ . The visualisation follows Quirk & Karni (Reference Quirk and Karni1996) and Meng & Colonius (Reference Meng and Colonius2015). Solid red lines indicate the liquid-column interface identified from a density level; ISW: incident shock wave; LC: liquid column; TSW: transmitted shock wave; RSW: reflected shock wave; MS: Mach stem; R-TW: retransmitted wave; REW: reflected expansion wave.

Figure 17. Long description

The diagram illustrates the interaction of shock waves with a liquid column at three different stages. In the first stage, an incident shock wave (ISW) approaches a liquid column (LC). In the second stage, the incident shock wave interacts with the liquid column, creating a transmitted shock wave (TSW) and a reflected shock wave (RSW). In the third stage, the interaction becomes more complex with the formation of a Mach stem (MS), a retransmitted wave (R-TW), and a reflected expansion wave (REW). Solid red lines indicate the liquid-column interface identified from a density level.

and the overbar in (4.37) indicates a Gaussian filter. Furthermore (Fiorina & Lele Reference Fiorina and Lele2005)

(4.39) \begin{equation} k^{ {num}} = \rho c_s C_k \delta x^{r+1}\overline {\lvert \boldsymbol{\nabla }^{r-1}\lvert \boldsymbol{\nabla }s\rvert \rvert }, \end{equation}

with the entropy gradient obtained from (2.33) as

(4.40) \begin{equation} {\boldsymbol{\nabla }s=\frac {c_v}{T}\boldsymbol{\nabla }T-\frac {R}{\rho (1-b\rho )}\boldsymbol{\nabla }\rho .} \end{equation}

Here, to minimise the numerical dissipation we have set $r=5$ . The flow evolution, represented by a Schlieren image for ${\textit{Ma}}_s=1.47$ is shown in figure 17. The wave structures that arise during the initial stages of the shock–liquid-column interaction are commonly used to validate numerical schemes. In the present study, representative wave patterns are extracted and illustrated in figure 17. Only the early-stage interaction between a planar shock wave and a cylindrical liquid column is considered here. As the incident shock wave travels from left to right across the liquid column, both a transmitted wave and a reflected shock wave are generated. The reflected shock wave propagates upstream into the surrounding vapour, while the transmitted wave moves downstream within the liquid column. Notably, the transmitted shock wave moves faster than the incident shock wave, as the speed of sound in the liquid is larger than that in the vapour phase. Upon reaching the downstream interface of the column, the transmitted wave re-emerges into the downstream vapour. Simultaneously, expansion waves reflect repeatedly within the liquid column. At the upper lateral edge of the liquid column, the incident shock wave, the Mach stem and the slip line intersect to form a triple point. These wave structures along with the two-phase interface represent characteristic features of early-stage shock–liquid-column interaction. They appear as discontinuities of varying intensity, posing significant challenges for numerical modelling. The liquid column subsequently flattens in the flow direction and expanding in the transverse direction. As further quantitative validation, figure 18 shows the evolution of the width of the column $W$ along the centreline, compared with experiments and numerical simulations as reported in Reyhanian et al. (Reference Reyhanian, Dorschner and Karlin2020) for three Mach numbers. The results are in good agreement with both the experiments and simulations showing that the deformation of the liquid column was accurately captured by the proposed scheme.

Figure 18.

Evolution of the column width $W$ along $x$ -axis over time for three different Mach numbers. Simulations: (solid black line) ${ Ma}_s=1.47$ , (red dashed line) Ma = 1.3 and (blue dotted line) ${\textit{Ma}}_s=1.18$ . Experiments (Igra & Takayama Reference Igra and Takayama2001): (black filled circular markers) Ma = 1.47, (red filled triangle markers) ${\textit{Ma}}=1.3$ and (blue filled square markers) ${\textit{Ma}}=1.18$ . Numerical results from Reyhanian et al. (Reference Reyhanian, Dorschner and Karlin2020): (black unfilled circular markers) ${\textit{Ma}}=1.47$ , (red unfilled triangle markers) ${\textit{Ma}}=1.3$ and (blue unfilled square markers) ${\textit{Ma}}=1.18$ .

Figure 18. Long description

The line graph illustrates the evolution of the column width along the x-axis over time for three different Mach numbers. The x-axis represents the normalized time (t/t0), and the y-axis represents the normalized column width (W/D0). The graph features three solid lines representing simulation data: a solid black line for Ma = 1.47, a red dashed line for Ma = 1.3, and a blue dotted line for Ma 1.18. Additionally, the graph includes experimental data from Igra & Takayama 2001, represented by filled markers: black circles for Ma 1.47, red triangles for Ma 1.3, and blue squares for Ma 1.18. Numerical results from Reyhanian et al. 2020 are shown with unfilled markers: black circles for Ma = 1.47, red triangles for Ma 1.3, and blue squares for Ma 1.18. The data points and lines indicate how the column width changes over time for each Mach number, with all values approximated.

In addition to the flow features discussed above, we examine the thermodynamic states sampled by the solution. Unlike previously documented configurations, the present case evolves near the critical point. Figure 19 shows the thermodynamic state at $t/t_0=0.3$ . Along the centreline, from left to right, the fluid first occupies the post-shock vapour state labelled PSV, which lies above both the critical pressure and the critical-temperature isotherm and is therefore supercritical. Behind the reflected shock, the pressure and temperature increase sharply, with the pressure reaching nearly twice the critical value. Across the liquid-column interface, the pressure remains approximately constant while the temperature decreases. This brings the fluid below the critical-temperature isotherm to the state labelled PSL. Inside the liquid column, the shock reduces both pressure and temperature, driving the fluid toward the saturated liquid state. Overall, the shock interaction heats the initially saturated liquid column to temperatures close to, but slightly below, the critical temperature, while the surrounding vapour is driven well above it. The resulting configuration is a supercritical vapour surrounding in contact with a marginally subcritical liquid column.

Figure 19.

Shock–liquid-column interaction for ${\textit{Ma}}=1.47$ at $t/t_0=0.3$ . (Bottom left) numerical Schlieren image. (a) reduced pressure distribution along the $x$ -axis centreline. (b) Fundamental derivative iso-contours in $P_r$ – $\rho _r$ plane. Red square symbols represent the state of the domain. The blue dash lined is the critical-temperature isotherm. Solid black line is the co-existence curve; PSV: post-shock vapour, PRSV: post-reflected shock vapour, PSL: post-shock liquid, SL: saturated liquid and SV: saturated vapour.

Figure 19. Long description

The image contains three graphs illustrating shock liquid column interaction for a specific fluid. The first graph (a) on the top left shows a line graph of reduced pressure distribution along the x-axis centerline. The second graph (b) on the bottom left displays a numerical Schlieren image with fundamental derivative iso-contours in the plane. Red square symbols represent the state of the domain, the blue dashed line indicates the critical temperature isotherm, and the solid black line denotes the co-existence curve. The third graph (c) on the right is a phase diagram with pressure on the y-axis and density on the x-axis. It includes regions labeled PSV (post-shock vapor), PRSV (post-reflected shock vapor), PSL (post-shock liquid), SL (saturated liquid), and SV (saturated vapor). The color gradient represents the value of the fundamental derivative. All values are approximated.