# Galilean Invariant Preconditioned Central Moment Lattice Boltzmann Method without Cubic Velocity Errors for Efficient Steady Flow Simulations

###### Abstract

Lattice Boltzmann (LB) models used for the computation of fluid flows represented by the Navier-Stokes (NS) equations on standard lattices can lead to non Galilean invariant (GI) viscous stress involving cubic velocity errors. This arises from the dependence of their third order diagonal moments on the first order moments for standard lattices, and strategies have recently been introduced to restore GI without such errors using a modified collision operator involving either corrections to the relaxation times or to the moment equilibria. Convergence acceleration in the simulation of steady flows can be achieved by solving the preconditioned NS equations, which contain a preconditioning parameter that can be used to tune the effective sound speed, and thereby alleviating the numerical stiffness. In the present study, we present a GI formulation of the preconditioned cascaded central moment LB method used to solve the preconditioned NS equations, which is free of cubic velocity errors on a standard lattice, for steady flows. A Chapman-Enskog analysis reveals the structure of the spurious non-GI defect terms and it is demonstrated that the anisotropy of the resulting viscous stress is dependent on the preconditioning parameter, in addition to the fluid velocity. It is shown that partial correction to eliminate the cubic velocity defects is achieved by scaling the cubic velocity terms in the off-diagonal third-order moment equilibria with the square of the preconditioning parameter. Furthermore, we develop additional corrections based on the extended moment equilibria involving gradient terms with coefficients dependent locally on the fluid velocity and the preconditioning parameter. Such parameter dependent corrections eliminate the remaining truncation errors arising from the degeneracy of the diagonal third-order moments and fully restores GI without cubic defects for the preconditioned LB scheme on a standard lattice. Several conclusions are drawn from the analysis of the structure of the non-GI errors and the associated corrections, with particular emphasis on their dependence on the preconditioning parameter. The new GI preconditioned central moment LB method is validated for a number of complex flow benchmark problems and its effectiveness to achieve convergence acceleration and improvement in accuracy is demonstrated.

###### pacs:

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

^{†}preprint: PREPRINT

## I Introduction

The lattice Boltzmann (LB) method has now been established as a powerful kinetic scheme based computational fluid dynamics approach (Succi (2001), Guo and Shu (2013), Kruger et al. (2016)). It is a mesoscopic method based on local conservation and discrete symmetry principles, and may be derived as a special discretization of the Boltzmann equation. Algorithmically, it involves the streaming of the particle distribution functions as a perfect shift advection step along the lattice directions and followed by a local collision step as a relaxation process towards an equilibria, and accompanied by special strategies for the implementation of impressed forces. The hydrodynamic fields characterizing the fluid motion are then obtained via the various kinetic moments of the evolving distribution functions and its consistency to the Navier-Stokes (NS) equations may be established by a Chapman-Enskog expansion or Taylor series expansions under appropriate scaling between the discrete space step and time step. As such, the LB method has been applied for the computation of a wide range of complex flows including turbulence, multiphase and multicomponent flows, particulate flows and microflows (Chen and Doolen (1998), Aidun and Clausen (2010)). Its various appealing features, including its inherent parallelism, natural framework to incorporate kinetic models for complex flows and the ease of boundary conditions has made it an unique and efficient approach for computational fluid dynamics (CFD). During the last decade, many efforts were made to further improve its numerical stability, accuracy and efficiency. In particular, sophisticated collision models based on multiple relaxation times and involving raw moments, central moments or cumulants, and entropic formulations have significantly expanded the capabilities of the LB method. The significant achievements of these developments and their applications to a variety of complex flow problems have been discussed, for example, in Aidun and Clausen (2010); Luo et al. (2010); Geller et al. (2013); Karlin et al. (2014); Frapolli et al. (2015); Succi (2015); Mazloomi et al. (2015); Ning et al. (2015); Krafczyk et al. (2015); Yang et al. (2016); Dorschner et al. (2016); Hajabdollahi and Premnath (2018).

There exist various additional aspects in the LB approach that require further attention and present scope for improvements. In particular, the finiteness of the lattice can introduce certain truncation errors that manifest as non-Galilean invariant viscous stress, i.e. fluid velocity dependent viscosity. This lack of Galilean invariance (GI) arises due to the fact that the diagonal terms in the third-order moments are not independently supported by the standard tensor-product lattices (i.e. D2Q9 and D3Q27). More specifically, for example,

Here, and in the following, the primed quantities denote raw moments. In other words, there is a degeneracy of the third-order diagonal (longitudinal) moments that results in a deviation between the emergent macroscopic equations derived by the Chapman-Enskog expansion and the Navier-Stokes (NS) equations. Such cubic-velocity truncation errors are grid independent and persist in finer grids especially under high shear and flow velocity. Moreover, such emergent anisotropic viscous stress may have a negative impact on numerical stability as a result of a negative dependence of the emergent viscosity on the fluid velocity. In order to overcome this shortcoming, various attempts have been made.

One possibility is to consider a lattice with a larger particle velocity set, such as the D2Q17 lattice in two-dimensions Qian and Zhou (1998), which was pursued after Qian and Orszag (1993) pointed out nonlinear, cubic-velocity deviations of the emergent equations of the LB models with standard lattice sets from the NS equations. This involved the use of higher order velocity terms in the equilibrium distribution. However, Hazi and Jimenez (2006) showed that the specific equilibria adopted in Qian and Zhou (1998) does not fully eliminate the cubic-velocity errors. Moreover, the use of non-standard lattice stencils with larger number of particle velocities increases the computational cost and propensity of the numerical instability at grid scales. On the other hand, it was shown more recently by various others (Wagner and Yeomans (1999), Hazi and Jimenez (2006), Keating et al. (2007)) that partial corrections to the GI errors on the standard lattice (i.e. D2Q9 lattice) may be achieved by adopting special forms of the off-diagonal, third-order moments in the equilibria. That is,

Here, is the speed of sound and the particular choices of the cubic-velocity terms that are underlined are crucial to partially restore GI for the above identified moments. Here, we also point out that the above forms of the off-diagonal, third-order raw moment equilibria that allow such partial GI corrections naturally arise in the central moment LB formulations, when the equilibrium central moment components are set to zero and and then rewritten in terms of their corresponding raw moments. However, since and due to the degeneracy of the third-order longitudinal moments, which is inherent to the standard tensor-product lattices, additional corrections are required to restore GI completely free of cubic-velocity errors. In this regard, in order to compensate the terms which violate GI on standard lattices, LB schemes with single relaxation time models were augmented with finite difference expressions Prasianakis and Karlin (2007, 2008); Prasianakis et al. (2009). On the other hand, more recently, Dellar (2014) introduced small intentional anisotropies into a matrix collision operator that corrects the anisotropy in the resulting viscous stress tensor thereby addressing the above mentioned issue. In addition, independently, Geier et al. (2015) introduced additional corrections involving velocity gradients to the equilibria that achieved equivalent results. These two studies provided strategies to represent the Navier-Stokes equations in LB models on standard lattices completely free of cubic-velocity errors. In addition, Dubois et al. (2017) presented finite difference based corrections to the method proposed in Asinari et al. (2012) to reduce the resulting spurious velocity dependent viscosity effects on standard lattices.

While the LB schemes have found applications to a wide range of fluid flow problems, there has also been considerable interest to an important class of problems related to low Reynolds number steady state flows. They include analysis and design optimization of a variety of Stokes flows through capillaries, porous media flows, heat transfer problems under stationary conditions, and since the LB methods are explicit marching in nature, efficient solution techniques need to devised to accelerate their convergence (see e.g. Verberg and Ladd (1999); Bernaschi et al. (2001); Tölke et al. (2002); Guo et al. (2004); Mavriplis (2006); Pingen et al. (2008); Izquierdo and Fueyo (2008); Premnath et al. (2009); Hübner and Turek (2010); Patil et al. (2014); Hu et al. (2015); Atanasov et al. (2016); F.Hajabdollahi and Premnath (); De Rosis (2017); Hajabdollahi and Premnath (2017)). A recent review of the literature in the LB approach for such problems can be found in F.Hajabdollahi and Premnath (); De Rosis (2017). Generally, multigrid and preconditioning techniques can be devised to improve the steady state convergence of the LB scheme. A comparison of a multigrid LB formulation with the conventional solvers showed significant improvement in efficiency Patil et al. (2014). At low Mach numbers, the convergence can be further accelerated by means of preconditioning for both the traditional single grid LB methods Guo et al. (2004); Izquierdo and Fueyo (2008); Premnath et al. (2009); F.Hajabdollahi and Premnath (); De Rosis (2017) and multigrid LB scheme Hajabdollahi and Premnath (2017). The present work addresses a further refinement to the LB techniques for steady state flows, viz., improving the accuracy of the acceleration strategy based on the preconditioned LB formulation without GI cubic velocity and parameter dependent errors.

Thus, it is clear that another aspect of the LB method, similar to certain schemes based on the classical CFD, is its slow convergence to steady state at low Mach numbers. In such conditions, there is a relatively large disparity between the sound speed and the convection speed of the fluid motion resulting in higher eigenvalue stiffness and larger number of iterations for convergence. This stiffness can be alleviated and convergence can be significantly improved by preconditioning. Reference Guo et al. (2004) presented a preconditioned LB method based on a single relaxation time model by modifying the equilibrium distribution function by using a preconditioning parameter. Then, Izquierdo and Fueyo (2008) and Premnath et al. (2009) presented preconditioned LB formulations based on multiple relaxation times. More recently, F.Hajabdollahi and Premnath () presented a preconditioned scheme for the central moment based cascaded LB method Geier et al. (2006) in the presence of forcing terms Premnath and Banerjee (2009) and demonstrated significant convergence acceleration.

In general, such preconditioned LB schemes are intended to solve the preconditioned NS equations, which can be written as (Choi and Merkle (1993), Turkel (1999))

(1a) | |||

(1b) |

where , and are the pressure, strain rate tensor and the impressed force, respectively. Here, is the preconditioning parameter, which can be used to tune the pseudo-sound speed, thereby alleviating the eigenvalue stiffness and improving convergence acceleration (e.g. F.Hajabdollahi and Premnath ()). However, the existing LB models for the preconditioned NS equations are not Galilean invariant and are expected to involve both velocity- and parameter-dependent anisotropic form of the viscous stress tensor. Development of the Galilean invariant preconditioned central moment based LB method without cubic-velocity defects and parameter free truncation errors for steady flow simulations is the main focus of this study. It may be noted that the preconditioned NS equations may be considered as a specific example of what may be called as the generalized NS equations containing a free parameter. In the present case, such a parameter is imposed by numerics due to preconditioning. On the other hand, such generalized NS equations arise in other contexts such as in the simulation of the fluid saturated variable porous media flows represented by the Brinkman-Forchheimer-Darcy equation. In such cases, the free parameter appearing in the generalized NS equations is imposed by physics, viz., the porosity. Thus, our present investigation on the development of the Galilean invariant LB models for the preconditioned NS equations on standard lattices without cubic-velocity and parameter dependent errors also has wider implications in other contexts.

In order to first identify such truncation errors, we perform a Chapman-Enskog analysis of the preconditioned central moment LB formulation and isolate various cubic-velocity and parameter dependent errors at various moment orders. It will be seen that the anisotropy of the stress tensor depends not just on the cubic-velocity terms (like in the previous studies), but also on the preconditioning parameter . Furthermore, we will also demonstrate that even to achieve partial corrections for the GI defects on the standard lattice, the cubic velocity terms appearing in the off-diagonal components of the third-order moment equilibria need to be appropriately scaled by (e.g. ). In general, the various truncation error terms that arise due to the degeneracy of the third-order diagonal elements will be seen to have complex dependence on both the velocities and the preconditioning parameter. Once such GI defect terms are identified, new corrections are derived for the preconditioned central moment LB formulation based on the extended moment equilibria. This results in a GI central moment LB method for the preconditioned NS equations without cubic-velocity and parameter based defects on standard lattices. The present scheme is targeted towards efficient and accurate low Reynolds number steady state laminar flows by a preconditioned LB formulation without the discrete cubic velocity and parameter dependent effects via corrections to the moment equilibria. On the other hand, for high Reynolds number turbulent flow simulations, higher-order lattice based LB methods such as that presented in Chikatamarla et al. (2010) appears as one of the attractive approaches.

This paper is organized as follows. In the next section (Sec. 2), our previous central moment based preconditioned LBM with forcing terms on the D2Q9 lattice is summarized first. Section 3 performs a more refined analysis based on the Chapman-Enskog expansion and identifies various cubic-velocity and parameter dependent GI defect errors. Then, Sec. 4 derives new corrections based on the extended moment equilibria and Sec. 5 presents a GI preconditioned central moment LB method free of cubic-velocity and parameter dependent errors. Numerical results are presented in Sec. 6, which compares our numerical results for a variety of benchmark problems, including the lid-driven cavity flow, flow over a square cylinder, backward facing step flow, the Hartmann flow and the four-roll mills flow problem for the purpose of validation. In addition, convergence acceleration due to preconditioning and improvement in accuracy due to the GI corrected LB scheme are also illustrated. Finally, the main findings of our study are summarized in Sec. 7.

## Ii Preconditioned Cascaded Central Moment Lattice Boltzmann Method: Non-Galilean Invariant Formulation

In our previous work, we presented a modified cascaded central moment lattice Boltzmann method (LBM) with forcing terms for the computation of preconditioned NS equations F.Hajabdollahi and Premnath (). However, this preconditioned LBM formulation is not Galilean invariant on standard lattices. This is because it results in grid-independent cubic-velocity errors that are sensitive to the preconditioning parameter. In fact, the derivation of the precise expression for the non-GI truncation errors will be derived in the next section. It may be noted that all other prior preconditioned LB schemes are also not Galilean invariant. However, the choice of central moments here partially corrects parts of the cubic velocity defects in the off-diagonal third order moments naturally (Sec. III) and simplifies derivation of the correction terms to completely restore GI free of cubic velocity errors on standard lattice (Sec. IV). Here, we summarize our previous preconditioned central moment LB model setting the stage for further development in the following.

The preconditioned cascaded central moment LBM with forcing terms may be written as F.Hajabdollahi and Premnath ()

(2a) | |||

(2b) |

where a variable transformation is introduced to maintain second order accuracy in the presence of forcing terms. In the above, is the orthogonal transformation matrix and is the collision operator. In order to list the expressions for the collision kernel for the standard two-dimensional, nine particle velocity (D2Q9) lattice, we first define various sets of raw moments as follows on which it is based:

(3) |

The preconditioned collision kernel set for the orthogonal moment basis using the preconditioning parameter can be written as F.Hajabdollahi and Premnath ()

(4) | |||

For further details, and including the choice of the collision matrix and source raw moments , see F.Hajabdollahi and Premnath (). This scheme results in a tunable pseudo-sound speed , where , and the emergent viscosity is given by , . While this scheme is intended to simulate the preconditioned NS equations given in Eq. (1), as will be shown via a consistency analysis based on the Chapman-Enskog expansion in the next section that it leads to velocity-and preconditioning parameter-dependent non-GI truncation errors. In particular, it will be seen that the components of the non-equilibrium parts of the second order moments, which contribute to the viscous stress tensor, depends on cubic velocity truncation errors and modulated by the preconditioning parameter .

## Iii Derivation of Non-Galilean Invariant Spurious Terms in the Preconditioned Cascaded Central Moment LB Method: Chapman-Enskog Analysis

In order to facilitate the Chapman-Enskog analysis, the central moment LB formulation can be equivalently rewritten in terms of a collision process involving relaxation to a generalized equilibria in the lattice or rest frame of reference F.Hajabdollahi and Premnath (). This strategy is considered in this work to further investigate the structure of the cubic velocity non-GI truncation errors for our preconditioned LB method. In this regard, it is convenient to define the non-orthogonal transformation matrix which is the basis to obtain the orthogonal collision matrix used in the previous section and on which the subsequent analysis follows:

(5) |

where the usual bra-ket notation is used to represent the raw and column vectors in the q-dimensional space () for the D2Q9 lattice. Then, the relation between the various sets of the raw moments and their corresponding states in the velocity space can be defined via this nominal, non-orthogonal transformation matrix as

(6) |

where

are the various quantities in the velocity space, and

(7a) | |||

(7b) | |||

(7c) | |||

(7d) |

are the corresponding states in the moment space.

To facilitate the Chapman-Enskog analysis, we can rewrite the preconditioned LB model presented in Eq. (2a) and Eq. (2b) in terms of the raw moment space given in Eq. (6) as (Premnath and Banerjee (2009), F.Hajabdollahi and Premnath ())

(8) |

where the diagonal relaxation time matrix is defined as

(9) |

The preconditioned raw moments of the equilibrium distribution and source terms can be represented as

(10) |

and

(11) |

The following comments are in order here. Up to the second order moments, the above expressions coincide with those presented in our previous work F.Hajabdollahi and Premnath ()). In other words, terms in the moment equilibria are preconditioned by , while the first and second order moment terms, i.e. and are preconditioned by and , respectively. As a first new element towards a LB scheme with an improved GI, we precondition the third-order moment equilibria terms terms by (see the terms inside boxes in Eq. (10)). This partially restores GI without cubic velocity defects for the preconditioned LB model for the off-diagonal components of the third-order moments. In fact, as will be shown later in this section, in order to remove the spurious cross-velocity derivative terms appearing in the equivalent macroscopic equations of our preconditioned LB scheme (e.g. and ), such a scaling of the cubic velocity terms in the third order moment equilibria is essential. Then, applying the standard Chapman-Enskog multiscale expansion to Eq. (8), i.e.

(12) | |||||

(13) |

where is small bookkeeping perturbation parameter, and also using a Taylor expansion to simplify the streaming operator in Eq. (8), i.e.

(14) |

After converting all the resulting terms into the moment space using Eq. (6), we get the following moment equations at consecutive order in :

(15a) | |||

(15b) | |||

(15c) |

where . The relevant components of the first-order equations Eq. (15b), i.e. up to the second order in moment space needed for deriving the preconditioned macroscopic hydrodynamics equations are given as

(16a) | |||

(16b) | |||

(16c) | |||

(16d) | |||

(16e) | |||

(16f) |

Similarly, the leading order moment equations at can be obtained from Eq. (15c) as

(17a) | |||

(17b) | |||

(17c) |

In the above equations, the second-order, non-equilibrium moments , and (corresponding to, , and , respectively) are unknowns. Ideally, they should only be related to the strain rate tensor components to recover the correct physics related to the viscous stress. However, as will be should below, on the standard D2Q9 lattice there will be non-GI contributions dependent on the preconditioning parameter . In what follows, , and will be obtained from Eq. (16d), Eq. (16e) and Eq. (16f), respectively.

Now, from Eq. (16d), the non-equilibrium moment can be written as

(18) |

In order to simplify Eq. (18) further, one needs to obtain expressions, in particular, for , , and . It follows from Eq. (16b) that

(19) |

Rearranging as

Using Eq. (19) and Eq. (16a) to replace the time derivative in the first and second terms respectively, on the right hand side of the above equation, we get.

(20) | |||||

Similarly, we may write

(21) | |||||

Thus, the time derivative can be replaced with the spatial derivative. Also , it readily follows that

(22a) | |||

(22b) |

Rearranging Eq. (20) and simplifying it further by retaining all cubic velocity terms and neglecting all others higher order terms in velocity (e.g. fifth order and higher) we get

(23) |

Similarly, it follows from Eq. (21) that

(24) |

Now, to obtain an expression for , we group all the higher order terms given in Eqs. (22a), (22b), (23) and (24). It follows that owing to the choice of the off-diagonal third-order equilibrium moments with the cubic velocity terms scaled by (i.e. at the outset following Eq. (9) earlier, all the cross-derivative spurious terms, i.e. and cancel. Then, simplifying the grouping of all the remaining higher order terms in Eq. (22a), Eq. (22b), Eq. (23) and Eq. (24) and retaining all cubic velocity terms and neglecting terms of negligible higher orders and after considerable rearrangement, we get