Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

51, 2, pp. 447-462, Warsaw 2013

THE EFFECTS OF INCLINATION ANGLE AND PRANDTL NUMBER ON

THE MIXED CONVECTION IN THE INCLINED LID DRIVEN CAVITY

USING LATTICE BOLTZMANN METHOD

Arash Karimipour

Department of Mechanical Engineering, Najafabad Branch, Islamic Azad University, Isfahan, Iran

e-mail: arashkarimipour@gmail.com

Alireza Hossein Nezhad

University of Sistan and Baluchestan, Department of Mechanical Engineering, Zahedan, Iran

e-mail: nezhadd@hamoon.usb.ac.ir

Annunziata D’Orazio

Dipartimento di Ingegneria Astronautica, Elettrica ed Energetica, Sapienza Università di Roma, Rome, Italy

e-mail: annunziata.dorazio@uniroma1.it

Ebrahim Shirani

Foolad Institute of Technology, Fooladshahr, Esfahan, Iran

e-mail: eshirani@ictp.it

The laminarmixed convection in a two-dimensional rectangular inclined cavitywithmoving
top lid is investigatedusing thedoublepopulation thermal latticeBoltzmannmethod (LBM)
at different values of the Richardson number, inclination angle and the Prandtl number. In
this problem, velocity components are changed by both buoyancy forces and the inclina-
tion angle of the cavity. Comparison of the present results with other available data shows
good agreement. As the results, the velocity and temperature profiles, the Nusselt num-
ber, streamlines and isotherms are presented and discussed. It is shown that the increase of
Prandtl number enhances the heat transfer rate, especially at higher values of inclination
angle and Richardson number. Moreover, the average Nusselt number at the upper limit
of the considered range of the Richardson and Prandtl numbers variability increases by a
factor of 9.

Key words: LBM, inclination angle, Prandtl number, mixed convection

Nomenclature

AR – cavity aspect ratio (L/H)
cs,e – lattice speed of sound and internal energy
f,g – momentum and internal energy functions
fe,ge – equilibrium distribution functions
g – gravity vector
Gr,Ma,Pr – Grashof, Mach and Prandtl number, respectively
H,L – height and length of the cavity
k – thermal conductivity
Num – average Nusselt number
q – heat flux
R – constant of gas
Ra,Re,Ri – Rayleigh, Reynolds and Richardson number, respectively
t,T – time and temperature, respectively
Tc,Th – cold and hot wall temperature



448 A. Karimipour et al.

u – macroscopic flow velocity vector, u= [u,v]
U0 – top lid velocity
(U,V ) – dimensionless flow velocity, (U,V )= (u/U0,v/U0)
Uw,Vw – velocity components of the cavity walls
x – dimensional Cartesian coordinate vector, x=(x,y)
(X,Y ) – dimensionless coordinates, (X,Y )= (x/H,y/H)
Z – viscous heating term

Greek symbols
α – thermal diffusivity ν – kinematic viscosity
β – volumetric expansion coefficient θ – dimensionless temperature,
ρ – density θ=(T −Tc)/(Th−Tc)
γ – cavity inclination angle τf,τg – relaxation times

1. Introduction

The lattice Boltzmann method is a particle based approach being used for the numerical simu-
lation of fluid flow and heat transfer. The particle characteristic of this method has increased
its application in a wide range of fluid flow and heat transfer problems, so that in addition to
the simulation of macroflows (Grucelski and Pozorski, 2012; Kefayati et al., 2011; Nemati et al.,
2010; Yang andLai, 2011), it is used for the simulation ofmicro andnanoflows (Kandlikar et al.,
2006;Karimipour et al., 2012;Niu et al., 2007; Tian et al., 2010).Moreover, LBMhas foundwide
application in micro-electro-mechanical-systems (MEMS) and nano-electro-mechanical systems
(NEMS).Compared to the conventional numericalmethods and other particle based simulations
such as molecular dynamics simulation and direct simulation Monte Carlo, LBM is more ap-
propriate for parallel processing. Moreover, using LBM, the pressure field is directly calculated
without the need for solving another system of equations, multiphase and complex flows can be
solved easier, and less computational memory and time are needed (Chen et al., 1992; Chen and
Doolen, 1998; Oran et al., 1998).Moreover, LBMconsists of only first-order PDEs, whichmakes
discretization and programming simpler than Navier-Stokes equations which are second-order
PDEs. Moreover, the nonlinear convective term in Navier-Stokes equations is written simpler
in LBM (Tallavajhula et al., 2011). These advantages give incentives to researchers to study
the application of the LBM to solve more realistic problems by improving and innovating the
LBM models and related boundary conditions. However, there are some difficulties and draw-
backs in LBM: it is a compressible model for ideal gas, and theoretically always simulate the
compressible Navier-Stokes equation. However, the incompressible Navier-Stokes equations can
be derived from the LBM through theChapman-Enskog expansion at the nearly incompressible
limit. Itmeans LBMcan simulate an incompressible flowunder lowMach number (Ma< 0.15).
The compressible nature of LBM produces a compressibility error, which at low values of the
Mach number will be negligible (order of Ma2) (Buick and Greated, 2000; He and Luo, 1997;
Mohamad, 2011; Shi et al., 2006). Moreover, the LBMmulti phase model can not simulate the
systems with large viscosity ratio fluids (Kuzmin and Mohamad, 2009). In addition, using the
regular square grids is another difficulty of LBM for simulation of the curved boundaries. Some
researchers have dealt with curved boundaries and using unstructured meshes in LBM (Cheng
and Hung, 2002; Kao and Yang, 2008; Peng et al., 1999). Kao and Yang (2008) applied an
interpolation-based approach (under a uniform Cartesian mesh) to track the position of boun-
dary for solving the distribution functions near the curved boundary. This method results in a
loss of mass conservation and reduces the accuracy at the boundary. Ubertini and Succi (2008)
used non-uniform or unstructured meshes for LBM to improve both stability and accuracy.
However, they reported that further improvements are necessary to obtain accurate results at



The effects of inclination angle and Prandtl number... 449

different flow conditions. To simulate heat transfer, different lattice Boltzmann methods have
been proposed such as themulti-speed, passive scalar, and doubled populations internal energy
method. The last method has been widely used to simulate natural convection problems (He et
al., 1998). Guo et al. (2007) used thermal LBM for solving lowMach number thermal flowswith
viscous dissipation and compression work. They obtained a lattice Boltzmann equation model
from a kinetic model for the decoupled hydrodynamic and energy equations. Their model was
tested by simulating the thermal Poiseuille flow in a planer channel and natural convection in a
square cavity. Natural convection in the inclined cavity using LBMhas been reported in various
articles in the recent years (Jafari et al., 2011; Mezrhab et al., 2006). Numerous investigations
have been conducted in the past on the lid-driven cavity flow and heat transfer, considering
various combinations of the imposed temperature gradients and cavity configurations. Sharif
(2007) studied numerically two-dimensional shallow rectangular driven cavities of aspect ratio
AR=10 for Ra ranging from 105 to 107, keeping the Reynolds number fixed at Re= 408.21.
Basak et al. (2009) performed finite element simulations to investigate the influence of linearly
heated side wall(s) or cooled right wall on mixed convection lid-driven flows in a square cavi-
ty. Sivasankaran et al. (2010) performed a numerical study, with the finite volume method, on
mixed convection in a lid-driven cavity with vertical sidewalls maintained with sinusoidal tem-
perature distribution and top and bottomwall adiabatic; the results were analyzed over a range
of Ri, amplitude ratios and phase deviations. The amplitude ratio was defined as the ratio of
the amplitude of temperature oscillations of the right wall to that of the left wall, and the phase
deviation was defined as the phase difference of temperature oscillations between the right and
left walls. The effects of Prandtl numbers (0.7<Pr< 70) on natural convection in the cavity
using LBM were investigated by Kao and Yang (2007). Satisfying the nearly incompressible
flow (Ma< 0.1), they determined different characteristic velocities at each Pr. Their approach
showed good performance at Pr = 0.7 and Pr = 7; however, it needed more time at higher
Pr values. Their method is usually applied for Pr¬ 7 (Nemati et al., 2010). Parmigiani et al.
(2009) used two supplementary methods for higher values of Pr and Ra (10 < Pr < 104 and
104 <Ra< 109) to simulate natural convection usingLBM. In their firstmethod, the timescales
of thermal and density-momentum distribution functions are separated at higher Pr values. In
their second method, a smaller grid size is used for the thermal distribution function than for
the density-momentumdistribution function. Simulation ofmixed convection using LBMat dif-
ferent conditions has been one of the interesting topics for researchers in the recent years (Guo
et al., 2010; Du et al., 2011; Fattahi et al., 2011). Among them, Rosdzimin et al. (2010) studied
the effects of a heated square inside the lid driven cavity, using the nine-velocity model for the
velocity field and the four-velocity model for the thermal field.

The mentioned review shows that mixed convection in an inclined lid driven cavity has
not been investigated by LBM. So, in this work, for a laminar mixed convection, the effects of
gamma and Pr at different Ri values on the thermal and hydrodynamic fluid properties inside
a two-dimensional inclined enclosure with hot moving top lid are studied.

2. Problem statement

The fluid mixed convection inside an enclosure shown in Fig. 1 (L/H = 3) is studied utilizing
LBM. The upper lid temperature is larger than that of the lower wall. U0 is the upper wall
velocity and the side walls are insulated. A computer program in Fortran language is developed
to simulate the fluid parameters for Re=U0H/ν =200 and Ri=Gr/Re

2 =0.1,1,10, in which
Gr= gβ∆TH3/ν2. Theeffects of γ=0◦, 30◦, 60◦, 90◦ and Pr= να=0.07,0.7,7 are studied on
heat transfer and fluid flow; ν and α are the kinematic viscosity and thermal diffusivity. LBM is
applied in near-incompressible regimes. Thus, in theMach number definition, Ma=U∗/cs ≪ 1,



450 A. Karimipour et al.

the characteristic velocity of the flow for both natural, U∗ =
√
gβ∆TH, and forced convection,

U∗ = νRe/H, must be small compared with the speed of sound. The compressibility errors are
proportional to Ma2, soat lowvalues of theMachnumber,densityvariationswill benegligible. In
the presentwork, Ma is assumedas 0.1; therefore, the compressibility errorswould be negligible.

Fig. 1. Inclined cavity geometry and coordinates axis

3. Formulation

3.1. Lattice Boltzmann method

The dimensionless lattice Boltzmann equations for hydrodynamic and thermal models are
as follows (He et al., 1998)

∂fi
∂t
+ ciα

∂fi
∂xα
=−fi−f

e
i

τf

∂gi
∂t
+ ciα

∂gi
∂xα
=−gi−g

e
i

τg
−fiZi (3.1)

where fi is the discretisedmomentum distribution function and denotes the probability density
of particles having velocity around ciα at an infinitesimal volume element centered at xα. g is
called the internal energy density distribution function. Indices i and α are lattice velocity
directions and x-y geometry components, respectively, τf and τg arehydrodynamicand thermal
relaxation times, respectively, fe and ge arehydrodynamicand thermal equilibriumdistribution
functions, respectively, and ci representsmicroscopicparticles velocity.UsingD2Q9 lattice (Qian
et al., 1992), shown in Fig. 2, the subscript i varies from 1 to 9. So, the microscopic particle
velocities are calculated as follows

ci=0 = [0,0] ci=1,2,3,4 =
[
cos
i−1
2
π,sin

i−1
2
π
]
c

ci=5,6,7,8 =
√
2
[
cos
((i−5)π
2

+
π

4

)
,sin
((i−5)π
2

+
π

4

)]
c

(3.2)

Z is the heat dissipation term defined as

Zi =(ciα−uα)
(δuα
δt
+ ciα

∂uα
∂xα

)
(3.3)

Theequilibriumdistribution functions fedescribe the equilibriumstate of f, andarewritten
as (He et al., 1998)

fei=0,1,...,8 =ωiρ
(
1+
3ci ·u
c2
+
9(ci ·u)2
2c4

− 3u
2

2c2

)

ω0 =
4

9
ω1,2,3,4 =

1

9
ω5,6,7,8 =

1

36

(3.4)

and



The effects of inclination angle and Prandtl number... 451

ge0 =−
2

3
ρe
(u2

c2

)
ge1,2,3,4 =

1

9
ρe
(3
2
+
3

2

c1,2,3,4 ·u
c2

+
9

2

(c1,2,3,4 ·u)2
c4

− 3
2

u2

c2

)

ge5,6,7,8 =
1

36
ρe
(
3+6

c5,6,7,8 ·u
c2

+
9

2

(c5,6,7,8 ·u)2
c4

−
3

2

u2

c2

) (3.5)

where c2 =3RT and R is the constant of gas. The discretized form of Eq. (3.1)1 is written as

fi(x+ci∆t,t+∆t)−fi(x, t)=−
∆t

τf
[fi(x, t)−fei (x, t)] (3.6)

Fig. 2. D2Q9 lattice

The constant value of BGK collision operator results in the second-order truncation error in
the lattice Boltzmann equation. This error is absorbed into the physical viscous term by using
ν =(τf−0.5∆t)RT for isothermalmodels.But, for thermalmodels, theviscosity appears inboth
momentum and energy equations; however the viscosity in the viscous heat dissipation term of
energy equation must be considered as ν = τfRT , which is inconsistent with the viscosity form
in isothermalmodels. To solve this problem, a second-order strategy to integrate theBoltzmann
equation is applied, which leads to the following equations –more details can be found in He et
al. (1998) and Peng et al. (2003)

− fi(x, t)−f
e
i (x, t)

τf
=−∆t
2τf
[fi(x+ci∆t,t+∆t)−fei (x+ci∆t,t+∆t)]

− ∆t
2τf
[fi(x, t)−fei (x, t)]

fi(x+ci∆t,t+∆t)−fi(x, t)=−
∆t

2τf
[fi(x+ci∆t,t+∆t)−fei (x+ci∆t,t+∆t)]

−
∆t

2τf
[fi(x, t)−fei (x, t)]

(3.7)

Using the sameprocedure done onEqs. (3.6) and (3.7) forEq. (3.1)1, the following equations
are obtained for Eq. (3.1)2

− gi(x, t)−g
e
i(x, t)

τg
−fi(x, t)Zi(x, t)=−

∆t

2τg
[gi(x+ci∆t,t+∆t)−gei(x+ci∆t,t+∆t)]

−
∆t

2
fi(x+ci∆t,t+∆t)Zi(x+ci∆t,t+∆t)−

∆t

2τg
[gi(x, t)−gei(x, t)]

−
∆t

2
fi(x, t)Zi(x, t)

gi(x+ci∆t,t+∆t)−gi(x, t)=−
∆t

2τg
[gi(x+ci∆t,t+∆t)−gei(x+ci∆t,t+∆t)]

−
∆t

2
fi(x+ci∆t,t+∆t)Zi(x+ci∆t,t+∆t)−

∆t

2τg
[gi(x, t)−gei(x, t)]

−∆t
2
fi(x, t)Zi(x, t)

(3.8)



452 A. Karimipour et al.

To solve the difficulty of implicitness of Eqs. (3.7)2 and (3.8)2, two distribution functions f̃i
and g̃i are defined. The symbols f̃ and g̃ are defined for the numerical purpose. However, they
indicate the momentum and internal energy distribution functions like f and g, respectively
(He et al., 1998)

f̃i = fi+
∆t

2τf
(fi−fei )−

∆t

2
F g̃i = gi+

∆t

2τg
(gi−gei)+

∆t

2
fiZi (3.9)

In each time step, the collision and propagation stages are performed sequentially between
particles. In BGK model, these stages are stated as follows

f̃i(x+ci∆t,t+∆t)− f̃i(x, t)=−
∆t

τf +0.5∆t
[f̃i(x, t)−fei (x, t)]

g̃i(x+ci∆t,t+∆t)− g̃i(x, t)=−
∆t

τg+0.5∆t
[g̃i(x, t)−gei(x, t)]−

τg∆t

τg +0.5∆t
fiZi

(3.10)

Finally, using f̃i and g̃i, the hydrodynamic and thermal variables are calculated as

ρ=
∑

i

f̃i ρe=
∑

i

g̃i−
∆t

2

∑

i

fiZi ρu=
∑

i

cif̃i

q=
τg

τg+0.5∆t

(∑

i

cig̃i−ρeu−
1

2
∆t
∑

i

cifiZi

) (3.11)

where e=RT is the internal energy and q= [qx,qy] is the heat flux vector. ν and α are stated
as

ν = τfRT α=2τgRT (3.12)

Using the Chapman-Enskog expansion, the continuity and Navier-Stokes equations can be
obtained from LBM-BGK equation (Cercignani, 1998). The characteristic velocity is assumed
to be 0.1 to keep Ma< 1, so the kinematic viscosity is estimated as ν =Re/U0H. Considering
both RT = 1/3 and Eq. (3.12), the value of τf is determined. Now, the thermal diffusivity is
calculated as α= ν/Pr, and finally by Eq. (3.12) it is found that τg =α/2RT .

3.2. Gravity effects

The Boussinesq approximation is used as ρ= ρ[1−β(T −T)] where ρ, T are the reference
fluid density and temperature. The buoyancy force per unit mass is defined as G=βg(T −T)
and F =G ·(c−u)fe/RT in Eq. (3.13) refers to the buoyancy force effects in this problem (He
et al., 1998; Kuznik et al., 2007)

∂tf+(c ·∇)f =−
f−fe
τf
+F (3.13)

Now the discretized Boltzmann equation is written as follows (D’Orazio et al., 2004)

∂tfi+(ci ·∇)fi =−
fi−fei
τf

+
G · (ci−u)
RT

fei (3.14)

where G= [Gx,Gy], Gx = β|g|(T −T)sinγ =Gsinγ, and Gy = β|g|(T −T)cosγ =Gcosγ.
Using the same method as in Section 3.1 and by substituting u= [u,v] and ci = [cix,ciy], the
following equations are obtained

f̃i(x+ci∆t,t+∆t)− f̃i(x, t)=−
∆t

τf +0.5∆t
(f̃i−fei )

+
( ∆tτf
τf +0.5∆t

3G(cix−u)
c2

fei

)
sinγ+

( ∆tτf
τf +0.5∆t

3G(ciy −v)
c2

fei

)
cosγ

(3.15)



The effects of inclination angle and Prandtl number... 453

Equation (3.9), which includes the external force term F, is written as the following (He et
al., 1998)

f̃i = fi+
∆t

2τf
(fi−fei )−

∆t

2
F ⇒ fi =

τff̃i+0.5∆tf
e
i

τf +0.5∆t
+
0.5∆tτf
τf +0.5∆t

F

fi =
τff̃i+0.5∆tf

e
i

τf +0.5∆t
+
0.5∆tτf
τf +0.5∆t

G · (ci−u)
RT

fei

(3.16)

Equation (3.16)2 in the x and y directions is stated as

fi =
τff̃i+0.5∆tf

e
i

τf +0.5∆t
+
( 0.5∆tτf
τf +0.5∆t

G(cix−u)
RT

fei

)
sinγ+

( 0.5∆tτf
τf +0.5∆t

G(ciy −v)
RT

fei

)
cosγ

(3.17)

Using the sameprocedureand considering the effects of gravity and γ, the following formulae
are derived to calculate the macroscopic hydrodynamic variables

ρ=
∑

i

f̃i u=
1

ρ

∑

i

f̃icix+
∆t

2
Gsinγ v=

1

ρ

∑

i

f̃iciy +
∆t

2
Gcosγ (3.18)

3.3. Hydrodynamic boundary conditions

The non-equilibrium bounce back model is used to simulate the no-slip boundary condition
on the walls (Zou and He, 1997). This model improves accuracy compared to the bounce back
boundarycondition and satisfies the zeromassflowrate atnodeson thewall.The collision occurs
on the nodes located at the solid-fluid boundaries. The distribution functions are reflected in
suitable directions, satisfying the equilibrium conditions (He et al., 1998). As an example for
the west wall, the known populations are f̃0, f̃2, f̃3, f̃4, f̃6, f̃7, and after collision to the wall
nodes, the unknown populations will be f̃1, f̃5, f̃8. Using Eqs. (3.18), the following equations
are obtained

ρ=
∑

i

f̃i ⇒ f̃1+ f̃5+ f̃8 = ρw− (f̃0+ f̃2+ f̃3+ f̃4+ f̃6+ f̃7)

u=
1

ρ

∑

i

f̃icix+
∆t

2
Gsinγ ⇒ f̃1+ f̃5+ f̃8 = ρwUw+(f̃3+ f̃6+ f̃7)−

∆t

2
ρwGsinγ

(3.19)

v=
1

ρ

∑

i

f̃iciy +
∆t

2
Gcosγ ⇒ f̃5− f̃8 = ρwVw+(−f̃2+ f̃4− f̃6+ f̃7)−

∆t

2
ρwGcosγ

where ρw and Uw are density and velocity at thewall nodes, respectively. Using the bounceback
rule for the non-equilibrium part of the particle distribution normal to the boundary and Eqs.
(3.4) (Zou and He, 1997) gives

f̃1− f̃e1 = f̃3− f̃e3 ⇒ f̃1 = f̃3+
2

3
ρwUw (3.20)

substitutingEq. (3.20) inEqs. (3.19)2,3, and then adding and subtracting the resulting equations
leads to

f̃8 = f̃6−
f̃4− f̃2
2
+
1

6
ρwUw−

1

2
ρwVw+

∆t

4
ρwG(cosγ− sinγ)

f̃5 = f̃7+
f̃4− f̃2
2
+
1

6
ρwUw+

1

2
ρwVw−

∆t

4
ρwG(cosγ+sinγ)

(3.21)

Equations (3.20) and (3.21) are presented to include the effects of gravity and inclination an-
gle for theno sliphydrodynamicboundarycondition on thewestwall, and the rest corresponding
equations are written similarly for other walls and corners.



454 A. Karimipour et al.

3.4. Thermal boundary conditions

The top lid and bottom wall of the cavity are maintained at constant but different tempe-
ratures, Th and Tc, respectively, and the sidewalls are insulated. The general purpose thermal
boundary condition (GPTBC) is used to implement the constant thermal boundary condition
on the top and bottom walls. This model was developed by D’Orazio et al. (2003) and (2004),
based on the non-equilibrium bounce back boundary condition of Zou and He (1997) and He
et al. (1998). In this model, the unknown thermal populations are assumed to be equilibrium
distribution functions with a counter slip thermal energy density ρe′, which is determined so
that suitable constraints are verified. For example, for the topmoving lid

g̃4,7,8 = ρ(e+e
′)
ge4,7,8
ρe

ρe′ =2ρe+
3

2
∆t
∑

i

fiZi−3K (3.22)

K represents the sum of the six known thermal distribution functions of the neighboring nodes
and e is the imposed thermal energy density at the wall. The GPTBC shows suitable stability
and accuracy for different boundary conditions. Finally, the unknown distribution functions
reflecting from the top wall g̃4, g̃7 and g̃8 are selected as follows (D’Orazio et al., 2003, 2004)

g̃7 =

(
3ρe+

3

2
∆t
∑

i

fiZi−3(g̃0+ g̃1+ g̃2+ g̃3+ g̃5+ g̃6)
)
(3.0−6U0+3U20)

1

36

g̃4 =

(
3ρe+

3

2
∆t
∑

i

fiZi−3(g̃0+ g̃1+ g̃2+ g̃3+ g̃5+ g̃6)
)(3
2
−
3

2
U20

)1
9

g̃8 =

(
3ρe+

3

2
∆t
∑

i

fiZi−3(g̃0+ g̃1+ g̃2+ g̃3+ g̃5+ g̃6)
)
(3.0−6U0+3.0U20)

1

36

(3.23)

Similar procedure is performed for the cold bottom wall. It is noted that for the implemen-
tation of the adiabatic boundary condition on the side walls, the condition qx = 0 should be
substituted for q in Eq. (3.11). For example for thewestwall, the following equation is obtained
(D’Orazio et al., 2003, 2004)

∑

i

cixg̃i =
1

2
∆t
∑

i

cixfiZi (3.24)

Using Eqs. (3.22) and (3.24) for i=1,5,8, the unknown distributions are selected as

ρe′ =3(g̃6+ g̃3+ g̃7)+
3

2
∆t
∑

i

cix
c
fiZi−ρe (3.25)

The Nusselt numbers along the top and bottom walls are estimated as follows

NuX =−
(∂θ
∂Y

)

Y=0,1
Num =

1

AR

AR∫

0

NuX dX (3.26)

4. Results

The effects of γ and Pr on the fluid inside the inclined cavity shown in Fig. 1 with Re= 200
at different Ri are studied using the lattice Boltzmann method. The top and bottom walls are
at Y =1, 0 respectively, whereas the side walls are considered at X =0, 3.



The effects of inclination angle and Prandtl number... 455

4.1. Grid independency and validation

In order to obtain a grid independent solution, a grid refinement study is performed for a
horizontal cavity (γ = 0). Grid independence of the results has been established in terms of
Num on the lid and dimensionless values of x-velocity U, y-velocity V , and temperature θ at
X =1.5 and Y =0.5 (cavity centre) for three different grid sizes, namely 300×100, 450×150
and 600× 200 lattice nodes. In Table 1, the results are reported as obtained for Ri = 0.1,
Re= 200 and Pr=0.7; because of small differences, the 450×150 grid is selected to continue
the calculations.

To validate the computer code, three cases are examined. The first one is a benchmark
numerical solutionof a free convection square cavityflowwith sidewalls atdifferent temperatures
and horizontal walls adiabatic, obtained by Davis (1983). The results for different Ra ranging
from 104 to 105, and Pr = 0.7, are reported in Table 2 in which V ∗ = ν/PrH represents the
diffusion velocity. Table 2 shows the maximum horizontal velocity umax/V

∗ at x/L=0.5, the
maximumvertical velocity vmax/V

∗ at y/L=0.5 and their locations. Num at the hot-wall also
is reported.The second case is amixed convection problem, investigated by Iwatsu et al. (1993);
it concerns a square cavity, heated from the topmoving wall and cooled from the bottom, with
adiabatic sidewalls. Comparisons of U and T profiles along the vertical centreline for Gr=102

and Re=400 are shown in Fig. 3. As the last case for validation, themixed convection of fluid
flowandheat transfer in a vertical channel (x-direction) is studied, and the results are compared
with those of Habchi and Acharya (1986). The right wall temperature is Th (hot temperature)
at y = 0 and the left wall is assumed adiabatic at y = 1. A hot block with length L is also
attached to the right wall. The inlet fluid temperature is Tc (cold temperature). Dimensionless
temperature profiles at the channel cross-section x/L = 0.77 for Pr = 0.7, Ra = 105 and
different values of Ri, are presented in Fig. 4, and good agreement is seen.

Table 1.Grid study for Ri= 0.1, Re=200, Pr=0.7 at X =1.5 and Y =0.5

Parameters
Mesh

300×100 450×150 600×200
U −0.197 −0.195 −0.194
V 0.063 0.066 0.067

θ 0.560 0.564 0.567

Num 2.331 2.367 2.382

Table 2.Comparisonof themaximumhorizontal velocity umax/V
∗ at x/L=0.5, themaximum

vertical velocity vmax/V
∗ at y/L=0.5, and their locations obtained from present results with

those of Davis (1983)

Ra
umax/V

∗ and (y/L) vmax/V
∗ and (x/L) Num

Present Davis error % Present Davis error % Present Davis error %

104
15.951 16.178 −1.403 19.338 19.617 −1.422

2.210 2.243 −1.471
(0.817) (0.823) (−0.729) (0.123) (0.119) (3.361)

105
34.239 34.730 −1.414 67.501 68.590 −1.588

4.456 4.519 −1.394
(0.851) (0.855) (−0.468) (0.067) (0.066) (1.515)

4.2. Effects of cavity inclination angle

In order to show the effects of γ on the flowfield and heat transfer, in Fig. 5 streamlines and
isotherms are reported at different γ = 0◦, 30◦, 60◦ and 90◦ for the case Ri = 1 at Pr = 0.7,



456 A. Karimipour et al.

Fig. 3. Comparison of U and T along the cavity vertical centreline with Iwatsu et al. (1993)

Fig. 4. Dimensionless temperature profiles with those of Habchi and Acharya (1986)

respectively. A clockwise rotational cell is produced in the fluid flowbecause of cavity lidmotion
which transports the hot fluid to the lower parts of the cavity space. The desired pressure
gradient in the vertical direction is made due to such hot fluid motion and, consequently, the
buoyancy forces are generated to push the hot fluid to the upper parts again. Combination of
these twomechanisms of heat transfer, resulting from lidmotion and buoyancy forces, is named
as mixed convection. The Richardson number is defined as Ri = Gr/Re2; it means that for
Ri≪ 1 and Ri≫ 1, the forced and free convection are the dominant heat transfermechanisms,
respectively, and for Ri 1 the mixed convection is considered.

At forced convection domination (Ri = 0.1), there will be only one powerful cell which
will covere almost all the cavity space, and the increasing γ leads to a slight increase in the
power of such a cell, and no more other important effects are seen. Figure 5 shows the two
cells affecting the fluid flow as Ri = 1, however the upper one is grater than the lower one.
At a larger inclination angle, these two cells merge, so that for γ = 90◦, there will be a large
strong cell which covers the whole cavity space. The inclination angle has significant effects on
the thermal and hydrodynamic fluid parameters when natural convection dominates the cavity
space. The straight isotherms, which are almost perpendicular to the sidewalls, can be seen in
this case as γ =0◦, especially in the lower half of the enclosure spacewhich shows the conduction
heat transfer in this region. Figure 6 shows the dimensionless horizontal velocity profile U and
dimensionless temperatureprofile θ along the cavity vertical centreline at x/H =1.5 for Ri= 1,
Pr=0.7 and different inclination angles.

At Y =0, U is zero andat Y =1, it approaches the lid velocity. A larger γ corresponds to a
larger absolute value of U at 0<Y < 0.3. The larger γ corresponds to the larger absolute value
of U at 0 <Y < 0.3; it means faster fluid movement in the lower part of the cavity with the



The effects of inclination angle and Prandtl number... 457

Fig. 5. Streamlines (left) and isotherms (right) for Ri=1, at γ=0◦, 30◦, 60◦ and 90◦, Pr=0.7

Fig. 6. Profiles of U and θ at x/H =1.5 for Ri=1 and Pr=0.7

increasing inclination angle. At γ =0, the temperature profile changes almost linearly from zero
(cold wall temperature at Y = 0) to the hot wall temperature at Y = 1. A larger inclination
angle leads to less temperature differences in the central region of the cavity, which shows the
core of the rotational cell. Thin thermal boundary layers along the top and bottom walls can
be seen due to high temperature variations close to these walls. At γ = 60◦ and γ = 90◦, the
temperature values at Y =0.3 are higher than the corresponding values at Y =0.75; however
the bottom cold wall is closer to the region at Y = 0.3. This physical phenomenon illustrates
the desired temperature gradient to generate the buoyancy forces in these regions. At higher Ri,
the larger γ increases the absolute value of U adjacent to the upper and lower walls. Thus, in
spite of previous articles (Iwatsu et al., (1993), the Umax value can be higher than the moving
lid velocity (at Ri = 10 and γ = 90◦). This phenomenon occurs because of the effects of both
buoyancy forces and lid motion.



458 A. Karimipour et al.

4.3. Effects of Prandtl number

The effects of Pr = 0.7 were studied in Section 4.2, and in this section, the effects of
Pr=0.07, 7 are studied. Figure 7 shows the profiles of U and θ for Ri= 0.1 at Pr=7.

Fig. 7. Profiles of U and θ at x/H =1.5 for Ri=0.1 and Pr=7

It is seen that the largest temperature variations occur in the lower part of the cavity
(Y < 0.25), which is different from the related results in the previous section. This cavity
space is far from the top lid, so its properties depend more on the buoyancy forces than the
lid motions. The increasing γ leads to more significant effects on temperature profiles in this
region. Figure 8 shows profiles of U and θ for Ri= 10 at Pr=7. A larger Pr corresponds to a
larger Umax in the state of natural convection dominance. Figure 9 shows Num on the hot wall
as a function of γ for Ri= 0.1, 1, 10 and Pr= 0.07, 0.7, 7. It has to be noted that for γ =0,
Num decreases when Ri increases, it implies the weak contribution of natural convection in the
horizontal configuration. With regard to γ 6= 0, for Ri = 0.1, Num increases with γ slightly;
and for Ri­ 1 it would increase more intensively.

Fig. 8. Profiles of U and θ at x/H =1.5 for Ri= 10 and Pr=7

A larger Pr corresponds to a larger Num, especially at higher values of Ri and γ. Num at
Ri= 10 and Pr=7 increasesmore stronglywith the increasing of γ=0◦ to γ=90◦.Moreover,
for engineering applications and showing the physical effects of the parameters, the following
single non-linear correlation is obtained to estimate the average Nusselt number as a function
of Ri, Pr and γ, used for 0.1¬Ri¬ 10, 0.07¬Pr¬ 7, 0¬ γ¬π/2 (γ in radians), Re=200
and AR=3. The average deviation of this correlation is 4.5%

Num =0.2364+2.957Ri
0.2522Pr0.3277γ0.5074+1.637Ri−0.1629Pr0.3344−0.8402γ0.4872 (4.1)



The effects of inclination angle and Prandtl number... 459

Fig. 9. Nu
m
on the hot wall as a function of γ ? for different Ri and Pr

5. Conclusion

As an alternativemethod, a thermal LBM-BGKwas developed to study the fluidmixed convec-
tion in an inclined enclosure. In this problem, gravity effects and the inclination angle changed
the velocity components resulted from the moving top lid. In order to use LBM, the collision
term of the Boltzmann equation and the calculation procedure of the macroscopic properties
and hydrodynamic boundary conditions weremodified so that the buoyancy forces and the inc-
lination angle could be incorporated properly into the solution process. Moreover, a correlation
was introduced for the average Nusselt number as a function of Ri, Pr and γ.

The results show that with increasing γ and for Pr = 0.7, the thermal and hydrodynamic
flow parameters change more. For Ri = 10, γ = 60◦ and γ = 90◦, the absolute U adjacent
to the upper and lower walls could be more than the U0. At natural convection dominance,
the variations of the inclination angle affect more the fluid flow and heat transfer rate. For the
inclined cavity, at γ = 0◦, the maximum value of Num is obtained at Ri = 0.1, but at larger
inclination angles it occurs at Ri = 10. At larger Pr values, the heat transfer rate is more
sensitive to variations of γ at natural convection dominance. Num increases with the increase
of Pr, especially at higher values of Ri and γ. Num at Ri = 10 and Pr = 7 increases more
strongly (by a factor of 9) as the inclination angle increases from γ =0◦ to γ =90◦. A higher
heat transfer rate occurs at larger Pr values; however, to obtain the higher value of Num at
γ = 0◦ (horizontal cavity) and γ ­ 30◦, the dominant mechanism of heat transfer must be
forced and natural convection, respectively.

To increase the heat transfer rate in the horizontal cavity, forced convection must be the
dominant heat transfer mechanism. In this state, increasing the lid velocity results in an appro-
priate heat transfer growth. However, for the inclined and vertical cavity, the maximum values
of heat transfer occur at the state of natural convection dominance. Thus, in this state, the
increase of buoyancy forces leads to the increase of heat trasfer. To achieve an increased heat
transfer rate, physical geometry and fluid Pr valuemust be changed.A higher heat transfer rate
occurs at larger Pr values. This phenomenon is significant at larger values of the inclination
angles for the state of natural convection dominance. It is recommended to use the cavity in the
vertical position at the state of natural convection dominance to obtain the larger heat transfer
rate.

References

1. Basak T., Roy S., Sharma P.K., Pop I., 2009, Analysis of mixed convection flows within a
square cavity with linearly heated side wall(s), Int. J. of Heat and Mass Transfer, 52, 2224-2242



460 A. Karimipour et al.

2. BhatnagarP.L.,Gross E.P., KrookM., 1954,Amodel for collision process in gases. I. Small
amplitude processes in charged and neutral one-component system,Phys. Rev., 94, 511-1954

3. Bogusławski L., 2007, Estimation of the influence of inflow turbulence on heat convection from
a sphere surface, Journal of Theoretical and Applied Mechanics, 45, 505-511

4. Buick J.M., Greated C.A., 2000, Gravity in a lattice Boltzmann model, Phys. Review E, 61,
5307-5319

5. CercignaniC., 1988,TheBoltzmannEquationand itsApplications,AppliedMathematicalScien-
ces 67, Springer-Verlag, NewYork

6. Chen H., Chen S., Mathaaeus W.M., 1992, Recovery of the Navier-Stokes equations using a
lattice-gas Boltzmannmethod,Physical Review A, 45, 5339-5342

7. Chen S., Doolen G., 1998, Lattice Boltzmannmethod for fluid flows,Annual Rev. Fluid Mech.,
30, 329-364

8. Cheng M., Hung K.C., 2002, Lattice Boltzmannmethod on nonuniformmesh,Recent Advances
in Computational Science And Engineering, 196-199

9. Davis G.V., 1983, Natural convection of air in a square cavity: a benchmark numerical solution,
Int. J. Numer. Methods Fluids, 3, 249-264

10. D’Orazio A., Corcione M.,Celata G., 2004, Application to natural convection enclosed flows
of a lattice Boltzmann BGK model coupled with a general purpose thermal boundary condition,
Int. J. of Thermal Science, 43, 575-586

11. D’Orazio A., Succi S., Arrighetti C., 2003, Lattice Boltzmann simulation of open flowswith
heat transfer,Phys. of Fluids, 15, 2778-2780

12. Du H.Y., Chai Z.H., Shi B.C., 2011, Lattice Boltzmann study of mixed convection in a cubic
cavity,Commun. Theor. Phys., 56, 144-150

13. Fattahi E., Farhadi M., Sedighi K., 2011, Lattice Boltzmann simulation of mixed convection
heat transfer in eccentric annulus, Int. Communications in Heat andMass Transfer, 38, 1135-1141

14. Grucelski A., Pozorski J., 2012, Lattice Boltzmann simulation of fluid flow in porous media
of temperature-affected geometry, Journal of Theoretical and Applied Mechanics, 50, 193-214

15. GuoY., BennacerR., Shen S., Ameziani D., BouzidiM., 2010, Simulation ofmixed convec-
tion in slender rectangular cavity with lattice Boltzmann method, Int. J. of Numerical Methods
for Heat and Fluid Flow, 20, 130-148

16. Guo Z., Zheng C., Shi B., Zhao T.S., 2007, Thermal lattice Boltzmann equation for lowMach
number flows: Decoupling model,Phys. Rev. E, 75 (036704) 1-15

17. Habchi S., Acharya S., 1986, Laminarmixed convection in a partially blocked, vertical channel,
Int. J. Heat Mass Transfer, 29, 1711-1722

18. He X., Luo L.S., 1997, Lattice BoltzmannModel for the Incompressible Navier-Stokes Equation,
J. of Statistical Phys., 88, 927-944

19. He X., Chen S., Doolen G., 1998, A novel thermal model for the lattice Boltzmannmethod in
incompressible limit, J. of Comp. Phys., 146, 282-300

20. Iwatsu R., Hyun J.M., Kuwahara K., 1993,Mixed convection in a driven cavity with a stable
vertical temperature gradient, Int. J. Heat Mass Transfer, 36, 1601-1608

21. JafariM.,NaysariA.,BodaghiK., 2011,LatticeBoltzmannSimulation ofNaturalConvection
Heat Transfer in an Inclined Open Ended Cavity, World Academy of Science Engineering and
Technology, 78, 493-498

22. Kandlikar S.G., Garimella S., Li D., Colin S., King M., 2006, Heat Transfer and Fluid
Flow in Minichannels and Microchannels, First ed., Britain: Elsevier

23. Kao P.H., Yang R.J., 2007, Simulating oscillatory flows in Rayleigh-Bénard convection using
the lattice Boltzmannmethod, Int. J. of Heat and Mass Transfer, 50, 3315-3328



The effects of inclination angle and Prandtl number... 461

24. Kao P.H., Yang R.J., 2008, An investigation into curved and moving boundary treatments in
the lattice Boltzmannmethod, Journal of Computational Physics, 227, 5671-5690

25. Karimipour A., Hossein Nezhad A., D’Orazio A., Shirani E., 2012, Investigation of the
gravity effects on the mixed convection heat transfer in a microchannel using lattice Boltzmann
method, Int. Journal of Thermal Sciences, 54, 142-152

26. Kefayati GH.R., Hosseinizadeh S.F., Gorji M., Sajjadi H., 2011, Lattice Boltzmann simu-
lation of natural convection in tall enclosures using water/SiO2 nanofluid, Int. Communications in
Heat and Mass Transfer, 38, 798-805

27. Kuzmin A., Mohamad A., 2009,Multiphase Simulations with Lattice Boltzmann Scheme, Ph.D.
thesis, University of Calgary, Alberta

28. Kuznik F., Vareilles J., Rusaouen G., Krauss G., 2007, A double-population lattice Bolt-
zmannmethod with non-uniformmesh for the simulation of natural convection in a square cavity,
International Journal of Heat and Fluid Flow, 28, 862-870

29. Mezrhab A., Jami M., Abid C., Bouzidi M., Lallemand P., 2006, Lattice-Boltzmann mo-
delling of natural convection in an inclined square enclosure with partitions attached to its cold
wall, Int. J. of Heat and Fluid Flow, 27, 456-465

30. Mohamad A.A., 2011, Lattice Boltzmann Method Fundamentals and Engineering Applications
with Computer Codes, Springer, Canada

31. Nemati H., Farhadi M., Sedighi K., Fattahi E., Darzi A.A.R., 2010, Lattice Boltzmann
simulation of nanofluid in lid-driven cavity, Int. Communications in Heat and Mass Transfer, 37,
1528-1534

32. NiuX.D., ShuC., ChewY.T., 2007,A thermal lattice Boltzmannmodel with diffuse scattering
boundary condition for micro thermal flows,Comput. Fluids, 36, 273-281

33. Oran E.S., Oh C.K., Cybyk B.Z., 1998, Direct simulation Monte Carlo: recent advances and
applications,Ann. Rev. Fluid Mech., 30, 403-441

34. Parmigiani A., Huber C., Chopard B., Latt J., Bachmann O., 2009, Application of the
multi distribution function lattice Boltzmann approach to thermal flows, Eur. Phys. J. Special
Topics, 171, 37-43

35. Peng, Y., Shu, C., Chew, Y.T., 2003, Simplified thermal lattice Boltzmann model for incom-
pressible thermal flows,Physical Review E, 68, 026701-1-8

36. PengG., Xi H., DuncanC., Chou S.H., 1999, Finite volume scheme for the lattice Boltzmann
method on unstructuredmeshes,Phys. Rev. E, 59, 4675-4682

37. Qian Y., Humières D., Lallemand P., 1992, Lattice BGKmodels for Navier-Stokes equation,
Europhys. Lett., 17, 479-484

38. Rosdzimin A.R.M., Zuhairi S.M., Azwadi C.S.N., 2010, Simulation of mixed convective heat
transfer using lattice Boltzmann method, Int. J. of Automotive and Mechanical Engineering, 2,
130-143

39. Sharif M.A.R., 2007, Laminar mixed convection in shallow inclined driven cavities with hot
moving lid on top and cooled from bottom,Applied Thermal Engineering, 27, 1036-1042

40. Shi Y., Zhao T.S., Guo Z.L., 2006, Lattice Boltzmann method for incompressible flows with
large pressure gradients,Phys. Review E, 73, 026704-1-11

41. Sivasankaran S., Sivakumar V., Prakash P., 2010, Numerical study onmixed convection in
a lid-driven cavity with non-uniform heating on both sidewalls, Int. J. of Heat andMass Transfer,
53, 4304-4315

42. TallavajhulaA.,Kharagpur I.,RuedeU.,BartuschatD., 2011, Introduction to the Lattice
Boltzmann Method, 10th Indo-GermanWinter Academy

43. Tian Z., Chen S., Zheng C.G., 2010, Lattice Boltzmann simulation of gaseous finite-Knudsen
microflows, Int. J. Mod. Phys., 21, 769-783



462 A. Karimipour et al.

44. Ubertini S., Succi S., 2008, A generalised lattice Boltzmann equation on unstructured grids,
Communications in Computational Physics, 3, 342-356

45. YangY.T., Lai F.H., 2011,Numerical study of flow and heat transfer characteristics of alumina-
water nanofluids in a microchannel using the lattice Boltzmann method, Int. Communications in
Heat and Mass Transfer, 38, 607-614

46. Zou Q., He X., 1997, On pressure and velocity boundary conditions for the lattice Boltzmann
BGKmodel,Phys. Fluids, 9, 1591-1598

Analiza za pomocą siatki Boltzmanna wpływu kąta pochylenia oraz liczby Prandtla na

mieszaną konwekcję w ukośnej szczelinie domkniętej ruchomą pokrywą

Streszczenie

W pracy zajęto się problemem mieszanej konwekcji laminarnej w dwuwymiarowej, prostokątnej i
ukośnie usytuowanej szczelinie domkniętej od góry ruchomą pokrywą.Wbadaniach zastosowanometodę
siatki termicznej Boltzmanna (LBM) podwójnej populacji, uwzględniając różne wartości liczby Richard-
sona, kąta pochylenia szczeliny oraz liczby Prandtla. W rozważanym zagadnieniu, składowe prędkości
zostały poddane zmianom indukowanym siłami wyporu oraz kątem pochylenia szczeliny. Porównanie
otrzymanychwyników analizy z dostępnymi w literaturze danymi wykazało dobrą zgodność. Rezultatem
badańw pracy są także profile rozkładu prędkości i temperatury, liczbaNusselta, linie prądu oraz izoter-
my, które szczegółowo przedyskutowano. Pokazano, że wzrost liczby Prandtla zwiększa transfer ciepła,
zwłaszcza dla wyższych wartości kąta pochylenia szczelin i liczby Richardsona. Co więcej, średnia liczba
Nusselta przy górnych wartościach przyjętego zakresu zmienności liczb Richardsona i Prandtla wzrasta
9-krotnie.

Manuscript received December 22, 2011; accepted for print September 26, 2012