Transient and steady bouyancy-driven flow in a horizontal cylinder

Transient and steady buoyancy-driven flow in a horizontal cylinder J. L. Xia,*t B. L. Smith,1 and G. Yadigaroglut -l-Nuclear Engineering IThermal Hydraulics

Laboratory, Laboratory,

Swiss Federal Institute of Technology, Paul Scherrer Institute, Switzerland

Zurich,

Switzerland

The ASTEC code is used to simulate both transient and steady two-dimensional, compressible, buoyancydriven laminarjow and heat transfer in a horizontal cylinder full of air. Calculations cover Rayleigh numbers from IO4 to 10’. The local and average Nusselt numbers, the velocity and temperature projiles along the vertical and horizontal diameters, and the velocity vector plots and temperature jields are presented for both the transient and steady-state cases. Results show that a near-stagnant, thermally stratiJied core exists only at moderately high Rayleigh numbers (> I@), the dimensionless temperature difSerence vertically across the core decreases with increasing Rayleigh number, and the time evolution of the flow and thermal fields is quite complicated. Numerical predictions are in good agreement with available experimental data.

Keywords: natural circulation, horizontal cylinder, numerical, transient

Introduction Buoyancy-driven natural convection in an enclosure has many practical applications that include, for example, cooling of nuclear reactors, thermal insulation problems, and others. Furthermore, solar energy collection, simulation of processes with computer codes provides suitable comparison of the performance of various numerical models. An increasing amount of research on this subject has been undertaken, but the configuration most studied is a rectangular enclosure. Though there have been benchmark solutions for this geometry, relatively sparse information is available for other basic geometries such as cylindrical. For the horizontal cylinder geometry, the first theoretical work was by Ostrach,’ who studied steady, laminar flow inside a horizontal cylinder with a cosine temperature distribution on the surface. He assumed that the core would be in isothermal, solid-body rotation. However, this assumption is not consistent with the experimental data of Martini and Churchill’ for a horizontal cylinder with the two halves maintained at different, uniform temperatures. Although this particular

Address reprint requests to Dr. Xia at the Nuclear Engineering Laboratory, Swiss Federal Institute of Technology, Clausiusstr.33, ETHZentrum, CH-8092 Zuerich, Switzerland. * Permanent address: Institute of Engineering Thermophysics, Chongqing University, Chongqing 630044, China. Received

10 June 1993; revised

‘0 1994 Butterworth-Heinemann

11 April 1994; accepted

8 June 1994

boundary condition is different from the cosine distribution of Ostrach, the similarity of the two cases is obvious. In fact, the experimental results showed that velocity and thermal boundary layers exist only within a narrow circulating band of fluid next to the cylinder surface, but the core was found to be thermally stratified and relatively stagnant. In later years, numerous theoretical analyses (based on Oseen linearization) and experiments have been carried out3-6; these are reviewed by Ostrach and Hantman6 and Ostrach.7 To the authors’ knowledge, there has been no systematic numerical analysis for natural convection in a horizontal cylinder, including transient behavior. The first numerical solution for this problem appears to be by Hellums and Churchill,8 which covers Ra up to 3.5 x lo6 and is the companion work to the Martini and Churchill experiments.2 Reasonably good agreement was obtained between numerical predictions and experimental data with guidance from the experiments. Leong and Davis’ numerically studied steady-state, natural convection in a half-heated, half-cooled horizontal cylinder and reported the results, mainly at Ra = 105. Other numerical calculations include those of Robillard et al.,” who examined the influence of time-dependent boundary conditions, and Takeuchi and Cheng,” who studied laminar convective heat transfer in horizontal cylinders with wall temperature decreasing at a constant rate. Detailed numerical simulations at higher Ra numbers and a computational benchmark solution for this basic geometry are not yet available. From an

Appl.

Math.

Modelling,

1994,

Vol. 18, December

691

Buoyancy-driven

flow in a horizontal

cylinder:

J. L. Xia et al.

engineering standpoint, it is also important to study transient processes. For example, the understanding of the system transient behavior during start-up and other time-dependent processes could be important. This paper reports the use of the ASTEC code” to simulate laminar, natural convection in a horizontal cylinder, and has the following objectives: (a) to supply details of the transient behavior; (b) to make a systematic examination of the influence of Ra on how and heat transfer at steady state; (c) to extend the range to Ra N 108; (d) to test the capability of the commercial code ASTEC in modelling natural convection by comparison with experimental data.

Governing equations and solution procedures The physical problem to be considered here is schematically depicted in Figure Z(a). The configuration studied consists of a horizontal cylinder with two vertical halves maintained at temperatures Th and T,, except for two narrow sectors of 7.5” at the top and bottom where adiabatic, or constant-heat-flux, conditions are applied. This arrangement is similar to the one studied by Hellums and Churchill8 and Leong and Davis,’ but with the narrow adiabatic regions at the top and bottom defined to avoid temperature discontinuities. The flow and heat transfer are assumed to be two-dimensional. Initially, the fluid is at rest and at temperature T,. The set of governing partial differential equations in tensor notation is

ap

apuj

(1)

z+c=O

XJ

aui BUi P~+P”jaxj=-~

ap

with R the gas constant for air. The boundary conditions are no-slip conditions on all the rigid walls, isothermal on the right and left side surfaces at their respective temperatures, and adiabatic for the narrow regions at the top and bottom of the cylinder. The solution of the governing equations (l)-(4), subject to the given boundary conditions, is carried out using the ASTEC code, which has been developed for modelling fluid dynamics and associated physical phenomena in complex geometries. This code uses an unstructured computational mesh, and the equations are solved using a finite-volume technique. For the present problem, we have structured the mesh as shown in Figure I(b). To obtain node-independent solutions, we conducted a series of trial calculations with 22, 33, 48, and 96 nonuniformly-spaced nodes along the diameter and found that the maximum deviations between 48 and 96 nodes are very small. For instance, at Ra = 3.55 x 106, the maximum deviation for the nondimensional temperature 4 is 0.57%, for velocity, l.O%, and for heat flux, 0.73%. Thus we chose 48 nonuniform nodes along the diameter (the minimum node spacing is 0.0014D) and 48 nodes along the azimuthal angle for all the ASTEC runs, as a compromise between numerical accuracy and computational efficiency. The transient behavior is different for different initial conditions. To examine the influence of Ra on the transient flow and heat transfer, we start calculating from zero-flow, constant-temperature (T = T,) conditions and set the same convergence criteria for all Ra values. Several convergence criteria are used. For example, the convergence limits are set to 10e6 J/kg for internal energy and to 10m6 m/s for velocity. The time step is adjusted to satisfy the above criteria (typical time step is 10m4 s). To validate the numerical solutions, comparison with available experimental data’ is made.

Results and discussion

a2T

aT

PC,ir+PC,uj~,=kar,+,+uj~

ap

J

(3)

+$(~+$)]

The density of air can be adequately perfect gas law,

represented

by the

p=pRT

(4)

Th

(b)

(4 Figure

1.

692

Ad.

Physical

Math.

model and refined finite-element

Modellina,

1994,

mesh.

Vol. 18, December

To characterize the time evolution of the flow, and the thermal behavior, a complete computation is carried out for Ra = 3.55 x 106. Figure 2 shows the velocity vector plots and the isotherm plots at various instants (the interval between isotherms is 0.1 (Th - T,)). Soon after initialization, the fluid in a thin layer adjacent to the hot surface is first driven upward; the fluid on the cold side and in the core is nearly motionless, Figure 2(a). That is, in the initial phase, the thermal and flow boundary layers are confined to a thin region adjacent to the hot surface. Later, the flow boundary layer gradually progresses from top to bottom along the cold side of the cylinder, and the thermal diffusion front evolves similarly (Figure 2(b) and 2(c)). At z z 0.025, a full flow ring is established adjacent to the cylinder surface (Figure 2(6)). In order to study the effect of Rayleigh number on transient flow and heat transfer, a complete solution is also computed for Ra = 108. Figure 3 shows the velocity vector plots and isotherms at various times for this case. The evolution is basically similar to that for Ra = 3.55 x 106, though at the higher Ra the flow is stronger,

flow in a horizontal cylinder: J. L. Xia et al.

Buoyancy-driven

(b)

7=

3.1 x 10-3

=2.8

(d) 7=2.5

(d) 7=3.4

x IO-*

Temporal evolution of velocity vectors and isotherms Figure 2. for Ra = 3.55 x 1 O6 (interval between isotherms is 0.1 ( Th- T,)).

and the thermal and velocity boundary layers are thinner. Almost at once, the flow field becomes noticeably unsteady, leading to a localized recirculation zone in the bottom region of the cylinder, as shown in Figure 3(u). At later times, Figures 3(b) and 3(c), the velocity vector plots show several asymmetric eddies growing in the core but, later still (Figure 3(d)), these dissipate, resulting in a nearly motionless core at steady state. Figures 4 and 5 give the velocity and temperature profiles along the vertical and horizontal diameters at different instants for both Ra = 3.55 x lo6 and Ra = 108. The change in the velocity with time is quite complicated and the behavior is more complex with increasing Ra. Before a full circulating band of fluid is formed near the cylinder surface, there exists an obvious flow in the core, though this is weak compared with the flow close to the cylinder surface. When steady state is reached, the velocities in the core are close to zero. At late times, close to steady state, the temperature and velocity profiles (solid curves in Figures 4(a), 4(c), 5(a), and 5(c)) show inverse symmetry.

Figure various

-4000

3. Velocity instants.

-.5

Figure horizontal

Appl.

’ : ’

’ : ’ : !

vector

’

1. ’ : 1

X0 4.

and

!

1 5

x 1O-3

x 10m3 isotherm

-4000

-.I

plots

’

!

’

I

for

’

’ . ’ : ’

Modelling,

1994,

! ’

Y

Temporal velocity profiles along the diameters for Ra = 3.55 x lo6 and lOa.

Math.

Ra = IO* at

vertical

Vol. 18, December

’

’ :

.?

and

693

Buoyancy-driven

0

-.-

3.06E-3

-

4.8X-2 Ra-3

flow in a horizontal

cylinder:

5X+6

Figure 7. History of the average cold surfaces.

-.-

--

X Figure 5. horizontal

J. L. Xia et al.

-

,.0,6E-3 ,.839E-3 4.676E-3 2.427E-2

Y

Temporal temperature diameters for Ra = 3.55

profiles along the vertical x 1 O6 and IO’.

and

Nusselt

number

on the hot and

multiple local extreme values, indicating more complex transient characteristics for the higher Ra value. Note also that the dimensionless time interval needed to reach the steady state, zS, becomes shorter with increasing Rayleigh number (7, z 0.048 for Ra = 3.55 x 106, and rS z 0.024 for Ra = log) as a consequence of the increased buoyancy effect. Figure 8 shows the velocity vector plots and the isotherm profiles at steady state for various Ra values. With Ra = lo4 (Figure 8(a)), no core is formed, the fluid circulates counterclockwise around the axis of the

The local Nusselt number along the cylinder surface at various times is presented in Figure 6. At early times, part of the cold surface does not contribute to the heat transfer because the front of the thermal diffusion has not yet reached that region. The maximum local Nusselt number occurs at the edge of the adiabatic section where the heating (or cooling) begins. This is because the largest temperature differences appear there. The local Nusselt number changes drastically with angle and a sharp change occurs near the top of the cold surface and near the bottom of the hot surface. The larger the Ra value, the more complicated the local Nusselt number profile becomes and the more dramatic the change of the local Nusselt number with time. The average Nusselt numbers of the hot and cold surface as a function of time are shown in Figure 7. At Ra = 3.55 x lo6 the average Nusselt number for the hot surface first decreases, then increases slightly, and finally decreases to a steady value; that of the cold surface increases monotonically with time. At Ra = 108, the situation is quite different. The variation of the average NM on the hot and cold surfaces with time now portrays

(4 !?a

-.-

3.08E-3 LZOE-2 2.07E-2 ,.wc:-2 Ra-3.55E.6

--

Figure 6. various

694

Local instants.

Awl.

-

Nusselt

Math.

-

number

Modelling,

along

the

1994,

cylinder

surface

at

Vol. 18, December

Figure 8.

Velocity

(c) Ra

vectors

and isotherms

at steady state

Buoyancy-driven cylinder; convection is rather weak and heat conduction has a clear influence. With increasing Ra, the flow becomes stronger. At Ra = lo’, two weak eddies are formed in the central region (Figure 8(b)). At still larger Ra values, fully developed thermal and velocity boundary layers are formed with the core relatively stagnant and thermally stratified. For Ra = 3.55 x 106, an eddy is formed in the upper hot surface region. The eddy moves toward the upper cold surface region at Ra = 10’. The temperature field in the core consists of a series of isothermal, nearly horizontal contours with the temperature increasing vertically (as shown in Figures 8(c) and 8(d)). These results are in good agreement with the experimental findings of Martini and Churchill’ for a heating-from-the-side configuration, as shown in Figure 9. They also qualitatively agree with those of Sabzevari and Ostrach’ for the cosine surface temperature distribution. With increasing Ra the velocity and thermal boundary layers become thinner. For example, the average boundary layer thickness is about 0.170 for Ra = 3.55 x lo6 and O.lD for Ra = 10’. Moreover, the dimensionless temperature difference vertically across the core decreases with increasing Ra value. With Ra = 3.55 x 106, the dimensionless temperature drop in the core is about 0.5, but this reduces to about 0.2 at Ra = lo*. The local Nusselt numbers along the cylinder surface at steady state are plotted in Figure 10 for various Ra values. Inversely symmetrical profiles exist with the air flow upward on the hot wall and downward on the cold wall. The average Nusselt numbers at various Ra values are given below; with increasing Ra, the Nusselt number increases due to the stronger convective flow. Ra Nu

lo4

lo5

2.03

3.79

Comparison with experimental

3.55 x 106 10.4

108 49.3

(4

cylinder:

J. L. Xia et al.

at various

Rayleigh

J

Local Figure 10. at steady state.

Nusselt

numbers

numbers

profiles along the vertical diameter (Figure Zl(b)), and the local Nusselt number along the cylinder surface (Figure II(c)), the numerical calculation predicts higher values in the top and bottom regions of the cylinder. This is explained by the differences in the boundary condition in these areas. In a further calculation, we have assumed a heat loss (59.5 W/m”) at the top, and a heat gain (53.0 W/m”) at the bottom for Ra = 3.55 x lo6 (the values are taken from Martini and Churchil12) to more closely model the actual conditions of the experiments. The results are also plotted in Figure 11 (dashed lines). It is seen that the predicted temperature profiles along the horizontal and vertical diameters now virtually coincide with the experimental data. For the Nusselt number profiles around the cylinder surface there remain some discrepancies. These may also be caused by the uncertainty in determining the heat transfer coefficient from the temperature gradient in the experiments. In fact, the thermocouple placement in the experiment is probably not able to capture the sharp change of local Nusselt numbers near the top of the hot surface, and

data

Martini and Churchill2 carried out experiments on natural convection inside horizontal cylinders. The boundary conditions are similar to those used in the present study, except that some heat transfer occurred at the top and bottom surface regions of the cylinder in their experiments, whereas these are assumed to be adiabatic in the present case. Comparison of the present results with the experimental data for Ra = 3.55 x lo6 is given in Figures 9 and 11. For the temperature profiles along the horizontal diameter, excellent agreement is achieved, as shown in Figure I I (a). For the temperature

Figure 9. Ra = 3.55

flow in a horizontal

lb)

Comparison of predicted isotherms with experiment* x 106. (a) experimental; (b) predicted.

at

---

with heat losses

~I!l:I!l:pI!l.lL X

(a)

(b)

Figure 11. Comparison of numerical experimental data*. (a) horizontal diameter; (c) cylinder surface.

Appt. Math. Modelling,

.5

Y

prediction with the (b) vertical diameter;

1994, Vol. 18, December

695

Buoyancy-driven

flow in a horizontal cylinder:

J. L. Xia et al. UiD

near the bottom of the cold surface, as predicted by the numerical simulations. A closer comparison is therefore impossible.

dimensionless velocity component, ~ % Cartesian coordinates dimensionless Cartesian coordinates, 2

Conclusions At steady-state, a near-stagnant flow region is not always formed in the horizontal cylinder. At low Ra values, the entire fluid circulates around the cylinder axis, whereas at moderately high Ra values, the fluid circulates rapidly in a thin layer adjacent to the cylinder surface and a relatively stagnant, thermally stratified core is formed. With increasing Rayleigh number, the flow and thermal boundary layer thicknesses decrease (resulting in an increase in the Nusselt number) and the dimensionless temperature difference vertically across the core also decreases. Transient characteristics in the horizontal cylinder are quite complex. At higher Rayleigh numbers there appear to be somewhat periodic, unsteady flow changes. Single or multiple eddies may occur in the core during the transient. With low Rayleigh numbers, the Nusselt number first decreases, then increases and finally decreases to a steady-state value. With high Rayleigh numbers, several peaks occur. Comparison with the available experimental data was excellent and has validated the ASTEC predictions.

thermal diffusivity thermal expansion coefficient T-T,

dimensionless temperature, ~ K---T, azimuthal angle measured counterclockwise the top dynamic viscosity density dimensionless time, s Subscripts

cold wall hot wall reference state at temperature

: 0

c, D

thermal conductivity specific heat at constant pressure diameter of the cylinder

Nu

local Nusselt number, -

7% P r

average Nusselt number pressure radius

R

dimensionless radius, i

Ra

Rayleigh number, “(’

5

w

- ‘)03 vll %

t T

time temperature

To

average temperature,

%

velocity component

696

Appl.

Math.

ml+ T,)

____ 2

Modelling,

1994,

Vol. 18, December

TO

References

Nomenclature k

from

Ostrach, S. A boundary layer problem in the theory of free convection. Ph.D. Thesis, Brown University, Providence, RI, 1950 Martini, W. R. and Churchill, S. W. Natural convection inside a horizontal cylinder. AIChE J. 1960, 6, 251 Weinbaum, S. Natural convection in a horizontal cylinder. J. Fluid Mech. 1964, 18, 409 Ostrach, S. Completely confined natural convection, development in mechanics. Proceedings of the 10th Midwestern Mechs. Conference, 1968, 4, 53 Sabzevari, A. and Ostrach, S. Experimental studies of natural convection in a horizontal cylinder. Proceedings of the 5th International Heat Transfer Conference, 1974, 3, 100 Ostrach, S. and Hantman, R. G. Natural convection inside a horizontal cylinder. Chem. Eng. Comm. 1981, 9, 213 Ostrach, S. Natural convection in enclosures. J. Heat Transfer 1988, 110, 1175 Hellurns, J. D. and Churchill, S. W. Transient and steady state free and natural convection, numerical solutions, Part II-the region inside a horizontal cylinder. AIChE J. 1962, 8, 692 Leong, S. S. and de Vahl Davis, G. Natural Convection in a Horizontal Cylinder, Numerical Methods in Thermal Problems, eds. R. W. Lewis and K. Morgan, Pineridge Press, 1979, p. 287 Robillard, L., Vasseur, P., and Nguyen, H. T. Thermal stratification induced by a periodic time-dependent temperature imposed on the boundary. J. Heat Transfer 1987, 109, 525 Takeuchi, M. and Cheng, K. Transient natural convection in horizontal cylinders with constant cooling rate. Waerme und Stoffiebertragung 1976, 9, 215 Lonsdale, R. D. ASTEC release 3.0, User Manual, 1991