New thermal MRT lattice Boltzmann method for simulations of convective flows

International Journal of Thermal Sciences 100 (2016) 98e107

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

New thermal MRT lattice Boltzmann method for simulations of convective flows Mohammed Jami a, Fayçal Moufekkir a, Ahmed Mezrhab a, *, Jean Pierre Fontaine b, M'hamed Bouzidi c Laboratoire de M ecanique & Energ etique, D epartement de Physique, Facult e des Sciences, Universit e Mohammed 1er, 60000 Oujda, Morocco Clermont Universit e, Universit e Blaise Pascal, Institut Pascal e axe GePEB (UMR-CNRS 6602), BP 20206, F-63174 Aubi ere Cedex, France c Clermont Universit e, Universit e Blaise Pascal, Institut Pascal e axe MMS (UMR-CNRS 6602), IUT de Montluçon, av. A. Briand, CS 82235, F-03101 Montluçon Cedex, France a

b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 25 February 2015 Received in revised form 1 July 2015 Accepted 16 September 2015 Available online xxx

The purpose of this paper is to develop a numerical scheme based on the MRT-D2Q9 model for fluid velocity and the MRT-D2Q5 cross model for the temperature to study the natural convection phenomenon in a differentially heated square enclosure containing air. The present work aims to prove the accuracy and the capacity of this model to predict the thermodynamics phenomena. In fact, the main goal consists to demonstrate the efficiency of the double population approach based on new MRT-LBM scheme in the assessment of laminar and transitional flows in natural convection case for different Rayleigh numbers (up to 108). The numerical results of the present work were successfully checked for accuracy and are in satisfactory agreement with those available in literature. The results analysis shows that the present method has proven its ability to handle similar problem. © 2015 Published by Elsevier Masson SAS.

Keywords: New thermal MRT lattice Boltzmann method Convective flows

1. Introduction During last years, the Lattice Boltzmann Equation (LBE), which was derived from Lattice Gas Automata (LGA), has received considerable attention of the research community in reason to the success achieved in simulating hydro-thermodynamic problems. It is considered as potentially efficient for solving a large class of problems. In fact, it has been successfully used to simulate fluid flows and heat transfer. Because of its mesoscopic origin, the LBM appears as a flexible method with many advantages, compared to conventional numerical methods, which are based on the discretization of macroscopic governing equations. Actually, the advantages of LBM include simplicity, efficient implementation for parallel computing and easy handling for complex geometries or boundary conditions. Another merit of this method is its high performance in terms of stability and accuracy. Frisch et al. [1] was among the first authors who conducted twodimensional simulation of hydrodynamic problems using a Lattice

* Corresponding author. E-mail address: [email protected] (A. Mezrhab). http://dx.doi.org/10.1016/j.ijthermalsci.2015.09.011 1290-0729/© 2015 Published by Elsevier Masson SAS.

Boltzmann method. Then, McNamara and Zanetti [2] used the real variables instead of the Boolean quantities in reason to remove the back draws of LGA models related to the Galilean invariance, isotropy and background noise due to the use of discrete quantities. In order to meet the computation needs, many authors [3,4] interested to the linearization of the collision expression and the improvement of its term form, which was fit continuously. The most sophisticated release of the collision expression named BGK operator was adopted by Chen and Qian [5,6]. Various numerical simulations connected to thermal and fluid flows were performed to demonstrate the reliability and capability of the LBE in modeling complex phenomena. A review on the use of LBM in the simulation of more complex phenomena such as compressible flows, multiphase and porous media flows can be found in the references [7e16]. Several models of thermal lattice Boltzmann have been developed to check their ability to simulate more physical phenomena. The first one called thermal LBE model with passive scalar consists of the cleavage of the thermal and dynamical fields in the numerical computations [17,18]. The stability and efficiency in computations are increased under suitable modeling of the thermal field, which neglect the thermal viscous dissipation and compression. The

M. Jami et al. / International Journal of Thermal Sciences 100 (2016) 98e107

Nomenclature A g H L Nu Pr Ra RD T T Tm u, v U, V x, y X, Y

aspect ratio, H/L ¼ 1 gravitational acceleration, m s2 cavity height, m length of the cavity, m average Nusselt number Prandtl number, y/a Rayleigh number relative difference time, s temperature, K average temperature, (Th þ Tc)/2, K velocity components, m s1 horizontal and vertical velocity Cartesian coordinates, m dimensionless coordinates, (x/H, y/H)

multi-speed approach [19,20] is known in the literature as the second thermal LBE model and is characterized by adding the density distribution function to the internal energy in their terms. The hybrid method is considered as a sophisticated approach of the last model, since it allows improving the drawback in terms of numerical stability due to the use of limited thermal intervals in computations. The last model is recognized by it high compatibility when coupled with various CFD tools such as the finite difference method [21,22], which explain the improvement of computation efficiency. The next model, named the double population distribution function DDF approach, uses the entropic Boltzmann equation and the kinetic principles to model the Knudsen effects. This method is unconditionally stable and producing faster computations, when the H-theorem is taken into account. This approach affords to avoid the negative entropy generation and fixes the instabilities occurring related to the ELBE method [23,24]. Despite these advantages, there are some deficiencies related to the LBGK method used in most works based on a DDF approach, such as the inaccuracy of implementing boundary conditions and the numerical instability that can sometimes occur [25,26]. These short-comings can be overcome by using the multiple relaxation time lattice Boltzmann method (MRT-LBM) introduced by res [27]. Because of its advantages compared to the LBMd'Humie BGK commonly called the single relaxation time lattice Boltzmann method (SRT-LBM), the MRT-LBM was designed to produce more stable solution than its LBM-BGK counterpart [28e31]. The main idea is that the collision process is accomplished in the momentum space by a linear transformation and the advection is still performed in the velocity space. Recently, the LBM was considered as an alternative CFD tool to other numerical methods when dealing with different topics related to the addition of force term, the turbulent flows and unsteady flows of nanofluids. Indeed, the numerical assessment of adding force term to LBM with BGK method was applied to the natural convection case by comparing with results provided by usual CFD tools [32]. In fact, the incorporation of dynamic subgrid scale (SGS) models in the lattice-Boltzmann method (LBM) for large-eddy simulation (LES) of turbulent flows was studied using the MRT-LBM [33]. Furthermore, the LBM was proven its ability to simulate the laminar mixed convection of watereCu nanofluid in an inclined shallow driven cavity [34]. More details about the LBM applications previously discussed can be found in the literature such as the analytical investigation of the heat transfer for the

99

Greek symbols thermal diffusivity, m2 s1 volumetric expansion coefficient, K1 dynamic viscosity of fluid, kg m1 s1 kinematic viscosity, m2 s1 density of fluid, kg m3 relaxation time dimensionless temperature, (T  Tm)/(Th  Tc) stream function

a b m n r t q j

Subscripts c cold eq equilibrium h hot max maximum Superscripts þ post collision state

microchannel heat sink cooled by different nanofluids conducted by Ref. [35]. In the same trend, the last analytical model was investigated for unsteady flow of a nanofluid squeezing between two parallel plates by the authors [36]. Similar study was carried out for unsteady motion of a rigid spherical particle in a quiescent shear-thinning power-law fluid by the same authors [37]. Most methods described in the literature produces similar results, despite that some authors claim that their scheme is more accurate than the other schemes. The main purpose of the present work focuses on modeling natural convection in a square cavity using a novel DDF approach, which provides accurate solutions for 2-D convective flows. To carry out this, the adopted numerical approach uses the multiple relaxation time-lattice Boltzmann equation (MRT-LBE) with D2Q9 lattice model to obtain the velocity distribution, and the MRT-LBE with D2Q5 cross lattice model to obtain the temperature distribution. The novelty of this letter, compared to the existing works on LBM schemes, is the use of new double population approach based on cross MRT-LBM scheme (with view to evaluate its effectiveness) to compute thermal field in simulation of benchmark problems related to laminar and transitional convective flows shows. The paper is organized as follows: in the next section the MRTLBE associated to the D2Q9 lattice model and to the D2Q5 cross lattice model are presented, the validation test is reported in the third section, then the results of numerical simulations of the 2-D incompressible flow inside a differentially heated square cavity are presented in Section 2. Finally conclusions of numerical study are done in Section 3. 2. Application of the MRT-LBE to a convective flow in a square cavity 2.1. D2Q9-MRT LBE model Unlike the traditional approaches depending on the numerical solution of the Navier-Stokes equations (macroscopic continuum equation), the LBM divides time and space into steps to form a lattice and discretize the fluid as particles. These particles are placed on the nodes of this grid called lattice sites (or cells) and they move at successive times from nodes where they are located toward the neighboring node. In the LBM, a simulation region is regarded as a lattice system, and fluid particles are made to move from site to site in the lattice

100

M. Jami et al. / International Journal of Thermal Sciences 100 (2016) 98e107

system. Numerous lattice models have already been introduced in the modeling of thermal flow phenomena and hydrodynamics systems. In the present study, the D2Q9 model is employed to compute the dynamical field as shown in Fig. 1. In this model, the fluid particles have 8 possibilities for moving to the neighboring sites, and including the quiescent state, they have 9 velocities as shown in Fig. 1. If an arbitrary lattice point is denoted by r and the particle distribution function in the i direction at time t is denoted by fi(rj, tn), then the particle distributions fi(rj þ eidt, tn þ dt) at rj þ eidt after the time dt can be obtained from the following equation:



 

 f r j þ ei dt ; tn þ dt  f i rj ; tn ¼ S f rj ; tn  f eq r j ; tn

(1)

in which feq(rj, tn) is the discrete equilibrium distribution function and S is the collision matrix which can be expressed in the moment space [22,27] and gives the following equation:



 

 f r j þ ei dt ; tn þ dt  f rj ; tn ¼ M1 S* m rj ; tn  meq rj ; tn (2) where M is a transformation matrix projecting the discrete distribution function f onto the moment space m ¼ Mf, and S* is the diagonal matrix of relaxation rates {sa, 0  a  8}. The bold face symbols such as f stand for nine-component vectors, m and meq are nine-component vectors of moments and their equilibria. Nine being the number of discrete velocities, as follows: T

f ¼ ðf0 ; f1 ; …; f8 Þ ;

(4)

where m0 ¼ r is the fluid density, m1 ¼ e is related to the total energy, m2 ¼ 3 is related to the energy square, (m3, m5) ¼ (jx, jy) is the flow momentum, (m4, m6) ¼ (qx, qy) is related to the energy flux, and (m7, m8) ¼ (pxx, pyy) are related to the diagonal and off-diagonal components of the symmetric stress tensor, respectively. In addition, in Eq. (2), the equilibrium moments meq(r, j) are [25]:

  eq m1 ¼ eeq ¼ 2r þ 3 j2x þ j2y .  eq r m2 ¼ 3eq ¼ r  3 j2x þ j2y eq meq 4 ¼ qx ¼ jx eq eq m6 ¼ qy ¼ jy .  eq 2 2 rm m7 ¼ peq xx ¼ jx  jy . eq eq m8 ¼ pxy ¼ jx jy r

(5)

rm is the mean density and is taken equal to one in our study. The non-conserved moments relax linearly towards their equilibrium values as suggested in Ref. [25]. According to the references [22,25,27,31], the collision process of the MRT model for the D2Q9 lattice is achieved in step:

mþ ¼ m  S  ðm  meq Þ

(6)

where mþ is the moment vector in the post-collision state. The post-collision vector fþ is carried out as follows

T

m ¼ ðm0 ; m1 ; …; m8 Þ

where the superscript T denotes transpose operator. It should be noted that the nine-velocity square lattice Boltzmann (D2Q9) model has widely and successfully been used for simulating two-dimensional thermal flows. In the D2Q9 model, ea denotes the discrete velocity set, namely

8 < ð0; 0Þ; ðcos½ða  1Þp=2; sin½ða  1Þp=2Þc; ea ¼ p : ffiffiffi 2ðcos½ð2a  9Þp=4; sin½ð2a  9Þp=4Þc;

T  m ¼ r; e; 3; jx ; qx; jy ; qy ; pxx ; pxy

f þ ¼ f  M 1 S* Mðf  f eq Þ

Note that the transformation matrix M is an orthogonal one and it can be constructed via the GrameSchmidt orthogonalization process [22,25,27,31,38]. From these references, and for the D2Q9 model, the matrix M can be written as:

0

a¼0 1a4 5a8 (3)

where dx and dt are the spacing and the time step, respectively (these quantities are fixed at the unity (dx ¼ dt ¼ 1)), c ¼ dx/dt denotes the particle velocity. In the following, all quantities are normalized by the lattice spacing dx and time increment dt and defined in dimensionless formulation. The moments are arranged in the following order:

(7)

B B B B B B M¼B B B B B B @

1 4 4 0 0 0 0 0 0

1 1 2 1 2 0 0 1 0

1 1 2 0 0 1 2 1 0

1 1 2 1 2 0 0 1 0

1 1 2 0 0 1 2 1 0

1 2 1 1 1 1 1 0 1

1 2 1 1 1 1 1 0 1

1 2 1 1 1 1 1 0 1

1 2 1 1 1 1 1 0 1

1 C C C C C C C C C C C C A (8)

It is worth noting that this matrix is such that the matrix product MMT is a non-normalized diagonal matrix. This facilitates the calculation of the inverse matrix M1 using the formula [29]:

 1 M1 ¼ M T MM T

(9)

It should be noted that the relaxation S* is given by:


 S* ¼ diag s0 ; s1 ; s2 ; s3 ; s4 ; s5 ; s6 ; s7; s8

(10)

where s0 ¼ s3 ¼ s5 ¼ 0 enforces mass and momentum conservation before and after collision (see Refs. [22,25,27,31,39,40] for example), and the relaxation rates s7 ¼ s8 ¼ 1/t determine the dimensionless kinematic viscosity of model:

n ¼ ðt  0:5Þ=3 Fig. 1. Lattice structure (unit cell) for the D2Q9 model in 2D LBM.

(11)

It should be stressed that the relaxation parameters sa can be

M. Jami et al. / International Journal of Thermal Sciences 100 (2016) 98e107

101

determined by a linear stability analysis [25], and in the present simulation, they are chosen as the following: s1 ¼ s2 ¼ 1.4, s4 ¼ s6 ¼ 1.2 and s7 ¼ s8 ¼ 2/(1 þ 2n/3) being the kinematic viscosity of the fluid. What should be mentioned here is that it is possible to recover the SRT-LBE, also called LBGK after its authors BhatnagareGrosseKrook [39], by setting s1 ¼ s2 ¼ s4 ¼ s6 ¼ s7 ¼ s8 ¼ 1/t. In the other hand, the macroscopic fluid variables, density r and momentum ru on each lattice node are determined by the following formulas:

r¼

8 X

fi

(12)

i¼0 8   X fi ei þ G=2 J jx ; jy ¼ ru ¼

Fig. 2. Lattice structure for the D2Q5 model.

(13)

i¼1

where G is the G ¼ rbg(T  Tm).

external

force

acting

per

unit

mass:

2.2. D2Q5-MRT LBE model for thermal problem In contrast to our previous works [41e43], in which the energy equation was solved by finite difference approach, the temperature field in the present work is computed using the double-distribution function (DDF) model [31], in reason to demonstrate the capability of the MRT-D2Q5 cross model to deal with similar problems. A validation test comparing the D2Q5 model with finite-difference Modified LaxeWendroff (MLW) scheme was conducted by Rasin et al. [44]. The authors concluded that the kinetic LB scheme is significantly faster than MLW in terms of computational cost and they found that their numerical results are in satisfactory agreement with the analytical solution. The 2-D transient conduction and radiation heat transfer problem was addressed by Mishra et al. [45] using alternatively a finite volume method (FVM) and a lattice Boltzmann method to compute the temperature field. It is drawn that the LBM present many advantages for complex geometries and multi-dimensional flows, such as fast convergence to the steady state and accuracy. In the same trend, the present work aims to highlight the effectiveness of the double MRT-TLBE to simulate convective flows using a new powerful method instead of the classical CFD tools such as finite difference methods (FDM) or finite volume methods (FVM). The D2Q5 LBE cross model is introduced to solve the advectionediffusion equation for the temperature by using a second set of velocity distribution functions fga 2Vð¼ R5 Þ 0  a  4g defined on each node of the regular lattice parameterized by a space step dx, rj 2 (dxZ)2. c ¼ dx/dt denotes the celerity and dt is defined as the time step of the evolution of LBE. The velocities ci, i 2 (0, … , 4) are chosen such that ci ¼ eidx/dt, where ei are the vectors connecting neighboring nodes of the regular lattice. Fig. 2 illustrates the D2Q5 (five-velocity) cross model used in this work. The discrete velocities are defined as follows:

ea ¼

ð0; pffiffiffi0Þ; 2ðcos½ð2a  1Þp=4; sin½ð2a  1Þp=4Þc;

a¼0 1a4 (14)

A LBE with multiple relaxation times for distribution functions gj can be expressed as [46]:

gj ðri þ ci dt; t þ dtÞ ¼ gj* ðri ; tÞ

(15)

In LBM, two fundamental processes (the collision substep and

the advection substep) take place during each time step dt. Here the superscript (*) describes the post-collision state. In the advection substep, the “particles” move from a lattice node ri to either itself with the velocity c0 ¼ 0, or one of the four nearest neighbors with the velocity ci, i 2 (0, … , 4). As for the collision substep, it consists in the redistribution of the populations gj at each node ri, and it is modeled by the operator g*j (ri, t) (see Eq. (15)). The collision substep is not easy to carry out in the velocity space. It becomes convenient when this substep is carried out in the moment space as shown in Ref. [27]. To this end, one defines some moments mk based on the distribution gj through a linear transformation,

mk ¼

4 X

Mkj gj ;

0k4

(16)

j¼0

The transformation matrix M (invertible and orthogonal) for this model can be expressed by:

2

1 6 0 6 M¼6 6 0 4 4 0

1 1 1 1 1

1 1 1 1 1

1 1 1 1 1

3 1 1 7 7 1 7 7 1 5 1

(17)

Only the first moment is conserved, in the simulation of the diffusion problems,

m0 ¼ T ¼

4 X

gj

(18)

j¼0

It is assumed that non-conserved moments, in the collision substep, relax towards equilibrium moments meq that are nonlinear k functions of conserved quantities according to the following macroscopic scalar equation:

m*k ¼ ð1  sk Þmk þ sk meq ; k

0k4

(19)

where sk ¼ dt/tk is a relaxation rate which satisfy the constraint 0 < sk < 2 to get a numerically stable scheme. The relaxation rates sk are not necessarily identical as in the so-called BGK method [6,39]. eq eq eq eq Choosing m1 ¼ 0, m2 ¼ 0, m3 ¼ aT, m4 ¼ bT and using Taylor expansion [46,47], we find the equivalent anisotropic diffusion equation up to order three in dt:

102

M. Jami et al. / International Journal of Thermal Sciences 100 (2016) 98e107

 

fj ðxs ; t þ 1Þ ¼ fi xs þ ej ; t þ 1 and gj ðxs ; t þ 1Þ ¼ gi xs þ ej ; t þ 1

  2 vT c2 dt 1 1 v2 T 1 1 v T 2  ð4 þ aÞ   c dtb þ  1 vt 5 s1 2 vx2 s1 s2 vxvy    c2 dt 1 1 v2 T ð4 þ aÞ   ¼ 0 dt 3 2 5 s2 2 vy

(23)

(20) Eq. (20) reduces to the standard isotropic diffusion equation for b ¼ 0, s1 ¼ s2 ¼ s, dx ¼ dt ¼ 1, a ¼ 2 and whose diffusion coefficient is

a¼

 ð4 þ aÞ 1 1  5 s 2

(21)

With a given velocity field (u, v), if we take meq 1 ¼ uT and meq 2 ¼ vT, the LBE scheme describes the following advectionediffusion equation:

vT v2 T v2 T  uVT  a þ vt vx2 vy2

!

  ¼ 0 dt 2

(22)

2.3. Boundary conditions The usual bounce-back condition [31,40] is applied at all solid cavity walls that are located with a good approximation between the last ‘fluid’ node and the first ‘solid’ point. Indeed, this condition is applied at all walls for the fluid velocity and on the adiabatic horizontal walls for the temperature. The bounce-back boundary condition is considered as appropriate in terms of accuracy for simple boundary geometries [48,49]. As shown in Fig. 3, this con-

   1 1 1 u N þ ; y ¼ 0; v N þ ; y ¼ 0; q N þ ; y ¼ qc ¼ 0:5 2 2 2    1 1 1 u ; y ¼ 0; v ; y ¼ 0; q ; y ¼ qh ¼ 0:5 2 2 2    1 1 vq 1 u x; N þ ¼ 0; v x; N þ ¼ 0; x; N þ ¼0 2 2 vy 2    1 1 vq 1 ¼ 0; v x; ¼ 0; x; ¼0 u x; 2 2 vy 2

where fj and gj are the distribution functions of the velocity and the temperature, respectively and j is such that ej ≡ ei. For the isothermal condition applied at the left and right walls of the cavity, we adopt the following condition: The distribution function g can be written as a function of the inverse moment matrix M1 as:

gj ¼

4  X k¼0

M 1

 jk

mk

(24)

and, g1 þ g3 ¼ 2T(1/5 þ a/20), g2 þ g4 ¼ 2T(1/5 þ a/20) with a ¼ 2. At the top adiabatic wall (x, N þ 1/2): g3 ¼ g1, g4 ¼ g2. At the bottom adiabatic wall (x, 1/2): g1 ¼ g3, g2 ¼ g4. Fig. 4 illustrates this concept applied at the bottom wall. It is assumed that the node (x, 0) from which g2 and g1 come has distributions mirroring those of the node (x, 1) from which g4 and g3 come (see Fig. 4). Consequently, the temperature at the node (x, 0) is the same as that of the node (x, 1). This assumption implies that applying a second order central difference approximation to vq/vy at the wall boundary node will lead to the following condition at the macroscopic level, vq/vy ¼ 0. At the left wall boundary (1/2, y), g1 ¼ g3 þ 2/5qh(1 þ a/4), g4 ¼ g2 þ 2/5qh(1 þ a/4) At the right wall boundary (N þ 1/2, y), g3 ¼ g1 þ 2/5qc(1 þ a/ 4), g2 ¼ g4 þ 2/5qc(1 þ a/4) This is the mesoscopic form of the boundary conditions used. The macroscopic form is given by:

for

1 1 yNþ 2 2

for

1 1 yNþ 2 2

for

1 1 xNþ 2 2

for

1 1 xNþ 2 2

dition defines the position of the physical wall at the half grid spacing beyond the last fluid node xf. Note that the particle moves from xf toward xw and comes back to its location after being reflected by the stationary wall. This is expressed by the following equations:

We note that the computational domain is a square grid of size Nx*Ny ¼ N2 for Nx ¼ Ny ¼ N. According to the bounce-back boundary rule, the effective boundaries will be at 1/2 and N þ 1/2 as shown in Fig. 5.

Fig. 3. Bounce-back boundary condition.

Fig. 4. Illustration of the boundary wall.

M. Jami et al. / International Journal of Thermal Sciences 100 (2016) 98e107

103

Fig. 5. Computational domain.

addition, all thermo-physical properties of the fluid are assumed constant except the density which varies in so far as the Boussinesq approximation is respected. The dimensionless parameters governing this problem are the Prandtl number Pr and the Rayleigh number Ra defined by

Pr ¼ y=a

and

. Ra ¼ bgðTh  Tc ÞH3 Pr n2 :

In our present numerical test, we describe the main results that characterize the flow and heat transfer. Therefore, these results are presented in terms of maximum horizontal velocity Umax, maximum vertical velocity Vmax, maximum stream function jmax and average Nusselt number Nu. The stream function is defined as follows: Fig. 6. Schematic geometry of the physical problem.

2.3.1. Example of simulation: square differentially heated cavity The physical problem under study and displayed in Fig. 6 concerns a differentially heated square cavity filled with air (Pr ¼ 0.71). Indeed, the vertical active walls (left and right) are maintained isothermal at constant temperatures Th and Tc, respectively and the horizontal walls are perfectly insulated (adiabatics). The fluid is Newtonian, incompressible and viscous dissipation is neglected. In

61*61 81*81 101*101 121*121 141*141 161*161 181*181 201*201 221*221 241*241

umax

vmax

jmax

8.615 8.706 8.769 8.794 8.811 8.813 8.814 8.817 8.819 8.819

64.345 64.343 64.472 64.378 64.313 64.264 64.638 64.617 64.576 64.576

219.821 218.020 220.351 220.824 220.320 219.896 221.093 221.104 220.807 220.807

16.5447 16.5751 16.6231 16.6259 16.6284 16.6305 16.7059 16.7088 16.7122 16.7122

(25)

Table 2 Comparison of laminar flow for Ra ¼ 106. Ra

Ra ¼ 106 Nu

vv vu  vx vy

The heat exchange between the fluid and the enclosure is

103

Table 1 Grid dependence study. N2 ¼ Nx*Ny

V2 jðx; yÞ ¼

104

105

106

Nu umax vmax

jmax Nu umax vmax

jmax Nu umax vmax

jmax Nu umax vmax

jmax

De Vahl Davis [50] R D (%)

Shu et al. [51] R D (%)

Hortmann et al. [52] R D (%)

Present results

1.117 3.649 3.697 1.174 2.238 16.178 19.617 5.071 4.509 34.73 68.59 9.612 8.817 64.63 219.36 16.75

1.117 3.646 3.694 No data 2.244 16.158 19.676 No data 4.520 34.301 68.188 9.631 8.804 63.590 218.320 16.815

No data No data No data No data 2.244 16.180 19.629 No data 4.521 34.740 68.639 No data 8.825 64.837 220.461 No data

1.115 3.658 3.712 1.173 2.246 16.149 19.684 5.043 4.525 34.557 68.624 9.558 8.825 64.549 220.915 16.715

0.17 0.24 0.37 0.08 0.35 0.17 0.34 0.55 0.35 0.49 0.04 0.56 0.09 0.12 0.70 0.20

0.17 0.32 0.48 0.08 0.05 0.04 0.11 0.74 0.63 0.75 0.23 1.50 1.18 0.59

0.08 0.19 0.28 0.08 0.52 0.02 0.00 0.44 0.20

104

M. Jami et al. / International Journal of Thermal Sciences 100 (2016) 98e107

Table 3 Comparison of transitional flow with previous works for the square heated cavity. re  [53] Le Que R D (%)

Ra 107

Nu umax vmax

16.523 148.580 999.236 30.170 30.225 321.876 2222.390 53.760

jmax

108

Nu umax vmax

jmax

Dixit and Babu [23] R D (%) 0.06 0.62 30.09 0.60 0.08 1.65 0.34 0.42

16.790 164.236 701.922 No data 30.506 389.877 2241.374 No data

characterized by the average Nusselt number (Nu) that is defined as follows:

1 Nu ¼  Th  Tc

ZH 0

vT vx

dy

(26)

wall

The averaged Nusselt number along the Y-axis is expressed by:

Nu ¼

1 A

ZA 0

 vq  dY vX wall

Kuznik et al. [24] R D (%) 1.53 10.09 0.48 0.83 18.80 1.19

16.408 148.768 702.029 No data 29.819 321.457 2243.360 No data

Present results 0.76 0.75 0.49 1.44 1.52 1.27

16.533 147.651 698.532 29.986 30.25 316.562 2214.692 53.531

In reason to determine the optimum grid that insures a good compromise between accuracy and computational costs, various simulations are conducted in the present work. Prandtl and Rayleigh numbers were kept constant at 0.71 and 106, respectively. The average Nusselt number, the maximum horizontal velocity, the maximum vertical velocity and the maximal stream function are presented in Table 1 for Ra ¼ 106 as a function of the grid size N2. From this Table, it can be observed that the solutions become grid independent at N2 ¼ 221*221. Consequently, this grid seems to be sufficiently fine for this numerical calculation and it is sufficient to capture the transient flow.

Fig. 7. Streamlines for the laminar flow.

M. Jami et al. / International Journal of Thermal Sciences 100 (2016) 98e107

The transition from the motionless conduction dominated regime to the convection dominated regime takes place after Ra ¼ 103. For Ra  107, the fluid motion takes place in the boundary layer near the walls. In our paper, we present these cases as examples of simulation in order to compare the physical phenomena behavior obtained for the Benchmark problems and to show that our code is efficient for the simulation of the different flow cases. The results obtained are presented in Tables 2 and 3, and are in agreement with those found in the literature [23,24,50e53]. 2.3.2. Laminar flow In the first part of the present numerical work, the different values of Ra (103  Ra  106) number considered are summarized in Table 2, where a comparison is made with results available elsewhere in the literature [50e52]. As can be seen from the relative differences concerning the average Nusselt number Nu, maximum horizontal velocity Umax, maximum vertical velocity Vmax and maximum stream function jmax our results are in good agreement with those listed above. Figs. 7 and 8 show respectively the isotherms and streamlines for the laminar flow case. In all the streamline and isotherm plots, the contour lines correspond to equispaced absolute values of the dimensionless stream function j/a and the dimensionless temperature q.

105

As the Rayleigh number increases, the fluid motion mainly takes place near the differentially heated walls and the flow in the core of the cavity becomes quasi-motionless: these flow features are well captured by the proposed numerical method. One may conclude that the isotherms are more stratified with the increase of Rayleigh number, which means that the heat transfer at lower Rayleigh numbers (Ra ¼ 103) is dominated by the conductive effects (Fig. 7). These effects disappear with the increase of the Rayleigh number and the heat transfer is mainly characterized by the presence of convective effects at higher Rayleigh numbers (Ra ¼ 106). Regarding the streamlines (Fig. 8), the flow is more accelerated with increase of Rayleigh number and bicellular structure of streamlines appears at higher Rayleigh number (Ra ¼ 106). It should be noted that the heat transfer increases as function Rayleigh number. 2.3.3. Transitional flow One of the main interests of this paper is to prove that the transitional flow can be simulated using Multiple Relaxation TimeLattice Boltzmann Equation with D2Q9 model for the fluid velocity variables and D2Q5 for the temperature. Table 3 summarizes the numerical values of Nu, Umax, Vmax and jmax compared with literature results. From this table, the numerical values are in good agreement with those from literature [23,24,53].

Fig. 8. Isotherms for the laminar flow.

106

M. Jami et al. / International Journal of Thermal Sciences 100 (2016) 98e107

Fig. 9. Streamlines for the transitional flow.

Fig. 10. Isotherms for the transitional flow.

The isotherms and streamlines for the transitional flow are shown in Figs. 9 and 10. The temperature field is becoming more and more stratified, with horizontal streamlines. Moreover, reversed flows occur at the upper-left and bottom-right corners, destabilizing the laminar flow: the transitional flow features are well captured by our model. In other words, the most denser isotherms are found near the thermally active walls by increasing the Rayleigh number. This densification of iso-lines explains the appearance of the boundary layer regime (Fig. 9). It is found that the dynamical flow is more accelerated versus the increase of Rayleigh number and multicellular structure appears at large Rayleigh numbers (Fig. 10). 3. Conclusion The aim of this work is to present numerical studies concerning air flow and heat transfer in a two-dimensional square cavity. This study is based on Relaxation Time-Thermal Lattice Boltzmann Equation with the D2Q9 model for fluid velocity variables and the D2Q5 for the temperature. The main conclusions concern the validation and the importance of the double population approach needed to solve the heat transfer problem in the case of the heated

square cavity. The results obtained for different Rayleigh numbers show that both laminar and transitional flow features are well captured by the proposed numerical method and the results provide a good comparison with those available in previous papers. This shows that the developed numerical method is reasonably accurate and robust to solve the heat transfer problem in the case of the heated square cavity for Ra numbers up to 108.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]

U. Frisch, B. Hasslacher, Y. Pomeau, Phys. Rev. Lett. 56 (14) (1986) 1505. G.R. McNamara, G. Zanetti, Phys. Rev. Lett. 61 (1988) 2332. nez, Europhys. Lett. 9 (7) (1989) 663. F.J. Higuera, J. Jime F.J. Higuera, S. Succi, R. Benzi, Europhys. Lett. 9 (4) (1989) 345. H. Chen, S. Chen, W.H. Matthaeus, Recovery Phys. Rev. A 45 (8) (1992) R5339. res, P. Lallemand, Europhys. Lett. 17 (1992) 479. Y.H. Qian, D. d'Humie E. Semma, M. El Ganaoui, R. Bennacer, A.A. Mohamad, Int. J. Therm. Sci. 47 (3) (2008) 201. K.N. Premnath, J. Abraham, J. Comput. Phys. 224 (2) (2007) 539. S. Succi, E. Foti, F. Higuera, Europhys. Lett. 10 (1989) 433. A. Cali, S. Succi, A. Cancelliere, R. Benzi, M. Gramignani, Phys. Rev. A 45 (8) (1992) 5771. A. Ladd, J. Fluids Mech. 271 (1994) 311. O. Filippova, D. Hanel, Comput. Fluids 26 (1997) 697. G.W. Yan, Y.S. Chen, S.X. Hu, Phys. Rev. E 59 (1999) 454.

M. Jami et al. / International Journal of Thermal Sciences 100 (2016) 98e107 [14] G.W. Yan, J.Y. Zhang, Y.H. Liu, Y.F. Dong, Int. J. Numer. Methods Fluids 55 (1) (2007) 41. [15] R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: theory and applications, Phys. Rep. 222 (3) (1992) 145. [16] S.P. Dawson, S.Y. Chen, G.D. Doolen, J. Chem. Phys. 98 (1993) 1514. [17] S. Succi, G. Bella, F. Papetti, J. Sci. Comput. 12 (4) (1997) 395. [18] R.G.M. van der Sman, M.H. Ernst, A.C. Berkenbosch, Int. J. Heat Mass Transf. 43 (4) (2000) 577. [19] Y. Gan, A. Xu, G. Zhang, X. Yu, Y. Li, Phys. A Stat. Mech. Appl. 387 (8e9) (2008) 1721. [20] X.F. Pan, A. Xu, G. Zhang, S. Jiang, Int. J. Mod. Phys. C 18 (11) (2007) 1747. [21] A. Mezrhab, M. Bouzidi, P. Lallemand, Comput. Fluids 33 (4) (2004) 623. [22] P. Lallemand, L.-S. Luo, Phys. Rev. E 68 (2003) 036706. [23] H.N. Dixit, V. Babu, Int. J. Heat Mass Transf. 49 (2006) 727. [24] F. Kuznik, J. Vareilles, G. Rusaouen, G. Krauss, Int. J. Heat Fluid Flow 28 (2007) 862. [25] P. Lallemand, L.-S. Luo, Phys. Rev. E 61 (2000) 6546. res, Phys. Rev. E 68 (2003) 066614. [26] I. Ginzburg, D. d'Humie res, Progress in Astronautics and Aeronautics, vol. 159, AIAA, [27] D. d'Humie Washington, DC, 1992, p. 450. [28] R. Du, B. Shi, X. Chen, Phys. Lett. A 359 (2006) 564. [29] X.D. Niu, C. Shu, Y.T. Chew, Y. Peng, Phys. Lett. A 354 (2006) 173. [30] L. Zheng, B. Shi, Z. Guo, Phys. Rev. E 78 (2008) 026705. [31] A. Mezrhab, M.A. Moussaoui, M. Jami, H. Naji, M. Bouzidi, Phys. Lett. A 374 (2010) 3499. [32] A. Mohamad, A. Kuzmin, A critical evaluation of force term in lattice Boltzmann method, natural convection problem, Int. J. Heat Mass Transf. 53 (2010) 990. [33] K.N. Premnath, M.J. Pattison, S. Banerjee, Dynamic subgrid scale modeling of turbulent flows using lattice-Boltzmann method, Physica A 388 (2009) 2640. [34] A. Karimipour, M.H. Esfe, M.R. Safaei, D.T. Semiromi, S. Jafari, S.N. Kazi, Mixed convection of copper-water nanofluid in a shallow inclined lid driven cavity using the lattice Boltzmann method, Physica A 402 (2014) 150. [35] M. Rahimi-Gorji, O. Pourmehran, M. Hatami, D. Domiri ganji, Statistical optimization of microchannel heat sink (MCHS) geometry cooled by different nanofluids using RSM analysis, Eur. Phys. J. Plus 130 (22) (2015) 1.

107

[36] O. Pourmehran, M. Rahimi-Gorji, M. Gorji-Bandpy, D.D. Ganji, Analytical investigation of squeezing unsteady nanofluid flow between parallel plates by LSM and CM, Alex. Eng. J. 54 (2014) 17. [37] M. Rahimi-Gorji, O. Pourmehran, M. Gorji-Bandpy, D.D. Ganji, An analytical investigation on unsteady motion of vertically falling spherical particles in non-Newtonian fluid by collocation method, Ain Shams Eng. J. 6 (2) (2015) 531e540. http://dx.doi.org/10.1016/j.asej.2014.10.016. res, P. Lallemand, L.-S. Luo, J. Comput. Phys. 172 (2) [38] M. Bouzidi, D. d'Humie (2001) 704. [39] P.L. Bhatnagar, E.P. Gross, M. Krook, Phys. Rev. 94 (3) (1954) 511. [40] M. Bouzidi, M. Firdaouss, P. Lallemand, Phys. Fluids 13 (11) (2001) 3452. [41] M. Jami, A. Mezrhab, M. Bouzidi, P. Lallemand, Int. J. Therm. Sci. 46 (1) (2007) 38. [42] A. Mezrhab, M. Jami, M. Bouzidi, P. Lallemand, Comput. Fluids 36 (2) (2007) 423. [43] M.A. Moussaoui, M. Jami, A. Mezrhab, H. Naji, Heat Mass Transf. 45 (11) (2009) 1373. [44] I. Rasin, S. Succi, W. Miller, J. Comput. Phys. 206 (2005) 453. [45] S.C. Mishra, A. Lankadasu, K.N. Beronov, Int. J. Heat Mass Transf. 48 (2005) 3648. [46] F. Dubois, Equivalent partial differential equations of a lattice Boltzmann scheme, Comput. Math. Appl. 55 (2008) 1441e1449. res, B. Hasslacher, P. Lallemand, Y. Pomeau, J.P. Rivet, [47] U. Frisch, D. d'Humie Complex Syst. 1 (1987) 649. [48] I. Ginzbourg, P.M. Adler, Boundary flow condition analysis for the threedimensional lattice Boltzmann model, J. Phys. Paris. II 4 (1994) 191. [49] P. Lallemand, L.S. Luo, Lattice Boltzmann method for moving boundaries, J. Comput. Phys. 184 (2003) 406. [50] G. De Vahl Davis, Int. J. Numer. Methods Fluids 3 (1983) 249. [51] C. Shu, Y. Peng, Y.T. Chew, Int. J. Mod. Phys. C 13 (2002) 1399. [52] M. Hortmann, M. Peric, G. Scheuerer, Int. J. Numer. Methods Fluids 11 (1990) 189. re , Comput. Fluids 20 (1991) 29. [53] L. Le Que