Extending the LBM to 2D and 3D
August 2024
In this article I’ll show the extension to two and three dimensions of the Lattice Boltzmann Method, as detailed in the previous article on the Lattice Boltzmann Method in 1D.
There are really only two things which have to change;
- We need a multidimensional definition of the approximate equilibrium distribution.
- We need a suitable Gauss-Hermite velocity set which covers 2D and 3D.
For the first item, we will need a definition of the Hermite polynomials for arbitrary dimensions, so that is our starting point. Then for the second we will see how to combine D1Q3 with itself to get the higher dimensional velocity sets.
Hermite Polynomials
The Hermite series expansion was key to our derivation of the 1D Lattice Boltzmann Method. Wikipedia has a nice introduction to the 1D Hermite polynomials . To extend our derivation to 2 and 3 dimensions we need a definition of the Hermite polynomials in higher dimensions.
One definition of the (probabilist’s) Hermite polynomials in 1D is
$$H_n(q) = (-1)^n \frac{1}{w(q)} \frac{d^n}{dq^n} w(q) $$where
$$w(q) = \frac{1}{\sqrt{2\pi}}\exp\Big(-\frac{q^2}{2}\Big)$$is the pdf of a standard Gaussian.
An obvious extension of this to multiple dimensions - $\mathbf q\in \mathbb R^d$ - is to consider the $n$-th order partial derivative instead
$$\partial^n_x = \underbrace{\partial_{x}\partial_{x}\dots\partial_{x}}_{n\text{ times}}$$However now, since we have more than 1 dimension, at each position in the chain of partial derivatives we can choose any dimension to differentiate wrt to. This gives us a rank $n$ tensor with each dimension being of size $d$.
Denote the label of each dimension as $\alpha \in \{1,\dots d\}$. In 2D this would give us $\alpha=1$ for the $x$ dimension and $\alpha=2$ for the $y$ dimension. The $\alpha$-th dimension of the velocity is then $(\mathbf q)_{\alpha} = q_\alpha$. For our $n$-th order tensor we have $n$ such labels $\mathbf a = \alpha_1,\dots \alpha_j\dots \alpha_n$. This notation looks cumbersome (and potentially opaque) but it lets us write the fully generic case of $n$-th order partial differentiation as
$$\prod_{j=1}^n\frac{\partial}{\partial q_{\alpha_j}}$$… okay, so definitely a little opaque. An example will help elucidate things; say $n=3$ and $d=2$ and we have $\alpha_1=1$, $\alpha_2=2$ and $\alpha_3=1$, this gives us (being deliberately very explicit)
$$\prod_{j=1}^3\frac{\partial}{\partial {q_{\alpha_j}}} = \frac{\partial}{\partial {q_{\alpha_1}}}\frac{\partial}{\partial {q_{\alpha_2}}}\frac{\partial}{\partial {q_{\alpha_3}}} = \frac{\partial}{\partial {q_x}}\frac{\partial}{\partial {q_y}}\frac{\partial}{\partial {q_x}}=\partial_x\partial_y\partial_x$$
Note that I’ve shown the examples here using the (micro) velocity $q$ since that is ultimately what we want to expand. But it has another advantage!
Often expositions on vector calculus topics use the position, denoted by $\mathbf x$, as their domain of interest. This can lead to confusion since $x_\alpha$ is the position in the $\alpha$-th dimension, it has nothing to do with the $x$-dimension, which is just the case $\alpha=1$. This is just an unfortunate use of $x$ both to represent the position, and as a label for the first dimension. In fact it is even worse since the actual value of the first component of position is often denoted $x$, so we have ($\mathbf x)_x = x_x = x$. The reason for the conflation is obvious, but it still means we need to be a bit careful.
Discussions written denoting the position with $\mathbf r$ side-step the first of these issues but not necessarily the second.
where inside $w$ we have $q^2 = |\mathbf q|^2 = \sum_{\alpha=1}^d q_\alpha^2$.
Examples - 2D
In 2D we get the following for the first few orders;
$$H_0(\mathbf q) = 1$$ $$\begin{align} \big(\mathbf H_1(\mathbf q)\big)_{x} = -\frac{1}{w(\mathbf q)}\partial_xw(\mathbf q) = q_x \qquad \big(\mathbf H_1(\mathbf q)\big)_{y} = -\frac{1}{w(\mathbf q)} \partial_y w(\mathbf q) = q_y \end{align}$$which is typically denoted in vector notation as
$$\mathbf H_1(\mathbf q) = -\frac{1}{w(\mathbf q)}\nabla w = \left[\begin{matrix}q_x \\ q_y\end{matrix}\right]$$then second order gives a tensor
$$\begin{align} \big(\mathbf H_2(\mathbf q)\big)_{xx} &= \frac{1}{w(\mathbf q)}\partial_{xx} w(\mathbf q) = q_x^2 - 1 \qquad &\big(\mathbf H_2(\mathbf q)\big)_{xy} = \frac{1}{w(\mathbf q)} \partial_{xy} w(\mathbf q) = q_xq_y \notag \\ \big(\mathbf H_2(\mathbf q)\big)_{yx} &= \frac{1}{w(\mathbf q)}\partial_{yx} w(\mathbf q) = q_xq_y\qquad &\big(\mathbf H_2(\mathbf q)\big)_{yy} = \frac{1}{w(\mathbf q)} \partial_{yy} w(\mathbf q) =q_y^2 - 1 \notag \end{align}$$which is related to the Hessian matrix of $w$ (I’m not going to write that step because too many $H$’s.) In component notation we can write this as $\big(\mathbf H_2(\mathbf q)\big)_{ij} = q_i q_j - \delta_{ij}$ or in vector notation
$$\mathbf H_2(\mathbf q) = \left[\begin{matrix} q_x^2 - 1 & q_xq_y \\ q_xq_y & q_y^2 - 1 \end{matrix}\right] = (\mathbf q\otimes \mathbf q - \mathbf I)$$Orthogonality
Like their 1-dimensional counterparts, the $d$-dimensional Hermite polynomials are orthogonal to each other wrt $w(\mathbf q)$
$$\int_{\mathbb R^d} \big(\mathbf H_n(\mathbf q)\big)_\mathbf a\,\big(\mathbf H_m(\mathbf q)\big)_\mathbf b\,w(\mathbf q)\,d\mathbf q = \delta_{nm}\delta_{\mathbf a\mathbf b}$$where $\delta_{\mathbf a \mathbf b}$ is 1 if the set of indices $\mathbf a$ is a permutation of the set of indices $\mathbf b$, and 0 otherwise. For example $\mathbf a = (x, y, x)$ and $\mathbf b=(x, x, y)$ would give $\delta_{\mathbf a \mathbf b} = 1$. This is because partial derivatives commute.
Hermite Series
Using the higher dimensional Hermite polynomials we can, again, build a series expansion for certain functions in $\mathbb R^d$; namely those which satisfy
$$\int_{\mathbb R^d} w(\mathbf q) f(\mathbf q)^2 < \infty$$The $n$-th order coefficient is now an $n$ dimensional tensor but otherwise the expansion looks much as in 1D;
$$f(\mathbf q) = w(\mathbf q)\sum_{n=0}^\infty \frac{1}{n!}\mathbf a_n\cdot\mathbf H_n(\mathbf q) $$with the coefficients given by
$$\mathbf a_n = \int f(\mathbf q) \mathbf H_n(\mathbf q)\,d\mathbf q$$Equilibrium Distribution
Approximation & Discretization
In 1D we used the Hermite series to approximate our equilibrium distribution. The approximation doesn’t change greatly in higher dimensions. We already noted in the previous article that the equilibrium distribution is the Maxwell distribution
$$\mathbf q^\mathrm{eq} \sim \mathcal{N}\left(\mathbf u, \,\frac{k_B T}{m}\mathbf I\right)$$where we defined $\theta = k_B T / m$ for convenience. The pdf of such a multivariate Gaussian is
$$f^\mathrm{eq}(\mathbf q) = \frac{\rho}{\sqrt{(2\pi\,\theta)^d}}\exp\left(-\frac{|\mathbf q - \mathbf u|^2}{2\theta}\right) = \frac{\rho}{\sqrt\theta} w\left(\frac{\mathbf q - \mathbf u}{\sqrt\theta}\right)$$If we evaluate the integrals above we get the following coefficients
$$\begin{align} a_0 &= \rho \\ (\mathbf a_1)_\alpha &= \rho u_\alpha \\ (\mathbf a_2)_{\alpha\beta} &= \rho(u_\alpha u_\beta + (\theta - 1)\delta_{\alpha\beta}) \end{align}$$In vector notation these would be written
$$\begin{align} a_0 &= \rho \\ \mathbf a_1 &= \rho \mathbf u \\ \mathbf a_2 &= \rho(\mathbf u\otimes\mathbf u + (\theta - 1)\mathbf I) \end{align}$$We can see the coefficients of the 1D expansion are just the first component of each of these coefficients. Using these coefficients gives us our approximate equilibrium distribution
$$\begin{align} f^\mathrm{eq}(\mathbf q) &\approx w(\mathbf q)\rho\left(1 + \mathbf u\cdot\mathbf q + \frac{1}{2}\big(\mathbf u\otimes\mathbf u + (\theta -1)\mathbf I\big)\cdot\big(\mathbf q\otimes \mathbf q - \mathbf I\big)\right) \\ &= w(\mathbf q)\rho\left(1 + u_\alpha q_\alpha + \frac{1}{2}(u_\alpha u_\beta + (\theta-1)\delta_{\alpha\beta})(q_\alpha q_\beta - \delta_{\alpha\beta})\right) \end{align}$$with Einstein summation convention used for the latter.
The next step was to discretize the velocity using Gauss-Hermite quadrature, so let’s do that again. It gives us
$$f^\mathbf{eq}_i = w_i\rho\left(1 + \mathbf u\cdot \mathbf q_i + \frac{1}{2}\big(\mathbf u\otimes\mathbf u + (\theta - 1)\mathbf I\big)\cdot(\mathbf q_i\otimes \mathbf q_i - \mathbf I)\right)$$Velocity Sets
To actually use this method we need a suitable $d$-dimensional velocity set, analogous to D1Q3. Fortunately we can simply take a tensor product of D1Q3 $d$ times to get the most commonly used velocity sets for lattice Boltzmann simulations.
D2Q9
“Taking the tensor product” of our velocity set D1Q3 means that we form the 9 vectors whose 2 components are each combination of the 3 velocities in D1Q3. This is most natural if you think about the new set of velocities and weights living on a grid indexed by $(\alpha, \beta)$ where each dimension is an index into D1Q3
$$\mathbf c^\prime_{\alpha\beta}= (c_{\alpha}, c_{\beta})$$where $\mathbf c^\prime$ is the D2Q9 velocity and $c$ is the D1Q3 velocity. The weights are determined by the (renormalized) product of the D1Q3 weights
$$w^\prime_{\alpha\beta} = \frac{w_{\alpha} w_{\beta}}{W},\qquad W = \sum_{\alpha,\beta} w_\alpha w_\beta$$This give us the following velocities and weights for D2Q9;
| $i\vphantom{\frac{1}{6}}$ | $0$ | $1$-$4$ | $5$-$8$ |
|---|---|---|---|
| $w_i$ | $\frac{4}{9}$ | $\frac 1 9$ | $\frac 1 {36}$ |
| $c_i\vphantom{\frac{1}{6}}$ | $(0, 0)$ | $(\pm1, 0)$ $(0, \pm 1)$ |
$(\pm1, \pm1)$ |
The exact enumeration of our new velocity set does not matter, so long as each unique $(\alpha,\beta)$ pair maps to a unique index in our higher dimensional velocity set.
D3Q27
To get a 3-dimensional velocity set we can similarly take the tensor product of D1Q3 with our new D2Q9 velocity set. The procedure is exactly the same and yields
| $i\vphantom{\frac{1}{6}}$ | $0$ | $1$-$6$ | $7$-$18$ | $19$-$26$ |
|---|---|---|---|---|
| $w_i$ | $\frac{8}{27}$ | $\frac 2 {27}$ | $\frac 1 {54}$ | $\frac{1}{216}$ |
| $c_i\vphantom{\frac{1}{6}}$ | $(0, 0, 0)$ | $(\pm1, 0, 0)$ $(0, \pm1, 0)$ $(0, 0, \pm1)$ |
$(\pm1, \pm1, 0)$ $(\pm1, 0, \pm1)$ $(0, \pm1, \pm1)$ |
$(\pm1, \pm1, \pm1)$ |
Isothermal assumption
In 1-dimension we made the isothermal assumption because without it our problem was over-constrained and the collision operator would have had no effect. People also make the assumption that $\theta=1$ in 2- and 3-dimensional LBM. I am not sure if this is necessary in these cases, but it does reduce the number of macroscopic moments we need to calculate and simplifies the form of $f_i^\mathrm{eq}$. Perhaps it also has implications for the stability of the method. In effect we are breaking the conservation the second moment (i.e. energy) so does this lead to dissipation? Something to ponder.
Simulation
The actual steps involved for our simulation do not change with the number of dimensions, but for completeness I’ll restate them here:
- Advect / Stream - For point $(\mathbf x, \mathbf q)$ shift the value of $f$ to $(\mathbf x+\mathbf q, \mathbf q)$.
- Calculate - Calculate the new macroscopic moments at each point as simple finite sums.
- Infer - Infer the local equilibrium distribution using our $f_i^\mathrm{eq}$ in lattice-units.
- Collide - Perform the BGK collision operation.
The resulting simulation can be seen in the next article Lattice Boltzmann in 2D - Interactive.