Milestone 03

Author

Andreas Greiner, Lars Pastewka

Collision operator

Learning goals

The student will…

  • …understand the BGK (relaxation time) approximation of the collision operator.
  • …be able to compute the equilibrium distribution function from the local density and velocity.
  • …combine streaming and collision into a full lattice Boltzmann time step.

Introduction

The Boltzmann transport equation (BTE) was given in the lecture:

\[\begin{equation} \frac{\partial f\left(\mathbf{r},\mathbf{v},t\right)}{\partial t}+\mathbf{v}\nabla_{\mathbf{r}} f\left(\mathbf{r},\mathbf{v},t\right) +\mathbf{a}\nabla_{\mathbf{v}} f\left(\mathbf{r},\mathbf{v},t\right)=C(f). \end{equation}\]

The l.h.s. of this equation is called the streaming term. The collision term on the r.h.s. represents the interactions between particles and usually is a complicated two particle scattering integral. We approximate this term by a relaxation time approximation. This approximation implies that the distribution function \(f\left(\mathbf{r},\mathbf{v},t\right)\) locally relaxes to an equilibrium distribution \(f^\mathrm{eq}\left(\mathbf{r},\mathbf{v},t\right)\). Viewing the streaming term as the total time derivative of the distribution function we can write the BTE as \[ \frac{d}{dt} f\left(\mathbf{r},\mathbf{v},t\right) = -\frac{f\left(\mathbf{r},\mathbf{v},t\right)- f^\mathrm{eq}\left(\mathbf{r},\mathbf{v},t\right)}{\tau} \tag{1}\] where \(\tau\) is the relaxation time constant.

What does this relaxation process look like? For the sake of clarity we imagine for a moment the distribution function to be independent of \(\mathbf{r}(t)\) and \(\mathbf{v}(t)\). Let us take the discrete form to read \[ f_i(\mathbf{r}+\mathbf{c}_i\Delta t,t+\Delta t) = f_i(\mathbf{r},t)+\omega\left(f_i^{eq}(\mathbf{r},t)-f_i(\mathbf{r},t)\right) \tag{2}\] The remaining question is how does \(f_i^{eq}(\mathbf{r})\) look like. Observe that the equilibrium is a local one, i.e. it depends on local variables. As given in the lecture these are the local density \(\rho(\mathbf{r})\) and the local average velocity \(\mathbf{u}(\mathbf{r})\). To this end you first have to calculate these two local quantities

\[ \rho(\mathbf{r})=\sum_i f_i \tag{3}\]

\[ \mathbf{u}(\mathbf{r})= \frac{1}{\rho(\mathbf{r})}\sum_i \mathbf{c}_i f_i(\mathbf{r}) \tag{4}\]

And the equilibrium distribution function is given as \[ f_i^\mathrm{eq}(\rho,\mathbf{u})=w_i\rho(\mathbf{r}) \left[1 +3\mathbf{c}_i\cdot\mathbf{u} +\frac{9}{2}\left(\mathbf{c}_i\cdot\mathbf{u}\right)^2 -\frac{3}{2}|\mathbf{u}|^2 \right] \tag{5}\] where \(w_i\) for the D2Q9 lattice is as follows: \[ \mathbf{w} = \left(\dfrac{4}{9}, \dfrac{1}{9}, \dfrac{1}{9}, \dfrac{1}{9}, \dfrac{1}{9}, \dfrac{1}{36}, \dfrac{1}{36}, \dfrac{1}{36}, \dfrac{1}{36}\right) \tag{6}\]

Figure 1: One full lattice Boltzmann time step: collision (relaxation towards the local equilibrium) and streaming (propagation to the neighboring nodes) alternate.

Tasks

So how do you run dynamic simulations in practice? As illustrated in Figure 1, every time step consists of a collision step followed by a streaming step:

  1. After the streaming step recalculate \(\rho\) from Equation 3.
  2. Calculate \(\mathbf{u}\) from Equation 4.
  3. Compute the new \(f_i(\mathbf{r},t)\) using Equation 2. What does \(\omega\) stand for? Discretize Equation 1 using a forward Euler step in time (with time step \(\Delta t\)) to arrive at Equation 2 and identify \(\omega = \Delta t / \tau\).

Let us run two tests:

  • Create a uniform density on your grid and set the density to a slightly higher value at the center. What happens? You should see the density bump spread outwards as (damped) sound waves that travel across the periodic domain; in the long run the density becomes uniform again.
  • Choose an initial distribution of \(\rho(\mathbf{r})\) and \(\mathbf{u}(\mathbf{r})\) at \(t=0\). Observe what happens dynamically as well as in the long time limit \(t\rightarrow\infty\). The dynamics depends on your choice, but for \(t\rightarrow\infty\) the collisions dissipate all gradients: the system relaxes towards a uniform density with a spatially constant velocity (zero in the co-moving frame), i.e. the global equilibrium.

Validation

Use the following checks (ideally as unit tests) to convince yourself that your collision operator is correct:

  • Mass conservation: the total mass \(\sum_{\mathbf{r},i} f_i(\mathbf{r})\) must be conserved to machine precision at every time step. The BGK operator conserves \(\rho\) locally because \(\sum_i f_i^\mathrm{eq} = \rho\).
  • Momentum conservation: without external forcing, the collision step must not change the local momentum \(\rho\mathbf{u} = \sum_i \mathbf{c}_i f_i\) at any node, since \(\sum_i \mathbf{c}_i f_i^\mathrm{eq} = \rho\mathbf{u}\). Check this before and after a collision step.
  • Equilibrium is a fixed point: if you initialize \(f_i = f_i^\mathrm{eq}(\rho, \mathbf{u})\) with uniform \(\rho\) and \(\mathbf{u}\), nothing must change — the collision term vanishes identically and streaming only shifts identical values. The fields \(\rho\) and \(\mathbf{u}\) must stay constant in time.
  • Relaxation: a small density perturbation on top of a uniform background must relax back towards uniform density, with the excess density radiating away as damped sound waves (your first test above).

Notes

  • Choose \(\omega\) with care, i.e. \(0<\omega<2\).
  • Choose \(0<\rho<1\) and \(|\mathbf{u}|<0.1\).
  • You need to set \(\rho\) and \(\mathbf{u}\) as initial conditions before the first streaming step.