Alternative kinetic theory based lattice Boltzmann model for incompressible axisymmetric flows

Computers and Mathematics with Applications 72 (2016) 2751–2772

Contents lists available at ScienceDirect

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

Alternative kinetic theory based lattice Boltzmann model for incompressible axisymmetric flows Liangqi Zhang a,b , Shiliang Yang a , Zhong Zeng b , Liping Yao c , Jia Wei Chew a,d,∗ a

School of Chemical and Biomedical Engineering, Nanyang Technological University, Singapore 637459, Singapore

b

Department of Engineering Mechanics, College of Aerospace Engineering, Chongqing University, Chongqing, 400044, PR China

c

College of Engineering and Technology, Southwest University, Chongqing, 400044, PR China

d

Singapore Membrane Technology Center, Nanyang Environment and Water Research Institute, Nanyang Technological University, Singapore 637141, Singapore

article

info

Article history: Received 6 March 2016 Received in revised form 11 July 2016 Accepted 30 September 2016 Available online 24 October 2016 Keywords: Axisymmetric Boltzmann equation Kinetic theory based model Incompressible flows

abstract Based on the axisymmetric Boltzmann equation, an incompressible lattice Boltzmann model for axisymmetric flows is proposed within the framework of the kinetic theory based model developed by Guo et al. (2009). While retaining the advantages of the Guo et al. (2009) model in terms of the solid physics basis and simple source terms involving no gradient calculations, the present model further improves the numerical stability, and reduces compressibility errors and computational requirements. Armed with the assumption that the fluid density is a constant and thus the fluid pressure has no direct relation with the density, the incompressibility conditions are realized by applying the Hermite expansion. Then, the present model employs a novel way of calculating the fluid pressure which is derived from the modified second-order moment equation. Additionally, based on the regularized lattice BGK (RLBGK) model, an extra relaxation parameter pertaining to the ghost mode is introduced to enhance the numerical stability of the present model. The accuracy and applicability of the present model are verified by both the Chapman–Enskog theoretical analysis and numerical validations. It is demonstrated via well-acknowledged test cases that the present model is accurate and reliable for incompressible axisymmetric flows, and is able to effectively reduce the compressibility errors vis-à-vis the Guo et al. (2009). © 2016 Elsevier Ltd. All rights reserved.

1. Introduction The lattice Boltzmann (LB) equation, originating from the lattice gas micro-dynamics, provides us with a convenient alternative for the conventional computational method based on the direct discretization of the Navier–Stokes (N–S) equations [1–5]. The resulting LB method has advantageous features such as a simple evolution procedure, intrinsic parallel nature and the easy treatment of boundary conditions for systems with complex geometry. Besides, the LB equation can also be regarded as the discretization of the Boltzmann equation [6]. From this perspective, the LB method is claimed to be a second-order explicit finite difference approximation of the relevant macroscopic equations, with the adopted finite difference stencils determined by the relevant lattice models [7,8]. Moreover, despite the relevance between the LB method and the macroscopic N–S equations, the applications of the LB method have been successfully extended to some more

∗

Corresponding author at: School of Chemical and Biomedical Engineering, Nanyang Technological University, Singapore 637459, Singapore. E-mail address: [email protected] (J.W. Chew).

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

2752

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

complicated systems beyond the N–S solutions, such as the interfacial dynamics [9,10], flows at small scales [11–13] and non-equilibrium flows [14,15]. Axisymmetric fluid flows, and heat and mass transfer in cylindrical systems are widely encountered in engineering practices and can be directly resolved with the three-dimensional (3D) LB method [16–18]. However, in such applications, the axisymmetric conditions are neglected, otherwise the 3D axisymmetric problems can be reduced to quasi-2D ones in the meridian plane and the computation efficiency is highly enhanced. Two strategies are available for designing the axisymmetric LB models [19]. First, at the macroscopic level, the axisymmetric N–S equations in cylindrical systems are represented in a pseudo-Cartesian system with the extra terms regarded as source terms which are then incorporated into the standard LB equation. Second, the other strategy is realized at the distribution function level, including the modifications of the LB equation and the relevant equilibrium distribution functions. Harnessing the first strategy, the first axisymmetric LB model was proposed by Halliday et al. [20]. It was reported that the recovered momentum equations were incorrect with a missing term related to the radial velocity [21–23]. Such mistakes were removed by the following models proposed by Lee et al. [21,24] and Reis and Phillips [22,23], which were able to predict the accurate radial velocity distributions. The model proposed by Halliday et al. [20] achieved widespread applications, for example, in axisymmetric multiphase flows [25–27] and rotating thermal flows, whereby a hybrid method involving the second-order finite difference method was used for describing the azimuthal velocity and temperature evolution [28]. However, as stated by Zhou [29], the source terms in these earlier models were rather complex; for example, in the method of Halliday et al. [20], as many as five terms were involved in the second forcing term (which is related to the axisymmetric momentum equation), and the resulting forcing term expression contains more than ten items when the missing term (which improves the accuracy of the model) was added. In order to simplify the source terms, Chen et al. [30,31] derived an axisymmetric LB model from the vorticity stream equation, but the drawbacks of the Chen et al. [30,31] model were that the vorticity at the boundaries was hard to determine, and that a Poisson equation had to be solved at each time step, which significantly decreased the computation efficiency. As another means to simplify the source terms, Zhou [29] developed a LB model based on the standard LB equation, in which the source terms were simply the extra term in the transformed pseudo-Cartesian macroscopic equations and a centered scheme was adopted to remove the discrete lattice effect of the forcing terms. Despite the efforts in simplifying the source terms, the velocity gradient terms in the source terms cannot be removed by above models. In the axisymmetric thermal LB model proposed by Li et al. [32], an additional collision term, in which the relaxation parameter was correlated to the discrete particle velocity vector, was introduced and the temperature gradient terms were avoided in the source terms. Subsequently, an improved axisymmetric LB model was developed by the same authors [33], who showed that the gradient calculation in the source terms is realized by the non-equilibrium part of the distribution. This idea [33] was later introduced into the Zhou [34] model, in which velocity gradients in the force terms were eliminated. It should be noted that above models have all developed based on the idea from the first strategy stated above, in which the standard LB equation for the N–S equations in Cartesian systems was applied. Regarding the second strategy implemented at the distribution function level, an axisymmetric LB model was proposed by Guo et al. [35] from the axisymmetric continuous Boltzmann equation. Guo et al. [35] described the radial, axial and azimuthal velocity components in the same fashion, and the source terms naturally contained no velocity gradients. Notably, the kinetic theory origin implies a solid physics basis, which is a clear advantage. Subsequently, the numerical stability of the Guo et al. [35] model was enhanced by replacing the BGK collision operator with the multi-relaxation-time operator [36], after which the application was extended to axisymmetric thermal flows [37,38] and multiphase flows [39]. It has been demonstrated via the Chapman–Enskog analysis that the derived equations from the Guo et al. model can be transformed to the standard axisymmetric N–S equations by applying the incompressible continuity equation [35]. However, two shortcomings in the Guo et al. [35] model lies in: (i) the ideal gas state equation adopted leads to the compressible continuity equation, which in turn leads to compressibility errors, and (ii) the source terms were not expressed in a determinate way, and only one set of possible expressions was presented. To address the shortcomings, in this study, the moment equations of the distribution functions are firstly modified in accordance with the incompressibility conditions, i.e., the fluid density is assumed to be a constant such that the fluid pressure evolves independently. The Chapman–Enskog analysis in the Appendix demonstrates that the standard axisymmetric N–S equations for incompressible flows are obtained from the axisymmetric Boltzmann equation with the modified moments. Then, the expressions of the equilibrium distribution functions are derived from the relevant modified moments by applying the moment method developed by Grad [40,41]. Similarly, the source terms of the present model are determined with the moment constraints derived from the Chapman–Enskog analysis. Thereafter, the discretization in phase space for the present model is implemented via different means. Whereas the evolution equations for the two reduced distribution functions are both discretized by the D2Q9 lattice in the Guo et al. [35] model, the present model is such that the distribution function for the axial and radial velocity components evolves on the D2Q9 lattice, while the distribution function for azimuthal velocity is discretized by the D2Q4 lattice model. The overall accuracy of the present model is not only unaffected by the simplified phase space discretization, but the computation efficiency is also highly enhanced. Lastly, a new set of computing formulas for the independent macroscopic variables, i.e., the fluid velocity and pressure, emanates from the modified moment equations in the present model. Furthermore, motivated by the concept of the regularized lattice BGK (RLBGK) model [42,43], the evolution of the ghost modes in the present model is suppressed by adding an extra relaxation parameter pertaining to the non-hydrodynamic

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

2753

variable J , which is proven to be responsible for the instability of the LB evolution induced by the coupling between the hydrodynamic and non-hydrodynamic modes [44]. Hence, the present incompressible model becomes numerically more stable. The applicability and accuracy of the present model has been validated by several well-acknowledged benchmark tests, from which the accuracy of the present model is ascertained. Moreover, when a lower mesh resolution is used in the cylindrical cavity flow, the present model is more accurate than the Guo et al. [35] model particularly in the region near the top lid, the overall mass conservation of the present model is numerically verified (which validates the constant fluid density assumption), and the compressibility errors obtained are lesser than that of Guo et al. [35]. 2. Incompressible axisymmetric LB model 2.1. Review of the kinetic theory based model by Guo et al. [35] Generally, the Boltzmann equation in the cylindrical coordinates (r , θ , z ), where r, θ and z denotes respectively the radial, azimuthal and axial coordinates, is defined as

 ∂f ∂f ξθ ∂ f ∂f ξ 2 ∂f ξr ξθ ∂ f 1  + ξr + + ξz + θ − =− f − f (eq) (1) ∂t ∂r r ∂θ ∂z r ∂ξr r ∂ξθ τf   where f x0 , ξ 0 , t = f (r , θ , z , ξr , ξθ , ξz , t ) is the distribution function at position x0 = (r , θ , z ) and time t, and ξ 0 = (ξr , ξθ , ξz ) is the fluid particle velocity. Besides, the BGK collision model [45] is adopted in the present work with the single relaxation time τf , and the Maxwellian equilibrium distribution is defined as   2  ξ 0 − u0 ρ (eq) (2) exp − f = 2RT (2π RT )3/2 where ρ , u0 = (ur , uθ , uz ) and T are respectively the fluid density, velocity and temperature, and R is the gas constant. For ∂f the axisymmetric condition, i.e., ∂θ = 0, Eq. (1) is simplified as

 ∂f ∂f ∂f ξ 2 ∂f ξr ξθ ∂ f 1  + ξr + ξz + θ − =− f − f (eq) . ∂t ∂r ∂z r ∂ξr r ∂ξθ τf

(3)

Then, with the reduced distribution functions f˜ and f¯ defined as follows f˜ (x, ξ) = f¯ (x, ξ) =

 

fdξθ

(4a)

f ξθ dξθ

(4b)

where x = (r , z ) and ξ = (ξr , ξz ) are the corresponding reduced position and velocity vectors, we have

 ∂ f˜ ∂ f˜ ξr 1 ∂ φ˜ 1  ∂ f˜ + ξr + ξz + f˜ + =− f˜ − f˜ (eq) ∂t ∂r ∂z r r ∂ξr τf

(5a)

 ∂ f¯ ∂ f¯ ∂ f¯ 2ξr 1 ∂ φ¯ 1  + ξr + ξz + f¯ + =− f¯ − f¯ (eq) (5b) ∂t ∂r ∂z r r ∂ξr τf   where φ˜ (x, ξ) = f ξθ2 dξθ and φ¯ (x, ξ) = f ξθ3 dξθ . Moreover, the equilibrium reduced distribution functions are defined

as

˜ (eq)

f

  ρ (ξ − u)2 exp − = 2π RT 2RT

f¯ (eq) = uθ f˜ (eq)

˜ (eq)

φ

φ¯ (eq)



2

(6a) (6b)



˜ (eq)

= uθ + RT f   = u2θ + 3RT f¯ (eq)

(6c) (6d)

where u = (ur , uz ) is the fluid velocity pertaining to the particle velocity ξ = (ξr , ξz ) in the meridian plane. Additionally, the moment constraints for φ˜ and φ¯ are derived from the Chapman–Enskog analysis as follows



φ˜ (1) dξ = −

2pur r

(7a)

2754

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

∂ φ¯ (n) dξ = 0 ∂ξr   ∂ φ˜ (n) ξi dξ = −δir φ˜ (n) dξ ∂ξr



∂ φ˜ (n) dξ = ∂ξr



(7b)

(7c)

which lead to 1 ∂ φ˜

˜ = Φ

r ∂ξr 1 ∂ φ¯

¯ = Φ

r ∂ξr

u2θ

 =− 1+  =−

u2θ RT

−

RT

2τf ur



ξr − ur ˜ (eq)

r

r

f

(8a)

 ξr − ur ¯ (eq) +3 f .

(8b)

r

As for the numerical implementation of the kinetic theory based model, the phase space and time discretization of Eqs. (5a) and (5b) are required for completing the axisymmetric LB model. Before that, two new sets of distribution functions ξ 2ξ are introduced to remove the terms rr f˜ and r r f¯ .

 ∂g ∂g 1  ∂g + ξr + ξz =− g − g (eq) + G ∂t ∂r ∂z τf  ∂h ∂h ∂h 1  + ξr + ξz =− h − h(eq) + H ∂t ∂r ∂z τf

(9a) (9b)

˜ , h = r 2 f¯ , h(eq) = r 2 f¯ (eq) and H = −r 2 Φ ¯ . The macroscopic variables are related to g where g = r f˜ , g (eq) = r f˜ (eq) , G = −r Φ and h by ρ=

1



r

gdξ,

ρu =

1



g ξ dξ,

r

ρ uθ =

1 r2



hdξ.

(10)

Firstly, by applying the D2Q9 lattice Model [46]

 α=0 (0, 0) α = 1, 2, 3, 4 ξ α = c√(cos[(α − 1)π /2], sin[(α − 1)π /2])  2c (cos[(2α − 1)π /4], sin[(2α − 1)π /4]) α = 5, 6, 7, 8 √ with c = 3RT being the lattice speed, the discretization in the phase space for Eqs. (9a) and (9b) is realized as  ∂ gα ∂ gα ∂ gα 1  + ξαr + ξαz =− gα − gα(eq) + Gα ∂t ∂r ∂z τf  ∂ hα ∂ hα ∂ hα 1  + ξαr + ξαz =− hα − h(αeq) + Hα ∂t ∂r ∂z τf

(11)

(12a) (12b)

where gα(eq) and h(αeq) are the discrete equilibrium distribution functions

 gα(eq)

= r ρwα 1 +

ξα • u RT

 2 ξα • u u2 + − 2(RT )2 2RT 

(13a)

h(αeq) = ruθ gα(eq)

(13b)

with w0 = 4/9, w1,2,3,4 = 1/9, w5,6,7,8 = 1/36 being the weight coefficients of the D2Q9 lattice model. Besides, Gα and Hα are the discrete source terms

 Gα =

RT + u2θ −

Hα = r u2θ + 3RT



2τf RTur



ξαr − ur ˜ (eq)

r

RT

 ξαr − ur (eq) ¯ RT

fα

fα

(14a)

.

(14b)

Secondly, the trapezium rule is introduced for the time discretization, and the LB equation in the Guo et al. [35] model is obtained as follows



gˆα x + ξ α δ t , t + δ t − gˆα (x, t ) = −ω gˆα (x, t ) − gα(eq) (x, t ) + δ t 1 −

ω

Gα (x, t ) 2      ω hˆ α x + ξ α δ t , t + δ t − hˆ α (x, t ) = −ω hˆ α (x, t ) − h(αeq) (x, t ) + δ t 1 − Hα (x, t )











2

(15a) (15b)

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

2755

where gˆα and hˆ α are defined for eliminating the implicitness induced by the time discretization and related to the original distribution functions by

 δt  g Ωα + Gα ,

gˆα = gα −

Ωαg = −

2

δt 

hˆ α = hα −

2

 Ωαh + Hα ,

Ωαh = −

 1  gα − gα(eq)

(16a)

 1  hα − h(αeq)

(16b)

τf

τf

and the relaxation parameter is expressed as ω = 2δ t / 2τf + δ t and related to the kinematic viscosity by



ν=



1

ω

−

1



2



δ tRT .

(17)

Furthermore, the computing formulas for the hydrodynamic variables are derived from Eqs. (10) and (16) as

ρ=

1 r

ρ ui =

α

gˆα ,

ρ uθ =  

r r 2 + νδ t δir

α

1  r2 α

hˆ α

(18a)

  δt  2 gˆα ξα i + ρ RT + uθ δir .

(18b)

2

It is demonstrated from above derivations that the Guo et al. [35] model is directly derived from the axisymmetric Boltzmann equation, and therefore is characterized by the solid kinetic theory basis, which represents a significant advantage. The second advantage which distinguishes Guo et al. [35] from other existing axisymmetric LB models is that the velocity gradient terms in the source terms are automatically avoided without any special treatments. However, as discussed in Guo et al. [35], the incompressible continuity equation is necessary for obtaining the standard axisymmetric N–S equations, but the recovered mass conservation equation is compressible. Besides, although the source terms are critically important for obtaining the final hydrodynamic equations, only their moment constraints are determined by the Chapman–Enskog analysis and no explicit way of expressing φ˜ and φ¯ is provided. The present incompressible axisymmetric model has the potential to overcome the above-mentioned drawbacks of the Guo et al. [35] model, as described in the following. 2.2. Derivations of the present model In this section, to reduce the compressibility errors in the Guo et al. [35] kinetic theory based model, the incompressibility condition is incorporated by using the moment system put forth by Grad [40,41], and thus an incompressible LB model for axisymmetric rotating flows is developed. Besides, the moment expansion based on the Hermite tensorial polynomials provides us a means for determining the expressions of the source terms from their moment constraints. First, by assuming the fluid density is constant and the pressure is independent of density, the moments of the reduced distribution functions are given as

ρ0 =



ρ 0 uθ =

f˜ dξ =





f¯ dξ =



f˜ (eq) dξ,



f¯ (eq) dξ,

ρ0 u =



f˜ ξ dξ =

ρ 0 uθ u =



  φ˜ (eq) dξ, ρ0 u u2θ + RT =    ρ0 uθ u2θ + 3RT = φ¯ (eq) dξ ρ0 u2θ + p =



f˜ (eq) ξ dξ,

ρ0 ui uj + pδij =

f¯ (eq) ξ dξ



φ˜ (eq) ξ dξ



f˜ (eq) ξi ξj dξ

(19a) (19b) (19c) (19d)

where ρ0 is the constant fluid density. With the Chapman–Enskog analysis in Eqs. (A.1a) and (A.1b), the incompressible axisymmetric N–S equations are obtained as (see Appendix for details) ur ∂ uj + =0 ∂ xj r   ∂ ui ∂ ui u2 ∂p ∂ 2 ui µ ∂ ui µ ur ρ0 + uj − θ δir = − +µ + − 2 δir ∂t ∂ xj r ∂ xi ∂ xj ∂ xj r ∂r r     2 ∂ uθ ∂ uθ uθ ur ∂ uθ µ uθ ∂ uθ ρ0 + uj + =µ − − . ∂t ∂ xj r ∂ xj ∂ xj r r ∂r

(20a)

(20b)

(20c)

2756

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

Subsequently, the Hermite tensorial polynomials are introduced to derive the expressions of the equilibrium reduced distribution functions and the source terms. The dimensionless distribution function fˆ and particle velocity ξˆ are defined as fˆ =

RT (eq) f˜ ,

ξˆ = ξ

ρ0



√

RT .

(21) (n)

ˆ and the Hermite tensorial polynomials H (ξ) ˆ are provided as The definitions of the expansion weight function ω(ξ) i ˆ = ω(ξ)

1 2π

(n)

ˆ = Hi (ξ)

e−

2 ξˆ

(22)

2

(−1)n (n) ˆ ∇i ω(ξ) ˆ ω(ξ)

(23)

and the first few Hermite tensorial polynomials are

ˆ =1 H (0) (ξ)

(24a)

(1)

ˆ = ξˆ i Hi (ξ)

(24b)

(2) ˆ Hij (ξ) = ξˆ i ξˆ j − δij .

(24c)

The dimensionless equilibrium distribution function fˆ is expanded as

ˆ x, t ) = ω(ξ) ˆ fˆ (ξ,

∞  1 (n) (n) ˆ ai (x, t )Hi (ξ) n ! n =0

(25)

where the expansion coefficients are related to the moment constraints in Eq. (19) by (n)

M (n) =





(n)

1



(n)

ˆ dξ f˜ (eq) Hi (ξ)

(26a)

f˜ (eq) ξ dξ = ρ0 (RT )n/2 (a(n) + δ a(n−2) ).

(26b)

ai (x, t ) =

ˆ dξˆ = fˆ Hi (ξ)

ρ0

n

Thus, the expansion coefficients for the equilibrium reduced distribution functions are defined from the relevant moment equations in Eqs. (19a) and (19b) as f˜ (eq) : a(0) = 1,

(1)

ai

√ = uˆ i = ui / RT ,



f˜ (eq) =

ρ0 ξ exp − 2π RT 2RT 2

(1)

f¯ (eq) : a(0) = uθ ,

ai



p

 − 1 δij

ρ0 RT    2 u•ξ p ξ u2 (u • ξ) 1+ + + −1 −1 − RT 2RT ρ0 RT 2RT 2 (RT )2



2

= uθ uˆ i 

ρ0 uθ ξ2 f¯ (eq) = exp − 2π RT 2RT



(2)

aij = uˆ i uˆ j +

1+

u•ξ RT



.

(27a)

(27b)

Similarly, making use of the moment constraints in Eqs. (A.9), (19c) and (19d), the source terms are also derived by applying the moment expansion as

    p 2 (1) φ˜ : a(0) = u2θ + − τf RTur , ai = uˆ i u2θ + RT ρ0 r       2   u•ξ ρ0 ξ p 2 φ˜ = exp − u2θ + − ν ur + u2θ + RT 2π RT 2RT ρ0 r RT  2  ( 0 ) φ¯ : a = uθ uθ + 3RT     ρ0 ξ2 ¯φ = exp − uθ u2θ + 3RT 2π RT 2RT       1 ∂ φ˜ ρ0 ξ2 p 2 ξr u2 + RT u • ξξr ˜ = Φ =− exp − u2θ + − ν ur + θ − ur r ∂ξr 2π rRT 2RT ρ0 r RT RT RT

(28a)

(28b)

(28c)

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

¯ = Φ

1 ∂ φ¯ r ∂ξr



=−

ρ0 ξ2 exp − 2π rRT 2RT

 uθ u2θ + 3RT



 ξr RT

2757

.

(28d)

Subsequently, the discrete equilibrium distribution functions and the source terms are obtained as

   2   2 u • ξα ξα p u2 ˜fα(eq) = ρ0 wα 1 + u • ξ α + + −1 −1 − RT 2RT ρ0 RT 2RT 2 (RT )2     ρ0 wα 2 p ξαr u2θ + RT u • ξ α 2 ˜ Φα = − − ν ur uθ + + ξαr − ur r ρ0 r RT RT RT   u • ξk (eq) = ρ 0 w k uθ 1 + f¯k ′

(29a)

(29b)

(29c)

RT

¯k = − Φ

wk ρ0 uθ 

u2θ + 3RT ′

r

 ξkr RT ′

.

(29d)

We note that the equilibrium distribution function and source term for radial and axial velocity components are discretized by the D2Q9 model, since the nine discrete velocity vectors are necessary for maintaining the accuracy when evaluating the second order moments [6]. However, contrary to the Guo et al. [35] model, the discretization of f are accomplished by the following D2Q4 model

(eq)

and H

ξ k = c (cos[(k − 1)π /2], sin[(k − 1)π /2]) k = 1, 2, 3, 4 (30) √ with weight coefficients wk = 1/4 and c = 2RT ′ . The accuracy requirements are well satisfied by the four discrete velocity vectors since only the zeroth and first order moments are involved in the evolution of f . The LB equation for the present model is given as



gˆα x + ξ α δ t , t + δ t − gˆα (x, t ) = −ωg gˆα (x, t ) − gα(eq) (x, t ) + δ t 1 −





hˆ k x + ξ k δ t , t + δ t − hˆ k (x, t ) = −ωh







ωg 

Gα (x, t ) 2    ωh (eq) Hk (x, t ) hˆ k (x, t ) − hk (x, t ) + δ t 1 −





2

(31a) (31b)

where δ t is the discrete time step and related to the corresponding spatial step δ x through the lattice speed c = δ x/δ t. We note that the time and spatial steps in the present work are assumed to be dimensionless constant 1, and therefore c = 1 and RT = 1/3, RT ′ = 1/2. Thus, the relaxation parameters ωg and ωh are respectively related to the kinematic viscosity by

ν= ν=

1



ωg

3 1 2

1



1

ωh

1

−

 (32a)

2 1

−

 (32b)

2 4ω

which implies the relationship between ω, ωg and ωh of ωg = ω = 6−ωh . h The equilibrium distribution functions and the source terms are defined as

 gα(eq)

= r ρ 0 wα 1 +

Gα = ρ0 wα (eq)

hk



u2θ +

u • ξα

+

RT p

ρ0



= r ρ 0 w k uθ 1 + 2



u • ξα

2 (RT )2

2

− ν ur



r

u • ξk

ξαr

+

RT

−

u2 2RT

 +

u2θ + RT



p

ρ0 RT u • ξα

RT

RT

 −1

ξ 2α 2RT

ξαr − ur

 −1

(33a)

 (33b)

 (33c)

RT ′

Hk = r ρ0 uθ wk u2θ + 3RT ′



2

 ξkr RT ′

.

(33d)

Under the incompressibility limit, the macroscopic variables in the present model are the fluid velocity u and pressure p. The fluid velocity is calculated as

ρ 0 ui =

r r2

+ νδ t δir

  α

gˆα ξα i +

δt  2



p + ρ0 uθ δir 2



,

ρ 0 uθ =

1  r2

k

hˆ k

(34)

2758

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

which is straightforwardly derived from the modified moment equations in Eq. (19) and the definitions of gˆα and hˆ α in Eq. (16). Moreover, an alternative way of describing the fluid pressure, instead of the ideal gas state equation in Guo et al. [35], is proposed in the present model by applying the second-order moment equations



gα ξα i ξα j =

α

gα(eq) + gα(neq) ξα i ξα j = r ρ0 ui uj + rpδij − r ρ0 ν

 α





∂ uj ∂ ui + ∂ xj ∂ xi

 (35)

where gα(neq) = gα − gα(eq) denotes the non-equilibrium part of the distribution function gα . The second-order moment of gα(neq) is derived from the Chapman–Enskog analysis in the Appendix. Substituting the incompressible continuity equation in Eq. (20a) into Eq. (35), the formula for computing fluid pressure is obtained from the diagonal part of Eq. (35) as p=

1

 

2

r

1

  α

gˆα ξ α − 2ρ0 ν ur 2

 − ρ0 u

2

.

(36)

The present incompressible LB model is completed by the above formula for computing fluid pressure. Conclusions drawn from the detailed derivations of the present model include the following. First, the distinct advantages of the Guo et al. [35] model are retained in the present model, including the kinetic theory basis and the simple source terms. Second, the derived macroscopic equations from the present model accord well with the standard incompressible axisymmetric N–S equations, and therefore the present model is expected to perform better in avoiding the compressibility errors. Third, the D2Q4 lattice model is adopted for the evolution of the azimuthal velocity distribution function, which makes the present model more computationally efficient than Guo et al. [35]. Lastly, in contrast to Guo et al. [35] whereby the source terms were not explicitly derived, the source terms in the present model are systematically derived from their relevant moment constraints by applying the Hermite tensorial polynomials. 2.3. Stabilization of the present model In this work, the first three tensor Hermite polynomials, namely, 1, ξ α and ξ α ξ α − 13 I, are applied for constructing the equilibrium reduced distribution functions and the source terms, and the corresponding components compose the six polynomials, namely, 1, ξα r , ξα z , ξα2r − 31 , ξα r ξα z and, ξα2z − 13 . For the D2Q9 lattice model, each of these polynomials may be regarded as a nine-dimensional lattice vector

|1⟩ = (1, 1, 1, 1, 1, 1, 1, 1, 1)T

(37a)

|ξαr ⟩ = (0, 1, 0, −1, 0, 1, −1, −1, 1)

T

(37b)

|ξαz ⟩ = (0, 0, 1, 0, −1, 1, 1, −1, −1)    2 ξ − 1 = 1 (−1, 2, −1, 2, −1, 2, 2, 2, 2)T  αr 3 3 T

(37c) (37d)

|ξαr ξαz ⟩ = (0, 0, 0, 0, 0, 1, −1, 1, −1)T    2 ξ − 1 = 1 (−1, −1, 2, −1, 2, 2, 2, 2, 2)T .  αz 3 3

(37e) (37f)

These polynomials are orthogonal with respect to the weighted inner product ⟨p, q⟩ = α wα pα qα , which is the direct analogy of the orthogonality of the continuous tensor Hermite polynomials [40,41]. The above six orthogonal lattice vectors in Eq. (37) can form an R9 orthogonal basis by including the following ‘‘ghost vectors’’ gαhost and gαhost ξα [1,44,47,48]



   host  9  4 15  2 2 g ξ α  − ξ α  + = = (1, −2, −2, −2, −2, 4, 4, 4, 4)T α 2

9

9

 host  g ξαr = (0, −2, 0, 2, 0, 4, −4, −4, 4)T α  host  g ξαz = (0, 0, −2, 0, 2, 4, 4, −4, −4)T . α

(38a) (38b) (38c)

The ghost variables pertaining to the ‘‘ghost vectors’’ are defined as N =

 α

J =

 α

gˆα gαhost

(39)

gˆα gαhost ξ α

(40)

which are simply the higher-order moments relative to the hydrodynamic moments, and have no effect on the evolution of the hydrodynamic variables but do affect the numerical stability of the LB scheme. In particular, as claimed in Dellar [44], it

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

2759

is the coupling between the hydrodynamic and ghost modes that leads to instability of the LB algorithms (refer to Dellar [44] for detailed stability analysis of the LB method). Therefore, many efforts have been devoted to suppress the effects of these non-hydrodynamic modes (i.e., ghost modes). With the orthogonal basis defined in Eqs. (37) and (38), the BGK collision operator in Eq. (31a) is generalized as

Ωα = −wα



9 2

(neq) s3 Mrr(neq) Qα rr + 9s4 Mrz Qα rz +

9 2

(neq) s5 Mzz Qα zz + gαhost



1 4

s6 N +

3 8

s7 Jr ξα r +

3 8

s8 Jz ξα z



.

(41)

With Qα = ξ α ξ α − 31 I, and M (neq) is the corresponding non-equilibrium part of the second-order hydrodynamic mode M (neq) =

   1 (ˆgα − gα(eq) ) ξ α ξ α − I .

(42)

3

α

If all the relaxation parameters si are prescribed to be the same constant, Eq. (41) is reduced to the BGK model. Particularly, by choosing the relaxation parameters as s3 = s4 = s5 = ω,

s6 = s7 = s8 = 1

(43)

the influences from the ghost modes are completely removed, which leads to the regularized lattice BGK (RLBGK) scheme [42,43] gˆα x + ξ α δ t , t + δ t =





gα(eq)

9 1 − ωg wα



(x, t ) +



2



Qα : M (neq) + δ t 1 −

ωg  2

Gα (x, t ) .

(44)

In this work, for the purpose of enhancing the stability of the present model, an additional relaxation parameter ωJ pertaining to the ghost mode J is added to the BGK operator in Eq. (31a) as

    3 gˆα x + ξ α δ t , t + δ t − gˆα (x, t ) = −ωg gˆα (x, t ) − gα(eq) (x, t ) − (ωJ − ωg ) wα gαhost ξ α • J 8  ωg  Gα (x, t ) . + δt 1 −

(45)

2

With ωJ = 1, the non-hydrodynamic mode J would be damped completely. Since J is mainly responsible for the coupling between the hydrodynamic and non-hydrodynamic modes, as demonstrated in Dellar [44], the numerical stability of the present model is significantly enhanced by Eq. (45). 3. Numerical results In this section, four benchmark tests are adopted for validating the proposed incompressible axisymmetric LB model, namely, the Hagen–Poiseuille flow with analytical solution, the cylindrical cavity flow, the natural convection in an annulus between two coaxial vertical cylinders, and the swirling thermal flows in a cylindrical container with counter-rotating end disks. For the last two thermal tests, the LB equation describing temperature convection–diffusion proposed by Zheng et al. [37] is applied for the temperature field



(eq)

Tk x + ξ k δ t , t + δ t − Tk (x, t ) = −ωT Tk (x, t ) − Tk





  ωT  source Tk (x, t ) + δ t 1 − (x, t ) 2

(46)

which evolves on the D2Q4 lattice model defined in Eq. (30). The relaxation parameter ωT is related to the thermal diffusion coefficient α by

α=

1



2

1

ωT

−

1

 (47)

2

and the source term, as well as the equilibrium distribution, for the temperature distribution function are expressed as Tksource = ρ0 wk T ξkr (eq)

Tk



= r ρ0 wk T 1 +

(48a) u • ξk



RT ′

.

(48b)

The temperature is calculated as T =

1  r ρ0

Tk .

(49)

k

As for the boundary conditions, the non-equilibrium extrapolation scheme (NEES) [49–51] is applied for the physical boundary constraints, but not for the periodic boundary condition in the z-direction for the Hagen–Poiseuille flow. In

2760

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

particular, the symmetry boundary condition is also accurately realized by the NEES instead of the symmetry scheme in Guo et al. [35] and the other special treatments for avoiding the singularity at r = 0 [35,19,22,26,33,34,52]. The present implementation for the symmetry boundary is concisely introduced as follows. The physical constraint for the symmetry condition is defined as r =0:

∂ϕ = 0, ∂r

∀ϕ.

(50)

First, the unknown macroscopic variables at the symmetry axis are obtained from the nearest two neighboring nodes in the r-direction by using the gradient constraints. We then have the corresponding equilibrium distribution functions. Second, after the collision process at the inner fluid nodes, the non-equilibrium part of the distribution function at the boundary nodes is extrapolated from the neighboring nodes. The approximated equilibrium and non-equilibrium parts constitute the post-collision distribution function at r = 0 with a second-order accuracy. Moreover, the singularity is totally avoided since the collision and the calculation of the macroscopic variables are implemented at the inner fluid nodes. In this work, the idea is similarly applied to the temperature and velocity Dirichlet boundaries and insulated (Neumann) boundaries. 3.1. Hagen–Poiseuille flow The simulation of Hagen–Poiseuille flow in a pipe, driven by a constant body force in the axial direction az , is carried out. For the present test case, an analytical solution is available for the axial velocity uz ( r ) = u0

 1−

r2

 (51)

R2

where u0 = az R2 /4ν is the maximum axial velocity at the symmetry axis. The radius of the pipe R and the maximum axial velocity u0 lead to the definition of the Reynolds number, Re = 2Ru0 /ν . In addition to the symmetry conditions, the periodic conditions in the z-direction are applied at the inlet and outlet of the pipe and the boundary constraints at other boundaries are described as r = R : ur = uz = 0 .

(52)

A Nz × Nr = 33 × 17 lattice, where Nz and Nr are the number of lattice nodes respectively in the z and r directions, is adopted for the simulation of the present test case of ν = 0.2 and Re = 10. Fig. 1(a) presents the axial velocity (uz /u0 ) with respect to radial position (r /R, where r is the radial position assessed and R is the radius of the pipe; r /R = 0 and 1 refer respectively to the pipe radial center and wall). As demonstrated in Fig. 1(a), the velocity profile obtained with the present model agrees exactly with the analytical solutions and the results from Guo et al. [35]. Furthermore, the convergence rate of the present model with respect to mesh resolution is also investigated using the following global error definition

   present  − uanalytical  uz  z 2   . E (u) =  analytical   uz 

(53)

2

The global errors of the present model and that of Guo et al. [35] for different mesh resolutions are plotted in Fig. 1(b). It is demonstrated that the global errors pertaining to the proposed model converge at the same rate as that of Guo et al. [35], and more importantly that the present model confers slightly less errors. Therefore, the NEES-based boundary treatment for the symmetry condition is equivalent to the symmetry scheme in terms of accuracy as well as convergence rate, which confirms the reliability of the present symmetry boundary treatment. 3.2. Cylindrical cavity flow The cylindrical cavity flow, driven by the top lid which rotates around the cylindrical axis with a constant angular velocity Ω , has been widely adopted to evaluate different numerical schemes [35,33,34,53–55]. The flow structure is well-acknowledged to depend on two key parameters, namely, the aspect ratio RA = H /R (where H and R denotes respectively the height and radius of the cylindrical cavity) and Reynolds number Re = Ω R2 /ν . The boundary conditions are ur = uz = uθ = 0 at the bottom z = 0 and the outer wall r = R, and ur = uz = 0 and uθ = Ω r for the top lid. For the purpose of validating the present model, four test cases are investigated in this work with two mesh resolutions (namely, R = 50 and R = 100): (A) RA = 1.5, Re = 990, (B) RA = 1.5, Re = 1290, (C) RA = 2.5, Re = 1010, and (D) RA = 2.5, Re = 2200. The flow structure of the cylindrical cavity flow is characterized by the recirculation region (i.e., the breakdown bubble) along the symmetry axis. Accordingly, the streamlines of the swirling flow resulting from the present model for the four cases are depicted in Fig. 2. For the two cases with relatively lower Re, namely, cases (A) in Fig. 2(a) and (C) in Fig. 2(c), no breakdown bubble is found. With the increase of Re, one single breakdown bubble appears in case (B) and two bubbles develop in case (D), as demonstrated respectively in Fig. 2(b) and (d). The results presented in Fig. 2 agree well with those from previous numerical and experimental studies [35,33,34,53–56].

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

2761

Fig. 1. (a) Axial velocity (uz /u0 ) with respect to radial position (r /R), and (b) global error with respect to mesh resolution for Hagen–Poiseuille flow with ν = 0.2 and Re = 10.

Fig. 2. Streamlines of the cylindrical cavity flow pertaining to the R = 100 lattice: (a) case (A) whereby RA = 1.5 and Re = 990; (b) case (B) whereby RA = 1.5 and Re = 1290; (c) case (C) whereby RA = 2.5 and Re = 1010; and (d) case (D) whereby RA = 2.5 and Re = 2200.

Further quantitative comparisons are also carried out between the present model and the results obtained via the spectral element method (SEM), which has a high degree of accuracy, and also by Guo et al. [35]. In particular for SEM, it should be noted that the results in the literature deviate from one another and no consistent results are available [35,53–55]; however, the results from the present LB simulations has been proven to converge to the reference SEM solutions when a higher mesh resolution is applied. Fig. 3 exhibits the cylinder axial position (z /H, where z is the height assessed and H is the total height of the cylinder) with respect to the axial velocity on the symmetry axis (uz /u0 ) for the four cases obtained using the mesh resolution of R = 100. At this mesh resolution, excellent agreement is found for all four cases between the results from

2762

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

Fig. 3. Axial position (z /H) with respect to axial velocity on the symmetry axis (uz /u0 ) obtained with a mesh resolution of R = 100 (i.e., the radius R was divided into 100 lattices): (a) case (A) whereby RA = 1.5 and Re = 990; (b) case (B) whereby RA = 1.5 and Re = 1290; (c) case (C) whereby RA = 2.5 and Re = 1010; and (d) case (D) whereby RA = 2.5 and Re = 2200.

the present model, Guo et al. [35] and that obtained from SEM. Analogous to Fig. 3, Fig. 4 presents the results for a lower mesh resolution of R = 50. The agreement among all three results are reasonable for Fig. 4(a)–(c), but Guo et al. [35] gave a significant deviation in Fig. 4(d) at the higher axial positions. In particular for Fig. 4(d), the highest Re investigated of 2200 resulted in the axial velocity (uz /u0 ) from the Guo et al. model [35] becomes a lot lower than the SEM results, which indicates that the numerical evolution is prone to instability at this condition. Notably, results from the present model are consistent with the SEM results under all four conditions, as demonstrated in Fig. 4. Two more quantitative checks on the accuracy of the present model are carried out. The first is regarding one of the underlying assumptions that the fluid density ρ0 is a constant set to be 1. In the following, the mass conservation of the present model is investigated based on the current test cases. First, an average fluid density is defined as follows

 ρ¯ =

i ,j

|M0 (i, j)| (54)

Nr Nz

ˆα is the fluid density at the lattice nodes, and the mass difference, determined as 1ρ = ρ¯ − 1, is where M0 = 1r αg tabulated in Table 1. We find that the mass differences for all four test cases and both mesh resolutions are much less than 1.0e−8 and far smaller than the global errors in Fig. 1(b), which indicates that the assumption of ρ0 is sound in terms of having negligible effect on the evolution of the distribution functions. The second is regarding the accuracy of maintaining the incompressibility condition, which is investigated with the following definition of compressibility error 

  ∂(rur )  ∂r + E compress =

i ,j

Nr Nz



∂(ruz )  ∂z 

(55)

which is supposed to be 0 in accordance with the mass conservation for incompressible axisymmetric flows. Table 2 provides a quantitative comparison between the present model and the Guo et al. model [35] in terms of the compressibility errors, which reflect the accuracy in maintaining the mass conservation of each model. It is found that the compressibility errors

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

2763

Fig. 4. Axial position (z /H) with respect to axial velocity on the symmetry axis (uz /u0 ) obtained with a mesh resolution of R = 50 (i.e., the radius R was divided into 50 lattices): (a) case (A) whereby RA = 1.5 and Re = 990; (b) case (B) whereby RA = 1.5 and Re = 1290; (c) case (C) whereby RA = 2.5 and Re = 1010; and (d) case (D) whereby RA = 2.5 and Re = 2200. Table 1 Mass differences caused by the present model in the cylindrical cavity flow test. Grid

R = 50 R = 100

Re Re = 990

Re = 1010

Re = 1290

Re = 2200

6.71E−10 2.59E−10

7.70E−09 2.28E−09

9.66E−11 1.70E−09

3.06E−09 3.59E−09

Table 2 Compressibility errors caused by the present LB simulations of the cylindrical cavity flow. Grid

R = 50 R = 100

Present model

Guo et al. model [35]

Re = 990

Re = 1010

Re = 1290

Re = 2200

Re = 990

Re = 1010

Re = 1290

Re = 2200

2.77E−04 1.06E−04

1.61E−04 6.19E−05

2.98E−04 1.18E−04

2.02E−04 8.21E−05

3.33E−04 1.29E−04

1.92E−04 7.49E−05

3.60E−04 1.44E−04

2.40E−04 1.01E−04

from the Guo et al. model [35] are nearly 20% larger than that of the present model for each case. Because the present model is also based on the kinetic theory, the relationship between density and pressure is not completely eliminated. Thereby, the compressibility errors of the present LB simulations are improved but not significantly vis-à-vis Guo et al. [35], as shown in Table 2. Collectively, Tables 1 and 2 prove that the constant fluid density assumption is reasonable and the compressibility errors associated with the Guo et al. model [35] can be effectively reduced by the present model. 3.3. Natural convection in an annulus between two coaxial vertical cylinders Two thermal test cases are employed in the validation of the present model, whereby the D2Q4 base thermal LB equation in Eq. (46) is adopted for the temperature fields. Fig. 5 shows the schematic of the setup, whereby Ri , Ro and H represent

2764

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

Fig. 5. Schematic of the coaxial vertical cylinders.

Table 3 Nusselt numbers for the natural convection between coaxial vertical cylinders at different test cases. Ra

Present model

Guo et al. model [35]

Ref. [19]

Ref. [57]

103 104 105

1.697 3.228 5.810

1.697 3.228 5.804

1.692 3.215 5.787

– 3.163 5.882

respectively the inner cylinder radius, outer cylinder radius and height of the cylinders, and the ratio of the dimensions are such that Ro /Ri = H / (Ro − Ri ) = 2. Firstly, the simulation of natural convection in an annulus without swirling is performed. The no-slip boundary conditions are applied for the velocities on the four boundaries, namely, the inner and outer cylindrical surfaces, and the top and bottom walls. Constant temperatures of Ti and To (whereby Ti > To ) are imposed at the inner and outer cylindrical surfaces, respectively, while insulation condition is imposed on the bottom and top walls, i.e., ∂ T /∂ z = 0. The flow is driven by the temperature field through invoking the Boussinesq approximation, and the relevant body force ρ az = −ρ g β (T − Tm ) is added to the momentum equation, where Tm = (Ti + To ) /2 is the reference temperature, g is the gravitational acceleration and β is the thermal-expansion coefficient. In addition, the flow patterns depend on the dimensionless Prandtl number Pr = ν/α and the Rayleigh number Ra = g β (Ti − To ) (Ro − Ri )3 /αν . In this study, the mesh resolution corresponding to Ro − Ri = 100 is applied for the cases with Pr = 0.7 and Ra = 103 , 104 and 105 . Figs. 6 and 7 are respectively the streamlines and isotherms for the three cases. As Ra increases, the magnitude of the stream function becomes larger, thus the driving force for convection becomes larger and the steep changes in the flow variables gets more confined to the vicinity of the cylindrical surfaces. As demonstrated in Fig. 6, with the increase of Ra, the vortex pattern in the central region for Ra = 105 in Fig. 6(c) is much more complex than the other two cases (Fig. 6(a)–(b)). The corresponding temperature distributions in Fig. 7 demonstrate similar trends as the flow structure: with the increase of Ra, the temperature gradients near the cylindrical walls become larger. It should be noted that both streamlines and the isotherms from the present model are consistent with the reference results in the literature [35,19,57]. Subsequently, quantitative validations of the numerical results are performed by comparing the average Nusselt numbers between the present model and others in the literature [35,19,57]. Specifically, the average Nusselt number is defined as Nu = (Nui + Nuo ) /2, whereby Nui and Nuo represent the surface-averaged Nusselt numbers respectively on the inner and outer cylindrical walls:

Nui = −

Ri H (Ti − To )

Nuo = −

Ro H (Ti − To )

H

 0

H

 0

 ∂ T  dz ∂ r i  ∂ T  dz . ∂ r o

(56a)

(56b)

Table 3 lists the Nu values of all four models for the three test cases. Excellent agreement between the results from the present model and other published results are obtained, which further reflects the reliability of the present model.

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

2765

Fig. 6. Streamlines of the natural convection between the coaxial vertical cylinders at Pr = 0.7 and (a) Ra = 103 , (b) Ra = 104 and (c) Ra = 105 .

3.4. Swirling thermal flows in a cylindrical container Lastly, the present model is verified by the swirling flows in a cylindrical container with counter-rotating end disks (Fig. 8) under a constant temperature difference. The aspect ratio between the height and the radius of the container H /R equals 2. As for the boundary conditions, constant temperatures of Th and Tc (whereby Th > Tc ) are designated on the top and bottom wall, respectively, and the cylindrical surface is insulated. Besides, a no-slip boundary condition is imposed on the solid walls, and the top and bottom disks are rotating with a constant angular speed of Ωt and Ωb , respectively. The swirling flows are characterized by the following dimensionless parameters: Pr = ν/α , Re = Ωb R2 /ν , s = Ωt /Ωb and the Richardson number Ri = g β (Th − Tc ) /RΩb2 . The first three parameters are fixed as Pr = 1.0, Re = 1000 and s = −1. The Boussinesq assumption is valid in the present test for Ri = 1, and the force expression as well as the reference temperature definition is the same as the former test. The converged results for cases with Ri = 0 and 1 are obtained from the present model with respective mesh resolutions of R = 100 and R = 150. Fig. 9 provides the flow structures in terms of the contours of the azimuthal velocity and the temperature distribution for case of Ri = 0. Analogous to Fig. 9, Fig. 10 presents the same plots for the case of Ri = 1. As mentioned above, when Ri switches from 0 to 1, the Boussinesq assumption is invoked, and the resulting body forces, ρ az > 0 for z > H /2 and ρ az < 0 for z < H /2, causes the azimuthal velocity and the temperature distribution to be more confined near the top and bottom walls. The same phenomenon is also observed in earlier reports [38,58]. Furthermore,

2766

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

Fig. 7. Isotherms for the natural convection between coaxial vertical cylinders at Pr = 0.7 and (a) Ra = 103 , (b) Ra = 104 and (c) Ra = 105 .

profiles of each velocity component, namely, ur , uz and uθ , and the temperature distribution along the vertical line of r /R = 0.8 are given in Fig. 11 for the Ri = 0 case. Analogous to Fig. 11, Fig. 12 presents the same plots for the Ri = 1 case. Figs. 11 and 12 show that all the velocity and temperature profiles predicted by the present model agree excellently with Guo et al. [35] and Omi and Iwatsu [58], which again attests to the accuracy of the present model. Additionally, the quantitative profiles of the velocity components and the dimensionless temperature T ∗ = (T − Tm ) / (Th − Tc ) (whereby Tm = (Th + Tc ) /2) in Figs. 11 and 12 reflect the same trends as in Figs. 9 and 10, in that the macroscopic variables distribute linearly in the vicinity of the mid-height plane (z /H = 0.5) when Ri increases from 0 to 1. Additionally, the compressibility errors arising from the present LB simulations are compared in Table 4, based on the average compressibility error definition in Eq. (55). Again, the advantage of the proposed incompressible model is well demonstrated: the compressibility errors from the Guo et al. model [35] for various mesh resolutions are nearly 20% and 10% larger than those from the proposed model for the test cases of Ri = 0 and Ri = 1, respectively. We note that the introduction of the buoyancy force based on the Boussinesq assumption generally √ increases the compressibility errors of the involved LB simulations. However, as a smaller Mach number (Ma = u0 3/c) is applied in the case of Ri = 1 and

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

2767

Fig. 8. Schematic of the cylindrical container with counter-rotating end disks.

Fig. 9. (a) Azimuthal velocity contours and (b) isotherms for the case of Ri = 0.

R = 150 by choosing u0 = Ωb R = 0.05 (note that u0 = 0.1 in other cases), the compressibility error is therefore smaller than that of the Ri = 0 and R = 150 case.

2768

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

Fig. 10. (a) Azimuthal velocity contours and (b) isotherms for the case of Ri = 1.

Fig. 11. Axial profiles of velocity components, namely, ur , uz and uθ , and dimensionless temperature at r /R = 0.8 for case of Ri = 0.

4. Conclusion In this study, an alternative LB model for incompressible axisymmetric flows is developed, which advances the kinetic theory-based model developed by Guo et al. [35]. Summarily, the moment equations are first revised by introducing

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

2769

Fig. 12. Axial profiles of velocity components, namely, ur , uz and uθ , and dimensionless temperature at r /R = 0.8 for case of Ri = 1.

Table 4 Compressibility errors caused by the present LB simulations of the swirling thermal flows. R = 50

Guo et al. model [35] Present model

R = 100

R = 150

Ri = 0

Ri = 1

Ri = 0

Ri = 1

Ri = 0

Ri = 1

4.77E−04 4.02E−04

5.07E−04 4.63E−04

1.86E−04 1.52E−04

2.12E−04 1.94E−04

9.74E−05 8.06E−05

5.20E−05 4.48E−05

the incompressibility conditions, following which the moment method based on the Hermite tensorial polynomials expansion [40,41] is employed to derive the expressions of the equilibrium reduced distribution functions and the source terms. The present model is then completed by a novel formula for computing the fluid pressure. It should be noted that the present model is derived from the continuous axisymmetric Boltzmann equation, and therefore the advantages of the Guo et al. [35] model over other existing axisymmetric LB models, including the solid physics basis and the simple source terms containing no velocity gradients, are retained in the present model. Notably, since the fluid density in the present model is assumed to be a constant and the fluid pressure evolves independently of density, the macroscopic equations derived from the present model are the standard incompressible axisymmetric N–S equations. Furthermore, an additional relaxation parameter pertaining to the ghost mode is introduced for further improving the numerical stability of the present model in accordance with the idea of the RLBGK scheme [42,43]. Numerical validations of the proposed model are carried out with typical test cases. Comparisons with the reference solutions demonstrate three key highlights. First and foremost, accurate results are obtained from the present model for all the test cases investigated. Second, the convergence rates of the global errors with respect to the analytical solutions reflect the second-order accuracy of the present model. Third, for the cylindrical cavity flow, (i) the present model is slightly more accurate near the top lid when a lower mesh resolution is applied, (ii) the mass difference of the present model is far less than the overall errors, which affirms the constant density assumption, and (iii) the compressibility errors of Guo et al. [35] model are about 20% larger than the present model, which indicates that the present model is more in line with the incompressibility conditions.

2770

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

Acknowledgments We would like to thank for the Financial support from the National Research Foundation (NRF), Prime Minister’s Office, Singapore under its Campus for Research Excellence and Technological Enterprise (CREATE) program. This work is also financially supported by the National Natural Science Foundation of China (No. 11572062 and No. 11402215) and the Program for Changjiang Scholars and Innovative Research Team in University (No. IRT13043). Appendix In this section, the Chapman–Enskog analysis is conducted to investigate if the macroscopic equations pertaining to the modified moment equations in Eq. (19) of the present model agree well with the standard incompressible axisymmetric N–S equations. Substituting the following multiscale expansions ∞  ∂ ∂ = τfn , ∂t ∂ tn n=0

φ˜ =

∞ 

τfn φ˜ (n) ,

n =0

∞ 

f =

τfn f (n) ,

f˜ =

n =0

φ¯ =

∞ 

∞ 

τfn f˜ (n) ,

f¯ =

n =0

∞ 

τfn f¯ (n)

(A.1a)

n =0

τfn φ¯ (n)

(A.1b)

n =0

into the axisymmetric Boltzmann equation (Eq. (3)) and the derived evolution equations for the reduced distribution function (i.e. Eqs. (5a) and (5b)), the following set of equations in the successive order of the adopted expansion parameter τf is obtained as

τf0 : f (0) = f (eq) , τf1 :

f˜ (0) = f˜ (eq) ,

f¯ (0) = f¯ (eq)

(A.2)

∂ ξ 2 ∂ f (0) ξr ξθ ∂ f (0) ∂ (0) f + ξj f (0) + θ − = −f (1) ∂ t0 ∂ xj r ∂ξr r ∂ξθ

(A.3a)

∂ ˜ (0) ∂ ξr 1 ∂ φ˜ (0) = −f˜ (1) f + ξj f˜ (0) + f˜ (0) + ∂ t0 ∂ xj r r ∂ξr

(A.3b)

∂ ¯ (0) ∂ 2ξr (0) 1 ∂ φ¯ (0) f + ξj f¯ (0) + f¯ + = −f¯ (1) ∂ t0 ∂ xj r r ∂ξr

(A.3c)

τf2 :

∂ ˜ (1) ∂ ξr 1 ∂ φ˜ (1) ∂ ˜ (0) f + f + ξj f˜ (1) + f˜ (1) + = −f˜ (2) ∂ t1 ∂ t0 ∂ xj r r ∂ξr

∂ ¯ (0) ∂ ¯ (1) ∂ 2ξr (1) 1 ∂ φ¯ (1) = −f¯ (2) . f + f + ξj f¯ (1) + f¯ + ∂ t1 ∂ t0 ∂ xj r r ∂ξr

(A.4a)

(A.4b)

With the definition of the moment equations in Eq. (26b) and the reduced distribution functions expressed in Eqs. (27) and (28), the moments equations for the present model are determined as

ρ0 =



 f (0) dξ,

ρ0 u =



 f (0) ξ dξ,

ρ0 ui uj + pδij =



 f (0) ξi ξj dξ

   ρ0 RT ui δjk + uj δik + uk δij =  f (0) ξi ξj ξk dξ    ρ0 uθ = f¯ (0) dξ, ρ0 uθ u = f¯ (0) ξ dξ, ρ0 uθ RT δij = f¯ (0) ξi ξj dξ    2  2 (0) ˜ ρ0 uθ + p = φ dξ, ρ0 u uθ + RT = φ˜ (0) ξ dξ    ρ0 uθ u2θ + 3RT = φ¯ (0) dξ.

(A.5a) (A.5b) (A.5c) (A.5d) (A.5e)

Taking the moments of Eq. (A.3b), we have

 ρ0 ur ∂  ∂ρ0 + ρ0 uj + =0 ∂ t0 ∂ xj r

(A.6a)

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

 ρ 0 ui ur ∂ (ρ0 ui ) ∂  δir + ρ0 ui uj + pδij + − ρ0 u2θ = 0 ∂ t0 ∂ xj r r    ∂  ∂  ρ0 ui uj + pδij + ρ0 RT ui δjk + uj δik + uk δij ∂ t0 ∂ xk  1    ρ0 RT  (1) + ui δjr + uj δir + ur δij − ρ0 u2θ + RT ui δjr + uj δir = −Pij r

(1)

where Pij

(1)

Pij

=

r



2771

(A.6b)

(A.6c)

f˜ (1) ξi ξj dξ , and further derivations lead to

 = −ρ0 RT

∂ uj ∂ ui + ∂ xj ∂ xi



.

(A.6d)

Similarly, moments of Eq. (A.3c) are

 2ρ0 uθ ur ∂ (ρ0 uθ ) ∂  + ρ0 uθ uj + =0 ∂ t0 ∂ xj r  2ρ0 RTuθ  ∂ ∂  ρ 0 uθ  2 (1) δir − uθ + 3RT δir = −Qi ρ0 RTuθ δij + (ρ0 uθ ui ) + ∂ t0 ∂ xj r r  (1) where Qi = f¯ (1) ξi dξ is derived from Eq. (A.7b) as (1)

Qi

=

ρ0 RT r

uθ δir − ρ0 RT

∂ uθ . ∂ xi

(A.7a) (A.7b)

(A.7c)

The leading-order moments of Eqs. (A.4a) and (A.4b) are then provided as

∂ρ0 =0 ∂ t1  ∂ (ρ0 ui ) ∂ (1) 1 (1) δir (1) + Pij + Pir − φ dξ = 0 ∂ t1 ∂ xj r r ∂ (1) 2 (1) ∂ (ρ0 uθ ) + Q + Qr = 0 ∂ t1 ∂ xj j r 

(A.8a) (A.8b) (A.8c)

φ˜ (1) dξ , can be derived by taking the second-order moments of Eq. (A.3a)     ∂ ∂ 3 − φ˜ (1) dξ = φ˜ (0) dξ + φ˜ (0) ξj dξ + φ˜ (0) ξr dξ ∂ t0 ∂ xj r

where the unknown term in Eq. (A.8b),

=

2 r

ρ0 RTur .

(A.9) (1)

(1)

Substituting the definitions of Pij , Qi

and



φ˜ (1) dξ , Eqs. (A.8b) and (A.8c) are rewritten as

∂ 2 ui ρ0 RT ∂ ui ur ∂ (ρ0 ui ) − ρ0 RT − + ρ0 RT 2 δir = 0 ∂ t1 ∂ xj ∂ xj r ∂r r   2 ∂ (ρ0 uθ ) ∂ uθ ρ0 RT uθ ∂ uθ − ρ0 RT + − = 0. ∂ t1 ∂ xj ∂ xj r r ∂r

(A.10a)

(A.10b)

Combining Eqs. (A.10), (A.8a), (A.6) and (A.7), the macroscopic equations pertaining to the present model are obtained as in Eq. (20) with the dynamic viscosity µ = ρ0 τf RT . It is hence demonstrated that the macroscopic equations recovered from the present model agree well with the standard incompressible axisymmetric N–S equations. References [1] [2] [3] [4] [5] [6]

R. Benzi, S. Succi, M. Vergassola, The Lattice Boltzmann-equation - Theory and applications, Phys. Rep. 222 (1992) 145–197. F.J. Higuera, J. Jimenez, Boltzmann approach to lattice gas simulations, Europhys. Lett. 9 (1989) 663–668. F.J. Higuera, S. Succi, Simulating the flow around a circular-cylinder with a Lattice Boltzmann-equation, Europhys. Lett. 8 (1989) 517–521. F.J. Higuera, S. Succi, R. Benzi, Lattice gas-dynamics with enhanced collisions, Europhys. Lett. 9 (1989) 345–349. G.R. Mcnamara, G. Zanetti, Use of the Boltzmann-equation to simulate lattice-gas automata, Phys. Rev. Lett. 61 (1988) 2332–2335. X.Y. He, L.S. Luo, Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation, Phys. Rev. E 56 (1997) 6811–6817. [7] M. Junk, A finite difference interpretation of the lattice Boltzmann method, Numer. Methods Partial Differential Equations 17 (2001) 383–402.

2772

L. Zhang et al. / Computers and Mathematics with Applications 72 (2016) 2751–2772

[8] M. Junk, A. Klar, Discretizations for the incompressible Navier–Stokes equations based on the lattice Boltzmann method, SIAM J. Sci. Comput. 22 (2000) 1–19. [9] X. Shan, H. Chen, Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation, Phys. Rev. E 49 (1994) 2941–2948. [10] X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815–1819. [11] S. Tao, Z.L. Guo, Boundary condition for lattice Boltzmann modeling of microscale gas flows with curved walls in the slip regime, Phys. Rev. E 91 (2015) 12. [12] X.J. Yue, Z.H. Wu, Y.S. Ba, Y.J. Lu, Z.P. Zhu, D.C. Bo, Lattice-Boltzmann simulation for pressure driven microscale gas flows in transition regime, Internat. J. Modern Phys. C 26 (2015) 9. [13] J. Zhang, Lattice Boltzmann method for microfluidics: models and applications, Microfluid. Nanofluid. 10 (2011) 1–28. [14] C. Lin, A. Xu, G. Zhang, Y. Li, S. Succi, Polar-coordinate lattice Boltzmann modeling of compressible flows, Phys. Rev. E (3) 89 (2014) 013307. [15] G. Yanbiao, X. Aiguo, Z. Guangcai, Y. Yang, Lattice BGK kinetic model for high-speed compressible flows: Hydrodynamic and nonequilibrium behaviors, Europhys. Lett. 103 (2013) 24003. [16] A.M. Artoli, A.G. Hoekstra, P.M.A. Sloot, 3D pulsatile flow with the lattice Boltzmann BGK method, Internat. J. Modern Phys. C 13 (2002) 1119–1134. [17] R.S. Maier, R.S. Bernard, D.W. Grunau, Boundary conditions for the lattice Boltzmann method, Phys. Fluids 8 (1996) 1788–1801. [18] R. Mei, W. Shyy, D. Yu, L.-S. Luo, Lattice Boltzmann method for 3-D flows with curved boundary, J. Comput. Phys. 161 (2000) 680–699. [19] L. Li, R. Mei, J.F. Klausner, Multiple-relaxation-time lattice Boltzmann model for the axisymmetric convection diffusion equation, Int. J. Heat Mass Transfer 67 (2013) 338–351. [20] I. Halliday, L.A. Hammond, C.M. Care, K. Good, A. Stevens, Lattice Boltzmann equation hydrodynamics, Phys. Rev. E 64 (2001) 011208. [21] T.S. Lee, H. Huang, C. Shu, An axisymmetric incompressible lattice Boltzmann model for pipe flow, Internat. J. Modern Phys. C 17 (2006) 645–661. [22] T. Reis, T.N. Phillips, Numerical validation of a consistent axisymmetric lattice Boltzmann model, Phys. Rev. E 77 (2008) 026703. [23] T. Reis, T.N. Phillips, Modified lattice Boltzmann model for axisymmetric flows, Phys. Rev. E 75 (2007) 056703. [24] T.S. Lee, H. Huang, C. Shu, An axisymmetric incompressible lattice BGK model for simulation of the pulsatile flow in a circular pipe, Internat. J. Numer. Methods Fluids 49 (2005) 99–116. [25] S. Mukherjee, J. Abraham, Lattice Boltzmann simulations of two-phase flow with high density ratio in axially symmetric geometry, Phys. Rev. E 75 (2007) 026701. [26] K.N. Premnath, J. Abraham, Lattice Boltzmann model for axisymmetric multiphase flows, Phys. Rev. E 71 (2005) 056706. [27] H. Huang, M. Sukop, X. Lu, Multiphase Lattice Boltzmann Methods: Theory and Application, John Wiley & Sons, 2015. [28] Y. Peng, C. Shu, Y.T. Chew, J. Qiu, Numerical investigation of flows in Czochralski crystal growth by an axisymmetric lattice Boltzmann method, J. Comput. Phys. 186 (2003) 295–307. [29] J.G. Zhou, Axisymmetric lattice Boltzmann method, Phys. Rev. E 78 (2008) 036701. [30] S. Chen, J. Tölke, S. Geller, M. Krafczyk, Lattice Boltzmann model for incompressible axisymmetric flows, Phys. Rev. E 78 (2008) 046703. [31] S. Chen, J. Tölke, M. Krafczyk, Simulation of buoyancy-driven flows in a vertical cylinder using a simple lattice Boltzmann model, Phys. Rev. E 79 (2009) 016704. [32] Q. Li, Y.L. He, G.H. Tang, W.Q. Tao, Lattice Boltzmann model for axisymmetric thermal flows, Phys. Rev. E 80 (2009) 037702. [33] Q. Li, Y.L. He, G.H. Tang, W.Q. Tao, Improved axisymmetric lattice Boltzmann scheme, Phys. Rev. E 81 (2010) 056707. [34] J.G. Zhou, Axisymmetric lattice Boltzmann method revised, Phys. Rev. E 84 (2011) 036704. [35] Z. Guo, H. Han, B. Shi, C. Zheng, Theory of the lattice Boltzmann equation: Lattice Boltzmann model for axisymmetric flows, Phys. Rev. E 79 (2009) 046708. [36] L. Wang, Z. Guo, C. Zheng, Multi-relaxation-time lattice Boltzmann model for axisymmetric flows, Comput. & Fluids 39 (2010) 1542–1548. [37] L. Zheng, B. Shi, Z. Guo, C. Zheng, Lattice Boltzmann equation for axisymmetric thermal flows, Comput. & Fluids 39 (2010) 945–952. [38] L. Zheng, Z. Guo, B. Shi, C. Zheng, Kinetic theory based lattice Boltzmann equation with viscous dissipation and pressure work for axisymmetric thermal flows, J. Comput. Phys. 229 (2010) 5843–5856. [39] H. Liang, Z.H. Chai, B.C. Shi, Z.L. Guo, T. Zhang, Phase-field-based lattice Boltzmann model for axisymmetric multiphase flows, Phys. Rev. E 90 (2014) 063311. [40] H. Grad, Note on N-dimensional hermite polynomials, Comm. Pure Appl. Math. 2 (1949) 325–330. [41] H. Grad, On the kinetic theory of rarefied gases, Comm. Pure Appl. Math. 2 (1949) 331–407. [42] A. Montessori, G. Falcucci, P. Prestininzi, M. La Rocca, S. Succi, Regularized lattice Bhatnagar-Gross-Krook model for two- and three-dimensional cavity flow simulations, Phys. Rev. E 89 (2014) 053317. [43] A. Montessori, M. La Rocca, G. Falcucci, S. Succi, Regularized lattice BGK versus highly accurate spectral methods for cavity flow simulations, Internat. J. Modern Phys. C 25 (2014). [44] P.J. Dellar, Nonhydrodynamic modes and a priori construction of shallow water lattice Boltzmann equations, Phys. Rev. E (3) 65 (2002) 036309. [45] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511. [46] Y.H. Qian, D. Dhumieres, P. Lallemand, Lattice Bgk models for Navier–Stokes equation, Europhys. Lett. 17 (1992) 479–484. [47] R. Benzi, S. Succi, M. Vergassola, Turbulence modelling by nonhydrodynamic variables, Europhys. Lett. 13 (1990) 727. [48] P.J. Dellar, Non-hydrodynamic modes and general equations of state in lattice Boltzmann equations, Physica A 362 (2006) 132–138. [49] Z.L. Guo, C.G. Zheng, B.C. Shi, An extrapolation method for boundary conditions in lattice Boltzmann method, Phys. Fluids 14 (2002) 2007–2010. [50] G. Zhao-Li, Z. Chu-Guang, S. Bao-Chang, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chin. Phys. 11 (2002) 366. [51] H.B. Huang, X.Y. Lu, M.C. Sukop, Numerical study of lattice Boltzmann methods for a convection–diffusion equation coupled with Navier–Stokes equations, J. Phys. A 44 (2011) 055001. [52] W. Xie, An axisymmetric multiple-relaxation-time lattice Boltzmann scheme, J. Comput. Phys. 281 (2015) 55–66. [53] S.K. Bhaumik, K.N. Lakshmisha, Lattice Boltzmann simulation of lid-driven swirling flow in confined cylindrical cavity, Comput. & Fluids 36 (2007) 1163–1173. [54] Y. Wang, C. Shu, C.J. Teo, A fractional step axisymmetric lattice Boltzmann flux solver for incompressible swirling and rotating flows, Comput. & Fluids 96 (2014) 204–214. [55] T. Zhang, B. Shi, Z. Chai, F. Rong, Lattice BGK model for incompressible axisymmetric flows, Commun. Comput. Phys. 11 (2012) 1569–1590. [56] K. Fujimura, H.S. Koyama, J.M. Hyun, Time-Dependent vortex breakdown in a cylinder with a rotating lid, J. Fluids Eng. 119 (1997) 450–453. [57] M. Venkatachalappa, M. Sankar, A.A. Natarajan, Natural convection in an annulus between two rotating vertical cylinders, Acta Mech. 147 (2001) 173–196. [58] Y. Omi, R. Iwatsu, Numerical study of swirling flows in a cylindrical container with co-/counter-rotating end disks under stable temperature difference, Int. J. Heat Mass Transfer 48 (2005) 4854–4866.