Improving the low Mach number steady state convergence of the cascaded lattice Boltzmann method by preconditioning

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Improving the low Mach number steady state convergence of the cascaded lattice Boltzmann method by preconditioning Farzaneh Hajabdollahi, Kannan N. Premnath ∗ Department of Mechanical Engineering, University of Colorado Denver, 1200 Larimer street, Denver, CO 80217, USA

article

info

Article history: Available online xxxx Keywords: Lattice Boltzmann method Central moments Preconditioning Convergence acceleration

abstract Cascaded lattice Boltzmann method (LBM) involves the use of central moments in a multiple relaxation time formulation in prescribing the collision step. When the goal is to simulate low Mach number stationary flows, the greater the disparity between the magnitude of fluid speed and the sound speed, the higher is the numerical stiffness, which results in relatively large number of time steps for convergence of the LBM. One way to improve the steady state convergence of the scheme is to precondition the cascaded LBM, which reduces the disparities between the characteristic speeds or equivalently those between the eigenvalues of the system. In this paper we present a new preconditioned, two-dimensional cascaded LBM, where a preconditioning parameter is introduced into the equilibrium moments as well as the forcing terms. Particular focus is given to preconditioning differently the moments due to forcing at first and second orders so as to avoid any spurious effects in the emergent macroscopic equations. A Chapman–Enskog analysis performed on this approach shows consistency to the preconditioned Navier–Stokes equations. This modified central moment based scheme is then validated for accuracy by comparison against prior analytical or numerical results for certain benchmark problems, including those involving spatially variable body forces. Finally, significant steady state convergence acceleration of the preconditioned cascaded LBM is demonstrated for a set of characteristic parameters. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Over the last two decades, the lattice Boltzmann method has attracted considerable attention as a novel method for computational fluid dynamics (CFD) [1–4]. It has been successfully applied to a variety of problems, such as multiphase flows [5,6], turbulence [7,8], micro flows [9], multicomponent flows [10], complex hydrodynamic and particulate flows [11]. Its success stem from its inherent parallelism, ease of implementation of boundary conditions and natural framework to incorporate kinetic models for complex flows. During the last decade or so, the LBM has undergone a number of refinements, especially with the development of more sophisticated collision models as well as in adopting advanced strategies from classical computational fluid dynamics (CFD) to further improve its numerical stability, accuracy and efficiency. Like other explicit time-marching methods in classical CFD, the LBM suffers from slow convergence rate at low Mach numbers (Ma). In fact, this difficulty in solving the LBM to steady state at low Ma is due to the relatively large disparity between the acoustic wave speed and the convection speed of the fluid. In other words, low Ma simulations are subjected to high eigenvalue stiffness, or, equivalently, larger condition numbers which can cause inefficiencies. Therefore to maintain numerical stability, relatively small values of the time step to march to steady state are required. Moreover, to limit the

∗

Corresponding author. E-mail addresses: [email protected] (F. Hajabdollahi), [email protected] (K.N. Premnath).

http://dx.doi.org/10.1016/j.camwa.2016.12.034 0898-1221/© 2017 Elsevier Ltd. All rights reserved.

2

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

compressibility errors in the LBM so as to improve accuracy, the Ma has to be kept small in the simulations which can lead to a slower convergence rate. To accelerate lattice Boltzmann (LB) simulations various different approaches that exploit its inherently data-parallel nature have been presented in the literature. For example a nonlinear multigrid (MG) solution approach using a direct discretization of the stationary Boltzmann equation by an implicit second-order finite difference scheme was introduced by Tolke et al. [12]. This was redesigned and MG formulations that exploit the stream-and-collide nature of the LB algorithms were developed by Mavriplis [13] and Patil et al. [14]. In order to alleviate the numerical stiffness at low Mach numbers for compressible and incompressible Navier–Stokes (NS) equations, preconditioning techniques were developed by many authors in the context of classical CFD [15–25]. The first use of preconditioning in the LBM using a single relaxation time formulation for steady incompressible flows was proposed by Guo et al. [26]. They modified the equilibrium distribution function in such a way that the preconditioned NS equations were obtained by using a Chapman–Enskog expansion. They also showed that the steady state macroscopic equations derived from this LB model are the same as those obtained from the standard LB model, but with an improved eigenvalues of the system. Izquierdo and Fueyo [27] applied this idea to the generalized lattice Boltzmann equation (GLBE) based on multiple relaxation times. They also obtained optimal values of the preconditioning parameter [28]. A preconditioned GLBE with forcing term including an extension for a lattice kinetic scheme based on a vector distribution function for fast simulation of magnetohydrodynamics was presented by Premnath et al. [29]. The use of the cascaded formulation for the collision step of the mesoscopic based LBM represents one of its more recent developments. Introduced by Geier et al. [30], it is based on the use of central moments to construct the collision operator. The equilibrium moments at successively higher orders represent a cascaded structure as they are constructed for a moving frame of reference, which were shown to be equivalent to considering a generalized local equilibria in the rest or lattice frame of reference by Asinari [31]. Premnath and Banerjee [32] constructed forcing terms based on the central moments using a variable transformation to maintain second order accuracy within this approach, and performed a Chapman–Enskog analysis to demonstrate its consistency to the NS equations. As the central moments are obtained by shifting the particle velocity with the local fluid velocity, it provides a natural setting to maintain Galilean invariance of the chosen set of moments, as well as enhancing the numerical stability by tuning its relaxation times. Recently, a detailed study of the numerical properties of the cascaded LBM and its comparison with various other collision models was performed by Ning et al. [33]. In this work, we develop a simple and efficient preconditioned formulation of the cascaded LBM with forcing terms to accelerate its convergence to steady state at low Mach numbers. This is achieved by modifying the equilibrium moments and the source moments representing the forcing terms by introducing a preconditioning parameter. Special care is taken in constructing the preconditioned source moments of various orders to incorporate the general case of spatially and/or temporally varying body forces, which are important in practical applications, without introducing any spurious lattice effects into the resulting preconditioned NS equations. This will be theoretically demonstrated by a Chapman–Enskog analysis and then verified numerically. This article is organized as follows: First, we present a derivation of the preconditioned cascaded lattice Boltzmann equation (LBE) with forcing term in Section 2. In Section 3, we perform a Chapman–Enskog analysis of this new preconditioned cascaded LBE and demonstrate its consistency to the NS equations without any spurious effects in the presence of variable body forces. Section 4 deals with the implementation of our preconditioned cascaded formulation, where the numerical results obtained from computational experiments are compared with the standard benchmark problems for the purpose of validation and then illustrates its convergence acceleration. Conclusions and outlook of this work are given in Section 5. 2. Preconditioned cascaded lattice Boltzmann method In general, to accelerate the steady state convergence in simulations of compressible, low Mach number flows preconditioning techniques have been introduced. In order to reduce the eigenvalue stiffness and accelerate convergence, one can introduce a preconditioning parameter so as to reduce the condition number while maintaining steady state solution accuracy [23,24]. In this regard, for the purpose of illustration, a standard two-dimensional, nine velocity (D2Q9) lattice is employed, while the overall procedure is also equally applicable to various other models. As the cascaded LBM is a moment n based approach, we start with a set of the following nine nonorthogonal basis vectors obtained using the monomials em α x eα y in an ascending order:

80 = (1, 1, 1, 1, 1, 1, 1, 1, 1)Ď , 81 = ex = (0, 1, 0, −1, 0, 1, −1, −1, 1)Ď , 82 = ey = (0, 0, 1, 0, −1, 1, 1, −1, −1)Ď , 83 = ex ◦ ex + ey ◦ ey , 84 = ex ◦ ex − ey ◦ ey , 85 = ex ◦ ey , 86 = ex ◦ ex ◦ ey , 87 = ex ◦ ey ◦ ey , 88 = ex ◦ ex ◦ ey ◦ ey ,

(1)

where Ď is the transpose operator. In the above, for any two q-dimensional vectors (q = 9 here) a and b, we define the elementwise vector multiplication (Hadamard product) by a ◦ b. That is, (a ◦ b)α = aα bα , where α represents a component of the vector and the implicit summation convention is not assumed here. To facilitate the presentation in the following, we

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

3

q

also define a standard scalar inner product of any two such vectors as ⟨a, b⟩. That is, ⟨a, b⟩ = α=0 aα bα . By applying the Gram–Schmidt orthogonalization method they result into the following equivalent set of orthogonal basis vectors: K0 = 80 ,

K1 = 81 ,

K3 = 383 − 480 ,

K2 = 82 ,

K4 = 84 ,

K6 = −386 + 282 ,

K5 = 85 ,

K7 = −387 + 281 ,

K8 = 988 − 683 + 480 .

Collecting them as an orthogonal matrix K, we get K = [K0 , K1 , K2 , K3 , K4 , K5 , K6 , K7 , K8 ] .

(2)

See [30–32], which enumerates the details of the matrix K. Now we define the continuous equilibrium central moments as [30]:

xMm yn = Π



∞



∞

f M (ξx − ux )m (ξy − uy )n dξx dξy

(3)

−∞

−∞

where f M is the local Maxwell–Boltzmann distribution in the continuous particle velocity space ξ = (ξx , ξy ) and is given  

by: f M ≡ f M (ρ, u, ξ) =

ρ

(ξ−u)2

. Here, ρ is the density and u = (ux , uy ) is the macroscopic fluid velocity. By 2cs2 M  evaluating Πxm yn in the increasing order of moments for the D2Q9 lattice we obtain 2π cs2

exp −

0M = ρ, Π

xM = 0, Π

yM = 0, Π

xxM = cs2 ρ, Π

M xxy Π = 0,

M xyy Π = 0,

M xxyy Π = cs4 ρ.

M yy Π = cs2 ρ,

M xy Π = 0,

Similarly, we define the central moments of the sources of order (m + n) due to a force field F = (Fx , Fy ) as

xFm yn = Γ



∞

−∞



∞

1f F (ξx − ux )m (ξy − uy )n dξx dξy

(4)

−∞

where 1f F is the change in the distribution function due to force fields. Again we can evaluate Eq. (4) as

0F = 0, Γ F xxy Γ = cs2 Fy ,

xF = Fx , Γ

yF = Fy , Γ

F xyy Γ = cs2 Fx ,

xxF = 0, Γ

yyF = 0, Γ

xyF = 0, Γ

F xxyy Γ = 0.

Based on the above continuous central moments and using the trapezoidal rule for evaluating the source term to retain second order accuracy, the elements of the cascaded LBE can be formulated as follows [30,32]: 1

[Sα (x, t ) + Sα (x + eα , t + 1)] . 2 In the above equation, the collision term can be modeled by fα (x + eα , t + 1) = fα (x, t ) + Ωαc (x, t ) +

Ωαc ≡ Ωαc (f, g) = (K · g)α ,

(5)

(6)

where f = (f0 , f1 , f2 , . . . , f8 )Ď is the vector of distribution functions and  g = ( g0 , g1 , g2 , . . . , g8 )Ď is the vector of unknown collision kernel to be determined later. It should be noticed the discrete form of the source term Sα in the cascaded LBE represents the influence of the force field (Fx , Fy ) in the velocity space as: S = (S0 , S1 , S2 , . . . , S8 )Ď . As Eq. (5) is semiimplicit, by using the standard variable transformation by He et al. [34,35], and, in particular, by He et al. [36] f α = fα − 21 Sα , so that the implicitness is removed and a second-order accuracy is maintained, we get f α (x + eα , t + 1) = f α (x, t ) + Ωαc (x, t ) + Sα (x, t ).

(7)

In the above transformation, if necessary, one may also include the collision term, as we had already shown in the context of the cascaded formulation in our earlier work [32], which could be useful for problems such as those involving multicomponent flows. However, the above simpler version of the transformation with appropriate corrections for the higher order source moments is adequate for problems considered in the present study, as will be confirmed through a Chapman–Enskog analysis later. Such an approach is similar to the well-known method of Guo et al. [37] for incorporating body forces. Alternatively, as shown recently by Dellar [38], one can derive the LBE with forcing term using a Strang splitting approach [39]. In order to determine the structure of the cascaded collision operator and the source terms in the presence of general spatially and/or temporally variable body forces, we define the following set of central moments of order (m + n), respectively, as

   κˆ xm yn fα κˆ xeqm yn   fαeq  m n   S  (eαx − ux ) (eαy − uy ) . σˆ xm yn  = α α f¯α κˆ¯ xm yn 

(8)

4

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

By equating the discrete central moments for both the distribution functions and source terms with the corresponding  Mm n  Fm n , we get continuous central moments  κxeqm yn = Π σxm yn = Γ x y x y

 κ0eq = ρ,

 κxeq = 0,

eq  κxxy = 0,  σ0 = 0,

eq  κxyy = 0,  σx = Fx ,

 σxxy = 0,

 σxyy = 0,

 κyeq = 0,

eq  κxx = cs2 ρ,

eq  κxxyy = cs4 ρ.  σy = Fy ,  σxx = 0,

eq  κyy = cs2 ρ,

eq  κxy = 0,

(9)

 σyy = 0,

 σxy = 0,

 σxxyy = 0.

(10)

Here, we apply the effect of the variable body forces only on the first order central moments of the sources to obtain consistent macroscopic equations [32]. Similarly, for the transformed central moments we have 1 eq  κ x = − Fx ,

eq  κ 0 = ρ,

1 eq  κ y = − Fy ,

2

2

2

2

c eq  κ xxy = − s Fy ,

c eq  κ xyy = − s Fx ,

eq  κ xx = cs2 ρ,

eq  κ yy = cs2 ρ,

eq  κ xy = 0,

eq

 κ xxyy = cs4 ρ. (11) 2 2 On the other hand, the actual calculations in the cascaded formulation are carried out in the terms of the raw moments. Hence, we define the following set of raw moments (designated here with the prime (′ ) symbol), which will be used later: 

κˆ x′m yn



   eq′   fαeq κˆ xm yn  fα  m n   S  eαx eαy . σˆ ′m n  = α  x y  α f¯α ′ κ¯ˆ xm yn

(12)

We can readily convert the central moments into a combination of the raw moments by using the binomial theorem. Now, in order to accelerate its steady state convergence, we are in a position to apply a preconditioning parameter, designated here as γ , for the various elements of the cascaded LBE, such as the equilibria and the source terms written in terms of the raw moments. Since terms only up to the second order moments explicitly influence the consistency to the macroscopic equations, we limit the application of the preconditioning parameter to only those orders. The actual form of the preconditioned raw moments for various orders is obtained by applying the Chapman–Enskog analysis, which will be discussed in the next section. Here, we first summarize the following set of raw moments for the source terms modified by preconditioning:

 σ0′ = 0,  σxx′ =

Fx

 σx′ =

2Fx ux

γ2

,

γ

,

 σyy′ =

 σy′ = 2Fy uy

γ2

′  σxxy = Fy u2x + 2Fx ux uy ,

Fy

γ

,

,

 σxy′ =

Fx uy + Fy ux

γ2

,

′  σxyy = Fx u2y + 2Fy uy ux ,

′  σxxyy = 2Fx ux u2y + 2Fy uy u2x .

Notice the presence of γ 2 in the second order terms above, which will be elaborated in the next section. Similarly, the sβ = ⟨Kβ , S⟩ are obtained as preconditioned source moments projected to the orthogonal moment space m

s0 = 0, m

s1 = m

Fx

γ

,

s2 = m

Fy

γ

,

(Fx uy + Fy ux ) , γ2     1 1 s6 = s7 = m − 3u2x Fy − 6Fx ux uy , m − 3u2y Fx − 6Fy uy ux , γ γ      2 2 s 2 2  m8 = 3 6uy − 2 Fx ux + 6ux − 2 Fy uy . γ γ s3 = m

6F · u

γ2

,

s4 = m

2(Fx ux − Fy uy )

γ2

,

s5 = m

s1 , m s2 , . . . , m s8 )Ď . Now, inverting the Equivalently, the above equations can be written in a matrix form as K · S = ( ms0 , m above equation, we can obtain explicit expressions for Sα in terms of the general variable body force F and fluid velocity u in the particle velocity space as S0 =

1

 s8 , − ms3 + m

9  1  s s3 + 9m s4 + 6m s7 − 2m s8 , 1 − m S1 = 6m 36

(13a) (13b)

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

S2 = S3 = S4 = S5 = S6 = S7 = S8 =

)

–

 1  s 2 − m s3 − 9m s4 + 6m s6 − 2m s8 , 6m 36  1  s1 − m s3 + 9m s4 − 6m s7 − 2m s8 , −6 m 36  1  s2 − m s3 − 9m s4 − 6m s6 − 2m s8 , −6 m 36  1  s  1 + 6m s2 + 2m s3 + 9m s5 − 3m s6 − 3m s7 + m s8 , 6m 36  1  s1 + 6m s2 + 2m s3 − 9m s5 − 3m s6 + 3m s7 + m s8 , −6 m 36  1  s1 − 6m s2 + 2m s3 + 9m s5 + 3m s6 + 3m s7 + m s8 , −6 m 36  1  s  1 − 6m s2 + 2m s3 − 9m s5 + 3m s6 − 3m s7 + m s8 . 6m 36

5

(13c) (13d) (13e) (13f) (13g) (13h) (13i)

Finally, to determine the structure of the preconditioned cascaded collision operator in the presence of forcing terms we start from the lowest-order nonconservative postcollision central moments and successively set them equal to their corresponding equilibrium states. In this regard, we need to carefully precondition the equilibria of the resulting raw moments such that the emergent macroscopic equations lead to the preconditioned NS equations (see next section for details). In particular, we need to precondition terms involving the quadratic products of the fluid velocity as they relate to the convective terms in the emergent equations. Thus, the preconditioned collision kernel in the presence of the forcing terms is introduced as follows:

 g0 =  g1 =  g2 = 0,   ρ(u2x + u2y ) ω3 2 1 ′ ′ ′ ′  g3 = κ xx +  κ yy ) − ( ρ+ − ( σxx +  σyy ) , 12 3 γ 2   1 ′ ω4 ρ(u2x − u2y ) ′ ′  κ xx −  κ yy ) − ( − ( σxx −  σyy′ ) , g4 = 4 γ 2   1 ′ ω 5 ρ ux uy  ′  − κ xy −  σxy , g5 = 4 γ 2   ω6 1 1 ′ ′ ′ 2     g6 = 2ρ ux uy + κ xxy − 2ux κ xy − uy κ xx −  σxxy − uy (3 g3 +  g4 ) − 2ux g5 , 4 2 2   ω7 1 1 ′ ′ ′  g7 = 2ρ ux u2y +  κ xyy − 2uy κ xy − ux κ yy −  σxyy − ux (3 g3 −  g4 ) − 2uy g5 , 4 2 2   ′ ω8 1 ′ ′ ′ ′  κ xxyy − 2ux κ xyy − 2uy κ xxy + u2x κ yy + u2y κ xx g8 = ρ + 3ρ u2x u2y −  4 9   1 1 1 ′ ′ + 4ux uy κ xy −  σxxyy − 2 g3 − u2y (3 g3 +  g4 ) − u2x (3 g3 −  g4 ) − 4ux uy g5 − 2uy g6 − 2ux g7 . 2

2

2

(14a) (14b)

(14c)

(14d) (14e) (14f)

(14g)

In the above equations, we limited preconditioning by applying the preconditioning parameter γ only to the second order collision kernel terms  g3 ,  g4 and  g5 as they are just adequate to recover preconditioned NS equations and higher order terms are not preconditioned out of stability and efficiency considerations. The preconditioned cascaded lattice Boltzmann equation with forcing term can then be presented in terms of following collision and streaming steps, respectively:

 f α (x, t ) = f α (x, t ) + ΩαC + Sα (x, t )  f (x + e , t + 1) = f (x, t ) α

α

(15) (16)

α

where, ΩαC is obtained using Eq. (6) and Eqs. (14a)–(14g). Finally, the explicit form of the post-collision distribution function in the velocity space can be written as follows:

 f 0 = f 0 + [ g0 − 4( g3 −  g8 )] + S0 ,  g1 −  g3 +  g4 + 2( g7 −  g8 )] + S1 , f 1 = f 1 + [ g0 +   g + g − g − g + 2( g − g )] + S , f = f + [ 2

2

0

2

3

4

6

8

2

 f 3 = f 3 + [ g0 −  g1 −  g3 +  g4 − 2( g7 +  g8 )] + S3 ,

6

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

 f 4 = f 4 + [ g0 −  g2 −  g3 −  g4 − 2( g6 +  g8 )] + S4 ,  f 5 = f 5 + [ g0 +  g1 +  g2 + 2 g3 +  g5 −  g6 −  g7 +  g8 ] + S5 ,  f = f + [ g − g + g + 2 g − g − g + g + g ]+S , 6

0

6

1

2

3

5

6

7

8

6

 f 7 = f 7 + [ g0 −  g1 −  g2 + 2 g3 +  g5 +  g6 +  g7 +  g8 ] + S7 ,  f = f + [ g + g − g + 2 g − g + g − g + g ]+S . 8

0

8

1

2

3

5

6

7

8

8

3. Chapman–Enskog analysis of the preconditioned cascaded lattice Boltzmann method In this section, in order to confirm that the preconditioned cascaded LBM with forcing term recovers the preconditioned NS equations, a Chapman–Enskog analysis will now be performed. The central moment based LBM formulation can be equivalently rewritten in terms of a process involving the relaxation to a generalized equilibria in the lattice or rest frame of reference which facilitates analysis. To establish consistency to the macroscopic equations, it is sufficient to retain only terms up to second order in the Mach number in such an equivalent formulation in the lattice frame of Refs. [31,32], which is adopted in the present work. First, we define the following nominal, non-orthogonal transformation matrix [32]: T = [80 , 81 , 82 , 83 , 84 , 85 , 86 , 87 , 88 ] ,

(17)

where the nominal basis vectors 8β are obtained from Eq. (1). It is convenient in calculations to obtain the raw moments and apply the binomial theorem to compute them from central moments. In the following, the superscript ‘‘prime’’ symbol is used to display raw moments (see Eq. (12)). Conversely, the central moments (see Eq. (8)) can be obtained from their corresponding raw moments as follows: ′  κ xx =  κ xx − ρ u2x + Fx ux , ′  κ yy =  κ yy − ρ u2y + Fy uy ,

1 ′  κ xy =  κ xy − ρ ux uy + (Fx uy + Fy ux ), 2

1 ′ ′ ′  κ xxy =  κ xxy − 2ux κ xy − uy κ xx + 2ρ u2x uy − Fy u2x − Fx ux uy , 2

1 ′ ′ ′  κ xyy =  κ xyy − 2uy κ xy − ux κ yy + 2ρ ux u2y − Fx u2y − Fy uy ux , 2

′ ′ ′ ′ ′ ′  κ xxyy =  κ xxyy − 2ux κ xyy − 2uy κ xxy + u2x κ yy + u2y κ xx + 4ux uy κ xy − 3ρ u2x u2y + Fx ux u2y + Fy uy u2x .

The preconditioned raw moments of the equilibrium distribution and source terms can be written as follows: ′

′

 κ0eq = ρ,

ρ ux uy , γ

′

eq  κxy = ′

eq =  κxxyy

1 9

′

 κxeq = ρ ux , ′

eq  κxxy =

′

 κyeq = ρ uy , 1 3

ρ uy + ρ u2x uy ,

eq  κxx = ′

eq  κxyy =

1 3 1 3

ρ+

ρ u2x , γ

′

eq  κyy =

ρ ux + ρ ux u2y ,

1 3

ρ+

ρ u2y γ

, (18)

1

ρ + ρ(u2x + u2y ) + ρ u2x u2y , 3

and

 σ0′ = 0,  σxy′ =

 σx′ =

Fx uy + Fy ux

γ2

Fx

γ ,

,

 σy′ =

Fy

γ

,

 σxx′ =

′  σxxy = Fy u2x + 2Fx ux uy ,

2Fx ux

γ

2

,

 σyy′ =

2Fy uy

γ2

,

(19)

′  σxyy = Fx u2y + 2Fy uy ux ,

′  σxxyy = 2(Fx ux u2y + Fy uy u2x ),

respectively. As can be seen, only terms up to the second order in the components of the raw moments are preconditioned by the parameter γ , since only those are expected to contribute to the corresponding macroscopic equations. The higher order moments are left un-preconditioned out of numerical stability considerations as well as to reduce its effect on introducing additional floating point operations. Clearly, ui uj terms are preconditioned by γ , while Fi and Fi uj by γ and γ 2 , respectively, which are crucial to recover the required emergent behavior.

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

7

According to Eq. (6), and using Eq. (15), the cascaded collision term can be written as follows:

 (K · g) = (f − f) + S.

(20)

By applying the central moment operator ⟨(ex − ux ) (ey − uy ) , ·⟩ on both sides of the above equation we get m

n

  ⟨(ex − ux )m (ey − uy )n , Kβ ⟩ gβ = ( κ xm yn −  κ xm yn ) +  σxm yn .

(21)

β

eq

For the specific case when the post-collision state is in the ‘‘equilibrium state’’ we can set  κ xm yn =  κ xm yn ,  σxm yn = 0, for ∗ which we represent the collision kernel as  g = g . For this special case, the collision operator, Eq. (20), can also be written eq ∗ as K · g = f − f = feq − f − 21 S, which can be inverted to yield



  1 ∗  g = K−1 feq − f − S .

(22)

2

Now by using a relaxation-time matrix 3, we propose to ‘‘over-relax’’ the above system as [32] ∗  g = 3 g ,

(23)

and by combining Eqs. (22) and (23) we get

 ∗ f = f + K · g + S = f + K3 g +S   1 −1 eq = f + K3K f − f − S + S.

(24)

2

Let 3∗ = K3K−1 . Hence,

    1  f = f + 3∗ feq − f + I − 3∗ S

(25)

2

where I is the identity matrix. For convenience, the relation between the raw moments and their corresponding states in the velocity space is defined using the transformation matrix T as [32]

 f = Tf,

eq  f = Tfeq ,

 f = Tf,

 S = TS ,

(26)

where

   ′ ′ ′ ′ Ď ′ ′ ′ ′ ′ ′ ′      Ď f = f 0, f 1, f 2, . . . , f 8 =  κ 0,  κ x,  κ y,  κ xx +  κ yy ,  κ xx −  κ yy ,  κ xy ,  κ xxy ,  κ xyy ,  κ xxyy , Ď  Ď  ′ ′ ′ ′ ′ ′ ′  f=  f0 ,  f1 , f2 , . . . , f8 =  κ0 , κx , κy , κxx +  κyy′ , κxx′ −  κyy′ , κxy′ , κxxy , κxyy , κxxyy ,  Ď   ′ ′ ′ eq eq eq eq eq Ď eq ′ eq ′ eq ′ eq ′ eq ′ eq ′ eq ′ eq ′  f =  f0 ,  f1 , f2 , . . . ,  f8 =  κ0eq , κxeq , κyeq , κxx + κyy , κxx − κyy , κxy , κxxy , κxyy , κxxyy ,  Ď  ′ ′ ′ ′  Ď ′ ′ ′  S=  S0 ,  S1 ,  S2 , . . . ,  S8 =  σ0 ,  σx ,  σy ,  σxx +  σyy′ ,  σxx′ −  σyy′ ,  σxy′ ,  σxxy , σxyy , σxxyy . Now, we can rewrite Eq. (25) in terms of the raw moment space using Eq. (26) as

      1 eq  + I−  f = f + T−1 − 3  f − f 3  S , 2

(27)

where  3 is a diagonal collision relaxation time matrix which is given as

 3 = T3∗ T−1 = diag(0, 0, 0, ω3 , ω4 , ω5 , ω6 , ω7 , ω8 ).

(28)

Let us now apply a Chapman–Enskog multiscale expansion to Eq. (27). First, we introduce the expansions

 f=

∞ 

(n)

ϵ n f ,

∂t =

n=0

∞ 

ϵ n ∂tn

(29)

n=0

where ϵ is a small bookkeeping perturbation parameter. By using a Taylor expansion for the representation of the streaming operator, we get f(x + eϵ, t + ϵ) =

n  ϵn n =0

n!

(∂t + e · ∇)f(x, t ).

(30)

8

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

By substituting all the above three expansions into the preconditioned cascaded LBE, the following equations at consecutive orders of the parameter ϵ in moment space are obtained: O(ϵ 0 ) :

(0) eq  f = f ,

O(ϵ 1 ) : O(ϵ 2 ) :

(31a) (0)

Ei ∂i ) f (∂t0 + 

(1)

= − 3 f + S,   1 (0) (1) (2) Ei ∂i ) I −  f + (∂t0 +  ∂t1 3  f = − 3 f ,

(31b) (31c)

2

where  Ei = T(ei I)T−1 , i ∈ {x, y}. The components of the first-order equations in moment space, i.e. Eq. (31b) in the presence of preconditioning can be written up to the third order moments, which are relevant for deriving the macroscopic hydrodynamic equations, as

∂t0 ρ + ∇ · (ρ u) = 0,   ρ  F ρ uux x ∂t0 (ρ ux ) + ∇ · = −∂x + , γ 3 γ   ρ  F   ρ uuy y = −∂y ∂t0 ρ uy + ∇ · + , γ 3 γ       2 ρu · u 4 4 2F · u (1) ∂t0 ρ+ + ∂x ρ ux + ρ ux u2y + ∂y ρ uy + ρ u2x uy = −ω3 , f3 + 3 γ 3 3 γ2       ρ(u2x − u2y ) 2 2 2(Fx ux − Fy uy ) (1) ∂t0 + ∂x ρ ux − ρ ux u2y + ∂y − ρ uy + ρ u2x uy = −ω4 , f4 + γ 3 3 γ2       1 1 Fx uy + Fy ux ρ ux uy (1) ∂t0 + ∂x ρ uy + ρ u2x uy + ∂y ρ ux + ρ ux u2y = −ω5 f5 + . γ 3 3 γ2

(32a) (32b)

(32c)

(32d)

(32e)

(32f)

Similarly, the second-order moment equations in the presence of preconditioning can be derived from Eq. (31c), which can be written up to the first order moments that are relevant for deriving the macroscopic hydrodynamic equations, as

∂t0 ρ = 0,

(33a)



1

1 (1) (1) f3 + (1 − ω4 /2) f4 (1 − ω3 /2)



  (1) + ∂y (1 − ω5 /2) f5 = 0,

(33b)

      1 1 (1) (1) (1) ∂t1 ρ uy + ∂x (1 − ω5 /2) f5 + ∂y f3 − (1 − ω4 /2) f4 = 0. (1 − ω3 /2)

(33c)

∂t1 (ρ ux ) + ∂x

2

2

2

2

Combining Eqs. (32a)–(32c), with ϵ times Eqs. (33a)–(33c), respectively, and using ∂t = ∂t0 + ϵ∂t1 , we get the preconditioned dynamical equations for the conserved or hydrodynamic moments after setting the bookkeeping parameter ϵ to unity. That is,

∂t ρ + ∇ · (ρ u) = 0,     ρ    F ρ uux ω3(1) ω4(1) x (1) = −∂x − ∂x f3 + f4 − ∂y ω5 f5 + , ∂t (ρ ux ) + ∇ · γ 3 2 2 γ     ρ    ρ uuy ω3(1) ω4(1) Fy (1) ∂t (ρ uy ) + ∇ · = −∂y − ∂x ω5 f5 − ∂y f3 − f4 + γ 3 2 2 γ

(34a) (34b)

(34c)

where, for the purpose of clarity, we have used ω3 = (1 − ω3 /2), ω4 = (1 − ω4 /2) and ω5 = (1 − ω5 /2). The non(1) (1) (1) equilibrium, second-order raw moments  f3 ,  f4 and  f5 in Eqs. (34a)–(34c), are unknowns and can be determined from Eqs. (32d)–(32f), respectively. Thus, 1 (1)  f3 = ω3 1 (1)  f4 = ω4 1 (1)  f5 = ω5

        2 ρu · u 4 4 2F · u −∂t0 − ∂x + , ρ+ ρ ux + ρ ux u2y − ∂y ρ uy + ρ u2x uy 3 γ 3 3 γ2         ρ(u2x − u2y ) 2 2 2(Fx ux − Fy uy ) 2 2 −∂t0 − ∂x ρ ux − ρ ux uy − ∂y − ρ uy + ρ ux uy + , γ 3 3 γ2         ρ ux uy 1 1 Fx uy + Fy ux −∂t0 − ∂x ρ uy + ρ u2x uy − ∂y ρ ux + ρ ux u2y + . γ 3 3 γ2

(35a)

(35b)

(35c)

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

9

In order to simplify the above equations, when used with Eqs. (32a)–(32c), we neglect terms of O(u3 ) or higher and get ∂t0 (ρ u2x ) ≈ 2Fx ux /γ 2 , ∂t0 (ρ u2y ) ≈ 2Fy uy /γ 2 , ∂t0 (ρ ux uy ) ≈ (Fx uy + Fy ux )/γ 2 . By substituting the above relations in Eqs. (35a)–(35c), we get

  2 1  2  (1) (1) (1)    f3 ≈ − f5 ≈ − (36) ∇ · j, f4 ≈ − ∂x jx − ∂y jy , ∂x jy + ∂y jx , 3ω3 3ω4 3ω5 where jx = ρ ux , jy = ρ uy . Note that the proper choice of preconditioning the second order moments of the source terms, (1) (1) (1) which corresponds to terms such as Fi uj by γ 2 allows the correct forms for  f , f and  f without spurious lattice effects. 3

4

5

By substituting Eq. (36) into the conserved moment equations, Eqs. (34a)–(34c), we finally obtain:

∂t ρ + ∇ · j = 0,       jjx p∗ ϑ4 ϑ3 ϑ5 Fx ∂t jx + ∇ · = −∂x + ∂x (2∂x jx − ∇ · j ) + ∇ · j + ∂y (∂x jy + ∂y jx ) + , γ γ γ γ γ γ       p∗ ϑ5 ϑ4 ϑ3 Fy jjy = −∂y + ∂x (∂x jy + ∂y jx ) + ∂y (2∂y jy − ∇ · j ) + ∇ · j + , ∂t jy + ∇ · γ γ γ γ γ γ where j = (jx , jy ) and       γ 1 1 γ 1 1 γ 1 1 ϑ3 = − , ϑ4 = − , ϑ5 = − , 3 ω3 2 3 ω4 2 3 ω5 2 γ p∗ = ρ.

(37a) (37b)

(37c)

(38) 3 To maintain the isotropy of the viscous stress tensor ϑ4 = ϑ5 , as these terms determine the shear kinematic viscosity of the fluid and ϑ3 controls its bulk viscous behavior which are influenced by the preconditioning parameter γ and so is the definition of the pressure p∗ . Therefore, the preconditioned cascaded lattice Boltzmann equation is consistent with the weakly compressible preconditioned Navier–Stokes equations. Based on the above considerations, the following comments are now in order to better clarify the explicit relation between the present preconditioned cascaded LBM and the corresponding macroscopic equations so as to interpret and identify the features of the scheme and the role of the preconditioning parameter γ that enable convergence acceleration. First, we rearrange Eqs. (37a)–(37c) more compactly as follows:

∂t ρ + ∇ · (ρ u) = 0,

(39a)

γ ∂t (ρ u) + ∇ · (ρ uu) = −∇ p + ∇ · [ρν (2S − I∇ · u) + ρζ I∇ · u] + F , (39b) 1 Ď where S = 2 (∇ u + (∇ u) ) and I is the identity tensor. Also,   ζ = (γ /3) (1/ω3 − 1/2) , ν = (γ /3) 1/ωβ − 1/2 , β = 4, 5, which represent the bulk and shear viscosities, respectively, and F = (Fx , Fy ). In other words, by collecting the factor γ , ∗

which originally appears in the denominator in the convective term, pressure term, viscous term and the body form term into the single time derivative term (see Eq. (39b)), we arrive at the time-derivative preconditioning of the NS equations as the equivalence of our numerical scheme. Thus, it is clear then that the preconditioning parameter γ controls a pseudotime and by appropriately adjusting it, as discussed below, one can significantly improve convergence by sacrificing the temporal accuracy. At steady state, this free parameter γ is eliminated and Eqs. (39a) and (39b) reduce to the standard steady NS equations, and thereby the scheme would retain the correct stationary state physics. While it may appear that as usual by neglecting the time derivative term in the continuity equation to approximately satisfy the divergence free velocity field and rescaling the time variable in the above momentum equation, one may try to build a time accurate preconditioned NS solver, care needs to be exercised in such an interpretation. It may be noted that as the time derivative term after preconditioning, which is referred to as the pseudo-time derivative, becomes equal to γ times the true time derivative, this for a discrete time marching scheme implies taking larger time steps. This is also implied by the reduced pseudo sound speed with preconditioning when the same grid spacing is employed. Such an effective larger time step after preconditioning is purely dependent on the choice of the numerical parameter γ . On the other hand, the LBM is an explicit scheme and to resolve time dependent flow physics, its time step should be limited by the CFL condition. Since, in general, γ ≪ 1, the effective time step after preconditioning may be one or more orders of magnitude larger than the time step of the original scheme, which could result in a loss of temporal accuracy as the CFL condition may not be satisfied any longer for the explicit scheme after preconditioning. Hence, in general cases, our present preconditioned LBM scheme is restricted for convergence acceleration of the computation of steady state flows. To further analyze the above system (Eqs. (39a) and (39b)), as is usual in the CFD literature (e.g. [24]), and as exemplified in the case of the multiple relaxation time LBM formulation by [27], we further express it as a PDE system in vector form given by P

∂M ∂Y ∂N ∂Y ∂Y + + = R (Y ), ∂t ∂Y ∂x ∂Y ∂y

(40)

10

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

where Y = (ρ, ρ ux , ρ uy )Ď , M = (ρ ux , ρ u2x + p∗ , ρ ux uy )Ď , N = (ρ uy , ρ ux uy , ρ u2y + p∗ )Ď and R (Y ) lumps all the viscous and

√

body force terms. Here, the modified pressure in Eq. (38) can be rewritten as p∗ = γ cs2 ρ , where cs = 1/ 3 is the speed of sound of the original (unpreconditioned) cascaded LBM. The time derivative preconditioning matrix P in Eq. (40) follows as P = diag(1, γ , γ ).

To quantify the role of preconditioning in the evolution of the macroscopic conserved variables Y , it is crucial to identify the eigenvalues λi of the following matrices appearing in Eq. (40): P−1

∂M , ∂Y

P−1

∂N . ∂Y

, we have For example, evaluating the corresponding matrix in the x-direction, i.e. P−1 ∂∂M Y



0 1/(3γ ) 0

0

ux /γ 0

0 0



ux /γ

and hence the corresponding 3 eigenvalues are

  1 ∂M = (ux , ux − cs , ux + cs ) , (41) λ P−1 ∂Y γ  1/2 where  cs = ux 1 − γ + (γ cs /ux )2 is the effective or pseudo speed of sound. Notice that in the case of low-Mach number flows, which is the typical application range of the LBM, as ux → 0, the pseudo-sound speed reduces to the following  cs = cs∗ = γ 1/2 cs . This may also be obtained directly by writing our effective pressure p∗ in Eq. (38) as p∗ = (cs∗ )2 ρ and comparing it with p∗ = (γ /3)ρ , where cs2 = 1/3. From the above discussion, it then follows that for any characteristic fluid velocity U, the Mach number for the original (unpreconditioned) cascaded LBM, Ma and effective (preconditioned) Mach number for the preconditioned cascaded LBM, Ma∗ are given by Ma = U /cs and Ma∗ = U /cs∗ , respectively, and the two related by Ma∗ = γ −1/2 Ma.

(42)

Clearly, while Ma is fixed by the flow problem statement, the preconditioned Ma∗ is tunable and dependent on γ . In particular, by changing (usually decreasing) the preconditioning parameter (γ = 1 corresponds to the original scheme), the differences between the characteristic wave speeds of the pseudo-sound propagation cs∗ and the local fluid advection (e.g. ux ) can be reduced and the convergence to stationary states of the solution can be accelerated. This may also be quantified in terms of the condition number (CN), which measures the numerical stiffness in terms of the maximum disparities in the eigenvalues of the system. Thus, the condition number for the original cascaded LBM is given by CN = max λi (∂ M /∂ Y ) λj (∂ M /∂ Y ) ∼ 1/Ma,







and that for the preconditioned cascaded LBM follows as CN∗ = max λi P−1 ∂ M /∂ Y

 

   −1  λj P ∂ M /∂ Y  ∼ 1/Ma∗ .

By appropriately tuning (i.e. decreasing) γ , one can ensure Ma∗ ≫ Ma and hence CN∗ ≪ CN. Thus, our preconditioned solver can effectively remove the eigenvalue stiffness in the standard cascaded LBM that arise at low-Ma flows and can accelerate the steady state convergence. We now compare our preconditioning formalism to that carried out for the other existing LBM models. Izquierdo and Fueyo [27] modified the raw moments of equilibria by introducing the parameter γ in the multiple relaxation time formulation, while earlier, Guo et al. [26] modified the equilibrium distribution function in the velocity space to obtain a pseudo-sound speed to alleviate the eigenvalue stiffness. In contrast, in our present work, we obtain a tunable pseudosound speed by modifying the terms in the moment equilibria appearing in the collision kernel in the cascaded LBM, which is designed a priori based on central moments. Also, the earlier studies did not present a detailed consideration or analysis of preconditioning in the presence of variable body forces, which is addressed in the present investigation. The general aim, in each case, however, remains the same, which is to devise a mesoscopic scheme that is equivalent to the time-derivative preconditioning of the NS equations. To derive an efficient scheme, we explicitly limit preconditioning the moments of the equilibria and the source terms appearing in the collision kernel only up to the second order, which is adequate to recover the desired preconditioned macroscopic equations. Nevertheless, due to the cascaded structure of our formulation, the modified second order moment terms influence the third and higher order moment terms in the collision kernel. Central moments were used in the original model development primarily to improve the numerical stability. This is equivalent to considering additional terms in the higher order (third and beyond) in the generalized equilibria, when the model is rewritten in terms of raw moments. While these higher order additional terms do not play a role in adjusting the pseudo-sound speed they help to improve the numerical stability of the scheme (see e.g., [40,33]). Thus, the overall goal of our approach is to present a technique which

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

1.5

×10-3

1.5

1

0.5

0.5

–

11

×10-4

u(z)

1

)

0 –50

0

50

0 –50

0

50

z Fig. 1. Comparison of the computed velocity profiles using the preconditioned cascaded LBM (γ = 0.1) with the analytical solution for Poiseuille flow at Ma = 0.01 and Ma = 0.001, when Re = 10.

modifies the standard cascaded LBM with significantly improved convergence and stability characteristics for accurate and efficient solution of steady state flows, i.e. to solve the time-derivative preconditioned NS equations and not necessarily designed to preserve the nature of the central moments at the kinetic level after preconditioning. It may be noted that, while central moments were introduced with the additional intent of yielding a Galilean invariant model, due to the finiteness of the lattice, the resulting scheme does not fully satisfy this desired property and a posteriori additional corrections are still needed in this regard to match the flow physics prescribed by the NS equations [40]. Our present approach is in line with such a philosophy of making any additional changes necessary to the cascaded LBM to improve its numerical characteristics, such as the convergence acceleration in the present case and not necessarily to maintain its original characteristics imposed at the kinetic level. As such, work is in progress to develop a preconditioned technique for the cascaded LBM with additional corrections to obtain a fully Galilean invariant model to represent the equivalent preconditioned NS equations, i.e. for accurate and efficient simulation of steady state flows without the gridindependent cubic velocity truncation errors, which will be reported in future. 4. Computational experiments 4.1. Poiseuille flow In this section to demonstrate the performance of our proposed preconditioned cascaded LBM with forcing term, some numerical simulations are carried out. We first apply our preconditioned model to compute the steady plane Poiseuille flow where, the flow is steady and driven by a pressure gradient. In this case, a spatially and temporally constant body force Fx is applied to represent the uniform pressure gradient. In the simulations, we apply periodic boundary conditions at the inlet and outlet and the half-way bounce-back scheme is utilized for the no-slip boundary conditions at the wall. This flow problem has the well-known analytical solution given by the parabolic velocity profile u(z ) = Umax [1 − (z /L)2 ], where Umax = Fx L2 /(2ρν). Here, L, ρ and ν are the channel half-width, fluid density and kinematic viscosity, respectively. Simulations are performed at Re = 10 and two different Ma of 0.01 and 0.001 with the preconditioning parameter γ set to 0.1. To set the dimensionless parameters, the shear viscosity is adjusted by varying the relaxation times for the second order moments, while the relaxation times for other higher moments are to be set to unity for simplicity. Fig. 1 shows comparisons between the computed velocity profiles using the preconditioned cascaded LBM and the analytical solutions. Excellent agreement is seen. Next, convergence histories for various values of the preconditioning parameter γ are presented in Fig. 2. It should be pointed out that we analyze the convergence to the steady state by the second norm of the residual error of the velocity field ∥u(t + 10) − u(t )∥2 . It can be observed that significant convergence acceleration occurs as the preconditioning parameters γ are decreased. In particular, at least one order of magnitude improvement in acceleration to steady state can be seen for γ = 0.05, when compared to the cascaded LBM without preconditioning γ = 1. Furthermore, it is seen that, in general, the lower Ma cases take larger number of iterations for steady state convergence due to the higher associated eigenvalue stiffness or the condition number. The preconditioned cascaded LBM helps to significantly alleviate such numerical stiffness and evidently the benefit of preconditioning is seen to be greater at lower Ma. In order to demonstrate that preconditioning has an advantage over choosing a larger Mach number for the simulation of the same incompressible flow problem, we carried out the following numerical case study. For conditions similar to Fig. 2, i.e. considering Poiseuille flow with the same Reynolds number and grid resolution as above, we performed an error analysis by comparing the computed velocity profiles against the analytical solution for different combinations of the preconditioning parameter γ and Mach number Ma. The error between the converged computed velocity field uc and the analytical solution ua is defined under a second norm ∥uc − ua ∥2 and the results obtained are summarized in Table 1.

12

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

Fig. 2. Convergence histories of the preconditioned cascaded LBM and the standard cascaded LBM (γ = 1) for the computation of Poiseuille flow at Ma = 0.01 and Ma = 0.001, when Re = 10. Table 1 Error analysis with and without preconditioning at various Mach numbers for the simulation of Poiseuille flow. Error between the computed uc and the analytical ua velocity profiles is obtained by taking a second norm. The first row corresponds to the results with preconditioning, with the rest obtained without preconditioning. Preconditioning parameter,

Mach number, Ma

Error = ∥uc − ua ∥2

0.1 1.0 1.0 1.0

0.001 0.01 0.02 0.1

4.1500 × 10−5 4.1182 × 10−4 8.1965 × 10−4 4.0049 × 10−3

γ

It is evident that, in general, using lower Mach numbers together with preconditioning (i.e. γ = 0.1) gives much lower errors than using higher Mach numbers but without preconditioning (i.e. γ = 1.0) to simulate the same incompressible flow problem at a given Reynolds number. That is, the significant increase in the compressibility errors negates any reduction in the numerical stiffness achieved when the Mach number is increased. The major reduction in accuracy simply cannot compensate for any gain achieved in terms of efficiency if the Mach number were to be increased. Thus, to simulate and represent incompressible flows using a (weakly) compressible solver at finite Mach numbers with a chosen threshold for the desired error, using a preconditioned scheme at lower Mach numbers would be much more accurate and efficient than trying to solve it using higher Mach numbers but without preconditioning. 4.2. Hartmann flow As the next benchmark case, we consider the application of our preconditioned cascaded formulation to a linear magnetohydrodynamics problem, viz., the Hartmann flow in which a uniform pressure gradient drives the flow between two parallel plates in the presence of its interaction with an imposed uniform magnetic field Bz = B0 in the direction normal to the flow. A characteristic feature of this benchmark problem is that the flow is effectively driven by a spatially varying body force. Hence, it represents an important test problem for the preconditioned cascaded LBM where the force terms for various moments need to be carefully preconditioned in different ways so as to avoid any spurious effects in the emergent macroscopic equations involving such a variable body force. The interaction between the imposed magnetic field B0 and the fluid motion results in an induced magnetic field Bx along the channel given by Bx = (Fb L/B0 ) [sinh (Haz /L) / sinh(Ha) − z /L], where L and Fb are the channel half-width and the uniform driving force due to an imposed pressure gradient, respectively. Here, Ha is the Hartmann number, which is the ratio of the electromagnetic force to the viscous force. The interaction of this induced field and the flow results in a variable retarding Lorentz force components Fmx = B0 (dBx /dz ) and Fmz = B0 (dBx /dy). As a result, the effective variable body force components imposed in the simulation are given as F√ x = Fb + Fmx and Fz = Fmz . The analytical solution for the Hartmann flow is given by the velocity profile ux = (Fb L/B0 ) η/ν coth(Ha) [1 − cosh (Haz /L) / cosh(Ha)], where η is the magnetic resistivity given by η = B0 2 L2 /Ha2 ν . Fig. 3 presents comparisons of the computed velocity profile using the preconditioned cascaded LBM with γ = 0.1, and Ma = 0.01 and Ma = 0.001 against the analytical solution for different values of Ha. It can be noticed that the preconditioned cascaded formulation is able to reproduce the benchmark solution very well. In particular as the Ha is increased, the higher magnitudes of the Lorentz force cause significant flattering of the velocity profiles which can be seen to be reproduced by the preconditioned approach with good accuracy. Next, we analyze the influence of the preconditioning parameter on the steady state convergence of this problem. Fig. 4 presents the convergence histories for Ha = 5 and

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

13

Fig. 3. Comparison of the computed velocity profile using the preconditioned cascaded LBM (γ = 0.1) with the analytical solution for Hartmann flow for various values of Ha at Ma = 0.01 and Ma = 0.001. The lines indicate analytical results, and the symbols are the solutions obtained by the preconditioned cascaded LBM.

10-8

10-10

10-12

10-14

10-16 0 10

102

104

106

Fig. 4. Convergence histories of the preconditioned cascaded LBM and the standard cascaded LBM (γ = 1) for the computation of Hartmann flow with Ma = 0.01 and Ma = 0.001 when Ha = 5.

at two different values of Ma of 0.01 and 0.001. It can be seen that when compared to the usual cascaded LBM without preconditioning (γ = 1), the preconditioned formulation (e.g. for γ < 0.1) is able converge to steady state significantly faster, with the residual error being reduced to the machine round off error at a rate at least 10 times more rapidly. Thus the preconditioned cascaded LBM provides an effective approach to handle problems with variable body force efficiently and with good accuracy. 4.3. Lid-driven cavity flow As the last test problem, the preconditioned cascaded LBM is applied for the simulation of two-dimensional, steady flow within a square cavity driven by the motion of its top lid. For the purpose of illustration, two different Reynolds numbers of 100 and 1000 are considered in the simulation, which are resolved by computational meshes of 100 × 100 and 200 × 200, respectively. To implement the moving top wall at a velocity Up , the momentum augmented half-way bounce back scheme is used. To validate the preconditioned cascaded formulation for this complex flow problem, the computed dimensionless horizontal and vertical velocity profiles along the vertical and horizontal center lines, respectively, for Re = 400 and 1000 at Ma = 0.01 and γ = 0.05 are presented together with benchmark solutions of Ghia et al. [41] in Fig. 5. It is clear that the velocity profiles for all the cases agree very well with the prior numerical data. Finally the computed streamlines at the steady state for Re = 100 and 1000 are shown in Fig. 6. For Re = 100 a vortex is formed at lower right corners of the cavity. As the Re increases to 1000 a secondary vortex can also be clearly seen at the lower left corner, which is consistent with features observed in prior benchmark results. Let us now investigate how the steady state convergence histories are influenced by the use of preconditioning within the cascaded LBM formulation for this benchmark problem. Fig. 7 presents the convergence histories for Re = 100 and Re = 1000, where for each case, two sets of Ma of 0.01 and 0.005 are considered by varying the preconditioning parameter γ . Clearly, the use of preconditioning accelerates the steady state convergence by at least one order of magnitude. Moreover, for Re = 100, it can be seen that when compared to the usual case without preconditioning (γ = 1), the preconditioned cascaded LBM is about 27

14

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

Fig. 5. Comparison of the computed horizontal velocity u/Up and vertical velocity v/Up profiles along the geometric centerlines of the cavity using the preconditioned cascaded LBM with the benchmark results of Ghia et al. [41] (symbols) for Re = 100 and 1000 at Ma = 0.01 and γ = 0.05. 1

0.8

y

y

0.6

0.4

0.2

0

0.2

0.4

0.6

0.8

0

0

x (a) Re = 100.

0.2

0.4

0.6

0.8

1

x (b) Re = 1000.

Fig. 6. Streamlines contours computed using the preconditioned cascaded LBM for Re = 100 and Re = 1000 when γ = 0.1.

times faster for Ma = 0.01 and about 50 times faster for Ma = 0.005. The greater benefit of preconditioning at lower Ma observed numerically is consistent with its role in reducing the eigenvalue stiffness when there is greater disparities in the characteristic flow speed and the sound speed. On the other hand, as the Re is increased for a given Ma and the preconditioning parameter γ , it takes longer for steady state convergence due to reduced viscous damping effects. It is also noticed that the greater benefit of preconditioning at lower Ma mentioned above seems somewhat diminished at higher Re. It may be noted that other acceleration strategics such as the multigrid formulation has also observed similar behavior in steady state convergence at higher Re [13]. 5. Summary and conclusions Simulation of fluid flows at relatively low Mach numbers is characterized by large condition numbers due to higher disparities in flow speeds and the sound speed. This results in slow convergence of numerical schemes to steady state. Many efforts have been attempted to overcome this issue and one of the simplest and efficient among them involves preconditioning. This effectively avoids such disparities in the characteristic speeds in the system by sacrificing the temporal accuracy while preserving the accuracy of the steady state results. Cascaded lattice Boltzmann method (LBM) is one of the more recent developments in computational fluid dynamics. It is based on central moments. In particular, the equilibrium moments of different orders and the forcing terms can also be selfconsistently constructed using central moments. In this work, we focused on developing a new preconditioned cascaded formulation of the LBM to improve the steady state convergence acceleration at lower Mach numbers. Such a scheme is derived by introducing a preconditioning parameter γ into the equilibrium moments in the collision kernel and source moments representing external forces. We limit applying the preconditioning parameter up to the second order moments for both equilibria and sources out of stability and efficiency considerations, while preserving accuracy. Particular attention

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

15

Fig. 7. Convergence histories of the preconditioned cascaded LBM and the standard cascaded LBM (γ = 1) for lid-driven cavity flow with Ma = 0.01 and Ma = 0.005 for Re = 100 (top) and Re = 1000 (bottom).

is paid to preconditioning the first and second order moments of sources differently (by applying γ and γ 2 , respectively) to avoid any spurious lattice effects in the presence of spatially and/or temporally varying body forces in the resulting macroscopic equations. The consistency of the proposed preconditioned cascaded LBM with forcing terms to the equivalent time-derivative preconditioned Navier–Stokes equations is demonstrated by using the Chapman–Enskog analysis. The use of the preconditioning technique results in an adjustable pseudo-sound speed depending on γ and hence a tunable preconditioned Mach number in the cascaded formulation that effectively enables acceleration of steady state convergence. The resulting numerical scheme is first validated by comparison against prior analytical or numerical benchmark results with good accuracy for Poiseuille flow, Hartmann flow and the lid driven cavity flow. Furthermore, significant convergence acceleration was found for all the test problems considered, with at least one or more orders of magnitude improvements in efficiency over those without preconditioning. As expected, the numerical results show the benefit of preconditioning is greater at lower Mach numbers. For example an improvement by a factor of about 50 was noticed for the steady state simulation of lid driven cavity flow at a Reynolds number of 100 and a Mach number of 0.005. The present formulation is naturally applicable to other central moment approaches such as the factorized LBM or cumulant LBM [40] and can be extended to three dimensions to accelerate convergence. To further enhance the efficiency, our future work will focus on combining the preconditioning formulation within a multigrid strategy to provide a robust, local and fast approach to handle flow problems under broader range of conditions. Acknowledgment The authors F. Hajabdollahi and K.N. Premnath would like to acknowledge the travel support from US National Science Foundation (NSF) to attend ICMMES-2015, held in CSRC (www.csrc.ac.cn), Beijing, July 20–24, 2015 under the Grant CBET1549614.

16

F. Hajabdollahi, K.N. Premnath / Computers and Mathematics with Applications (

)

–

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41]

R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: theory and applications, Phys. Rep. 222 (1992) 145–197. Y. Qian, S. Succi, S. Orszag, Recent advances in lattice Boltzmann computing, Annual Rev. Comput. Phys. 3 (1996) 195–242. S. Chen, G. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. L.-S. Luo, The lattice-gas and lattice Boltzmann methods: past, present, and future, in: International Conference on Applied Computational Fluid Dymamics, Beijing, China, 2000, pp. 52–83. L.-S. Luo, Theory of the lattice Boltzmann method: Lattice Boltzmann models for nonideal gases, Phys. Rev. E 62 (2000) 4982–4996. X. He, G. Doolen, Thermodynamic foundations of kinetic theory and lattice Boltzmann models for multiphase flows, J. Stat. Phys. 107 (2002) 1572–4996. H. Chen, S. Orszag, R. Shock, S. Succi, V. Yakhot, Extended Boltzmann kinetic equation for turbulent flows, Science 301 (2003) 633–636. K. Premnath, M. Pattison, S. Banerjee, Dynamic subgrid scale modeling of turbulent flows using lattice-Boltzmann method, Physica A 388 (2009) 2640–2658. J. Zhang, Lattice Boltzmann method for microfluidics: models and applications, Microfluid. Nanofluid. 10 (2011) 1–28. P. Asinari, Semi-implicit-linearized multiple-relaxation-time formulation of lattice Boltzmann schemes for mixture modeling, Phys. Rev. E 73 (2006) 056705–24. C. Aidun, J. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439–472. J. Tolke, M. Krafczyk, E. Rank, A multigrid-solver for the discrete Boltzmann equation, J. Stat. Phys. 107 (2002) 573–591. D. Mavriplis, Multigrid solution of the steady-state lattice Boltzmann equation, Comput. Fluids 35 (2006) 793–804. D. Patil, K. Premnath, S. Banerjee, Multigrid lattice Boltzmann method for accelerated solution of elliptic equations, J. Comput. Phys. 265 (2014) 172–194. A. Godfrey, R. Walters, B. van Leer, Preconditioning for the Navier-Stokes equations with finite rate chemistry, AIAA Paper No. 0535, Reno, NV, 1993. B. Koren, B. van Leer, Analysis of preconditioning and multigrid for Eluer flows with low-subsonic regions, Adv. Comput. Math. 4 (1995) 127–144. D. Lee, B. van Leer, J. Lynn, A local Navier-Stokes preconditioner for all Mach and cell Reynolds numbers, AIAA Paper No. 2024, Washington DC, 1997. W. Lee, Local preconditioning of the Euler equations (Ph.D. thesis), University of Michigan, 1991. M. Liou, C.J. Steffen, A new flux splitting scheme, J. Comput. Phys. 107 (1993) 23–39. D.J. Mavriplis, Directional agglomeration multigrid techniques for high-Reynolds number viscous flows, J. Comput. Phys. 145 (1997) 141–165. Y.-H. Choi, C. Merkle, The application of preconditioning in viscous flows, J. Comput. Phys. 105 (1993) 207–223. E. Morgan, J. Periaux, F. Thomasset, Analysis of laminar flow over a backward facing step, in: Notes on Numerical Fluid Mechanics, vol. 9, 1984, pp. 245–267. E. Turkel, Review of preconditioning methods for fluid dynamics, Appl. Numer. Math. 12 (1993) 257–284. E. Turkel, Preconditioning techniques in computational fluid dynamics, Rev. Fluid Mech. 31 (1999) 385–416. J. Weiss, W. Smith, Preconditioning applied to variable and constant density flows, AIAA J. 33 (1995) 2050–2057. Z. Guo, T.S. Zhao, Y. Shi, Preconditioned lattice-Boltzmann method for steady flows, Phys. Rev. E 70 (2004) 066706–8. S. Izquierdo, N. Fueyo, Preconditioned Navier-Stokes schemes from the generalized lattice Boltzmann equation, Prog. Comput. Fluid Dyn. 8 (2008) 189–196. S. Izquierdo, N. Fueyo, Optimal preconditioning of lattice Boltzmann methods, J. Comput. Phys. 228 (2009) 6479–6495. K. Premnath, M. Pattison, S. Banerjee, Steady state convergence acceleration of the generalized lattice Boltzmann equation with forcing term through preconditioning, J. Comput. Phys. 228 (2009) 746–769. M. Geier, J. Greiner, F. Korvink, Cascaded digital lattice Boltzmann automata for high Reynolds number flow, Phys. Rev. E 73 (2006) 066705. P. Asinari, Generalized local equilibrium in the cascaded lattice Boltzmann method, Phys. Rev. E 78 (2008) 016701. K. Premnath, S. Banerjee, Incorporating forcing terms in cascaded lattice Boltzmann approach by method of central moments, Phys. Rev. E 80 (2009) 036702. Y. Ning, K. Premnath, D. Patil, Numerical study of the properties of the central moment lattice Boltzmann method, Internat. J. Numer. Methods Fluids 82 (2015) 59–90. X. He, X. Shan, G. Doolen, Discrete Boltzmann equation model for nonideal gases, Phys. Rev. E 57 (1998) R13–R16. X. He, S. Chen, G. Doolen, A novel thermal model of the lattice Boltzmann method in incompressible limit, J. Comput. Phys. 146 (1998) 282–300. X. He, S. Chen, R. Zhang, A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability, J. Comput. Phys. 152 (1999) 642–663. Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E 65 (2002) 046308. P. Dellar, An interpretation and derivation of the lattice Boltzmann method using Strang splitting, Comput. Math. Appl. 65 (2013) 129–141. G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 506–517. M. Geier, M. Schonherr, A. Pasquali, M. Krafczyk, The cumulant lattice Boltzmann equation in three dimensions: Theory and validation, Comput. Math. Appl. 704 (2015) 507–547. U. Ghia, K. Ghia, J. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387–411.