NASA Contractor Report 178027

leASE REPORT NO. 86-6

ICASE

NASA-CR-178027

19860014382

SIMULATION OF THE TURBULENT RAYLEIGH-BENARD PROBLEM

USING A SPECTRAL/FINITE DIFFERENCE TECHNIQUE

T. M. Eidson

M. Y. Hussaini

T. A. Zang

Contract Nos. NASl-17070 and NASl-18l07

February 1986

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING

NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

NI\SI\

National Aeronautics and

Space Administration

Langley Research Center

Hampton,Virginia 23665

AD P n r; h I ~J6 {""".f 1\ ~LI I.:JU

JJ.r~SLtY I,ES'>'li ;r::NTER

L1SR,l..RY. i\IASA

/~."~.:~TCN. VIRGINIA

SIMULATION OF THE TURBULENT RAYLEIGH-BENARD PROBLEM

USING A SPECTRAL/FINITE DIFFERENCE TECHNIQUE

T. M. Eidson

Georgia Institute of Technology

M. Y. Hussaini

Institute for Computer Applications in Science and Engineering

T. A. Zang

NASA Langley Research Center

ABSTRACT

The three-dlmensional, incompressible Navler-Stokes and energy equations

with the Bousslnesq assumption have been directly simulated at a Rayleigh

number of 3.8 x 105 and a Prandtl number of 0.76. In the vertical direction,

wall boundaries were used and in the horizontal, periodic boundary conditions

were applied. A spectral/finite difference numerical method was used to

simulate the flow. At these conditions the flow is turbulent, and a

sufficiently fine mesh was used to capture all relevant flow scales. The

results of the simulation are compared to experimental data to justify the

conclusion that the small scale motion was adequately resolved.

Research was supported by the National Aeronautics and Space

Administration under NASA Contract Nos. NASI-17070 and NASI-18107 while the

second author was in residence at ICASE, NASA Langley Research Center,

Hampton, VA 23665-5225.

I° INTRODUCTION

Direct simulation of turbulent fluid flows is now possible with the large

vector computers that have become available [1,2]. Prediction of low-order

flow statistics is definitely within current capabilities, and some results

have already been published which show predictions of small scale turbulent

features which are consistent with experimental observations [3-5]. The

current study was undertaken to explore the quality of information that can be

extracted from a direct flow simulation (DFS) of turbulence on a sufficiently

fine mesh.

The turbulent Rayleigh-Benard problem (natural convection)was chosen for

study since it is a simple turbulent flow for which a good body of

experimental measurements exists. Moreover, some DFS and large-eddy

simulations (LES) of this problem have been published albeit on coarser

meshes. While experimental data do exist, measurements of velocity, where no

mean flow exists, are difficult. Hence, there is much to be learned about

turbulent natural convection from an accurate simulation.

The two requirements for conducting such a study are a hlgh-speed computer

and an efficient, accurate flow simulation code. The CYBER-205 computer with

a 16 mega-word memory provides sufficient computation power. This current

code has been extensively tested, and various versions of it have been used to

study transition in channel flow [6]. The version used in this study includes

the addition of the energy equation and a modified vertical momentum equation

that includes buoyancy consistent with the Boussinesq assumption.

A simulation of a turbulent flow was then conducted, and these data as

well as a discussion of the code will be presented in this paper. Though the

overall goal of this work is an in-depth examination of the quality and type
of information that can be extracted from such a simulation, the purpose of

this paper is to document the basic simulation. The simulation results will

be compared with experimental mean measurements as well as previous DFS

results. The increased resolution of this work over previous DFS resulted in

an improvement in the prediction of the Nusselt number; it was sufficiently

close to experimental results to suggest that in addition to a good prediction

of the large-scale flow, the small-scale features are accurately represented.

Grotzbach discusses this connection extensively [7,8]. Comparisons with

experimental data which are more dependent on the small scale components of

the flow will also be presented to justify further this conclusion.

II. RAYLEIGH-BENARD PROBLEM

The Rayleigh-Benard problem is a simple geometry, laboratory-type problem

used to study natural convection (Figure I). Chandrasekhar [9] and Busse [I0]

have described the basic problem and discuss both the stability analysis and

some experimental results. Krishnamurti [II,12] has summarized much

experimental data and developed a map showing the qualitative flow at

different values of Ra and Pr, the principal independent problem parameters

(defined below). For Pr = 0.76 (air) and Ra = 3.8 × 105 the motion is

turbulent, although it should possibly be qualified as low Reynolds number

turbulence. Several experimental studies [13-18] and numerical simulations

[7,19,20] (both LES and DFS) have been completed in the qualitatively similar

Pr-Ra region. The flow at these values of Pr and Ra consists of a core

flow (a horizontal layer in the middle 80% of the fluid layer) and a boundary

region near each plate. The turbulence is statistically homogeneous in the
horizontal directions for both layers. In the core the vertical variation of

most statistical quantities is small. In the boundary layer there is a

transition from molecular dominated physical processes near the wall to the

fully turbulent core flow.

This flow is described by the incompressible Navier-Stokes equations

modified to include the effect of temperature-induced density variations on

the buoyancy force (Boussinesq assumption) plus the temperature equation.

These equations, when non-dimensionalized by _, h, and AT, are

aui a(uiu j ) aP a2ui

--+ = ---+ Pr----_+ Pr Ra T_i3 , (la)

at ax.3 axi ax.3

aT a(Tuj) a2T

--+ - 2 + u3 ' (Ib)

at ax. ax.

3 J

and

au.

.__i= O. (lc)

ax.

3

The temperature and pressure in Equations (la,b,c) are the difference between

the actual temperature, Ta, and pressure, Pa' and the values due to the static

temperature gradient only. These are defined as follows:

Ta (x,t) _ = TO - x3 + T(x,t)

= T + T (x,t)

o r --

and
4

_Pa(_ 't) = _ gh3 _P(!, t)

_x3 ---_-c Pr Ra x3+ _x3 .

The dependent variables in the problem are the velocity components, ui or

(u,v,w), the temperature, T, and the fluid pressure, P. The independent

variables are the spacial coordinates, xi or __x= (x,y,z), and time, t. The

indices i = 1 and i = 2 signify the horizontal directions, and i = 3

denotes the vertical direction. The problem parameters are the thermal

dlffuslvlty, c; the kinematic viscosity, _; the acceleration of gravity,

g; the reference fluid density, Po; the coefficient of thermal expansion, B;

a reference temperature (the temperature of the lower plate), To; and the

temperature difference between the two plates, AT. The dependent and

independent variables have been non-dimensionalized by c/h (velocity), AT

(temperature), Po c2/h2 (pressure), h (coordinates), and h2/c (time). The

Raylelgh number Ra - gBATh3 and Prandlt number, Pr = v/c, are the

principal non-dlmensional problem parameters.

III. NUMERICAL METHOD

The Fourier finite difference algorithm developed by Moin and Kim [21]

for their large-eddy simulations of turbulent channel flow has been applied,

with modifications, to the present direct simulation of turbulent RaylelghBenard flow. This Raylelgh-Benard Fourier-finite difference Method (RBFFDM)

is an unspllt method on a grid staggered (for the pressure variable and the

continuity equation) in the vertical (z) direction. Fourier collocation is

used for the spaclal discretizatlon in the x and y directions whereas in

the vertical direction second-order flnlte-differences are employed on the
5

non-uniform grid,

zk = (I - cos(k_/Nz))/2 , k = 0,1,-.-,N z.

The time discretization is Crank-Nicolson for the viscous and conductive

terms and backward Euler for the pressure gradient term. The advection and

buoyancy terms are handled by a third-order Adams-Bashforth method.

The implicit part of the algorithm requires, for each pair of horizontal

Fourier wavenumbers, the solution of 2 real, block tridiagonal systems

(involving the velocities and pressure) and, independently, 2 real, scalar

tridiagonal systems (for the temperature). The block-tridiagonal equations

were scaled as described by Zang & Hussaini [6] for their Fourler-Chebyshev

version of the corresponding channel flow algorithm. Pivoting has proven to

be unnecessary for this system. The block-trldiagonal solution algorithm

takes advantage of the many zero elements which occur. Vectorization of this

phase of the algorithm is achieved by solving for many pairs of Fourier

wavenumbers at the same time.

Equations (la,b,c) were solved on the region, 0 < x < Ax, 0 _ y < Ay,

0 _ z < I. Under the present scaling of the vertical direction, the lengths

A x and Ay correspond to the aspect ratios of the two horizontal directions

to the vertical one. The boundary conditions at the lower and upper walls,

z = 0 and z = I, are the conventional no-slip and no temperature jump

conditions. In the horizontal directions, periodic boundary conditions are

assumed. These aspect ratios must be large for reasonable correspondence with

experiments. The aspect ratios are related to the resolution by A = N Ax x X

and A = N Ay. The computer memory limitation places an upper bound on

Y Y
6

N = Nx Ny Nz; therefore, to have a large Axl along with sufficiently small

Axi, some compromise is necessary. The values chosen, Ax = 4 and _ = 2,

allow the available computer memory to be used for better small scale

resolution. This will be further discussed in Section IV.

The RBFFDM code has been implemented on a CDC Cyber 205 with 2 pipes and

16 million 64-blt words of main memory. For each grid point 13 variables were

stored. Additional storage equivalent to 7 variables per grid point was used

to facilitate vectorlzation. A total of II million words was used for the

simulation on the 128 x 64 x 64 grid. Vector lengths for the explicit portion

of the algorithm were between 4 and 40 thousand. Typically, one-fourth of the

implicit equations were solved together. The vector lengths here were roughly

1,000. The linked triad feature was heavily exploited. A single tlme-step

required 6.8 seconds of CPU time and no I/0 time since the job was run

entirely within the central memory. The sustained speed of the calculation

was I00 MFLOPS.

The start-up phase of the calculation took 2100 time-steps and the data

collection an additional 5600 steps. A total of 12 hours of CPU time was

required for the data collection. This includes the time for some preliminary

diagnostics. The Courant number, defined as the maximum over the grid of the

quantity

ranged between 0.19 and 0.26 and averaged 0.23 for the data collection phase.
IV. COMPARISON WITH PREVIOUS WORK

The turbulent flow of the Rayleigh-Benard problem is assumed to be

homogeneous in both horizontal directions as well as statistically steady in

time (after a start-up period). Experimental data are usually presented using

some combination of a long-time average as well as a spacial average in one or

both horizontal directions. The simulation results presented as horizontal

averages, < >, have been averaged over both the x and y directions (except

for the I-D x-spectra which were averaged in y only). In addition, they

have been time averaged over a time period equal to 10/Wc, where

W c = (Nu Pr Ra) I/3. (See Deardorff [22] for a discussion of Wc, the scaling

velocity for the large eddies.) This period, which should consist of several

large eddy turn-over times, was found adequate by Eidson [19]. Volume

averages were tlme averaged as well. All the simulation results presented

below are both horizontally and temporally averaged unless otherwise

specified.

Prediction of the average vertical heat flux, the Nusselt number in nondimensional form, is an important result of any natural convection study.

Previous simulations have predicted values of the Nusselt number, Nu, which

are slightly higher than those measured experimentally. Grotzbach [8] has

discussed extensively this discrepancy and has shown that inadequate

resolution is partially responsible. In Table I, the results from both

simulations and experiments are shown. The prediction of the present study

lies at the the upper range of experimental measurements and below the Nusselt

numbers predicted by previous, coarser-grid simulations. In the present

simulation, Nu was calculated at each z level. The average value is

reported in Table I. The variation with z was small (approximately ± 0.I)

except very near the lower wall where Nu increased to 7.0.
The aspect ratio of the horizontal to vertical boundary lengths also is

known to affect Nu [8,15]. The values of Ax = 4 and Ay = 2 are smaller

than the values of 4 to 7 suggested by experimentalists as the minimum for

removing significant side boundary effects. Although the aspect ratio effect

is not negligible, especially on Nusselt number [8], a large aspect ratio was

foregone in favor of better small scale resolution in view of the goal of this

study: to resolve eddies down to nearly the dissipation scales.

The relative temperature, Tr, vertical profile is compared to data of

_Tr

studies by Chu and Goldstein [16] in Figure 2. A line with slope, B--_- = Nu,

is drawn in this figure. From this one can estimate the conductive layer

thickness, _ . This will be used later in a more extensive examination of C

the temperature data near the wall.

In Figures 3 and 4, the vertical dependence of the velocity and

temperature RMS values are compared to the simulation of Grotzbaeh [7] and the

measurements of Deardorff and Willis [13] (slightly larger Ra). In Table 1 a

comparison for the centerline values is presented for a wide range of studies.

The present improved resolution reduces slightly the predicted levels of the

RMS temperature over previous simulations. However, no systematic change in

the RMS velocity levels with the improved spacial resolution was observed.

Both trends are consistent with Grotzbach's results. Considering the

variation in experimental values and the uncertainty in aspect ratio effects,

the results appear quite satisfactory.

The wT correlation coefficient, C(wT), is constant for 0.2 < z < 0.8,

giving a value of 0.71 (Figure 5). Both Grotzbach [7] and Eidson [20]

previously obtained a value of approximately 0.67. Deardorff and Willis

measured approximately 0.60 for a slightly higher Ra, but C(wT) should
decrease with increasing Ra. Near the wall all four studies differ.

Deardorff and Willis measured a significant drop as the wall was approached,

but warned that since the numerator and denominator of C(wT) become small,

their results are uncertain. The increase in C(wT) near the wall found in

the current study was not observed in either previous simulation, but neither

of these had sufficient resolution in this region. Note that only the data of

Grotzbach's case 9 with 16 vertical grid points was available for

comparison. Also, horizontal averages of the uv, uw and vw correlations

were calculated. They were all approximately zero as would be expected for

turbulence homogeneous in the horizontal directions.

In Figure 6, the several terms in the horizontally averaged kinetic

energy equation are plotted versus z and compared with experimental data at

a slightly larger value of Ra. These terms are

a<E> > _ a

_-E - Pr Ra <wT a _-_ <w(E + Pa )>

Production Diffusion

aui aui a2<E>

- Pr <_-_]_-__> + Pr

.] .] _ '

Dissipation Molecular Transfer

where

I

E = _ ui ui .

Since the production term in the core is equal to (Pr Ra Nu), the

experimental data of Deardorff & Willis [13] at Ra = 6.3 × 105 and Nu = 5.7

would be expected to be about 40% higher than the current simulation. For
10

clarity, only a few key values from Deardorff and Willis are shown to

demonstrate the general agreement between experiment and simulation. The four

terms are all essentially constant and approximately 40% below the

experimental data for the core region. In the boundary region the variation

with z is similar for both experiment and simulation. The diffusion term

(for the experiment and simulation) peaks at a z value of 0.035, and the

molecular diffusion term of the simulation peaks at a slightly higher z. The

experimental point plotted at a low z value is at these peaks for these two

terms. Actually, one would expect the peaks for the lower Ra to be at a

slightly larger z since the boundary region thins with increasing Ra. From

a closer inspection of the data, the molecular dissipation also can be seen to

have a small positive value between z = 0.1 and z = 0.2 as was found by

Deardorff and Willis. The production and dissipation terms have the same

general shape in both studies. The volume averages of the production and

dissipation terms were 1.641 x 106 and 1.637 × 106 , respectively. Also, the

volume average of the molecular transfer term was 9.0 x 103 .

One advantage of the simulation technique is that quantities which are

experimentally difficult to measure can be easily calculated. Figures 7 to 9

give three examples. In Figure 7 the absolute value of each of the three

components of vorticity is shown. As expected, the x and y components are

large near the wall due to the creation there of a boundary layer by the large

eddies. Near the center the flow is more isotropic. The ratio of the volume

average of the horizontal to vertical vorticity was 5.0 and 4.3 for the x

and y directions, respectively. This is higher by a factor of 2 than in the

previous DFS by Eidson [19]. Figure 8 shows that the only significant

velocity skewness is for the w component near the wall. A negative value
ii

for the w component is reasonable since the fluid particles with negative

velocity (near the lower wall) come from the core where there is more velocity

variation due to the turbulent cascading process. The particles with a

positive velocity originate near the wall where the motion is damped and they

have a more uniform value. Recall that the horizontal average of the w

velocity equals zero. The flatness profile of the velocities is shown in

Figure 9. For the horizontal components, these approach 3 in the core, a

value which is similar to that in channel flow turbulence away from the

wall [21]. The larger flatness factor near the wall suggests that the flow is

more intermittent in this region.

The wT cross-spectra for several z values are plotted in Figures i0

to 13. These are I-D spectra in the x-direction. They have been normalized

so that they integrate to <wT>. The wavenumber, <, is defined such that 2_x

is the argument of the periodic expansion functions. Notice that the spectra

decrease with < except near the wall where they are roughly constant with

<. This trend was found also in the u, w and T spectra but not in v.

Near the center the most energy would be expected in the largest available

scales for the small values of A in this study. Near the wall, smaller

horizontal structures would be expected, and thus significant energy at a

wavenumber greater than the minimum is a reasonable result. Since this is low

Ra turbulence, an equilibrium region would not be expected, and indeed none

was found here. The drop-off with wavenumber is rapid due to the importance

of viscous effects at all the scales.

To evaluate the simulation further, a more sensitive measure of the

small scale features of the flow field, especially in the wall region, was

sought. Carroll [14] has made direct instantaneous measurements of the
12

temperature and the vertical temperature gradient and calculated the mean,

standard deviation, and skewness. Carroll's data are compared with the

results of this simulation in Figures 14 to 20. The data of Carroll in these

figures are a mean curve (drawn by Carroll) through the data points. The bars

show the range of data scatter. The mean of the relative temperature

(simulation data shown in Figure 2) and its vertical derivative (Figure 18)

show good agreement with Carroll's data. The higher order statistics give a

variation with z which is similar to the experimental results although the

magnitudes differ in some cases. In Table I, a significant variation in

experimental measurements of the RMS of Tr, R(Tr), between different

experimental studies is apparent. Carroll used a stationary probe, and this

may have reduced the fluctuations coming from the direction of the probe and

lowered the statistical levels. Also Carroll measured the same RMS

temperature level for Nu ranging from 5.5 to 14.0. Both experiments [13,15]

and simulations [20] have shown that R(Tr) increases with decreasing Nu

(or Ra). Since the magnitude of the data is in question in some cases, the

conclusions drawn from the simulation/experiment comparison are based mainly

on changes of slope with vertical distance from the wall. These variations

correspond to differences in the dominant physics in the different layers, and

the simulation was able to predict these layers in agreement with the

experimental results.

A four-layer (or three-layer with the middle divided into 2 sub-layers)

model is consistent with both the experimental and simulation'results. These

are :
13

(i) a conduction dominated layer 0 < Nu*z < 61

0.i < 6" < 0.3

1

(ii) an inner boundary layer 61 < Numz < 62

6 2 _ 0.5

(iii) an outer boundary layer 62 < Nu*z < 63

63 _ 1.5

(iv) the core 63 < Nu*z < Nu/2 •

These layers are all well resolved by the simulationas they contain 8, 4, 9,

and ii grid points respectively (assuming 61 = 0.2).

Determination of 61 is uncertain. As was mentioned previously (Figure

2), divergence of the temperature profile from a straight line with a

normalized slope equal to Nu (the temperature gradient at the wall) gives

61 = _c = 0.26. Another measure is that the wT correlation rises from I%

to 10% of its core value near Nu*z _ 0.i.

The inner boundary layer is characterized by the transition from

molecular to turbulent dominated processes. Both experiment and simulation

8Tr

show that 8T has been reduced by 50% between Nu*z of 0.4 and 0.5 (Figure

18). In the inner boundary layer, R(Tr) reaches a maximum (Figure 14), the

skewness factor for Tr, S(Tr) , changes sign (Figure 16), and the flatness

factor for Tr, F(Tr) , reaches a minimum (Figure 17). The skewness, S(Tr) ,

continues to increase reaching a positive maximum in the outer boundary layer.

The above trends, shown in both studies, suggest significant changes in the

flow character in these layers.
14

In the outer layer, the flow transitions to a core region where the

large eddies (the most efficient heat transfer agents) carry the heated fluid

_Tr

from the bottom to the top. At Nu*z = 1.5, the temperature gradient, _ ,

has been reduced to I% of its maximum value (Figure 18). The skewness, S(Tr) ,

decreases in the outer layer and approaches zero near the center of the core

(Figure 16). Also note, in the outer layer R(T r) varies as z-I/3 compared

with z+I in the inner layer, a trend which is found in both studies (Figure

15). The simulation predicts a transition to a constant level for R(T r) in

the center in agreement with Deardorff and Willis [13] (Figure 4). However,

Carroll's data suggest that R(T r) varies as z-I/3 in the core as well as

the outer layer.

A possible discrepancy between the simulation results and Carroll's data

_Tr

is seen in the R(z_i ) plot (Figure 19). At values of Nu*z between 0.4

and 1.0 the simulation shows a hump. However, Carroll's choice of an "average

curve" has significant scatter only on top of the curve. Re-drawing the data

fit through this scatter would result in a similar hump in the outer boundary

layer region. Another example, where the simulation results suggest that the

"averaged curve" through sparse, scattered data is not correct, is Carroll's

choice for the temperature skewness profile for Nu*z < 0.i (Figure 16). He

suggests that the curve should continue decreasing to S(Tr) = -1.5 and then

rise to zero at Nu*z = 0.01. However, only a few data points which gave

significant scatter were measured in this region. From an examination of the

data in Carroll's paper, it is not inconsistent to draw the experimental data

fit to turn between Nu*z = 0.1 and 0.2 and be parallel to the simulation

results somewhere between the dashed lines shown in Figure 16. One

discrepancy that cannot be explained by data scatter is the comparison between
15

8T r 8Tr

studies of S(_--_-) away from the wall. The simulation predicts that S(_--_--)

returns to zero in the outer layer and remains zero in the core (Figure 20).

However, the experiments show a return to negative skewness in the core.

The final comparison is for the temperature gradient versus Nu*z near

the wall (Figure 18). Various power laws have been hypothesized which predict

a linear slope for the data plotted in log-log form. Carroll [14], Businger

[23], and Monin and Yaglom [24] discuss the various theories. Carroll's data

show that only for the range 0.4 < Nu*z < 1.0 does a simple power law seem

reasonable and even then the data have a slow change in slope in this

region. The current simulation fits Carroll's data nicely and thus also

predicts a slope of -2 in this range of z. The value of -2 is predicted by

the theory of Malkus [25].

From these comparisons, it can be seen that the simulation data

represent a turbulent realization which has a good agreement with

experimental, averaged measurements. The agreement, at least in prediction of

sub-layers using the variance and skewness data, shows that the small scale

turbulent features are accurately represented near the wall. This is

important because a complete understanding of the Rayleigh-Benard convection

requires studying the thin boundary layer which exists as a result of the

large scale motion in the core. From isothermal, flat plate boundary layer

studies, it is known that significant small scale events occur in this type of

region which are important to the global flow.

Some of the features that we have observed in this calculation were also

found in a reference pointed out to us by a referee [26]. However, we have

not yet performed a detailed comparison with this data.
16

V° CONCLUDING REMARKS

From the comparison with the gross properties of the experimental data,

it can be concluded that the completed computer simulation is a reasonably

accurate realization of a turbulent flow even down to the dissipation scales.

Moreover, the detailed database has enabled us to identify a feature in the

skewness of the vertical temperature derivative which has been overlooked in

the experiment due to data scatter.

Obviously, a far more detailed analysis of the simulation results is in

order. This simulation is of sufficient accuracy to study the small scale

structure of a natural convection flow in a similar manner to the work

previously done for channel flows [4,21]. Moreover, it can be used to help

interpret experimental measurements which are typically limited to a few

spacial points and fewer than all three velocity components. Finally, we feel

this database can be used to evaluate turbulent models and other theoretical

ideas.

Acknowledgment

The authors would llke to express their appreciation to Control Data

Corporation for the use of a CYBER-205 at the Arden Hills Plant in Minnesota

during the initial phase of this study, and to J. Lamblotte and G. Tennille

for their assistance in running the code.
17

References

I. Rogallo, R. S. and Moin, P. [1984] Numerical simulation of turbulent

flows. Ann. Rev. FI. Mech., 16, 99-137.

2. Schumann, U., Grotzbach, G. and Klelser, L. [1980] Direct numerical

simulation of turbulence. Prediction Methods for Turbulent Flows,

Hemisphere Publ. Corp., (W. Kollman, ed.), 124.

3. Brachet, M. E., et al., [1983] Small-scale structure of the TaylorGreen vortex. J. Fluld Mech., 130, 411-453.

4. Moin, P. and Kim, J. [1985] The structure of vorticlty field in

turbulent channel flow. Part 2. Study of ensemble-averaged fields. J.

Fluld Mech. To appear.

5. Moser, R. D. and Moin, P. [1984] Direct numerical simulation of curved

turbulent channel flow. NASA TM 85974. Also Report TF-20, Department

of Mechanical Engineering, Stanford University.

6. Zang, T. A., and Hussaini, M. Y. [1985] Numerical experiments on

subcritical transition mechanisms. AIAA Paper 85-0296.

7. Grotzbach, G. [1982] Direct numerical simulation of laminar and

turbulent Benard convection. J. Fluid Mech. I19, 27.
18

8. Grotzbach, G. [1983] Spaclal resolution for direct numerical simulation

of the Rayleigh-Benard convection. J. Comp. Phys. 49, 241.

9. Chandrasekhar, S. [1961] Hydrodynamic and Hydromagnetlc Stability.

Clarendon Press, Oxford.

i0. Busse, F. H. 1981 Transition to turbulence in Rayleigh-Benard

convection. Topics in Applied Physics. Vol. 45: Hydrodynamic

Instabilities and the Transition to Turbulence. (H. L. Swinney &

J. P. Gollub, eds.), Springer-Verlag.

II. Krishnamurti, R. [1970a] On the transition to turbulent convection.

Part i. The transition from two- to three-dimensional flow. J. Fluid

Mech. 42, 295.

12. Krishnamurti, R. [1970b] On the transition to turbulent convection.

Part 2. The transition to time-dependent flow. J. Fluid Mech. 42, 309.

13. Deardorff, J. W., and G. E. Willis [1967] Investigation of turbulent

thermal convection between horizontal plates. J. FluldMech. 28, 675.

14. Carroll, J. J. [1976] The thermal structure of turbulent convection.

J° Atmos. Sc. 33, 642.

15. Fitzjarrald, D. E. [1976] An experimental study of turbulent

convection in air. J. Fluid Mech. 73, 693.
19

16. Chu, T. Y. and R. J. Goldstein [1973] Turbulent convection in a

horizontal layer of water. J. FluldMech. 60, I, 141.

17. Goldstein, R. J. and T. Y. Chu [1969] Thermal convection in a

horizontal layer of air. Prog. Heat and Mass Transfer. 2, 55.

18. Denton, R. A. and I. R. Wood [1979] Turbulent convection between two

horizontal plates. Int. J. Heat Mass Transfer. 22, 1339.

19. Eidson, T.M. [1982] Numerical simulation of the turbulent RayleighBenard problem using subgrid modeling. Ph.D. Thesis, University of

Michigan.

20. Eidson, T.M. [1985] Numerical simulation of the turbulent RayleighBenard problem using subgrid modeling. J. Fluld Mech. 158, 245.

21. Moin, P. and Kim, J. [1982] Numerical investigation of turbulent channel

flow. J. Fluid Mech. 118, 341.

22. Deardorff, J. W. [1970] Convection velocity and temperature scales for

the unstable planetary boundary layer and for Rayleigh convection. J.

Atmos. Scl. 27, 1211.

23. Businger, J.A. [1973] Turbulent transfer in the atmospheric surface

layer. Workshop in Micrometeorology, American Meteorology Society, 67.
20

[24] Monin, A. S. and Yaglom, A. M. [1971] Statistical Fluid Mechanics, 1

MIT Press.

[25] Malkus, W. V. R. [1954] The heat transport of thermal turbulence.

Proc. Roy. Soc. (London), A. 225 196.

Grotzbach, G. [1980] Uber das raumllche Auflosungsvermogen numerlscher [26]

Simulatlonen von turbulenter Benard-Konvektlon. Kernforschungszentrum

Karlsruhe. 2981B
21

Table I. Horizontal Averages at the Vertical Center

Direct Simulation Grid Ra Nu URMS WRMS TrRMS

Current Study 128×64x64 3.8xi05 6.6 52 116 0.08

Grotzbach* 3.8xi05

(case 7) 16x16x16 7.8 -- 109 --

(case 9) 32×16x16 7.4 50 107 0.08

(case 14) 64x32x32 6.9 -- 109 0.08

Eidson 64x64x16 3.8xi05 8.1 70 112 0.II

Experiment

Deardorff & Willis 6.8xi05 5.8 65 120 0.08

Carroll 3.8×105 5.9* .... 0.05*

Fitzjarrold 3.8xi05 6.2* 80* 120" 0.05*

Chu & Goldsteln 3.8xi05 6.5* ......

Goldstein & Chu 3.8xi05 5.4. ......

*Calculated from curve fit of data over a range of Ra.
22

s_-To-AT

_I//I / II IIIIIII I/I//II // /_._/// / /i/ //i I/// I ii

I

g

z Fluid: vl _1 13_Pot Cp h

_//I/// IIII/II/II//II/II/'I /II ////

L. _ To_/

Figure 1. Geometry and Parameters of the Rayleigh-Benard Problem
23

o

o

o

o Chu & Goldstein. Ref. 16

O=1/2·Nu

0c~ 0.04

Nu *Oc~ 0.26

o

Oc

N

oI Current Simulation

o·

I

IJ)

o

I

o·

I

to

O-+--------r----.------r-----.------, ·

I 0.0 0.1 0.2 0.3 0.4 0.5

z

Figure 2. Estimation of the region for which the relative

temperature varies linearly with z.
a.

a

N

- .....,

..'

o ........... ...... .... ..'

o

o o

-----------.

.'

.' .' o

.Ii

lY.:'

o

o

'6',<>', <>,' " I ,

I ' , " , , , , . , 0 '/ 0

: : " I "

'0 , 0, " , : , , "

: 0 '" I " , ' ,

I

,

I

I

I

I

I

I

I

,

a

a

o

a

a

OJ

o.

o

"t"

o -.J •

w@

>

lY. :' .-............ Current Simulation

o

a

N

o o

o

6

o

Grotzbach, Case 9, ReI. 7

Deardorff & Willis, ReI. 13

fj,:

a :

o +------r---

0.0 0.1

----,Ir--------"------TI-------,,

0.2 0.3 0.4 O.S

z

Figure 3. Comparison with experiments and previous

simulations of the RMS of the velocity components
C3

om

o

Current Simulation

O Grotzbach,Case 9, Fief. 7

u') O Deardorlf & Willis, Fief.13

_-_ c_-

c_ D

o

o

I I I I I

0.0 0.i 0.2 0.3 0.4 0.5

Z

Figure 4. Comparison with experiments and previous simulations

of the RMSof relative temperature.
2f5

a

en.

a

co

a

CD

a

l/)

.

0-t------..,.-----.,...-----~----___r----____.

0.0 0.1 0.2

z

0.3 0.4 0.5

Figure 5. wT correlation versus height.
27

c)

u-;-

Deardorff

&

Current

O Simulation Willis, Ref,13

"-' r7 Production

----- O Dissipation

i

A .......... A Diffusion

k

o I ...... O Molecular Transfer I

I._ I °°"

I o" ",

£0 _" ".

0 :1 "", [3

-- "l •

" I ".

__° ..... _=4 _: .....

LU l .._c___._ .

\ _ 0

I

I'I t

o_;- /°

,/

°

c; /

I

o

I I I I I I

0.0 0.1 0.2 0.3 0.4 0.5

Z

Figure 6. Terms in the kinetic energy equation versus height.
28

C:)

c)

c)-

('x]

o

x component

c)

(:D

y component

('kl

_" .......... z component I--

t,J I'

I-- o

O9

rn

co

... .... ° ..................................................

.........-*'°*'" °..*

• ..-

I I I I I

0.0 0.I 0.2 0.3 0.4 0.5

Z

Figure 7. RMS of the three vorticity components versus height.
29

=i'm

°°°° ..... o°o°°"

°°o..O.°°''°''"

°°

o..°

°°°

U_ °°

°

€5-

| °°

.o

.o

°

." <u 3 >/<u2>3 / 2

o°

.o

r,1 _ • °.. <v3>/<v2>3/2

_T' "°

[.1.-] °° o'

_,/ .. ............ <w3>/<w2> 3/2

0"-') *"

,"

o"

L_ ,°

o'

o°

°

,°

0 °'

I I I I I I

0.0 0.I 0.2 0.3 0.4 0.5

Z

Figure 8. Skewness factor of the three velocity components versus height,
30

c)

<u 4>/<u 2>2

c3 . <v4>/<v2>2

............ <w4>/<w2> 2

u'?

U') "

[.t._lo "

7 • "

_um-

(3Z

_]

L

• a

I | I I I

0.0 0.1 0.2 0.3 0.4 0.5

Z

Figure9. Flatness factorof the three velocitycomponentsversus height.
31

10 -_

u-z = 0.005

10-

a:::

I-- (O

IJJ

13_

0")

10-

0 -1 10 o 10 1

WAVENUMBER

Figure 10. One-dimensional wT cross spectra (averaged in y only)

versus x wavenumber.
32

10 oNu°z = 0.318

10":

iv"

I--

LLI

Q--10 "--

O3

I--

10-.

10- ' ' ' ' ' ' ' '1 ' ' ' ' ' _ r_l

10 -1 10 o 10

WAVENUMBER

Figure 11. One-dimensional wT cross spectra (averaged in y only)

versus x wavenumber.
33

10 o-

= . 7

10-

0")

10

10" ' ' " ' ' ' i"' I ' ' i ' ' ' ' I I

0 -_ 10 o 10 1

WAVENUMBER

Figure12. One-dimensional wT cross spectra (averagedin y only)

versus x wavenumber.
34

lO°-

Nu_z = 2.O4 '

10-

or"

I--

J.J

O_ 10-

10-

10- ' €' • ' ' ' ' '1 , t i" , , , i , i

10 -1 10 o 10 1

WAVENUMBER

Figure 13. One-dimensional wT cross spectra (averaged in y only)

versus x wavenumber.
35

o

N

o

Current Simulation

-e-- Carroll, Ref. 14

o

o

lJ')

o

o·

·O+-----r-.---.---;:,.T1rT-------r-.---..--rTTlrr---r-,-,--,rrTTrr----r-r-T1rTrrn

lJ')

.....

o·

,......... 0

0::::- E-to

0::::

Figure 14. Comparison with experiment of the RMS relative

temperature versus normalized height.
36

+|

rYE-. -

n_,_

CD_

_- Current Simulation

- j

!

o

I i i I I Ill I i I i I I I I i i I i i i I i I I i I i I I ! i i i I

_i0-3 I0-2 I0-_ I0° I0!

NU_Z

Figure 15. RMSrelative temperature versus normalized height.

The experimentally determined slopes of Carroll are shown

for comparison.
O

0

0

c_

-- Current Simulation

Carroll, Flel. 14

0

.2- _

I

I I i i i t III I i i I i I iii I i I I I i ill i I i i t i i"li i

i0-3 10-2 I0-I i0° iOt

NU_Z

Figure 16. Comparison with experiment of the skewness of the

reTative temperature versus n0rmalized height.
u_

Current Simulation

p-,

r,

Oo

U'3

#k]

o

I i I i I ii[J i ! i i I Ill I l i i J i ill i I i I i i fill

0-3 10-2 10-_ 18o 10_

NU_Z

Figure 17. Flatness fdctor of the relative temperature

versus normalized height,
-2

"0

_D

T" Current Simulation

(_)_

0 Carroll, I_et. 14

0

J i i i J J _i i J l _ a I I i I I i I i | I I li"J 0

10-3 [I]-2 I0-I 10

HlJ_Z

Figure 18. Comparison with experiment of the norma]ized relative

temperature derivative versus normalized height. The

s]ope shown is from the theory of Malkus.
O3

0

Current Simulation

_r. -0- Carroll, Per. 14

Figure lq. Comparison with expuri,,eflt or" the RHSof the

normalized relative temperature derivative.
u_

o

0

o

u_

I

IZ°

L_

I

-- Current Simulation

c::) /' --0--- Carroll, Fief. 14

r4

I

t._

, --, , , ,,,,,j , , ,,,,,,I , , , ,,,,,I , ' , ,'"',

0-3 ! 0-2 10-' 10° 101

NUxZ

Figure _0. Comparison with experiment of the skewness factor

for the normalized relative temperature derivative.
1. Report No. NASA CR-178027 12. Government Accession No. 3. Recipient's Catalog No.

ICASE Report No. 86-6

4. Title and Subtitle 5. Report Date

SIMULATION OF THE TURBULENT RAYLEIGH-BENARD PROBLEM February 1986

USING A SPECTRAL/FINITE DIFFERENCE TECHNIQUE 6. Performing Organization Code

7. Author{s) 8. Performing Organization Report No. T. M. Eidson, M. Y. Hussaini, and T. A. Zang 86-6

9. Prffwt1~l2'i!atfWioe6\itI%t~<t Alt'ifpi'fi cat ions

10. Work Unit No.

in Science

and Engineering

Mail Stop 132C, NASA Langley Research Center 11. C~Sf~~7~1a~N~ASl-18107

Hampton, VA 23665-5225

13. Type of Report and Period Covered 12. Sponsoring Agency Name and Address

Contractor Report

National Aeronautics and Space Administration 14. Sponsoring Agency Code

Washington, D.C. 20546 505-31-83-01

15. Supplementary Notes

Langley Technical Monitor: Proc. EUROMECH 199 Direct and Large

J. C. South Eddy Simulation of Turbulent Flows,

Munich, Germany, 9/30 - 10/2, 1985.

Final Report

16. Abstract

The three-dimensional, incompressible Navier-Stokes and energy equations

with the Bousines~ assumption have been directly simulated at a Rayleigh

number of 3.8 x 10 and a Prandtl number of 0.76. In the vertical direction,

wall boundaries were used and in the horizontal, periodic boundary conditions

were used. A spectral/finite difference numerical method was used to simulate

the flow. The flow at these conditions is turbulent and a sufficiently fine

mesh was used to capture all relevant flow scales. The results of the

simulation are compared to experimental data to justify the conclusion that

the small scale motion is adequately resolved.

17. Key Words (Suggested by Authors{s)) 18. Distribution Statement

Raleigh-Benard problem, 34 - Fluid Mechanics & Heat Transfer

direct turbulence simulation

Unclassified - unlimited

19. Security Classif.{of this report) 120. Security Classif.~of this page) 21. N

3

· of Pages 122. Price Unclassified Unclassifie 4 A03

For sale by the National Technical Information Service, Springfield, Virginia 22161

NASA-Langley, 1986