...

Madera2011-ConvectiveHT-PorousChannel.pdf

by user

on
Category: Documents
2

views

Report

Comments

Transcript

Madera2011-ConvectiveHT-PorousChannel.pdf
International Journal of Thermal Sciences 50 (2011) 1355e1368
Contents lists available at ScienceDirect
International Journal of Thermal Sciences
journal homepage: www.elsevier.com/locate/ijts
Convective heat transfer in a channel partially filled with a porous medium
Carlos G. Aguilar-Madera a, *, Francisco J. Valdés-Parada a, Benoît Goyeau b, J. Alberto Ochoa-Tapia a, *
a
b
División de Ciencias Básicas e Ingeniería, Universidad Autónoma Metropolitana-Iztapalapa, 6 Apartado Postal 55-534, 09340 Mexico D.F., Mexico
Laboratoire EM2C, UPR-CNRS 288, Ecole Centrale Paris, Grande Voie des Vignes, 92295 Châtenay-Malabry, Cedex, France
a r t i c l e i n f o
a b s t r a c t
Article history:
Received 20 October 2010
Received in revised form
8 March 2011
Accepted 8 March 2011
Available online 14 April 2011
In this work we solve effective-medium equations for modeling momentum and heat transfer in
a parallel-plate channel partially filled with a porous insert. In order to avoid specifying the boundary
conditions at the fluideporous boundary, we solve equations involving position-dependent coefficients
(i.e., a one-domain approach). The solution of the momentum-transport problem is carried out using
implicit integral equation formulations based on Green’s functions, whereas the heat transfer problem
was solved numerically using the finite element method. The simulations were performed in terms of
several values of the porosity, the Péclet number, the size of porous insert and the thermal conductivities
ratio. In agreement with previous works, it was found that the thermal performance is improved by
either increasing the size of the porous insert or by favoring mixing inside the channel. A drawback of
this approach is the high computational demand associated to modeling transport in the vicinity of the
porous medium and the fluid. In this way, the extents and limitations about the use of a one-domain
formulation are exposed in a practical application. The results from this work should serve as motivation
for more experimental and theoretical research; in particular, the derivation and application of jump
boundary conditions.
2011 Elsevier Masson SAS. All rights reserved.
Keywords:
One-domain approach
Local thermal equilibrium model
Porous insert
Parallel-plate channel
1. Introduction
The study of porous inserts for improving heat transfer in
modified equipments has recently acquired substantial interest
[1e6]. This is triggered mainly due to the wide range of applications
in engineering systems. Some examples are: drying processes [7],
heterogeneous reactors [8], filtering, cooling and heating processes
[9], geothermal energy management [10], among others.
The improvement of heat transport by using porous substrates
in different types of heat exchangers has been emphasized in
several works in the literature. Alkam and Al-Nimr studied the
problem of transient developing forced convective flow in
concentric pipes [4] and in circular channels [1] partially filled with
porous inserts. They showed the role played by the thickness of the
porous medium, the Darcy and Forchheimer numbers, upon the
hydrodynamic and thermal system performances. These researchers
reported that heat transfer can be enhanced up to 8 times and that
the Nusselt number may improve up to 1200% by the inclusion of
a porous substrate. In subsequent works [3,11], these authors suggested that the equipment efficiency can increase for 3e32% in
a modified solar collector, but at the expense of increasing the
* Corresponding authors. Tel.: þ52 55 5804 4648; fax: þ52 55 5804 4900.
E-mail address: [email protected] (J.A. Ochoa-Tapia).
1290-0729/$ e see front matter 2011 Elsevier Masson SAS. All rights reserved.
doi:10.1016/j.ijthermalsci.2011.03.005
pressure drop up to 30 times. This situation happens regardless
whether parallel or counter-current arrangements are considered
in the system [2]. In addition, it was reported that a critical value of
the porous insert thickness can be reached for which there is no
significant improvement in the thermal performance. Recent
studies have been performed in 2D parallel-plate channels [12] and
in 3D configurations [13]. These works evidence that convection is
not the only heat transport mechanism that is improved. In fact,
Chen and Sutton [14] demonstrated that the heat transfer by
radiation can be augmented about 105% by using ceramic porous
inserts in circular ducts. For turbulent flow, Yang and Hwnag [15]
recently investigated the heat transfer enhancement in a tube
with a porous core for a broad range of the Reynolds number values
(5000e15,000). Moreover, these authors reported that the
optimum porous-radius ratio is around 0.8. Recently, Lukisha and
Prisnyakov [16] carried out a detailed analysis about the conditions
that favor the thermal efficiency using Robin-type boundary
conditions in round channels. They summarize that the optimal
conditions are found by: decreasing the diameter of the channel,
decreasing the Biot number, increasing the porosity and increasing
the Reynolds number.
Pavel and Mohamad [17] experimentally investigated the effect
of metallic porous materials on the heat transfer rate and the
pressure drop. Recently, Huang et al. [6] carried out experiments
inserting a porous medium with high porosity (0.951e0.975) in the
1356
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
Nomenclature
Abs
bb
bs
ci
Cp
(cP)b
e
g
G
H
I
Kb
kb
ks
Kb,u
K*eff
Keff
;u
[b
[s
[cell
li
L
Ltotal
nbs
^
P
hpb ib
hPi0
hPiL
Q
r
r0
t*
ti
hTi
location of the fluidesolid interface, dimensionless
closure variable associated to the fluid phase, m
closure variable associated to the solid phase, m
fit parameter in Eq. (4), i ¼ 0,...,3, dimensionless
mass fraction weighted heat capacity, J/(Kg K)
specific heat capacity at pressure constant of the fluid,
J/(Kg K)
thickness of the porous insert, m
gravity constant, m/s2
Green’s function, dimensionless
channel height, m
identity tensor, dimensionless
permeability tensor, m2
thermal conductivity of the fluid, W/(m K)
thermal conductivity of the solid matrix, W/(m K)
permeability at the bulk of the porous insert, m2
*
total thermal disKeff
persion tensor, W/(m K)
;u
total thermal dispersion in the bulk of the porous
insert, W/(m K)
characteristic length of the pores, m
characteristic length of the solid particles, m
length of the unit cell to reproduce the entire porous
insert, m
lattice vectors for periodic conditions, m
characteristic length of the macroscale, m
channel length, m
unit normal vector pointing from the fluid to the solid,
dimensionless
non-dimensional pressure, dimensionless
intrinsic average pressure, Pa
average pressure at the entrance of the channel, Pa
average pressure at the exit of the channel, Pa
average volumetric flow rate, m3/s
position vector, m
radius of the averaging domain in the upscaling
process, m
characteristic time of heat transport at the pore level, s
fit parameter in Eq. (4), i ¼ 0,...,3, dimensionless
surface average temperature, K
core of a tube. They studied the flow resistance and heat transfer
characteristics experimentally and numerically. In these and other
applications, one of the main problems is to characterize the rapid
variations of transport properties taking place in the vicinity of the
fluideporous dividing surface. Commonly, the boundary conditions
applying at the fluideporous interface are postulated from the
macroscale (normally imposing continuity conditions), which can
lead to physical inconsistencies with the phenomenon taking place
at the microscale [18]. Generally, continuity conditions are
acceptable when the excess (surface) quantities are negligible. But,
evidently, this course of action is not infallible, and it is required to
theoretically build the jump conditions accordingly with the
transport phenomena at the microscale. According to Goyeau et al.
[19], modeling of transport phenomena between a porous medium
and a fluid can be carried out by the following alternatives:
1. The one-domain approach (ODA). The whole system, i.e., the free
fluid zone and the porous insert, are considered as a single
domain. The rapid variations of properties are accounted by
either solving associated closure problems in representative
unit cells or from experimental data. The portion of the system
hTi0
hTiwall
ub
V
vb
Vb
hvb i
^b
v
~v
b
hvb ibu
x
^
x
w
W
average temperature of the fluid input, K
average temperature of the channel walls, K
x-component of the velocity, m/s
volume of the averaging domain, m3
y-component of the velocity, m/s
location of the fluid within the averaging domain,
dimensionless
average velocity, m/s
non-dimensional average velocity, dimensionless
spatial deviations of the velocity, m/s
norm of the intrinsic average velocity in the bulk of the
porous insert, m/s
position vector locating the centroid of the averaging
domain, m
non-dimensional position vector, dimensionless
weight parameter for iterative scheme, dimensionless
channel width, m
Greek letters
thermal diffusivity of the fluid phase, m2/s
thermal diffusivity of the solid phase, m2/s
slip coefficient, dimensionless
fluid phase, dimensionless
Dirac’s delta function, dimensionless
porosity, dimensionless
volumetric fraction of the solid, dimensionless
porosity at the bulk of the porous insert, dimensionless
jump coefficient for heat flux boundary condition,
W/(m2 K)
F
heterogeneous term in the dimensionless momentum
governing equation, dimensionless
k
ratio of thermal conductivities, dimensionless
l
characteristic length, m
h
domain of the free fluid region, dimensionless
mb
fluid viscosity, Pa s
spatial average density, Kg/m3
hri
rb
fluid density, Kg/m3
Q
non-dimensional average temperature, dimensionless
Qm
local mean temperature, dimensionless
x
thickness of the inter-region, m
u
domain of the porous insert, dimensionless
ab
as
aT
b
d
3b
3s
3b,u
f
where the average coefficients exhibit spatial variations is
called the inter-region. In this way, the transport models consist
of effective-medium equations involving position-dependent
coefficients.
2. The two-domain approach (TDA). Here the free fluid zone and
the porous insert are considered as two different homogeneous
domains. In this context, we use the word homogeneous in the
sense of Quintard and Whitaker [20], i.e., a region is said to be
homogeneous when the effective properties (for example
permeability, porosity and effective conductivity) are independent of position. In this approach, it is required to derive the
corresponding jump boundary conditions at the dividing
surface in order to account for transport phenomena in the
inter-region [21]. Furthermore, according to Valdés-Parada
et al. [22] in order to compute the associated jump coefficients
it is necessary to account for the spatial variations of the
effective transport coefficients involved in the ODA.
In a TDA context, the derivation of complete boundary conditions still remains a challenge. Since the pioneer experimental and
theoretical work by Beavers and Joseph [23] about flow in channel
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
partially filled with a porous medium, several efforts have been
devoted to this subject [21,22,24e28]. For the same configuration
used by Beavers and Joseph, Alazmi and Vafai [29] summarized and
compared several boundary conditions reported in the literature
for momentum and heat transfer. They found that momentumtransport jump conditions are more sensitive to variations between
models than those for heat transfer. As a matter of fact, the interfacial boundary conditions in the TDA reported in the literature, are
mainly combinations of continuity or discontinuity of the flux and/
or the fields. In order to clearly point out the variety of reported
jump conditions, we reproduce the information originally reported
by Alazmi and Vafai [29] concerning the boundary conditions for
the heat transport problem in Table 1. There aT y f are unknown
jump coefficients, l is a characteristic length, and hTiu and hTih are
the volume-averaged temperatures corresponding to the porous
medium and the free fluid, respectively.
Within the TDA there are some analytical solutions reported in
the literature. Vafai and Thiyagaraja [30] derived an approximate
analytical expression, based on the matched asymptotic expansion
method, for flow and heat transfer through three types of interface
regions of a porous medium: the porouseporous interface, the
fluideporous interface and the impermeable walleporous interface. Subsequently, Vafai and Kim [31] presented the exact solution
for momentum transport in the fluideporous boundary. Assuming
continuity of the velocity and using the stress jump boundary
condition derived by Ochoa-Tapia and Whitaker [21,32], Kuznetsov
[33] reported the analytical solution for the velocity field, whereas
Silva and Lemos [34] numerically investigated the fluid mechanics
under turbulent conditions. Concerning the forced convection heat
transfer problem, Kuznetsov [35] utilized boundary layer solutions
and continuity of temperatures and heat fluxes to develop the
corresponding analytical solutions. Also, analytical solutions for
unsteady fluid flow in parallel-plates channels were derived by AlNimr and Alkam [36] using Green’s function method.
Regarding the modeling options, and without an in-depth
analysis, we stress that the main differences between the ODA and
TDA contexts are: 1) The inter-region in the ODA is replaced by
a dividing surface in the TDA; 2) The derivation of the jump
conditions requires that the ODA is satisfied on the average in the
inter-region and 3) In the ODA it is necessary to account for the
thickness of the inter-region, whereas in the TDA it is required to fix
the location of the dividing surface. The latter is still an open
problem in the literature [37,38].
Thus, the main objectives of this work are:
To model convective heat transfer in a system composed of two
homogeneous regions (fluid and porous insert), under an ODA
formulation. This requires accounting for the spatial dependence
of the permeability, the porosity and the total effective thermal
conductivity appearing in the effective-medium equations.
To solve the corresponding variable-coefficients partial differential equations, using robust numerical schemes.
To show the extents and limitations of the ODA under several
transport conditions.
1357
We stress that in the ODA context, handled in this work, the
transport properties undergo rapid and continuous changes through
the inter-region [19], which does not necessarily imply the use of
continuity conditions in the equivalent TDA formulation [22]. In
this regarding, Jamet et al. [39] showed that the TDA with continuity boundary conditions is exactly the discontinuous ODA, i.e.,
with stepwise changes in transport properties. Nevertheless,
previous studies devoted to transport in the inter-region [40] have
shown that the effective properties do not, in general, obey a stepwise functionality. These results arise from the statement and
solution of the corresponding closure problems in representative
unit cells for the inter-region. Hence, in this work, we use these
results instead of imposing a stepwise spatial functionality of the
effective properties.
2. System description
In Fig. 1 we provide a sketch of the problem under consideration,
which consists of a parallel-plate channel with a homogeneous
porous insert adhered to the inferior wall. The forced convection is
caused by a pressure drop hPiL hPi0 over the horizontal axis and,
for simplicity, we assume that the temperature of both walls are
fixed at hTiwall .
We commence our modeling effort by adopting the Darcye
Brinkman model for the volume-averaged velocity, under an ODA
context (i.e., involving position-dependent coefficients). This
equation was derived by discarding the inertial terms using the
volume averaging method by Ochoa-Tapia and Whitaker [21] and
more recently by Valdés-Parada et al. [27] and, therefore, it is valid
for low values of the particle Reynolds number (i.e., Rep < 1),
b
2
0 ¼ V pb þrb g þ mb 31
b ðxÞV vb
1
mb K1
mb 31
b ðxÞ$ vb
b ðxÞV3b ðxÞ$V 3b ðxÞ vb
Here 3b(x) is the porosity evaluated at position x, Kb is the
permeability tensor, mb and rb are the fluid viscosity and density,
respectively, whereas hpb ib and hvb i are the volume-averaged
pressure and velocity, respectively. It should be noted that in Eq. (1)
the dependent variables (pressure and velocity) are defined
everywhere in the system, i.e., in the homogeneous and heterogeneous regions.
In the context of the method of volume averaging [41], the
predictions of effective-medium coefficients is carried out by
solving associated closure problems in representative zones of the
microstructure (i.e., unit cells). To have a clear idea about how we
computed the spatial dependence of the average coefficients
involved in Eq. (1), we sketch several unit cells at different positions
in the boundary region in Fig. 2. Here, the microstructure of the
Table 1
Types of jump boundary conditions for heat transfer at the fluideporous interface.
Model
Temperature
Temperature gradient
I
hTiu ¼ hTih
*
Keff
;u
dhTih
dhTiu
¼ kb
dy
dy
II
hTiu ¼ hTih
*
Keff
;u
dhTih
dhTiu
¼ kb
þf
dy
dy
III
a
dhTiu
¼ T ðhTiu hTih Þ
l
dy
*
Keff
;u
dhTih
dhTiu
¼ kb
dy
dy
(1)
Fig. 1. Parallel-plate channel with a porous insert and characteristic lengths.
1358
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
porous insert is represented as a periodic model of non-touching
squared cylinders and cubes. In this way, according to the centroid
of each unit cell, we can specify the position of each average
coefficient.
For the heat transport process, in the absence of thermal sources, the governing macroscale equation for the average temperature hTi, assuming local thermal equilibrium, is given by [25]
hriCP
vhTi
þ ðrcP Þb V$ vb hTi ¼ V$½K*eff ðxÞ$VhTi
vt
(2)
where hriCP ¼ 3b ðxÞðrcP Þb þ ½1 3b ðxÞðrcP Þs is the averaged
volumetric heat capacity at constant pressure with rs being the
density of the solid composing the porous insert while (cP)b and
(cP)s are the specific heat capacities for the fluid and the solid,
respectively; in addition K*eff is the total thermal dispersion tensor
introduced by Quintard and Whitaker [42]. Once again, Eq. (2) has
been previously derived from the rigorous upscaling of the governing equations at the pore level via the method of volume
averaging by Ochoa-Tapia and Whitaker [25]. As mentioned
8 3 ; y=r < 1
0
b;u
>
2 <
y
3b ðy=r0 Þ ¼ 1 3b;u þ 1 þ 1 3b;u 1 y
3 ;
>
4
r0
r0
:2
1; y=r0 > þ1
above, in order to use Eq. (2), it is necessary to account for the
spatial dependence of the effective parameter, K*eff , as well as the
averaged velocity field. Nevertheless, in this case, the velocity is
independent of the heat transport problem, therefore the velocity
is treated as input data in Eq. (2). In the following sections we
discuss the spatial dependence of the coefficients involved in Eqs.
(1) and (2) to later present the solution of the macroscale models.
3. Implicit integral solution of the average velocity
3.1. Spatial dependence of porosity and permeability
The upscaled momentum-transport equation (Eq. (1)) is
expressed in terms of effective-medium coefficients that depend, in
general, of all the spatial coordinates. However, for the sake of
simplicity we assume that the most important variations take place
along the vertical direction (i.e., the y-direction). For a periodic
model of spheres representing the porous insert, Valdés-Parada
et al. [27] provide the following expression for the spatial dependence of the porosity
(3)
1 y=r0 þ1
where 3b,u is the porosity evaluated in the bulk of the porous insert
and r0 is the radius of the averaging domain in the upscaling
process. Here we utilized the normalized coordinate y/r0 to set the
interval of major changes to be [1, þ1]. It is worth stressing that, in
the context of the one-domain approach, the interval of major
changes is the so-called inter-region and generally its size, x, is close
to 2r0 [21] (see Fig. 2). Notice that in Eq. (3), the coordinate axis is
located at the fluideporous medium boundary.
Recently, Valdés-Parada et al. [22] computed the spatial
dependence of the permeability for a periodic arrangement of nontouching squared cylinders. This situation corresponds to solving
the closure problem that defines the permeability (see Appendix A
in [22]) in 2D unit cells as those shown in Fig. 2. However, we notice
from our simulations that the permeability strongly depends of the
microscale geometry in which the closure problem is solved. With
this in mind, in order to obtain more realistic predictions, the
closure problem presented in ref. [22] is solved in 3D unit cells as
those shown in Fig. 2. In this way, the spatial dependence of the
horizontal component of the permeability was fitted (R2 > 0.999)
according to the expression
8
1; y=r0 < 1
>
>
<
h y
3
Kb;u
P
ci exp t0 ti1 ; 1 y=r0 þ1
¼ c0 þ
r0
Kb ðy=r0 Þ >
>
i¼1
:
0; y=r0 > þ1
(4)
Table 2
Values of coefficients in Eq. (4).
Fig. 2. Unit cells for predicting the effective coefficients at different positions in the
inter-region.
3b,u
c0
t0
c1
t1
c2
t2
c3
t3
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0.175
0.124
0.144
0.231
0.284
0.292
0.557
1.006
1.000
1.000
1.000
1.000
1.045
1.074
0.189
0.109
0.362
0.392
0.388
0.445
0.539
0.012
0.100
0.282
1.364
0.744
1.360
1.956
0.674
0.516
0.387
0.436
0.388
0.445
0.539
0.340
0.405
1.177
0.568
0.858
1.360
1.956
0.402
0.498
0.395
0.405
0.506
0.445
0.539
2.375
1.411
1.186
1.754
2.491
1.360
1.956
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
1359
3.2. Implicit solution with Green’s functions
To conclude this section, we solve the governing equation for
momentum transport. To this end, it is convenient to introduce the
following non-dimensional variables and parameters
^
x¼
b
K
v
x
yþe
^ ¼ pb H
^b ¼ b ; Da ¼ b ; P
^¼
; y
; v
H
H
H2
nb bu
mb nb bu
(6)
where hnb ibu is the characteristic velocity in the bulk of the porous
insert and Da is the Darcy number. In addition, notice that the
origin of the y-axis is located at the walleporous medium
boundary. In this way, Eq. (1) takes the form
ˇ
Fig. 3. Spatial dependence of permeability by solving the closure problem in 3D unit
cells as a function of the porosity.
where Kb,u is the permeability evaluated in the bulk of the porous
insert and the values of coefficients ci and ti are given in Table 2 as
functions of the porosity. Moreover, in Fig. 3 we show the spatial
dependence of Kb for several porosity values, noticing that the
profiles exhibit a similar dependence as the one shown in Fig. A-2
in ref. [22]. It is stressed that in this work we only present simulations for 3b,u ˛ [0.3, 0.9] due to the prohibitive computational
demands for 3b,u < 0.3.
Finally, in order to use Eq. (4) we need the bulk permeability,
Kb,u, which can be computed from the solution of the closure
problem stated in [22] within a representative unit cell of the
porous insert bulk. For a periodic arrangement of non-touching
square cylinders and cubes the permeability as a function of the
porosity is presented in Fig. 4. These results were satisfactorily
fitted (R2 > 0.999) according to the expression,
Kb;u =[2cell ¼ 2:430 103 þ 1:360
106 exp 13:9003b;u ;
Kb;u =[2cell
¼ 1:606 10
5
10
^Þ
3 ðx
^ þ V3
^P
^ b ðx
^ 31 ðx
^ 2v
^b
^ ÞV
^ Þ$V
^ Þv
^b þ b
^ b ¼ 3b ð x
I$v
V
b
^Þ
Daðx
were it has been assumed that the permeability tensor is isotropic,
i.e., Kb ¼ KbI, where I is the identity tensor.
The system configuration suggests imposing the following
simplifications
^nb ¼ 0;
(5a)
^
vP
¼ 0;
^
vy
^b
d2 u
^
dy
(9)
F u^ b ; y^
"
2 #
^
dP
d
1
1
^b
^Þ
^ÞDa ðy
^ Þ 3b ð y
^ Þ 3b ð y
^Þ
u
¼ 3b ð y
þ 3b ð y
^
d^
x
dy
^
du
d
3b ðy^Þ b
^
^
dy
dy
(10)
(5b)
here [cell is the length of the unit cell used to solve the closure
problem in the porous medium bulk.
(8)
where
^
þ 31
b ðyÞ
3D unit cells
v 3b
¼ 0
v^
x
^
^b ; y
¼ F u
2
4
þ 2:047
exp 8:3973b;u ;
^b
vu
¼ 0;
v^
x
^ b and ^nb are the x- and y-components of the velocity field.
where u
Taking these simplifications into account, we can write Eq. (7) as
the following scalar equation for the x-component of the velocity
2D unit cells
(7)
which is subject to the following boundary conditions:
^ ¼ 0; 1
at y
^b ¼ 0
u
(11)
In its current form, it is not easy to find an analytical solution to
Eq. (9) because of the non-trivial spatial dependence of the
permeability and the porosity. With this in mind, it is convenient to
regard the terms on the right-hand side of Eq. (9) as a single
heterogeneous term. This is done in order to find an implicit
analytical solution using integral equation formulations in terms of
Green’s functions as described in ref. [43]. The associated Green’s
function solves the following boundary-value problem,
^; y
^0 Þ
d2 Gðy
2
^
dy
At
Fig. 4. Bulk permeability as a function of the porosity for 2D and 3D unit cells.
^y
^0 Þ
¼ dðy
^ ¼ 0; 1;
y
^; y
^0 Þ ¼ 0
Gðy
(12a)
(12b)
where d represents the Dirac’s delta function. The solution of this
problem is straightforward and can be expressed as follows
1360
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
( bb
y y0 1 ;
^; y
^0 Þ ¼
Gðy
b
y1 ;
y0 b
b
y<b
y0
(13)
b
y>b
y0
Now, let us define the linear operator
Lð Þ ¼
d2
^
dy
2
ðÞ
(14)
that F can be simply modified in order to explore different alternatives of the spatial dependency of the coefficients. This is one of
the main advantages of using integral formulations over other more
classical schemes based on discretization of the differential equation. The evaluation of the solution given in Eq. (17) is provided in
Section 5.
4. Determination of the average temperature field
thus, we have the equalities
^b ; y
^b ¼ F u
^
L u
(15a)
The macroscale heat transport problem consists of the differential Eq. (2), subject to the boundary conditions
^y
^0 Þ
LðGÞ ¼ dðy
(15b)
hTi ¼ hTiwall ;
Subsequently, the Green’s formula, adapted to our problem, is
given by [44]
^Z¼ 1
y
^¼0
y
h
i
^ b y^¼1
du
dG
^ b LðGÞ GL u
^ ¼ u
^
^ b dy
u
G
^
^ y^¼0
dy
b dy
(16)
Eventually, after taking into account the boundary conditions
^ b and G, along with the filtration property of the Dirac’s delta
for u
function, we have
^ b ðy
^Þ ¼
u
^0 ¼ 1
y
Z
F u^ b ; y^0 Gðy^; y^0 Þdy^0
(17)
^0 ¼ 0
y
This result represents the formal solution to the momentumtransport problem. This is an implicit integral solution as it is
required to integrate the source term throughout the entire domain
evaluated with the same field solution. Therefore, it is necessary to
implement an iterative scheme to find the velocity field. This type
of implicit integral formulations have been utilized, for example, in
solving transport-reaction mass problems with linear and nonlinear kinetics in basic geometries [43,45,46].
In Eq. (17), the Green’s function represents the influence of the
heterogeneous term at the position y0 over the velocity field
located at y. Also, the opposite case is true because of the symmetry
property of the Green’s function, i.e., G(y,y0) ¼ G(y0,y). Also, note
that the non-homogeneous term can be as complicated or as simple
so that the only difficulty appearing, is related to the stability of the
iterative scheme used to find the solution. This means, for example,
ˇ
vhTi
*
ðrcP Þb ub hTi0 ¼ ðrcP Þb ub hTi Keff
;
xx vx
vhTi
¼ 0;
vx
(18a)
At x ¼ 0;
At x ¼ Ltotal
(18b)
(18c)
where (Keff*)xx represents the xx-component of the thermal tensor.
The boundary condition in Eq. (18b) is known as Danckwerts-type
[47], which states that at the boundary x ¼ 0 the heat entering by
the left side, ðrcP Þb ub hTi0 , is the same escaping by the right side,
* Þ vhTi=vx. Notice that this condition reduces to
ðrcP Þb ub hTi ðKeff
xx
a Dirichlet-type for sufficiently high flow rates, i.e., hTizhTi0 .
For the sake of brevity we omit the details concerning the
determination of K*eff ðxÞ which is predicted through the solution of
the corresponding closure problems provided in the Appendix.
Under purely conductive conditions, the numerical results are
straight lines as shown in Fig. 5 in [40] and obey the following
relation
1
y
*
*
k
1
þ
I;
K*eff ¼ Keff
þ
K
;u
eff ;u
2 b
r0
Pecell < 1
(19)
*
where kb is the thermal conductivity of the fluid and Keff
repre;u
sents the total effective thermal conductivity in the bulk of the
porous insert. In Eq. (19), we have introduced the cell Péclet
number as
ˇ
ˇ
ˇ
ˇ ˇ
a
At y ¼ e; H e
Pecell ¼
D Eb
ðrcP Þb nb [cell
kb
u
(20)
b
Fig. 5. Velocity profiles for: a) several values of the porosity taking the permeability from 3D unit cells and b) different permeabilities computed in 2D and 3D unit cells taking
3b,u ¼ 0.5. In all cases was used e/H ¼ 0.5 and x/H ¼ 0.01.
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
where [cell represents the minimum length of the periodic unit cell
to completely reproduce the homogeneous porous insert. Generally, [cell is smaller than the macroscale characteristic lengths H or
Ltotal, therefore, we can affirm that
Pecell < PeH
(21)
where PeH is the Péclet number based on the channel height, i.e.,
PeH h HPecell/[cell.
In order to use Eq. (19), it is necessary to account for Keff, u. As
explained by Ochoa-Tapia et al. [48], an analytical expression for
this coefficient can be obtained using Chang’s unit cell under
diffusive regime. The results for two-dimensional arrays of nontouching cylinders are given by
2k 3b;u ðk 1Þ
Keff ;u
¼
2 þ 3b;u ðk 1Þ
kb
5. Results and discussion
5.1. Hydrodynamic performance
As explained in Section 3, in order to obtain the volume-averaged velocity fields, we need to implement an iterative scheme. In
this case, we use the recursive formula
^b
u
jþ1
^b
¼ w u
j1
^b
þð1 wÞ u
a
where j refers to the number of iteration and w is the weight
parameter, which ranges from 0 to 1, and it strongly influences the
scheme convergence. In our computer code, the iterative scheme
was forced to satisfy an error tolerance based on the L2 e norm, i.e.,
8
91
^
>
< yZ¼ 1 h i2 >
=2
^b u
^b
^
u
error ¼
dy
>
>
j
j1
:
;
j
(23)
(24)
^¼0
y
in our simulations we stopped the iteration process whenever
error 104. The pressure drop was fixed as
D
D pb
Eb
Ltotal
(22)
For situations in which conduction and convection compete, the
spatial dependence of K*eff is no longer lineal and its behavior is
similar to the results shown in Fig. 6 in [40]. It should be stressed
*
under convective regime,
that in order to compute the value of Keff
;u
we cannot use Eq. (19) and we must solve the corresponding closure
problem in representative cells of the homogeneous porous insert.
In other words, we require to solve Eq. (34aeh) and then substitute
the closure variable field, bb, in Eq. (32) to compute the total thermal
tensor (see Section 6.1 for details).
Notice that, contrary to the momentum-transport problem, in
this case there are serious difficulties related to the dispersive term
in order to find an analytical solution based on Green’s functions.
For this reason, the solution of the heat transfer problem is carried
out by numerical methods. The results for hydrodynamic and
thermal performances are provided in the next section.
1361
¼ 1 Pa m1
(25)
Thus, we are able to compute the characteristic velocity hnb ibu,
when the channel is fully filled with the porous insert. In this way,
we have not fixed the channel Reynolds number,
D Eb
ReH ¼
rb nb
mb
u
H
(26)
because the characteristic velocity depends, in general, of the
porosity. Of course, an alternative procedure is to prescribe the
value of the channel Reynolds number and then find the corresponding pressure drop. Indeed, both approaches are equivalent.
In Fig. 5 we plot the spatial dependence of the velocity.
According to Fig. 5a, the velocity is strongly influenced by the
porous medium porosity, which is in agreement with previous
theoretical and experimental works [19,23,32]. Notice that, as the
porosity approaches 1, the resistance exerted by the insert
decreases and the velocity profile is similar to that one found
between two parallel plates (i.e., Poiseuille flow). For the opposite
case, i.e., 3b,u / 0, eventually the insert acts like an impermeable
wall. It is worth noting that approximately at the half of the porous
^z0:25Þ the non-dimensional velocity equals the porosity
insert ðy
due to the definition adopted in Eq. (6).
Furthermore, in Fig. 5b, we show the importance of the unit cell
dimensionality in the predictions of the velocity profiles. We
observe that the largest differences take place at the peak of the
^ b;max ). In fact, for the case under
velocity profiles (i.e., u
b
Fig. 6. Velocity profiles for: a) three sizes of the porous insert and taking x/H ¼ 0.01, b) three sizes of the channel height and taking e/H ¼ 0.5. In all cases was used 3b,u ¼ 0.5 and the
permeability from 3D unit cells.
1362
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
^ b;max corresponding to 2D
consideration, we have that the value of u
unit cells is about one third of the one obtained using 3D unit cells.
To the best of our knowledge, this influence of the unit cell
dimensionality has not been reported in the literature and it is
crucial in macroscale modeling.
Following the experimental observations by Beavers and Joseph
[23], the velocity profiles are also influenced by the configuration of
the channel partially filled with the porous medium. In Fig. 6a we
present the velocity profiles for three different sizes of the porous
insert (e/H ¼ 0.25, 0.5, 0.75). As expected, the velocity increases as
long as the porosity approaches to 1 or the porous insert thickness
decreases.
One of the primordial requirements in the upscaling process, is
to identify the scaling laws that allow discarding redundant information from the pore scale [49]. In other words, we need to identify
the disparity of length-scales in order for the macroscale equations
to be valid. Usually, this disparity is expressed in the form of the
inequality [b r0 L, with [b and L being the characteristic
lengths associated to the microscale and the macroscale, respectively. While a constraint in the above form is certainly informative,
it may not be sufficient for practical applications. For this reason, it
is necessary to determine what exactly is meant, by “r0 L”, for
example. With this aim, we investigate the influence of the ratio x/H
over the velocity profiles in Fig. 6b. We observe that the profiles
collapse in a single curve when
x=H 0:01
(27)
This means that, the averaging domain used to upscale (or to
compare with experimental data or direct numerical simulations)
must be, at least, 1% of the height channel. It is inferred that the
characteristic length of the pores must be smaller than that of the
average volume (i.e., [b r0 ). Nevertheless, an exhaustive study
involving all the design variables as the properties and shape of the
materials composing the walls and the porous matrix, etc., is
indeed in order but it lies beyond of the objectives of this work and
will be pursued in a future study.
To finalize this section, as a quantitative measurement about the
influence of the porous insert over the hydrodynamic performance
of the channel, let us introduce the average volumetric flow rate per
unit of channel width as
Q
¼
W
y¼
Z He
D
E
ub dy
(28)
y ¼ e
a
In Fig. 7a we present the variation of Q with the size of the
porous insert. Notice that the flow rate is divided with the one
corresponding when there is no porous insert within the channel.
We observe an exponential decrease of the flow rate as long as the
porous insert fills the parallel-plate channel, and an asymptotic
value, Q/Qe¼0 ¼ 0.0164, is reached for e/H ¼ 1. It is worth noting
that, it is possible to compute the volumetric flow rate directly from
Darcy’s law when the porous insert completely saturates the
channel. In addition, in Fig. 7b the influence of the porous medium
porosity over Q is presented for three values of the porous insert
thickness. Contrary to Fig. 7a, here we used the permeabilities from
2D unit cells since this allows us to simulate for porosities close to
zero. A minimum value is observed in all the curves and its corresponding porosity changes with the porous medium thickness.
From our simulations, we noticed the absence of a minimum value
when the channel is completely filled with the porous medium. In
that case, the behavior of the curve is similar to those shown in
Fig. 7b. In a general way, if the porosity is less than 0.6, the flow rate
moderately decreases.
5.2. Thermal performance
With the velocity field satisfying the iterative scheme of Eq. (23),
this is supplied as input data for the heat transport problem. This
problem was solved using the finite element solver Comsol Multiphysics 3.5a, which uses adaptive mesh refinements to ensure that
the results are not dependent upon the number of mesh nodes. In
addition, we used the UMFPACK routine, provided in Comsol, to
carry out the factorized matrices inversion. The main difficulty in
the numerical solution lies on the correct quantification of the
spatial changes of the total thermal tensor K*eff ðxÞ. From the above
developments we notice that there are several degrees of freedom
for the solution of the heat transport problem, here we only consider
the following: the channel Péclet number, PeH; the ratio of thermal
conductivities in the porous medium, k; the porosity in the bulk of
the porous medium, 3b,u and the size of the porous insert, e/H. In
each case we hold all the other degrees of freedom constant while
variating one at the time; our observations for each numerical
experiment are provided below.
In Fig. 8 we present the non-dimensional temperature fields,
defined as
Q¼
hTi hTi0
hTiwall hTi0
(29)
b
Fig. 7. a) Effect of the size of the porous insert on Q taking the permeability from 3D unit cells, 3b,u ¼ 0.5 and x/H ¼ 0.01. b) Effect of the porosity on Q taking the permeability from
2D unit cells, for e/H ¼ 0.25, 0.5 and 0.75 taking x/H ¼ 0.01.
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
1363
Fig. 8. Effect of the Péclet number on the non-dimensional temperature field in the channel: a) PeH ¼ 0.1, b) 1, c) 10 and d) 100. Here was taken k ¼ 100, e/H ¼ 0.5, x/H ¼ 0.01,
3b,u ¼ 0.5 and the permeability computed in 3D unit cells.
for different channel Péclet numbers, PeH. We observe that, for the
conditions here considered, under moderate flow rates (i.e.
PeH < 1), practically all the fluid passing through the channel
acquires the wall temperature. In the opposite case, for PeH > 1,
convective heat transport predominates over axial conduction, and,
in this situation, only the fluid near the walls is at hTiwall . Moreover,
we notice that the portion of the system occupied by the porous
insert has a stronger tendency to preserve the wall temperature;
because of this we notice drastic variations of the temperature near
the fluideporous medium boundary, especially for PeH > 1. In
addition, under these conditions, the fluid outside the porous
medium conserves practically the entrance temperature. It is worth
noting that, these results might change slightly for larger (or
shorter) channel lengths, because of the augmentation (or
decrease) of the contact time between the fluid and the channel
walls.
The effect of using porous inserts with different values of the
thermal conductivity for the solid matrix is shown in Fig. 9. As
expected, the more conductive is the insert, the more heat that is
transferred from the channel walls to the fluid (and the opposite is
true). Notice that for k ¼ 0.01 the fluid does not reach the
temperature of the walls at the exit channel, whereas for k ¼ 100
practically all the fluid leaves the channel at hTiwall . These results
suggest that one way to improve heat exchangers performance is by
using porous substrates with high thermal conductivity and with
large interfacial area available for heat exchange between the
interstitial fluid and the solid composing the substrate.
Regarding the influence of porosity, in Fig. 10 we show the
channel temperature profiles for 3b,u ¼ 0.4, 0.7 and 0.9. It should be
stressed that the channel Péclet number changes with the porosity,
however the applied pressure drop was the same in all cases and it
corresponds to the one for PeH ¼ 10 and 3b,u ¼ 0.5. For the sake of
illustration, a channel fully filled with the porous insert was
utilized. For low values of the porosity, 3b,u < 0.3, the temperature
profile is similar to the one found in a solid rod with same
dimensions as the channel. The opposite case occurs for porosity
values close to 1. In this case, the temperature and velocity profiles
resemble those found in parallel plates separated by a gap of size H;
this is due to the negligible resistance offered by the porous matrix.
In this way, it is expected that low porosities favor the heat
exchange.
Finally, the influence of the size of the porous insert is presented
in Fig. 11. As more substrate fills the channel, the dimensionless
temperature profiles are more uniform and close to 1. It is interesting to note that only the fluid near the walls increases its
temperature in all the cases studied in this section. In fact, if no
Fig. 9. Effect of the ratio of thermal conductivities on the non-dimensional temperature field in the channel: a) k ¼ 0.01 and b) 100. Here was taken PeH ¼ 10, e/H ¼ 0.5, x/H ¼ 0.01,
3b,u ¼ 0.5 and the permeability computed in 3D unit cells.
1364
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
Fig. 10. Effect of the porosity on the non-dimensional temperature field in the channel: a) 3b,u ¼ 0.4, b) 0.7 and c) 0.9. Here was taken e/H ¼ 1, x/H ¼ 0.01, k ¼ 10 and the pressure
drop corresponding to PeH ¼ 10 for 3b,u ¼ 0.5 and using the permeability computed in 3D unit cells.
porous insert is present, the temperature field is the typical one
encountered in most heat exchangers. Even though the usage of
porous substrates leads to increasing pressure drops [11], their
utility can be noted clearly in the cases plotted in Fig. 11.
Moreover, in order to quantify the thermal performance of the
heat exchanger let us introduce the local Nusselt number referred
to the lower wall [11,12] as,
3
2 * K
1 4 eff yy dQ5
Nux ¼
^ ^
kb
dy
1 Qm
y¼0
(30)
where a mean temperature weighted with the Darcy velocity has
been introduced as
Z
Qm ¼
^¼1
y
^¼0
y
Z
Qu^ dy^
b
^¼1
y
^ b dy
^
u
(31)
^¼0
y
We begin our analysis of the local Nusselt number by varying
the ratio of thermal conductivities as shown in Fig. 12a. For the
values of parameters used, we found that the velocity and
temperature profiles are practically developed at x ¼ H. For this
reason, we only present the results for ^
x˛½0; 1:5. In all cases, at the
entrance channel the Nusselt number is higher than in the rest of
the channel, since in that location the mean temperature difference
(1Qm) is the largest. Notice that the most important changes
occur for k > 1. It is worth noting that the Nusselt number, as
defined in Eq. (30), is a direct measurement of the heat transfer rate
of the thermal equipment. Because of this, high values of the
Nusselt number are expected for highly conductive porous inserts.
This point is in agreement with the tendency of the temperature
fields shown in Fig. 9.
The influence of the porosity over the thermal performance of
the channel is summarized in Fig. 12b. We observe that the heat
transfer rate is improved as the porous medium porosity is
decreased. In addition, an extra phenomenon occurs inside the
porous insert in real situations. This is the turbulence, which is
originated by the numerous and complex pore networks [50].
However, a thorough investigation of turbulence is beyond the
scope of this work.
In Fig. 13a we present the effect of the flow intensity over the
local Nusselt number. For high Péclet number values (i.e., PeH > 10)
the heat transfer rate is greater at the entrance, but this situation is
reversed as the fluid passes through the channel. This phenomenon
is due to the fact that in the case of conductive transport (i.e.,
PeH < 1) most of the fluid acquires the wall temperature in a shorter
distance within the channel. In the opposite case, for PeH > 1, the
Fig. 11. Effect of the size of the porous insert on the non-dimensional temperature field in the channel: a) e/H ¼ 0.0, b) 0.5 y c) 1.0. Here was taken PeH ¼ 10, k ¼ 100, 3b,u ¼ 0.5, x/
H ¼ 0.01 and the permeability computed in 3D unit cells.
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
a
1365
b
Fig. 12. Influence on the local Nusselt number of: a) the ratio of thermal conductivities taking 3b,u ¼ 0.5 and PeH ¼ 1; b) the porosity taking k ¼ 10 and the pressure drop satisfying
PeH ¼ 1 for 3b,u ¼ 0.5 and e/H ¼ 1. In all cases was taken e/H ¼ 1, x/H ¼ 0.01 and the permeability computed in 3D unit cells.
a
b
Fig. 13. Influence on the local Nusselt number of: a) the Péclet number taking e/H ¼ 0.5; b) the size of the porous insert taking PeH ¼ 1. In all cases was taken k ¼ 10, 3b,u ¼ 0.5, x/
H ¼ 0.01 and the permeability computed in 3D unit cells.
fluid needs to travel a larger distance in order to reach such
temperature. It is worth stressing that these observations are only
applicable to the channel configuration here considered. For
example, it is not hard to imagine that if the length of the channel is
increased so that the maximum value of ^
x[1, the heat performance of the system will be favored at low channel Péclet number
values (i.e., PeH 1).
In Fig. 13b we plot the results varying the thickness of the porous
insert. Our results are in agreement with those presented in Fig. 11,
i.e., an increment in the size of the porous insert leads to a better
heat transfer rate. However, this is not directly proportional to the
quantity of porous material within the channel, i.e. there is a certain
thickness (close to 25% of the channel height) for which the Nusselt
number reaches a minimum value. Nevertheless, it can be
concluded that better thermal performance corresponds for
a channel fully filled with the porous insert. It is interesting to
notice that there are some cases in which the Nusselt number is
a decreasing function of the axial position and there are other cases
where the opposite happens. Directing the attention to the results
shown in Fig. 13b we observe that if the channel is fully filled with
the porous medium or the fluid, the Nusselt number is maximum at
the entrance and exponentially decreases with position. However,
if 0.5 e/H < 1 the porous medium enhances the thermal flux at
the lower channel wall. This is attributed to the fact that the ratio of
thermal conductivities is 10 in this case, which means that the solid
matrix is a better heat conductor than the fluid phase. In this way,
the portion of the channel that is only occupied by the fluid behaves
like a thermal sink that generates considerable temperature
gradients (see Fig. 11b) for the e/H ¼ 0.5 case. For this reason, we
observe that the Nusselt number increases with the axial position,
reaching a maximum and constant value sufficiently far from the
entrance (i.e., ^
x > 2).
Summing up, in order to improve the heat transfer rate, it is
recommendable to:
Use a highly conductive porous insert.
Use a porous insert with low porosity.
Maintain, on the one hand, channel Péclet number values
smaller than 1 for large channels (i.e., ^
x[1) and on the other
hand, for small channels (i.e., ^
x < 1), it is convenient that
PeH [1.
Use a channel fully filled with the porous insert.
It is clear that a more detailed study, for example including costs
due to pressure drops and porous materials, needs to be done in
order to assess the viability of the above recommendations in
practical application in order to optimize the heat transfer rate.
6. Conclusions
In this work we present the modeling of momentum and heat
transfer in parallel-plate channels partially filled with a porous insert
1366
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
using a one-domain approach. Thus, volume-averaged equations are
utilized for momentum and heat problem that contain positiondependent effective properties, namely: the porosity, permeability
and total thermal dispersion tensor. These effective coefficients
undergo rapid changes near the inter-region separating the homogeneous porous insert and the zone of free fluid. With this methodology, we avoid specifying the boundary conditions at the
fluideporous insert boundary at the cost of solving more complicated
transport models.
We showed that, for the DarcyeBrinkman model, an integral
formulation based on standard Green’s functions leads to an implicit
solution for the velocity field. In this way, an iterative scheme was
implemented in order to set the velocity field within the channel. We
found that significant differences in the velocity profiles arise when
the permeability is computed from 2D and 3D unit cells. Thus, we
emphasize the importance to account for accurate methods to
predict effective properties. The hydrodynamic performance (which
encompasses the velocity profiles and the volumetric flow rate) of
the porous channel was analyzed in terms of pertinent parameters
such as the porosity and the size of the porous insert. In addition, it
was found that the minimum requirement to satisfy disparity of
length-scales in the system (needed for the upscaling process), was
that the size of the averaging domain being lesser than 1% of channel
height.
Finally, concerning the heat transport problem, we presented
the temperature fields and the local Nusselt number by varying the
flow intensity, the ratio of thermal conductivities, the thickness of
the porous insert and the porosity. As expected, and also in
agreement with previous works, we found that the heat transfer
rate between the fluid and the channel walls can be increased by
either: using highly conductive porous inserts, decreasing the
porous medium porosity, set large or low values of the channel
Péclet number depending upon the channel length and saturate as
much as possible the channel with the porous insert.
The results from this work provide useful guidelines to design
future experimental and theoretical studies of momentum and heat
transfer in channels partially filled with a porous medium by
supplying information about disparity of length-scales and values
of the non-dimensional numbers governing the heat transport. It is
stressed that, despite the fact that we were able to solve the ODA in
this work, this does not necessarily mean that this alternative
would be optimal in practical situations. For example, the considerable computational demands involved in the solution of partial
differential equations with variable coefficients may turn out to be
impractical in continuous monitoring applications, as those
required in control systems. For this reason, it is certainly of interest
to compare the results from this work with those arising from TDA
formulations, which involve solving more (at least twice as much as
those in the ODA) partial differential equations, with the difference
that they are expressed in terms of constant effective-medium
coefficients. Nevertheless, the predictions from the TDA are driven
by the jump boundary conditions between the porous medium and
the fluid. As mentioned in the introduction, in order to compute the
coefficients involved in the jump conditions, it is necessary to
account for the spatial variations of the effective transport coefficients from the ODA computed here. For this reason, the material
presented in this work should be useful in future studies involving
the derivation, solution and comparison of the one- and twodomain approaches.
dispersion tensor definition arising from using the method of
volume averaging [41],
h
i
kb ks
K*eff ðxÞ ¼ kb 3b ðxÞ þ ks 3s ðxÞ I þ
V
ðrcP Þb
V
Z
Z
nbs bb dA
Abs ðxÞ
~ b dV
v
b b
(32)
Vb ðxÞ
where 3s is the volumetric fraction of the solid, bb is the so-called
closure variable, nbs is the unit normal vector pointing from the
fluid to the solid (nbs ¼ nsb), V is the volume of the averaging
domain, Abs represents the spatial location of the fluidesolid
interface, Vb is the portion of the averaging domain occupied by
~ represents the spatial deviations of the velocity.
fluid and v
b
The closure variable maps the gradient of the average temperature field onto the spatial deviations of the fluid temperature field
at the pore level, according to the expression
T~ b ¼ bb $VhTi
(33)
Moreover, the closure variable, bb, is determined by solving the
following boundary-value problem,
2
1
~
ðrcP Þb vb $Vbb þ v
b ¼ kb V bb ðrcP Þb 3b c b ;
in the fluid
(34a)
0 ¼ ks V2 bs ðrcP Þs 31
s cs ;
in the solid
(34b)
nbs $ kb Vbb ¼nbs $ðks Vbs Þþnbs ks kb ; at the fluid
solid interface
bb ¼ bs ;
at the fluid solid interface
(34c)
(34d)
bb ðr þ li Þ ¼ bb ðrÞ; bs ðr þ li Þ
¼ bs ðrÞ; periodicity at the vertical boundaries
D E
bb ¼ hbs i ¼ 0;
cb ¼
ab Z
V
restriction on the average
(34e)
(34f)
nbs $Vbb dA
(34g)
nsb $Vbs dA
(34h)
Abs
cs ¼
as
Z
V
Abs
where ab and as are the thermal diffusivities of the fluid and the
solid, respectively, r is a position vector, li are the lattice vectors for
periodicity and bs represents the closure variable associated to the
solid phase. In this case, bs maps the gradient of the average
temperature field onto the spatial deviations of the point-wise solid
temperature field according to the expression,
Appendix A. Closure problem for heat transport
T~ s ¼ bs $VhTi
In this section, we present the definition of the total thermal
dispersion tensor and the closure problem associated to the heat
transport problem. Let us begin by providing the total thermal
In order to predict the spatial dependence of the total thermal
dispersion tensor, the boundary-value problem defined by Eq. (34)
is solved in representative unit cells of the inter-region as those
(35)
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
sketched in Fig. 2. In obtaining the closure problem given by Eq. (34),
and the representation of the spatial deviations of the temperature,
Eqs. (33) and (35), some constraints have been imposed which,
generally, are satisfied due to the disparity of length-scales that
multiphase systems present. These constraints are:
r0 L;
r02 L2 ;
[b ; [s L
(36)
where [s is the characteristic length of the solid particles. In addition, in order to consider the closure problem as one in quasisteady state, we impose the time constraint
[2b =kb
t*
D Eb
nb [b
u þ1
kb
(37)
where t* represent the characteristic time of heat transport at the
pore level. We stress that the presentation of the numerous details
involved in the solution of the closure problem are beyond the
scopes of this work. Thus, for a better comprehension we refer the
reader to [40].
References
[1] M.K. Alkam, M.A. Al-Nimr, Transient non-darcian forced convection flow in
a pipe partially filled with a porous material, International Journal of Heat and
Mass Transfer 41 (2) (1998) 347e356.
[2] M.K. Alkam, M.A. Al-Nimr, Improving the performance of double-pipe heat
exchangers by using porous substrates, International Journal of Heat and Mass
Transfer 42 (1999) 3609e3618.
[3] M.K. Alkam, M.A. Al-Nimr, Solar collector with tubes partially filled with
porous substrate, Journal of Solar Energy Engineering 121 (1) (1999) 20e24.
[4] M.A. Al-Nimr, M.K. Alkam, Unsteady non-darcian forced convection analysis in
an annulus partially filled with a porous material, Journal of Heat Transfer 119
(4) (1997) 799e804.
[5] S.Y. Byun, S.T. Ro, J.Y. Shin, Y.S. Son, D.Y. Lee, Transient thermal behavior of
porous media under oscillating flow condition, International Journal of Heat
and Mass Transfer 49 (2006) 5081e5085.
[6] Z.F. Huang, A. Nakayama, K. Yang, C. Yang, W. Liu, Enhancing heat transfer in
the core flow by using porous medium inserts in a tube, International Journal
of Heat and Mass Transfer 53 (2010) 1164e1174.
[7] K.N. Shukla, Diffusion Processes During Drying of Solids, Vol. 11 of Series in
Theoretical and Applied Mechanics. World Scientific Publishing Co., 1999.
[8] O. Levenspiel, Chemical Reaction Engineering, third ed. John Wiley & Sons,
1999.
[9] W.L. McCabe, J.C. Smith, P. Harriot, Unit Operations of Chemical Engineering,
seventh ed. McGraw Hill, 2005.
[10] M.H. Dickson, M. Fanelli, Geothermal Energy: Utilization and Technology.
Earthscan Publications, 2005.
[11] M.A. Al-Nimr, M.K. Alkam, A modified tubeless solar collector partially filled
with porous substrate, Renewable Energy 13 (2) (1998) 165e173.
[12] M.K. Alkam, M.A. Al-Nimr, M.O. Hamdan, Enhancing heat transfer in parallelplate channels by using porous inserts, International Journal of Heat and Mass
Transfer 44 (2001) 931e938.
[13] T.C. Jen, T.Z. Yan, Developing fluid flow and heat transfer in a channel partially
filled with porous medium, International Journal of Heat and Mass Transfer 48
(2005) 3995e4009.
[14] X. Chen, W.H. Sutton, Enhancement of heat transfer: combined convection
and radiation in the entrance region of circular ducts with porous inserts,
International Journal of Heat and Mass Transfer 48 (2005) 5460e5474.
[15] Y.T. Yang, M.L. Hwang, Numerical simulation of turbulent fluid flow and heat
transfer characteristics in heat exchangers fitted with porous media, International Journal of Heat and Mass Transfer 52 (2009) 2956e2965.
[16] A.P. Lukisha, V.F. Prisnyakov, The efficiency of round channels fitted with
porous, highly heat-conducting insert in a laminar fluid coolant flow at
boundary conditions of the third kind, International Journal of Heat and Mass
Transfer 53 (2010) 2469e2476.
[17] B.I. Pavel, A.A. Mohamad, An experimental and numerical study on heat
transfer enhancement for gas heat exchangers fitted with porous media,
International Journal of Heat and Mass Transfer 47 (2004) 4939e4952.
[18] M. Prat, On the boundary conditions at the macroscopic level, Transport in
Porous Media 4 (1989) 259e280.
[19] B. Goyeau, D. Lhuillier, D. Gobin, M.G. Velarde, Momentum transport at
a fluid-porous interface, International Journal of Heat and Mass Transfer 46
(2003) 4071e4081.
1367
[20] M. Quintard, S. Whitaker, Écoulement monophasique en milieu poreux: Effet
des hétérogénéités locales, Journal Mécanique Théorique et Appliquée 6
(1987) 691e726.
[21] J.A. Ochoa-Tapia, S. Whitaker, Momentum transfer at the boundary
between a porous medium and a homogeneous fluid-I e theoretical
development, International Journal of Heat and Mass Transfer 38 (14)
(1995) 2635e2646.
[22] F.J. Valdés-Parada, J. Alvarez-Ramírez, B. Goyeau, J.A. Ochoa-Tapia, Computation of jump coefficients for momentum transfer between a porous medium
and a fluid using a closed generalized transfer equation, Transport in Porous
Media 78 (2009) 439e457.
[23] G.S. Beavers, D.D. Joseph, Boundary conditions at a naturally permeable wall,
Journal of Fluid Mechanics 30 (1967) 197e207.
[24] M. Sahraoui, M. Kaviany, Slip and no-slip temperature boundary conditions at
the interface of porous, plain media: convection, International Journal of Heat
and Mass Transfer 37 (6) (1994) 1029e1044.
[25] J.A. Ochoa-Tapia, S. Whitaker, Heat transfer at the boundary between a porous
medium and a homogeneous fluid: the one-equation model, Journal of Porous
Media 1 (1) (1998) 31e46.
[26] B.D. Wood, M. Quintard, S. Whitaker, Jump conditions at a non-uniform
boundaries: the catalytic surface, Chemical Engineering Science 55 (2000)
5231e5245.
[27] F.J. Valdés-Parada, B. Goyeau, J.A. Ochoa-Tapia, Jump momentum boundary
condition at a fluid-porous dividing surface: derivation of the closure
problem, Chemical Engineering Science 62 (2007) 4025e4039.
[28] M. Chandesris, D. Jamet, Jump conditions and surface-excess quantities at
a fluid/porous interface: a multiscale approach, Transport in Porous Media 78
(2008) 419e438.
[29] B. Alazmi, K. Vafai, Analysis of fluid flow and heat transfer interfacial conditions between a porous medium and a fluid layer, International Journal of
Heat and Mass Transfer 44 (2001) 1735e1749.
[30] K. Vafai, R. Thiyagaraja, Analysis of flow and heat transfer at the interface
region of a porous medium, International Journal of Heat and Mass Transfer
30 (1987) 1391e1405.
[31] K. Vafai, S.J. Kim, Fluid mechanics of the interface region between a porous
medium and a fluid layer e an exact solution, International Journal of Heat
and Fluid Flow 11 (3) (1990) 254e256.
[32] J.A. Ochoa-Tapia, S. Whitaker, Momentum transfer at the boundary between
a porous medium and a homogeneous fluid-II. Comparison with experiment, International Journal of Heat and Mass Transfer 38 (14) (1995)
2647e2655.
[33] A.V. Kuznetsov, Analytical investigation of the fluid flow in the interface
region between a porous medium and a clear fluid in channels partially filled
with a porous medium, Applied Scientific Research 56 (1996) 53e67.
[34] R.A. Silva, M.J.S. Lemos, Turbulent flow in a channel occupied by a porous
layer considering the stress jump at the interface, International Journal of
Heat and Mass Transfer 46 (2003) 5113e5121.
[35] A.V. Kuznetsov, Analytical study of fluid flow and heat transfer during forced
convection in a composite channel partially filled with a BrinkmaneForchheimer porous medium, Flow, Turbulence and Combustion 60
(1998) 173e192.
[36] M.A. Al-Nimr, M.K. Alkam, Unsteady non-darcian fluid flow in parallel plates
channels partially filled with porous materials, Heat and Mass Transfer 33 (4)
(1998) 315e318.
[37] D. Jamet, M. Chandesris, On the intrinsic nature of jump coefficients at the
interface between a porous medium and a free fluid region, International
Journal of Heat and Mass Transfer 52 (2009) 289e300.
[38] S. Veran, Y. Aspa, M. Quintard, Effective boundary conditions for rough
reactive walls in laminar boundary layers, International Journal of Heat and
Mass Transfer 52 (2009) 3712e3725.
[39] D. Jamet, M. Chandesris, B. Goyeau, On the equivalence of the discontinuous
one- and two-domain approaches for modeling of transport phenomena at
a fluid-porous interface, Transport in Porous Media 78 (2009) 403e418.
[40] C.G. Aguilar-Madera, F.J. Valdés-Parada, B. Goyeau, J.A. Ochoa-Tapia, Onedomain approach for heat transfer between a porous medium and a fluid,
International Journal of Heat and Mass Transfer 54 (2011) 2089e2099.
[41] S. Whitaker, The Method of Volume Averaging. Kluwer Academic Publishers,
Netherlands, 1999.
[42] M. Quintard, M. Kaviany, S. Whitaker, Two-medium treatment of heat transfer
in porous media: numerical results for effective properties, Advances in Water
Resources 20 (2e3) (1997) 77e94.
[43] F.J. Valdés-Parada, J. Alvarez-Ramírez, J.A. Ochoa-Tapia, Analysis of mass
transport and reaction problems using Green’s functions, Revista mexicana de
ingeniería química 6 (3) (2007) 283e294.
[44] R. Haberman, Applied Partial Differential Equations with Fourier Series and
Boundary Value Problems, fourth ed. Pearson Prentice Hall, New Jersey, 2004.
[45] F.J. Valdés-Parada, A.M. Sales-Cruz, J.A. Ochoa-Tapia, J. Alvarez-Ramírez, An
integral equation formulation for solving reaction-diffusion-convection
boundary-value problems, International Journal of Chemical Reactor Engineering 6 (2008) 1e20.
[46] F.J. Valdés-Parada, M. Sales-Cruz, J.A. Ochoa-Tapia, J. Alvarez-Ramírez, On
Green’s functions methods to solve nonlinear reaction-diffusion systems,
Computers & Chemical Engineering 32 (2008) 503e511.
[47] P.V. Danckwerts, Continuous flow systems. Distribution of residence times,
Chemical Engineering Science 2 (1986) 219e230.
1368
C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368
[48] J.A. Ochoa-Tapia, P. Stroeve, S. Whitaker, Diffusive transport in two
phase-media: spatially periodic models and Maxwell’s theory for isotropic
and anisotropic systems, Chemical Engineering Science 49 (1994)
709e726.
[49] B.D. Wood, The role of scaling laws in upscaling, Advances in Water Resources
32 (2009) 723e736.
[50] M. de Lemos, Turbulence in Porous Media: Modeling and Applications.
Elsevier Science, 2006.
Fly UP