## Abstract

A lattice Boltzmann (LB) approach is presented for solving scalar transport equations. In addition to the standard LB for fluid flow, a second set of distribution functions is introduced for transport scalars. This LB approach fully recovers the macroscopic scalar transport equation satisfying an exact conservation law. It is numerically stable and scalar diffusivity does not have a Courant–Friedrichs–Lewy-like stability upper limit. With a sufficient lattice isotropy, numerical solutions are independent of grid orientations. A generalized boundary condition for scalars on arbitrary geometry is also realized by a precise control of surface scalar flux. Numerical results of various benchmarks are presented to demonstrate the accuracy, efficiency and robustness of the approach.

## 1. Introduction

The lattice Boltzmann method (LBM) has now become a promising alternative numerical method to conventional computational fluid dynamics (CFD) for simulating various fluid flows. Unlike conventional CFD based on directly solving macroscopic continuum equations, the LBM starts from a mesoscopic kinetic (i.e. Boltzmann) equation to determine macroscopic fluid dynamics [1]. Its kinetic nature brings many advantages over the conventional methods, such as easy handling of complex geometries, parallel computation, low numerical dissipation and efficient simulations of multi-phase flows with thermodynamic phase transitions (cf. [1–4]). It has been applied extensively in research and real engineering applications.

In many fluid flow simulations, additional scalar quantities usually need to be solved; for example, temperature variation in heat transfer problems. Because of lack of a global H-theorem [5], most LB approaches with full energy conservation encounter challenging numerical instability issues unless high-order LB models are applied [6]. Simultaneous realization of mass, momentum and temperature boundary conditions (BCs) is also extremely challenging. To overcome these difficulties, one solution is to solve an additional scalar equation for temperature [4].

Generally, there are two types of methods to solve scalar transport equations along with a basic LB fluid solver, one is the finite difference (FD) method and the other is an LB scalar solver. Because of the Lagrangian nature of the background LB fluid solver, explicit FD schemes are commonly applied on the same cubic lattice grids for convenience to solve the scalar equation. There are many explicit FD schemes in the literature. Nonetheless, owing to the very strict Courant–Friedrichs–Lewy (CFL) constraint, there exists a stability upper bound for the normalized value of diffusivity, and also a mixture of upwind and central schemes for scalar advection is usually necessary [7]. A direct consequence of such a treatment is the increased numerical diffusion and loss of exact Galilean invariance. Symmetry of the schemes could also be sacrificed, which causes grid orientation dependence of numerical solutions. Furthermore, it is difficult to accurately calculate scalar gradient and scalar advection in regions containing arbitrary geometries due to missing neighbour information. Thus, near wall numerical accuracy could be significantly compromised.

There are not many works on LB scalar solvers in the literature. In 1987, the concept of using multi-component formulation was proposed in lattice gas framework [8]. In 1997, Shan [9] used a separate distribution function to solve the temperature equation and conducted Rayleigh–Benard simulations in two and three dimensions. Later on He *et al*. [10] proposed a similar approach to solve internal energy distribution function. Peng *et al.* [11] further extended the internal energy model for three-dimensional incompressible thermal simulations. One of the advantages of these LB scalar solvers is numerical stability. Unlike explicit FD, LB scalar solvers naturally satisfy the CFL condition. The scalar advection is exact and there is no need to calculate scalar gradient explicitly. Scalar diffusivity does not have an upper limit.

Although the existing LB scalar solvers possess many advantages, they all fail on a fundamental invariance condition. That is, if there is no external scalar source, a uniform scalar distribution should be maintained at all times no matter the background fluid flow.

In this paper, a new LB scalar approach is proposed that addresses the above problem. In addition, a generalized BC for the scalar on arbitrary geometry, which is largely missing in previous studies, is also presented. Theoretical details are described in §2 followed by results for simple but representative benchmark problems to demonstrate the accuracy, stability and symmetry of the new LB scalar solver. Conclusions are given in §4.

## 2. New lattice Boltzmann scalar solver

### (a) The fluid solver

Various types of LBM could be applied for solving fluid flows, which serve as the background carrier for scalar transport. For simplicity, the standard D3Q19 Bhatnagar–Gross–Krook (BGK) model [1] is used in this study:
2.1
Here, *f*_{i}(**x**,*t*)(*i*=1,…,19) is the particle distribution function, *τ* is the single relaxation time and is the equilibrium distribution function with a third-order expansion in fluid velocity,
2.2
where *T*_{0}=1/3. The discrete lattice velocities **c**_{i} are
2.3
with *w*_{0}=1/3 for rest particle, *w*_{i}=1/18 for states of Cartesian directions and *w*_{i}=1/36 for states of bi-diagonal directions. *g*_{i}(**x**,*t*) is the external body force term [12]. The hydrodynamic quantities *ρ* and **u** are moments of the particle distribution function
2.4

### (b) The scalar solver

In addition to the fluid solver, a separate set of distribution functions, *T*_{i}, is introduced for scalar transport:
2.5
2.6
and
2.7
*T*_{i} is the scalar distribution function and *T* is the scalar being solved for. *τ*_{T} is the relaxation time corresponding to scalar diffusivity. The *f*_{i}, *ρ*, *T*_{0} and **u** are defined in equations (2.1), (2.2) and (2.4), respectively. In general, the lattice speed set for scalar distribution does not need to be the same as the lattice speed set for fluid distribution because the current scalar solver is an additional system attached to the basic fluid solver. Different lattice speed set for scalar could be applied as long as the scalar lattice speed set is a subset of the fluid lattice speed set. For example, a six-speed LB model may be used for scalar simulation when the 19-speed LB model is used for fluid simulation. Since the 19-speed LB model has a higher order lattice symmetry than the six-speed LB model, the same 19-speed lattice model for scalar is used here.

Unlike the regular BGK, the so-called BGK regularized/filtered collision operator form is used here [13,14]. Such a form ensures that only the first-order non-equilibrium moment contributes to scalar diffusion in the hydrodynamic range. All non-equilibrium moments of higher orders are filtered out by this collision process. Its direct benefits include elimination of numerical noise exhibited in BGK and improved robustness. The scalar *T* serves as its own equilibrium and no complicated expression of scalar equilibrium distribution function is needed. The overall calculation of collision operator *Φ*_{i} is rather efficient.

It can be shown easily that the collision in equation (2.5) obeys the scalar conservation law. Multiplying *f*′_{i}(**x**,*t*)=*f*_{i}(**x**+**c**_{i},*t*+1) on both sides of equation (2.5) and noticing
2.8
we have
2.9
where *T*′_{i}(**x**,*t*) denotes the right-hand side of equation (2.5). Hence, the scalar collision operator conserves local *ρT*, which implies realization of local energy conservation if the scalar is considered as temperature. Since *T*_{i} propagates along with *f*_{i}, the energy distribution *E*_{i}=*f*_{i}*T*_{i} is fully maintained during advection. The global conservation of *ρT* is therefore achieved. Furthermore, and most notably, this scheme maintains the exact invariance on uniformity of scalar *T*. It is straightforward to see that if everywhere, then *Φ*_{i}(**x**,*t*)≡0, and everywhere at all later times, regardless of the background flow field. This fundamental property is not demonstrated in any previous LB scalar models.

Using Chapman–Enskog expansion, it can be shown that equation (2.5) recovers the following second-order macroscopic scalar transport equation:
2.10
with *κ*=(*τ*_{T}−1/2)*T*_{0}. The uniformity invariance condition ensures that *ρ* is outside of ∇*T*.

### (c) Boundary condition

One substantial advantage of the LBM is the capability of dealing with complex geometry [1]. In 1998, Chen *et al*. [15] formulated a general volumetric LB surface algorithm to achieve frictionless/frictional BCs on arbitrary geometry. In their scheme, mass is conserved, both tangential and normal momentum fluxes on the boundary are precisely realized. The local detailed balance is fully satisfied. Many benchmark studies have well demonstrated the accuracy and robustness of this scheme [2,3,16]. An adiabatic (zero scalar flux) BC on arbitrary geometry for scalar can be derived as a direct extension of this scheme. Once the adiabatic BC is realized, a prescribed finite flux BC can be accomplished.

Unlike other point-wise LB, BCs are conducted on a discretized set of surface elements. These piece-wise flat surface elements together represent a curved geometry. During particle advection, each surface element collects incoming particles from its neighbouring fluid cells. The incoming distributions are weighted by volume overlapping of parallelepipeds from the underling surface element with cells in particle moving directions (cf. [15]). After receiving the incoming quantities, the following surface scalar algorithm is applied to determine the outgoing distributions from the surface:
2.11
where
2.12
Here, **n** is the surface normal pointing towards the fluid domain and **c**_{i}⋅**n**<0, **c**_{i*}=−**c**_{i}. *P*_{i} (≡|**n**⋅**c**_{i}|*A*) is the volume of parallelepiped in particle direction **c**_{i} associated with the surface normal **n** and area *A* of a given surface element, and obviously *P*_{i}=*P*_{i*}. Finally, the outgoing distributions are propagated back from the surface element to fluid cells according to the same surface advection process as in Chen *et al*. [15].

It is not difficult to show that the above surface scalar collision achieves exact zero surface scalar flux. Taking summation over outgoing directions, the outgoing scalar flux is
2.13
Note, *P*_{i}=*P*_{i*} and the definition of *T*^{in} in equation (2.12), the second summation term on the right-hand side is zero. In addition, because of the mass flux conservation , the total outgoing scalar flux is the same as the total incoming scalar flux
2.14
Therefore, zero net surface flux (adiabatic) BC is fully satisfied on arbitrary geometry.

If an external scalar source *Q*(*t*) is specified on the surface, a source term can be directly added to equation (2.11)
2.15
If the BC has a prescribed scalar quantity *T*_{w} (for example, surface temperature), surface heat flux can be calculated accordingly:
2.16

## 3. Numerical verification

Four sets of simulation results are presented to demonstrate the capability of the LB scalar solver regarding its numerical accuracy, stability, Galilean invariance, grid orientation independence, etc. Results using two different second-order FD schemes, van Leer type of flux limiter scheme and direct mixing scheme (mixture of central and first-order upwind schemes), are also presented as comparisons.

### (a) Shearwave decay

The first test case is a temperature shearwave decay carried by a constant uniform fluid flow. The initial temperature distribution is a uniform one plus a spatial sinusoidal variation with lattice wavelength *L*=16 and magnitude . *T*_{A} is a constant. The velocity of background mean flow is 0.2 and the thermal diffusivity *κ* is 0.002. With such a low resolution and *κ*, the numerical stability and accuracy can be well validated. For temperature decay without background flow, both the LB scalar solver and the FD methods show excellent agreements with theory. With non-zero background mean flow, the LB scalar solver is still able to accurately compare with theory. However, the FD results show noticeable numerical errors. The temperature profiles at lattice time step 81 are plotted in figure 1. Numerical diffusion is seen clearly for the flux-limiter FD scheme, while neither the correct temperature profile nor its location can be maintained by the mixing FD scheme.

### (b) Inclined channel with volume heat source

The second test case is a simulation of temperature distribution in a channel flow with different lattice orientations. The channel walls are free-slip, and the fluid flow stays uniform with *U*_{0}=0.2 as a result. The channel width is 50 (lattice spacing) and the flow *Re* is 2000. The thermal diffusivity *κ* is 0.005. The temperature on the wall is fixed at *T*_{w}=1/3. A constant volume heat source *q*=5×10^{−6} is applied in the bulk fluid domain. The flow has periodic BC in the streamwise direction, which is easy to realize in lattice aligned situation. When the channel (light colour) is tilted as shown in figure 2, the in and out channel boundaries are matched perfectly in coordinate directions so that the periodic BC is once again realized in the streamwise direction. In order to demonstrate lattice independence, we choose the tilted angle to be 26.56^{°}. Like the first test case, the temperature distributions using the LB scalar solver and the two FD schemes match the analytic solution very well when the channel is lattice aligned. However, the results from the FD schemes depart significantly from theory when the channel is tilted. The simulation results of temperature distribution across the tilted channel are shown in figure 2. The LB results are clearly shown to be independent of lattice orientations. The errors from the FD methods are also originated from their fundamental difficulty in dealing with gradient calculation on a tilted boundary orientation, so that additional numerical artefacts are introduced. Since the LB scalar particle advection is exact with the BC presented here, it is thus able to achieve a lattice orientation independent scalar evolution.

### (c) Temperature propagation in an inclined channel

Owing to lack of neighbour information in a non-lattice aligned near wall region on Cartesian grids, it is extremely difficult to get accurate estimation of local gradients, which is essential for FD based methods. Furthermore, because of strong dependency on upwind information, the calculation of scalar advection can be further compromised for FD methods. In contrast, the boundary treatment of the LB scalar solver achieves exact local scalar flux conservation as discussed above. The scalar advection in such a near wall region can be computed accurately.

High temperature convection in a channel tilted by 30^{°} is conducted as a demonstration. The free-slip and adiabatic BCs are enforced at solid walls and fluid velocity is constant *U*_{0}=0.0909 along the channel. The thermal diffusivity *κ* is 0.002. Initially, the temperature is 1/3 everywhere except for *T*=4/9 at the inlet. Then this temperature front should be convected by the uniform background fluid flow to downstream locations at later times without distortion. The computed temperature front distributions across the channel at lattice time step 2000 are shown in figure 3. The temperature front of the LB scalar solver maintains a nearly straight profile. On the other hand, the temperature fronts from the two FD schemes show substantial distortions in near wall regions. It is also worth mentioning that the LB scalar result shows the thinnest temperature front which implies that the LB scalar solver has a smaller numerical diffusion.

### (d) Rayleigh–Bénard natural convection

Rayleigh–Bénard (RB) natural convection is a classical benchmark for accuracy verification of numerical solvers. It has a simple case set-up but complex physics phenomena. When Rayleigh number *Ra* exceeds a certain critical value, the system experiences a transition from no-flow to convection. Current study is carried out under the Boussinesq approximation, in that the buoyancy force acting on the fluid is defined as *αρ***g**(*T*−*T*_{m}) where *α* is the thermal expansion rate, **g** is gravity, and *T*_{m} is the average temperature value of the top and bottom boundaries. Since the most unstable wave number is *k*_{c}=3.117 when *Ra* exceeds the critical value *Ra*_{c}=1707.762, the resolution 101×50 is used in the study. *Pr* used here is 0.71.

A set of simulations with various *Ra* are conducted. Owing to page limitation, only Nusselt number *Nu* versus *Ra* is presented here. When RB convection is established, the heat transfer between two plates is greatly enhanced. The enhancement of the heat transfer is described by *Nu*=1+〈*u*_{y}*T*〉*H*/*κ*Δ*T*. Our results are compared with those of Clever & Busse [17] in figure 4. As shown, the agreements are very good for *Ra* less than 20 000. Our results slightly underestimate the heat transfer at higher *Ra*. Other LB methods also showed a similar trend [9,10].

## 4. Conclusions

A new LB scalar solver including a generalized BC on arbitrary geometry is presented in this paper. The algorithm is simple and efficient. The fundamental flaw associated with lack of invariance of uniform temperature in the previous LB scalar models is fully removed. Its numerical stability, accuracy, symmetry and robustness are demonstrated by some representative benchmark cases. Finally we would like to point out that this LB scalar solver can be trivially extended to simulate multiple scalars.

## Footnotes

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

- This journal is © 2011 The Royal Society