Natural convection heat transfer in shallow horizontal rectangular enclosures uniformly heated from the side and filled with non-Newtonian power law fluids

Energy Conversion and Management 47 (2006) 2535–2551

www.elsevier.com/locate/enconman

Natural convection heat transfer in shallow horizontalrectangular enclosures uniformly heated from the side and filled

with non-Newtonian power law fluids

M. Lamsaadi a, M. Naımi a,*, M. Hasnaoui b

a Universite Cadi Ayyad, Faculte des Sciences et Techniques de Beni-Mellal, Departement de Physique, UFR de Chimie Appliquee et Sciences

de l’Environnement, Equipe de Modelisation des Ecoulements et des Transferts (EMET), BP 523 Beni-Mellal, Moroccob Universite Cadi Ayyad, Faculte des Sciences Semlalia, Departement de Physique, UFR de Thermique et Mecanique des Fluides,

Laboratoire de Mecanique des Fluides et d’Energetique (LMFE), BP 2390 Marrakech, Morocco

Received 10 May 2005; accepted 30 October 2005Available online 13 December 2005

Abstract

A combined analytical and numerical study is conducted for two dimensional, steady state, buoyancy driven flows ofnon-Newtonian power law fluids confined in a shallow rectangular cavity submitted to uniform fluxes of heat along bothits short vertical sides, while its long horizontal walls are considered adiabatic. The effect of the non-Newtonian behavioron the fluid flow and heat transfer characteristics is examined. An approximate theoretical solution is developed on thebasis of the parallel flow assumption and validated numerically by solving the full governing equations.� 2005 Elsevier Ltd. All rights reserved.

Keywords: Heat transfer; Natural convection; Non-Newtonian fluids; Rectangular enclosures

1. Introduction

Natural convection in a fluid filled enclosure is of great importance because this geometry is present in along list of engineering and geophysical systems. Because of the non-linear nature of the governing equations,analytical solutions are possible in some specific situations related to the geometry and also the nature of theboundary conditions imposed for a given system. According to Cormack et al. [1], analytical calculations arepossible in situations where the length of the cavity is much larger than its height. Thus, using matched asymp-totic expansions, in the case of a shallow cavity with differentially heated ends, these authors showed that theflow is parallel to the long sides of the cavity in the core region and presents a non-parallel structure near itsends. An analytical solution was found possible in the core region, and the limits of the parallel flow structurevalidity were discussed. However, for this case, the authors reported that the determination of the constant

0196-8904/$ - see front matter � 2005 Elsevier Ltd. All rights reserved.

doi:10.1016/j.enconman.2005.10.028

* Corresponding author. Tel.: +212 2348 5112/5122/5182; fax: +212 2348 5201.E-mail addresses: [email protected], [email protected] (M. Naımi).

mailto:[email protected]

mailto:[email protected]

Nomenclature

A aspect ratio of rectangular enclosure, Eq. (10)C dimensionless temperature gradient in horizontal direction (x-direction)g acceleration due to gravity (m/s2)H 0 height of enclosure (m)k consistency index for power law fluid at reference temperature (Pa sn)L 0 length of rectangular enclosure (m)n flow behavior index for power law fluid at reference temperatureNu local Nusselt number, Eqs. (11), (12) and (29)Nu mean Nusselt number, Eqs. (13) and (29)Pr generalized Prandtl number, Eq. (10)q 0 constant density of heat flux (W/m2)Ra generalized Rayleigh number, Eq. (10)T dimensionless temperature, ½¼ ðT 0 � T 0cÞ=DT ��T 0c reference temperature at geometric centre of enclosure (K)DT* characteristic temperature (=q 0H 0/k) (K)(u,v) dimensionless axial and transverse velocities [=(u 0,v 0)/(a/H 0)](x,y) dimensionless axial and transverse co-ordinates [=(x 0,y 0)/H 0]

Greek symbols

a thermal diffusivity of fluid at reference temperature (m2/s)b thermal expansion coefficient of fluid at reference temperature (1/K)k thermal conductivity of fluid at reference temperature (W/m C)l dynamic viscosity for Newtonian fluid at reference temperature (Pa s)la dimensionless apparent viscosity of fluid, Eq. (6)q density of fluid at reference temperature (kg/m3)X dimensionless vorticity, ½¼ X0=ða=H 02Þ�w dimensionless stream function, (=w 0/a)

Superscript0 dimensional variables

Subscriptsc value relative to centre of enclosure (x,y) = (A/2,1/2)max maximum value

2536 M. Lamsaadi et al. / Energy Conversion and Management 47 (2006) 2535–2551

horizontal temperature gradient in the parallel flow region is not obvious and requires a combination of ana-lytical and numerical methods, which leads to onerous computation time. On the other hand, when uniformheat fluxes are imposed instead on the vertical short sides of the enclosure, the calculation of the constant axialtemperature gradient in the parallel flow region becomes easy by performing an energy balance at any verticalsection of the enclosure (see for instance Bejan [2]). Later, the parallel flow approximation was used by severalauthors for both porous [3–5] and clear fluid [6,7] media submitted to horizontal uniform heat flux. The pres-ent study is motivated by the lack of analytical works related to non-Newtonian fluids confined in shallowcavities. Compared to the Newtonian case, studies concerning non-Newtonian flows in clear fluid mediaare relatively less numerous despite their importance and presence in many industrial applications, such aspaper making, oil drilling, slurry transporting, food processing, polymer engineering and many others. Someof these applications are discussed in Jaluria [8]. To the best of the authors� knowledge, the numerical studyconducted by Ozoe and Churchill [9] on natural convection in the case of confined non-Newtonian fluidsheated from below with a constant temperature, counts among the first studies that considered the rheological

M. Lamsaadi et al. / Energy Conversion and Management 47 (2006) 2535–2551 2537

behavior of fluids. It was mainly a question of determining the threshold for the onset of Rayleigh–Benardconvection in power law fluids. The critical Rayleigh number was found to increase with the flow behaviorindex, but the values obtained overestimated those obtained experimentally and theoretically in an earlierstudy by Tien et al. [10] relative to thermal instability in a horizontal layer of a non-Newtonian power lawfluid heated from below with a constant temperature. In the case of a rectangular cavity differentially heatedand filled with non-Newtonian fluids, a numerical modeling of natural convection was conducted by Turki[11]. According to the numerical results reported by this author, the non-Newtonian behavior affects consid-erably the flow structure and heat transfer. On the other hand, when the fluid viscosity is important, the Archi-medes forces are generally insufficient to induce a motion with a significant effect of convection. Hence, itbecomes necessary to promote the fluid circulation by putting either one or two moving horizontal boundariesof the cavity (the model of the dragged cavity), which leads to important modifications in the flow structureand heat transfer.

In the present study, the attention is focused on the natural convection heat transfer problem within a twodimensional horizontal rectangular enclosure filled with a non-Newtonian fluid. The cavity is heated andcooled with uniform heat fluxes from its short vertical sides, while its long horizontal boundaries are insulated.The power law model, suggested originally by Ostwald–de Waele, is used to characterize the non-Newtonianfluid behavior. A numerical solution of the full governing equations is obtained for wide ranges of the con-trolling parameters. The results, in terms of velocity and temperature profiles and heat transfer rates, are pre-sented for different values of Rayleigh number and power law index. Useful correlations of heat transferresults are also obtained in wide ranges of n and Ra.

2. Mathematical formulation

The studied configuration, sketched in Fig. 1, is a rectangular enclosure of height H 0 and length L 0. Thelong horizontal walls are thermally insulated, while the vertical short sides are submitted to constant heat flux,q 0. The non-Newtonian fluids considered here are those whose rheological behavior can be described by thepower law model (also known as the Ostwald–de Waele, or two parameters model), which leads to a relation-ship between the shear stress and shear rate. In terms of the laminar apparent viscosity, this relation can bewritten, in the present case, as follows:

l0a ¼ k 2ou0

ox0

� �2

þ ov0

oy 0

� �2 !

þ ou0

oy 0þ ov0

ox0

� �2 !n�1

2

ð1Þ

where n is a dimensionless constant called the power law index, and k is an empirical coefficient known as theconsistency factor, which is an indicator of the degree of the fluid�s viscosity. Note that for n = 1, the powerlaw model reduces to Newton�s law by setting k = l. Thus, the deviation of n from unity characterizes the de-gree of non-Newtonian behavior of the fluid. Specifically, when n is in the range 0 < n < 1, the fluid is said tobe pseudo-plastic (or shear thinning), and the viscosity is found to decrease by increasing the rate of strain. Onthe other hand, when n > 1, the fluid is said to be dilatant (or shear thickening), and the viscosity increases byincreasing the rate of strain. Dilatant fluids are generally much less common than pseudo-plastic ones. Asexamples of pseudo-plastic and dilatant fluids, we mention the 4% paper pulp in water (k = 20 Pa sn,

Fig. 1. Sketch of the enclosure and co-ordinates system.


n = 0.575) and 41% of corn starch suspensions in an aqueous sucrose solution of 50% concentration by weight(k = 0.639 Pa sn, n = 1.397), respectively. Though the Ostwald–de Waele model does not converge to aNewtonian behavior in the limit of zero and maximum rates of shear, it presents, however, the advantageof being simple and mathematically tractable. In addition, the rheological behavior of many substances canbe adequately represented by this model for a relatively large range of shear rates (or shear stresses), whichmakes it useful, at least for engineering purposes, and justifies its use in most theoretical investigations of fluidshaving pseudo-plastic or dilatant behaviors. The main assumptions made here are those commonly used, i.e.,

• The fluid velocities are small enough to consider the flow as laminar. In most buoyancy driven motions, thefluid circulation is slow due to moderate temperature gradients [12].

• The fluid is incompressible. For pressures close to atmospheric, liquids can be considered as incompressiblefluids with a good approximation.

• The viscous dissipation is negligible. According to Turki [11], this approximation is valid for polymer solu-tions, weakly or moderately concentrated, such as aqueous solutions of carboxymethylcellulose (CMC).For polymer solutions highly concentrated or polymer melts, the viscous frictions are so significant thatit becomes not obvious to generate the fluid motion simply by buoyancy effects.

• The physical properties are considered independent of temperature except for the density in the buoyancyterm, which obeys the Boussinesq approximation, whose systematic derivation from the equations of com-pressible fluid dynamics is quite tedious and requires careful ordering of several limiting processes. Moredetails are given in the paper of Gray and Giorgini [13].

• The third dimension of the cavity is large enough so that the problem is two dimensional. This assumptionis generally relatively well satisfied and provides insight into the more complicated three dimensional flows[12].

Then, the dimensionless governing equations, written in terms of vorticity, X, temperature, T, and streamfunction, w, are as follows:

oXotþ oðuXÞ

oxþ oðvXÞ

oy¼ Pr la

o2Xox2þ o2X

oy2

� �þ 2

ola

oxoXoxþ ola

oyoXoy

� �� þ SX ð2Þ

oTotþ oðuT Þ

oxþ oðvT Þ

oy¼ o

2Tox2þ o

2Toy2

ð3Þ

and

o2wox2þ o2w

oy2¼ �X ð4Þ

where

u ¼ owoy

; v ¼ � owox

ð5Þ

la ¼ 2ouox

� �2

þ ovoy

� �2" #

þ ouoyþ ov

ox

� �2" #n�1

2

ð6Þ

and

SX ¼ Pro2la

ox2� o2la

oy2

� �ouoyþ ov

ox

� �� 2

o2la

oxoyouox� ov

oy

� �� þ PrRa

oTox

ð7Þ

This formulation presents the advantage of reducing the number of equations by eliminating the pressure asa variable, though it can be calculated after the solution is obtained if need be. Such an approach is generallyadvantageous, as compared to the primitive variable approach, for two dimensional and axisymmetric flows.The latter approach is more appropriate for three dimensional circumstances.


The dimensionless variables are obtained by using the characteristic scales H 0, H 02=a, a/H 0, a=H 02

, q 0H 0/kand a corresponding to length, time, velocity, vorticity, characteristic temperature and stream function,respectively.

For all the above equations, derivations from basic principles are sketched in the Refs. [9,11,14].To complete the problem formulation, the following non-dimensional boundary conditions are used:

TableConve

Grids

n

0.6

0.8

1.0

1.2

1.4

u ¼ v ¼ w ¼ oT =oxþ 1 ¼ 0 for x ¼ 0 and A ð8Þu ¼ v ¼ w ¼ oT =oy ¼ 0 for y ¼ 0 and 1 ð9Þ

In addition to the power law index, n, the present problem is governed by three other dimensionless param-eters: the aspect ratio of the enclosure, A, the generalized Prandtl number, Pr, and the generalized Rayleighnumber, Ra, which are, respectively, defined as

A ¼ L0

H 0; Pr ¼ ðk=qÞH

02�2n

a2�nand Ra ¼ gbH 02nþ2q0

ðk=qÞankð10Þ

It is easy to verify that the classical expressions of Pr and Ra, corresponding to Newtonian fluids, are recov-ered for n = 1, provided that k is replaced by the Newtonian viscosity l.

3. Numerics

The two dimensional governing equations are solved by using the well known second order central finitedifference method with a regular mesh size. The integration of the vorticity and energy equations, Eqs. (2)and (3), is performed with the alternating direction implicit method (ADI). This method, frequently usedfor Newtonian fluids, was also successfully used for non-Newtonian power law fluids by Ozoe and Churchill[9], Turki [11] and Naımi [14]. To satisfy conservation of mass, the stream function equation, Eq. (4), wassolved by a point successive over relaxation method (PSOR) with an optimum relaxation factor calculatedby the Franckel formula, which is given in the book by Roache [15]. The calculations were performed on aPentium 4 (Processor: Intel · 86 Family 15 Model 1 Stepping 2–1.6 GHz, Memory: 128 MB RAM). As indi-cated in connection with Table 1 concerning the computation time versus accuracy, Table 2 summarizes theCPU (Central Processing Unit) time against grid size. As it can be seen, while starting the numerical code withthe same initial conditions, an accurate steady state solution requires relatively excessive CPU time. The pro-cedure followed consists of refining the grid size until the analytical results are recovered with very satisfactory

1rgence tests for A = 8 and various values of Ra and n

(161 · 41) (241 · 41) (201 · 41) (201 · 33) (201 · 49)

Ra wmax Nu wmax Nu wmax Nu wmax Nu wmax Nu

103 �2.923 4.541 �2.961 4.523 �2.956 4.526 �2.864 4.421 �2.942 4.543104 �7.629 25.885 �7.659 25.868 �7.653 25.878 �7.553 24.446 �7.634 25.98105 �17.234 142.04 �17.362 143.12 �17.31 142.75 �16.442 141.23 �17.256 142.88

103 �2.016 2.635 �2.029 2.621 �2.024 2.627 �2.035 2.523 �2.026 2.623104 �5.200 12.298 �5.215 12.271 �5.210 12.2919 �5.243 12.142 �5.225 12.284105 �11.801 62.112 �11.856 61.548 �11.874 61.694 �12.014 60.445 �11.899 61.651

103 �1.447 1.849 �1.461 1.828 �1.450 1.839 �1.441 1.820 �1.451 1.838104 �3.772 6.904 �3.781 6.888 �3.778 6.893 �3.768 6.689 �3.778 6.892105 �8.440 30.713 �8.456 30.723 �8.446 30.692 �8.412 30.225 �8.444 30.701

103 �1.012 1.609 �1.019 1.595 �1.017 1.599 �1.008 1.415 �1.017 1.586104 �2.821 4.276 �2.825 4.267 �2.824 4.269 �2.805 4.096 �2.824 4.256105 �6.118 16.864 �6.128 16.798 �6.124 16.823 �6.101 16.325 �6.125 16.817

103 �0.753 1.216 �0.781 1.236 �0.779 1.258 �0.742 1.189 �0.776 1.251104 �2.161 2.938 �2.207 2.952 �2.205 2.971 �2.173 2.758 �2.195 2.962105 �4.581 9.975 �4.669 10.032 �4.671 10.036 �4.443 9.958 �4.662 10.041

Table 2CPU time versus accuracy

Grids (161 · 41) (241 · 41) (201 · 41) (201 · 33) (201 · 49)

N CPU time (s)

0.6 692.44 1061.22 857.87 649.88 1070.131.0 284.95 544.75 400.07 268.43 551.441.4 4542.66 6864.12 5654.05 4288.33 12474.62


accuracy. Hence, a uniform grid of 201 · 41 was selected for A = 8 (value used for the numerical computa-tions). It was estimated to be sufficient to model accurately the fluid flow and temperature distribution insidethe cavity. To satisfy the continuity equation, the convergence criterion

Pi;jjw

kþ1i;j � wk

i;jj < 10�4Pi;jjw

kþ1i;j j was

adopted, where wki;j is the value of the stream function at the kth iteration level. The time step size, dt, was

varied in the range 10�76 dt 6 10�4, depending on the values of the governing parameters. More precisely,

the smallest values of dt were used for high values of n and Ra. In the present study, the relation of Woodswas selected among the various formulae proposed in the literature to evaluate the vorticity on the rigidboundaries (see for instance [15] for more details).

With the Ostwald power law model, the dimensionless viscosity given by Eq. (6) tends toward infinity at theimmediate vicinity of the corners of the cavity where the velocity gradients are very small, which renders directnumerical computations impossible. This numerical difficulty is, however, surpassed by using average valuesfor the corner viscosity, which makes the computations possible and stable.

The local heat transfer through the fluid layer filling the cavity can be expressed in terms of local Nusseltnumber, defined as

NuðyÞ ¼ q0=ðkDT 0=L0Þ ¼ A=DT ¼ 1=ðDT=AÞ ð11Þ
where DT = T(0, y) � T(A,y) is the side to side dimensionless local temperature difference. Because of the edgeeffects, this definition is notoriously inaccurate owing to the uncertainty of the temperature values evaluated atthe two vertical walls. Instead, the Nusselt number is calculated on the basis of a temperature difference be-tween two vertical sections, far from the end sides. Thus, by analogy with Eq. (11) and considering two infin-itesimally close sections, the local Nusselt number can be defined by
NuðyÞ ¼ limdx!0

dx=dT ¼ limdx!0

1=ðdT=dxÞ ¼ �1=ðoT=oxÞx¼A=2 ð12Þ

where dx is the distance between two symmetrical sections with respect to the central one (section at x = A/2).The mean Nusselt number is calculated at different locations as follows:

Nu ¼Z 1

0

NuðyÞdy ð13Þ

The validity of the above remarks, in connection with Eqs. (11) and (12), is confirmed by the results pre-sented in Table 3. In fact, as it can be seen from this table, while moving away from the side edges, Nu tends tobe constant, which justifies the approach adopted for evaluation of this quantity. The values of Nu, calculatedby integrating Eq. (11) (dx = A) and Eq. (12) (dx! 0), differ notably, and the corresponding relative differ-ences increase by increasing Ra and decreasing n.

As an additional check of the accuracy of the results, an energy balance was systematically verified for thesystem at each numerical code running. Thus, the overall heat transfer through each vertical plane was eval-uated and compared with the quantity of heat furnished to the system at x = 0. For the results reported here,the energy balance was satisfied within 2% as a maximum difference.

Typical numerical results, in terms of streamlines and isotherms, are presented in Fig. 2 for A = 8 and dif-ferent values of n and Ra. It is seen from this figure that the flow is parallel to the horizontal boundaries in thecentral part of the enclosure and the temperature is linearly stratified in the horizontal direction for the con-sidered values of n and Ra. These remarks are at the origin of the assumptions allowing appropriate simpli-fications useful to develop an approximate analytical solution.

Table 3Effect of the side containment on the mean Nusselt number for A = 8 and various values of n and Ra

n Ra Nu

dx = A dx = 3A/4 dx = 2A/3 dx = A/2 dx = A/3 dx! 0

0.6 103 3.757 4.519 4.514 4.522 4.523 4.526104 14.217 25.305 25.427 25.683 25.808 25.878105 42.109 143.352 143.478 144.779 144.900 142.752

0.8 103 2.335 2.633 2.624 2.623 2.622 2.627104 8.273 12.271 12.259 12.284 12.294 12.291105 24.796 61.141 61.317 61.749 61.761 61.694

1.0 103 1.703 1.839 1.832 1.832 1.833 1.839104 5.243 6.905 6.888 6.892 6.893 6.893105 15.538 30.455 30.483 30.608 30.666 30.692

1.2 103 1.463 1.561 1.563 1.572 1.581 1.599104 3.537 4.270 4.279 4.269 4.249 4.269105 10.157 16.689 16.792 16.804 16.785 16.823

1.4 103 1.221 1.256 1.259 1.257 1.254 1.258104 2.603 2.973 2.978 2.973 2.959 2.971105 6.921 9.993 10.033 10.016 9.992 10.036


4. Approximate analytical solution

On the basis of the results presented in Fig. 2, the following simplifications are considered:

uðx; yÞ ¼ uðyÞ; vðx; yÞ ¼ 0; wðx; yÞ ¼ wðyÞ and T ðx; yÞ ¼ Cðx� A=2Þ þ hðyÞ ð14Þ
where C is the unknown, but constant, dimensionless horizontal temperature gradient in the core region.Using these approximations, the simplified resulting non-dimensional governing equations become
d2

dy2

dudy

��n�1

dudy

" #¼ CRa ð15Þ

Cu ¼ o2Toy2¼ d2h

dy2ð16Þ

with the following boundary conditions:

u ¼ dh=dy ¼ 0 for y ¼ 0 and 1 ð17Þ
and the return flow condition Z 1
0

uðyÞdy ¼ 0 ð18Þ

In the past, the above concept has been successfully used by many authors to predict the flow behavior inthe case of a pure fluid (see for instance Cormack et al. [1]) and fluid saturated porous media (see for instanceBejan [2]).

The integration of Eqs. (15) and (16), coupled with conditions (17) and (18), leads to analytical expressionsfor velocity and temperature. However, it is mentioned that the integration operation is complicated due to thenature of the governing equations and requires a special numerical treatment. Indeed, the non-linearity of thefluid behavior and the change of the velocity gradient sign due to the return flow imposes that the velocityexpressions are different depending on whether 0 6 y 6 y0, y0 6 y 6 y1 or y1 6 y 6 1, where y0 and y1

(y1 = 1 � y0, because of the centro-symmetry of the core flow) are the vertical co-ordinate values at whichthe vertical velocity gradient is zero. They are derived from Eq. (18), which is numerically solved by usinga combination of the Regula–Falsi and Wegstein iteration methods [16] and the Gauss–Legendre integrationmethod [17]. To reduce the velocity and temperature expressions, the function f(y) = (y2 � y + y0y1)/2 is intro-duced. Thus, for 0 6 y 6 y0, one obtains

(a)

(b)

(c)

(d)

(e)---------(1)---------

(a)

(b)

(c)

(d)

(e)---------(2) ---------

(a)

(b)

(c)

(d)

(e)---------(3)---------

Fig. 2. Streamlines (left) and isotherms (right) for A = 8 and various values of n ((a) n = 0.6, (b) n = 0.8, (c) n = 1.0, (d) n = 1.2 and (e)n = 1.4) and Ra ((1) Ra = 103, (2) Ra = 104 and (3) Ra = 105).



uðyÞ ¼ C1=nRa1=n

Z y

0

½f ðyÞ�1=n dy� �

ð19Þ

hðyÞ ¼ C1þ1=nRa1=n

Z y

0

Z y

0

Z y

0


dy� �

dy� �

þ hð0Þ ð20Þ

For y0 6 y 6 y1, the obtained expressions are

uðyÞ ¼ C1=nRa1=n

Z y0

0

½f ðyÞ�1=n dy þZ y0

y½�f ðyÞ�1=n dy

� �ð21Þ

hðyÞ ¼ C1þ1=nRa1=n ðy � y0Þ2

2

Z y0

0

½f ðyÞ�1=n dy þZ y

y0

Z y

y0

Z y0


� �dy

" #dy

"

þðy � y0ÞZ y0

0

Z y

0


dy þZ y0

0

Z y

0

Z y

0


dy� �

dy

#þ hð0Þ ð22Þ

Finally, for y1 6 y 6 1, the expressions are as follows:

uðyÞ ¼ C1=nRa1=n

Z y0

0


y1

½�f ðyÞ�1=n dy þZ y

y1

½f ðyÞ�1=n dy

" #ð23Þ

hðyÞ ¼ C1þ1=nRa1=n 1

2ðy � y1Þðy þ y1 � 2Þ

Z y0

0


y1

½�f ðyÞ�1=n dy

" #"

þZ y

y1

Z y

1

Z y

y1

½f ðyÞ�1=n dy

" #dy

" #dy þ

Z y1

y0

Z y

y0

Z y0


� �dy

" #dy

þ 1

2ðy1 � y0Þ

2

Z y0

0

½f ðyÞ�1=n dy þ ðy1 � y0ÞZ y0

0

Z y

0


dy

þZ y0

0

Z y

0

Z y

0


dy� �

dy

#þ hð0Þ ð24Þ

The expression for w(y) can be deduced from that of u(y) by integration of Eq. (5), taking into account the asso-ciated boundary conditions. As a result, the stream function at the center of the enclosure can be expressed by

wc ¼ wðA=2; 1=2Þ

¼ C1=nRa1=n ð1=2� y0ÞZ y0

0


0

Z y

0


dy þZ 1=2

y0

Z y0


� �dy

" #

ð25Þ
The constant h(0), which is determined while exploiting the centro-symmetry of the core temperature field,
is given by the following expression:

hð0Þ ¼ �C1þ1=nRa1=n ð1=2� y0Þ2

2

Z y0

0

½f ðyÞ�1=n dy þZ 1=2

y0

Z y

y0

Z y0


� �dy

" #dy

"

þð1=2� y0ÞZ y0

0

Z y

0


dy þZ y0

0

Z y

0

Z y

0


dy� �

dy

#ð26Þ

The constant C characterizes the axial temperature gradient. Its value is obtained from the thermal bound-ary conditions imposed on the end walls. Because of the turning flow at the end regions of the fluid layer, theboundary conditions in the x-direction (Eq. (8)) can not be satisfied by the parallel flow approximation.Instead, the expression for the constant C was determined by matching the core solution (Eq. (14)) to the inte-gral solution for the end regions (integration of Eq. (3), together with boundary conditions (8) and (9), by con-sidering the arbitrary control volume of Fig. 1). This yields the following equation:

C þ 1 ¼Z 1

0

uðyÞhðyÞdy ð27Þ

Table 4Dependence of y0 and An on n

n y0 An

0.6 0.199 �0.485 · 10�7

0.8 0.206 �0.599 · 10�6

1.0 0.211 �0.276 · 10�5

1.2 0.216 �0.768 · 10�5

1.4 0.219 �0.160 · 10�4


By substituting the expressions of u(y) and h(y) into Eq. (27), we obtain:

C þ 1 ¼ AnRa2=nC1þ2=n ð28Þ
where the coefficient An, which depends only on n, is calculated with the Gauss–Legendre integration method,and its values are presented with those of y0 in Table 4. It should be pointed out, from this table, that y0 is anincreasing function of n, which means that the velocity maximum is displaced away from the lower wall byincreasing n. On the other hand, for a given n, there is no dependence between y0 and Ra, and the latter seemsnot to affect the shape of the velocity profile.
To determine the value of C for given values of n and Ra, Eq. (28) was firstly solved by the Regula–Falsiiteration method. Then, the resulting solution is used as an initial estimation to solve the same equationwith the Wegstein iteration method (which is of higher order than the former) in order to improve theaccuracy.

Taking into account Eqs. (12) and (13), the Nusselt number is constant, and it can be expressed as

Nu ¼ �1=C ¼ Nu ð29Þ

Thus, for Ra = 0, which corresponds to the pure conductive regime, Eq. (28) gives C = (oT/ox)x=0,A = �1and, according to Eq. (29), Nu ¼ Nu ¼ 1.

5. Results and discussion

The thermal boundary conditions imposed (uniform heat flux) lead to flow characteristics independent ofthe aspect ratio A when this parameter is large enough. The approximate solution, developed in the precedingsection on the basis of the parallel flow assumption, thus is valid asymptotically in the limit of a shallow cavity(A� 1). Therefore, numerical tests are performed to determine the smallest value of A leading to results rea-sonably close to those of the large aspect ratio approximation. The natural convection flow developed insidethe enclosure is controlled by the Rayleigh number, Ra, and the power law index n, which is varied here from0.6 to 1.4 to cover shear thinning (0 < n < 1), Newtonian (n = 1) and shear thickening (n > 1) fluids. It is nec-essary to point out that, contrary to the case where the cavity is heated from below for which the onset ofconvection occurs when some threshold of Ra is reached [7], the convection phenomenon is present in the caseof lateral heating at any value of Ra different from zero. For the non-Newtonian fluids considered here, thePrandtl number is large enough (Pr!1) that the convection heat transfer becomes insensitive to any changeof the large values of this parameter [14,18]. This finding is also confirmed by the present analytical solution,which is independent of Pr in its range of validity. A parametric study is, thus, required to determine the Pr

values corresponding to the fluids concerned by the present study and, at the same time, leading to significantreductions in terms of computation time.

5.1. Determination of the asymptotic Pr as a large non-Newtonian Prandtl number value

The generalized Prandtl number, Pr, considered here does not represent an intrinsic property of the fluidbecause of its dependence on H 0 and n. However, this parameter is reduced to its conventional expression cor-responding to a Newtonian fluid (n = 1). Unfortunately, it turns out that, for essentially all non-Newtonianfluids of interest, Pr is much larger than unity which requires an excessive calculation time since, as reported

Table 5Effect of the Prandtl number on the flow intensity and heat transfer rate forA = 8 and various values of n and Ra

N Pr Ra = 103 Ra = 104 Ra = 105

wmax Nu wmax Nu wmax Nu

0.6 1 �2.948 4.627 �7.842 27.199 �20.889 325.79910 �2.939 4.535 �7.646 25.868 �16.783 144.396100 �2.950 4.524 �7.651 25.879 �17.267 142.3491 �2.956 4.526 �7.653 25.878 �17.310 142.752

0.8 1 �2.006 2.646 �5.216 12.269 �11.783 64.04610 �2.021 2.621 �5.203 12.304 �11.851 61.754100 �2.021 2.622 �5.200 12.3106 �11.861 61.7191 �2.024 2.627 �5.210 12.291 �11.874 61.694

1.0 1 �1.442 1.832 �3.782 6.891 �8.423 30.71610 �1.444 1.834 �3.779 6.895 �8.426 30.720100 �1.446 1.836 �3.776 6.892 �8.430 30.7041 �1.450 1.839 �3.778 6.893 �8.446 30.692

1.2 1 �1.049 1.437 �2.820 4.269 �6.163 16.69310 �1.019 1.591 �2.825 4.271 �6.121 16.839100 �1.018 1.595 �2.826 4.269 �6.122 16.8381 �1.017 1.599 �2.824 4.269 �6.124 16.823

1.4 1 �0.727 1.208 �2.199 2.934 �4.698 9.87910 �0.787 1.258 �2.201 2.964 �4.657 10.024100 �0.784 1.260 �2.200 2.968 �4.661 10.0291 �0.779 1.258 �2.205 2.971 �4.671 10.036


by Ozoe and Churchill [9], the system reaction is a decreasing function of Pr. On this basis, trial calculationswere performed for different values of this parameter (Pr = 1, 10, 100 and Pr!1) to determine the Pr valuethat will avoid excessive time of computation and, at the same time, represent fluids with large Prandtl num-bers. The results presented in Table 5 show that the increase of Pr from 1 to 10 is followed by an importantdecrease in terms of wmax and Nu for small n (n = 0.6) and large Ra (Ra = 105), but any increase of Pr from 10leads only to a limited effect on the fluid flow and overall heat transfer characteristics (Table 5). This can beeasily deduced from Eq. (2) in which the convective terms (those on the left hand side) vanish by increasing Pr

when compared to the diffusive ones (terms on the right hand side of the same equation). Hence, the governingparameters of the study are restricted to Ra and n, and the results obtained are valid in the limit Pr!1. Asindicated in Table 6, this presents the advantage of making the computations faster than those conducted withfinite values of Pr (i.e., for Pr 6 100) for which the convective terms (also called inertial terms) are present.

5.2. Search for A satisfying the large aspect ratio approximation

In this paragraph, the objective consists of determining the lower value of the aspect ratio A that lead toconvection heat transfer results in agreement with those obtained analytically. Hence, the effect of A on wc

and Nu, evaluated numerically, is illustrated in Fig. 3a and b for Ra = 104 and various values of n. It is seenfrom these figures that an increase of A leads to an asymptotic behavior of wc and Nu. For all the values of n

considered, the asymptotic analytical limits are largely reached for A = 8. The numerical values of wc and Nu,

Table 6Example of the dependence of the CPU time and the time step on the Prandtl number

Pr 1 10 100 1N dt CPU time (s) dt CPU time (s) dt CPU time (s) dt CPU time (s)

0.6 10�4 935.78 10�5 2283.28 10�6 13195.87 10�4 857.871.0 10�4 592.37 10�5 1771.07 10�6 10629.03 10�4 400.071.4 10�5 5900.18 10�6 17944.92 10�7 68351.70 10�5 5654.05

-7.8

-5.2

-2.6

0.0

0.8

1.0

1.2

1.4

n = 0.6

(a) Parallel flow solution Numerical solution

Parallel flow solution Numerical solution

ψc

1 2 3 4 5 6 7 8 9 100

15

30

45

0.8

1.4

1.21.0

n = 0.6

(b)

Nu

A

1 2 3 4 5 6 7 8 9 10A

Fig. 3. Evolution of the stream function (a) and the Nusselt number (b) in the central part of the cavity with the aspect ratio for Ra = 104

and various values of n.


obtained with this value of A and depicted in Fig. 6, are seen to be in good agreement with the analyticalresults. Such a limit of A appears to be sufficient for the values of Ra shown in these figures and justifiesits use, at least for the numerical calculations performed in the present study. On the other hand, the peaksobserved in Fig. 3b for n = 0.6 and 0.8 are not attributed to a change in the flow structure. In fact, the exam-ination of the evolutions of the streamlines and isotherms with A (not presented here) showed that the flowstructure remains unicellular. According to Eq. (11), these maxima are probably related to the fact that theshear thinning behavior acts in such a way that the increase of the overall temperature gradient is slower thanthat of A in the case of significant end-side walls effect (1 6 A 6 1.4) and more rapid just beyond A = 1.4(value corresponding to the peaks). This tendency is not observed for n = 1.2 and 1.4 due to the reducing effectof the shear thickening behavior on natural convection heat transfer.

5.3. Validation of the analytical solution

To validate the approximate analytical solution, the numerical results (full circles), obtained by solving thefull governing equations, are compared against those obtained analytically (solid lines) for A = 8, Ra = 104

and different values of n. Thus, Fig. 4 shows an excellent agreement between the numerical and analyticalresults presented in terms of the profiles of velocity (a), stream function (b) and temperature (c), calculatedat the mid-length of the cavity along the vertical direction. In addition, the profiles of stream function andtemperature along the horizontal direction, obtained with both approaches at the mid-height of the cavity,are presented, respectively, in Fig. 5a and b. The excellent agreement observed in the range 1 < x < 7, i.e.,

-24

-12

0

12

24

0.0 0.2 0.4 0.6 0.8 1.0

0.8

1.21.0

1.4

n = 0.6

(a)

(b)

(c)

Parallel flow solutionNumerical solution



y

0.0 0.2 0.4 0.6 0.8 1.0y

0.0 0.2 0.4 0.6 0.8 1.0y

u

-8

-6

-4

-2

0

1.0

0.8

1.2

1.4

n = 0.6

-0.210

-0.105

0.000

0.105

0.2101.4

1.21.0

0.8n = 0.6

T

ψ

Fig. 4. Horizontal velocity (a), stream function (b) and temperature (c) profiles at mid-length of the cavity, along the vertical co-ordinatefor A = 8, Ra = 104 and various values of n.


far from the vertical confining sides, confirms the validity of the parallel flow hypothesis for the consideredvalues of the parameter n. Finally, computed and calculated values of the stream function at the center, wc,and mean Nusselt number, Nu, presented in Fig. 6a and b, show also good agreement between the analyticaland numerical results for a wide range of Rayleigh number, Ra, and various values of the power law index, n.

5.4. Combined effects of the power law index and Rayleigh number

Typical numerical results in terms of streamlines (left) and isotherms (right) are presented in Fig. 2 forA = 8 and various values of n and Ra. Their examination gives an idea about the conjugate effect of the powerlaw index, n, and the Rayleigh number, Ra, on the flow structure and temperature distribution inside thecavity. The unicellular nature of the flow is always preserved for n and Ra varying in their respective ranges.Globally, if we except the end regions, the flow remains parallel to the long walls, and the distribution of thetemperature inside the cavity is characterized by a linear stratification in the horizontal direction. Theseobservations allow important simplifications that were used to transform the governing partial differential

-8

-6

-4

-2

0

0.8

1.4

1.2

n = 0.6

1.0

0 2 4 6-1.70

-0.85

0.00

0.85

1.70

0.8

n = 0.6

1.4

1.2

1.0

(b)

(a)



T

ψ

x8

0 2 4 6x

8

Fig. 5. Stream function (a) and temperature (b) profiles at the mid-height of the cavity, along the horizontal co-ordinate for A = 8,Ra = 104 and various values of n.


equations to ordinary ones. Qualitatively, the isotherms seem to be more sensitive to the variations of Ra or nthan the streamlines. Their inclination with respect to the vertical direction (that of the conductive regime) ismore and more marked by enhancing the fluid circulation (increase of Ra or a decrease of n). These qualitativeobservations are consistent with the flow intensity and heat transfer rates, which are presented in Table 1 interms of wmax and Nu. Indeed, for a given value of Ra, the decrease of these quantities (in absolute value)accompanying the decrease of n is mainly attributed to the diminution of the generalized shear rate. Accordingto Eq. (6), this leads to an increase of the apparent viscosity whose slowing down role of the fluid motion iswell known. In terms of computing time, the retarding character of the non-Newtonian behavior to achieveconvergence towards the steady state solution is illustrated in Table 2. It is seen from this table that the devi-ation from the Newtonian behavior is accompanied by an increase of the CPU time, and this increase is moreimportant in the case of the shear thickening fluids (n > 1). The non-linearity of the vorticity transport equa-tion (Eq. (3)) augments in the case of fluids having non-Newtonian behavior, which makes the calculation timelonger before reaching the steady state. On the other hand, for a given value of n, both wmax and Nu increasewith Ra, due to the obvious contribution of the buoyancy effects in promoting the convection. In this respect,Table 7 indicates that the CPU time is a decreasing function of Ra, whereas the time step is invariable. Theexamination of Fig. 6a and b, presenting, respectively, the evolution of wc and Nu with Ra for various values ofn, shows clearly that the Rayleigh number and the power law index have opposite effects on the flow charac-teristics and heat transfer. In fact, for a given n, an increase of Ra leads to enhancing the convection effect, anda similar behavior is induced by decreasing n for a given Ra. The convection is, however, more sensitive to the

102 103 104 105 106

102 103 104 105 106

-50

-40

-30

-20

-10

0

(a)

(b)

0.8

n = 0.6

1.0

1.4

1.2

Parallel flow solution Numerical solution

ψc

Ra

100

101

102

103

1.4

1.2

0.8

1.0

n = 0.6 Parallel flow solution Numerical solution Correlating relation Eq. (30)

Nu

Ra

Fig. 6. Evolution of the stream function (a) and the Nusselt number (b) in the central part of the cavity with the Rayleigh number forA = 8 and various values of n.

Table 7Example of the dependence of the CPU time and the time step on the Rayleigh number

Ra 103 104 105

n dt CPU time (s) dt CPU time (s) dt CPU time (s)

1.0 10�4 1900.13 10�4 400.07 10�4 287.56


change of Ra for the shear thinning behavior (0 < n < 1) compared to the shear thickening one (n > 1). Similarobservations concerning the influences of n and Ra on the convection heat transfer were reported in the pastby Turki [11] while investigating natural convection in a power law non-Newtonian fluid filled cavity heatedfrom the vertical side with a constant temperature.

To render the obtained results more useful, the heat transfer results are correlated in terms of Nu versusadequate combinations, including n and Ra, in the following mathematical form:

Nu ¼ aðnÞRabðnÞ ð30Þ
The expressions for the coefficient a(n) and the exponent b(n) given in Table 8 are found to depend on theranges of Ra (ranges of validity) given also in this table. The results obtained with Eq. (30) are presentedin Fig. 6b with dashed lines. They are seen to be in good agreement with both the parallel flow solutionand the numerical results with a maximum difference not exceeding 3%. Note that the correlation given by

Table 8Expressions of a(n) and b(n) for different ranges of Ra

1026 Ra < 103 103

6 Ra < 104 1046 Ra 6 106

a(n) 0.271n2.791 0.093n2 � 0.132n + 0.070 �0.009ln(n) + 0.018b(n) 0.511n2 � 1.673n + 1.425 �0.037n2 � 0.366n + 0.985 0.092n2 � 0.429n + 0.986


Eq. (30) could be of some interest in thermal processing applications involving non-Newtonian shear thinningor shear thickening fluids in the presence of the natural convection phenomenon.

6. Conclusion

The present paper is devoted to numerical and analytical experiments on natural convection in a twodimensional horizontal shallow enclosure (A� 1) filled with non-Newtonian power law fluids in the casewhere both short vertical sides are submitted to uniform heat flux (Neumann type condition), while the hor-izontal boundaries are insulated.

The full partial differential equations governing the problem are integrated numerically using a finite dif-ference method. The computations were conducted for various values of A (1 6 A 6 10), n (0.6 6 n 6 1.4),Pr (Pr = 1, 10, 100 and1) and Ra (0 6 Ra 6 106). The main conclusions of the present investigation are sum-marized on the following points:

• The fluid flow and heat transfer characteristics are found to be insensitive to any increase of the aspectratio, A, and the Prandtl number, Pr, when these parameters are large enough (A P 8 and Pr P 100). Thisallows concluding that for the non-Newtonian fluids considered in this study, which are characterized byhigh values of Pr, natural convection in shallow cavities is mainly controlled by the flow behavior index,n, and the Rayleigh number, Ra.

• The fluid flow and heat transfer characteristics appear to be rather sensitive to the non-Newtonian powerlaw behavior. Compared to Newtonian fluids (n = 1), a shear thinning behavior (0 < n < 1) enhances thefluid circulation and the convection heat transfer, while the shear thickening behavior (n > 1) producesan opposite effect.

• The approximate analytical solution developed on the basis of the parallel flow concept is found to agreewell (in the core region of the cavity) with the numerical solution obtained by solving the full governingequations.

References

[1] Cormack DE, Leal LG, Imberger J. Natural convection in a shallow cavity with differentially heated end walls. Part 1. Asymptotictheory. J Fluid Mech 1974;65:209–29.

[2] Bejan A. The boundary layer regime in a porous layer with uniform heat flux from the side. Int J Heat Mass Transfer1983;26:1339–46.

[3] Bian W, Vasseur P, Bilgen E. Natural convection of non-Newtonian fluids in an inclined porous layer. Chem Eng Commun1994;129:79–97.

[4] Amari B, Vasseur P, Bilgen E. Natural convection of non-Newtonian fluids in a horizontal porous layer. Warme-undStoffubertragung 1994;29:185–93.

[5] Amahmid A, Hasnaoui M, Mamou M, Vasseur P. Double-diffusive parallel flow induced in a horizontal Brinkman porous layersubjected to constant heat and mass fluxes: analytical and numerical studies. Heat Mass Transfer 1999;35:409–21.

[6] Mamou M, Vasseur P, Bilgen E. Analytical and numerical study of double-diffusive convection in a vertical enclosure. Heat MassTransfer 1996;32:115–25.

[7] Mamou M, Vasseur P, Hasnaoui M. On numerical stability analysis of double-diffusive convection in confined enclosures. J FluidMech 2001;433:209–50.

[8] Jaluria Y. Thermal processing of materials: from basic research to engineering. J Heat Transfer 2003;125:957–79.[9] Ozoe H, Churchill SW. Hydrodynamic stability and natural convection in Ostwald–De Waele and Ellis fluids: the development of a

numerical solution. AIChE J 1972;18:1196–207.


[10] Tien C, Tsuei HS, Sun ZS. Thermal instability of a horizontal layer of non-Newtonian fluid heated from below. Int J Heat MassTransfer 1969;12:1173–8.

[11] Turki S. Contribution to numerical study of natural and mixed convection heat transfers in confined non-Newtonian fluids, Ph.D.Thesis, CNAM, Paris, France, 1990.

[12] Siginer DA, Valenzuela-Rendon A. On the laminar free convection and stability of grade fluids in enclosures. Int J Heat MassTransfer 2000;43:3391–405.

[13] Gray DD, Giorgini A. The validity of the Boussinesq approximation for liquids and gases. Int J Heat Mass Transfer 1976;19:545–51.[14] Naımi M. Contribution to the thermal Marangoni effect: study in non-Newtonian fluids under gravity and micro-gravity conditions,

State Doctorate Thesis, Cadi Ayyad University, Marrakech, Morocco, 2001.[15] Roache PJ. Computational fluid dynamics. Albuquerque, New Mexico: Hermosa Publishers; 1982.[16] Gourdin A, Boumahrat M. Applied numerical methods, technique and documentation. Paris: Lavoisier; 1989.[17] Sibony M, Mardon J-CI. Numerical analysis II, approximation and differential equations. Paris: Hermann; 1982.[18] Ng ML, Hartnett JP. Natural convection in power law fluids. Int Commun Heat Mass Transfer 1986;13:115–20.

Documents

Natural convection heat transfer in shallow horizontal rectangular enclosures uniformly heated from the side and filled with non-Newtonian power law fluids