08 JPM 1278revised2 PDF

Generated by Foxit PDF Creator © Foxit Software
http://www.foxitsoftware.com For evaluation only.

1-32
Numerical Analysis of Natural Convection in Porous Cavity with
Partial Convective Cooling Condition
W. Pakdee1 and P. Rattanadecho

Department of Mechanical Engineering, Faculty of Engineering,
Thammasat University, Rangsit Campus,
Klong Luang, Pathumtani 12120, Thailand
Tel: 011 662 5643001-5, Fax: 011 662 5643010
Abstract
The transient natural convection flow through a fluid-saturated porous medium in a

square enclosure with a partially cooling surface condition was investigated using
Brinkmann-extended Darcy model. Physical problem consists of a rectangular cavity filled
with porous medium. The cavity is insulated except the top wall that is partially exposed to an
outside ambient. The exposed surface allows convective transport through the porous
medium, generating a thermal stratification and flow circulations. The formulation of
differential equations is nondimensionalized and then solved numerically under appropriate
initial and boundary conditions using the finite difference method. The finite difference
equations handling the convection boundary condition of the open top surface are derived for
cooling condition. In addition to the negative density gradient in the direction of gravitation, a
lateral temperature gradient in the region close to the top wall induces the buoyancy force
under an unstable condition. The two-dimensional flow is characterized mainly by the
clockwise and anti-clockwise symmetrical vortices driven by the effect of buoyancy. The
directions of vortex rotation generated under the cooling condition are in the opposite
direction as compared to the heating condition. Unsteady effects of associated parameters
were examined. The modified Nusselt number (Nu) is systematically derived. This newly
developed form of Nu captures the heat transfer behaviors reasonably accurately. It was found
that the heat transfer coefficient, Rayleigh number, Darcy number as well as flow direction
strongly influenced characteristics of flow and heat transfer mechanisms.
Keywords: convective transport, numerical analysis, saturated porous media, convection
boundary condition, partial cooling
1
Author to whom correspondence should be addressed; E-mail address: pwatit@engr.tu.ac.th, and
wpele95@yahoo.com
2-32
Nomenclature
cp specific heat capacity [J/kgK]
Da Darcy number [-]

g gravitational constant [m/s 2 ]
H cavity length [m]
h convective heat transfer coefficient [W/m2K]
k thermal conductivity of the porous medium [W/mK]
p pressure [Pa]
Pr Prandtl number [-]
Ra Rayleigh number [-]
T temperature [C]
t time [s]
u, v velocity component [m/s]
x, y Cartesian coordinates
X, Y dimensionless Cartesian coordinates
W cavity width [m]
Greek symbols
κ permeability of porous medium[m 2 ]
α thermal diffusivity [m 2 /s]
β coefficient of thermal expansion [1/K]
ε porosity [-]
µ dynamic viscosity [Pa/s]
ν kinematics viscosity [m 2 /s]

ρf fluid density [kg/m 3 ]
τ dimensionless time
θ dimensionless temperature
ω vorticity [s-1]
ψ stream function [m2/s]
ς dimensionless vorticity
Ψ dimensionless stream function
3-32
Subscripts
∞ ambient condition
i initial condition and index for a number of points in x direction
j index for a number of points in y-direction
e effective
1. Introduction
The convective heating or cooling that causes heat and fluid flows inside cavity is
found in various applications including lakes and geothermal reservoirs, underground water
flow, solar collector etc. (Bergman et al., 1986). Associated industrial applications include
secondary and tertiary oil recovery, growth of crystals (Imberger and Hamblin, 1982), heating
and drying process (Stanish et al., 1986; Rattanadecho et al., 2001, 2002), electronic device
cooling, solidification of casting, sterilization etc. Natural or free convection in a porous
medium has been studied extensively. Cheng (1978) provides a comprehensive review of the
literature on free convection in fluid-saturated porous media with a focus on geothermal
systems. In the framework of porous media models, Darcy proposed the phenomenological
relation between the pressure drop across a saturated porous medium and the flow rate. The
Darcy model has been employed in the recent investigations. Bradean et al. (1997) assumed
Darcy’s law and used Boussinesq approximation to numerically simulate the free convection
flow in a porous media adjacent to vertical or horizontal flat surface. The surface is suddenly
heated and cooled sinusoidally along its length. The Darcy law with the Boussinesq
approximation was also employed by Bilgen and Mbaye (2001) to study the development of
Be’n-ard cell in fluid-saturated porous cavity whose lateral walls are cooled. It was found that
the existence of two convective solution branches is related to the Darcy-Rayleigh and Biot
numbers. Recently, a numerical study was conducted to solve the problem of thermosolutal
convection within a rectangular enclosure (Bera and Khalili, 2002). The results revealed that
anisotropy causes significant changes in Nusselt and Sherwood numbers. Many works of flow
4-32
in porous media, such as ones addressed above, have used the Darcy law. Although the Darcy
law is applicable to slow flows, it does not account for initial and boundary effects. In the
situation when the flow is strong, and solid boundary effect and viscous effect are not
negligible, these effects termed non-Darcy effects, become important (Khanafer and
Chamkha, 1998). Bera et al. (1998) considered double diffusive convection due to constant
heating and cooling on the two vertical walls, based on a non-Darcy model inclined
permeability tensor. Two distinguished modifications of Darcy’ law are the Brinkmann’s and
the Forchheimer’s extensions which treats the viscous stresses at the bounding walls and the
non-linear drag effect due to the solid matrix respectively (Nithiarasu et al., 1997). The
Darcy-Forchheimer- Brinkman model was used to represent the fluid transport within the
porous medium in the investigation of a convective flow through a channel (Marafie and
Vafai, 2001). In this work, the two-equation model was used to describe energy transport for
solid and fluid phase. The Brinkman-extended Darcy model has been considered in a
literature (Tong and Subramanian, 1985; Laurat and Prasad, 1987; Kim et al., 2001; Pakdee
and Rattanadecho, 2006). Darcy-Forchheimer model has been used in a number of published
works (Beckermann et al., 1985; Lauriat and Prasad, 1989; Basak et al., 2006). In the study of
effects of various thermal boundary conditions applied to saturated porous cavity, the
conduction dominant regime is within Da ≤ 10−5 . Nithiarasu et al. (1998) examined effects of
applied heat transfer coefficient on the cold wall of the cavity upon flow and heat transfer
inside a porous medium. The differences between the Darcy and non-Darcy flow regime are
clearly investigated for different Darcy, Rayleigh and Biot numbers and aspect ratio.
Variations in Darcy, Rayleigh and Biot numbers and aspect ratio significantly affect natural
flow convective pattern.
Natural convection flows with a variety of configurations were investigated for different
aspects. Oosthuizen and Patrick (1995) performed numerical studies of natural convection in
5-32
an inclined square enclosure with part of one wall heated to a uniform temperature and with
the opposite wall uniformly cooled to a lower temperature and with the remaining wall
portions. The enclosure is partially filled with a fluid and partly filled with a porous medium,
which is saturated with the same fluid. The main results considered were the mean heat
transfer rate across the enclosure. Nithiarasu et al. (1997) examined effects of variable
porosity on convective flow patterns inside a porous cavity. The flow is triggered by
sustaining a temperature gradient between isothermal lateral walls. The variation in porosity
significantly affects natural flow convective pattern. Khanafer and Chamkha (1998)
performed numerical study of mixed convection flow in a lid-driven cavity filled with a fluid-
saturated porous media. In this study, the influences of the Richardson number, Darcy number
and the Rayleigh number play an important role on mixed convection flow inside a square
cavity filled with a fluid-saturated porous media. Recently, Al-Amiri (2000) performed
numerical studies of momentum and energy transfer in a lid-driven cavity filled with a
saturated porous medium. In this study, the force convection is induced by sliding the top
constant-temperature wall. It was found that the increase in Darcy number induces flow
activities causing an increase in the fraction of energy transport by means of convection. With
similar description of the domain configuration, Khanafer and Vafai (2002) extended the
investigation to mass transport in the medium. The buoyancy effects that create the flow are
induced by both temperature and concentration gradients. It was concluded that the influences
of the Darcy number, Lewis number and buoyancy ratio on thermal and flow behaviors were
significant. Furthermore, the state of art regarding porous medium models has been
summarized in the recently published books (Nield and Bejan, 1999; Vafai, 2000; Pop and
Ingham, 2001; Basak et al., 2006).
Previous investigations have merely focused on momentum and energy transfer in cavity
filled with a saturated porous medium subjected to prescribed temperature and prescribed wall
6-32
heat flux conditions. However, only a very limited amount of numerical and experimental
work on momentum and energy transfer in a cavity filled with a saturated porous medium
subjected to heat transfer coefficient boundary condition at the exposed portion of the top wall
has been reported. Moreover, only very few published work is pertinent to partially heated or
cooled porous media although they are found in a number of applications such as in flush
mounted electrical heater or buildings (Desai et al, 1997; Al-Amiri, 2002; Oztop, 2007). The
very recent work of Oztop (2007) investigated natural convection in partially cooled and
inclined porous enclosures. His study presented the steady state results within the enclosure of
isothermal heated and cooled walls. In our study, the surface is partially cooled under the
convective boundary condition, allowing the surface temperature to change with time. The
convective cooling condition or so-called condition of the third kind is systematically derived.
While the focus of the present study is on the cooling effect, our recently published work
(Pakdee and Rattanadecho, 2006) studied the influence of partially heated surface on
thermal/flow behaviors. In this previous work, although the results were qualitatively
discussed in detail, no quantitative description of heat transfer in terms of Nusselt number
(Nu) was reported. Therefore, in order to gain better insights into the analysis, our present
study proposes a new formulation of Nu employed to analyze the heat transfer behaviors.
Moreover, to the best knowledge of the authors, no attention has been paid to transient
convection due to surface partial convective cooling.
In the present study, the quantitative study in terms of Nu is taken into account. The new
formulation of Nu is developed to correctly capture heat transfer behaviors. The study of heat
transfer due to cooling condition has been carried out for transient natural convective flow in
a fluid-saturated porous medium filled in a square cavity. In contrast to the heating condition,
the cooling condition changes the direction of the induced flows. The top surface is partially
open to the ambient, allowing the surface temperature to vary, depending on the influence of
7-32
convection heat transfer mechanism. Computed results are depicted using temperature, flow
distributions and heat transfer rates in terms of local and average Nusselt numbers. The
influences of associated parameters such as Rayleigh number and Darcy number on the flow
and thermal configurations were examined.
2. Problem Description
The computational domain, depicted in figure. 1 is a rectangular cavity of size W×H
filled with a fluid-saturated porous medium. Aspect ratio of unity (A=1) is used in the present
study. The domain boundary is insulated except the top wall, which is partially exposed to an
ambient air. The initial and boundary conditions corresponding to the problem are of the
following forms.
u = v = 0, T=Ti for t=0 (1)
u = v = 0 at x = 0, W 0 ≤ y ≤ H
 (2)
u = v = 0 at y = 0, H 0 ≤ x ≤ W
∂T 
= 0 at x = 0,W 0≤ y≤H 
∂x 
∂T 
= 0 at y = 0 0≤ x≤ W  (3)
∂y 
∂T 
= 0 at y = H 0 ≤ x ≤ L and W-L ≤ x ≤ W 
∂y 
The boundary condition at the exposed portion of the top wall is defined as
∂T
−k = h[T − T∞ ] at y = H L ≤ x ≤ W-L, (4)
∂y
where k and h are effective thermal conductivity and convection heat transfer coefficient. This
type of condition corresponds to the existence of convective heat transfer at the surface and is
obtained from the surface energy balance.
The porous medium is assumed to be homogeneous and thermally isotropic. The
saturated fluid within the medium is in a local thermodynamic equilibrium (LTE) with the
8-32
solid matrix (El-Refaee et al., 1998; Nield and Bejan, 1999; Al-Amiri, 2002). The validity
regime of local thermal equilibrium assumption has been established (Mohammad, 2000;
Marafie and Vafai, 2001). The porous porosity is uniform. The fluid flow is unsteady, laminar
and incompressible. The pressure work and viscous dissipation are all assumed negligible.
The thermophysical properties of the porous medium are taken to be constant. However, the
Boussinesq approximation takes into account of the effect of density variation on the
buoyancy force, in which the fluid density is assumed constant except in the buoancy term of
the equation of motion. Furthermore, the solid matrix is made of spherical particles, while the
porosity and permeability of the medium are assumed to be uniform throughout the
rectangular cavity. Using standard symbols, the governing equations describing the heat
transfer phenomenon are given by
∂u ∂v
+ =0 (5)
∂x ∂y
1 ∂u u ∂u v ∂u 1 ∂P υ  ∂ 2u ∂ 2u 
+ 2 + 2 =− +  + 
ε ∂t ε ∂x ε ∂y ερ f ∂x ε  ∂x 2 ∂y 2 
µu
− (6)
ρfκ
1 ∂v u ∂v v ∂v 1 ∂P υ  ∂ 2 v ∂ 2 v 
+ 2 + 2 =− +  + 
ε ∂t ε ∂x ε ∂y ερ f ∂y ε  ∂x 2 ∂y 2 
µv
+ g β (T − T∞ ) − (7)
ρfκ
∂T ∂T ∂T  ∂ 2T ∂ 2T 
σ +u +v =α  2 + 2  (8)
∂t ∂x ∂y  ∂x ∂y 
[ε ( ρ c p ) f + (1 − ε )( ρ c p ) s ]
σ= , (9)
(ρcp ) f
where κ is medium permeability, β is thermal expansion coefficient, α is effective thermal
diffusivity of the porous medium, μ and υ are viscosity and kinematic viscosity of the fluid
respectively. Symbols ε and ν denotes porosity of porous medium and fluid viscosity,
respectively. In the present study, the heat capacity ratio σ is taken to be unity since the
thermal properties of the solid matrix and the fluid are assumed identical (Bergman et al.,
9-32
1986; Khanafer and Vafai, 2002). The momentum equation consists of the Brinkmann term,
which describes viscous effects due to the presence of solid body (Brinkmann, 1947). This
form of momentum equation is known as Brinkmann-extended Darcy model. Lauriat and
Prasad (1987) employed the Brinkmann-extended Darcy formulation to investigate the
buoyancy effects on natural convection in a vertical enclosure. Although the viscous
boundary layer in the porous medium is very thin for most engineering applications, inclusion
of this term is essential for heat transfer calculations (Al-Amiri, 2000). However, the inertial
effect was neglected, as the natural convection flow was studied (Basak et al., 2006).
The variables are transformed into the dimensionless quantities defined as,
x y tα uH vH 
X= ,Y = , τ = 2 , U = ,V=
H H H α α 
, (10)
ωH 2 ψ T − Tl 
ς= , Ψ= ,θ =
α α Th − Tl 
where ω and ψ represent dimensional vorticity and stream function, respectively. Symbol α
denotes thermal diffusivity. Temperatures Tl and Th change their values according to the
problem type. In the heating case, Tl is initial temperature of a medium, and Th is an ambient
temperature. In the other case of cooling, Th is set to be an initial temperature of the medium,
while Tl is an ambient temperature instead. The governing equations are transformed into a
vorticity −stream function formulation. Thus the dimensionless form of the governing
equations can be written as
∂ 2Ψ ∂ 2Ψ
+ = −ς (11)
∂X 2 ∂Y 2
∂ς ∂ς ∂ς  ∂ 2ς ∂ 2ς   ∂θ  ε Pr
2
ε +U +V = ε Pr  2 + 2  + ε 2 Ra Pr   − ς (12)
∂τ ∂X ∂Y  ∂X ∂Y   ∂X  Da
∂θ ∂θ ∂θ  ∂ 2θ ∂ 2θ 
σ +U +V = α  2 + 2 , (13)
∂τ ∂X ∂Y  ∂X ∂Y 
∂Ψ ∂Ψ
U= , V =− (14)
∂Y ∂X
10-32
where the Darcy number, Da is defined as κ / H 2 , and Pr= ν/α is Prandtl number, where α =
ke/(ρcp)f is the thermal diffusivity. The Rayleigh number Ra, which gives the relative
magnitude of buoyancy and viscous forces, is defined as Ra = g β (Ti − T∞ )H 3 /(υα ) .
3. Numerical Procedure
The thermal properties of the porous medium are taken to be constant. Specific heat ratio
of unity is assumed. The effective thermal conductivity of the porous medium considered is
10 W/m·K.
In the present study, the iterative finite difference method was used to solve the transient
dimensionless governing equations (Eqs. (10)−(12)) subject to their corresponding initial and
boundary conditions given by Eqs. (1)−(4). Central−difference formulae were used for all
spatial derivatives. The transient transport equations, Eqs (12)−(13), were solved explicitly.
Successive over relaxation method (SOR) was utilized to solve for the flow kinematics
relation given by Eq. (11). The velocity components, U and V, were computed according to
Eq. (14). Approximation of convective terms is based on a second−order upwind finite
differencing scheme, which correctly represents the directional influence of a disturbance. A
uniform grid resolution of 61 × 61 was found to be sufficient for all smooth computations and
computational time required in achieving steady-state conditions. Finer grids did not provide
a noticeable change in the computed results.
3.1 Convective cooling boundary condition
The finite difference form of boundary condition at the open part of the top surface is
systematically derived, based on energy conservation principle. The boundary values of
dimensionless temperature of a node i, j θi,j in the heating case are expressed as

11-32
h∆y
2θij −1 + θ i −1 j + θi +1 j + 2
θij = k , (15)
h
2( ∆Y + 2)
k
where ∆Y is the mesh size in y-direction.
In the different case of cooling phenomenon, the expression is given by
2θij −1 + θi −1 j + θi +1 j
θij = (16)
h
2( ∆Y + 2)
k
It can be noticed that both the equations (15) and (16) are independent of an ambient
temperature T ∞ as it has been eliminated during the derivation. This feature is attractive
since the solutions can be obtained regardless of a value of T ∞ .
3.1 Corrected formulation of Nusselt number
The local Nusselt number (Nu) at the cooled horizontal surface is used as a tool to
determine the ratio of convection heat transfer to conduction heat transfer within the porous
enclosure. The accurate derivation of Nu is extremely important from the standpoint of
determining the rate of heat transfer occurring at a surface. Based on the concept of energy
balance at the surface for the cooling case,
dT
−k = h(TH − T∞ ) , (17)
dy y= H
where H is indicated in figure 1, and with the definition of Nu,
hH H dT
Nu = =− (18)
k (TH − T∞ ) dy Y =H
In terms of the dimensionless quantities θ and Y defined in the preceding equation (10), Nu
will take the form,
1 dθ
Nu = − , (19)
θ H dY Y =1
12-32
where θ H is the dimensionless temperature at the top surface.
The new formulation of Nu in the present work has not yet been found in the literature. This
modified form of Nu takes into account of temperature variation at the cooled surface. The
average Nusselt Number, Nu is computed according to
W −L
Nu( x)dx
Nu = ∫ , (20)
L l
where l is the length of the gap at the top wall.
In order to verify the accuracy of the present numerical study, the results obtained by
the present numerical model were validated against the Benchmark solutions for natural
convection in a cubic cavity (Wakashima and Saitoh, 2004). The comparisons tabulated in
table 1 reveal an excellent agreement within 1.5 percent difference. Also present computed
results were compared with those obtained by Aydin (2000) for a free convection flow in a
cavity, with side-heated isothermal wall, filled with pure air (Pr = 0.7) for Rayleigh number of
104. It was found that the solutions have good agreement with the previously published work.
The results of selected tests are given in table 2 that shows a good agreement of the maximum
value of the stream function and the maximum values of the horizontal and vertical velocity
components between the present solution and that of Aydin. Moreover, the results from the
present numerical model were compared with the solution of Nithiarasu et al. (1997) in the
presence of porous medium for additional source of confidence, as shown in figure 2 for
streamlines and isotherms for which the compared contours have the same range of contour
levels. The values of Ra = 104, Da = 0.01 and ε = 0.6 were chosen. Table 3 clearly shows a
good agreement of the maximum values of the stream function and vertical velocity
component between the present solution and that of Nithiarasu et al (1997). All of these
favorable comparisons lend confidence in the accuracy of the present numerical model.
13-32
4. Results and Discussion
The following discussions include the numerical results from the present study, which
focuses on transient flow and thermal behaviors. Initial values of θ for an entire domain are
set 1, based on equation (10) as the ambient temperature is lower than temperature of the
medium in cavity. The investigations were conducted for a range of controlling parameters,
which are Darcy number (Da) Rayleigh number (Ra) and convective heat transfer coefficient
(h). The uniform porosity ε of 0.8 and unity aspect ratio (A=1) were considered throughout in
the present study. In order to assess global effects of these parameters, the streamlines and
isotherm distributions inside the entire cavity are presented. All the figures have the same
range of contour levels to facilitate direct comparisons.
The resulting computational fields were extracted at the time adequately long to ensure
sufficient energy transferred throughout the domain. Figure 3 displays instantaneous images
of the contour plots during the thermal and flow evolution. The Rayleigh number of 5x104,
Da = 0.1, Pr = 1.0, h = 60 w/m2K, and ε =0.8 are considered. The two columns represent
contours of temperature and stream function respectively from left to right. With the same
contour levels, comparisons can be made directly. The four snapshots from top to bottom in
each column are results taken at the dimensionless times τ = 0.013, 0.088, 0.168, and 0.245.
The vertical temperature stratification is observed. The streamline contours exhibit circulation
patterns, which are characterized by the two symmetrical vortices. The fluid flows as it is
driven by the effect of buoyancy. This effect is distributed from the top wall of cavity where
the fluid is cooled through the partially open surface, causing lower temperature near the top
boundary. The existence of the non-uniform temperature along the top surface, and a decrease
of density in the direction of gravitational force lead to an unstable condition. Thus the
buoyancy effect is associated with the lateral temperature gradients at locations near the top
surface. High temperature portions of fluid become lighter than the lower temperature
14-32
portions at the middle where the wall is open. Theses light portions from two sides then
expand laterally towards the center, compressing the lower temperature portions, which are
heavier. As a result, the downward flows along the vertical centerline are originated, while the
lighter fluid will rise, cooling as it moves. Consequently, the circulation flow pattern is
generated. The clockwise and counter-clockwise circulations are located respectively on the
left side and right side within the enclosure. The circulations get larger and expand downward
with time. An increase in strength of the vortices develops fast during early simulation times,
and its maximum magnitude reaches 6.0. Subsequently the vortices are weakened. Similarly,
temperature distribution progressively evolves relatively fast in the early times. Slow
evolution is observed after that. This result corresponds well with the decrease in strength of
flow circulations.
The resulting computational fields of the heating scenario were demonstrated in figure 4.
For a purpose of comparison the parameter set remains unchanged. Similarly, the two
columns represent temperature and stream function taken at the dimensionless times τ =
0.013, 0.088, and 0.168. The vertical temperature stratification is observed. The streamline
contours exhibit circulation patterns, which are characterized by the two symmetrical vortices.
The fluid flows as it is driven by the effect of buoyancy. This effect is distributed from the top
wall of cavity where the fluid is heated through the partially open area. Unlike the cooling
case, in which a presence of negative density mainly causes an unstable condition, in the
heating case the lateral density gradient near the top surface is the only cause to the unstable
condition that actually leads to the buoyancy force. This reason explains why the heated
circulations are weaker than the cooled circulations presented earlier. Heated portions of the
fluid become lighter than the rest of fluid, and are expanded laterally away from the center to
the sides then flow down along the two vertical walls, leading to the counter-clockwise and
clockwise flow circulations. These results suggest that the buoyancy forces are able to
15-32
overcome the retarding influence of viscous forces. It should be noted that directions of
circulations are opposite to those under cooling condition. An increase in strength of the
vortices develops fast during early simulation times, and its maximum magnitude reaches
0.25, which is considerably small. Therefore, profiles of temperature contours look similar to
those for a stationary fluid, in which the heat transfer is caused by conduction. Similarly,
temperature distribution progressively evolves relatively fast in the early times. This result
corresponds to the decrease in strength of flow circulations. In the remaining area, the fluid is
nearly stagnant suggesting that conduction is dominant due to minimal flow activities. This is
because of prevailing viscous effects. It is evident from figures 3 and 4 that the cooling case
provides a considerably faster thermal evolution thereby greater convection rate. Furthermore,
heat transfer in the vertical direction is much greater than that in the span wise direction. The
reader is directed to (Pakdee and Rattanadecho, 2006) for more detailed discussions of
heating configuration.
Figure 5 shows the roles of Rayleigh number on heat transfer mechanism. The
computed data was extracted at τ = 0.155. Various Rayleigh numbers (Ra = 5×103, 104, 5×104
and 105) are examined whereas the Darcy number of 0.1, porosity of 0.8, and h of 60 w/m2 K
are fixed. The Rayleigh number provides the ratio of buoyancy forces to change in viscous
forces. As Rayleigh number increases, the buoyancy-driven circulations inside the enclosure
become stronger as seen from greater magnitudes of stream function. For large Ra (Ra =
5×104 and 105), contour lines of temperature penetrate faster relative to the low Ra case
especially near the central locations. The result is more pronounced for larger Ra. This
incident results from strong flow in the downward direction around the central domain. The
downward flows assist heat to transfer towards the bottom of the enclosure. In contrast, near
the vertical walls where the upward flows are present, the thermal propagation is hindered.
Effects of the Darcy number on the fluid flow and temperature inside the rectangular
16-32
cavity are depicted in figure 6. The contour of isotherms and streamlines at τ = 0.155 are
plotted for different Darcy numbers while ε , Pr and h are kept at 0.8, 1.0 and 60 w/m2K
respectively. Relatively high Ra of 5×104 is chosen. The Darcy number, which is directly
proportional to the permeability of the porous medium, was set to 0.0001, 0.001 and 0.1. The
case in which the porous medium is absent corresponds to infinite Darcy number. The
presence of a porous medium within rectangular enclosure results in a force opposite to the
flow direction which tends to resist the flow which corresponds to suppress in the thermal
currents of the flow as compared to a medium with no porous (infinite Darcy number). It is
evident that the increase in Da enhances the streamline intensities thereby assisting downward
flow penetration, which causes the streamline lines, i.e., two symmetrical vortices to stretch
further away from the top surface. This results in expanding the region for which the
convection significantly influences an overall heat transfer process. Further, the evolution
results reveal faster rate of vertical temperature distribution than lateral rate. The results are
consistent with the thermal behaviors observed in figure 5 for the same reasoning, which
confirms how a flow direction impacts the convection heat transfer. On the other hand, as the
Darcy number decreases, the flow circulations as well as thermal penetration are
progressively retarded due to the reduced permeability of the medium. Figure 6d (Da =
0.0001) indicates that as Darcy number approaches zero, the two circulations confined within
the top domain appear very weak. In the remaining area, the fluid is nearly stagnant with very
small temperature gradient suggesting that conduction is dominant due to minimal flow
activities.
Figure 7 presents how the average Nusselt number changes with time for a variety of
Rayleigh numbers. The local Nu at the open portion on the top boundary is computed
according to equation (19). The average Nusselt number Nu is then obtained based on
equation (20). Initially, the value of Nu decreases rapidly for all cases of Rayleigh numbers,
17-32
clearly due to the fast reduction of temperature gradients. In the case of low Rayleigh number
of 2x104, Nu progressively decreases with time. While for higher Ra (5x104, 8x104), Nu
values become greater and reach peak values after some time. Further increasing Ra (8x104),
higher maximum Nu is reached more quickly due to greater flow intensities. At late
simulation times when stable state is approached, the values of Nu continually decrease and
essentially level off at late times, thereby diminishing heat transfer by means of heat
convection. It can be expected that Nu will continue to decrease with time as the steady state
is reached.
To gain insights into the observation made, the local values of the corresponding
thermal and flow behaviors were traced for Ra of 8x104. The data are extracted and depicted
in figure 8 at τ of 0.02, 0.1 and 0.16. The streamlines and isotherms are illustrated in figure
8a-c at τ = 0.02, 0.1 and 0.16 respectively. At τ = 0.02, the averaged Nu is small due to
minimal flow activities. Then Nu gets higher as the flows gets stronger, which can be seen in
figure 8(b) at τ = 0.1. The effect of the rigorous flows overcomes the continual reduction of
temperature gradient, resulting in the increase in Nu . At the subsequent times, the viscous
effect increasingly weakens the flows as shown in figure 8(c). As a result, the reduction of
temperature gradient prevails, causing Nu to decrease. These results correspond well with the
variation with time of the averaged Nu, depicted in figure 7. The results confirm the
validation of the proposed formulation of Nu.
To better understand the effects of Darcy number on the heat transfer behavior, variations
of Nu with time for different Darcy number are shown in figure 9. The resulting plots show
an interesting evidence of similar variations of Nu on Da and those on Ra, which was
observed previously in figure 7. Average Nu correlates with Ra in a way similar to correlation
of Nu with Da. Further increasing values of Da (0.05 and 0.1) cause larger Nu variations.
18-32
Locations of the peak values are altered relative to Da value. A peak of profile is reached
more quickly for higher Da. Greater Da gives higher Nu , suggesting that the higher overall
heat transfer rate is due to more energetic vortices. However, Nu substantially reduces at late
times.
5. Conclusions
Numerical simulations of natural convection flow through a fluid-saturated porous
medium in a rectangular cavity due to cooling convection at top surface were performed.
Transient effects of associated controlling parameters were examined. The two-dimensional
flow is characterized mainly by two symmetrical eddies that are initiated by the presence of
buoyancy effect. In the cooling case, the buoyancy effect is associated not only with the
lateral temperature gradient at locations near the top surface, but also with the condition that
the density gradient is negative in the direction of gravitational force. On the other hand, the
buoyancy force is induced solely by the lateral temperature gradient in the heating case. The
cooling and heating flow directions are opposite. Cooling flows are much stronger due to
greater buoyancy effects, indicating higher overall convection rate. The heat transfer
mechanism is analyzed using the newly derived formulation of Nu. Heat transfer rate is faster
around vertical symmetric line relative to the near-wall regions. Large values of Rayleigh
number increase streamline intensities, thus enhancing the downward flow penetration. The
temperature stratification penetrates deeper toward the bottom wall, and temperature range
within the domain is extended. Therefore it enlarges the region where convection mode is
significant. Small values of Darcy number hinder the flow circulations. Therefore the heat
transfer by convection is considerably suppressed. Furthermore, the new formulation of Nu
captures the heat transfer behaviors reasonably correctly. Interestingly, the dependences of Nu
on Da and on Ra are found to have the same trends.

19-32
Acknowledgments
This research was supported by the Thailand Research Fund (TRF) 2008 under contract
number: MRG5180238
References
[1] Al-Amiri, A.M., Analysis of Momentum and Energy Transfer in a Lid-Driven Cavity
Filled with a Porous Medium. International Journal of Heat and Mass Transfer, 43, pp.3513-
3527, 2000.
[2] Al-Amiri, A.M., Natural convection in porous enclosures: The application of the two-
energy equation model, Numerical Heat Transfer Part A, 41, pp.817-834, 2002.
[3] Aydin, O., Determination of Optimum Air-Layer Thickness in Double-Pane Windows.

Energy and Buildings, 32, pp.303-308, 2000.
[4] Basak, T., Roy, S., Paul, T., and Pop, I., Natural convection in a square cavity filled with
porous medium: Effects of various thermal boundary conditions. International Journal of Heat
and Mass Transfer, 49, pp.1430-1441, 2006.
[5] Beckermann, C., Viskanta, R., and Ramadhyani, S., A numerical study of non-Darcian
natural convection in a vertical enclosure filled with a porous medium. Numerical Heat
Transfer, 10, pp.557-570, 1986.
[6] Bera, P., Eswaran, V., and Singh, P., Numerical study of heat and mass transfer in an
anisotropic porous enclosure due to constant heating and cooling. Numerical Heat Transfer
Part A, 34, pp.887-905, 1998.
[7] Bera, P. and Khalili, A., Double-diffusive natural convection in an anisotropic porous
cavity with opposing buoyancy forces: multi-solutions and oscillations. International Journal
of Heat and Mass Transfer, 45, pp.3205-3217, 2002.
[8] Bergman, T.L., Incropera, F.P., and Viskanta, R., Correlation of Mixed Layers Growth in
Double-Diffusive, Salt-Stratified System Heated from Below. Journal of Heat Transfer, 108,
pp.206-211, 1986.
[9] Bilgen, E. and Mbaye, M., Be’nard cells in fluid-saturated pourous enclosures with lateral
cooling. International Journal of Heat and Fluid Flow, 22, pp.561-570, 2001.
[10] Brinkmann, H.C., On the permeability of media consisting of closely packed porous
particles. Applied Sciences Research, 1, pp.81-86, 1947.
[11] Bradean, R., Ingham, D.B., Heggs P.J., and Pop, I., The unsteady penetration of free
convection flows caused by heating and cooling flat surfaces in porous media. International
Journal of Heat and Mass Transfer, 40, pp.665-687, 1997.
[12] Cheng, P., Heat transfer in geothermal systems. Advances in Heat Transfer, 4, pp.1-105,
1978.
20-32
[13] Desai, C.P., Vafai, K., and Keyhani, M., On the natural-convection in a cavity with a
cooled top wall and multiple protruding heaters, J. Electronic Packaging, 117, pp.34-45, 1995.
[14] El-Refaee, M. M., Elsayed, M.M., Al-Najem, N.M., and Noor, A.A., Natural convection
in partially cooled tilted cavities, Int. J. Numerical Methods, 28, pp.477-499, 1998.
[15] Imberger, J. and Hamblin, P.F., Dynamics of Lakes, Reservoirs, and Cooling Ponds.
Annual Review of Fluid Mechanics, 14, pp.153-187, 1982.
[16] Khanafer, K.M., Chamkha, A.J., (1998). Mixed Convection Flow in a Lid-Driven
Enclosure Filled with a Fluid-Saturated Porous Medium. International Journal of Heat and
Mass Transfer 42, 2465-2481.
[17] Khanafer, K. and Vafai, K., Double-Diffusive Mixed Convection in a Lid-Driven

Enclosure Filled with a Fluid-Saturated Porous Medium. Numerical Heat Transfer Part A, 42,
pp.465-486, 2002.
[18] Kim, S.J., Kim, D., and Lee, D.Y., On the local thermal equilibrium in microchannel heat
sinks, Int. J. Heat Mass Transfer, 43, pp.1735-1748, 2000.
[19] Lauriat, G., and Prasad, V., Natural convection in a vertical porous cavity: a numerical
study for Brinkmann-extended Darcy formulation. ASME Journal of Heat Transfer, 109,
pp.295-320, 1987.
[20] Lauriat, G. and Prasad, V., Non-Darcian effects on natural convection in a vertical
porous enclosure. International Journal of Heat and Mass Transfer, 32, pp.2135-2148, 1989.
[21] Marafie, A and Vafai, K., Analysis of non-Darcian effects on temperature differentials in
porous media, Int. J. Heat Mass Transfer, 44, pp.4401-4411, 2001.
[22] Mohammad, A.A., Nonequilibrium natural convection in a differentially heated cavity

filled with a saturated porous matrix. Journal of Heat Transfer, 122, pp.380-384, 2000.
[23] Nield, D.A., Bejan, A., Convection in Porous Media, Springer, New York, 1999.
[24] Nithiarasu, P., Seetharamu, K.N., and Sundararajan, T., Natural convective heat transfer
in a Fluid Saturated variable porosity medium, Int. J. of Heat Mass Transfer, 40, pp.3955-
3967, 1997.
[25] Nithiarasu, P., Seetharamu, K.N., and Sundararajan, T., Numerical Investigation of
Buoyancy Driven Flow in a Fluid Saturated non-Darcian Porous Medium. International
Journal of Heat and Mass Transfer, 42, pp.1205-1215, 1998.
[26] Oosthuizen, P.H. and Patrick, H., Natural Convection in an Inclined Square Enclosure
Partly Filled with a Porous Medium and with a Partially Heated Wall, American Society of
Mechanical Engineers. Heat Transfer Division, (Publication) HTD 302, pp.29-42, 1995.
[27] Oztop, H.F., Natural convection in partially cooled and inclined porous rectangular
enclosure. International Journal of Thermal Sciences, 46, pp.149-156, 2007.
21-32
[28] Pakdee, W. and Rattanadecho, P., Unsteady effects on natural convective heat transfer
through porous media in cavity due to top surface partial convection Applied Thermal
Engineering. 26, pp.2316-2326, 2006.
[29] Pop, I. and Ingham, D.B., Convective Heat Transfer, Mathematical and Computational
Modeling of Viscous Fluids and Porous Media, Pergamon, Oxford, 2001.
[30] Ratanadecho, P., Aoki, K., and Akahori, M., A Numerical and Experimental
Investigation of the Modeling of Microwave Drying Using a Rectangular Wave Guide.
Journal of Drying Technology, 19, pp.2209-2234, 2001.
[31] Ratanadecho, P., Aoki, K., and Akahori, M., Influence of Irradiation Time, Particle Sizes
and Initial Moisture Content During Microwave Drying of Multi-Layered Capillary Porous
Materials. ASME Journal of Heat Transfer, 124, pp.151-161, 2002.
[32] Wakashima, S. and Saitoh, T.S., Benchmark solutions for natural convection in a cubic
cavity using the high-order time-space method, Int. J. Heat Mass Transfer, Vol. 47, pp.853-
864, 2004.
[32] Stanish, M.A., Schager, G.S., and Kayihan, F., A., Mathematical Model of Drying for
Hygroscopic Porous Media. AIChE Journal, 32, pp.1301-1311, 1986.
[33] Tong, T.W. and Subramanian, E., A boundary-layer analysis for natural convection in
vertical porous enclosures-use of the Brinkman-extended Darcy model, ASME Journal of
Heat Transfer, 28, pp.563-571, 1985.
[34] Vafai, K., (2000). Handbook of Porous Media, Marcel Dekker, New York.
22-32
TABLE
Table 1. Comparison of the results obtained in the present study with those of the benchmark
solutions for natural convection of air (Wakashima and Saitoh, 2004).
Ra ω center Difference Umax Difference Vmax Difference

(%) (%) (%)
104 present 1.111 0.82 0.202 1.51 0.220 0.2
previous work 1.102 0.199 0.222
105 present 0.262 1.55 0.144 1.41 0.249 1.22

previous work 0.258 0.142 0.246
Table 2. Comparison of the results obtained in the present study with those of Aydin (2000).
Present Published Difference

work work (%)
ψ max 5.070 5.087 0.33

Umax 16.300 16.225 0.46
Vmax 19.730 19.645 0.43
Table 3. Comparison of the results obtained in the present study with those of Nithiarasu et al.
(1997). (Da=0.01, Ra = 104, porosity = 0.6)
Present Published Difference

work work (%)
ψ max 2.53 2.56 1.17
Vmax 9.49 9.34 1.60

23-32
LIST of FIGURE CAPTIONS
Figure 1. Schematic representation of the computational domain.
Figure 2. Test results for validation purpose: a) Nithiarasu et al. (1997): Non-Darcian model
(including inertial and boundary effect) b) present simulation: Brinkman-extended Darcy
model, which accounts for viscous effects.
Figure 3. Sequential files with the cooling boundary for contours of temperature and
streamlines at times τ = (a) 0.013, (b) 0.088, (c) 0.168, and (d) 0.245. (Ra = 5x104, Da = 0.1,
Pr =1.0, ε = 0.8, and h = 60 W/m2K)
Figure 4. Sequential files with the heating boundary for contours of temperature and
streamlines at times τ = (a) 0.013, (b) 0.088, and (c) 0.168. (Ra = 5x104, Da = 0.1, Pr =1.0, ε
= 0.8, and h = 60 W/m2 K)
Figure 5. Contours of temperature and streamlines for the cooling case (a) Ra = 5x103 (b) Ra
= 104 (c) Ra = 5x104 (d) Ra = 105. (Da = 0.1, h = 60 W/m2 K, Pr = 1.0, and ε = 0.8)
Figure 6. Contours of temperature and streamlines for the cooling case (a) Da = infinity (b)
Da = 0.01 (c) Da = 0.001 (d) Da = 0.0001. (Ra = 5x104, h = 60 W/m2K, Pr = 1.0, and ε =
0.8)
Figure 7. Variations of the average Nusselt number with time for different Rayleigh numbers.
(Da = 0.01, h = 60 W/m2K, Pr = 1.0, and ε = 0.8)
Figure 8. (a)-(c) temperature contours overlaid by velocity vectors at τ = 0.02, 0.1 and 0.16
respectively. Data is taken from that of figure 7 for Ra = 8x104.
Figure 9. Variations of the average Nusselt number with time for different Darcy numbers.
(Ra = 5x104, h = 60 W/m2K, Pr = 1.0, and ε = 0.8)

24-32
FIGURES
Figure 1
25-32
(a)
(b)
Figure 2
26-32
Figure 3
27-32
Figure 4
28-32
Figure 5
29-32
Figure 6
30-32
Figure 7
31-32
Figure 8
32-32
Figure 9

08 JPM 1278revised2 PDF

Uploaded by

Copyright:

Available Formats

08 JPM 1278revised2 PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

08 JPM 1278revised2 PDF

Uploaded by

Copyright:

Available Formats

Generated by Foxit PDF Creator © Foxit Software

http://www.foxitsoftware.com For evaluation only.

Numerical Analysis of Natural Convection in Porous Cavity with

Partial Convective Cooling Condition

W. Pakdee1 and P. Rattanadecho

The transient natural convection flow through a fluid-saturated porous medium in a

boundary condition, partial cooling

cp specific heat capacity [J/kgK]

Da Darcy number [-]

ν kinematics viscosity [m 2 /s]

cooling, solidification of casting, sterilization etc. Natural or free convection in a porous

literature on free convection in fluid-saturated porous media with a focus on geothermal

flow convective pattern.

Ingham, 2001; Basak et al., 2006).

discussed in detail, no quantitative description of heat transfer in terms of Nusselt number

convection due to surface partial convective cooling.

and thermal configurations were examined.

The computational domain, depicted in figure. 1 is a rectangular cavity of size W×H

u = v = 0, T=Ti for t=0 (1)

obtained from the surface energy balance.

The porous medium is assumed to be homogeneous and thermally isotropic. The

transfer phenomenon are given by

where κ is medium permeability, β is thermal expansion coefficient, α is effective thermal

form of momentum equation is known as Brinkmann-extended Darcy model. Lauriat and

Prasad (1987) employed the Brinkmann-extended Darcy formulation to investigate the

buoyancy effects on natural convection in a vertical enclosure. Although the viscous

equations can be written as

magnitude of buoyancy and viscous forces, is defined as Ra = g β (Ti − T∞ )H 3 /(υα ) .

Eq. (14). Approximation of convective terms is based on a second−order upwind finite

differencing scheme, which correctly represents the directional influence of a disturbance. A

a noticeable change in the computed results.

3.1 Convective cooling boundary condition

systematically derived, based on energy conservation principle. The boundary values of

dimensionless temperature of a node i, j θi,j in the heating case are expressed as

where ∆Y is the mesh size in y-direction.

In the different case of cooling phenomenon, the expression is given by

since the solutions can be obtained regardless of a value of T ∞ .

3.1 Corrected formulation of Nusselt number

enclosure. The accurate derivation of Nu is extremely important from the standpoint of

balance at the surface for the cooling case,

where H is indicated in figure 1, and with the definition of Nu,

will take the form,

where θ H is the dimensionless temperature at the top surface.

average Nusselt Number, Nu is computed according to

4. Results and Discussion

range of contour levels to facilitate direct comparisons.

validation of the proposed formulation of Nu.

an interesting evidence of similar variations of Nu on Da and those on Ra, which was

observed previously in figure 7. Average Nu correlates with Ra in a way similar to correlation

Numerical simulations of natural convection flow through a fluid-saturated porous

Transient effects of associated controlling parameters were examined. The two-dimensional

transfer by convection is considerably suppressed. Furthermore, the new formulation of Nu

on Da and on Ra are found to have the same trends.

[3] Aydin, O., Determination of Optimum Air-Layer Thickness in Double-Pane Windows.

[17] Khanafer, K. and Vafai, K., Double-Diffusive Mixed Convection in a Lid-Driven

[22] Mohammad, A.A., Nonequilibrium natural convection in a differentially heated cavity

Ra ω center Difference Umax Difference Vmax Difference

105 present 0.262 1.55 0.144 1.41 0.249 1.22

Present Published Difference