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 ﬁlled 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 ﬁlled with a porous insert. In order to avoid specifying the boundary conditions at the ﬂuideporous boundary, we solve equations involving position-dependent coefﬁcients (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 ﬁnite 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 ﬂuid. 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 modiﬁed 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 , heterogeneous reactors , ﬁltering, cooling and heating processes , geothermal energy management , 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 ﬂow in concentric pipes  and in circular channels  partially ﬁlled 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 efﬁciency can increase for 3e32% in a modiﬁed 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 . In addition, it was reported that a critical value of the porous insert thickness can be reached for which there is no signiﬁcant improvement in the thermal performance. Recent studies have been performed in 2D parallel-plate channels  and in 3D conﬁgurations . These works evidence that convection is not the only heat transport mechanism that is improved. In fact, Chen and Sutton  demonstrated that the heat transfer by radiation can be augmented about 105% by using ceramic porous inserts in circular ducts. For turbulent ﬂow, Yang and Hwnag  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  carried out a detailed analysis about the conditions that favor the thermal efﬁciency 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  experimentally investigated the effect of metallic porous materials on the heat transfer rate and the pressure drop. Recently, Huang et al.  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 ﬂuidesolid interface, dimensionless closure variable associated to the ﬂuid phase, m closure variable associated to the solid phase, m ﬁt parameter in Eq. (4), i ¼ 0,...,3, dimensionless mass fraction weighted heat capacity, J/(Kg K) speciﬁc heat capacity at pressure constant of the ﬂuid, 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 ﬂuid, 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 ﬂuid 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 ﬂow 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 ﬁt parameter in Eq. (4), i ¼ 0,...,3, dimensionless surface average temperature, K core of a tube. They studied the ﬂow 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 ﬂuideporous dividing surface. Commonly, the boundary conditions applying at the ﬂuideporous interface are postulated from the macroscale (normally imposing continuity conditions), which can lead to physical inconsistencies with the phenomenon taking place at the microscale . 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. , modeling of transport phenomena between a porous medium and a ﬂuid can be carried out by the following alternatives: 1. The one-domain approach (ODA). The whole system, i.e., the free ﬂuid 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 ﬂuid 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 ﬂuid 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 ﬂuid phase, m2/s thermal diffusivity of the solid phase, m2/s slip coefﬁcient, dimensionless ﬂuid 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 coefﬁcient for heat ﬂux 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 ﬂuid region, dimensionless mb ﬂuid viscosity, Pa s spatial average density, Kg/m3 hri rb ﬂuid 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 coefﬁcients exhibit spatial variations is called the inter-region. In this way, the transport models consist of effective-medium equations involving position-dependent coefﬁcients. 2. The two-domain approach (TDA). Here the free ﬂuid 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 , 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 . Furthermore, according to Valdés-Parada et al.  in order to compute the associated jump coefﬁcients it is necessary to account for the spatial variations of the effective transport coefﬁcients 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  about ﬂow in channel C.G. Aguilar-Madera et al. / International Journal of Thermal Sciences 50 (2011) 1355e1368 partially ﬁlled with a porous medium, several efforts have been devoted to this subject [21,22,24e28]. For the same conﬁguration used by Beavers and Joseph, Alazmi and Vafai  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 ﬂux and/ or the ﬁelds. In order to clearly point out the variety of reported jump conditions, we reproduce the information originally reported by Alazmi and Vafai  concerning the boundary conditions for the heat transport problem in Table 1. There aT y f are unknown jump coefﬁcients, l is a characteristic length, and hTiu and hTih are the volume-averaged temperatures corresponding to the porous medium and the free ﬂuid, respectively. Within the TDA there are some analytical solutions reported in the literature. Vafai and Thiyagaraja  derived an approximate analytical expression, based on the matched asymptotic expansion method, for ﬂow and heat transfer through three types of interface regions of a porous medium: the porouseporous interface, the ﬂuideporous interface and the impermeable walleporous interface. Subsequently, Vafai and Kim  presented the exact solution for momentum transport in the ﬂuideporous boundary. Assuming continuity of the velocity and using the stress jump boundary condition derived by Ochoa-Tapia and Whitaker [21,32], Kuznetsov  reported the analytical solution for the velocity ﬁeld, whereas Silva and Lemos  numerically investigated the ﬂuid mechanics under turbulent conditions. Concerning the forced convection heat transfer problem, Kuznetsov  utilized boundary layer solutions and continuity of temperatures and heat ﬂuxes to develop the corresponding analytical solutions. Also, analytical solutions for unsteady ﬂuid ﬂow in parallel-plates channels were derived by AlNimr and Alkam  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 satisﬁed 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 ﬁx 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 (ﬂuid 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-coefﬁcients 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 , which does not necessarily imply the use of continuity conditions in the equivalent TDA formulation . In this regarding, Jamet et al.  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  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 ﬁxed 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 coefﬁcients). This equation was derived by discarding the inertial terms using the volume averaging method by Ochoa-Tapia and Whitaker  and more recently by Valdés-Parada et al.  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 ﬂuid 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 deﬁned everywhere in the system, i.e., in the homogeneous and heterogeneous regions. In the context of the method of volume averaging , the predictions of effective-medium coefﬁcients 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 coefﬁcients 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 ﬂuideporous 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 coefﬁcient. 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  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 speciﬁc heat capacities for the ﬂuid and the solid, respectively; in addition K*eff is the total thermal dispersion tensor introduced by Quintard and Whitaker . 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 . 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 ﬁeld. 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 coefﬁcients 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 coefﬁcients 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.  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  (see Fig. 2). Notice that in Eq. (3), the coordinate axis is located at the ﬂuideporous medium boundary. Recently, Valdés-Parada et al.  computed the spatial dependence of the permeability for a periodic arrangement of nontouching squared cylinders. This situation corresponds to solving the closure problem that deﬁnes the permeability (see Appendix A in ) 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.  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 ﬁtted (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 coefﬁcients in Eq. (4). Fig. 2. Unit cells for predicting the effective coefﬁcients 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 coefﬁcients 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 proﬁles exhibit a similar dependence as the one shown in Fig. A-2 in ref. . 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  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 ﬁtted (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 conﬁguration suggests imposing the following simpliﬁcations ^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 ﬁeld. where u Taking these simpliﬁcations 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 ﬁnd 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 ﬁnd an implicit analytical solution using integral equation formulations in terms of Green’s functions as described in ref. . 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 deﬁne the linear operator Lð Þ ¼ d2 ^ dy 2 ðÞ (14) that F can be simply modiﬁed in order to explore different alternatives of the spatial dependency of the coefﬁcients. 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 ﬁeld 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  ^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 ﬁltration 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 ﬁeld solution. Therefore, it is necessary to implement an iterative scheme to ﬁnd the velocity ﬁeld. 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 inﬂuence of the heterogeneous term at the position y0 over the velocity ﬁeld 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 difﬁculty appearing, is related to the stability of the iterative scheme used to ﬁnd 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 , 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 sufﬁciently high ﬂow 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  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 ﬂuid 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 proﬁles 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 afﬁrm 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. , an analytical expression for this coefﬁcient 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 ﬁelds, 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 inﬂuences 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 ﬁxed 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 . 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 ﬁeld, 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 difﬁculties related to the dispersive term in order to ﬁnd 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 ﬁlled with the porous insert. In this way, we have not ﬁxed 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 ﬁnd 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 inﬂuenced 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 proﬁle is similar to that one found between two parallel plates (i.e., Poiseuille ﬂow). 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 deﬁnition adopted in Eq. (6). Furthermore, in Fig. 5b, we show the importance of the unit cell dimensionality in the predictions of the velocity proﬁles. We observe that the largest differences take place at the peak of the ^ b;max ). In fact, for the case under velocity proﬁles (i.e., u b Fig. 6. Velocity proﬁles 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 inﬂuence 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 , the velocity proﬁles are also inﬂuenced by the conﬁguration of the channel partially ﬁlled with the porous medium. In Fig. 6a we present the velocity proﬁles 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 . 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 sufﬁcient 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 inﬂuence of the ratio x/H over the velocity proﬁles in Fig. 6b. We observe that the proﬁles 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 ﬁnalize this section, as a quantitative measurement about the inﬂuence of the porous insert over the hydrodynamic performance of the channel, let us introduce the average volumetric ﬂow 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 ﬂow rate is divided with the one corresponding when there is no porous insert within the channel. We observe an exponential decrease of the ﬂow rate as long as the porous insert ﬁlls 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 ﬂow rate directly from Darcy’s law when the porous insert completely saturates the channel. In addition, in Fig. 7b the inﬂuence 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 ﬁlled 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 ﬂow rate moderately decreases. 5.2. Thermal performance With the velocity ﬁeld satisfying the iterative scheme of Eq. (23), this is supplied as input data for the heat transport problem. This problem was solved using the ﬁnite element solver Comsol Multiphysics 3.5a, which uses adaptive mesh reﬁnements 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 difﬁculty in the numerical solution lies on the correct quantiﬁcation 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 ﬁelds, deﬁned 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 ﬁeld 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 ﬂow rates (i.e. PeH < 1), practically all the ﬂuid 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 ﬂuid 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 ﬂuideporous medium boundary, especially for PeH > 1. In addition, under these conditions, the ﬂuid 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 ﬂuid 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 ﬂuid (and the opposite is true). Notice that for k ¼ 0.01 the ﬂuid does not reach the temperature of the walls at the exit channel, whereas for k ¼ 100 practically all the ﬂuid 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 ﬂuid and the solid composing the substrate. Regarding the inﬂuence of porosity, in Fig. 10 we show the channel temperature proﬁles 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 ﬁlled with the porous insert was utilized. For low values of the porosity, 3b,u < 0.3, the temperature proﬁle 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 proﬁles 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 inﬂuence of the size of the porous insert is presented in Fig. 11. As more substrate ﬁlls the channel, the dimensionless temperature proﬁles are more uniform and close to 1. It is interesting to note that only the ﬂuid 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 ﬁeld 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 ﬁeld 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 ﬁeld is the typical one encountered in most heat exchangers. Even though the usage of porous substrates leads to increasing pressure drops , 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 proﬁles 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 deﬁned 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 ﬁelds shown in Fig. 9. The inﬂuence 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 . However, a thorough investigation of turbulence is beyond the scope of this work. In Fig. 13a we present the effect of the ﬂow 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 ﬂuid 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 ﬂuid 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 ﬁeld 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. Inﬂuence 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. Inﬂuence 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. ﬂuid 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 conﬁguration 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 ﬁlled 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 ﬁlled with the porous medium or the ﬂuid, 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 ﬂux 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 ﬂuid phase. In this way, the portion of the channel that is only occupied by the ﬂuid 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 sufﬁciently 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 ﬁlled 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 ﬁlled 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 coefﬁcients undergo rapid changes near the inter-region separating the homogeneous porous insert and the zone of free ﬂuid. With this methodology, we avoid specifying the boundary conditions at the ﬂuideporous 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 ﬁeld. In this way, an iterative scheme was implemented in order to set the velocity ﬁeld within the channel. We found that signiﬁcant differences in the velocity proﬁles 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 proﬁles and the volumetric ﬂow 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 ﬁelds and the local Nusselt number by varying the ﬂow 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 ﬂuid 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 ﬁlled 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 coefﬁcients 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 coefﬁcients. Nevertheless, the predictions from the TDA are driven by the jump boundary conditions between the porous medium and the ﬂuid. As mentioned in the introduction, in order to compute the coefﬁcients involved in the jump conditions, it is necessary to account for the spatial variations of the effective transport coefﬁcients 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 deﬁnition arising from using the method of volume averaging , 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 ﬂuid to the solid (nbs ¼ nsb), V is the volume of the averaging domain, Abs represents the spatial location of the ﬂuidesolid interface, Vb is the portion of the averaging domain occupied by ~ represents the spatial deviations of the velocity. ﬂuid and v b The closure variable maps the gradient of the average temperature ﬁeld onto the spatial deviations of the ﬂuid temperature ﬁeld 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 ﬂuid 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 ﬁeld onto the spatial deviations of the point-wise solid temperature ﬁeld according to the expression, Appendix A. Closure problem for heat transport T~ s ¼ bs $VhTi In this section, we present the deﬁnition 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 deﬁned 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 satisﬁed 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 . References  M.K. Alkam, M.A. Al-Nimr, Transient non-darcian forced convection ﬂow in a pipe partially ﬁlled with a porous material, International Journal of Heat and Mass Transfer 41 (2) (1998) 347e356.  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.  M.K. Alkam, M.A. Al-Nimr, Solar collector with tubes partially ﬁlled with porous substrate, Journal of Solar Energy Engineering 121 (1) (1999) 20e24.  M.A. Al-Nimr, M.K. Alkam, Unsteady non-darcian forced convection analysis in an annulus partially ﬁlled with a porous material, Journal of Heat Transfer 119 (4) (1997) 799e804.  S.Y. Byun, S.T. Ro, J.Y. Shin, Y.S. Son, D.Y. Lee, Transient thermal behavior of porous media under oscillating ﬂow condition, International Journal of Heat and Mass Transfer 49 (2006) 5081e5085.  Z.F. Huang, A. Nakayama, K. Yang, C. Yang, W. Liu, Enhancing heat transfer in the core ﬂow by using porous medium inserts in a tube, International Journal of Heat and Mass Transfer 53 (2010) 1164e1174.  K.N. Shukla, Diffusion Processes During Drying of Solids, Vol. 11 of Series in Theoretical and Applied Mechanics. World Scientiﬁc Publishing Co., 1999.  O. Levenspiel, Chemical Reaction Engineering, third ed. John Wiley & Sons, 1999.  W.L. McCabe, J.C. Smith, P. Harriot, Unit Operations of Chemical Engineering, seventh ed. McGraw Hill, 2005.  M.H. Dickson, M. Fanelli, Geothermal Energy: Utilization and Technology. Earthscan Publications, 2005.  M.A. Al-Nimr, M.K. Alkam, A modiﬁed tubeless solar collector partially ﬁlled with porous substrate, Renewable Energy 13 (2) (1998) 165e173.  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.  T.C. Jen, T.Z. Yan, Developing ﬂuid ﬂow and heat transfer in a channel partially ﬁlled with porous medium, International Journal of Heat and Mass Transfer 48 (2005) 3995e4009.  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.  Y.T. Yang, M.L. Hwang, Numerical simulation of turbulent ﬂuid ﬂow and heat transfer characteristics in heat exchangers ﬁtted with porous media, International Journal of Heat and Mass Transfer 52 (2009) 2956e2965.  A.P. Lukisha, V.F. Prisnyakov, The efﬁciency of round channels ﬁtted with porous, highly heat-conducting insert in a laminar ﬂuid coolant ﬂow at boundary conditions of the third kind, International Journal of Heat and Mass Transfer 53 (2010) 2469e2476.  B.I. Pavel, A.A. Mohamad, An experimental and numerical study on heat transfer enhancement for gas heat exchangers ﬁtted with porous media, International Journal of Heat and Mass Transfer 47 (2004) 4939e4952.  M. Prat, On the boundary conditions at the macroscopic level, Transport in Porous Media 4 (1989) 259e280.  B. Goyeau, D. Lhuillier, D. Gobin, M.G. Velarde, Momentum transport at a ﬂuid-porous interface, International Journal of Heat and Mass Transfer 46 (2003) 4071e4081. 1367  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.  J.A. Ochoa-Tapia, S. Whitaker, Momentum transfer at the boundary between a porous medium and a homogeneous ﬂuid-I e theoretical development, International Journal of Heat and Mass Transfer 38 (14) (1995) 2635e2646.  F.J. Valdés-Parada, J. Alvarez-Ramírez, B. Goyeau, J.A. Ochoa-Tapia, Computation of jump coefﬁcients for momentum transfer between a porous medium and a ﬂuid using a closed generalized transfer equation, Transport in Porous Media 78 (2009) 439e457.  G.S. Beavers, D.D. Joseph, Boundary conditions at a naturally permeable wall, Journal of Fluid Mechanics 30 (1967) 197e207.  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.  J.A. Ochoa-Tapia, S. Whitaker, Heat transfer at the boundary between a porous medium and a homogeneous ﬂuid: the one-equation model, Journal of Porous Media 1 (1) (1998) 31e46.  B.D. Wood, M. Quintard, S. Whitaker, Jump conditions at a non-uniform boundaries: the catalytic surface, Chemical Engineering Science 55 (2000) 5231e5245.  F.J. Valdés-Parada, B. Goyeau, J.A. Ochoa-Tapia, Jump momentum boundary condition at a ﬂuid-porous dividing surface: derivation of the closure problem, Chemical Engineering Science 62 (2007) 4025e4039.  M. Chandesris, D. Jamet, Jump conditions and surface-excess quantities at a ﬂuid/porous interface: a multiscale approach, Transport in Porous Media 78 (2008) 419e438.  B. Alazmi, K. Vafai, Analysis of ﬂuid ﬂow and heat transfer interfacial conditions between a porous medium and a ﬂuid layer, International Journal of Heat and Mass Transfer 44 (2001) 1735e1749.  K. Vafai, R. Thiyagaraja, Analysis of ﬂow and heat transfer at the interface region of a porous medium, International Journal of Heat and Mass Transfer 30 (1987) 1391e1405.  K. Vafai, S.J. Kim, Fluid mechanics of the interface region between a porous medium and a ﬂuid layer e an exact solution, International Journal of Heat and Fluid Flow 11 (3) (1990) 254e256.  J.A. Ochoa-Tapia, S. Whitaker, Momentum transfer at the boundary between a porous medium and a homogeneous ﬂuid-II. Comparison with experiment, International Journal of Heat and Mass Transfer 38 (14) (1995) 2647e2655.  A.V. Kuznetsov, Analytical investigation of the ﬂuid ﬂow in the interface region between a porous medium and a clear ﬂuid in channels partially ﬁlled with a porous medium, Applied Scientiﬁc Research 56 (1996) 53e67.  R.A. Silva, M.J.S. Lemos, Turbulent ﬂow 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.  A.V. Kuznetsov, Analytical study of ﬂuid ﬂow and heat transfer during forced convection in a composite channel partially ﬁlled with a BrinkmaneForchheimer porous medium, Flow, Turbulence and Combustion 60 (1998) 173e192.  M.A. Al-Nimr, M.K. Alkam, Unsteady non-darcian ﬂuid ﬂow in parallel plates channels partially ﬁlled with porous materials, Heat and Mass Transfer 33 (4) (1998) 315e318.  D. Jamet, M. Chandesris, On the intrinsic nature of jump coefﬁcients at the interface between a porous medium and a free ﬂuid region, International Journal of Heat and Mass Transfer 52 (2009) 289e300.  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.  D. Jamet, M. Chandesris, B. Goyeau, On the equivalence of the discontinuous one- and two-domain approaches for modeling of transport phenomena at a ﬂuid-porous interface, Transport in Porous Media 78 (2009) 403e418.  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 ﬂuid, International Journal of Heat and Mass Transfer 54 (2011) 2089e2099.  S. Whitaker, The Method of Volume Averaging. Kluwer Academic Publishers, Netherlands, 1999.  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.  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.  R. Haberman, Applied Partial Differential Equations with Fourier Series and Boundary Value Problems, fourth ed. Pearson Prentice Hall, New Jersey, 2004.  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.  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.  P.V. Danckwerts, Continuous ﬂow 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  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.  B.D. Wood, The role of scaling laws in upscaling, Advances in Water Resources 32 (2009) 723e736.  M. de Lemos, Turbulence in Porous Media: Modeling and Applications. Elsevier Science, 2006.