1 Introduction↩︎

Non-ideal compressible fluid dynamics is a novel and rapidly developing branch of fluid mechanics, mainly due to the emergence of methods and technologies operating in the near-, trans- and super-critical regimes. It is, in part, concerned with gas-dynamics of single phase fluids in non-ideal thermodynamic states, i.e. states where compressibility factor differs from unity [1]. These states are illustrated in a P-T diagram for \({\rm CO}_2\) in Figure 1.

In addition it also includes fluids of higher molecular complexity with negative fundamental gas dynamics derivative [2], also known as Bethe-Zel’dovich-Thompson (BZT) fluids [3], [4]. The rapidly growing interest in the energy industry for fluids operating in those thermodynamic states along with their radically different behavior as compared to fluids in ideal states are two important factors demonstrating the need for systematic studies of such flows. However, as noted in the literature, while considerable effort has been put in developing experimental set-ups for non-ideal fluid dynamics in recent years [5]–[8], experimental data remains scarce and is complicated to acquire [1]. This means that development of consistent and efficient numerical tools for the simulation of such flows is of the utmost importance in the current state of development of both our understanding of non-ideal compressible fluid dynamics physics and technologies such as organic Rankine cycles or supercritical \({\rm CO}_2\) turbines [1], [9], [10]. Recent studies on dense vapor effects on fluid structures in homogeneous isotropic turbulence have shown that while large-scale flow structures are mostly affected by molecular complexity, smaller-scale flow structures are impacted by local variations in sound speed [11], [12]. In addition, the dense vapor has been shown to display modified shocklet structures, with considerably reduced jumps in pressure, density, and entropy in compression shocklets, and emergence of expansion shocklets [13], [14]. Similar importance has been reported for wall bounded turbulence where substantially different – from ideal gas – dynamics at smaller scales have been reported [15] and strong normal gradients in density and viscosity –and therefore local speed of sound and Mach number– are primary elements affecting turbulence [16], [17]. Reliable direct numerical simulation studies of non-ideal compressible fluid dynamics can fill existing gaps in understanding of physics of such flows. A reliable numerical scheme pre-supposes a physical model valid and consistently covering all regime of interest here, i.e. super-, trans- and near-critical flows involving pronounced non-ideal and compressibility effects.
LBM solvers are conservative and low-dissipative schemes with excellent spectral properties –both dissipation and dispersion– as compared to conventional solvers of the same order [18]–[23], especially for normal propagation modes, i.e. acoustic modes [24]–[26]. Along with incompressible ideal fluid dynamics, two-phase flow simulation has witnessed major success and growth in popularity in LBM, starting with formulations such as the color-gradient [27], pseudo-potential [28] and free energy [29] in the early 90’s. The latter two are of special interest in the context of non-ideal fluid dynamics as in effect they solve a form of the NSK. The large majority of the literature that focuses or uses this class of models has been tailored towards becoming essentially highly efficient interface tracking models for multi-phase flows [30]–[34]. Non-ideal compressible thermodynamics have long been neglected in LBM applications. The free energy model for instance, which is the only one out of the previously-listed approaches which results from minimization of a free energy functional under constraint on total mass, can readily be shown to be related to mean-field kinetic models such as the Boltzmann-Enskog-Vlasov [32], [35] equations and recover, at the hydrodynamic limit and under specific scaling, the NSK equations [35] and mean-field van der Waals fluid thermodynamics. This means that it can not only model two-phase flow dynamics and phase separation, but also properly recover all mean-field, near-, trans- and super-critical behaviors associated to non-ideal fluids, as demonstrated for instance through studies of properties such as the Tollman length  [32], [35]–[40].
Use and extension of this class of thermodynamically consistent model to compressible non-ideal fluid dynamics is a largely under-explored area that can considerably impact research on non-ideal compressible fluid dynamics. First attempts at proposing LBM for compressible non-ideal flows were documented in [41]. Since then, thermal multi-phase models based on the pseudo-potential approach have witnessed steady growth. However, the vast majority of the models and studies in the literature have been tailored to boiling applications, see for instance [42]–[46]. To the author’s knowledge, only documented attempts at modeling compressible non-ideal flows beyond evaporation are [47] where a hybrid lattice Boltzmann finite differences scheme –finite differences for the energy equation– was proposed and [36], [37] where the authors proposed a numerical model based on the Particles-on-Demand realization of the lattice Boltzmann method [48]. We propose to address this gap in the literature with a lattice Boltzmann model for non-ideal fluids in the compressible regime. To retain the main advantages of the lattice Boltzmann method the proposed model relies on classical first neighbor lattices, taking advantage of a second distribution function for the energy balance equations [49]–[52]. In addition, the model removes the Galilean-variant errors in the diagonal components of the third-order equilibrium moments tensor and provides for independent control over the bulk viscosity; The latter is a critical point as the bulk viscosity dictated by the BGK structure, for equilibria carrying a pressure other than the ideal gas pressure, can take on negative values in the hydordynamic limit, see [35]. Together with a consistent second-order discretization and treatment of source terms, the model will be shown to correctly recover the target hydrodynamic limit.
The article starts with a brief introduction of the target hydrodynamic limit followed by presentation of the proposed lattice Boltzmann model. The proposed model is then validated with a list of configurations of increasing complexity.

2 Balance equations for compressible non-ideal fluid↩︎

We begin with a brief overview of the compressible non-ideal fluid system.

\[\label{eq:de} {\color{black}{de=c_vdT+\left[T\left(\dfrac{\partial P}{\partial T}\right)-P\right]{dv},}}\tag{1}\]

The \(D\)-dimensional mass, momentum and bulk energy balance equations are, \[\begin{align} & \partial_t \rho + \boldsymbol{\nabla}\cdot \rho\boldsymbol{u} = 0, \tag{2} \\ &\partial_t (\rho\boldsymbol{u}) + \boldsymbol{\nabla}\cdot \rho\boldsymbol{u}\otimes\boldsymbol{u} + \boldsymbol{\nabla}P + \boldsymbol{\nabla}\cdot \boldsymbol{T}_{\rm NS} + \boldsymbol{\nabla}\cdot \boldsymbol{T}_{\rm K}= 0, \tag{3}\\ &{\color{black}{\partial_t (\rho E) + \boldsymbol{\nabla}\cdot \left(\rho E\boldsymbol{u}+{P}\boldsymbol{u}\right) + \boldsymbol{u}\cdot\left(\boldsymbol{\nabla}\cdot \boldsymbol{T}_{\rm K}\right) + \boldsymbol{\nabla}\cdot(\boldsymbol{u}\cdot\boldsymbol{T}_{\rm NS}) + \boldsymbol{\nabla}\cdot\boldsymbol{q}= 0.}} \tag{4} \end{align}\] The viscous stress tensor is defined as, \[\label{eq:TNS} \boldsymbol{T}_{\rm NS} = {-}\mu\left(\boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^\dagger - \frac{2}{D}(\boldsymbol{\nabla}\cdot\boldsymbol{u})\boldsymbol{I}\right) - \eta (\boldsymbol{\nabla}\cdot\boldsymbol{u})\boldsymbol{I},\tag{5}\] Furthermore, the Korteweg surface tension tensor is defined as, \[\boldsymbol{T}_{\rm K} = \kappa \boldsymbol{\nabla}\rho\otimes\boldsymbol{\nabla}\rho - \kappa\left(\rho \boldsymbol{\nabla}^2\rho+\frac{1}{2}{\lvert \boldsymbol{\nabla}\rho\lvert}^2\right)\boldsymbol{I},\] Finally, the Fourier heat flux is defined as, \[\boldsymbol{q} = -k\boldsymbol{\nabla}T,\]

Consequently, and without loss of generality, we use the van der Waals equation of state in the numerical examples below, \[\label{eq:Pvdw} P(\rho,T) = \frac{\rho R T }{1 - b\rho} - a\rho^2.\tag{6}\] The excluded volume parameter \(b\) and the long-range molecular attraction parameter \(a\) are defined in terms of the critical state thermodynamic data, critical point density \(\rho_c\), critical temperature \(T_c\) and critical pressure \(P_c\), as follows: \(a ={27 R^2 T_c^2}/{64 P_c}\), \(b = {R T_c}/{8 P_c}\). The differentials of the specific internal energy [eq:de95rho] and of the specific entropy \(s\) for the van der Waals fluid are, respectively, \[\begin{align} \tag{7} &{\color{black}{de={c_vdT}-ad\rho,}}\\ \tag{8} &{\color{black}{ds=\frac{c_vdT}{T}-\frac{Rd\rho}{\rho(1-b\rho)}.}} \end{align}\] In the next section we shall introduce a lattice Boltzmann model that recovers the above system of balance equations 2 , 3 and 4 in the hydrodynamic limit.

3 Lattice Boltzmann model for non-ideal compressible flows↩︎

We consider the standard \(D3Q27\) discrete velocity set \(\boldsymbol{v}_i=c\boldsymbol{c}_i\) in \(D=3\) dimensions, with \(Q=27\) velocities, \[\label{eq:D2Q9c} \boldsymbol{c}_i=(c_{ix},c_{iy},c_{iz}),\;c_{i\alpha}\in\{-1,0,1\}.\tag{9}\] The \(D3Q27\) lattice 9 is characterized with the lattice speed of sound, \[\label{eq:cs} \varsigma=\frac{1}{\sqrt{3}}c.\tag{10}\] Below, we consider lattice units by setting \(c=1\). Motivated by the generic approach initialized in [35], [53], we consider the lattice Boltzmann equations for the two sets of populations \(f_i(\boldsymbol{x},t)\) and \(g_i(\boldsymbol{x},t)\), \(i=1,\dots, Q\), \[\label{eq:beta} {\color{black}{\beta = \frac{P\delta t}{{2\mu}+P\delta t}.}}\tag{11}\] The last term in [eq:fLBGK] and [eq:gLBGK] is a generalized forcing, characterized by shifted-equilibrium populations \(f_i^*\) and \(g_i^*\), respectively.

To define the discrete equilibrium and shifted equilibrium, we first introduce functions in two variables, \(\xi_{\alpha}\) and \(\zeta_{\alpha\alpha}\), \[\label{eq:phi} \Psi_{{i\alpha}}(\xi_{\alpha},\zeta_{\alpha\alpha})=1-c_{i\alpha}^2+ \frac{1}{2}\left[(3c_{i\alpha}^2-2)\zeta_{\alpha\alpha}+c_{i\alpha}\xi_{\alpha}\right].\tag{12}\] The equilibrium populations \(f_i^{\rm eq}\) are defined by setting the parameters as follows: \[\begin{align} &\xi_{\alpha}^{\rm eq}=u_{\alpha},\tag{13}\\ &\zeta_{\alpha\alpha}^{\rm eq}=\frac{P}{\rho}+u_{\alpha}^2.\tag{14} \end{align}\] With the definitions 13 and 14 in the functions 12 , the local equilibrium populations are represented with a product-form, \[\label{eq:LBMeq} f_i^{\rm eq}= \rho\prod_{\alpha}\Psi_{{i\alpha}}\left(u_\alpha,\frac{P}{\rho}+u_{\alpha}^2\right).\tag{15}\] For the shifted-equilibrium populations \(f_i^*\), the parameters \(\xi_\alpha\) and \(\zeta_{\alpha\alpha}\) in the functions 12 are set as follows:

\[\begin{align} &\xi_{\alpha}^{*} = u_{\alpha}+{\delta t}\frac{F_{\alpha}}{\rho},\tag{16} \\ &\zeta_{\alpha\alpha}^{*} = \frac{P}{\rho} +\left(u_{\alpha}+{\delta t}\frac{F_{\alpha}}{\rho}\right)^2 + {\delta t}\frac{\Phi_{\alpha\alpha}}{\rho}.\tag{17} \end{align}\] When compared to their equilibrium counterparts, the added terms in [eq:xistar] and [eq:zetastar] are, respectively, \[\begin{align} & {F}_\alpha = -\kappa \rho \partial_\alpha \boldsymbol{\nabla}\cdot\boldsymbol{\nabla} \rho,\tag{18}\\ \tag{19} & \Phi_{\alpha\alpha} = -\partial_{\alpha}\left(\rho u_{\alpha} \left(u_{\alpha}^2 + \frac{3P}{\rho}-3\varsigma^2\right)\right) + {\Phi'},\\ &\Phi'={P\left(\frac{D+2}{D}-\frac{\rho c_s^2}{P} - \frac{\eta}{\mu}\right) (\boldsymbol{\nabla}\cdot\boldsymbol{u})}, \end{align}\] where \(c_s\) is the speed of sound [eq:speed95of95sound]. We finalize the construction of the shifted-equilibrium populations by using again the products of functions 12 but with the parameters set according to [eq:xistar] and [eq:zetastar], \[\label{eq:LBMstar} f_i^*=\rho\prod_{\alpha}\Psi_{{i\alpha}}\left(u_\alpha+{\delta t}\frac{F_{\alpha}}{\rho},\frac{P}{\rho}+{\left(u_{\alpha}+{\delta t}\frac{F_{\alpha}}{\rho}\right)}^2+{\delta t}\frac{\Phi_{\alpha\alpha}}{\rho}\right),\tag{20}\] which allows us to define the forcing term in the lattice Boltzmann equation [eq:gLBGK]. Before moving on to the second distribution function a few comments are in order:

• Eq. 19 is made up of two contributions. The first term is introduced to correct errors associated with diagonal components of the viscous stress tensor stemming from the \(c_i^3=c_i\) bias of the \(D3Q27\) lattice. This correction fixes the bulk viscosity to \(\eta'=\left(\frac{D+2}{D}-\frac{\rho c_s^2}{P}\right)\mu\).

• The above-mentioned bulk viscosity is only positive-definite under \(\frac{D+2}{D}\geq\frac{\rho c_s^2}{P}\) which is not guaranteed to hold. The second term in Eq. 19 is introduced to set the hydrodynamic dissipation rate of normal modes to an independent value of \(\eta\) guaranteeing positive-definiteness.

Moving onto the \(g\)-populations, we consider a generating function (bulk energy [eq:rhoE] per unit mass), \[\label{eq:Egen} E(\rho,\boldsymbol{u},T) = e(\rho,T)+\frac{{u}^2}{2},\tag{21}\] and introduce operators \(\mathcal{O}_\alpha\) acting on the generating function \(E\) as, \[\label{eq:Oa} \mathcal{O}_\alpha E = \left(\frac{P}{\rho}\right) \dfrac{\partial E}{\partial u_\alpha}+ {u}_{\alpha} E.\tag{22}\] The energy moments can be computed by repeated application of the operators 22 on the generating function 21 . For clarity, let us compute the first two moments explicitly, \[\begin{align} \mathcal{O}_{\alpha} E &= \left(\frac{P}{\rho}+E\right) {u}_{\alpha},\\ \mathcal{O}_{\alpha}^2 E &= u_\alpha^2\left(\frac{2P}{\rho} + E\right) + \frac{P}{\rho}\left(\frac{P}{\rho} + E\right). \end{align}\] \[g_i^*(\rho,\boldsymbol{u},T) = g_i^{\rm eq}\left(\rho,\boldsymbol{u}^*,T^*\right) + \left\{\begin{align} & \frac{1}{2}\boldsymbol{c}_i\cdot\boldsymbol{q}', &\text{ if } c_i^2=1, & \\ &0, & \text{otherwise}.&\\ \end{align}\right.\]

where the non-equilibrium energy flux correction is the vector \(\boldsymbol{q}'\) with the components, \[\label{eq:correction95q} {q}_{\alpha}' = \delta t \left[P\partial_{\alpha}\left(e+\frac{P}{\rho}\right) - \left(\frac{k}{\mu}\right)T\partial_\alpha T+ {u}_{\alpha}{\Phi}'\right].\tag{23}\] The moments of orders zero and one of \(g^*_i\) are as follows: \[\begin{align} & \sum_{i=1}^Q g_i^*= \rho E\left(\rho,\boldsymbol{u}^*, T^*\right),\tag{24}\\ & \sum_{i=1}^Q \boldsymbol{c}_i g_i^*= \boldsymbol{u}^*\left(P(\rho,T^*) + \rho E(\rho,\boldsymbol{u}^*,T^*)\right) + \boldsymbol{q}'.\tag{25} \end{align}\] A few comments about \(g^*_i\) are in order:

• The equilibrium of Eq. 21 would have naturally resulted in a Navier-Stokes level heat flux that is a function of specific enthalpy, i.e. \(\propto \partial_\alpha (e+P/\rho)\). In the case of an ideal gas equation of \(\partial_\alpha (e+P/\rho)=(c_v+R)\partial_\alpha T\) since both internal energy and pressure are only functions of temperature. For a non-ideal equation of state however, since the latter are functions of temperature and density this would lead to an additional flux driven by density gradients. The first two terms in Eq. 23 restore the correct form of the Fourier heat flux to the model while the last term accounts for viscous heating brought in by the bulk viscosity in the \(f\)-population quasi-equilibrium, see Eq. 19 .

• Looking at \(g_i^*\) as an attractor towards a shifted velocity \(\boldsymbol{u}^*\) and temperature \(T^*\), we note that the shifted temperature is systematically lower than \(T\), see Eq. [eq:Tstar]. Developing the moment of order zero as, \[\sum_{i=1}^Q g_i^* = \frac{1}{2}\rho\boldsymbol{u}^2 + \delta t \boldsymbol{u}\cdot\boldsymbol{F} + \frac{\delta t^2}{2\rho}\boldsymbol{F}\cdot\boldsymbol{F} + \rho e(\rho,T^*),\] with, \[e(\rho,T^*) = e(\rho,T) +\left.\frac{\partial e(\rho, T)}{\partial T}\right|_{T} \delta T + \mathcal{O}\left(\delta T^2\right),\] where we introduced \(\delta T=T^*-T\) for readability. To correctly recover fluxes up to order \(\delta t^2\) one must have, \[\left. \frac{\partial e(\rho, T)}{\partial T}\right|_{T} \delta T = -\frac{\delta t^2}{2\rho^2}\boldsymbol{F}\cdot\boldsymbol{F},\] which explains Eq. [eq:Tstar] and why \(T^*\) is systematically lower than \(T\). Note that while higher order terms appear for non-linear-in-temperature equations of state, such as Peng–Robinson [54], for linear equations such as van der Waals there are no higher order terms, see Eq. [eq:eos95general] and [eq:cv95zero].

Finally, the locally conserved fields elsewhere in the equilibrium and shifted-equilibrium populations, the density \(\rho\), the momentum \(\rho\boldsymbol{u}\) and the bulk energy \(\rho E\), are consistently defined via the zeroth- and the first-order moments of the populations, \[\begin{align} &\rho=\sum_{i=1}^Q f_i,\tag{26}\\ &\rho\boldsymbol{u}=\sum_{i=1}^Q \boldsymbol{c}_i f_i + \frac{\delta t}{2} \boldsymbol{F},\tag{27}\\ &\rho E =\sum_{i=1}^Q g_i + \frac{\delta t}{2}\boldsymbol{u}\cdot\boldsymbol{F}.\tag{28} \end{align}\]

4 Numerical applications↩︎

4.1 Consistency: Dispersion and dissipation of hydrodynamic eigen-modes↩︎

As a first step, we shall probe the dispersion and dissipation properties of hydrodynamic modes comprising the shear, the normal and the entropic modes, in the limit of resolved flows simulation. This benchmark will involve: (a) speed of sound, (b) dissipation of shear waves, (c) shear stress viscous heating and entropic mode dissipation, and (d) normal mode dissipation rate. For all cases, for the sake of readability and without loss of generality, we consider Nitrogen with properties listed in Table 1.

The setup to investigate the speed of sound consists of a one-dimensional domain of size \(L_x = 0.1~[m]\) discretized with the spacing \(\delta x=10 [\mu m]\). The initial conditions are, \[\begin{align} P(x,t) &=& \begin{cases} P_0,\;\forall x\leq L_x/2\\ P_0 + \delta P,\;\forall x> L_x/2 \end{cases},\\ u_x(x,t) &=& u_{x0},\\ T(x,t) &=& T_0. \end{align}\] Once the initial conditions are set, the system is left to evolve, resulting in two opposite-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 for different initial conditions on both liquid and vapor branches of the saturation curve for \(T_r\in[0.7,1]\) were run and compared to the analytical speed of sound [eq:speed95of95sound] for the van der Waals equation of state 6 , \[\label{eq:vdW95speed95of95sound} c_s = \sqrt{\frac{R T(1+R/c_v)}{{\left(1-b\rho\right)}^2} - 2 a\rho}.\tag{29}\] The results are shown in Figure 2 and point to excellent agreement.

The next configuration is set to measure the effective shear viscosity. To that end, we set up a one-dimensional domain of size \(L_x=0.1~[m]\) with periodic boundary conditions in both \(x-\) and \(y-\)directions and discretized with \(N_x=100\) grid points. Initial conditions are set as, \[\begin{align} u_y(x) &=& {\rm Ma}_y c_s(\rho_0, T_0) + \delta {\rm Ma}_y c_s(\rho_0, T_0) \sin{\left(\frac{2\pi x}{L_x}\right)},\\ u_x(x) &=& 0,\\ \rho(x) &=& \rho_0,\\ T(x) &=& T_0. \end{align}\] The maximum \(u_y-{\rm Ma}_y c_s(\rho_0, T_0)\) in the domain is monitored throughout the simulation and its evolution in time is fitted with an exponential decay function, \[u_{y}^{\rm max}(t) \propto \exp{\left(-{\frac{4\pi^2}{L_x^2}} \frac{\mu}{\rho} t\right)}.\] The viscosity measured from simulations is compared to those predicted from the multi-scale analysis. The results are shown in Fig. 3 and point to excellent agreement and Galilean invariance of the measured viscosity.

As the next case, to validate both viscous heating and entropic modes dissipation rate, we 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, \[\begin{align} \{u_x, u_y, T\}(x=0) &=& \{0, 0, T_w\},\\ \{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, \[\begin{align} u_y(x) &=& U_w \frac{x}{L_x},\\ 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~[{\rm mm}]\) discretized with 100 grid-points. Simulations were run for \({\rm Pr}\in\{0.6, 1.2, 4.9\}\) and \({\rm Ma}\in\{0.8, 1.2, 1.6\}\). Results are displayed in Fig. 4 and show excellent agreement with analytical solutions.

Finally, we move on to probing dissipation rate of normal modes, i.e. acoustics. To that end we set up a 1-D domain of size \(L_x\) with periodic boundary conditions in both \(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\), \[\begin{align} u_x(x) &=& {\rm Ma}_{x} c_s(\rho_0, T_0),\\ u_y(x) &=& 0,\\ P(x) &=& P_0 + \delta P. \end{align}\] The density and temperature fields can be readily computed using isentropic relations for van der Waals fluid [55], [56], \[\begin{align} \frac{P_0}{\rho_0^{\gamma_{P\rho}^0}} &=& \frac{P}{\rho^{\gamma_{P\rho}}},\\ \frac{T_0}{\rho_0^{\gamma_{T\rho}^0-1}} &=& \frac{T}{\rho^{\gamma_{T\rho}-1}}, \end{align}\] where, \[\begin{align} \gamma_{P\rho} &=& \frac{\rho}{P}\frac{c_p}{c_v} \left(\frac{\partial P}{\partial \rho}\right)_T,\\ \gamma_{T\rho} &=& \frac{1}{\rho c_v}\left(\frac{\partial P}{\partial T}\right)_\rho+1. \end{align}\] Here \(c_p\) is, \[c_p = c_v + T{\left(\frac{\partial P}{\partial T}\right)}_v {\left(\frac{\partial v}{\partial T}\right)}_P,\] which for the van der Waals equation of state leads to, \[c_p = c_v + \frac{R^2T}{RT-2a\rho{(1-b\rho)}^2}.\] One can readily show that for \(P\rightarrow \rho R T\), i.e. the limit of an ideal gas, the following holds \((\gamma_{P\rho},\gamma_{T\rho})\rightarrow\frac{c_p}{c_v}\). We leave the system to evolve over time and monitor the acoustic energy [57], \[E_{\rm 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] dx,\] where \(\rho'=\rho-\rho_0\). It can readily be shown that the decay for a propagating plane wave is proportional to, \[E_{\rm acoustic} \propto \exp{\left(-{\frac{4\pi^2}{L_x^2}} \alpha t\right)}.\] with [57], Simulations were conducted for different initial values of velocity and effective dissipation rates measured. Results are shown in Fig. 5.

The data points to very good agreement with analytical predictions.

4.2 Multi-phase regime↩︎

4.2.1 Liquid-vapour coexistence↩︎

Here we look into the multiphase regime; as a standard validation of thermodynamic and mechanical consistency we first probe the co-existence densities. For that purpose we consider the van der Waals equation of state fitted to the \(N_2\) critical point. See table 1 for properties. Simulations are conducted in a 1-D domain of size \(L_x=0.4~[{\rm mm}]\) with periodic boundary conditions. The domain is filled with saturated vapor with a column of saturated liquid at the center. Simulations are left to evolve until convergence of the density field.

Simulations were run for \(T_r\in[0.3\,0.99]\). Results are shown in Fig. 6 and show excellent agreement with theory. Theoretical results are obtained using the Maxwell equal area construction.

4.2.2 Interface consistency and convergence↩︎

To illustrate the consistency and convergence of the proposed solver to the target hydrodynamic system we compare the density profile of Nitrogen liquid vapor interface. Nitrogen system critical properties are listed in table 1. The capillarity coefficient is set to \(\kappa=10^{-10}~[{\rm m}^7/{\rm kg} .{\rm s}^2]\). The simulation consists of a periodic 1-D domain of size \(L_x=0.5[{\rm mm}]\) with a liquid column at the center surrounded with vapor. We consider a system at \(T_r=0.9\). The simulations are initialized with the corresponding co-existence densities in each phase, i.e. \((\rho^l,\rho^v)=(399.8,102.71)~[{\rm kg}/{\rm m}^3]\).

The results obtained from simulations are compared to data from a high resolution iterative finite differences solver for, \[\partial_x P = \kappa \rho \partial_x^3 \rho.\] The results are shown in Fig. 7.

Results point to excellent agreement with reference solution and convergent behavior as resolution increases.

4.3 Compressible configurations↩︎

4.3.1 1-D non-ideal shock tubes↩︎

Until the early 1980s, shock tube experiments were limited to gases exhibiting classical wave behavior. [58] first reported a shock tube experiment aimed at investigating the non-classical wave phenomena in a dense gas, i.e. near the thermodynamic critical point. The nonlinear dynamics of gases are fundamentally governed by the fundamental derivative of gas dynamics [59]: \[\Gamma = 1 + \frac{\rho}{c_s}\left(\frac{\partial c_s}{\partial \rho}\right)_s.\] For simple waves, \(\Gamma\) represents the rate of change of the convected sound speed with respect to density. When \(\Gamma>0\), the flow exhibits positive nonlinearity, i.e. disturbances steepen forward to form compression shocks. Conversely, when \(\Gamma<0\), negative nonlinearity occurs and disturbances steepen backward, leading to expansion shocks. In regions of negative nonlinearity, gases display distinct non-classical phenomena.
In all three cases, simulations consist of a 1-D domain of size \(L_x=1[{\rm m}]\), initially divided into right and left halves. The initial conditions set for each half are listed in table 3.

In all presented cases the grid-size is set to \(\delta x = 0.001~[m]\). Other simulation-specific parameters are listed in table 4.

In all cases the shock front is initially located on the half-length of the domain. Results for case I are shown and compared to reference data in Fig. 8.

This configuration is a typical example of non-classical gas dynamics as both left and right half-domains lead to \(\Gamma<0\). In Fig. 8 one can clearly observe a rarefaction shock moving from right to left, i.e. low to high pressure. Additionally, one also observes a compression waves propagating into the low pressure region.
The second configuration does not fall in the non-classical regime. As shown in Fig. 9, 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.

Finally, the last tube is different from previous cases in the sense that it dynamically crosses the \(\Gamma=0\) line during the simulations. Both sides of the initial conditions are in the classical region. The crossing of \(\Gamma=0\) results in the formation of mixed rarefaction wave composed of a rarefaction shock connected to a rarefaction fan. Additionally, as shown in Fig. 10 a compression shock is traveling toward the low pressure area.

All three configurations show excellent agreement with reference data and show that the model properly captures non-classical compressible gas dynamic behavior.

4.4 Shock-liquid column interaction↩︎

As our final configuration, to further showcase the consistency and suitability of the model for compressible regimes we consider the case of a circular liquid column interacting with a planar shock-wave. The case consists of a 2-D domain of size \(L_x\times L_y\), here resolved with \(800\times800\) grid-points, divided into subdomains via a shock positioned at \(x_s\). On the right hand side of the shock the pre-shock state \((\rho_1,T_1,u_{x,1})\) is set to that of a saturated vapor at \(T_r=0.9\), following [36]. The post-shock state, to the left of the shock, \((\rho_2,T_2,u_{x,2})\) is derived using the Rankine-Hugoniot conditions. Furthermore, a saturated liquid column at \(T_r=0.9\) of radius \(R\), here resolved with \(65\) grid-points, is placed at \((x_c,y_c)\) in the pre-shock domain. Times are non-dimensionalized using the characteristic time \(t_0 = \frac{2R_0 \varsigma^v}{u_s \varsigma^l}\sqrt{\frac{\rho^l}{\rho^v}}\) where \(u_s\) is the shock speed, \(\rho^{v,l}\) are the vapor/liquid densities and \(\varsigma^{v,l}\) are the saturated vapor/liquid speed of sound. The shock speed is defined via the shock Mach number \({\rm Ma_s}\) as \(u_s = {\rm Ma}_s \varsigma^{v}\).
To further stabilize simulations especially near sharp fronts, we use a non-linear numerical viscosity scheme as devised in [61], [62]. This consists in adding a numerical contribution to dissipation coefficients as, \[\{\mu,\eta,k\}^{\rm eff} = \{\mu,\eta,k\} + \{\mu,\eta,k\}^{\rm num},\] where the numerical contributions are defined as, \[\label{eq:shock95cap} \{\mu,\eta\}^{\rm num} = C_{\{\mu,\eta,k\}} \delta x^{r+1}\overline{\lvert \boldsymbol{\nabla}^{r-1}S\lvert},\tag{30}\] where, \[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)},\] and the overbar in Eq. 30 indicates a Gaussian filter. Furthermore [62], \[\label{eq:shock95cap95k} k^{\rm num} = \rho \varsigma C_k \delta x^{r+1}\overline{\lvert \boldsymbol{\nabla}^{r-1}\lvert\boldsymbol{\nabla}s\lvert\lvert},\tag{31}\] \[{\color{black}{\boldsymbol{\nabla}s=\frac{c_v}{T}\boldsymbol{\nabla}T-\frac{R}{\rho(1-b\rho)}\boldsymbol{\nabla}\rho.}}\] Here to minimize the numerical dissipation we have set \(r=5\). The evolution of the field, represented via a Schlieren image for \({\rm Ma}_s=1.47\) is shown in Fig. 11.

The wave structures that arise during the initial stages of shock-droplet interaction are commonly used to validate numerical schemes. In the present study, representative wave patterns are extracted and illustrated in Figure 11. Only the early-phase interaction between a planar shock wave and a cylindrical water 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 vapor, while the transmitted wave moves downstream within the liquid column. Note that the transmitted shock wave moves faster than the incident shock wave, as speed of sound in the liquid is larger than in the vapor phase. Upon reaching the downstream wall of the column, the transmitted wave re-emerges into the downstream vapor. Simultaneously, expansion waves reflect repeatedly within the liquid column. At the upper lateral edge of the water column, the incident shock wave, Mach stem, and slip line intersect to form a triple point. These wave structures along with the two-phase interface represent characteristic features of early-stage shock-droplet interaction. They appear as discontinuities of varying intensity, posing significant challenges for numerical modeling. The liquid column then starts flattening in the flow direction and expanding in the radial direction. As further quantitative validation Figure 12 represents the evolution of the width of the column, \(W\) on the center-line along the \(x\)-axis as compared to experiments and numerical simulations as reported in [36] for three different Mach numbers.

The results are in good agreement with both the experiments and simulations showing that the deformation of the droplet was accurately captured by the proposed scheme.

5 Conclusion↩︎

Development and validation of lattice Boltzmann models for non-ideal fluids in the compressible regime is an area that has been left under-developed. Given the growing interest in technologies involving flows in compressible regime where non-ideal effects are quite relevant, development of numerical tools able to capture non-ideal effects in the compressible regime is of critical importance. In the present contribution, we presented a lattice Boltzmann model for such flows. The model uses two double distribution functions, allowing it to be consistently realized on a classical first-neighbor lattice and providing for the possibility to freely tune the specific heat capacity at constant volume. In addition, to maintain positive-definiteness the bulk viscosity is made independent from the relaxation coefficient and isentropic sound speed [35], [66]. Surface tension is introduced into the system via a consistent realization of the force contribution, in both distribution functions. All of these ingredients led to a numerical model that properly captures the target physics, as evidenced from various studied test-cases, and allow for high speed simulation within the efficient framework of the lattice Boltzmann method. In the present work a simple single relaxation time formulation was used to carry out simulations. Extension of the solver to more robust collision models better suited to carry out turbulent configurations will be the topic of future publications.

Acknowledgment↩︎

This work was supported by European Research Council (ERC) Advanced Grant No. 834763-PonD and by the Swiss National Science Foundation (SNSF) Grants 200021-228065 and 200021-236715. Computational resources at the Swiss National Super Computing Center (CSCS) were provided under Grants No. s1286, sm101 and s1327. Authors would like to thank Patrick Jenny for his support and fruitful discussions.

Declaration of interests↩︎

The authors report that they do not have a conflict of interest.

Data Availability Statement↩︎

The data that support the findings of this study are available from the corresponding author upon request.

6 Multi-scale analysis of lattice Boltzmann model for compressible non-ideal flows↩︎

The first step in the multi-scale analysis is a Taylor expansion of the lattice Boltzmann equations, \[\begin{gather} \{f_i,g_i\}(\boldsymbol{x}+\boldsymbol{c}_i\delta t, t+\delta t)= \{f_i,g_i\}(\boldsymbol{x}, t) + 2\beta\left(\{f_i^{\rm eq},g_i^{\rm eq}\}(\boldsymbol{x},t) - \{f_i,g_i\}(\boldsymbol{x},t)\right) \\+ \left(1-\beta\right) \left(\{f_i^{*},g_i^{*}\}(\boldsymbol{x},t) - \{f_i^{\rm eq},g_i^{\rm eq}\}(\boldsymbol{x},t)\right), \end{gather}\] around \((\boldsymbol{x},t)\), leading to the following space and time-evolution equations, \[\begin{gather} \delta t\mathcal{D}_t \{f_i,g_i\} + \frac{\delta t^2}{2}{\mathcal{D}_t}^2 \{f_i,g_i\} + \mathcal{O}\left(\delta t^3\right) = 2\beta\left(\{f_i^{\rm eq},g_i^{\rm eq}\} - \{f_i,g_i\}\right) \\ + \left(1-\beta\right)\left(\{f_i^*,g_i^*\} - \{f_i^{\rm eq},g_i^{\rm eq}\}\right). \end{gather}\] Introducing the flow characteristic size and time, \(\mathcal{L}\) and \(\mathcal{T}\) the equations are made non-dimensional as, \[\begin{gather} \frac{\delta x}{\mathcal{L}}\mathcal{D}'_t \{f_i,g_i\} + \frac{\delta x^2}{2\mathcal{L}^2}{\mathcal{D}'_t}^2 \{f_i,g_i\} = 2\beta\left(\{f_i^{\rm eq},g_i^{\rm eq}\} - \{f_i,g_i\}\right) \\ + \left(1-\beta\right)\left(\{f_i^{*},g_i^{*}\} - \{f_i^{\rm eq},g_i^{\rm eq}\}\right), \end{gather}\] where, \[\mathcal{D}'_t = \frac{\mathcal{L}/\mathcal{T}}{\delta x/\delta t}\left(\partial'_t + \boldsymbol{c}'_i\cdot\boldsymbol{\nabla}'\right).\] Assuming acoustic scaling and hydrodynamic scaling, \(\varepsilon\sim\delta x/\mathcal{L} \sim\delta t/\mathcal{T}\) and dropping the primes, \[\varepsilon\mathcal{D}_t \{f_i,g_i\} + \frac{\varepsilon^2}{2}{\mathcal{D}_t}^2\{f_i,g_i\} = 2\beta\left(\{f_i^{\rm eq},g_i^{\rm eq}\} - \{f_i,g_i\}\right) + \left(1-\beta\right)\left(\{f_i^{*},g_i^{*}\} - \{f_i^{\rm eq},g_i^{\rm eq}\}\right),\] \[\begin{align} \varepsilon^0 &: \{f_i^{(0)},g_i^{(0)}\} = \{f_i^{\rm eq},g_i^{\rm eq}\},\\ \varepsilon &: \mathcal{D}_{t}^{(1)} \{f_i^{(0)},g_i^{(0)}\} = -2\beta \{f_i^{(1)},g_i^{(1)}\} + \left(1-\beta\right)\{{f^*}_i^{(1)},{g^*}_i^{(1)}\},\tag{32}\\ \varepsilon^2 &: \partial_t^{(2)}\{f_i^{(0)},g_i^{(0)}\} + \mathcal{D}_{t}^{(1)}(1-\beta) \left(\{f_i^{(1)},g_i^{(1)}\} + {\color{black}{\frac{1}{2}}} \{{f^*}_i^{(1)},{g^*}_i^{(1)}\}\right) = \nonumber \\ &-2\beta \{f_i^{(2)},g_i^{(2)}\} + \left(1-\beta\right)\{{f^*}_i^{(2)},{g^*}_i^{(2)}\}. \end{align} \tag{33}\]

The following solvability conditions apply, \[\begin{align} \sum_{i=1}^Q {f}_i^{(k)} = 0,\;\forall k>0, \tag{34}\\ \sum_{i=1}^Q \boldsymbol{c}_i {f}_i^{(1)} + \frac{1}{2}\sum_{i=1}^Q \boldsymbol{c}_i {f^*}_i^{(1)} = 0, \tag{35}\\ \sum_{i=1}^Q \boldsymbol{c}_i {f}_i^{(k)} = 0,\;\forall k>1, \tag{36}\\ \sum_{i=1}^Q {g}_i^{(1)} + \frac{1}{2} \sum_{i=1}^Q {g^*}_i^{(1)} = 0, \tag{37}\\ \sum_{i=1}^Q {g}_i^{(k)} = 0,\;\forall k>0. \tag{38} \end{align}\] Taking the zeroth-order moment of Eq. 32 , \[\partial_t^{(1)}\rho + \boldsymbol{\nabla}\cdot\rho\boldsymbol{u} = 0,\] where we have used, \[\sum_{i=1}^{Q} {f_i^*}^{(1)} = 0,\] and solvability condition 34 . For the first-order moment of 32 , \[\partial_t^{(1)}\rho\boldsymbol{u} + \boldsymbol{\nabla}\cdot\rho\boldsymbol{u}\otimes\boldsymbol{u} + \boldsymbol{\nabla} P = -2\beta\overbrace{\left(\sum_{i=1}^Q \boldsymbol{c}_if_i^{(1)}+\frac{1}{2}\boldsymbol{c}_i{f_i^*}^{(1)}\right)}^{=0} + \underbrace{\sum_{i=1}^Q \boldsymbol{c}_i{f_i^*}^{(1)}}_{\boldsymbol{F}},\] where we used 35 and further set, \[\boldsymbol{F} = -\boldsymbol{\nabla}\cdot T_K,\] where, \[T_K = \kappa \boldsymbol{\nabla}\rho\otimes\boldsymbol{\nabla}\rho - \kappa\left(\rho \boldsymbol{\nabla}^2\rho+\frac{1}{2}{\lvert \boldsymbol{\nabla}\rho\lvert}^2\right)\boldsymbol{I}.\] The force can also be shown to simplify to, \[\boldsymbol{F} = \kappa \rho \boldsymbol{\nabla} \boldsymbol{\nabla}^2\rho.\] Finally taking the zeroth-order moment for \(g_i\), \[\partial_t^{(1)}\rho E + \boldsymbol{\nabla}\cdot \boldsymbol{u}\left(\rho E + P\right) = -2\beta\overbrace{\left(\sum_{i=1}^Q g_i^{(1)}+\frac{1}{2}{g_i^*}^{(1)}\right)}^{=0} + \underbrace{\sum_{i=1}^Q {g_i^*}^{(1)}}_{\boldsymbol{u}\cdot\boldsymbol{F}}.\] Summing up balance equations at order \(\varepsilon\), \[\begin{align} \partial_t^{(1)}\rho + \boldsymbol{\nabla}\cdot \rho \boldsymbol{u} &=& 0,\tag{39}\\ \partial_t^{(1)}\rho \boldsymbol{u} + \boldsymbol{\nabla}\cdot \left(\rho \boldsymbol{u}\otimes\boldsymbol{u} + P\boldsymbol{I} \right) - \boldsymbol{F} &=& 0,\tag{40}\\ \partial_t^{(1)}\rho E + \boldsymbol{\nabla}\cdot \rho \boldsymbol{u}\left(E + P/\rho\right) - \boldsymbol{u}\cdot\boldsymbol{F} &=& 0. \tag{41} \end{align}\] The last equation Eq. 41 can be transformed into a balance equation for internal energy, using , \[\partial_t^{(1)} \mathcal{K} + \boldsymbol{\nabla}\cdot\boldsymbol{u} \mathcal{K} + \boldsymbol{u}\cdot\boldsymbol{\nabla}P - \boldsymbol{u}\cdot\boldsymbol{F} = 0,\] as, Furthermore, using \[de = c_v dT - \left(T \left(\frac{\partial P}{\partial T}\right)_\rho- P\right)\frac{d\rho}{\rho^2},\] and Eq. 39 a balance equation for temperature can be derived as, \[\partial_t^{(1)}T + \boldsymbol{u}\cdot\boldsymbol{\nabla}T + \frac{T}{\rho c_v}\left(\frac{\partial P}{\partial T}\right)_\rho\boldsymbol{\nabla}\cdot\boldsymbol{u} = 0.\] Finally, using \[dP = \left(\frac{\partial P}{\partial \rho}\right)_T d\rho + \left(\frac{\partial P}{\partial T}\right)_\rho dT,\] we can also write a balance equation for pressure as, \[\partial_t^{(1)} P + \boldsymbol{u}\cdot\boldsymbol{\nabla}P + \rho c_s^2\boldsymbol{\nabla}\cdot\boldsymbol{u} = 0,\] where, \[c_s^2 = \left(\frac{\partial P}{\partial \rho}\right)_T + \frac{T}{c_v\rho^2} \left(\frac{\partial P}{\partial T}\right)_\rho^2.\] At order \(\varepsilon^2\), the zeroth order moment leads to, \[\partial_t^{(2)} \rho = 0,\] the first order moments, \[\partial_t^{(2)}\rho\boldsymbol{u} + \boldsymbol{\nabla}\cdot\left(1-\beta\right)\left[\left(\sum_{i=1}^Q \boldsymbol{c}_i\otimes\boldsymbol{c}_i {f_i}^{(1)}\right) + \frac{1}{2}\left(\sum_{i=1}^Q \boldsymbol{c}_i\otimes\boldsymbol{c}_i {f^*}_i^{(1)}\right)\right] = 0.\] Here we can use Eq. 32 to obtain, \[\partial_t^{(2)}\rho\boldsymbol{u} + \boldsymbol{\nabla}\cdot\left(\frac{1}{2}-\frac{1}{2\beta}\right)\left(\partial_t^{(1)}\Pi_2(f^{(0)}) + \boldsymbol{\nabla}\Pi_3(f^{(0)}) - \sum_{i=1}^Q \boldsymbol{c}_i\otimes\boldsymbol{c}_i {f^*}_i^{(1)}\right) = 0.\label{eq:epsilon295mom951951}\tag{42}\] Here, \(\Pi^{(0)}_2(f)\) and \(\Pi^{(0)}_3(f)\) are the second and third order moments tensor of the equilibrium distribution function, \[\begin{align} \Pi_2(f^{(0)}) &=& \rho\boldsymbol{u}\otimes\boldsymbol{u} + P\boldsymbol{I},\\ \Pi_3(f^{(0)}) &=& \rho\boldsymbol{u}\otimes(\boldsymbol{u}\otimes\boldsymbol{u}+P\boldsymbol{I})\circ(\boldsymbol{1}-\boldsymbol{J}) + 3\rho \varsigma^2 \boldsymbol{u}\otimes\boldsymbol{I}\circ\boldsymbol{J}. \end{align}\] Using Eqs. 39 and 40 we can write, \[\begin{gather} \partial^{(1)}_t \left(\rho\boldsymbol{u}\otimes\boldsymbol{u} + P\boldsymbol{I}\right) = \boldsymbol{u}\otimes\boldsymbol{F} + \boldsymbol{F}\otimes\boldsymbol{u} -\boldsymbol{\nabla}\cdot\rho\boldsymbol{u}\otimes\boldsymbol{u}\otimes\boldsymbol{u} - \boldsymbol{\nabla}P\boldsymbol{u} - (\boldsymbol{\nabla}P\boldsymbol{u})^\dagger \\ P\left(\boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^\dagger\right) + \partial^{(1)}_t P \boldsymbol{I}. \end{gather}\] For the last term, i.e. \(\partial^{(1)}_t P\), we write a balance equation for \(P\), \[\partial_t^{(1)} P = \left(P - \rho c_s^2\right)\boldsymbol{\nabla}\cdot\boldsymbol{u} - \boldsymbol{\nabla}\cdot P\boldsymbol{u},\] while, \[\boldsymbol{\nabla}\cdot\Pi_3(f^{(0)}) = \left[\boldsymbol{\nabla}\cdot\rho\boldsymbol{u}\otimes\boldsymbol{u}\otimes\boldsymbol{u} + \boldsymbol{\nabla}P\boldsymbol{u} + \boldsymbol{\nabla}P\boldsymbol{u}^\dagger\right] + \Psi +\boldsymbol{I}\boldsymbol{\nabla}\cdot P\boldsymbol{u},\] where, \[\Psi_{\alpha\alpha} = -\partial_\alpha u_\alpha\left(u_\alpha^2 + 3(P-\rho\varsigma^2)\right).\] Adding up all terms, \[\begin{gather} \partial_t^{(1)}\Pi_2(f^{(0)}) + \boldsymbol{\nabla}\cdot\Pi_3(f^{(0)}) = \boldsymbol{u}\otimes\boldsymbol{F} + {\color{black}{\boldsymbol{F}\otimes\boldsymbol{u}}} + P\left(\boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^\dagger\right) \\+ \left(P - \rho c_s^2\right)\boldsymbol{\nabla}\cdot\boldsymbol{u}\boldsymbol{I} - \Psi. \end{gather}\] Setting, \[\sum_{i=1}^Q \boldsymbol{c}_i\otimes\boldsymbol{c}_i {f^*}_i^{(1)} = \left(\boldsymbol{u}\otimes\boldsymbol{F} + \boldsymbol{F}\otimes\boldsymbol{u} + \Psi + P\left(\frac{D+2}{D}-\frac{\rho c_s^2}{P} -\frac{\eta}{\mu}\right)\boldsymbol{\nabla}\cdot\boldsymbol{u}\boldsymbol{I}\right)\] and plugging it back into Eq. 42 , \[\partial_t^{(2)}\rho\boldsymbol{u} - \boldsymbol{\nabla}\cdot\left[ \mu\left(\boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^\dagger - \frac{2}{D}\boldsymbol{\nabla}\cdot\boldsymbol{u} \boldsymbol{I}\right) + \eta \boldsymbol{\nabla}\cdot\boldsymbol{u}\boldsymbol{I}\right] = 0,\] Where we used, \[\mu = \left(\frac{1}{2\beta}-\frac{1}{2}\right) P \delta t.\] For the second population, at order \(\varepsilon^2\), \[\partial_t^{(2)}\rho E + \boldsymbol{\nabla}\cdot\left(\frac{1}{2}-\frac{1}{2\beta}\right)\left( \partial_t^{(1)}\Pi_1(g_i^{(0)}) + \boldsymbol{\nabla}\cdot\Pi_2(g_i^{(0)}) - \sum_{i=1}^Q \boldsymbol{c}_i {g_i^*}^{(1)}\right) = 0,\] where, \[\begin{align} \Pi_1(g_i^{(0)}) &=& \boldsymbol{u}(\rho E + P),\\ \Pi_2(g_i^{(0)}) &=& \boldsymbol{u}\otimes\boldsymbol{u}\left(\rho E+2P\right) + \left(E + P/\rho\right)P\boldsymbol{I}. \end{align}\] Here we can use, \[\partial_t^{(1)}P\boldsymbol{u} + \boldsymbol{\nabla}\cdot P\boldsymbol{u}\otimes\boldsymbol{u} + \frac{P}{\rho}\boldsymbol{\nabla}P -\frac{P}{\rho}\boldsymbol{F} + \boldsymbol{u}\left(\rho c_s^2 - P\right)\boldsymbol{\nabla}\cdot\boldsymbol{u} = 0,\] and \[\partial_t^{(1)}\rho E\boldsymbol{u} + \boldsymbol{\nabla}\cdot (\rho E + P)\boldsymbol{u}\otimes\boldsymbol{u} + E\boldsymbol{\nabla}P - E\boldsymbol{F} - P\boldsymbol{u}\cdot\boldsymbol{\nabla}\boldsymbol{u} - \boldsymbol{u}(\boldsymbol{u}\cdot\boldsymbol{F})= 0,\] Adding up both contributions, \[\begin{gather} \partial_t^{(1)}\Pi_1(g_i^{(0)}) + \boldsymbol{\nabla}\cdot\Pi_2(g_i^{(0)}) = \underbrace{-2\boldsymbol{\nabla}\cdot P \boldsymbol{u}\otimes\boldsymbol{u} + 2\boldsymbol{\nabla}\cdot P \boldsymbol{u}\otimes\boldsymbol{u}}_{=0} \underbrace{+ \boldsymbol{\nabla}\cdot \rho E \boldsymbol{u}\otimes\boldsymbol{u} - \boldsymbol{\nabla}\cdot \rho E \boldsymbol{u}\otimes\boldsymbol{u}}_{=0} \\ \underbrace{+ \boldsymbol{\nabla} P(P/\rho + E) - (E+P/\rho)\boldsymbol{\nabla}P}_{P\boldsymbol{\nabla}(E+P/\rho)} + P\boldsymbol{u}\cdot\boldsymbol{\nabla}\boldsymbol{u} - \boldsymbol{u}\left(\rho c_s^2 - P\right)\boldsymbol{\nabla}\cdot\boldsymbol{u} \\+(P/\rho+E)\boldsymbol{F} + \boldsymbol{u}(\boldsymbol{u}\cdot\boldsymbol{F}). \end{gather}\] Further expanding, \[\begin{gather} \partial_t^{(1)}\Pi_1(g_i^{(0)}) + \boldsymbol{\nabla}\cdot\Pi_2(g_i^{(0)}) = P\boldsymbol{\nabla}h \overbrace{+ P\boldsymbol{\nabla}(\boldsymbol{u}^2/2) + P\boldsymbol{u}\cdot\boldsymbol{\nabla}\boldsymbol{u}}^{=P\boldsymbol{u}\cdot(\boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^\dagger)} + \boldsymbol{u}\left(P - \rho c_s^2\right)\boldsymbol{\nabla}\cdot\boldsymbol{u} \\ +(P/\rho+E)\boldsymbol{F} + \boldsymbol{u}(\boldsymbol{u}\cdot\boldsymbol{F}), \end{gather}\] where \(h = H - K\). Plugging this back into the balance equation, \[\begin{gather} \partial_t^{(2)}\rho E - \boldsymbol{\nabla}\cdot\underbrace{\boldsymbol{u}\cdot\left[\mu\left(\boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u} - \frac{2}{D}\boldsymbol{\nabla}\cdot\boldsymbol{u}\boldsymbol{I}\right) + \eta\boldsymbol{\nabla}\cdot\boldsymbol{u}\boldsymbol{I}\right]}_{-\boldsymbol{u}\cdot\boldsymbol{T}_{\rm NS}} \\ + \boldsymbol{\nabla}\cdot\left(\frac{1}{2}-\frac{1}{2\beta}\right)\left(P\boldsymbol{\nabla}h + \boldsymbol{u}\cdot P\left(\frac{D+2}{D}-\frac{\rho c_s^2}{P}-\frac{\eta }{\mu}\right)\boldsymbol{\nabla}\cdot\boldsymbol{u}\boldsymbol{I} \right. \\ \left. + (P/\rho+E)\boldsymbol{F} + \boldsymbol{u}(\boldsymbol{u}\cdot\boldsymbol{F}) - \sum_{i=1}^Q \boldsymbol{c}_i {g_i^{*}}^{(1)}\right) = 0. \end{gather}\] Using the definition of the shifted equilibrium, \[\begin{align} \Pi_0(g_i^{*(1)}) &=& \boldsymbol{u}\cdot\boldsymbol{F},\\ \Pi_1(g_i^{*(1)}) &=& \left(\boldsymbol{u}(\boldsymbol{u}\cdot\boldsymbol{F}) + \boldsymbol{F} \left(\frac{P}{\rho} + E\right) \right. \nonumber \\ & & \left. + P\left(\boldsymbol{\nabla}h - \frac{\kappa}{\mu}\boldsymbol{\nabla}T + \boldsymbol{u} \cdot\left(\frac{D+2}{D}-\frac{\rho c_s^2}{P}-\frac{\eta}{\mu}\right)\boldsymbol{\nabla}\cdot\boldsymbol{u}\boldsymbol{I}\right)\right). \end{align}\] we recover, \[\partial_t^{(2)}\rho E + \boldsymbol{\nabla}\cdot \boldsymbol{u}\cdot\boldsymbol{T}_{\rm NS} - \boldsymbol{\nabla}\cdot k\boldsymbol{\nabla}T = 0.\]

References↩︎