Steady Rayleigh–Bénard convection between no-slip boundaries

Rayleigh–Bénard convection (RBC) is the buoyancy-driven flow in a fluid layer heated from below and cooled from above in the presence of gravity. The emergent convective flow enhances heat flux from the warm bottom boundary to the cool top boundary beyond the conductive flux from diffusion alone. This dimensionless enhancement factor—the ratio of bulk-averaged vertical heat flux from both conduction and convection to the flux from conduction alone—defines the Nusselt number $Nu$. In Rayleigh’s mathematical model (Rayleigh, 1916) $Nu$ depends on several dimensionless quantities characterizing the problem at hand: (i) what we now call the Rayleigh number $Ra$, which is proportional to the imposed temperature drop across the layer, (ii) the fluid’s Prandtl number $Pr$, which is the ratio of kinematic viscosity to thermal diffusivity, and (iii) details of the spatial domain, often captured by an aspect ratio $\Gamma$ that is a ratio of a horizontal length scale to the vertical layer height.

Convection is coherent at $Ra$ values not too far above the critical value $Ra_{c}$ beyond which the conductive no-flow state is linearly unstable. By coherent we mean flows with few scales present; spatial scales might include a horizontal period and the vertical thickness of boundary layers, and temporally the flow may be steady or time-periodic. Meanwhile, convection is turbulent at the large $Ra$ values pertinent to many engineering and scientific applications. Turbulent flows are complex and contain a range of spatial and temporal scales and, in the present context, have thermal and viscous boundary layers at the top and bottom boundaries from which thermal plumes emerge and mix the bulk. In a given domain it is expected that a scaling of $Nu$ with respect to both $Pr$ and $Ra$ will emerge in the $Ra\to\infty$ limit (Kadanoff, 2001).

After nearly a century of increasingly sophisticated mathematical analysis, increasingly resolved direct numerical simulations (DNS), and increasingly refined laboratory experiments, two quantitatively distinct conjectures remain in contention for the heat transport scaling law at large $Ra$ (Chillà & Schumacher, 2012; Doering, 2020). The two conjectures follow from heuristic physical arguments that both seem plausible but give incompatible predictions: the ‘classical’ scaling $Nu\sim Pr^{0}Ra^{1/3}$ and the ‘ultimate’ scaling $Nu\sim Pr^{1/2}Ra^{1/2}$, with the latter sometimes including logarithmic-in-$Ra$ modifications.

For RBC between flat, no-slip, isothermal boundaries, rigorous analysis of the governing equations has yielded upper bounds of the form $Nu\leq\mathcal{O}(Ra^{1/2})$ uniformly in $Pr$ and $\Gamma$ (Howard, 1963; Doering & Constantin, 1996), but this still allows for either classical or ultimate scaling. Upper bounds that rule out ultimate scaling by being asymptotically smaller than ${\cal O}(Ra^{1/2})$ have been derived in the limit of infinite $\Pran$ (Doering et al., 2006; Otto & Seis, 2011; Whitehead & Doering, 2012) and for two-dimensional convection between stress-free boundaries (Whitehead & Doering, 2011). For the no-slip boundaries relevant to experiments, however, it remains an open question whether an upper bound asymptotically smaller than $Ra^{1/2}$ is possible.

In view of the problem’s stubbornness, a new strategy is called for to determine—or at least to bound—$Nu$ as a function of $Ra,Pr$, and $\Gamma$. Toward that end we have undertaken an indirect approach consisting of two parts. The first part is to study coherent flows for which one can reasonably hope to determine asymptotic heat transport, and the second part is to investigate how transport by those coherent flows compares with transport by turbulent convection. The simplest coherent flows are steady—i.e., time-independent—solutions of the equations of motion. Many such states exist, although they are generally unstable at large $Ra$. We focus on what might be called the simplest type of steady states: two-dimensional convection rolls like the counter-rotating pairs shown in figure 1($a$, $b$). In horizontally periodic or infinite domains in two or three dimensions, such rolls bifurcate supercritically from the conductive state in the linear instability identified by Rayleigh (1916). A roll pair of any width-to-height aspect ratio $\Gamma$ admitted by the domain exists for sufficiently large $Ra$.

For steady rolls, the dependence of $Nu$ on the parameters $(\Gamma,Pr,Ra)$ at asymptotically large $Ra$ is accessible to computation. As for whether heat transport by steady rolls can be connected to transport by turbulence, there are several reasons for optimism. Relationships between turbulent attractors and the unstable coherent states embedded therein have been established in models of wall-bounded shear flows (Graham & Floryan, 2021), where particular steady states, traveling waves, and time-periodic states have been found that closely reflect turbulent flows in terms of integral quantities as well as particular flow structures. Analogous study of RBC began only recently but indeed suggests that certain steady states capture qualitative aspects of turbulent convection (Waleffe et al., 2015; Sondak et al., 2015; Kooloth et al., 2021; Motoki et al., 2021). Our findings add to this evidence. The desire to understand and perhaps strengthen the mathematical bound $Nu\leq\mathcal{O}(Ra^{1/2})$ is further motivation for studying unstable states since bounds apply to all solutions of the governing equations regardless of stability. It is an open question whether any solutions can achieve ultimate scaling, let alone turbulent solutions.

Here we report numerical computations of steady convection rolls for a $Pr=1$ fluid contained between no-slip isothermal top and bottom boundaries. We reach sufficiently large $Ra$ values to convincingly reveal several asymptotic scalings of $Nu$, depending on the horizontal periods of the rolls. These are the first clearly asymptotic scalings found for any type of flow—steady, turbulent, or otherwise—for RBC in the no-slip case. Notably, the largest heat transport among steady rolls of all horizontal periods displays the classical $Nu\sim Ra^{1/3}$ scaling. We further observe that $Nu$ for these steady rolls is larger than turbulent $Nu$ from all laboratory experiments and two- or three-dimensional (2D or 3D) simulations at comparable parameters. This observation supports the conjecture that steady states maximize $Nu$ among all stable or unstable flows, as was recently verified for a truncated model of RBC (Olson et al., 2021) using methods that are not yet applicable to the full governing equations. If steady-roll transport continues to dominate turbulent transport as $Ra\to\infty$, then our finding of classical scaling for steady rolls would rule out ultimate scaling of turbulent convection.

The asymptotic scaling of steady rolls is already known in the case of stress-free velocity conditions at the top and bottom boundaries, which were considered for mathematical convenience in Rayleigh’s original work. In that case $Nu\sim Ra^{1/3}$ as $Ra\to\infty$ at fixed $Pr$ and $\Gamma$, and the aspect ratio of the roll pair maximizing $Nu$ at each $Ra$ and $Pr$ approaches $\Gamma\approx 1.9$ (Chini & Cox, 2009; Wen et al., 2020). Recent computations of steady rolls in the no-slip case for pre-asymptotic $Ra$ values up to $10^{9}$ revealed significant differences from the stress-free problem (Waleffe et al., 2015; Sondak et al., 2015). The dependence $Nu(\Gamma)$ for no-slip rolls at fixed $Ra$ and $Pr$ can have multiple local maxima, as shown in figure 1($c$), and the aspect ratio $\Gamma^{*}$ that globally maximizes $Nu(\Gamma)$ approaches zero rather than a constant as $Ra\to\infty$. Steady rolls of $Nu$-maximizing aspect ratios $\Gamma^{*}$ were reported in Sondak et al. (2015) for $Ra\in[5\times 10^{6},3\times 10^{8}]$ at $Pr=1$, yielding fits of $\Gamma^{*}\sim Ra^{-0.217}$ and $Nu(\Gamma^{*})\sim Ra^{0.31}$. This heat transport scaling is faster than with $\Gamma$ fixed: computations in Waleffe et al. (2015) for $Ra\in[5\times 10^{5},5\times 10^{6}]$ at $Pr=7$ with $\Gamma=2$ fixed yield the fit $Nu\sim Ra^{0.28}$. These best-fit scaling exponents are, however, not asymptotic.

Steady convection rolls are dynamically unstable at large $Ra$ and cannot be found by standard time integration, so we employed a purpose-written code that iteratively solves the time-independent equations. We computed rolls with $\Gamma=2$ fixed for $Ra\lesssim 2\times 10^{10}$ and with the parameter-dependent aspect ratios $\Gamma^{*}$ and $\Gamma^{*}_{loc}$ (cf. figure 1) that globally and locally maximize $Nu(\Gamma)$, respectively, for $Ra\leq 10^{14}$. These $Ra$ values are evidently large enough to reach asymptotia: the results reported below strongly suggest that fixed-$\Gamma$ rolls asymptotically transport heat like $Nu\sim Ra^{1/4}$ while the ever-narrowing rolls of aspect ratio $\Gamma^{*}$ achieve the classical $Nu\sim Ra^{1/3}$ scaling.

Following Rayleigh (1916), we model RBC using the Boussinesq approximation to the Navier–Stokes equations with constant kinematic viscosity $\nu$, thermal diffusivity $\kappa$, and coefficient of thermal expansion $\alpha$. We nondimensionalize lengths by the layer height $h$, temperatures by the fixed difference $\Delta$ between the boundaries, velocities by the free-fall scale $U_{f}=\sqrt{g\alpha h\Delta}$, and time by the free-fall time $h/U_{f}$. Calling the horizontal coordinate $x$ and the vertical coordinate $z$, the gravitational acceleration of magnitude $g$ is in the $-\mathbf{\hat{z}}$ direction. The evolution equations governing the dimensionless velocity vector $\mathbf{u}=(u,w)$, temperature $T$, and pressure $p$ are then

where

The dimensionless spatial domain is $(x,z)\in[0,\Gamma]\times[0,1]$, and all variables are horizontally periodic. The top and bottom boundaries are isothermal with $T=0$ and $T=1$, respectively, while no-slip conditions require $\mathbf{u}$ to vanish on both boundaries. The conductive state $(\mathbf{u},T)=(\mathbf{0},1-z)$ becomes unstable when $Ra$ increases past the critical value $Ra_{c}\approx 1708$ (Jeffreys, 1928), at which a roll pair with horizontal period $\Gamma\approx 2.016$ bifurcates supercritically. As $Ra\to\infty$ the horizontal period of the narrowest marginally stable roll pair decreases as $\mathcal{O}(Ra^{-1/4})$, while the horizontal period of the fastest-growing linearly unstable mode decreases more slowly as $\mathcal{O}(Ra^{-1/8})$.

In terms of the dimensionless solutions to 1, the Nusselt number is

where $\langle\cdot\rangle$ denotes an average over the spatial domain and infinite time. For steady states no time average is needed.

To compute rolls at $Ra$ values large enough to reach the asymptotic regime we developed a numerical scheme by adapting the approach of Wen et al. (2020) and Wen & Chini (2018) to the case of no-slip boundary conditions. In these numerics the temperature is represented using the deviation $\theta$ from the conductive profile, meaning $T=1-z+\theta$, and the velocity is represented using a stream function $\psi$, where $\mathbf{u}=\partial_{z}\psi\mathbf{\hat{x}}-\partial_{x}\psi\mathbf{\hat{z}}$ so that the (negative) scalar vorticity is $\omega=\partial_{x}w-\partial_{z}u=-\nabla^{2}\psi$. In terms of these variables, steady ($\partial_{t}=0$) solutions of 1 satisfy

with fixed-temperature and no-slip boundary conditions,

To compute solutions of the time-independent equations 4 and 5 by an iterative method, we do not need to impose all boundary conditions precisely on each iteration—the conditions need to hold only for the converged solution. Thus we do not impose (5$c$) exactly, instead using approximate boundary conditions on $\omega$ for equation 4a. These are derived by Taylor expanding $\psi$ about the top and bottom boundaries to find

where $\delta>0$ is small. Combining equations (4b) with (5$b$,$c$) and neglecting $\mathcal{O}(\delta^{4})$ terms in 6 give the approximate boundary conditions,

In computations we set $\delta$ to be the distance between the boundary and the first interior mesh point.

The time-independent equations (4) are solved numerically subject to boundary conditions (5$a$,$b$) and (7) using a Newton–GMRES (generalized minimal residual) iterative scheme. The spatial discretization is spectral, using a Fourier series in $x$ and a Chebyshev collocation method in $z$ (Trefethen, 2000). All of our computations had at least 20 collocation points in the viscous and thermal boundary layers. At $Ra$ just above the linear instability, iterations starting from the unstable eigenmode converge to the steady rolls we seek. At larger $Ra$, already-computed steady rolls from nearby $Ra$ and $\Gamma$ values were used as the initial iterate. Every 2 to 4 Newton iterations, we change the boundary values of the iterate to match the $\partial_{z}\psi=0$ boundary condition exactly. Prior to convergence this makes the boundary values slightly inconsistent with the governing equations, but the converged solutions satisfy the equations and the no-slip boundary conditions to high precision. Newton iterations were carried out until the Lebesgue $L^{2}$-norm of the residual of the governing steady equations had a relative magnitude less than $10^{-10}$. To accurately locate $\Gamma^{*}$ and $\Gamma^{*}_{loc}$, rolls were computed at several nearby $\Gamma$, and then $Nu(\Gamma)$ was interpolated with cubic splines like those in figure 1($c$). Details of computational results, including resolutions used, are included in the supplementary material.

We computed steady $Pr=1$ rolls for aspect ratios $\Gamma$ encompassing the three distinguished values indicated by figure 1($c$): the fixed value $\Gamma=2$ and the $Ra$-dependent values $\Gamma^{*}$ and $\Gamma^{*}_{loc}$ that globally and locally maximize $Nu$ over $\Gamma$. As previously observed by Sondak et al. (2015), the $Nu(\Gamma)$ curve has a single maximum when $Ra$ is small and develops a second local maximum at smaller $\Gamma$ when $Ra$ increases past roughly $2\times 10^{5}$. The value of $Nu$ at this second local maximum remains less than the value at the first, so the picture remains as in figure 1($c$) with $\Gamma^{*}$ on the left and $\Gamma^{*}_{loc}$ on the right, in contrast to the $Pr=10$ and 100 cases (Sondak et al., 2015). For most $Ra$ values we did not compute rolls over a full sweep through $\Gamma$ as in figure 1($c$), instead searching over $\Gamma$ only as needed to locate $\Gamma^{*}$ and $\Gamma^{*}_{loc}$. The rest of this section reports Nusselt number and Reynolds number scalings for the computed steady rolls, and the supplementary material provides tabulated data.

To compare heat transport by steady rolls with that by turbulent thermal convection, we compiled Nusselt number data from high-$Ra$ DNS with $Pr=1$ or 0.7 and laboratory experiments where the estimated $Pr$ is between 0.7 and 1.3. Figure 5 shows these $Nu$ values compensated by $Ra^{1/3}$, along with $Nu$ values of steady convection rolls at the $Nu$-maximizing aspect ratios $\Gamma^{*}$. Strikingly, heat transport by the $Nu$-maximizing 2D steady rolls is larger than transport by turbulent convection in all cases.

The turbulent data shown in figure 5, as detailed in the figure caption, include DNS in horizontally periodic 2D and 3D domains, wherein 2D steady rolls solve the equations of motion, as well as 3D DNS and laboratory experiments in cylinders that do not admit 2D rolls. Values of $Nu$ for steady rolls with $\Gamma=2$ fixed are omitted from figure 5 for clarity, but they lie below all turbulent values once $Ra$ approximately exceeds $2\times 10^{9}$ (cf. figure 2), and this gap would only widen at larger $Ra$ if their $Nu\sim Ra^{1/4}$ scaling persists. The laboratory data sets in figure 5 have unavoidably varying $Pr$ values that can be hard to estimate, as well as non-Oberbeck–Boussinesq effects (Urban et al., 2011, 2012, 2014). The figure includes only a narrow range of estimated $Pr$ values in order to avoid significant non-Oberbeck–Boussinesq effects. When data over a wider range of estimated $Pr$ is included, a few data points from the experiments of Chavanne et al. (2001) lie above the $Nu(\Gamma^{*})$ values of steady rolls, as shown in the supplementary material.

Our finding that steady rolls of $Nu$-maximizing aspect ratios apparently display classical $Nu\sim Ra^{1/3}$ asymptotic scaling does not ineluctably imply anything about turbulent convection. Taking a dynamical systems point of view, however, steady solutions admitted by the domain are fixed points of 1, so they and their unstable manifolds are part of the global attractor. Turbulent trajectories may linger near these fixed points and so inherit some quantitative features (Kooloth et al., 2021), as has been found for unstable coherent states in shear flows (Nagata, 1990; Waleffe, 1998; Wedin & Kerswell, 2004; Gibson et al., 2008; Suri et al., 2020; Graham & Floryan, 2021). Indeed, figure 5 shows scaling similarities between steady and turbulent convection. Further exploration of the global attractor calls for study of 3D steady flows. Recently computed ‘multi-scale’ 3D steady states (Motoki et al., 2021) give larger $Nu$ values than all 2D rolls at moderate $Ra$, but their scaling at large $Ra$ is unknown. Simpler 3D steady convection patterns remain to be computed as well. Analytically, it is an open challenge to construct approximations of 2D or 3D steady flows that are asymptotically accurate as $Ra\to\infty$, as has been done for 2D rolls between stress-free boundaries (Chini & Cox, 2009; Wen et al., 2020). Such constructions could be used to verify that $Nu\sim Ra^{1/3}$ is indeed the exact asymptotic scaling for the $Nu$-maximizing rolls we have computed, as well as to determine the precise $Re$–$Ra$ scaling relations for rolls of both $Nu$- and $Re$-maximizing aspect ratios.

More generally, figure 5 highlights the absence of reproducible evidence for ultimate $Nu\sim Ra^{1/2}$ scaling, and it raises the intriguing possibility that steady rolls with $Nu\sim Ra^{1/3}$ might transport more heat than turbulent convection as $Ra\to\infty$. We know of no counterexamples to this hypothesis, including in the case of stress-free boundaries (Wen et al., 2020). Heat transport by solutions of 1 with no-slip isothermal boundaries has been mathematically proved to be limited by $Nu\leq\mathcal{O}(Ra^{1/2})$ (Howard, 1963; Doering & Constantin, 1996), but it remains unknown whether any solutions attain the ultimate scaling of this upper bound. One avenue for pursuing a stronger mathematical statement is to study two conjectures suggested by our computations: that steady convection maximizes $Nu$ among all solutions of 1 regardless of their stability or time-dependence, and that steady solutions of 1 are subject to an upper bound of the form $Nu\leq\mathcal{O}(Ra^{1/3})$. Therefore, although numerically computed flows can never determine $Ra\to\infty$ scaling definitively, our results suggest a new mathematical approach that may be able to finally resolve the question of asymptotic $Nu$ scaling in turbulent convection.

After this manuscript was written our senior author, Charles Doering, passed away too soon. Beyond his many contributions to the present study, we are forever indebted to him for his mentorship, to say nothing of his many lasting contributions to the field of fluid dynamics. He will be deeply missed by us and many others. We also want to acknowledge helpful discussions about the present work with L.M. Smith, D. Sondak, and F. Waleffe. This work was supported by US National Science Foundation awards (DMS-1515161, DMS-1813003), Canadian NSERC Discovery Grants Program awards (RGPIN-2018-04263, RGPAS-2018-522657, DGECR-2018-00371), and computational resources provided by Advanced Research Computing at the University of Michigan.

The authors report no conflict of interest.