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

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.


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][3][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 poloidaltoroidal 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.

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 ℎ 2 /, the fluid layer half depth ℎ = /2, temperature difference between the rigid plates Δ, and the dimensionless form of the equations take the form: where 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Ω 2 ).
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 u = (, V, ) and the temperature , which represents the deviation from the conductive profile.
The no-slip boundary conditions at the upper and lower rigid walls are considered: at the rigid plates ( = ±1) in the -direction where e  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].

Considerations.
The assumption of periodicity allows the use of Fourier representation in the and -directions for the velocity and temperature fields in the form: ( (1)   V (1)   () +  (2)   V (2)   ()) , for the wave numbers where   =   /ℎ,   =   /ℎ are the scaled periods in and -directions, respectively, 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 using the orthogonality of Fourier expansion for each wavenumber pair (  ,   ), where g belongs to the physical space and f 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 V ()  (), the key consideration is the elimination of the pressure term from the resulting equations.This is accomplished by the conditions: where n 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 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.
Here, D 2 =  2 / 2 stands for a part of the Laplacian ∇ 2 .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, V (1)   , are lacking in their vertical velocity components while the second bases, V (2)   , are lacking in their vertical vorticity components.According to this classification, the expansion coefficient  (1)   is associated with the toroidal component and  (2)   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: Case 2.   = 0 and   = 0: where () =   (), () = (1 −  2 )  (), and ℎ() = (1− 2 ) 2   () 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 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.

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  (1) ,  (2) , and  in (19).The resulting linear system of ordinary differential equations reduces to a generalized eigenvalue problem for the eigenvalues  when a time dependence in the form [ (1) ;  (2) ; ] = E  (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.

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 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].

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.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.
Coriolis term depends linearly on the horizontal velocity.Hence, increasing rotation rate magnifies the toroidal energy (that is associated with V (1)   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.

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 ()  (),   () are the expansion coefficients, and {V ()  (),   ()} are the set of basis functions.The set of basis functions are chosen to satisfy

Figure 2 :
Figure 2: Marginal stability curves as effected by different rotation rates.

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