International Journal of Heat and Mass Transfer 115 (2017) 854–870

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer

journal homepage: www.elsevier.com/locate/ijhmt

Mixed convective heat transfer from a permeable square cylinder:

A lattice Boltzmann analysis
T.R. Vijaybabu, K. Anirudh, S. Dhinakaran ⇑
The Centre for Fluid Dynamics, Discipline of Mechanical Engineering, Indian Institute of Technology Indore, Khandwa Road, Simrol, Indore, Madhya Pradesh 453 552, India

a r t i c l e i n f o a b s t r a c t

Article history: The flow and mixed convection heat transfer from a two-dimensional porous square cylinder under the
Received 10 February 2017 influence of aiding buoyancy in an infinite stream are analysed employing a mesoscopic approach.
Received in revised form 28 June 2017 Reynolds number (Re) and Darcy number (Da) considered in this study vary from 2 to 40 and 10
6 to
Accepted 11 August 2017
10
2 , respectively. The flow and heat transfer characteristics at the Prandtl number value of 0.71 (air)
is compared for three different values of Richardson number (Ri) i.e. 0, 0.5 and 1. The numerical exper-
iments in this study are carried out by using Lattice Boltzmann technique with two distribution func-
Keywords:
tions. The BGK collision operator with Darcy-Forchheimer and Boussinesq force terms are added to the
Square porous cylinder
Mixed convection
LB collision equation. Mach number annealing process is also carried out to accelerate the simulations.
Steady flow Flow and heat transfer characteristics are found to be a function of non-dimensional permeability (Da),
Darcy-Forchheimer force term buoyancy condition and Reynolds number. It is observed that a monotonous reduction occurs in the wake
LBM length and drag coefficient values at higher permeability levels. Whereas, aiding buoyancy depicts a pro-
D2Q9 lattice model nounced reduction in wake length and an increment in drag coefficient values. The heat transfer
enhancement ratio for all surfaces of the cylinder and mean Nusselt number were calculated to compare
the thermal behaviour at various Ri and Da values. A significant augmentation in heat dissipation is
reported for increasing values of Ri and/or Da. The percentage increment in mean Nusselt number at
Re = 40, Da ¼ 10
2 is found to be 18% and 34% for Ri = 0.5 and 1, respectively with reference to the forced
convection case. Also, heat transfer is maximum at Da ¼ 10
2 and Ri = 1 for the flow regime considered in
this study. Correlations for mean Nusselt number, valid for the range of parameters considered in the pre-
sent study, are also provided. The key results obtained from this study can be helpful for further research
in different realms of engineering sciences, especially thermal engineering, aided by porous media mod-
eling approach.
Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction providing more surface area for heat dissipation in the same vol-
ume compared to solid bodies. A vast amount of books and articles,
Diverse engineering applications can be modeled or approxi- depicting typical case studies for application and numerical tools
mated as transport through porous media. Such instances can be based on porous medium, are available in the literature [1–11].
found in packed bed heat exchangers, drying technology, catalytic Additionally, recent research works related to porous medium
reactors, thermal insulation, petroleum industries and electronic [12,13] have shown its ability on heat transfer magnification.
cooling. These modeling techniques can be an invaluable tool for Another prime advantage with the porous media modeling tech-
the designers in different fields of engineering. Transport phenom- nique is that, it can be used to simplify a complex system with
ena in porous media has been predominantly used in the field of numerous elements into a single porous object. This empowers
geothermal engineering such as underground flow and filtration. the engineer to discover experimentally unapproachable problems
However, more recently, attention is being paid for the application with significant computational economy. A relevant example
of porous media theory in thermal engineering, specifically for the would be of a nuclear reactor core, wherein the arrangement of
purpose of heat transfer enhancement in compact devices. This fuel and control rods resemble a porous square system, as shown
occurrence accounts due to the inherent property of porous bodies, in Fig. 1. Such arrangements can be tuned to receive required flow
and heat transfer characteristics. Thus, numerous research
attempts are in progress to capture enhancement in the heat trans-
⇑ Corresponding author.
fer performance by using porous media.
E-mail address: [email protected] (S. Dhinakaran).

http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.08.033
0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.
 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870 855

Nomenclature

gb4hH
Notations Ri Richardson number, v 20
AB front face of the cylinder T dimensional temperature, [ C]
BC right face of the cylinder u dimensional x-component velocity, [m s
1 ]
CD rear face of the cylinder v dimensional y-component, [m s
1 ]
DA left face of the cylinder u non-dimensional x-component velocity, [m s
1 ]
CD coefficient of drag, FD
0:5qv 20
v non-dimensional y-component velocity, [m s
1 ]
V auxiliary velocity, [m s
1 ]
cF non-dimensional Forchheimer form-drag coefficient
wi weighing factor in direction i
cs speed of sound, [m s
1 ]
x ; y dimensional horizontal & vertical coordinate
K
Da Darcy number, H2 x; y non-dimensional horizontal & vertical coordinate
dp particle diameter, [m]
ei discrete lattice velocity in direction i; dx
dt
i
Greek symbols
Fb Boussinesq force term, [N] q fluid density, [kg m
3 ]
fi particle density distribution function along ith link s dimensionless relaxation time for density
direction s0 dimensionless relaxation time for temperature
eq
fi equilibrium distribution function of density along ith 4x lattice space
link direction 4t time step
F body force due to the presence of the porous medium, m fluid kinematic viscosity, [m2 s
1 ]
[N] / porosity
Fi total force term due to porous medium, [N] h dimensionless temperature, TT
T 1
w
T 1
g gravitational acceleration, [m s
2 ]
g eq equilibrium distribution function of temperature along b thermal expansion coefficient, [ C
1 ]
i
ith link direction st non-dimensional time,tmH1
gi temperature distribution function along ith link direc-
tion Subscripts
G body force due to gravity, [N] 1 far field value
H height of the cylinder, [m] o inlet value
K permeability of the material, [m2 ] M mean value
LD downstream length, [m] i lattice link direction
LU upstream length, [m] w wall
Ma Mach number, cms f front face
N number of lattices on the cylinder r rear face
@h
Nu local Nusselt number, @n s1 left face
Pr Prandtl number, a m
s2 right face
Re Reynolds number, m1mH

Fig. 1. An example of a scenario of flow around and through square shaped porous body - Fuel and control rods of a nuclear reactor mimic this situation.
856 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870

Porous bluff body research is also in vogue amongst researchers cylinder to channel ratio 0.1, 0.3 and 0.5 at Re = 250. Their results
due to primarily two main reasons: (a) Its simplicity in terms of depict that although at higher channel to cylinder ratios (B/
studying flow and heat transfer from it; (b) Representative model- H = 0.3 & B/H = 0.5) the buoyancy effect is not clear, it still has rich
ing of staggered real-life complex structures, which are numeri- heat transfer characteristics. However, at B/H = 0.1 the laminar
cally expensive to model. These porous bodies, of different flow motion is only noticeable near the porous square cylinder
shapes and at permeability levels, can seriously affect the flow and there is not much effect on cylinder wall. Guerroudj and
phenomena i:e. delay in occurrence and later, suppression of wake Kahalerras [30] numerically examined the influence of porous
and vortex shedding. And they can also provide significant aug- block shape and height on a partially heated lower plate under dif-
mentation in thermal dissipation [14–19]. At very low permeabil- ferent cross-flow buoyancy conditions. Their results show that the
ity levels (i.e. Da  10
6 ), the porous bluff body mimics a solid intensity of cross-flow buoyancy and Reynolds number increases
body and the heat transfer performance is same as that of an the strength of vorticity irrespective to the size and shape of por-
impermeable cylinder. Increase in permeability levels, allows fluid ous body. Furthermore, it has been reported that the heat transfer
to flow through the permeable cylinder, resulting in increased heat rate is maximum with buoyancy parameter (i.e. Ri) and it is found
transfer rate. Most of the studies [18–23] have concentrated on to be high for triangular shaped cylinder even at low permeability
forced convective heat transfer from a permeable body rather than values. Later, Guerroudj and Kahalerras [31] studied the flow and
practically occurring mixed convection. Consequently, it becomes heat transfer in an inclined channel provided with heated porous
necessary to look into the conjugative effects of mixed convection block on its lower plate. They have found that maximum and min-
heat transfer phenomena and porous media. Hwang and Chao [24] imum heat dissipation occurs when channel angle approaches
 
investigated the effects of permeability on laminar mixed convec- þ90 (aiding buoyancy) and
90 (opposing buoyancy), respec-
tion in a fully developed square channel filled with porous med- tively. However, an optimum angle for maximum Nusselt number

ium. They revealed that the nonuniform wall temperature has been reported as 30 for low permeability values. It can be
distribution greatly influences the flow and heat transfer charac- understood from the above discussion that mixed convection heat
teristics. Bae et al. [25] studied the mixed convection heat transfer transfer from an isolated porous square cylinder has not been stud-
in a channel containing multiple porous blocks heated from bot- ied yet, although studies on forced convection exist. Moreover,
tom. They used Brinkman-Forchheimer-extended Darcy model there is no information available in the literature on the drag and
and two-equation energy model to characterize the flow and tem- Nusselt number of the porous square cylinder which is important
perature inside the porous blocks in which higher heat transfer to design heat exchangers, LED backlight module cooling system
performances is reported. With reference to the experiments, Chou under the mixed convective conditions.
et al. [26] investigated the fully developed non-Darcian mixed con-
Complex applications usually require appropriate and, in most
vection in horizontal packed sphere channels experimentally. Their
cases, more sophisticated mathematical and numerical modeling.
analysis has shown that the buoyancy effect significantly affects
Since the advent of Boltzmann method for solving transport equa-
the flow structure and heat transfer rate at low Peclet number.
tions, LBM has shown its effectiveness against conventional
Kurtbas and Celik [27] experimentally investigated the mixed con-
approaches, especially while dealing with complex flow domain,
vection flow through a horizontal rectangular channel with open-
flow through porous media and multiphase flow. In 1988, LBM
cell metal foams of different pore densities. It has been revealed
was introduced by McNamara and Zanetti to overcome the draw-
that the buoyancy driven secondary flow and higher pore density
backs of the lattice gas cellular automata [32]. Since then research-
supplements the augmentation of the mean Nusselt number. From
ers have been started smearing this mesoscopic technique in
the aforementioned studies, it can be seen that higher heat transfer
various domains. LBM can be employed in wide range of applica-
rates could be attained by using porous media under the influence
tions in the field of Micro Electro Mechanical Systems (MEMS)
of buoyancy. Mixed convection flow through a porous medium has
and Nano Electro Mechanical Systems (NEMS) [33]. Also, compared
been applicable to electronic component cooling system, light
to the convectional numerical methods LBM is more appropriate
emitting diode (LED) backlight module system, LED street lamp
for parallel processing and it solves multiphase flows easily. Recent
cooling equipment, heat exchangers, thermal insulations and dry-
research articles [34–41] indicate that the flow modeling for non-
ing process [28]. A thorough knowledge of the flow and heat trans-
Newtonian fluids and nano-fluids is quite simple which carries an
fer parameters like drag, lift, wake length and Nusselt number is
addition of external force term with the LB flow evolution equa-
crucial to design aforementioned applications.
tion. Numerical simulations of unsteady flow around a square
On the parallel lines of the above discussion, work has also been
cylinder located between two parallel plates using FVM and LBM
done on mixed convection from porous bluff bodies for different
were investigated by Breuer et al. [42]. The vortex shedding phe-
operating conditions. It is well-known that a porous cylinder at
nomenon behind a porous square cylinder is probed by Babu and
lower values of permeability, behaves like a solid one. Hence,
Narasimhan [43] at Re = 100 and 200. They have employed
mixed convection from such a porous cylinder would obviously
D2Q9i lattice model in which entire collision equation with force
exhibit heat transfer trends similar to that from a solid cylinder.
term is written in terms of pressure. The main advantage of this
Such an instance can be visualised from the numerical study per-
approach is Navier-Stokes equations can be recovered from the
formed by Sharma et al. [29], wherein they showed the influence
of aiding buoyancy on a 2D solid square cylinder in vertical uncon- LB equations within OðMa2 ) [43]. Also, a force term F i which
fined configuration. They have calculated the local and average accounts for drag effects due to the presence of porous medium
Nusselt number for Reynolds number and Richardson number is added explicitly with the collision equation to characterize the
ranges from 1–40 and 0–1, respectively. It was found that the drag flow through the porous medium. Most of the other studies
coefficient increases at higher values of Ri, compared to the forced [44–49] imply that employing LBM for porous medium is quite
convection case. It was further emphasised that the average Nus- simple and efficient. By considering these advantages, we have
selt number is weakly dependent on the cross-buoyancy condition applied LB technique in our study. The objective of the present
in steady regime for an unconfined flow domain. They have also study is, therefore to investigate in detail the flow and thermal
formulated simple correlations for the wake length and average field around an isolated porous square cylinder and to explore
Nusselt number for the range of conditions embraced. Wu and the influence of permeability and aiding buoyancy on heat transfer
Wang [28] numerically investigated mixed convection effect on a from it. A LB model for the flow through porous media proposed by
porous square cylinder with the heated cylinder bottom for Guo and Zhao [45] is included.
 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870 857

2. Mathematical formulation cylinder when it is placed in an infinite stream, for making compu-
tationally feasible, the domain was limited to rectangular configu-
2.1. Problem definition ration. The following assumptions were made in order to make the
problem amenable to numerical simulations:
Consider a square cylinder of height H, maintained at a constant
temperature T w , and exposed to a constant free stream velocity v 1  The problem under consideration is steady, laminar, two-
(flowing in the positive y-direction) which is at an ambient tem- dimensional, incompressible and the fluid is Newtonian. The
perature T 1 as shown in Fig. 2. Although the objective is to simu- fluid flows around and through an isotropic, homogenous por-
late the flow and mixed convective heat transfer around the ous matrix with constant porosity and permeability.

Fig. 2. Computational domain used in the analysis of mixed convection flow and heat transfer from a porous square cylinder.
858 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870

 !
 All the fluid properties are considered to be constant and the @v 1 @v @v @p K @ 2 v @ 2 v
body forces are negligible, except for the body force terms in þ u þv ¼
/ þ þ
@ st / @x @y @y Re @x2 @y2
the collision equation (Boussinesq approximation in y direction)
1 1:75 1
in this study.
C1/ v
C 2 pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi
 There is no heat generation in the porous cylinder and viscous ReDa 150 Da
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
dissipation is negligible. Further, the fluid phase temperature u2 þ v 2
is equal to that of the porous phase (i.e. Local Thermal Equilib-  pffiffiffiffi v þ /Rih: ð7Þ
/
rium - LTE). This boundary condition should not be considered
for rapid cooling or heating process or problems with signifi-
 !
@h 1 @h @h Rk @2h @2h
cant temperature variation across the porous media [50]. þ u þv ¼ þ : ð8Þ
However, the temperature difference between the fluid and
@ st / @x @y Re  Pr @x2 @y2
porous phase in the present study is not significant and so
Here, C 1 and C 2 are the binary constants which hold the value of 0 at
the assumption is valid. The radiation from the permeable
clear fluid zone and 1 at the porous zone. The dimensionless vari-
cylinder is negligible.
ables are defined as: Pr ¼ lC p =k is the Prandtl number, K ¼ le =l,

2.2. Conservation equations denotes the viscosity ratio, Ri ¼ Gr=Re2 is the Richardson number.
The viscosity ratio and thermal conductivity ratio are assumed to
Considering above assumptions, the following governing equa- be unity in the present study.
tions are used to represent flow and heat transfer phenomena [51].
Continuity equation: 2.3. Lattice Boltzmann Method [LBM]

@u @ v  LBM is a kinetic based mesoscopic approach that can be viewed

þ ¼ 0: ð1Þ
@x @y as a special discrete scheme for the Boltzmann equation with dis-
crete velocities. In LBM, a single particle distribution function is
Momentum equations:
used to model the fluid. It exhibits many advantages over the con-

 !
1 @u u @u v  @u @p l @ 2 u @ 2 u ventional computational fluid dynamics. In LB equation the con-
q þ 2 þ 2  ¼
þ e þ vective term (i.e. streaming equation) is linear. Whereas, it is
/ @t / @x / @y @x / @x2 @y2
nonlinear in case of Navier stokes equations. Complex boundary
l qF s ~  conditions can be easily implemented by using particle reflection

u
pffiffiffiffi jUju : ð2Þ
K K rule. The conventional numerical methods usually needs more
! computations to solve the Poisson equation for the pressure field


1 @ v  u @ v  v  @ v  @p l @ 2 v  @ 2 v  l of the incompressible flows while in LB method the fluid pressure
q þ 2 þ ¼
þ e þ
v can be easily formulated with an equation of state [32].
/ @t / @x /2 @y @y / @x2 @y2 K
qF s ~ 

pffiffiffiffi jUjv þ qgbðT
T 1 Þ: ð3Þ 2.3.1. LB equation for velocity field
K In the traditional lattice Boltzmann method, the particles at
where le is the effective dynamic viscosity, p is the dimensional each node undergo collision followed by streaming. The extended
ffiffiffiffiffiffi 3=2
pressure and F s ¼ p1:75 1
is the inertial factor defined by Fu et al. LB collision equation for the porous medium [43] can be expressed
150 /
as
[52].
Energy equation: 1 eq 
f i ðx þ ei Dt; t þ DtÞ
f i ðx; tÞ ¼
f i ðx; tÞ
f i ðx; tÞ þ DtF i

 ! s
@T qC p @T v  @T @2T @2T þ DtF b : ð9Þ
qC p þ u  þ ¼ k e þ : ð4Þ
@t / @x @y @x2 @y2
eq
In Eq. (9), f i is the instantaneous particle density at a link, f i is the
The above governing equations can be written in the non- corresponding equilibrium density and ei are the velocity direction
dimensional form by introducing the following characteristic vectors of the particles residing in the lattice. In this study, D2Q9
scales: lattice arrangement which has eight particles around the central
particle has been used. The speeds of the particles can be written as
x y t1 p u 8
x¼ ; y¼ ; st ¼ ; p¼ ; u¼ ;
H H H qm21 m1 < ð0; 0Þ;
> i¼0
v T
T1 ei ¼ ðcos½ði
1Þp=2 ; sin½ði
1Þp=2 Þe; i¼1
4 ð10Þ
v¼ ; h¼ : >
: pffiffiffiffiffiffi
v1 TW
T1 ðcos½ð2i
9Þp=4 ; sin½ð2i
9Þp=4 Þ 2e; i ¼ 5
8:

Thus, the non-dimensional governing equations of macroscopic Also, the right side term of this collision equation is the BGK col-
scale can be written as follows: lision operator or Single Relaxation Time (SRT). In SRT model, the
population of particles resides in a single lattice are relaxed
@u @ v
þ ¼ 0: ð5Þ towards equilibrium with the same relaxation time. The relaxation
@x @y factor ‘s’ relaxes the particles towards to the equilibrium value

 ! during the collision of particles. The Chapman-Enskog relation
@u 1 @u @u @p K @ 2 u @ 2 u
þv
relates the LB fluid viscosity and the non-dimensional relaxation
þ u ¼
/ þ þ
@ st / @x @y @x Re @x2 @y2 time and it can be written as


1 1:75 1 1

C1/ u
C 2 pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi m¼ s
Dtc2s ð11Þ
ReDa 150 Da 2
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
u2 þ v 2 A force term F i which accounts the drag effects due to the presence
 pffiffiffiffi u: ð6Þ
/ of porous medium is directly coupled with the BGK collision
 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870 859

operator to recover fluid flow through the porous region. The force distribution functions f ðx; tÞ are then used in the temperature dis-
term that includes the effect of porosity (/) of the porous medium tribution function gðx; tÞ. The LB equation that describes the tem-
[43,45] can be written as perature filed can be written as

 
1 9 3 1 
F i ¼ wi q 1
3ðei :FÞ þ ðei :UÞðei :FÞ
ðU:FÞ : ð12Þ g i ðx þ ei Dt; t þ DtÞ
g i ðx; tÞ ¼
g i ðx; tÞ
g eq
i ðx; tÞ : ð20Þ
2s / / s0
Here, F is the body force due to the presence of porous medium Here, s0 is the dimensionless relaxation time and it is evaluated
which is otherwise termed as Darcy-Forchheimer force term [43]. from thermal diffusion coefficient (a ¼ m=Pr). The Chapman-
This can be written as Enskog relationship for thermal relaxation time can be written as


/m /cF 1
F¼
U
pffiffiffiffi jUjU þ /G: ð13Þ a ¼ s0
Dtc2s : ð21Þ
K K 2

where m is the kinematic viscosity of the fluid, K is the permeability The temperature equilibrium distribution function g eq
i in the
of the porous medium, cF is the non-dimensional Forchheimer Eq. (20) can be written as
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
form-drag coefficient or inertial factor, jUj ¼ u2 þ v 2 in which u g eq
i ¼ wi T½1 þ 3ðei :UÞ : ð22Þ
and v are the components of U in the x- and y-directions, and G is
the body force due to gravity. The equilibrium distribution function In the above equation, wi is the same weighing factors which has
which appears on the right side of Eq. (9) should consider the poros- been defined in the flow equilibrium distribution function and T
ity of the porous medium. It can be written as follows: is the temperature of the fluid. The temperature T can be calculated
  from the distribution function gðx; tÞ as
eq 9 3
f i ¼ wi q 1 þ 3ðei  UÞ þ ðei  UÞ2
ðU  UÞ : ð14Þ X
8
2/ 2/
T¼ gi : ð23Þ
The weighing factors wi are an important part in the calculation of i¼0

equilibrium distribution functions. Each particle residing in the lat-

tice will have its own weighing factor. In case of D2Q9, w0 ¼ 4=9, 2.3.3. Boundary conditions
wi ¼ 1=9 for i = 1–4 and wi ¼ 1=36 for i = 5–8. After performing In practical cases, the boundary conditions are applied in terms
the simulation in the mesoscopic level, the macroscopic properties of macroscopic variables like q and U. This Lattice Boltzmann sim-
needs to be calculated. The addition of distribution function of each ulation has the information of particle distribution functions
lattice represents density of the fluid. This relationship between the before calculating the macroscopic values. Therefore, boundary
distribution function and density is very helpful in finding out the conditions also need to be specified in terms of distribution func-
macroscopic properties such as density and velocity of the fluid. tions only. Various boundary conditions used in this study are
Density q and velocity u can be evaluated from the distribution described below.
function as follows: Inlet boundary condition: The flow simulated by LB method is
X
8 weekly compressible and hence, the inlet velocity has to be chosen
q¼ f i: ð15Þ as Mach number (Ma), less than 0.1 [32]. A free uniform velocity
i¼0 (v 1 ), which has less compressibility effects is given to the inlet
The auxiliary velocity V can be defined as and it can be indicated as v 0 . The distribution functions f 2 ; f 5
and f 6 fall inside the domain at the inlet region (bottom of the
X8
Dt computational domain). The Zou and He boundary condition [53]
qV ¼ ei f i þ qF: ð16Þ
is applied to inlet of the domain to find out the unknown distribu-
i¼1
2
tion functions and it can be written as follows:
The actual velocity components of porous medium can be calcu-
1
lated from the following equation: q0 ¼ ½f þ f 1 þ f 3 þ 2ðf 4 þ f 7 þ f 8 Þ : ð24Þ
1
v0 0
V
U¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : ð17Þ
2 1 1
c0 þ c20 þ c1 jVj f 2 ¼ f 4 þ q0 v 0 ; f 5 ¼ f 7 þ q0 v 0 ; f 6 ¼ f 8 þ q0 v 0 : ð25Þ
3 6 6
The two parameters c0 and c1 in Eq. (17) are given by In the above equation q0 and v 0 are the initial density and velocity

 (in y-direction) respectively.
1 Dt m Dt c F
c0 ¼ 1þ/ ; c1 ¼ / pffiffiffiffi : ð18Þ Left and Right boundary: The study aims to simulate flow around
2 2 K 2 K
and through the square porous cylinder in the unconfined domain.
To analyse forced convection, a separate distribution function is Artificial left and right boundaries are placed far away from the
used and to show the effects of natural convection, an extra force cylinder in order to avoid the shear stress developed by the walls.
term (F b ) is added to the LB collision equation, because the flow Accordingly, the slip boundary condition needs to be adopted to
in this case is driven by temperature gradient. The Boussinesq satisfy this requirement. In case of left boundary, the unknown dis-
approximation in LB can be written as follows: tribution functions are f 1 ; f 5 and f 8 . As per the slip condition [54],
for the entire length of left boundary the following conditions are
F b ¼ 3wi gbheiy : ð19Þ
used. f ð1;xÞ ¼ f ð1;x
1Þ ; f ð5;xÞ ¼ f ð5;x
1Þ and f ð8;xÞ ¼ f ð8;x
1Þ .
Here, eiy is the velocity vectors acting in the y-direction. Outlet boundary condition: General practice of evaluating the
unknown distribution functions at the outlet is through extrapola-
2.3.2. LB equation for temperature field tion [32]. At the outlet of the domain the unknown distributions
To evaluate the temperature physiognomies, one more distribu- are f 4 ; f 7 and f 8 which can be evaluated by using second order
tion function called gðx; tÞ is used. In both forced and mixed con- polynomial.
vection, the energy equation can be solved after flow field is The porous square region is kept at the non-dimensional tem-
calculated. The macroscopic values obtained from the flow field perature value of 1 (i.e. h = 1).
860 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870

2.4. Dimensionless parameters 3. Numerical procedure and validation

Some of the dimensionless parameters used in this study are 3.1. Numerical method
defined as follows:
The dimensionless parameters Reynolds number (Re) and Darcy A FORTRAN code has been written based on lattice Boltzmann
number (Da) are written in terms of lattice numbers. method to perform the present numerical simulations. The colli-
sion equation with the external force term (F i ) (Eq. (9)) is applied
v 1H v1N
Re ¼ ¼ : ð26Þ to the entire computational domain. A constant value of 1 is mul-
m m tiplied with this force term in the porous square region. In all other
where the height of the cylinder ‘H’ is the characteristic length and places the porous force term (F i ) is nullified to get the clear fluid
the free stream velocity ‘v 1 ’ is the characteristic velocity. The region flow. However, the buoyancy force term (F b ) is applied in
porosity (/) and Darcy number (Da) could be related through the the entire domain. In this study, BGK collision operator with single
Carman-Kozeny relation [4] given by, non-dimensional relaxation factor (s) has been used. The initializa-
tion process carries the formulation of initial velocity (v 0 ) and LB
2
1 /3 dp K K kinematic viscosity based on the number of lattices on characteris-
Da ¼ ¼ ¼ : ð27Þ tic length (N). Through Chapman-Enskog relationship [56] the
180 H ð1
/Þ2 H2 N2
2
relaxation factors (s&s0 ) are calculated. Initially, the velocities
In LBM, the numerical stability and accuracy of the results are (i.e. u and v) at all nodes are assumed to be zero and the density
mainly dependent on the number of lattices on characteristic length value q is kept 1. The Boussinesq force term which appears in
and in present analysis; N represents the number of lattices that are Eq. (9) has gb and this can be calculated from Richardson number.
present on the front face of the cylinder. In this study we have taken It should be noted that the velocity scale is proportional to
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
the particle diameter (dp ) as 1 cm (for H = 1 m) [55]. gb 4 TH. In order to ensure the problem in incompressible
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
The dimensionless stream function w is defined as regime, the value of gb 4 TH should be less than 0.1 [32]. It is
@w @w seen that the numerical algorithm of lattice Boltzmann is relatively
u¼ and v ¼
: ð28Þ simpler than the conventional CFD methods. In the present work,
@y @x
the LBM computations are carried out until the following conver-
The curl of the flow velocity or vorticity (x) for two-dimensional gence is satisfied.
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
i2ffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
i2ffi
flows can be evaluated as P h kþ1 k P h kþ1 k
i;j U i;j
U i;j i;j T i;j
T i;j
@ m @u rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 6 1  10
6 ; rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 1  10
6 :
x¼

@x @y
: ð29Þ P h kþ1 i2 P h kþ1 i2
i;j U i;j i;j T i;j
Heat transfer from a permeable body to the surrounding fluid is
@h
ð33Þ
evaluated from the local Nusselt number is given by Nu ¼
Rk @n .
where U i;j is the fluid velocity, T i;j is the fluid temperature and k is
Here, n is the direction normal to the body surface and Rk is the
the iteration level.
ratio of thermal conductivity of the porous material and fluid,
which is assumed to be 1 in this study. Average heat transfer at
any face of the permeable body is obtained by integrating the local 3.2. Mach number annealing
Nusselt number along that face. The surface averaged Nusselt
number on the front (f), rear (r), left (s1) and right (s2) face of In lattice Boltzmann, the compressibility errors can be reduced
the cylinder (see Fig. 2) is calculated as significantly by reducing Mach number (Ma). This unfortunately
blows up the computational time needed for convergence [57].
Z B Z D
Hence, in case of low Re flows, for the same values of characteristic
Nuf ¼
Nu  dx; Nur ¼
Nu  dx; length (H), the initial velocity (v 0 ) has to be reduced. It implies that
A C
Z A Z C
ð30Þ the Ma(=v/cs ) is actually getting reduced. To accelerate the simula-
Nus1 ¼
Nu  dy and Nus2 ¼
Nu  dy: tions even at low velocities, Artoli et al. [58] have proposed an idea
D B
which is based on stepwise reduction of Ma after the simulation
The mean Nusselt number value (NuM ), which represents the converges at higher value of Mach number. In order to perform
total heat transfer from the permeable cylinder is evaluated from simulation at low Ma, the velocity v 0 and consequently the viscos-
the average Nusselt number of all surfaces as shown below ity (m) have to be reduced to produce the same Re. Since Ma is
reduced in steps, this technique is called Mach number annealing.
1X
NuM ¼ Nuf ;r;s1;s2 : ð31Þ The annealing factor n can be define as
4 AB
Ma v m
n¼ ¼ ¼ : ð34Þ
The heat transfer enhancement obtained by permeable cylinder Ma0 v 0 m0
with respect to solid cylinder (or permeable cylinder with Da = This equation implies that all the velocity (in terms of driving force)
10
6 ) is expressed in terms of Heat Transfer Enhancement Ratio and the viscosity values are to be reduced n times. This annealing
(E). It is given by the ratio of Nusselt number of any permeable can be of single (direct) and multilevel. In one level annealing strat-
cylinder to Nusselt number of permeable cylinder with Da = 10
6 . egy, after the simulation converges with a higher Mach number, the
viscosity and driving force are reduced n times in a single step and
ðNuÞat any value of Da the simulation converges to the final solution. The gradual reduc-
E¼ : ð32Þ tion of these parameters towards n is followed in case of multilevel
ðNuÞat Da¼10
6
annealing. This level of annealing is mainly dependent on the
The present study elaborates on the variation of E at different faces stability and tolerance constraints. In this study, we have used
of the cylinder and, hence, heat transfer enhancement ratio at dif- direct or single step annealing to accelerate the simulation. The
ferent faces is represented as Ef ; Er ; Es1 & Es2 for front, rear, left non-dimensional relaxation values for hydrodynamic (s) and ther-
and right faces, respectively. mal (s0 ) are taken as 0.56 and 0.5845, respectively.
 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870 861

3.3. Grid sensitivity analysis 11

Present study (Da = 10-6 Ri =0)
Dhinakarn and Ponmozhi [51]
10
A lattice sensitivity study is carried out at Re = 40, Da = 10
6 Sharma and Eswaran [59] Solid & Ri = 0
Paliwal et al. [60]
and Ri = 1 to check the dependence of flow and heat transfer -2
Present study (Da = 10 & Ri =0)
9 Dhinakarn and Ponmozhi [51], (Da = 10-2 & Ri =0)
parameters such as C D and NuM on the number of lattice units -6
Present study (Da = 10 , Ri =1)
Sharma et al.(Da = 10-6 & Ri =1) [29]
and downstream length (LD ). Based on number of lattice units 8
on characteristic length, the relaxation parameter (s) that
appears in the collision equation changes and so the stability 7
and accuracy of the solution. Although more number of lattices

CD
could deliver accurate results, the convergence rate of numerical 6
solution slows down due to the decrement of initial velocity
5
value. A systematic test simulations have been performed with
20–40 lattice units on front face of the cylinder and is presented 4
in Table 1. The simulation running time at LD ¼ 30 is noticed to
be 37.35, 44.25, 51.45, 60.3, 72.4 and 84.2 h for N = 20, 24, 28, 3
32, 36 and 40, respectively. Also, the results obtained at this
downstream length has shown less variations with LD ¼ 40. 2
Hence, the study suggest that, the grid size of 1312  960 (32
(a) 1
lattices on characteristic length) with downstream length of 5 10 15 20 25 30 35 40
30D will produce results independent of lattice units with less Re
computational time.
3.5
3.4. Code validation
3.2
The numerical results obtained from present LBM code with air
as a working fluid (Pr = 0.71) is compared with the studies avail- 2.9

able in the literature at 5 6 Re 6 40, Da = 10
3 and 10
2 and

Ri = 0 and 1. To consider the variation of both permeability and 2.6

buoyancy conditions, different aforementioned parameters have

been chosen to compare with the available literatures. In Fig. 3 2.3
NuM

(a), the obtained drag coefficient from present code is compared

with the results of Dhinakaran and Ponmozhi [51], Sharma and 2

Eswaran [59], Paliwal et al. [60] and Sharma et al. [29]. Fig. 3(b)
1.7
graphically showing the comparison of mean Nusselt number
-6
Present study (Da = 10 , Ri = 0)
obtained from the present code with above literatures at different Dhinakaran and Ponmozhi [51]
1.4 Solid & Ri =0
permeability levels and buoyancy condition. The graphical illustra- Sharma and Eswaran [59]
-3
Present study (Da = 10 & Ri =0)
tion clearly represents that the present numerical results for C D Dhinakaran and Ponmozhi [51], (Da = 10-3& Ri =0)
1.1 -6
Present study (Da = 10 , Ri =1)
and NuM are in good agreement with the results available in the lit- Sharma et al.(Da = 10-6 & Ri =1) [29]
erature. A maximum error of 1.4% and 3.4% is obtained in the case
(b) 0.8
of drag coefficient and mean Nusselt number, respectively, and this 5 10 15 20 25 30 35 40
deviation is due to the numerical approach and discretization Re
schemes. It is noteworthy that the literatures considered for
numerical result comparison have used finite volume approach; Fig. 3. Comparison of present study results with literature for Pr = 0.71: (a) Drag
whereas, the present results are simulated from lattice Boltzmann coefficient (C D ) and (b) Mean Nusselt number (NuM ) at different values of Da and Ri.

method.

4. Numerical results and discussions  Reynolds number, Re = 2–40.

 Richardson number, Ri = 0, 0.5 and 1.
The 2-D numerical simulations are carried out for flow across a  Darcy number, Da = 10
6 ; 10
4 ; 10
3 and 10
2 .
permeable square cylinder in an unconfined vertical flow regime,  Porosity, / = 0.629, 0.935, 0.977, 0.993.
under the influence of aiding buoyancy for the following range of  Prandtl number, Pr = 0.71(air).
physical parameters:

Table 1
Analysis on the effects of number of lattice units (N) and downstream length on drag coefficient (C D ) and mean Nusselt number (NuM ) at Re ¼ 40; Da ¼ 10
6 ; Pr ¼ 0:71 and
Ri ¼ 1.

N Grid (m  n) Drag coefficient (C D ) Mean Nusselt number (NuM )

LD = 20 LD = 30 LD = 40 LD = 20 $L_D = 30 LD = 40
20 820  600 4.175 4.054 4.055 3.205 3.241 3.241
24 989  720 3.612 3.587 3.520 3.225 3.289 3.289
28 1148  840 3.701 3.562 3.562 3.255 3.302 3.302
32a 1312  960 3.570 3.534 3.534 3.266 3.305 3.306
36 1476  1080 3.570 3.534 3.534 3.266 3.305 3.306
40 1640  1200 3.570 3.534 3.534 3.266 3.305 3.306
a
Number of lattice units placed on front face of the cylinder in this study.
862 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870

Da = 10−6 Da = 10−4 Da = 10−3 Da = 10−2

Ri = 0
2
Re = 20

0
X

-1

-2
2
Re = 40

0
X

-1

-2 -2 0 2 4 6 -2 0 2 4 6 -2 0 2 4 6
-2 0 2 Y 4 6 Y Y Y

Ri = 0.5
2
Re = 20

0
X

-1

-2
2
Re = 40

0
X

-1

-2 -2 0 2 4 6 -2 0 2 4 6 -2 0 2 4 6
-2 0 2 Y 4 6 Y Y Y

Ri = 1
2
Re = 20

0
X

-1

-2
2
Re = 40

0
X

-1

-2 -2 0 2 4 6 -2 0 2 4 6 -2 0 2 4 6
-2 0 2 Y 4 6 Y Y Y

Fig. 4. Streamlines at Re ¼ 20 and 40; Ri ¼ 0; 0:5 and 1 at various values of Darcy number for the flow around and through the porous square cylinder.

This numerical study will serve us to gauge the effects of buoy- that at very low Re (62), the flow is completely attached to the
ancy and permeability on flow and thermal patterns. The results cylinder due to excess viscous forces on its surface. Upon increas-
are labelled into flow and heat transfer characteristics and are pre- ing fluid velocity, the inertial forces become dominant and hence,
sented as follows: flow separation takes places at the trailing edge of the cylinder.
This phenomenon is common in the case of flow over a solid square
4.1. Flow characteristics obstacle and the same pattern can be observed at Da = 10
6 with
no buoyancy factor (i.e. Ri = 0). Furthermore, the wake length
A detailed discussion on variation of flow behaviour with Rey- decreases with an increase of Ri for all values of Da (except
nolds number (Re), Darcy number (Da) and Richardson number Da = 10
2 ) at constant Re. At higher buoyancy force, the velocity
(Ri) is presented in this section. gradient is more and it results in reduction of pressure over the
cylinder. Owing to this, reduction of wake zone is noticed with
4.1.1. Streamline pattern an increase in Ri. The increment of permeability level gradually
In Fig. 4 we present the representative flow patterns in the allows fluid into the cylinder with different deviation levels. For
vicinity of permeable square cylinder through streamline contours example, the deviation of flow through the cylinder at Da = 10
4
at different values of Reynolds number (Re = 20 and 40), Richard- is greater due to the excess viscous drag term when compared with
son number (Ri = 0, 0.5 and 1) and Darcy number Da = 10
3 . A closer view of streamline pattern through the perme-

6
4
3
2
(Da = 10 ; 10 ; 10 and 10 ). It is clear that, streamlines are able cylinder is presented in Fig. 5, which shows the variation of
a function of Reynolds number and are also altered with the intro- fluid deviation inside the porous zone with permeability and buoy-
duction of different permeability and buoyancy levels. It is seen ancy. The flow deviation level inside the cylinder also varies due to
 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870 863

Da = 10−4 Da = 10−2

(a)

(b)

(c)
Fig. 5. A closer view of streamline pattern through the permeable cylinder at Re = 40 for two different Da values and (a) Ri = 0, (b) Ri = 0.5 and (c) Ri = 1.

buoyancy effects. For a fixed Re and Da, upon increasing Ri the flow cylinder. Therefore, this streamline flow pattern concludes that
deviation starts reducing and this is due to the dominance of buoy- the aiding buoyancy and permeability prove to be influential on
ancy over viscous drag term of porous region. To emphasis this, we drag and inertial forces inside the cylinder.
present the average normal velocity (U av g ) variation over leeward,
middle and right/left surface of the cylinder in Fig. 6. The variation 4.1.2. Vorticity pattern
of average normal velocity over different surfaces clearly exhibits To further provide an insight on the behaviour of fluid flow,
the influence of permeability and buoyancy. For a constant value especially near the porous square zone, the vorticity patterns are
of Da and Ri; U av g is higher at the left/right surface followed by also presented in Fig. 7. The solid and dotted lines in this contour
middle and leeward surface of the cylinder. This clearly indicates frames indicate positive and negative values of vorticity. The flow
that, due to resistance created by the porous medium, the fluid phenomenon for different Re has been well elucidated earlier in
tends to deviate; consequently, average normal velocity is higher Section 4.1.1. In this section, effects of Ri and Da on vorticity
at these surfaces. The velocity of fluid on different surfaces strength are discussed, and same is presented in Fig. 7 for
increases with an increase in permeability level and this increment Re = 40. The intensity of local rotation increases with Re for a con-
is drastic when the permeability level is shifted from 10
3 to 10
2 . stant value of Ri, due to the increase of inertial forces and it
Hence, more streamlines pass through the cylinder at this level of decreases with the increasing Ri for a fixed value of Re. The similar
Da. It can also be seen from Fig. 6 that the average normal velocity behaviour of vorticity pattern with respect to Re and Ri is observed
increases with increasing buoyancy levels for same value of Da. It is for a long square bar under the influence of buoyancy by Sharma
evident that increasing buoyancy can stretch the fluid inside the et al. [29]. Considering the effect of Da, the vorticity strength
porous medium with relatively higher velocity. However, for all changes due to the variation of fluid penetration. At Ri ¼ 0, the
Re and buoyancy levels the flow is observed to be more elongated changes in permeability at constant Re shows the similar patterns
at Da = 10
2 and there is no wake formation downstream of the as derived by Dhinakaran and Ponmozhi [51]. At Da = 10
6 , the
864 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870

0.7 4.1.3. Drag coefficient ðC D Þ

To convey the quantitative effects of Da and Ri at various Re on
Ri = 0
0.6 Ri = 0.5 Leeward surface
flow behaviour around and through the permeable cylinder, the
Ri = 1 coefficient of drag (C D ) variations have been presented in Fig. 8.
Ri = 0
Ri = 0.5 Middle surface The drag force experienced by the cylinder is calculated using
0.5 Ri = 1 momentum exchange method [61], in which the exchange of
Ri = 0
Ri = 0.5 Left/Right surface momentum between two opposing directions of the neighbouring
0.4
Ri = 1 lattices are used to evaluate the drag and lift forces. It is a well
known fact that, the drag value decreases with the increment in
Uavg

Re due to the depletion of viscous layer around the cylinder. The

0.3 Left
low fluid resistance inside the porous zone at higher permeability
Leeward
Middle

reduces the pressure gradient across the cylinder. This occurrence

0.2 results in decrement of drag coefficients for all the ranges of buoy-
Right ancy conditions considered. In short, the drag coefficient decreases
0.1
with the increase of Da irrespective of Re and/or Ri. At Re P 25,
drag values of the cylinder at Da ¼ 10
3 approaches lower perme-
ability drag values for all ranges of Ri. At Da ¼ 10
2 , fluid pene-
0
trates with very less deviation inside the cylinder; the cylinder
10
-6
10
-5
10
-4
10
-3
10
-2
experiences lower drag. The non-linear drag trend, with variation
Da of Re is found to be consistent for all values of Ri. The viscous dom-
inance due to aiding buoyancy increases the drag value of the
Fig. 6. Average normal velocity (U Av g ) along the leeward, middle and right/left cylinder for constant Re and Da. Thus, it can be concluded that,
surface of the porous square cylinder for different Da and Ri values at Re = 40.
the drag coefficient increases with an increase in Ri substantially
and the same decreases monotonously while increasing the per-
meability level.
vorticity patterns are symmetric along a vertical line and it is com-
pletely attached to the cylinder surface. Due to the reduction of vis- 4.1.4. Recirculation length (Lr )
cous resistance and increase in inertial forces, the fluid flow inside From the streamline patterns it can be observed that due to
the cylinder increases drastically. Therefore, the fluid diffusion dominant inertial forces with the increase in Re, two counter rotat-
inside the cylinder increases with an increase in permeability. It ing eddies are formed in downstream of the cylinder. The stream-
can be concluded from vorticity frames that both buoyancy effects wise distance between top surface and wake bubble stagnation
and permeability of the cylinder reduces the strength of vorticity. point is termed as wake/recirculation length. Fig. 9 neatly presents
Buoyancy force degrades the pressure gradient outside the cylinder the variation of wake length with Reynolds number at
whereas, more permeable cylinder reduces the viscous drag inside Da = 10
6 ; 10
4 and 10
3 and at Ri = 0, 0.5 and 1. The symmetric
porous zone. vortex pair behind the cylinder grows with the increase of Re

Da = 10−6 Da = 10−4 Da = 10−3 Da = 10−2

Ri = 0
-0.2
1 -0.6 -0.4 -0.8 -0.4 -0.8 -0.4
-0.8
-0.6 -0.6 -0.6
-0.4

0
X

0.7 0.7 0.7

0.9 0.5 0.9 0.5
-1 0.8 0.6 0.5
0.4 0.3

-0.5 0.5 1.5 Y

2.5 3.5 4.5 -0.5 0.5 1.5 Y
2.5 3.5 4.5 -0.5 0.5 1.5 Y
2.5 3.5 4.5 -0.5 0.5 1.5 Y
2.5 3.5 4.5

Ri = 0.5
1 0.1
-0.3 -0.1 -0.1
-0. -0.3 -0 . -0.1
-0

5 6 -0.4 4
-0 . -0.7
.7

-0.7
-0.2
0
X

0.2
0.7
1

0.3
5

0.4 0.5
0.

0.
6 0.6 0.4 0.2
0.

0.2
-1

-0.5 0.5 1.5 Y

2.5 3.5 4.5 -0.5 0.5 1.5 Y
2.5 3.5 4.5 -0.5 0.5 1.5 Y
2.5 3.5 4.5 -0.5 0.5 1.5 Y
2.5 3.5 4.5

Ri = 1
0.1
1 0.1 0.1 0.1
0.3
0.2 0.5
-0.3 -0.1 -0 . -0 -0. 0.3 -0
3 .1 3 .3
-0.1

0
X

2 -0.2
1 0.
6

0.5 0.2 -0 . 0.4 0.2 1 0.4 -0.3 3

0.

-0 . 0.1 -0 . -0.5
-0.1
-1 -0.1

-0.5 0.5 1.5 Y

2.5 3.5 4.5 -0.5 0.5 1.5 Y
2.5 3.5 4.5 -0.5 0.5 1.5 Y
2.5 3.5 4.5 -0.5 0.5 1.5 Y
2.5 3.5 4.5

Fig. 7. Vorticity contours at Ri ¼ 0; 0:5 and 1; Re = 40 at various values of Darcy number for the flow around and through the porous square cylinder.
 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870 865

and it decreases slightly with the introduction of higher thermal boundary layer in the region due to fluid recirculation. The
permeability. At Re P 35 the stagnation point of wake moves value of h for Da ¼ 10
6 along x-axis is almost same as that of
downstream when compared with the lower permeable cylinder Da ¼ 10
4 . The temperature profile expands further with the incre-
(i.e. Da ¼ 10
6 ) for the case of forced convection (i.e. Ri = 0). ment in Da and it becomes bell-like curve at Da ¼ 10
2 . This sug-
However, the aiding buoyancy cases (Ri = 0.5 and 1) have shown gests the presence of a thick thermal boundary layer at leeward
inverse trend with the increase in Da for all values of Re. From surface of the cylinder. From Fig. 11(b), one can see the variation
Fig. 9 it can be understood that buoyancy is more influential on of temperature profile with buoyancy. The lateral width of the
wake length variation than permeability. For fixed values of Re curve shrinks with the increment in Richardson number. Hence,
and Da, the wake stagnation point move towards the top surface it is clear that the thermal boundary layer around the cylinder
of the cylinder. The pressure reduction over the cylinder surfaces reduces while increasing Richardson number. At higher Ri, more
due to aiding buoyancy is the cause for the streamwise compres- fluid sticks over the cylinder which results in rich heat transfer
sion of wake bubble. The wake length reduction is found to be rate. It is to be noted that the shrinkage of temperature profile is
largest when the Richardson number increases from 0 to 0.5. The prominent while Ri shifts from 0 to 0.5 and it is less sensitive while
percentage decrement in wake length at Re = 40 with Da = 10
3 is it varies from 0.5 to 1.
494.24%; whereas, the corresponding difference is observed
88.8%. Thus, it can be concluded that buoyancy effect is prominent
on wake length reduction than the permeability levels of the 11
porous zone.
10 Da = 10-4
-3
Da = 10 Ri = 0
4.2. Heat transfer characteristics Da = 10-2
9 -4
Da = 10
Heat transfer characteristics are a function of Reynolds number Da = 10-3 Ri = 0.5
8 Da = 10
-2
(Re), non-dimensional permeability (Da) and type of buoyancy
Da = 10-4
condition on the porous cylinder surface. A detailed study on vari- 7
-3
Da = 10 Ri = 1
-2
ation of heat transfer performance with aforementioned parame- Da = 10
CD
ters is presented in this section. 6

4.2.1. Isotherms 5
The effect of physical parameters such as Re; Da and Ri on tem-
4
perature distributions near the porous square zone for the buoy-
ancy aided flow is presented in Fig. 10 for a Prandtl number 3
value of 0.71. The thermal pattern around the cylinder is greatly
affected by the Reynolds number, Darcy number and Richardson 2
number. As Re increases, the isotherm pattern spread along the
streamwise direction with a subsequent reduction of the lateral 1
5 10 15 20 25 30 35 40
width. This is because at higher Re, the convective transport
becomes dominant and thermal boundary layer thickness reduces. Re
It can also be seen from Fig. 10 that the size of the thermal plume Fig. 8. Variation of drag coefficient (C D ) with Re for different values of Ri and Da, for
increases with an increase in Da for a fixed value of Re and Ri. For the flow through and around a porous square cylinder.
example, at Re = 40, looking at different values of Da in Fig. 10, one
can observe the largest thermal plume at Da = 10
2 for all values of
Ri. As the higher permeable cylinder allows more fluids into it,
more heat is dissipated from the cylinder and hence, thermal field 2.8
is stretched at higher permeability than lower permeability case. Da = 10
-6

-4
Da = 10 Ri = 0
Similar patterns have also been reported by Dhinakaran and Pon- -3
2.4 Da = 10
mozhi [51]. Coming to the effect of buoyancy, the isotherms are Da = 10-6
Da = 10-4 Ri = 0.5
more clustered near the surfaces of porous square cylinder with Da = 10-3
the increase in Richardson number from 0 to 0.5. Also, the thermal 2 Da = 10
-6

-4
Da = 10 Ri = 1
plume stretch is large and contract for a constant Reynolds number
Wake length, Lr

-3
Da = 10
and Darcy number due to this increment. Therefore, the combined
1.6
effect of Da and aiding buoyancy has the largest and narrow ther-
mal field which implies higher heat transfer rate. Thus, the reduc-
tion of thermal boundary layer around the porous zone due to the 1.2
higher velocity gradients around the cylinder and dominant iner-
tial forces inside the porous zone results in rich heat transfer
performance. 0.8

4.2.2. Temperature profile 0.4

The influence of Darcy number and aiding buoyancy forces on
non-dimensional temperature (h) profile along the x-direction at
y = 11.5H for Re = 40 is shown in Fig. 11(a) and (b). In the vicinity 0
10 15 20 25 30 35 40
of the hot porous cylinder, the non-dimensional temperature value
Re
increases and it reaches its maximum value near the centre point,
followed by a minor dip at the centre for lower values of Da and at Fig. 9. Wake length (Lr ) variation at different values of Re; Da (10
6 - 10
3 ) and Ri
Ri = 0. This dip in the middle of domain indicates the reduction of (0, 0.5 and 1) for flow around and through porous square cylinder.
866 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870

Da = 10−6 Da = 10−4 Da = 10−3 Da = 10−2

Ri = 0
1.5

0.15 0.15
0.25
0.2 0.35
Re = 20

0.25 0.35 0.25 0.45

0.3 0.4 0.45
5
0.5

0.45 0.35
X
-0.5

0.4
0.25 0.4 0.3 0.4 0.4
0.3
0.2 0.2 0.3
0.15 0.2
1.5 -1.5

0.15
0.15 0.15 0.25
0.15 0.25 0.35
0.45
Re = 40

0.25 0.45 0.35 0.25 0.35

0.4 0.3 0.45
5 5
0.5
X
-0.5

0.4 0.3 0.4 0.3 0.4

0.2 0.2 0.3 0.4
0.1 0.2 0.3
0.1 0.1 0.2
-1.5

-1 0 1 2 Y 3 4 5 -1 0 1 2 Y 3 4 5 -1 0 1 2 Y 3 4 5 -1 0 1 2 Y 3 4 5

Ri = 0.5
1.5

0.15 0.15 0.15

0.15
Re = 20

0.4 0.25 0.25

5 0.35 0.4 0.35 0.25 0.4 0.25 0.45 0.35
0.5

5 5 0.35
X
-0.5

0.4 0.3 0.4 0.3 0.4

0.4 0.3 0.2 0.3
0.2 0.2 0.2
0.1 0.1 0.1 0.1
1.5 -1.5
Re = 40

0.15 0.15 0.15 0.25

0.5

0.4 0.35 0.25 0.4 0.35 0.25 0.4 0.25 0.35

5 5 5 0.35 0.45
X
-0.5

0.4 0.3 0.4

0.4 0.3 0.2 0.4 0.3 0.3
0.2 0.2 0.2
0.1 0.1 0.1 0.1
-1.5

-1 0 1 2 Y 3 4 5 -1 0 1 2 Y 3 4 5 -1 0 1 2 Y 3 4 5 -1 0 1 2 Y 3 4 5

Ri = 1
1.5

0.15
Re = 20

0.15 0.15 0.15

0.25
0.5

0.4 0.35 0.25 0.4 0.25 0.45 0.25 0.35

5 5 0.35 0.35 0.45
X
-0.5

0.4 0.3 0.4 0.3 0.4

0.2 0.2 0.4 0.3 0.2 0.3
0.2
0.1 0.1 0.1 0.1
1.5 -1.5
Re = 40

0.15 0.25
0.5

0.4 0.15 0.15

5 0.25 0.4 0.25 0.45 0.35 0.25
0.35 5 0.35 0.45
X

0.4 0.3 0.4 0.3

-0.5

0.4 0.3 0.2 0.2 0.2 0.3

0.1 0.1 0.1 0.1
-1.5

-1 0 1 2 Y 3 4 5 -1 0 1 2 Y 3 4 5 -1 0 1 2 Y 3 4 5 -1 0 1 2 Y 3 4 5

Fig. 10. Isotherms at Re ¼ 20 and 40; Ri ¼ 0; 0:5 and 1 at various Darcy numbers for the flow around and through the porous square cylinder (shaded contour levels between
1 6 h 6 0:5).

4.2.3. Heat transfer enhancement for different faces causes heat transfer reduction at all surfaces except front face of
Fig. 12(a–c) depict the representative variation of the heat the cylinder. The heat transfer enhancement ratio at left/right sur-
transfer enhancement ratio on each surface of the porous square face (Es1 or Es2 ) increases at lower values of Re with the increase of
cylinder for different values of Re; Ri and Da. Since it is well known Da for all values of Ri and later it reduces with the increase of Da.
that the heat transfer performance increases with Re, we have pre- Looking into the effect of aiding buoyancy with permeability on
sented the effects of aiding buoyancy and Da in terms of heat trans- enhancement ratio (E), the top surface is inversely proportional
fer enhancement ratio (E). The front face of the square cylinder is in and at other surfaces it is almost directly proportional to the Darcy
direct contact with fluid, and so produces thermal enrichment at numbers. However, it is seen that without the permeability effect,
this surface compared to other surfaces. It can be observed from the heat transfer performance at all surfaces increase with the
Fig. 12(a–c) that at Da ¼ 10
4 , the heat transfer enhancement ratio increase of Ri irrespective of Re.
is slightly more than 1 for all faces of the cylinder which implies
that, at this permeability level the heat transfer augmentation is 4.2.4. Mean Nusselt number (NuM )
not much significant. The chance given by front surface to enter The effect of non-dimensional permeability on mean Nusselt
more fluid at higher permeability results in rich heat transfer at number at various values of Re and Ri is illustrated in Fig. 13(a–
this surface. Upon increasing Da, the front surface heat transfer c). NuM increases almost linearly with the increase of Re at all val-
enhancement ratio (Ef ) increases and it is more drastic when the ues of Ri and Da. Although different surfaces of the cylinder has
non-dimensional permeability switches from 10
3 to 10
2 . The shown diverse effects on heat transfer performance with the vari-
inverse effect on heat transfer enhancement can be observed at ations of Darcy number, the mean Nusselt number increases with
rear, left and right surfaces of the cylinder with the increase in the increase of the Da. At low values of Re (2 6 Re 6 10), there is
Da. At higher permeability levels, from the front surface of the no significant variation in NuM for all values of Da; in this range
cylinder, more fluid penetrates and carries heat away from it; this of Re, only less fluid enters into the cylinder even at higher Da
 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870 867

1 3.3 -4
Da = 10 s1
Da = 10 -6 Da = 10-3 Ri = 0
-2
-4 Da = 10 f r
Da = 10 3 -4
Da = 10
Da = 10-3 Da = 10-3 Ri = 0.5
0.8 -2 s2
Da = 10-2 Da = 10
2.7 Da = 10
-4

Da = 10-3 Ri = 1
-2
Da = 10
2.4
0.6

Ef
2.1
θ

0.4 1.8

1.5
0.2
1.2

(a) (a) 0.9

0 5 10 15 20 25 30 35 40
0 240 480 720 960
X Re
0.8

Ri = 0 1
Ri = 0.5
Ri = 1
0.6
0.8
Er

0.4 0.6

0.4
0.2

0.2
(b) (b)
0 5 10 15 20 25 30 35 40
0 240 480 720 960
Re
X
1.05
Fig. 11. Influence of (a) Darcy number at Ri = 0 and (b) Richardson number at
Da ¼ 10
3 , on the non-dimensional temperature (h) profile along the direction of x
at y = 11.5H for Re = 40. 1

0.95
values due to low momentum of the fluid particles. Heat transfer
rate is more prominent at higher values of Re for variation of Da 0.9
Es1 = Es2

and the same has been reported by Dhinakaran and Ponmozhi

[51]. One can observe from Fig. 13(a–c) that the NuM value of
0.85
Da ¼ 10
4 is close to Da ¼ 10
6 at Ri = 0 and the gap between these
two Da levels increases with the increase in Richardson number.
For example, the percentage increase in mean Nusselt number at 0.8

4
Re = 40, Da = 10 at Ri = 0,0.5 and 1 is 2.04%, 2.2% and 3.03%
respectively, compared to solid cylinder. The reduction in thermal 0.75
boundary layer due to higher velocity gradients at higher Ri is the
cause for such heat transfer enhancement. From this section, we
0.7
can conclude that, the heat transfer performance can be aug-
(c)
mented by increasing permeability as well as aiding buoyancy. 5 10 15 20 25 30 35 40
The following nonlinear regression correlations are established Re
for the calculation of Nusselt number for the intermediate values
Fig. 12. Effect of Darcy number on heat transfer enhancement ratio at different
of Reynolds number (2 6 Re 6 40) and Darcy number
Reynolds number and Richardson number on (a) front face (b) rear face and (c) left/

6
2
(10 6 Da 6 10 ). right face of the heated porous square cylinder in an unbounded domain.
868 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870

pffiffiffiffiffiffi
4.5 NuM ¼ 1:258  Re0:696  Da þ 0:116 þ 0:761 for Ri ¼ 0:
Ri = 0 ð35Þ
4 Sharma et al. (solid cylinder) [29]
Da = 10-6
Da = 10-4 pffiffiffiffiffiffi
Dhinakaran and Ponmozhi (Da = 10-3) [51] NuM ¼ 1:260  Re0:739  Da þ 0:118 þ 0:835 for Ri ¼ 0:5:
3.5 Da = 10-3
Da = 10-2
ð36Þ

3 pffiffiffiffiffiffi
NuM ¼ 1:305  Re0:784  Da þ 0:103 þ 0:902 for Ri ¼ 1:
NuM

2.5 ð37Þ
The above correlations (Eqs. (35)–(37)) have the correlation
2 coefficient (R2 ) values of 0.995, 0.995 and 0.993 with the computed
numerical data for the Richardson number of 0, 0.5 and 1, respec-
tively. However, a common expression for different values of Ri
1.5
with maximum deviation of 8.77% is also been formulated as

(a) 0:444
pffiffiffiffiffiffi
1 NuM ¼ 1:062  Re0:746  ½1 þ Ri  Da þ 0:111 þ 0:839:
5 10 15 20 25 30 35 40
Re ð38Þ
5.5

Ri = 0.5 5. Conclusion
5

4.5
The flow and heat transfer characteristics of a 2-D porous
square cylinder placed in a uniform flow has been analysed in
4 detail under the influence of aiding buoyancy using lattice Boltz-
mann method. Darcy-Forchheimer terms that corresponds to vis-
3.5 cous and inertial effects due to porous medium, is coupled with
NuM

the SRT collision equation of LBM for porous zone simulation.

3 The outcomes of this numerical investigation can be summarized
as follows:
2.5
 At lower values of Da (6 10
6 ), the permeable cylinder behaves
2 like a solid cylinder as the flow restriction is more. Rich perme-
able cylinder reduces the recirculation zone that forms behind
1.5 the cylinder and at higher Da (=10
2 ) no wake forms behind
(b) the cylinder. The inertial forces inside the porous cylinder
1 increases as Da increases; which makes the flow less deviated
5 10 15 20 25 30 35 40
Re inside the porous zone at higher values of Da. The vortex diffu-
sion is more inside the higher permeable cylinder due to the
6
decrement of resistance to the flow and higher inertial force
5.5
Ri = 1 of the porous zone.
 For the same level of Da and Re, the deviation of fluid is found to
5 be more elongated inside the cylinder at higher values of Ri > 0.
The higher velocity gradients due to the aiding buoyancy,
4.5 stretches the fluid and this enables more fluid to penetrate
inside the cylinder for a constant values of Da and Re. The signif-
4
icance of permeability variation on drag coefficient and wake
length is less than that of buoyancy variations that is considered
NuM

3.5
in this study.
3  The mean Nusselt number increases with the increase of Da and
this increment is found to be more significant between 10
3 and
10
2 irrespective to Re and buoyancy condition. The stretched
2.5

2
and narrow thermal field behind the porous cylinder corre-
sponds for the reduction of thermal boundary layer at higher
1.5 Da and buoyancy levels and hence, heat transfer rate is found
(c) to be maximum at this levels. At low Re, there is no significant
1 change in heat transfer rate with Da and Ri variations.
5 10 15 20 25 30 35 40
 Finally, to quantify the mean Nusselt number at different buoy-
Re
ancy levels, three correlations as a function of Re and Da are
Fig. 13. Effect of Darcy number on mean Nusselt number (NuM ) at different obtained. Also, a common correlation that accounts the varia-
Reynolds numbers for flow through and around a porous square cylinder in an tion of Ri is formulated with correlation coefficient and maxi-
unbounded domain at (a) Ri = 0, (b) Ri = 0:5 and (c) Ri =1. mum deviation being 0.994 and 8.77%, respectively.
 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870 869

From this study, it is clear that the choice of permeability value [21] I.V. Morenko, V.L. Fedyaev, E.R. Galimov, Cross-flow and heat transfer
of porous permeable cylinder, IOP Conf. Ser.: Mater. Sci. Eng. 86 (2015)
and aiding buoyancy condition plays an important role to tune the
012012.
flow and heat transfer characteristics for practical applications. The [22] S. Chikh, A. Boumedien, K. Bouhadef, G. Lauriat, Analysis of fluid flow and heat
heat transfer magnification is more pronounced in case of channel transfer in a channel with intermittent heated porous blocks, Heat Mass
flow [18], cross buoyancy flow and it is expected to be more signif- Transfer 33 (5–6) (1998) 405–413.
[23] S. Rashidi, M. Bovand, I. Pop, M.S. Valipour, Numerical simulation of forced
icant for a completely heated permeable square cylinder with dif- convective heat transfer past a square diamond-shaped porous cylinder,
ferent incident angle of flow. Hence, a detailed study on cross Transport Porous Media 102 (2) (2014) 207–225.
buoyancy heat transfer characteristics around a porous square [24] G. Hwang, C. Chao, Effects of wall conduction and Darcy number on laminar
mixed convection in a horizontal square porous channel, J. Heat Transfer 114
cylinder in a confined channel with various incident angle of flow (3) (1992) 614–621.
can be carried out in future. [25] J.H. Bae, J.M. Hyun, J.W. Kim, Mixed convection in a channel with porous
multiblocks under imposed thermal modulation, Numer. Heat Transfer, Part A:
Appl. 46 (9) (2004) 891–908.
Conflict of interest [26] F. Chou, C. Cheng, W. Lien, Analysis and experiment of non-Darcian convection
in horizontal square packed-sphere channels-II. Mixed convection, Int. J. Heat
Mass Transfer 35 (5) (1992) 1197–1207.
We wish to confirm that there are no known conflicts of [27] I. Kurtbas, N. Celik, Experimental investigation of forced and mixed convection
interest. heat transfer in a foam-filled horizontal rectangular channel, Int. J. Heat Mass
Transfer 52 (5) (2009) 1313–1325.
[28] H.-W. Wu, R.-H. Wang, Mixed convective heat transfer past a heated square
Acknowledgement porous cylinder in a horizontal channel with varying channel height, J. Heat
Transfer 133 (2) (2011) 022503.
[29] N. Sharma, A.K. Dhiman, S. Kumar, Mixed convection flow and heat transfer
One of the authors, S. Dhinakaran, gratefully acknowledges the
across a square cylinder under the influence of aiding buoyancy at low
fund received from Council of Scientific & Industrial Research Reynolds numbers, Int. J. Heat Mass Transfer 55 (9) (2012) 2601–2614.
(CSIR), Government of India, through a project grant (Project Refer- [30] N. Guerroudj, H. Kahalerras, Mixed convection in a channel provided with
ence No. 22(0642)/13/EMR-II) for carrying out this work. heated porous blocks of various shapes, Energy Convers. Manage. 51 (3) (2010)
505–517.
[31] N. Guerroudj, H. Kahalerras, Mixed convection in an inclined channel with
References heated porous blocks, Int. J. Numer. Meth. Heat Fluid Flow 22 (7) (2012) 839–
861.
[32] A.A. Mohamad, Lattice Boltzmann Method: Fundamentals and Engineering
[1] D.B. Ingham, I. Pop, Transport Phenomena in Porous Media, Elsevier, UK, 1998.
Applications with Computer Codes, Springer Science & Business Media,
[2] I. Pop, D.B. Ingham, Transport Phenomena in Porous Media II, Elsevier, UK,
London Dordrecht Heidelberg New York, 2011.
2002.
[33] A. Karimipour, A.H. Nezhad, A. D’Orazio, E. Shirani, The effects of inclination
[3] D. Das, S. Hassanizadeh, Upscaling Multiphase Flow in Porous Media, vol. 58,
angle and Prandtl number on the mixed convection in the inclined lid driven
Springer, The Netherlands, 2005.
cavity using lattice Boltzmann method, J. Theor. Appl. Mech. 51 (2) (2013)
[4] D.A. Nield, A. Bejan, Convection in Porous Media, Springer Science & Business
447–462.
Media, New York Heidelberg Dordrecht London, 2006.
[34] G. Kefayati, Mesoscopic simulation of double-diffusive mixed convection of
[5] P. Vadász, Emerging Topics in Heat and Mass Transfer in Porous Media: From
Pseudoplastic Fluids in an enclosure with sinusoidal boundary conditions,
Bioengineering and Microelectronics to Nanotechnology, vol. 22, Springer
Comput. Fluids 97 (2014) 94–109.
Science & Business Media, USA, 2008.
[35] G. Kefayati, Mesoscopic simulation of mixed convection on non-Newtonian
[6] M.J. De Lemos, Turbulence in Porous Media: Modeling and Applications,
nanofluids in a two sided lid-driven enclosure, Adv. Powder Technol. 26 (2)
second ed., Elsevier, USA, 2012.
(2015) 576–588.
[7] K. Vafai, Handbook of Porous Media, second ed., Crc Press, Taylor & Francis
[36] G.R. Kefayati, FDLBM simulation of double-diffusive mixed convection of
Group, Boca Raton, 2015.
shear-thinning fluids between two-square concentric duct annuli, Heat Mass
[8] M. Dehghan, M.S. Valipour, S. Saedodin, Temperature-dependent conductivity
Transfer 51 (11) (2015) 1505–1521.
in forced convection of heat exchangers filled with porous media: a
[37] A. Karimipour, M.H. Esfe, M.R. Safaei, D.T. Semiromi, S. Jafari, S. Kazi, Mixed
perturbation solution, Energy Convers. Manage. 91 (2015) 259–266.
convection of copper–water nanofluid in a shallow inclined lid driven cavity
[9] J. Ma, Y. Sun, B. Li, H. Chen, Spectral collocation method for radiative–
using the lattice Boltzmann method, Phys. A: Stat. Mech. Appl. 402 (2014)
conductive porous fin with temperature dependent properties, Energy
150–168.
Convers. Manage. 111 (2016) 279–288.
[38] A. Karimipour, A.H. Nezhad, A. D’Orazio, M.H. Esfe, M.R. Safaei, E. Shirani,
[10] B. Kundu, K.-S. Lee, A proper analytical analysis of annular step porous fins for
Simulation of copper–water nanofluid in a microchannel in slip flow regime
determining maximum heat transfer, Energy Convers. Manage. 110 (2016)
using the lattice Boltzmann method, Euro. J. Mech. - B/Fluids 49 (2015) 89–99.
469–480.
[39] A. Karimipour, A.H. Nezhad, A. D’Orazio, E. Shirani, Investigation of the gravity
[11] B. Wang, Y. Hong, L. Wang, X. Fang, P. Wang, Z. Xu, Development and
effects on the mixed convection heat transfer in a microchannel using lattice
numerical investigation of novel gradient-porous heat sinks, Energy Convers.
Boltzmann method, Int. J. Therm. Sci. 54 (2012) 142–152.
Manage. 106 (2015) 1370–1378.
[40] Z. Nikkhah, A. Karimipour, M.R. Safaei, P. Forghani-Tehrani, M. Goodarzi, M.
[12] G. Kefayati, Heat transfer and entropy generation of natural convection on
Dahari, S. Wongwises, Forced convective heat transfer of water/functionalized
non-Newtonian nanofluids in a porous cavity, Powder Technol. 299 (2016)
multi-walled carbon nanotube nanofluids in a microchannel with oscillating
127–149.
heat flux and slip boundary condition, Int. Commun. Heat Mass Transfer 68
[13] G. Kefayati, Simulation of double diffusive natural convection and entropy
(2015) 69–77.
generation of power-law fluids in an inclined porous cavity with Soret and
[41] A. Karimipour, New correlation for Nusselt number of nanofluid with Ag/
Dufour effects (Part I: Study of fluid flow, heat and mass transfer), Int. J. Heat
Al2O3/Cu nanoparticles in a microchannel considering slip velocity and
Mass Transfer 94 (2016) 539–581.
temperature jump by using lattice Boltzmann method, Int. J. Therm. Sci. 91
[14] X. Chen, P. Yu, S. Winoto, H.-T. Low, Numerical analysis for the flow past a
(2015) 146–156.
porous square cylinder based on the stress-jump interfacial-conditions, Int. J.
[42] M. Breuer, J. Bernsdorf, T. Zeiser, F. Durst, Accurate computations of the
Numer. Meth. Heat Fluid Flow 18 (5) (2008) 635–655.
laminar flow past a square cylinder based on two different methods: lattice-
[15] T.-C. Jue, Numerical analysis of vortex shedding behind a porous square
Boltzmann and finite-volume, Int. J. Heat Fluid Flow 21 (2) (2000) 186–196.
cylinder, Int. J. Numer. Meth. Heat Fluid Flow 14 (5) (2004) 649–663.
[43] V. Babu, A. Narasimhan, Investigation of vortex shedding behind a porous
[16] P. Yu, Y. Zeng, T. Lee, H. Bai, H. Low, Wake structure for flow past and through a
square cylinder using lattice Boltzmann method, Phys. Fluids 22 (5) (2010)
porous square cylinder, Int. J. Numer. Meth. Heat Fluid Flow 31 (2) (2010)
053605 (1994-present).
141–153.
[44] Q. Liu, Y.-L. He, Q. Li, W.-Q. Tao, A multiple-relaxation-time lattice Boltzmann
[17] B. Çuhadaroğlu, O. Turan, Numerical simulation of turbulent flow around a
model for convection heat transfer in porous media, Int. J. Heat Mass transfer
square cylinder with uniform injection or suction and heat transfer, Numer.
73 (2014) 761–775.
Heat Transfer, Part A: Appl. 55 (2) (2009) 163–184.
[45] Z. Guo, T. Zhao, Lattice Boltzmann model for incompressible flows through
[18] S.-W. Perng, H.-W. Wu, R.-H. Wang, T.-C. Jue, Unsteady convection heat
porous media, Phys. Rev. E 66 (3) (2002) 036304.
transfer for a porous square cylinder varying cylinder-to-channel height ratio,
[46] C.-Y. Zhao, L. Dai, G. Tang, Z. Qu, Z. Li, Numerical study of natural convection in
Int. J. Therm. Sci. 50 (10) (2011) 2006–2015.
porous media (metals) using Lattice Boltzmann Method (LBM), Int. J. Heat
[19] Y.-L. Yen, P.-C. Huang, C.-F. Yang, Y.-J. Chen, Numerical study of heat transfer of
Fluid Flow 31 (5) (2010) 925–934.
a porous-block-mounted heat source subjected to pulsating channel flow,
[47] Z. Guo, T. Zhao, A lattice Boltzmann model for convection heat transfer in
Numer. Heat Transfer, Part A: Appl. 54 (4) (2008) 426–449.
porous media, Numer. Heat Transfer, Part B 47 (2) (2005) 157–177.
[20] S.-W. Perng, H.-W. Wu, T.-C. Jue, Numerical investigation of heat transfer
[48] T. Seta, E. Takegoshi, K. Okui, Lattice Boltzmann simulation of natural
enhancement on a porous vortex-generator applied to a block-heated channel,
convection in porous media, Math. Comput. Simul. 72 (2) (2006) 195–200.
Int. J. Heat Mass Transfer 55 (11–12) (2012) 3121–3137.
870 T.R. Vijaybabu et al. / International Journal of Heat and Mass Transfer 115 (2017) 854–870

[49] P.A.K. Lam, K.A. Prakash, A numerical investigation of heat transfer and [56] X. He, L.-S. Luo, Theory of the lattice Boltzmann method: from the
entropy generation during jet impingement cooling of protruding heat sources Boltzmann equation to the lattice Boltzmann equation, Phys. Rev. E 56 (6)
without and with porous medium, Energy Convers. Manage. 89 (2015) 626– (1997) 6811.
643. [57] A. Artoli, A. Hoekstra, P. Sloot, 3D pulsatile flow with the lattice Boltzmann
[50] W. Minkowycz, A. Haji-Sheikh, K. Vafai, On departure from local thermal BGK method, Int. J. Modern Phys. C 13 (08) (2002) 1119–1134.
equilibrium in porous media due to a rapidly changing heat source: the [58] A. Artoli, A. Hoekstra, P. Sloot, Accelerated Lattice BGK method for unsteady
Sparrow number, Int. J. Heat Mass Transfer 42 (18) (1999) 3373–3385. simulations through Mach number annealing, Int. J. Modern Phys. C 14 (06)
[51] S. Dhinakaran, J. Ponmozhi, Heat transfer from a permeable square cylinder to (2003) 835–845.
a flowing fluid, Energy Convers. Manage. 52 (5) (2011) 2170–2182. [59] A. Sharma, V. Eswaran, Heat and fluid flow across a square cylinder in the two-
[52] W.-S. Fu, H.-C. Huang, W.-Y. Liou, Thermal enhancement in laminar channel dimensional laminar flow regime, Numer. Heat Transfer, Part A: Appl. 45 (3)
flow with a porous block, Int. J. Heat Mass Transfer 39 (10) (1996) 2165–2175. (2004) 247–269.
[53] Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice [60] B. Paliwal, A. Sharma, R. Chhabra, V. Eswaran, Power law fluid flow past a
Boltzmann BGK model, Phys. Fluids 9 (6) (1997) 1591–1598 (1994-present). square cylinder: momentum and heat transfer characteristics, Chem. Eng. Sci.
[54] Y. Zhang, R. Qin, Y. Sun, R.W. Barber, D. Emerson, Gas flow in microchannels–a 58 (23) (2003) 5315–5329.
lattice Boltzmann method approach, J. Stat. Phys. 121 (1–2) (2005) 257–267. [61] D. Yu, R. Mei, L.-S. Luo, W. Shyy, Viscous flow computations with the
[55] P.D. Noymer, L.R. Glicksman, A. Devendran, Drag on a permeable cylinder in method of lattice Boltzmann equation, Prog. Aerospace Sci. 39 (5) (2003)
steady flow at moderate Reynolds numbers, Chem. Eng. Sci. 53 (16) (1998) 329–367.
2859–2869.