Two-dimensional Rayleigh–Bénard convection without boundaries

Abstract

We study the effects of Prandtl number $\mathit {Pr}$ and Rayleigh number $\mathit {Ra}$ in two-dimensional Rayleigh–Bénard convection without boundaries, i.e. with periodic boundary conditions. For Prandtl numbers in the range $10^{-3} \leqslant \mathit {Pr} \leqslant 10^2$, the viscous dissipation scales as $\epsilon _\nu \propto \mathit {Pr}^{1/2}\mathit {Ra}^{-1/4}$, which is based on the observation that enstrophy $\langle {\omega ^2}\rangle \propto \mathit {Pr}^0 \mathit {Ra}^{1/4}$, and the Nusselt number tends to follow the ‘ultimate’ scaling $\mathit {Nu} \propto \mathit {Pr}^{1/2}\mathit {Ra}^{1/2}$ for all values of $\mathit {Pr}$ considered. The inverse cascade of kinetic energy forms the power-law spectrum $\hat {E}_u(k) \propto k^{-2.3}$, which is close to $k^{-11/5}$ proposed by the Bolgiano–Obukhov (BO) scaling. The potential energy flux is not constant, in contrast to one of the main assumptions underlying the BO phenomenology. So, the direct cascade of potential energy forms the power-law spectrum $\hat {E}_\theta (k) \propto k^{-1.2}$, which deviates from the expected $k^{-7/5}$. Finally, at $\mathit {Pr} \to 0$ and $\infty$, we find that the dynamics is dominated by vertically oriented elevator modes that grow without bound, even at high Rayleigh numbers and with large-scale dissipation present.

1. Introduction

The fundamental challenge of understanding thermally driven flow in the strongly nonlinear regime has puzzled the fluid dynamics community for more than a century (Doering 2019). The typical set-up that is studied theoretically consists of a fluid confined between two horizontal plates that is heated from below and cooled at the top. This is the so-called Rayleigh–Bénard convection (RBC) problem after Rayleigh's proposed model (Rayleigh 1916) for Bénard's experiment on buoyancy-driven thermal convection (Bénard 1900, 1927). The classical problem depends on two dimensionless parameters: the Rayleigh number $\mathit {Ra}$, which measures the driving effect of buoyancy relative to the stabilising effects of viscosity and thermal diffusivity; and the Prandtl number $\mathit {Pr}$, which is the ratio of kinematic viscosity to thermal diffusivity. The global flow properties may be characterised by the Nusselt number $\mathit {Nu}$, a further dimensionless parameter that measures the total heat flux relative to the purely conductive heat flux. A way to progress our understanding on RBC is to understand the dependence of the dynamics on $\mathit {Ra}$, $\mathit {Pr}$ and the boundary conditions.

Bolgiano (1959) and Obukhov (1959) proposed the so-called Bolgiano–Obukhov (BO) phenomenology for stably stratified turbulence, in which buoyancy balances inertia in the momentum equation and the potential energy flux is approximately constant for length scales larger than the Bolgiano scale. These assumptions lead to kinetic and potential energy spectra of the form $\hat {E}_u(k) \propto k^{-11/5}$ and $\hat {E}_\theta (k) \propto k^{-7/5}$, respectively, where $k$ is the wavenumber. This scaling law is considered to be universal under certain conditions, irrespective of the specific details of the flow geometry, boundary conditions or forcing mechanisms. However, deviations from universality can occur. Despite originally being proposed for homogeneous stably stratified flows, BO scaling has been partly identified in experimental and numerical studies of inhomogeneous (i.e. with no-slip or free-slip boundary conditions on the top and bottom plates) three-dimensional (3-D) RBC in different geometries (Wu et al. 1990; Calzavarini, Toschi & Tripiccione 2002; Mishra & Verma 2010; Kaczorowski & Xia 2013). However, the existence of the BO scaling is still debatable (Lohse & Xia 2010; Kunnen & Clercx 2014), with some studies of inhomogeneous and homogeneous (i.e. with periodic boundary conditions) 3-D RBC reporting (Kolmogorov 1941) power-law scaling $k^{-5/3}$ for the kinetic and potential energy spectra (Borue & Orszag 1997; Verma, Kumar & Pandey 2017; Verma 2018). Similar controversy remains in simulations of two-dimensional (2-D) RBC, with some studies partly identifying BO scaling (Toh & Suzuki 1994; Mazzino 2017; Xie & Huang 2022; Samuel & Verma 2024), while some argue against the validity of the BO scaling (Biskamp & Schwarz 1997; Biskamp, Hallatschek & Schwarz 2001; Celani, Mazzino & Vergassola 2001; Celani et al. 2002).

For the strongly nonlinear regime of thermal convection, which is of paramount importance for geophysical and astrophysical applications, there are two competing theories for the behaviour of $\mathit {Nu}$ as $\mathit {Ra}$ tends to infinity for arbitrary $\mathit {Pr}$. These two proposed asymptotic scaling laws are the ‘classical’ theory $\mathit {Nu} \propto \mathit {Pr}^0\mathit {Ra}^{1/3}$ by Malkus (1954) and the ‘ultimate’ theory $\mathit {Nu} \propto \mathit {Pr}^{1/2}\mathit {Ra}^{1/2}$, at least for $\mathit {Pr} \leqslant 1$, by Kraichnan (1962) and Spiegel (1971). The Rayleigh number at which the transition to the ultimate scaling is presumed to occur is not known, and laboratory experiments and numerical simulations (Niemela et al. 2000; Chavanne et al. 2001; Niemela & Sreenivasan 2003; Urban, Musilová & Skrbek 2011; He et al. 2012; Zhu et al. 2018; Doering, Toppaladoddi & Wettlaufer 2019) have reported different power laws for a wide range of Rayleigh numbers.

In the analogy of RBC with geophysical phenomena, the top and bottom boundaries are often absent, particularly when the focus is on understanding the dynamics of the bulk flow. In this study, we choose the most obvious theoretical approach to bypass boundary layer effects by considering a fully periodic domain (Toh & Suzuki 1994; Borue & Orszag 1997; Biskamp et al. 2001; Celani et al. 2002; Ng et al. 2018) for the Rayleigh–Bénard problem, with an imposed constant vertical temperature gradient. For this homogeneous RBC set-up it has been claimed that $Nu \propto Ra^{1/2}$ (Lohse & Toschi 2003; Calzavarini et al. 2005), in line with the ultimate theory, and this scaling is also suggested from simulations of axially periodic RBC in a vertical cylinder (Schmidt et al. 2012).

Homogeneous RBC, however, exhibits exponentially growing solutions in the form of axially uniform vertical jets, called ‘elevator modes’ (Baines & Gill 1969; Batchelor & Nitsche 1991; Calzavarini et al. 2006; Liu & Zikanov 2015). As soon as these modes grow to a significant amplitude, they experience secondary instabilities, ultimately leading to statistically stationary solutions (Calzavarini et al. 2006; Schmidt et al. 2012). In recent numerical simulations, the elevator modes were suppressed by introducing an artificial horizontal buoyancy field (Xie & Huang 2022) or large-scale friction (Barral & Dubrulle 2023). The inverse cascade that is observed in 2-D homogeneous RBC (Toh & Suzuki 1994; Verma 2018; Xie & Huang 2022) is another source of energy for the large-scale modes that grow to extreme values, forming a condensate whose amplitude saturates when the viscous dissipation at the largest scale balances the energy injection (Chertkov et al. 2007; Boffetta & Ecke 2012). So, in this study, to avoid the unbounded growth of energy we include large-scale dissipation to mimic the effect of friction when there are boundaries and to be able to reach statistically stationary solutions.

The aim of this paper is to bring more insight to the effect of the control parameters and the boundaries on the dynamics of 2-D RBC. Most of the attention on 2-D RBC focuses on the Rayleigh number dependence of the dynamics, with only a few studies considering the effects of the Prandtl number (Calzavarini et al. 2005; Zhang, Zhou & Sun 2017; He et al. 2021). So, we extensively study the effects of the Prandtl number using numerical simulations of 2-D RBC driven by a constant temperature gradient in a periodic domain, to effectively remove the boundaries. We also consider high Rayleigh number simulations with normal viscosity and hyperviscosity to permit large scale separation, which is crucial for the analysis of the multi-scale dynamics, that will also shed light on the long standing debate of the power-law scalings of the spectra.

The paper is structured as follows. Section 2 contains the dynamical equations, numerical methods and the definitions of the global and spectral observables under study. The results of our simulations are presented in § 3, where we analyse how the global observables scale with the Prandtl and Rayleigh numbers and then investigate the effects of these dimensionless parameters on the spectral dynamics. In Appendix A we discuss the behaviour of the system at zero and infinite Prandtl number and in § 4 we summarise our conclusions.

2. Problem description

2.1. Governing equations

We consider 2-D RBC of a fluid heated from below in a periodic square cell $(x,y)\in [0,L]^2$. The temperature $T(x,y,t)$ is decomposed as $T=-\Delta T y/L+\theta$, where $\Delta T/L$ is the constant imposed temperature gradient and the temperature perturbation $\theta (x,y,t)$ satisfies periodic boundary conditions. As usual, for simplicity we employ the Oberbeck–Boussinesq approximation (Oberbeck 1879; Boussinesq 1903; Tritton 2012), in which the kinematic viscosity $\nu$ and the thermal diffusivity $\kappa$ are taken to be constant while the temperature dependence of the fluid density $\rho$ is neglected, except in the buoyancy term of the momentum equation.

The governing equations of the problem in two dimensions can be written in terms of $\theta (x,y,t)$ and the streamfunction $\psi (x, y, t)$ as follows:

(2.1a)$$\begin{gather} \partial_t \nabla^2 \psi + \{\psi, \nabla^2 \psi\} = \alpha g \partial_x \theta + \nu {(-\nabla^2)^{n+1}} \psi + \mu \psi, \end{gather}$$

(2.1b)$$\begin{gather}\partial_t \theta + \{\psi, \theta \} = \frac{\Delta T}{L} \partial_x \psi - \kappa {(-\nabla^2)^n} \theta, \end{gather}$$

where $\{A, B\} = \partial _x A \partial _y B - \partial _y A \partial _x B$ is the standard Poisson bracket, $\alpha$ is the thermal expansion coefficient and $g$ is the gravitational acceleration. To prevent the formation of a large-scale condensate in the presence of an inverse cascade, to avoid the elevator modes and to reach a turbulent stationary regime, we supplement our system with a large-scale dissipative term $\mu \psi$ that is responsible for saturating the inverse cascade. We consider both normal and hyperviscosity by raising the Laplacian to the powers of $n=1$ and $n=4$, respectively. The hyperviscous case, albeit not physically realisable, gives a wider inertial range, as diffusive and viscous terms kick in abruptly at much smaller scales compared with the normal viscosity case. In the ideal limit of $\nu \to 0$, $\kappa \to 0$ and $\mu \to 0$, the quantity that is conserved is $E_u(t) - ({agL}/{\Delta T})E_\theta (t)$, where the instantaneous kinetic energy $E_u(t)$ and the potential energy $E_\theta (t)$ are defined by

(2.2a,b)\begin{equation} E_u(t) =\tfrac{1}{2}\langle{|\nabla\psi|^2}\rangle_\textbf{x}, \quad E_\theta(t) = \tfrac{1}{2}\langle{\theta^2}\rangle_\textbf{x}, \end{equation}

with the angle brackets $\langle {\cdot }\rangle _\textbf{x}$ here denoting average in space. Spatio-temporal averages will be denoted as $\langle {\cdot }\rangle$.

Equations (2.1) depend on three dimensionless parameters, namely

(2.3a–c)\begin{equation} \mathit{Pr} = \frac{\nu}{\kappa}, \quad \mathit{Ra} = \frac{\alpha g\Delta T L^{4n-1}}{\nu \kappa}, \quad \mathit{Rh} = \mu \left(\frac{L^5}{\alpha g \Delta T}\right)^{1/2}, \end{equation}

which are the Prandtl number, Rayleigh number and friction Reynolds number, respectively.

We perform direct numerical simulations of (2.1) using the pseudospectral method (Orszag & Gottlieb 1977). We decompose the streamfunction into basis functions with Fourier modes in both the $x$ and $y$ directions, viz.

(2.4)\begin{equation} \psi({\boldsymbol{x}},t) = \sum^{N/2}_{{\boldsymbol{k}} ={-}N/2} \hat{\psi}_{\boldsymbol{k}}(t) {\rm e}^{{\rm i} {\boldsymbol{k}} \boldsymbol{\cdot}\boldsymbol{x}}, \end{equation}

where $\hat { \psi }_{\boldsymbol {k}}$ is the amplitude of the ${\boldsymbol {k}} = (k_x, k_y)$ mode of $\psi$, and $N$ denotes the number of aliased modes in the $x$- and $y$-directions. We decompose $\theta$ in the same way. A third-order Runge–Kutta scheme is used for time advancement and the aliasing errors are removed with the two-thirds dealiasing rule (Gómez, Mininni & Dmitruk 2005). In both the normal and hyperviscous simulations, we find that $\mathit {Rh}=(2{\rm \pi} )^{5/2}$ yields a saturated turbulent state that dissipates enough kinetic energy at large scales such that the kinetic energy spectrum peaks at $k = 2$ without overdamping the system. So, we fix $\mathit {Rh} = (2{\rm \pi} )^{5/2} \simeq 100$ while varying the Rayleigh and Prandtl numbers in the ranges $6.2 \times 10^7 \leqslant \mathit {Ra} \leqslant 6.2 \times 10^{11}$ and $10^{-3} \leqslant \mathit {Pr} \leqslant 10^2$. To model the large Rayleigh number dynamics, we set $\mathit {Ra} = 9.4 \times 10^{49}$ in our hyperviscous simulations. Figure 1 shows the parameter values simulated in the $(\mathit {Ra},\mathit {Pr})$-plane as well as the resolution, $N$, used in each case. The criterion we used to ensure that our runs are well resolved was to verify that the peaks of the potential and kinetic energy dissipation spectra are at least one order of magnitude above the tail of the spectra. Time-averaged quantities are computed after the system has reached a statistically stationary regime. We have accumulated statistics over long time integrations to ensure that our statistics are well converged.

Figure 1. Parameter values simulated in the $(\mathit {Ra},\mathit {Pr})$-plane with the resolution, $N$, used at each instant colour coded in the legend. The black dashed line separates runs with normal viscosity (below line) from those with hyperviscosity (above line).

2.2. Global and spectral observables

Next, we briefly outline the global and spectral observables that will be explored in our numerical simulations below. The energy spectra of the velocity field $\hat {E}_u(k,t)$ and the temperature field $\hat {E}_{\theta }(k,t)$, referred to as the kinetic energy and potential energy spectra, are defined as

(2.5a)$$\begin{gather} \hat{E}_u(k,t) = \frac{1}{2}\sum_{k \leqslant |{\boldsymbol{k}}| < k+\Delta k} k^2\left| \hat{\psi}_{\boldsymbol{k}}(t) \right|^2, \end{gather}$$

(2.5b)$$\begin{gather}\hat{E}_{\theta}(k,t) = \frac{1}{2}\sum_{k \leqslant |{\boldsymbol{k}}| < k+\Delta k} \left| \hat{\theta}_{\boldsymbol{k}}(t) \right|^2, \end{gather}$$

where the sum is performed over the Fourier modes with wavenumber amplitude $k = |{\boldsymbol {k}}| = \sqrt {k_x^2 + k_y^2}$ in a shell of width $\Delta k = 2{\rm \pi} /L$. Using the Fourier transform, one can derive the evolution equations of kinetic and potential energy spectra from (2.1), namely

(2.6a)$$\begin{gather} \partial_t \hat{E}_u(k,t) ={-} \partial_k \varPi_u(k,t) - D_\nu(k,t) - D_\mu(k,t) + \alpha g F_B(k,t), \end{gather}$$

(2.6b)$$\begin{gather}\partial_t \hat{E}_{\theta}(k,t) ={-}\partial_k \varPi_\theta(k,t) - D_\kappa(k,t) + \frac{\Delta T}{L}F_B(k,t). \end{gather}$$

The energy flux $\varPi$ is a measure of the nonlinear cascades in turbulence (Alexakis & Biferale 2018). The energy flux for a circle of radius $k$ in the 2-D wavenumber space is the total energy transferred from the modes within the circle to the modes outside the circle. Consequently, we define the flux of kinetic energy $\varPi _u(k,t)$ and potential energy $\varPi _{\theta }(k,t)$ as

(2.7a)$$\begin{gather} \varPi_u(k,t) = \sum_{k' \leqslant k} T_u(k',t), \end{gather}$$

(2.7b)$$\begin{gather}\varPi_{\theta}(k,t) = \sum_{k' \leqslant k} T_\theta(k',t), \end{gather}$$

where $T_u(k,t)$ and $T_\theta (k,t)$ are the nonlinear kinetic and potential energy transfer across $k$

(2.8a)$$\begin{gather} T_u(k,t) ={-}\sum_{k \leqslant |{\boldsymbol{k}}| < k+\Delta k} \hat{\psi}^*_{\boldsymbol{k}}(t) \widehat{\{\psi,\nabla^2 \psi \}}_{\boldsymbol{k}}(t), \end{gather}$$

(2.8b)$$\begin{gather}T_\theta(k,t) = \sum_{k \leqslant |{\boldsymbol{k}}| < k+\Delta k} \hat{\theta}^*_{\boldsymbol{k}}(t) \widehat{\{\psi,\theta \}}_{\boldsymbol{k}}(t). \end{gather}$$

The notation $\widehat {\{ . \}}_{\boldsymbol {k}}$ represents the Fourier mode of the Poisson bracket expanded using (2.4), and the asterisk denotes complex conjugation.

The spectra of the small-scale viscous dissipation $D_\nu (k,t)$, the large-scale frictional dissipation $D_\mu (k,t)$ and the thermal dissipation $D_\kappa (k,t)$ are defined as

(2.9a)$$\begin{gather} D_\nu(k,t) = 2 \nu k^{2n}\hat{E}_u(k,t), \end{gather}$$

(2.9b)$$\begin{gather}D_\mu(k,t) = 2 \mu k^{{-}2} \hat{E}_u(k,t), \end{gather}$$

(2.9c)$$\begin{gather}D_\kappa(k,t) = 2 \kappa k^{2n} \hat{E}_\theta(k,t), \end{gather}$$

and the spectrum of the buoyancy term is given by

(2.10)\begin{equation} F_B(k,t) = \sum_{k \leqslant |{\boldsymbol{k}}| < k+\Delta k} {\rm i} k_x \widehat{\psi^*}_{\boldsymbol{k}}(t) \hat{\theta}_{\boldsymbol{k}}(t). \end{equation}

The Nusselt number is a dimensionless measure of the averaged vertical heat flux relative to the conductive heat flux. The form of the conductive heat flux is determined by the thermal diffusive term. Hence, the Nusselt number is defined mathematically by

(2.11)\begin{equation} \mathit{Nu} = 1 + \frac{\langle{\partial_x \psi \theta}\rangle}{\kappa \Delta T / L^{2n-1}}. \end{equation}

Using the above definition, one can derive the following exact relations for the kinetic and potential energy balances in the statistically stationary regime, as in Shraiman & Siggia (1990):

(2.12a)$$\begin{gather} \epsilon_u = \epsilon_\nu + \epsilon_\mu = \frac{\nu^3}{L^{6n-2}}(\mathit{Nu} - 1)\frac{\mathit{Ra}}{\mathit{Pr}^2}, \end{gather}$$

(2.12b)$$\begin{gather}\epsilon_\theta = \epsilon_\kappa = \frac{\kappa \Delta T^2}{L^{2n}} (\mathit{Nu} - 1), \end{gather}$$

where $\epsilon _u = \alpha g \langle {\psi \partial _x \theta }\rangle$ is the injection rate of kinetic energy due to buoyancy, $\epsilon _\nu = \nu {\langle \psi \nabla ^{2(n+1)} \psi \rangle }$ is the viscous dissipation rate, $\epsilon _\mu = \mu {\langle \psi ^2\rangle }$ is the large-scale dissipation rate, $\epsilon _\theta = ({\Delta T}/{L}) \langle {\partial _x \psi \,\theta }\rangle$ is the injection rate of potential energy due to buoyancy and $\epsilon _\kappa = \kappa {\langle \theta \nabla ^{2n} \theta \rangle }$ is the thermal dissipation rate.

2.3. Elevator modes

Upon linearising (2.1) about the conductive state ($\psi =\theta =0$), we find that infinitesimal solutions with $\psi (\boldsymbol {x},t)=\exp ({{\rm i} \boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x}+\sigma t})$ are possible, provided the normalised linear growth rate $\sigma$ satisfies the relation

(2.13)\begin{equation} \left(\sigma+k^{2n}\right) \left(\sigma+\mathit{Pr} k^{2n}+\frac{\mathit{Rh}\sqrt{\mathit{Ra}\mathit{Pr}}}{k^2}\right) =\frac{\mathit{Ra}\mathit{Pr} (2{\rm \pi} k_x)^2}{k^2}, \end{equation}

which has two real roots for $\sigma$. From these two roots, one is positive if and only if

(2.14)\begin{equation} \mathit{Rh}\sqrt{\frac{\mathit{Ra}}{\mathit{Pr}}}<\frac{\mathit{Ra}(2{\rm \pi} k_x)^2}{k^{2n}}-k^{2n+2}. \end{equation}

One can show that $\sigma$ is a monotonic decreasing function of $k_y$, so the most dangerous modes are independent of $y$, with

(2.15)\begin{equation} \psi(\boldsymbol{x},t)=\exp({{\rm i}k_x x+\sigma t}),\quad k_x\in\mathbb{Z}_{{>}0}. \end{equation}

Indeed, such a unidirectional mode satisfies the nonlinear governing equations (2.1) exactly, because the nonlinear Poisson bracket terms are identically zero.

Although the maximum growth rate does not necessarily occur at the minimum wavenumber $k_x=1$, it is straightforward to show that (2.14) is satisfied for some $(k_x,k_y)$ if and only if it is satisfied at $(k_x,k_y)=(1,0)$. In other words, exact solutions of the problem (2.1) that grow exponentially without bound exist whenever

(2.16)\begin{equation} \mathit{Rh}\sqrt{\frac{\mathit{Ra}}{\mathit{Pr}}}< \frac{\mathit{Ra}}{(2{\rm \pi})^{2(n-1)}}-(2{\rm \pi})^{2n+2}. \end{equation}

In practice, unless the initial conditions are equal to an elevator mode, we find that elevator modes do not dominate the dynamics at finite Prandtl number if $\mathit {Rh}$ is large enough. So, we fix $\mathit {Rh}=(2{\rm \pi} )^{5/2}$, which is sufficient to prevent elevator modes from growing and to allow the flow to reach turbulent saturated states without overdamping the system. For zero and infinite Prandtl number convection, elevator modes are seemingly more attractive, even for randomly generated initial conditions (see Appendix A). Thus, our analysis focuses on finite values of Prandtl number.

3. Results

3.1. Global observables

In figures 2(a)–2(c) we show how the global observables in the kinetic and potential energy balances (2.12) vary with $\mathit {Pr}$ keeping $\mathit {Ra} = 6.2 \times 10^{11}$, while in (d–f) we show how the same quantities vary with $\mathit {Ra}$ while keeping $\mathit {Pr}=1$. In both cases, we use normal viscosity ($n=1$) and keep $\mathit {Rh}=(2{\rm \pi} )^{5/2}$ constant. Firstly, figures 2(a) and 2(d) show that the kinetic and potential energies remain virtually constant while $\mathit {Pr}$ and $\mathit {Ra}$ vary by at least four orders of magnitude. Secondly, figures 2(b) and 2(e) show that $\epsilon _u \approx \epsilon _\mu$, which indicates that the majority of the kinetic energy, injected by buoyancy in the flow, is dissipated at large scales. This effect is caused by the inverse cascade of kinetic energy, which will be investigated in more detail in § 3.2. Note that $\epsilon _u$ and $\epsilon _\mu$ are almost independent of both $\mathit {Pr}$ and $\mathit {Ra}$, while the viscous dissipation scales like

(3.1)\begin{equation} \epsilon_\nu \propto \mathit{Pr}^{1/2} \mathit{Ra}^{{-}1/4}. \end{equation}

Finally, in figures 2(c) and 2(f), we see that $\epsilon _\theta =\epsilon _\kappa$, as required by the potential energy balance (2.12b), and both observables, are also virtually independent of both $\mathit {Pr}$ and $\mathit {Ra}$.

Figure 2. Spatio-temporal averages of the kinetic and potential energies and the terms in the kinetic and potential energy balances (2.12) as functions of $\mathit {Pr}$ with $\mathit {Ra} = 6.2 \times 10^{11}$ fixed (a–c) and as functions of $\mathit {Ra}$ with $\mathit {Pr} = 1$ fixed (d–f). We use normal viscosity ($n=1$) and keep $\mathit {Rh} = (2{\rm \pi} )^{5/2}$ fixed in all plots.

To explain the scaling of the viscous dissipation rate (3.1), we now look at how the enstrophy $\langle {\omega ^2}\rangle$ varies with the Prandtl and Rayleigh numbers, where $\omega = \nabla ^2 \psi$ is the vorticity of the flow. In figure 3(a) we plot enstrophy as a function of $\mathit {Pr}$ with $\mathit {Ra} = 6.2 \times 10^{11}$ fixed, and in figure 3(b) as a function of $\mathit {Ra}$ with $\mathit {Pr} = 1$ fixed. We keep $\mathit {Rh} = (2{\rm \pi} )^{5/2}$ fixed throughout. From figure 3(a) we observe that enstrophy can be considered approximately independent of $\mathit {Pr}$, because it varies by less than a factor of two over the five decades of Prandtl numbers considered, while figure 3(b) demonstrates that enstrophy scales like $\mathit {Ra}^{1/4}$, i.e.

(3.2)\begin{equation} \langle{\omega^2}\rangle \propto \mathit{Pr}^0 \mathit{Ra}^{1/4}. \end{equation}

With normal viscosity, by definition we have $\epsilon _\nu = \nu \langle {\psi \nabla ^4 \psi }\rangle = \nu \langle {\omega ^2}\rangle$ and, writing $\nu \propto (\mathit {Pr}/\mathit {Ra})^{1/2}$, we thus obtain the scaling of (3.1) that is observed in figure 2.

Figure 3. Enstrophy $\langle {\omega ^2}\rangle$ as a function of $\mathit {Pr}$ with $\mathit {Ra} = 6.2 \times 10^{11}$ fixed (a), and as a function of $\mathit {Ra}$ with $\mathit {Pr} = 1$ fixed (b). We keep $n=1$ (normal viscosity) and $\mathit {Rh} = (2{\rm \pi} )^{5/2}$ fixed in all plots.

According to figure 2, the energy injection rates $\epsilon _u$ and $\epsilon _\theta$ are both approximately constant in the range of $\mathit {Pr}$ and $\mathit {Ra}$ we considered. With $\mathit {Ra}\mathit {Pr} \gg 1$ and $n = 1$, the net energy balances (2.12) thus produce the Nusselt number scaling

(3.3)\begin{equation} \mathit{Nu} \propto \mathit{Pr}^{1/2}\mathit{Ra}^{1/2}. \end{equation}

This relation is in agreement with the ultimate scaling. In figure 4, we plot the Nusselt number compensated by the classical scaling, $\mathit {Nu} \propto \mathit {Pr}^{0}\mathit {Ra}^{1/3}$, and by the ultimate scaling, $\mathit {Nu} \propto \mathit {Pr}^{1/2}\mathit {Ra}^{1/2}$. The ultimate scaling provides a much more convincing collapse of the data, with the fit becoming progressively more accurate as $\mathit {Ra}$ increases for all values of $Pr$ considered. Indeed, the ultimate scaling in terms of Rayleigh number dependence was expected to be exhibited by our simulations, as we have effectively removed the boundary layers by applying periodic boundary conditions, similarly to previous studies (Lohse & Toschi 2003; Calzavarini et al. 2005). The ultimate scaling asserts (3.3) at least for $\mathit {Pr} \leqslant 1$ and this is what is expected in 3-D RBC with boundaries. Nevertheless, figure 4 demonstrates that in homogeneous 2-D RBC the Nusselt number depends on all the values of Prandtl number we considered, even for $\mathit {Pr} > 1$.

Figure 4. The Nusselt number compensated with the ‘classical’ and the ‘ultimate’ scaling in (a,b), respectively. The Prandtl numbers are colour coded in the legend.

3.2. Spectra and cascades

In this section, we examine the multi-scale dynamics of the flows using spectra. To eliminate finite Rayleigh number effects and to have large enough scale separation, we focus on the hyperviscous simulations (i.e. $n=4$ in (2.1)) at $\mathit {Ra} = 9.4 \times 10^{49}$ and $\mathit {Rh} = (2{\rm \pi} )^{5/2}$. Results with normal viscosity (i.e. $n=1$ in (2.1)) at $\mathit {Ra} = 6.2 \times 10^{11}$ are also presented for comparison. According to the BO scaling (Bolgiano 1959; Obukhov 1959), the ratio of the kinetic to the potential energy spectra scales as $\hat {E}_u(k) / \hat {E}_\theta (k) \propto k^{-4/5}$, since $\hat {E}_u(k) \propto k^{-11/5}$ and $\hat {E}_\theta (k) \propto k^{-7/5}$. In figure 5 we plot $\hat {E}_u(k) / \hat {E}_\theta (k)$ compensated by $k^{4/5}$ for the hyperviscous simulations and, in the inset, for the runs with normal viscosity. Instead of finding a wavenumber range where this scaling is valid, we observe a $k^{-0.3}$ power law for all Prandtl numbers considered, suggesting that BO scaling is not followed in our simulations. For the normal viscosity simulations we find the $k^{-0.3}$ power law again to be followed, but within a narrower wavenumber range (see inset of figure 5). Figure 5 also shows that, for $\mathit {Pr} \ll 1$, the kinetic energy is much larger than the potential energy at large wavenumbers. This is expected as the small scales are dominated by thermal diffusivity when $\mathit {Pr} \ll 1$, and so the potential energy is dissipated much more effectively than the kinetic energy. The opposite is true for $\mathit {Pr} \gg 1$.

Figure 5. The ratio of the kinetic energy to the potential energy spectra $\hat {E}_u(k)/\hat {E}_\theta (k)$ compensated by $k^{4/5}$ for the hyperviscous simulations with $\mathit {Ra} = 9.4 \times 10^{49}$, $\mathit {Rh} = (2{\rm \pi} )^{5/2}$ and different Prandtl numbers as colour coded in the legend. The inset shows the same ratio but for runs with normal viscosity at $\mathit {Ra} = 6.2 \times 10^{11}$.

In figures 6(a) and 6(b), we plot the kinetic and potential energy spectra compensated by the power laws proposed by the BO phenomenology. The hyperviscous runs are shown in the main plots, and normal viscosity runs in the insets. From the hyperviscous runs we observe that the compensated kinetic energy spectrum $k^{11/5}\hat {E}_u(k)$ is flat within a range of wavenumbers and hence it is close to the BO scaling, with the best fit being $\hat {E}_u(k) \propto k^{-2.3}$. However, the BO scaling is not evident from the inset showing the runs with normal viscosity due finite Rayleigh number effects. On the other hand, it is clear that the compensated potential energy spectrum $k^{7/5}\hat {E}_\theta (k)$ does not follow the BO scaling for either the hyperviscous and normal viscosity runs, with the best fit being $\hat {E}_\theta (k)\propto k^{-1.2}$. The power-law exponents of the spectra we observe differ from those in 3-D RBC with periodic boundary conditions, where the kinetic and potential energy exhibit $k^{-5/3}$ power-law spectra (Borue & Orszag 1997), similar to those observed in passive scalar turbulence (Warhaft 2000).

Figure 6. The time-averaged energy spectra compensated by best-fit power laws (a,b), and spectral fluxes normalised by the time-averaged dissipation rates (c,d) for the runs with hyperviscosity at $\mathit {Ra} = 9.4 \times 10^{49}$, $\mathit {Rh} = (2{\rm \pi} )^{5/2}$ and different Prandtl numbers as colour coded in the legend. The insets show the same quantities but for runs with normal viscosity at $\mathit {Ra} = 6.2 \times 10^{11}$.

To understand the dynamics of the turbulent cascades, we plot the associated kinetic and potential energy fluxes normalised by the respective time-averaged injection rates of energy due to buoyancy in figures 6(c) and 6(d). The positive potential energy flux suggests a strong direct cascade, while the weak negative kinetic energy flux suggests an inverse cascade. For $\mathit {Pr} = 1$ these types of cascades are in agreement with Xie & Huang (2022). The negative kinetic energy flux peaks at low wavenumbers, while the potential energy flux peaks at high wavenumbers. The inverse cascade of kinetic energy is not affected by the Prandtl number; however, the direct cascade of potential energy moves to higher wavenumbers along with the peak of $\varPi _\theta (k)$ as $\mathit {Pr}$ increases. We emphasise the wavenumber dependence of kinetic and potential energy fluxes, with $\partial _k \varPi _u(k) > 0$ and $\partial _k \varPi _\theta (k) > 0$ in the intermediate range of wavenumbers for all of the Prandtl numbers considered.

In figure 7, we present the time-averaged spectra of the magnitudes of the terms in the kinetic and potential energy balances (2.6) for runs with hyperviscosity and three different values of $\mathit {Pr}\in \{10^{-2},1,10^2\}$. The corresponding plots with normal viscosity and $\mathit {Ra}=6.2\times 10^{11}$ are shown in figure 8. In both figures, the red dots indicate where the inertial flux terms $\partial _k\varPi _u$ and $\partial _k\varPi _\theta$ become negative. In the kinetic energy balance, we identify three distinct wavenumber ranges, labelled I to III in the plots. In region I, the kinetic energy injected by buoyancy is dissipated by the large-scale friction. In region II, the inertial term balances buoyancy, which is positive for all wavenumbers in this region, i.e.

(3.4)\begin{equation} \partial_k \varPi_u(k) \approx \alpha g F_B(k) > 0, \end{equation}

with $F_B(k) \propto k^{-1.5}$. Equation (3.4) implies the $k$ dependence of the kinetic energy flux we see in figure 6(c), which is expected based on the BO phenomenology. Note that region II is largest for small Prandtl numbers, especially in the runs with normal viscosity, as evident from the results presented in figures 8(a)–8(c). In region III, the balances between terms depend on the Prandtl number. For $\mathit {Pr} = 10^{-2}$ buoyancy decays rapidly and so small-scale viscous dissipation is balanced by the inertial term, which is negative in this range of wavenumbers. As $\mathit {Pr}$ increases, buoyancy becomes more significant in the balance of region III between the small-scale viscous dissipation and the inertial term. For $\mathit {Pr} \gtrsim 10^2$, small-scale viscous dissipation seems to be balanced by buoyancy rather than by the inertial term. This effect is shown more clearly in the runs with normal viscosity shown in figure 8(c), where the small-scale dissipation range is much larger.

Figure 7. Time-averaged spectra of the terms in the kinetic energy balance (2.6a) (a–c); and the potential energy balance (2.6b) (d–f) for the runs with hyperviscosity, ${\mathit {Ra} = 9.4 \times 10^{49}}$, $\mathit {Rh} = (2{\rm \pi} )^{5/2}$ with $\mathit {Pr} = 10^{-2}$ (a,d), $\mathit {Pr} = 1$ (b,e) and $\mathit {Pr} = 10^{2}$ (c,f). The terms displayed in the legends are defined in (2.7), (2.9) and (2.10). In (a–c) and (d–f), we observe three and two distinct dominant balances, respectively, annotated I–III and A–B. The red dots indicate where the inertial terms become negative.

Figure 8. Time-averaged spectra of the terms in the kinetic energy balance (2.6a) (a–c); and the potential energy balance (2.6b) (d–f) for the runs with normal viscosity, ${\mathit {Ra} = 6.2 \times 10^{11}}$, $\mathit {Rh} = (2{\rm \pi} )^{5/2}$ with $\mathit {Pr} = 10^{-2}$ (a,d), $\mathit {Pr} = 1$ (b,e) and $\mathit {Pr} = 10^{2}$ (c,f). The terms displayed in the legend are defined in (2.7), (2.9) and (2.10). In (a–c) and (d–f), we observe three and two distinct dominant balances, respectively, annotated I–III and A–B. The red dots indicate where the inertial terms become negative.

In the potential energy balance, we identify two distinct wavenumber ranges, labelled A and B in the plots (see (d–f) in figures 7 and 8). In these plots we observe the inertial term to be balanced by buoyancy in region A and by small-scale thermal dissipation in region B. In other words, the potential energy injected by buoyancy in region A is cascaded to larger wavenumbers, where it is dissipated by thermal diffusivity.

We recall that the red dots in figures 7 and 8 show where the inertial terms of the kinetic and potential energy become negative. This sign change occurs primarily in region III for $\partial _k \varPi _u(k)$ and region B for $\partial _k \varPi _\theta (k)$, the latter corresponding to the large negative gradient of $\varPi _\theta (k)$ observed in figure 6 at large wavenumbers. Near the boundary between regions A and B in figures 7 and 8(d–f), we observe that $\partial _k \varPi _\theta (k)$ exhibits fluctuations between positive and negative values over these wavenumbers. However, for the majority of the wavenumbers in region A we find that

(3.5)\begin{equation} \partial_k \varPi_\theta(k) \approx \frac{\Delta T}{L} F_B(k) > 0. \end{equation}

This relation implies the $k$ dependence of the potential energy flux we observe in figure 6(d) for all the Prandtl numbers we studied. Note that in the wavenumber range where region II and A overlap, we have $\partial _k \varPi _u(k)/\alpha g \approx ({L}/{\Delta T}) \partial _k \varPi _\theta (k) \approx F_B(k)$, such that $\varPi _u(k) - ({\alpha g L}/{\Delta T}) \varPi _\theta (k)$ is approximately constant in the intermediate range of wavenumbers as $F_B(k)$ cancels out. This is the inertial range of scales, where the viscous, diffusive and friction effects can be neglected and so the energy flux is constant. This is expected for the ideal invariant $E_u(t) - ({agL}/{\Delta T} )E_\theta (t)$, which is conserved in the limit of $\nu \to 0$, $\kappa \to 0$ and $\mu \to 0$. In Appendix B we plot the absolute value of the time-averaged spectra $|\hat {E}_u(k) - ({agL}/{\Delta T}) \hat {E}_\theta (k)|$ and the flux $\varPi _u(k) - ({\alpha g L}/{\Delta T}) \varPi _\theta (k)$ for the hyperviscous runs at $Ra = 9.4 \times 10^{49}$, $Rh = (2{\rm \pi} )^{5/2}$ and different Prandtl numbers as colour coded in the legend (see figure 10). The spectra $|\hat {E}_u(k) - ({agL}/{\Delta T} )\hat {E}_\theta (k)|$ scale like $k^{-1}$ and the flux $\varPi _u(k) - ({\alpha g L}/{\Delta T} )\varPi _\theta (k)$ is constant as expected for all values of $Pr$ considered.

4. Conclusions

In this paper we study the effects of varying the Prandtl and Rayleigh numbers on the dynamics of 2-D RBC without boundaries, i.e. with periodic boundary conditions. For zero and infinite Prandtl number convection, we find that large-scale elevator modes dominate the dynamics, unless large-scale dissipation is made so strong as to suppress turbulent convection. Such elevator modes have long been known (Baines & Gill 1969). In all parameter values simulated, the inequality (2.16) holds, implying that the system admits exact single-mode solutions that grow exponentially without bound. Whether or not the solution blows up must depend on the initial conditions, in a way that is not currently understood. Instead, we focus on finite Prandtl numbers, where we find that nonlinear interactions between modes allow the system to reach turbulent stationary states.

Examining the Prandtl and Rayleigh number dependences of the terms in the kinetic and potential energy balances, we find that the enstrophy scales as $\langle {\omega ^2}\rangle \propto \mathit {Pr}^0 \mathit {Ra}^{1/4}$ and hence the small-scale viscous dissipation scales as $\epsilon _\nu = \nu \langle {\omega ^2}\rangle \propto \mathit {Pr}^{1/2}\mathit {Ra}^{-1/4}$. On the other hand, we observe that the injection rate of kinetic energy $\epsilon _u$ due to buoyancy is effectively independent of both the Prandtl and the Rayleigh number. Using this observation, we find that $\mathit {Nu} \propto \mathit {Pr}^{1/2}\mathit {Ra}^{1/2}$, which agrees with the so-called ultimate scaling. The collapse of the data using the ultimate scaling becomes progressively more accurate as $\mathit {Ra}$ increases for all values of $\mathit {Pr}$ considered.

Analysing the kinetic and potential energy spectral fluxes, we find an inverse cascade of kinetic energy and a direct cascade of potential energy. The inverse cascade is independent of the Prandtl number, while the peak of the potential energy flux moves to higher wavenumbers as $\mathit {Pr}$ increases. The kinetic and potential energy fluxes, $\varPi _u(k)$ and $\varPi _\theta (k)$, are not constant because both of their derivatives in $k$ are balanced by the buoyancy term $F_B(k)$, which is predominantly positive in the intermediate range of wavenumbers. These two balances imply a positive slope in $k$ for the fluxes. On the other hand, the potential energy flux was recently claimed to be constant in the inhomogeneous set-up of 2-D RBC (Samuel & Verma 2024). Nevertheless, their observations of the potential energy spectrum did not follow the BO scaling. Although we observe no range of wavenumbers where either $\varPi _u(k)$ or $\varPi _\theta (k)$ is constant, we find that they are connected by the relation $\varPi _u(k) - ({\alpha g L}/{\Delta T}) \varPi _\theta (k) \approx \text {constant}$ in the overlap between regions II and A.

The kinetic energy spectra scale as $\hat {E}_u(k) \propto k^{-2.3}$, which is close to $k^{-11/5}$ behaviour of the BO phenomenology. However, the potential energy spectra scale as $\hat {E}_\theta (k) \propto k^{-1.2}$, which deviates from the $k^{-7/5}$ scaling predicted by the BO arguments. The deviation from the BO phenomenology is also clear when we test the scaling $\hat {E}_u(k) / \hat {E}_\theta (k) k^{4/5} \approx \text {const.}$, which is clearly not followed by our spectra, which scale as $\hat {E}_u(k) / \hat {E}_\theta (k) k^{4/5} \propto k^{-0.3}$. For the hyperviscous simulations, the observed power laws in the inertial range of the kinetic and potential energy spectra do not show any dependence on the Prandtl number. The only dependence we observe is at the dissipative range of wavenumbers, where the viscosity and the thermal diffusivity dominate the dynamics. In the spectra from the normal viscosity simulations, the effects of the Prandtl number are more significant due to the comparatively low scale separation. Hence, the wavenumber range over which a power-law behaviour can be observed is truncated.

This study clearly demonstrates the necessity for large scale separation to be able to make clearer conclusions on the spectral dynamics and the power-law exponents of 2-D RBC. This requirement makes similar studies in three dimensions more challenging. The development of a phenomenology where buoyancy acts as a broadband spectral forcing is required to interpret the current observations. Numerical simulations at Prandtl and Rayleigh numbers outside the ranges we investigated, i.e. $\mathit {Pr} < 10^{-3}$, $\mathit {Pr} > 10^2$ and $\mathit {Ra} > 10^{12}$, are challenging but would be of great interest to see if they agree with the hyperviscous simulations we have performed and to provide a more complete picture of the asymptotic regime of buoyancy-driven turbulent convection.