Prediction of natural convection in horizontal annular enclosures at high Rayleigh numbers

INT. COMM. HEAT MASS T R A N S F ~ 0735-1933/88 $3.00 + .00 Vol. 15, pp. 785-797, 1988 ~Pergamc~ Press plc Printed in the United States

PREDICTION OF NATURAL CONVECTION IN HORIZONTAL ANNULAR ENCLOSURES AT HIGH R A Y L E I ( ~ H ~

Ashley M. Thornhill, E. Kobla Glakpe Department of Mec_hanical Engineering Howard University Washington, D.C. 20059 and Derej e Agorkafer Data Systems Division IH~4 Poughkeepsie, N.Y. 12602

(Cx~l,iLlnicated by J.P. Hartnett and W.J. Minkcwycz)

ABSTRACT This paper presents ntw0erical results of natural cenvection at high Rayleigh numbers in ar~uli formed by concentric and vertically eccentric horizontal cylinders. Soluticns are obtained in a rmm~rically generated bot~dary-fitted coordinate syst~n in which coordinate lines are coincident with the enclosure ~ i e s . Results obtained reveal the effects of eccentricity and Rayleigh mLmlber on convective heat transfer. Results are presented in the form of vector plots of velocity, local heat transfer data in addition to contours for t e r ~ r a t u r ~ turbulent kinetic energy, and turbulent viscosity. All results presented are for a fluid of Prandtl r~xdoer 0.7 and Rayleigh numbers in the range of 105 to 2xlO 7.

Introduction Turbulent natural convection in enclosed cavities is a p h ~ mendous practical interest.

of tre-

/~mmng the r~x~rous practical applications of tur-

bulent natural convection in enclosures are concentrating solar collector receiver designs, ventilation of rooms, nuclear reactor designs, and electronic packaging.

One area of current interest is the design of shipping casks for

the transportation of spexlt nuclear fuel assemblies from nuclear reactor sites to reprocessing facilities. Relatively little work has been done experimentally and theoretically in the area of tumbulent natural convection.

Kuehn and Goldstein [1 ], and other

investigators cited therein, studied convective heat transfer in ar~nular enclo785

786

A.M. Thornhill, E.K. Glakpe and D. Agonafer

sures.

Vol. 15, No. 6

In particular, Kuehn and Goldstein [1 ] investigated the influence of

Rayleigh number on natural convection heat transfer through a fluid bounded b? horizontal isothermal cylinders.

This study was performed using pressur'-

ized nitrogen to obtain results up to a Rayleigh number of 7.7xi07 for a concentric configuration.

The effect of eccentricity on heat transfer was also

studied by varying the position of the inner cylinder with the Rayleigh nt~nber fixed at 5xlO &.

In all cases, the Prandtl number and the radius ratio were

kept constant at 0.7 and 2.6 respectively. Recently Farouk [2] employed the k-E n~)del of turbulence in a vorticity based f o l i a t i o n

to study turbulent natural convection in the annulus between

horizontal concentric cylinders with isothermal boundaries.

Results for the

turbulent regime were generated using air with Prandtl ntm~oer of 0.7, a radius ratio of 2.6 and for Rayleigh nt~3ers ranging from 106 to 107 . In the present study, the k-~ model equations as en]01oyed by Farouk [2], and by Ozoe et al [3], were transformed and solved with the Navier-Stokes and energy equations simultaneously in a boundary-fitted coordinate system.

Un-

like these investigators, however, the primitive variable form of the goveming

equations is retained.

The boundary-fitted coordinate system approach

was first introduced by Thompson [4 ] and subsequently employed by other investigators [5,6] in obtaining solutions to fluid dynamic problems in irregular gec~etries. A major contribution of this paper is to demonstrate the effectiveness of the boundary-fitted coordinate method in obtaining solutions for problems involving natural convection at high Rayleigh ntrabers. Results generated are for concentric and vertically eccentric a m u l i between horizontal cylinders. Although Kuehn and Goldstein [i ] obtained results for a few eccentric configurations, no heat transfer data exists in the literature for Rayleigh nt~doers beyond 5x]O 4.

Thus, an additional contribution of the paper is to investigate

the influence of vertical eccentricity on natural convection in horizontal annuli for high Rayleigh ntr~oers.

To facilitate comparisons with previous re-

sults for the concentric case, the enclosure is filled with air having a Prandtl ntm%ber of 0.7 and has a diameter ratio of 2.6. Governing Equations In the present study, the BODYFIT-IFE ccn~outer code developed by Chen et al. [5], and modified by Glakpe et al. [6], has been further modified to account for turbulent natural convection.

The use of the k-c model of turbu-

Vol. 15, No. 6

~ I C ~

IN ANNULAR~NKILOSURES

787

lence requires that two additional equations, one for the turbulent kinetic energy k, and one for its rate of dissipation ~, also be solved in the nt~nerically generated coordinate system. Using the following dimensionless variables (an overbar denotes a dJxne_nsional quantity), t = Vo{/Lo,

P= P/Po' p = P/Po' xi = xi/Uo ' gi = gi/go ' k = ~ V o 2,

ui/Vo, T

ui =

: T/To, and E : L ° $/VoS,

the k and ~ equations can be written in a generic non-dimensionalized form as P ~¢ ~

PT ~ui ~ui ~u + P uj ~x--] ~¢ = AI ~-6 [8--~j ( ~xj-- + ~ l)] + A3 p E + ToSPT

8T

1

+ A2 ReF2OT (~--~j g j) + Re

8

[(UT~=

~¢ ] $ + ~) ~ 3

~xj

in which the generic variables and constants assume the following values: A1 = A2 = %

= i for

= k

and A

n

= C

n

E/k (n : 1,2,5) for ¢ : e.

The above equation can be cast in a conservative form and integrated over a control vol~ue to yield the integral form of the equations.

As a stabili-

zing measure, the dependent variable is multiplied by the continuity equation and added to the integral form of the equations.

In addition to the trans-

formed Navier-Stokes and energy equations [5], the above equations describe the turbulent flow structure in the enclosures of interest.

Additional detail.<

can be found in the thesis of Thornhill [7]. Boundary and Initial Conditions Because of the symmetric configuratic~Is considered, the problem is solved in half of the physical domain.

No-flux conditions are inposed on the symme-

tric boundaries of the enclosures. walls.

The no-slip condition is used on the solid

Turbuler~e kinetic energy k is zero at the solid walls, while its rate

of dissipation E becomes undefined.

The e-equation is therefore solved in a

reduced cc~putational dcmain which excludes the wall. for c is defined at a point close to the wall as

ep

Cp 3/4 k3/2 < ~r

The boundary condition

788

A.M. Thornhill, E.K. Glakpe and D. Agonafer

Vol. 15, No. 6

in which Cp is an empirical constant set equal to 0.09; < is the Von Karman constant equal to 0.42 and Ar is the distance between the computational node and the solid wall. Solution Method The modifiedBODYFIT-iFEcode uses a staggered grid arrangement in which density, temperature, and pressure are stored at the center of the control volume, and velocities are stored at the intersection of the grid lines.

In

the present study, k, c and PT are also stored at the center of the control volume. The convective terms are treated using the upwind differencing scheme, and partial derivatives in the transformed domain are evaluated using central differencing and some averaging procedures.

The above procedures yield a

general finite-difference equation for k and ~ in the form: ApCp = ZAi¢ i + S¢o + S¢I + S¢2 + S¢3 The buoyancy term S¢2 can be either positive or negative depending on the physical location of any point within the flow domain.

Large negative values

of S¢2 could occur and cause negative k values in the computational scheme. To ensure a stable computational scheme, the k and E equations were reformulated in a manner similar to that done by Farouk [2]. The reformulated equationbecomes ~Ai()i + (S¢o + S¢5) ~p: -S¢1

(-S~2,,

In the solution scheme the converged results for the laminar case is used as guess values to rt~ the turbulent case. ties k, ~ and ~r are set as follows.

Guess values for turbulent quanti-

The laminar viscosity was used as an

initial guess for the turbulent viscosity.

It was observed that the stability

of the scheme was dependent on the initial value of k. of k yielded negative values for k.

Small initial values

An initial guess value of l0 S was found

to yield satisfactory results at a Rayleigh ntr~ber of 105 .

The guess values

of k and PT were used to set the initial value of ~. Results and Discussion Most of the calculations were performed using a 45 (angular) by 40 (radial) grid system on the IBM 3033 (model S) on the Howard University computing

Vol. 15, No. 6 facility.

~ I G N

IN ANNUIARI~kKZ/3SURES

789

The remaining calculations were done on an IH~4 3090 which is about

two to five times faster than the IBM 303S. ent Rayleigh ~ r s ,

In solving the problem for differ-

the k-E model of turbulence was introduced for Rayleigh

numbers greater than or equal to 105 .

In order to achieve ccnver~ence, it was

necessary to tunder-r~lax all the equations.

Due to the relatively large vari-

ations during initial iterations for k and E, relaxation factors used for k and c were smaller than that used for velocities and enthalpy.

Typical values

of relaxation parameters at high Rayleigh numbers were about O. 1 for k and ~, and O. 3 for enthalpy and velocities. Being interested in the heat transfer results, the convergence criteria were set relative to the change in enthalpy.

Convergence is achieved when

the maxinlzn absolute change of enthalpy became less than 10 -7 .

At this level,

relative changes in velocities, k and ~ were about 0.08%, 0.05% and 0.09% respectively. To prove that the cede can generate reliable results, test calculations were done for the concentric annulus problem with isothermal boundaries. Fig. 1 shows the distribution of local equivalent conductivity on the inner and outer cylinders, for a Rayleigh nLw~ber of 106 in cc~parison with results of Farouk [2].

The present results c c ~

well with those of Farouk [2]

even though the grid spacing that he employed was finer.

The performance of

the coarser grid system can he attributed to the use of the boundary-fitted coordinate system in which bot~K~ary conditions are better represented nt~nerically.

Since results for the concentric configuration has been presented by

Farouk [2] and other investigators, the rest of the presentation in this paper will focus on the eccentric configurations.

Eccentricity defines the dis-

30, r l Inn~ boundlry

/

• outer boundary

20 '

/

lg

0

r

-

,

100

-

200

FIG. 1 Distribution of localequivalent conductivity for concentric configuration: Ra = i0~

790

A.M. Thornhill, E.K. Glakpe and D. Agonafer

placement of the outer cylinder relative to the inner one.

Vol. 15, No. 6 Positive eccentri-

city is defined as the upward vertical displacement of the outer cylinder rela rive to the inner cylinder.

In the following sections, results are presented

for both positive and negative eccentricity configurations, with a stmrmary of the overall effect of eccentricity on natural convection.

Results are pre-

sented for two eccentricities, 0.62Z and -0.652. Positive Eccentricity The general flow pattern in the annulus is such that hot fluid near the inner cylinder rises while the cold fluid near to the outer cylinder falls. The velocity field is shown in Fig. 2 for Rayleigh nu~abers of 2x105 and 2xlO 7 . The velocity vector plots show a general flow pattem, with a region of stagnant fluid near the top of the annulus around which the fluid moves. fluid in the bottom of the annulus is also stagnant.

The

The configuration in

this region resembles that of horizontal parallel flat plates with the upper one heated.

As expected, the magnitude of the velocities increase with the

Rayleigh number. The isotherms shown in Fig. 3 are equally spaced with the one on the outer cylinder being at the lowest t e ~ r a t u r e at the highest.

and the one on the inner cylinder

Isotherms indicate the presence of a thermal plmme above the

inner cylinder where the boundary layer separates.

The t e ~ r a t u r e

profile

near the bottom of the annulus indicates the dc~ninance of conduction heat transfer in that region.

l

,

7 T rI

,,~ '

.

~,..,t

\ -~\

II

,tl..... •

(a)

(b) FIG. 2 Velocity distribution in positive eccentric geometry: (a) Ra = 2xlO 5, (b) Ra = 2xlO 7

Vol. 15, No. 6

OONVECTION IN ANNULARENCLOSURES

(a)

FIG. 3

791

(b)

Distribution of isotherms in positive eccentric geometry: (a) Ra = 2xlO 5, (b) Ra : 2xlO 7

(a)

(b) FIG. 4

Isoplots of turbulence viscosity: e = 0.623; (a) Ra = 2xlO 5, (b) Ra -- 2xlO 7 The profiles for turbulent viscosity ~T shown in Fig. 4, display two regions of significant turbulent viscosity in the upper half of the annulus. The isoplots are equally spaced between the maxint~ and minim~n values of ~T" As the Rayleigh number is increased, a third region of significant turbulent viscosity appears close to the inner cylinder. Isoplots for selected values of turbulent kinetic energy k are shown in Fig. 5 for Rayleigh ~ r s low Rayleigh ~ r s ,

of 2xlO 5 and 2xlO 7 .

The profiles indicate that at

only a small region of the annulus is characterized by

regions of large turbulent kinetic energy.

As the Rayleigh ntm~er is in-

792

A.M. Thornhill, E.K. Glakpe and D. Agonafer

(a)

Vol. 15, No. 6

(b) FIG. 5 Isoplots of turbulence kinetic energy: e = 0.625; (a) Ra = 2x105, (b) Ra = 2x107

-20

¸

/

1~,,,,2E07

/

FIm

w loo

200

FIG. 6 Distribution of local equivalent conductivity on inner and outer cylinders creased, there is a significant increase in the turbulent kinetic energy with most of the annulus being turbulent.

As expected, the bottom of the annulus

is characterized by much lower values of turbulent kinetic energy. Results for local equivalent conductivity on the outer and inner cylinders for Rayleigh nt~nbers of 2x105 and 2x107 are presented in Fig. 6.

Results

for the outer cylinder indicate that keqo is m i n i ~ n at the bottom of the annulus (0 = 0), and m a x i n ~ at the top (8 = 180).

As the Rayleigh ntm~oer is in-

creased, it is observed that the region of fairly constant keqo , near the bottom of the cylinder begins to get smaller.

The ~xirmxn value of keqo is

also observed to increase at a slower rate with increasing Rayleigh numbers. On the inner cylinder, the local equivalent conductivity keqi is fairly constant near the bottcm of the annulus.

Closer to the top of the annulus,

keqi increases to a n~xirman and then decreases to a local mininun at the top

Vol. 15, No. 6

CONVECTION IN A N N U I A R ~ C L O S U R E S

793

of the annulus.

As the Rayleigh number is increased, the initial value of keq i

also increases.

The region of constant keqi also becomes smaller and the max-

value of keqi increases, with the position of the m a x ~

remaining about

the same. Negative Eccentricity The vector plots for the negative eccentric case (Fig. 7) reveal the same general flow pattern as the positive eccentric case.

In this case, however,

there is an enlarged stagnant region at the bottom of the annulus.

The stag-

nant region about which the fluid moves does not occur as high in the annulus as is the case for the positive eccentric case.

The velocity increases with

increasing Rayleigh ntl~Iber. The isotherms shown in Fig. 8 display the temperature inversion pattern characteristic of convection. creasing Rayleigh numbers.

The inversion becomes more pronounced with in-

The steep gradients near the bottom of the inner

cylinder and near the top of the outer cylinder are characteristic of the high levels of heat transfer in these regions.

The thermal phx~e at the top of the

inner cylinder is not pronounced in this case because of the narrow gap between the two cylinders in this region.

The profiles near the bottom of the

annulus resemble that of conduction. In the negative eccentric case, there is only one region of significant turbulent viscosity.

Fig. 9 shows this region to be close to the top, and it

is observed that the position of the m a x i m ~ value of ~T does not change appreciably as the Rayleigh ntrd0er is increased.

The turbulent kinetic energy

profiles also show one region of high turbulent kinetic energy.

At a Rayleigh

ntmi~r of 2xlO 5 the turbulent kinetic energy does not exceed i0,000 which indicates that the annulus is still mostly laminar.

At a Rayleigh nt~ber of 106

there is a region of significant turbulent kinetic energy in the upper part of the annulus (Fig. i0).

The distributions of local equivalent conductivity on

the outer and inner cylinders for Rayleigh numbers of 2x105 and 106 are shown in Fig. Ii.

The local equivalent conductivity on the outer cylinder keqo is

m i n i m ~ at the bottcm of the annulus and increases to a m a x i n ~ at about 0 = 120 ° .

Above 120 ° , keqo begins to drop indicating the reduced effect due

to convection as the gap between the two cylinders becomes smaller.

As the

Rayleigh ntm~oer is increased keqo also increases but the angular position of the maxJxm~n remains about the same.

It is also observed that as the Rayleigh

number is increased, keqo close to the top of the annulus does not decrease as fast as at lower Rayleigh nt~nbers. This is because of the increased effect of

794

A.M. Thornhill,

E.K. Glakpe and D. Agonafer

Vol. 15, No. 6

! "''!! ! !'*'.~R ~,, ~ °~ : 1 ' . - . " , ,,,

: . ....-.

,

,.

(a)

i~

(b) FIG. 7 Velocity distribution in negative eccentric configuration: (a) Ra = 2xlO 5, (b) Ra = 10 6

(a) FIG. 8 Distribution of i s o t h e ~ s in negative eccentric geometry: (a) Ha : 2x10 , (b) Ha : 10 6

(a)

r±t/. 9 Isoplots of turbulence viscosity for and (a) Ra = 2xlO 5, (b) Ra = 106

e~-'0.652

Vol. 15, NO. 6

(I)NVECYION IN A N ~ k K I L O S U R E S

795

buoyancy on convection. Local equivalent conductivity on the inner cylinder keqi is maxirm~ at the bottom of the annulus and minirsmn near the top.

Initially, there were some

regions close to the top and bottcrn of the annulus where keg i either decreased relatively slowly or remained constant. leigh n t ~ r

was incr£~sed.

These regions disappeared as the Ray-

The minim~n value of keg i also occured closer and

closer to the top of the annulus.

The m a x ~

value of keqi was also observed

to increase as the Rayleigh ntrdoer was increased. Effect of Eccentricity Results obtained indicate that the position of the inner cylinder has a significant effect on the heat transfer rates in horizontal annuli. lationship between the average Nusselt ~ r ous eccentricities is shown in Fig. 12.

The re-

and the Rayleigh ntm%ber for vari-

One observes that the heat transfer

rate is highest for the positive eccentric case making this configuration most suitable for cooling applications involving high rates of convective heat transfer.

The small gap at the top of the annulus in the nega£ive eccentric

case restricts the movement of the fluid and causes lower levels of convective heat transfer. A comparison of the distributions of the local equivalent conductivity on the inner cylinder keqi for various eccentricities is presented in Fig. 13 for a Rayleigh nt~ber of 106 . At the bottc~ of the annulus, the heat transfer rate is lowest for the positive eccentric case and highest for the negative eccentric case.

The opposite is true at the top of the annulus (e = 180 ° ).

This observation clearly illustrates the effect of the distance between the two cylinders on convective heat transfer.

At the bottom of the annulus, the

cylinders are closer together for the positive eccentric case than they are for the negative eccentric case.

One would therefore expect that the heat

transfer rate at the bottom of the annulus would be higher for the negative eccentric case when compared with the positive eccentric configuration.

For

the top of the annulus, the cylinders are further apart for the positive eccentric case and a higher rate of heat transfer results in this case. Conclusions The ntxnerical results obtained for natural convection of air in horizontal concentric and eccentric annuli with isothermal boundaries have been presented. The positive eccentric case exhibits stronger convective effects cenpared to the negative eccentric case.

Increasing the Rayleigh ntraber causes higher

796

A.M. Thornhill,

E.K. Glakpe and D. Agonafer

Vol. 15, No. {~

J

(a)

(b)

FIG. I0 Isoplots of turbulence kinetic energy: e = -0.652 and (a) Ra : 2xlO 5, (b) Ra : 106 lO

8

~r

6

4

2

0

lOO

200

FIG. ii Distribution of local equivalent conductivity on outer and inner cylinders lO-

10-

5

o

.....-'':

2

I~ ,o-

s

4

5

•

+ vl

Kc~i.cily

•

- vl

,cc,nt,~Jly

-:.: ..... "'--.... "\ 2

2

,

1

. . . .

,,

,

,

,,,,

,

i

J

,,,,

0

i

L

I

I

J

L

a

I

l

I

IE4

Rayleigh

Number

FIG. 12 Average Nusselt number distribution for various configurations

Angle (degrees)

FIG. iZ Local equivalent conductivity distributions on inner cylinder for various eccentricities and Ra = 106

Vol. 15, No. 6

CONVECTION IN ANNULAR ~%K2I/DSURES

797

levels of heat transfer to be obtained in all cases; the effect on the positive eccentric case being greatest.

Most of the positive eccentric annulus becomes

turbulent at high Rayleigh nt~nbers, while only parts of the negative eccentric case exhibit turbulent behavior at higher Rayleigh numbers.

In the future, the

modified code will be exercised in obtaining results for three-dimensional geometries of practical interest. Acknowledgement The support of the U.S. Deparhaent of Energy t~nder Grant Number DE-FGO5840H21495 is gratefully acknowledged. References i.

T.H. Kuehn and R.J. Goldstein, J. Heat Transfer, i00, 655 (1978)

2.

B. Farouk, Ph.D. Dissertation, Univ. of Delaware (1981)

3.

H. Ozoe, A. Mouri, M. OhTuro, S.W. Churchill and N. Lior, Int. J. Heat Mass Transfer, 28, 125 (1985)

4.

J.F. Thc~son, ~ r i c a l

5.

B.C.J. Chen, W.T. Sha, M.L. Doria and C. Schmitt, NRC Report # NUREG/CR1874 (1980)

6.

E.K. Glakpe, C.B. Watkins and J.N. Cannon, Nt~n. Heat Transfer, iO, 279 (1986)

7.

A.M. Thornhill, Master's Thesis, Howard University (1987)

Grid Generation, Elsevier Science Publ. Co. (1982)