So as showed by Frisch et al. (a), the 2D Euler equation $$v_{t}+ v\cdot \nabla v=\mu \Delta v$$ can be derived by the Hexagonal-placed automaton (for low velocity).

I am curious about the existence of similar derivations for the following equations:

  1. the passive scalar with white noise and known velocity field v: $$\theta_{t}+ v\cdot \nabla \theta=\mu \Delta \theta+\phi$$

  2. the vorticity form with $\omega=\nabla \times v$: $$\omega_{t}+ v\cdot \nabla \omega=\mu \Delta \omega$$

  3. The vorticity form with white noise: $$\omega_{t}+ v\cdot \nabla \omega=\mu \Delta \omega+\phi$$

Did I miss any obvious references that contain these derivations? Thank you.

Remark: Once we have a derivation for $\omega_{t}+ v\cdot \nabla \omega=\mu \Delta \omega$, then the force can be represented by adding a momentum at each site (see Rothman (b) for details). So (3) should follow from (2).

Update: So nluigi resolved the passive scalar problem (1) and similar approach will guide (2). Specifically, for (2) I found this paper from 2011 "A lattice Boltzmann model for the eddy–stream equations in two-dimensional incompressible flows" (c), which seems to tackle this problem. But I still would like to see any other more detailed treatments or extensions because I don't want to rely on a single paper.


(a) Uriel Frisch, Brosl Hasslacher and Yves Pomeau. Lattice-gas automata for the Navier-Stokes equation. Physical Review Letters 56, 14 (1986) 1505.

(b) Daniel H. Rothman and Stiphane Zaleski. Lattice-gas cellular automata: simple models of complex hydrodynamics. Vol. 5. Cambridge University Press, 2004.

(c)Yan, Bo, et al. "A lattice Boltzmann model for the eddy–stream equations in two-dimensional incompressible flows." Applied Mathematical Modelling 35.5 (2011): 2358-2365.



  1. I hope you are ok with a LB derivation (your title suggests you are).
  2. I am going to assume you are familiar with the fundamentals of LB: discretization, BGK collision, weights $w_i$, velocity set $\boldsymbol{e}_i$, symmetry conditions on moments, etc.
  3. I will be deriving the convection-diffusion equation using a distribution function $g_i$ assuming the velocity field follows from a hydrodynamic equation often described by a distribution function $f_i$. This means that velocity field does not follow from $g_i$ and that any moments of $g_i$ higher than zeroth are undefined.
  4. I try to keep this short because including all details will cost a lot time. However if anything is unclear let me know and i will try to explain by making edits.


The LB equation describes the evolution of a distribution function $g_i$ according to the following special discretization of the Boltzmann equation: $$\frac{g_i(\boldsymbol{x}+\boldsymbol{e}_i\Delta t, t+\Delta t)-g_i\left(\boldsymbol{x},t\right)}{\Delta t}=\Omega_i+w_i\phi\qquad\Omega_i=-\frac{g_i-g_i^{(eq)}}{\tau}$$

where $\Omega_i$ is the BGK collision operator, $\phi$ is a source/sink term, $\tau$ is a relaxation constant and $g_i^{(eq)}$ is an equilibrium distribution which we define as: $$g_i^{(eq)}=w_i\rho\theta\left[1+\frac{\boldsymbol{e}_i\cdot\boldsymbol{v}}{c_s^2}+\frac{1}{2}\frac{\left(\boldsymbol{e}_i\boldsymbol{e}_i-c_s^2\boldsymbol{I}\right):\boldsymbol{v}\boldsymbol{v}}{c_s^4}\right]$$

in terms of macroscopic quantities; density $\rho$, scalar variable $\theta$ and velocity $\boldsymbol{v}$. Taking moments of $g_i^{eq}$ yields:

$$\sum_i g_i^{(eq)}=\rho\theta \quad \sum_i \boldsymbol{e}_i g_i^{(eq)}=\rho\theta\boldsymbol{v} \quad \sum_i \left(\boldsymbol{e}_i\boldsymbol{e}_i-c_s^2\boldsymbol{I}\right)g_i^{(eq)}=\rho\theta\boldsymbol{v}\boldsymbol{v}$$

Multi-scale analysis

Next we will apply a multi-scale analysis to determine the continuum equations from the discretized equation. We first convert the discretized equation to a continuous LB equation using a Taylor expansion up to second-order terms:

$$g_i(\boldsymbol{x}+\boldsymbol{e}_i\Delta t, t+\Delta t)\approx g_i\left(\boldsymbol{x},t\right)+Dg_i\left(\boldsymbol{x},t\right)\Delta t+\frac{1}{2}D^2g_i\left(\boldsymbol{x},t\right)\Delta t^2$$ $$Dg_i+\frac{1}{2}D^2g_i\Delta t\approx\Omega_i+w_i\phi$$

where $D=\partial_t+\boldsymbol{e}_i\cdot\boldsymbol{\nabla}$ is the material derivative.

Next, a smallness parameter $\epsilon$ is introduced which signifies the degree to which we are away from equilibrium:

$$g_i=g_i^{(eq)}+\epsilon g_i^{(neq)}=g_i^{(eq)}+\epsilon g_i^{(1)}+\epsilon^2 g_i^{(2)}$$

Here $(neq)$ refers to 'non-equilibrium', i.e. all terms as a result of the system being out of equilibrium are contained in here. These non-equilibrium terms are further seperated into a convective scale $g_i^{(1)}$ and a diffusive scale $g_i^{(2)}$.

Note: A system at low velocities guarentees we are always near equilibrium and is required for stability in BGK. Sometimes $\epsilon$ is refered to the Knudsen number and/or made equivalent to $\Delta t$, either way it should be small.

The same applies to the temporal and spatial scale: $$\partial_t=\epsilon\partial_t^{(1)}+\epsilon^2\partial_t^{(2)}\qquad\boldsymbol{\nabla}=\epsilon\boldsymbol{\nabla}^{(1)}$$

The assumption made here is that the diffusive time scale is much shorter than the convective time scale. On the other hand, convective and diffusive transport are assumed to occur on the same spatial scale, hence only first order terms are retained. The material derivative becomes: $$D=\epsilon D^{(1)}+\epsilon^2 D^{(2)} \qquad D^{(1)}=\partial_t^{(1)}+\boldsymbol{e}_i\cdot\boldsymbol{\nabla}^{(1)}\qquad D^{(2)}=\partial_t^{(2)}$$

Likewise for simplicity the collision operator and source term are expanded as follows: $$\Omega_i=\epsilon\Omega_i^{(1)}+\epsilon^2\Omega_i^{(2)}\qquad \Omega_i^{(k)}=-\frac{g_i^{(k)}}{\tau}\qquad \phi=\epsilon^2 \phi^{(2)}$$

We now substitute into the continous LB equation where only terms up to $O\left(\epsilon^2\right)$ are retained: $$\left(\epsilon D^{(1)} + \epsilon^2 D^{(2)}\right)\left(g_i^{(eq)}+\epsilon g_i^{(1)}\right) + \frac{1}{2}\epsilon^2D^{(1)2}g_i^{(eq)}\Delta t\approx\epsilon\Omega_i^{(1)}+\epsilon^2\Omega_i^{(2)}+\epsilon^2 w_i\phi^{(2)}$$

Collecting equivalent orders in $\epsilon$ yields, respectively: $$O\left(\epsilon\right):\quad D^{(1)}g_i^{(eq)} = \Omega_i^{(1)}$$ $$O\left(\epsilon^2\right):\quad D^{(2)}g_i^{(eq)} + D^{(1)}g_i^{(1)} + \frac{1}{2}D^{(1)2}g_i^{(eq)}\Delta t= \Omega_i^{(2)}+w_i\phi^{(2)}$$

The latter equation may be 'simplified' using the former: $$D^{(2)}g_i^{(eq)} + D^{(1)}\left[g_i^{(1)} + \frac{\Delta t}{2}\Omega_i^{(1)}\right] = \Omega_i^{(2)}+w_i\phi^{(2)}$$

Next, we take the zeroth moment of the $O\left(\epsilon\right)$ equation: $$\partial_t^{(1)}\sum_i g_i^{(eq)} + \nabla^{(1)}\cdot\sum_i\boldsymbol{e}_i g_i^{(eq)} = \sum_i\Omega_i^{(1)}$$ $$\partial_t^{(1)}\rho\theta + \boldsymbol{\nabla}^{(1)}\cdot\rho\theta\boldsymbol{v} = 0$$

Now, we take the zeroth moment of the $O\left(\epsilon\right)$ equation: $$\partial_t^{(2)}\sum_i g_i^{(eq)} + \partial_t^{(1)}\sum_i\left[g_i^{(1)} + \frac{\Delta t}{2}\Omega_i^{(1)}\right] + \boldsymbol{\nabla}^{(1)}\cdot\sum_i\boldsymbol{e}_i\left[g_i^{(1)} + \frac{\Delta t}{2}\Omega_i^{(1)}\right]= \sum_i\Omega_i^{(2)}+\sum_i w_i\phi^{(2)}$$ $$\partial_t^{(2)}\rho\theta + \boldsymbol{\nabla}^{(1)}\cdot \boldsymbol{j}^{(1)} = \phi^{(2)}$$

where the diffusive flux is identified as:

$$\boldsymbol{j}^{(1)} = \sum_i\boldsymbol{e}_i\left[g_i^{(1)} + \frac{\Delta t}{2}\Omega_i^{(1)}\right]$$

This is further simplified to:

$$\begin{aligned} \boldsymbol{j}^{(1)} &= -\Delta t\left(\frac{\tau}{\Delta t}-\frac{1}{2}\right)\sum_{i}\boldsymbol{e}_{i}\Omega_{i}^{(1)} \\ &= -\Delta t\left(\frac{\tau}{\Delta t}-\frac{1}{2}\right)\sum_{i}\boldsymbol{e}_{i}D^{(1)}g_i^{(eq)} \\ &= -\Delta t\left(\frac{\tau}{\Delta t}-\frac{1}{2}\right)\left[\partial_{t}^{(1)}\sum_{i}\boldsymbol{e}_{i}g_{i}^{(eq)}+\boldsymbol{\nabla}^{(1)}\cdot\sum_{i}\boldsymbol{e}_{i}\boldsymbol{e}_{i}g_{i}^{(eq)}\right] \\ &= -\Delta t\left(\frac{\tau}{\Delta t}-\frac{1}{2}\right)\left[\partial_{t}^{(1)}\rho\theta\boldsymbol{v}+\boldsymbol{\nabla}^{(1)}\cdot\left(\rho\theta\boldsymbol{v}\boldsymbol{v}+\rho\theta c_{s}^{2}\boldsymbol{I}\right)\right] \\ &= -\Delta t\left(\frac{\tau}{\Delta t}-\frac{1}{2}\right)\left[\rho\theta\left(\partial_{t}^{(1)}\boldsymbol{v}+\boldsymbol{v}\cdot\boldsymbol{\nabla}^{(1)}\boldsymbol{v}\right)+\boldsymbol{v}\left(\partial_{t}^{(1)}\rho\theta+\boldsymbol{\nabla}^{(1)}\cdot\rho\theta\boldsymbol{v}\right)+\boldsymbol{\nabla}^{(1)}\rho\theta c_{s}^{2}\right] \\ &= -\rho \mathcal{D}\boldsymbol{\nabla}^{(1)}\theta\end{aligned}$$

which may be identified as a Fick/Fourier law of diffusion with a diffusivity defined as: $$\mathcal{D}=c_{s}^{2}\Delta t\left(\frac{\tau}{\Delta t}-\frac{1}{2}\right)$$

Note: From any canonical LB literature reference it can be found that: $$\partial_t^{(1)}\rho+\boldsymbol{\nabla}^{(1)}\cdot\rho\boldsymbol{v}=0$$ $$\partial_t^{(1)}\rho\boldsymbol{v}+\boldsymbol{\nabla}^{(1)}\cdot\rho\boldsymbol{v}\boldsymbol{v}=-\boldsymbol{\nabla}\rho c_s^2$$ This is used in the above derivation to simplify the second to last line using: $$\rho\left(\partial_{t}^{(1)}\boldsymbol{v}+\boldsymbol{v}\cdot\boldsymbol{\nabla}^{(1)}\boldsymbol{v}\right)=-\boldsymbol{\nabla}\rho c_s^2$$ To show it would require another analysis which is similar to this one in context and length, which doesn't add any new information other than this result. The analysis is therefore ommitted.


We now combine the $O\left(\epsilon\right)$ and $O\left(\epsilon^2\right)$ scales: $$\epsilon\partial_t^{(1)}\rho\theta + \epsilon^2\partial_t^{(2)}\rho\theta + \epsilon \boldsymbol{\nabla}^{(1)}\cdot\rho\theta\boldsymbol{v} + \epsilon^2 \boldsymbol{\nabla}^{(1)} \cdot \boldsymbol{j}^{(1)} = \epsilon^2 \phi^{(2)}$$ to give: $$\partial_t\rho\theta + \boldsymbol{\nabla}\cdot\rho\theta\boldsymbol{v} = -\boldsymbol{\nabla}\cdot\left(\rho \mathcal{D}\boldsymbol{\nabla}\theta\right) + \phi $$

which is exactly the advection-diffusion equation with a source term. Combining it with the continuity equation yields: $$\rho\left[\partial_t\theta + \boldsymbol{v}\cdot\boldsymbol{\nabla}\theta\right] = -\boldsymbol{\nabla}\cdot\left(\rho \mathcal{D}\boldsymbol{\nabla}\theta\right) + \phi $$

which is as far as i can see the equation you are looking for. This ofcourse does not include the case of the vorticity because this is a scalar equation. However, if you limit yourself to 2D incompressible flows then the vorticity is a scalar defined by the z-component of the vorticity. In that case this approach may be used.

  • $\begingroup$ Thank you @nluigi, I will go through this asap. I am curious about your comment on 2D vorticity. Is this a very recent development? "A lattice Boltzmann model for the eddy–stream equations in two-dimensional incompressible flows" was in 2011. Do you know any other papers/books that attack this problem in more detail? $\endgroup$ – user133100 Sep 1 '16 at 3:00
  • $\begingroup$ @user133100 - In an incompressible 2D flow by definition the voriticity is $\boldsymbol{w}=\boldsymbol{\nabla}\times\boldsymbol{v}=\boldsymbol{k}w_z$, i.e. it suffices to only look at the z-component of the vorticity vector. You should be able to set $\theta=w_z$ in the analysis above. I think the above analysis should atleast be able to guide you through the derivation in the article you link. I am unfamiliar with canonical literature that tackle this problem; perhaps check Succi or Wolf-Gladrow and see if they have a section on it? $\endgroup$ – nluigi Sep 1 '16 at 15:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.