Discrete Velocity Boltzmann Model for Quasi-Incompressible Hydrodynamics

: In this paper, we consider the development of the two-dimensional discrete velocity Boltzmann model on a nine-velocity lattice. Compared to the conventional lattice Boltzmann approach for the present model, the collision rules for the interacting particles are formulated explicitly. The collisions are tailored in such a way that mass, momentum and energy are conserved and the H -theorem is fulﬁlled. By applying the Chapman–Enskog expansion, we show that the model recovers quasi-incompressible hydrodynamic equations for small Mach number limit and we derive the closed expression for the viscosity, depending on the collision cross-sections. In addition, the numerical implementation of the model with the on-lattice streaming and local collision step is proposed. As test problems, the shear wave decay and Taylor–Green vortex are considered, and a comparison of the numerical simulations with the analytical solutions is presented.


Introduction
In the kinetic theory, the distribution function of a rarefied gaseous system is governed by the Boltzmann equation or its models [1]. In the applications, the discretization of these equations in the velocity (and physical) space is usually performed. One of the most popular discretization approaches is the Lattice-Boltzmann (LB) method [2][3][4][5] which was initially developed as an alternative to the continuum fluid methods like Navier-Stokes equations [6]; furthermore, the method has been extended to the rarefied flows modeling [7][8][9][10][11][12][13][14][15][16][17][18][19]. The conventional LB model has the following form where f i (t, x) is the distribution function related to the particles with the velocity c i , i = 1 . . . N, τ is the relaxation time, f eq i is the local equilibrium, N is the number of the discrete velocities, d dt = ∂ ∂t + c i ∂ ∂r , r is the spatial variable. In this approach, the collisions between the particles are described in a phenomenological way, i.e., it is postulated that, due to the collisions, the distribution function tends to the local equilibrium state at a rate proportional to f eq i − f i . For LB models, the local equilibrium is usually taken as a finite-order polynomial on the bulk velocity, and the conservation laws for mass and momentum are satisfied by construction. On the other hand, for this form of the local equilibrium, the H-theorem does not exist [20][21][22]. To overcome this issue, models with non-polynomial equilibria have been proposed [23][24][25][26].
Another possible discretization technique is the discrete velocity (DV) Boltzmann method [27][28][29][30], the general DV Boltzmann model reads as where A ij kl = A kl ij ≥ 0 are the transition probabilities.
Compared to the LB method, the DV models have some attractive properties. Similarly to the Boltzmann equation, the binary collisions are described explicitly. Moreover, by construction, the H-theorem is valid for these models [28], i.e., Moreover, the local equilibrium for DV kinetic Boltzmann models can be obtained as an exponential function of the macroscopic variables. The DV Boltzmann approach attracted the attention of many researchers several decades ago, but, at present, is significantly less popular than the LB method. For instance, the well-known four velocity Broadwell equation in two dimensions has been investigated thoroughly [31][32][33][34][35][36], this model has correct collision invariants, but its discrete velocity set is too small and lacks isotropy [37]; therefore, the correct description of the hydrodynamics is impossible in the framework of this model. In addition, another subtle feature should be mentioned: for the discrete velocity models, the molecular chaos hypothesis can be violated, i.e., the particles can be correlated before the collision [38]. This is undesirable, but the influence of this effect on the flow properties in applications is not clear. Furthermore, one should construct the DV Boltzmann models in such a way that the only conserved variables are mass, momentum and energy. The equilibrium state is obtained as minimum of the H-function under the constraint that these variables are not changed by collisions. The presence of other conserved quantities (spurious invariants) changes the form of the local equilibrium state, this, in turn, leads to a distortion of the hydrodynamic equations. The construction of DV Boltzmann models without excessive invariants is a non-trivial procedure [39][40][41][42].
In this paper we consider a DV Boltzmann model on a nine velocity, two-dimensional lattice. As a starting point, we consider the local equilibrium for the general DV Boltzmann model and its expansion at the vicinity of the absolute Maxwell distribution. Next, the Chapman-Enskog expansion for the DV Boltzmann model is performed in order to derive the hydrodynamic equations. In addition, we show that the model does not have invariants without physical meaning. The considered model has four different possible transition probabilities. In terms of the LB theory, this model can be considered as a scheme with multiple relaxation times. For viscosity, we obtain a closed expression depending on the values of the transition probabilities. If the viscosity is fixed, we obtain a constraint on the transition probabilities, but three of them can be chosen as free parameters; for instance, they can be adjusted to improve stability properties. As benchmark problems, we consider the shear wave decay and Taylor-Green vortex. The numerical experiments show excellent agreement between the numerical simulation results and analytical solutions.

Equilibrium for DV Boltzmann Kinetic Model and the Euler Equations
The local equilibrium of the model (1) is obtained as a minimum of the H functional with the constraints corresponding to the conservation laws; it has the following form (Formula (5) in [43]) where the coefficients a i , b i , d i depend on the density, flow velocityand temperature ρ, u, θ and "·" defines scalar product. In this paragraph, we consider the particle's dynamics in D spatial dimensions. We assume that the local equilibrium is close to the absolute equilibrium with the density ρ 0 = 1 flow velocity u 0 = 0 and the temperature θ 0 ; then, one can write down the absolute equilibrium denoted as w i in the form where a 0 = a 0 (ρ 0 , θ 0 ), d 0 = d 0 (ρ 0 , θ 0 ), we also term w i as lattice weights. The conservation laws for mass, momentum, energy yield the following equations for the lattice velocities and the lattice weights note that c i c i is a tensor with elements c i,α c i,β , α, β = 1 . . . D, and that D is the number of spatial dimensions. For the coefficients, we have the expression where ∆a, ∆b, ∆d are small quantities. Similarly to the previous studies [43] we expand the expression (2) on ∆a, ∆b, ∆d, and one has where o(∆ 2 ) stands for o(∆a 2 ), o(∆b 2 ), o(∆d 2 ), also the operator ":" is tensor convolution. Next, we assume that where ∆ρ, ∆u, ∆θ are small. Taking (4) into account, for the first local equilibrium moments, we derive the following equations ∑ i f eq i = 1 + ∆ρ = 1 + ∆a + Dθ 0 ∆d + and we omitted third-order terms and used the definitions We assume that ∆ρ, ∆u, ∆θ are of the same order of smallness, which we define as O(∆). Then, we seek solutions to the Equations (5)-(7) in the form ∆a = ∆a lin + ∆a nonl , ∆b = ∆b lin + ∆b nonl , ∆d = ∆d lin + ∆d nonl , where the terms ∆a lin , ∆b lin , ∆d lin are solutions to the linearized Equations (5)-(7) of order O(∆) and ∆a nonl , ∆b nonl , ∆d nonl are nonlinear corrections of order O(∆ 2 ). First, from the linearized Equations (5)-(7), one has ∆ρ = ∆a lin + Dθ 0 ∆d lin , ∆u = θ 0 ∆b lin , Dθ 0 ∆ρ + D∆θ = Dθ 0 ∆a lin + m 4 ∆d lin , these equations have the solutions Next, we find the nonlinear corrections ∆a nonl , ∆b nonl , ∆d nonl from the Equations (5)- (7). It would be convenient to start with the Equation (6), which can be rewritten as ∆ρ∆u = θ 0 ∆a lin ∆b lin + D −1 m 4 ∆d lin ∆b lin + θ 0 ∆b nonl , from the last equation, we immediately obtain Consideration of the Equations (5) and (7) yields by applying (9) we get the solutions The combination of (9) and (10)- (12) leads to the following expression for f eq i Proposition 1. The DV local equilibrium f eq i in the form (2) can be expressed as where and ∆ρ, ∆u ∆θ are small density, flow velocity and temperature variations of order O(∆); moreover, the moments m 4 and m 6 are defined by (8), and the absolute equilibrium w i satisfies the conditions where θ 0 is the reference temperature.
By applying (13)- (19), the pressure tensor P with the components P αβ , α, β = 1 . . . D at the local equilibrium can be derived where ρθ = (1 + ∆ρ)(θ 0 + ∆θ) and Now, let us assume that R is isotropic tensor; in such a case, its components can be written in the form (Formula (69) in [37]) one can see that the tensor P equalling ρθδ Compared to P for the local Maxwell distribution, in the DV approach, the error O(∆ 3 ) is observed; therefore, we can conclude that the hydrodynamics (mass and momentum equations) at the Euler level of accuracy is recovered with the errors of order O(∆ 3 ) if the conditions (3), (20) and (21) are satisfied.
Finally, we consider the heat flow q at the level of the Euler equations one can see that the terms proportional ∆u, ∆ρ∆u are eliminated if m 4 satisfies (21), in addition, the second term can be removed if m 6 = D(D + 2)(D + 4)θ 3 0 or we are restricted by the isothermal flows ∆θ = 0; in such a case, the heat flow (which equals zero for the Euler equations) is only of order O(∆ 3 ). In the present study, we assume that the temperature variations are negligible, ∆θ = 0.

Navier-Stokes Equations
In order to obtain the Navier-Stokes equations, one needs to find the corrections to the pressure tensor corresponding to the viscous terms. This can be performed by applying the Chapman-Enskog expansion for DV Boltzmann model [29]. Then, following the previous results [29], we assume that the solution to (1) can be expressed in the form and Kn is the Knudsen number.
At the limit of small Mach numbers, the equations for f (1) i read as (Equation (19) in [43]) one can see from (22) (22) can be obtained as (Formula (22) in [43]) where Q i is a second-order tensor whose exact form we will discuss further, and a i , b i are numerical coefficients.

Spurious Invariants
For a collision, in which the particles with the velocities c i , c j turn into the particles with the velocities c k , c l , we introduce the following reaction vector [40,41,44] where the entries denoted by dots equal zero. Assume that we have p linearly independent reaction vectors e s , s = 1 . . . p . We denote a matrix consisting of all reaction row vectors e s as the collision matrix C = (e 1 ; e 2 ; . . . ; e p ) ∈ R p × R N .
Note that the collision invariants ϕ(c 1 , . . . c N ) ∈ R N are defined by the relation [28] this condition can also be rewritten in the following form [44,45] ϕ · e s = 0, s = 1 . . . p, i.e., the linear subspace spanned by the invariants is orthogonal to the subspace spanned by the reaction vectors. The condition (24) can be applied for the detection of spurious invariants: Proposition 2. Assume that, for some DV Boltzmann models, the number of linearly independent physical collision invariants equals q, then additional invariants do not exist if [40,41,44,45] where N is the number of the discrete velocities.

Nine Velocity DV Boltzmann Model for D2Q9 Lattice
We consider the DV Boltzmann model on a nine-velocity lattice (Figure 1). This lattice is popular in LB theory [3], since the corresponding LB model recovers hydrodynamics at small Mach numbers limit and, in addition, its numerical implementation is very simple.  It is well-known that these lattice velocities and weights satisfy the conditions (3), (20) and (21); therefore, if it is possible to construct the collisions in such a way that the mass and momentum are conserved then the Euler equations are satisfied. We mention that the lattices and collision rules for DV Boltzmann models, which can potentially recover the hydrodynamics, have been considered previously [46,47]-for instance, the model with single-relaxation time describing Navier-Stokes equations has been proposed [46,47]. In here we consider only the collisions for the nine-bit lattice in a more detailed way; the considered model is of the multiple-relaxation-time type: a. Broadwell type collision is the reaction between the particles 1 and −1, which turn into the particles 2 and −2 ( Figure 1); schematically, we can denote this reaction as (1, −1) −→ (2, −2). The contribution of this collision to right side of (1) denoted as J 0 is as follows b. the collisions linking all three different energy states, they define transitions between the particle's states with different kinetic energies, and evidently can not be excluded. We have four different reactions (1, 2) −→ (0, 3 c. Broadwell type collision between the particles with the velocity magnitudes √ 2c is defined by the reaction (3, −3) −→ (4, −4), the contributions to the collision kernel are The collisions (25)-(28) conserve mass, momentum and energy; the corresponding D2Q9 DV Boltzmann model reads as where α, β, λ, γ in (29)-(37) are positive transition probabilities. Now, we can consider the analogs of the Navier-Stokes equations for the model (29)-(37).

Proposition 3. The Equations (29)-(37) lead to Navier-Stokes equations for nearly incompressible flows with errors of order O(∆ 3 ) if
where ∂ ∂r = ( ∂ ∂x , ∂ ∂y ). According (23) we can try add the terms proportional div(∆u), but they equal zero for the incompressible limit; then, we seek the solution in the form where the coefficients a i are equal for the indexes i corresponding to the discrete velocities c i with the same kinetic energy. The substitution of (41) into (29)-(37) leads to three algebraic equations for the coefficients a i , (29)-(32) yield the first equation from which we obtain a 1 = − 3 4α , (33)-(36) yield the second equation The third equation, which can be obtained from (37) where σ, η equal x or y. For instance, we require 4α = γ + 4β + 4λ, then by applying div(∆u) = 0 we finally obtain For the model (29)- (37), there are ten collisions. If we consider all reaction vectors and the corresponding collision matrix, one can convince that rank(C) = 5, the number of the discrete velocities N = 9. This means that we do not have any collision invariants except mass, momentum, energy (Proposition 2). We can exclude up to five reactions from the model; for instance, we can keep only the Broadwell collisions (type a.) and the collisions of type b., i.e., we set γ = λ = 0. On the other side, the numerical simulations show that the addition of the collisions from the group c. or d. enhances the stability properties.
Finally, we emphasize that, for the model (29)-(37), all the collisions conserve energy (elastic). Generally speaking, this is not necessary because we are focused on the correct reproduction of the mass and momentum equations. For instance, it is possible to construct the model of DV Boltzmann type in one spatial dimension with inelastic collisions [26] (quasi-chemical model with three discrete velocities) which leads to the correct Navier-Stokes equation at small Mach limit.

Numerical Implementation and Test Problems
The model is implemented similarly to the conventional LB D2Q9 model [3]. Firstly, we perform the collision step, then post-collision distribution functions are streamed at appropriate directions. It is well-known from the LB theory that the discretization of space-time affects the viscosity. The DV Boltzmann model discretized in a similar form as LB model reads as by applying the Taylor expansion this equation can be rewritten as therefore, we can conclude that the scheme (42) led to the hydrodynamic equations in which the contributions from − 1 The additional terms for the Navier-Stokes equations can be obtained with the application of the Chapman-Enskog expansion [3]. Note that the terms O(δt 2 ) do not affect the Navier-Stokes equations, since they contain third order derivatives, which, in the Chapman-Enskog multiple-scale expansion, enter the equations for the moments at the Burnett level. For the Navier-Stokes equation, the additional viscosity terms result from − δt dr d dt f eq i , remembering that d dt f eq i can be expressed by (40), we eventually obtain δt 6 ( ∂ 2 ∂x 2 + ∂ 2 ∂y 2 )u. Then, for the DV D2Q9 Boltzmann model in the form (42), the viscosity is given by In the simulations, the parameters are taken as follows i.e., we have six different collisions.
To validate the second-order convergence of the presented scheme, we estimate the simulation error defined as where z denotes the spatial variable, u m , u bench are the modeled variable (velocity) and the benchmark solution, respectively. The convergence rate is evaluated by fitting the values of log(error) for the various log(h) = log( 1 N ) (N is the number of the lattice nodes, h is proportional to the lattice spacing) using the linear regression, the second-order convergence is achieved if the regression slope coefficient is close to 2.
Compared to LB D2Q9 model, the scheme (29)-(37) differs only in the collision term and the expression for the viscosity. This means that the computation time for (29)- (37) implemented in the form (42) is approximately the same as for LB D2Q9 model.

Shear Wave Decay
We consider the dynamics in terms of the time of the sinusoidal velocity wave in a square domain. The initial flow velocity in x direction is dependent on y coordinate and is given by where L is the length of the domain equals N lattice nodes and U 0 = 0.01. The periodic boundary conditions are applied for the present problem. This problem has the following analytical solution u x (x, y, t) = U 0 sin(ky)e −νk 2 t .
In the present case, we consider ν = 0.001 and N = 101, the time step δt = 1. We compare the analytical solutions with the velocity profiles obtained by the application of the model (29)-(37) (implemented in the form (42)). The peak velocity time history and the velocity profiles for the different moments of time are plotted, Figure 2. One can see that the simulation results are very similar to the analytical profiles. It is worth mentioning that it is possible to shorten the model and take γ = λ = 0, in this case α = β, and we have only five different collisions. The numerical experiments show that this model becomes unstable for ν < 0.1, while the setting (43) allows to model the flow with small viscosity and no instabilities are observed.

Taylor-Green Vortex
Similarly to the previous problem, we consider a square domain, and the initial velocity field is given by the formula u x (x, y, t = 0) = −U 0 cos(kx) sin(ky), u y (x, y, t = 0) = U 0 sin(kx) cos(ky), where the size of the domain is L × L (or N × N in lattice units, where N is the number of the lattice nodes) and k = 2π L . The periodic boundary conditions are applied. For the present problem we set U 0 = 0.01, ν = 0.001, N = 51, the time step δt = 1. The analytical solution to the problem is as follows u x (x, y, t) = −U 0 cos(kx) sin(ky)e −2νk 2 t , u y (x, y, t) = U 0 sin(kx) cos(ky)e −2νk 2 t , one can see that the initial structure of the velocity field persists in time, and only uniform decay of the velocity amplitudes is observed. The numerical simulations for the model (29)-(37) (implemented in the form (42)) show that the form of the velocity field does not change. We also present the behavior of the velocity u x (x, y = L/2, t) over time, obtained analytically and numerically for three different moments of time; obviously, both approaches give very similar profiles ( Figure 3).
Finally, we consider the convergence rates of the numerical simulation results to the benchmark solutions. This can be performed by considering the logarithms of the simulation errors (44) for the different values of log(h) = log(1/N). In the present case, we take N = 25, 49, 73, 101. In Figure 4, the logarithms of the errors of the velocities are presented for DV and the conventional LB D2Q9 models; the results are very similar for both models. One can see that the estimated slope values are close to 2; this indicates that the proposed scheme is accurate in the second-order.   . Convergence rates for the shear wave decay and Taylor-Green vortex problems are shown. The results are obtained by applying DV and the conventional LB D2Q9 models. In the (first slide) (shear wave decay), the logarithms of the errors (44) for the velocity u x (y, t) computed at the moment of time t = 1/(νk 2 ) are presented; in the (second slide) (Taylor-Green vortex), the logarithms of the errors of the velocity u x (x, y = L/2, t) computed at the moment of time t = 1/(2νk 2 ) are presented, where the variable h is proportional to the lattice spacing. The slope estimates are obtained by fitting the values of log(error) using the linear regression.

Results and Discussion
In this paper, we have considered the DV Boltzmann model applicable to the modeling of viscous quasi-incompressible flows at a small Mach number limit. The presented model has the same discrete velocity structure and absolute equilibrium as LB D2Q9, but the collision rules for the particles are postulated exactly. There are four types of collision and ten possible different collisions; the unique transition probability corresponds to all possible reactions in the group. Moreover, these collisions conserve only mass, momentum and energy (spurious invariants do not exist). In terms of LB theory, this model can be considered as a scheme with multiple relaxation times. Note that the H-theorem is valid for the model by construction (at least for the continuous space-time variables).
We have demonstrated that DV Boltzmann equations can be a viable tool in modeling of hydrodynamic flows. The shear wave decay and Taylor-Green vortex have been considered as benchmark problems. The comparison of the simulation results with the analytical solutions has shown good accuracy.
One of the most intriguing problems is the evaluation of the stability properties of the presented DV Boltzmann system and the optimal choice of transition probabilities. One can expect that the DV Boltzmann model for D2Q9 lattice has a better stability than the conventional LB D2Q9 model, since the H-theorem is satisfied. In order to elucidate this issue, one can consider additional problems like Sod shock tube, double shear layer and lid-driven cavity. These problems are left for future study.

Data Availability Statement:
The simulation code that supports the findings of this study is available from the author upon reasonable request.