2.1. Equations, computational domain and scaling

The computation domain $\varOmega$ in this study is a cylinder with height $H$ and diameter $D$. We can express $\varOmega$ in cylindrical coordinates that are normalized by $H$ and symmetrized about the mid-plane $(H/2 \rightarrow z=0)$ such that the normalized $\varOmega$ is defined as

(2.1)\begin{equation} \{\varOmega(r,\theta,z) \,|\, r\in[0,\varGamma/2],\ \theta \in [0,2{\rm \pi}), \ z \in[-0.5,0.5]\}, \end{equation}

where $\varGamma$ is the aspect ratio $(\varGamma =D/H)$ of the cylinder. Here $\varOmega$ is also aligned with the gravitational vector $(\boldsymbol {g})$ such that ${\boldsymbol {g}}/{|\boldsymbol {g}|}=-\hat {e}_z$, where $\hat {e}_z$ is the unit normal in the $z$-direction.

Velocity and temporal units are normalized by the ‘free-fall’ velocity $(w_f=\sqrt {\beta g {\rm \Delta} T H})$ and time $(t_f=H/w_f)$, where $\beta$ is the coefficient of thermal expansion, $g$ is the gravitational constant and ${\rm \Delta} T$ is the temperature difference between the top and bottom plates of the convection cell. Here ${\rm \Delta} T$ is used to normalize the temperature field and the Boussinesq approximation is applied to the incompressible Navier–Stokes equations for computation and analysis. The reference temperature for the Boussinesq approximation is taken to be the average mid-plane temperature such that

(2.2)\begin{equation} \vartheta=\frac{T-T_{ref}}{{\rm \Delta} T} \in [-0.5,0.5]. \end{equation}

Utilizing these scales, the non-dimensional form of the Boussinesq equations for RBC can be expressed as

(2.3)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}

(2.4)\begin{gather}\frac{\partial{\boldsymbol{u}}}{\partial t}+(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u}=-\boldsymbol{\nabla} p + \sqrt{\frac{Pr}{Ra}} \nabla^2 \boldsymbol{u} + \vartheta \hat{e}_z, \end{gather}

(2.5)\begin{gather}\frac{\partial \vartheta}{\partial t}+(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \vartheta=\frac{1}{\sqrt{Ra\, Pr}} \nabla^2 \vartheta, \end{gather}

where $\boldsymbol {u}, p$ and $\vartheta$ are the dimensionless velocity, pressure and temperature. The Rayleigh $(Ra)$ and Prandtl $(Pr)$ numbers in (2.4) and (2.5) are defined as

(2.6)\begin{equation} Ra=\frac{\beta g {\rm \Delta} T {{H}}^3}{\alpha \nu}, \end{equation}

(2.7)\begin{equation}Pr=\frac{\nu}{\alpha}, \end{equation}

where $\alpha$ and $\nu$ are the thermal diffusivity and kinematic viscosity, respectively.

In this study $\varGamma =6.3$, $Ra=9.6\times 10^7$, $Pr=6.7$. When presenting the results, the velocity field is expressed in terms of cylindrical coordinates such that $\boldsymbol {u}=\{u_r,u_\theta ,u_z\}$ and $u_i$ represents any of the components $\{u_r,u_\theta ,u_z\}$.

Two primary time scales are used in this work. The free-fall time scale $(t_f)$, and an eddy turnover $(t_\epsilon )$. Here $t_\epsilon$ approximates the mean up and down times for large-scale motions that span the depth of the domain. It is defined as

(2.8)\begin{equation} t_\epsilon = \frac{2H}{\sqrt{\langle u_z^2\rangle_{V, t}}}, \end{equation}

where $\langle \rangle _V$ indicates a spatial average over the entire domain volume $V$ and $\langle \rangle _t$ is an average in time.

Eddy-turnover time can be thought of as the average time it will take for a particle to cross the layer depth twice driven by a turbulent diffusion. For reference, $t_\epsilon$ is approximately $31t_f$ in this study.

2.2. Numerical method

The data in this study is obtained by DNS using the open source spectral element code Nek5000. Nek5000 is a thoroughly validated research code that has been used extensively in scientific literature (see Fischer, Lottes & Kerkemeier Reference Fischer, Lottes and Kerkemeier2008). The spectral element method uses high-order polynomial approximation of the solution within each element, while local element solutions are assembled globally through gather-scatter operations (Deville, Fischer & Mund Reference Deville, Fischer and Mund2002). A $P_n-P_{n-2}$ (Fischer Reference Fischer1997) spectral element formulation is employed herein. Temporal integration is of second order, utilizing an implicit backward-differentiation formula for the viscous terms, and an explicit extrapolation for nonlinear and Boussinesque terms. Pressure decoupling is performed using the operator splitting method (Karniadakis, Israeli & Orszag Reference Karniadakis, Israeli and Orszag1991), and a resulting pressure Poisson equation is solved by GMRES (Saad & Schultz Reference Saad and Schultz1986) with a multigrid preconditioning. The computational grid is discretized with hexahedral elements and a marginal amount of biasing toward the upper and lower plates is applied to the element distribution. The spectral element method (SEM) used in this simulation also applies a Gauss–Lobatto–Legendre (GLL) quadrature within each element which clusters points toward the boundaries of each element and greatly improves resolution at the walls. Ninth-order polynomials were used for the quadrature resulting in roughly 44.6 million grid points. A schematic of the computational domain and the spectral element grid employed is provided in figure 1(a). In the experimental work of Fernandes (Reference Fernandes2001) it is reported that the Kolmogorov length for RBC at this $\varGamma$, $Ra$ and $Pr$ is approximately $1.2\times 10^{-2}H$, and our simulation's grid had five points within this range at the wall. We also determined that this grid satisfies the spatial resolution criteria of Grötzbach (Reference Grötzbach1983). The temporal resolution for each time step was approximately $t_\epsilon \times 10^{-4}$ with a corresponding Courant–Friedrichs–Lewy range of ${\sim }0.6\text {--}0.7$. In our prior work we conducted a a-posteriori analysis to evaluate the resolution of our results utilizing the techniques outlined by Scheel, Emran & Schumacher (Reference Scheel, Emran and Schumacher2013) and confirmed that our simulation met the requirements for dissipation continuity across elements, as well as the resolution with respect to the height dependent Kolmogorov and Batchelor scales (see Sakievich et al. Reference Sakievich, Peet and Adrian2016). Boundary conditions are specified as follows: no-slip conditions for velocity are set at all bottom, top and side walls. Temperature is set to a constant value $T=T_{hot}$ at the bottom wall, and $T=T_{cold}$ at the top wall, with the temperature difference ${\rm \Delta} T=T_{hot}-T_{cold}$ used for non-dimensionalization in (2.2) and (2.6), while adiabatic boundary conditions are used for the temperature at the side walls. Additional details regarding resolution, convergence and comparison with experiments for the specific computations in this work can be found in the prior work of Sakievich et al. (Reference Sakievich, Peet and Adrian2016).

Figure 1. Images of the computational grid using spectral elements (44.6 million grid points with a ninth-order polynomial basis in the spectral elements) (a) and an example of the cylindrical coordinates post-processing grid with a resolution of [60, 512, 32] points in the $r$, $\theta$ and $z$ directions, respectively (b). The sampling resolution in (b) is reduced from the actual resolution used for analysis [160, 2048, 64] to improve image quality in the figure.

The total run time of the simulations in the current work is 3054 free-fall time units, or close to 100 eddy-turnover times, which makes it one of the longest DNS studies of turbulent Rayleigh–Bénard convection up to date (Sakievich et al. Reference Sakievich, Peet and Adrian2016; Pandey et al. Reference Pandey, Scheel and Schumacher2018), and perhaps the longest for the considered Rayleigh number. Note that while the current study undoubtedly pushes the limit in terms of the simulation time, it still only covers less that one viscous time unit $t_v=\sqrt {Ra/Pr}\, t_f=3800t_f$, and only $10\,\%$ of a thermal diffusion time unit $t_d=\sqrt {Ra\, Pr}\,t_f=25\,460t_f$ for the current $Pr$ number, once again highlighting the challenges of accessing the longest possible time scales in the RBC studies at high $Ra$ number. While the total span of the simulations is more than enough to capture the superstructures (Pandey et al. Reference Pandey, Scheel and Schumacher2018) and their dominant dynamics, it might not be enough to capture the effect of slow diffusive processes on their evolution.

The data for calculating statistical quantities was sampled every three free-fall times. This choice of a sampling rate is a compromise between fine resolution in time and a storage requirement for a very long temporal study, such as the one performed in the current work. The sampling rate is more than adequate for studying the long time scales and slow energetic processes associated with the evolution of the large-scale structures, which is the focus of the current work. For post-processing, each snapshot is projected onto cylindrical coordinates using spectral interpolation routines native to Nek5000, and the velocity components are transformed from Cartesian to cylindrical. This was previously done on a smaller scale (Sakievich et al. Reference Sakievich, Peet and Adrian2016), but in this work the transformation has been extended to the entire domain. Cylindrical coordinates are the logical choice for analysing the current dataset and facilitate operations along the domain's periodic, azimuthal direction. The DNS snapshots are resampled with [160, 2048, 64] points in $r,\theta$ and $z$, respectively, to generate the cylindrical grids used for analysis. Non-uniform, Gauss–Legendre (GL) quadrature is used to sample in the $r$ and $z$ directions, but the $\theta$ direction uses equispaced sampling points to facilitate Fourier transforms. Gauss–Legendre quadrature does not include the endpoints and is defined on the standard interval $x\in (-1,1)$. Gauss–Legendre quadrature is selected to facilitate high accuracy numerical integration and to remove unnecessary sampling at the walls of the cell where the system is well defined. Boundaries in the $z$ direction are constrained with Dirichlet boundary conditions, so that sampling on them for Fourier transforms is trivial. Points along the central axis $(r=0)$ are at a spatial singularity in the cylindrical coordinate representation and provide no additional data when Fourier transforms in $\theta$ and integration over the $r$–$z$ plane are applied. The points along $r=\varGamma /2$ have Neumann boundary conditions in the temperature field but virtually no information is lost since the gradient at the wall is zero (adiabatic), and the GL quadrature samples very close to the boundaries. An example image of a post-processing grid in cylindrical coordinates is shown in figure 1(b).

The post-processing grid has a little less than half the number of grid points when compared with the computational grid, but the change in coordinate system and quadratures leads to non-uniform sampling ratios with respect to the original domain $\varOmega$. In the vertical direction this causes the post-processing grid to have approximately twice as many points in the boundary layer region, and about half as many in the bulk region. The horizontal sampling is not as straightforward to compare because the coordinate systems are fundamentally different. However, it is definitely the case that the centre of the post-processing grid is sampled more finely than the computational grid. The spectral interpolation used during resampling ensured that no wavenumber information was lost on a post-processing grid according to its resolution, but neither is gained, compared to the original computational grid. Figure 1 provides an example of the two different grids (computational and post-processing), however, the number of depicted grid points is reduced in figure 1(b) with respect to the actual post-processing grid used for analysis to make the visualization comprehensible. The spectral element grid in figure 1(a) does not require a reduction in sampling for visualization purposes because the GLL points within each element can be given a different contrast.

2.3. Fourier decomposition

Fourier decomposition in the azimuthal direction renders insights into the structure of the flow field. Fourier modes are a natural choice because the azimuthal direction is geometrically periodic. Fourier decomposition provides additional benefits in this study that extend beyond the mathematical significance of the modes. For example, azimuthal motions for RBC in cylinders tend to evolve on extremely long time scales, and the azimuthal velocity signals are relatively weak (Brown et al. Reference Brown, Nikolaenko and Ahlers2005; Mishra et al. Reference Mishra, De, Verma and Eswaran2011). Performing an analytical decomposition such as Fourier analysis allows the azimuthal evolution of the flow to be studied in a well understood format with precise measurements.

For a general signal $u(r,\theta ,z,t)$, azimuthally periodic with a period $2{\rm \pi}$, the Fourier series decomposition is given by

(2.9)\begin{equation} u(r,\theta,z,t)=\sum_{k=0}^{\infty} \hat{u}_k(r,z,t)\,\textrm{e}^{jk\theta}, \end{equation}

where $j=\sqrt {-1}$ and $k$ is the integer Fourier mode number. Fourier coefficients are given by the inverse operation

(2.10)\begin{equation} \hat{u}_k(r,z,t)=\frac{1}{2{\rm \pi}}\int_{0}^{2{\rm \pi}} u(r,\theta,z,t) \,\textrm{e}^{-jk\theta}\,\textrm{d}\theta. \end{equation}

The Fourier series representation can be approximated by a finite number of modes, $N$, and a uniform convergence of a truncated series to the original signal $u(x)$ is guaranteed given the signal is (a) smooth, and (b) periodic. A discrete Fourier transform in this case can be defined as

(2.11)\begin{equation} \hat{u}_k(r,z,t)=\frac{1}{N}\sum_{i=0}^{N} u(r,\theta_i,z,t) \,\textrm{e}^{-j k \theta_i}. \end{equation}

Note that in this formulation, $k$ can be understood as the integer azimuthal mode number, and, at each particular radius $r$, the angular wavenumber $\tilde {k}$ could be defined as

(2.12)\begin{equation} \tilde{k}=\frac{2{\rm \pi}}{2{\rm \pi} r}k=\frac{k}{r}, \end{equation}

so that the wavelength, that gives a relation to a physical length of the structures, for each particular azimuthal mode $k$ at each particular radius is given as

(2.13)\begin{equation} \lambda_k(r)=\frac{2{\rm \pi}}{\tilde{k}}=\frac{2{\rm \pi} r}{k}, \end{equation}

where the dependence on $r$ denotes that the wavelength for a certain $k$ mode has been calculated at a given radius.

Throughout this work, Fourier coefficients are indicated by the $\hat {u}$ accent, and the Fourier operator (computed with its discrete representation (2.11)) is indicated by $\mathcal{F}[u]$. For each complex Fourier coefficient $\hat {u}_k$, its amplitude

(2.14)\begin{equation} |\hat{u}_k|=\sqrt{\textrm{Re}(\hat{u}_k)^2+\textrm{Im}(\hat{u}_k)^2} \end{equation}

and phase

(2.15)\begin{equation} \varPhi=\tan^{-1}\frac{\textrm{Im}(\hat{u}_k)}{\textrm{Re}(\hat{u}_k)} \end{equation}

can be defined, where $\textrm {Re}(\hat {u}_k)$ and $\textrm {Im}(\hat {u}_k)$ correspond to the real and imaginary parts of $\hat {u}_k$, respectively. All averages are noted by the brackets $\langle \rangle$ and subscripts are listed by the order in which the averaging operations are applied. For instance, $\langle u_z(r,\theta ,z,t)\rangle_{\theta ,t} (r,z)$ is the vertical profile of the vertical velocity field after averaging in the azimuthal direction and in time.

Throughout the paper, we also look at the scaled volume integrated values of the Fourier coefficients defined as

(2.16)\begin{equation} \{\hat{u}_k\}_{V/2{\rm \pi}}(t)=\frac{1}{2{\rm \pi}}\iiint_{\varOmega} \hat{u}_k(r,z,t)\,\textrm{d} V. \end{equation}

Note that the integral defined in (2.16) is also equivalent to

(2.17)\begin{equation} \{\hat{u}_k\}_{V/2{\rm \pi}}(t)=\int_z\int_r \hat{u}_k(r,z,t)r\,\textrm{d} r\,\textrm{d} z. \end{equation}

5.1. Temporal evolution of the flow field

In the previous section it was shown that the large-scale flow transitioned from a structure dominated by a $k=3$ Fourier mode to one dominated by a $k=2$ Fourier mode. In this section the temporal evolution of the transition will be investigated in greater detail using visualization of the unfiltered mid-plane temperature field.

Figure 8 presents snapshots of $\vartheta (r,\theta ,z=0,t)$ at times that bracket the transition each separated by $90t_f$. In figure 8(a) the flow is dominated by the $k=3$ mode, albeit with the warm spokes located at two o'clock being rather weaker than the other two spokes. In figure 8(b) the one o'clock spoke weakens more as the upward plume at the centre of the cylinder becomes stronger. The process continues in figure 8(c) wherein the warm plume at one o'clock has almost vanished, and the central updraft plume has moved toward the spoke at ten o'clock. Finally, in figure 8(d) the central upward plume has merged with the ten o'clock spoke and the one o'clock spoke has vanished completely, leaving a large structure that is clearly dominated by the $k=2$ mode. Compare figure 8(d) to the pure $k=2$ mode in figure 7(c) and the mixture of $k=0,1,2$ modes in figure 7(g). The process illustrated by the snapshots in figure 8 requires $270t_f$, or approximately 8.7 eddy-turnover times.

Figure 8. Instantaneous snapshots of temperature at the mid-plane during the transition of the global pattern. Snapshots are spaced $90t_f$ apart covering approximately 8.7 eddy-turnover time units.

A more detailed investigation of the transition is performed by plotting the scaled volume integrated Fourier coefficients for a given mode defined by (2.16). Volume integration removes the localized spatial variations of the mode and allows the temporal evolution to be investigated from a macro perspective. Even though the volume integrated Fourier coefficients only depend on wavenumber and time, they are still complex variables. The phase and amplitude of the volume integrated coefficient can simultaneously change with time.

Figures 9 and 10 display the temporal evolution of $\hat {\boldsymbol {u}}_r$ and $\hat {\vartheta }$ fluctuations for the first five ($k=0:4$) volume integrated Fourier modes in terms of their phase and amplitude. Note that the $k=0$ modes represent an azimuthal mean of the corresponding physical quantities, and their further integration with (2.17) results in volume averaging in physical space, which is identically zero for $\hat {u}_z$ and $\hat {u}_r$ due to mass conservation. The reader is cautioned that while interpreting figures 9 and 10 to remember that any phase jumps of $2{\rm \pi}$ in the plots are continuous in phase space. In general, the modes in figures 9 and 10 share some common behaviour. The variables and plots are divided into two highly correlated groups, $\hat {u}_z$ with $\hat {\vartheta }$ in figure 9 and $\hat {u}_r$ with $\hat {u}_\theta$ in figure 10. A strong degree of correlation between $\hat {\vartheta }$ and $\hat {u}_z$ dynamics testifies that temperature and vertical velocity are the signatures of the same large-scale structure, consistent with the findings of Krug et al. (Reference Krug, Lohse and Stevens2019). The $\hat {\vartheta }$ and $\hat {u}_z$ coefficients for the low-order modes except for $k=2$ show larger mean amplitudes and smaller fluctuations in phase during the initial part of the simulation. During the transition from a $k=3$ to $k=2$ in the interval $500<t/t_f<1000$, the amplitude of the dominant structures tends to decrease, and the phase fluctuations tend to increase.

Figure 9. Temporal evolution of the scaled volume integrated Fourier coefficients of $\hat {u}_z$ and $\hat {\vartheta }$ for $k=0:4$ plotted in terms of amplitude (a,c,e,g,i) and phase (b,d,f,h,j). Each row of figures corresponds to a separate wavenumber with the top row corresponding to $k=0$ and the bottom row corresponding to $k=4$. The units of $\varPhi$ are in radians. The phase plots have been rescaled to cover a range of $4 {\rm \pi}$ to highlight the low frequency cycles that occur in the temporal evolution of these modes. $(\textit{a},\!\textit{b}) \ k= 0, \ (\textit{c},\!\textit{d}) \ k= 1, \ (\textit{e},\!\textit{f}) \ k= 2, \ (\textit{g},\!\textit{h}) \ k= 3, \ (\textit{i},\!\textit{j}) \ k= 4.$

Figure 10. Temporal evolution of the scaled volume integrated Fourier coefficients of $\hat {u}_r$ and $\hat {u}_\theta$ for $k=0:4$ plotted in terms of amplitude (a,c,e,g,i) and phase (b,d,f,h,j). Each row of figures corresponds to a separate wavenumber with the top row corresponding to $k=0$ and the bottom row corresponding to $k=4$. The units of $\varPhi$ are in radians. The phase plots have been rescaled to cover a range of $4 {\rm \pi}$ to highlight the low frequency cycles that occur in the temporal evolution of these modes. $(\textit{a},\!\textit{b}) \ k= 0, \ (\textit{c},\!\textit{d}) \ k= 1, \ (\textit{e},\!\textit{f}) \ k= 2, \ (\textit{g},\!\textit{h}) \ k= 3, \ (\textit{i},\!\textit{j}) \ k= 4.$

While the signature of the transition from a $k=3$ to $k=2$ dominant structure around $500\text {--}1000t_f$ can be seen in each of the low wavenumber modes, the effects are most clearly displayed in the $k=2$ and $k=3$ modes. Inspection of the amplitude and phase for $k=2$ and $k=3$ shows that the amplitude of $k=3$ starts out at a relatively large value and, beginning at around $500t_f$, decays steadily to a negligible value over another $500t_f$, completing the process around $1000t_f$. Furthermore, while the amplitude is large, the phase remains relatively constant, but once the amplitude dies down the phase fluctuations increase. The opposite effects are seen in the $k=2$ mode where the amplitude is initially low but then gradually increases over the same period that $k=3$ decays.

A simple analysis of the time scales involved with the dynamical processes occurring in Rayleigh–Bénard convection can be presented by estimating a ‘mode turnover time’, which is the time it takes a particle to travel one full circuit in a particular mode $k$. The full circulation length for each mode (normalized with $H$) can be defined as

(5.1)\begin{equation} L=2 + {\varLambda_k},\end{equation}

where $\varLambda _k$ is the $k$th mode wavelength given by (4.1), which yields

(5.2)\begin{equation} L=2\left(1+\frac{{\rm \pi} \varGamma}{2k}\right). \end{equation}

The ‘mode turnover’ time therefore can be defined through a classical eddy-turnover time as

(5.3)\begin{equation} t_{k}=\left(1+\frac{{\rm \pi}\varGamma}{2k}\right) t_{\epsilon}. \end{equation}

Note that in the eddy-turnover time definition, (2.8), the variance of the vertical velocity fluctuations is used as a velocity scale. In the present study the magnitude of the horizontal and vertical velocity fluctuations is of the same order, consistent with the observations in Stevens et al. (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018), so the same velocity scale can be used in the mode turnover time definition. It can be seen that both the mode number and the aspect ratio play a role in the calculation of a mode turnover time. The presented time scale definition is similar in spirit to the filtering time scale defined in Pandey et al. (Reference Pandey, Scheel and Schumacher2018), albeit they use the statistically averaged structure size in their definition, while we adapt the mode circulation length as the length scale, which allows us to distinguish between the temporal scales of the different modes that contribute to the overall dynamics of the large-scale structure. Also note that Pandey et al. (Reference Pandey, Scheel and Schumacher2018) introduce an empirical factor of three into their time scale definition, which ‘accounts for the fact that an individual parcel is not perfectly circulating around in such a roll when the flow is turbulent’. We choose not to introduce any empirical factors, but admit that the proposed time scale is only an order of magnitude estimate, and it might take several of such time scales for a particular event to happen, as discussed below.

Since the mode turnover time reflects the time it takes an information (disturbance) to propagate across the entire mode, it should be representative of the time scales associated with the mode destabilization. The time scale analysis with (5.3) predicts, for an aspect ratio $\varGamma =1$, where a single $k=1$ mode dominates, a time scale of $t_1^{(1)}=2.57 t_{\epsilon }$. For an aspect ratio $\varGamma =1$, Mishra et al. (Reference Mishra, De, Verma and Eswaran2011) identified a time scale of $10 t_{\epsilon }$ as the time scale associated with the LSC reversal, while Brown et al. (Reference Brown, Nikolaenko and Ahlers2005) reported $30t_{\epsilon }$ as the average time between reorientations. However, these time scales refer to the time between the reorientation events, and not to the duration of the events themselves. If one closely examines the data presented in, for example, Brown et al. (Reference Brown, Nikolaenko and Ahlers2005), Mishra et al. (Reference Mishra, De, Verma and Eswaran2011) and Zürner et al. (Reference Zürner, Schindler, Vogt, Eckert and Schumacher2019) one can see that the time scale associated with the process of transition of the global structure into a state with a new LSC reorientation, i.e. from the beginning of the mode destabilization to the time when it stabilizes again, is on the order of $2 t_{\epsilon }$ in Brown et al. (Reference Brown, Nikolaenko and Ahlers2005) and Zürner et al. (Reference Zürner, Schindler, Vogt, Eckert and Schumacher2019), and on the order of $2.63 t_{\epsilon }$ in Mishra et al. (Reference Mishra, De, Verma and Eswaran2011), well in line with the predicted $t_1^{(1)}=2.57 t_{\epsilon }$ from (5.3). Thus, while cessations in unit aspect-ratio cells are often perceived as instantaneous events, they are, in fact, events of a finite, albeit short, duration. While a phase reversal happens almost instantaneously during the LSC reorientation, the other accompanied events, such as a decrease of the $k=1$ mode amplitude preceding the phase reversal (or partial reversal), and a subsequent increase back to its original value, happen more gradually (Brown et al. Reference Brown, Nikolaenko and Ahlers2005; Brown & Ahlers Reference Brown and Ahlers2007). For a current aspect ratio of $\varGamma =6.3$, the time scales associated with the destabilization of the modes are about six times larger for a given $k$. The corresponding time scales for the first three modes $k=1, 2, 3$ would be $t_{1}^{(6.3)}=10.90t_{\epsilon }$, $t_{2}^{(6.3)}=6.00t_{\epsilon }$ and $t_{3}^{(6.3)}=4.30t_{\epsilon }$, which, if scaled with the free-fall time units, are $t_{1}^{(6.3)}\sim 300t_{f}$, $t_{2}^{(6.3)}\sim 200t_{f}$ and $t_{3}^{(6.3)}\sim 130t_{f}$. Note that the duration of the observed 3 to 2 mode transition correlates well with these time scales.

One can observe that as the mode number $k\rightarrow \infty$, the mode turnover time converges to a constant value of $t_{\epsilon }$, irrespective of $\varGamma$. This perhaps can explain a self-similarity of higher-order modes, and an apparent lack of dependence on the problem geometry, as will be illustrated below. This also shows that the time scale of an eddy turnover is still a valid measure to describe the fast processes in the RBC cell.

The 3 to 2 mode transition in the global structure in the current $\varGamma =6.3$ domain occurs on a time scale of approximately $500t_f$, or 16 eddy turnovers, but it involves not only modes $k=2$ and $k=3$, but the other low-order modes, for example, $k=1$. A careful observation of the $k=1$ mode behaviour in figure 9(c) shows a rather rapid decrease in its amplitude around the time $t=500t_f$ followed by a relatively rapid increase back to its original value, with another event of a vanishing amplitude for $k=1$ at $t=1000t_f$. A rapid decay in amplitude remarkably resembles the cessations observed in Mishra et al. (Reference Mishra, De, Verma and Eswaran2011) where it was shown that the $k=1$ mode amplitude rapidly vanishes during cessation. A close inspection of the phase diagram in figure 9(d) reveals that this rapid decrease is indeed accompanied by a phase shift close to $180^{\circ }$, pointing to a reorientation in a mode 1. These mode 1 cessations happen to fall onto the interval corresponding to the mode $k=3$ to 2 transition, which suggests that the two events might be associated with each other. Interestingly, the duration between the consecutive cessations correlates well with $t_1^{(6.3)}$ predicted by (5.3), which might signify that they are associated with the same $k=1$ destabilization event and might resemble double cessations identified in $\varGamma =1$ cells (Xi et al. Reference Xi, Zhou and Xia2006; Mishra et al. Reference Mishra, De, Verma and Eswaran2011), albeit on longer time scales here than in unit aspect-ratio domains, commensurate with the 6.3 times difference in the aspect ratio. Note that mode 2 also exhibits an event similar to a cessation at a time of approximately 300$t_f$ where its amplitude essentially vanishes and its orientation rapidly changes. While such cessations in a dominant mode would lead to a complete reorientation of the structure, they do not result in observable changes when the modes are not dominant. This once again testifies of a relative complexity of a wide aspect-ratio system, where a global reorganization is a more complex process that involves the mode interaction.

In general, the $k=2$ and $k=3$ modes exhibit similar characteristics during the time periods when each respective mode is the most dominant mode in the flow field. However, there is one notable difference between the temporal evolution of the phases for these two modes. While the phase of $k=3$ remains virtually constant during its period of dominance (albeit showing signs of low-amplitude, low-frequency fluctuations), the $k=2$ mode's phase changes at a relatively constant rate. This indicates a steady rotation event for the $k=2$ mode that becomes strikingly clear when the scaled volume integrated coefficients are plotted on the complex plane. Plotting on the complex plane (see figures 11 and 12) is another intuitive way to view the changes in amplitude and phase for a given wavenumber.

Figure 11. Temporal evolution of the scaled volume integrated Fourier coefficients plotted on the complex plane for $k=2$ (a,c) and $k=3$ (b,d), where (a,b), $\circ =\hat {u}_r$, and (c,d), $\unicode{x2B20} ={\hat {u}_\theta }$. Markers are colour-coded according to their time, from blue (start of the simulations) to red (finish).

Figure 12. Temporal evolution of the scaled volume integrated Fourier coefficients plotted on the complex plane for $k=2$ (a,c) and $k=3$ (b,d), where (a,b), $\Box =\hat {u}_z$, and (c,d), $\triangle =\hat {\vartheta }$. Markers are colour-coded according to their time, from blue (start of the simulations) to red (finish).

The trajectory of the $k=2$ mode for the temperature $\hat {\vartheta }$ in figure 12(c) begins near the origin and as time progresses it tracks up along the imaginary axis and then begins to drift into quadrant two of the real-complex plane. The rate of rotation manifested by a steady increase in the phase angle $\varPhi$ is measured to be 1.1 degrees per eddy turnover by employing a least squares fitting of the $k=2$ phase in figure 9(f) from a time of $1800t_f$ to the end of the simulation. Careful inspection of figure 5 also shows that a very slow clockwise rotation is starting to occur in the large-scale structure. However, it is hard to discern by just looking at visualizations of the flow field because the individual lobes of the large-scale structure modulate and shift in size. Figures 12(a) and 12(c) give a more clear indication that rotation is indeed occurring in the large-scale structure of the flow.

The trajectory of $\hat {\vartheta }$ for $k=3$ in figure 12(d) begins in quadrant three of the real-complex plane. Low-amplitude fluctuating azimuthal motions of this mode manifested in a swinging shift in phase are noticeable in the early parts of the trajectory. As time progresses, the $\hat {\vartheta }$ coefficient moves toward the origin and upon arrival it begins to fluctuate about the origin.

As a general observation, when the mode is the dominant mode in the large-scale structure, its complex Fourier coefficient drifts away from the origin, and when it loses its dominance it becomes centred around the origin with a rapidly changing amplitude scattered between zero and some threshold value. The same can be said about the behaviour of $\hat {u}_r$ and $\hat {u}_{\theta }$ coefficients, even in the modes 2 and 3, shown in figure 11, which reflects the fact that these variables do not play an active role in the global structure dynamics. The key observation here is that active components of the dominant modes display persistence in terms of phase and amplitude of their Fourier coefficients, while the supporting modes display chaotic behaviour with zero mean in all their Fourier coefficients.

Even though there appears to be no net rotation when the $k=3$ mode is dominant on a time scale while it persisted, a rotation event is occurring in the $k=1$ mode. Figure 9(d) shows a steady change in phase for the first $1500t_f$ of the simulation, after which it suddenly begins to see large fluctuations in phase like the $k=3$ mode.

During the first $500t_f$, the $k=1$ mode rotates approximately 60$^{\circ }$ with a least squares fit providing a rotation speed of approximately 3 degrees per eddy turnover. The phase jumps up and down during the period 500–1000 $t_f$ between the two cessations (marking the period of mode 3 to 2 transition), when at the time of $1000 t_f$, a sudden acceleration occurs, and a higher value of rotation velocity persists until $1500 t_f$, at which time the phase randomizes. Interestingly, the moment of sudden acceleration in $k=1$ rotation coincides with a second cessation of this mode, and a completion of $k=3$ to $k=2$ transition. Brown & Ahlers (Reference Brown and Ahlers2007) reported a similar phenomenon of a rapid acceleration of motion in the azimuthal direction accompanying a cessation predicted by their stochastic LSC model which was explained by the fact that when the amplitude becomes small, the azimuthal motion becomes fast, due to a reduction in the LSC angular momentum. The above observations indicate that there may be a connection between the dynamics of the $k=1$ mode, i.e. rotation and cessation, and the transition of the flow's global structure. The $k=1$ mode rotation and cessation supports our observation that the transition is marked by the movement of the central column, while the phase fluctuations of the $k=3$ mode can explain the azimuthal shifts of the large-scale structure lobes along the edge of the domain (see figure 8).

The rotation rate of 3 degrees per eddy turnover corresponds to a time scale of rotation of approximately 40–60 eddy turnovers (per 1/2 revolution). A rotation rate of 1.1 degree per eddy turnover observed for mode 2 corresponds to the time scale of 160 eddy turnovers (per 1/2 revolution). Interestingly, very different time scales associated with the rotation were also observed in Brown et al. (Reference Brown, Nikolaenko and Ahlers2005). Fast rotation typically classified as a reorientation seemed to scale with the time of about 10 eddy turnovers for $\varGamma =1$, which is approximately $3.6t_1^{(1)}$, while much slower drift that was not associated with a reorientation event would occur on time scales of an order of magnitude larger. Similar time scale discrepancy in a duration between different reorientation events was observed with $\varGamma =1$ in Sreenivasan, Bershadskii & Niemela (Reference Sreenivasan, Bershadskii and Niemela2002). The current $\varGamma =6.3$ cell is shown to exhibit similar dynamics, where the rotation of mode 1 accompanied by a mode transition scales on a mode turnover time $3.6\text {--}5.5t_1^{(6.3)}$, while a slow rotation in a persistent mode 2, not undergoing a transition, is approximately $26.6t_2^{(6.3)}$. Brown & Ahlers (Reference Brown and Ahlers2006, Reference Brown and Ahlers2007) attributed slow rotations, or drifts, to diffusive processes. Viscous time scale is given by $t_v=\sqrt {Ra/Pr}\,t_f$ (Pandey et al. Reference Pandey, Scheel and Schumacher2018), which, for given flow parameters, is equal to $3800t_f$, or $120t_{\epsilon }$, while a thermal diffusion time scale is $t_d=\sqrt {Ra\,Pr}\,t_f=25\,460 t_f=804t_{\epsilon }$ for the current $Pr=6.7$ case. It is seen that the slow mode rotation events (‘azimuthal jitter’) are expected to occur on very long time scales in the current case, and a slow rotation observed in mode 2 indeed falls in between these time scales.

So far virtually all of the discussion for the volume integrated Fourier coefficients has centred on the $\hat {\vartheta }$ and $\hat {u}_z$ fields. These two fields are highly correlated, and they tend to describe the events that are oriented in the inhomogeneous, vertical direction, and provide a profound amount of information regarding the time scales and dynamics that are associated with the large-scale structure's transition.

Figure 10 shows that the amplitude of $\hat {u}_r$ changes rapidly in time, but the magnitude of these fluctuations stays within a specific band for each wavenumber. Meanwhile, the phase of $\hat {u}_r$ for each wavenumber tends to evolve at a slower pace with several low frequency components. These time scales are measured by performing a fast Fourier transform (FFT) analysis of the $\hat {u}_r$ and $\hat {u}_\theta$ signals in figure 10 (using time as the abscissa). This analysis allows us to identify the temporal frequencies and their accompanying periods that are associated with the consistent phase oscillations in these fields. Since these temporal signals are not strictly periodic, a Blackman window function is multiplied with the signal prior to taking the FFT to improve peak detection by limiting spectral leakage (Blackman & Tukey Reference Blackman and Tukey1958). The measured peaks are not consistent between wavenumbers, but the general trend is that the $\hat {u}_r$ and $\hat {u}_\theta$ amplitudes fluctuate with time scales on the order of 1 eddy-turnover time, while the phases change at a much slower rate with time scales of $O(10)$ eddy turnovers. This is consistent with the behaviour that is exhibited by $\hat {u}_z$ and $\hat {\vartheta }$ in the non-dominant modes. Non-dominant modes are characterized by the mode numbers $k$ below 2, and $k$ greater than 3. The $k=2$ mode is non-dominant at a time below $500t_f$, and the $k=3$ mode is non-dominant at a time above $1000t_f$. The primary difference between $\{\hat {u}_z, \hat {\vartheta }\}$ and $\{\hat {u}_r, \hat {u}_\theta \}$ groups of variables is that the $\{\hat {u}_r, \hat {u}_\theta \}$ coefficients of the non-dominant as well as the dominant modes show no noticeable change in the behaviour of their phase and amplitude throughout the course of the simulation, and, thus, seem to be unaffected by the large-scale mode transition, as was previously observed in figure 11.

The distinction between the coefficients of the horizontal velocity components and $\hat {u}_z$ and $\hat {\vartheta }$ becomes less clear at higher wavenumbers. The temporal evolution of two additional modes ($k=5$ and $k=10$) is plotted in figure 13 to show how the phase and amplitude of the volume integrated coefficients are affected by increasing Fourier wavenumber. Figure 13 shows that as the wavenumber increases, so do the frequencies associated with the temporal evolution of the coefficients, consistent with the $t_k$ prediction in (5.3). The rapid oscillations and lack of distinction between the behaviour of the individual variables indicates that the high-wavenumber fields are becoming increasingly random in nature, and less descriptive of the physical system as a whole. Also, it is worth noting that the amplitudes of the high-frequency modes are higher for the vertical and radial velocities as opposed to the temperature and the azimuthal velocity. Stronger high-frequency contributions to the energy content for the vertical velocity rather than the temperature were also previously observed by Krug et al. (Reference Krug, Lohse and Stevens2019).

Figure 13. Temporal evolution of the scaled volume integrated Fourier coefficients of $\hat {u}_z$, $\hat {\vartheta }$ (left pane) and $\hat {u}_r, \hat {u}_{\theta }$ (right pane) for $k=5$ (a–d), $k=10$ (e–g) plotted by amplitude (left) and phase (right). The phase plots have been rescaled to cover a period of $4 {\rm \pi}$ to highlight the low frequency cycles that occur in the temporal evolution of these modes.

5.2. Integral time scale

In this section we attempt to quantify the temporal behaviour of the modes by measuring the classical integral time scale, $\mathcal {T}$. Normally $\mathcal {T}$ is interpreted to be the coherence time of the flow field. It is defined in terms of the autocorrelation function,

(5.4)\begin{equation} R_{ii}(\varOmega,\tau)=\langle \psi_i(\varOmega,t+\tau)\psi_i(\varOmega,t)\rangle_{t},\end{equation}

where $\psi$ is a vector containing variables of interest and $\tau$ is the temporal offset between the two instances of the flow field, or snapshots. Usually $\psi$ is taken to be the turbulent velocity field $(\psi =\{u_r',u_\theta ',u_z'\})$, where the prime indicates the fluctuating portion of the velocity. An autocorrelation based on this particular vector will determine a correlation based on the turbulent kinetic energy. A summation over the indices $i=1\ldots d$, where $d$ is the dimension of the vector, is then implied in the definition of the autocorrelation function in (5.4).

The interest of this work includes the global correlation times as well as the time scales of the individual Fourier modes, since the large-scale structures have been shown to contain a relatively small number of Fourier modes. In terms of the Fourier coefficients

(5.5)\begin{equation} R_{ii}(r,\theta,z,\tau)=\left\langle\sum_{k=-\infty}^\infty \sum_{k'=-\infty}^\infty \hat{\psi}_i(r,k,z,t+\tau)\hat{\psi}_i(r,k',z,t)\,\textrm{e}^{\,j(k+k')\theta}\right\rangle_t \end{equation}

averaging (5.5) over $\theta$ gives us

(5.6)\begin{align} \langle R_{ii}(r,\theta,z,\tau)\rangle_\theta&=R_{ii}(r,z,\tau) \nonumber\\ &=\sum_{k=-\infty}^\infty \sum_{k'=-\infty}^\infty \langle\langle\hat{\psi}_i(r,k,z,t+\tau)\hat{\psi}_i(r,k',z,t)\,\textrm{e}^{j(k+k')\theta}\rangle_t\rangle_\theta \nonumber\\ &=\sum_{k=-\infty}^\infty \sum_{k'=-\infty}^\infty \langle\hat{\psi}_i(r,k,z,t+\tau)\hat{\psi}_i(r,k',z,t)\langle \,\textrm{e}^{j(k+k')\theta}\rangle_\theta\rangle_t \nonumber\\ &=\sum_{k=-\infty}^\infty \sum_{k'=-\infty}^\infty \langle\hat{\psi}_i(r,k,z,t+\tau)\hat{\psi}_i(r,k',z,t)\delta_{k,-k'}\rangle t, \end{align}

where $\delta _{k,-k'}$ is the Kronecker delta function, in which statistical stationarity in $\theta$ requires $k=-k'$.

Since all the flow variables are real signals, the negative Fourier mode $-k$ can be expressed as the complex conjugate of the positive Fourier mode $k$. Therefore, the Dirac delta function in (5.6) shows that all wavenumbers will contribute to the correlation when multiplied by their complex conjugates. This also ensures that the correlation will be comprised entirely of real numbers which is required since the dependent variables are all defined in real space. The discrete representation of (5.6) is

(5.7)\begin{equation} R_{ii}(r,z,\tau) =\sum_{k=-N_\theta/2-1}^{N_\theta/2-1}\langle\hat{\psi}_i(r,k,z,t+\tau) \hat{\psi}^*_i(r,k,z,t)\rangle_{t} , \end{equation}

where $N_\theta$ is the number of samples for the Fourier transform in the $\theta$ direction and $^*$ indicates the complex conjugate. Defining an autocorrelation function per wavenumber

(5.8)\begin{equation} R_{ii}(r,k,z,\tau)=\langle\hat{\psi}_i(r,k,z,t+\tau)\hat{\psi}^*_i(r,k,z,t)\rangle_{t} \end{equation}

and recognizing that

(5.9)\begin{equation} R_{ii}(r,-k,z,\tau)=R_{ii}^*(r,k,z,\tau), \end{equation}

which follows from the corresponding properties of the Fourier coefficients (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988), one can see that (5.7) can be rewritten as

(5.10)\begin{align} R_{ii}(r,z,\tau) &= R_{ii}(r,0,z,\tau)+\sum_{k=1}^{N_\theta/2-1} R_{ii}(r,k,z,\tau)+R_{ii}(r,-k,z,\tau) \nonumber\\ &=R_{ii}(r,0,z,\tau)+\sum_{k=1}^{N_\theta/2-1} R_{ii}(r,k,z,\tau)+R_{ii}^*(r,k,z,\tau) \nonumber\\ &=R_{ii}(r,0,z,\tau)+\sum_{k=1}^{N_\theta/2-1} 2\textrm{Re}\{R_{ii}(r,k,z,\tau)\}, \end{align}

where $\textrm {Re}\{R_{ii}(r,k,z,\tau )\}$ is the real part of the autocorrelation (5.8).

Utilizing (5.7) and (5.8), the total integral time scale $(\mathcal {T})$ and the integral time scales for each wavenumber $(\mathcal {T}_k)$ are

(5.11)\begin{equation} \mathcal{T}(r,z) = \int_0^{\infty} \frac{R_{ii}(r,z,\tau)}{R_{ii}(r,z,0)}\textrm{d}\tau, \end{equation}

(5.12)\begin{equation}\mathcal{T}_k(r,z) = \int_0^{\infty} \frac{R_{ii}(r,k,z,\tau)+R_{ii}(r,-k,z,\tau)}{R_{ii}(r,k,z,0)+R_{ii}(r,-k,z,0)}\textrm{d}\tau,\quad k\ge 0, \end{equation}

which is equivalent to

(5.13)\begin{equation} \mathcal{T}_k(r,z)= \int_0^{\infty} \frac{\textrm{Re}\{R_{ii}(r,k,z,\tau)\}}{\textrm{Re}\{R_{ii}(r,k,z,0)\}}\textrm{d}\tau,\quad k\ge 0.\end{equation}

The quantity under the integral in (5.11) can also be referred to as a normalized autocorrelation,

(5.14)\begin{equation} \bar{R}_{ii}(r,z,\tau)=\frac{R_{ii}(r,z,\tau)}{R_{ii}(r,z,0)}.\end{equation}

These integral time scales depend on the field $\psi$ that is chosen. $\psi =\{u_r', u_{\theta }',u_z'\}$, We used $\psi =\{\vartheta '\}$ and $\psi =\{u_r', u_{\theta }',u_z',\vartheta '\}$ to perform proper orthogonal decomposition in the study of Bailon-Cuba et al. (Reference Bailon-Cuba, Emran and Schumacher2010) on data from turbulent RBC simulations at various $\varGamma$. We shall also, likewise, be referring to the norms of these fields as the turbulent kinetic energy, thermal energy and total energy, respectively. Thus, in the definition of the autocorrelation functions in (5.4)–(5.14), the values $i=1:3$ are assumed for the turbulent kinetic energy, $i=4$ for the thermal energy and $i=1:4$ for the total energy.

Figure 14(a–c) shows $\mathcal {T}$ for the entire field when $\psi$ is defined as the turbulent kinetic energy, the turbulent thermal energy and the total turbulent energy, respectively. The $r$–$z$ plots of $\mathcal {T}$ also give insight into the structure of the flow field by indicating which regions of the flow field have longer correlation times, and also by how much the correlation times vary.

Figure 14. Spatially varying, azimuthally averaged integral time scales $\mathcal {T}(r,z)$ based on the kinetic energy (a) (min: $12.9 t_f$, max: $617 t_f$); temperature fluctuations (b) (min: $10.5 t_f$, max: $663 t_f$); and total turbulent energy (c) (min: $13.8 t_f$, max: $569 t_f$). White boxes $(\square )$ and circles $(\circ )$ indicate the locations of minimum and maximum integral scales, respectively. Here $\tau$ is in multiples of $t_f$. All length scales are normalized with $H$.

Figures 14(a) and 14(b) show very different behaviour between the correlation of the turbulent kinetic and thermal energy fields. The two fields have little overlap between regions with very long correlation times. The kinetic energy field has its longest correlation times in the boundary layer while the turbulent thermal energy field has its longest correlation times in the bulk region. Regions where the fluctuations change sign frequently tend to have lower correlation times while regions where fluctuations maintain the same sign for long periods of time have longer correlations.

Further understanding of these correlation times can be gained by reviewing the variance of the various components that comprise the correlation metrics. Variance of the velocity and temperature fluctuations is plotted in figure 15 to show where the strongest fluctuations most frequently appear. Variance is defined as azimuthally and temporally averaged turbulence fluctuations $\sigma _{u_i}(r,z)=\langle u_i'(r,\theta ,z,t)^2\rangle _{\theta , t}$. The peak variance of the velocity components are aligned with the shape of the large-scale structures. The horizontal components ($u_r$ and $u_\theta$) have a strong variance in the boundaries where the large-scale roll cells will have the strongest horizontal velocity contributions. Likewise the vertical velocity ($u_z$) has a larger variance in the bulk region where the large-scale structures are principally updrafts or downdrafts. This is where the opposing plumes pass one another as they cross the layer depth (see figure 15d), and this appears to be the factor that is driving the correlation times of the turbulent kinetic energy field down in the bulk region, due to a prevalence of the opposite sign fluctuations in this region.

Figure 15. Variance $\sigma = \langle u'^2(r,\theta ,z,t)\rangle _{\theta ,t}$ of $\vartheta$ (a) (min: $5.48\times 10^{-3}$, max: $2.29\times 10^{-1}$), $u_r$ (b) (min: $0.0$, max: $6.80\times 10^{-3}$), $u_\theta$ (c) (min: $0.0$, max: $9.82\times 10^{-3}$), $u_z$ (d) (min: $0.0$, max: $8.12\times 10^{-3}$), turbulent kinetic energy (e) (min: $0.0$, max: $1.68\times 10^{-2}$), and total energy (f) (min: $2.33\times 10^{-3}$, max: $2.29\times 10^{-1}$) in the $r$–$z$ plane. White boxes $(\square )$ and circles $(\circ )$ indicate the locations of minimum and maximum values, respectively. Also, note that the variance for the velocity components and turbulent kinetic energy is analytically $0.0$ at the walls. All length scales are normalized with $H$.

The variance of $u_z$ also shows a peak in the central region near $r=0$ which can be attributed to the time period when $k=3$ was the dominant structure, and the uncertainty as the lobes of the $k=2$ structure dance back and forth across the centre of the convection cell.

The peak variance of the temperature near the boundaries in figures 15(a) and 16(b) highlight how well mixed the thermals are in the bulk region. Note that the variance of the total energy (figure 15f) has its peak values in the thermal boundary layers, the same as the variance of temperature, showing a noticeable contribution of the thermal fluctuations into a variance of the total energy field. However, in the bulk region, the variance of the total turbulent energy is similar to the variance in kinetic energy. This also explains while integral time scales shown in figure 14 are similar for the total and kinetic turbulent energy fields in the bulk region, but different from the thermal energy field: the large integral time scales for the temperature in the bulk are generated by relatively small temperature fluctuations divided by a small variance (the denominator of (5.11)), which account for a negligible contribution to the time scales when added to the velocity field with larger values of the fluctuations and the larger variance in the bulk. However, an increase of integral times scales for the total turbulent energy on the side walls due to a thermal field contribution is pronounced.

Figure 16. (a) Skewness of temperature, (5.15), in the $r$–$z$ plane; line plots of variance (b) and skewness (c) of temperature at select radial locations. All length scales are normalized with $H$.

Figure 15 also shows the skewness of the temperature field defined as

(5.15)\begin{equation} \mathcal{S}_{\vartheta}(r,z)=\frac{\langle\vartheta'(r,\theta,z,t)^3\rangle_{\theta,t}}{\sigma_{\vartheta}(r,z)^{3/2}}, \end{equation}

which is a measure of asymmetry of the probability distribution of a fluctuating temperature field about its mean. It can be seen that the skewness of temperature peaks just as the variance begins to increase, and the two quantities are anti-correlated as they approach the wall (see figure 16b,c). The skewness profile shows that even though the fluctuations about the mean are very small in the bulk region, the lower half of the domain is still biased toward warmer fluid and the upper half toward colder.

When the integral time scales are averaged over the volume, the differences between the three metrics narrows (see table 1). Following (2.16), volume integrated integral time scales (scaled with $2 {\rm \pi}$) are defined as

(5.16)\begin{equation} \{\mathcal{T}_k\}_{V/2{\rm \pi}} =\frac{1}{2{\rm \pi}}\iiint_{\varOmega} \mathcal{T}_k(r,z)\,\textrm{d} V=\int_z\int_r\mathcal{T}_k(r,z)r\,\textrm{d} r\,\textrm{d} z. \end{equation}

The turbulent kinetic energy shows the correlation time averaged across all Fourier modes that is about 50 % longer than the turbulent thermal energy, but the total turbulent energy shows about the same correlation time as the turbulent thermal energy. This means that there are regions where correlation among the velocity components cancels out correlation from the thermal region. Table 1 and figure 17 further demonstrate the effect of the dominant Fourier modes where the scaled volume integrated integral time scales are shown per mode number, $\{\mathcal {T}_k\}_{V/2{\rm \pi} }$. Here it can be seen that correlation times for the low-order modes, and specifically the $k=2$ mode increase in magnitude in all the presented metrics. It is also noteworthy to see that the correlation time of the $k=2$ mode is dramatically larger than the other modes with the second longest mode ($k=0$ for thermal energy and $k=3$ for kinetic and total energy) being almost five times shorter.

Figure 17. Scaled volume integrated integral time scales (in multiples of $t_f$) for each mode number, $\{\mathcal {T}_k\}_{V/2{\rm \pi} }$. Modes are plotted versus $k+1$ to make the $k=0$ mode visible on the log-scale plot and $u_i$ is a shorthand reference for all three velocity components. This data includes the values listed in table 1, and all the additional wavenumbers that were not included in table 1.

Table 1. Scaled volume integrated integral time scales (in multiples of $t_f$) for the total field $\{\mathcal {T}\}_{V/2{\rm \pi} }$ and a selection of Fourier modes $\{\mathcal {T}_k\}_{V/2{\rm \pi} }$.

Table 1 and figure 17 clearly show that the low frequency modes are responsible for the large integral time scales, and that the majority of the correlation time can be attributed to the $k=2$ mode. Further evidence of this can be seen by comparing the spatial distribution of the integral time scale for $k=2$ (figure 19c) and the entire system (figure 14c). Figure 17 also shows that $\mathcal {T}_k$ decays to a value of approximately $3t_f$ for all three energy vectors after the first 12 Fourier modes. This is the minimum limit that can be obtained with this dataset since the snapshots were sampled $3t_f$ apart. Shorter $\mathcal {T}_k$'s are probable for the higher wavenumbers.

6.1. Spatial variability of the integral time scales

In the previous section we presented the integral time scales for three different quantities. Moving forward we will restrict our analysis to the integral time scale of the total turbulent energy since it is an aggregate of the other two. The normalized autocorrelation for the total turbulent energy field, $\bar {R}_{ii}(r,z,\tau ),\ i=1:4$, is plotted versus snapshot spacing in figure 18 at a selection of points in the $r$–$z$ plane. While the entire simulation time is over 3000$t_f$, only about half of that period is used to calculate the autocorrelations, due to the large values of the maximum time separations.

Figure 18. Temporal normalized autocorrelation based on total turbulent energy at select points throughout the domain. Subplot (a) shows the normalized autocorrelation $\bar {R}_{ii}(r,z,\tau ),\ i=1:4$, (5.14), and subplot (b) marks $(r,z)$ points where the temporal normalized autocorrelations are calculated. $\tau$, $\mathcal {T}$ are in multiples of $t_f$. All length scales are normalized with $H$.

By definition (see (5.11)), $\mathcal {T}(r,z)$ is equal to the area under each of the plots in figure 18(a). Figure 18(a) shows that while fluctuations at points $C$ and $D$ in the highly correlated viscous boundary layer monotonically decay, they remain correlated over the entire dataset. Linear extrapolation of the lines in figure 18 can be used to estimate the time it will take for $\bar {R}_{ii}(r,z,\tau )$ at points $C$ and $D$ to reach a value of zero. This time turns out to be approximately $90$–$100$ eddy turnovers $({\approx }2700\text {--}3000 t_f)$. Based on the measured rotation rate of $1.1$ degrees per eddy turnover, a rotation during the decorrelation time is ${\approx }90^{\circ }$, corresponding to the $k=2$ mode changing from positive to negative vertical velocity. The updraft and downdraft of the $k=2$ mode will exactly cancel one another. This suggests that decorrelation of the most persistent structures in this flow is caused by the observed, global rotation. The other two probes are taken at the mid-plane. The probe at point $B$ is at a local minimum in $\mathcal {T}$ and shows sufficient decay in $\bar {R}_{ii}(r,z,\tau )$ to indicate that the values become uncorrelated during this computation. The other probe at point $A$ is near a local maxima in $\mathcal {T}$. It shows signs of a weak long-lived transient as the correlation decays to zero with a separation time of approximately $800 t_f$, but then begins to grow again. Some possible sources for this transient include the low-amplitude, low-frequency transient in the magnitude of $\hat {\vartheta }'$ for $k=0$ (see figure 9a), which can potentially be associated with a time scale of updraft and downdraft reversals (Sakievich et al. Reference Sakievich, Peet and Adrian2016).

Note also the higher-frequency oscillations in the correlation function at all the probes $A, B, C$ and $D$. The time scale of these higher-frequency processes stays very close between all the probes and varies between 84–94 $t_f$, which corresponds to roughly 3 eddy turnovers. This time correlates well with the mode turnover time for the high-order modes that converges to a scale of an eddy turnover. Since the correlation functions here include the representation of all the modes, it is probable that higher-order modes are responsible for these high-frequency oscillations. Similar oscillations in a correlation function have been observed in Xi et al. (Reference Xi, Zhou and Xia2006) and Mishra et al. (Reference Mishra, De, Verma and Eswaran2011) in unit aspect-ratio cells with a time scale close to an eddy-turnover time. This suggests a similar origin of these oscillations in different aspect-ratio cells, coming from the contribution of high-order modes which are independent of geometry.

Additional insight into the spatial variance of $\mathcal {T}_k(r,z)$ can be found by investigating the contribution from the individual Fourier modes. Recall that $R_{ii}(r,z,\tau )$ can be defined as a summation of $R_{ii}(r,k,z,\tau )$ over all $k$'s. A $\mathcal {T}_k(r,z)$ field (5.13) can be calculated for each Fourier mode giving an indication as to how the individual modes contribute in the total correlation. Plots of the total turbulent energy based $\mathcal {T}_k(r,z)$ for a selection of Fourier modes is provided in figure 19. Note that the values of the integral time scales for some modes can be negative, which explains why the modal values at some locations exceed their cumulative value in figure 14(c). Only low wavenumber plots are included in figure 19 due to the short correlation times of modes $k>10$ (see figure 17 and table 1).

Figure 19. Spatially varying integral time scale based on total turbulent energy for modes $k=$ 0 (a) (min: $-0.939 t_f$, max: $397 t_f$), $k=1$ (b) (min: $-2.05 t_f$, max: $77.0 t_f$), $k=2$ (c) (min: $0.083 t_f$, max: $1260 t_f$) , $k=3$ (d) (min: $0.015 t_f$, max: $392 t_f$), $k=4$ (e) (min: $-1.02 t_f$, max: $1050 t_f$) and 10 (f) (min: $-7.07 t_f$, max: $45.3 t_f$), where white boxes $(\square )$ and circles $(\circ )$ indicate the locations of minimum and maximum integral scales, respectively. Time scale $\mathcal {T}$ is in multiples of $t_f$. Note the difference in maximum $\mathcal {T}$ in each subfigure. All length scales are normalized with $H$.

One observation of the subplots in figure 19 is that the globally dominant, highly correlated modes (subplots (c) and (d), and also subplot (a)) show a high level of symmetry about the mid-plane with long correlation times covering a substantial portion of the $r$–$z$ plane, but the other modes do not. Subplots (b), (e) and (f) also show much smaller peak values for $\mathcal {T}_k$. The long correlation times in the dominant modes support our previous observations that the $k=2$ and $k=3$ modes are the primary sources for the long correlation times, and the high level of symmetry indicates that the shapes of these modes do not see much spatial variation over time. It is interesting to see that the peak in the correlation time $\mathcal {T}_k$ for the mode $k=0$ (subplot (a)) is the same as for the mode $k=3$, but the spatial locations of the highly correlated events are different. This supports an earlier observation that the modes $k=0$ and $k=3$ coexisted when the mode $k=3$ was strong, where mode $k=0$ was dominant in the central region and the mode $k=3$ was responsible for the formation of the side plumes. The lack of spatial symmetry and smaller range of $\mathcal {T}_k$ in the non-dominant modes indicate that these modes see a large amount of variation in their spatial structure and could contain rare energetic events in the flow field which have life spans much less than the length of the simulation, but much longer than the sampling rate of $3t_f$.

6.2. Statistics of the Fourier modes and their spatial variability

In this section the temporally averaged statistics is presented for the Fourier modes which allows one to judge about the spatial variability of modes and the contribution of different length scales into the overall flow structure. This is achieved by evaluating the time-averaged energy spectra of the modes at different $r$–$z$ locations. Up to this point in the paper all data has been presented with respect to the azimuthal Fourier modes, which, as discussed in § 2.3, are simply the integer mode indicators and are not directly related to the structure sizes. Therefore, seven different locations is provided iexamining Fourier coefficients at different radii corresponds to different physical length scales and energy densities per unit length. A more consistent way to compare the flow structure at various locations in the flow field is to normalize the energy spectra and frequency with respect to a geometric length scale, which is the wavelength $\lambda _k(r)=2 {\rm \pi}r/ k$ defined in (2.13). This is done by premultiplying the energy spectra with the radial location and plotting against the inverse of the wavelength $1/\lambda _k(r)$. A sampling of the spectra at seven different locations is provided in figure 20. These locations are at various points within the boundary layers (bottom plate and side walls), and bulk region of the flow field; regions where different physical phenomena dominate. Here $z=-0.45$ and $r=3.1$ are within the viscous boundary layers for the bottom and side walls, respectively, while $z=-0.4$ is just outside the viscous boundary layer in the vertical direction.

Figure 20. Time-averaged energy spectra for each of the components in the total turbulent energy vector at various locations in the flow field. Subplots (a) and (b) are for the temperature field, (c) and (d) are for the radial velocity component, (e) and (f) are for the azimuthal velocity component and (g) and (h) are for the vertical velocity component. Subplots (a), (c), (e) and (g) are at a fixed height of $z=-0.4$, and various radii. Subplots (b), (d), (f) and (h) are at a fixed radius $r=2.0$ and various vertical locations. All length scales are normalized with $H$.

6.2.1. Variations in radial location

Spectra of $\vartheta$, $u_\theta$ and $u_z$ evaluated at various radial locations with a fixed height $z=-0.4$ (figure 20a, c, e and g) show a good collapse across virtually all length scales. For length scales that are greater than $\varGamma =6.3$ ($k/2{\rm \pi} r\lessapprox 2\cdot 10^{-1}$), $\vartheta$ and $u_\theta$ collapse less well, and this behaviour is also seen in the $u_r$ plot. From the plots, it can be concluded that the effect of the side wall boundary conditions is rather small and is mostly prominent in the $u_r$ variable. It is not unexpected, since $u_r$ analytically must tend to zero at the wall. In fact, the same can be said about the centreline location, where $u_r$ is, again, analytically zero, and a lack of collapse in $u_r$ is again observed at low wavenumbers. The same behaviour is manifested in $u_{\theta }$. Other than $u_r$, the side walls influence the large length scales in all the variables but $u_z$, which shows an almost perfect collapse across all radii.

Failure to collapse in the larger length scales can be attributed to the dominance of low-order Fourier modes that describe the flow field's large-scale structure. Since Fourier mode $k=2$ contains a large amount of energy throughout the entire domain it will disrupt the collapse of the spectra by affecting different length scales at each radii. In fact, if the spectra were to collapse across all length scales for all variables then it would be horizontally homogenous as in the canonical form of RBC with infinite $\varGamma$. In a sense the side walls of the convection cell act as a high pass filter because they limit the size of the largest length scales that can be observed in the flow. The fact that the $k=2$ mode dominates the energy spectra at multiple length scales indicates that the underlying structure has a modal nature, and that it is the principle cause for radial inhomogeneity. This is most likely due to the confining, geometric effects of the cylinder. Table 2 illustrates this point by providing the time-averaged energy values for the first several Fourier coefficients at several different radial locations at $z=-0.4$. For a reference, the scaled volume integrated value of $|\vartheta _k|^2$ corresponding to the values in figure 3 is also provided. This data complements the data in figure 20 for $z=-0.4$. Table 2 shows that the $k=2$ mode is indeed energetically dominant at two of the three tabulated radii, and is responsible for the peaks in energy observed in figure 20. The $k=0$ mode is approximately 40 % more energetic at $r=1.0$, $z=-0.4$. This is most likely due to the residual effects of the central updraft early in the time series combined with geometric effects. Careful observation of figure 20(a) shows that as the radius decreases, the strength of the energy peak associated with the $k=2$ mode decays. Due to a singular nature of approaching $r=0$, an azimuthal alignment of the most dominant energetic structures might break, contributing its energy to an azimuthally invariant $k=0$ mode. However, it is worth noting that the $k=2$ mode is still larger than the $k=1$ and $k=3$ modes at $r=1.0$, $z=-4.0$, thus breaking a monotonic decay of the energy spectrum and displaying a level of dominance as the second most energetic mode.

Table 2. Time-averaged energy of Fourier coefficients from the temperature field with variations in the radial location for the first 16 modes at $z = -0.4$.