Numerical study on mixed convection in a horizontal flow past a square porous cavity using UNIFAES scheme

Numerical Study on Mixed Convection in a Horizontal Flow Past a Square Porous Cavity using UNIFAES Scheme

Jorge Llagostera

José R. Figueiredo

UNICAMP - Universidade Estadual de Campinas. Faculdade de Engenharia Mecânica

Departamento de Energia. C.P. 6122

13083-970 Campinas. SP. Brazil

[e-mail], [e-mail]

Mixed convection on the flow past a heated length and past a porous cavity located in a horizontal wall bounding a saturated porous medium is numerically simulated. The cavity is heated from below. The steady-state regime is studied for several intensities of the buoyancy effects due to temperature variations. The influences of Péclet and Rayleigh numbers on the flow pattern and the temperature distributions are examined. Local and global Nusselt numbers are reported for the heated surface. The convective-diffusive fluxes at the volume boundaries are represented using the UNIFAES, Unified Finite Approach Exponential-type Scheme, with the Power-Law approximation to reduce the computing time. The conditions established by Rivas for the quadratic order of accuracy of the central differencing to be maintained in irregular grids are shown to be extensible to other quadratic schemes, including UNIFAES, so that accuracy estimates could be obtained.

Keywords: UNIFAES scheme, porous medium, mixed convection, open cavity, irregular grids

Introduction

Mixed convection in saturated porous media is important in fields such as geothermal resources, high performance thermal insulation systems, contaminants dispersion, packed beds, and the exploration of petroleum and gas fields. Horizontal layers heated from below and rectangular cavities with vertical temperature gradients have been studied by several authors, since the classical studies on flow stability on this configuration due to Horton and Rogers (1945), and Lapwood (1948), as can be seen in the reviews by Cheng (1978), Nield and Bejan (1992), and Bejan (1995).

The combined natural and forced convection in porous media was studied experimentally by Combarnous and Bia (1971), focusing the influence of an imposed pressure gradient. Similarity solutions for thermal boundary layers in porous media involving mixed convection near heated surfaces were obtained by Cheng (1977). Prasad et al. (1988) studied numerically the mixed convection in horizontal porous layers locally heated from below. Lai and Kulacki (1991a) obtained experimental and numerical results for a similar configuration with uniform heat flux on the heated section. In other paper Lai and Kulacki (1991b) report numerical results for oscillatory mixed convection in horizontal porous layers locally heated from below . Brouwers (1994) studied the heat transfer process between a saturated porous medium and a permeable wall with fluid injection or withdrawal. Manole et al. (1994) considered theoretically the convection induced by thermal and solutal gradients in a horizontal layer. Vynnycky and Pop (1997) investigated analytically and numerically the forced and free convective flow past a heated, or cooled, body in a semi-infinite porous medium. Dufour and Neel (1998) obtained numerically the two-dimensional convective flow patterns in a long horizontal porous layer heated from below, where a cross flow is imposed.

This paper applies a new discretization scheme for convective and diffusive transport equations, UNIFAES, to the analysis of combined natural and forced convection in a two dimensional, saturated porous layer, with a cavity on its bottom wall. The influences of Péclet and Rayleigh numbers on the Nusselt number and the flow pattern are investigated.

Figure 1 presents a schematic view of the porous configuration studied. The fluid temperature upstream of the cavity is T₀ and the bottom of the cavity is isothermally heated (T₁). The top region of the layer is isothermally cooled (T₀) and the other walls are adiabatic. In order to simulate infinite flow domains the total height of the calculation domain (H) as well as the upstream (L_u) and downstream (L_d) regions sizes are long compared to the heated length L. The cavity depth (C/L) values considered were 0 and 1.

The smooth behavior of the solution functions far from the cavity allows the use of coarse grids. Geometrically varying cell sizes were employed, increasing with the distance from the cavity or the heated surface. The consequences of the variable size of the cells on the accuracy are evaluated on the basis of the theoretical analysis presented by Rivas (1972), which established certain conditions for the order of spatial convergence of the scheme to be maintained.

Nomenclature

Latin Symbols

B = layer thickness, Fig. 1

C = cavity depth, Fig. 1

g = gravitational acceleration

i, j = node indexes

K = source term of the interpolating curve's generating equation

H = calculation domain height, H = B + C

L = cavity length, Fig. 1.

M = real parameter in exact solution of Burgers' equation.

n = total number of intervals in discretized Burgers' equation

p = cell Péclet number, Eqs. (6) and (7)

P = global modified Péclet number, P = U L / a

R = geometric ratio, Eq. (11)

Ra = modified Rayleigh number:

r = geometric ratio between successive intervals in irregular grids

T = temperature

T_o, T₁ = cavity boundary temperatures, Fig. 1

u,v = velocity components, normalized by U

U = characteristic velocity

x,y = spatial coordinates, normalized by L

Greek Symbols

a = porous medium effective thermal diffusivity

b = isobaric thermal expansion coefficient of the fluid

c = function of the cell Péclet number, Eq. (10)

dx, dy = cell dimensions, Fig. 2

Dx, Dy = node distances, Fig. 2

f = non-dimensional temperature, f = (T - T_o)/(T₁- T_o)

G = thermal conductivity

K = permeability of the saturated porous medium

m = fluid dynamic viscosity

p = function of the cell Péclet number, Eq. (8)

P = influence coefficient in generalized Allen and Southwell scheme, Eqs. (13) and (14)

r = fluid density

r_o = fluid density at T_o

s = porous medium thermal capacity ratio

t = time normalized by sL /U

y = stream function, normalized by UL

Y = extra convective-diffusive discretization term, Eq. (9)
w = parameter in geometrically varying irregular grid,
w = n 1n(r).
x = regular, transformed space, x ( x _i ) = i/n

Subscripts

e, w = cell faces, Fig. 2

n, s = cell faces, Fig. 2

E, W = neighbor nodes, Fig.2

N, S = neighbor nodes, Fig. 2

i = node index

P = central node, Fig. 2

u, d = upstream and downstream, Fig. 1

x = x-direction

Superscripts

+ = positive side of central node

- = negative side of central node

Formulation

Impermeable walls bound the porous layer. The flow is incompressible, steady state and two-dimensional. Constant fluid properties are assumed, and the porous matrix is assumed as rigid, homogeneous and isotropic. Boussinesq hypothesis about the buoyancy effects and the applicability of Darcy law for the relation between velocity and pressure gradient are assumed. The fluid and the porous matrix are in local thermal equilibrium.

Adopting the stream function formulation expressed by Eq. (1), the mass conservation and the momentum equations are represented by Eq. (2).

Equation (3) corresponds to the energy conservation, where a transient term is put for the numerical approach to the steady solution:

The differential equations were approximated numerically using a computer code developed to simulate problems dealing with porous media (Llagostera, 1990 and 1994). The grid was defined by dividing the domain in rectangles representing the finite volumes, with the nodal points located at their centers. The equations were discretized by integration over finite volumes. The Poisson type Eq. (2) was discretized by finite-volume central differencing, with linear interpolation. For the energy transport, Eq. (3), the convective and diffusive fluxes on the volume boundaries were represented using the UNIFAES, Unified Finite Approach Exponential-type Scheme, summarized in the next section, which has shown excellent accuracy and stability records even at very high Péclet numbers (Figueiredo, 1997; Figueiredo and Llagostera, 1999).

Discretization Scheme

UNIFAES belongs to the family of conservative-form exponential-type schemes, whose interpolating curves in the x-direction, for instance, are given by the exact solution of the ordinary equation:

This interpolating curve's generating equation approximates the partial differential equation (3) assuming locally constant values for the velocity 'u' as well as for the transient and cross-derivative terms, globally represented by the source term K_x.

The exact solution to Eq. (4) has two integration constants determined by fitting to the nodes P and E for the 'e' cell boundary, etc., allowing the combined convective-diffusive fluxes in each cell to be found. Using the continuity equation, the difference equation of transport becomes:

where

The terms K_e, etc., are obtained by means of the Allen and Southwell (1955) scheme generalized to irregular grids, which is also based on the interpolating curve (4), yielding:

where

In the case Dx ⁺ =Dx ^- =Dx this expression reduces to the Allen and Southwell's scheme:

The source term K_e for the node (i, j) is found by linear interpolation from the generalized Allen's estimates of K_x at nodes (i, j) and (i+1, j), and so on for the other cell boundary source terms. At cell boundaries neighbor to the domain frontiers the source term is linearly extrapolated from the internal nodes.

The computational cost of the exponential function in Eqs. (8), (10), (13) and (14) is minimized using Patankar's (1980) Power-Law approximation to function p (Figueiredo and Llagostera, 1999).

Grid Irregularity

Rivas (1972) demonstrated by Taylor series analysis that the quadratic spatial convergence of the central differencing is maintained for irregular grids provided there is a function x(x) mapping the regularly spaced transformed domain x into the irregularly spaced physical domain x, with bounded derivatives different from zero, with the possible exception of the origin.

A simple one-dimensional test case will show that Rivas' conclusion applies not only to the central differencing schemes (the conservative form finite-volume central differencing with linear interpolation, and the non-conservative form finite-difference central differencing with parabolic interpolation), but also to the exponential type schemes, including the above generalization of the Allen and Southwell non-conservative form exponential scheme, the conservative form simple Exponential Scheme, that coincides with LOADS in this problem, and UNIFAES.

The physical domain 0 £ x £ 1 is discretized according to a geometric series with ratio 'r'. The variable grid index 'i' is defined from 0 to 'n', and the node coordinates are:

The equivalence between the two forms above is established by setting w =nln(r) and x =i/n. The second form presents a transformation rule that suits Rivas' conditions if w is kept constant with refinement, i.e., if the refinement path is such that the geometric ratio 'r' of the variable cell size decreases with the number of nodes 'n' in such a way that nln(r) is constant.

The investigation is performed with the steady-state non-linear Burgers' equation, which can be written in both conservative and non-conservative forms:

with boundary conditions u(0)=1 and u(1)=0. Its exact solution is:

where M=1.0 (to precision 10^-10) for P=50.0

Figure 3 shows results for the case P =50, n =10 and r=0.8. All the schemes required a simple iterative procedure to cope with the non-linearity. Most numerical results are visually coincident with the exact solution, except the non conservative central differencing that shows some wiggles.

Figures 4 and 5 provide the percentage root mean square error against the arithmetic mean grid spacing, Dx_m = 1/n, employing distinct refinement strategies. In Fig. 4 the parameter w is constant with refinement, w = n ln(r) = 10 ln(0.8) = - 2.2314355, so the geometric ratio 'r' decreases with the number of nodes 'n'; the quadratic asymptotic behavior of all curves indicates that Rivas's conclusion applies to all schemes used. To emphasize the importance of Riva's requirement, in Fig. 5 the geometric ratio is kept constant with refinement (r = 0.8) instead of the parameter w , so the space transformation is increasingly abrupt with refinement; all convergence tendencies are ultimately sub-linear.

Results for Mixed Convection

The system of discretized equations was solved using the line-to-line process, with alternating direction. For each iteration the domain was swept 16 times for stream function calculation and 16 times for temperature calculation, to accelerate the convergence process. The initial guesses for stream function and temperature values were null distributions. The stopping criteria were that the relative variation of both the stream function and the temperature between successive iterations should be below 10^-5, sufficient to obtain stabilized results.

The calculation domain, according to Fig. 1, has a non-dimensional length of 16.554, both the upstream length (L_u/L) and the downstream length (L_d/L) being 7.777. The thickness of the porous layer (B/L) is 11.466 for the cases where the cavity depth (C/L) is 1. For the configuration with C/L=0 the thickness B/L was equal to 24.743, to allow for more significant temperature gradients caused by the direct exposition of the heated length to the forced flow inside the horizontal layer.

For the configuration with C/L=1, the region of the square cavity is divided in 50 equal intervals, both in horizontal and vertical directions. The same interval is maintained, in the horizontal direction, in two regions of width 0.2 L, upstream and downstream of the cavity. From there, in both upstream and downstream directions, irregular grids with geometric ratio 1.04 are employed. Therefore, the horizontal grid is composed by five regions, with 70, 10, 50, 10 and 70 intervals. In the vertical direction the grid is uniform inside the cavity (50 intervals); above the cavity an irregular grid with ratio 1.04 is employed, corresponding to 80 intervals. The total number of finite control volumes in the domain for C/L=1 is then equal to 19,300. For the configuration with C/L=0, the horizontal grid is the same which is used for the previous case and the vertical grid uses 100 intervals with ratio 1.04, the domain being composed by a total of 21,000 finite control volumes.

The global Nusselt values at the heated surface, calculated according to Eq. (18), are presented in Table 1, for various Péclet and Rayleigh numbers, for the geometric configurations with C/L=0 and C/L=1. Some results obtained by Prasad et al. (1988) for a configuration with B/L=1 and C/L=0 are presented in Table 1 for comparison. The left corner of the heated length corresponds to x=0, and the vertical position of its surface is given by y=0.

Figure 6 presents the variation of global Nusselt number for the heated surface as a function of Rayleigh number for Péclet values 1, 10 and 100, with C/L=0 and C/L=1. The curves show that at low Rayleigh numbers the influence of forced convection is dominant, and for high Rayleigh values the influence of buoyancy is significant, mainly for P=1. For P=100, in the configuration with C/L=1, a decrease in Nusselt number is observed between Rayleigh numbers 200 and 300, caused by the formation of a circulation cell inside the square cavity.

For low Peclet (P=1) the flow pattern is more sensible to the buoyancy effects. For the configuration with C/L=1, the values of the stream function indicate the formation of a circulating cell for Rayleigh numbers equal or above 50. For Ra=20 or less the maximum stream function value is Y_max =11.466. For Rayleigh number equal to 50, 60, 70 and 80, the results for Y_max are 13.505, 15.118, 16.549 and 17.764, respectively, indicating the occurrence of circulation. All the values obtained for the minimum stream function are null. For the configuration with C/L=0 no circulation was detected for the studied range. The streamlines and constant temperature lines near the heated surface and inside the open square cavity for P=1 and Ra=50 are presented in Fig. 7, using Dy =1 and Df =0.1. The influence of Rayleigh number on the flow pattern and temperature gradients inside the cavity can be observed.

For P=10 the Nusselt number shows an important increase for low Rayleigh numbers resulting from the strong influence of the pressure gradients imposed on the horizontal porous layer, overcoming the influence of buoyancy effects. For P=10 and higher values of Rayleigh number no circulation is observed in the range studied. The streamlines and isothermal lines near the heated wall and inside the open square cavity for P=10, with Ra=50 and Ra=150, are shown in Fig. 7, adopting Dy =0.1 and Df =0.1. The influence of Rayleigh number on the flow pattern and temperature gradients inside the square cavity for P=10 is clearly visualized.

For Péclet equal to 100 no circulation is observed in the horizontal layer, and a more intense forced flow pattern is established until Ra=300, when a small circulation cell is observed inside the cavity, indicated by the value of Y_min equal to -0.0653. For all the other cases studied Y_min is null. The streamlines and isothermal lines near the heated surface and inside the open cavity for P=100, with Ra=50 and Ra=300, are shown also in Fig. 7, using Dy =0.05 and Df =0.1. The influence of Rayleigh number on the flow pattern and temperature gradients are shown, and for C/L=1, with Ra=300 the circulation cell formed near the cavity bottom is also observed.

The maps for P=1 show a more intense influence of buoyancy effects on fluid flow. The maps for P=100 reveal the strong domination of the pressure gradient in the horizontal direction, presenting a pattern similar to the observed in typical cases of forced convection. The isothermal lines near the heated surface are grouped more tightly for the case with P=100, indicating a more intense heat transfer.

Local Nusselt distributions on the heated surface for P=100 are presented in Fig. 8 for both geometric configurations. For C/L=0 the values of local Nusselt are very similar for the Rayleigh numbers considered, because the heated surface is directly exposed to the forced flow and the influence of Rayleigh number is less important. For C/L=1 the distributions for Ra=0 and Ra=50 are similar, and show a forced convection pattern. For Ra=300 one notices an increase of the local Nusselt number in the right region of the heated surface indicating a clear change in the flow pattern inside the cavity, caused by the formation of a circulation cell (clockwise).

Discussion on Accuracy

Two main features are involved in the accuracy considerations for this problem: the grid refinement and the finite size of the modeled infinite domain.

The case defined by Ra=200 and P=100 is employed for grid refinement testing. The cavity length previously divided in 50 intervals, was divided in 60 intervals, and the others parts of the domain were divided proportionally. According to the procedure described in section 4, the geometric ratio must decrease from 1.04 to 1.033224 so that quadratic decrease of the error with refinement can be expected. Results for the global Nusselt number are given in Table 3. The last line shows the value obtained with Richardson's extrapolation assuming quadratic behavior. The extrapolated value indicates the high accuracy of UNIFAES. The global Nusselt number obtained with the grid employed differs from the extrapolated value by about 0.2%.

The size of the domain was tested by parts. First, the influence of the layer thickness on the heat transfer results from the heated surface was studied. For the cases P=1, Ra=50 and P=10, Ra=100, eight different values of H/L were used, corresponding to 30, 40, 50, 60, 70, 80, 90 and 100 intervals in the layer. The results obtained for the global Nusselt number on the heated (y=0) and on the cooled surfaces (y=H/L) are reported in Table 2, and illustrated in Fig. 9.

For P=1 and Ra=50, where the buoyancy forces are important, H/L=12.466 is sufficient to make the influence of the top wall on the heat transfer from the heated surface negligible. Figure 10 presents the distribution of local Nusselt values on the heated surface (y=0) where the decrease of H/L influence on the results is evident. The distribution of local Nusselt values on the cooled surface (y=H/L) for the same cases is shown in Figure 11. For the highest values of H/L the heat transfer through the cooled surface becomes negligible, because nearly all the energy received from the heated wall is carried away by forced convection.

With P=10 and Ra=100, where the forced convection influence is stronger, H/L=5.950 is enough to make the influence of the upper wall negligible, as can be seen in Figure 9.

In both cases the global Nusselt number in the heated surface varied less than 0,03% when the total height varied from the standard value 12.4659 to the greatest value 26.7426.

Finally the influences of the inlet and outlet lengths are considered through the case P=10, Ra=100. When those lengths are given the values 7.777, 17.422 and 25.943, the Nusselt number assumes the values 3.4334, 3.4333 and 3.4333 respectively, with a variation much below 0.01%.

Concluding Remarks

The approximated form of UNIFAES was successfully applied to the problem of mixed convection on a cavity in the bottom wall of a porous layer. A comprehensive set of numerical simulations was presented, covering several cases of mixed convection, represented by different Péclet and Rayleigh numbers.

Two main flow conditions have been found: the pressure gradient dominated flows, at low Rayleigh numbers, and the flows influenced by buoyancy effects at high Rayleigh numbers, particularly for low Péclet flows. When the heated length is exposed directly to the forced flow the heat transfer is more intense and the influence of the Rayleigh number is less significant.

Inspired by Rivas's analysis, it has been shown that varying irregular meshes may be employed without loss of the quadratic spatial convergence if the geometric ratio 'r' and the number of nodes 'n' are such that 'n ln(r)' is constant. This allowed the accuracy of the present computations to be investigated by grid refinement with Richardson's extrapolation for a case having intense buoyancy and pressure gradients, indicating that the global Nusselt number was measured with error about 0.2% with respect to spatial refinement. The domain size was shown to introduce even smaller errors.

Acknowledgements

The authors are grateful to the "Centro Nacional de Processamento de Alto Desempenho em São Paulo (CENAPAD-SP)" for the use of its computing facilities.

Manuscript received: May 1999. Technical Editor: Angela Ourívio Nieckele.

• Allen, D.N. de G., and Southwell, R.V., 1955, "Relaxation Methods Applied to Determine the Motion, in Two Dimensions, of a Viscous Flow Past a Fixed Cylinder", Quart. J. Mech. and Applied Math., Vol.8, pp. 129-145.
• Bejan, A., 1995, "Convection Heat Transfer", 2nd ed., John Wiley, New York
• Brouwers, H.J.H., 1994, "Heat-Transfer Between a Fluid Saturated Porous-Medium and a Permeable Wall with Fluid Injection or Withdrawal", Int. J. Heat Mass Transfer, Vol. 37, pp. 989-996.
• Combarnous, M.A. and Bia, P., 1971, "Combined Free and Forced Convection in Porous Media", Soc. Petrol. Engng. Journal, Vol.11, pp.399-405.
• Cheng, P., 1977, "Similarity Solutions for Mixed Convection from Horizontal Impermeable Surfaces in Saturated Porous Media", Int. J. Heat Mass Transfer, Vol.20, pp.893-898.
• Cheng, P., 1978, "Heat Transfer in Geothermal Systems", Adv. Heat Transfer, Vol.14, pp.1-105.
• Dufour, F. and Neel, M.C., 1998, "Numerical Study of Instability in a Horizontal Porous Channel with Bottom Heating and Forced Horizontal Flow", Physics of Fluids, Vol.10, pp.2198-2207.
• Figueiredo, J.R., 1997, "A Unified Finite-Volume Finite-Differencing Exponential-Type Scheme for Convective-Diffusive Fluid Transport Equations", J. Brazilian Society Mech. Sciences, Vol.19, No.3, pp.371-391.
• Figueiredo, J.R., and Llagostera, J., 1999, "Comparative Study of the Unified Finite Approach Exponential-Type Scheme (UNIFAES) and its Application to Natural Convection in Porous Cavity", Numerical Heat Transfer, Part B., Vol.35, pp. 347-367.
• Horton, C.W. and Rogers R.T., 1945, "Convection Currents in a Porous Medium", J. Appl. Phys., Vol.16, pp.367-370.
• Lai, F.C. and Kulacki, F.A., 1991a, "Experimental Study of Free and Mixed Convection in Horizontal Porous Layers Locally Heated from Below", Int. J. Heat Mass Transfer, Vol.34, pp.525-541.
• Lai, F.C. and Kulacki, F.A., 1991b, "Oscillatory Mixed Convection in Horizontal Porous Layers Locally Heated from Below", Int. J. Heat Mass Transfer, Vol.34, pp.887-890.
• Lapwood, E.R., 1948, "Convection of a Fluid in a Porous Medium", Proceedings of Cambridge Philosophical Society, Vol.44, pp.508-521.
• Llagostera, J., 1990, "Mixed Convection Past a Porous Cavity", Doctoral Thesis, Universidade Estadual de Campinas (UNICAMP), Campinas, SP, Brazil. (In Portuguese)
• Llagostera, J., 1994, "Natural Convection in an Irregular Porous Cavity: Effects of Thermal Boundary Conditions", in Proceedings of the 10th International Heat Transfer Conference, Brighton-UK, edited by Hewitt, G.F., Taylor & Francis, Vol.7, pp. 103-108.
• Manole, D.M., Lage, J.L., and Nield, D.A., 1994, "Convection Induced by Inclined Thermal and Solutal Gradients, with Horizontal Mass-Flow, in a Shallow Horizontal Layer of a Porous Medium", Int. J. Heat Mass Transfer, Vol.37, pp.2047-2057.
• Nield, D.A., and Bejan, A., 1992, "Convection in Porous Media", Springer-Verlag, Berlin.
• Patankar, S.V., 1980, "Numerical Heat Transfer and Fluid Flow", McGraw-Hill, New York.
• Prasad, V., Lai, F.C., and Kulacki, F.A., 1988, "Mixed Convection in Horizontal Porous Layers Heated from Below", J. Heat Transfer, Vol.110, pp. 395-402.
• Rivas, E.K. de, 1972, "On the Use of Non-Uniform Grids in Finite-Difference Equations", J. Computational Physics, Vol.10, pp.202-210.
• Vynnycky, M. and Pop, I., 1997, "Mixed Convection Due to a Finite Horizontal Flat Plate Embedded in a porous medium", J. Fluid Mechanics, Vol.351, pp.359-378.

Publication Dates

• Publication in this collection
  14 Sept 2001
• Date of issue
  2000

History

• Received
  May 1999