A two-field model for natural convection in a porous annulus at high Rayleigh numbers

PII:

Chemical Engineering Science, Vol. 53, No. 8, pp. 1505—1516, 1998 ( 1998 Published by Elsevier Science Ltd. All rights reserved Printed in Great Britain S0009–2509(97)00421–1 0009—2509/98 $19.00#0.00

A two-field model for natural convection in a porous annulus at high Rayleigh numbers Julio A. Deiber* and Rau´l A. Bortolozzi Instituto de Desarrollo Tecnolo´gico para la Industria Quı´ mica, INTEC, Gu¨emes 3450, 3000 Santa Fe, Argentina (Received 13 March 1997; accepted 12 November 1997) Abstract—This work studies the natural convection in a vertical porous annulus for Rayleigh numbers up to 8]104. A model is deduced within a simple context of the mixture theory in which the following effects are included: (1) The solid and the fluid may be at different temperatures depending upon the flow regime and the permeability of the porous matrix. (2) Density, heat conductivity and viscosity of the fluid depend on temperature. (3) The porosity of the rigid solid matrix varies near the cavity walls. (4) Heat conduction in the fluid phase is accompanied with thermal dispersion. (5) Forchheimer inertial, Brinkman viscous and convective terms are considered in the momentum balance of the fluid phase. The resulting model is used to evaluate the heat transfer in a porous medium, the matrix of which is composed of equal spheres of well-specified diameters. Inner and outer walls of concentric cylinders are maintained at different temperatures and the two horizontal walls are adiabiatic. Numerical predictions of the Nusselt number are compared with available experimental data. For high Rayleigh and Darcy numbers (porous media of high permeability) all these effects have to be considered in the mathematical model. Thermal dispersion has a significant effect on the heat transfer through the solid—fluid composite. ( 1998 Published by Elsevier Science Ltd. All rights reserved. Keywords: Natural convection; porous media; thermal dispersion; non-Darcian effects; twofield model; mixture theory.

INTRODUCTION

Transport phenomena in porous media is a subject of wide technological interest. It has many applications which may range from natural underground reservoirs and biological systems to man-made packed-bed reactors, where chemical and physical transformations take place. In particular, natural convection in porous media is a complex phenomenon that involves momentum and heat transport in two irregular phases; one is the fixed solid matrix, in general assumed rigid, and the other is the fluid — usually Newtonian — flowing through complex interstices. At present, most of the transport mechanisms associated to this problem have been described and discussed in the literature. Therefore, in the last years, special emphasis has been placed on the generation of a full model that put together all the phenomena that affect the heat transfer through the solid—fluid composite. In this sense, the works of Amiri and Vafai (1994), and Vafai and Kim (1995) are relevant. This work studies the natural convection in a vertical porous annulus (see Fig. 1) mainly at high Rayleigh numbers (greater than 103), where most of the possible mechanisms of momentum and heat

* Corresponding author.

transfer are present. To classify the flow regimes of natural convection in this particular system, the values taken by the Rayleigh number, defined here as Ra"oo gbK D*¹/ko ao , are required. In this f,r = f, r m, r number *¹"¹ !¹ is the difference between hot h c (¹ ) and cold (¹ ) wall temperatures, oo and ko are h c f, r f, r the density and the viscosity of the fluid at the reference temperature ¹ "(¹ #¹ )/2, respectively; ao r h c m, r is the thermal diffusivity of the solid—fluid system, b is the fluid isobaric thermal expansion coefficient, K is = the permeability of the porous matrix far from impermeable walls and D is the width of the porous cavity. For low values of Ra (less than 500), the model that describes this phenomenon is basically the classical Darcy equation for momentum transport and the energy balance involving convection and conduction of heat (see, for example, Prasad and Kulacki, 1985). Unfortunately, this simple model is not appropriate at high Rayleigh numbers, because additional mechanisms of momentum and heat transport are enhanced to a point that they cannot be ignored in the calculations (Lauriat and Prasad, 1989). This has been evident when experimental data are compared with typical models of natural convection in porous media (see, for example, Prasad et al., 1985; Kladias and Prasad, 1991). Still more interesting is the fact that these mechanisms may compensate each other when

1505

1506

J. A. Deiber and R. A. Bortolozzi

Fig. 1. Vertical annulus, coordinate system and scales.

their effects on the Nusselt number Nu"hD/ko are m, r evaluated. Here, h is the average heat-transfer coefficient and ko is the stagnant thermal conductivity of m, r the saturated porous medium. Therefore, at high Ra, one may consider the effect of variable porosity near the impermeable walls to improve the prediction of the theoretical Nusselt number in relation to the experimental one (Vafai, 1984; Hong et al., 1987; David et al., 1989). In addition, the Forchheimer inertial term (see, for example, Vafai and Tien, 1981; Poulikakos and Bejan, 1985) and the Brinkman viscous term (Vafai and Tien, 1981; Lauriat and Prasad, 1987) are relevant. The effect of thermal dispersion has been considered by Georgiadis and Catton (1988) (see also Combarnous and Bories, 1975). In this way, one gets several mechanisms in the full model with some undetermined parameters and with a fair fitting curve of experimental results at high Rayleigh numbers reported as Nu vs Ra for different Darcy numbers Da"K /D2. = The purpose of this work has two clear targets. One is concerned with the proposal of a model for natural convection in porous media, where most of the mechanisms of momentum and heat transfer between solid and fluid are included and evaluated. Only at relatively high Ra these mechanisms can be detected, and hence, discussed separately. Otherwise, at low Ra, most of them become negligible and are screened by the Darcian term. The other aspect is concerned with providing a more rigorous model to many practical applications, where high Ra can be obtained as, for example, in packed-sphere beds, grain storage systems, high permeability geothermal and oil reservoirs, building thermal insulation, energy storage units and solar receiver devices. The model presented here has been deduced within a simple context of the mixture theory (see, for example, Truesdell, 1969; Bowen, 1976) where the following effects are included: (1) The solid and the fluid may be at different temperatures depending upon the flow regime and the permeability of the porous matrix, i.e. the values of Rayleigh and Darcy numbers. (2) Density, heat conductivity and viscosity of the fluid depend on temperature. (3) The porosity of the rigid solid matrix varies with the perpendicular distance from impermeable walls. (4) Heat conduction

in the fluid phase is accompanied with thermal dispersion. (5) Forchheimer inertial, Brinkman viscous and convection terms are considered in the momentum balance of the fluid phase. It is clear that a theoretical framework shall be used, like the mixture theory, for the incorporation of all these effects into a mathematical model, without any arbitrary consideration, mainly in what concerns to spatial variations and gradients of the dependent variables. The resulting model is applied to natural convection in a porous medium, the matrix of which is composed of equal spheres of well-specified diameters and is confined between two vertical concentric cylinders. The inner one is being heated at ¹ and the outer h is at a colder temperature ¹ . Upper and bottom c surfaces are adiabatic (see Fig. 1). This geometry is chosen because there exist in the literature well documented experimental values of Nu as function of Ra for different sphere diameters (Prasad et al., 1985; David et al., 1989). This allows one a confrontation of numerical predictions with experimental points as well as a rigorous validation of the emerging theories on the subject. Therefore, the proposed model framed in the mixture theory is presented next. A section considering the numerical method is also addressed here, because the quality of numerical results may be affected when momentum and thermal boundary layers appear at high Rayleigh numbers. It is thus found necessary to evaluate the Nusselt number at both constant temperature walls (¹ and ¹ ) of the h c porous cavity to control the consistency of numerical calculations. The Results and Discussion section compares the theoretical predictions with experimental data expressed as Nu vs Ra, while Da is parametrically fixed. Therefore, the experiments considered are those of the packed beds studied and reported by Prasad et al. (1985) and David et al. (1989), where the spheres are 3.0, 6.0 and 22.3 mm in diameter. The convection cell described above has an aspect ratio ¸/D"1. In this section, the validity of the Boussinesq approximation for the experimental data used is also analyzed. FORMULATION OF THE MODEL WITHIN THE MIXTURE THEORY

In the context of the mixture theory, N species are considered, the properties of which are designated with a subscript a"1 . . . N (Truesdell, 1969; Bowen, 1976). Therefore, mass, momentum and internal energy of species satisfy local balances as follows: d a o #o $ ) v "0 a a dt a

(1)

d a v "!$ ) T #o g#m a dt a a a a

(2)

d a º "!$ ) q !T : $v !m ) v #et a dt a a a a a a a

(3)

o o

where the species velocity v , internal energy º and a a stress tensor T are considered. This last tensor is a assumed to be symmetric here. The only volumetric

A two-field model at high Rayleigh numbers

1507

Table 1. Parameters of the water—glass system

Da Fs Pr* f j* l* c e = A B

d "3]10~3 m p

d "6]10~3 m p

d "22.3]10~3 m p

3.36]10~7 3.73]10~4 4.11—5.46 0.79—0.81 1.10—1.11 0.0242 0.351 0.25 7.50

1.66]10~6 7.68]10~4 4.37—6.17 0.79—0.81 1.11—1.12 0.0484 0.370 0.25 7.50

4.42]10~5 3.16]10~3 4.20—6.01 0.87—0.88 1.10—1.12 0.1796 0.431 0.30 7.50

* Since Pr , j and l depend on the reference temperature, ranges of variation are f indicated in each case.

force is gravity g. In eqs (1)—(3), (d /dt)( ) )" a (L/Lt)( ) )#v ) $( ) ) is the species time derivative. This a theory considers the exchanges of momentum m and a total energy (kinetic plus internal energies) et between a species, with the constraints + N m "0 and a/1 a + N et "0, so that, the classical local equations of a/1 a the mixture are recovered when eqs (1)—(3) are summed on a. For the purposes here, we point out that o"+ o and v"+ o v /o. Further details are dea a a a a scribed in the original works. In eq. (3), q is the flux by a heat conduction of species a. In writing eqs (1)—(3), it is considered that chemical reaction and volumetric heat transfer do not exist in the system. These equations can be applied to the flow and heat transfer of Newtonian fluid in an undeformable porous medium by taking a"f, s, where f indicates fluid and s stands for solid. Thus, o "oo e and f f o "oo (1!e). Throughout this work superscript (o) s s indicates that properties correspond to pure species. Several additional definitions are next required. First, density variations are accounted for in the model (Boussinesq approximation is not considered) by using the following expression for the fluid density oo"oo [1!b(¹ !¹ )]. Since the fluid is f f,r f r Newtonian, the stress tensor is T "po ed!eko f f f ($v #$v T ). The porosity is allowed to vary near the f f cavity walls according to (Benenati and Brosilow, 1962) e"e

=

C

A BD

Bn 1#A exp ! d p

(4)

where n is the perpendicular distance from any wall of the porous cavity. Following Vafai (1984), in eq. (4) we take into account only the exponential decay of porosity from the wall, neglecting any spatial oscillations, which are considered to be a secondary effect. A and B are empirical constants and the values assigned to them in this work are reported in Table 1. In this aspect, it is appropriate to mention here that the quantification of parameters A and B is comprised in a relatively wide range of variation, according to data obtained from the literature (0.2—1.7 for A and 2—8 for B) (Vafai, 1984; Cheng and Hsu, 1986; Hong et al.,

1987; David et al., 1989; Nield and Bejan, 1991; Kladias and Prasad, 1991; Chen et al., 1992; Amiri and Vafai, 1994). The reason for finding this wide range of variation of parameters A and B is that they have been obtained by comparing experimental data with different models that include variable porosity. In fact, it is in general found that the consideration of the effect of isotropic thermal dispersion in the energy balance yields lower values of A (see, for example, Kladias and Prasad, 1991; Amiri and Vafai, 1994). From the numerical calculations in this work (see Table 1 and Results and Discussion section) we find then enough evidence to sustain this conclusion. Also, ko depends on temperature according to, f ko"ko exp[!C (¹ !¹ )] (5) f f,r 1 f r The mixture heat fluxes are defined as q "!ko (1!e)$¹ (6) s s s q "!(ko ed#k ) ) $¹ (7) f f d f where ko and ko are the thermal conductivities of solid s f and fluid, respectively, while k is the thermal disperd sion tensor generated by the fluid convective fluctuations in the interstices. This tensor is expressed (see Mercer et al., 1982; Amiri and Vafai, 1994; Howle and Georgiadis, 1994) as e(v v ) (8) k "oo Co d l(n)a : f f . d f f p Dv D f Since the porous media is considered isotropic, the thermal hydrodynamic dispersivity tensor a has symmetric properties as described by Scheidegger (1974); thus, a "a , a "a , a "a "1 (a !a ), for iiii l iijj t ijij ijji 2 l t all possible subscript permutations when i"r, z and j"r, z in a cylindrical coordinate system (see Fig. 1). Otherwise, components of a are zero. In eq. (8), Co is f the heat capacity of the fluid. Also, l(n) is the Van Driest function (Cheng and Hsu, 1986; Cheng and Zhu, 1987) which considers the damping of fluid fluctuations near the walls and is expressed as l(n)"1!exp

A B n ud

p

(9)

1508

J. A. Deiber and R. A. Bortolozzi

where u is an empirical constant taking values between 1 and 2.5. The fluid conductivity ko depends on temperature f according to ko "ko [1#C (¹ !¹ )]. (10) f f, r 2 f r The internal energy exchange between species, which is e "et !m ) v , is defined (see, for example, a a a a Combarnous and Bories, 1975; Amiri and Vafai, 1994) as e "!e "h (v )a (¹ !¹ ) (11) s f sf f v f s where a is the specific surface of the porous medium v and is evaluated from, a "6(1!e)/d . This last exv p pression is deduced from geometric considerations (Bird et al., 1960). In eq. (11) the heat-transfer coefficient h depends, sf as a first approximation, on velocity and properties of the fluid, according to (Wakao et al., 1979) ko h (v )" f,r [2#1.1 Pr1@3 Re0>6] (12) sf f f p d p where Re "oo eDv Dd /ko is the particle Reynolds p f,r f p f,r number and Pr "Co ko /ko is the fluid Prandtl f f f,r f,r number. The momentum exchange between species is deduced here on the base of the terms accounted for in previous models of natural convection in porous media (see, for example, Vafai and Tien, 1981; Lauriat and Prasad, 1989). Thus, Darcy and Forchheimer terms are interpreted in this work as a part of the momentum exchange between solid and fluid as follows: boo ko m "!m "! f e2 v ! f e3Dv Dv #po $e f s f f f K K f (13) where the permeability K"e3d2 /180(1!e)2 and the p empirical Forchheimer factor b"1.8d /180(1!e) p have been included (Prasad et al., 1985; David et al., 1989). Combining eqs (2) and (13), and the constitutive equation for T , the momentum balance of the fluid f species is obtained in terms of the dimensionless superficial velocity V"ev D/ao and dimensionless f m, r fluid pressure po"po K /ko ao , where ao " f = f, r m, r m, r ko /oo Co . Additionally, the energy balance results m, r f f from eqs (3), (6), (7) and (11). Therefore, after intensive algebraic work, the following dimensionless model for the steady-state regime is obtained. Continuity: $ ) (/ V)"0. 0 Momentum balance of fluid:

(14)

/ Da 0 V ) $(V/e)"!e$po#(h !")Raek f jPr f / / / Fs ! 1 e 1# 0 3 DVD V / / jPr 2 1 f

A

B

#/ eDaM+2(V/e) 1 #$[$ ) (V/e)]#[$(V/e) #$(V/e) T] ) $ ln(/ e)N. 1 Energy balance of fluid: / V ) $h "$ ) [(je/ d#k* ) ) $h ] 0 f 4 d f !H(h !h ). f s Energy balance of solid:

(15)

(16)

0"$ ) [l(1!e)$h ]#H(h !h ). (17) s f s In eqs (14)—(17), the following dimensionless functions, which depend on position, are defined: oo / " f , 0 oo f, r b / " , 3 b =

ko K / " f , / " , 1 ko 2 K f, r = ko k / " f , k*" d , 4 ko d ko f, r m, r 6(1!e) h D2 sf . H" d ko p m, r The relation Fs"b /D is the Forchheimer num= ber. Also, the constant ratios ""1/b*¹, j"ko / f, r ko and l"ko/ko are defined in eqs (15)—(17), where m, r s m, r ko "e ko #(1!e )ko . Temperatures and coorm, r = f, r = s dinates are scaled with *¹ and D, respectively. NUMERICAL METHOD

Equations (14) and (15) are written using the vorticity-stream function scheme by following a similar procedure presented by Peirotti et al. (1987), where the definition of the stream function includes the dimensionless variable density / so that the continuity 0 equation is satisfied. The resulting model [eqs (14)—(17)] is expressed in finite differences as follows: (1) Convective terms are discretized with the second upwind technique (see Roache, 1972). (2) Central differences are used for second derivatives. (3) Heat conduction terms with variable effective conductivities on the right-hand side of eq. (16) are discretized according to the procedure described by Paeceman (1977). The discrete equations are solved using the relaxation method through successive iterations as it was described in the work of Peirotti et al. (1987). The iteration procedure is continued until a convergence criterion is satisfied as follows: + DSk`1!Sk D

ij

ij

ij )p

(18) ij ij where S is the dependent variable that is being nuij merically calculated (vorticity, stream function and temperatures of fluid and solid). The subscripts i and j imply spatial position in the computational mesh, and k refers to the iteration number. The values of p are assigned between 10~5 and 10~8, as discussed later. Since high values of Ra generate thermal and momentum boundary layers near the cavity walls, it is + DSk D

A two-field model at high Rayleigh numbers

1509

Fig. 2. Numerical Nusselt number Nu as function of grid interval *R near the walls of the porous cavity. h Solid circles (d) are numerical results for equally spaced grids. Open circles (s) are numerical results for non-uniform grids [eqs (23) and (24)]. Dashed line (— — —) is the matching between calculations carried out with and without coordinate transformation.

highly recommended to carry out a coordinate transformation so that grid intervals are very small (of the order of 10~3) near these walls. On the other hand, this transformation yields a coarse mesh at the cavity center. This interplay between mesh sizes and cavity regions requires a careful verification of consistency between numerical evaluations and the macroscopic energy balance. In fact, it is simple to show that the Nusselt numbers

PA

Nu "! h

0 r at R" i D and,

(19)

PA

Nu "! c

1

0 r at R" o D

B

Lh Lh f#l(1!e) s dZ je/ 4 LR LR

1

B

Lh Lh f#l(1!e) s dZ je/ 4 LR LR (20)

which are evaluated at hot and cold walls, respectively, must satisfy the relation Nu "Nu i (21) h c where i"r /r is the radius ratio (see also Prasad and o i Kulacki, 1984). Therefore, in order to control the consistency of computational calculations, a numerical parameter is defined as follows: Nu h . s" (22) Nu i c Thus, the numerical procedure is giving appropriate results when sP1. It is found that the requirement sP1 is a very severe test for the accuracy of results, and this limit is only attained when p and the grid intervals are made small enough for each value of Ra and Da.

In this work, the following coordinate transformation is used (Ka´lnay de Rivas, 1972):

AB AB

nm 1 R" #sin2 (i!1) 2 Z"sin2

ng 2

(23) (24)

where m and g are the new coordinates, so that constant grid intervals *m and *g yield a fine mesh near walls when *R and *Z are considered. A criterion is required in order to decide when the coordinate transformation must be used. With this purpose, a study was carried out for a typical value Ra"300, with a fixed p"10~5. Computations were performed to evaluate Nu for constant mesh sizes (i.e. h using *R"*Z, because the grid is square and the aspect ratio ¸/D"1) first with the grid of 40]40 and ending with a grid of 100]100. These results are shown in Fig. 2, where Nu vs *R is plotted. Since h a stable value of Nu shall be expected as *R deh creases, and a grid higher than 100]100 (*R"0.01) is becoming numerically unsuitable because the computation time in substantially increased, the decrease of *R near walls must be achieved with the coordinate transformation [eqs (23) and (24)] once the uniform grid of 100]100 is reached. Thus, for *R(0.0015 near the wall, it is found that Nu gets a small stable h region as Fig. 2 shows through the squared portion at the upper left. Since this curve is dependent on the average error of the convergence criterion, more details about the influence of p on the final values of Nu h are presented in Fig. 3. It is also appropriate to point out that further decreases of *R, i.e. for *R(10~3, differences between the Nusselt numbers for this typical run do not exceed 0.16% when p"10~6. Therefore, for grids of around 50]50 or greater, and p"10~6, and using the above coordinate transformation, a stable value of Nu can be computed h

1510

J. A. Deiber and R. A. Bortolozzi

Fig. 3. Numerical Nusselt numer Nu as function of grid interval *R near the walls of porous cavity when h the convergence criterion p is varied between 10~5 and 10~8. Open circles are numerical results for non-uniform grids. The connecting lines between points show the trends schematically.

throughout the range of Ra studied here without appreciable differences. RESULTS AND DISCUSSION

Once the thermophysical properties of the fluid and solid are fixed (see Table 1), the mathematical model [eqs (14)—(17)] shows that the Nusselt number can be plotted as function of the Rayleigh number when the Darcy number is kept parametrically constant. From this point onwards, whenever Nusselt number is mentioned, we are referring to the value of this parameter as evaluated at the hot wall, according to eq. (19). The Darcy number can be varied by choosing different particle diameters. Here three cases are modeled, according to data reported previously (see Prasad et al., 1985; David et al., 1989, and the Appendix): Case a: d "3]10~3 m at low Ra and Da" p 0.34]10~6. Case b: d "6]10~3 m at intermediate Ra and p Da"0.17]10~5. Case c: d "22.3]10~3 m at high Ra and Da" p 0.44]10~4. The model is applied to the system described in Fig. 1, which is composed of glass spheres and water as saturating fluid. The thermophysical properties of fluid vary with temperature according to eqs (5) and (10), where C "0.02 and C "0.0015. Also, b is 1 2 taken between 0.00025 and 0.00040 K~1, depending upon the value of the reference temperature of experimental points reported in Table A1. It is appropriate to point out here that most of the results published in the literature concerning this porous cavity presented the numerical evaluation of Nu for Ra(1.7]104. More recently, Marpu (1995) published results of Nu vs Ra by computing a model without variable porosity and without thermal dispersion for 500)Ra)5000. These calculations were done by David et al. (1989) for Ra(1.7]104 includ-

ing variable porosity in their model. Marpu compared his numerical results with experimental data published by Prasad et al. (1985) but, unfortunately, this comparison is only approximate. In fact, Prasad et al. defined effective Nusselt and Rayleigh numbers (Fig. 12 of reference) that are not the same as those of Marpu’s analysis. In the present work, it is found that in case a, numerical results fit well experimental data. Thus, Nu h calculated with the Darcy model (Fs"0, A"B"0, C "0, C "0, a "a "0 and neglecting mo1 2 l t mentum convective terms) and with the model proposed in this work takes almost the same values. In other words, the effects on the Nusselt number of Forchheimer, Brinkman and momentum convective terms as well as thermal dispersion and variable porosity are negligible. This is illustrated in Fig. 4. As reported by David et al. (1989), we also find a lesser agreement at low Rayleigh numbers due to larger errors associated with the experimental data at small temperature differences. For case b, the Darcy model underpredicts the experimental values of Nusselt number, when Rayleigh numbers are greater than 500. When all the effects discussed in this work are included, the numerical predictions are substantially better. This situation is illustrated in Fig. 5. For case c, at high Ra, the Darcy model overpredicts significantly the experimental values of Nusselt number (see Fig. 6, curve 1). When Forchheimer, Brinkman and momentum convective terms are added to the model, the effect is to underpredict the Nusselt number (Fig. 6, curve 2). Furthermore, when the variable porosity is considered in this flow regime [see eq. (4)] the slope of the curve Nu(Ra) calculated numerically is not the same as that obtained experimentally (Fig. 6, curve 3). This difference is enhanced for high Rayleigh numers (Ra'103) where numerical results start to be progressively below the experimental values. In this work, it is found that thermal

A two-field model at high Rayleigh numbers

1511

Fig. 4. Comparison between numerical and experimental Nusselt numbers as function of Rayleigh number for the water—glass system, when d "3]10~3 m. Parameters are presented in Table 1. p

Fig. 5. Comparison between numerical and experimental Nusselt numbers as function of Rayleigh number for the water—glass system, when d "6]10~3 m. Parameters are presented in Table 1. p

dispersion corrects this abnormal prediction of the model in relation to experimental data at high Rayleigh numbers (Fig. 6, curve 4). This aspect of the problem was analyzed previously by Georgiadis and Catton (1988), who proposed an expression for the thermal dispersion terms deduced from the statistical fluctuations of temperature and velocity fields of the fluid in the interstices of the porous matrix. These authors observed that their model overpredicted the values of Nusselt numbers, mainly when the fluid velocity was high. The calculations with the present model show that, when the thermal dispersion [see eq. (8)] is considered to be isotropic, the coefficient a "a l t reported by Georgiadis and Catton (1988) and also by Chen et al. (1992) makes the model overpredicts the Nusselt number for high values of Ra. On the other hand, when the thermal dispersion is considered anisotropic and homogeneous, i.e. a Oa so that l t 0.25(a /a (0.40, and the damping effect of the Van t l

Driest function is also included [see eq. (9)] the numerical prediction of Nusselt numbers is very close to the experimental values. This improvement of the model predictions was also observed by Howle and Georgiadis (1994) when a saturated porous cavity heated from below was studied and results were reported until Ra+500. These authors also assumed the slip velocity boundary condition on the cavity walls, i.e. the Brinkman term was neglected. Figure 7 shows our results for the three cases studied in this work. As Ra is increased, there is almost a matching between the two curves of Nu(Ra) for cases a and b. This is a consequence of the relatively low values of d (3 and 6]10~3 m) and porosities, p giving Darcy numbers less than 10~5. Thus, for Da(10~5, the Nusselt number is almost insensitive to this dimensionless number in the cavity of Fig. 1. However, for higher values of Darcy number (Da'10~5), which in general are found in the high

1512

J. A. Deiber and R. A. Bortolozzi

Fig. 6. Comparison between numerical and experimental Nusselt numbers as function of Rayleigh number for the water—glass system, when d "22.3]10~3 m. Different lines show numerical results obtained p through different models, as explained in the text. Parameters are presented in Table 1.

Fig. 7. Nusselt number as function of Rayleigh number for the three cases studied in this work. Full lines are the model predictions. Parameters are the same as in Figs. 4—6.

Rayleigh number regime, our numerical results show that there is a clear branch from the curve of low and intermediate Rayleigh number regimes, in a consistent agreement with experimental data reported in the literature. Furthermore, at high Rayleigh numbers, solid and fluid temperature profiles cannot be considered equal (see Fig. 8) and the use of the two-field model is highly recommended. In fact, in Fig. 8 we find that for high permeability values (when d is large p and hence Da is high) there are zones of the annular saturated cavity where the solid temperature differs substantially from the fluid temperature (see also Amiri and Vafai, 1994). Finally, although the Boussinesq approximation cannot be considered for drastic thermal conditions in which the fluid density varies substantially (see, for example, Gartling and Hickox, 1985; Peirotti et al., 1987) we found that, in particular, this approximation is good for the experimental data used in our work. In

fact, by making / "1 in eqs (14)—(16) of our model, 0 we found that for the Rayleigh numbers involving the higher temperature differences (runs 3, 7 and 14) the Nusselt numbers calculated with and without the Boussinesq approximation do not differ in more than 0.4% with each other. This result is consistent with the fact that the experimental data at high values of Ra in the porous cavity of Fig. 1 are obtained with high values of Da, instead of using high values of *¹ (see also the Appendix).

CONCLUSIONS

From the comparison between experimental data of Nu vs Ra for three different values of Da and the h numerical predictions of the model proposed in this work, it is clearly concluded that for Ra'500 and relatively high permeability porous media, the following mechanisms of momentum and heat transfer shall

A two-field model at high Rayleigh numbers

1513

Fig. 8. Solid and fluid dimensionless temperatures at midheight of the porous cavity for different Rayleigh numbers and d "22.3]10~3 m: (a) near the hot wall, (b) near the cold wall. Other parameters are p presented in Table 1.

be considered apart from Darcian momentum exchange and convective—conductive heat transfer: (1) Viscous effects, through the Brinkman term. (2) Inertial effects, evaluated through the Forchheimer term. (3) Convection of momentum. (4) Fluid channeling due to the increase of porosity near the walls for relatively high particle diameter. (5) Anisotropic thermal dispersion in the fluid phase. (6) Energy exchange between phases (two temperatures fields) when the state of local thermal equilibrium between fluid and solid is not strictly valid. The inclusion of thermal dispersion in the energy balance of the fluid phase was found to be crucial in order to predict the correct slope of Nu as function of Ra for high particle diameters in a logarithmic plot. On the other hand, the consideration of the Boussinesq approximation and the neglect of terms involving porosity gradients and variations of thermophysical properties with temperature are excellent hypothesis for the porous cavity described in Fig. 1, although they are not used in the present study. Any model that fits experimental data by neglecting one or more of the mechanisms mentioned above, might be placing emphasis on those effects remaining in the theory and the resulting model is somewhat arbitrary for high Rayleigh and Darcy numbers. From previous works and the results presented here, there are enough evidence that for this regime the full model suggested in this work should be computed. Nevertheless, several theoretical aspects still need further elucidation, mainly concerning parameter values like a , a , u, l t A and B.

Acknowledgements The authors wish to thank the financial aid received from CONICET (Consejo Nacional de Investigaciones Cientı´ ficas y Te´cnicas de Argentina) and from Secretarı´ a de Ciencia y Te´cnica de la UNL (Universidad Nacional del Litoral, Argentina) — Programacio´n CAI#D : 94-0810-005033/06/ID. NOTATION

a a l a t a v A b B C 1 C 2 Co f d p D Da e a et a Fs g h

thermal hydrodynamic dispersivity tensor of components a ijkl longitudinal component of the dispersivity tensor transverse component of the dispersivity tensor specific surface of the porous medium, m~1 empirical constant used in eq. (4) Forchheimer inertial constant, m empirical constant used in eq. (4) viscosity constant, K~1 thermal conductivity constant, K~1 fluid heat capacity, J/kg K particle diameter, m porous cavity width, m Darcy number ("K /D2) = exchange of internal energy (a"s, f ), W m~3 exchange of total energy (a"s, f ), W m~3 Forchheimer number ("b/D) acceleration due to gravity, m s~2 heat-transfer coefficient at the cavity walls, W m~2K~1

1514

h (v ) sf f H k ko a ko m k d K l(n) ¸ m a n Nu i po f Pr f q a r r i r o R Ra Ra f Re p t ¹ a ¹ i T f v f V z Z

J. A. Deiber and R. A. Bortolozzi

solid—fluid heat-transfer coefficient, W m~2 K~1 dimensionless heat-transfer coefficient unit vector in z-direction thermal conductivity of pure species, W m~1 K~1 (a"s, f ) thermal conductivity of the saturated porous medium, W m~1 K~1 thermal dispersion tensor, W m~1 K~1 permeability of the porous medium, m2 Van Driest function height of the porous cavity, m exchange of momentum, Pa m~1 (a"s, f ) perpendicular distance from any wall, m overall Nusselt number defined in eqs (19) and (20) (i"h, c) fluid pressure, Pa fluid Prandtl number ("Coko /ko ) f f, r f, r mixture heat flux, W m~2(a"s, f ) radial coordinate, m inner radius of porous cavity, m outer radius of porous cavity, m ("r/D) Rayleigh number ("oo gbK D*¹/ f, r = ko ao ) f, r m, r fluid Rayleigh number ("Ra/Daj) particle Reynolds number ("eDv Doo d / f f, r p ko ) f, r time, s species temperature, K (a"s, f ) wall temperature, K (i"h, c) fluid stress tensor, Pa fluid velocity, m s~1 dimensionless superficial velocity axial coordinate, m ("z/D)

Greek letters ao thermal diffusivity of saturated porous mem dium, m2 s~1 b isobaric thermal expansion coefficient, K~1 c ("d /D) p d unity tensor e porosity g transformed axial coordinate h ["(¹ !¹ )/(¹ !¹ )](a"s, f ) a a r h c i ("r /r ) o i j ("ko /ko ) f, r m, r " ("1/b*¹) ko fluid viscosity, Pa s f l ("ko/ko ) s m, r m transformed radial coordinate oo fluid density, kg m~3 f p numerical average error defined by eq. (18) / functions used in eqs (14)—(16) ( j"0—4) j s numerical parameter defined by eq. (22) u Van Driest empirical constant Subscripts c cold wall f fluid phase

h m r s R

hot wall properties of composite (solid and fluid) reference temperature solid phase any property far from cavity walls

Superscripts o properties of pure species t total property REFERENCES

Amiri, A. and Vafai, K. (1994) Analysis of dispersion effects and non-thermal equilibrium, non-Darcian, variable porosity imcompressible flow through porous media. Int. J. Heat Mass ¹ransfer 6, 939—954. Benenati, R. F. and Brosilow, C. B. (1962) Void fraction distribution in beds of spheres. A.I.Ch.E. J. 8 (3) 359—361. Bird, R. B., Stewart, W. E. and Lightfoot, E. N. (1960) ¹ransport Phenomena. p. 197. Wiley, New York, U.S.A. Bowen, R. M., (1976) Theory of Mixtures. In Continuum Physics, ed. A. C. Eringen, Vol. III, pp. 1—127. Academic Press, New York, U.S.A. Chen, C.-K., Chen, C.-H., Minkowycz, W. J. and Gill, U. S. (1992) Non-Darcian effects on mixed convection about a vertical cylinder embedded in a saturated porous medium. Int. J. Heat Mass ¹ransfer 35(11), 3041—3046. Cheng, P. and Hsu, C. T. (1986) Applications of Van Driest’s mixing theory to transverse thermal dispersion in forced convective flow through a packed bed. Int. Comm. Heat Mass ¹ransfer 13, 613—625. Cheng, P. and Zhu, H. (1987) Effects of radial thermal dispersion on fully developed forced convection in cylindrical packed bed. Int. J. Heat Mass ¹ransfer 30(11), 2373—2383. Combarnous, M. A. and Bories, S. A. (1975) Hydrothermal convection in saturated porous media. Adv. Hydrosci. 1, 231—307. David, E., Lauriat, G. and Prasad, V. (1989) NonDarcy natural convection in packed-sphere beds between concentric vertical cylinders. A.I.Ch.E. Sym. Ser. 85(269) 90—95. Gartling, D. K. and Hickox, C. E. (1985) A numerical study of the applicability of the Boussinesq approximation for a fluid-saturated porous medium. Int. J. Numer. Methods Fluids 5, 995—1013. Georgiadis, J. G. and Catton, I. (1988) An effective equation governing convective transport in porous media. ASME J. Heat ¹ransfer 110, 635—641. Hong, J., Yamada, Y. and Tien, C. L. (1987) Effects of non-Darcian and nonuniform porosity on verticalplate natural convection in porous media. ASME J. Heat ¹ransfer 109, 356—362. Howle, L. E. and Georgiadis, J. G. (1994) Natural convection in porous media with anisotropic dispersive thermal conductivity. Int. J. Heat Mass ¹ransfer 37(7), 1081—1094. Ka´lnay de Rivas, E. (1972) On the use of nonuniform grids in finite difference equations. J. of Comput. Phys. 10, 202—210. Kladias, N. and Prasad, V. (1991) Experimental verification of Darcy—Brinkman—Forchheimer flow model for natural convection in porous media. J. ¹hermophys. 5, 560—576.

A two-field model at high Rayleigh numbers

Lauriat, G. and Prasad, V. (1987) Natural convection in a vertical porous cavity: a numerical study for Brinkman-extended Darcy formulation. ASME J. Heat ¹ransfer 109, 688—696. Lauriat, G. and Prasad, V. (1989) Non-Darcian effects on natural convection in a vertical porous enclosure. Int. J. Heat Mass ¹ransfer 32(11), 2135—2148. Marpu, D. R. (1995) Forchheimer and Brinkman extended Darcy flow model on natural convection in a vertical cylindrical porous annulus. Acta Mechanica 109, 41—48. Mercer, J. W., Faust, C. R., Miller, W. J. and Pearson, F. J. (1982) Review of simulation techniques for aquifer thermal energy storage. Adv. Hydrosci. 13, 1—129. Nield, D. A. and Bejan, A. (1991) Convection in Porous Media. Springer, New York, U.S.A. Peaceman, D. W. (1977) Fundamentals of Numerical Reservoir Simulation, pp. 37—41. Elsevier, Amsterdam, The Netherlands. Peirotti, M. B., Giavedoni, M. D. and Deiber, J. A. (1987) Natural convective heat transfer in a rectangular porous cavity with variable fluid properties—Validity of Boussinesq approximation. Int. J. Heat Mass ¹ransfer 30(12), 2571—2581. Poulikakos, D. and Bejan, A. (1985) The departure from Darcy flow in natural convection in a vertical porous layer. Phys. Fluids. 28(12), 3477—3484. Prasad, V. and Kulacki, F. A. (1984) Natural convection in a vertical porous annulus. Int. J. Heat Mass ¹ransfer 27(2), 207—219. Prasad, V. and Kulacki, F. A. (1985) Natural convection in porous media bounded by short concentric vertical cylinders. ASME J. Heat ¹ransfer 107, 147—154. Prasad, V., Kulacki, F. A. and Keyhani, M. (1985) Natural convection in porous media. J. Fluid Mech. 150, 89—119. Roache, P. J. (1972) Computational Fluid Dynamics, p. 73. Hermosa Publishers, Albuquerque, NM, U.S.A. Scheidegger, A. E. (1974) ¹he Physics of Flow through Porous Media, pp. 197—198. Univ. Toronto Press, Toronto.

1515

Truesdell, C. (1969) Rational ¹hermodynamics, pp. 81—109, McGraw Hill, New York, U.S.A. Vafai, K. (1984) Convective flow and heat transfer in variable-porosity media. J. Fluid Mech. 147, 233—259. Vafai, K. and Kim, S. J. (1995) On the limitations of the Brinkman—Forchheimer—extended Darcy equation. Int. J. Heat and Fluid Flow 16, 11—15. Vafai, K. and Tien, C. L. (1981) Boundary and inertia effects on flow and heat transfer in porous media. Int. J. Heat Mass ¹ransfer 24(1), 195—203. Wakao, N., Kaguei, S. and Funazkri, T. (1979) Effect of fluid dispersion coefficients on particle-to-fluid heat transfer coefficients in packed beds. Chem. Engng Sci. 34, 325—336.

APPENDIX

The purpose of this appendix is to report Table A1, where the values of several parameters required by the numerical algorithm for the 14 numerical runs are presented, consistently with the experimental data reported by Prasad et al. (1985) and David et al. (1989). These data are also useful for the Results and Discussion section. Additionally, one has to consider D"0.1239 m and r /r "5.338. In Table A1, numerical runs 1—3 correspond to o i d "3]10~3 m, runs 4—7 are for d "6]10~3 m; the other p p runs take d "22.3]10~3 m. p The following procedure must be carried out in order to obtain the experimental information from David et al.’s work (runs 1—11): (1) From Table 2 of the reference, the average temperature can be obtained for all the experimental data available by using Pr as entry value (the relation between Pr and the f f temperature is uniquely determined). (2) With the average temperature ¹ , all the thermophysir cal properties of the fluid (water) can be evaluated. Then, for a given Ra "Ra/Daj of Table 2 of the reference and with f the value of D and d (Table 1 of the reference reports d /D p p and d ) the value of *¹ is finally obtained. p (3) Additionally, the Darcy number is also given in Table 1 of the reference and j can be readily calculated for

Table A1. Dimensionless numbers and thermophysical properties for numerical runs Run no.

Ra (—)

Da (—)

j (—)

Pr f (—)

¹ r (K)

oo f, r (kg m~3)

b (K~1)

ko f, r (Pa s)

1 2 3 4 5 6 7 8 9 10 11 12 13 14

173 511 818 187 476 1489 3187 1806 5647 10 558 16 707 30 862 53 371 81 101

3.36]10~7 3.36]10~7 3.36]10~7 1.66]10~6 1.66]10~6 1.66]10~6 1.66]10~6 4.42]10~5 4.42]10~5 4.42]10~5 4.42]10~5 4.42]10~5 4.42]10~5 4.42]10~5

0.792 0.800 0.804 0.792 0.795 0.801 0.806 0.873 0.873 0.874 0.875 0.875 0.876 0.877

5.46 4.51 4.11 6.17 5.72 4.98 4.37 6.01 5.86 5.65 5.46 4.98 4.60 4.20

304 311 316 297 301 307 312 298 300 302 304 307 310 314

995 993 991 997 996 994 993 997 996 996 995 994 993 992

3.1]10~4 3.7]10~4 4.0]10~4 2.5]10~4 2.8]10~4 3.4]10~4 3.8]10~4 2.6]10~4 2.8]10~4 3.0]10~4 3.1]10~4 3.4]10~4 3.6]10~4 3.9]10~4

7.9]10~4 6.8]10~4 6.2]10~4 9.0]10~4 8.4]10~4 7.4]10~4 6.6]10~4 8.8]10~4 8.6]10~4 8.3]10~4 7.9]10~4 7.4]10~4 6.9]10~4 6.4]10~4

ko *¹ f, r (Wm~1K~1) (K) 0.619 0.630 0.635 0.609 0.615 0.624 0.631 0.610 0.613 0.617 0.619 0.624 0.628 0.633

13.3 28.6 38.9 4.2 8.5 20.0 33.9 1.3 3.5 5.9 8.8 14.3 21.9 28.4

1516

J. A. Deiber and R. A. Bortolozzi

each temperature. Therefore, with these calculations, the Rayleigh number is determined for each experimental point. It is important to mention here that, these values can be verified at the extremes of each experimental interval, because they are reported in this table. The values of parameters reported in Table A1 for runs 12—14, are not available in Table 2 of the reference,

however, a similar procedure to that described above can be followed with the Fig. 7 reported by Prasad et al. (1985). Also, it should be mentioned here that the 14 experimental data of Table A1 correspond to the same experimental program.