## Abstract

This paper presents an analysis of the simultaneous incorporation of force and mass source terms into the multi-relaxation-time (MRT) collision operator. MRT force incorporation was obtained through Chapman–Enskog analysis. The numerical scheme was tested on different benchmark problems, including the decay of a shear wave with different bulk and kinematic viscosities and axisymmetric flow.

## 1. Introduction

The lattice Boltzmann equation (LBE) is a new technique for solving fluid dynamics and transport problems. Unlike conventional computational fluid dynamics, the LBE is based on the kinetic theory. The kinetic nature brings LBE many advantages and enables it to be an appealing tool for modelling and simulating fluid flows. In the past two decades, LBE has been successfully applied for a variety of hydrodynamics problems [1]. In many fluid flows, the working fluid is usually exposed to an external force field (e.g. gravity) or internal interactions (e.g. van der Waals forces). In order to capture the physics in such systems, it is critical to treat the force correctly in the LBE framework.

It is well understood that LBE models can be classified into two main types in terms of the collision operator employed, i.e. single-relaxation-time (SRT) or Bhatnagar–Gross–Krook (BGK) model and multi-relaxation-time (MRT) model. BGK–LBE models are widely used and have gained much success in solving a variety of fluid flow problems. On the other hand, MRT–LBE models have attracted much attention recently owing to some noticeable features, such as enhanced numerical stability, exact realization of boundary conditions, etc. [2].

There have been successful implementations of a body force in the literature for lattice BGK models [3–5]. However, MRT force implementation has not yet been addressed thoroughly, except for a few publications [6–8]. Furthermore, most of the previous studies on force implementation, for either BGK or MRT–LBE models, did not consider the cases when a source exists in the continuity equation. The first related work on this subject is due to Halliday *et al.* [9], who treated the axisymmetric Navier–Stokes equations in cylindrical coordinates as two-dimensional equations with mass and force terms in pseudo-Cartesian coordinates, which was solved by a lattice BGK model with a forcing term that can replicate the sources in the continuity and momentum equations. Some other versions were also developed later with similar ideas. Another work considering both mass source and body force in LBE is attributed to Ginzburg *et al.* [8], which is developed within the framework of a two-relaxation-time (TRT) model. To the best of the authors’ knowledge, there are no reports so far on implementation of a forcing term accounting for both mass source and body force in the full MRT–LBE.

In this work, we present an analysis of the MRT–LBE with a forcing term that can recover the mass source and body force through Chapman–Enskog expansion. The expressions of the forcing term in moment space are found explicitly. It is found that the force implementation moments depend on the eigenvalues of the collision operator. When all the eigenvalues coincide, then the MRT implementation agrees with that of Guo *et al.* [4]. Some numerical tests, including axisymmetric flow in a concentric pipe, and a decaying wave with different bulk and shear viscosities, are then carried out to validate the proposed algorithm.

The paper is organized as follows. We briefly mention the lattice Boltzmann equation with the force term in §2. The Chapman–Enskog analysis is performed in §3 for the most spread matrix for the two-dimensional D2Q9 model. In §4, the results are compared with analytical solutions for axisymmetric flow and wave decay with different bulk and kinematic viscosities.

## 2. Lattice Boltzmann method

The MRT lattice Boltzmann approach is as follows:
2.1
This can be rewritten for the sake of clarity and simplicity in the matrix model with normalization factors and relaxation times as follows:
2.2
Here the collision matrix *M*_{ij} consists of the normalized eigenvectors and relaxation times:
2.3
and *F*_{i} is the force population responsible for the force inclusion. In this work, we took the two-dimensional nine-velocity (D2Q9) model as an example, where the transform matrix *M* can be found in the work by Lallemand & Luo [2]. The moments related to the matrix are expressed as
2.4
which are known for the equilibrium distribution function as . The algorithm consists of two processes—collision and propagation. The collision is performed in moment space. The propagation is done in the distribution population space. Note that, in comparison with the one relaxation parameter for BGK, the MRT approach operates with nine collision rates corresponding to moments *ω*_{ρ}, *ω*_{e}, *ω*_{ε}, *ω*_{jx}, *ω*_{qx}, *ω*_{jy}, *ω*_{qy}, *ω*_{pxx} and *ω*_{pxy}. If all the relaxation rates equal each other, then the MRT–LBE is equivalent to the most spread BGK–LBE.

Matrix *M*_{ij} has certain features. The most used feature in this work is when the vector *M*^{l} multiplied on matrix *M*_{ij} and summed up by index *i* produces the eigenvector multiplied on the corresponding eigenvalue:
2.5

## 3. The Chapman–Enskog analysis

In the following paragraphs, the equilibrium distribution function and the force populations *F*_{i} will be derived to restore the Navier–Stokes equation with the mass and force terms:
3.1
The moments for the D2Q9 model (2.4) imply the macroscopic quantities as and . However, in the presence of the source terms, the macroscopic parameters are not changed [4,8] in a unique way. One possible choice is to take macroscopic parameters with half shifted source terms [8]:
3.2
where *ϵ* is the Knudsen parameter, and index m stands for the macroscopic parameters. The Chapman–Enskog analysis is done with the expansion of the distribution function in terms of the Knudsen number as . The time derivative is expanded as ∂_{t}=∂_{t0}+*ϵ*∂_{t1}+⋯ . From system (3.2) one can obtain the moments for (note that higher-order terms with *n*>2 are eliminated for the conservation laws to be fulfilled):
3.3
After the substitution of the series expanded distribution function into the LBE (2.2), the infinite consecutive series of equations (Chapman–Enskog system) can be obtained:
3.4
By multiplying equations (3.4) with the eigenvector *M*^{ρ}, one can restore the continuity equation:
3.5

After substitution of the known moments of (3.3) and specification of some of the force moments ( and ), system (3.5) is considerably simplified:
3.6
which restores the continuity equation by taking into account that ∂_{t}=∂_{t0}+*ϵ*∂_{t1}. The next iteration, the Chapman–Enskog system summed up with *M*^{jα}, gives
3.7
where already known force moments are substituted and used. We require terms containing *c*_{iα}*c*_{iβ} to be expressed through the eigenvectors of matrix *M* to perform the summation :
3.8
Therefore, the term can be represented through the sum of terms that have the form . The moments in the expression above can be found from the summation of the second equation of the Chapman–Enskog system (3.4) multiplied with the corresponding eigenvectors :
3.9
After lengthy algebra using the decomposition of the terms as through the matrix eigenvectors, the moments of can be obtained:
3.10
The expressions for and contain partial derivatives ∂_{t0}*ρ*^{m} and . These terms can be obtained through the continuity equations:
3.11
The third-order terms such as and can be neglected:
3.12
After the moments, calculations of (equation (3.10)) and substitution into the system (3.7), the force moments should have the form
3.13
in order to restore the full system of Navier–Stokes equation with a force and the continuity equation with a mass source term (3.1), where *ω*_{η}=*ω*_{pxx}=*ω*_{pxy}, *ω*_{ζ}=*ω*_{e} and the viscosities are . When all the eigenvalues coincide with each other, the force representation coincides with the Guo force term representation [4].

## 4. Numerical results

### (a) Flow in a concentric annular pipe

A flow in a concentric annular pipe driven by a constant pressure gradient is simulated. The radius of the inner and outer cylinders are *r*_{1} and *r*_{2}, respectively. With the transformation (*z*,*r*)→(*x*,*y*) and (*v*_{z},*v*_{r})→(*u*,*v*), the governing equations of the flow can be expressed in pseudo-Cartesian coordinates as [9]:
4.1
where *S*=*ρ*−*v*/*y*, *F*_{x}=(*ρν*/*y*)∂_{y}*u* and *F*_{y}=(*ρν*/*y*)(∂_{y}*v*−*v*/*y*); is the Laplacian operator in the pseudo-Cartesian coordinates. The analytical solution of this flow at steady state is
4.2
where *G*=−∂_{x}*p*/*ρ* is a constant, and *α*=*r*_{2}/*r*_{1} is the radius ratio.

In the test, we set *r*_{1}=100 and *r*_{2}=132 or 164, and the length of the pipe is set to 8. The driven force is 10^{−5}, and periodic boundary conditions are applied to the inlet and outlet of the pipe, and the half-way bounce-back scheme is applied to the surfaces of the inner and outer cylinders. The gradients in *F*_{x} and *F*_{y} are discretized using the second-order central finite-difference scheme. In figure 1, the velocity profiles predicted by the MRT–LBE with the present force formulation are compared with the analytical solutions at different values of *ω*_{η} with *ω*_{ζ}=0.8. Good agreement between the numerical results and the analytical solutions is observed, which confirms the validity of the force scheme.

### (b) Wave decay with different bulk and shear viscosities

The system in the region of the one-dimensional sound waves with the mass term in the advection–diffusion equation can be represented in terms of pressure *p*=*c*^{2}_{s}*ρ* and the flux *j*_{x}=*ρ*_{0}*u*_{x}:
4.3
For simplification, the initial conditions along with the mass and force terms are assumed as
4.4
The overall solution satisfying the initial conditions (4.3) can be written as
4.5
where *λ*_{real}=−*k*^{2}(*η*+*ζ*)/2 and . The numerical simulation is performed on a domain of 40×1 with the parameters *A*=0.001,*B*=0.002 and *C*=0.003. A few examples for *ω*=1.0 and *ω*_{bulk}=0.8,1.2 along with the analytical profiles are given in figure 2. We examined the range of parameters *ω* and *ω*_{bulk} with the given force and mass source representations. The relative error is introduced as *e*=|*j*_{x,analytical}−*j*_{x,simulation}|/*j*_{x,amp}, where *j*_{x,amp} is the amplitude of the signal that bounds the oscillating solution. One of the possible measurements of the error can be a least-squares fit of the signal and simulation to calculate *λ*_{real} and *λ*_{imag}.The results are shown in figure 3.

## 5. Conclusion

This paper demonstrated the simultaneous incorporation of the mass and force terms in the frame of the MRT collision operator. A few benchmark problems are solved, and predicted results were compared with analytical solutions. Results show that the approach can be used for the simulation of fluid flow with source terms. However, more tests need to be done. For example, the results are useful for researchers working in the LBE combustion and reacting flows.

## Acknowledgements

A.K. wants to thank Alberta Ingenuity Fund for their financial support. Z.L.G. is supported by the National Natural Science Foundation of China (10972087). A.A.M. is supported by NSERC, Canada.

## Footnotes

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

- This journal is © 2011 The Royal Society