Preface
- I hope you are ok with a LB derivation (your title suggests you are).
- 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.
- 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.
- 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.
Introduction
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.
Conclusion
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.
