J . Fluid Mech. (1973), uol. 58, part 2 , p p .

289-312 289
Printed in Great Britain

Two-dimensional Rayleigh-Benard convection

By D. R. MOORE A N D N.0.WEISS
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

Department of Applied Mathematics and Theoretical Physics,

University of Cambridge

(Received 7 July 1972 and in revised form 1 December 1972)

Two-dimensional convection in a Boussinesq fluid confined between free

boundaries is studied in a series of numerical experiments. Earlier calculations
by Fromm and Veronis were limited to a maximum Rayleigh number R 50 times
the critical value R, for linear instability. This range is extended to IOOOR,.
Convection in water, with a Prandtl number p = 6.8, is systematically investi-
gated, together with other models for Prandtl numbers between 0.01 and infinity.
Two different modes of nonlinear behaviour are distinguished. For Prandtl
numbers greater than unity there is a viscous regime in which the Nusselt
number N M 2 ( B / R C ) * ,independently of p . The heat flux is a maximum for cells
whose width is between 1.2 and 1-4 times the layer depth. This regime is found
when 5 < R/Rc 5 p i . At higher Rayleigh numbers advection of vorticity becomes
important and N cc B0.365. When p = 6.8 the heat flux is a maximum for square
cells; steady convection is impossible for wider cells and finite amplitude oscilla-
tions appear instead, with periodic fluctuations of temperature and velocity in
the layer. For p < 1 it is also found that N cc R0.365, with a constant of propor-
tionality equal to 1-90when p < 1 and decreasing slowly as p is increased. The
physical behaviour in these regimes is analysed and related to astrophysical
convection.

1. Introduction
The characteristic features of astrophysical convection are not easily repro-
duced in the laboratory. They can, however, be investigated through numerical
experiments, and more physical understanding may be gained from studying
a planned sequence of idealized models than from attempting to describe every
detail of a star’s convective zone. The simplest relevant problem, that of con-
vection in a Boussinesq fluid, confined between slippery horizontal planes and
heated from below, wa8 first defined by Rayleigh (1916) and has been treated
by many others since (Chandrasekhar 1961; Brindley 1967; Spiegel 1971b). For
computation it is convenient to restrict the flow to two dimensions and this
https://doi.org/10.1017/S0022112073002600

idealized fcrm of Rayleigh-BBnard convection was investigated by Fromm (1965)

and Veronis (1966). We first looked at this problem in order to check a general
program for two-dimensional convection and assess the limitations of our
numerical approximation by comparison with published results. To our surprise
we found that the subject had not yet been exhausted.
Steady convection can be described by expressing the Nusselt number N as
19 FLM 58
 290 D. R. Hoore and N . 0.Weiss
a function of the Prandtl number p and the ratio of the Rayleigh number R to
its critical value R, for linear instability. This new investigation helps to provide
a physical understanding of two-dimensional convection between free boundaries
for the rangep > 0.01 and R 6 1000 R,,limited by the computer time available.
The maximum Rayleigh number reliably examined has been extended from
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

Fromm’s limit of 18.5R, and Veronis’ of 5 0 8 , to a value of lOOOR,. This has

clarified the existence of two distinct regimes, with different dependence of N
on R. A one-third power law was found by Herring (1963), using the mean field
approximation (Spiegel 1967), and by Fromm; it is implicit also in the work of
Veronis. For viscous convection (when advection of vorticity is negligible) we
have confirmed that
N M 2(R/R,)*, (1)
and is independent of p , but this holds only for 5 6 R/R, 5 p*. At high Reynolds
numbers, when advection dominates diffusion of vorticity,

AT M A ( p )(R/R,)0*3s5,
where A(p) E 1.95p-0.05 for p 6.8 and tends to a limit of 1.90 as p +- 0.
Steady two-dimensional convection can be found at any Rayleigh number for
suitably chosen cell widths; its nature is determined by the balance between
advection and diffusion of vorticity. It has sometimes been conjectured that
two-dimensional rectangular cells must always settle down to a steady state.
However, we have found that time-dependent behaviour persists in flattened
cells for p > I , R 2 lOOR,. The coupling between temperature and velocity
maintains periodic finite amplitude oscillations, as in Welander’s (1967) simple
model.
I n the next section we formulate the mathematical system to be solved and
express the equations in dimensionless form. The numerical techniques are then
described in order t o discuss the accuracy of results obtained on a given mesh.
It was necessary to develop finite-difference methods on a rectangular grid so
as to resolve boundary layers at high Rayleigh numbers. The results of the
numerical experiments are presented and discussed in $4.Their physical in-
terpretation is considered in the following section. (Note that the results in § 4
are given in dimensionless units though the discussion of $ 5 uses physical,
dimensional quantities.) Finally, we indicate how the work might be extended
and relate it to other, more complicated problems.

2. Formulation of the problem

I n a Boussinesq fluid the velocity u satisfies the Navier-Stokes equation
https://doi.org/10.1017/S0022112073002600

[aU 1
po z + ( u . V ) u = -vP+pg+p,vv~u (3)

and the incompressibility condition

v . u = 0, (4)
 Two-dimensional Rayleigh-Bdnard convection 29 1
where P is the pressure, g the gravitational acceleration and v the kinematic
viscosity; the equation of state is

P(T)= POP - 4 T - G)I, (5)

Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

where the density p has the value po when the temperature T is equal to Toand a
is the coefficient of thermal expansion. Pressure can be eliminated by taking the
curl of (3); then
awlat = v A (uA W ) -aVT A g + V V ~ W , (6)
where the vorticity w = V A u. Finally, the heat flow equation gives
a q a t = - v .( T U ) +K V ~ T , (7)
where K is the thermometric conductivity.
For two-dimensional convection we choose Cartesian axes with x vertical so
that all flow is confined to the x,z plane and independent of y. Then

u = ( U )0 )W ) ) W = ( 0 )w , 0)) (8)
while, from (4))there exists a stream function @ such that

u = - az/r/az, w = a@/ax, w = - v2$. (9)

Equation (6) then simplifies to
a@ = -v.(wu)-ga-+vv2w)
- 8T
at ax
which, apart from the buoyancy term, is similar to (7).
We restrict the calculation to the rectangular region 0 < x < L,0 < z < d and
adopt the following boundary conditions.
(i) On z = 0, d the normal velocity and tangential stress both vanish, while
the temperatures are fixed : hence

?+9= 0) w = 0 (z = 0,d) (11)

and T = To (X = d), T = To+AT (Z = 0). (12)
(ii) The stream function is assumed to be an odd function of x and periodic,
with period 2L, SO that
$ = 0, w = 0, aT/ax= 0 ( x = 0,L). (13)
I n the absence of convection there would be a uniform vertical temperature
gradient and a thermometric flux Po = tdT/d. When convection occurs, the
https://doi.org/10.1017/S0022112073002600

horizontally averaged thermometric flux is

F = wT-KaT/az (14)
and the normalized flux is given by their ratio, the Nusselt number

(15)
19-2
 2 92 D . R.Moore and N . 0. Weiss
It is usual to express the equations in dimensionless form. We find it con-
venient to use a time unit based on the buoyancy force (analogousto the Brunt-
Vaisala period) so that the dimensionless (primed) quantities become

(x',z') = d-l(x,z), -
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

t' =
AT
I n terms of the Rayleigh number
R =~ uAT~~~Kv

and the Prandtl number P = V/K, (18)

the dimensionless velocity, conductivity and viscosity are then

U' = (pR)-&dK-'U, K' = (pR)-*, V' = (p/R)t. (19)

Finally, suppressing primes, we can rewrite (7), (10) and (9) in the form
aTlat = - V. (Tu)+ K V ~ T ,
_
aw - aT
at - - v . ( w u ) - -+
ax VV2-w

and VZ$ = --w. (22)

Equations (20), (21) and (22) have to be solved in the region 0 < x < A , 0 < z < 1,
where the normalized cell width
h = L/d, (23)
subject to boundary conditions
T=1 (X = 0)) T=0 (Z = l), aT/ax = 0 (Z= 0,h) (24)
and w = 0, $= 0 (x = 0 , A ; z = 0 , l ) . (25)
The effectiveness of convection is now measured by a Nusselt number

while the vigour of the motion is indicated by the average, dimensionless kinetic
energy density

The boundary conditions (24) and (25) already imply that the partial dif-
ferential equations (20)-(22) need only be solved over half of the full convection
cell with 0 < x < 2h; the solutions show mirror symmetry about the planes x = 0,
x = A. We can also take advantage of a further symmetry property: symmetry
https://doi.org/10.1017/S0022112073002600

of $ and T about the midpoint x = +A, z = 4)such that

$(x, 2) = $(A -x, 1- 2 )
and T ( x , z ) + T ( h - x ,1 - 2 ) = 1
is preserved by the equations. This can be verified by Fourier analysis (Veronis
1966). Thus the equations have to be solved only over the quarter-cell 0 < x < A,
 Two-dimensional Rayleigh-Bdnard convection 293
0 < z < 4, with appropriate boundary conditions on z = 4. This simplification
is particularly convenient for computation.
Our boundary conditions allow no horizontal velocity at x = 0. Greater
generality could be achieved by letting $be periodic in x without insisting that
+(O,z) = 0. This would permit lateral sloshing (Deardoa & Willis 1965) and
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

also allow the development of cells with widths hlm, where m is any integer.
We have preferred to impose more stringent boundary conditions so as to
facilitate investigations at high Rayleigh numbers. Thus (24) and (25) inhibit
the transition from a configuration with odd m t'o cells with m even (and vice
versa). The conditions (28) and (29) are even stronger and prohibit the growth of
cells unless m is odd. It is easy to .find particular cases where the solution may
depend upon the boundary conditions that have been adopted, but the overall
behaviour of convection is generally unaffected.

3. Numerical methods
The nonlinear partial differential equations (20)-( 22) were solved by finite-
difference methods; these will be described in detail elsewhere (Moore,Peckover &
Weiss 1973) and need only be summarized here. The variables T , w and @ are
represented on a rectangular grid with spacing Ax, hz; let x j = jAx, z k = kAz
and tn = nAt, wherej = 0, 1, ...,N, and k = 0, 1, ...,N,, and put T2k = T(xj,xk;
tn)
etc. Equations (20) and (21) are advanced by a second-order leapfrog scheme,
centred in time and space (Roberts & Weiss 1966). With Ax = Az, equation (21)
becomes

- (@??l!k+I- @Y:?k-l) ( ! @ $ ? k + l - @ ~ ~ ~ k + l ) W a+$

i,k+l
f (@y$$k-l- @T?:k-l) wzf!l- 2 A x ( T y 2 t k -T?$k)

+ 4 ~ ( w ~ $ ~ k + ~ ~ ?j ,~k +kl f+W E ; ! 1 - 4 W E k ) ] . (30)

This, with a similar equation for TZZl, is solved explicitly for all points withj + k
even, while values of o and T at points withj + k odd are subsequently calculated
a t tn+Q. The stream function @ is obtained by solving Poisson's equation in the
form
@ j + l , k + l + @ j + l , k - l + @j-l, k f l f @ j - l , k - 1 - 4@j.k = - 2Az2wj,
k, (31)
using fast Fourier analysis in the x direction and tridiagonal elimination in the
z direction. Unfortunately, (31) gives @ at the same points as w , which is not
what is wanted for (30), so the required values must be obtained by fourth-order
interpolation from the formula
https://doi.org/10.1017/S0022112073002600

@j, k = % [ ( @ j + l , k -k @j-1, k f @j, k + l f @j, k-1)

f %Az2(wj+l, k f @ j - l , k f wj, k + l +wj,k-l)]- (32)

This method was applied over an entire half-cell (0 < x < A, 0 c x c d ) , with
a maximum of 48 intervals in the range 0 < x < A. The boundary values at
z = 0 , l remain unaltered throughout the calculation and the temperatures at
x = 0 , h are easily computed.
 294 D . R. Moore and N . 0.Weiss
It is essential to establish the accuracy of any difference scheme and, for
convection problems, the maximum Rayleigh number that can be studied on
a given mesh. At high Rayleigh numbers the temperature variation is con-
centrated into narrow boundary layers (see figure 6) of thickness 6 s (2N)-1.
To resolve such a layer adequately it is necessary to have at least three intervals
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

across it (cf. the discussion of convection with fixed boundaries by Schneck &
Veronis 1967). Hence we need to have N, > 6 N ; conversely, for a mesh with
fixed N, we are, from (1)) restricted to Rayleigh numbers such that
R/Rc 6 (33)
For example, a model computed with 18 intervals is accurate only for N < 3
or R/R,6 4;with N, = 34, we need N 6 6 and AIR, 6 25. Comparison of actual
results (see table 1 below) confirms that this is a correct, if somewhat cautious,
(z
criterion: for N, = 18 the error in the Nusselt number rises from 1 at N = 3
to 8 yoat N = 6, though a global quantity like N may provide an over-optimistic
measure of accuracy.
The horizontal length scale associated with variations of temperature and
vorticity is greater than the vertical. So it is possible to use a grid with Ax > Az,
and appropriate difference schemes were therefore developed in order to resolve
the thermal boundary layers at high Rayleigh numbers. The modifications to
(30) and ( 3 2 ) when Ax + Ax are straightforward but the left-hand side of (31)
has to be replaced by a seven-point representation of the Laplacian, requiring
pentadiagonal elimination in the z direction. This version was implemented over
the quarter-cell 0 < x < A, 0 < z 6 4, with Ax > Az (Mooreet al. 1973). The mesh
only covers half the layer depth, so N, = (2Az)-1 and values of $j,-&,z+l etc. are
derived from (38)and (29). This procedure halves the number of points used and
runs were made with Ax 2 iAx on meshes with up to 48 x 100 points (the largest
practicable on the machine). The correctness of the method was checked by
comparison with cases run using both half-cells and the full cell, 0 < x < 2A,
O<z<l.
For R 6 200R, the vertical plumes were adequately described with N, = 24;
this was confirmed by cases run with N, = 48. The width of the cell could be set
to any desired values and the number of vertical intervals chcsen to resolve the
structure at the expected Nusselt number. There is only one boundary layer in
the quarter-cell and we allowed four intervals across it, giving N, s 4 N . The
adequacy of this approximation was again confirmed by comparison with cases
run with N , = 8N. With 48 intervals horizontally and 100 vertically, Nusselt
numbers up to 25 could be examined; this corresponds t o R 2 1700Rc 2 lo6.I n
fact, limited computer time restricted our investigation to R = 1000R,;this case
was run for 4000 times steps, and took 5 h of computing time on an TBN 360144.
https://doi.org/10.1017/S0022112073002600

The various arrays were monitored during eaoh run, and two averaged quan-
tities, the Nusselt number and the kinetic energy density, were calculated from
approximations t o (26) and (27). The Nusselt number was evaluated in the middle
of each vertical interval by calculating the average temperature gradient and
the mean of the convected heat fluxes at the two adjacent levels. The alternative
procedure, of averaging the gradient vertically, is less accurate.
 Two-dimensional Rayleigh-Btnard convection 295
For R < 40R,, cases were run until N converged t o the limit set by machine
accuracy (1 in lo7);the Nusselt number was constant from level to level to one
part in lo4.For R > 40Rcthe solutions oscillated and some cases did not converge
to steady states. The oscillations in converging cases were quite distinct from
those that persisted at a constant amplitude. The former were computed until
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

their amplitude had fallen by at least a factor of l/e (typically 4-8 cycles). The
latter were followed until the oscillations were constant for at least, 5 cycles.
Checks were made to confirm the persistence of unsteady behaviour: for the
cases R/R, = 100, 200, p = 6.8, h = 4 2 runs were made on meshes with 24 x 34
and 48 x 68 intervals and the results agreed to within 2 yo. Furthermore, the
range of variation in N actually increased slightly with higher resolution, as
shown in table 3. I n general it was found that inadequate vertical resolution
suppressed unsteady convection when it was present. On the other hand, when
the solutions were converging to a steady state the high resolution cases con-
verged more rapidly than the coarse ones (cf. Schneck & Veronis 1967).
The finite-differencemethod described above differs slightly from that adopted
by Fromm (1964),who calculated all points at each time level, thereby simplifying
the solution of Poisson’sequation. For the convection calculations he had N, < 26.
Deardorff (1964) used a similar method, with N, = 3 2 ; at R x lOOOR, he found
a Nusselt number of 15, which is 30 yolower than the value we predict. Veronis
(1966) devised an entirely different approach: he expanded T,w and $ in
truncated Fourier series and solved a set of coupled nonlinear ordinary differential
equations for the Fourier coefficients. For a given number of points the Fourier
representation is more accurate (Orszag 1971), though the nonlinear interactions
may be more difficult to compute. Veronis’s expansion had the form

y? = C amnsin (mnxlh)sinnm, (34)

m n

summed over positive integral values of m and n such that m + n < p. His
published results have q 6 10, which enabled him to treat Nusselt numbers up to
7 with errors of less than 1yo.

4. Results
4.1. High Rayleigh number convection in water
We have systematically investigated the effect of varying the Rayleigh number
at a fixed Prandtl number p = 6-8, corresponding to water at room temperature:
Veronis (1966) solved the equations for R < 50R,;we have repeated his calcula-
tions and extended them to R = lOOOR,. This has enabled us to identify different
regimes of behaviour.
https://doi.org/10.1017/S0022112073002600

Viscous regime. Linear theory (Chandrasekhar 1961, chap. 2) predicts in-

stability, at R = R, = 3 3 ~ %4 657.5, to convection in rolls with h = 4 2 . Veronis
accurately treated cases with this cell width for R 6 40R,. Our investigation for
1.1 < R/R,< 40 was done primarily to confirm the accuracy of the computer
program by reproducing his results. Values of the Nusselt number are compared
in table 1 : each value of R was run on a 24 x 18 mesh over the half-cell 0 < x < $,
 296 D.R.Moore and N . 0.Weiss

Veronis (1900), Low resolution, High resolution,

R/Re q = 10 2 4 x 18 48 x 34
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

1 1 -
1.1 1.18 1.16
1.2 1.34 1.32
1.4 1.61 1.58 -
2 2-14 2-12 2.14
3 2.68 2.65 -
4 3.04 3.01
6 3.55 3-52
10 4.24 4-19 -
15 4-85 4-78 -
20 5.33 5.24 5-30(5.32)
30 6.08 5.97 -
40 6.68 6.56 6.67
50 7.16 7.05 - (7.34)

TABLE1. Heat flux as a function of Rayleigh number; p = 6.8, R < 50Rc.Bracketed

figures obtained with Ax + Az, see table 2.

- 0
0
0
I I I I I1111111111ll I I I I I1111111111111 I I I I I I I~1l111111111
1 10 102 103
RlRc
FIGURE 1. Hcat flux as a function of Rayleigh number, p = 6.8. Logarithmic plot of N
against R/RC: - - -, N = 1*96(R/Rc)~; -, N = 1.8(R/R,)0.S65;0 , h = 42; 0,h = 1.

0 < z < 1;some values were repeated on a 48 x 34 mesh. I n models with Ax = Az

it was not possible to obtain the exact width h = 4 2 . However, this seems of
little consequence for low Rayleigh numbers: at R = 10R,Veronis showed that
https://doi.org/10.1017/S0022112073002600

a 10 yovariation in h produced a variation of only 0.5 :(in,N .

The agreement betweeii our finite-difference calculations and Veronis's in-
tegration in Fourier space is very satisfactory, provided enough points are used,
and our results for R = 2022,' p = 6.8 and h = 4 2 are identical with those in
figures 3 and 5 of his paper. Contour plots of T,w and $ at R/RG= 3 , 5 , 7, 10'15,
18.5 have already been published b) Fromm (1965).
 Two-dimensional Rayleigh-Bdnard convection 297
The Reynolds number Re = w/v remains small throughout this range and
convection settles down to a steady state. The variation of N with R is shown in
figure 1. At the onset of instability, R = R, and N = 1. I n the range 1 < R/R, < 5
N rises steeply and then flattens out, with a logarithmic slope decreasing from
2 to an asymptotic value of 4. Por 5 6 RIR, < 40 there is a simple one-third
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

power law:
N = 1.96(R/RC)0333, (35)
with an error of less than 1 yo in the exponent. Herring (1963) predicted this
behaviour, though his mean field calculations overestimated the flux by one-
third. The relation (35) was obtained by Fromm and is also implicit in Veronis’s
results, though he expressed them in the form
[(RIR,) N - 11oc [(R/R,) - 1]1‘26.
The behaviour of the Nusselt number for steady convection with R > 50R, will
be described below and the distinction between the viscous and advective
regimes will be analysed in 5 5.
Dependence on cell width: oscillatory convection.At the critical Rayleigh number
convection first sets in with h = 4 2 . For higher Rayleigh numbers it is frequently
asserted that fluid in an infinite layer will choose the horizontal form that
maximizes the Nusselt number (Malkus 1954). This hypothesis provides a con-
venient criterion for adopting a particular cell width in two-dimensional calcula-
tions, though Foster (1969) and Ogura (1971) have shown that steady-state
solutions are affected by the initial conditions and not determined uniquely by
R and p . Nevertheless, we shall assume that the favoured cell width is indicated
by a maximum of N .
The variation of N with h has been determined at R/R, = 2 and 4 (forp = 0.71)
by Ogura (1971)) at R = lOR, by Veronis and at R = 18.5R, (for p = 1) by
Fromm. The values of h for which N is a maximum lie in the range 2/2 > h 3 1.25
and appear to decrease slightly with increasing R. Table 2 shows the effect of
varying the cell width at higher Rayleigh numbers. At R = SOR, there are
transient oscillations about a state of steady convection for all h 6 2 ; but the
Nusselt number is now a maximum for square cells, with h = 1. When R = IOOR,
the form of the solutions has changed: steady convection is no longer possible
for h = 4 2 . Instead there are finite amplitude oscillations, with a variation of
5 % in N . I n square cells steady convection still occurs and the heat flux is
greater than the average over an oscillation at h = 4 2 . This behaviour persists
when R is increased, as can be seen from figure 2 , where our results are compared
with Fromm’s and with those obtained by Herring (1963) and Huppert (1973)
https://doi.org/10.1017/S0022112073002600

using the mean field approximation. We shall first examine the oscillatory solu-
tions and then return to steady convection with h = 1.
The extreme values of the Nusselt number are listed in table 3. I n a series of
runs on a 24 x 34 mesh with h = 4 2 , the oscillations persisted as the Rayleigh
number was increased to lOOOR, (though the actual numerical values are not very
meaningful for R > 400R,).The boundary conditions prohibit the lateral sloshing
investigated by Deardorff & Willis (1965), and the direction of the flow does not
 298 D. R.Moore and N . 0.Weiss

20 - 5.32 5.35 - - 24 34
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

50 8-74 7.34 7.45 6.61 - 24 34

100 - (9.5) 9.61 - - 24 34
200 (9.8) (11.7) 12.42 11.25 9.75 24 34
12.35 48 48
400 - - 16.22 - - 24 64
15.90 48 64
1000 - - 22.16 - - 48 100
TABLE2. Heat flux as a function of Rayleigh number; p = 6.8, R 2 20R,. Bracketed
figures are average values for oscillatory convection; the mesh intervals are Ax = h/N,,
Az = (2NZ)-'.

0 1 2
h
FIGURE 2. Heat flux as a function of cell width, p = 6.8. Nusselt number plotted against h
for RlR, = 18.5, 50, 100, 200. Points show steady solutions; for oscillatory convection the
average value of N is shown together with the range of variation. - - -, maximum: Nusselt
number as predicted by the mean field equations (Huppert 1973) ; 0-0, Fromm (1965);
0-0, present results.
https://doi.org/10.1017/S0022112073002600

reverse during an oscillation. The variations correspond to fluctuations in

temperature, et c., with a period that is approximately one-third of the circulation
time for the cell.
Figure 3 shows a seqnence of vorticity profiles, streamlines and isotherms
during one oscillation at R = 200Rc. The changes in vorticity are striking-the
maximum value of IwI varies between 0.16 and 0.24in the course of the oscillation-
 Two-dimensional Rayleigh-Bknard convection 299

N
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

100 - - 9.28 9.65 24 34

9-30 9-68 24 68
9.17 9.56 48 68
150 - - 10.4 11.5 24 34
200 8-2 11.4 10.5 13.0 24 34
9-8 12.9 48 68
300 - - 10.1 14.6 24 34
400 - - 10.8 16.3 24 34
TABLE3. Oscillatory convection, p = 6.8.

but the form of the oscillation appears most clearly from the temperature. Three
pairs of hot and cold regions circulate round the cell, affecting the generation
and dissipation of vorticity. The approximately sinusoidal behaviour of the
velocities and heat fluxes, together with their phase relationships, can be seen
from figure 4.Although N varies by only 20 yoa t z = 0,1, the heat flux at z = 4
changes by an order of magnitude. Outside the boundary layer, heat transport
is dominated by rising thermals and their sinking counterparts. The mechanism
producing this behaviour will be discussed in § 5.
Advective regime. When R 2 50R, the heat flux is a maximum for square cells.
VaIues of the Nusselt number calculated with h = 1 for 50 < R/R, 6 1000
are given in table 2 and also plotted as a function of the Rayleigh number in
figure 1. The one-third law no longer holds over this range; instead, the heat flux
follows a power law of the form
N = 1*8(R/R,)0’36s. (37)
(The exponent in (37) has an error of less than 1 yo;the slope of the curve
diminishes slightly towards the highest value of R but this decrease is not
significant.) This value of N is greater than that from equation (35) for R > 15Rc.
At high Rayleigh numbers the favoured cell size is narrower and convection in
square cells becomes more efficient than the viscous regime described by (35).
Vorticity, instead of being dissipated locally in the vertical plumes, is advected
into the thermal boundary layers and lost there.
The detailed structure of the flow for R = lOOOR,, the highest Rayleigh
number studied, is depicted in figure 5 . Comparison of the streamlines and iso-
https://doi.org/10.1017/S0022112073002600

therms with those for R = 20Rc shows that the noillinear features have become
even more pronounced. The streamlines are, as always, similar but the thermal
boundary layers at z = 0 , l are very narrow, while the rising plume has become
attenuated and curves sharply round. The central isothermal region occupies
almost all of the cell. The mean temperature gradient for R = 20R, is stably
stratified in this central region (Veronis 1966) but this reversal can no longer be
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.
https://doi.org/10.1017/S0022112073002600

300

FIGURE

plots is one-sixth of the period.

D.R.Moore and N . 0.Weiss

3. Oscillatory convection, R = ZOOR,, p = 6.8,h = 42. Vorticity profiles, stream-

lines and isotherms during EL single oscillation. The time increment between successive
 Two-dimensionalRayleigTcBdnard convection 301
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

FIGURE 4. Oscillatory convection: variation with time of the velocities and heat flux.
(a)Maximum values of u and w, the horizontal and vertical velocities. (a) Nusselt number
N , evaluated at z = 0 and z = &. Times are given in terms of the oscillation period
(3.6 dimensionless units).
https://doi.org/10.1017/S0022112073002600

5 . Convection in the advective regime, for R = lOOOR,, p = 6.8, h = 1. (a)Stream-

FIGURE
lines; ( b ) vertical profile of mean temperature; (c) isotherms.
 302 D. R.Moore and N . 0.Weiss

I -
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

0.5 Jc=*

0 0.5 1.0 0.8 0.4 0

T W RIB,
FIGURE 6. (a)Vertical profiles for R = lOOOR,, p = 6.8, h = 1 . Temperature: T ( i ,z ) and
T ( 0 ,z ) ;normalized vertical velocity: ~ ( 0z)., ( 6 ) Advective regime: variation of velocities
and boundary-layer thickness with Rayleigh number for p = 6.8. Logarithmic plots (with
slopes shown in brackets) of the thickness of the horizontal boundary layer SIT( - 0.32), the
width of the vertical plume Sv ( - 0.29) and then of the maximum horizontal and vertical
velocities (0.22 and 0.18 respectively), as functions of RIR,.

1.o

w o

- 1.0
0 0.5 1 .o
https://doi.org/10.1017/S0022112073002600

xlh
FIGURE 7. Horizontal profiles for R = 20R, ( A = 4 2 ) and R = 1000R, ( A = 1) with
p = 6.8. The temperature T(z/h,4) and the normalized vertical velocity w(z/A,$) are
plotted as functions of x/A. T : -, R = lOOOR,; W-W, R = 2013,. W : .-a,
R = 1000R,; A-A, R = 20R,.
 Two-dimensional Rayleigh-Bdnard convection 303
discerned a t R = lOOOR,. The temperature in the central core is virtually con-
stant, though the kinks in the mean temperature profile just before the boundary
layers are real and caused by the plumes spreading out horizontally as they
approach the boundaries. This behaviour appears clearly in vertical profiles of
T a t x = 0, i, shown in figure 6 (a).The vertical velocity in the plumes varies
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

approximately sinusoidally with x : acceleration by the buoyancy force is counter-

acted by an adverse pressure gradient.
Figure 7 shows horizontal profiles of temperature and vertical velocity at
z = 8. The vorticity in the isothermal core is nearly constant and the central
region therefore rotates with a uniform angular velocity, accelerating slightly
in the plumes. This becomes more marked at high Rayleigh numbers and the
temperature profiles show that convection of heat is increasingly concentrated
into plumes. Although the central core becomes more nearly isothermal there
is always a region with a small reversed horizontal temperature gradient between
the core and the plumes. This is a consequence of advecting the temperature
(cf. the magnetic fields described by Weiss 1966) and causes a local decrease in
vorticity.
The variation of the convective velocities and boundary-layer thicknesses with
R is plotted in figure 6 ( b ) .For the horizontal boundary layers the thickness is
defined by the points where T = 0.5; for the plumes it is calculated similarly
from T ( z ,4). The width of the vertical plumes is then about 1.7 times the thermal-
boundary-layer thickness. For the one-third power law of (35) we wouId expect
the boundary-layer thickness 6, the dimensionless velocity IuI and the kinetic
energy density 8 to vary as powers of R:
6ccR?+, \ U I c c B + , 8CcH (38)
(see Q 5 ) . The actual thicknesses of boundary layers and plumes are consistent
with (38) but the r.m.s. velocity increases more rapidly, as WlS. Thus the
enhanced heat flux of equation (37) is apparently achieved by increasing the
velocity without significantly changing the temperature fieId.

4.2. Ee
’ct ofvarying the Prandtl number
Variation of heat .flux. Veronis calculated the Nusselt number N ( R , p ,A) for
models with h = 4 2 , 100 2 p 3 0.005 and R < 20R,. We have computed heat
fluxes for h = 1, p 2 0.01 and R < lOOOR,. These results are shown in table 4:
values of N are accurate to within 1 % except for runs with R = IOOOR,, where
the error may be as high as 2.5 yo.Models with infinite Prandtl number required
a modified version of the program, developed for geophysical problems (McKenzie,
Roberts & Weiss 1973). The Nusselt numbers for p = 00 agree with those com-
puted by Straus (1972) for R < 60& and are consistent with an asymptotic
https://doi.org/10.1017/S0022112073002600

dependence of the form N = 2-OO(R/Rc)). We have also recomputed N(20Rc,

0.01, 4 2 ) for comparison with Veronis’s result and they agree to within 1 % .
The variation of N with h is somewhat greater for the advective than for the
viscous regime.
Veronis observed that the values in table 4 show a definite pattern. For p < 1
the heat flux is independent of the Prandtl number; as p increases the Nusselt
 304 D. R.Moore and N . 0.Weiss

P
A
7 7
R/Rc 0.01 0.1 1.0 6.8 16 100 03

6 3.61 3.61 3-58 3.55 3.54 3.55 3.55

Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

10 4-39 4.39 4-38 4.24 4.22 4-24 4-24

15 5.11 5.11 5.06 4.85 4.82 4.85 4.83
20 5.68 5.67 5.62 5.33 5.29 5.33 5.29
5-62 5.62 5-53 5-35 5.35 5-36 5.37
50 7.91 7.92 7.81 7.45 7.31 7.30 7.32
100 - io.18 10.09 9.61 9.32 - 9.25
200 - 13-01 12.93 12.35 11.84 - 11.52
400 - - 16.55 15.90 - - 14.45
1000 - - 22-93 22.16 - - 19.47
TABLE4. Nusselt number as a function of R and p . Values in bold-face for h = 1. Other
values for h = 42 (results for p = 0.01, 0.1, 1.0, 6.8, 100 from Veronis 1966).

'06 i 1
I
I
/
1
I
I
111
1
I
1

I
1 1 I I I J
1 o2
1 104

23
FIGURE8. Regimes of convection in the R, p plane. I, low Rayleigh number convection;
11, viscous regime ; 111, advective regime ; the dashed line separates high and low Prandtl
number behaviour.

number declines slightly between p = 1 and p = 6.8. Thereafter, the heat flux
is once again independent ofp for Rayleigh numbers less than 20Rc.The variation
of N with R provides more information: for p = 0.1 and p = I we find that
N cc R0.362 over the range 20 < RIB, 6 200; while, from Veronis's results,
d In N/d In R z 0.365 when R 20R, and p < 0,025. Thus convection at low
N
https://doi.org/10.1017/S0022112073002600

Prandtl numbers shows the same power-law dependence as we 'found for the
advective regime at p = 6.8.
We can summarize this behaviour in terms of a general power law
fl = 4 P )( R I R 2 . (39)
The R,p plane can be divided into three regions (figure 8). In region I the Nusselt
number rises steeply from unity at R = Reand the heat flux can be calculated
 Two-dimensional Rayleigh-Be'nard convection 305
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

9. The effect of varying Prandtl number, with R = 50R,,h = 1. Vorticity profiles

FIUURE
and isotherms (a) for p = 0.01, (b) for p = 100.

by high-order expansion methods (Malkus & Veronis 1958; Kuo 1961). Region I1
is the viscous regime discussed above, with N given by (35) SO that
p = 0.333 f 0.002
and A is independent of p in (39). This region only exists for Prandtl numbers
greater than unity. Region I11 is the advective regime. Here /3 = 0.365 k 0.003
and A is a slowly decreasing function of p : A(16) = 1.7, A(6.8) = 1-8 and
A @ )+ 1.90 as p -+0. As the Prandtl number increases, the efficiency of this
advective regime is impaired and the viscous regime extends to higher Rayleigh
https://doi.org/10.1017/S0022112073002600

numbers. Ultimately, when p is infinite, the advective term in (10) vanishes for
all finite R. The boundary between regions I1 and I11 is poorly determined from
our data; for p = 6-8, 16 it is given by RIRc w p1.593. This is consistent with the
expression
R/R,= (1.1 f O.l)p+. (40)
From (35), (37) and (40)it follows that A ( p ) ,N 1.95p-oo5 for p 2-6.8.
20 FLM 58
 306 D . R. Moore and N . 0.Weiss
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

FIGURE 10. The dependence of the vorticity on the Rayleigh number, for p = 6.8.
Vorticity profiles: ( a ) for R = 20R,, h = 4 2 , (b) for R = 100R,, h = 1 and (c) for
R = 1000R,, h = 1.

General behaviour. For a fixed value of R the Nusselt number changes by only
I % when p is altered by lo4. Figure 9 shows the isotherms for R = 50R, with
p = 0.01 a n d p = 100. They too are scarcely affected, though the velocities vary
in order to maintain a constant flux and the dimensionless kinetic energy density
Q is inversely proportional to p .
The vorticity shows more striking differences and provides the best means of
distinguishing the different regimes. Its behaviour at low Rayleigh numbers is
shown in figure 5 of Fromm’s paper; for R < 5R,the vorticity resembles the
https://doi.org/10.1017/S0022112073002600

linear solution with w cc sin nxlh sin m. Profiles for p = 6.8 with R/Rc = 20, 100
and 1000 are plotted in figure 10. When R = 20R, the central peak has separated
into two humps, centred on the regions with horizontal temperature gradients,
where vorticity is generated. I n the viscous regime advection is unimportant
and vorticity is dissipated where it is created, in the plumes. As R is increased,
vorticity is carried towards the horizontal boundaries and dissipated by viscosity
 Two-dimensional Rayleigh-Be'nard convection 307
near the corners of the cell, while w scarcely varies in the central core; this
behaviour characterizes the advective regime.
The effect of changing the Prandtl number is shown by the vorticity profiles
for a fixed Rayleigh number, R = 5OR,,at p = 0.01 and p = 100, plotted in
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

figure 9. For p = 100, the vorticity is concentrated on the sides of the thermal
plumes and w is approximately harmonic in the centre, with a saddle-point at
x = z = 4.Whenp = 0.01, vorticity is advected round, producing a band of high
vorticity with a flattish minimum in the centre, and the angular velocity in-
creases to a maximum at the edges of the plumes. Thus convection a t high and
low Prandtl numbers is distinguished by different distributions of vorticity.

5. Physical interpretation
5.1. Steady convection
The simplest behaviour occurs for the viscous regime in region I1 of figure 8,
where N cc R). This viscous regime can easily be understood in terms of a simple
boundary-layer model (Turcotte & Oxburgh 1967; Robinson 1967). We shall
henceforth revert to physical (dimensional) units. Consider a roll with width
L d. Let w be the vorticity and w the velocity in the plumes; then at the
N

perimeter of the isothermal core

w wd. N

I n the horizontal boundary layers,

whence w N (i)", (43)

if we assume that u w and thus that horizontal boundary layers and vertical
N

plumes have the same thickness 6.I n the plumes, vorticity is created and dis-
sipated locally, so that
vv2w ga aTpx (44)

or, from (41), (45)

. .
Then, from (17), (43) and (45),
(S/d) N R-t, w N R*K/d,
while the Nusselt number
N N (d/S) cc Ri.
These predictions are consistent with the numerical results summarized in (38),
if the units are transfclrmed according to equation (19). The overall Reynolds
number
https://doi.org/10.1017/S0022112073002600

Re = wd/v cc p-IR* (48)

and so the transition from viscous to nonlinear flow occurs a t a Rayleigh number
R cc pQ.Equation (40) shows that this is consistent with the numerical results
and the one-third law applies only for Re 5 10. Hence the vorticity in the core
is approximately harmonic and the steady 'inviscid' solution with w constant in
the central region (Batchelor 1956; Robinson 1967) is not attained in this regime.
20-2
 308 D. R. Moore and N . 0. Weiss
When advection of vorticity becomes important vorticity is dissipated around
the stagnation points where the plumes impinge on the boundaries. It is no
longer feasible to construct a simple model but region I11 of figure 8 can still
be subdivided into two regimes. At large Prandtl numbers an element of fluid
gains vorticity as it rises or falls but then loses what was gained through friction
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

as it traverses the thermal boundary layers. Thus the flow is essentially com-
posed of rising and sinking plumes. The vorticity w, in the central core is approxi-
mately constant and

(49)

The advective regime for p < 1 is different (Veronis 1966). Viscosity is too
feeble for the fluid to lose more than a small fraction of its vorticity as it passes
along the horizontal boundaries. So the fluid gains vorticity as it circulates round
the cell and this process continues until the total vorticity is so great that viscous
losses can balance what is gained. This value of u is much greater than wo, the
vorticity gained in one passage through a plume: for example, with p = 0.01
and R = 50R,, the central vmticity is 9000 times the value given by (49). Con-
vection is no longer dominated by the plumes: instead the flow resembles an
inertial flywheel driven by a fixed torque, with friction proportional to its
angular velocity. I n the isothermal core the Reynolds number is high and so the
vorticity is virtually constant in a steady state (Batchelor 1956). Th’is con-
figuration is only slowly achieved and the Nusselt number creeps exponentially
to its final value at a rate determined by the viscous time scale for the central
core. Since the sign of aT/ax reverses between the plumes and the central core,
as shown in figure 7, w diminishes before settling down to a constant value. For
p = 0.01. and R = 50R,, w falls by 14 yo of its maximum value and the profiles
in figure I0 outline a circular volcano. At much higher Rayleigh numbers we
would expect an elevated racetrack enclosing a flat interior basin.
Oscillations. I n the advective regime for Prandtl numbers greater than unity
there are periodic variations in all quantities. For square cells these fluctuations
slowly decay as the convection converges t o a steady state. Finite amplitude
oscillations persist in wider cells. A linearized analysis might show that steady
convection is either stable to small perturbations (so the oscillatory modes are
damped) or unstable (so that oscillationsgrow until they are limited by nonlinear
effects), depending on the width. We shall attempt to describe the mechanism
that maintains the steady oscillations. The same phenomena appear in a simple
model put forward by Welander (1967).
The oscillations themselves are thermal in origin. Suppose that an abnormally
cold parcel of fluid falls from the upper surface: then it will remain colder as it
https://doi.org/10.1017/S0022112073002600

circulates round the cell. On a fixed Eulerian mesh there will be variations of
temperature, repeated each time the cold blob passes by. This behaviour is
visible in the isotherms of figure 3. The cold plume bulges and the hot rising
plume is pinched off as colder fluid enters it. At the same time, symmetry requires
a parcel of hotter fluid, following the cold fluid around, and in general we might
find m such pairs. The symmetry assumed for high resolution runs allows only
 Two-dimensional Rayleigh-Be'nard convection 309
modes with m odd. However, no even modes were found in models computed for
a full half-cell (0 Q x Q L, 0 < z < d).
Thermal diffusion normally eliminates these fluctuations: for example, at
R/Rc = 400 and 1000 with h = 1 the oscillations showed periodicities corre-
sponding to m = 1 and to m = 3 which decayed exponentially on the thermal
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

time scale. The oscillations are maintained only if the velocity is coupled to the
thermal field, so that cold fluid moves rapidly past the lower (hot) boundary and
slowly past the upper one (Welander 1967). This coupling can just be detected
in the streamlines of k u r e 3 and it is clearly seen in the vorticity profiles. Now
vorticity generated by horizontal temperature gradients in the plumes is advected
into the horizontal boundary layers. The vorticity near the boundaries depends
also on Iocal temperature gradients and viscous dissipation. Colder fluid reaches
the lower boundary with more vorticity and therefore travels faster, losing less
heat as it goes across. As this fluid crosses the cell the horizontal temperature
gradient steepens and vorticity continues to be created, reducing the effect of
viscous losses. Behind the cold blob, however, the isotherms are drawn out and
friction retards the flow until the temperature rises and isotherms reconnect.
There is thus an adverse temperature gradient which further reduces the speed.
Conduction of heat from the boundary then causes the formation of a hot spot.
This follows at a distance of about *L behind the cold region, corresponding t o
an oscillation with m = 3, the only mode found in the numerical experiments.
This process requires that vorticity be advected into the boundary layer, where
viscous diffusion must be significant. For p < 1 the overall vorticity cannot be
altered sufficiently in an oscillation, so large Rayleigh and Prandtl numbers are
necessary for oscillations to occur. Increasing the cell width enhances the effect
of viscosity, as can be seen by comparing the vorticities in figure 10. Eddies with
>
h 1 might be expected to split in three. However, we have found no indication
of iission for cells with R = 200Rc,h < 4. The behaviour of an infinite layer is not
obvious from our numerical experiments on a single cell. It is probable that
infinitesimalperturbations would developinto steady convection,though suitable
initial conditions would allow oscillatory solutions to persist. Once some pattern
of cellular convection is established, the total number of cells can change only
by the consumption of small cells or the fission of larger ones, while the cell widths
grow more uniform on a time scale that is determined by diffusion.

6 . Conclusion
These numerical experiments have enlarged our physical understanding of
two-dimensional convection between free boundaries and it is not obvious that
extensions to higher Rayleigh numbers, or to smaller Prandtl numbers, will
https://doi.org/10.1017/S0022112073002600

introduce any new phenomena. Some points still need to be clarified. The transi-
tion to the advective regime - the curve separating regions I1 and 111of figure 8 -
might be more precisely determined. The exponent, /3 = 0.365, in (37) appears to
be constant over the range 50 < R/Rc 6 1000, but much greater values of R
would have to be investigated in order to determine whether the heat flux has
the form N 0~ fi[h(WW', (50)
 310 D.R.Moore and N . 0.Weiss
, where y and R,are constants of order unity. Similarly, we could find no evidence
that the favoured cell width changes for R > 50Rc,though it might decrease with
increasing p .
Our results provide the first clear evidence for a power law with an exponent
greater than 8. Howard (1963)used integrals of (3)and (7) to obtain an upper
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

bound of 4 for the exponent in the limit of large R.For a single horizontal mode
satisfying (4) and the same integral constraints, the exponent is reduced t o 8
(Howard 1963; Straus 1973); this provides an upper bound for all mean field
and single mode calculations. It is interesting that the value obtained from
numerical experiments lies between & and 8. However, when many modes are
allowed (corresponding to fine structure in thermal boundary layers and the
vertical plumes) the upper bound has an exponent of 4 for R % 1 (Busse 1969).
The computed value of 0.365 lies well below this limit.
Gough, Spiegel & Toomre (1973)have recently studied three-dimensional con-
vection between free and fixed boundaries, obtaining asymptotic results for
a cellular model together with numerical solutions of the modal equations up to
R = 1025.For two-dimensional rolls their model is identical to Herring’s; for
three-dimensional motion with free boundaries they find that N cc Rfwhen p is
infinite but that N cc [R In R]*when p is of order unity. Numerical experiments
confirm their prediction that for high Prandtl number the RQ law holds at much
lower Rayleigh numbers for free than for fixed boundaries.
Two-dimensional rolls are found in convection experiments at low Rayleigh
numbers, but it is doubtful whether our results are relevant to convection be-
tween fixed boundaries, where the thermal boundary layer itself becomes un-
stable. Busse (1967)showed that when p % 1 rolls become unstable to three-
dimensional disturbances for R 2 13R, and the resulting bimodal convection has
been investigated experimentally by Busse & Whitehead (1971). Although con-
vection becomes time-dependent a t higher Rayleigh numbers, or when p < 1
(Rossby 1969),it is uncertain whether the oscillations correspond to circulating
hot and cold spots (Krishnamurti 1970)or to displacements of the cells (Willis &
Deardorff 1970).
The onset of convection between free boundaries has been studied in the
laboratory (Goldstein & Graham 1969) but there is no experimental evidence
on the persistence of two-dimensional rolls. The thermal boundary layers remain
stableifdlnN/dlnR 2 &(Busse1967)andStraus(1972)hasconhmed,forp > 1,
that rolls are stable to the three-dimensional perturbations that lead to bimodal
convection. The viscous regime described above should therefore be a stable
solution to the full three-dimensional problem. At low Prandtl numbers, however,
rolls are unstable to oscillatory three-dimensional disturbances which generate
vertical vorticity and so produce wavelike distortions of the rolls (Busse 1972).
https://doi.org/10.1017/S0022112073002600

These probably correspond to the oscillations found in experiments with fixed

boundaries. Thus it seems likely that the two-dimensional advective regime
described here will become unstable to wavelike disturbances and that these will
eventually lead to fully three-dimensional cellular convection.
Convection between fixed plates is dominated by viscous boundary layers but
stress-free boundaries are more appropriate in astrophysics. Apart from its
 Two-dimensional RayleigTGBdnard convection 311
intrinsic interest, the idealized model that we have studied is a necessary pre-
liminary to more complicated problems, such as penetrative convection and the
effects of magnetic fields or compressibility. I n stars the radiative conductivity
is high and the effective Prandtl number is small. It has been conjectured that
the heat flux should be independent of both v and d, so that
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

N cc (pR)). (51)
I n calculating models of convective stellar atmospheres the depth d is replaced
by a mixing length; it is generally assumed that fluid loses its kinetic energy
after rising through a mixing length, so that

(IPR)) (52)
>
for p R 1 and the heat flux is independent of both K and v (Spiegel 1971a, b ) .
Our two-dimensional models with p R of order unity show no such variation of
N with p . Stellar convection requires a study of three-dimensional models a t
low Prandtl numbers, with p R 2 lo3.

We wish to thank Dr D. 0.Gough for many discussions which have improved

our understanding of convection. The program was developed in collaboration
with Dr R. S. Peckover and most of the computing was done on the IBM 360/44
at the Institute of Theoretical Astronomy; we are grateful to Mr N. J.Butler for
his co-operation and assistance in producing the diagrams. Further computations
were carried out at the Max-Planck-Institut fur Astrophysik, Munich, and we are
grateful to Dr K. von Sengbusch for his assistance and advice. This paper was
written mainly at the Max-Planck-Institut fiir Astrophysik and improved as a
result of comments by Dr 3’. H. Busse, Dr Gough and Dr H. E. Huppert. D. R. M.
is grateful to I.B.M. (United Kingdom) Ltd. for a graduate fellowship held while
this research was being done.

REFERENCES
BATCHELOR, G. K. 1956 J . Fluid Mech. 1, 177.
BRINDLEY,J. 1967 J . In&. Math. Applics. 3, 313.
BUSSE, F.H. 1967 J . Math. & Phys. 46, 140.
BUSSE, F.H. 1969 J . Fluid Mech. 37, 457.
BUSSE, F.H. 1972 J . Fluid Mech. 52, 97.
BUSSE,F.H. & WHITEHEAD,J. A. 1971 J . Fluid Mech. 47, 305.
CHANDRASEKHAR, S. 196 1 Hydrodynamic and Hydromagnetic Stability. Clarendon Press.
DEARDORFF, J. W. 1964 J . Atmos. Sci. 21, 419.
DEARDORFF, J. W. & WILLIS, G. E. 1965 J . FZuid Mech. 23, 337.
https://doi.org/10.1017/S0022112073002600

FOSTER,T. D. 1969 J . Fluid Mech. 37, 81.

FROMM,5. E. 1964 Methods in Comp. Phys. 3, 346.
FROMM,J. E. 1965 Phys. Fluids, 8, 1757.
GOLDSTEIN,R. J. & GRAHAM,D. J. 1969 Phys. Fluids, 12, 1133.
GOUUH,D. O., SPIEGEL,E. A. & TOOMRE, J. 1973 Model equations for turbulent con-
vection. To be published.
HERRING,J. R. 1963 J . Atmos. Sci. 30, 325.
 312 D. R.Moore and N . 0.Weiss
HOWARD, L. N. 1963 J . l%sMe&. 17,405.
HUPPERT,H. 1973 The asymptotic solution of the mean field equations of thermal
convection. I n preparation.
KRISHNAMURTI, R. 1970 J. Fluid Mech. 42, 309.
KUO,H. L. 1961 J. Fluid Mech. 10,611.
Downloaded from https:/www.cambridge.org/core. Cornell University Library, on 12 Jul 2017 at 08:47:46, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

MCKENWE,D. P., ROBERTS, J. M. & WEISS, N. 0. 1973 Numerical experiments on

convection in the earth's mantle. I n preparation.
MALKUS, W. V. R. 1954 Proc. Roy. SOC.A 225, 185.
MALKUS,W. V. R. & VERONIS,G. 1958 J. Fluid Mech. 4, 225.
MOORE,D. R., PECKOVER, R. S. t WEISS, N. 0. 1973 Difference methods for time-
dependent two-dimensionalconvection. I n preparation.
OGURA,Y. 1971 J . Atmos. Sci. 28, 709.
ORSZAG, S. 1971 J. Fluid Mech. 49, 7 5 .
RAYLEIGH, LORD 1916 Phil. Mag. 32, 529.
ROBERTS, K. V. & WEISS, N. 0. 1966 Math. Gomp. 20, 272.
ROBINSON, J. L. 1967 J. Fluid Mech. 30, 577.
ROSSBY,H. T. 1969 J. Fluid Mech. 36, 309.
SCHNECK, P. & VERONIS, G. 1967 Phys. Fluids, 10,927.
SPIEGEL, E. A. 1967 I n Aerodynamic Phenomena i n Stellar Atmospheres (ed.R. N. Thomas),
p. 348. Academic.
SPIEGEL, E. A. 1971a Comm. Astrophys. & Space Phys. 3, 53.
SPIEGEL,E. A. 1971b Ann. Rev. Astro. Astrophys. 9, 323.
STRAUS, J. M. 1972 J. Pluid Mech. 56, 353.
STRAus, J. M. 1973 I n preparation.
TURCOTTE, D. L. & OXBURGH, E. R. 1967 J. Fluid Mech. 28, 29.
VERONIS,G. 1966 J . Fluid Mech. 26, 49.
WEISS,N. 0. 1966 Proc. Roy. Soc. A293, 310.
WELANDER, P. 1967 J. Fluid Mech. 29, 17.
WILLIS,G. E. & DEARDORFF, J. W. 1970 J. PZuid Mech. 44, 661.
https://doi.org/10.1017/S0022112073002600