1. Introduction
Studying stability of lattice Boltzmann schemes is a non-trivial problem. Classically for this purpose, the scheme is linearized around a constant state and a Fourier analysis is performed. We refer to the work of Lallemand and Luo [9] for the D2Q9 scheme applied to hydrodynamics. Note also the work of Ginzburg et al. [8] extending the Fourier analysis to a wide variety of different two and three dimensional lattice Boltzmann schemes. Instabilities and their interpretattion in terms of bulk viscosity has been proposed by Dellar [2]. But no mathematical analysis has been performed.
A new way of improving stability is proposed by Geier [7], who proposed a new generalized lattice Boltzmann scheme with the approach of relative velocities and utilized it for hydrodynamics applications [4]. The tentative of analysis of this method for a two-dimensional scalar linear equation has been also proposed [5].
Even if it is a difficult task, it is well known that Fourier analysis is not the best method for analyzing nonlinear hyperbolic equations. Total variation diminishing schemes, developed for suppressing oscillations in higher order CFD algorithms, provide an alternative nonlinear stability analysis tool for analysing the schemes for nonlinear wave propagation. The convergence of such schemes is well established [10]. The underlying stability notion concerns the maximum principle. This notion can be extended to nonlinear cases and a first attempt has been proposed in [1] for lattice Boltzmann schemes for the D1Q2 scheme for scalar nonlinear hyperbolic equations.
In this contribution, we propose to investigate the stability in the maximum sense, of a linear mono-dimensional lattice Boltzmann scheme with three velocities. More precisely, we look to a positivity constraint for a particle distribution function in the context of a relative velocities.
In Section 2, we describe the scheme and the underlying advection model. More precisely, the local relaxation step is written as a linear operator on the particle distribution functions. If all the coefficients of the underlying matrix are nonnegative, the non-negativity of the distribution is maintained during this step. Because the transport step is just a change of locus, the non-negativity is maintained for the whole time step of the scheme. The question is then to find appropriate conditions to handle this property. In Section 3, a necessary and sufficient condition is derived on the parameters in order to ensure that the scheme has the stability property. In Section 4, we completely describe the classical case where the relative velocity is reduced to zero. In Section 5, the general case is presented. With an analytical study for necessary conditions and a numerical one for a complete description of the stability zones. In Section 6, numerical experiments show the correlation of the positivity constraint for a particle distribution and the presence of oscillations for discontinuous profiles.
2. Description of the framework
2.1. Description of the scheme
In this contribution, we investigate a mono-dimensional three velocities linear lattice Boltzmann scheme with relative velocity [4]. Denoting the spatial step, the time step, and the lattive velocity, this scheme can be described in a generalized D. d’Humière’s framework [3]:
-
[(1)]
-
the 3 velocities , , and ;
-
the 3 associated distributions , , and ;
-
the equilibrium value of the 3 moments
where and are given scalars (without loose of generality, we assume that );
-
the 2 relaxation parameters and such that the relaxation phase reads
In this formalism, the moments are defined as polynomial functions of the discrete velocities and the discrete distribution functions. Indeed, introducting , , and , the three moments read
The equilibrium values are chosen such that the equilibrium distributions do not depend on the relative velocity . Indeed, we have:
Note that this scheme can be used (see e.g. [6]) to simulate a scalar transport equation with constant velocity given by
Denoting , , the discret spatial points and , , the discret time points, one time step of the scheme reads
2.2. Matrix notation for the relaxation step
In this contribution, we are concerned with the non-negativity of the particle distribution functions. This property can be viewed as refering to a weak maximum principle for linear schemes. Indeed, it is always possible, by adding a constant, to assume that all the particle distribution functions are initialy nonnegative. If the scheme ensures that this property of non-negativity remains as time marches, each particle distribution function is then bounded, as their total sum is conserved. Moreover, as the transport step consists simply in exchanging the position of the particle distribution functions, we focus on the relaxation step.
We use a matrix notation for the relaxation step as it can be read as a multiplication by a matrix. As this step is local in space, we omit the dependancy on time and space. We define the vector of the distribution functions
One relaxation step then reads
where the matrix is defined by
with
The coefficients of the matrix are obtained by the relations , and those of the matrix by the change of basis formula: is the coefficient of the -element in the definition of according to
The matrix is then the change of basis that transforms the vector into the vector :
The matrix can then be viewed as the change of basis matrix from the classical moments without relative velocity toward the moments with relative velocity.
2.3. Remark on the choice of the moments
Note that the last moment that is chosen in this contribution is not the energy but a moment that is orthogonal to the two first ones, and . In this section, we show that all the results of the contribution would be identical by choosing the last moment as the energy: the relaxation matrix would still be the same.
We consider two schemes with two different choices of polynomials: the moments of the first scheme are defined by while the moments of the second scheme by . The first moment is the same in both schemes to be able to simulate the same transport equation. We then have . We define the change of basis matrix associated to the tranformation into :
The first line of is then .
We assume that the equilibrium values of the distribution functions are the same, that is , and that the relaxation parameters are the same, that is . Then, we have for all iff
Proof.
First, we immediately obtain the following relations by identifying the coefficients of and of :
We deduce that
Moreover, as the second and the third columns of are zero (the equilibrium values depend only on the first moment ), and as the first line of is then , we have
We have
Then
As the matrices , , and are invertible, the condition is equivalent to . Denoting
a straightforward calculation yields
Then the property for all values of and is equivalent to , that ends the proof. ∎
2.4. Positivity of the matrix
The velocity being fixed, we propose to give a full description of the sets
Indeed, the non-negativity of the matrix imposes that all the distributions , , remain non-negative if they are so at the initial time. These sets are first described by a set of nine inequalities that can be joined into just one. Numerical illustrations are then given to visualize it in the characteristic cases including SRT, MRT, and relative velocity scheme.
3. Positivity of the iterative matrix
The nine inequalities obtained from the matrix can be combined neatly into one formula.
The inequalities are
(1) |
We prove now that the previous nine inequalities can be written in a much more lucid way.
We introduce the reduced parameters and according to
(2) |
Then the nine previous inequalities displayed in (1) are equivalent to
(3) |
Proof.
Consider first the two inequalities associated with and :
They can be synthetized in the following form
and we can write this relation as
(4) |
The last inequality can be written as and this inequality is equivalent to . In other terms,
(8) |
4. The particular case
In this section, we suppose that the relative velocity is reduced to zero. Then the necessary and sufficient conditions (3) for the stability can be written as
(9) |
To fix the ideas, we suppose that the advection velocity is non-negative:
(10) |
The case follows directly. When , the reduced stability conditions
(11) |
are equivalent to the following conditions for the relaxation parameters
(12) |
joined with a natural Courant type condition for explicit schemes on the advection velocity
(13) |
Of course, the conditions (9) have still to be imposed for the equilibrium parameter when the pair is given. In particular,
(14) |
and
(15) |

Proof.
We first observe that . Then . Secondly, we have and because both and are positive, we have . We have also and (13) is established.
Moreover, and implies . Due to the positivity of the parameter , we deduce from (9) a new set of inequalities: . Then . Moreover, because .
Conversely, if the relations (12) and (13) are satisfied, we have , and . Then . Moreover, , because and . Thus .
Finally the inequalities (11) are established and the proposition is proven. ∎
5. The general case
We analyse the general case of non-zero in this section, with suitable illustrations of the stability region for various ranges of the parameters.
5.1. Necessary conditions for stability
In this subsection, we prove the following proposition.
We suppose that the necessary and sufficient stability conditions (3) are satisfied:
with the notations (2):
We suppose also that the advection velicity is positive. Then if the scheme is stable, the point satisfies the following inequalities
(16) |
With these necessary stability conditions, the parameter has been eliminated. We have also, necessarily
(17) |
and
(18) |
Proof.
We start from the inequalities (3). Then we have and . With the same type of inequalities, and .
If we remark that , we deduce that and the inequality (18) is proven.
We have the triangular inequality . Then from the genaral stability conditions (3), we deduce and . In a similar way, , from (3) and finally . We put together the two inequalities and we have .
In consequence, we have and the third point is proven. Since we made the choice of then . The first inequality of the two first points are established. From we have and the relation (17) is true. Moreover, and the last inequality of (16) is true.
Consider now the inequalities . We deduce and . Due to the positvity of the absolute values, we have also and . A part of the fifth inequality of (16) is proven.
We illustrate the zones of necessary stability in the Figure 2 for five particular velocities: .
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
5.2. Numerical study of necessary and sufficient conditions for stability
We now illustrate the necessary and sufficient stability regions for , , , , , and for various ranges of in the Figure 3.
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
For each value of the velocity , the necessary and sufficient stability regions given by the relations (3) are displayed for several values of the relative velocity: . In dotted line, the necessary stability region of the Proposition 5.1 is added for remind. The particular value is enhanced by filling the region in gray.
Some analysis can be drawn from these figures:
-
the stability region changes with the relative velocity;
-
the maximal value of the first relaxation parameter is obtained (not only) for ;
-
the stability region is not clearly more interessant (larger or including greater value of ) for ;
-
the point is always in the stability region.
To conclude this section, this notion of stability allows a large set of values for the relaxation parameters. If the scheme is used to simulate the hyperbolic advection equation without diffusive term, we can try to minimize the numerical diffusion while maintaining this stability property. This task is complicated and out of the scope of the paper as the second-order term reads as a non-linear formula that links all the parameters , , and .
6. Numerical illustrations
In this section, we illustrate the stability property with numerical simulations involving the D1Q3 model with and without relative velocity used to simulate the linear advection equation. The stability is demonstrated for the parameters chosen according to the analysis presented in the previous sections. Oscillations are seen whenever the parameters go beyond the stability limits presented, as highlighted in the results.
The parameters chosen for the simulations are the following:
left (stable) | 0.25 | 0.0 | 1.6 | 1.3 | 0.3076923076923076 |
0.25 | 0.25 | 1.6 | 1.3 | -0.17548076923076938 | |
right (unstable) | 0.25 | 0.0 | 1.9 | 1.4 | 0.14285714285714302 |
0.25 | 0.25 | 1.9 | 1.4 | -0.10491071428571441 |
For the left figures, the parameters are chosen in order to satisfy the stability property. We observe numerically also a maximum principle, even if it is not formally the stability notion that we investigate. For the right figures, the parameters are chosen in order not to break the stability property.
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
If the profile is smooth, we have observed that no numerical oscillations occur even if the non-negativity property of the matrix is not satisfied. If the profile is just continuous, small negative values of the macroscopic quantity are observed when our non-negativity property is not satisfied. Last but not least, classical oscillations are visible for discontinuous profiles if our non-negativity property is not satisfied. These oscillations are eliminated when the non-negativity property of the matrix is realized.
7. Conclusion
In this contribution, we have investigated a stability property for a classical mono-dimensional linear three velocities lattice Boltzmann scheme with relative velocity. This property ensures that the particle distribution functions remain nonnegative if they are at initial time. We then give a necessary and sufficient condition to describe the stability region. The case without relative velocity is completely describe and simplest necessary conditions are given for the general case. We finally propose some numerical simulations that illustrate the stability property: even if the stability notion that we investigate is not exactly a cosntraint of convexity, a numerical maximum principle is observed if the parameters are in the stability region whereas numerical oscillations appear (in particular for non smooth profils) if the parameters are not there.
Moreover, relative velocities modify the stability array in a nontrivial manner. For instance, intuition might have suggested that the stability region for the relative velocity equal to the advection velocity contains all the others but it is realy not the case. For a given advection velocity, relative velocities cannot be used to increase the value of the first relaxation parameter, the one involved in the numerical diffusion.
The non-negativity of the relaxation matrix could be extended to nonliear schemes. The theoretical study will be much more technical and has not been performed. Nevertheless, numerical experiments for the Burgers equation show that the behavior of the D1Q3 scheme, and in particular the possibility or not of oscillations, is analogous to the linear case.
References
- [1] F. Caetano, F. Dubois, and B. Graille. A result of convergence for a mono-dimensional two-velocities lattice Boltzmann scheme. arXiv:1905.12393 [math], May 2019. arXiv: 1905.12393.
- [2] P. J. Dellar. Nonhydrodynamic modes and a priori construction of shallow water lattice Boltzmann equations. Physical Review E, 65(3):036309, February 2002.
- [3] D. d’Humières. Generalized Lattice-Boltzmann equations. In Rarefied Gas Dynamics: Theory and simulation, volume 159, pages 450–458. AIAA Progress in astronomics and aeronautics, 1992.
- [4] F. Dubois, T. Février, and B. Graille. Lattice Boltzmann schemes with relative velocities. Communications in Computational Physics, 17(4):1088–1112, 2015.
- [5] F. Dubois, T. Février, and B. Graille. On the stability of a relative velocity lattice Boltzmann scheme for compressible Navier–Stokes equations. Comptes Rendus Mécanique, 343(10-11):599–610, 2015.
- [6] F. Dubois and P. Lallemand. Towards higher order lattice Boltzmann schemes. Journal of Statistical Mechanics: Theory and Experiment, 2009:P06006, 2009.
- [7] M. Geier, A. Greiner, and J. C. Korvink. Cascaded digital lattice Boltzmann automata for high reynolds number flow. Physical Review E, 73:066705, 2006.
-
[8]
I. Ginzburg.
Truncation Errors, Exact And Heuristic Stability Analysis Of Two-Relaxation-Times Lattice Boltzmann Schemes For Anisotropic Advection-Diffusion Equation.
Communications in Computational Physics, 11(5):1439–1502, May 2012. - [9] P. Lallemand and L.-S. Luo. Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability. Physical Review E, 61:6546–6562, 2000.
- [10] S. Osher. Convergence of generalized muscl schemes. SIAM Journal on Numerical Analysis, 22(5):947–961, 1985.
Comments
There are no comments yet.