This chapter contains a basic introduction to convective instabilities. We give a brief sketch of the main results for Rayleigh–Bénard convection in cylindrical containers and we focus on phenomena reported by Hof, Lucas & Mullin (1999) and Wanschura, Kuhlmann & Rath (1996). These two studies define the two sets of control parameters we used in our simulations.
An instability occurs in a dynamical system when a slight perturbation to a solution does not die, but grows, causing an evolution of the whole system towards another solution. As in the real world, slight perturbations always exist, the stable solution is that observed experimentally. Usually the stability of a given solution depends on several control parameters. This is the case of Rayleigh–Bénard convection, where the stability of the flow structure is determined by the properties of the fluid and its container.
Rayleigh–Bénard instability occurs when a fluid layer is heate d from below and cooled from above. A vertical temperature gradient is then created, with upper layers heavier than lower lay-ers. The situation in which the fluid remains motionless and i ts density stratified, is called the conductive state. This state can become unstable in the presence of gravity. If the balance be-tween the stabilising and destabilising forces is favourable to the latter, the onset of convection is observed. Hot zones begin to ascend and then are cooled as they approach the upper plate; cold zones descend, and are heated. Such cycles are at the origin of convective motions, which take the form of rolls (see figure 1.1). A rich variety of conve ctive flows has been reported: convective rolls, straight or curved, can be parallel or concentric, forming targets, spirals, and more complex combinations.
The Rayleigh–Bénard problem is one of several widely-studied convective systems, in which a density variation in the presence of gravity engenders fluid motions. A similar prob-lem is Bénard–Marangoni convection, where the rigid upper plate is replaced by a free surface; surface tension must be then taken into account. In convection in mixtures the instability is driven by a concentration gradient, either with or without the presence of a temperature gra-dient. In natural convection, the instability occurs in a cavity with differentially heated lateral walls. We should mention here also forced convection, where the flow influencing the temper-ature field is engendered mainly by forces other than buoyanc y. Finally, more complex models are investigated in industrial research.
Convection is an important phenomenon in nature; its comprehension is necessary in un-derstanding the dynamics of such large-scale systems as the atmosphere, the oceans or the interiors of stars and planets (Busse, 1978; Getling, 1998; Marcus, 2004). At smaller scales, the interest in studying convection comes from multiple applications, such as designing heating systems of buildings, controlling crystal growth or constructing cooling systems for electronic devices. Another motivation for studying Rayleigh–Bénard convection comes from the fact that this is a classic pattern-forming system. Formation of patterns and their dynamics are a com-mon challenge occurring in many nonlinear problems, not only in fluid mechanics, but also in condensed-matter physics, chemistry and biology. As patterns observed in different systems can be very similar, there exist attempts to find unifying concep ts (e.g. Cross & Hohenberg, 1993); results for one particular configuration can then elucidate the dynamics of another system or contribute to a generalising theory.
Rayleigh number and convection onset
Lord Rayleigh (1916) showed that the balance between effects favourable and unfavourable to convection onset can be described by a dimensionless number, called the Rayleigh number. For a fluid layer of depth d, thermal diffusivity, kinematic viscosity and thermal expansion co-efficient (at constant pressure), confined between two parallel plates of temperature differ ence T , this number is defined by gd3 Ra ≡ T, where g is gravitational acceleration. Lord Rayleigh found analytically that for plates of infinite extent, convection occurs above the critical number Rac = 1708. For finite systems, the critical Rayleigh number and the convective flow structure depend on th e form of the lateral boundaries.
Secondary flows and Prandtl number
The flow which appears at the convection onset is often referr ed to as the primary flow. As the Rayleigh number is increased, this flow can lose its stabil ity to a secondary flow. This transition is governed not only by the Rayleigh number and container form, but also by the Prandtl number, defined as the ratio of the kinematic viscosity.
The Prandtl number varies widely: it is of order of 10−2 for liquid metals, 1 for gases and up to 105 for oils. An extensive study of the stability of parallel rolls in function of their horizontal wavenumber, Rayleigh number and Prandtl number was carried out by Busse and co-workers (Busse, 1978). The stability zone in the three-dimensional parameter space Ra-Pr-k (where k is the wavenumber of the rolls) is therefore called the Busse balloon.
Higher Rayleigh numbers
As the Rayleigh number is increased, the convective patterns become more complicated and more irregular, their symmetries and periodicities being successively broken. If we continue to increase Ra, spatio-temporal chaos appears and finally, the flow becomes turbulent. In the present work we focus on lower Rayleigh numbers; the description of turbulent regimes is a separate branch in fluid dynamics research.
Since Bénard’s experiments in the beginning of the twentiethcentury (Bénard, 1901), there has been constant progress in experimental techniques, allowing researchers to vary largely the Prandtl number, to work with quasi-Boussinesq fluids, and to o bserve the patterns not only in fluids but also in gases (Croquette, 1989). The experiments pr ovide either two-dimensional views of patterns (shadowgraphy technique) or approximate velocity fields (Particle Image Ve-locimetry), and usually also the Nusselt number, which measures the global heat transport. A great variety of convective patterns has been observed, especially for large rectangular, cylin-drical and annular containers (Koschmieder, 1993; Croquette, 1989; Bodenschatz et al., 2000).
Theoretical analysis of convective instabilities, with purely mathematical tools, could be per-formed only for the simplest cases because of the complexity of the problem – the laws gov-erning the system are described by partial differential equations, which are, in general, three-dimensional, nonlinear and time-dependent. There are some common simplifications: one of them is treating the fluid as incompressible and using the Bous sinesq approximation. The latter assumes that the fluid properties do not vary with temperatur e, except for the density, which varies linearly in the buoyancy term. Sometimes free-slip boundary conditions are used, which substantially facilitates calculations (Rosenblat, 1982). Convective patterns are also modelled by several kinds of amplitude equations, with simpler nonlinearities and vertical dependence maximally reduced (Manneville, 1990).
Linear analysis is a useful method in stability problems. Its principle is to add a small perturbation to a given base solution, linearise the equations (retaining only terms linear in the perturbation), and checking whether the perturbation grows or decreases. If the perturbation decreases, then the base solution is stable. If, increasing a control parameter, we cross a limit above which the base solution is unstable, the initial perturbation grows, evolving into a form called the critical eigenvector. However, as this formalism does not take into account nonlinear effects, linear analysis does not provide quantitative results concerning the new stable solution.
The structure of the critical eigenvector can, at best (in the case of a supercritical bifurcation), give an idea of the symmetry of the new stable solution, which is the sum of the base solution and the eigenvector, modified by nonlinear effects.
A new chapter in investigating convective instabilities began with the appearance of comput-ers. At first, numerical computations were used in the last st ages of theoretical analysis, or to perform linear analysis for perturbations with different wavenumbers, then slightly nonlinear or simplified models were studied numerically. Finally, fast i ncrease of computer performance and development of numerical methods made feasible three-dimensional, nonlinear, high-resolution simulations of the governing equations.
The incompressibility and Boussinesq approaches are widely used in simulations. This lim-its the possible range of investigations; incompressible models work only for flow velocities far below the speed of sound and the Boussinesq approximation is adequate only for small temper-ature gradients. There are also approaches with free-slip boundary conditions (Siggers, 2003), but contrary to the incompressibility and Boussinesq approximation, these models usually give rather unrealistic results (see Wanschura et al., 1996).
Several methods for numerical integration of the governing equations are in use. In the fi-nite difference method, discrete analogs of derivatives are calculated using local values. This method, the simplest in implementation, is still used in some models (Leong, 2002). Finite vol-ume methods, where the equations are integrated over control volumes corresponding to grid-points, have good local conservation properties and are therefore used even at high Rayleigh numbers, where turbulent regimes appear. Finite element methods, with their solid mathemati-cal foundations and a great deal of generality, are appreciated by engineers, as they are suitable for complicated geometries. They are widely used, but rather expensive numerically. Finally spectral methods, based on projecting unknown quantities on series of functions, are the most precise for a fixed number of points, but convenient only for s imple geometries. In our simula-tions we used a spectral code; we describe the projection further, in § 3.1.
Features of the system
In the present work we focus on convective flows in cylindrica l containers with aspect ratios G ≡ radius/height between 1.5 and 2.0 and with lateral boundaries that are either perfectly con-ducting or insulating. The advantage of cylindrical containers in comparison with finite rectan-gular boxes is the azimuthal symmetry, which permits us to use trigonometric basis functions; The system is invariant under rotation around the vertical axis, which influences the symme-try of the flow patterns. The disadvantage of the cylindrical geometry is the singularity at the cylinder axis.
As the azimuthal direction is periodic, a Fourier decomposition is a convenient representation : f (r, z, q) = åm fˆ(r, z) exp(imq). Because of the influence of the form of lateral boundaries on the flow geometry, for small aspect ratios there is often only one Fourier mode exp(imq) visibly dominating in a given convective pattern. The corresponding wavenumber is then used to describe the periodicity of the flow structure.
Large aspect ratios
A summary covering the developments since the mid 1980s for convective systems with large aspect ratio ( ≫ 1) can be found in Bodenschatz, Pesch & Ahlers (2000). In such domains a great variety of patterns was reported: “Pan Am” patterns ( arches with several centres of curvature, see Ahlers, Cannell & Steinberg, 1985), straight parallel rolls (Croquette, 1989; Cro-quette, Le Gal & Pocheau, 1986), concentric rolls (targets, see Koschmieder & Pallas, 1974; Croquette, Mory & Schosseler, 1983), one- and several- armed rotating spirals (Plapp, Egolf, Bodenschatz & Pesch, 1998), targets with dislocated centre (Croquette, 1989), hexagonal cells (Ciliberto, Pampaloni & Pérez-García, 1988) and spiral-defect chaos (Morris et al., 1993). A large overview on convective phenomena observed experimentally before this time can also be found in Koschmieder (1993).
Small aspect ratios
We focus on cylinders with moderate aspect ratio ∼ 1. 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), Sto rk & Müller (1975), Buell & Catton (1983). Critical Rayleigh numbers Rac are about 2000 for ∼ 1, increasing steeply for lower and decreasing asymptotically towards Rac = 1708 for →.
Convective state stability
Charlson & Sani (1970) estimated by a numerical variational technique the onset of axisymmet-ric convection in cylinders of aspect ratios between 0.5 and 8, with insulating and conducting sidewalls. They found the critical Rayleigh numbers (Rac = 2545 for = 1, decreasing for higher) and the corresponding number of rolls. They then generalised this analysis (Charl-son & Sani, 1971), including non-axisymmetric modes and predicting Rac and corresponding critical azimuthal wavenumbers.
Stork & Müller (1975) observed experimentally convective p atterns in annuli and cylinders of aspect ratio 0.7 ≤ ≤ 3.2, 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 free-slip boundary conditions, using a severely truncated expansion in a small number of eigenmodes. He de-scribed non-axisymmetric 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 perform ing linear analysis for the aspect ratio range 0 < ≤ 4. They determined the critical Rayleigh number and azimuthal wavenum-ber 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. The flow succeeding the conductive state is three-dim ensional over large ranges of aspect ratios, contrary to the expectations of Koschmieder (1993).
The stability of the first convective state, depending on bot h 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. Müller, Neumann & Weber (1984) investigated convective flows experimentally and th eoretically for aspect ratios 0.1 ≤ ≤ 10. For = 1 and Prantl numbers 0.02 and 6.7, they observed that the axisymetric primaryflow was stable up to Ra ∼ 10Rac.
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 m = 2 for = 1, Pr = 6.7 and Rac2 = 2430.
A very complete numerical study of secondary convective instabilities for moderate aspect ratio cylinders and adiabatic sidewalls was performed by Wanschura, Kuhlmann & Rath (1996). They focused on aspect ratios 0.9 < < 1.57, where the primary flow is axisymmetric. They calculated the nonlinear two-dimensional solutions and conducted linear analysis for Prandtl numbers 0.02 and 1. They thus found Rac2 and the azimuthal wavenumbers of corresponding critical eigenvectors.
Touihri, Ben Hadid & Henry (1999) numerically investigated the stability of the conductive state for aspect ratios = 0.5 and = 1. They described the main critical modes and established a diagram of primary bifurcations, including unstable branches. For = 1 they also found a secondary bifurcation point Rac2, at which the axisymmetric flow becomes unstable towards a two-roll flow and calculated Rac2 for 0 < Pr < 1.
Wanschura et al. (1996), in their study, predicted the flows succeeding the pr imary axisymmetric state to be steady, except over a narrow aspect ratio range 1.45 ≤ ≤ 1.57 at Pr = 1. For this configuration they found a Hopf bifurcation. The correspond ing eigenvectors had azimuthal wavenumbers m = 3 and m = 4 (see table 1.1). Thus the secondary flows should be oscillat ory, and according to the theory of Hopf bifurcation in systems with this symmetry, should have the form of either standing or travelling waves. As the study was linear, it was not possible to verify which of these two solutions succeeds the axisymmetric flow.
Multiplicity of convective states
Far above the convection threshold, a description in which only one stable state exists, with its form depending only on Rayleigh number, is insufficient. For c ylinders with small aspect ratios this was shown by the experimental study carried out by Hof, Lucas & Mullin (1999). Varying the Rayleigh number through different sequences of values, for fixed parameters = 2.0 and Pr = 6.7, they obtained several different stable patterns for the same final Rayleigh number. For 14 200 they observed two, three and four parallel rolls, a “me rcedes” pattern with three spokes of ascending or descending fluid and even an axisymmetric sta te. They also reported a transition from an axisymmetric steady state towards azimuthal waves.
More recently, convective patterns were investigated by Rüd iger & Feudel (2000), who used a spectral simulation. For aspect ratio = 4, they found stability ranges of several patterns: parallel rolls, target and spirals. The stability zones had common sections – depending on the initial conditions, different stable flows were observed at the same Rayleigh number.
Leong (2002) used a finite difference method to simulate conv ective flows for Prandtl num-ber 7 in cylinders of aspect ratios 2 and 4 and adiabatic lateral walls. He observed several steady patterns: parallel rolls, three-spoke flow and axisymmetri c state, all of which were stable in the range 6250 ≤ Ra ≤ 37500. He also calculated the heat transfer for each pattern.
The first part of our numerical study of Rayleigh–Bénard convec tion was stimulated by the re-sults of Wanschura 5 Oscillatory states are common in binary fluid convection and in rotating convection, but not as much so in ordinary Rayleigh–Bénard convection, particularly in small aspect ratios. We wished to study in detail this configuratio n in which an instance of time-periodic Rayleigh–Bénard convection was suggested. Hence we have simulated numerically the loss of stability of the first convective axisymmetric so lution undergoing an oscillatory bi-furcation for 1.45 ≤ ≤ 1.57 and Pr = 1. In chapter 5 we describe the results of linear stability analysis and nonlinear simulations, which identify the scenario in terms of bifurcation theory in systems with symmetries.
We wished also to explore the multiplicity of stable states reported by Hof et al. (1999). This phenomenon is known in other hydrodynamic systems and in large cylinders, but rarely reported for Rayleigh-Bénard convection in cylinders of smal aspect ratio. We ran a series of nonlinear simulations for Prandtl number Pr = 6.7, aspect ratio = 2.0 and two types of thermal boundary conditions. Chapter 6 presents the results for perfectly insulating sidewalls and chapter 7 for perfectly conducting sidewalls.
The methods used in our study are described in the three chapters which follow. In chapter 2 we present the nonlinear equations governing the system and we explain the principles of linear analysis. We define also the symmetries characterisi ng the problem and introduce and derive an appropriate representation describing linear and nonlinear dynamics of the system. Chapter 3 describes the most important aspects of the numerical methods used for integrating the differential equations. This comprises discretisation in space and time, imposing boundary conditions and finding the leading eigenvectors and eigenva lues for linearised equations. We also give some details concerning the optimisation of the code and post-processing of the re-sults. The choice of the simulation parameters and some of the tests performed on the code are discussed in chapter 4.
Table of contents :
1 Rayleigh–Bénard instability
1.1 Rayleigh–Bénard instability
1.1.2 Rayleigh–Bénard convection
1.1.3 Rayleigh number and convection onset
1.1.4 Secondary flows and Prandtl number
1.1.5 Higher Rayleigh numbers
1.1.6 Investigation tools
1.2 Cylindrical system
1.2.1 Features of the system
1.2.2 Azimuthal wavenumber
1.2.3 Large aspect ratios
1.2.4 Small aspect ratios
2 Governing equations
2.1 Navier–Stokes equations
2.2 Boundary conditions
2.3 Linear analysis
2.3.1 Linearised equations
2.3.2 Power method
2.3.3 Complex eigenvalues
2.4.1 Symmetries of Boussinesq equations
2.4.2 Z2 symmetry and inverse patterns
2.4.3 O(2) symmetry
2.4.4 Representations of complex eigenvectors in O(2) symmetry
2.5 Amplitude equations and normal form
3 Numerical integration
3.1 Spectral discretisation – Galerkin projection
3.1.1 Galerkin decomposition
3.1.2 Fourier series in q
3.1.3 Chebyshev polynomials in z
3.1.4 Chebyshev polynomials in r
3.2 Time discretisation
3.3 Influence matrix
3.3.1 Tau method
3.3.2 Velocity–pressure decoupling
3.3.3 Influence matrix
3.3.4 Tau correction
3.4 Linear simulation
3.4.1 Power method
3.4.2 Arnoldi–Krylov method
3.4.3 Calculating the base state
3.5 Code vectorisation
3.6 Program output
4 Code validation
4.1 Choice of resolution
4.2 Choice of timestep
4.3 Choice of parameters for linear runs
4.4 Convergence criteria
4.5 Test cases
5 Standing and travelling waves
5.1 Oscillatory bifurcation of the axisymmetric state — linear analysis
5.1.1 Steady axisymmetric state
5.1.2 Eigenvalues and eigenvectors
5.2 Nonlinear simulation of time-dependent states
5.2.1 Weakly unstable standing waves
5.2.2 Stable travelling waves
5.2.3 Amplitudes and frequencies
5.2.4 Normal form coefficients
6 Convective patterns – insulating sidewalls
6.1 Simulation of the experiment
6.2 Evolution from perturbed conductive state
6.3 Three-roll patterns
6.4 Evolution from four-rolls
6.5 Evolution from pizza pattern
6.6 Evolution from two rolls
6.7 Axisymmetric flows
6.8 Evolution from mercedes pattern
6.9 Evolution from dipole pattern
6.10 Dipole-shaped perturbation
6.11 Summary diagram
7 Convective patterns – conducting sidewalls
7.1 Start from perturbed conductive state
7.2 Three rolls
7.3 Four rolls
7.4 Evolution from dipole pattern
7.5 Evolution from hourglass pattern
7.6 Evolution from star pattern
7.7 Evolution from da Vinci pattern
7.8 Evolution from Y pattern
7.9 Summary diagram