Session 15
NATURAL CONVECTION III
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, 21945970, 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^{1}.
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^{2}.
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 nonlinear 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^{6} and used in reference^{7}. 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 nonexpanded 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 ^{3} by using a FORTRAN code. More details and a variety of solved examples
are given in the forthcoming book^{8}.
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^{9}.
The Laplace equation is solved here by generating finite difference templates on a non
uniform grid. The rule, reported in ^{10}, that derives finite difference scheme is extended to include the
first and second kind boundary conditions.
In order to obtain optimal finite difference formulae^{11}, 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 ^{10}. 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
 S. Wolfram, The Mathematica Book, 3rd ed., Wolfram Media/Cambridge University Press, 1996.
 M. D. Mikhailov and M. N. Ozisik, Unified Analysis and Solutions of Heat and Mass Diffusion, John Wiley, NY 1984, Dover,1994.
 R. M. Cotta, Integral Transforms in Computational Heat and Fluid Flow, CRC Press, Boca Raton, Florida (1993).
 R. M. Cotta, The Integral Transform Method in Computational Heat and Fluid Flow, Special Keynote Lecture, Proc. of the 10th Int. Heat Transfer Conf., Brighton, UK, SK3, Vol.l, pp.4360, August 1994.
 R. M. Cotta, Integral transforms in transient convection: Benchmarks and Engineering Simulations, ICHMT, International Symposium on Transient Convection, 1923 August 1996, Cesme, Turkey.
 R. M. Cotta and M. D. Mikhailov, Heat Conduction Lumped Analysis, Integral Transforms, Symbolic Computation, John Wily,1997 (in press).
 M. D. Mikhailov, Mixed Coinputation in Transient Heat Convection, ICHMT, International Symposium on Transient Heat Transfer,1923 August 1996, Cesme, Turkey.
 M. D. Mikhailov and R. M. Cotta, Integral Transform Method with Mathematica (in preparation).
 M. N. Ozisik, Finite Differer:ce Methods in Heat Transfer, CRC Press,1994.
 M. D. Mikhailov, Finite difference method by using Mathematica , Int. J Heat Mass Transfer, 37, Sutipl.l, pp.375379,1994.
 R. Manohar and J. W. Stephenson, Optimal Finite Analytic Methods, Journal of Heat Transfer, Transactions of the ASME,104, 432437,1982.
THE STRUCTURE OF BUOYANCYDRIVEN MIXING OF TWO FLUIDS IN A CLOSED CAVITY
Mahmut D. MAT and Olusegun J. ILEGBUSI
Department of Mechanical, Industrial and Manufacturing Engineering
Northeastern University, Boston, USA
ABSTRACT
Buoyancydriven 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 twofluid 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 twofluid formulation which derives from the twophase 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 buoyancydriven flow by thermal and solutal gradients is rather extensive, these studies have been limited to perfectly premixed 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.^{1} has been extended to nonisothermal 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 breakup of the interface an produce chaotic mixing depending on the level of buoyancy. A twofluids 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.
REFERENCES
 Ilegbusi, O.J., Mat, M.D., and Andrews, M., The deformation and kinematic mixing of a collapsing interface, Appl. Math. Modelling, Vol. 21, pp 6676, 1997.
MEAN AND SECONDARY FLOW CONTRIBUTION TO THE HEAT TRANSFER IN AN INCLINED SLENDER BOX WITH AXIAL HEATING
R. D. Buscalioni^{*}, E. Crespo del Arco^{*}, P. Bontoux^{**}, J. Ouazzani^{***}.
^{*Depto. Fisica Fundamental, UNED, Madrid (Spain).
** I.M.F.M., C.N.R.S., Marseille, (France).
*** ArcoFluid., Aix en Provence, (France).
}
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 endtoend temperature difference have received an increasing attention since the last decade due to its relevance in a variety of technological and fundamental areas^{1}. 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^{23}. Natural convection in tilted fluid layers is also important in many geophysical situations^{45}. An interesting application resides in the transport and rate of spread of passive contaminants (for instance, radioactive material) in long tilted liquidfilled rock fractures. Woods et al. found that the contribution of the base convective flow to the transport rates is larger than that of diffusion^{5}. 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.
MEAN FLOW
Inside the inclined box, the gravity vector forms a certain angle with the imposed temperature gradient, so for nonvanishing 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^{6}.
The heat transfer is analized by calculating the Nusselt number along the crosssections 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 crossstream 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.
SECONDARY FLOW
The onset of secondary flow was investigated by studying the stability of the mean flow to small perturbations. The resulting OrrSommerfeld and heat equations for the amplitude of the perturbations were solved by a TauChebyshev method. The stability of the Hadley circulation in limiting horizontal configuration was studied by Hart^{7}, Laure et al.^{8} and Kuo et al.^{9}. 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 nonhorizontal tilts, the secondary motion arises as thermal long transversal oscillatory rolls^{10} (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 aPr 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 cutoff value for the wavelength of the secondary rolls^{11} and also as it regulates the influence of the recirculation eddies at the end regions on the flow at the core^{12}. 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^{12}.
NUMERICAL SIMULATION OF THE FLOW
The numerical part of the calculations were done by a collocation pseudospectral method. To study the twodimensional base flow and the transversal instability we solved the unsteady NavierStokes and heat equations in vorticitystream 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 streamfunction 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 crossstream 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 AdamBashforth, second order Backward Euler scheme^{13}. For each time cycle, the equation for the temperature at the next time step consists on a Helmholtztype equation which is solved by means of a double diagonalization procedure for the algebraic system that arise from the collocation method^{14}. The corresponding equations for the streamfunction and vorticity consist on a Stokestype problem. This is solved by using the Influence Matrix technique^{14} that leads to the solution of several Helmholtz equations with Dirichlet boundary conditions. REFERENCES
 Daniels, P. G. , Wang P., On the evolution of thermally driven shallow cavity flows, J. Fluid Mech., Vol. 259, pp. 107124, 1994.
 Markham, B. L. and Rosenberger, F., Diffusiveconvective vapor transport across horizontal and inclined rectangular enclosures, J. Crystal Growth Vol. 67, pp. 241254, 1984.
 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 threedimensional buoyancy driven flows in tilted cylinders, in Advances in Space Research, Pergamon, 1986, pp. 155160.
 Woods, A. W. and Lintz S. J., Natural convection and dispersion in a tilted fracture. J. Fluid Mech., Vol. 241, pp. 5974, 1992.
 Cessi, P. and Young, W. R., Fixedflux convection in a tilted slot, J. Fluid Mech., Vol. 237, pp. 5771, 1992.
 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.
 Hart, J. E., A note on the stability of lowPrandtlnumber Hadley circulations, J. Fluid Mech, Vol 132, p. 271, 1983.
 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.
 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
 Buscalioni R. D., Ouazzani, J., Bontoux P. and Crespo del Arco, E., Unicellular and multicellular flows in inclined cavities driven by an endtoend temperature difference, Proceedings of Second European Symposium, Fluid in Space. Naples, April 2226, 1996, pp. 376382.
 Janssen, R. J. A. and Henkens, R. A. W. M., Instabilities in threedimensional diffentially heated cavities with adiabatic horizontal walls. Phys. Fluids, Vol. 8, p. 62, 1996.
 Drummond, E. and Korpela, S. A., Natural convection in a shallow cavity. J. Fluid Mech Vol., 182, p. 543, 1987.
 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. 115128, 1984.
 Vanel, J. M., Peyret, R. and Bontoux, P., A pseudospectral olution of vorticitystream function equations using the influence matrix technique, Numerical Methods for Fluid Dynamics II, K. W. Morton and M. J. Baines, Eds., 463475, Clarendom Press, Oxford, 1986.
A NUMERICAL STUDY OF STABILITY OF SQUARES IN RAYLEIGHBENARD 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, RayleighBenard convection represents the simplest fluid system
which exhibits a sequence of transitions from twodimensional laminar to more complicated
threedimensional and finally turbulent convection. It has been shown experimentally (White^{1}) and
theoretically (Busse & Frick^{2} and Christensen & Harder^{3}) that square convective patterns become
stable when there is variation of viscosity.
The term "square" has been used originally for threedimensional 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 GlycerolWater 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 vorticityvector
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^{4}. Numerical solutions are presented for Rayleigh
numbers above the onset of convection up to the limit of timedependence and viscosity ratios
between 1 and 50. The results are a small part of a detail numerical and experimental study of
RayleighBenard 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^{1} and Christensen & Harder^{3}). 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
GlycerolWater mixture with 80% Glycerol and for Ra_{avt} = 2250. The wavenumber of the squares
shown in Figure 2, according to the definition ( , is the nondimensional 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 viscositytemperature relation. When the
convecting fluid has an exponential viscositytemperature 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 viscositytemperature 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 midplane 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
 White, D.B., The planform and onset of convection with a temperaturedependent viscosity. J. Fluid Mech., Vol.191, pp 247286,1988.
 Busse, F.H. & Frick, H., Square pattern convection in fluids with strongly temperature dependent viscosity. J. Fluid Mech. Vol.150, pp 451465,1985.
 Christensen, U. & Harder, H., 3D convection with variable viscosity. Geophys. J. Int., Vol. 104, pp 213226,1991.
 Khorasanizadeh, H., A Numerical and Experimental Study of RayleighBenard 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 wellknown 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 xdirection
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.
Computational domain  L_{z} x L_{y} x L_{x} =12.0 x 6.0 x 1.0
 Number of grid points  N_{z} x N_{y} x N_{X} = 432 x 216 x 96
 Resolution  dX_{min} =0.00044 h and dX_{max} =0.0131 h
 Rayleigh numbers computed  Ra = 5.4e5, 8.23e5, 2.0e6, 5.0e6

To give an impression of our results, we show in fig.2 the average velocity
and temperature profiles for the four Rayleigh numbers that we have computed.
The results agree quite well with the few experimental data that are available.
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.
^{1}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.
