Analysis of Discrete Velocity Models for Lattice Boltzmann Simulations of Compressible Flows at Arbitrary Specific Heat Ratio

Abstract

This paper is devoted to the comparison of discrete velocity models used for simulation of compressible flows with arbitrary specific heat ratios in the lattice Boltzmann method. The stability of the governing equations is analyzed for the steady flow regime. A technique for the construction of stability domains in parametric space based on the analysis of eigenvalues is proposed. A comparison of stability domains for different models is performed. It is demonstrated that the maximum value of macrovelocity, which defines instability initiation, is dependent on the values of relaxation time, and plots of this dependence are constructed. For double-distribution-function models, it is demonstrated that the value of the Prantdl number does not seriously affect stability. The off-lattice parametric finite-difference scheme is proposed for the practical realization of the considered kinetic models. The Riemann problems and the problem of Kelvin–Helmholtz instability simulation are numerically solved. It is demonstrated that different models lead to close numerical results. The proposed technique of stability investigation can be used as an effective tool for the theoretical comparison of different kinetic models used in applications of the lattice Boltzmann method.

1. Introduction

In recent years, the lattice Boltzmann method (LBM) has been widely used for numerical simulation of different physical phenomena [1,2]. Unlike the conventional methods, which are based on discretization of the equations for macrovariables (density, velocity, energy, etc.), in the LBM, the systems of fully discrete kinetic equations (so-called lattice Boltzmann equations, LBEs) are considered.

Based on LBEs, the computational schemes deal with the intrinsic coupling of time and space steps. It leads to the fixed value of the Courant–Friedrichs–Lewi (CFL) number. For LBEs, the CFL number is equal to unity, and it is a constraint stated on the values of grid steps. It brings a limitation to the standard LBM, because it leads to the impossibility of using variable time and space steps. This problem is relevant, for example, for the simulation of gas flows, where shock and rarefication waves take place [3]. Such phenomena can lead to the necessity of the application of unstructured and adaptive grids and variable time steps. To overcome this problem, so-called off-lattice schemes, which are constructed using finite-difference [4], finite-volume [5] and finite-element [6] methods, are proposed. These approaches are based on the independent discretization of continuous kinetic equations with discrete velocities on time and space variables.

Most existing schemes, based on LBEs, are restricted to weakly compressible simulations. However, for highly compressible flows, such schemes lead to some difficulties. One of the main reasons is that the equilibrium distribution functions used in the standard LBM are derived from the truncated Taylor expansions of Maxwellian distributions with the assumption of a low or moderate Mach number.

The straightforward extension of the LBM to the simulation of compressible flows can be achieved by increasing the number of discrete velocities and using appropriate equilibrium distributions. This leads to the commonly known multi-speed models [7]. The standard approach for LBEs, based on collide-and-stream techniques, is less computationally expensive than off-lattice schemes. However, its application to multi-speed lattices leads to numerical instabilities [8]. To avoid this problem, for many multi-speed models, the off-lattice approach is used [9,10].

Different compressible models are constructed for the opportunities to vary specific heat ratio, and Mach and Prandtl numbers. For the case of the fixed specific heat ratio, Yan et al. [11] proposed a compressible model with more than one rest energy level. In [12], the LB model for the compressible Navier–Stokes system with an adjustable specific heat ratio is constructed by adding a new molecular velocity. Watari [13] constructs the kinetic equations with discrete velocities and additional degrees of freedom for 2D Navier–Stokes and 3D Euler models with flexible specific heat ratios and Prandtl numbers. Nie et al. [14] propose an approach for treating the internal degrees of freedom of the polyatomic gas molecules. As it is mentioned above, traditional equilibrium distribution functions are based on the truncated Taylor expansion of the Maxwellian distribution, limited to low and moderate Mach number [3]. Authors of [15,16] proposed the general approach equilibrium distribution functions constructing, based on the expansions on generalized Hermite polynomials presented in [17]. Qu et al. [18] obtained equilibrium distributions based on the circular function, which is not related to the Maxwellian distribution. In [19], a general technique for the construction of equilibrium distributions for kinetic models of high-speed compressible flows is proposed. The approach is based on the solution of a linear algebraic system based on the equations for moments of distribution functions for the selected lattice and velocity set. In [20], the control of the Prandtl number is realized by the modification of the collision term. In [21], the multiple-relaxation-time scheme with a flexible value of this number is constructed. Zhang et al. [22] tried to control the Prandtl number in the case of Burnett equations using an ellipsoidal statistical model.

Another class of discrete velocity models used for the computation of compressible flows is the double distribution function (DDF) models. In addition to a set of distribution functions for the computation of hydrodynamic variables (density, velocity), a set of functions for energy (or temperature) is introduced. One of the first models where functions for temperature are introduced is proposed in the paper of He et al. [23]. Li et al. [24] propose a coupled DDF-LB model for the compressible Navier–Stokes equations in order to achieve a flexible specific heat ratio. In [25], this model is applied for the simulation of high-speed compressible viscous flows with a boundary layer. In [26], this model is used for the simulation of shock wave propagation and its interactions with the boundary layer in a shock tube. As for the multi-speed models, in DDF models at high Mach numbers, the traditional collide-and-stream method, typical for the standard LBM, is not used due to the numerical instabilities (for example, see works [27,28,29], where the off-lattice schemes are used). Saadat et al. [30] propose a DDF model for the D2Q9 lattice with a correction of the collision term and a shifted lattice, proposed in [31] for the simulation of flows with high Mach numbers. Qiu et al. [32] constructed a DDF model based on the D2Q17 lattice for the simulation of viscous compressible flows.

In the literature, alternative approaches for compressible flow simulations, distinguished from multi-speed and DDF models, are used. The typical examples are the hybrid LBM, which uses both kinetic equations and equations for macrovariables [33]; the pressure-based regularized LBM, which uses explicit correction terms for the pressure [34]; hybrid recursive regularized versions of the LBM [35,36,37]; and multiple-relaxation-time (MRT) models (e.g., see [38]).

The presented paper is devoted to the analysis of discrete velocity models used in the LBM for the simulation of compressible flows and for construction of off-lattice computational schemes for practical simulations. In the theoretical aspect, attention is focused on the stability analysis of the solutions associated with the steady flow regimes. Multi-speed and DDF models, widely used in literature, are considered. The analysis is based on the stability criterion, related to the negativity of the real parts of the eigenvalues of matrices of systems, obtained after linearization of the governing equations and separation of variables. The stability domains in the parameter space are constructed and analyzed. It is demonstrated that some models can have larger stability domains than others. Selected multi-speed and DDF models in practice are compared on the numerical solutions of nonlinear problems, such as the Riemann problem and the Kelvin–Helmholtz instability (KHI) simulation. For the numerical solution, the parametric off-lattice finite-difference scheme is used. It is demonstrated that the numerical solutions obtained by different models are close to each other.

The following new results are obtained in the presented work:

1. The numerical procedure for stability analysis, based on the construction of stability domains in parametric space, is proposed.

2. The comparison of stability domains, constructed by the proposed procedure, is performed.

3. It is demonstrated that the value of the Prandtl number does not seriously affect the stability of DDF models.

4. The models are compared on the nonlinear problems of gas dynamics: Riemann problems and simulation of KHI. It is demonstrated that models lead to close numerical solutions.

The paper has the following structure. Section 2 presents the discrete velocity models considered in the paper. Section 3 is devoted to the stability analysis of the governing equations of these models. The approach of constructing stability domains in parametric space is proposed. The domains, constructed for different models, are compared. Section 4 is devoted to the numerical experiments. A parametric off-lattice finite-difference scheme for the numerical simulations is proposed. Selected nonlinear problems in gas dynamics are solved numerically. Some concluding remarks are made in Section 5.

2. Kinetic Models with Discrete Velocities

The Boltzmann equation, which describes gas dynamics in the case of an arbitrary value of the Knudsen number, is written as follows:

$$\frac{\partial f}{\partial t}+\mathbf{v}·\nabla f+\frac{\mathbf{F}}{m}·{\nabla }_{\mathbf{v}}f=I\left(f\right),$$

(1)

where term $I\left(f\right)$ can be treated as the model of particle interactions.

The macrovariables $\rho (t,\mathbf{r})$ and $\mathbf{u}(t,\mathbf{r})$ are computed as moments of the distribution function:

$$\rho (t,\mathbf{r})=\int f(t,\mathbf{r},\mathbf{v})d\mathbf{v},\mathbf{u}(t,\mathbf{r})=\frac{1}{\rho }\int \mathbf{v}f(t,\mathbf{r},\mathbf{v})d\mathbf{v}.$$

The exact solution of Equation (1), known as a Maxwellian equilibrium distribution function, is written as [39]

$$f^{\left(eq\right)}(t,\mathbf{r},\mathbf{v})=\frac{\rho }{m}{\left(\frac{m}{2\pi \mathrm{k}T}\right)}^{\frac{3}{2}}\exp \left(-\frac{m{(\mathbf{v}-\mathbf{u})}^2}{2\mathrm{k}T}\right).$$

(2)

In the presented paper, Equation (1) is considered in the absence of body force action ($\mathbf{F}=\mathbf{0}$) and the Bhatnagar–Gross-Krook model of collision term (see [40]) is considered:

$$\frac{\partial f}{\partial t}+\mathbf{v}·\nabla f=-\frac{f-f^{\left(eq\right)}}{\lambda }.$$

(3)

One of the main methods of numerical solution of problems for the Boltzmann equation is the so-called discrete velocities method. According to this approach, a finite set of velocities is introduced: ${\mathbf{v}}_i=v{\mathbf{e}}_i$, $i=\overline{1,n}$, where $v=l/\delta t$ and ${\mathbf{e}}_i$ are the given vectors. After considering Equation (3) on this set of velocities, the following system of kinetic equations with discrete velocities is obtained:

$$\frac{\partial f_i}{\partial t}+{\mathbf{v}}_i·\nabla f_i=-\frac{f_i-f_i^{\left(eq\right)}}{\lambda },$$

(4)

where $f_i=f_i(t,\mathbf{r})=f(t,\mathbf{r},{\mathbf{v}}_i)$.

The specific model with discrete velocities for simulation of compressible gas flow is defined by the specific velocity set and by the expression, which approximates Equation (2). In the presented paper, we consider four models widely used in the literature to simulate highly compressible flows. These models are constructed for the following equation-of-state: $p=(\gamma -1)\rho e$. We focus our attention on the models constructed for 2D flows. We chose two typical examples from the class of multi-speed models and two from the class of DDF models. Considered models are widely used in many works on LB simulation of compressible flows.

2.1. Multispeed Model, Based on D2Q16 Lattice

This model is presented by the system of kinetic Equation (4). The D2Q16 lattice, represented by the following velocity set (see Figure 1, a for schematic representation of velocity set) is used:

$${\mathbf{v}}_1=v\left(1,0\right),{\mathbf{v}}_2=v\left(0,1\right),{\mathbf{v}}_3=v\left(-1,0\right),{\mathbf{v}}_4=v\left(0,-1\right),$$

$${\mathbf{v}}_5=v\left(1,1\right),{\mathbf{v}}_6=v\left(-1,1\right),{\mathbf{v}}_7=v\left(-1,-1\right),{\mathbf{v}}_8=v\left(1,-1\right),$$

$${\mathbf{v}}_9=2v\left(1,0\right),{\mathbf{v}}_{10}=2v\left(0,1\right),{\mathbf{v}}_{11}=2v\left(-1,0\right),{\mathbf{v}}_{12}=2v\left(0,-1\right),$$

$${\mathbf{v}}_{13}=2v\left(1,1\right),{\mathbf{v}}_{14}=2v\left(-1,1\right),{\mathbf{v}}_{15}=2v\left(-1,-1\right),{\mathbf{v}}_{16}=2v\left(1,-1\right).$$

The macrovariables are computed as follows:

$$\rho =\sum _i{f}_i,\rho \mathbf{u}=\sum _i{f}_i{\mathbf{v}}_i,\rho \left(e+\frac{u^2}{2}\right)=\sum _i{f}_i\left(\frac{{\mathbf{v}}_i^2}{2}+\frac{{\eta }_i^2}{2}\right),$$

where ${\eta }_i$ are the free parameters introduced to describe the extra degrees of freedom, corresponding to the molecular rotation and vibration: ${\eta }_i={\eta }_0\ne 0$, $i=\overline{1,4}$, ${\eta }_i=0$, $i=\overline{5,16}$.

The vector of equilibrium distributions is computed as: ${\mathbf{f}}^{\mathbf{eq}}={\mathbf{C}}^{-1}\mathbf{M}$, where $\mathbf{M}$ is the vector of moments and $\mathbf{C}$ is a square matrix, presented in

Table 1. The expressions for component of vector $\mathbf{M}$ and matrix $\mathbf{C}$ [19].

$M_1$	$\rho$	$C_{i1}$	1
$M_2$	$\rho u_x$	$C_{i2}$	${\mathrm{v}}_{ix}$
$M_3$	$\rho u_y$	$C_{i3}$	${\mathrm{v}}_{iy}$
$M_4$	$\rho \left(bT+u_x^2+u_y^2\right)$	$C_{i4}$	${\mathrm{v}}_{ix}^2+{\mathrm{v}}_{iy}^2+{\eta }_i^2$
$M_5$	$\rho u_x{u}_y$	$C_{i5}$	${\mathrm{v}}_{ix}{\mathrm{v}}_{iy}$
$M_6$	$\rho u_x^2+P$	$C_{i6}$	${\mathrm{v}}_{ix}^2$
$M_7$	$\rho u_y^2+P$	$C_{i7}$	${\mathrm{v}}_{iy}^2$
$M_8$	$\rho u_x\left(\left(b+2\right)T+u_x^2+u_y^2\right)$	$C_{i8}$	${\mathrm{v}}_{ix}\left({\mathrm{v}}_{ix}^2+{\mathrm{v}}_{iy}^2+{\eta }_i^2\right)$
$M_9$	$\rho u_y\left(\left(b+2\right)T+u_x^2+u_y^2\right)$	$C_{i9}$	${\mathrm{v}}_{iy}\left({\mathrm{v}}_{ix}^2+{\mathrm{v}}_{iy}^2+{\eta }_i^2\right)$
$M_{10}$	$\rho u_x\left(3T+u_x^2\right)$	$C_{i10}$	${\mathrm{v}}_{ix}^3$
$M_{11}$	$\rho u_y\left(3T+u_y^2\right)$	$C_{i11}$	${\mathrm{v}}_{iy}^3$
$M_{12}$	$\rho u_y\left(T+u_x^2\right)$	$C_{i12}$	${\mathrm{v}}_{ix}^2{\mathrm{v}}_{iy}$
$M_{13}$	$\rho u_x\left(T+u_y^2\right)$	$C_{i13}$	${\mathrm{v}}_{ix}{\mathrm{v}}_{iy}^2$
$M_{14}$	$$\begin{gathered} \left(b+4\right)\rho T{u}_x{u}_y+ \\ \rho u_x{u}_y\left(u_x^2+u_y^2\right) \end{gathered}$$	$C_{i14}$	${\mathrm{v}}_{ix}{\mathrm{v}}_{iy}\left({\mathrm{v}}_{ix}^2+{\mathrm{v}}_{iy}^2+{\eta }_i^2\right)$
$M_{15}$	$\begin{array}{c}\left(b+2\right)\rho T^2+\\ \left(\left(b+4\right)u_x^2+u_x^2+u_y^2\right)\rho T+\\ \rho u_x^2\left(u_x^2+u_y^2\right)\end{array}$	$C_{i15}$	${\mathrm{v}}_{ix}^2\left({\mathrm{v}}_{ix}^2+{\mathrm{v}}_{iy}^2+{\eta }_i^2\right)$
$M_{16}$	$\begin{array}{c}\left(b+2\right)\rho T^2+\\ \left(\left(b+4\right)u_y^2+u_x^2+u_y^2\right)\rho T+\\ \rho u_y^2\left(u_x^2+u_y^2\right)\end{array}$	$C_{i16}$	${\mathrm{v}}_{iy}^2\left({\mathrm{v}}_{ix}^2+{\mathrm{v}}_{iy}^2+{\eta }_i^2\right)$

. The value of $\gamma$ is incorporated into the model in the following way: the parameter $b=2/(\gamma -1)$ is substituted in the expressions for $f_i^{\left(eq\right)}$.

2.2. Multispeed Model, Based on D2Q8 Lattice

This model is based on the following system:

$$\frac{\partial f_{ki}}{\partial t}+{\mathbf{v}}_{ki}·\nabla f_{ki}=-\frac{1}{\lambda }(f_{ki}-f_{ki}^{\left(eq\right)}),$$

(5)

where in the 2D case index $k=\overline{0,8}$, is related to the groups of particles with extra ($k=\overline{1,4}$) and without degrees of freedom ($k=\overline{5,8}$), with corresponding velocities, $i=1$ for $k=0$, $i=\overline{1,8}$ if $k>0$. So, the D2Q8 lattice is used for the model (see Figure 1b for schematic representation of velocity set):

$${\mathbf{v}}_0=0,{\mathbf{v}}_{ki}=\left({\mathrm{v}}_{kix},{\mathrm{v}}_{kiy}\right)={\mathit{v}}_k\left[cos\left(\frac{2\pi i}{N}\right),sin\left(\frac{2\pi i}{N}\right)\right],i=\overline{1,8},$$

where $v_5=v_1$, $v_6=v_2$, $v_7=v_3$, $v_8=v_4$, $v_1=v$, $v_2=2v_1$, $v_3=3v_1$ and $v_4=4v_1$. The system (4) consists of 65 equations. The values of ${\eta }_k$ are defined as for the model from [19], but ${\eta }_k\ne 0$, $k=\overline{1,4}$, ${\eta }_k=0$ and $k=\overline{5,8}$.

The equilibrium distributions are computed as [13]

$$\begin{array}{c}\hfill f_{ki}^{\left(eq\right)}=\rho F_k\left(\left(1-\frac{u^2}{2e}+\frac{u^4}{8e^2}\right)+\left(1-\frac{u^2}{2e}\right)\frac{{\mathrm{v}}_{ki\alpha }{u}_{\alpha }}{e}+\left(1-\frac{u^2}{2e}\right)\frac{{\mathrm{v}}_{ki\alpha }{\mathrm{v}}_{ki\beta }{u}_{\alpha }{u}_{\beta }}{2e^2}+\right.\\ \hfill \left.\frac{{\mathrm{v}}_{ki\alpha }{\mathrm{v}}_{ki\beta }{\mathrm{v}}_{ki\gamma }{u}_{\alpha }{u}_{\beta }{u}_{\gamma }}{6e^3}+\frac{{\mathrm{v}}_{ki\alpha }{\mathrm{v}}_{ki\beta }{\mathrm{v}}_{ki\gamma }{\mathrm{v}}_{ki\xi }{u}_{\alpha }{u}_{\beta }{u}_{\gamma }{u}_{\xi }}{24e^4}\right),\end{array}$$

where the weights $F_k$ are obtained in such a way, that the moment of these functions are equal to the moments of distribution functions:

$$F_0=1-8\left(F_1+F_2+F_3+F_4+F_5+F_6+F_7+F_8\right),$$

$$F_1=\frac{n}{{\eta }_1^2}\frac{6e^4-\left(v_2^2+v_3^2+v_4^2\right)e^3+\frac{1}{4}\left(v_2^2{v}_3^2+v_3^2{v}_4^2+v_4^2{v}_2^2\right)e^2-\frac{1}{8}{v}_2^2{v}_3^2{v}_4^{2}e}{\left(v_1^2-v_3^2\right)\left(v_1^2-v_3^2\right)\left(v_1^2-v_4^2\right)},$$

$$F_2=\frac{n}{{\eta }_2^2}\frac{6e^4-\left(v_3^2+v_4^2+v_1^2\right)e^3+\frac{1}{4}\left(v_3^2{v}_4^2+v_4^2{v}_1^2+v_1^2{v}_3^2\right)e^2-\frac{1}{8}{v}_3^2{v}_4^2{v}_1^{2}e}{\left(v_2^2-v_3^2\right)\left(v_2^2-v_4^2\right)\left(v_2^2-v_1^2\right)},$$

$$F_3=\frac{n}{{\eta }_3^2}\frac{6e^4-\left(v_4^2+v_1^2+v_2^2\right)e^3+\frac{1}{4}\left(v_4^2{v}_1^2+v_1^2{v}_2^2+v_2^2{v}_4^2\right)e^2-\frac{1}{8}{v}_4^2{v}_1^2{v}_2^{2}e}{\left(v_3^2-v_4^2\right)\left(v_3^2-v_1^2\right)\left(v_3^2-v_2^2\right)},$$

$$F_4=\frac{n}{{\eta }_4^2}\frac{6e^4-\left(v_1^2+v_2^2+v_3^2\right)e^3+\frac{1}{4}\left(v_1^2{v}_2^2+v_2^2{v}_3^2+v_3^2{v}_1^2\right)e^2-\frac{1}{8}{v}_3^2{v}_2^2{v}_3^{2}e}{\left(v_4^2-v_1^2\right)\left(v_4^2-v_2^2\right)\left(v_4^2-v_3^2\right)},$$

$$F_5=-F_1+\frac{48e^4-6\left(v_2^2+v_3^2+v_4^2\right)e^3+\left(v_2^2{v}_3^2+v_3^2{v}_4^2+v_4^2{v}_2^2\right)e^2-\frac{1}{8}{v}_2^2{v}_3^2{v}_4^{2}e}{v_1^2\left(v_1^2-v_2^2\right)\left(v_1^2-v_3^2\right)\left(v_1^2-v_4^2\right)},$$

$$F_6=-F_2+\frac{48e^4-6\left(v_3^2+v_4^2+v_1^2\right)e^3+\left(v_3^2{v}_4^2+v_4^2{v}_1^2+v_1^2{v}_3^2\right)e^2-\frac{1}{4}{v}_3^2{v}_4^2{v}_1^{2}e}{v_2^2\left(v_2^2-v_3^2\right)\left(v_2^2-v_4^2\right)\left(v_2^2-v_1^2\right)},$$

$$F_7=-F_3+\frac{48e^4-6\left(v_4^2+v_1^2+v_2^2\right)e^3+\left(v_4^2{v}_1^2+v_1^2{v}_2^2+v_2^2{v}_4^2\right)e^2-\frac{1}{4}{v}_4^2{v}_1^2{v}_2^{2}e}{v_3^2\left(v_3^2-v_4^2\right)\left(v_3^2-v_1^2\right)\left(v_3^2-v_2^2\right)},$$

$$F_8=-F_4+\frac{48e^4-6\left(v_1^2+v_2^2+v_3^2\right)e^3+\left(v_1^2{v}_2^2+v_2^2{v}_3^2+v_3^2{v}_1^2\right)e^2-\frac{1}{4}{v}_1^2{v}_2^2{v}_3^{2}e}{v_4^2\left(v_4^2-v_1^2\right)\left(v_4^2-v_2^2\right)\left(v_4^2-v_3^2\right)}.$$

The macrovariables are computed as follows:

$$\rho =\sum _{ki}{f}_{ki},\rho \mathbf{u}=\sum _{ki}{f}_{ki}{\mathbf{v}}_{ki},\rho \left(\frac{D+\sigma }{\sigma }e+\frac{u^2}{2}\right)=\sum _{ki}{f}_{ki}\left(\frac{v_{ki}^2}{2}+\frac{{\eta }_k^2}{2}\right),$$

where $D=2$, the energy level, which corresponds to the extra degree of freedom, is defined as ${\eta }_k^2/2$.

Specific heat ratio is computed as follows:

$$\gamma =\frac{D+\sigma +2}{D+\sigma }.$$

2.3. DDF Model, Based on D2Q13 Lattice

The model is presented by the following system:

$$\frac{\partial f_i}{\partial t}+{\mathbf{v}}_i·\nabla f_i=-\frac{1}{{\lambda }_f}\left(f_i-f_i^{eq}\right),$$

(6)

$$\frac{\partial h_i}{\partial t}+{\mathbf{v}}_i·\nabla h_i=-\frac{1}{{\lambda }_h}\left(h_i-h_i^{eq}\right)+\frac{1}{{\lambda }_{hf}}\left({\mathbf{v}}_i·\mathbf{u}\right)\left(f_i-f_i^{\left(eq\right)}\right),i=\overline{0,n},$$

(7)

where $f_i$ are the density distribution functions, $1/{\lambda }_{hf}=1/{\lambda }_h-1/{\lambda }_f$. For this model, the D2Q13 lattice ($n=12$) is used (see Figure 1c for schematic representation of velocity set):

$${\mathbf{v}}_0=(0,0),{\mathbf{v}}_1=v\left(1,0\right),{\mathbf{v}}_2=v\left(0,1\right),{\mathbf{v}}_3=v\left(-1,0\right),{\mathbf{v}}_4=v\left(0,-1\right),$$

$${\mathbf{v}}_5=v\left(1,1\right),{\mathbf{v}}_6=v\left(-1,1\right),{\mathbf{v}}_7=v\left(-1,-1\right),{\mathbf{v}}_8=v\left(1,-1\right),$$

$${\mathbf{v}}_9=2v\left(1,0\right),{\mathbf{v}}_{10}=2v\left(0,1\right),{\mathbf{v}}_{11}=2v\left(-1,0\right),{\mathbf{v}}_{12}=2v\left(0,-1\right).$$

Equilibrium distributions $f_i^{\left(eq\right)}$ are computed as follows:

$$f_0^{\left(eq\right)}\left(u_x,u_y\right)=\frac{\rho }{4}\left(u_x^4+5p^2-10p+4+4u_x^2{u}_y^2+u_y^4+\left(10p-5\right)\left(u_x^2+u_y^2\right)\right),$$

$$f_1^{\left(eq\right)}\left(u_x,u_y\right)=-\frac{\rho }{6}\left(-4u_x^2+3p^2+u_x^4-4p+3u_x{u}_y^2+9p{u}_x^2+6p{u}_x+3u_x^2{u}_y^2-u_x+u_x^3\right),$$

$$f_5^{\left(eq\right)}\left(u_x,u_y\right)=\frac{\rho }{4}\left(u_x{u}_y^2+u_x{u}_y+p{u}_x+p{u}_y+u_x^2{u}_y+0.5p^2+u_x^2{u}_y^2+p{u}_x^2+p{u}_y^2\right),$$

$$f_9^{\left(eq\right)}\left(u_x,u_y\right)=\frac{\rho }{24}\left(-2u_x+u_x^4-p-u_x^2+6p{u}_x^2+1.5p^2+2u_x^3+6p{u}_x\right),$$

$$f_2^{\left(eq\right)}\left(u_x,u_y\right)=f_1^{\left(eq\right)}\left(u_y,u_x\right),f_6^{\left(eq\right)}\left(u_x,u_y\right)=f_5^{\left(eq\right)}\left(-u_x,u_y\right),$$

$$f_3^{\left(eq\right)}\left(u_x,u_y\right)=f_1^{\left(eq\right)}\left(-u_x,u_y\right),f_7^{\left(eq\right)}\left(u_x,u_y\right)=f_5^{\left(eq\right)}\left(-u_x,-u_y\right),$$

$$f_4^{\left(eq\right)}\left(u_x,u_y\right)=f_1^{\left(eq\right)}\left(-u_y,u_x\right),f_8^{\left(eq\right)}\left(u_x,u_y\right)=f_5^{\left(eq\right)}\left(u_x,-u_y\right),$$

$$f_{10}^{\left(eq\right)}\left(u_x,u_y\right)=f_9^{\left(eq\right)}\left(u_y,u_x\right),f_{11}^{\left(eq\right)}\left(u_x,u_y\right)=f_9^{\left(eq\right)}\left(-u_x,u_y\right),$$

$$f_{12}^{\left(eq\right)}\left(u_x,u_y\right)=f_9^{\left(eq\right)}\left(-u_y,u_x\right).$$

Equilibrium distributions for total energy $h_i^{\left(eq\right)}$ are computed as follows:

$$h_i^{\left(eq\right)}=\left(E+\left({\mathbf{v}}_i-\mathbf{u}\right)·\mathbf{u}\right)f_i^{\left(eq\right)}+{\omega }_i\frac{p}{v_1^2}RT,$$

where ${\omega }_0=0$, ${\omega }_i=-\frac{1}{3}$$i=\overline{1,4}$, ${\omega }_i=\frac{1}{4}$$i=\overline{5,8}$, ${\omega }_i=\frac{1}{12}$$i=\overline{9,12}$ and $E=e+\frac{1}{2}{u}^2$ is the total energy.

The system (6)–(7) consists of 26 equations. The macrovariables are computed as follows:

$$\rho =\sum _i{f}_i,\rho \mathbf{u}=\sum _i{f}_i{\mathbf{v}}_i,\rho \left(e+\frac{u^2}{2}\right)=\sum _i{h}_i\frac{v_i^2}{2}.$$

2.4. DDF Model, Based on D2Q17 Lattice

The model is based on system (6)–(7) in the case of D2Q17 lattice (see Figure 1d for schematic representation of velocity set):

$${\mathbf{v}}_0=\left(0,0\right),{\mathbf{v}}_1=\frac{2}{3}v\left(1,0\right),{\mathbf{v}}_2=\frac{2}{3}v\left(0,1\right),{\mathbf{v}}_3=\frac{2}{3}v\left(-1,0\right),$$

$${\mathbf{v}}_4=\frac{2}{3}v\left(0,-1\right),{\mathbf{v}}_5=\frac{2}{3}v\left(2,1\right),{\mathbf{v}}_6=\frac{2}{3}v\left(1,2\right),{\mathbf{v}}_7=\frac{2}{3}v\left(-1,2\right),$$

$${\mathbf{v}}_8=\frac{2}{3}v\left(-2,1\right),{\mathbf{v}}_9=\frac{2}{3}v\left(-2,-1\right),{\mathbf{v}}_{10}=\frac{2}{3}v\left(-1,-2\right),{\mathbf{v}}_{11}=\frac{2}{3}v\left(1,-2\right),$$

$${\mathbf{v}}_{12}=\frac{2}{3}v\left(2,-1\right),{\mathbf{v}}_{13}=\frac{2}{3}v\left(3,0\right),$$

$${\mathbf{v}}_{14}=\frac{2}{3}v\left(0,3\right),{\mathbf{v}}_{15}=\frac{2}{3}v\left(-3,0\right),{\mathbf{v}}_{16}=\frac{2}{3}v\left(0,-3\right).$$

The presented model consists of 34 equations.

Equilibrium distributions $f_i^{\left(eq\right)}$ are computed as follows:

$$f_0^{\left(eq\right)}\left(u_x,u_y\right)=\frac{1}{16}\rho (6u_x^4+33u_x^4{u}_y^2+6u_y^2{u}_x+69\theta u_x^2+69\theta u_y^2+69{\theta }^2-$$

$$28u_x^2-28u_y^2-56\theta +16),$$

$$f_1^{\left(eq\right)}\left(u_x,u_y\right)=\frac{3}{1280}\rho (27u_x^5+27u_y^5+270\theta u_x^3+270\theta u_y^3-60u_x^4-180u_x^3{u}_y-420u_x^2{u}_y^2-$$

$$180u_x{u}_y^3-60u_y^4+405{\theta }^2{u}_x+405{\theta }^5{u}_y-780\theta u_x^2-1080\theta u_x{u}_y-780\theta u_y^2-180u_x^3-$$

$$120u_x^2{u}_y-120u_x{u}_y^2-180u_y^3-780{\theta }^2-660\theta u_x-660\theta u_y+240u_x^2+400u_x{u}_y+240u_y^2+$$

$$480\theta +288u_x+288u_y),$$

$$\begin{array}{c}\hfill f_5^{\left(eq\right)}\left(u_x,u_y\right)=-\frac{3}{5120}\rho \left(108u_x^5-27u_y^5+1080\theta u_x^3-270\theta u_y^3+60u_x^4-360u_x^3{u}_y-480u_x^2{u}_y^2-\right.\\ \hfill 120u_y^4+1620{\theta }^2{u}_x-405{\theta }^2{u}_y-120\theta u_x^2-1080\theta u_x{u}_y-1200\theta u_y^2-480u_x^3-480u_x^2{u}_y+\\ \hfill \left.60u_y^3-660{\theta }^2-1440\theta u_x-300\theta u_y-240u_x^2+160u_x{u}_y+480u_y^2+240\theta +192u_x+192u_y\right),\end{array}$$

$$\begin{array}{c}\hfill f_{13}^{\left(eq\right)}\left(u_x,u_y\right)=\frac{1}{2560}\rho \left(81u_x^5+810\theta u_x^3+120u_x^4-240u_x^2{u}_y^2-60u_y^4+1215{\theta }^2{u}_x+480\theta u_x^2-\right.\\ \hfill \left.600\theta u_y^2-180u_x^3-60{\theta }^2-540\theta u_x-160u_x^2+240u_y^2+80\theta +64u_x\right),\end{array}$$

$$f_2^{\left(eq\right)}\left(u_x,u_y\right)=f_1^{\left(eq\right)}\left(-u_x,u_y\right),f_6^{\left(eq\right)}\left(u_x,u_y\right)=f_5^{\left(eq\right)}\left(u_y,u_x\right),$$

$$f_3^{\left(eq\right)}\left(u_x,u_y\right)=f_1^{\left(eq\right)}\left(-u_x,-u_y\right),f_7^{\left(eq\right)}\left(u_x,u_y\right)=f_5^{\left(eq\right)}\left(u_y,-u_x\right),$$

$$f_4^{\left(eq\right)}\left(u_x,u_y\right)=f_1^{\left(eq\right)}\left(u_x,-u_y\right),f_8^{\left(eq\right)}\left(u_x,u_y\right)=f_5^{\left(eq\right)}\left(-u_x,u_y\right),$$

$$f_{14}^{\left(eq\right)}\left(u_x,u_y\right)=f_{13}^{\left(eq\right)}\left(u_y,u_x\right),f_9^{\left(eq\right)}\left(u_x,u_y\right)=f_5^{\left(eq\right)}\left(-u_x,-u_y\right),$$

$$f_{15}^{\left(eq\right)}\left(u_x,u_y\right)=f_{13}^{\left(eq\right)}\left(-u_x,u_y\right),f_{10}^{\left(eq\right)}\left(u_x,u_y\right)=f_5^{\left(eq\right)}\left(-u_y,-u_x\right),$$

$$f_{16}^{\left(eq\right)}\left(u_x,u_y\right)=f_{13}^{\left(eq\right)}\left(-u_y,u_x\right),f_{11}^{\left(eq\right)}\left(u_x,u_y\right)=f_5^{\left(eq\right)}\left(-u_y,u_x\right),$$

$$f_{12}^{\left(eq\right)}\left(u_x,u_y\right)=f_5^{\left(eq\right)}\left(u_x,-u_y\right),$$

where $\theta =T/T_c$. Functions $h_i^{\left(eq\right)}$ are presented as follows:

$$h_i^{\left(eq\right)}=\left(E+\left({\mathbf{v}}_i-\mathbf{u}\right)·\mathbf{u}\right)f_i^{\left(eq\right)}+{\omega }_i\frac{p}{v^2}RT,$$

where ${\omega }_0=0$, ${\omega }_i=-\frac{17}{56}$$i=\overline{1,4}$, ${\omega }_i=\frac{1}{8}$$i=\overline{5,12}$ and ${\omega }_i=\frac{3}{56}$$i=\overline{13,16}$. Macrovariables are computed in the same way as the model from [24].

3. Stability Analysis

Let the following dimensionless variables be introduced:

$$\tilde{t}=\frac{t}{\delta t},\tilde{\mathbf{r}}=\frac{\mathbf{r}}{l},{\tilde{f}}_i=\frac{f_i}{\mathsf{\Phi }},{\tilde{h}}_i=\frac{h_i}{H},$$

(8)

where $\mathsf{\Phi }$ and H are the typical values of $f_i$ and $h_i$.

After the substitution of (8) into (4), the dimensionless system is obtained:

$$\frac{\partial {\tilde{f}}_i}{\partial \tilde{t}}+{\mathbf{e}}_i·\nabla {\tilde{f}}_i=-\frac{{\tilde{f}}_i-f_i^{\left(eq\right)}\left(\tilde{\mathbf{f}}\right)}{\tau },$$

(9)

where $\tau =\lambda /\delta t$. The tilde sign will be ignored in the text below.

As in previous works (e.g., see [8,41]) the stability only on the initial conditions is investigated, so the flows in the unbouded domain are considered. The stability of the steady-state solutions of (9), associated with the steady flow regimes, is analyzed. Let $\widehat{f_i}$ correspond to the steady flow in the unbounded domain, defined by the constant values of dimensionless macrovariables: $\rho ={\rho }_0$, $\mathbf{u}={\mathbf{U}}_0$ and $e=e_0$. These solutions are derived from $f_i^{\left(eq\right)}$, as $\widehat{f_i}=f_i^{\left(eq\right)}({\rho }_0,{\mathbf{U}}_0,e_0)$.

For the stability analysis of this steady flow regime, the solutions of (9) are presented:

$$f_i(t,\mathbf{r})={\widehat{f}}_i+\delta f_i(t,\mathbf{r}),$$

(10)

where $\left|\right|\delta \mathbf{f}\left|\right|\ll 1$ is a vector of small perturbations.

Let us consider the dynamics of $\delta f_i$ at $t\to \infty$. For this purpose, (10) is substituted into (9):

$$\frac{\partial }{\partial t}\left({\widehat{f}}_i+\delta f_i\right)+{\mathbf{e}}_i·\nabla \left({\widehat{f}}_i+\delta f_i\right)=-\frac{1}{\tau }\left({\widehat{f}}_i+\delta f_i-f_i^{\left(eq\right)}\left({\widehat{\mathbf{f}}}_i+\delta {\mathbf{f}}_i\right)\right).$$

(11)

According to the proposition on the small values of perturbations, the Taylor expansion can be used:

$$f_i^{\left(eq\right)}({\widehat{\mathbf{f}}}_i+\delta {\mathbf{f}}_i)\approx f_i^{\left(eq\right)}\left(\widehat{\mathbf{f}}\right)+\sum _{s=1}^n\frac{\partial f_i^{\left(eq\right)}}{\partial f_s}\left(\widehat{\mathbf{f}}\right){\tilde{f}}_s.$$

After the substitution of this expression into (11), the following linearized system is obtained for the perturbations $\delta f_i$:

$$\frac{\partial \left(\delta f_i\right)}{\partial t}+{\mathbf{e}}_i·\nabla \delta f_i=-\frac{1}{\tau }\left(\delta f_i-\sum _{s=1}^n{A}_{is}\delta f_s\right),$$

(12)

where the components $A_{is}$ are computed as

$$A_{is}=\frac{\partial f_i^{\left(eq\right)}}{\partial f_s}\left(\widehat{\mathbf{f}}\right).$$

According to the Fourier method, the solutions of (12) can be presented as follows:

$$\delta f_i(t,\mathbf{r})=F_i\left(t\right)exp(-j\mathbf{k}\mathbf{r}),$$

(13)

where $j^2=-1$ and $\mathbf{k}$ is the wave vector, $k_x,k_y\in [-\pi ,\pi ]$. The structure of the solutions, defined by (13), corresponds to the idea of separation of variables. After the substitution of (13) into (12), the following system is obtained for $F_i\left(t\right)$:

$$\dot{\mathbf{F}}=\mathbf{G}\mathbf{F},$$

(14)

where the components of $\mathbf{G}$ are written as follows:

$$G_{is}=\left\{\begin{array}{cc}j(e_{ix}{k}_x+e_{iy}{k}_y)-\frac{1}{\tau }(1-A_{is}),\hfill & i=s,\hfill \\ \frac{1}{\tau }{A}_{is},\hfill & i\ne s.\hfill \end{array}\right.$$

For the model from [19], the dimension of $\mathbf{G}$ is $16\times 16$; for the model from [13], it is equal to $65\times 65$.

For the DDF models presented by system (6)–(7), the steady-state solutions, associated with the same steady flow regime, are computed as ${\widehat{f}}_i=f_i^{\left(eq\right)}({\rho }_0,{\mathbf{U}}_0,e_0)$, ${\widehat{h}}_i=h_i^{\left(eq\right)}({\rho }_0,{\mathbf{U}}_0,e_0)$. System (6)–(7) in dimensionless form is written as

$$\frac{\partial f_i}{\partial t}+{\mathbf{e}}_i·\nabla f_i=-\frac{1}{{\tau }_f}(f_i-f_i^{\left(eq\right)}),$$

(15)

$$\frac{\partial h_i}{\partial t}+{\mathbf{e}}_i·\nabla h_i=-\frac{1}{{\tau }_h}\left(h_i-h_i^{eq}\right)+\frac{1}{{\tau }_{hf}}\left({\mathbf{e}}_i·\mathbf{u}\right)\left(f_i-f_i^{\left(eq\right)}\right),$$

(16)

${\tau }_f=\frac{{\lambda }_f}{\delta t}$, ${\tau }_h=\frac{{\lambda }_h}{\delta t}$, $1/{\tau }_{hf}=1/{\tau }_h-1/{\tau }_f$.

The solution of (15)–(16) is presented in the same form as the previous case:

$$f_i(t,\mathbf{r})={\widehat{f}}_i+\delta f_i(t,\mathbf{r}),h_i(t,\mathbf{r})={\widehat{h}}_i+\delta h_i(t,\mathbf{r}),$$

where it is proposed, that $\parallel \delta \mathbf{f}\parallel \ll 1$, $\parallel \delta \mathbf{h}\parallel \ll 1$.

The dynamics of perturbations $\delta f_i$ and $\delta h_i$ are described by the following system of differential equations:

$$\frac{\partial \left(\delta f_i\right)}{\partial t}+{\mathbf{e}}_i·\nabla \left(\delta f_i\right)=-\frac{1}{{\tau }_f}\left(\delta f_i-\sum _{s=1}^n{B}_{is}^1\delta f_s-\sum _{s=1}^n{B}_{is}^2\delta h_s\right),$$

$$\frac{\partial \left(\delta h_i\right)}{\partial t}+{\mathbf{e}}_i·\nabla \left(\delta h_i\right)=-\frac{1}{{\tau }_h}\left(\delta h_i-\sum _{s=1}^n{B}_{is}^3\delta f_s-\sum _{s=1}^n{B}_{is}^4\delta h_s\right)+\left(\sum _{s=1}^n{D}_{is}^1\delta f_s+\sum _{s=1}^n{D}_{is}^2\delta h_s\right),$$

where

$$B_{is}^1=\frac{\partial f_i^{\left(eq\right)}}{\partial f_s}\left(\widehat{\mathbf{f}},\widehat{\mathbf{h}}\right),B_{is}^2=\frac{\partial f_i^{\left(eq\right)}}{\partial h_s}\left(\widehat{\mathbf{f}},\widehat{\mathbf{h}}\right),B_{is}^3=\frac{\partial h_i^{\left(eq\right)}}{\partial f_s}\left(\widehat{\mathbf{f}},\widehat{\mathbf{h}}\right),B_{is}^4=\frac{\partial h_i^{\left(eq\right)}}{\partial h_s}\left(\widehat{\mathbf{f}},\widehat{\mathbf{h}}\right),$$

$$D_{is}^1=\frac{\partial C_i}{\partial f_s}\left(\widehat{\mathbf{f}},\widehat{\mathbf{h}}\right),D_{is}^2=\frac{\partial C_i}{\partial h_s}\left(\widehat{\mathbf{f}},\widehat{\mathbf{h}}\right),C_i\left(\mathbf{f},\mathbf{h}\right)=\frac{1}{{\tau }_{hf}}\left({\mathbf{e}}_i·\mathbf{u}\right)\left(f_i-f_i^{\left(eq\right)}\left(\mathbf{f},\mathbf{h}\right)\right).$$

The solution of this system can be presented as follows:

$$\delta f_i(t,\mathbf{r})=F_i\left(t\right)exp(-j\mathbf{k}\mathbf{r}),\delta h_i(t,\mathbf{r})=H_i\left(t\right)exp(-j\mathbf{k}\mathbf{r}),$$

where functions $F_i\left(t\right)$ and $H_i\left(t\right)$ are the solutions of the following system:

$${\dot{F}}_i\left(t\right)=\sum _{s=1}^n{G}_{is}^1{F}_s\left(t\right)+\sum _{s=1}^n{G}_{is}^2{H}_s\left(t\right),$$

$${\dot{H}}_i\left(t\right)=\sum _{s=1}^n{G}_{is}^3{F}_s\left(t\right)+\sum _{s=1}^n{G}_{is}^4{H}_s\left(t\right),$$

where

$$G_{1is}=\left\{\begin{array}{cc}j(e_{ix}{k}_x+e_{iy}{k}_y)-\frac{1}{{\tau }_f}\left(1-B_{is}^1\right),\hfill & i=s,\hfill \\ -\frac{1}{{\tau }_f}{B}_{is}^1,\hfill & i\ne s,\hfill \end{array}\right.G_{2is}=\frac{1}{{\tau }_f}{B}_{is}^2,$$

$$G_{3is}=\left\{\begin{array}{cc}j(e_{ix}{k}_x+e_{iy}{k}_y)-\frac{1}{{\tau }_h}\left(1-B_{is}^3\right)+D_{is}^1,\hfill & i=s,\hfill \\ \frac{1}{{\tau }_h}{B}_{is}^3+D_{is}^3,\hfill & i\ne s,\hfill \end{array}\right.G_{4is}=\frac{1}{{\tau }_h}{B}_{is}^4+D_{is}^2.$$

If the extended vector $\mathbf{S}={(F_1,\dots ,F_n,H_1,\dots ,H_n)}^T$ is introduced, this system can be rewritten as

$$\dot{\mathbf{S}}=\mathbf{G}\mathbf{S},$$

(17)

where

$$\mathbf{G}=\left(\begin{array}{cc}{\mathbf{G}}_1& {\mathbf{G}}_2\\ {\mathbf{G}}_3& {\mathbf{G}}_4\end{array}\right).$$

By this approach, the problems of stability analysis of steady-state solutions are reduced to the problems of stability analysis of zero solutions of linear systems (14) and (17). These solutions are asymptotically stable if and only if the real parts of all eigenvalues of matrix $\mathbf{G}$ are negative. These eigenvalues are dependent on the following parameters: $\tau$, $U_{0x}$, $U_{0y}$, ${\rho }_0$, $e_0$, $k_x$ and $k_y$ (the models from [24,32] are analyzed at a fixed value of Prandtl number $Pr={\tau }_f/{\tau }_h$, so only the value ${\tau }_f=\tau$ can be varied if the value of $Pr$ is inputted), where only the components of $\mathbf{k}$ are considered as the internal parameters. For the simplification of analysis, the values of two dimensionless macrovariables are fixed: ${\rho }_0=e_0=1$.

The stability is analyzed for the case of the following flow regime: $U_{0x}=U$, $U_{0y}=0$, as it is realized for other LB models [41,42]. The wave vector $\mathbf{k}$ can rotate relatively to ${\mathbf{U}}_0$ in all 2D space, as it is realized in [8]. So, we restrict the parametric dependence of eigenvalues to the dependence of two internal parameters, $k_x$ and $k_y$, and two input parameters $\tau$ and U. The values of $\tau$ are chosen in small intervals near zero because the presented models in most problems are applied to the simulation of inviscid flows (when $\tau =0$). The values of U are varied in a bounded interval including zero, in order to check the maximum value, which defines the stability interval on this parameter.

The stability is investigated by the construction and analysis of the stability domains in the 2D space of the input parameters $(\tau ,U)$. These domains are constructed by the following algorithm:

(1) The grids with nodes ${\tau }_k$, $U_l$, $(k_{xp},k_{yp})$ are constructed in the intervals $(0,T]$, $(0,\mathcal{U}]$ and rectangle $[-\pi ,\pi ]\times [-\pi ,\pi ]$, respectively. The values of T and $\mathcal{U}$ are selected by the investigator.

(2) At every fixed node (${\tau }_k$,$U_l$), the eigenvalue problems for matrix $\mathbf{G}$ are solved at every node $(k_{xp},k_{yp})$ and eigenvalues ${\lambda }_m({\tau }_k,U_l,k_{xp},k_{xq})$ are obtained at every p and q. If at some node $(k_{x\overline{p}},k_{x\overline{q}})$ the following inequality is realized: $\underset{m}{min}Re\left({\lambda }_m({\tau }_k,U_l,k_{x\overline{p}},k_{x\overline{q}})\right)\ge 0$, so the solutions of the systems (14) or (17) are treated as not asymptotically stable at $\tau ={\tau }_k$ and $U=U_l$ and node $({\tau }_k,U_l)$ is not included into the stability domain in parametric space.

The eigenvalue problems are solved numerically using the QR algorithm, realized on FORTRAN90 in EISPACK subroutines [43]. We perform the computations for the monoatomic gas with 3 translational degrees of freedom per atom, which is typical for the air flow. For this physical condition, the value of $\gamma$ is equal to 5/3. The models presented in [13,19] are additionally dependent on the values of ${\eta }_0$ [13] and ${\eta }_i$, $i=\overline{1,4}$. As it is mentioned in [19,44,45], the value of this parameter can affect the numerical stability. So, we decided to investigate the influence of these parameters on the solutions. For the model from [13] the values of ${\eta }_i$ are chosen as ${\eta }_i=i{\eta }_0$. It is inspired by the same approach of choosing ${\eta }_i$ in numerical examples, presented in [13]. In Figure 2 and Figure 3, the boundaries of the stability domains for different values of ${\eta }_0$ are presented. As is demonstrated, for both models, the acceptable value of ${\eta }_0$ can be chosen as ${\eta }_0>0.25$, and the greatest values of U take place for the values $\tau <0.05$. As can be seen, the largest areas of the stability domains correspond to the model from [13]. Moreover, its area does not change after the value $\tau \approx 0.05$ (see Figure 3).

For the DDF models, the role of the external parameter at fixed $\gamma$ is played by the Prandtl number $Pr$. For the air, we consider the following values of $Pr$: 0.674 (corresponds to 273.15 K), 0.683 (corresponds to 303.15 K), 0.703 (corresponds to 343.15 K) and 0.724 (corresponds to 433.15 K). In Figure 4, the plots of the boundaries of stability domains are presented. As can be seen, the values of $Pr$ do not seriously affect the area of the domain. So, the physically correct temperature values do not seriously affect the stability of such DDF models; it depends only on values of U and relaxation time. It must be noted that the smallest areas occur for the model from [32], which is based on the D2Q17 lattice. For both DDF models, the areas of domains decrease with the increase of $\tau$. As for the multi-speed models, the largest values of U correspond to the interval $\tau <0.05$. It must be noted that the constructed models correspond to the case of viscous flows, which is described by the Navier–Stokes system, and in the expression for viscosity obtained by the application of Chapman–Enskog expansions, the dynamic viscosity is directly proportional to the value of relaxation time. So, the case of small values of $\tau$ corresponds to a small value of viscosity. So, the largest values of U correspond to the case of inviscid flows. In most of the works, the models are constructed for the simulation of this type of flow in order to simulate many interesting effects, such as shock and rarefication waves, instabilities, etc., and the best stability properties are also observed for the conditions typical for these flows.

4. Numerical Experiments

In the presented section, we try to compare the considered models, applied to the numerical solution of nonlinear 2D problems: the Riemann problems and KHI simulation. All parameter values from the stated conditions are chosen from the literature where these problems are stated.

The computations are performed using the finite-difference off-lattice scheme with the use of the technique based on the splitting of physical processes (e.g., see [46,47,48,49,50,51]). Let us consider a time grid constructed with the step $\Delta t$ and nodes $t_k$. According to the splitting conception, on one time step $[t_k,t_{k+1}]$ the two processes are considered separately. First, the collision of particles, which is described by the following system:

$$\frac{\partial {\overline{f}}_i}{\partial t}=-\frac{{\overline{f}}_i-f_i^{\left(eq\right)}\left(\overline{\mathbf{f}}\right)}{\lambda },t\in (t_k,t_{k+1}],$$

(18)

with the initial condition ${\overline{f}}_i(t_k,\mathbf{r})=f_i(t_k,\mathbf{r})$. Second, their free streaming after the collision, which is described by the following system of linear advection equations:

$$\frac{\partial f_i}{\partial t}+{\mathbf{v}}_i·\nabla f_i=0,t\in (t_k,t_{k+1}],$$

(19)

which is solved with the initial condition $f_i(t_k,\mathbf{r})={\overline{f}}_i(t_{k+1},\mathbf{r})$. By applying this technique, we can use different schemes for the implementation of both stages in order to improve the stability and accuracy of the computational procedure.

Let us consider the spatial grid with nodes $(x_j,y_s)$, constructed with step h on both variables. For the discretization of the linear system (19) we use the following parametric scheme, constructed as a combination of schemes with second-order central differences and first-order upwind differences:

$$\begin{array}{c}{f}_{i,js}^{k+1}=f_{i,js}^k-\kappa \frac{\Delta t}{2h}\left(v_{ix}\left(f_{i,j+1s}^k-f_{i,j-1s}^k\right)+v_{iy}\left(f_{i,js+1}^k-f_{i,js-1}^k\right)\right)-\\ (1-\kappa )\frac{\Delta t}{h}\left(v_{ix}{R}_{ix}+v_{iy}{R}_{iy}\right),\end{array}$$

(20)

where $f_{i,js}^k\approx f_i(t_k,x_j,y_s)$ and

$$R_{ix}=\left\{\begin{array}{c}{f}_{i,js}^k-f_{i,j-1s}^k,if{v}_{ix}\ge 0,\hfill \\ f_{i,j+1s}^k-f_{i,js}^k,if{v}_{ix}<0,\hfill \end{array}\right.R_{iy}=\left\{\begin{array}{c}{f}_{i,js}^k-f_{i,js-1}^k,if{v}_{iy}\ge 0,\hfill \\ f_{i,js+1}^k-f_{i,js}^k,if{v}_{iy}<0,\hfill \end{array}\right.$$

and $\kappa \in [0,1]$ is the dimensionless parameter. By choosing the proper value of this parameter, we can compensate the effects of numerical dispersion, typical for schemes with central differences, and the effect of numerical dissipation, typical for schemes with upwind differences [52]. The same approach is used for the equations of the DDF model.

For the implementation of (18) the explicit Euler method is used:

$$\frac{{\overline{f}}_{i,js}^{k+1}-{\overline{f}}_{i,js}^k}{\Delta t}=-\frac{1}{\lambda }\left({\overline{f}}_{i,js}^k-f_i^{eq}\left({\overline{\mathbf{f}}}_{i,js}^k\right)\right).$$

The computations are performed in the case of small value of CFL number ($v\Delta t/h=0.05$) in order to exclude any numerical instabilities.

As an example of the multi-speed model, we consider the model from [19], and as a DDF model, the model from [24] is used. The model from [13] is not considered due to the large number of equations; the model from [32] is not considered due to the small areas of the stability domain (see Figure 4).

The parameter values are selected in such a way that the solutions of kinetic equations will be stable (the parameter values correspond to the internal points of the stability domains). For considered problems, the following parameter values are used: For the model from [19], $\delta t=1$, $l=1$, $\gamma =5/3$ and ${\eta }_0=25$. For the model from [24], $\delta t$, l and $\gamma$ are the same. For all problems, $\lambda =1\times {10}^{-5}$ is used as the value of relaxation time in order to reproduce the inviscid flow because, as is mentioned above, the dynamic viscosity for this model is proportional to $\lambda$ [19]. For the DDF model from [24], we set $Pr=0.71$ and ${\lambda }_f=1\times {10}^{-5}$, ${\lambda }_h$ is computed as ${\lambda }_f/Pr$. For all problems, the value $\sigma =0.75$ is used.

The initial-boundary-value problems considered in this section are stated in the unbounded spatial domain (2D for the KHI problem and 1D for Riemann problems). In order to solve such problems numerically, the computational domain is defined as a rectangle. The effects of infinite boundaries are simulated by periodic or extrapolated boundary conditions.

4.1. Riemann Problems

Let us consider three problems, stated in the infinite interval $x\in (-\infty ,+\infty )$ with the following initial conditions for macrovariables:

(a) Problem with discontinuous velocity:

$$\rho (0,x)=1,p(0,x)=1,u_x(0,x)=\left\{\begin{array}{c}0,x\le 0,\hfill \\ 1,x>0.\hfill \end{array}\right.$$

(b) Problem with discontinuous density:

$$\rho (0,x)=\left\{\begin{array}{c}0.25,x\le 0,\hfill \\ 1,x>0.\hfill \end{array}\right.T(0,x)=1,u_x(0,x)=1,$$

where temperature T is related with the internal energy according to the used equation-of-state: $T=(\gamma -1)e$.

(c) Sod’s problem, which is widely used as a test for solvers in gas dynamics [53,54]:

$$\rho (0,x)=\left\{\begin{array}{c}1,x\le 0,\hfill \\ 0.125,x>0.\hfill \end{array}\right.,p(0,x)=\left\{\begin{array}{c}1,x\le 0,\hfill \\ 0.1,x>0.\hfill \end{array}\right.,u_x(0,x)=0.$$

This initial condition leads to the solution with a right-propagating shock, a right-propagating contact wave and a left-propagating rarefication wave.

Stated problems can be solved exactly [55], so they can be considered as powerful tools for testing numerical schemes and their implementations.

For numerical realization, we consider the finite interval on variable x: $x\in [-L,L]$, where the grid with step h and N nodes is constructed. We perform computations in the rectangular domain, where the interval on y is defined as $y\in [0,10h]$. On the upper and lower boundaries, the periodic boundary conditions are stated. At left and right boundaries, the extrapolated conditions are stated: $f_i(t,-L,y)=f_i(t,-L+h,y)$, $f_i(t,L,y)=f_i(t,L-h,y)$, and the same conditions are stated for $h_i$. The initial values of $f_i$ and $h_i$ are computed from the values of equilibrium distributions computed from the initial conditions stated for macrovariables.

The results demonstrated in Figure 5, Figure 6 and Figure 7 are obtained in the case of $L=1$ for the spatial grid with $N=500$ at the endpoints of time intervals. For problem 1, the time interval $[0,0.4]$ is considered, and for problems 2 and 3, the intervals $[0,0.3]$ and $[0,0.25]$ are used. Some minor numerical artifacts, associated with small-amplitude osclillations near the contact line, are observed according to the presence of the term with a second-order central difference in the finite-difference scheme (20). As can be seen, both models lead to very close results.

4.2. Kelvin–Helmholtz Instability

The KHI occurs in a system when two immiscible fluids flow at different tangential velocities (e.g., in opposite directions) and there are exist small perturbations near the interface (Figure 8). Such perturbations between two fluids grow in time and may evolve into turbulent mixing through nonlinear interaction. The KHI is observed in nature and industry and plays an important role in such processes as supernova dynamics, interaction of the solar wind with the Earth’s magnetosphere, cloud dynamics, etc. [56,57,58]. Furthermore, this type of instability is a typical example of non-equilibrium flow, and its non-equilibrium effect is significant near the interface, where the instabilities grow, which is why the kinetic equations are intensively used in the analysis of its evolution [59].

For the numerical simulation, the system, formed by two immiscible fluids with different initial densities and tangential velocities, is considered (see Figure 8 for schematic representation). The interface in the initial time moment is represented by the straight line, which corresponds to $x=0$.

The considered models are not adopted for multiphase simulations or simulations with interfaces between flows with different physical characteristics. So, the effect of the interface is simulated by the initial condition with the steep profile [60]. The solution corresponding to this condition evolves as a model that reproduces the dynamics of a system with immiscible fluids. So, LBM in this situation is considered as shock-capturing method.

The following initial conditions are stated [59,60]:

$$\rho \left(0,x,y\right)=\frac{{\rho }_L+{\rho }_R}{2}-\frac{{\rho }_L-{\rho }_R}{2}\tanh \left(\frac{x}{D_{\rho }}\right),$$

$$u_y\left(0,x,y\right)=\frac{u_L+u_R}{2}-\frac{u_L-u_R}{2}\tanh \left(\frac{x}{D_v}\right),$$

$$u_x(0,x,y)=u_{0}sin\left(ky\right)\exp \left(-k\left|x\right|\right),p_L(0,x,y)=p_R(0,x,y)=p_0,$$

where ${\rho }_L$, ${\rho }_R$ are constants corresponding to densities of fluids, $u_L$ and $u_R$ are corresponding tangential velocities, $D_{\rho }$ and $D_v$ are positive constants, which characterize the width of density and velocity transitional layers, $k>0$ is the frequency of the perturbations of the interface, $u_0$ characterizes its amplitude and $p_0>0$ is the initial pressure.

For the computations, we consider the rectangular domain $x\in \left[-\frac{L}{2},\frac{L}{2}\right]$, $y\in [0,L]$ in the case of $L=0.2$. The following parameter values are used: ${\rho }_L=5$, $u_L=0.5$, ${\rho }_R=2$, $u_R=-0.5$, $u_0=0.05$ and $p_0=2.5$. On all boundaries, the zeroth-order extrapolation conditions, as for the Riemann problems, are stated. The distribution functions in the initial time moment are computed as the equilibrium distributions, computed on the initial conditions, stated for macrovariables. The computations are performed for the time interval $[0,1]$. A spatial grid consists of $500\times 500$ nodes.

In Figure 9, the contour plots of $\rho$ in moment $t=0.7$ for values $k=\pi$, $D_{\rho }=0.004$ and $D_v=0.008$ are presented. As can be seen, the results obtained by both models are close to each other, and the coordinates of the vortex center are similarly predicted by both models.

In Figure 10, the plots of density profiles corresponding to the fixed value $y=0.5L$ are presented for the values $D_{\rho }=0.002$ and $D_v=0.008$ in the cases of various values of k, which characterize the frequency of interface perturbation. As it can be seen, the solutions obtained by both models lead to similar profiles with weak discrepancies in values of $\rho$ (not larger than 0.2). The profiles, obtained for various values of $D_v$, are presented in Figure 11 for the case of $k=10\pi$, $D_{\rho }=0.002$. As can be seen, both models lead to similar profiles (discrepancies on values of $\rho$ are not larger than 0.25 ($\approx 5$ %)). It should be noted that the complete investigation of parameters influence on the solution of this problem is performed and discussed in detail in [60]. In the presented paper, we focus attention only on the comparison of numerical solutions obtained by different models based on kinetic equations.

For the estimation of $\rho$, obtained by kinetic models in a whole spatial domain, the averaged value can be used:

$$\overline{\rho }\left(t\right)=\frac{1}{L^2}\underset{-\frac{L}{2}}{\overset{\frac{L}{2}}{\int }}\underset{0}{\overset{L}{\int }}\rho (t,x,y)dydx.$$

The computation of this characteristic is performed by the following quadrature formula:

$$\overline{\rho }\left(t_s\right)\approx \frac{h^2}{L^2}\sum _{i,j}{\rho }_{i,j}^s,$$

where ${\rho }_{i,j}^s\approx \rho (t_s,x_i,y_j)$ are computed from the values of distribution functions. In Figure 12, the plots of $\overline{\rho }$, obtained for $D_v=0.008$ and $D_{\rho }=0.004$ for different values of k, are presented. In Figure 13, the plots for the case of $k=10\pi$, $D_{\rho }=0.004$ and different values of $D_v$ are presented. As can be seen, the values of the averaged density obtained by both models are close to each other. The highest deviations of values of $\overline{\rho }$ are not larger than 0.25.

5. Conclusions

In the presented paper, the analysis of discrete velocity kinetic models used in the LBM simulation of compressible flows in the LBM for cases of varied specific heat ratios in the equation-of-state is performed. The popular multi-speed and DDF models, presented in the literature, are compared. In the theoretical part, the stability analysis of steady-state solutions associated with steady flow regimes is performed. The stability analysis problem is reduced to the analysis of zero solutions of linear systems of ordinary differential equations. These systems are obtained after the linearization of kinetic equations near the steady-state solution. The analysis is based on the construction of stability domains in parametric space. The results obtained for different models are compared. In the practical part, the selected models are applied to the simulation of compressible flows. The Riemann problems and the simulation of KHI are considered. For the numerical simulations, the parametric finite-difference scheme is used. As is demonstrated, different kinetic models lead to close numerical solutions for both types of problems.

As the main results, the following conclusions can be formulated:

1. The procedure for the linear stability analysis of compressible discrete kinetic models is proposed, and stability domains are constructed. The approach is based on the considering of steady-state solutions associated with steady flow regimes. According to the proposed procedure, we reduce the analysis of steady-state solutions of nonlinear kinetic equations to the analysis of zero solutions of linear systems of ordinary differential equations. According to the stability criterion, these solutions are asymptotically stable if and only if the real parts of the eigenvalues of the matrices of the corresponding systems have negative real parts. With the use of this criterion, we can construct the stability domains in the parametric space. Only the stability of solutions of kinetic equations, which formed the model, is analyzed. In practice, the numerical instabilities associated with the used computational schemes also take place.

The proposed approach to analysis can be used as an effective tool for the theoretical comparison of different kinetic models used in the LBM in other fields, for example, for multiphase and magnetohydrodynamic simulations. It allows for the possibility to define the range of free parameters presented in the models.

2. The effects of the input parameters on the boundaries of the stability domains are analyzed. For the models from [13,19], the areas of the domains can be enlarged by the increase of ${\eta }_0$ (see Figure 2, Figure 3 and Figure 4). For the model from [19], the upper values of U, which provide stability, decrease with the increase of $\tau$ (Figure 2). However, for the model from [13], this value for $\tau >0.05$ can be approximated by the constant 2.4 (see Figure 3). For the DDF models, the value of $Pr$ very slightly affects the boundaries of the stability domains (see Figure 4).

For both types of models, the largest value of U occurs for the small values of $\tau$ ($\tau <0.05$). This means that the models provide the most stable results if the semi-inviscid flows are simulated. Despite the fact that the considered models are developed in order to mathematically reproduce Navier–Stokes equations, which describe viscous flows, in most works they are used for the simulation of inviscid flows, which are described by the Euler equations. The expression for viscosity is obtained by applying the Chapman–Enskog asymptotic expansion method. As it is demonstrated, the dynamic viscosity $\mu$ is proportional to $\lambda$ (that is why the inviscid flows are simulated by the use of a small relaxation time). For example, for the model from [13] the following expression takes place:

$$\mu =\frac{2}{D}\rho e\lambda .$$

The case of the small value of $\tau$ corresponds to the small value of $\lambda$. So, in the presented work, it is demonstrated that these models are better suited to be used in practice for the simulation of inviscid flows because, in this case, they lead to larger stability intervals on values of U.

3. As it is obtained after numerical experiments for the case of the air flow, the smallest areas of the domains are observed for the model from [32] and for the model from [19] at small values of ${\eta }_0$. The best stability properties (largest areas) occur for the model from [13]. However, this model consists of 65 equations, so it is not optimal for practical computations.

4. As is demonstrated in the numerical experiments on the simulation of compressible flows, both types of models (multi-speed and DDF) provide close values for macrovariables. So, in practice, any of the constructed models can be used. However, their parameters should be defined based on the stability conditions obtained after the analysis of stability domains.

5. The presented study can be considered as the first stage of the analysis of models used in the LBM for compressible flow simulations. The presented approach to linear stability analysis can be used for different models of different physical phenomena. In the LBM over the last decades, multiparametric methods, which improve stability, have been proposed; for example, MRT methods [38,61,62] and cascaded models [63,64,65]. The stability of such schemes and their variations will be analyzed in future studies.