# Lattice Boltzmann Method: How is shear flow handled in D2Q5?

I've implemented 2-dimensional, incompressible, high-reynolds fluid-flow using the Lattice Boltzmann Method on a D2Q9 lattice. My main goal is just visual plausibility, not quantitative accuracy. The simulation is initialised in equilibrium and in a small portion of space in the center a constant force to the right is exerted onto the fluid (to simulate a small "pump" to test the flow). The flow nicely moves rightward with emerging vorticies at the top and bottom of the flow as well as turbulence due to shear forces, as expected.

In order to speed up the simulation I experimented with using a D2Q5 lattice instead, but visual plausibility breaks down as no vorticies or any sort of turbulence emerges, just a boring narrow laminar flow to the right. I'd assume the reason is that in D2Q9 shear forces are considered due to the exchange with diagonally neighboring cells which are missing in D2Q5. Therefore the latter does not capture shear forces which are essential for turbulence, vorticies and so on. Since D2Q5 seems to be a practical choice in a bunch of papers I strongly assume I'm missing something here, as completely neglecting shear forces can not result in a physically sound (and not even visually pleasing) flow. See here for comparison:

Can somebody explain to me how D2Q5 is supposed to capture shear forces ? I really don't understand it conceptually.

• Physicists won't be able to answer this easily, consider asking an altered version of this question on scicomp.stackexchange.com. – whn Oct 10 '19 at 14:36
• also looks like your equilibrium step needs to be changed for d2Q5, ie, the probability of water moving up and colliding with other particles causing perpendicular flow is not zero, where it looks like there is zero perpendicular flow happening in your D2Q5 version. – whn Oct 10 '19 at 14:38
• Hey Lenny, the answer whn has given you is wrong. You simply can't simulate fluid flow with a D2Q5 lattice, but you can use it to simulate advection-diffusion of scalars. In LBM the number of velocities can't be chosen arbitrarily - it is connected to the conserved quantities: A simple diffusion equation requires the amount of substance to be conserved, an incompressible simulation requires the conservation of mass and momentum and a thermal/compressible simulation additionally conservation of energy. – 2b-t Aug 1 '20 at 1:09
• Therefore the requirements for the lattice discretisation increase, a D2Q5 suffices for scalar advection-diffusion, but for incompressible fluid flow you will need a D2Q9 and for fully compressible simulations you generally will need a very high number of discretisation speeds (multi-speed lattices such as a D2Q37). If you are still interested I can take my time and try to give you a more complete answer. Sorry I did not see your question some month earlier. – 2b-t Aug 1 '20 at 1:15
• @Lenny Don't worry, I will give you a longer and more detailed answer tomorrow. I would not recommend reading papers for getting familiar with any new topic: Publications are intended to be more like an open discussion between people that more or less familiar with the key aspects rather than introducing people unfamiliar with the basics. In my opinion there is only a single good book available on the topic by Krüger et al.. – 2b-t Aug 3 '20 at 20:00

Lattice-Boltzmann algorithms have several applications ranging from the traditional fluid flow to other topics such as scalar transport (such as species or energy equation) or radiation. Depending on the application you will find a lot of different discretisations:

• with/without rest node
• with a single speed for the non-rest node velocities or with multiple velocities (multi-speed lattices)

and of course with a different number of speeds.

At the first glance one might think that it is possible to use discretisations arbitrarily but this is not the case. The discretisation is strictly coupled with the conserved moments and thus a certain differential equation demands a certain isotropy of the lattice. The D2Q5 lattice can only be used for applications where only the first moment has to be conserved such as scalar transport (the scalar itself is preserved), while it lacks the isotropy for simulating incompressible fluid flow where the first two momenta (density and momentum) have to be conserved, in such a case you have to switch to the D2Q9 lattice. If you furthermore want to simulate compressible flows you might end up with as many as 37 distinct velocities in 2D.

The lattice-Boltzmann algorithms have their background in kinetic theory of gases, precisely the Boltzmann equation. I have already written some posts on this in case you are not familiar with it:

This equation describes the propagation of the single-particle probability distribution, an extended concept of density, due to the free motion and collisions

$$f(\vec x, \vec \xi, t) := \frac{d N_p}{d \vec x \, d \vec \xi}$$

on a microscopic scale. Contrary to the macroscopic description with density, velocity and pressure, one looks at probability density function that describes the probability of finding a certain particle with a particular momentum at a particular point in space. Macroscopic variables like density and velocity are not directly accessible, they emerge as expected values from the distribution. One basically averages over the underlying particles to obtain a single macroscopic value at a particular point in space. The velocity space in this case is continuous, you can have virtually any velocity magnitude and any direction.

Boltzmann published this equation in 1872 at a time where the atomistic view was not widely accepted. People were using the continuum-based approaches based on the Navier-Stokes equations (I have written posts on how to derive the stress tensor and their interpretation in case you are interested). Boltzmann was not able to convince the scientific community of his theory and he died struggling with depression pretty tragically by suicide in 1906. Only after his death his theory would gain wide acceptance but even when it did, one question still remained open. How can one derive a macroscopic description from this equation? The equation tells you what happens after a collision, how the distribution is re-destributed but what happens on much larger time and spatial scales?

Only several years later in 1916 and 1917 respectively two mathematicians, Sydney Chapman and David Enskog came up with a solution. They found out that - as continuum-based flows are assumed to be close to equilibrium - the macroscopic behaviour for dense fluids ($$Kn := l_{mfp} / L \to 0$$) can be recovered by applying an asymptotic perturbation analysis ($$\epsilon \propto Kn$$),

$$f = \sum_{n=0}^{\infty} \epsilon^n f^{(n)}$$

named after them as Chapman-Enskog expansion, around the equilibrium $$f^{(0)} = f^{(eq)}$$ and recombining the conserved quantities. Surprisingly this simplified explicit equation can be found to preserve the full Navier-Stokes equations for compressible flow: The equilibrium part preserves the Euler equations and the viscous contribution emerges from the terms of first and second order $$f^{(1)}$$ and $$f^{(2)}$$, while the meaning of higher-order terms, the Burnett and super-Burnett approximations, is not yet completely clear.

Motivated by the success of cellular automata, simplified molecular systems that mimic the macroscopic behaviour of complex system by prescribing rules on a microscale (maybe you have hear of "Conway's Game of Life"), researchers tried to come up with simple numerical systems that describe the complex flow of fluids. The lattice-gas automata were born, first from a purely empiric stand-point (trail and error) and at the end of the 1980s the framework of the Boltzmann equation was used to formalise the findings and construct a physically consistent fictional gas with convenient numerical properties by including statistical mechanics into it. From the previous experiments people knew that some sort of grid where particles collide with their neighbours on a regular grid were sufficient to to mimic the behaviour of a fluid. The problem with the derivation is that you have a distribution in space, velocity and time and you have not only to discretise the space and time but also quantise the velocity space and then do a discrete Chapman-Enskog expansion to recover its macroscopic equations. Clearly if you want to break down the algorithm in free-motion and collision each discrete velocity of a lattice has to move with such a velocity that it ends up at another node in the next time step. So the first approach was to choose the time $$\Delta_t$$ and space step $$\Delta_x$$ to one ("lattice units") respectively. The populations in this lattice-Boltzmann model (LBM) now travel with a single speed $$c = \Delta_x / \Delta_t$$ on a discrete equidistant grid and collide with populations from neighbouring nodes.

Already from kinetic theory it was known that for this purpose simpler collision models that describe the effect, the relaxation of a distribution towards equilibrium with a relaxation time $$\tau$$, rather than the collision itself, suffice as they preserve the leading moments linked to the hydrodynamic quantities. The simplest collision operator is given by the single relaxation time BGK collision operator

$$\Omega = \frac{\Delta_t}{\tau} (f_\alpha^{(eq)} - f_\alpha).$$

The derivation starts by discretising the time of this equation and integrating along a characteristics,

$$\frac{\partial f}{\partial t} + \vec \xi \boldsymbol{\cdot} \vec \nabla f = \frac{1}{\tau} (f^{(eq)} - f)$$

$$f(\vec x + \Delta_t \, \vec \xi, \vec \xi, t+ \Delta_t) - f(\vec x, \vec \xi, t) = \frac{\Delta_t}{\tau} \Big[ f^{(eq)}(\vec x, \vec \xi, t) - f(\vec x, \vec \xi, t) \Big] + \mathcal{O}(\Delta_x^2) + \mathcal{O}(\Delta_t^2)$$

in order to have an discretised equation to work on. Then it expands the equilibrium distribution for a certain macroscopic velocity (low Mach number expansion) in order to remove the exponential from the equilibrium distribution

$$f^{(eq)}=\frac{\rho}{\left(2 \pi R_m T\right)^{\frac{D}{2}}} \, e^{-\frac{(\vec \xi - \vec u)^2}{2 R_m T}} = \frac{\rho}{\left(2 \pi R_m T\right)^{\frac{D}{2}}} \, e^{-\frac{\vec \xi}{2 R_m T}} e^{-\frac{\vec u^2 - 2 \vec u \vec \xi}{2 R_m T}}$$

$$f^{(eq)}=\frac{\rho}{(2 \pi R_m T)^{\frac{D}{2}}} e^{- \frac{\vec \xi^2}{2 R_m T}} \Big[ 1 + \frac{\vec \xi \vec u}{R_m T} + \frac{(\vec \xi \vec u)^2}{2 (R_m T)^2} - \frac{\vec u^2}{2 R_m T} + \cdots \Big]$$

and evaluates the moments of the distribution, which are connected to the macroscopic properties,

$$\int \vec \xi^m f^{(eq)} d \vec \xi$$

with a numerical quadrature. From this quadrature with a particular number of integration points with their corresponding velocities emerge the weights $$w_\alpha$$ you are using as well as their velocities (velocity subsets, so called lattices with discrete directions $$\alpha$$). These weights must sum up to unity, meaning that you can't combine weights of different quadratures!

So far we could have used any discretisation with an arbitrary number of velocities but now comes the Chapman-Enskog expansion in which the macroscopic equations are recovered. For the incompressible algorithm this means the two moments

$$\sum_{\alpha} f_{\alpha}^{(0)} \begin{bmatrix} 1 \\ \vec c_{\alpha} \\ \end{bmatrix} = \begin{bmatrix} \rho \\ \rho \vec u \\ \end{bmatrix}$$

have to be considered.

In the process tensors of the form

$$E^{(n)} = \sum_{\alpha \neq 0} w_{\alpha} c_{\alpha i_1} c_{\alpha i_2} \cdots c_{\alpha i_n}$$

have to be evaluated. If your lattice has not a sufficient amount of velocities some of these moments will not evaluate correctly, certain terms in the resulting differential equation will not be present. This is be the case for the D2Q5 lattice. In other applications where only the sum of the populations has to be conserved such as scalar transport this is not problematic, and one might use discretisations with less velocities. For the full compressible Navier-Stokes equations on the other hand even the energy has to be conserved and an even higher amount of velocities is required: Single-speed lattices lack the isotropy in order to preserve the fully compressible conservation equations. Therefore the most common single-speed models apply a low Mach-number expansion of the Maxwell-Boltzmann equilibrium, an asymptotic approximation to the incompressible conservation equations, and are thus restricted to weakly compressible flow. Generally such schemes are used to accelerate the solution of incompressible fluid flows by allowing controlled compressibility.

First note that LBM doesn't actually work for in-compressible fluids unless you replace the density in your D2Q5 vectors with pressure (or energy etc.. something that doesn't relate to the number of particles per unit volume/area), as density is constant in in-compressible fluids (and I'll be from this point onward considering your LBM implementation in terms of pressure, not density). Second I would probably need to see the code, but it appears your equilibrium equation is wrong.

It looks like your equilibrium equation results in collisions with little to no perpendicular flow, and hence why it just looks like a straight shot across your grid.

Here is visual explanation:

imagine the following adjacent grid points, where blue represents non zero pressure.

On the next step, this is what your simulation looks like it is doing:

This is what it should be doing:

Collisions should result in perpendicular pressure flow to some degree. Think of billiard balls, if you hit one to the side, the ball doesn't just go straight, it hits off at an angle. Now you may be thinking "even with the billard ball situation, I can't get it to go exactly perpendicular to my collision!" and there are two factors to this:

• First LBM equilibrium collision equations happen at the particle level, and the equilibrium equation is supposed to not handle just one collision, but multiple collisions in a time-step. You in effect have multiple "billiard ball" collisions to account for (example):

• The quantized 4 directions + stationary (the Q5 in D2Q5) in LBM are not just supposed to be for the aligned particle movement to that quantized direction, but rather all particles angles more closely aligned with that angle than others. Here is what I'm talking about:

Each color should correspond to all particles flowing in that solid angle range, not just the particles that happen to flow aligned to the axis. Even in the single collision in the billard ball situation, the ball could still end up going in the "left general direction" despite not moving directly left, and would be advected leftwards.

I'm not sure what your equilibrium equation is, but I would suggest making it more isotropic, making sure that (if for example you use energy) you conserve energy.

• Very nice explanation and illustration. Thank you very much. Though, at least to me it seems that the equliibrium distribution function is independent of the lattice and number of neighbors used. I'm using the one stated in wikipedia (en.wikipedia.org/wiki/…). – Lenny Oct 18 '19 at 8:23
• @Lenny collision is $1/tf * (f_i^{eq} - f_i)$, the only thing you really can change in lattice boltzman is the $f_i^{eq}$, so it's pretty much guarnteed the issue is there. Additionally nowhere does it state that equilibrium is independent of lattice, it even shows two different types for D2Q9 and D3Q19, $w_i$ changes between those two, D2Q9 $w_0$ != D3Q19 $w_0$ ($4/9$ vs $1/3$), and so will also be different for D2Q5. – whn Oct 18 '19 at 15:55
• @Lenny Notice how $w_i$'s add up to 1, you can't just steal the values from D2Q9 for D2Q5 – whn Oct 18 '19 at 15:58
• you are correct. i did adjust the w's to 1/3 (center) and 1/6 for N,S,E,W – Lenny Oct 19 '19 at 17:37