What is the Lattice Boltzmann Method?

Fluid Mechanics: Standard Continuum Modeling

The equations governing fluid mechanics have been known since the 1840's with the development of the Navier-Stokes equations. The theoretical underpinning of fluid mechanics is surprisingly simple; one must simply evoke conservation of momentum on a small deformable volume, stating that the change in momentum of the element in time equals the applied force. In addition to conservation of momentum, if one assumes that the fluid is incompressible (like water at standard temperatures and pressures) and that the stresses in the fluid are proportional to the fluid's rate of deformation (strain rate), one obtains the incompressible Navier-Stokes equation

where u(x,y, t) is the velocity of the fluid at each point x, y at time t in space (it is a vector-field; it has a magnitude and direction at each point in space), ρ is the density of the fluid (it is a constant as the fluid is incompressible), P(x, y, t) is the pressure field imposed on the fluid, and ν is the kinematic viscosity which is essentially a measure of how hard it is to shear the fluid. For example, honey has a larger kinematic viscosity ν than water.

Many numerical schemes such as finite volume, finite element, and finite difference focus on discretizing this partial differential equation in space and solving for the time evolution of the velocity field u. These methods are generally intuitive but suffer from a couple of problems:

1. Depending on the technique used, the underlying numerical scheme may not obey all physical laws and symmetries (such as conservation of mass or invariance of the laws of physics under a rotation), leading to erroneous results.
2. It is difficult to choose how to discretize a complex geometry. A poor choice may give incorrect results.
3. To calculate spatial derivatives, one must look at the changes in u and P(x,y, t) in space. This requires many nonlocal-memory lookups to update a single grid point. Furthermore, it is often difficult to choose an appropriate scheme to calculate the derivative of a field; instabilities and other bizarre behavior may occur. Also note that the nonlinearities in the Navier-Stokes equations are particularly sensitive to ones choice of spatial discretization and differentiation, causing additional problems.

Mesoscopic View: The Lattice Boltzmann Technique

Introduction to the Technique

Motivated by the kinetic theory of fluids, a new technique to simulate fluids was developed in the 1970s: the Lattice Boltzmann technique.

At microscopic scales, fluids are composed of individual atoms bumping around and moving in a collective manner. When the behavior of enough atoms are averaged, we obtain the macroscopic Navier Stokes equations as discussed above. In between the microscopic and macroscopic scale lies the mesoscopic scale. The mesoscopic scale describes the probability density f(x,y, u, v t) of finding a particle at the position x, y, with horizontal velocity u and vertical velocity v at time t. Note that f explicitly depends on space and velocity now. At large enough scales, the fluctuations due to the probability density f are averaged out and we obtain the macroscopic Navier-Stokes equation. At smaller scales, we must track the motion of each atom individually which is a computationally infeasible task.

By prescribing physical rules that govern the evolution of f(x,y, u, v,t), it can be shown that, when coarse grained, the behavior of f reproduces the Navier Stokes equations. Deriving the equations that govern f are highly technical; we will not discuss them in detail here. The basic idea is that the evolution of f must obey all physical laws, such as conservation of momentum, mass, and energy, which constrains how f must evolve in time. For a good discussion of deriving the equations governing f, see the Lattice Boltzmann Scholarpedia page. Ultimately, the prescription that works is

This equation describes fields of particle f_i that move with a constant velocity c_i (i.e. the ith type of particle moves with velocity c_i). If the particles are on a lattice, they will hop one lattice site per time step in the c_i direction. The lattice we used was called "D2Q9," the 2 standing for two dimensions and 9 standing for 9 types of jumpers. The different directions of the jumpers can be seen in the figure below. Note that jumper 0 does not move.

In the equation governing the evolution of f, feq is a local equilibrium that depends on the current hydrodynamic variables, i.e. density ρ, pressure P, and velocity u and v. w is the relaxation frequency towards the local equilibrium field feq. It can be shown that the viscosity of a fluid ν is related to w by

where cs is the local speed of sound on the lattice. For the lattice we used (D2Q9), c_s=1/sqrt(3).

Coding Lattice Boltzmann

Now that we have discussed the basics of the Lattice Boltzmann method, let us write some pseudo code to help the reader imagine how a simulation would work. Note that our actual code is structured like the below (see it here).

# Begin by creating dimensionless parameters from input parameters...discussed more in "Simulation Verification"

...

# Initialize the grid based on the input resolution

...

#### System is setup; run lattice Boltzmann! ####

# Based on the initial conditions, initialize the hydrodynamic fields, like density and velocity

init_hydro()

# Based on the hydrodynamic fields, create the local equilibrium feq that the jumpers f will relax to

update_feq()

# Based on feq, create the initial population of jumpers. The population should be close in value to feq initially.

init_pop()

# Step the simulation forwards in time!

for cur_iteration in range(num_iterations):

# Moving the jumpers on the boundaries must be handled independently of the normal move step

move_bcs()

# Move all jumpers other than those on the boundary

move()

# Based on the new positions of the jumpers, update the hydrodynamic variables

update_hydro()

# Based on the new values of the hydrodynamic variables, update the equilibrium fields

update_feq()

# Relax the jumping fields towards the equilibrium fields. Depends on the value of w (the relaxation frequency)

collide_particles()

Although the discussion above may have been intimidating, it is actually relatively straightforwards to computationally implement the Lattice Boltzmann technique; all one has to do is write each of these functions for a desired lattice (in our case D2Q9). Note that the structure of this pseudo-code and our eventual final code was loosely based on Sauro Succi's Lattice Boltzmann code written in Fortran (the professor teaching AC274: Computational Fluid Dynamics).

Advantages of the Lattice Boltzmann Technique

The Lattice Boltzmann technique is clearly non-intuitive in comparison to a naive discretization of the Navier-Stokes equations. Why do people use it? There are several reasons:

1. All non-local spatial derivatives in the Navier Stokes Equation are replaced by a local memory lookup; this eliminates instabilities and other complexities when numerically solving the Navier Stokes equation.
2. The numerical scheme obeys all conservation laws (by definition), making simulations more reliable
3. It is easy to handle complex boundary conditions and geometries; no sophisticated meshing techniques are needed.
4. Since the Lattice Boltzmann method is a mesoscale technique, it is sometimes easier to incorporate complex physics missed by the Navier-Stokes equation.
5. Perhaps most importantly for the purposes of CS205, since most of the update step is independent of neighboring lattice sites, the lattice boltzmann technique is embarrassingly parallelizable. Very few memory lookups outside of a single lattice site are needed per iteration.