Effects of buoyancy ratio on unsteady double-diffusive natural convection in a cavity filled with porous medium with non-uniform boundary conditions

International Journal of Heat and Mass Transfer 85 (2015) 401–413

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Effects of buoyancy ratio on unsteady double-diffusive natural convection in a cavity filled with porous medium with non-uniform boundary conditions Sabyasachi Mondal ⇑, Precious Sibanda School of Mathematics, Statistics and Computer Science, University of KwaZulu-Natal, Scottsville, Pietermaritzburg, South Africa

a r t i c l e

i n f o

Article history: Received 8 September 2014 Received in revised form 24 January 2015 Accepted 27 January 2015

Keywords: Double-diffusive convection Finite-difference method Non-uniformly heated walls Non-uniformly concentrated walls Porous media

a b s t r a c t The effects of buoyancy ratio on unsteady double-diffusive natural convection in a cavity filled with porous medium with uniform and non-uniform boundary conditions are analyzed in this paper. It is assumed that the left vertical wall and bottom wall are heated and concentrated (uniformly and nonuniformly), while the right vertical wall is maintained at a constant cold temperature, and the top wall is well insulated. The governing equations are solved numerically using a staggered grid finite-difference method to determine the streamlines, isotherms, isoconcentrations, local Nusselt number, local Sherwood number, average Nusselt number and average Sherwood number for various values of buoyancy ratio and Rayleigh number. The change of flow patterns with respect to time depicted and described here. The results are compared with previously published work and excellent agreement has been obtained. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Studies of flow through porous medium have attracted considerable research attention in recent years because of their several important applications notably in the flow through packed beds, extraction of energy from geothermal regions, filtration of solids from liquids, flow of liquids through ion-exchange beds, the evaluation of the capability of heat removal from nuclear fuel debris that may result from an accident in a nuclear reactor and in chemical reactors for the separation or purification of mixtures [1–4]. Fluid flow, heat and mass transfer induced by double-diffusive natural convection in fluid saturated porous media have practical importance in many engineering applications [5]. This aspect of fluid dynamics has gained considerable attention in recent years among researchers. The migration of moisture in fibrous insulation, drying processes, chemical reactors, transport of contaminants in saturated soils and electro-chemical processes are some examples of double-diffusive natural convection phenomena. Double-diffusion occurs in a wide range of scientific fields such as oceanography, astrophysics, geology, biology and chemical processes. In the recent past, a significant number of researchers have shown a keen interest in the study of heat and mass transfer in ⇑ Corresponding author. E-mail addresses: [email protected] (S. Mondal), sibandap@ukzn. ac.za (P. Sibanda). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.01.129 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.

enclosures and cavities. Double-diffusive natural convection in cavities has been subject to an intensive research due to its importance in engineering and geophysical problems. This includes nuclear reactors, solar ponds, geothermal reservoirs, solar collectors, crystal growth, electronic cooling and chemical processing equipments. The numerical investigation of natural convection in porous trapezoidal enclosures has been performed for uniformly or non-uniformly heated bottom wall by Basak et al. [6]. Sathiyamoorthy et al. [7] studied non-Darcy buoyancy flow in a square cavity filled with porous medium for various temperature difference aspect ratios. Deng et al. [8] investigated fluid, heat and contaminant transport structures of laminar double diffusive mixed convection in a two-dimensional ventilated enclosure numerically. Roy et al. [9] performed finite element simulation on natural convection flow in a triangular enclosure due to uniform and non-uniform bottom heating. Later, Alimi et al. [10] studied the buoyancy effects on mixed convection heat and mass transfer in an inclined duct preceded with a double step expansion. Brown and Lai [11] numerically examined combined heat and mass transfer from a horizontal channel with an open cavity heated from below numerically. Teamah [12] studied double-diffusive convective flow in a rectangular enclosure with the upper and lower surfaces being insulated and impermeable by imposing constant temperature and concentration along the left and right walls of the enclosure and a uniform magnetic field was applied in a horizontal direction. Saha et al.

402

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413

[13] investigated the new characteristics of the airflow and heat/contaminant transport mechanism inside a vented cavity in terms of streamlines, isotherms and isoconcentration lines. Minkowycz et al. [14] showed that the discontinuity can be avoided by choosing a non-uniform temperature distribution along the walls (i.e. non-uniformly heated walls). Roy and Basak [15] solved the nonlinear coupled partial differential equations for flow and temperature fields with both uniform and non-uniform temperature distributions prescribed at the bottom wall and at one vertical wall. Teamah et al. [16] studied the effect of the heater length, Rayleigh number, Prandtl number and buoyancy ratio on both average Nusselt and Sherwood number with uniform heating at left vertical wall. Karimi-Fard et al. [17] studied double diffusive natural convection in a cavity filled with a porous medium. Patil et al. [18] investigated double diffusive mixed convection flow over a vertical plate. Nithiarasu et al. [19] studied the development of the variable porosity model in natural convection heat transfer in detail. Recently, Mahapatra et al. [20] investigated the effects of buoyancy ratio and the thermal Rayleigh number on double diffusive natural convection in a cavity when the boundaries are uniformly and non-uniformly heated and concentrated. The aim of this investigation is to study the effects of buoyancy ratio and Rayleigh number on the heated and concentrated walls in terms of streamlines, isotherms, isoconcentrations, local Nusselt number, average Nusselt number, local Sherwood number and average Sherwood number when the bottom wall and left vertical wall are heated and concentrated (uniformly and non-uniformly), right vertical wall is cooled by means of a constant temperature and top wall is well insulated. The thermal and mass exchanges generated in the case of co-operating thermal and concentration buoyancy effects with uniform and non-uniform boundary conditions have been analyzed. 2. Governing equations and boundary conditions An unsteady-state two-dimensional square cavity of height L as shown in Fig. 1 is considered. It is assumed that the top wall is considered to be adiabatic. The bottom wall and left vertical wall are heated and concentrated (uniformly and non-uniformly) and right vertical wall is cooled by means of a constant temperature. The thermophysical properties of the fluid are assumed to be constant except the density variation in the buoyancy force, which is approximated according to the Boussinesq approximation. This variation, due to both temperature and concentration gradients, can be described by the following equation:

! @T @T @T @2T @2T ; þU þV ¼a þ @t0 @X @Y @X 2 @Y 2

ð5Þ

! @C @C @C @2C @2C : þU þV ¼D þ @t 0 @X @Y @X 2 @Y 2

ð6Þ

The associated boundary conditions are when t0 ¼ 0 for 0 6 X; Y 6 L:

UðX; YÞ ¼ 0 ¼ VðX; YÞ;

ð7Þ

TðX; YÞ ¼ T c ;

ð8Þ

when t > 0 for 0 6 X; Y 6 L:

UðX; LÞ ¼ UðX; 0Þ ¼ Uð0; YÞ ¼ UðL; YÞ ¼ 0;

ð9Þ

VðX; 0Þ ¼ VðX; LÞ ¼ Vð0; YÞ ¼ VðL; YÞ ¼ 0;

ð10Þ

ð1Þ

where bT and bS are the thermal and concentration expansion coefficients, respectively. In the cartesian coordinate system, the fundamental governing equations are as follows (see Marcondes et al. [21], Mahapatra et al. [22] and Mahapatra et al. [23]):

@U @V þ ¼ 0; @X @Y

ð2Þ

@U U @U V @U 1 @P @2U @2U þ ¼ þm þ 0 þ @t  @X  @Y q @X @X 2 @Y 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1:75 U 2 þ V 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U;  150K  @V U @V V @V 1 @P þ þ ¼ þm @t0  @X  @Y q @Y

! 

m K

U

ð3Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi @ V @ V m 1:75 U 2 þ V 2  V  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V þ K 150K  @X 2 @Y 2 2

2

!

þ gbðT  T c Þ þ gbS ðC  C c Þ;

ð4Þ

@T ðX;LÞ ¼ 0; @Y

TðX;0Þ ¼ T h or TðX;0Þ ¼ ðT h  T c ÞsinðpX=LÞ þ T c ;

ð11Þ Tð0;YÞ ¼ T h or Tð0; YÞ ¼ ðT h  T c ÞsinðpY=LÞ þ T c ; TðL;YÞ ¼ T c ; ð12Þ CðX; 0Þ ¼ C h or CðX; 0Þ ¼ ðC h  C c ÞsinðpX=LÞ þ C c ;

@C ðX; LÞ ¼ 0; @Y ð13Þ

Cð0; YÞ ¼ C h or Cð0; YÞ ¼ ðC h  C c ÞsinðpY=LÞ þ C c ; CðL; YÞ ¼ C c : ð14Þ where, X and Y are the distances measured along the horizontal and vertical directions respectively; U and V are velocity components in the X- and Y- directions respectively; T and C denote the temperature and concentration respectively; m; a and D are kinematic viscosity, thermal diffusivity and mass diffusivity respectively; P is the pressure and q is the density; T h and T c are the temperatures at the hot and cold walls respectively; C h and C c are the concentrations at the hot and cold walls respectively; L is the side of the square cavity. We now introduce dimensionless variables given as follows:

t¼

h¼

q ¼ q0 ½1  bT ðT  T c Þ  bS ðC  C c Þ

CðX; YÞ ¼ C c ;

0

at0 L2

;

x¼

T  Tc ; Th  Tc

X ; L

y¼

Y ; L

S¼

C  Cc : Ch  Cc

u¼

UL

a

;

v¼

VL

a

;

p¼

PL2

qa2

;

ð15Þ

ð16Þ

Here x and y are dimensionless coordinates along the horizontal and vertical directions respectively; u and v are dimensionless velocity components in the x- and y- directions respectively; h and S denote the dimensionless temperature and concentration respectively; p is the dimensionless pressure parameter. Using these dimensionless variables, we obtain the following dimensionless governing equations from the Eqs. (2)–(6):

@u @ v þ ¼ 0; @x @y

ð17Þ

! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   @u @p @2u @2u 1 @u2 @uv Pr  1:75 u2 þ v 2 p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u; ¼  þ Pr þ þ u     @x @y Da @t @x @x2 @y2 150Da

ð18Þ

!   @v @p @2v @2v 1 @ v 2 @uv Pr  ¼  þ Pr þ þ v  @y @t @x2 @y2  @y @x Da pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1:75 u2 þ v 2  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v þ PrRaðh þ NSÞ; 150Da

ð19Þ

403

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413

Fig. 1. Schematic diagram of the physical system.

@h ¼ @t

!   @2h @2h @uh @ v h  ; þ þ @x2 @y2 @x @y

@S 1 @2S @2S ¼ þ @t Le @x2 @y2

!

 

ð20Þ

 @uS @ v S : þ @x @y

ð21Þ

The dimensionless boundary conditions are when t ¼ 0 for 0 6 x; y 6 1:

uðx; yÞ ¼ 0 ¼ v ðx; yÞ; hðx; yÞ ¼ 0;

Here N; Le; Pr; GrC ; GrT ; Ra and Da are buoyancy ratio, Lewis number, Prandtl number, solutal Grashof number, thermal Grashof number, Rayleigh number and Darcy number respectively. The heat transfer coefficient in terms of the local Nusselt number ðNuÞ is Nu ¼ ð@hÞ=ð@nÞ, where n denotes the normal direction to a plane. The local Nusselt number at the left vertical wall and bottom walls are defined as Nul ¼ ð@hÞ=ð@xÞjx¼0 and Nub ¼ ð@hÞ=ð@yÞjy¼0 . The average Nusselt numbers at the bottom wall and left vertical wall are given by,

ð22Þ

Sðx; yÞ ¼ 0;

ð23Þ

when t > 0 for 0 6 x; y 6 1:

uðx; 1Þ ¼ uðx; 0Þ ¼ uð0; yÞ ¼ uð1; yÞ ¼ 0;

ð24Þ

v ðx; 0Þ ¼ v ðx; 1Þ ¼ v ð0; yÞ ¼ v ð1; yÞ ¼ 0;

ð25Þ

NuH jy¼0 ¼

Z

Z

1

0

Nub dx and NuH jx¼0 ¼

1

Nul dy:

ð32Þ

0

The local Sherwood number at the left vertical wall and bottom walls are defined as Shl ¼ ð@SÞ=ð@xÞjx¼0 and Shb ¼ ð@SÞ=ð@yÞjy¼0 . The average Sherwood numbers at he bottom wall and left vertical wall are given by,

ShH jy¼0 ¼

Z

1

Shb dx and ShH jx¼0 ¼

Z

1

Shl dy:

ð33Þ

hðx; 0Þ ¼ 1 or hðx; 0Þ ¼ sinðpxÞ;

@h ðx; 1Þ ¼ 0; @y

ð26Þ

hð0; yÞ ¼ 1 or hð0; yÞ ¼ sinðpyÞ;

hð1; yÞ ¼ 0;

ð27Þ

3. Method of solution

Sðx; 0Þ ¼ 1 or Sðx; 0Þ ¼ sinðpxÞ;

@S ðx; 1Þ ¼ 0; @y

ð28Þ

Sð0; yÞ ¼ 1 or Sð0; yÞ ¼ sinðpyÞ;

Sð1; yÞ ¼ 0:

ð29Þ

Control-volume based finite-difference discretization of the continuity, momentum, temperature and concentration equations are carried out using staggered grid, popularly known as MAC cell method which was also studied by Mahapatra et al. [24]. The u and v velocity components are evaluated at different locations of the control volume whereas the pressure, temperature and concentration are evaluated at same location of the control volume as shown in Fig. 2. It can be noted from the Fig. 2 that u; v velocity components are stored at the mid point of the vertical and horizontal faces respectively whereas pressure, temperature and concentration

Where

GrT ¼ gbT ðT h  T c ÞL3 =m2 ; N¼

GrC ; GrT

Pr ¼

GrC ¼ gbS ðC h  C c ÞL3 =m2 ;

m K ; Ra ¼ Gr T Pr; Da ¼ 2 : a L

ð30Þ ð31Þ

0

0

404

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413

values are stored at the centre of the cells. In MAC method we use different cells to discretized different equations. While discretizing continuity, temperature and concentration equations, we use the cell ABCD as can be seen in Fig. 2. The cells EFGH and A0 B0 C 0 D0 are used to discretize u-momentum and v-momentum equations. Now, the non-dimensional Equations can be written in the following form:

@u @p ¼  þ UDPC; @t @x

ð34Þ

@v @p ¼  þ VDPC; @y @t

ð35Þ

The convective terms in the momentum equations are differenced with a hybrid formula consisting of central differencing and second-order upwinding. The diffusive terms are differenced by second order accurate three point central difference formula in both the cases. The source terms are evaluated at the centre of the corresponding control volumes. The finite difference forms of p; u; v ; h; S with t ¼ ndt; x ¼ idx; y ¼ jdy are

pðx; y; tÞ ¼ pðidx; jdy; ndtÞ ¼ pnij ;

ð42Þ

uðx; y; tÞ ¼ uðidx; jdy; ndtÞ ¼ unij ;

ð43Þ

v ðx; y; tÞ ¼ v ðidx; jdy; ndtÞ ¼ v

ð44Þ

n ij ;

hðx; y; tÞ ¼ hðidx; jdy; ndtÞ ¼ hnij ;

ð45Þ

Snij ;

ð46Þ

@h ¼ TDPC; @t

ð36Þ

Sðx; y; tÞ ¼ Sðidx; jdy; ndtÞ ¼

@S ¼ CDPC: @t

ð37Þ

where superscript n refers to the time direction, subscripts i and j refer to the spatial directions, dt is the time increment and dx; dy are the length and width of the control volume. Now discretization of the continuity equation at ði; jÞ cell takes the following form,

where,

!   @2u @2u 1 @u2 @uv Pr   UDPC ¼ Pr þ þ u @x2 @y2  @x @y Da pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1:75 u2 þ v 2  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 150Da !   @2v @2v 1 @ v 2 @uv Pr   VDPC ¼ Pr þ þ v @x2 @y2  @y @x Da pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1:75 u2 þ v 2  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v þ PrRaðh þ NSÞ 150Da !

  @uh @ v h  ; þ @x @y

TDPC ¼

@2h @2h þ @x2 @y2

CDPC ¼

1 @2S @2S þ Le @x2 @y2

!

 

 @uS @ v S : þ @x @y

unij  uni1j dx ð38Þ

þ

v nij  v nij1 dy

¼ 0:

ð47Þ

The finite-difference formulation of the momentum equation in xdirection using uniform grid-spacing is

unþ1  unij ij dt

¼

pniþ1j  pnij dx

þ ðUDPCÞnij ;

ð48Þ

where

ð39Þ

n

ðUDPCÞnij ¼ PrDiffuij 

Conunij

Pr  n  u sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Da ij

0 1 1:75unij ðv t þ v b Þ2 A 2 n @ ; ðuij Þ þ  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 150Da

ð40Þ

ð49Þ

n

ð41Þ

where, Diffuij ; Conunij are diffusive and convective terms of the u-momentum equation at the nth time level at ði; jÞth cell and v t ¼ 0:5ðv nij þ v niþ1j Þ; v b ¼ 0:5ðv nij1 þ v niþ1j1 Þ. In the y-direction, Xthe finite-difference formulation of the momentum equation is

v nþ1  v nij ij dt

¼

pnijþ1  pnij dy n

ðVDPCÞnij ¼ PrDiff v ij 

þ ðVDPCÞnij ;

ð50Þ

Conv nij

Pr  n  v sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Da ij ffi

0 1 2 1:75v nij ðu þ u Þ t b 2 n  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi @ þ ðv ij Þ A 4 150Da ! hnij þ hnijþ1 Snij þ Snijþ1 ; þ 0:5RaPr þN 2 2

ð51Þ

n

where, Diff v ij ; Conv nij are diffusive and convective terms of the v-momentum equation at the nth time level at ði; jÞth cell and ut ¼ 0:5ðunijþ1 þ uni1jþ1 Þ; ub ¼ 0:5ðunij þ uni1j Þ. Now, the discritized form of energy equation

hnþ1  hnij ij ¼ ðTDPCÞnij ; dt

ð52Þ

where

ðTDPCÞnij ¼ Diff hnij  Conhnij : Diff hnij ; Conhnij

Fig. 2. Control volume for u-momentum, v-momentum, temperature and concentration equations.

ð53Þ

Here, are diffusive and convective terms of the energy equation at the nth time level at ði; jÞth cell. Now, the discretized form concentration equation

405

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413 n

Snþ1  Snij ij ¼ ðCDPCÞnij ; dt

ð54Þ

where

n

n

@ 2 h hijþ1  2hij þ hij1 ¼ þ Oðdy2 Þ; @y2 dy2

ð69Þ

and convective terms in temperature equation are discretized as

1 n ðCDPCÞnij ¼ DiffSij  ConSnij : Le

ð55Þ

un hn  unl hnl un /n  unl /nhl @uh ¼ ð1  bÞ r r þ b r hr ; @x dx dx

ð70Þ

v n hn  v nb hnb v nt /nht  v nb /nhb @v h ¼ ð1  bÞ t t þb : @y dy dy

ð71Þ

n

Here, DiffSij and ConSnij denote diffusive and convective terms of the concentration equation at the nth time level at ði; jÞth cell. The diffusive terms in u-momentum equation are discretized as n

n

n

n

Thus the momentum fluxes for temperature equation are given by

n

@ 2 u uiþ1j  2uij þ ui1j ¼ þ Oðdx2 Þ; @x2 dx2

ð56Þ

n

@ 2 u uijþ1  2uij þ uij1 ¼ þ Oðdy2 Þ; @y2 dy2

ð57Þ

and convective terms in u-momentum equation are discretized as

un un  unl unl un /n  unl /nul @u2 ¼ ð1  bÞ r r þ b r ur ; @x dx dx

ð58Þ

un v n  unb v nb v n /n  v nb /nub @uv ¼ ð1  bÞ t t þ b t ut : @y dy dy

ð59Þ

if

unr P 0;

/nur ¼ unij ;

ð60Þ

if

unr < 0;

/nur ¼ uniþ1j ;

ð61Þ

if

unl P 0;

/nul ¼ uni1j ;

ð62Þ

if

unl < 0;

/nul ¼ unij ;

ð63Þ

    where unr ¼ 0:5 unij þ uniþ1j ; unl ¼ 0:5 uni1j þ unij . Suffix u of / represents the quantity for u-momentum. Similarly,

if if if

v nt P 0; v nt < 0; v nb P 0; v nb < 0;

/nut ¼ unij ;

ð64Þ

/nut ¼ unijþ1 ;

ð65Þ

/nub ¼ unij1 ; /nub ¼ unij ;

ð66Þ ð67Þ

where unt ¼ 0:5ðunij þ unijþ1 Þ; unb ¼ 0:5ðunij1 þ unij Þ; v nt ¼ 0:5ðv nij þ v niþ1j Þ and v nb ¼ 0:5ðv nij1 þ v niþ1j1 Þ. The finite-difference discretization of the diffusive and convective terms of the v -momentum equation are similar to those of the u-momentum equation. The diffusive terms in temperature equation are discretized as n

n

unr P 0;

/nhr ¼ hnij ;

ð72Þ

if

unr

< 0;

/nhr ¼ hniþ1j ;

ð73Þ

if

unl unl

P 0;

/nhl ¼ hni1j ; /nhl ¼ hnij :

ð75Þ

if

if if if

ð68Þ

ð74Þ

v nt P 0; v nt < 0; v nb P 0; v nb < 0;

/nht ¼ hnij ;

ð76Þ

/nht ¼ hnijþ1 ;

ð77Þ

/nhb ¼ hnij1 ; /nhb ¼ hnij :

ð79Þ

ð78Þ

  The convected momentum fluxes unt ; unb ; unr ; unl ; v nt ; v nb ; v nr ; v nl at the interfaces are linearly interpolated for both u-momentum, v -momentum and temperature equations. Here, suffixes t; b; r; l denote corresponding value at top, bottom, right and left middle positions of the suitable cell faces. In the second-order upwinding scheme, the choice of taking the momentum flux / passing through the interface of the control volume depends on the sign of convecting velocities at that interface. Here b is the combination factor between the central-differencing scheme and second order upwind-differentiating scheme. When b ¼ 0, the convective terms are represented by the central-difference whereas it is second-order upwinding when b ¼ 1. A pressure-Poisson equation is derived combing the discretized momentum and continuity equations. Here the local dilation term nþ1

at ðn þ 1Þth level Div ij is set equal to zero. Thus the final form of the pressure-Poisson equation is n

n

2ða þ bÞpnij  apniþ1j  apni1j  bpijþ1  bpij1 " n # o 1 Dij 1 n ¼ þ ðUDPCÞnij  ðUDPCÞni1j þ fðVDPCÞnij  ðVDPCÞni1j g ; dy dt dx ð80Þ Dnij

where is the divergence of the velocity field at the cell ði; jÞ at nth time level, which is to be zero for the convergence of the flow in successive iterations of the method and

n

@ 2 h hiþ1j  2hij þ hi1j ¼ þ Oðdx2 Þ; @x2 dx2

< 0;

Suffix h of / represents the quantity for temperature equation. Similarly,

if

Here b, the combination factor the central-differencing and second order upwind-differencing, is determined from the numerical stability criteria. When b ¼ 0, the convective terms are discretized by the central-difference, whereas it is second-order upwind scheme when b ¼ 1. In the second-order upwind scheme, the choice of taking the momentum flux / passing through the interface of the control volume depends on the sign of convecting velocities at that interface. Thus the momentum fluxes for u-momentum are given by

if

if

a¼

 n  uij  uni1j v nij  v nij1 1 1 : ; b ¼ 2 ; Dnij ¼ þ 2 dx dy dx dy

ð81Þ

Table 1 Grid independence test when Pr ¼ 0:7; Le ¼ 2:0; N ¼ 1; Da ¼ 103 and Ra ¼ 103 . No. of grid points

20  20 40  40 80  80 160  160

Uniformly heated and concentrated

Non-uniformly heated and concentrated

Iteration

jwmin j

Iteration

jwmin j

18707 74163 295857 1051253

0.113576 0.111777 0.111222 0.111212

18066 71804 287014 1074605

0.08735 0.08595 0.08552 0.08551

406

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413

Table 2 Comparison of average Nusselt number NuH jx¼0 in absence of concentration equation, thermal radiation and heat generation of the generalized model in the non-Darcian regime with same boundary conditions and N ¼ 0; Pr ¼ 1; Da ¼ 106 ;  ¼ 0:4. NuH jx¼0 Ra

Nithiarasu et al. [19]

Medeiros et al. [25]

Mahapatra et al. [26]

Present work

8

2.97

3.05

3.07

3.07

109

11.46

12.10

12.23

12.23

5  109

23.09

24.74

24.43

24.43

10

Table 3 Average Nusselt number and average Sherwood number at different walls with different time when Pr ¼ 0:7; Le ¼ 2:0; N ¼ 1; Da ¼ 103 and Ra ¼ 103 . Time

0.0279 0.1395 0.5580 1.6015

Average Nusselt number at wall

Average Sherwood number at wall

Bottom wall

Left vertical wall

Bottom wall

Left vertical wall

2.1798 1.1741 1.0383 1.0376

1.8843 0.5882 0.3375 0.3360

3.0565 1.4762 1.0785 1.0560

2.8391 0.9916 0.3636 0.3195

Fig. 3. Contour plots for Pr ¼ 0:7; Le ¼ 2:0; Da ¼ 103 ; Ra ¼ 103 and (a) N ¼ 50 (b) N ¼ 1 (c) N ¼ 50 (uniformly heated and concentrated).

3.1. Solution procedure and numerical stability criteria Now the iteration process is described to obtain the solutions of the basic equations with appropriate boundary conditions. In the derivation of pressure-Poisson equation, the divergence term at

nth time level ðDnij Þ is retained and evaluated in the pressure-Poisson iteration. It is done because the discretized form of divergence of velocity field, i.e, Dnij is not guaranteed to be zero initially. The solution procedure starts with the initializing the velocity field. This is done either from the result of previous cycle or from the

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413

407

Fig. 4. Contour plots for Pr ¼ 0:7; Le ¼ 2:0; Da ¼ 103 ; Ra ¼ 103 and (a) N ¼ 50 (b) N ¼ 1 (c) N ¼ 50 (non-uniformly heated and concentrated).

prescribed initial and boundary conditions. Using this velocity field pressure-Poisson equation is solved using Bi-CG-Stab method. Knowing the pressure field, equation for u-momentum, v -momentum, temperature and concentration are updated to get u; v ; h and S at ðn þ 1Þth time level. Then using the values of u and v at ðn þ 1Þth time level, the value of the divergence of velocity field is obtained and checked for its limit. If its absolute value is less than 0:5  105 and steady state is reached then iteration process stops, otherwise pressure-Poisson equation is solved again for pressure. h i dx dy Linear stability of fluid flow gives dt1 6 Min juj ; jv j , which is related to the convection of fluid, i.e., fluid should not move more than one cell width per time step (Courant, Friedrichs and Lewy condition). Also, from the Hirt’s stability analysis, we have h i 2 dy2 1 dt 2 6 Min 2Pr  ðdxdx2 þdy 2 Þ . This condition roughly stated that momentum cannot diffuse more than one cell width per time step. The time step is determined from dt ¼ FCT  ½Minðdt 1 ; dt 2 Þ, where the factor FCT varies from 0:2 to 0:4. The upwinding parameter b is h i governed by the inequality condition 1 P b P Max j udt j; j vdydt j . As dx a rule of thumb, b is taken approximately 1:2 times larger than what is found from the above inequality condition.

4. Results and discussion Grid independence test is provided in Table 1 for various grid sizes when Pr ¼ 0:7; Le ¼ 2:0; N ¼ 1; Da ¼ 103 and Ra ¼ 103 . It is important to note that as the number of grid points are

increased, the number of iterations required for getting the converged results for jwmin j drastically increases. But when the number of grid points increases from 80  80 to 160  160, no significant change found in the value of jwmin j. Hence, all the results are computed taking 80  80 grid points. A Comparison of the average Nusselt number at the hot wall ðx ¼ 0Þ in absence of concentration equation, thermal radiation and heat generation of the generalized model in the non-Darcian regime with same boundary conditions and N ¼ 0; Pr ¼ 1; Da ¼ 106 ;  ¼ 0:4 is made with Nithiarasu et al. [19], Medeiros et al. [25] and Mahapatra et al. [26] for various values of Rayleigh number in Table 2. It is noted from this table that an average Nusselt number increases with increase in the value of Rayleigh number and from this table it is observed that a very good agreement has been obtained with the previously published results. From Table 3 it is seen that average Nusselt number and average Sherwood number at bottom wall and left vertical wall are decrease as time increases for both uniform and non-uniform boundary conditions. Eqs. (83)–(86)show a least square curve fitting of average Nusselt number and average Sherwood number with Rayleigh number at bottom and left vertical walls for uniform and nonuniform boundary conditions. The fitted model is of the form,

NuH

or

ShH ¼ coefficient  Rapower

ð82Þ

The values of ‘‘coefficient’’ and ‘‘power’’ with Standard Error of Least Square Curve fitting formula of average Nusselt number or Sherwood number with Rayleigh number at bottom and left vertical

408

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413

Fig. 5. Contour plots for Pr ¼ 0:7; Le ¼ 2:0; Da ¼ 103 ; N ¼ 1 and (a) Ra ¼ 103 (b) Ra ¼ 104 (c) Ra ¼ 105 (uniformly heated and concentrated).

walls for uniform and non-uniform boundary conditions when Pr ¼ 0:7; Le ¼ 2:0; Da ¼ 103 ; Ra ¼ 103 and N ¼ 1 is found from equations (83)–(86). The following relation are obtained for cases a and b (for uniform boundary conditions) and c and d (for non-uniform boundary conditions) as follows:

Case a NuH ¼ 1:55214  Ra0:11195

ðat bottom wallÞ

NuH ¼ 0:07278  Ra0:19718

ðat left vertical wallÞ

ð83Þ

Case b

ShH ¼ 1:18323  Ra0:14943

ðat bottom wallÞ

ShH ¼ 0:01063  Ra0:40102

ðat left vertical wallÞ

ð84Þ

Case c

NuH ¼ 0:12477  Ra0:25513

ðat bottom wallÞ

NuH ¼ 0:04013  Ra0:18772 ð at left vertical wallÞ

ð85Þ

Case d

ShH ¼ 0:10754  Ra0:29169 ðat bottom wallÞ ShH ¼ 0:00105  Ra0:53969 ðat left vertical wallÞ

ð86Þ

We observe from equations (83)–(86) that the standard error of estimation is very low, so the relation obtained for average Nusselt number and average Sherwood number as a function of the

Rayleigh number would lead to almost the same values as obtained from a numerical computation of the set of differential equations with appropriate boundary conditions. Thus it is interesting to note without further computations that equations (83)–(86) can be used to estimate the average Nusselt and Sherwood numbers for various values of Rayleigh number. The numerical results for streamlines, isotherms, isoconcentrations, local Nusselt number, average Nusselt number, local Sherwood number and average Sherwood number are presented in Figs. 3–9 for uniformly and non-uniformly heated and concentrated walls. The relative importance of the thermal and solutal buoyancy forces is denoted by the buoyancy ratio (N), and is defined as the ratio of the solutal buoyancy force to the thermal buoyancy force. This parameter is varied through the range 50 < N < 50; covering the concentration-dominated opposing flow (N ¼ 50), pure thermal convection dominated flow (N ¼ 0), and concentration-dominated aiding flow (N ¼ 50). When N is sufficiently small i.e, the mass buoyancy is greater than the thermal buoyancy, the buoyancy forces that drive the fluid motion are mainly due to the gradients of temperature. Negative values of N represent the opposing nature of two buoyancy forces due to the negative coefficient of concentration expansion. In this limit, the so-called heattransfer-driven flows, the distribution of constituent does not influence the flow pattern and the heat transfer rate. When N ¼ 1, the flow is steady; this is because in this case, the two buoyancy forces are equal and oppose each other. When N > 1, the flows driven by buoyancy due to solutal gradients where the flow are mainly due to gradients of solute concentration. The effects of N on streamline, isotherm and isoconcentration contours for Pr ¼ 0:7; Le ¼ 2:0 Da ¼ 103 and Ra ¼ 103 are dis-

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413

409

Fig. 6. Contour plots for Pr ¼ 0:7; Le ¼ 2:0; Da ¼ 103 ; N ¼ 1 and (a) Ra ¼ 103 (b) Ra ¼ 104 (c) Ra ¼ 105 (non-uniformly heated and concentrated).

played in Figs. 3 and 4 for uniform and non-uniform boundary conditions. Fig. 3 shows the effect of N on the streamlines, isotherms and as well as on the isoconcentration for wide range of variations in the buoyancy ratio (N) with uniformly heated and concentrated wall. As expected due to a hot vertical wall, fluids rise up along the side of the hot vertical wall and flow down along the cold vertical wall forming a roll with clockwise rotation inside the cavity. As N increases from 50 to 1, the values of stream function decrease i.e, the strength of the flow circulation decreases. Again, when N is further increased from 1 to 50 the values of the stream function decreases, i.e, the strength of the flow circulation decreases with an increase in N. The streamlines are concentrated near the cavity wall and the cavity core is empty. It is clear from the Fig. 3(a) that when N ¼ 50 the contours of both the isotherms and isoconcentrations are mostly located near the lower half of the right cold vertical wall and the upper half of the left hot vertical wall indicating that most of the heat transfer is due to heat conduction. This is due to an increase in the thermal boundary layer thickness. However, when N is increased to 1, they are dispersed to the cold wall (see Fig. 3(b)). Again when the N is further increased to 50, both the isotherms and isoconcentrations contours are mainly located near the cold vertical wall (see Fig. 3(c)). At low buoyancy ratio, the entire enclosure is influenced by the flow structure and as the buoyancy ratio increases the boundary

layer thickness becomes thinner. This change in the flow structure for high buoyancy ratio has a significant influence on the concentration field, which builds up a vertical stratification of the enclosure. Again it interesting to note that for N ¼ 50, the effect of solutal buoyancy force is in the opposite direction of thermal buoyancy force. Therefore, we note that the magnitude of the thermal buoyancy force is very small compared to the solutal buoyancy force. The isothermal and isoconcentration tend toward the left wall, while for N ¼ 1 the thermal and solutal buoyancy forces are equal. The isothermal and isoconcentration contours tend toward the right wall and for N ¼ 50, the effect of solutal buoyancy force is in the same direction as the thermal buoyancy force. In such cases, the isothermal and isoconcentration tend toward the right wall. The streamline, isotherm and isoconcentration contours for Pr ¼ 0:7; Le ¼ 2:0 Da ¼ 103 and Ra ¼ 103 with non-uniformly heated and concentrated wall are displayed in Fig. 4. Comparing Figs. 3 and 4, it can be said that the patterns of streamlines are almost similar for uniformly and non-uniformly heated and concentrated cases. However, significant differences are observed in the patterns of isotherms and isoconcentrations for uniformly and non-uniformly heated and concentrated cases. The effects of Ra on streamline, isotherm and isoconcentration contours for Pr ¼ 0:7; Le ¼ 2:0 Da ¼ 103 and N ¼ 1 are displayed in Figs. 5 and 6 for uniform and non-uniform boundary conditions.

410

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413

Fig. 7. Contour plots for Pr ¼ 0:7; Ra ¼ 105 ; Le ¼ 2:0; Da ¼ 103 ; N ¼ 1 and (a) t ¼ 0:0279 (b) t ¼ 0:1395 (c) t ¼ 1:6015 (non-uniformly heated and concentrated).

Fig. 8. Local Nusselt number and local Sherwood number for uniformly heated and concentrated on the left vertical wall when Pr ¼ 0:7; Le ¼ 2:0; Da ¼ 103 and Ra ¼ 103 .

Fig. 5 shows the effect of Ra on the streamlines, isotherms and as well as on the isoconcentration for wide range of variations in Ra with uniform boundary conditions. We observe that there is no change in the values of stream functions, isotherms and isoconcen-

from 104 to 105 , the values of the stream function increase due to

trations values for Ra ¼ 103 and Ra ¼ 104 . However as Ra increases

concentrations for Ra ¼ 103 and Ra ¼ 104 . Fig. 6 shows the effect

convection. We also observe in from Fig. 5(c) that when Ra ¼ 105 , the contours of both the isotherms and isoconcentrations are located near the cold vertical wall compared to the isotherms and iso-

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413

411

Fig. 9. Local Nusselt number and local Sherwood number for uniformly heated and concentrated on the bottom wall when Pr ¼ 0:7; Le ¼ 2:0; Da ¼ 103 and Ra ¼ 103 .

Fig. 10. Local Nusselt number and local Sherwood number for non-uniformly heated and concentrated on the left vertical wall when Pr ¼ 0:7; Le ¼ 2:0; Da ¼ 103 and Ra ¼ 103 .

Fig. 11. Local Nusselt number and local Sherwood number for non-uniformly heated and concentrated on the bottom wall when Pr ¼ 0:7; Le ¼ 2:0; Da ¼ 103 and Ra ¼ 103 .

412

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413

of Ra on the streamlines, isotherms as well as on the isoconcentration for a range of variations in Ra with non-uniform boundary condition. We observe that the values of the stream function increase as Ra increases. Fig. 6(c) that when Ra ¼ 105 , the contours of both the isotherms and isoconcentrations are mainly concentrated near the cold vertical wall but when Ra decreases to 104 (see Fig. 6(b)), they are dispersed. When Ra is further decreased to 103 (see Fig. 6(a)), the contours are more dispersed towards the left vertical wall. Fig. 7 shows the changes in the streamline, isotherm and isoconcentration patterns with respect to time when the bottom and left walls are non-uniformly heated and concentrated when Ra ¼ 105 , Pr ¼ 0:7; Le ¼ 2:0; N ¼ 1 and Da ¼ 103 . The stream lines are very much concentrated near the left wall of the cavity when the non-dimensional time t ¼ 0:0279. Beyond this time, the stream lines are dispersed to the center of the cavity when non-dimensional time t ¼ 0:1395. Ultimately, the flow inside the cavity reaches a steady state and the stream lines are mostly concentrated near the right wall when t ¼ 1:6015. The isotherm pattern changes significantly with time. Initially when t ¼ 0:0279, the isotherm contours are concentrated near the left and lower walls. After that, as time progresses isotherms spread through out the cavity. The similar thing happens to the pattern of isoconcentrations also. Fig. 8 displays the variation of local Nusselt number and local Sherwood number at the heated and concentrated left vertical wall for various values of N due to uniform temperature and uniform concentration boundary conditions when Ra ¼ 103 ; Pr ¼ 0:7; Le ¼ 2:0 and Da ¼ 103 . It is interesting to note that in this figure, the local Nusselt and Sherwood numbers are zero at the bottom edge of left vertical wall. However, as the value of N increases, the local Nusselt number and the local Sherwood number decreases at the top edge of the left vertical wall. Fig. 9 shows the local Nusselt number and local Sherwood number at the heated and concentrated bottom wall for various values of N due to uniform temperature and uniform concentration boundary conditions when Ra ¼ 103 ; Pr ¼ 0:7; Le ¼ 2:0 and Da ¼ 103 . The local Nusselt and Sherwood numbers are zero at the left edge of bottom wall. Also, the local Nusselt and Sherwood numbers attain a maximum value near the right edge of the bottom wall. Figs. 10 and 11 show the distribution of the local Nusselt and Sherwood numbers at the bottom wall and left vertical wall when Ra ¼ 103 ; Pr ¼ 0:7; Le ¼ 2:0 and Da ¼ 103 due to non-uniformly heated and concentrated cavity walls. All these results show that heat transfer and mass transfer rates are high near the central region of the walls due to a discontinuity in temperature and concentration boundary conditions at these walls. 5. Conclusion The main objective of the current investigation was to study the influence of a non-uniform boundary conditions on doublediffusive natural convection in a cavity filled with porous media.  As, the buoyancy ratio increases the boundary layer thickness becomes thinner. The change in the flow structure for high buoyancy ratio has a significant influence on the concentration field, which builds up a vertical stratification of the enclosure.  When buoyancy ratio N ¼ 50, the effect of solutal buoyancy force is in the opposite direction of the thermal buoyancy force. The magnitude of the thermal buoyancy force is very small compared with the solutal buoyancy force due to which isotherm and isoconcentration tend to move to left wall whereas, when the buoyancy ratio N ¼ 1 both thermal and solutal

buoyancy forces are equal for which isothermal and isoconcentration tend to move to the right wall. For N ¼ 50, the effect of solutal buoyancy force is in the same direction of thermal buoyancy force. In such case the isothermal and isoconcentration tend toward the right wall.  The change of flow patterns with respect to time from t ¼ 0:0279 to t ¼ 1:6015 depicted and described here. Also, average Nusselt number and average Sherwood number at bottom wall and left vertical wall decrease as time increases from t ¼ 0:0279 to t ¼ 1:6015 for both uniform and non-uniform boundary conditions.  It is concluded that the heat and mass transfer rate are very high at the right edge of the bottom wall for uniformly heated and concentrated boundary conditions.  Finally, least square curves are fitted of average Nusselt number and average Sherwood number with Rayleigh number when Pr ¼ 0:7; Le ¼ 2:0; Da ¼ 103 ; Ra ¼ 103 and N ¼ 1 at the bottom and left vertical walls for uniform and non-uniform boundary conditions and the standard error are calculated. Conflict of interest I wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. Also I confirm that the manuscript has been read and approved by both authors and that there are no other persons who satisfied the criteria for authorship but are not listed. I further confirm that the order of authors listed in the manuscript has been approved by us. References [1] D.A. Nield, A. Bejan, Convection in Porous Media, fourth ed., Springer, New York, 2013. [2] D.B. Ingham, I. Pop (Eds.), Transport Phenomena in Porous Media III, Elsevier, Oxford, 2005. [3] K. Vafai (Ed.), Handbook of Porous Media, second ed., Taylor and Francis, New York, 2005. [4] I. Pop, D.B. Ingham, Convective Heat Transfer: Mathematical and Computational Modeling of Viscous Fluids and Porous Media, Pergamon, Oxford, 2001. [5] P. Vadasz (Ed.), Emerging Topics in Heat and Mass Transfer in Porous Media, Springer, New York, 2008. [6] T. Basak, S. Roy, A. Matta, I. Pop, Analysis of heatline for natural convection within porous trapezoidal enclosures: effect of uniform and non-uniform heating of bottom wall, Int. J. Heat Mass Transfer 53 (2010) 5947–5961. [7] M. Sathiyamoorthy, T. Basak, S. Roy, Non-Darcy buoyancy flow in a square cavity filled with porous medium for various temperature difference aspect ratios, J. Porous Media 14 (2011) 649–657. [8] Q.H. Deng, Jiemin, C. Mei, Y.M. Shen, Fluid, heat and contaminant transport structures of laminar double diffusive mixed convection in a two dimensional ventilated enclosure, Int. J. Heat Mass Transfer 47 (2004) 5257–5269. [9] S. Roy, T. Basak, Ch. Thirumalesha, C.M. Krishna, Finite element simulation on natural convection flow in a triangular enclosure due to uniform and nonuniform bottom heating, ASME Trans. J. Heat Transfer 130 (2008). Article Number 032501. [10] S.E. Alimi, J. Orfi, S.B. Nasrallah, Buoyancy effects on mixed convection heat and mass transfer in a duct with sudden expansions, Heat Mass Transfer 41 (6) (2005) 559–567. [11] N. Brown, F. Lai, Correlations for combined heat and mass transfer from an open cavity in a horizontal channel, Int. Commun. Heat Mass Transfer 32 (8) (2005) 1000–1008. [12] M.A. Teamah, Numerical simulation of double diffusive natural convection in rectangular enclosure in the presences of magnetic field and heat source, Int. J. Therm. Sci. 47 (2008) 237–248. [13] S. Saha, M.N. Hasan, I.A. Khan, Double diffusive mixed convection heat transfer inside a vented square cavity, Chem. Eng. Res. Bull. 13 (2009) 17–24. [14] W.J. Minkowycz, E.M. Sparrow, G.E. Schneider, R.H. Fletcher, Handbook of Numerical Heat Transfer, Wiley, New York, 1988. [15] S. Roy, T. Basak, Finite element analysis of natural convection flows in a square cavity with non-uniformly heated wall(s), Int. J. Eng. Sci. 43 (2005) 668–680. [16] M.A. Teamah, M.M. Khairat Dawood, W.M. El-Maghlany, Double diffusive natural convection in a square cavity with segmental heat sources, Eur. J. Sci. Res. 54 (2) (2011) 287–301.

S. Mondal, P. Sibanda / International Journal of Heat and Mass Transfer 85 (2015) 401–413 [17] M. Karimi-Fard, M.C. Charrier-Mojtabi, K. Vafai, Non-Darcian effects on double-diffusive convection within a porous medium, Numer. Heat Transfer Part A Appl. 31 (1997) 837–852. [18] P.M. Patil, S. Roy, A.J. Chamkha, Double diffusive mixed convection flow over a moving vertical plate in the presence of internal heat generation and a chemical reaction, Turk. J. Eng. Environ. Sci. 33 (2009) 193–205. [19] P. Nithiarasu, K.N. Seetharamu, T. Sundarararjan, Natural convective heat transfer in a fluid saturated variable porosity medium, Int. J. Heat Mass Transfer 40 (1997) 3955–3967. [20] T.R. Mahapatra, D. Pal, S. Mondal, Effects of buoyancy ratio on double-diffusive natural convection in a lid-driven cavity, Int. J. Heat Mass Transfer 57 (2013) 771–785. [21] F. Marcondes, J.M. Mederios, J.M. Gurgel, Numerical analysis of natural convection in cavities with variable porosity, Numer. Heat Transfer Part A 40 (2001) 403–420.

413

[22] T.R. Mahapatra, D. Pal, S. Mondal, Natural convection in a lid-driven square cavity filled with Darcy–Forchheimer porous medium in the presence of thermal radiation, Int. J. Nonlinear Sci. 11 (3) (2011) 366–379. [23] T.R. Mahapatra, D. Pal, S. Mondal, Mixed convection in a lid-driven square cavity filled with porous medium in the presence of thermal radiation and non-uniform heating, Int. J. Appl. Math. Mech. 9 (13) (2013) 23–51. [24] T.R. Mahapatra, G.C. Layek, M.K. Maiti, Unsteady laminar separated flow through constricted channel, Int. J. Non-linear Mech. 37 (2002) 171–186. [25] J.M. Medeiros, F. Marcondes, J.M. Gurgel, Natural convection in a porous cavity using the generalized model with uniform porosity, XV Brazilian Congress of Mechanical Engineering, in CD-ROM, Aguas de Lindoia, Sao Paulo, Brazil, 1999. [26] T.R. Mahapatra, D. Pal, S. Mondal, Non-Darcian natural convection and radiation in a cavity filled with fluid saturated porous Medium, Nonlinear Anal. Modell. Control 17 (2) (2012) 223–237.