Applied Scientific Research 47: 195-220, 1990. © 1990 Kluwer Academic Publishers. Printed in the Netherlands.

195

On the stability of the natural convection flow in a square cavity heated from the side R.A.W.M. HENKES & C.J. H O O G E N D O O R N Delft University of Technology, Department of Applied Physics, P.O. Box 5046, 2600 GA Delft, The Netherlands Received 5 April 1989; accepted in revised form 28 August 1989

Abstract. The stability of the steady laminar natural-convection flow of air (Prandtl number

0.71) and water (Prandtl number 7.0) in a square cavity is calculated by numerically solving the unsteady, two-dimensional Navier Stokes equations. The cavity has a hot and cold vertical wall and either conducting or adiabatic horizontal walls. The flow looses its stability at a lower Rayleigh number in the case of conducting horizontal walls than in the case of adiabatic horizontal walls. The flow of water is more stable than the flow of air. Directly above the critical Rayleigh number the unsteady flow shows a single-frequency oscillation. Air in the case of adiabatic horizontal walls is an exception and shows two frequencies. The instabilities in the cavity seem to be related to well-known elementary instability mechanisms. In the case of conducting and adiabatic horizontal walls the instability seems to be related to a Rayleigh/ B6nard and a Tollmien-Schlichting instability respectively. The second instability for air in the case of adiabatic horizontal walls seems to be related to an instability after a hydraulic jump.

1. Introduction

The natural-convection flow in an enclosure undergoes a transition from a laminar to a turbulent state when the temperature difference (or dimensionless: the Rayleigh number) between the two vertical walls exceeds a critical value. This transition can be calculated by numerically solving the unsteady, three-dimensional Navier-Stokes equations. Because of the appearance of very fine eddy structures in the final phase of the transition, this calculation would require an enormous computational effort. Therefore in the present study we restrict ourselves to the calculation of the initial phase of the transition: we calculate how the steady, two-dimensional laminar flow looses its stability by solving the unsteady, two-dimensional Navier-Stokes equations. The stability can be examined by following the evolution of disturbances on a steady laminar flow. Beyond a critical value these disturbances no longer die out and a bifurcation to an unsteady flow is found. For very simple geometries some classical results exist, referring to elementary instability

196

R.A.W.M. Henkes & C.J. Hoogendoorn

mechanisms. For example we know (i) the Rayleigh/B6nard instability above an infinite hot horizontal plate, (ii) the instability after a hydraulic-jump in a horizontal moving layer and (iii) the Tollmien-Schlichting instability in a boundary layer. In the present approach three-dimensional disturbances are excluded. Three-dimensional effects are certainly important in the fully turbulent flow. Whether the three-dimensional disturbances lead to an earlier instability of the steady, two-dimensional flow than the two-dimensional disturbances remains to be investigated. The unsteady, two-dimensional Navier-Stokes equations have been solved numerically for a square cavity that is heated differentially over the vertical walls. The stability is calculated for both the case of conducting and adiabatic horizontal walls. Both air and water are considered. The stability for air in the square cavity in the case of conducting horizontal walls has already been calculated by Winters [1] and by Le Qu6r6 & Alziary de Roquefort [2, 3]. Further, the latter authors [4] also calculated the stability for air in a cavity with adiabatic horizontal walls and a height-width ratio of 4 and larger. We do not know of stability calculations for air in the square cavity with adiabatic horizontal walls, nor of any stability calculations for water in the cavity. Although the cavity configuration is more complicated than any of the mentioned simple configurations, the present calculations suggest that mechanisms leading to instabilities in the cavity are closely related to the Rayleigh/B6nard instability, the instability after a hydraulic-jump and the Tollmien-Schlichting instability. A better knowledge of the instability mechanisms is important, because they form the roots of turbulence: it can help us to better model the turbulence. 2. Mathematical formulation

The flow in the square cavity, as depicted in Fig. 1, is described by the unsteady, two-dimensional, incompressible Navier-Stokes equations Du 073 + -- 0, 8x 82 8u 8u 8u 8-7 -}- bl ~X "-~ 7d 8y

8v

8v

_

8v _

8-7 + u Uxx + 73 8y

a-t- + U ~ x + v 8y -

18p {82u 8Zu) p 8X -[- F \ OX2 -[- Oy~ ) ,

1@ p a y + g/ (r -

P r k S x 2 + 8y2] '

(82v r,) + v \ ax 2 +

82v~ ],

(1)

Stability oJ natural convection flow i,,l

197

H

///

/ //

/ / / /

/

./ /

/

H

/

/ /

r~ u

Th

/ / /

/

~£///////// Fig.1.

/

~x

Geometry of the problem.

(t time; x, y coordinates; u, v velocities; T temperature; p pressure; p density; v kinematic viscosity; fi coefficient of thermal expansion; Pr Prandtl number; T,- reference temperature; g gravitational acceleration). The Boussinesq aproximation has been applied. The cavity has a hot left vertical wall (T~,) and a cold right vertical wall (T~.). All the walls are rigid. Two cases are considered: the case with conducting horizontal walls (i.e. these walls have a linear temperature profile from the hot to the cold vertical wall) and the case with adiabatic horizontal walls (c?T/@ = 0). The variables are non-dimensionalized with the length scale x, = H, the velocity scale u, = x/--gfiATH, the time scale t, = H2/v, the temperature scale T,, = T,. and the temperature difference AT = T~, - T,.. This gives

{u

v T-

u,' u~'

T,

AT

p } =

'Ou~

f{~ "

x y_,_,Ra, Pr}.

(2)

'x~ x,

Two dimensionless numbers appear in this problem: the Prandtl number Pr and the Rayleigh number Ra = gfiATH3Pr/v 2. The evolution of the following quantities will be examined u.... =

m a x i m u m of lu(y)l at x =

H/2,

v ..... =

m a x i m u m of Iv(x) l at y =

H/2,

\

~?x /.,-=o

(3)

R.A.W.M. Henkes & C.J. Hoogendoorn

198

Here Nu (Nusselt number) is the non-dimensionalized averaged heat transfer through the hot vertical wall. An oscillation in a quantity q~ can be characterized by its maximum ~bmax, its amplitude 6 ( = ~bmax - qSm~n)and by its frequency f

3. Numerical method The spatial derivatives in the Navier-Stokes equations are discretized with the finite-volume method on a staggered grid. The solution domain is covered with finite volumes, having a grid point for the pressure and the temperature in the centre of each volume, gridpoints for the u-velocity in the middle of the west and east side, and gridpoints for the v-velocity in the middle of the north and south side of each volume. The equations are integrated over each volume, after which mass, momentum and heat fluxes are discretized with finite differences. The finite-volume method has the advantage that it gives a conservative discretization: numerical mass, momentum and heat fluxes can be indicated that are conserved over the domain. All spatial derivatives are discretized such that the familiar secondorder accurate central scheme is obtained in the case that the grid is equidistant. The use of the central scheme for the convection, instead of the first-order upwind scheme, is essential in the present stability calculations. The numerical diffusion in the upwind scheme damps oscillations and was often found to give an evolution to a steady state in cases where the central scheme predicted an oscillating final state. In the present calculations, the appearance of thin boundary layers along the vertical walls requires a non-equidistant grid. The u-gridpoints are positioned in the x-direction according to

H

-

/

l

/max

2re

i = 0, l , . . . ,

/max"

(4)

The same spacing is used for the v-grid points in the y-direction. The time dependence is treated implicitly; the spatial derivatives are evaluated at the new time level n ÷ 1, and the time derivatives are approximated with three time levels (B3 scheme) "+l

(0~b) ~-]

=

3~bn+l - 4~b" + q~"-~ 2At + O(At2).

(5)

The discretized system at the new time level n + 1 consists of nonlinear, algebraic equations, which are solved iteratively with a line Gauss-Seidel

Stability of natural convectionflow

199

method. Alternating sweeps are made from the west to east side and from the east to west side of the computational domain. The updating of the iterative solution is such that only tri-diagonal matrices have to be solved at a line. In a Gauss-Seidel sweep the variables u, v and Tare updated one after the other at a line. After each sweep the pressure is updated with the SIMPLE pressure-correction method (Patankar & Spalding [5]). This method determines a pressure correction as long as the continuity equation in a cell is not satisfied: the pressure is increased/decreased as long as there is a net inflow/outflow of mass. In the limit At --+ 0, the equation for the pressure correction p' is nothing but a Poisson equation

632p,

•2p,

~c?x + @2

3p(c3u -

2At

Ov)

~xx + ~yy "

(6)

This pressure-correction equation is solved with a direct solver. Discretization of (6) gives a (symmetric) band matrix for the Laplace operator, having fixed coefficients which only depend on the geometry. The LU decomposition of this matrix has to be determined only once and can be stored. The sweep process at the new time level is stopped when the pressure correction in each cell and the heat flux through the boundaries are below a certain criterion. Typically 5 to 10 sweeps at each time level are sufficient.

4. Stability and bifurcation Although the boundary conditions are steady, an unsteady solution of the Navier-Stokes equations does not necessarily converge to a steady solution at large time. Its large-time behaviour is studied in the qualitative theory of partial differential equations. This theory distinguishes the following solutions to which an initial solution can be attracted for large time: (i) steady attractor, (ii) periodic attractor, (iii) quasi-periodic attractor, (iv) strange or chaotic attractor. All these attractors, except the first, show an essentially unsteady behaviour at large time. The attractors (ii) and (iii) give a solution with only discrete frequencies (spikes) in the Fourier-spectrum. The periodic attractor (ii) gives a single-frequency oscillation (or so called limit cycle). The quasi-periodic attractor (iii) consists of several incommensurate frequencies. The solution of the chaotic attractor is non-periodic, showing a broadband spectrum. The attractor might depend on the initial solution at t = 0; multiple steady or

200

R.A.W.M. Henkes & C.J. Hoogendoorn

unsteady solutions of the unsteady Navier-Stokes equations at t ---, oo can exist. Therefore, variation of the Rayleigh number can give different branches of solutions. These branches can be either stable or unstable against very small disturbances. The phenomenon that a solution remains on its branch as long as disturbances are small enough (linear stability), and only jumps to another branch for larger disturbances (nonlinear instability), is referred to as hysteresis. We can study the large-time behaviour of the Navier-Stokes equations by increasing the Rayleigh number for a fixed Prandtl number. The solution can proved to be steady, and unique, up to a critical Rayleigh number. Above Racr the steady solution becomes unstable, i.e. disturbances no longer decay for increasing time. Linear stability theory predicts an exponential growth of the unstable disturbances. However, during the growth of the disturbances nonlinear effects become important and the solution ends up in a new solution that is stable against small disturbances. The Rayleigh number at which the solution branches to another one (Racr) is a bifurcation point of the equations. Different types of bifurcation exist; two well-known types are the pitchfork bifurcation (branching from one steady to another steady solution) and the Hopf bifurcation (branching from a steady solution to an unsteady solution with a single-frequency oscillation). The Hopf bifurcation is related to the periodic attractor (ii). Winters [1] explains how the bifurcation in a steady solution can numerically be determined. The unsteady Navier-Stokes equations are linearized around a steady discrete solution. After discretization a system of linear algebraic equations results:

+ J¢

=

0.

(7)

J is the Jacobian matrix of the steady part. The steady solution ~bis linearly stable if all eigenvalues o- of J have Re(o-) > 0. If Re(a) = Im(o-) = 0 (i.e. d e t ( J ) = 0), the steady solution shows the pitchfork bifurcation. If Re(a) = 0 and Im(o-) = _+ ie), the steady solution shows a Hopf bifurcation, giving the single frequency f = co/(2rc). Winters' method has the advantage that it can locate the position of the Hopt bifurcation without the calculation of the unsteady oscillating solution. If, on the other hand, one wants to determine the bifurcation by the calculation of the unsteady solution, the position of the Hopf bifurcation can be extrapolated from the supercritical oscillating solution with the help of the series expansion of

Stability of natural convectionflow

201

Joseph & Sattinger [6] 6 + x/Ra - Racr f-for

(8)

+ Ra -- Racr.

Winters' method of examining the eigenvalues to determine a bifurcation can also be generalized to determine the linear stabililty of a singlefiequency oscillation; in that case the Floquet multipliers of the system have to be determined (see Guckenheimer & Holmes [7]). Determining the stability of less simple unsteady solutions by examining the eigenvalues is complicated. The use of the unsteady solution has the advantage that besides the Hopf bifurcation, it can also determine the bifurcation and (nonlinear) stability from any other steady or unsteady solution. In the present study the stability is examined by the calculation of the unsteady Navier-Stokes solution. Recently (Roux [8]) benchmark numerical results were obtained for the Hopf bifurcation in the Navier-Stokes solution at low Prandtl numbers in a shallow cavity with conducting horizontal walls and differentially heated vertical walls. Our computational code gave good results for this benchmark case (see Henkes & Hoogendoorn [9]). The Hopf bifurcation introduces an unsteady solution. Because its spectrum is discrete, we cannot denote this as a turbulent solution. However, the occurrence of a Hopf bifurcation can be seen as the initial phase of the transition to turbulence. The total path to turbulence seems to consist of different transitions between the attractors (ii), (iii) and (iv). If the Rayleigh number is increased far enough, the chaotic attractor is found, with a broadband spectrum that is characteristic for the turbulence.

5. Instability mechanisms The present study investigates the stability for air (Pr = 0.71) and water (Pr = 7.0) in a square cavity that is heated from the side. Both the cases of conducting and adiabatic horizontal walls are considered. For increasing Rayleigh number the steady solution shows boundary layers along the vertical and horizontal walls and a stagnant core with horizontal streamlines. Further, the core is thermally stratified (horizontal isotherms). Gill [10] has suggested that the vertical velocity in the vertical boundary layers scales with x~ff~ATH, whereas the boundary-layer thickness scales with H U a -1/4. The horizontal velocity in the core scales with

202

R . A . W . M . Henkes & C.J. Hoogendoorn

(gflATv)l/3Ra-~/12, and the characteristic lengthscale in the core is H. We found good agreement of these scalings when compared with steady solutions of the Navier-Stokes equations. With the help of these spatial scalings Patterson & Imberger [11] derived timescales (i.e. the ratio of the characteristic length and velocity) which are expected to be important in an unsteady evolution. The vertical lengthscale H and the vertical velocity scale (gflATH) 1/2 in the vertical boundary layer gives the timescale (HZ/v)Ra -1/2. The horizontal length scale H and the horizontal velocity scale (g~ATv) ~/3Ra-1/~2 in the core gives the timescale (H2/v)Ra - 1/4. For large Rayleigh numbers, the timescale in the core is largest. The present calculations confirm that the timescale in the core is indeed characteristic for the time required for the initial solution to reach its final state (see also Yewell et al. [12], Hyun [13]). The molecular timescale ts = (H2/v)Ra ° only enters the problem at moderate Rayleigh numbers, but not at the large Rayleigh numbers as considered here. The distinction of the steady solution in a stratified core region, horizontal boundary layers and vertical boundary layers gives at least three possible instability mechanisms: (i) a Rayleigh/B6nard instability in the core, (ii) an instability after the hydraulic-jump in the horizontal boundary layer and (iii) a Tollmien-Schlichting instability in the vertical boundary layer. 5.2. Rayleigh/BOnard instability

For large Rayleigh numbers the stratified core can be considered as inviscid. Stability analysis (see for example Drazin & Reid [14]) on the inviscid equations for the core shows that the stratification is stable if the density decreases with the height (or, if/3 is positive: if the temperature increases with the height). Moreover the analysis shows that if the steady solution in a stable stratification is disturbed, the flow returns to its original state in an oscillatory way, showing damped internal gravity waves with the BruntVfiisfil~ frequency ,,/ff~ATH/H as the characteristic frequency. Viscosity can delay the instability in the unstably stratified environment until a critical Rayleigh number is exceeded. Rayleigh found a pitchfork bifurcation at R a , = 1708 in the classical problem of an infinite horizontal layer that is heated from below. Here the Rayleigh number is based on the thickness of the layer and the temperature difference across the layer. The horizontal walls are rigid. Below Racr there is a linear thermal stratification and zero convection. Above R a , there is convection, showing the well-known B6nard cells.

Stability oJ natural convection,/tow

203

5.2. Instability after the hydraulic jump The hydraulic jump is the phenomenon that a thin high-speed horizontal layer almost discontinuously jumps to a broader low-speed horizontal layer. The hydraulic jump can occur in the corner of the cavity, where the vertical boundary layer hits the horizontal wall and bends to a horizontal layer. If the jump is too strong the flow after the jump becomes unstable and the freecoming jump energy is dissipated via an unsteady oscillating wave pattern. Ivey [15] has observed such a hydraulic jump during the time evolution to a steady state of water in the heated cavity at high Rayleigh numbers.

5.3. Tollmien-Schlichting instability The Tollmien-Schlichting instability of the natural-convection boundary layer along a heated, semi-infinite, vertical plate in an isothermal environment has been studied theoretically by Nachtsheim [16]. He numerically solved the Orr-Sommerfeld equations. These equations approximately describe the evolution of disturbances on the steady solution of the boundarylayer equations. The approximations are: linear disturbances, locally parallel flow (i.e. zero normal velocity) and certain assumptions about the shape of the disturbances (both in space and time). Some studies have made a similar stability analysis, but without the parallel-flow assumption (Haaland & Sparrow [17], Tzuoo et al. [18], Lee et al. [19]). The calculated critical Rayleigh number and resulting frequency and wavelength are summarized in Table 1. The Rayleigh number is based here on the y-coordinate along the vertical plate and on AT = 2(Th - T~), where T~ is the isothermal environment temperature. The factor 2 is introduced here to enable us to make a senseful comparison with the cavity geometry; replacing T~ by the centre temperature in the cavity (Th -T~.)/2 gives the earlier definition

AT=Th-T

c.

Because of the horizontal walls, the boundary layer along the vertical walls in the cavity is different from the one along a heated vertical plate in an isothermal environment. Henkes et al. [20] showed that the main influence of these horizontal walls on the development of the vertical boundary layer is via the thermal stratification in the core. The linear stability analysis of Gill & Davey [21] takes this stratification into account. It determines the stability against small disturbances of the boundary layer in an infinite vertical slot having a linear stably stratified core. The main difference with the present square cavity situation is that the steady basic flow of Gill & Davey is parallel (and therefore does not start to develop at the bottom), and that their wall temperature is not isothermal, but increases with the height in the

204

R.A.W.M. Henkes & C.J. Hoogendoorn

same linear way as the core temperature. The core stratification is described by the parameter ~ ( = temperature gradient nondimensionalized with AT and H). Their results in Table 1 show that Rac~ largely depends on the stratification parameter 3. Jaluria & Gebhart [22], Gebhart & Mahajan [23] and Gebhart [24] have experimentally studied the transition of the natural-convection boundary layer for water along a vertical plate with a constant wall-heat flux placed in an isothermal environment. The first unsteadiness in the experiments shows a single-frequency oscillation. Moreover the unsteadiness for larger Ra~. values still shows a single oscillation, in which the frequency is independent of y. This y-independence implies that the non-dimensional frequency fvl/~/(g~AT) 2/3 is independent of Ra~.. Gebhart & Mahajan [23] indicate this non-dimensional frequency as the "characteristic frequency" of the boundary layer. The Orr-Sommerfeld stability analyses do not find a single frequency for supercritical Ra~. values, but a continuous spectrum of unstable modes. Of course, one has to realize that the Orr-Sommerfeld equations only approximately describe the linear stability. However, Gebhart & Mahajan note that the frequency which is most strongly amplified in the Orr-Sommerfeld solution for supercritical Ra~. values closely follows a curve in whichfv~/3/(g~AT) 2/3 is independent of Ra.~.. This is in particular true for larger Ra, values, whereas the deviation for smaller values is a bit larger. As indicated in Table 1, the theoretical frequencyfv~/3/(g~AT) 2/3 for water, as found by Gebhart & Mahajan along this way, is close to the experimental result of Jaturia & Gebhart [22]. The agreement for air between this theoretical value and the experimental value of Eckert & Soehnghen [25] is a bit worse. If Ra). is further increased the experiments show that more frequencies enter the flow and the transition to the turbulent state is found. The precise transition is a complicated process. Jaluria & Gebhart [22] emphasize that the whole transition process cannot be described as a function of Ra~ (for Pr fixed) alone. In mathematical terms, they suggest that in the transition regime several (linearly) stable solutions exist; whether the experiments give the one or the other solution depends on the precise experimental situation. This can explain the large difference between the critical Rayleigh numbers in Table 1. Although the fully turbulent natural-convection boundary layer flow shows a broadband spectrum of frequencies, Kitamura et al. [26] experimentally observed coherent structures in it. The main characteristic of this coherent structure is a motion of large eddies, with a predominant frequency that can be fitted t o fvL/3/(g~AT) 2/3 = 0.0063 (see Table 1). Here we have given a short review of both experiments and theoretical stability analyses to determine the first unsteadiness in the natural-convection

Stability

of natural

convection f l o w

205

Table 1. Comparison of existing theoretical and experimental stability results for the heated vertical plate (also added are the present results in the cavity, see Section 7); critical Rayleigh number, frequency (f) and wavlength (2) (a) air (Pr = 0.71)

Study

f vU3

}~(gflA T )l/3

(g/~AT)2/3

v2/3

3.82 x 105 1.73 x 106 1.13 x 108 z3.

0.0162 0.0263 0.0351 0.0239

211 156

6 x 108 -

0.0415 0.0063

2 x 108*

0.0217

(R%,)or

fv t/3

(Ra),)c~

Theoretical: Nachtsheim (1963) Lee et al. (1987) Gebhart & Mahajan (1975) Gill & Davey (1969)

161

Experimental: Echert & Soehnghen (1951) Kitamura et al. (1985)

100

Present: adiabatic cavity (b) water (Pr = 7.0) Study

(g/~AT)2/3

2(g~AT) l/~ ~2/3

Theoretical: Nachtsheim (1963) Lee et al. (1987) Gebhart & Mahajan (1975) Gill & Davey (1969)

2.82 x 105 1.12 x 106 2.21 x 108 z3.

0.0228 0.0293 0.0251 0.0234

57 53

2 x 109 -

0.0248 0.0063

100

4 x 109*

0.0403

41

56

Experimental: Jaluria & Gebhart (1973) Kitamura et al. (1985)

Present: adiabatic cavity *y = H b o u n d a r y l a y e r a l o n g a h e a t e d v e r t i c a l p l a t e in a n i s o t h e r m a l o r s t r a t i f i e d e n v i r o n m e n t . F r o m T a b l e 1, w e c a n c o n c l u d e t h a t b o t h s t a b i l i t y a n a l y s e s a n d e x p e r i m e n t s m o r e o r less p r e d i c t t h e s a m e frequencyfvl/3/(g~AT)2/3 f o r t h e first u n s t e a d i n e s s . T h i s is in p a r t i c u l a r t r u e f o r w a t e r , b u t t h e a g r e e m e n t for air is less g o o d . A l l t h e t h e o r e t i c a l s t a b i l i t y a n a l y s e s m e n t i o n e d in this s e c t i o n h a v e in common that they are simplified linear stability analyses for simple conf i g u r a t i o n s . It is n o t b e f o r e h a n d c l e a r h o w these r e s u l t s r e l a t e to t h e s t a b i l i t y o f t h e flow in t h e s q u a r e c a v i t y h e a t e d f r o m t h e side. T h e r e f o r e in t h e n e x t t w o s e c t i o n s we give t h e s t a b i l i t y r e s u l t s f o r the c a v i t y flow, as c a l c u l a t e d from the unsteady, two-dimensional Navier-Stokes equations.

206

R.A.W.M. Henkes & C.J. Hoogendoorn

6. Conducting horizontal walls Firstly we reconsider the stability for air. Le Qu6r6 & Alziary de Roquefort [2, 3] already investigated this problem by a similar approach as used in the present study, namely by solving the unsteady, two-dimensional NavierStokes equations. Increasing the Rayleigh number from 106 up t o 107, they find transitions (showing hysteresis) from the steady solution to three unsteady oscillating branches, with the non-dimensionalized frequencies fH/x~flATH = 0.255, 0.29 and 0.32 respectively. The frequency of the oscillation for supercritical Rayleigh numbers is almost independent of the Rayleigh number if it is non-dimensionalized with the Brunt-V/iis/il/i freThe intermediate frequency corresponds to a centroquency g f l ~ / H . symmetry breaking solution. These frequencies were also experimentally found by Briggs & Jones [27]. Winters [1] detected five H o p f bifurcations in the steady solution, increasing the Rayleigh number from 2 x 106 up to 3 x 106. His non-dimensionalized frequencies are respectively 0.254, 0.215, 0.290, 0.173 and 0.325. The three highest frequencies seem to correspond to those found in the unsteady calculations and experiments. As summarized in Table 2, our results on the fine grid (60 x 60 space, AtxfgflATH/H = 1/8) agree with the values in the above-mentioned literature. Also the results on the coarse grid (30 x 30, Atx/-g~ATH/H = 1/8) do reasonably predict the bifurcation. Calculations on the fine grid were made in the Rayleigh number range 1 0 6 - - 3 x 106. The solution at Ra = 106 shows an evolution to the steady solution. The solution at Ra = 2.25 x 106 is attracted to an oscillating solution; all the characteristic quantities um~, 7,)max and Nu do oscillate with the same frequency. The oscillating solution at Ra = 3 x l06 has approximately the same non-dimensionalized frequency. Therefore, hysteresis effects seem to prevent jumps to the other unsteady oscillating branches in the present calculations. Extrapolating the amplitude of the oscillation of the solutions in the range Ra = 2 x 1 0 6 Table 2. Calculated frequency for air in the case of conducting horizontal walls (values at our coarse grid in parentheses)

Ra

Study

Type

f

106 2.25 x 106 3 x 106 5 x 108 2.2 x 106 2.1092 x 106 3.2 × 10 6

present present present present Le Qu~r~ & AdR (1986) Winters (1987) experimental (Briggs & Jones 1985)

steady oscill. oscill. oscill. oscill. bifurcation oscill.

0.248 (0.295) 0.253 (0.300) 0.238 0.255 0.254 0.248

Stability of natural convection flow

207

3 x 10 6 t o zero with eq. (8) gives Racr = 2.07 x 10 - 6 . Winters finds Racr = 2.11 x 10 6. The steamlines and isotherms at different Rayleigh numbers are shown in Fig. 2. We also calculated the much larger Rayleigh number 5 x 108 on the 60 x 60 spatial grid with At,,/-gflATH/H = 1/8. The solution is no longer periodic, but chaotic. However, the dominant frequency in the spectrum is fH/~/--g-~ATH = 0.238. This value is close to the frequency fH/,,J-g~ATH = 0.248 as found for Ra = 2.25 x 106. Secondly, we considered the stability for water with Pr = 7.0. The solution at Ra = 4 x 10 6 still converges to a steady state. An oscillation in all characteristic quantities is found at Ra = 5 x 10 6. The oscillation has a frequency fH/x~g-fiATH = 0.157 at the fine grid (60 x 60, Atx/-gfiATH/H = 1/8) and 0.146 on the coarser grid (30 x 30, Atx/gfiATH/H = 1/8). The evolution of Umaxis given in Fig. 3. The limit cycle at Ra = 6 × l06 is shown in Fig. 4. The cycle shows the clockwise circulation of eddies in the core, moving just along the horizontal and vertical boundary layers. Comparison of the stability for air and water in the configuration with conducting horizontal walls shows that an increase of the Prandtl number increases the critical Rayleigh number. In the case of adiabatic horizontal walls (see the next section) the instability is delayed until Ra~r = 2 × 108 for air. In order to interpret the nature of the instability in the case of conducting horizontal walls we have compared the streamlines and isotherms at Ra = 106 (just below Ra~r) as shown in Fig. 2a with the solution for adiabatic horizontal walls at the same Rayleigh number. Both solutions show boundary layers along the vertical walls and a stably stratified core. In those regions the solutions only slightly deviate. For example, the velocity maximum at half the cavity height v .... /x~fiATH is 0.304 in the conducting case and 0.263 in the adiabatic case, and (H/AT) c~T/?y in the cavity centre is 0.71 in the conducting case and 0.92 in the adiabatic case. The main difference occurs along the horizontal walls; in the conducting case the temperature decreases with the height in the layer along the horizontal walls. In Fig. 2a the broken line (with 3T/c~y = 0) indicates up to where this unstable layer extends. In the adiabatic case this line follows the horizontal wall. As explained in the previous section, viscosity can stabilize this unstable temperature up to a critical Rayleigh number. The maximum local Rayleigh number at the broken line (based on the local temperature difference with the T(x) temperature at the floor and based on the local height) is 6.1 x 10 3 for air (6.0 x 10 3 for water). This maximum is indicated with a dot in the figure. This value is in order of magnitude comparable to 1.7 × 10 3 found for the infinite horizontal layer heated from below. Therefore we can relate the instability in the cavity with the Rayleigh/ B6nard instability. For supercritical Rayleigh numbers the frequency

208

R.A.W.M. Henkes & C.J. Hoogendoorn

STREAMLINES

ISOTHERMS

(a)

(b)

\ (c) Fig. 2. Streamlines a n d i s o t h e r m s for air in the case o f c o n d u c t i n g horizontal walls; (a) steady solution at Ra = 106, (b) s n a p - s h o t f r o m the limit cycle at R a = 3 x 106, (c) time-averaged solution at R a = 5 x 108.

Stability of natural convection flow

209

0.045 ! Ra = 5 x 106

0.035 0.045

I

Ra = 4× 106

tv/H 2

0.2

Fig. 3. T i m e e v o l u t i o n of Urn.x for w a t e r in the case of c o n d u c t i n g h o r i z o n t a l walls.

in the cavity scales with the Brunt-V/iis/ilfi frequency. This indicates that, although the instability is introduced in the viscous horizontal layer, the resulting frequency is determined by the wave mode of the inviscid core.

7. Adiabatic horizontal walls

The stability for air in cavities with adiabatic horizontal walls and an aspect ratio ( = height-width ratio A) in the range 4-10 has been calculated by Le Qu6r6 & Alziary de Roquefort [4]. They calculated that the critical Rayleigh number (based on the cavity height) is almost independent of the aspect ratio, namely Rac~ ~ 108. The frequency increases from fH/x/-g~ATH = 0.478 at A = 4 to 0.847 at A = 10. Further, the frequency for supercritical Rayleigh numbers in a cavity with a given aspect ratio is almost independent of the Rayleigh number if it is scaled with the Brunt-Vfiisfilfi frequency. Le Qu~r6 & Penot [28] have verified the calculated stability for air in the A = 4 cavity experimentally. The boundary conditions along the horizontal walls in the experiment were not clearly known, but were expected to be somewhere between the conducting and adiabatic case. They measured Racr = 1.1 x 108 with fH/x/-g~ATH = 0.522, which reasonably agrees with the calculations (Rac~ = 4.5 x 107 withfH/x/-g~ATH = 0.578 for the conducting case, and Ra~ = 1.1 x 10~ withfH/x/-g~ATH = 0.478 for the adiabatic case).

210

R.A.W.M. Henkes & C.J. Hoogendoorn

i

Fig. 4. Limit cycle for water at Ra = 6 x 1 0 6 in the case of conducting horizontal walls. (8 equidistant time intervals; the cycle starts in the right upper corner where Ureax reaches its maximum.)

Stability of naturaI convectionflow

211

-.d ",d

(a)

(b)

(c)

(d)

Fig. 5. Streamlines a n d grid, Atgflx/TfiAfH/H = streamlines a t R a

i s o t h e r m s for air in the case of adiabatic horizontal walls (80 x 80 1/8); (a) streamlines at R a = 108, (b) i s o t h e r m s at Ra = 10~, (c) = 5 × 105, ( d ) streamlines at R a = 7 x 10L

The stability for air in the square cavity (A = 1) with adiabatic horizontal walls has been calculated in the present study. The flow at Ra = l0 s converges to a steady state. As illustrated in Fig. 5a/b, the flow shows thin boundary layers along the vertical walls, and an almost stagnant, thermally stratified core. In line with the theory of Patterson & Imberger [11], this high-Rayleigh number flow reaches its steady state in an oscillatory way. These damped oscillations are internal gravity waves in the core, having a frequency fH/x[-g-flATH = 0.107 (almost independent of the Rayleigh number). When the Rayleigh number is increased above l0 s, the flow no longer reaches a steady state at Racr = 2 x 10s. At this Rayleigh number the flow becomes quasi-periodic with two incommensurate frequencies: it has a lower frequency fH/.,~flATH = 0.038 and a higher frequency fH/,,/-gflATH = 0.555. These frequencies are almost independent of the Rayleigh number if they are non-dimensionalized with the Brunt-VS, is/ilfi

212

R.A.W.M. Henkes & C.J. Hoogendoorn

frequency (see also Fig. 7). A typical time evolution of the quantities Nu, vm,x and urn,x at a supercritical Rayleigh number is shown in Fig. 6. Nu and vm,Xshow both the lower and higher frequency, but Um~x(at x = 1t/2) only shows the lower frequency. Because Nu and vm,~are quantities of the vertical boundary layer and u~,x is not, this implies that the higher frequency is related to an instability in the vertical boundary layer. This is also clear from the streamline patterns in Fig. 5c/d; there is a short-wave instability at the edge of the vertical boundary layers. The lower frequency can be identified as the long-wave in the core. Because the left upper corner shows a hydraulic jump, with a recirculation region, the lower-frequency is expected to be related to an instability after this jump. In order to verify the accuracy of the calculations we have used different grids: 40 × 40, 60 × 60 and 80 × 80 gridpoints. The timestep was refined up to At,fg-~AT/H -- 1/16, which gives almost timestep-independent results. The 1/16 timestep gives approximately 421 timesteps in each lower-frequency oscillation and 29 timesteps in each higher-frequency oscillation. The results on the different spatial grids are given in Fig. 7. The lower frequency in Fig. 7a is already accurate at the coarse 40 × 40 grid. The higher frequency in Fig. 7b shows a slightly larger grid-dependence. The higher frequency on the 80 × 80 grid splits up into more branches at about Ra = 4 × 108; a new bifurcation seems to be found at this Rayleigh number. The upper branch gives the frequency for vm,~ and the lower branch gives the frequency for Nu. The value of Raor and the amplitude of the oscillation show a large grid-dependence. The amplitude in Vm,~ (Fig. 7c) is almost totally due to the higher-frequency oscillation. The 40 × 40 grid is much too coarse. The higher-frequency in ~Umax at the 60 × 60 grid appears at Racr = 3 × 108, but disappears between Ra = 4 × 108 and 6.25 × 108. At the 80 × 80 grid this frequency appears at Ra = 2 × 108. Its amplitude is very small up to Ra = 4 × 108, after which it considerably increases. We do not know of any stability calculations for water (Pr = 7.0) in the heated cavity with adiabatic horizontal walls. Unsteady experiments in the square cavity for water up to Ra = 1.2 × 109 are made by Ivey [15]. For these Rayleigh numbers he finds an evolution to a steady state. Ivey starts with a flow in rest and suddenly switches on the temperature difference over the vertical walls at t -- 0. The initial phase of the evolution shows (damped) oscillations. From the flow visualization it is clear that these oscillations are initiated in the left upper corner of the cavity, where the hot vertical boundary layer turns to a horizontal layer; during the evolution a hydraulic jump seems to occur in the corner. The released energy is dissipated in a (damped) oscillation. It is emphasized that the hydraulic jump for water at

Stability of natural convection flow

213

X

:gg e ¢,t

"~ II ""~

>

o

m

~a o < I D m

c5

0

Hzvr3X/v

0

-C}

c5

c5 S ~

Y

S

o

c5

o

o

.Hev~3~

¢-1

214

R.A.W.M. Henkes & C.J. Hoogendoorn o.lo I

7

r

t

f

g

GRID 80x80 60x60 40x40

• o + 0.05 .

10

- - -4-~ - ~ _ ~ - - -

~_~-~--~-

I

I

2

3

I

1

4 (a) i

o- -o-

I

I

I

5

6

7 X 10 . 8

Ra

I

i

i

_o

• ° +

I

•

o +

GRID 80x80 60x60 40x40

I

I

I

I

I

3

4

5

6

7

Ra x 10 -8

(b) 0.04

I

I

J

I - -

I

I

+~+

80x80 60x60 40x40

/

//i

0.02

I

/

'

3

4 (c)

5

6

7

Ra × 10 -8

Fig. 7. Grid dependence o f oscillating Vmax for air in the case o f adiabatic horizontal walls (Atg~/H = 1/16); (a) lower frequency, (b) higher frequency, (c) amplitude.

Ra = 1.2 x 109 Occurs only in the initial phase of the evolution, but is absent in the steady final state. We have calculated the stability for water. In contrast with Ivey's experiments, we initiated the time evolution for a certain Rayleigh number with

Stability of natural convection flow

215

(a)

I ------4

~

~

i

(b)

Fig. 8. Streamlines (a) and isotherms (b) for water in the case of adiabatic horizontal walls (Ra = 3 x 109, 120 x 120, Atg~x/~--T-H/H = 1/32; the right picture shows the hot vertical boundary layer a 20 times enlarged in the x direction).

the solution belonging to a slightly different Rayleigh number. Below Racr "~ 4 x 109 the flow converges to a steady state in an oscillatory way. The damped oscillations are internal gravity waves in the core with a frequency fH/x/-gflATH = 0.088. The numerical result that Racr is above 1.2 x 109 agrees with Ivey's experiments, who finds a steady solution up to at least that value. For supercritical Rayleigh numbers a periodic solution is found with the frequencyfH/x/-g-fiATH ~ 1.16. This frequency has twice the value for air. Its value is practically independent of the Rayleigh number if it is non-dimensionalized with the Brunt-V/iisfil/i frequency. This frequency is only visible in Vmax,but not in Nu or ureax. In contrast with air, the flow for water in Fig. 8a does not show a hydraulic jump in the corner; in the range of Rayleigh numbers considered, water does not show the lower-frequency as found for air. As can be seen in the same figure, the instability for water is related to a traveling wave in the vertical boundary layer.

216

R.A.W.M. Henkes & C.J. Hoogendoorn 1.5 y-grid: • non-equidistant o equidistant \

o\ 0.5

oo

I

O'025Ay/H at y=H/2 0.05

Fig. 9. Grid dependence of the frequency for water in the case of adiabatic horizontal walls. We investigated the accuracy by refinement of the timestep and the spatial grid. The timstep was refined up to At,fgflATH/H = 1/32, which gives almost timestep-independent results. The spatial grid was refined from 40 x 40 up to 160 x 160 points. Besides the non-equidistant grid in the y-direction (eq. (4)) we also used an equidistant grid in the y-direction. Figure 9 shows the frequency in Vm,xat half the cavity height for supercritical Rayleigh numbers. At half the cavity height the non-equidistant n x n grid gives a meshsize Ay = 2H/n, and the equidistant grid gives Ay = H/n. This Ay value is plotted on the horizontal axis in Fig. 8. The resulting frequency turns out to depend on this local meshsize, almost irrespective of the use of an equidistant or non-equidistant grid; the equidistant y-grid with n x n points predicts the frequency in Vm,~ as accurate as the non-equidistant 2n x 2n grid. However, even on the finest grids there still is some griddependence. If we assume that the spatial discretization leads to a secondorder global error in the solution on sufficiently fine grids, a Richardson extrapolation can be used to extrapolate the calculated frequencies at the two finest grids to fH/x/-g-flATH = 1.16 at a zero-grid. Besides the frequency, also Ra~r for water largely depends on the grid. It increases from R a , = 3 x 108 at the 40 x 40 non-equidistant grid to 4 x 109 at the 160 x 160 non-equidistant grid. The critical Rayleigh number at the 120 x 120 equidistant grid is 5 x 109. A remarkable difference between the results on the non-equidistant and equidistant grids is shown in Fig. 10: the non-equidistant grid gives a regular single oscillation, and the equidistant grid gives an intermittent solution. The frequency in both types of solutions is almost the same. Further grid refinement might exclude one of the types. But, the fact that all non-equidistant grids give a periodic solution and all the equidistant grids give a intermittent solution might suggest that above Racr two stable attractors can occur.

Stability of natural convection flow

217

0.105

0.075[ 0

0.088

0

tv/H2

(a)

tvlH 2

(b)

0.005

0.005

Fig. tO. Time evolution of vm~x for water in the case of adiabatic horizontal walls; (a) Ra = 4 × 109 on a non-equidistant y-grid (160 x 160, Atx/g[IATH/H = 1/32), (b) Ra = 5 x 109 on an equidistant y-grid (120 x 120, Atx/g~ATH/H

=

1/32).

It turns out that m u c h m o r e gridpoints are required for water t h a n for air. The n u m b e r o f gridpoints needed is determined by the critical wavlength; this wavelength is p r o b a b l y larger for water than for air. F r o m Fig. 8 we can derive that the traveling wave in the b o u n d a r y layer for water a p p r o x i m a t e l y has a wavelength 2/I-I = 0.05 at h a l f the cavity height. I f at least 6 gridpoints are used to describe such a wavelength, at least a 120 x 120 (equidistant) grid is required. This wave pattern was not clearly visible in the oscillating solutions for air, a n d the critical wavelength for air could not be derived f r o m it.

218

R.A.W.M. Henkes & C.J. Hoogendoorn

To summarize, above a critical Rayleigh number two instabilities are found for air and one for water. One instability for air seems to be related to the instability after the hydraulic-jump. The other instability for air and the instability for water is an instability in the vertical boundary layer. This boundary-layer instability relates to the Tollmien-Schlichting instability in the vertical boundary layer along the heated plate as described in Section 5. The calculated frequencies are added to Table 1, which also compares existing stability results for the heated vertical plate. This table shows a good agreement of the frequency for air, whereas the frequency for water is higher. A remarkable difference between the oscillations in the boundary layer along the heated cavity wall with the stratified core-environment and the oscillations in the boundary layer along the heated vertical plate in the isothermal environment (Gebhart & Mahajan [23]) is that the cavity oscillation scales with x/f/3A TH/H and the plate oscillation scales with (gilA T)2/3/ v1/3.The former is an inviscid scale, whereas the latter is viscous. A possible explanation of the inviscid scaling of the boundary-layer frequency is that although it seems to be initiated as an unstable Tollmien-Schlichting wave as for the plate configuration, the frequency at supercritical Rayleigh numbers results from interaction with the inviscid core.

8. Conclusion

By numerically solving the unsteady, two-dimensional Navier-Stokes equations, we have calculated the stability of the steady laminar naturalconvection flow of air (Pr = 0.71) and water (Pr = 7.0) in a square cavity that is differentially heated over the vertical walls. For conducting horizontal walls, the flow of air looses stability at Racr = 2 x 10 6 and water at 5 x 10 6. At supercritical Rayleigh numbers the unsteady flow is periodic with a frequencyfH/~/g~ATH = 0.248 for air and 0.157 for water. This frequency is almost independent of the Rayleigh number if it is scaled with the Brunt-Vfiisfil/i frequency x~flATH/H. The steady solution at subcritical Rayleigh numbers shows an unstable temperature profile (it decreases with the height) close to the horizontal walls. This indicates that the instability is related to the Rayleigh/Bdnard instability. The Brunt-V/iis/il/i scaling of the frequency refers to an interaction with the inviscid core. For adiabatic horizontal walls, the flow of air looses stability at Racr = 2 x 108. Below this value the evolution to the steady state shows damped internal gravity waves in the core with a frequencyfH/x~ATH = 0.107. Beyond Racr, the unsteady final solution shows two incommensurate frequencies,fH/x/g~ATH = 0.038 and 0.555. The lower frequency appears

Stability of natural convectionflow

219

most clearly in the core. Because the streamline pattern shows a hydraulic jump in the corner, where the vertical boundary layer is bent to a horizontal layer, this lower frequency seems to be related to an instability after the hydraulic-jump. The higher frequency appears in the natural-convection boundary layer along the vertical walls. Therefore, this instability seems to be related to the Tollmien-Schlichting instability, which is found in the boundary layer along a semi-infinite hot plate in an isothermal environment. However, the frequency for the plate is known to scale with (gflAT)2/3/v 1/3, whereas the cavity frequency scales with x/gflATH/H. This indicates that the cavity frequency results from interaction with the inviscid, stratified core. Below Racr ~ 4 x 109, the solution for water reaches a steady state after internal gravity waves with a frequency JH/x/gfiATH -- 0.088 have been damped. The solution beyond this Rayleigh number is periodic with a frequency fH/x/g-flATH ,,~ 1.16. The instability appears in the vertical boundary layer and seems to be related to the Tollmien-Schlichting instability. In contrast to air, water does not show the hydraulic jump and no lower frequency is found. The stability calculations for water in the case with adiabatic horizontal walls requires much more gridpoints than for air: at least 120 x 120 points are needed for water. In particular, the solution for water shows a periodic solution on a non-equidistant y-grid, whereas the oscillation is intermittent on an equidistant y-grid. Further grid refinement is required to check whether the different behaviour is physical or only of numerical nature. References 1. K.H. Winters, Hopf bifurcation in the double-glazing problem with conducting boundaries, J. Heat Transfer 109 (1987) 894-898. 2. P. Le Qu6r~ and T, Alziary de Roquefort, Transition to unsteady natural convection of air in vertical differentially heated cavities: influence of thermal boundary conditions on the horizontal walls. Proc. 8th Int. Heat Transfer Conf. San Francisco (1986) 1533-1538. 3. P. Le Qu~r+ and T. Alziary de Roquefort, On the existence of multiple periodic solutions of the Boussinesq equations. Mechanics of Fluids, C.R. Acad. Sci. Paris, 306, 1I (1988) 681 687. 4. P. Le Qu6r6 and T. Alziary de Roquefort, Transition to unsteady natural convection of air in differentially heated vertical cavities. 4th lnt. Conf. on Numerical Methods in Laminar and Turbulent Flow, Pineridge Press (1985) 841-852. 5. S.V. Patankar and D.B. Spalding, A calculation procedure for heat, mass and momentum transfer in three dimensional parabolic flows. Int. J. Heat Mass Tran,~'fbr 15 (1972) 1787-1806. 6. D.D. Joseph and D.H. Sattinger, Bifurcating time periodic solutions and their stability. Arch. Rational Mech. Anal. 45 (1972) 79-109. 7. J. Guckenheimer and P. Holmes, Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Applied Mathematical Sciences 42, Springer (1983). 8. B. Roux (ed), Proc. of the Gamin-workshop "Numerical simulation of oscillatory convection in low Pr fluids". To appear in the series '"Notes on Numerical Fluid Mechanics", Vieweg (1989).

220

R.A.W.M. Henkes & C.J. Hoogendoorn

9. R.A.W.M. Henkes and C.J. Hoogendoorn, Towards a benchmark solution for oscillating convection: Contribution of the heat-transfer group at Delft University. To appear in [8] (1989). 10. A.E. Gill, The boundary-layer regime for convection in a rectangular cavity. J. Fluid Mech. 26 (1966) 515-536. 11. J. Patterson and J. Imberger, Unsteady natural convection in a rectangular cavity. J. Fluid Mech. 100 (1980) 65-86. 12. R. Yewell, D. Poulikakos and A. Bejan, Transient natural convection experiments in shallow enclosures. J. Heat Tran,fer 104 (1982) 533-538. 13. J.M. Hyun, Transient buoyant convection of a contained fluid driven by the changes in the boundary temperatures. J. Appl. Mech. 52 (1985) 193-198. 14. P.G. Drazin and W.H. Reid, Hydrodynamic stability. Cambridge University Press (1981). 15. G.N. Ivey, Experiments on transient natural convection in a cavity. J. Fluid Mech. 144 (1984) 389-401. 16. P.R. Nachtsheim, Stability of free-convection boundary-layer flows. N A S A TN D-2089 (1963). 17. S.E. Haaland and E.M. Sparrow, Wave instability of natural convection on inclined surfaces accounting for nonparallelism of the.basic flow. J. Heat Transfer 96 (1973) 405-407. 18. K.L. Tzuoo, T.S. Chen and B.F. Armaly, Wave instability of natural convection flow on inclined surfaces. J. Heat Transfer 107 (1985) 107-111. 19. S.L. Lee, T.S. Chen and B.F. Armaly, Wave instability characteristics for the entire regime of mixed convection flow along vertical flat plates. Int. J. Heat Mass Transfer 30 (1987) 1743-1751. 20. R.A.W.M. Henkes, A.M. Lankhorst and C.J. Hoogendoorn, Structure of the laminar natural convection flow in a square cavity heated from the side for infinitely large Rayleigh number. Proc. A S M E Winter Annual Meeting, Natural Convection in Enclosures, HTD-99 (1988) 9-16. 21. A.E. Gill and A. Davey, Instabilities of a buoyancy-driven system. J. Fluid Mech. 35 (1969) 775-798. 22. Y. Jaluria and B. Gebhart, On transition mechanisms in vertical natural convection flow. J. Fluid Mech. 66 (1974) 309-337. 23. B. Gebhart and R. Mahajan, Characteristic disturbance frequency in vertical natural convection flow. Int. J. Heat Mass Tran~y[Or 18 (1975) 1143-1148. 24. B. Gebhart, Transient response and disturbance growth in vertical buoyancy-driven flows. J. Heat Transfer 110 (1988) 1166-1174. 25. E.R.G. Eckert and E. Soehnghen, lnterferometric studies on the stability and transition to turbulence of a free-convection boundary layer. Proe. Gen. Disc. Heat Transfer, London, (1951) 321 323. 26. K. Kitamura, M. Koike, I. Fukuoka and T. Saito, Large eddy structure and heat transfer of turbulent natural convection along a vertical flat plate. Int. J. Heat Mass Tran~fer 28 (1985) 837-850. 27. D.G. Briggs and D.N. Jones, Two-dimensional periodic natural convection in a rectangular enclosure of aspect ratio one. J. Heat Transfer 107 (1985) 850-854. 28. P. Le Qu6r6 and F. Penot, Numerical and experimental investigation of the transition to unsteady natural convection of air in a vertical differentially heated cavity. Proc. A S M E Winter Annual Meeting, Bifurcation Phenomena in Thermal Processes and Convection, HTD-94 (1987) 75-82.