Lattice Boltzmann Parameterization

February 2025

Fluid-Dynamics Lattice-Boltzmann Physics Simulation

The lattice Boltzmann method is defined in terms of ’lattice units’, where the microscopic velocities correspond to moving one lattice cell per iteration. How do these seemingly arbitrary units let us simulate real physical systems?

The trick is to use something called dynamic similarity to transform the simulation results into the desired physical system. To be able to do this, however, we must choose the correct parameters for our simulation (e.g. number of grid cells and collision timescale).

In this article I will show how to choose the parameters for a lattice Boltzmann simulation and how to map it onto the desired physical system.

Dynamic Similarity

Consider a physical system governed by the incompressible Navier-Stokes equation

$$\frac{\mathrm d \mathbf u}{\mathrm d t} = -\frac{\nabla P}{\rho} + \nu \nabla^2 \mathbf u + \mathbf g$$

This equation is “dimensionful”, with all terms having the units of acceleration; $LT^{-2}$. We have to choose a system of units with which to measure the various quantities. This choice obviously cannot change the behaviour we observe, so the NS equation must be valid for any system of units. This can be demonstrated easily.

Define characteristic length-, time- and mass-scales for the system; $\lambda$, $\tau$, $\beta$ respectively. For example for a flow around a cylinder $\lambda$ could be the cylinder size, and $\tau$ could be the time taken for the fluid to traverse one such distance at a typical velocity; $\tau = \lambda / V$. These characteristic scales are measured in our chosen system of units (SI, naturally).

Now rescale distances, times, and masses to be dimensionless by dividing by their respective characteristic scale

$$\tilde {\mathbf x} = \frac{\mathbf x}{\lambda} \qquad \tilde{t} = \frac{t}{\tau} \qquad \tilde{m} = \frac{m}{\beta}$$

This extends naturally to combined quantities, which we can simply divide by the right combination of $\lambda$, $\tau$ and $\beta$ to make them dimensionless. For example

$$\tilde{\mathbf u} = \mathbf u \frac{\tau}{\lambda} \qquad \tilde \rho = \rho \frac{\lambda^3}{\beta} \qquad \tilde P = P \frac{\lambda\tau^2}{\beta} \qquad \tilde\nu = \nu \frac{\tau}{\lambda^2} \qquad \text{etc}$$

Gradients introduce a factor reciprocal to the corresponding characteristic quantity. Given the definition of a derivative and the suggestive notation this is intuitive. But we can also see it clearly from the chain rule, for example

$$\partial_x = \partial_x \tilde x \, \partial_{\tilde x} = \frac{1}\lambda \partial_{\tilde x} $$

This is the opposite from the quantities themselves where we have $x = \lambda \tilde x$. This is effectively co- and contravariance of vectors under basis transformation.

We can invert these equations (e.g. $\mathbf x = \lambda \tilde{\mathbf x}$) and plug them into the Navier-Stokes equation. It should not be too surprising but every term just ends up with a prefactor $\beta\lambda^{-2}\tau^{-2}$, which corresponds to the dimensionality of a force density. We can cancel this prefactor to leave

$$\frac{\mathrm d \tilde{\mathbf u}}{\mathrm d \tilde t} = -\frac{\tilde \nabla \tilde P}{\tilde\rho} + \tilde \nu \tilde \nabla^2 \tilde{\mathbf u} + \tilde{\mathbf g}$$

where $\tilde\nabla = \partial_{\tilde{\mathbf x}}$ and similarly for $\tilde\nabla^2$. This equation is exactly the same as the original dimensionful equation but we’ve replaced every occurrence of a variable with its dimensionless counterpart.

The thing to notice here is that the characteristic scales are independent of the units we measure them with. If we had chosen a different system of units the numerical value of both $\mathbf x$ and $\lambda$ (for example) would be different by the same factor so we would get the same values for $\tilde{\mathbf x}$.

So far so unsurprising; units definitely shouldn’t affect things. However, this also tells us that different physical systems that have the same ratios between dimensional quantities and their characteristic scales will have the same dimensionless behaviour. They are said to be ’dynamically similar’ .

Rescaling all quantities in accordance with their dimensionality means that dimensionless combinations, such as the Reynold’s number , are maintained in the dimensionless system. For example

$$\mathrm{Re} = \frac{u L}{\nu} = \frac{\tilde u \frac\lambda \tau \, \tilde L\lambda}{\tilde\nu\,\frac{\lambda^2}{\tau} } = \frac{\tilde u\tilde L}{\tilde\nu} = \widetilde{\mathrm{Re}}$$

The flip-side of this is that if two systems have the same values for all relevent dimensionless combinations we know they are dynamically similar. We can convert between the systems using the ratio of their characteristic scales; for example distances can be converted with $\lambda / \lambda^\prime$:

$$\mathbf x = \lambda \tilde{\mathbf x} = \frac{\lambda}{\lambda^\prime} \mathbf x ^ \prime$$

Mapping to Physical Systems

Why is dynamic similarity useful? Because it lets us simulate things in the arbitrary units of the simulation then convert into our desired physical system. We only need to ensure that the Reynold’s number in lattice units matches the desired Reynold’s number in physical units.

In the following, any quantities that are measured in simulation units will be indicated with a prime; e.g. $\tau^\prime.$ Most of the time we will be working with simulation units so this is a bit clunky. However every other notation I tried was either worse or made things unclear at certain points. Oh well, let’s stick with the prime notation.

Unit Systems

Lets say that our simulation grid has a cell size of $\Delta x$ in physical units. By definition each cell corresponds to a size of 1 in lattice units. Similarly for the time-scale; one iteration of the simulation corresponds to a time interval of $\Delta t$ in physical units. The mass-scale, by contrast, is not directly constrained by the simulation method. Usually we are not interested in the value of $\rho$ directly but only through its effect on the kinematic viscosity $\nu$. For simplicity therefore we can assume a rest density of 1 in the simulation, which gives us a conversion factor of $\rho_0$. Thus we have

Units	Physical	Simulation	Conversion
Length	$\quad \Delta x \quad$	$\Delta x^\prime = 1$	$\Delta x/\Delta x^\prime = \Delta x$
Time	$\quad \Delta t \quad$	$\Delta t^\prime = 1$	$\Delta t/\Delta t^\prime = \Delta t$
Mass	$\quad \rho_0 \quad$	$\rho_0^\prime = 1$	$\rho_0 / \rho_0^\prime = \rho_0$

The table shows the conversion factor from simulation units into physical units. The conversion factor from physical into simulation units is given by the reciprocal of the appropriate factor. As above, the conversion factors for combined quantities can be simply derived from their dimensionality.

Recall that the speed-of-sound in lattice units is given by ${c_s^\prime} = 1/\sqrt{3}$. This is a dimensional quantity like any other, so requires conversion when it appears in an equation in physical units. In physical units it becomes

$$c_s = c_s^\prime \frac{\Delta x}{\Delta t} = \frac{1}{\sqrt{3}}\frac{\Delta x}{\Delta t}$$

The factor $c_s$ comes from the conversion from “raw units” to “lattice units”, but here we are directly converting from lattice units to physical units so we do not need to worry about this extra conversion from raw units. We only need to worry about the conversion of $c_s$ where it appears in equations relating lattice units to physical units, for example for the kinematic velocity (see below.)

Setting Viscosity

The direct conversion of the units is pretty trivial, but to obtain dynamic similarity we need to correctly specify the kinematic viscosity of the simulation.

We know the viscosity is related to the collision term in the Boltzmann equation. The most basic LBM collision operator is the Bhatnagar–Gross–Krook (BGK) operator, which exponentially relaxes the distribution towards equilibrium. In physical units this reads

$$\Omega = \frac{1}{\tau} (f^\mathrm{eq} - f)$$

where $\tau$ is a relaxation time, not the characteristic timescale of the system. $\Omega$ is an instantaneous rate; the Boltzmann equation reads $\mathrm{d} f/\mathrm{d} t = \Omega$. Therefore one time-step of our simulation corresponds to a change of $\Delta f = \Omega\,\Delta t$. This is equivalent to saying that in simulation units we have

$$\Delta f = \Omega^\prime = \frac{1}{\tau^\prime} (f^\mathrm{eq} - f)$$

with $\tau^\prime = \tau / \Delta t$ being the usual conversion of the relaxation time into simulation units. All of which is to say our code for collision will look something like

omega = dt / tau
f[i] += omega * (f_eq[i] - f[i])

But how do we choose the value of $\tau^\prime$?

A Chapman-Enskog type analysis of the time-discretized Lattice Boltzmann equation with BGK operator shows that the kinematic viscosity is given by

$$\nu = {c_s}^2 \Big(\tau - \frac{\Delta t}{2}\Big) \quad\Longleftrightarrow\quad \nu^\prime = \frac{1}{3}\Big(\tau^\prime - \frac{1}{2}\Big)$$

thus we can set the viscosity through our choice of $\tau^\prime$, or we can set $\tau^\prime$ through our choice of viscosity

$$\tau^\prime = 3\nu^\prime + \frac{1}{2} = 3\frac{\nu \Delta t}{\Delta x^2} + \frac{1}{2}$$

Choosing Parameters

Dynamic stability means we only care about the value of the Reynold’s number, not the specific values of $u^\prime$, $L^\prime$ and $\nu^\prime$. This gives us more freedom in how we parameterize our lattice Boltzmann simulation; we can change the simulation parameters so long as $\mathrm{Re}$ is maintained, then map back into physical units using the conversion factors.

Since we will be matching the Reynold’s number of the simulation to the desired physical system we have $\mathrm{Re} = \mathrm{Re}^\prime$, so we will just write $\mathrm{Re}$ henceforth.

We have three parameters we can tune;

Viscosity $\nu^\prime$ (via $\tau^\prime$).
Velocities $u^\prime$.
Simulation grid size $L^\prime$

Note that for the last we have $L^\prime = L / \Delta x$ being the number of grid cells for one characteristic length of our system. So increasing the simulation grid size corresponds to increasing resolution / decreasing $\Delta x$.

Parameter Bounds

Error Terms

Unfortunately we are not completely free to choose the parameters. Increasing the value of $\tau^\prime$ not only increases the viscosity, it also increases the magnitude of the error terms arising from the discretization of the Boltzmann equation. Therefore we must choose values for $\Delta x$, $\Delta t$ and $\nu$ which give us a sensible value for $\tau^\prime$ while still satisfying our needs. If we say we want a maximum value $\tau^\prime_\mathrm{max}$ then our parameterization must satisfy the following inequality:

$$\tau^\prime \leq \tau^\prime_\mathrm{max} \quad\implies\quad \nu^\prime \leq \frac 1 3 \Big(\tau^\prime_\mathrm{max} - \frac 1 2\Big)$$

CFL

Just like any numerical method, LBM has certain constraints that must be satisfied for the method to be stable. The CFL condition tells us the following inequality must be obeyed

$$M = \frac{u}{c_s} = \sqrt 3 \, u^\prime \leq M_\mathrm{max} \quad\implies\quad u^\prime \leq \frac{M_\mathrm{max}}{\sqrt{3}}$$

for constant $M_\mathrm{max}.$ That is, the flow must stay below some critical Mach number .

Note that because the lattice Boltzmann method is not Galilean invariant (the lattice & associated velocities pick out a preferred reference frame) the velocity here is absolute, not relative as it would be for other methods like SPH .

We want $M <\!\!< 1$ anyway because transient velocities approaching the speed of sound would lead to large compression and we are interested in (approximately) incompressible flows. Further the NS equation emerges from the LB equation only in the limit of small $\mathrm{M}$.

Collision

We want $\tau^\prime$ to have a value approaching $1/2$ to minimize error, but unfortunately we cannot set it arbitrarily close to $1/2$ because the simulations rapidly become unstable. In Niu, 2004 the authors investigate the stability of the BGK collision operator and empirically find the following constraint for the simulation to remain stable

$$u^\prime \leq 8\Big(\tau^\prime - \frac{1}{2}\Big)\qquad\text{for }\tau^\prime \leq 0.55$$

By plugging in the value for $\nu^\prime$ and rearranging we get an lower bound for $\nu^\prime$ in terms of $u^\prime$:

$$\nu^\prime \geq \frac{u^\prime}{24}$$

Note, this bound is specific to the BGK collision operator. Other collision operators such as ’two relaxation-time’ or ‘multi relaxation-time’ can have much lower values of $\tau^\prime$ and remain stable.

We may also just want to impose a minimum value for $\tau^\prime$ giving another bound

$$\tau^\prime \geq \tau^\prime_\mathrm{min} \quad\implies\quad \nu^\prime \geq \frac 1 3 \Big(\tau^\prime_\mathrm{min} - \frac 1 2\Big)$$

Fixed Resolution

It seems logical when setting up our simulations that we have a particular resolution in mind we would like to simulate. This means the value of $\Delta x$ is fixed, and thus $L^\prime$ is set. Since we are also targeting a particular value for the Reynold’s number this leaves us $u^\prime$ and $\nu^\prime$ to choose. Our choice must obey

$$\frac{\nu^\prime}{u^\prime} = \frac{L^\prime}{\mathrm{Re}}$$

This is easy to visualize in the $(u^\prime, \nu^\prime)$ plane as simply a line with gradient $\frac{L^\prime}{\mathrm{Re}}.$ Any point along this line will give us the desired Reynold’s number but obviously we have to be cognizant of the bounds discussed above. Valid parameterizations are given by the intersection of this line with the region bounded by the various inequalities.

Stable Region, Fixed L — Stable Region, fixed $L$

This plot makes it quite clear the choices we can make in our parameterization. Firstly, we can see that we cannot stably simulate systems where $L^\prime / \mathrm{Re} < 1/24$ for any value of $u^\prime.$ This means if we want to simulate high Reynold’s number flows we will need to increase our resolution to increase the gradient of this line.

Because $u^\prime$ and $\nu^\prime$ are linearly related, we can convert the limit on $u^\prime$ into a limit for $\nu^\prime$. Thus the various inequalities define an interval $\nu^\prime\in[\nu^\prime_\mathrm{min}, \nu^\prime_\mathrm{max}]$, whose bounds are

$$\begin{align} \nu^\prime_\mathrm{max} &= \min\Big\{\frac{1}{3}(\tau^\prime_\mathrm{max} - \frac{1}{2}), \frac{L^\prime}{\mathrm{Re}} \frac{M_\mathrm{max}}{\sqrt 3}\Big\} \\ \nu^\prime_\mathrm{min} &= \frac{1}{3}(\tau^\prime_\mathrm{min}-\frac{1}{2}) \end{align}$$

We can choose a value for either $\nu^\prime$ or $u^\prime$, and this immediately fixes the other. Once we have made our choice we can convert $\nu^\prime$ into $\tau^\prime.$

We could, of course, have worked directly with $\tau^\prime$ rather than $\nu^\prime$ since they are linearly related. It doesn’t change much and I decided to stick with $\nu^\prime$ since that’s what appears in the equation for the Reynold’s number.

As we increase $u^\prime$ it takes fewer time-steps for the fluid to traverse the characteristic length-scale $L^\prime$ and so the speed (in the computational sense) of our simulation will increase. Conversely if we choose a parameterization with a smaller $u^\prime$ then we must take more steps for the same characteristic time and our simulation will take longer. The benefit, however, is it will be more accurate. Note that this change in computational performance is linear in $u^\prime$.

If you you don’t particularly care about the accuracy of your simulation, and you just want it to run as fast as possible while staying stable, then we can forget about the bound $\tau_\mathrm{max}$. In that case you can set

$$\tau^\prime = \frac{\sqrt{3} M_\mathrm{max} L}{\mathrm{Re}} - \frac{1}{2}$$

which comes from converting the Mach number bound on $u^\prime$ first to be in terms of $\nu^\prime$ then to $\tau^\prime$

Fixed Viscosity

Another possible way to parameterize the simulation is to start by choosing a fixed value for $\tau^\prime$, and thus $\nu^\prime$. This gives us the two parameters $u^\prime$ and $L^\prime$ which must obey

$$u^\prime L^\prime = \nu^\prime \mathrm{Re}$$

Similar to above, we can visualize this in the $(L^\prime, u^\prime)$ plane and we get a hyperbola of valid choices. There will be a hard upper limit to $L^\prime$ simply because as we increase the resolution we use more memory, and at some point we will run out. For $u^\prime$ we have the same upper limits as above. It is a given that the choice of $\tau^\prime$ lies in the interval $[\tau^\prime_\mathrm{min}, \tau^\prime_\mathrm{max}]$ and so these bounds do not appear on this plot.

Stable Region, Fixed tau — Stable Region, fixed $\tau$

This plot shows us that increasing either $\nu^\prime$ or $\mathrm{Re}$ too much will push us outside the region of stable parameterizations (towards the top right of the plot). This is mainly due to the memory limit, however. Because if it weren’t for that we could always increase $L^\prime$ to bring ourselves below the Mach limit.

As above moving along the line of valid parameterizations changes the computational burden. However, unlike above things are more dire here. As we increase $L^\prime$ the number of lattice cells grows as ${L^\prime}^d$ and the required speed drops as $u^\prime \propto^{-1} {L^\prime}$ so we need to take more time-steps to simulate the same characteristic time. We can show that these combined effects mean the computation burden grows as ${L^\prime}^{d+2}!$

Proof

The characteristic time in simulation units is in effect the number of iterations per characteristic time in physical units. It is given by

$$T^\prime = \frac{L^\prime}{u^\prime} = \frac{{L^\prime}^2}{\nu^\prime\mathrm{Re}}$$

Where the second equality follows simply from rearranging the equation for $\mathrm{Re}$ and substituting for $u.$ The total number of cells in our simulation is $N \propto {L^\prime}^d$, and the total work we need to do is the number of cells multiplied by the number of iterations, so

$$NT^\prime \propto {L^\prime}^{d+2}$$

To fix the values we choose a specific value for either $L^\prime$ or $u^\prime$ and then infer the other. Choosing $L^\prime$ is equivalent to fixing the resolution and letting the number of time-steps per output vary, whereas choosing $u^\prime$ is equivalent to fixing the number of iterations per output and letting the resolution vary.

Recipe

A typical ‘recipe’ for our LBM simulations looks like this:

Decide on the physical system we are simulating.
- Dimensions, boundary-conditions, viscosity, etc.
Calculate the characteristic velocity, length- and time-scales for the system.
Use these to calculate the Reynold’s number.
Choose suitable simulation parameters in lattice units that match the Reynold’s number.
Run simulation in its ’native’ unit system.
Convert back into physical units by appropriately multiplying dimensional quantities by the correct powers of $\Delta x$ and $\Delta t$.

Below is a code snippet which I use to choose the parameters for my LBM simulations, where I am always assuming a fixed resolution (number of grid cells) and optionally fixed $\tau^\prime.$ If you specify $\tau^\prime$ it will just calculate the appropriate value for $u^\prime$ and verify the stability bounds. Alternatively if you leave $\tau^\prime$ unspecified it will choose the largest permissible value for maximum simulation speed. The velocity $u^\prime$ is always inferred.

Code Snippet

def calc_lbm_params(
    # params
    Re: float,
    L: float,
    tau: float | None = None,
    # bounds
    M_max: float = 0.1,
    tau_max: float = 1.0,
    tau_min: float = 0.5,
) -> tuple[float, float]:

    if tau is None:
        tau = sqrt(3) * M_max * L / Re + 0.5    # Maximum permissible by Mach condition
        tau = min(max(tau, tau_min), tau_max)   # Hard bounds on tau

    nu = (tau - 0.5) / 3
    u = Re * nu / L

    eps = 1e-8
    assert u - eps <= (M_max / np.sqrt(3)), "Mach condition violated."
    assert u - eps <= 24 * nu, "BGK stability condition violated."

    return (tau, u)

To convert back to physical units we need to know the conversion factors. Usually we will be setting $\Delta x$ explicitly and so is known, but if not it can be calculated from the characteristic length-scales:

$$\Delta x = \frac{L}{L^\prime}$$

Similarly the time-step can be calculated $\Delta t$ by using the characteristic velocities:

$$\Delta t = \frac{u^\prime}{u} \Delta x$$

Example - Plane Poiseuille

Let’s test out the parameterization and unit conversion with the example of plane Poiseuille flow. The characteristic length scale is the width of the pipe $L=W$ and the final max velocity is given by

$$u = \frac{1}{8}\frac{g}{\nu} W^2$$

This gives a Reynold’s number of

$$\mathrm{Re} = \frac{1}{8}\frac{g W^3}{\nu^2}$$

Some vaguely realistic numbers for these for water in a pipe would be $W=1 \text{ cm}$, $\nu=10^{-6}\text{ m}^2/\text{s}$ and $g=0.08\text{ m}/\text{s}^2$ (a pressure drop of $8\text{ Pa/m}$). These parameters correspond to a Reynold’s number of 1000.

Let’s try a few different LBM parameterizations to check our formulae all work

Param Set	$\tau^\prime$	$L^\prime$	$u^\prime$	$T^\prime$
#1	0.55	100	0.1666	600
#2	0.51	100	0.0333	3,000
#3	0.55	200	0.0833	2,400
#4	0.51	200	0.0166	12,000

“Simulation Results — Simulation Results

The figure shows the results for these simulations after running for the same length of physical time; 175 sec. There are clearly some numerical errors arising (for example at boundaries), but apart from that the results all look similar and in line with the analytical solution. This shows our parameterization is working. Phew!

Plane Poiseuille Flow

Dipolar Forces

Units	Physical	Simulation	Conversion
Length	\(\quad \Delta x \quad\)	\(\Delta x^\prime = 1\)	\(\Delta x/\Delta x^\prime = \Delta x\)
Time	\(\quad \Delta t \quad\)	\(\Delta t^\prime = 1\)	\(\Delta t/\Delta t^\prime = \Delta t\)
Mass	\(\quad \rho_0 \quad\)	\(\rho_0^\prime = 1\)	\(\rho_0 / \rho_0^\prime = \rho_0\)