Lattice Boltzmann
The Boltzmann transport equation (BTE) was introduced by Ludwig Boltzmann in the context of kinetic gas theory (Boltzmann 1995) and is a statistical model for the transport of molecular constituents during flow (Cercignani 1988). The Lattice Boltzmann method (LBM) is a numerical scheme based on a discretized version of the BTE, introduced by McNamara and Zanetti (1988). LBM models have been used for the last three decades to study the dynamics of fluids in many different applications, from multiphase flow (Gunstensen et al. 1991; Grunau et al. 1993), porous media (Aharonov and Rothman 1993; Gunstensen and Rothman 1993) to microfluidics (Zhang 2011). A thorough review of its applications in fluid dynamics and beyond can be found in Succi (2001).
Probability density
Given the positions \(\mathbf r_i\) and velocities \(\mathbf v_i\) of a huge number \(N\) of molecules, we would have to solve a huge number of equations of motion \(\mathbf {\dot r}_i\) and \(\mathbf {\dot v}_i\). Indeed, this is what one does in Molecular Dynamics simulations.
Instead of following the trajectories of all \(N\) particles in phase space, i.e. the space spanned by the \(2D \times N\) coordinates of positions and velocities, we take the average over an infinitesimal “volume” in phase space. Thus \(f(\mathbf v,\mathbf r,t)\) turns out to be a probability density for finding a molecule with velocity \(\mathbf{v}\) at position \(\mathbf{r}\) at time \(t\). The 0th, 1st and 2nd moments of this probability density define the mass density \(\rho(\mathbf{r}, t)\), the momentum density \(\rho(\mathbf{r}, t)\mathbf{u}(\mathbf{r}, t)\) and the temperature \(T(\mathbf{r}, t)\) respectively in \(D\)-dimensional space, \[ \begin{split} \rho(\mathbf{r}, t) &= m \int \mathrm{d}^Dv\, f(\mathbf{v},\mathbf{r},t) \\ \rho(\mathbf{r}, t)\mathbf{u}(\mathbf{r}, t) &=m \int \mathrm{d}^Dv\, \mathbf{v} f(\mathbf{v},\mathbf{r},t) \\ k_B T(\mathbf{r}, t) &= \frac{m^2}{D \rho(\mathbf{r}, t)} \int \mathrm{d}^Dv\, \left[\mathbf{v}-\mathbf{u}(\mathbf{r}, t)\right]^2 f(\mathbf{v},\mathbf{r},t), \end{split} \tag{1}\] where we have assumed a monoatomic system with molecules of mass \(m\). Here \(\mathbf{u}(\mathbf{r}, t)\) is the average velocity at position \(\mathbf{r}\) and \(k_B\) is the Boltzmann constant.
The Boltzmann transport equation
The BTE describes the temporal evolution of the probability density \(f(\mathbf{v},\mathbf{r},t)\), i.e. the total time-rate of change of this probability distribution, \(\mathrm{d}f/\mathrm{d}t\). We know that for large times \(t\) it must relax towards statistical equilibrium, given by the Maxwell velocity distribution function (Huang 1987), \[ f^\text{eq}(\mathbf{v}; \rho, \mathbf{u}, T) = \frac{\rho}{m} \left(\frac{m}{2\pi k_B T}\right)^{D/2} \exp\left\{-\frac{m(\mathbf{v}-\mathbf{u})^2}{2 k_B T}\right\}. \tag{2}\] The Maxwell distribution is a Gaussian centered on the local mean velocity \(\mathbf{u}\) whose width grows as \(\sqrt{T}\), as illustrated in Figure 1.
A common approximation is to assume relaxation of \(f\) towards \(f^\text{eq}\) with a single characteristic time \(\tau\), as suggested by Bhatnagar et al. (1954), \[ \frac{\mathrm{d}f(\mathbf{v}, \mathbf{r}, t)}{\mathrm{d}t} = -\frac{f(\mathbf{v},\mathbf{r},t)-f^\text{eq}(\mathbf{v}; \rho(\mathbf{r}, t), \mathbf{u}(\mathbf{r}, t), T(\mathbf{r}, t))}{\tau}. \tag{3}\] Equation 3 is the BGK-Boltzmann equation. Note that the total derivative of \(f\) is \[ \begin{split} \frac{\mathrm{d}f}{\mathrm{d}t} &=\mathbf{\dot v} \frac{\partial f(\mathbf{v},\mathbf{r},t)}{\partial \mathbf{v}}+ \mathbf{\dot r} \frac{\partial f(\mathbf{v},\mathbf{r},t)}{\partial \mathbf{r}} +\frac{\partial f(\mathbf{v},\mathbf{r},t)}{\partial t} \\ &=\frac{\mathbf{F}(\mathbf{r})}{m} \frac{\partial f(\mathbf{v},\mathbf{r},t)}{\partial \mathbf{v}}+ \mathbf{v}\frac{\partial f(\mathbf{v},\mathbf{r},t)}{\partial \mathbf{r}} +\frac{\partial f(\mathbf{v},\mathbf{r},t)}{\partial t}, \end{split} \] where \(\mathbf{F}(\mathbf{r})\) is an external force acting on the molecules. The left-hand side of Equation 3 is therefore a transport operator: it describes the streaming of the probability density through phase space (both in real space and in velocity space). The right-hand side represents the collision operator, which drives the distribution towards local equilibrium.
We will remain within the context of this single (BGK) relaxation time approximation, but modern developments of the method aim at introducing multiple relaxation times, for example by relaxing the cumulants of \(f\) individually (Geier et al. 2006; Geier et al. 2015).
Streaming
The left-hand side of the BTE is a transport operator. Think of the probability density as a quantity that is advected through phase space. In real space, the density is transported at velocity \(\mathbf{v}\); in velocity space, it is transported at the acceleration \(\mathbf{F}/m\). Together, these two processes constitute the streaming of phase space density along the characteristics of the BTE.
In the absence of external forces (\(\mathbf{F} = 0\)), the streaming reduces to free transport: \[ \frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla_{\mathbf{r}} f = \text{collision term}. \] This states that a packet of probability density at \((\mathbf{r}, \mathbf{v})\) moves to \((\mathbf{r} + \mathbf{v}\Delta t, \mathbf{v})\) after a time step \(\Delta t\).
Collision
The right-hand side of the BTE represents the effect of molecular collisions. The full Boltzmann collision integral involves a six-fold integration over all possible scattering partners and outcomes: \[\left(\frac{\partial f}{\partial t}\right)_\text{coll}= \int\int\int w(\mathbf{v}_1,\mathbf{v}_2,\mathbf{v}_3,\mathbf{v}) \left\{ f(\mathbf{r},\mathbf{v}_1,t) f(\mathbf{r},\mathbf{v}_2,t)- f(\mathbf{r},\mathbf{v}_3,t) f(\mathbf{r},\mathbf{v},t) \right\} \mathrm{d}\mathbf{v}_1\mathrm{d}\mathbf{v}_2\mathrm{d}\mathbf{v}_3 \] where \(w(\mathbf{v}_1,\mathbf{v}_2,\mathbf{v}_3,\mathbf{v})\) is the scattering kernel – the probability that two particles with incoming velocities \(\mathbf{v}_1\) and \(\mathbf{v}_2\) scatter into outgoing velocities \(\mathbf{v}_3\) and \(\mathbf{v}\). The two terms represent gain (scattering into velocity \(\mathbf{v}\)) and loss (scattering out of velocity \(\mathbf{v}\)).
This collision integral is nonlinear in \(f\) and extremely expensive to evaluate. In the BGK approximation (Equation 3), it is replaced by a simple relaxation towards the local equilibrium \(f^\text{eq}\) with a single time scale \(\tau\): \[\left(\frac{\partial f(\mathbf{r},\mathbf{v})}{\partial t}\right)_\text{coll}=-\frac{f(\mathbf{r},\mathbf{v})-f^\text{eq}(\mathbf{r},\mathbf{v})}{\tau}\] This means that the system is pushed towards the local equilibrium distribution \(f^\text{eq}(\mathbf{r},\mathbf{v})\) at a rate \(1/\tau\) (Figure 2). Although this is a drastic simplification, it preserves the conservation laws (mass, momentum, energy) and reproduces the correct Navier-Stokes equations in the macroscopic limit.
The Maxwell distribution from maximum entropy
Equation 2 states the equilibrium distribution as a result; we now sketch its derivation from the principle of maximum entropy. The fact that the equilibrium is a Gaussian, rather than any other shape, is not a postulate of kinetic theory – it is a consequence of choosing the distribution that has the largest entropy compatible with the macroscopic constraints (mass, momentum and energy). The same principle, applied to the discrete velocity set, will give us the discrete equilibrium used in the LBM (Equation 7).
The H-functional
Boltzmann’s H-functional is the negative of the entropy carried by the distribution \(f\), up to a factor of \(k_B\), \[ H[f] = \int \mathrm{d}^D v\, f(\mathbf{v}) \ln f(\mathbf{v}). \tag{4}\] The thermodynamic entropy density is \(s = -k_B H[f]\), and Boltzmann’s H-theorem states that under the full Boltzmann collision integral \(H\) is non-increasing in time. The equilibrium distribution \(f^\text{eq}\) is therefore the distribution that minimizes \(H\) (maximizes entropy) at fixed values of the conserved local moments.
Constrained variation
We minimize \(H[f]\) subject to the three local moment constraints from Equation 1 (rewritten as conditions on \(f\) itself): \[ \int \mathrm{d}^D v\, f = \frac{\rho}{m}, \qquad \int \mathrm{d}^D v\, \mathbf{v}\, f = \frac{\rho}{m}\, \mathbf{u}, \qquad \int \mathrm{d}^D v\, (\mathbf{v}-\mathbf{u})^2 f = \frac{D\, \rho\, k_B T}{m^2}. \] Introducing Lagrange multipliers \(\alpha\), \(\boldsymbol{\beta}\) and \(\gamma\) for each constraint and setting \(\delta \mathcal{L}/\delta f = 0\) in the Lagrangian \[ \mathcal{L}[f] = \int \mathrm{d}^D v\, \Big\{ f \ln f + \alpha\, f + \boldsymbol{\beta}\cdot\mathbf{v}\, f + \gamma (\mathbf{v}-\mathbf{u})^2 f \Big\} \] yields \[ \ln f^\text{eq} + 1 + \alpha + \boldsymbol{\beta}\cdot \mathbf{v} + \gamma (\mathbf{v}-\mathbf{u})^2 = 0, \] i.e. \[ f^\text{eq}(\mathbf{v}) = e^{-(1+\alpha)}\, \exp\!\Big[ -\boldsymbol{\beta}\cdot\mathbf{v} - \gamma(\mathbf{v}-\mathbf{u})^2 \Big]. \tag{5}\] Completing the square shows that the linear term \(\boldsymbol{\beta}\cdot\mathbf{v}\) simply shifts the centre of the Gaussian from \(\mathbf{u}\) to \(\mathbf{u} - \boldsymbol{\beta}/(2\gamma)\). The momentum constraint \(\int \mathbf{v}\, f^\text{eq}\, \mathrm{d}^D v = (\rho/m)\mathbf{u}\) therefore forces \(\boldsymbol{\beta} = \mathbf{0}\).
Identifying the multipliers
The remaining two multipliers are fixed by the normalization and variance constraints. Using the Gaussian integrals \[ \int \mathrm{d}^D v\, e^{-\gamma (\mathbf{v}-\mathbf{u})^2} = \left(\frac{\pi}{\gamma}\right)^{D/2}, \qquad \int \mathrm{d}^D v\, (\mathbf{v}-\mathbf{u})^2\, e^{-\gamma (\mathbf{v}-\mathbf{u})^2} = \frac{D}{2\gamma} \left(\frac{\pi}{\gamma}\right)^{D/2}, \] the normalization (mass) constraint gives \[ e^{-(1+\alpha)} = \frac{\rho}{m}\left(\frac{\gamma}{\pi}\right)^{D/2} \] and the variance (energy) constraint gives \[ \gamma = \frac{m}{2 k_B T}. \] Substituting both back into Equation 5 recovers exactly the Maxwell distribution Equation 2.
The identification \(\gamma = m/(2 k_B T)\) is what fixes the meaning of temperature in kinetic theory: it is the Lagrange multiplier dual to the energy, just as the chemical potential is dual to mass and the drift velocity is dual to momentum. This is why the equilibrium has the universal form \(f^\text{eq} \propto \exp[-(\text{energy})/k_B T]\) familiar from statistical mechanics.
The Lattice Boltzmann method in two dimensions
Discretization
The BTE is discretized in space, velocity and time. Discretizing Equation 3 on a regular (square) lattice in space is straightforward. However, we also require a suitable discretization of velocity space and the time step. Both are chosen such that the distance between interpolation points in velocity space multiplied by the time step equals the distance between points on the spatial lattice; in other words, a molecule traveling on the spatial lattice moves exactly between neighboring lattice points during one time step.
The particular realization of the velocities used by us is shown in Figure 3. The velocity set contains 9 directions. Each direction is assigned the index shown in Figure 3. Direction \(0\) describes the population of molecules at rest. This specific discretization in two-dimensional space (\(D=2\)) with nine directions is commonly denoted by D2Q9.
The full velocity set is listed in the following table. In lattice units (\(\Delta x = \Delta t = 1\), see below), the velocity vectors are simply the integer offsets to the neighboring sites. The last column gives the opposite direction \(\bar{\imath}\), defined by \(\mathbf{c}_{\bar{\imath}} = -\mathbf{c}_i\), which will be needed for the bounce-back boundary conditions.
| \(i\) | \(\mathbf{c}_i\) | \(w_i\) | \(\bar{\imath}\) (opposite) |
|---|---|---|---|
| 0 | \((0,0)\) | \(4/9\) | 0 |
| 1 | \((1,0)\) | \(1/9\) | 3 |
| 2 | \((0,1)\) | \(1/9\) | 4 |
| 3 | \((-1,0)\) | \(1/9\) | 1 |
| 4 | \((0,-1)\) | \(1/9\) | 2 |
| 5 | \((1,1)\) | \(1/36\) | 7 |
| 6 | \((-1,1)\) | \(1/36\) | 8 |
| 7 | \((-1,-1)\) | \(1/36\) | 5 |
| 8 | \((1,-1)\) | \(1/36\) | 6 |
The discrete velocities are given by the distances to the eight neighbors divided by the time step (the rest population has zero velocity). We have nine velocity vectors \(\mathbf{c}_i\) with \(i=0,\dots,8\) for each direction \(i\). The probability distribution \(f(\mathbf{r},\mathbf{v})\) is represented by nine discrete values \(f_i(\mathbf{x}_j, t)\) where \(\mathbf{x}_j\) is a discrete lattice point and \(i\) denotes the direction. The moments of the probability density (Equation 1) become \[\begin{align} \rho(\mathbf{x}_j, t) &= \sum_i f_i(\mathbf{x}_j, t) \\ \mathbf{u}(\mathbf{x}_j, t) &= \frac{1}{ \rho(\mathbf{x}_j, t)} \sum_i \mathbf{c}_i f_i(\mathbf{x}_j, t), \end{align}\] where \(\rho(\mathbf{x}_j, t)\) is now the number density, i.e. we assume unit molecular mass.
The discrete Boltzmann equation
To obtain the lattice Boltzmann equation (LBE), we discretize the BGK-Boltzmann equation (Equation 3) in time using a forward Euler step with time step \(\Delta t\). Along the characteristic \(\mathbf{r}(t+\Delta t) = \mathbf{r}(t) + \mathbf{c}_i \Delta t\) for direction \(i\), this gives: \[ \frac{f_i(\mathbf{x}_j + \mathbf{c}_i \Delta t, t+\Delta t) - f_i(\mathbf{x}_j, t)}{\Delta t} = -\frac{1}{\tau}\left[f_i(\mathbf{x}_j, t) - f_i^\text{eq}(\mathbf{x}_j, t)\right] \] Multiplying both sides by \(\Delta t\) and defining the dimensionless relaxation parameter \(\omega = \Delta t / \tau\), we obtain the LBE: \[ f_i(\mathbf{x}_j+\mathbf{c}_i\Delta t,t+\Delta t)=f_i(\mathbf{x}_j,t) -\omega\left[f_i(\mathbf{x}_j,t)-f_i^\text{eq}(\mathbf{x}_j,t)\right] \tag{6}\] The left-hand side represents streaming (transport to the neighboring lattice site) and the right-hand side contains the collision operator (relaxation towards equilibrium). The external forces \(\mathbf{F}(\mathbf{r})\) are set to zero. In lattice units, where \(\Delta t = 1\) (see below), the definition \(\omega = \Delta t/\tau\) reduces to \(\omega = 1/\tau\); combined with the viscosity relation Equation 12 this gives the convention used in the project, \(\tau = \nu/c_s^2 + 1/2\).
Equilibrium distribution
The expression for the discrete equilibrium distribution function is (Mohamad 2011; Wolf-Gladrow 2000) \[ f_i^\text{eq}(\mathbf{x}_j, t) = w_i \rho(\mathbf{x}_j, t) \left[ 1 + 3 \mathbf{c}_i \cdot \mathbf{u}(\mathbf{x}_j, t) + \frac{9}{2} \left(\mathbf{c}_i \cdot \mathbf{u}(\mathbf{x}_j, t)\right)^2 - \frac{3}{2} u^2(\mathbf{x}_j, t) \right] \tag{7}\] with weights \[ w_i = \begin{cases} 4/9 & \text{for } i=0 \\ 1/9 & \text{for } i \in \{1,2,3,4\} \\ 1/36 & \text{for } i \in \{5,6,7,8\} \end{cases} \] Like Equation 2 for the continuous distribution function, Equation 7 is obtained from the maximum entropy principle of Section 3, but now in the discrete nine-dimensional space of populations (Mohamad 2011; Wolf-Gladrow 2000). The derivation also explains why the numerical prefactors \(3\), \(9/2\) and \(3/2\) – which appear arbitrary in Equation 7 – are entirely fixed by the lattice geometry through the speed of sound \(c_s^2 = 1/3\).
Discrete maximum entropy problem
We maximize the discrete H-functional \[ H_d = \sum_i f_i \ln(f_i / w_i) \tag{8}\] subject to mass and momentum conservation only, \[ \sum_i f_i = \rho, \qquad \sum_i \mathbf{c}_i f_i = \rho\, \mathbf{u}. \] There is no energy constraint because the standard isothermal LBM does not transport temperature; the second moment of \(f_i\) is instead fixed by the lattice structure itself, as we will see below. The reference weights \(w_i\) appear inside the logarithm so that, at zero velocity, the entropy is maximized by \(f_i = w_i \rho\) – the lattice-symmetric rest state – rather than by the artificial uniform distribution \(f_i = \rho/9\) that would otherwise treat rest, cardinal and diagonal populations on the same footing.
Introducing Lagrange multipliers \(\alpha\) and \(\boldsymbol{\beta}\) for the two constraints and setting the variation of the Lagrangian to zero, \[ \delta \sum_i \Big[ f_i \ln(f_i/w_i) + \alpha\, f_i + \boldsymbol{\beta}\cdot \mathbf{c}_i\, f_i \Big] = 0, \] yields \[ \ln(f_i^\text{eq}/w_i) + 1 + \alpha + \boldsymbol{\beta}\cdot \mathbf{c}_i = 0, \] i.e. \[ f_i^\text{eq} = w_i\, \lambda\, e^{-\boldsymbol{\beta}\cdot\mathbf{c}_i}, \tag{9}\] with \(\lambda = e^{-1-\alpha}\). This is the exact maximum-entropy solution – the discrete analogue of Equation 5, but with the quadratic energy term in \(\mathbf{v}-\mathbf{u}\) replaced by a linear-in-\(\mathbf{c}_i\) momentum term because we have not constrained the energy.
Low-Mach expansion
Unlike the continuous Gaussian, the exponential in Equation 9 cannot be summed in closed form against only nine velocities. We therefore expand it in the Mach number \(|\mathbf{u}|/c_s\), which is small in the regime where the LBM reproduces Navier-Stokes (\(|\mathbf{u}| \ll c_s\)). To second order in \(\mathbf{u}\), the equilibrium has the polynomial form \[ f_i^\text{eq} = w_i\, \rho \left[ A + B\,(\mathbf{c}_i\cdot\mathbf{u}) + C\,(\mathbf{c}_i\cdot\mathbf{u})^2 + D\, u^2 \right], \tag{10}\] with four unknown coefficients \(A\), \(B\), \(C\) and \(D\). They are fixed by combining the conservation constraints with the moment identities of the D2Q9 velocity set, \[\begin{align} \sum_i w_i &= 1, \\ \sum_i w_i\, c_{i\alpha} &= 0, \\ \sum_i w_i\, c_{i\alpha} c_{i\beta} &= c_s^2\, \delta_{\alpha\beta}, \\ \sum_i w_i\, c_{i\alpha} c_{i\beta} c_{i\gamma} &= 0, \\ \sum_i w_i\, c_{i\alpha} c_{i\beta} c_{i\gamma} c_{i\delta} &= c_s^4\, \big(\delta_{\alpha\beta}\delta_{\gamma\delta} + \delta_{\alpha\gamma}\delta_{\beta\delta} + \delta_{\alpha\delta}\delta_{\beta\gamma}\big), \end{align}\] with \(c_s^2 = 1/3\). The weights \(w_i\) are in fact chosen precisely so that these moment identities hold up to fourth order, which is the deeper reason why the D2Q9 set has its specific values \(4/9\), \(1/9\) and \(1/36\).
Substituting Equation 10 into the mass constraint gives \[ A + (C c_s^2 + D)\, u^2 = 1 \quad\Longrightarrow\quad A = 1, \quad D = -C c_s^2, \] and into the momentum constraint gives \[ B\, c_s^2\, u_\beta = u_\beta \quad\Longrightarrow\quad B = \frac{1}{c_s^2}. \] This leaves one undetermined coefficient \(C\) (with \(D\) slaved to \(C\)). It is fixed by requiring the second moment of \(f_i^\text{eq}\) to reproduce the Euler-level momentum-flux tensor of the Maxwell distribution, \[ \sum_i c_{i\alpha} c_{i\beta} f_i^\text{eq} = \rho\, c_s^2\, \delta_{\alpha\beta} + \rho\, u_\alpha u_\beta. \tag{11}\] This is the condition that the LBM reproduce the inviscid pressure and convective momentum flux, and it is the bridge between the discrete kinetic model and the Euler/Navier-Stokes equations through the Chapman-Enskog analysis. Substituting Equation 10 into Equation 11 and using the fourth-order lattice moment yields \[ C = \frac{1}{2\, c_s^4}, \qquad D = -\frac{1}{2\, c_s^2}. \]
Specialising to D2Q9 (\(c_s^2 = 1/3\)), the four coefficients are \(A = 1\), \(B = 3\), \(C = 9/2\) and \(D = -3/2\), which recovers Equation 7 exactly. The numerical prefactors are therefore not arbitrary: they are entirely determined by the lattice speed of sound, which itself follows from the geometry of the D2Q9 velocity set and the values of the lattice weights.
The collide-stream algorithm
An algorithmic implementation of Equation 6 is conveniently split into a collision step and a streaming step (Figure 4). The algorithm for a single time step is:
Collision: at every lattice site \(\mathbf{x}_j\), compute the density \(\rho\) and velocity \(\mathbf{u}\) from the current \(f_i\), evaluate the equilibrium \(f_i^\text{eq}\), and relax towards it: \[ f_i^*(\mathbf{x}_j, t) = f_i(\mathbf{x}_j, t) - \omega\left[f_i(\mathbf{x}_j, t) - f_i^\text{eq}(\mathbf{x}_j, t)\right] \] where \(f_i^*\) denotes the post-collision distribution.
Streaming: propagate each post-collision population to the neighboring lattice site in the corresponding direction: \[ f_i(\mathbf{x}_j + \mathbf{c}_i \Delta t, t + \Delta t) = f_i^*(\mathbf{x}_j, t) \] Figure 5 illustrates this propagation on the lattice.
Boundary conditions: apply boundary conditions (periodic, bounce-back, etc.) to correct the populations at domain boundaries. See the boundary conditions lecture.
Repeat from step 1.
The collision step is purely local – each lattice site can be updated independently, making it trivially parallel. The streaming step involves nearest-neighbor data access, which requires communication between adjacent subdomains in a distributed-memory parallelization (see milestone 06).
Stability
The relaxation parameter \(\omega\) must satisfy \(0 < \omega < 2\) for the LBE to be stable. Recall from Equation 12 that the viscosity \(\nu = \frac{1}{3}(1/\omega - 1/2)\) decreases monotonically with increasing \(\omega\):
- \(\omega \to 0\): the collision has almost no effect; \(f_i\) relaxes only very slowly towards equilibrium and the viscosity diverges, \(\nu \to \infty\). This is the most diffusive (over-damped) limit.
- \(\omega = 1\): \(f_i\) relaxes fully to equilibrium in a single time step (\(f_i^* = f_i^\text{eq}\)), giving \(\nu = 1/6\) in lattice units.
- \(\omega \to 2\): the collision overshoots equilibrium (over-relaxation) and the viscosity vanishes, \(\nu \to 0\). This is the stability limit.
In practice, values of \(\omega\) close to 2 produce low viscosity but also numerical instabilities, particularly at higher Reynolds numbers.
Lattice units
In the LBM, all quantities are expressed in lattice units: the lattice spacing \(\Delta x = 1\), the time step \(\Delta t = 1\), and the molecular mass \(m = 1\). The lattice speed of sound is \(c_s = 1/\sqrt{3}\), giving \(c_s^2 = 1/3\).
To convert between lattice units and physical units, one must choose three reference scales – typically a reference length \(L_\text{ref}\) (e.g. the domain size in lattice units), a reference velocity \(u_\text{ref}\) (e.g. the lid velocity), and the density \(\rho_\text{ref}\). Physical quantities are then obtained by:
- Physical length: \(x_\text{phys} = x_\text{lattice} \cdot \Delta x_\text{phys}\) where \(\Delta x_\text{phys} = L_\text{phys} / L_\text{ref}\)
- Physical time: \(t_\text{phys} = t_\text{lattice} \cdot \Delta t_\text{phys}\) where \(\Delta t_\text{phys} = \Delta x_\text{phys} \cdot u_\text{lattice} / u_\text{phys}\)
- Physical viscosity: \(\nu_\text{phys} = \nu_\text{lattice} \cdot \Delta x_\text{phys}^2 / \Delta t_\text{phys}\)
The dimensionless Reynolds number \(\text{Re} = u L / \nu\) is the same in both lattice and physical units, providing a consistency check.