Central-moment-based discrete Boltzmann modeling of compressible flows

Compressible flows, characterized by both hydrodynamic and thermodynamic nonequilibrium effects, are ubiquitous in natural and engineering contexts, including applications such as inertial confinement fusion, gas pipelines, jet engines, and rocket motors [1]. Notably, modern high-speed airplanes and the jet engines that power them are wonderful examples of the application of compressible flows. Additionally, during the reentry of an aircraft into the atmosphere at high Mach numbers, the surrounding shock wave induces significant acceleration and heating of the air. Therefore, in the modern world of aerospace and mechanical engineering, a deep understanding of the principles of compressible flows is essential [2].

To predict the thermal process and nonequilibrium effects, there are three main categories of numerical methodologies. The first class involves modifying traditional macroscopic models, for example, the Navier-Stokes (NS) equations plus slip boundary conditions [3, 4, 5], the NS-type equations with multi-components and/or multi-temperatures [6, 7], Burnett and super-Burnett equations [8, 9], etc. Since the macroscopic models are constructed under the continuum assumption, the application scope is limited at small Knudsen numbers. The second type encompasses microscopic models, such as the molecular dynamics (MD) [10, 11]. The MD provides details of physical systems, but is only applicable to small temporal and spatial domains due to an excessive computational burden.

Moreover, to address the aforementioned challenges, the third category of methods encompasses mesoscopic approaches based on kinetic theory, such as the direct simulation Monte Carlo (DSMC) method [12], discrete velocity methods (DVM) [13, 14], (discrete) unified gas-kinetic schemes ((D)UGKS) [15, 16, 17, 18], the lattice Boltzmann method (LBM) [19, 20], and the discrete Boltzmann method (DBM) [21, 22], among others. These approaches act as a bridge connecting microscopic and macroscopic scales. The DSMC method, first introduced by G. A. Bird, has been extensively developed and widely applied in rarefied gas dynamics, particularly for high-speed flows [12]. However, DSMC is not well-suited to the simulation of low-speed flows, and has to balance the noise level and computational efficiency [23]. DVM represents another widely used approach, where the continuous particle velocity space is discretized into a finite set of velocity coordinate points, and numerical quadrature is employed to approximate the integration of moments [13, 14]. However, for high-speed compressible flows, especially in the near continuous flow region, the DVM exhibits limited computational efficiency [14]. Besides, in order to solve the discrete velocity Boltzmann equation, Xu et al. [15, 16] and Guo et al. [17, 18] proposed the unified gas-kinetic scheme (UGKS) and the discrete unified gas-kinetic scheme (DUGKS), respectively. In UGKS, the local integration of the discrete velocity Boltzmann equation is used to compute the flux of the distribution function, whereas in DUGKS, the flux is derived directly from the evolution equation.

Among the kinetic methodologies, the LBM stands out [19, 20]. The LBM has emerged as a competitive scheme for simulating various complex flows due to its canonical “collision-streaming” algorithm which disentangles non-linearity and non-locality and is easy to implement. Owing to these inherent advantages, LBM has been successfully applied to simulate various physical problems including multiphase [24], reactive [25], magnetohydrodynamic [26], nano- [27], biomechanics [28], and porous media flow [29]. Although the LBM has achieved significant success in simulating nearly incompressible complex flows, its application to compressible flows continues to present significant challenges.

Alexander et al. devised the first multi-speed lattice Boltzmann model containing 13 discrete velocities, representing the earliest application of the LBM to the compressible NS equations [30]. Shortly after the introduction of Alexander’s model, in 1993, Qian proposed a series of multi-speed models based on the DnQb thermal lattice Bhatnager–Gross–Krook (BGK) isothermal model for thermohydrodynamics, in which a proper internal energy is introduced and the energy equation is obtained [31]. In 1994, Chen et al. introduced a thermal lattice BGK model capable of recovering the standard compressible NS equations through a higher-order velocity expansion of the Maxwellian-type equilibrium distribution [32]. However, this model was constrained by assumptions of zero bulk viscosity and unit Prandtl number [32]. In 1997, Chen et al. constructed a two time relaxation parameters thermal lattice BGK model to achieve adjustable Prandtl number [33]. In the same year, McNamara et al. derived the distribution functions of thermal LBM using a series of moment conservation equations, enabling variability in the Prandtl number within this model [34]. Subsequently, compressible lattice Boltzmann models with a flexible specific heat ratio were introduced by Hu et al. in 1997 [35] and Yan et al. in 1999 [36]. With the aim of simulating compressible flows, Sun et al. developed an adaptive compressible LBM based on a simple delta function, where the lattice velocities vary with mean flow velocity and internal energy [37, 38, 39, 40, 41]. In 2007, Li et al. proposed a coupled double-distribution-function LBM for the NS equations with a flexible specific heat ratio and Prandtl number [42].

In recent years, Yang et al. developed a platform for constructing one-dimensional compressible LBM and subsequently extended the finite volume LBM to simulate compressible flows, including shock waves [43, 44]. Li et al. introduced a novel LBM model designed for fully compressible flows. Building upon a multi-speed model, an extra potential energy distribution function is introduced to recover the full compressible NS equations, featuring a flexible specific heat ratio and Prandtl number [45]. In 2018, Dorschner et al. constructed a particles-on-demand (PonD) method which is suitable for both incompressible and compressible flows [46]. In the PonD method, the discrete particle velocities are defined relative to a local reference frame, determined by the local flow velocity and temperature, which vary in both space and time. Subsequently, Kallikounis et al. further refined PonD method with a revised reference frame transformation using Grad’s projection to enhance stability and accuracy [47]. In 2024, Kallikounis et al. improved the PonD method for variable Prandtl number by employing the quasi-equilibrium relaxation, and the model was validated via simulations of high Mach compressible flows [48]. In the same year, Ji et al proposed an Eulerian realization of the PonD for hypersonic compressible flows, where a kinetic model formulated in the comoving reference frame [49].

Generally, lattice Boltzmann methods (LBMs) are broadly classified into the single-relaxation-time (SRT) model and the multiple-relaxation-time (MRT) model. The SRT or BGK operator is the simplest collision operator in the standard LBM, where all distribution functions relax to their local equilibrium values at a constant rate. However, the BGK-LBM may encounter stability issues in the zero-viscosity limit and/or for non-vanishing Mach numbers [20, 50, 51]. To mitigate this issue, numerous strategies have been proposed, including modifications to the numerical discretization, collision model, or both [50, 51]. In 2006, Geier et al. introduced a cascaded operator, conducting collisions in the central-moment space rather than that of raw moments as in the MRT-LBM[52]. Consequently, the corresponding method gradually interpreted as “central-moment-based” lattice Boltzmann method (CLBM) [53, 54, 55, 56, 57, 58]. However, the CLBM does not necessarily provide enhanced stability. For example, when relaxation in the central moment space is implemented using a uniform relaxation parameter for all moments, the approach reduces to a BGK collision operator with an extended equilibrium [59, 60]. To improve the stability of CLBMs, the relaxation frequencies of higher-order moments, including the trace of second-order moments, should be set to one or a value close to one [61, 62].

In this work, we constructed a central-moment-based discrete Boltzmann method (CDBM) tailored specifically for the simulation of compressible flows. The “central-moment-based” indicates that physical quantities and higher-order kinetic moments are computed through central moments. Unlike the majority of previously proposed CLBMs, which are constrained to incompressible flows and neglect the thermodynamic nonequilibrium effects inherent to the Boltzmann equation, the CDBM possesses significantly broader applicability. Based upon the Boltzmann equation, in addition to capturing the general hydrodynamic behaviors described by the NS model, the CDBM provides deeper insights into more detailed thermodynamic nonequilibrium behaviors in various complex fluids. Moreover, distinct from all existing DBMs constructed in raw-moment space [21, 22, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72], the CDBM, using the peculiar velocity, can quantify the nonequilibrium effects related to the thermal fluctuation directly. These distinctive features make the CDBM particularly valuable for investigating compressible flows that exhibit significant thermodynamic nonequilibrium effects.

The governing equations for the CDBM are expressed as follows:

where $f_{i}$ indicate the discrete distribution functions, with the subscript $i$ denoting the velocity index. The discrete velocities are given by:

where $v_{a}$, $v_{b}$, $v_{c}$ and $v_{d}$ are adjustable. Besides, in order to describe the vibrational and/or rotational energies, $I$ is introduced for extra degrees of freedom due to vibration and/or rotation, and the corresponding parameters for degrees of freedom, ${\eta}_{i}$, are defined as follows,

Figure 1 exhibits the schematic of discrete velocities.

It is worth noting that the values of these parameters are chosen based on the specific application to ensure numerical stability. In general, the magnitudes of $v_{a}$, $v_{b}$, $v_{c}$ and $v_{d}$ should approximate key physical quantities such as flow velocity $\mathbf{u}$, sound speed $v_{s}=\sqrt{\gamma T}$, and shock speed, etc. Similarly, the magnitudes of ${\eta}_{a}$, ${\eta}_{b}$, ${\eta}_{c}$ and ${\eta}_{d}$ should be around the value of $\sqrt{IT}$. This is because, in thermodynamic equilibrium, $m\bar{\eta}^{2}/2=IT/2$ where m = 1 is the particle mass, and $\bar{\eta}^{2}$ represents the average value of $\eta^{2}$, according to the equipartition of energy theorem. Hence, $\bar{\eta}=\sqrt{IT}$, and the values of ${\eta}_{a}$, ${\eta}_{b}$, ${\eta}_{c}$ and ${\eta}_{d}$ should be around $\bar{\eta}$.

The relaxation time $\tau$ is a constant which governs the relaxation speed of ${{{f}}_{i}}$ towards ${f}_{i}^{eq}$. ${{{f}}_{i}}$ and ${f}_{i}^{eq}$ represent the discrete velocity distributions and the corresponding equilibrium counterparts, respectively. The equilibrium discrete distribution functions, ${f}_{i}^{eq}$, are uniformly computed via a matrix inversion method [73]. In contrast to prior DBM constructions, this study employs central kinetic moments.

To recover the NS equations, through the Chapman-Enskog expansion, the discrete equilibrium distribution functions should satisfy the following seven central kinetic moments,

where $\Psi=1$, ${{\mathbf{v}}^{*}}$, ${{\mathbf{v}}^{*}}\cdot{{\mathbf{v}}^{*}}+{{\eta}^{2}}$, ${{\mathbf{v}}^{*}}\cdot{{\mathbf{v}}^{*}}$, $\left({{\mathbf{v}}^{*}}\cdot{{\mathbf{v}}^{*}}+{{\eta}^{2}}\right){{\mathbf{v}}^{*}}$, ${{\mathbf{v}}^{*}}{{\mathbf{v}}^{*}}{{\mathbf{v}}^{*}}$, $\left({{\mathbf{v}}^{*}}\cdot{{\mathbf{v}}^{*}}+{{\eta}^{2}}\right){{\mathbf{v}}^{*}}{{\mathbf{v}}^{*}}$, correspondingly, ${{\Psi}_{i}}=1$, $\mathbf{v}_{i}^{*}$, $\mathbf{v}_{i}^{*}\cdot\mathbf{v}_{i}^{*}+\eta_{i}^{2}$, $\mathbf{v}_{i}^{*}\cdot\mathbf{v}_{i}^{*}$, $\left(\mathbf{v}_{i}^{*}\cdot\mathbf{v}_{i}^{*}+\eta_{i}^{2}\right)\mathbf{v}_{i}^{*}$, $\mathbf{v}_{i}^{*}\mathbf{v}_{i}^{*}\mathbf{v}_{i}^{*}$, $\left(\mathbf{v}_{i}^{*}\cdot\mathbf{v}_{i}^{*}+\eta_{i}^{2}\right)\mathbf{v}_{i}^{*}\mathbf{v}_{i}^{*}$. Here ${{\mathbf{v}}^{*}}=\mathbf{v}-\mathbf{u}$ and $\mathbf{v}_{i}^{*}={{\mathbf{v}}_{i}}-\mathbf{u}$, with $\mathbf{u}$ the flow velocity. In addition, the equilibrium distribution function $f^{eq}$ is

where $D$ denotes the dimensional translational degree of freedom, $I$ stands for extra degrees of freedom due to vibration and/or rotation, and $\eta$ is used to describe the corresponding vibrational and/or rotational energies. The other parameters include $n$ as the particle number density, $T$ as the temperature, $m$ = 1 as the particle mass, and $\rho=nm$ as the mass density. Besides, the specific heat ratio is $\gamma=\left(D+I+2\right)/\left(D+I\right)$, the dynamic viscosity $\mu={\rho T}{\tau}$, kinematic viscosity $\nu={\mu}/{\rho}={T}{\tau}$, and bulk viscosity ${{\mu}_{B}}=\mu\left({2}/{D}-{2}/({D+I})\right)$.

The central kinetic moments can be expressed in a unified form

Here, ${{\bf{\bar{f}}}^{{\bf{eq}}}}={\left({\begin{matrix}{\bar{f}_{1}^{eq}}\>{\bar{f}_{2}^{eq}}\>\cdots\>{\bar{f}_{16}^{eq}}\end{matrix}}\right)^{\rm{T}}}$ and ${{\bf{f}}^{eq}}={\left({\begin{matrix}{f_{1}^{eq}}\>{f_{2}^{eq}}\>\cdots\>{f_{16}^{eq}}\end{matrix}}\right)^{\rm{T}}}$ represent the equilibrium velocity distribution functions in the central moment space and velocity space, respectively. The matrix ${\bf{M}}={{\left({{\mathbf{M}}_{1}}\;{{\mathbf{M}}_{2}}\;\cdots\;{{\mathbf{M}}_{16}}\right)}^{\mathrm{T}}}$ acts as the bridge for transforming the velocity distribution function between the moment space and the discrete velocity space, which contains the blocks, ${{\mathbf{M}}_{i}}=\left(\begin{matrix}{{M}_{i1}}\>{{M}_{i2}}\>\cdots\>{{M}_{i16}}\end{matrix}\right)$, with elements ${M}_{1i}=1$, ${M}_{2i}=v_{ix}^{*}$, ${M}_{3i}=v_{iy}^{*}$, ${M}_{4i}=v_{i}^{{*}2}+\eta_{i}^{2}$, ${M}_{5i}=v_{ix}^{{*}2}$, ${M}_{6i}=v_{ix}^{*}v_{iy}^{*}$, ${M}_{7i}=$ $v_{iy}^{{*}2}$, ${M}_{8i}=(v_{i}^{{*}2}+\eta_{i}^{2})v_{ix}^{*}$, ${M}_{9i}=(v_{i}^{{*}2}+\eta_{i}^{2})v_{iy}^{*}$, ${M}_{10i}=v_{ix}^{{*}3}$, ${M}_{11i}=v_{ix}^{{*}2}v_{iy}^{*}$, ${M}_{12i}=v_{ix}^{*}v_{iy}^{{*}2}$, ${M}_{13i}=v_{iy}^{{*}3}$, ${M}_{14i}=(v_{i}^{{*}2}+\eta_{i}^{2})v_{ix}^{{*}2}$, ${M}_{15i}=(v_{i}^{{*}2}+\eta_{i}^{2})v_{ix}^{*}v_{iy}^{*}$, ${M}_{16i}=(v_{i}^{{*}2}+\eta_{i}^{2})v_{iy}^{{*}2}$. Consequently, Eq. 6 can be expressed as

It should be mentioned that the mass, momentum and energy conservation are described by the first three kinetic moments in Eq. (4), where $f_{i}^{eq}$ can be replaced by $f_{i}$. In other words, the density, flow velocity and temperature are obtained from the kinetic moments of $f_{i}$. Replacing $f_{i}^{eq}$ with $f_{i}$ results in the imbalance in the last four kinetic moments. The difference between the results calculated by $f_{i}^{eq}$ and $f_{i}$ can be used to describe the deviation of the system from equilibrium state. Consequently, the CDBM contains the following nonequilibrium manifestations:

The second order tensor $\mathbf{\Delta}_{2}^{*}$ corresponds to the viscous stress tensor and twice the nonorganized energy. The vector $\mathbf{\Delta}_{3,1}^{*}$ is associated with the heat flux and twice the nonorganized energy flux. $\mathbf{\Delta}_{3}^{*}$ and $\mathbf{\Delta}_{4,2}^{*}$ are higher-order nonequilibrium quantities beyond traditional NS models [74]. It is crucial to note that, in comparison with previous DBMs constructed in raw-moment space, the CDBM can provide the nonequilibrium effects related to the thermal fluctuation directly.

In addition, in the basic lattice Boltzmann model, spatial and temporal discretizations are coupled, which restricts the choice of discrete velocities and also affects the construction of the discrete equilibrium distribution function. By contrast, the discrete Boltzmann method retains the use of discrete velocities but eliminates dependence on specific discretization schemes. Instead, it directly solves the continuous Boltzmann equation, allowing the CDBM to adopt spatial and temporal discretization schemes flexibly. In the subsequent simulations, the temporal derivatives are computed using the second-order Runge-Kutta scheme [75], while the spatial discretization employs the second-order non-oscillatory, non-free-parameter dissipation finite difference (NND) scheme [76]. The details are as follows,

The function ${\rm{minmod}}$ is also a type of flux limiter, defined as follows:

In addition,

where $ir$ denotes the $i$-th grid point in the $r$ direction. The NND scheme achieves second-order spatial accuracy. This scheme can suppress odd-even decoupling oscillations and effectively capture strong discontinuities.

In this section, let us verify that the CDBM is not only applicable for high-speed compressible flows, but can also capture the nonequilibrium effects accurately. To this end, five representative benchmarks are considered: the Sod shock tube and Lax shock tube validate the ability of the CDBM to capture discontinuities under varying specific heat ratios, the shock wave test demonstrates the suitability of the CDBM for hypervelocity compressible flows, the two-dimensional sound wave and decaying Taylor-Green vortex flow are used to verify the efficiency of CDBM in two-dimensional cases.

In summary, a CDBM is developed for kinetic modeling of compressible flows with both hydrodynamic and thermodynamic nonequilibrium effects. The central kinetic moments are employed for calculating the equilibrium discrete distribution functions. Conservation moments of the discrete velocity distribution function are utilized to derive macroscopic physical quantities, while higher-order central kinetic moments are employed to characterize nonequilibrium effects. Its ability to capture nonequilibrium effects is demonstrated through the simulation of the Sod shock tube. The results exhibit excellent agreement with exact solutions. Furthermore, the Lax shock tube benchmark confirms that the CDBM can accurately capture discontinuities under different specific heat ratios. Additionally, the simulation of shock waves at Mach number $\rm{Ma}=15$ verifies the capability of the CDBM for hypervelocity compressible flows. The results of two-dimensional sound wave propagation and decaying Taylor-Green vortex flow further showcase the efficiency of the CDBM in two-dimensional cases. Overall, the proposed CDBM proves to be a robust and promising tool for investigating complex compressible fluid systems, particularly those with significant thermodynamic nonequilibrium effects at the mesoscopic level.

This work is supported by Guangdong Basic and Applied Basic Research Foundation (under Grant No. 2022A1515012116), China Scholarship Council (Nos. 202306380179 and 202306380288), and Fundamental Research Funds for the Central Universities, Sun Yat-sen University (under Grant No. 24qnpy044). Support from the UK Engineering and Physical Sciences Research Council under the project “UK Consortium on Mesoscale Engineering Sciences (UKCOMES)” (Grant No. EP/X035875/1) is gratefully acknowledged. This work made use of computational support by CoSeC, the Computational Science Centre for Research Communities, through UKCOMES.