Have you heard about the Lattice Boltzmann Method? No?
…Or maybe you have, but you would like to know more details about it?
Then this article is for you!
My friend Nicolas Maquignon from France who is a PostDoctoral Researcher has massive experience in this area and he decided to write and share his knowledge about it with you all on this very blog.
I hope his article will bring new ideas to you!
Warning though… Nicolas made a big effort to keep the article simple, but it is still pretty advanced stuff…
so if you don’t understand fully, feel free to ask questions, comment or directly contact Nicolas (his contact is at the bottom of the page)
Here's a quick profile of Nicolas:
About the author:
Dr. Nicolas Maquignon
I am a PostDoctoral researcher who works in a French start-up located in Toulouse (south of France). My job is to develop physical models of explosions within a new FEA software using smoothed particle hydrodynamics methods (SPH)… and it’s really amazing!
Dr. Nicolas Maquignon,
Physical Model Development Specialist
1- Creating your own CFD solver with the Lattice Boltzmann Method, an alternate CFD model, from physics to computational physics
Before working on SPH, I developed equations and test cases simulation code in the field of the Lattice Boltzmann Method (LBM) for 5 years.
I used to study LBM equations for LNG (Liquefied Natural Gas) phase transition and scientific computing during my PhD and I could go through all this fascinating work thanks to a background in mechanical engineering.
I am always opened to discussions and debates in order to improve my understanding of numerical methods and also to help you to make your own idea... so if you have questions, don't hesitate!
Why the Lattice Boltzmann Method, you may ask... ?
The Lattice Boltzmann Method is a powerful method which is becoming famous and also a serious contender with traditional models used in CFD such as the finite element or finite volume methods…. and that’s why I decided to share all I know about it in this article today (Well, not ALL… obviously, but truly a lot)
LBM is a method which attracts many researchers and engineers, because it allows the incorporation of a lot of physics and because it has the advantage to be highly parallelisable.
Understand that LBM is very well suited for scientific CFD computing with Graphics Processing Units (GPU), which under certain circumstances can really accelerate the simulations.
Using this method, many fluid behaviour become accessible, like turbulence, acoustics, phase transition, thermal fluids… It is used in the field of CFD modeling for transport ( ie train-flight), meteorology, porous media… And many other topics.
2- Boltzmann and links between microscopic and macroscopic scales
Do you know the difference between a fluid seen at the human scale (the macroscopic scale) and the same fluid seen at a molecular scale (the microscopic scale)?
A fluid at the macroscopic scale is a continuous material with certain properties (density, speed) and its behaviour can be described by the famous Navier-Stokes equations... at the microscopic scale, the fluid is seen as particles and thus quantities such as density or speed take a whole other meaning.
Most CFD programs are solvers that numerically integrate the macroscopic Navier-Stokes equations, using finite difference, finite volume, or finite element methodologies (plus a graphical interface maybe).
Never heard about Navier-Stokes Equations?
The Navier Stokes equations describe the temporal and spatial evolution of density, momentum and energy.
They can be solved in a Lagrangian or Eulerian framework, depending on the problem to study.
In this equation(s), density and velocity are continuous macroscopic variables, which implies that:
- They can be described by continuous functions
- We can measure them by means of appropriate instruments at the human scale
Ex: A simple balance and a dosing glass make it possible to measure a mass, a volume, and thus a density. More sophisticated methods as particle tracking velocimetry allow us to study the velocity field of a fluid. Pressure sensors can inform us about the pressure at some strategic points.
In a general way, at a macroscopic scale, it is mostly possible to numerically describe the behavior of a fluid, with more or less moderate computational resources, depending on the geometry and the conditions of turbulence and compressibility.
However, even if the events seem to happen in a globally coherent way on our scale, things are quite different when we try to study this same fluid on a microscopic scale, that is to say in a dimension where one can discern its particles, and where the apparent continuity on our scale is broken.
...but why the heck should I care about the microscopic scale and the fluid's molecules?
Indeed, why not staying with the standard macroscopic methods developed dozens of years before that until now seemed to suffice?
First, a fluid seen at its microscopic scale IS an ensemble of molecules which have an unexpected behavior (and unpredictable behaviour)... we can't just ignore that!
Second, establishing links between micro an macro scales improves our understanding of the fluid behaviour, and thus we may become able to use this knowledge to improve current known methods.
If you think about it, many of the scientific and engineering breakthrough of the last century have been due to some scientists such as Einstein who have been willing to look closer at the real physical behaviour ...
If scientists weren't willing to do that, we would still probably believe that the earth is flat... ;-)
The LBM method is a very promising CFD model, which will play an important role in the next following years in the field of CFD (and It also is a very interesting topic for those interested in theoretical and applied physics...)
Understanding the LBM can also open the gates to a deeper understanding of the calculations done by the solver of your FEA software.
Moreover, I actually encourage you to develop your own CFD code, using any known method...
Never just trust blindly softwares made by others (Why? Read this article from Cyprien to understand)…you will truly learn a lot by writing your own code.
Now... What are the big problems with the microscopic scale?
On this microscopic scale, certain notions which are familiar to us disappear...
Just think about it:
- How can you define a density, while a large vacuum surrounds the molecules of the fluid ?
- How do you define a velocity, while in an infinitesimal volume, the particles move in a disordered way and in all directions ?
- And...what is the link between the movements of particles and the observable properties of a fluid?
Those are the questions the brilliant physicist Ludwig Boltzmann [1844-1906] tried to answer.
He is the one who established first all the links and equations that lead to the Lattice Boltzmann Method.
That's what we will talk about next.
3- Case of a gas at steady state and zero velocity
Consider the case of a volume of air at rest in a room, without apparent movement, and static.
From our human point of view, this volume of air is inert and in the absence of thermal or mechanical disturbances from the external environment, this volume of air can remain in this state for a very long time… Let’s now zoom in on an infinitesimal portion of this volume of air.
What do we observe? Well, we observe everything except a static state ;-)
The particles are in motion relative to the reference frame of the room, and this in a completely disordered and completely random manner, at least in appearance...but appearances are often misleading …
For the density of our volume of air, this is fairly easy to understand how to calculate it...
Because for a sufficiently large volume, defining a density actually means counting the number of particles:
mp is the mass of a particle and Np is the number of particles. The volume V corresponds simply to the volume of the room.
Here, shocks between particles are not taken into account. They are in any other case considered as hard, elastic and have no effect on the mass of the particles. Indeed, within the limit of Newtonian classical physics, the mass is preserved during collisions of particles.
Now ... how can we explain the velocity in a microscopic scale where velocities are so high and disordered that we see no global movement of the fluid at the macroscopic level?
4- A fundamental statement : Isotropy
First, we need to assume that the velocity vector of a particle can point in any direction. This statement seems quite obvious, as there isn’t any reason why an entity could not move into a certain direction.
This also means that there is no priviledged direction for the movement of a particle belonging to a fluid at rest. This is due to the fact that our geometrical space has the same properties in any direction, in other words it is isotropic. This property is responsible for the fact that the speed vector sum of the particles is zero.
This zero velocity vector sum is the speed we measure at the human scale, which explains why for a fluid at rest, we measure a zero velocity though a lot of disorder is present at the microscopic scale.
To summarise simply:
Isotropy = Same probability a particle moves to any direction = macroscopic velocity is zero
To demonstrate that, we need to do a bit of Maths ;-)
You can skip to the next section without much consequences if you’re not willing to go into some further mathematical details about it...
This property of isotropy is going to allow us to use a bit of math, to show that the vector sum of the particle velocities is zero. Let us assume an ensemble of particles « i » that have an individual speed norm and direction.
The speed vector sum V of this ensemble of particle can be written using Einstein summation convention :
In order to organise our algebraic development, we need to postulate that within a small range of speed norm [v-dv ;v+dv], every direction « k » is represented by a particle « j » owning the speed vj = v :
Thus, the analysis can be decomposed into a large number a small ranges [v-dv ;v+dv], correspondind to the spreading of speed magnitude. The calculations that are true for a certain speed range is also true for another speed range. The calculation can be simplified into :
Or with some simplified expression :
We have to suppose that a very large number of particles are present in our volume of study, in order to keep our approximations close to the continuum limit. This means that any direction or angle can be represented in a equi-probable proportion. This also means that for any a in [-1;1] :
Under the hypothesis that a very large number of particles are present, we could assume that the ensemble of particles behaves « almost » as in the continuum limit, and that the interval of vxk's is « dense ». Thus, we could make some strong assumption that we could arrange our sum of speed vector such as :
Thus, the vxk’s or vyk’s can be represented through an equi-probability distribution function on the [-1 ;1] interval. In the interval [-1 ;1], this distribution is simply :
Calculating the a mean value for vx, we have :
Calculating E(vx) we find E(vx)=0. This corresponds to the mean value of the vxk’s. Keeping in mind that :
So...that's it, we just demonstrated that the vector sum of particle velocities is 0 !
Though the vector sum is zero, this doesn’t mean that the mean speed magnitude value is also zero. The average speed is :
And the mean kinetic energy is :
We are now going to have a look at this kinetic energy and explain how it is related to temperature.
5- Air at 20°C and average kinetic energy
As I have just mentioned, though their directions are really disordered, the particles have a certain average magnitude velocity, which is non-zero.
This speed magnitude is around 400 m/s for air at a temperature of 20°C, as shown in the following figure showing speed distribution on a red curve.
I ploted this curve with Octave (GNU) to show the distribution of speed magnitude within a volume of air at rest at ambiant temperature. X axis reprensents a speed interval, as Y axis corresponds to a percentage of particles within a certain range of speed.
If the distribution function is called f(x), we can calculate the probability for a particle to own a certain velocity magnitude within a range [a ;b] :
Note that it was here necessary to introduce the notion of temperature when fixing an average velocity value of the particles. This is because the temperature is related to the velocity of the particles. On the following figure (green curve), I ploted the velocity distribution for air at 400°C. We can observe that the mean velocity is now around 600 m/s and that the distribution range is widened, compared to air at 20°C.
Knowing about this particular distribution speed leads us to the distribution of energy in the fluid.
First, we need to introduce the notion of kinetic energy. Kinetic energy is the energy that something have due to its mass and its speed.
Taking the example of a bullet going out of a rifle, it is not beyond imagination that this bullet has a certain (large) energy, regarding all the trouble it could cause when it encounters an obstacle. The energy of the bullet comes from the chemical energy that was embeded in the gun-powder that has been converted into speed for the bullet at the moment of detonation.
This « speed » energy is called kinetic energy. Boltzmann has established that there is in fact a link between the temperature and the kinetic energy of the particles, which is a function of the square of the velocity. Boltzmann tells us that for a monoatomic gas :
In this expression, the constant kb has appeared. This is a fundamental constant, called the Boltzmann’s constant.
The kinetic energy Ec that we get is in fact called the internal microscopic energy, as it contains the energy of the particles. If you try using this formula with some Argon, you’ll find a mean speed of 427 m/s.
6-Particule Position & Particle Velocity - Maxwell Distribution
The organization of particle velocities in disordered directions but with a fixed average magnitude, for a fixed temperature, can be further analyzed.
The mean energy may be fixed but nevertheless all the particles of an isothermal gas don't necessarily have the same kinetic energy.
In reality the velocities are distributed around a mean value, according to a certain organization. This is what we actually see in the previous red and green plots about speed magnitude distribution within a fluid at rest. This shows that the particles have a range of velocity magnitudes, distributed around a mean value which is located at the absciss of the maximum of the curve.
It is a priori not simple to find a probability law for a particle to have a certain velocity magnitude. However, for speed directions, the isotropy hypothesis, which leads to equiprobability of all directions, seems to hold, as we’ve just seen before.
To deal with the problem of speed magnitude distribution, we could postulate that a particle can have any kinetic energy, as long as the mean kinetic energy of the set of particles allows us to reconstruct a physically coherent macroscopic temperature, when summing over all the particles.
Only a fine observation of the properties of the particles, if it was possible, would teach us about the standard deviation around this mean kinetic energy (Appart from the averaged kinetic energy, the other statistical moments are also related to the macroscopic temperature).
Though, assuming that a particle can have any speed around an average value does not give us a precise idea about the distribution of velocities.
If we consider that each particle individually possesses a velocity according to an individually indeterminate probability law, the central limit theorem teaches us that for a large number of particles, the distribution law to be globally applied is the Gaussian normal law. Thus, the distribution of the speed magnitudes is done according to the normal law that follows :
This distribution law is the Maxwell equilibrium distribution, valid for a perfect gas.
The constant Rs is the specific gas constant, T is the temperature, e is the microscopic velocity of particles and u is the macroscopic velocity.
I used this formula to plot the red and green speed distribution curves, with u=0 for a fluid at rest.
For a perfect gas, the internal energy of the system is summed up in microscopic kinetic energy, because the potential inter-particle energy is considered as zero. A perfect gas is sufficiently dense for collisions to take place quite consistently.
7- Out-of-equilibrium distribution function and evolution
The notion of "phase space" is an elementary concept that needs to be understood in the field of Lattice Boltzmann methods.
First, you surely have heard about the usual geometric space, which is modeled with three coordinates (x,y,z) in the Cartesian framework.
Cartesian space (source)
This space isn’t sufficient when we come to the description of particles, because we’ve seen how important the speed of the particles is. Indeed, the macroscopic properties does not only depend on spatial distrbution of particles, but it also depends on the speed distribution.
If we are to study the variations of these spatial and velocity distributions through time, we need to include the time in our coordinates.
The phase space corresponds to the gathering of geometric and velocity spaces.
7-b-Phase space distribution function
At equilibrium, the knowledge of the microscopic properties of air, or more generally of a perfect gas, is now accessible to us, thanks to the Maxwell equilibrium distribution that we have introduced earlier in the article.
However, the study of a fluid rarely summarizes to its equilibrium properties.
Limiting ourselves to equilibrium properties would forbid the study of many cases.
We must therefore open ourselves to the study of properties that are out of equilibrium by simply introducing the distribution of particles out of equilibrium into the phase space.
This distribution of particles in phase space can be defined by a function f:
It corresponds to the probability for a particle at a position [x] to possess a velocity [e] at a time [t].
From here you can skip to the next section if you’re not willing to go into some further mathematical details about it:
Following the principle of inertia, we know that the velocity distribution is preserved over time, in the absence of an external force field and in the absence of interparticle collisions.
In the absence of compression or relaxation, the spatial distribution remains essentially unchanged over time.
At low Mach and in the absence of collision, the differential of f is conserved. We can thus deduce the following equality:
The transport expression of f can then be obtained with a Taylor development:
In absence of external forces:
In the case of collisions, this equality must be modified by taking into account the variations of the probability distribution in the phase space. This factor will be called a collision operator, and we will not deal with its complex expression in detail.
We’ve just written here the Lattice Boltzman Equation (LBE).The collision operator can nevertheless be seen as a rate of variation, corresponding to a relaxation towards the equilibrium distribution. Thus, a reasonable approximation of this collision factor can be established :
This collision operator is know in the LBM community as the "BGK" collision operator or « Bhatnagar-Gross-Krook ».
8- Lattice and discrete form - advection/collision
Having now a transport equation, a discretization can be envisaged, in order to achieve a numerical integration.
In this article, we will not go into the mathematical detail of this discretization, because it is tedious and is not really essential for an elementary understanding of the LBM. Interested readers can refer to articles in the scientific literature (you can contact me if you’re interested).
It is important to understand that the space to be discretized is not a classical four-dimensional space. Indeed, since the distribution function is defined as a function of position, speed and time, it is not possible to use a single spatial mesh on which time iterations are to be calculated. The mesh must correspond to a whole phase space discretisation.
This means that the mesh must include geometrical nodes, but also « speed nodes ». Thus, one can imagine the use of a spatial mesh, in which a set of discrete velocities are represented.
This adds a constraint on our spatial mesh, which can not be realized in a completely free way because of the presence of speeds.
In the same way that a finite number of discrete spatial nodes can discretise the geometry of the problem we want to solve, we must define a set of discrete velocities representative of the set of velocities a particle could take. In practice and for 2D cases with little complexity, a set of 9 speeds are sufficiently representative.
Indeed, 9 speeds provide sufficient isotropy to the model and it is sufficient for weakly compressible fluid applications.
In a Cartesian mesh, these velocities are taken such that a particle present on a node can remain on this node, or else go to one of the 8 nodes that are located in its direct neighborhood.
2D problem is discretized with the Lattice Boltzmann framework (Source)
This case is called the D2Q9 lattice: that means two spatial dimensions and nine discrete velocities. Other types of LBM networks can be established, in 2D and 3D. For LBM lattices, the DnQm denomination is usually used. A DnQm Lattice means that the lattice has n spatial dimensions and m discrete speeds.
The black dots represent spatial nodes that are occupied by particles. A site does not contain only one particle but a certain number that depends on the local density.
In a time step, a particle can only move to its nearest neighbors sites, with one of the discrete speeds that are here noted ei, with i from 0 to 8. We all together have 9 discrete velocities.
The fi’s represent the number of particles that move with a discrete speed ei. For a site containing N particles, fi*N particles move with a speed ei.
Once the ensemble particles have moved (advection) the one of the neighbouring site, they will collide with other ensemble of particles that arrive on the same site at the same time.
We have here identified the two fundamental steps of the LBM algorithm : advection and collision. And of course, the collision can be calculated with the collision operator that we have introduced earlier. The advection corresponds to the left had side of the LBM, and the collision corresponds to the right hand side of the LBE depicted before.
9- Discrete Equation and Algorithm
Since the different notions relating to the discretization of the phase space are now posed, it remains to describe the discrete form of the Boltzmann equation on a lattice of type DnQm.
Directly, this equation takes the following form, for a direction i of discrete speed:
The left term is a term of transport or advection. A Taylor development can also highlight a form of Lagrangian transport. The term on the right is the collision operator, it corresponds to the rule that applies when two particles arrive at the same time on a geometric node.
Based on this Boltzmann equation on lattice - LBM - an algorithm can be defined.
This is illustrated in the previous figure, and the steps are:
- Selection of a site
- Advection of incoming particles from direct neighbourhood.
- Collision of the new set of particles that have just been advected.
- Post-collision rearrangement
- Next site selection
10- An example : the pressure pulse
For this example, I wrote a D3Q19 LBM solver using C language under some ubuntu environment (If you’re interested in getting the code, please contact me. )
The pressure pulse is an elementary benchmark you can use to validate an implementation for some cases.
Here, the topic of pressure pulse is rather in the field of acoustics, as it can model very small pressure variations. This doesn’t mean that all we’ve seen before is only available for acoustics!
All the tools we’ve seen are quite powerful, and it can in general solve cases for low Mach laminar fluids. Some small modifications can bring the model up to soft turbulence modeling. Of course, if you’re interested in more complex flow, you should bring some further modifications.
The initialisation of the case is quite straightforward : A pressure peak is initialised at the center of the domain. I mean the density is initialised in a uniform manner, except at the center where a small pressure shoot is added to the average density, as it can be seen on the first 2D view of the following figures. The whole domain is initiaded with zero velocity. The boundary conditions are periodic. The mesh is cartesian and uniform has depicted earliear.
The spatial step is of 1.0 cm, and the time step is of 17.28 e-6 s. The speed of sound is set at 334 m/s, which is the speed of sound of air at 20°C. The physical simulation time is of 0.5 ms. I ploted illustration using Paraview (GNU) in order to depict initialisation and results :
Initialisation of density field :
Density field after 0.5 ms of physical simulation time :
We can observe that the initial density peak (linear to pressure) spreads to a slighly higher pressure wave, followed by a rarefaction wave : this is typical of that kind of test case, and it validates our implementation in a qualitative manner. The two following plots show a radial profile of the density (there is a spherical symmetry), at the initialisation and after the 0.5 ms of physical time:
We can observe the higher pressure peak followed by the rarefaction peak. And if we measure the distance that the high pressure speed has done during the simulation, and we divide it by the physical time, we find 334 m/s… Which is indeed the speed of sound that we set as parameter. This validates the code in a quantitative manner, though a comparison to analytical solution would still be interesting (for the readers who would like to go further). The two last figures show density and velocity fields at the end of the simulation in 3D.
That’s all basically for this article… and that's a LOT of knowledge to unpack ;-)
Take all the time you need to read and if you want to share your ideas, please just write what you think in the comments and let me know
Thank you so much for reading, I hope you learned a lot!