Structural stability of Lattice Boltzmann schemes

by   Claire David, et al.

The goal of this work is to determine classes of traveling solitary wave solutions for Lattice Boltzmann schemes by means of an hyperbolic ansatz. It is shown that spurious solitary waves can occur in finite-difference solutions of nonlinear wave equation. The occurence of such a spurious solitary wave, which exhibits a very long life time, results in a non-vanishing numerical error for arbitrary time in unbounded numerical domain. Such a behavior is referred here to have a structural instability of the scheme, since the space of solutions spanned by the numerical scheme encompasses types of solutions (solitary waves in the present case) that are not solutions of the original continuous equations. This paper extends our previous work about classical schemes to Lattice Boltzmann schemes.



There are no comments yet.


page 1

page 2

page 3

page 4


Finite Difference formulation of any lattice Boltzmann scheme

Lattice Boltzmann schemes rely on the enlargement of the size of the tar...

Lattice Boltzmann Method for wave propagation in elastic solids with a regular lattice: Theoretical analysis and validation

The von Neumann stability analysis along with a Chapman-Enskog analysis ...

The propagation of transient waves in two-dimensional square lattices

The aim of this article is to study the attenuation of transient low-fre...

Simulation of the Geometrically Exact Nonlinear String via Energy Quadratisation

String vibration represents an active field of research in acoustics. Sm...

Does the multiresolution lattice Boltzmann method allow to deal with waves passing through mesh jumps?

We consider an adaptive multiresolution-based lattice Boltzmann scheme, ...

On the numerical overshoots of shock-capturing schemes

This note introduces a simple metric for benchmarking shock-capturing sc...

Unexpected convergence of lattice Boltzmann schemes

In this work, we study numerically the convergence of the scalar D2Q9 la...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction:

The lattice Boltzmann method (LBM

) is used for the numerical simulation of physical phenomena, and serves as an alternative to classical solvers of partial differential equations. The primary domain of application is fluid dynamics; it is specially used to obtain the numerical solution of the incompressible, time-dependent Navier-Stokes equation.

The strength of the Lattice Boltzmann method is due to its ability to easily represent complex physical phenomena, ranging from multiphase flows to fluids with chemical reactions. The principle is to "mimic" at a discrete level the dynamics of the Boltzmann equation. Since it is based on a molecular description of a fluid, the knowledge of the microscopic physics can directly be used to formulate the best fitted numerical model.

This method can be regarded as either an extension of the lattice gas automaton (LGA) [5], [6],[7], or a special discrete form of the Boltzmann equation from kinetic theory. Although the connection between the gas kinetic theory and hydrodynamics has long been established, the Lattice Boltzmann method (LBM) needs additional special discretization of velocity space to recover the correct hydrodynamics. Due to the very same reason, the LBM works exactly opposite traditional CFD methods in deriving working schemes: LBM uses Navier-Stokes equations as its target while traditional CFD methods use Navier-Stokes equations as their starting point.

The LBM is applicable to the isothermal flow regime, i.e., the weakly compressible, low-Mach-number limit. This flow regime is traditionally treated as "incompressible", although there are CFD methods constructed to compute the Navier-Stokes equations in this regime. The argument for treating very low-Mach-number flows as incompressible is pragmatic rather than physical. The direct calculation of the isothermal Navier-Stokes equations requires time steps sufficiently small to resolve acoustic waves across a computational cell. This time step may be vastly smaller than the time scales of interest for the bulk fluid motion. Thus the computational cost of the many additional time steps required by an isothermal calculation may be vastly higher than the cost of an incompressible calculation. Of course, in reality there is no fluid or flow that is absolutely incompressible (i.e., with infinite acoustic velocity).

2 The Lattice Boltzmann method

The lattice Boltzmann (LB) method follows the same idea as its predecessor the Lattice Gas Automata (LGA

) when it also considers the fluid on a lattice with space and time discrete. Instead of directly describing the fluid by discrete particles and, thus Boolean variables, it describes the fictitious system in terms of the probabilities of presence of the fluid particles. A lattice Boltzmann numerical model simulates the time and space evolution of kinetic quantities, the particle distribution functions

, , .

The lattice Boltzmann equation is obtained by ensemble averaging the equation


where denotes the average number of particles at space position and time .
The system is supposed to satisfy the Boltzmann molecular chaos hypothesis, i.e. the fact that there is no correlation between particles entering a collision. Thus, the collision operator can be expressed as , which leads to the Lattice Boltzmann equation:


where, for , denotes the probability to have a fictitious fluid particle of velocity entering lattice site at time . The are also called the fluid fields, or the particle distribution functions.

The collision operator is normally a non-linear expression and requires a lot of computation time [13]. In a big lattice, e.g. 3D model, the computation becomes impossible even on a massively parallel computer. To overcome this problem, Higuera et al. [8], [9] proposed to linearize the collision operator around its local equilibrium solution to reduce the complexity of the operation. Using this idea, Bhatnager, Gross and Krook introduced the BGK lattice (LBGK) [10], in which the collision between particles is described in terms of the relaxation towards a local equilibrium distribution. The LBGK is considered to be one of the simplest forms of the Lattice Boltzmann equation and is mathematically expressed as


where is the single relaxation time, which is a free parameter of the model to determine the fluid viscosity, and , , denote the local equilibrium functions, which are functions of the density and the flow velocity .

In the lattice Boltzmann method, the space variable vacetor is supposed to live in a lattice included in an Euclidian space of dimension , .
The velocity belongs to a finite set composed by given velocities , , , chosen in such a way that


where denotes the time step.

The set of velocities is supposed to be invariant by space reflection, i.e.:


The numerical scheme is thus defined through the evolution of the population , with and towards a distribution at a new discrete time.
The scheme has two steps that take into account successively the left and right hand sides of the Boltzmann equation (1). The first step describes the relaxation of particle distribution towards the equilibrium. It is local in space and nonlinear in general. The second step is the collision process. Both steps have to be computed separately.

In general LB model, denoted by ,

, the actual macroscopic density, which is function of the space vector

, and the time variable , is obtained as a summation of the microscopic particle distribution functions:


The macroscopic velocity is calculated as the average of the microscopic velocities , , weighted by the related distribution functions:


In each time step, in each node, particles are streamed on to the neighbouring nodes, which thus lead to new distribution functions :


Also, once in each time step, the particles in each node collide, which is modeled as a relaxation of the distribution functions towards the equilibrium distributions



where is the lattice sound speed, such that:

and where the weights , , are determined by the velocity set. Of course, one has:

The lattice Boltzmann equation takes both steps into account:

The right-hand side gives thus the distribution of particles once collisions have occured, while the left-hand side represents particles appearing in neighbouring nodes in the next time step.

For our study, it is important to take into account the way the related algorithm is implemented:

  1. Initialization step of , , and the microscopic densities , , .

  2. Streaming step: densities are moved towards in the direction , , using (8).

  3. Computation step of the macroscopic variables and .

  4. Computation step of the microscopic variables , .

  5. Collision step, in order to obtain the updated densities , , using (9).

  6. Return to steps ii. to v.

3 Spurious lattice solitons

The discrete solution associated with the LB numerical scheme will admit spurious solitary waves, and therefore spurious local energy pile-up and local sudden growth of the error, if the discrete relation used to implement the scheme is satisfied by a solitary wave.

Following [18], [1], [2], [3], [4], and using [19], [20], [21], [22], [23], [24], [25], [26], we search solitary waves solution components under the form:



and where , , , , , are real constants.

The question one may ask is how to switch from the scheme, which is discrete, to this continuous expression ? In so far as if one can find a solution of the above form (10) that fits each step of the scheme, is the answer. This is how we will proceed.

Starting from the fact that once the discrete populations are known, the main fluid quantities can be obtained by simple linear summation upon the discrete speeds, it is legitimate, for those densities, to be of the above form (10), for the same reasons as previously.

3.1 A one-dimensional example: the D1Q3 model

In the following, we examine a one-dimensional case, given by the flow of an incompressible fluid, in a domain , , obeying the Navier-Stokes equation, in order to test the D1Q3 model, which has a zero velocity and two oppositely directed ones, moving thus the fluid particle to the left and right neighbour lattice sites.

Figure 1:

Discrete velocity vectors of the D1Q3 model

The velocities are:

while the weights are:

The conservation of mass quantity writes:


Since the density is a conserved quantity, we can suppose it to be constant, and normalized. In order to determine wether there exist, or not, solitary waves solutions, one requires to take into account related boundary conditions, in the case of an imposed velocity. Then, using (10), we search the densities , , under the form:

and the horizontal velocity as:

where, for , , , , , , are real constants to be determined.

By substituting those latter expressions in (11), one gets an equation of the form:

where we denote by and functions of , , , . By independance of the terms in sech and tanh, one gets, for :

i.e. a linear system of two equations, the unknowns of which are the , , , . It happens that in this D1Q3 model, the system is a very simple one:

which admits several sets of solutions. For sake of simplicity, we have choosen the following one:

It thus exhibits the existence of lattice solitons, related to the discrete numerical scheme, of the form


Figure 1 displays the related lattice solitary wave, as a function of the mesh points ; denotes the spatial number of mesh points in the -direction, the time one.

Figure 2:

A "lattice solitary wave",

as a function of the mesh points, in the case of the D1Q3 model

3.2 A two-dimensional example: The case of an isothermal Poiseuille flow between two parallel plates

The case of an isothermal Poiseuille flow, driven by a pressure gradient in the -direction, between two horizontal, parallel, plates is an interesting one: if the flow is a two-dimensional one, yet, since it is a longitudinal one, one only has to deal with the longitudinal velocity component, , which happens to depend only on the vertical space variable . In our case, that means that if there exists solitary waves, they also will be horizontal ones, and functions of . Also, since the exact analytical expression of the velocity is known (the flow develops a parabolic velocity profile), results can be tested profitably. In the following, we will denote by the distance between the two plates.

The velocity vectors of the D2Q9 model are:

Figure 3:

Discrete velocity vectors of the D2Q9 model

In the specific case of the D2Q9 model, one has:

The weights of the population are given by:

As in the above, in order to determine wether there exist, or not, solitary waves solutions, one requires to take into account related boundary conditions, in the case of an imposed velocity. We thus place ourselves in the case of the Zhou-He scheme [27]. The boundary between the fluid and the plates overlap the domain nodes. If one chooses to solve the undetermined population of the south boundary, , , , conservation of the mass and quantity of movement yield:


One can then, as in [ZhouHeScheme] use the bounceback of non-equilibrium in the direction normal to the wall:

which leads to:


Let us now study the existence of solitary waves solutions of (13) and (14). We assume the density to be constant, and normalized. Then, we search the densities , , and the longitudinal velocity component , as functions of the vertical space variable , and of the time one, , under the form:

where, for , , , , , , are real constants to be determined.

The boundary conditions on the plates require:


By substituting those latter expressions in (13) and (14), one gets a system of six equations of the form:

where, for , we denote by and functions of , , , . By independance of the terms in sech and tanh, one gets, for :

i.e. a linear system in 12 equations, the unknowns of which are the , , , .

Taking into account the boundary conditions (15), resolution with a symbolic calculus tool (Mathematica for instance) shows that this system admits several sets of solutions. For sake of simplicity, we have choosen the following one:

Figure 2 displays the related lattice solitary wave, as a function of the normalized vertical coordinate, which can be compared to the exact analytical solution.

Figure 4:

A "lattice solitary wave",

in the case of a Poiseuille flow between two parallel plates (in red),

compared with the analytical exact solution (dashed plot).

4 Concluding remarks

The existence of spurious numerical lattice solitary waves for Lattice Boltzmann schemes has been proved. Such lattice solitary waves, which are not solutions of the exact continuous original equation, nevertheless satisfy the numerical scheme, appearing as parasitic solutions of the correct one. Such schemes will be referred to as structurally instable ones. A solution to avoid such spurious solitary waves, would be to "lock" the scheme by adding a test in the main loop.
Such spurious solitary waves have constant energy, and therefore the numerical error norm does not vanish at arbitrary long integration times on unbounded numerical domains.

One may ask why do such solutions occur ? A partial answer can be found in the equivalent continuous equation, the principle of which is the following: for a small time step , and the associated space scale , a Taylor expansion in leads to establishing equivalent continuous equations as a formal limit. This matter has been explored by F. Dubois in [28], where it is proved that, at the first order in , the equivalent equation is the Euler equations, which are known to admit solitary waves solutions (see, for instance, ).


  • [1] Cl. David, P. Sagaut, Structural stability of discontinuous Galerkin schemes, Acta Applicandae, under press.
  • [2] Cl. David, P. Sagaut, Structural stability of finite dispersion-relation preserving schemes, Chaos, Solitons and Fractals, 41 (4), 2009, 2193-2199.
  • [3] Cl. David, P. Sagaut, Spurious solitons and structural stability of finite difference schemes for nonlinear wave equations, Chaos, Solitons and Fractals, 41 (2), 2009, 655-660.
  • [4] Cl. David, R. Fernando, Z. Feng, A note on "general solitary wave solutions of the Compound Burgers-Korteweg-de Vries Equation", Physica A: Statistical and Theoretical Physics, 375 (1), 2007, 44-50.
  • [5] J. Hardy, Y. Pomeau and O. de Pazzis, Time Evolution of a Two-Dimensional Classical Lattice System, Physical Review Letters, 31, 1973, 276-279.
  • [6] D. d’Humières, P. Lallemand and U. Frisch, Lattice gas models for 3D-hydrodynamics, Europhysics Letters, 2(4), 1986, 291-297.
  • [7] U. Frisch, B. Hasslacher and Y. Pomeau, Lattice gas automata for the Navier Stokes equation , Physical Review Letters, 56(14), 1986, 1505-1508, 1986.
  • [8] F. Higuera, J. Jimenez, and S. Succi. Boltzmann approach to lattice gas simulations, Europhys. Lett, 9, 1989, 663.
  • [9] F. Higuera, J. Jimenez, and S. Succi. Lattice Gas dynamics with enhanced collision. Europhys. Lett, 9,1989, 345.
  • [10] P. Bhatnagar, E. Gross and M. Krook. A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems , Physical Review, 94, 1954, 511-525.
  • [11] Y. Qian, D. d’Humières, and P. Lallemand, Lattice BGK models for Navier- Stokes equation, Europhys. Lett, 17, 1992, 470 84.
  • [12] D. d’Humières, Generalized Lattice-Boltzmann Equations , in Rarefied Gas Dynamics: Theory and Simulations, AIAA Progress in Astronautics and Astronautics, 1959, 1992, 450-458.
  • [13] B. Chopard, A. Dupuis, A. Masselot, and P. Luthi. Cellular Automata and Lattice Boltzmann techniques: An approach to model and simulate complex systems, Advances in Complex Systems, 5, 2002, 103 246.
  • [14] B. Chopard and M. Droz, Cellular Automata Modeling of Physical Systems, Cambridge University Press, 1998.
  • [15] M. B. Reider and J. D. Sterling, Accuracy of Discrete-Velocity BGK modor the simulation of the incompressible Navier-Stokes equations. Comput. fluids, 24(4), 1995, 459 467.
  • [16] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability, Physical Review E, 61, 2000, 6546-6562.
  • [17] Shan, Xiaowen and Yuan, Xue-Feng and Chen, Hudong, Kinetic theory representation of hydrodynamics: a way beyond the Navier Stokes equation, J. Fluid Mech., 550, 2006, 413-441.
  • [18] Z. Feng, G. Chen, Solitary Wave Solutions of the Compound Burgers-Korteweg-de Vries Equation, Physica A, 352, 2005, 419-435.
  • [19] B. Li, Y. Chen Y. and H.Q. Zhang, Explicit exact solutions for new general two-dimensional KdV-type and two-dimensional KdV Burgers-type equations with nonlinear terms of any order, J. Phys. A (Math. Gen.), 35, 2002, 8253 8265.
  • [20] G. B. Whitham, Linear and Nonlinear Waves, Wiley-Interscience, New York, 1974.
  • [21] M. J. Ablowitz, H. Segur, Solitons and the Inverse Scattering Transform, SIAM, Philadelphia, 1981.
  • [22] R. K. Dodd, J.C. Eilbeck, J. D. Gibbon, H.C. Morris, Solitons and Nonlinear Wave Equations, London Academic Press, London, 1983.
  • [23] R.S. Johnson, A Modern Introduction to the Mathematical Theory of Water Waves, Cambridge University Press, Cambridge, 1997.
  • [24]

    E.L. Ince, Ordinary Differential Equations, Dover Publications, New York, 1956.

  • [25] G. Birkhoff, G.C. Rota, Ordinary Differential Equations, Wiley, New York, 1989.
  • [26] A.D. Polyanin, V. F. Zaitsev, Handbook of Nonlinear Partial Differential Equations, Chapman and Hall/CRC, 2004.
  • [27] Q. Zou, and X. He, Pressure and velocity boundary conditions for the lattice Boltzmann, J. Phys. Fluids, 9, 1997, 1591-1598.
  • [28] F. Dubois, Equivalent partial differential equations of a lattice Boltzmann scheme, Computers and Mathematics with Applications, 55, 2008, 1441 1449.