# Numerical Study of the Properties of the Central Moment Lattice Boltzmann Method

## Abstract

Central moment lattice Boltzmann method (LBM) is one of the more recent developments among the lattice kinetic schemes for computational fluid dynamics. A key element in this approach is the use of *central* moments to specify collision process and forcing, and thereby naturally maintaining Galilean invariance, an important characteristic of fluid flows. When the different central moments are relaxed at different rates like in a standard multiple relaxation time (MRT) formulation based on *raw* moments, it is endowed with a number of desirable physical and numerical features. Since the collision operator exhibits a cascaded structure, this approach is also known as the cascaded LBM. While the cascaded LBM has been developed sometime ago, a systematic study of its numerical properties, such as accuracy, grid convergence and stability for well defined canonical problems is lacking and the present work is intended to fulfill this need. We perform a quantitative study of the performance of the cascaded LBM for a set of benchmark problems of differing complexity, viz., Poiseuille flow, decaying Taylor-Green vortex flow and lid-driven cavity flow. We first establish its grid convergence and demonstrate second order accuracy under diffusive scaling for both the velocity field and its derivatives, i.e. components of the strain rate tensor, as well. The method is shown to quantitatively reproduce steady/unsteady analytical solutions or other numerical results with excellent accuracy. The cascaded MRT LBM based on central moments is found to be of similar accuracy when compared with the standard MRT LBM based on raw moments, when detailed comparison of the flow fields are made, with both well reproducing even small scale vortical features. Numerical experiments further demonstrate that the central moment MRT LBM results in significant stability improvements when compared with certain existing collision models at moderate additional computational cost.

###### pacs:

47.11.Qr,05.20.Dd,47.27.-i^{1}

## I Introduction

Early developments in the area of computational fluid dynamics (CFD) have focused on the solution of the classical discretizations of the continuum description of fluid motion. During the last two decades, there has been much interest and effort in the development of schemes that derive their basis on a more smaller scale picture involving particle motion, which may be classified as mesoscopic methods. One of the most promising of such approaches is the lattice Boltzmann method (LBM) Chen and Doolen (1998); Succi (2001); Luo et al. (2010). Based on kinetic theory, it involves the solution of the lattice Boltzmann equation (LBE), which specifies the evolution of the particle populations along discrete directions, which comprise the lattice. This evolution involves a Lagrangian free streaming process along such lattice links and a local collision step specified as a relaxation process. Various elements involved in these two simple steps are constructed based on symmetry considerations, while obeying certain conservation constraints, in such a way that they recover the dynamics of fluid flow in the near incompressible limit. The resulting scheme has a number of desirable features. These include the ability to naturally represent complex fluid physics such as multiphase and multicomponent flows based on kinetic theory, amenability to parallelization due to the locality of the method and representation of flow through complex geometries. Furthermore, due to the exact conservation in the streaming step and machine round-off conservation in the collision process, it has considerably low numerical dissipation for a second-order numerical scheme Ubertini et al. (2010). Due to such competitive advantages, the LBM has found applications in the simulation of a wide range of fluid flow problems Chen and Doolen (1998); Succi (2001); Luo et al. (2010).

Since the LBM is usually developed by means of a bottom-up strategy, there is certain level of flexibility in the construction of its various elements to recover the macroscopic fluid motion. In particular, the choice of a suitable collision model can have profound influence on the fidelity as well as the stability of the approach. As such, the construction of the collision step has been the subject of considerable attention since the inception of the LBM. The simplest among these is the so-called single-relaxation-time (SRT) model Chen et al. (1992); Qian et al. (1992), which is based on the Bhatnagar-Gross-Krook (BGK) approximation Bhatnagar et al. (1954). While it is popular, it has limitations in the representation of certain flow problems and is generally prone to numerical instability, particularly at high Reynolds numbers. A major development to address these aspects is the moment approach d‘Humières (1992), which has been constructed based on multiple relaxation times (MRT) in particular to significantly improve the numerical stability Lallemand and Luo (2000). While it is related to its precursor involving a more general relaxation approximation Higuera and Jiménez (1989); Higuera et al. (1989), the characteristic difference being that it performs collision in an orthogonal moment space leading to an efficient and flexible numerical scheme. This moment approach, which is designated as the standard MRT formulation in this paper, has recently been studied and compared with some of the other collision models in detail Luo et al. (2011). A simpler version that is intermediate between the SRT and MRT model is the so-called two-relaxation-time (TRT) model Ginzburg (2005), in which the moments of even and odd orders are relaxed to their equilibrium at different rates. This, along with the MRT model, can be adjusted such that it results in a minimization of undesirable discrete kinetic effects near walls. Another significant development is the so-called entropic LBM Karlin et al. (1999). It involves an equilibria, which is based on a constrained minimization of a Lyapunov-type functional. By modulating the collision process through enforcing entropy involution locally, this approach aims to maintain non-linear stability. This approach has resulted in a number of simplified variants recently Asinari and Karlin (2009); Karlin et al. (2011).

An important physical feature of the fluid motion is that their description be independent of any inertial frame of reference (e.g. Pope (2000)). This invariance property, which is termed as the Galilean invariance, should be satisfied by any model or numerical scheme for its general
applicability. Furthermore, it has recently been shown that stabilization of classical schemes for compressible flow can be achieved when they are specifically constructed to respect this physical property Scovazzi (2007a, b); Hughes et al. (2010). Keeping these general notions in mind, Galilean invariance can be naturally prescribed in the LBM when its various elements are represented in terms of the *central* moments, i.e. moments obtained by shifting the particle velocity by the local fluid velocity. That is, any dynamical changes due to the collision process and impressed forces can be represented in terms of suitable variations of a set of such central moments. In particular, a collision model based on the relaxation of central moments was constructed recently Geier et al. (2006). The model exhibits a cascaded structure, which was later shown to be equivalent to considering a generalized equilibrium in the lattice or rest frame of reference Asinari (2008). These central moments can be relaxed at different rates during
collision leading to a cascaded MRT or central moment MRT formulation, whereas by contrast the standard MRT formulation considers *raw* moments. A systematic derivation of this approach by including the effect of impressed forces based on central moments was presented in Premnath and Banerjee (2009). This leads to considering generalized sources, analogous to the generalized equilibrium in the rest frame of reference. They also presented a detailed Chapman-Enskog analysis of the cascaded MRT LBM for its consistency with the macroscopic fluid dynamical equations of motion. This approach was further extended to various lattice models in three-dimensions in Premnath and Banerjee (2011), in the cylindrical coordinate system for axisymmetric flows in Premnath and Ning (2012) and for accounting of non-equilibrium effects in Premnath and Banerjee (2012).

Prior work on the cascaded LBM as discussed above have focused mainly on method developments or their mathematical analysis, with little attention towards their numerics except for few validation cases. In particular, a detailed numerical study of the properties of the cascaded LBM for established benchmark problems and also their performance against other LBM approaches is lacking. The focus of the present work is intended to fill this gap by presenting a systematic study of the numerical properties of the cascaded LBM, viz., grid convergence, accuracy and stability for various canonical problems of differing complexity in terms of flow features and temporal evolution. Establishing the reliability and merits of the method in quantitative terms could provide confidence in their extension and applications to various complex flow problems of interest. To study the numerics of the cascaded LBM, we consider the Poiseuille flow, decaying Taylor-Green vortex flow, and lid-driven cavity flow, for which either analytical solutions or detailed prior numerical results are available for comparison. Much of the literature on the LBM with other collision models on grid convergence studies have focused only on those for the velocity field. In this work, we present numerical results on the grid convergence of the cascaded LBM for the velocity field as well as its derivatives, i.e. the strain rate tensor. Furthermore, an advantage of the kinetic schemes such as the LBM is that the strain rate tensor can be computed locally in terms of non-equilibrium moments. In this work, we also present a direct comparison of the results obtained using the non-equilibrium moments of the cascaded LBM with those involving the finite differencing of the velocity field at various locations for the lid-driven cavity flow problem to assess their quantitative accuracy. It may be noted that a detailed comparison study of the SRT and the standard MRT models have recently been performed in Luo et al. (2011). Thus, in this work, we present a quantitative accuracy comparison between the standard MRT LBM and the cascaded or central moment MRT LBM for the lid-driven cavity flow. Finally, we will discuss the numerical stability performance of the various LBM schemes for the above benchmark problem.

The paper is organized as follows. Section II presents the details of the particular version of the cascaded MRT LBM used in this work. In Sec. III, the results of the grid convergence study of the cascaded MRT LBM together with the raw moment based standard MRT LBM for the three benchmark problems are discussed. Subsequently, the quantitative accuracy of the cascaded LBM is demonstrated by making detailed comparison with either analytical or other numerical solutions for the above problems in Sec. IV. In Sec. V, numerical stability test results are presented for the lid-driven cavity flow using the SRT LBM, standard MRT LBM and cascaded MRT LBM. Summary and conclusions of this work are given in Sec. VI.

## Ii Cascaded Lattice Boltzmann Method

We will now discuss the main features of the cascaded LBM. Similar to the standard MRT LBM, the cascaded MRT LBM also performs collisions in moment space, but these moments are obtained by shifting the particle velocity by the local fluid velocity, i.e. using central moments. As a result, the approach can naturally maintain Galilean invariance. Central moment relaxation process was specified in Geier et al. (2006), which was re-interpreted by considering generalized equilibrium in Asinari (2008). Its detailed mathematical consistency analysis in a MRT formulation with forcing was carried out in Premnath and Banerjee (2009). The computations of the cascaded LBM are actually performed after transforming the central moments into raw moments by means of a binomial formula. In this work, the specific formulation of the cascaded LBM given in Premnath and Banerjee (2009), whose details are somewhat different from that given in Geier et al. (2006), is used. This is briefly discussed in what follows.

In this work, the standard two-dimensional, nine velocity (D2Q9) lattice is employed. We consider the usual bra-ket notations in the description of the method as it provides a convenient representation. That is, we consider the depiction of vectors as and , where represents a row vector of of any state in the corresponding direction and represents a column vector . The inner product is then denoted by . As the cascaded LBM is a moment approach, we need a set of nine linearly independent moment basis vectors for its specification. The (raw) moments of the distribution function of different orders can be defined as . Here, is the discrete particle direction, and and are integers. Thus, a set of nine linearly independent nonorthogonal basis vectors obtained using the monomials in an ascending order can be written as

(1) |

This can be transformed by means of the Gram-Schmidt procedure into an equivalent set of *orthogonal* basis vectors, which provides a computationally more efficient and convenient setting for the description of the method. As a result, we have the following orthogonal set Premnath and Banerjee (2009):

(2) |

Collecting the above set of vectors as a matrix , it immediately follows that is a diagonal matrix, owing to orthogonality. This orthogonal matrix can be written in component form as

(3) |

To specify the collision step and forcing, we need the central moments of the local equilibrium and sources, which can be obtained as follows. First, the local Maxwell-Boltzmann distribution function in continuous particle velocity space is written as , where is the speed of sound. Typically, . Based on this, the continuous central moments of the equilibrium of order can be defined as , which yields

(4) |

Considering that the impressed forces only influence the fluid momentum, the central moments of the sources of order due to a force field defined by , where is the change in the distribution function due to force fields, can be simply written as Premnath and Banerjee (2009)

(5) |

Based on the above continuous central moments, the elements of the cascaded LBE can be formulated. Using the trepezoidal rule representation of the source term, the cascaded LBE can be written as Premnath and Banerjee (2009)

(6) |

Here, the collision term can be represented as , where is the vector of distribution functions and is the vector of unknown collision kernel to be obtained later. Owing to the cascaded nature of the central moment based approach, it satisfies the following functional relation . The discrete form of the source term in the cascaded LBE given above represents the influence of the force field in the velocity space and is defined as . Noting that Eq. (6) is semi-implicit, by using the standard variable transformation , its implicitness can be effectively removed. This yields

(7) |

The derivation of the collision term, i.e. the collision kernel and the source term involves matching the *discrete* central moments and the *continuous* central moments of equilibria and sources, which are specified above, of all orders supported by the lattice set. We designate this step as the *Galilean invariance matching principle*. First, the discrete central moments of the distribution functions and sources of order can be defined, respectively, as
and
. Also, in terms of the transformed distribution functions
we define , which satisfies
, and similarly for the local equilibria
. Then, the Galilean
invariance matching principle reads

(8) | |||

(9) |

This immediately specifies the various discrete central moments. Hence, we get

(10) |

(11) |

and

(12) |

The next important step is to transform all the above discrete central moments in terms of raw moments, which can be readily accomplished by means of the following binomial formula: , where . Thus, we obtain the following discrete raw moments of sources as

(13) |

Based on the above, we now obtain the source terms projected to the orthogonal moment basis vectors, i.e. , . This intermediate step is needed to obtain the source terms in the velocity space. It immediately follows that

Equivalently, this can be written in matrix form as . By exploiting the orthogonal property of , i.e. , where the diagonal matrix is , we exactly invert the above to obtain the source terms in velocity space as

(14) |

The discrete raw moments of the transformed distribution functions , which will be needed in the evaluation of the collision kernel, can be conveniently written as follows:

(15) |

where we have used , with , , as a compact summation operator for ease of presentation. Furthermore, the raw moments of the collision kernels are needed in its construction. Collision invariants of conserved moments imply . Exploiting the orthogonal property of the matrix , the non-conserved moments of at higher orders, i.e. can be obtained as follows Premnath and Banerjee (2009):

(16) | |||||

Using the above, the collision kernel of the cascaded collision operator can be obtained as follows. Starting from the lowest order central moments that are non-collisional invariants (i.e. and higher), they are successively set equal to their local attractors based on the transformed equilibria. This step provides tentative expressions for based on the equilibrium assumption. This is then modified to allow for relaxation process during collision. That is, they are multiplied with corresponding relaxation parameters Geier et al. (2006). In this step, care needs to be exercised to multiply the relaxation parameters only with those terms that are not yet in post-collision states (i.e. terms not involving ) for a given . See Premnath and Banerjee (2009) for various details involved in this procedure. Here, we summarize the final expressions of the non-conserved collision kernels, which are given as follows:

(17) | |||||

(18) | |||||

(19) | |||||

(20) | |||||

(21) | |||||

(22) | |||||