of 19

Standing and travelling waves in cylindrical Rayleigh Bénard convection

0 views
All materials on our website are shared by users. If you have any questions about copyright issues, please report us to resolve them. We are always happy to assist you.
Share
Description
J. Fluid Mech. (2006), vol. 559, pp c 2006 Cambridge University Press doi: /s Printed in the United Kingdom 279 Standing and travelling waves in cylindrical Rayleigh Bénard
Transcript
J. Fluid Mech. (2006), vol. 559, pp c 2006 Cambridge University Press doi: /s Printed in the United Kingdom 279 Standing and travelling waves in cylindrical Rayleigh Bénard convection By KATARZYNA BOROŃSKA AND LAURETTE S. TUCKERMAN Laboratoire d Informatique pour la Mécanique et les Sciences de l Ingénieur (LIMSI CNRS), BP 133, Orsay, France (Received 11 April 2005 and in revised form 23 December 2005) 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 time-dependent 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. 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 three-dimensional nonlinear high-resolution 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 Ra c, whose value depends on the aspect ratio Γ radius/height. Above Ra c, convective motions appear and form various roll structures. 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 rich variety of patterns was reported: Pan Am patterns (arches with several centres of curvature, see Ahlers, Cannell & Steinberg 1985), straight parallel rolls (Croquette, Le Gal & Pocheau 1986; Croquette 1989), concentric rolls (targets, see Koschmieder & Pallas 1974; Croquette, Mory & Schosseler 1983), oneand several-armed rotating spirals (Plapp et al. 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). We focus here 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), Stork&Müller (1975) and Buell & Catton (1983). Critical Rayleigh numbers Ra c are 280 K. Borońska and L. S. Tuckerman about 2000 for Γ 1, increasing steeply for lower Γ and decreasing asymptotically towards Ra c = 1708 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 (Ra c = 2545 for Γ = 1, decreasing for higher Γ ) and the corresponding number of rolls. They then generalized this analysis (Charlson & Sani 1971), including non-axisymmetric modes and predicting Ra c and corresponding critical azimuthal wavenumbers. Stork & Müller (1975) observed experimentally convective patterns 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 & Sani. Rosenblat (1982) investigated convective instabilities numerically for free-slip boundary conditions, using a severely truncated expansion in a small number of eigenmodes. He described 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 performing linear analysis for the aspect ratio range 0 Γ 4. 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 Marqués et al. (1993). The flow succeeding the conductive state is three-dimensional 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. Müller, Neumann & Weber (1984) investigated convective flows experimentally and theoretically. They observed axisymmetric flows for Γ =1 and non-axisymmetric flows for 0.1 Γ 0.5. 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 andra c2 = The most complete numerical study of secondary convective instabilities for moderate aspect ratio cylinders was performed by Wanschura, Kuhlmann & Rath (1996). For cylinders with insulating sidewalls and 0.9 Γ 1.57, the primary bifurcation to convection occurs at Ra c 2000 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 1.45 Γ 1.57 at Pr = 1, where they found oscillatory instabilities at Ra c towards flows with azimuthal wavenumbers m =3 and m =4. The primary aim of this paper is to provide a more detailed description of these bifurcations. 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. They also found a secondary bifurcation point Ra c2, at which the axisymmetric flow becomes unstable towards a two-roll flow and calculated Ra c2 for Γ =1 and 0 Pr 1. An interesting experimental study was carried out by Hof, Lucas & Mullin (1999). Varying the Rayleigh number through different sequences of values, for fixed parameters Γ =2.0 andpr =6.7, they obtained several different stable patterns for the Standing and travelling waves in cylindrical Rayleigh Bénard convection 281 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. Convective patterns were numerically investigated by Rüdiger & Feudel (2000) and by Leong (2002). Rüdiger & Feudel found stability ranges for multi-roll patterns, targets and spirals for Γ =4,Pr = 1. Leong observed several steady convective patterns for aspect ratios 2 and 4 and Prandtl number Pr = 7, all of which were stable in the range 6250 Ra , and calculated the heat transfer for each pattern. Convective systems often display oscillatory behaviour. In binary fluid or rotating convection, the primary bifurcation is usually to periodic states, while in Rayleigh Bénard convection, periodic behaviour occurs as a secondary bifurcation. The oscillatory and skew-varicose instabilities of long straight parallel rolls calculated in, for example, Clever & Busse (1974) and Busse & Clever (1979), are manifested as travelling waves along rolls and as periodic defect nucleation (Croquette et al. 1986; Croquette 1989; 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 behaviour 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 surface-tension 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, Kuhlmann & Rath (2000), 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 (buoyancy-driven) convection literature. We wished to study in detail the time-periodic non-axisymmetric 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 1.45 Γ 1.57 and Pr = 1. 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 bifurcation-theoretic 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 et al. 2003) to Bose Einstein condensation (Huepe et al. 2003). 2. Method 2.1. Governing equations We consider a fluid confined in a cylinder of depth d and radius R (figure 1). The aspect ratio is defined as Γ R/d. 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 T 0 T /2 and T 0 + T /2, respectively, leading to the linear conductive temperature profile 282 K. Borońska and L. S. Tuckerman z T 0 T 2 d 2 0 θ R r d 2 T 0 + T 2 Figure 1. Geometry and coordinate system. T (z)=t 0 z T /d. The lateral walls are insulating. The Rayleigh number Ra and the Prandtl number Pr are defined by T gγ d3 Ra, κν (2.1a) Pr ν κ. (2.1b) Using the units d 2 /κ, d, κ/d and νκ/γgd 3 for time, distance, velocity and temperature, we define u and h to be the non-dimensionalized velocity and deviation of the temperature from the basic vertical profile, respectively. We obtain the Boussinesq equations governing the system: Pr 1 ( t u +(u )u) = p + u + he z, t h + (u ) h = Ra u z + h, u =0. The boundary conditions for velocity are no-slip and no-penetration (2.2a) (2.2b) (2.2c) u =0 for r = Γ or z = ±1/2. (2.3) Since the horizontal plates are assumed to be perfectly conducting (Dirichlet condition for h) and the vertical walls are insulating (Neumann condition), the boundary conditions for the temperature are h =0 for z = ±1/2, (2.4a) h =0 for r = Γ. (2.4b) r 2.2. Symmetries Symmetries play an important role in the possible transitions undergone by this system. The Boussinesq equations (2.2) with boundary conditions (2.3) (2.4) have reflection symmetry in the vertical direction z, and rotational and reflection symmetry in the azimuthal direction θ. The reflection symmetry in z 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 O(2). Bifurcations in the presence of O(2) symmetry were studied and classified during the 1980s by a large number of workers (e.g. Bajaj 1982; Golubitsky & Stewart 1985; Knobloch 1986; van Gils & Standing and travelling waves in cylindrical Rayleigh Bénard convection 283 Mallet-Paret 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 reflection-symmetric 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 non-axisymmetric. If the critical eigenvalue is real, then the resulting bifurcation is a circle pitchfork, leading to a circle of steady states parameterized by phase. Each steady state is reflection symmetric in θ (about some value θ 0 ). 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 non-axisymmetric eigenvector, like that found by Wanschura et al. (1996) for parameters 1.45 Γ 1.57, Pr =1,Ra , leads to a Hopf bifurcation which engenders three nonlinear branches: standing waves, counterclockwise travelling waves, and clockwise travelling waves. The standing waves are reflection-symmetric in θ (again about some value θ 0 ), while the travelling waves break this symmetry. Our aim is to determine which of these types of waves is realized by our physical system Numerical integration We integrated the equations by a classical pseudospectral method (Gottlieb & Orszag 1977), in which each scalar f of the fields u and h is represented using Chebyshev polynomials in the radial and vertical direction and Fourier series in the azimuthal direction f (r, z, θ, t) = N r,n z,n θ j,k,m=0 ˆ f jkm (t)c j (r/γ )C k (2z)e imθ + c.c., (2.5) where the permitted combinations of (j,m) are restricted by the parity and regularity conditions described in Tuckerman (1989) for u r, u θ, u z and h. 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 implicit Euler formula. An influence matrix method was used to impose incompressibility (Tuckerman 1989). A resolution of N r +1 = 36, 2(N θ +1) = 80, N z + 1 = 18 gridpoints or modes was found to be sufficient for nonlinear simulations. All computations were performed on the NEC SX-5 vector supercomputer, with time step or , depending on Ra, with CPU time per time step per grid point of Linear stability analysis An important additional element in understanding the phenomena undergone by the system is linear stability analysis. The procedure, which we summarize below, is described in more detail in Mamun & Tuckerman (1995); Tuckerman & Barkley (2000) and references therein. We linearize the equations about a steady state (U,H): Pr 1 ( t u +(U )u +(u )U) = p + u + he z, t h +(U )h +(u )H = Ra u z + h, u =0. (2.6a) (2.6b) (2.6c) 284 K. Borońska and L. S. Tuckerman Equations (2.6) with boundary conditions (2.3) (2.4) are then integrated in time in the same way as the nonlinear equations (2.2). We abbreviate the linear evolution problem (2.6) by ( ) ( u u t = L. (2.7) h h) Temporal integration is equivalent to carrying out the power method on the approximate exponential operator, since ( ( ) u u (t + t) =e h) L t (t). (2.8) h 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.6) as follows. A small number of fields ( ( u u (0), (T ), h) h) ( ( u u (2T ),..., ((K 1)T ) (2.9) h) h) are calculated, by carrying out T/ t linearized timesteps. The Krylov space corresponding to initial vector (u,h) T and matrix e LT is the K-dimensional linear subspace consisting of all linear combinations of vectors in (2.9). These vectors are orthonormalized to one another to generate a set of vectors v 1,v 2,v 3,...,v K which form a basis for the Krylov space. The action of the operator on the Krylov space is represented by a small (K K) matrixm whose elements are M jk v j, e LT v k. (2.10) The small matrix M can be directly diagonalized. Its eigenvalues λ approximate a small number of the eigenvalues of the large matrix e LT : this is the essence of Arnoldi s method. The procedure of generating the Krylov space via repeated action of e LT selects preferentially the K dominant values (those of largest magnitude) of e LT,i.e.theK leading eigenvalues (those of largest real part) of L. The eigenvectors of M prescribe coefficients of the vectors v j which can be combined to form approximate eigenvectors φ of e LT. The accuracy of these approximate eigenpairs (λ,φ) is measured by the residue e LT φ λφ in the case of real eigenvalues or by the residues e LT φ R (λ R φ R λ I φ I ), e LT φ I (λ R φ I + λ I φ R ) in the case of complex eigenvalues. If the desired eigenvalues have sufficiently small residues, they are accepted; otherwise we continue integration of (2.6), replacing (2.9) by ( u h) (T ), ( u h) (2T ), ( u h) (3T ),..., ( u h) (KT) (2.11) and so on, until the residue is below the acceptance criterion. After integrating the axisymmetric version of the nonlinear equations (2.2) at a given Rayleigh number to create the nonlinear axisymmetric solution (U,H), we integrated the non-axisymmetric linearized equations (2.6) to evolve (u,h) from an arbitrary initial condition. To integrate (2.6), we used a time step of t =10 4 and a spatial resolution of N r =47,N z = 29 for each azimuthal mode. To construct the Krylov space (2.9) and approximate eigenpairs, we used K = 10 vectors, a time interval of T = 100 t =10 2, and an acceptance criterion of Complex eigenvectors and their representations The linear problem (2.6) for perturbations (u,h) about an axisymmetric convective state (U,H) can be divided into decoupled subproblems, each corresponding to Standing and travelling waves in cylindrical Rayleigh Bénard convection 285 a single azimuthal wavenumber m. The problem for wavenumber m can in turn be divided into two identical decoupled subproblems, corresponding to fields of the form and û r (r, z)cos(mθ), û θ (r, z) sin(mθ), û z (r, z)cos(mθ), ĥ(r, z)cos(mθ), (2.12a) û r (r, z) sin(mθ), û θ (r, z)cos(mθ), û z (r, z) sin(mθ), ĥ(r, z) sin(mθ). (2.12b) For simplicity, we will represent each of these types of vector fields by its temperature component ĥ(r, z) and leave the dependence on θ and on t to be written explicitly. We may write the linear evolution problem (2.7) restricted to fields with trigonometric dependence on mθ such as (2.12a) (2.12b) as t ĥ = ˆL m ĥ. (2.13) A real eigenvalue breaking azimuthal symmetry in an O(2) symmetric situation is associated with a two-dimensional eigenspace, consisting of linear combinations of vectors of type (2.12a) and (2.12b). Since α ĥ(r, z)cos(mθ)+β ĥ(r, z) sin(mθ) =C ĥ(r, z)cos(m(θ θ 0 )), (2.14a) where C = α 2 + β 2, mθ 0 =atan(β/α), (2.14b) all real eigenvectors have m nodal lines and reflection symmetry about some θ 0.Ifwe take C Ra Ra c2 and add (2.14a) 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 O(2) symmetric situation is associated with a four-dimensional eigenspace. Within each eigenvector class (2.12a) and (2.12b), the eigenspace is two-dimensional, spanned by two linearly independent eigenvectors ĥ R and
Related Search
We Need Your Support
Thank you for visiting our website and your interest in our free products and services. We are nonprofit website to share and download documents. To the running of this website, we need your help to support us.

Thanks to everyone for your continued support.

No, Thanks