[
Abstract
The Boussinesq equations for Rayleigh–Bénard convection are simulated for a cylindrical container with an aspect ratio near 1.5. The transition from an axisymmetric stationary flow to timedependent flows is studied using nonlinear simulations, linear stability analysis and bifurcation theory. At a Rayleigh number near , the axisymmetric flow becomes unstable to standing or travelling azimuthal waves. The standing waves are slightly unstable to travelling waves. This scenario is identified as a Hopf bifurcation in a system with O(2) symmetry.
Standing and travelling waves in cylindrical Rayleigh–Bénard convection]
Standing and travelling waves in cylindrical Rayleigh–Bénard convection
K. Borońska and L. S. Tuckerman]
Katarzyna Borońska
1 Introduction
Rayleigh–Bénard instability in a fluid layer heated from below in the presence of gravity is the classic prototype of pattern formation. A new chapter in its investigation began with the increase of computer performance that made feasible threedimensional, nonlinear, highresolution simulations of the Boussinesq equations governing this system.
We are interested in a fluid layer confined in a vertical cylinder whose upper and lower bounding surfaces are maintained at a temperature difference measured by the Rayleigh number. The conductive solution for this system is a motionless state with a uniform vertical temperature gradient. This solution is stable up to a critical Rayleigh number , whose value depends on the aspect ratio radius/height. Above , convective motions appear and form various roll structures.
A summary covering the developments since the mid 1980s for convective systems with large aspect ratio () can be found in ?BodPesAhl. In such domains a rich variety of patterns was reported: “Pan Am” patterns ?[arches with several centres of curvature, see][]AhlCanSte, straight parallel rolls ?[]Cro,CroGalPoc, concentric rolls ?[targets, see][]KosPal,CroMorSch, one and several armed rotating spirals ?[][]PlaEgoBodPes, targets with dislocated centre Croquette (1989), hexagonal cells ?[][]CilPamPer and spiraldefect chaos ?[][]MorBodCanAhl. A large overview on convective phenomena observed experimentally before this time can also be found in Koschmieder (1993).
We focus here on cylinders with moderate aspect ratio . The flow structure then depends strongly on system geometry. For this regime, the stability of the conductive state was well established in the 1970s–1980s by Charlson & Sani (1970), Stork & Müller (1975) and Buell & Catton (1983). Critical Rayleigh numbers are about for , increasing steeply for lower and decreasing asymptotically towards for . Charlson & Sani (1970) estimated by a numerical variational technique the onset of axisymmetric convection in cylinders of aspect ratios between 0.5 and 8, with insulating and conducting sidewalls. They found the critical Rayleigh numbers ( for , decreasing for higher ) and the corresponding number of rolls. They then generalised this analysis Charlson & Sani (1971), including nonaxisymmetric modes and predicting and corresponding critical azimuthal wavenumbers. Stork & Müller (1975) observed experimentally convective patterns in annuli and cylinders of aspect ratio , varying the sidewall insulation. Their critical Rayleigh numbers were in good agreement with those predicted by Charlson and Sani. Rosenblat (1982) investigated convective instabilities numerically for freeslip boundary conditions, using a severely truncated expansion in a small number of eigenmodes. He described nonaxisymmetric motions existing just above onset for aspect ratios between 0.5 and 2.0. Finally, Buell & Catton (1983) described how the onset of convection is influenced by the ratio of the fluid conductivity to that of the wall, by performing linear analysis for the aspect ratio range . They determined the critical Rayleigh number and azimuthal wavenumber as a function of both aspect ratio and sidewall conductivity, thus completing the results of the previous investigations, which considered either perfectly insulating or perfectly conducting walls. These results were confirmed by ?Marques. The flow succeeding the conductive state is threedimensional over large ranges of aspect ratios, contrary to the expectations of Koschmieder (1993).
The stability of the first convective state, depending on both aspect ratio and Prandtl number, has been investigated mainly for situations in which the primary flow is axisymmetric. Charlson & Sani (1975) attempted to predict numerically the stability of the primary axisymmetric flow, but the resolution available at that time was inadequate to the task. ?MulNeuWeb investigated convective flows experimentally and theoretically. They observed axisymmetric flows for and nonaxisymmetric flows for . Hardin & Sani (1993) calculated weakly nonlinear solutions to the Boussinesq equations for several moderate and small aspect ratios. They found a bifurcation from the axisymmetric state towards a mode with azimuthal wavenumber for , and .
The most complete numerical study of secondary convective instabilities for moderate aspect ratio cylinders was performed by ?WanKuhRat. For cylinders with insulating sidewalls and , the primary bifurcation to convection occurs at and leads to an axisymmetric flow whose stability was investigated for Prandtl numbers 0.02 and 1. Wanschura et al. predicted the succeeding flows to be steady, except over a narrow aspect ratio range at , where they found oscillatory instabilities at towards flows with azimuthal wavenumbers and . The primary aim of this paper is to provide a more detailed description of these bifurcations.
?TouHadHen numerically investigated the stability of the conductive state for aspect ratios and . They described the main critical modes and established a diagram of primary bifurcations, including unstable branches. They also found a secondary bifurcation point , at which the axisymmetric flow becomes unstable towards a tworoll flow and calculated for and .
An interesting experimental study was carried out by ?HofLucMul. Varying the Rayleigh number through different sequences of values, for fixed parameters and , they obtained several different stable patterns for the same final Rayleigh number. They also reported a transition from an axisymmetric steady state towards azimuthal waves. Our numerical simulations of this phenomenon are the subject of a separate investigation.
More recently convective patterns were numerically investigated by Rüdiger & Feudel (2000) and by Leong (2002). Rüdiger and Feudel found stability ranges for multiroll patterns, targets and spirals for , . Leong observed several steady convective patterns for aspect ratios 2 and 4 and Prandtl number , all of which were stable in the range , and calculated the heat transfer for each pattern.
Convective systems often display oscillatory behavior. In binary fluid or rotating convection, the primary bifurcation is usually to periodic states, while in RayleighBénard convection, periodic behavior occurs as a secondary bifurcation. The oscillatory and skewvaricose instabilities of long straight parallel rolls calculated in, e.g. Clever & Busse (1974) and Busse & Clever (1979), are manifested as travelling waves along rolls and as periodic defect nucleation Croquette (1989); Croquette et al. (1986); Rüdiger & Feudel (2000); rotating spirals were observed by the same investigators; and radially propagating patterns of concentric rolls were observed by Tuckerman & Barkley (1988). However, none of these manifestations of oscillatory behavior resemble the azimuthal waves we describe in this study.
Competition between standing and rotating azimuthal waves has been extensively studied in thermocapillary convection, driven by surfacetension gradients. For example, competition between rotating and standing waves is observed on the upper free surface of an open cylindrical container by Sim & Zebib (2002) and in the midplane of a cylindrical liquid bridge with free outer surface by ?Leypoldt, both of aspect ratio 1. These azimuthal waves are very similar to those we describe in this study; however, such flows are uncommon in the Rayleigh–Bénard (buoyancydriven) convection literature.
We wished to study in detail the timeperiodic nonaxisymmetric states in cylindrical Rayleigh–Bénard convection resulting from the bifurcation found by Wanschura et al. (1996). Hence we have simulated numerically the loss of stability of the first convective axisymmetric solution undergoing an oscillatory bifurcation for and . In this paper we describe the results of nonlinear simulations and linear stability analysis, which identify the scenario in terms of bifurcation theory in systems with symmetries.
In addition to obtaining results particular to cylindrical Rayleigh–Bénard convection with these parameter combinations, our purpose is to demonstrate how numerical and theoretical techniques can be combined in order to obtain a complete bifurcationtheoretic understanding of the oscillatory states produced by this secondary bifurcation. Such an approach can be applied to analyse transitions in a wide variety of other physical systems, ranging from flows driven by differentially rotating boundaries ?[]Nore to Bose–Einstein condensation ?[]Huepe.
2 Method
2.1 Governing equations
We consider a fluid confined in a cylinder of depth and radius (figure 1). The aspect ratio is defined as . The fluid has kinematic viscosity , density , thermal diffusivity and thermal expansion coefficient (at constant pressure) . The top and bottom temperatures of the cylinder are kept constant, at and , respectively, leading to the linear conductive temperature profile . The lateral walls are insulating.
The Rayleigh number and the Prandtl number are defined by
{subeqnarray}
Ra & ≡& ΔT g γd3κν,
Pr & ≡& νκ .
Using the units , , and
for time, distance, velocity and temperature,
we define and to be the nondimensionalised velocity and
deviation of the temperature from the basic vertical profile, respectively.
We obtain the Boussinesq equations governing the system:
{subeqnarray}
Pr^1(∂_tu+
(u⋅∇)u)
&=&  ∇p + Δu+he_z
∂_t h
+(u⋅∇)h &=& Ra u_z+Δh
∇⋅u&=&0.
The boundary conditions for velocity are noslip and nopenetration
(1) 
Since the horizontal plates are assumed to be perfectly conducting (Dirichlet condition for ) and the vertical walls are insulating (Neumann condition), the boundary conditions for the temperature are {subequations}
(2)  
(3) 
2.2 Symmetries
Symmetries play an important role in the possible transitions undergone by this system. The Boussinesq equations (1) with boundary conditions (1)–(1) have reflection symmetry in the vertical direction , and rotational and reflection symmetry in the azimuthal direction . The reflection symmetry in is broken by the first bifurcation to a convective state. If the first convective state consists of axisymmetric convective rolls, then its remaining symmetries are reflection and rotation in , which together comprise the symmetry group . Bifurcations in the presence of symmetry were studied and classified during the 1980s by a large number of researchers, e.g. Bajaj (1982); Golubitsky & Stewart (1985); Knobloch (1986); van Gils & MalletParet (1986); Kuznetsov (1998); Coullet & Iooss (1990). We give a brief summary of their results.
First, the critical eigenvector may be axisymmetric. This case may be further subdivided according to whether the eigenvector is reflectionsymmetric or antisymmetric in and whether the eigenvalue is real or complex. A reflectionsymmetric eigenvector can lead to a target pattern of radially propagating rolls, e.g. Tuckerman & Barkley (1988). The breaking of reflection symmetry is associated with azimuthal flow.
Secondly, the critical eigenvector may be nonaxisymmetric. If the critical eigenvalue is real, then the resulting bifurcation is a circle pitchfork, leading to a “circle” of steady states parametrised by phase. Each steady state is reflection symmetric in (about some value ). If reflection symmetry is broken by a subsequent bifurcation, the scenario is that of a drift pitchfork, leading to slow motion (“drift”) along the circle. A complex eigenvalue corresponding to a nonaxisymmetric eigenvector, like that found by Wanschura et al. (1996) for parameters , , , leads to a Hopf bifurcation which engenders three nonlinear branches: standing waves, counterclockwise travelling waves, and clockwise travelling waves. The standing waves are reflectionsymmetric in (again about some value ), while the travelling waves break this symmetry. Our aim is to determine which of these types of waves is realised by our physical system.
2.3 Numerical integration
We integrated the equations by a classical pseudospectral method Gottlieb & Orszag (1977), in which each scalar of the fields u and is represented using Chebyshev polynomials in the radial and vertical direction and Fourier series in the azimuthal direction
(4) 
where the permitted combinations of are restricted by the parity and regularity conditions described in Tuckerman (1989) for , , and . The nonlinear (advective) terms were calculated in physical space and integrated via the Adams–Bashforth formula, while the linear (diffusive) terms were calculated in spectral space and integrated via the Crank–Nicolson formula. An influence matrix method was used to impose incompressibility Tuckerman (1989). A resolution of , , gridpoints or modes was found to be sufficient for nonlinear simulations. All computations were performed on the NEC SX5 vector supercomputer, with time step or , depending on , with CPU time per time step per grid point of .
2.4 Linear stability analysis
An important additional element in understanding the phenomena undergone
by the system is linear stability analysis. The procedure, which we summarise
below, is described in more detail in Mamun & Tuckerman (1995); Tuckerman & Barkley (2000) and references therein.
We linearise the equations about a steady state
:
{subeqnarray}
Pr^1(∂_tu+
(U⋅∇)u+(u⋅∇)U
)
&=&  ∇p + Δu+he_z
∂_t h
+(U⋅∇)h
+(u⋅∇)H
&=& Ra u_z+Δh
∇⋅u&=&0.
Equations (2.4) with boundary conditions
(1)(1) are then integrated in time in the same way
as the nonlinear equations (1).
We abbreviate the linear evolution problem (2.4) by
(5) 
Temporal integration is equivalent to carrying out the power method on the approximate exponential operator, since
(6) 
In order to extract the leading real or complex eigenvalues (those of largest real part) and corresponding eigenvectors, we postprocess the results of integrating (2.4) as follows. A small number of fields
(7) 
are calculated, by carrying out linearised timesteps. The Krylov space corresponding to initial vector and matrix is the dimensional linear subspace consisting of all linear combinations of vectors in (7). These vectors are orthonormalised to one another to generate a set of vectors which form a basis for the Krylov space. The action of the operator on the Krylov space is represented by a small () matrix whose elements are
(8) 
The small matrix can be directly diagonalised. Its eigenvalues approximate a small number of the eigenvalues of the large matrix : this is the essence of Arnoldi’s method. The procedure of generating the Krylov space via repeated action of selects preferentially the dominant values (those of largest magnitude) of , i.e. the leading eigenvalues (those of largest real part) of .
The eigenvectors of prescribe coefficients of the vectors which can be combined to form approximate eigenvectors of . The accuracy of these approximate eigenpairs is measured by the residue in the case of real eigenvalues or by the residues , in the case of complex eigenvalues. If the desired eigenvalues have sufficiently small residues, they are accepted; otherwise we continue integration of (2.4), replacing (7) by
(9) 
and so on, until the residue is below the acceptance criterion.
After integrating the axisymmetric version of the nonlinear equations (1) at a given Rayleigh number to create the nonlinear axisymmetric solution , we integrated the nonaxisymmetric linearised equations (2.4) to evolve from an arbitrary initial condition. To integrate (2.4), we used a timestep of and a spatial resolution of for each azimuthal mode. To construct the Krylov space (7) and approximate eigenpairs, we used vectors, a time interval of , and an acceptance criterion of .
2.5 Complex eigenvectors and their representations
The linear problem (2.4) for perturbations about an axisymmetric convective state can be divided into decoupled subproblems, each corresponding to a single azimuthal wavenumber . The problem for wavenumber can in turn be divided into two identical decoupled subproblems, corresponding to fields of the form {subequations}
(10)  
and  
(11) 
For simplicity, we will represent each of these types of vector fields by its temperature component and leave the dependence on and on to be written explicitly. We may write the linear evolution problem (5) restricted to fields with trigonometric dependence on such as (10)–(11) as
(12) 
A real eigenvalue breaking azimuthal symmetry in an symmetric situation is associated with a twodimensional eigenspace, consisting of linear combinations of vectors of type (10) and (11). Since {subequations}
(13) 
where
(14) 
all real eigenvectors have nodal lines and reflection symmetry about some . If we take and add (13) to the basic axisymmetric state, we obtain the “circle” of steady states resulting from a circle pitchfork mentioned in § 2.2.
A complex eigenvalue in the symmetric situation is associated with a fourdimensional eigenspace. Within each eigenvector class (10) and (11), the eigenspace is twodimensional, spanned by two linearly independent eigenvectors and , which are transformed by as
(15) 
In (15), can be replaced by any linear combination of and , but once is selected, the choice of follows from (15). Although the components of equation (15) are the real and imaginary parts of the complex equation
(16) 
the customary designation of and as the real and the imaginary part of the eigenvector is arbitrary, as reflected by the fact that an eigenvector can be multiplied by any complex number.
To form eigenvectors of the full cylindrical problem belonging to the fourdimensional eigenspace, each of and is multiplied by a trigonometric function. This yields as a basis for the fourdimensional eigenspace: {subequations}
(17)  
(18)  
(19)  
(20) 
One choice for a complex eigenvector pair is (17)(18), since
(21) 
More generally, the trigonometric dependence can be taken as in (13), with the same trigonometric dependence for each of and , to form a complex conjugate eigenvector pair each of whose members has nodal lines and axes of reflection symmetry, including . The evolution in time under (12) for a field with an initial condition of this form is
(22) 
The subspace of fields with azimuthal dependence is invariant under linearised time evolution. (There also exists an invariant subspace under the nonlinear time evolution, which includes harmonics , with the same axes of reflection symmetry.) If we take and in (22), and add this to the basic axisymmetric solution, then we obtain to first order the standing wave solution mentioned in § 2.2.
Any combination of (17)(20) is also a member of a complex eigenvector pair. The calculation
(23) 
when compared with (15), shows that the two components of the vector in (23) form a complex conjugate pair of eigenvectors for the full cylindrical problem, as in (15). Because and have different functional forms in , these vectors, unlike those of (13), cannot be combined into a single trigonometric function. Neither of the two components of (23) has nodal lines or reflection symmetry about any axis if both and are nonzero. The evolution in time under (12) for a field whose initial condition is the first component of (23) is
(24) 
If , then (24) becomes
(25) 
where or may be replaced by or . If we take and in (25) and add the basic axisymmetric solution, then we obtain, to first order, the expression for clockwise or counterclockwise travelling waves mentioned in § 2.2.
2.6 Amplitude equations and normal form
The linearised evolution treated in the previous section permits any combinations of (17)–(20). The mathematical analysis of Hopf bifurcation in the presence of symmetry carried out by e.g. Bajaj (1982); Golubitsky & Stewart (1985); Knobloch (1986); van Gils & MalletParet (1986); Kuznetsov (1998) describes the effect of including generic nonlinear terms compatible with the symmetries. Following the formulation of these authors, we decompose the field into a sum of clockwise and counterclockwise travelling waves with complex amplitudes and , respectively. The four variables form another description of the fourdimensional space described in the previous section. The nonlinear evolution of near the bifurcation can be described by the following amplitude equations or normal form: {subequations}
(26)  
(27) 
We use the normal form to interpret the results of our full numerical simulations.
Separating (2.6) into equations for real amplitudes and phases leads to {subequations}
(28)  
(29)  
(30)  
(31) 
Periodic solutions to (2.6) must be either standing or travelling waves. Solutions to (2.6) and their properties are given in Table 1. This table shows that both standing and travelling wave solutions exist for if and are both negative. A positive growth rate from a solution indicates instability. Thus, the stability of the solutions depends on the sign of : if , then standing waves are stable and travelling waves unstable, and vice versa for . Figure 2 shows phase portraits for the amplitudes , for the cases in which all three branches coexist and either the standing or the travelling waves are stable.
3 Results
3.1 Conductive state
Figure 3 shows the linear stability limits of the conductive state to perturbations with azimuthal wavenumbers , 1, and 2 Borońska & Boroński (2001). These results, obtained with the linearised version of our code, agree very closely with those presented by Wanschura et al. (1996). Note that in the range , the primary instability is axisymmetric. Immediately below and above this range of aspect ratio, the first instability is to an eigenvector with azimuthal wavenumber . Instability of the conductive state is independent of . However, the resulting nonlinear states and their stability depend on ; in the remainder of the study we fix .
3.2 Steady axisymmetric state
We reproduced the primary flow for and , parameters for which, according to Wanschura et al. (1996) and figure 3, the conductive state is unstable only to axisymmetric perturbations. In a fully threedimensional simulation, starting the evolution from an arbitrary nonaxisymmetric perturbation about the conductive state, we obtained a flow consisting of one toroidal roll. While axisymmetric, this flow breaks the reflection symmetry in and thus two such states exist, with either upflow or downflow at the centre; these are illustrated in figure 4. We used the state with downflow at the centre as the initial condition for higher Rayleigh numbers. According to the calculations of Wanschura et al. (1996), the axisymmetric state first bifurcates towards a flow with azimuthal wavenumber for and with wavenumber for . The critical Rayleigh numbers at which this loss of stability occurs are given in table 2.
(a)  (b) 
3.3 Eigenvalues and eigenvectors
Using the methods described in § 2.4, we integrated the evolution equations (2.4) linearised about axisymmetric solutions for aspect ratios and several different Rayleigh numbers. The leading eigenpairs calculated for , are given in Table 3. For these parameter values, the critical eigenvectors are (in order of decreasing growth rate): two conjugate pairs with azimuthal wavelengths and , a real eigenvector with , and another conjugate pair with .
Figure 5 represents the dependence of the leading eigenvalues on Rayleigh number for aspect ratios and , along with the azimuthal wavenumbers of the corresponding eigenvectors. was calculated by determining the zero crossing of , the growth rate of the leading eigenvalue (that of largest real part), by linear interpolation. (Critical Rayleigh numbers calculated by introducing perturbations into nonlinear simulations at various values of , and fitting the initial evolution to an exponential to calculate growth or decay rates gave similar results.) We then calculated , also by linear interpolation. The values we obtained for two aspect ratios and , and the corresponding values published by Wanschura et al. (1996) are those given in Table 2. The critical wavenumbers are the same, and the errors in and in are less than 1%. In what follows, we will focus on the instability, since the transition is similar; the aspect ratio is unless otherwise specified.
present study  Wanschura et al.  error  

0.76%  
42.33  42.54  0.48%  
3  3  
0.70%  
45.26  45.47  0.45%  
4  4 
eigenvalue  eigenvector visualisation  wavenumber  error  

real part  imaginary part  







–  




(a)  (b) 
(c)  (d) 
(a)  (b) 
(c)  (d)  (e)  (f)  (g) 
We summarise here the differences between our numerical method and that of Wanschura et al. (1996). We linearised a timestepping code in order to, in effect, carry out the power method (supplemented by an Arnoldi decomposition) on the exponential of the Jacobian. Wanschura et al. constructed the Jacobian matrix and used inverse iteration to compute its eigenvalues. Our calculation was restricted to one of the two identical decoupled subproblems, corresponding to only one of the invariant subspaces of the form (10) or (11). As a result, the complex eigenfunctions we show in table 3 are all in the eigenspace corresponding to standing waves, with three axes of reflection symmetry. Basis vectors for the remainder of the fourdimensional eigenspace can be found by rotating the eigenvectors of table 3, i.e. multiplying by instead of . Wanschura et al., in contrast, used the travelling wave form as an initial condition or invariant subspace, as discussed below.
In figure 6, we show representative elements of the eigenspace associated with the complex eigenvector at . Figures 6 (a, b) show and , while figures 6 (c–g) are generated via
(32) 
a form equivalent to (24) after translation of and of . Clockwise travelling waves ensue for (c), counterclockwise travelling waves for (e), and standing waves at different temporal phases for (d,f). Thus, the angle is similar to that used in figure 2. An eigenvector which corresponds to neither travelling nor standing waves is shown in figure 6 (g). These are all depicted on the slice ; when we plot the field of figure 6(c) at , we recover the form shown by Wanschura et al. We emphasise, however, that the other fields depicted in figure 6 are all equally valid eigenvectors. In particular, a nonlinear analysis, such as the simulations presented below, is required to determine whether the resulting nonlinear flow near onset is a travelling or a standing wave.
3.4 Weakly unstable standing waves
Above the critical Rayleigh number , a slightly perturbed axisymmetric state evolved in our simulations towards a threedimensional timedependent state, presented in figures 9, 9 and 9. Figure 9 shows temperature contours on the midplane at six regularly spaced instants in time within one oscillation period. In contrast to the eigenvectors depicted previously, figure 9 displays full nonlinear temperature fields, which are dominated by a large axisymmetric component. There are six pulsing extrema, engendering oscillation between two triangular structures of opposite phases (figures 9 a and 9 d). At each instant, the flow is invariant under rotation in by . In addition, this flow is also symmetric with respect to three different axes of reflection. Figure 9 shows contours of azimuthal velocity at the same times as figure 9. Figure 9 shows the temperature dependence on the angle for fixed radius and height at different times. Six fixed nodes identify this state as a standing wave with azimuthal wavelength .
The standing wave state persists for such a long time that it might seem stable. However, a small reflectionsymmetry breaking imperfection develops that eventually leads to the transition to travelling waves. Figure 10 shows the temperature dependence on the angle for the same parameters as figure 9, but at a later time. The breaking of reflection symmetry can be observed when the amplitude of the standing wave is small. The standing waves can be stabilised by imposing reflection symmetry. When we did this, above a threshold , we discovered a new (unstable) standingwave solution, displayed in figure 11 for .
In order to study the transition from standing to travelling waves, we monitored the growth of antisymmetric components. When the standing wave is still dominant, the amplitude of the antisymmetric components behaves in time like , where is the growth rate from standing waves to travelling waves. The growth rate , shown on figure 12 as a function of , is about two thirds of , the growth rate from the axisymmetric state to an flow (denoted in the previous sections by ). The observed lifetime of the standing waves decreases as the Rayleigh number is increased, since the growth rate increases.
(a)  (b)  (c)  (d)  (e)  (f) 
(a)  (b)  (c)  (d)  (e)  (f) 
(a)  (b)  (c)  (d)  (e)  (f) 
3.5 Stable travelling waves
After the pattern has evolved sufficiently from the standing wave state, the fixed antinodes abruptly begin to rotate about the cylinder axis. The six pulsing spots change into three rotating spots, as the standing waves become travelling waves with the same azimuthal wavelength. Figure 15, 15 and 15 depict temperature profiles and contours of the temperature and the azimuthal velocity of the travelling waves at different times. The travelling waves, like the standing waves, have threefold rotational symmetry, but do not have reflection symmetry.
Travelling waves are the final state of the time evolution. The reason for which we obtained standing waves before travelling waves in our simulations is that our initial conditions were reflection symmetric and our numerical procedures introduce antisymmetric perturbations at a low rate. (This is also seen in the simulations of thermocapillary flow by Leypoldt et al. (2000).) When the Rayleigh number is decreased, travelling waves persist until reaches .
We conducted simulations for several values of in the range and observed weakly unstable standing waves and stable travelling waves for all of them. We believe that the same scenario also occurs for , but with azimuthal wavenumber instead of .
(a)  (b)  (c)  (d)  (e)  (f) 
(a)  (b)  (c)  (d)  (e)  (f) 
3.6 Amplitudes and frequencies
We calculated the energy of both types of waves by first defining a norm whose square is
(33) 
where denotes spatial integration; (33) is one of many possible choices for this system. We then simulated the nonlinear evolution equations and calculated as the difference between the threedimensional and the axisymmetric solution. We define to be the integral of (33) over one oscillation period.
The energies , and frequencies , as a function of are shown in figure 16. The energies and frequencies for the two types of waves are quite close. The frequency obtained from linear stability analysis is also reproduced from figure 5 (b) for comparison. For both types of waves, the frequencies near the threshold are close to the Hopf frequency and the energy satisfies . These are hallmarks of a supercritical Hopf bifurcation.
(a)  (b) 
3.7 Normal form coefficients
Using the growth rates, amplitudes and frequencies of the standing and travelling waves that we have presented in sections 3.3 and 3.6, it is possible to calculate the coefficients of the normal form (2.6) for our particular case. The bifurcation parameter and frequency vary linearly with , while the other coefficients , , , are constants.
From the data in figures 5 (a,b), we extract the fits {subequations}
(34)  
(35) 
From the data in figure 16 we extract the fits {subequations}
(36)  
(37)  
(38)  
(39) 
Equations (3.7) are used to determine the nonlinear coefficients as {subequations}
(40)  
(41)  
(42)  
(43) 
An additional equation is provided by the data in figure 12 showing the growth rate from standing to travelling waves:
(44) 
and provides a second determination of
(45) 
which differs by 6% from (41).
4 Conclusion
We have used both nonlinear simulations and linear stability analysis to elucidate the behaviour of Rayleigh–Bénard convection in the parameter region of , first studied by Wanschura et al. (1996). In this regime, the primary axisymmetric convective state loses stability to an perturbation via a Hopf bifurcation whose critical eigenspace is fourdimensional. We calculated representative eigenvectors and explained how these relate to those computed by Wanschura et al. The bifurcation scenario guarantees that branches of standing waves and of travelling waves are created at the bifurcation, but that at most one of these branches is stable. Our nonlinear simulations showed a supercritical bifurcation leading to longlived standing waves which were eventually succeeded by travelling waves, both as time progressed and as the Rayleigh number was increased. We explained this by showing that the rate of transition from standing waves to travelling waves, while positive, is nevertheless small. In the absence of longtime integration and of these analyses, it would be easy to conclude that the standing waves were stable. This underlines the importance of calculating growth rates, in addition to carrying out nonlinear simulations, and of using established bifurcation scenarios to interpret physical phenomena.
The numerical and theoretical techniques we have used can be generally applied to study transitions in hydrodynamic problems. Our main tool was direct numerical simulation of the governing Boussinesq equations using a pseudospectral semiimplicit timestepping code. We complemented this approach with several other techniques. To carry out stability analysis, we first linearised the code. This requires very little modification of the existing code, but yields results which are far more precise and robust than restricting integration to the time interval during which perturbations to the basic state are small. Integrating the linearised equations is, in effect, an implementation of the power method for finding the fastest growing eigenvalues and corresponding eigenvectors. Eigenvectors with different azimuthal wavenumbers can be found simultaneously, since the linearised evolution of each Fourier mode is independent of the others. For a single wavenumber, this use of the power method is rendered more accurate and more general by postprocessing the results of linearised time integration with the Arnoldi decomposition to extract several, possibly complex, eigenvectors. We also interpreted our results in light of known results concerning axisymmetrybreaking Hopf bifurcations in systems with symmetry. This framework allows us to generate the fourdimensional eigenspace by combining eigenvectors with different symmetries. Traditionally, eigenvectors corresponding to clockwise and counterclockwise travelling waves are combined to form standing waves; we used a complementary, but equivalent, approach of combining standing waves of different spatial phases to form travelling waves. Finally, we interpreted our results in terms of the four ordinary differential equations comprising the normal form for Hopf bifurcations in systems with symmetry. Using our nonlinear simulations of the governing Boussinesq equations, we were able to calculate the various coefficients in the normal form equations.
We have not sought to determine the limits of the range of this phenomenon, in aspect ratio and Prandtl number. As these ranges were given by Wanschura et al. only for , a future direction would be to determine the whole zone in the parameter space where the Hopf bifurcation occurs. It would be interesting also to examine more closely the pulsing pattern found by Hof et al. (1999) at , , , in order to determine whether this state, evolving from axisymmetric flow, is the result of a bifurcation similar to that described in the present paper.
Acknowledgements.
The computations were performed on the NEC SX5 of the IDRIS (Institut du Développement et des Ressources en Informatique Scientifique) supercomputer center of the CNRS (Centre National pour la Recherche Scientifique) under project 1119.Footnotes
 thanks: email: kasia@limsi.fr – Web page: http://www.limsi.fr/Individu/kasia/en
 thanks: email: laurette@limsi.fr – Web page: http://www.limsi.fr/Individu/laurette
References
 Ahlers, G., Cannell, D. & Steinberg, V. 1985 Time Dependence of Flow Patterns near the Convective Threshold in a Cylindrical Container. Phys. Rev. Lett. 54, 1373–1376.
 Bajaj, A. K. 1982 Bifurcating periodic solutions in rotationally symmetric systems. SIAM J. Appl. Math. 42, 1078.
 Bodenschatz, E., Pesch, W. & Ahlers, G. 2000 Recent Developments in RayleighBénard convection. Annu. Rev. Fluid Mech. 32, 709–778.
 Borońska, K. & Boroński, P. 2001 Unpublished results.
 Buell, J. C. & Catton, I. 1983 The effect of wall conduction on the stability of a fluid in a right circular cylinder heated from below. Journal of Heat Transfer 105, 255.
 Busse, F. & Clever, R. 1979 Instabilities of convection rolls in a fluid of moderate Prandtl number. J. Fluid Mech. 91, 319–335.
 Charlson, G. S. & Sani, R. L. 1970 Thermoconvective instability in a bounded cylindrical fluid layer. Int. Journal. Heat Mass Transfer 13, 1479–96.
 Charlson, G. S. & Sani, R. L. 1971 On thermoconvective instability in a bounded cylindrical fluid layer. Int. Journal. Heat Mass Transfer 14, 2157–60.
 Charlson, G. S. & Sani, R. L. 1975 Finite amplitude axisymmetric thermoconvective flows in a bounded cylindrical layer of fluid. J. Fluid Mech. 71, 209.
 Ciliberto, S., Pampaloni, E. & PérezGarcía, C. 1988 Competition between Different Symmetries in Convective Patterns. Phys. Rev. Lett. 61, 1198–1201.
 Clever, R. & Busse, F. 1974 Transition to timedependent convection. J. Fluid Mech. 65, 625–645.
 Coullet, P. & Iooss, G. 1990 Instabilities of onedimensional cellular patterns. Phys. Rev. Lett. 64, 866.
 Croquette, V. 1989 Convective pattern dynamics at low Prandtl number: Part II. Cont. Phys. 30, 153–171.
 Croquette, V., Le Gal, P. & Pocheau, A. 1986 Spatial features of the transition to chaos in an extended system. Physica Scripta T13, 135.
 Croquette, V., Mory, M. & Schosseler, F. 1983 RayleighBénard convective structures in a cylindrical container. J. Phys. 44, 293–301.
 van Gils, S. A. & MalletParet, J. 1986 Hopf bifurcation and symmetry: travelling and standing waves on the circle. Proc. Roy. Soc. Edinburgh 104A, 279.
 Golubitsky, M. & Stewart, I. 1985 Hopf bifurcation in the presence of symmetry. Arch. Rat. Mech. Anal. 87, 107.
 Gottlieb, D. & Orszag, S. A. 1977 Numerical Analysis of Spectral Methods: Theory and Applications. Philadelphia: SIAM.
 Hardin, G. R. & Sani, R. L. 1993 Buoyancydriven instability in a vertical cylinder: Binary fluids with with Soret effect. Part 2: Weakly nonlinear solutions. Intl J. Numer. Meth. Fluids 17, 755.
 Hof, B., Lucas, G. J. & Mullin, T. 1999 Flow state multiplicity in convection. Phys. Fluids 11, 2815–2817.
 Huepe, C., Tuckerman, L. S., Métens, S. & Brachet, M. E. 2003 Stability and decay rates of nonisotropic attractive boseeinstein condensates. Phys. Rev. A 68, 023609.
 Knobloch, E. 1986 Oscillatory convection in binary mixtures. Phys. Rev. A 34, 1538.
 Koschmieder, E. L. 1993 Bénard cells and Taylor vortices. Cambridge University Press.
 Koschmieder, E. L. & Pallas, S. G. 1974 Heat transfer through a shallow, horizontal convecting fluid layer. Int. Journal. Heat Mass Transfer 17, 991–1002.
 Kuznetsov, Y. 1998 Elements of Applied Bifurcation Theory. Springer.
 Leong, S. S. 2002 Numerical study of RayleighBénard convection in a cylinder. Numerical Heat Transfer, Part A 41, 673–683.
 Leypoldt, J., Kuhlmann, H. & Rath, H. 2000 Threedimensional numerical simulation of thermocapillary flows in cylindrical liquid bridges. J. Fluid Mech. 414, 285–314.
 Mamun, C. K. & Tuckerman, L. S. 1995 Asymmetry and Hopf bifurcation in spherical Couette flow. Phys. Fluids 7, 80–91.
 Marqués, F., Net, M., Massaguer, J. M. & Mercader, I. 1993 Thermal convection in vertical cylinders. A method based on potentials of velocity. Comput. Method. Appl. M. 110, 157–169.
 Morris, S. W., Bodenschatz, Cannell, D. S. & Ahlers, G. 1993 Spiral Defect Chaos in Large Aspect Ratio RayleighBénard Convection. Phys. Rev. Lett. 71, 2026–2029.
 Müller, G., Neumann, G. & Weber, W. 1984 Natural convection in vertical Bridgeman configurations. J. Cryst. Growth 70, 78–93.
 Nore, C., Tuckerman, L. S., Daube, O. & Xin, S. 2003 The 1:2 mode interaction in exactly counterrotating von Kármán swirling flow. J. Fluid Mech. 477, 51–88.
 Plapp, B. B., Egolf, D. A., Bodenschatz, E. & Pesch, W. 1998 Dynamics and Selection of Giant Spirals in RayleighBénard Convection. Phys. Rev. Lett. 81, 5334–5337.
 Rosenblat, S. 1982 Thermal convection in a vertical circular cylinder. J. Fluid Mech. 122, 395–410.
 Rüdiger, S. & Feudel, F. 2000 Pattern formation in Rayleigh Bénard convection in a cylindrical container. Phys. Rev. E 62, 4927–4931.
 Sim, B.C. & Zebib, A. 2002 Effect of free surface heat loss and rotation on transition to oscillatory thermocapillary convection. Phys. Fluids 14, 225–231.
 Stork, K. & Müller, U. 1975 Convection in boxes: An experimental investigation in vertical cylinders and annuli. J. Fluid Mech. 71, 231–240.
 Touihri, R., Ben Hadid, H. & Henry, D. 1999 On the onset of convective instabilities in cylindrical cavities heated from below. I. Pure thermal case. Phys. Fluids 11, 2078–2088.
 Tuckerman, L. S. 1989 Divergencefree velocity fields in nonperiodic geometries. J. Comput. Phys. 80, 403–441.
 Tuckerman, L. S. & Barkley, D. 1988 Global bifurcation to travelling waves in axisymmetric convection. Phys. Rev. Lett. 61, 408–411.
 Tuckerman, L. S. & Barkley, D. 2000 Bifurcation analysis for timesteppers, in Numerical Methods for Bifurcation Problems and LargeScale Dynamical Systems. Springer, New York, ed. by E. Doedel and L. S. Tuckerman.
 Wanschura, M., Kuhlmann, H. C. & Rath, H. J. 1996 Threedimensional instability of axisymmetric buoyant convection in cylinders heated from below. J. Fluid Mech. 326, 399–415.