Session 15

Chairman: R. L. Sani

COMPUTATIONAL HEAT TRANSFER WITH MATHEMATICA

M. D. Mikhailov
Mechanical Engineering Department
Laboratório de Transmissão et Tecnologia do Calor - LTTC
EE / COPPE / UFRJ, Universidade Federal do Rio de Janeiro
Cidade Universitária, Rio de Janeiro, RJ, 21945-970, BRASIL
mikhailov@serv.com.ufrj.br

Traditionally, computational methods are considered as analytical and numerical. In reality all methods are mixed - analytical and numerical, and becomes more powerful when using the Mathematica software system¹.

Then the results are presented in a notebook format. In contrast to standard text, a notebook is Web publishable interactive document. It can be used not only to be read, but also to execute the input statements, plot graphics and animate the results. The inputs can be modified to obtain new results.

The goal of this lecture is to show that Mathematica system is a natural environment for computational heat transfer. Some new perspectives for integral transform method (ITM) and finite difference method (FDM) are demonstrated.

The ITM transforms linear partial differential equations (PDE) into a system of algebraic equations (AE) or a system of ordinary differential equations (ODE). The resulting system is then solved analytically and an inversion formula is used to obtain the solution of the original problem².

During the last decade the potential of the ITM was considerably expanded by solving the resulting AE or ODE system numerically. This was an extremely important extension since it permitted non-linear PDE also to be treated. The reader should consult for the latest development in this generalized integral transform technique (GITT)^3,5.

One Mathematica version of GITT, when the inversion formula is expanded at the very beginning, is developed in the book⁶ and used in reference⁷. This version is easy to apply but it is convenient only for small number of eigenvalues.

Another more advanced integral transform technique is presented here. As an example the transient temperature distribution in a slab with temperature dependent thermal conductivity is found. The ODE system is derived in non-expanded form which permitted efficient generation of the truncated ODE system. The integrals that appear in ODE system are solved exactly by using a new approach. It consists of generating an algebraic system for the unknown integrals by using an integral identity and some properties of the eigenfunctions. The truncated ODE system is solved numerically and used in the inversion formulae. The numerical results for the temperature coincide with those obtained in ³ by using a FORTRAN code. More details and a variety of solved examples are given in the forthcoming book⁸.

A major drawback of the FDM appears to be in its inability to handle effectively the problems over arbitrarily shaped complex geometry, while the finite element method (FEM) has the flexibility in dealing with such tasks.

Only recently, with the advance of numerical grid generation approach, irregular geometries are transformed to regular computational domains in order to use FDM⁹.

The Laplace equation is solved here by generating finite difference templates on a non uniform grid. The rule, reported in ¹⁰, that derives finite difference scheme is extended to include the first and second kind boundary conditions.

In order to obtain optimal finite difference formulae¹¹, the templates are derived by using exact polynomial solutions of the Laplace equation. An efficient rule for generation of such polynomials is given.

The final Mathematica finite difference solution gives numerical values for any specified location in the solution domain.

This new finite difference technique is applicable to many other problems and extends the one presented in ¹⁰. With this approach the FDM becomes comparable to FEM in dealing with irregular geometry.

ACKNOWLEDGEMENTS

The kind hospitality of the Mechanical Engineering Department at EE/COPPE, Universidade Federal do Rio de Janeiro, and the financial support provided by FAPERJ and CNPq, both in Brazil, is gratefully acknowledged.

REFERENCES

1. S. Wolfram, The Mathematica Book, 3rd ed., Wolfram Media/Cambridge University Press, 1996.
2. M. D. Mikhailov and M. N. Ozisik, Unified Analysis and Solutions of Heat and Mass Diffusion, John Wiley, NY 1984, Dover,1994.
3. R. M. Cotta, Integral Transforms in Computational Heat and Fluid Flow, CRC Press, Boca Raton, Florida (1993).
4. R. M. Cotta, The Integral Transform Method in Computational Heat and Fluid Flow, Special Keynote Lecture, Proc. of the 10-th Int. Heat Transfer Conf., Brighton, UK, SK-3, Vol.l, pp.43-60, August 1994.
5. R. M. Cotta, Integral transforms in transient convection: Benchmarks and Engineering Simulations, ICHMT, International Symposium on Transient Convection, 19-23 August 1996, Cesme, Turkey.
6. R. M. Cotta and M. D. Mikhailov, Heat Conduction- Lumped Analysis, Integral Transforms, Symbolic Computation, John Wily,1997 (in press).
7. M. D. Mikhailov, Mixed Coinputation in Transient Heat Convection, ICHMT, International Symposium on Transient Heat Transfer,19-23 August 1996, Cesme, Turkey.
8. M. D. Mikhailov and R. M. Cotta, Integral Transform Method with Mathematica (in preparation).
9. M. N. Ozisik, Finite Differer:ce Methods in Heat Transfer, CRC Press,1994.
10. M. D. Mikhailov, Finite difference method by using Mathematica , Int. J Heat Mass Transfer, 37, Sutipl.l, pp.375-379,1994.
11. R. Manohar and J. W. Stephenson, Optimal Finite Analytic Methods, Journal of Heat Transfer, Transactions of the ASME,104, 432-437,1982.

Buoyancy-driven convection in a closed cavity has received considerable attention in recent years due to its practical relevance in many applications including materials processing (such as crystal growth), oceanography, and atmospheric transport phenomena. Most of previous theoretical studies have considered relatively small density variations due to the limitation of Boussiniseq approximation. The present study is concerned with the development and application of a novel two-fluid model that is not subject to such limitation. A model problem is considered, consist of two fluids which initially meet at a sharp density interface and mix under the combined effects of thermal and solutal buoyancy forces. This problem is closely related to thermosolutal convection in partitioned cavities. It is also relevant to material processing in microgravity environment. The mathematical model is based on a two-fluid formulation which derives from the two-phase flow theory.

Studies on thermosolutal convection in cavities subjected to horizontal temperature and concentration gradients have shown that a layered structure develops and at high thermal and solutal Grashof numbers, the flow is unsteady and transitional. the fluid flow, heat and mass transfer characteristics of such systems are determined by the following dimensionless groups: thermal Grashof number (Gr_T), solutal Grashof number (Gr_S), Lewis number (Le), Prandtl number (Pr), and aspect ratio (Ar). Although the literature on buoyancy-driven flow by thermal and solutal gradients is rather extensive, these studies have been limited to perfectly pre-mixed binary fluid mixture. The mixing of two fluids under effect of opposing buoyancy forces has not received much attention. This is the problem being addressed in the present study. The previous theoretical studies generally assume Boussinesq fluids and address mixing by considering the kinematics of a line or interfacial area due to a given flow field from a continuum mechanics point of view. Mixing has also been addressed with a Lagrangian approach in terms of iterated map models. The experimental studies have mainly considered mixing of fluids agitated through mechanical means and which involve motion of a boundary.

In the present study, the work of Ilegbusi et al.¹ has been extended to non-isothermal systems. The problem considered involves the mixing of two fluids initially separated by a sharp density interface in a closed cavity. The effects of Atwood number, buoyancy ratio and Prandtl number on interface behavior and flow field are studied. The evolution of the interface is quantified and it is shown that the break-up of the interface an produce chaotic mixing depending on the level of buoyancy. A two-fluids model is employed to calculate the zone averaged velocities, temperatures, and volume fraction of the fluids. These two fluids are in effect assumed to share occupancy of space but retain their identity throughout the mixing process. Transport equations are solved for flow the variables and allowance is made for interface transfer of momentum and energy. The calculated volume fraction (or existence probability of each fluid) coupled with the interface evolution provides a measure of the degree of mixing in the cavity. A van Leer scheme is used to resolve the sharp density interface between two fluids. It is well known that numerical methods with upwind schemes generally suffer from numerical diffusion which can act as an additional source terms to produce unphysical results. A higher order scheme can reduce this problem. However, due to the sharp density gradient in the present problem, a conventional high order scheme may cause numerical oscillations. For example, the volume fraction for each fluid should lie between 0.0 and 1.0 at every time step. Undershoots or overshoots of these values give unrealistic results and adversely affect the momentum equation since interface friction term involves a product of volume fractions. The van Leer scheme is used to prevent such undershoots and overshoots while maintaining high order numerical accuracy.

The preliminary results show that solutal buoyancy plays a dominant role in mixing process for the system considered. The material interface requires a longer time to develop for pure thermal convection compared to solutal convection. Interface elongates for a moderate buoyancy ratio two times more than pure solutal convection providing very large mixing area resulting an efficient mixing.

For small levels of Grashof numbers which may be achieved under microgravity conditions, the flow field and interface behavior are not much affected by buoyancy forces. Such situation would be ideal for crystal growth process in which flow control and minimal convection are desirable.

The analysis of the system dynamics shows that for high Grashof numbers, the flow field is aperiodic, and the phase space trajectory is irregular, and the power spectrum indicates that the flow field responds in a broadband frequency characteristic of chaotic flow.

1. Ilegbusi, O.J., Mat, M.D., and Andrews, M., The deformation and kinematic mixing of a collapsing interface, Appl. Math. Modelling, Vol. 21, pp 66-76, 1997.

MEAN AND SECONDARY FLOW CONTRIBUTION TO THE HEAT TRANSFER IN AN INCLINED SLENDER BOX WITH AXIAL HEATING

We present a numerical and theoretical study on the convective flow in an inclined box driven by a temperature gradient along their longest axis. Natural convection in shallow cavities driven by end-to-end temperature difference have received an increasing attention since the last decade due to its relevance in a variety of technological and fundamental areas¹. Most part of the published works on this subject have considered horizontal configurations although the inclined geometry is relevant in a variety of applications in industry and nature. For instance, in the process of crystal growth from vapor, larger transport rates are obtained by tilting the ampoule a certain angle respect to gravity^2-3. Natural convection in tilted fluid layers is also important in many geophysical situations^4-5. An interesting application resides in the transport and rate of spread of passive contaminants (for instance, radioactive material) in long tilted liquid-filled rock fractures. Woods et al. found that the contribution of the base convective flow to the transport rates is larger than that of diffusion⁵. They also asses the importance of the secondary flow in modifying the overall mass and heat flux, though they did not considered this part of the problem.

We consider an incompressible fluid filling a rectangular slender box whose dimension are HxDxL. The box is differentially heated along the longest side, L which is inclined an angle a respect to the gravity vector. The shape of the flow depends on the external Rayleigh number, R, (defined in terms of the overall temperature gradient), the Prandtl number, Pr, the inclination, a and the aspect ratios of the cavity.

Inside the inclined box, the gravity vector forms a certain angle with the imposed temperature gradient, so for non-vanishing Rayleigh number and inclination angle, a basic unicellular convective flow arises in response to the torque created by the buoyancy force. The shape of the mean velocity and temperature profiles are analytically obtained by assuming plane parallel flow in the core of the slender box. To calculate the amplitude of the basic flow in terms of the external Rayleigh number, the flow solution at the core is matched with constructed solutions at the end regions, near the closing walls, using a method similar to that proposed by Bejan and Tien⁶.

The heat transfer is analized by calculating the Nusselt number along the cross-sections of the cavity. In the conductive regime the Nusselt number increase almost linearly with the Rayleigh number and decrease with the inclination like sin(a). The advective regime appears with a sudden increase of the Nusselt number, for values of the Rayleigh number near a certain threshold value which essentially depends on the inclination. As long as the isothermals begin to be distorted, the flow intensity increases steeply due to a feedback mechanism that couples the mean velocity and temperature fields. The mean flow is amplified by the upslope buoyancy acting upon the cross-stream temperature differences while these, in turn, increase with the mean flow advection . Consequently to this fact in the advective regime, the largest heat transport rates along the upslope direction are found at intermediate inclinations (about 50 degrees) as shows the dependence of the Nusselt number with the inclination. As the external Rayleigh number further increases, the isothermals concentrates in a narrow region near the end walls, the local Rayleigh number at the core tends to a constant limiting value and the Nusselt number becomes almost independent on the inclination.

As a consequence of the form of the constructed solutions at the end regions, in the theoretical model, the Nusselt number is independent on the Prandtl number. The theoretical results were compared with results from the numerical calculation for Pr=0.7 (air) and Pr=6.7 (water). In the conductive and advective regimes a very good agreement was found. For larger Rayleigh number the Nusselt number presents a slight dependence with the Prandl number.

The onset of secondary flow was investigated by studying the stability of the mean flow to small perturbations. The resulting Orr-Sommerfeld and heat equations for the amplitude of the perturbations were solved by a Tau-Chebyshev method. The stability of the Hadley circulation in limiting horizontal configuration was studied by Hart⁷, Laure et al.⁸ and Kuo et al.⁹. The onset of secondary flow is due to transversal or longitudinal instabilities. For low enough Pr, transversal rolls appear in the form of shear stationary nearly squared rolls. For larger Pr and non-horizontal tilts, the secondary motion arises as thermal long transversal oscillatory rolls¹⁰ (with wavelength of about 10 H). In the case of perfectly conducting walls another oscillatory thermal mode with short wavelength (about 0.8 H) is found at large inclinations. The stability boundaries in the a-Pr space reveals that the critical transverse mode depends mostly on the Prandtl number. Longitudinal modes are driven by buoyancy though a part of their energy may come out from the mean flow. For almost all inclinations, the critical longitudinal instability consists on a stationary mode with very long wavelength (about 100 H). Near horizontal positions oscillatory longitudinal modes are present for moderate Pr and for larger values of Prandtl short wave stationary longitudinal rolls are dominant. The effect of secondary motions in the heat transport rate has been considered. The growth of the averaged Nusselt number with the Rayleigh number is reduced when the multicellular flow is established. The effect of aspect ratio in the onset of multcellular flow is important as it imposes a cut-off value for the wavelength of the secondary rolls¹¹ and also as it regulates the influence of the recirculation eddies at the end regions on the flow at the core¹². In this sense, we have found that in the case of transversal secondary rolls, the inclined closed geometry induces imperfect bifurcations to multicellular flow at Prandtl numbers one order of magnitude larger than in the case of the horizontal cavity¹².

The numerical part of the calculations were done by a collocation pseudospectral method. To study the two-dimensional base flow and the transversal instability we solved the unsteady Navier-Stokes and heat equations in vorticity-stream function variables. The longitudinal rolls were studied by assuming that the mean flow is independent on the upslope direction and solving for the upslope perturbation velocity, the temperature and the stream-function in the plane perpendicular to the base flow and the upslope direction. The equations of motion were solved for closed geometry and also for periodic boundary conditions on the cross-stream perpendicular direction. For the closed geometry, the spatial approximation in both directions was done by expanding the flow variables in truncated series of Chebyshev polynomials while in the case of periodic boundary conditions, truncated Fourier expansions were used for the periodic direction. The time discretization was obtained through an Adam-Bashforth, second order Backward Euler scheme¹³. For each time cycle, the equation for the temperature at the next time step consists on a Helmholtz-type equation which is solved by means of a double diagonalization procedure for the algebraic system that arise from the collocation method¹⁴. The corresponding equations for the stream-function and vorticity consist on a Stokes-type problem. This is solved by using the Influence Matrix technique¹⁴ that leads to the solution of several Helmholtz equations with Dirichlet boundary conditions.

1. Daniels, P. G. , Wang P., On the evolution of thermally driven shallow cavity flows, J. Fluid Mech., Vol. 259, pp. 107-124, 1994.
2. Markham, B. L. and Rosenberger, F., Diffusive-convective vapor transport across horizontal and inclined rectangular enclosures, J. Crystal Growth Vol. 67, pp. 241-254, 1984.
3. Bontoux, P., Smutek, C., Randriamampianina, A. , Roux B., Extremet G. P., Hurford, A. C., Rosemberger, F. and De Vahl Davis, Numerical solutions and experimental results for three-dimensional buoyancy driven flows in tilted cylinders, in Advances in Space Research, Pergamon, 1986, pp. 155-160.
4. Woods, A. W. and Lintz S. J., Natural convection and dispersion in a tilted fracture. J. Fluid Mech., Vol. 241, pp. 59-74, 1992.
5. Cessi, P. and Young, W. R., Fixed-flux convection in a tilted slot, J. Fluid Mech., Vol. 237, pp. 57-71, 1992.
6. Bejan, A. and Tien, C. L., Laminar natural convection heat transfer in a horizontal cavity with different end temperatures, J. Heat Transfer, Vol. 100, p. 641, 1978.
7. Hart, J. E., A note on the stability of low-Prandtl-number Hadley circulations, J. Fluid Mech, Vol 132, p. 271, 1983.
8. Laure, P. and Roux, B., Synthèse des résultats obtenus par l’étude de stabilité des mouvements de convection dans une cavité horizontale de grande extension., C. R. Acad. Sci. Paris, Vol., 305 Série II, 1137, 1987.
9. Kuo, P. and Korpela, S. A., Stability and finite amplitude natural convection in a shallow cavity with insulated top and bottom and heated from a side, Phys. Fluids Vol, 33, No. 1, p. 33, 1988
10. Buscalioni R. D., Ouazzani, J., Bontoux P. and Crespo del Arco, E., Unicellular and multicellular flows in inclined cavities driven by an end-to-end temperature difference, Proceedings of Second European Symposium, Fluid in Space. Naples, April 22-26, 1996, pp. 376-382.
11. Janssen, R. J. A. and Henkens, R. A. W. M., Instabilities in three-dimensional diffentially heated cavities with adiabatic horizontal walls. Phys. Fluids, Vol. 8, p. 62, 1996.
12. Drummond, E. and Korpela, S. A., Natural convection in a shallow cavity. J. Fluid Mech Vol., 182, p. 543, 1987.
13. Haldenwang, P., Labrosse, G., Abboudi, S. and Deville, M., Chebyshev 3D spectral and 2D pseudospectral solvers fort the Helmholtz equation, J. Comp. Phys., Vol. 55, pp. 115-128, 1984.
14. Vanel, J. M., Peyret, R. and Bontoux, P., A pseudospectral olution of vorticity-stream function equations using the influence matrix technique, Numerical Methods for Fluid Dynamics II, K. W. Morton and M. J. Baines, Eds., 463-475, Clarendom Press, Oxford, 1986.

A NUMERICAL STUDY OF STABILITY OF SQUARES IN RAYLEIGH-BENARD CONVECTION

Hossein Khorasanizadeh*, Eddie Leonardi* & John Reizes**
* School of Mechanical and Manufacturing Engineering, The University of New South Wales, Sydney, Australia 2052
** School of Mechanical Engineering, University of Technology, Sydney, 2007, Australia

ABSTRACT

In addition to having many practical thermal engineering applications such as thermal comfort, solidification processes and crystal growth, solar energy collectors, electronic equipment cooling and earth mantle convection, Rayleigh-Benard convection represents the simplest fluid system which exhibits a sequence of transitions from two-dimensional laminar to more complicated three-dimensional and finally turbulent convection. It has been shown experimentally (White¹) and theoretically (Busse & Frick² and Christensen & Harder³) that square convective patterns become stable when there is variation of viscosity.

The term "square" has been used originally for three-dimensional convective patterns with similar wavelengths in the two horizontal directions, even though this is not the general feature of the squares since they can exist with different horizontal wavelengths. In this paper stability of squares in both infinite and bounded domains are discussed.

A cavity with symmetry imposed on its lateral walls and a bounded cavity with horizontal aspect ratios of 9 x 3 (90 mm by 30 mm, see Figure 1 ) have been considered. The width of the Perspex lateral walls is 10 mm and the top and bottom walls are 2 mm thick glass. Experiments have been performed in the same cavity with Glycerol-Water mixtures. In order to investigate the effects of conduction through lateral walls, in the numerical study these walls have been considered to be (i) adiabatic and (ii) conducting.

The equations of motion of fluid and energy equations have been recast into a vorticity-vector potential formulation and then solved using a false transient method. A uniform mesh size of 0.0714 is used in all the three directions in the fluid part. This yields a mesh of 127x15x43 for the geometry shown in Figure 1. A similar mesh size in the solid lateral walls results in 15 mesh points and a mesh size of 0.0666 used in the top and bottom walls yields 4 mesh points. Detail of the governing equations, boundary and interface conditions, solution procedure, code validation and mesh refinement are given in Khorasanizadeh⁴. Numerical solutions are presented for Rayleigh numbers above the onset of convection up to the limit of time-dependence and viscosity ratios between 1 and 50. The results are a small part of a detail numerical and experimental study of Rayleigh-Benard convection in infinite and bounded domains.

The numerical results are presented in terms of: a Rayleigh number, Ra_avt, which is based on a thermally averaged viscosity, , viscosity ratio, , and a deviation ratio, .

Different from the structure of the rolls and bimodal patterns, the main feature of square patterns is that the central ascending region in each square is entirely surrounded by the descending region. For fluids with a viscosity that decreases with increasing temperature, the flow rises at the centre of each square and descends close to its edges (White¹ and Christensen & Harder³). The circular rising and the surrounding descending regions of the squares are indicated in Figure 2 which shows the isovels of the vertical velocity, v, on a horizontal plane at y = 0.43 in a solution which employed Glycerol-Water mixture with 80% Glycerol and for Ra_avt = 2250. The wavenumber of the squares shown in Figure 2, according to the definition ( , is the non-dimensional width of a complete square), is almost 2.8.

A study of the structure and stability of the squares for viscosity ratios between 1 and 50, 0.7 < r_d < 2.3 and Ra_avt between 2000 and 26000 revealed that they are strongly affected by the details of the viscosity variation. It has been shown that a minimum viscosity variation is required for the stability of the squares depending on the type of viscosity-temperature relation. When the convecting fluid has an exponential viscosity-temperature relation, a great deal of viscosity change occurs in the vicinity of the cold boundary and the rest of layer becomes almost isoviscous. However, in the case of an inverse or a linear viscosity-temperature relation the viscosity change is distributed across the layer and a stagnant high viscous and thick region close to the cold wall develops.

Squares with different horizontal wavelengths were obtained in the bounded cavity, when the lateral walls were assumed adiabatic (see Figure 3). The wavelength of the squares in the X direction is 2.25 while it is 3 in the Z direction. When the realistic boundary and temperature conditions were imposed (actual Perspex walls and convection to ambient) the squares became unstable and were transformed to rolls (see Figure 4). This is in agreement with the experiments in which squares were unstable. The temperature distribution on the X vertical mid-plane in the fluid and solid parts shown in Figure 5 is indicative of the effect of conduction through lateral walls that leads to the instability of squares.

REFERENCES

1. White, D.B., The planform and onset of convection with a temperature-dependent viscosity. J. Fluid Mech., Vol.191, pp 247-286,1988.
2. Busse, F.H. & Frick, H., Square pattern convection in fluids with strongly temperature dependent viscosity. J. Fluid Mech. Vol.150, pp 451-465,1985.
3. Christensen, U. & Harder, H., 3-D convection with variable viscosity. Geophys. J. Int., Vol. 104, pp 213-226,1991.
4. Khorasanizadeh, H., A Numerical and Experimental Study of Rayleigh-Benard Convection, PhD Thesis, University of N.S.W., Sydney, Australia,1997.

DNS OF NATURAL CONVECTION BETWEEN TWO VERTICAL DIFFERENCIALLY HEATED WALLS

T.A.M. Versteegh and F.T.M. Nieuwstadt
J.M. Burgers Centre
Delft University of Technology
Laboratory for Aero- and Hydrodynamics
Rotterdamseweg 145
2628 AL Delft
The Netherlands

The investigation of turbulence is most conveniently done in simple flow con- figurations where influence of various physical processes can be isolated and studied in detail. One of these well-known geometries is the flow between two infinite walls, also denoted as a channel. Here we consider this geometry in the form of two vertical, infinite walls which are kept different temperatures, i.e. one of the walls has a higher temperature with respect to the other. This geometry is illustrated in fig.l.

As a result of the temperature difference between the two walls a natural convection flow develops. If the temperature difference^l becomes sufficiently large, the flow is turbulent. This type of turbulent flow is interesting from a fundamental point of view, because the effect of shear and buoyancy in this case acts in the same direction, but also from an applied point of view, because this type of flow is found in many practical heat transfer applications. The objective of our investigation is to study this flow by means of direct numerical simulation (DNS).

As a result of the infinite extent in the vertical direction, one can consider the flow as being homogeneous in that direction. As a consequence, we can apply our DNS in a finite, but sufficiently large vertical domain by using periodic boundary conditions. Moreover, homogeneity implies that all statistical quantities are only a function of the distance x between the two walls.

The DNS code that we use, is based on a standard finite volume method with help of second order discretisation schemes. By applying these schemes on two grids, which difFer by approximately by a factor two, we can reach fourth order accuracy of the solution through application of the Richardson extrapolation technique. By means of testing on several grids, this technique has proven to be reliable for our type of flow. The grid size is variable along the x-direction with the smallest grid size near the wall and the largest grid size in the centre of the channel. The grid in the other two coordinate directions is chosen uniform. Some other computational details are collected in the table below.

The results of these simulations have been used to compute practically all separate components in the budget terms of various Reynolds stresses. These budgets can be used to validate closure models for natural convection flows. As an example we show a budget of the -stress, i.e. the budget of the turbulent kinetic energy in the vertical direction multiplied by two.

In our presentation we will discuss the features of these stress budgets and some of the implications for turbulence modelling. One of the results that we have found, is that near the wall the shear production term in the energy budget becomes negative. This invalidates the use of all turbulence models based on an eddy viscosity approach. Comparison with stress budgets in linear stability analysis shows that several features of the turbulent budget terms already are prominent in the first instabilities. Finally we will give outlines for scaling the budgets in the inner region of the flow.

¹A better parameter to characterise the flow is the Rayleigh number with b the thermal expansion coefñcient, n the viscosity and k the heat conduction coefficient.