A Spectral Solenoidal-Galerkin Method for Rotating Thermal Convection between Rigid Plates

Abstract

The problem of thermal convection between rotating rigid plates under the influence of gravity is treated numerically. The approach uses solenoidal basis functions and their duals which are divergence free. The representation in terms of the solenoidal bases provides ease in the implementation by a reduction in the number of dependent variables and equations. A Galerkin procedure onto the dual solenoidal bases is utilized in order to reduce the governing system of partial differential equations to a system of ordinary differential equations for subsequent parametric study. The Galerkin procedure results in the elimination of the pressure and is facilitated by the use of Fourier-Legendre spectral representation. Numerical experiments on the linear stability of rotating thermal convection and nonlinear simulations are performed and satisfactorily compared with the literature.

1. Introduction

Thermal convection in a fluid layer has been a cradle of nonlinear hydrodynamic stability studies. In particular, the classical Rayleigh-Bénard problem of thermal convection in a horizontal layer which is heated from below is the most studied problem of the convective flows. This is due to its stability behavior exhibiting a sequence of discrete steps from steady regime to periodic, quasiperiodic regimes and eventually to chaotic regime as well as the simplicity of its geometry. This geometry of infinite fluid layer confined between rigid plates has been approximated by a periodic horizontal extent in the numerical studies and by large-aspect-ratio containers in the experiments. The underlying nonlinear stability mechanism of transition has been extensively studied both experimentally [1] and numerically [2–4]. Rayleigh and Prandtl numbers as well as the wavelength of the imposed disturbances play important roles as control parameters in the sequence of transitions. Other control mechanisms studied may include externally imposed magnetic field and rotation [5].

In this work, the effect of rotation on thermal convection is studied using a direct numerical simulation. Early studies on this topic go back to the theoretical investigations on the linear stability with/without rotation by Chandrasekhar [5]. Another detailed linear stability analysis of rotating Rayleigh-Bénard flow was performed by Clever and Busse [6] for critical wavenumber and corresponding critical Rayleigh number values for a range of Coriolis parameter (Ω) values in a geometry of rigid plates rotating around its vertical axis. Some quantitative stability criterion is obtained depending on the varied rotation, and a stability diagram of the natural convection with respect to rotation is constructed. It is observed that coriolis force reduces the heat transport for low Rayleigh and high Prandtl number flows, while limited rotation enhances the heat transport for decreasing Prandtl number.

Early nonlinear studies include Veronis’ work [7] who worked on the two-dimensional system with stress-free boundary conditions. It was observed that finite amplitude instability occurs due to the nonlinear effects and that it can be damped by increasing rotation. Further, Veronis [8] theoretically examined the effects of rotation and viscosity on the cellular motion of convection for different boundary conditions. A clear explanation was provided on the energy releasing and dissipative mechanisms with regards to viscosity. The effects of rotation on convection in low and high Prandtl fluids confined between stress-free boundaries were also studied by Veronis [9]. A unique feature is observed to occur that steady finite amplitude convection can exist for Rayleigh number which is lower than the critical Rayleigh number for a limited range of rotation.

Theoretical studies include those of Küppers and Lortz [10] that investigated the stability behavior for the case of infinite Prandtl number and free-free boundaries. It was shown that there is no stable steady-state convection for Taylor number which is higher than the critical Taylor number, Ta=2285(Ω≈23.9). Küppers [11] extended the work to rigid boundaries and finite Prandtl numbers. Recently Clever and Busse [12] studied numerically two- and three-dimensional convections under the influence of rotation. Unusual dynamical features for low Prandtl number flow were the main attention in this work. Kurt et al. [13, 14] studied rotating cylindrical annulus with small gap approximation. They conducted stability analysis of convection influenced by rotation and magnetic field. They also observed some instabilities.

Experimentally, Rossby [15] investigated the natural convection in several fluids confined between rotating and stationary plates. The predictions of the stability theory on the onset of convection were successfully tested. In this experiment, the finite amplitude instability in mercury for a limited range of rotation was studied. Somerville and Lipps [16] repeated Rossby's work in a three-dimensional numerical simulation. Quasi-steady and quasi-two-dimensional flows with the same experimental parameters of Rossby's were observed. Knobloch and Clune [17] focused on the previous experimental work and obtained the results of linear stability and weakly nonlinear calculations under the same experimental conditions.

In this work, a solenoidal spectral representation for the flow field is used in a Galerkin approach. As a consequence, the incompressibility (divergence free) and boundary conditions on the flow field are strictly enforced. For the sensitive nature of the parametric study undertaken, this facilitates the accurate determination of the transitory dynamic picture of the flow. The incompressibility condition appears as a constraint in the governing system of equations and is an important source of difficulty in numerical simulations. There are schemes developed solely to satisfy the continuity constraints such as the fractional step scheme by Orzag and Kells [18] and the scheme by Kleiser and Schumann [19]. However, this can only be achieved to a limited degree of accuracy. Another issue is the numerical handling of the pressure variable that usually comes without any boundary conditions. In the current approach, the pressure term is eliminated in the Galerkin projection procedure onto a solenoidal dual space.

There have been various works utilizing solenoidal spectral expansions. Moser et al. [20] presented a spectral method to automatically satisfy the continuity equation and boundary conditions and tested their method on the channel flow and the flow between concentric cylinders. They expanded the vertical and horizontal extents with Chebyshev polynomials and Fourier series, respectively. Kessler [21] studied steady and oscillatory regimes of Rayleigh-Bénard convection with explicitly constructed solenoidal bases based on poloidal-toroidal decomposition. Trigonometric polynomials and the beam functions were used in the construction of the solenoidal bases satisfying the boundary conditions in a rectangular container. Gelfgat [22] carried out a parametric study for Rayleigh-Bénard convection in rectangular 2D and 3D boxes with divergence-free Galerkin method based on Chebyshev polynomials of the first and second types. Clever and Busse [12] used toroidal-poloidal expansion in their numerical approach satisfying the solenoidal condition exactly; however, the procedure for eliminating the pressure leads to higher order differentials. Puigjaner et al. [23] studied numerically stability and bifurcation in convective flow in air in a cubical cavity heated from below. They used a divergence-free Galerkin spectral method to discretize the system and a parameter continuation method to determine the different branches of solution. They used combination of trigonometric and hyperbolic functions instead of Jacobi polynomial family. Most recently Meseguer and Trefethen [24] proposed a spectral Petrov-Galerkin formulation based on divergence-free bases in terms of Chebyshev polynomials to study stability of pipe flow.

In this work, the application of solenoidal spectral Galerkin approach to thermal convection under rotation utilizes Legendre polynomials as the underlying spectral representation. The use of Legendre polynomials avails Gauss-Legendre-Lobatto quadrature integration for the accurate evaluation of the inner product integrals resulting from Galerkin projection onto the dual space. The dual space is spanned by dual bases that are also required to satisfy the solenoidal condition in a form that incorporates the associated weight arising from the inner product in the spectral representation. This is required in order to eliminate the pressure term in the projection procedure. However, the use of Legendre polynomials significantly simplifies the construction of the solenoidal dual bases due to associated unity weight.

2. Governing Equations

Boussinesq equations are used as a model of buoyancy driven thermal convection with a coriolis term arising due to rotation in a physical domain between rotating rigid plates as illustrated in Figure 1. Under the scaling of the respective physical variables based on the thermal diffusion time ℎ²/𝛼, the fluid layer half depth ℎ=𝐻/2, temperature difference between the rigid plates Δ𝑇, and the dimensionless form of the equations take the form: ∇⋅𝐮=0,𝜕𝐮+()𝐮𝜕𝑡𝐮⋅∇=−∇𝑝+PrRa𝜃𝐞_𝑧+Pr∇²𝐮−2PrΩ𝐞_𝑧×𝐮,𝜕𝜃+()𝜕𝑡𝐮⋅∇𝜃=𝐮⋅𝐞_𝑧2+∇²𝜃, (1)(2)(3) where ΩΩ=^𝑧ℎ²𝜈,Ra=𝑔𝛽Δ𝑇ℎ³𝜈𝜈𝛼,Pr=𝛼 (4) are the Coriolis parameter (Ω) and Rayleigh (Ra) and Prandtl (Pr) numbers, respectively. The use of fluid layer half depth as the spatial scale, so that −1≤𝑧≤1 for convenience, results in Ra and Ω based on the half depth. The square of the Coriolis parameter is known as Taylor number (Ta=4Ω²). The quantities involved are Ω_𝑧 as the rotation rate about the vertical axis, 𝛼 as the thermal diffusivity, 𝛽 as the thermal expansion coefficient, 𝑔 as acceleration of gravity, and 𝜈 as the kinematic viscosity. The convective motions are described by the velocity 𝐮=(𝑢,𝑣,𝑤) and the temperature 𝜃, which represents the deviation from the conductive profile.

Figure 1 

Convective box influenced by rotation.

The no-slip boundary conditions at the upper and lower rigid walls are considered: 𝐮=𝟎,𝜃=0, (5) at the rigid plates (𝑧=±1) in the 𝑧-direction where 𝐞_𝑧 is the direction vector opposite to gravity. The flow in the infinite horizontal extent is assumed to be periodic in the 𝑥- and 𝑦-directions. Our investigation is limited to those cases with small rotation rate so that centrifugal force is negligible in comparison with gravity [6].

3. Numerical Procedure and the Bases

3.1. Considerations

The assumption of periodicity allows the use of Fourier representation in the 𝑥- and 𝑦-directions for the velocity and temperature fields in the form: 𝐮(𝐱,𝑡)=_{|𝑚|≤𝑁𝑥/2|𝑛|≤𝑁𝑦/2}𝑒^{𝑖𝑘𝑥𝑥+𝑖𝑘𝑦𝑦}×^𝑄_𝑝=0𝑎_𝑝⁽¹⁾𝐕_𝑝⁽¹⁾(𝑧)+𝑎_𝑝⁽²⁾𝐕_𝑝⁽²⁾(𝑧,𝜃())𝐱,𝑡=_{|𝑚|≤𝑁𝑥/2|𝑛|≤𝑁𝑦/2}𝑒^{𝑖𝑘𝑥𝑥+𝑖𝑘𝑦𝑦𝑄}_𝑝=0𝑏_𝑝𝑇_𝑝(𝑧), (6)(7) for the wave numbers 𝑘_𝑥=2𝜋𝑛𝑠_𝑥𝑁,−^𝑥2𝑁≤𝑛≤^𝑥2,𝑘_𝑦=2𝜋𝑚𝑠_𝑦𝑁,−^𝑦2𝑁≤𝑛≤^𝑦2, (8) where 𝑠_𝑥=𝐿_𝑥/ℎ, 𝑠_𝑦=𝐿_𝑦/ℎ are the scaled periods in 𝑥- and 𝑦-directions, respectively, 𝑎_𝑝^(𝑗)(𝑡),𝑏_𝑝(𝑡) are the expansion coefficients, and {𝐕_𝑝^(𝑗)(𝑧),𝑇_𝑝(𝑧)} are the set of basis functions. The set of basis functions are chosen to satisfy ∇⋅𝑒^{𝑖𝑘𝑥𝑥+𝑖𝑘𝑦𝑦}𝐕_𝑝^(𝑗)(𝑧)=0,𝐕_𝑝^(𝑗)(𝑧=∓1)=𝟎, (9) for 𝑗=1,2 and 𝑇_𝑝(𝑧=∓1)=0. The velocity basis functions 𝐕_𝑝^(𝑗)(𝑧) come in pairs due to the fact that the continuity equation (1) reduces the degree of freedom in determining the velocity field to two by connecting the three components.

Substitution of the representations (6) and (7) into the governing (2) and (3) results in a system of ordinary differential equations for the coefficients 𝑎_𝑝^(𝑗)(𝑡),𝑏_𝑝(𝑡) after a projection procedure which is carried out by the inner product (𝐟,𝐠)=¹₋₁𝜔(𝑧)𝐟^∗⋅𝐠𝑑𝑧, (10) using the orthogonality of Fourier expansion for each wavenumber pair (𝑘_𝑥,𝑘_𝑦), where 𝐠 belongs to the physical space and 𝐟 to the dual space, ∗ denotes the conjugate transpose, and 𝜔(𝑧) is a suitable weight function. The projection step requires construction of a set of dual basis functions {𝐕_𝑝^(𝑗)(𝑧),𝑇_𝑝(𝑧)}. In the construction of the dual velocity basis functions 𝐕_𝑝^(𝑗)(𝑧), the key consideration is the elimination of the pressure term from the resulting equations. This is accomplished by the conditions: ∇⋅𝜔(𝑧)𝑒^{𝑖𝑘𝑥𝑥+𝑖𝑘𝑦𝑦}𝐕_𝑝^(𝑗)(𝑧)=0,𝐕_𝑝^(𝑗)(𝑧=∓1)⋅𝐧=0, (11) where 𝐧 is the unit outward normal at the plates. The dual basis functions 𝑇_𝑝(𝑧) is required to satisfy 𝑇_𝑝(𝑧=±1)=0.

Finally, an underlying spectral representation for the basis functions based on a class of Jacobi polynomials, namely, Legendre polynomials, is utilized. There are two considerations in this choice. First, the inherent weight function 𝜔(𝑧)=1 associated with Legendre polynomials makes the construction of the solenoidal dual bases easier. Second, the availability of the Gauss-Legendre-Lobatto (GLL) quadrature facilitates accurate evaluation of the projection integrals (𝐟,𝐠)=¹₋₁𝐟^∗⋅𝐠𝑑𝑧≈^𝑁𝑧_𝑗=0𝜔_𝑗𝐟^∗𝑧_𝑗⋅𝐠^∗𝑧_𝑗, (12) for the integrands of order 2𝑁_𝑧−1. Further, the denser distribution of the GLL nodes near the boundaries of the plates provides the desired resolution for the boundary layers.

3.2. Construction

The solenoidal condition (1) in Fourier representation becomes 𝑖𝑘_𝑥𝑢(𝑚,𝑛,𝑧,𝑡)+𝑖𝑘_𝑦𝑣(𝑚,𝑛,𝑧,𝑡)+𝔻𝑤(𝑚,𝑛,𝑧,𝑡)=0. (13) Here, 𝔻²=𝑑²/𝑑𝑧² stands for a part of the Laplacian ∇². This provides the main constraint for the construction of the solenoidal and dual bases. The degree of freedom in the representation of a solenoidal flow field is reduced to two as the three components of the velocity field are connected through (13). Solenoidal nature of the flow field also involves the usual toroidal-poloidal decomposition of the motion. Therefore, first bases, 𝐕_𝑝⁽¹⁾, are lacking in their vertical velocity components while the second bases, 𝐕_𝑝⁽²⁾, are lacking in their vertical vorticity components. According to this classification, the expansion coefficient 𝑎_𝑝⁽¹⁾ is associated with the toroidal component and 𝑎_𝑝⁽²⁾ is associated with the poloidal component of the velocity field. The subspace 𝑆 of the solenoidal flow field and the dual subspace 𝐷 as constrained by (9), (11), and (13) are constructed as follows.

Case 1. 𝑘_𝑥≠0 and/or 𝑘_𝑦≠0: −𝑘𝑆=span^𝑦𝑘_𝑥𝑔(𝑧)𝑔(𝑧)0,𝑖𝑘_𝑥𝔻ℎ(𝑧)𝑖𝑘_𝑦(𝑧)𝑘𝔻ℎ_2𝑥+𝑘²_𝑦ℎ(𝑧),−𝑘𝐷=span^𝑦𝑘_𝑥𝑓(𝑧)𝑓(𝑧)0,𝑖𝑘_𝑥𝔻𝑔(𝑧)𝑖𝑘_𝑦(𝑧)𝑘𝔻𝑔_2𝑥+𝑘²_𝑦𝑔(𝑧). (14)

Case 2. 𝑘_𝑥 = 0 and 𝑘_𝑦 = 0: (𝑧)00,0𝑔(𝑧)0,(𝑧)00,0𝑓(𝑧)0,𝑇𝑆=span−𝑔𝐷=span−𝑓_𝑝(𝑧)=𝑇_𝑝(𝑧)=𝑔(𝑧), (15) where 𝑓(𝑧)=𝐿_𝑝(𝑧), 𝑔(𝑧)=(1−𝑧²)𝐿_𝑝(𝑧), and ℎ(𝑧)=(1−𝑧²)²𝐿_𝑝(𝑧) are selected to satisfy the boundary conditions, and 𝐿_𝑝(𝑧) is the Legendre polynomial of order 𝑝. Besides satisfying the solenoidal and boundary conditions, these forms also satisfy ¹₋₁𝐕^(1)∗_𝑞⋅𝐕_𝑝⁽²⁾𝑑𝑧=¹₋₁𝐕^(2)∗_𝑞⋅𝐕_𝑝⁽¹⁾,𝑑𝑧=0,₁₋₁𝐕^(1)∗_𝑞⋅𝔻²𝐕_𝑝⁽²⁾𝑑𝑧=¹₋₁𝐕^(2)∗_𝑞⋅𝔻²𝐕_𝑝⁽¹⁾,𝑑𝑧=0. (16) that result in some desirable decoupling features in the resulting system of equations. It can be shown that, associated with these forms, the number of quadrature nodes 𝑁_𝑧 and the number of solenoidal bases 𝑄 should be related in the least by 𝑁_𝑧=𝑄+4 for the linear case and 𝑁_𝑧=3/2𝑄+5 for the nonlinear case to render the numerical quadrature exact in the inner product integrals.

4. Numerical Implementation

When a representation of the solenoidal flow field in the form of (6) and (7) is substituted into the governing equations (2) and (3), the residuals arise: 𝑅_𝐮=−𝜕𝐮−()𝜕𝑡𝐮⋅∇𝐮−∇𝑝+PrRa𝜃𝐞_𝑧+Pr∇²𝐮−2PrΩ𝐞_𝑧𝑅×𝐮,_𝜃=−𝜕𝜃−()𝜕𝑡𝐮⋅∇𝜃−𝐮⋅𝐞_𝑧2+∇²𝜃. (17) The projection of the residual onto the dual space spanned by 𝐕_𝑝⁽¹⁾, 𝐕_𝑝⁽²⁾, and 𝑇_𝑝 is annulled, 𝐕,𝑅_𝐮=0,𝑇,𝑅_𝜃=0, (18) in the Galerkin procedure to yield the weak form: 𝐕⁽¹⁾,𝐕⁽¹⁾000𝐕⁽²⁾,𝐕⁽²⁾000𝑇,_𝐀𝑇×𝑎̇⁽¹⁾𝑎̇⁽²⁾𝑏̇+𝑐⁽¹⁾𝑐⁽²⁾𝑑=_Pr^{𝐕(1),(∇2−2Ω𝑒𝑧)𝐕(1)𝐕(1),−2Ω𝑒𝑧𝐕(2)0𝐕(2),−2Ω𝑒𝑧𝐕(1)𝐕(2),(∇2−2Ω𝑒𝑧)𝐕(2)Ra𝐕(2),𝑇𝑒𝑧}_{0𝐕𝑇,(2)⋅𝑒𝑧2Pr𝑇,∇2𝑇Pr𝐁}×𝑎⁽¹⁾𝑎⁽²⁾𝑏, (19) where 𝑎⁽¹⁾, 𝑎⁽²⁾, and 𝑏 are the time dependent expansion coefficients, and 𝑐⁽¹⁾, 𝑐⁽²⁾, and 𝑑 are the projections of the nonlinear terms.

4.1. Linear Stability Analysis

The solenoidal bases and the projection procedure are tested on the linear stability of the conductive (no-motion) state leading to the critical values when the convective motion just sets in. This amounts to elimination of the terms 𝑐⁽¹⁾, 𝑐⁽²⁾, and 𝑑 in (19). The resulting linear system of ordinary differential equations reduces to a generalized eigenvalue problem 𝜆𝐀𝐄=𝐁𝐄, (20) for the eigenvalues 𝜆 when a time dependence in the form 𝑎⁽¹⁾;𝑎⁽²⁾;𝑏=𝐄𝑒^𝜆𝑡 (21) is assumed. As the rightmost eigenvalue in the complex plane crosses the imaginary axis as Rayleigh number is varied for a given wavenumber pair (𝑘_𝑥=2𝜋/𝑠_𝑥,𝑘_𝑦=0) or (𝑘_𝑥=0,𝑘_𝑦=2𝜋/𝑠_𝑦), the system becomes unstable to infinitesimal perturbations. The resulting critical Rayleigh number Ra_𝑐 and critical wavenumber 𝑘_𝑐 values form the marginal stability curve.

4.2. Nonlinear Regime

The computation of nonlinear terms consumes more time in Fourier-Legendre space than in real space. Thus, all nonlinear terms are computed in real space and then projected onto Fourier-Legendre space to obtain the projections of the nonlinear terms 𝑐^1,2=𝐕^(1,2),(𝐮⋅∇)𝐮,𝑑=𝑇,(𝐮⋅∇)𝜃. (22) The derivatives of velocity and temperature fields in horizontal directions are calculated basically by the Fast Fourier transform (FFT), while polynomial differentiation matrix based on GLL collocation points is used for the vertical direction. The details and the associated MATLAB files for the polynomial differentiation can be found in [25].

For time discretization, semi-implicit numerical integration is used. The nonlinear and driving (buoyancy) terms are integrated explicitly using the third-order Adams-Bashforth scheme, while diffusive terms are integrated implicitly by Adams-Moulton scheme [26].

5. Results

A Fortran code is developed for the implementation of the proposed numerical procedure. It is tested for possible bugs, and resolution independence is constantly monitored. A linear stability analysis is conducted first, and the results are compared with the literature. The critical wavenumber (𝑘_𝑐) and Rayleigh number (Ra_𝑐) values are listed in Table 1 for the rightmost eigenvalue just crossing the imaginary axis. Some most critical values, Ra^∗_𝑐 and 𝑘^∗_𝑐, corresponding to the least Ra_𝑐 value on each marginal stability curve as indicated by “+” markers in Figure 2, are shown in Table 2. These are obtained by the selection of (𝑛=0,𝑚=1) (or equivalently (𝑛=1,𝑚=0) due to degeneracy in the horizontal directions) as defined in (8). They are in agreement with the known linear stability results for the rotating Rayleigh-Bénard problem. In addition, the most critical wavenumber and the corresponding Rayleigh number variation with rotation, namely, 𝑘^∗_𝑐(Ω) and Ra^∗_𝑐(Ω), are presented in Figures 3 and 4, respectively. There is a good agreement between the computed critical values and the literature [5, 6, 27]. Linear stability analysis can be performed with high accuracy by means of the proposed method. While typical runs for the linear stability calculations take a few CPU seconds, nonlinear simulations require higher CPU time because of time integrations and Fourier-Legendre transforms to be performed in high spatial resolutions. For example, typical computing time is approximately 3 hours for integrating 100 nondimensional seconds a long with 𝑁_𝑥×𝑁_𝑦×𝑁_𝑧=16×16×20 resolution in 3.0 GHz CPU unit.

Table 1 

Critical wavenumber and rayleigh number values at which convection just sets in (Ω=0).

Table 2 

Comparison of the most critical values with Chandrasekhar's calculation [5] at 8×8×8 resolution.

Figure 2 

Marginal stability curves as effected by different rotation rates.

Figure 3 

The critical wavenumber (𝑘∗𝑐) variation with rotation: 𝑘∗𝑐(Ω).

Figure 4 

The critical Rayleigh number (Ra∗𝑐) variation with rotation: Ra∗𝑐(Ω).

Nonlinear simulations are performed for varying control parameters, and the Nusselt number is computed. It is observed in Figure 5(c) that limited rotation in moderate Prandtl fluids stabilizes the convection and suppresses the oscillations in high Rayleigh number flow. This phenomenon is also observed in [9]. Increased heat flux is caused by the stabilizing effect of limited rotation in the range Ω=0--30. Coriolis force can balance horizontal temperature gradients, and, hence, less potential energy is released by horizontal temperature gradients. On the other hand, increasing rotation rate destabilizes the system and introduces new oscillatory motions for Ra−Ra_𝑐=20000 [12]. Increased heat flux with limited rotation is also observed for low Rayleigh number flows in Figure 5((a) and (b)), but this increment is very low because convection is still two-dimensional at these parameter values. Heat flux starts decreasing at Ω=10 for Ra−Ra_𝑐=2000,Ω=30 for Ra−Ra_𝑐=10000, and Ω=30 for Ra−Ra_𝑐=20000 because horizontal velocities destabilize the system with increasing rotation. In [9], it is indicated that rotational constraint balances the nonlinear processes for up to Ta≤10^3.6(Ω≈31) for convection between free boundaries. More recent work conducted by [12] reported that rolls are unstable for Pr>1 beyond Ω≈27.

Figure 5 

Nusselt number variation with rotation, Pr=0.71, (a)Ra−Ra𝑐=2000, (b)Ra−Ra𝑐=10000, and (c)Ra−Ra𝑐=20000.

Coriolis term depends linearly on the horizontal velocity. Hence, increasing rotation rate magnifies the toroidal energy (that is associated with 𝐕_𝑝⁽¹⁾ component of the solenoidal flow field) as shown in Figure 6. For all Rayleigh number values, increasing rotation absorbs the poloidal energy and stimulates the toroidal energy. Increase in toroidal energy is very rapid for Ra−Ra_𝑐=2000 and Ra−Ra_𝑐=10000 because both cases correspond to steady two-dimensional convection regime in nonrotating system. Limited Coriolis force, for example Ω=10, does not affect the roll structure but introduces three-dimensional motions as it is seen in Figures 7 and 8. This motion appears in oblique angle to the roll direction as also observed by [8]. On the other hand, Ra−Ra_𝑐=20000 case exhibits a different behavior. Since it corresponds to periodic motion regime in nonrotating system, it contains small-scale motions and vertical vorticity component. The motion tends towards two-dimensional convection with increasing rotation and toroidal energy sharply decreases, while poloidal energy conversely increases for Ω≈15. It is caused by the stabilizing effect of limited rotation as explained above. This effect also increases the heat transport because of increasing poloidal component. As vertical shear increases, Ω≥15 for low Rayleigh and Ω≥30 for high Rayleigh numbers, toroidal motion dominates and brings unstable character.

(a)

(b)

(c)

Figure 6 

Toroidal and poloidal kinetic energy variation with rotation, Pr=0.71.

Figure 7 

v-w streamline (solid line) and temperature contour (dash line) on Y-Z plane at 𝑥=𝐿𝑥/4, Pr=0.71, Ra−Ra𝑐=2000, and Ω=10.

Figure 8 

u-v streamline (solid line) and temperature contour (dash line) on X-Y plane at𝑧=0.4861,Pr=0.71, andRa−Ra𝑐=2000, andΩ=10.

6. Conclusions

In this work, solenoidal spectral approach is used to study rotating Rayleigh-Bénard convection numerically for a range of parameter values. The expansion of the velocity field in terms of the solenoidal bases and the subsequent projection onto the solenoidal dual space provide the automatic satisfaction of the divergence-free condition and the elimination of the pressure. Removal of these common algorithmic difficulties helps to focus on the extraction of the dynamics hidden within the model dynamical system equations. The current formulation provides a robust numerical tool to seamlessly investigate the linear and nonlinear regimes in rotating Rayleigh-Bénard convection. Form of the resulting dynamical system facilitates the use of numerical tools from the dynamical system and bifurcation theory that is a natural extension of this work.

Acknowledgments

The authors are grateful to (The Scientific and Technological Research Council of Turkey) TÜBİTAK, for providing financial support through research project 109M435 and Dr. Baran Aydın, for his suggestions.