Scaling in Rayleigh–Bénard convection

Abstract

We consider the Nusselt–Rayleigh number problem of Rayleigh–Bénard convection and make the hypothesis that the velocity and thermal boundary layer widths, $\delta _u$ and $\delta _T$, in the absence of a strong mean flow are controlled by the dissipation scales of the turbulence outside the boundary layers and, therefore, are of the order of the Kolmogorov and Batchelor scales, respectively. Under this assumption, we derive $Nu \sim Ra^{1/3}$ in the high $Ra$ limit, independent of the Prandtl number, $\delta _T/L \sim Ra^{-1/3}$ and $\delta _u/L \sim Ra^{-1/3} Pr^{1/2}$, where $L$ is the height of the convection cell. The scaling relations are valid as long as the Prandtl number is not too far from unity. For $Pr \sim 1$, we make a more general ansatz, $\delta _u \sim \nu ^{\alpha }$, where $\nu$ is the kinematic viscosity and assume that the dissipation scales as $\sim u^3/L$, where $u$ is a characteristic turbulent velocity. Under these assumptions we show that $Nu \sim Ra^{\alpha /(3-\alpha )}$, implying that $Nu \sim Ra^{1/5}$ if $\delta _u$ were scaling as in a Blasius boundary layer and $Nu \sim Ra^{1/2}$ (with some logarithmic correction) if it were scaling as in a standard turbulent shear boundary layer. It is argued that the boundary layers will retain the intermediate scaling $\alpha = 3/4$ in the limit of high $Ra$.

1. Introduction

An ideal Rayleigh–Bénard convection experiment would be to let a fluid be enclosed by two infinite horizontal walls separated by a vertical distance $L$, keeping both the walls at constant temperature, with the lower wall at a higher temperature, and measure the heat flux through the walls in a stationary state. Given certain constraints the flow in such a configuration can be described by the Navier–Stokes equations under the Boussinesq approximation (Tritton 1988)

(1.1)$$\begin{gather} \frac{{\mbox{D}} \boldsymbol{u}}{{\mbox{D}} t} ={-} \frac{1}{\rho} \boldsymbol{\nabla} p + g \theta\boldsymbol{e}_z + \nu \nabla^2 \boldsymbol{u}, \end{gather}$$

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

(1.3)$$\begin{gather}\frac{{\mbox{D}} \theta }{{\mbox{D}} t} = \kappa \nabla^2 \theta, \end{gather}$$

where $\rho, p, g, \nu$ and $\kappa$ are density, pressure, acceleration due to gravity, kinematic viscosity and diffusivity, ${\boldsymbol {e} }_z$ is the vertical unit vector and $\theta = \beta T$ is non-dimensional temperature, with $\beta$ being the thermal expansion coefficient. We let $\theta _0$ and $\theta _1$ be the temperatures at the lower and upper walls and $\Delta \theta = \theta _0 - \theta _1$. If $\Delta \theta \ll \theta _0$, the input of the experiment can be specified by the Rayleigh and Prandtl numbers

(1.4a,b)\begin{equation} Ra = \frac{g \Delta \theta L^3}{\nu \kappa},\quad Pr = \frac{\nu}{\kappa}, \end{equation}

while the output can be specified by the Nusselt number

(1.5)\begin{equation} Nu ={-}{\frac{{\mbox{d}} \langle \theta \rangle }{{\mbox{d}} z} |_{z=0} } / {( \Delta \theta/L )}, \end{equation}

where $\langle \rangle$ is the horizontal mean. The Nusselt number is the ratio between the actual heat flux and the heat flux in the case the fluid is static. The Nusselt–Rayleigh–Prandtl number problem, in the following just called the Nusselt–Rayleigh number problem, is to predict the output given the input, that is, find $F$ such that

(1.6)\begin{equation} Nu = F(Ra, Pr). \end{equation}

Within the format of a brief introduction it is impossible to give a fair account of all theoretical and experimental efforts that have been devoted to the problem. Therefore, we will only consider some of the previous work that are of special relevance to the present paper. For more extensive reviews, the reader is referred to Siggia (1994), Ahlers, Grossmann & Lohse (2009) and Chillà & Schumacher (2012).

Building on a mixing length theory by Priestley (1959), Kraichnan (1962) developed an analysis predicting two types of scaling laws with respect to $Ra$, each divided into two scaling laws with respect to $Pr$. Putting $Nu \sim Ra^{\gamma } Pr^{\tau }$, he predicted $\gamma = 1/3$ for moderately large $Ra$, with $\tau = 0$ for $Pr > 0.1$ and $\tau = 1/3$ for $Pr < 0.1$, and $\gamma = 1/2$ for very high $Ra$, with $\tau = 1/2$ for $Pr < 0.15$ and $\tau = -1/4$ for $Pr > 0.15$. The expression for the large $Ra$ regime also included a logarithmic correction $(\ln (Ra))^{-3/2}$. The reason why there is a transition from $\gamma = 1/3$ to $\gamma = 1/2$ (with a logarithmic correction) is that the boundary layers are supposed to undergo a transition at some high $Ra$. The width of the viscous layer of a turbulent boundary layer scales almost linearly with the Reynolds number with a deviation that can be expressed with a logarithmic factor. This is the reason why there is a corresponding logarithmic factor in the high $Ra$ expression. An interesting essay on Kraichnan's hypothesis, including some biographical information on its genesis, is written by Doering (2020a).

Castaing et al. (1989) developed a model where the convection cell is divided into three zones: the boundary layers, a well mixed middle zone and a mixing zone separating the boundary layers from the middle zone. Using this model they managed to deduce $\gamma = 2/7$, which is close to the experimental value extracted from several experiments. Shraiman & Siggia (1990) deduced the same value based on a very different model in which the central assumption is that the boundary layers are fully turbulent. Grossmann & Lohse (2000) argued that there are eight different regimes of $Nu \sim Ra^{\gamma } Pr^{\tau }$ with different combinations of $(\gamma, \tau )$ depending on three different conditions: (1) whether kinetic energy dissipation is mainly taking place in the boundary layers or in ‘the bulk’, (2) whether thermal dissipation is mainly taking place in the boundary layers or in ‘the bulk’, and (3) whether $\delta _T < \delta _u$ or $\delta _T > \delta _u$, where $\delta _T$ and $\delta _u$ are the thermal and velocity boundary layer widths, respectively.

A central issue in all proposed solutions to the Rayleigh–Nusselt number problem is the nature of the boundary layers. Castaing et al. (1989) assume that the boundary layers resemble a laminar plane Couette flow and that they are stabilised by a ‘shearing wind’. Shraiman & Siggia (1990) make the assumption that there exists ‘a persistent mean flow, which sufficiently close to the horizontal plates can be approximated by a linear profile’ and that the boundary layers associated with the mean flow are turbulent. Grossmann & Lohse (2000) assume that the boundary layers are associated with ‘the mean large-scale velocity near the boundaries of the cell’ and that they start as laminar at low $Ra$ and undergo a transition to become turbulent at high $Ra$. The analogy with the fully developed turbulent shear boundary layer in the limit of high $Ra$ was introduced by Kraichnan (1962). In the following section we will take a closer look at this analogy.

An interesting case for $\gamma = 1/3$ is derived from the marginal stability argument of Howard (1966) that sometimes is attributed to Malkus (1954). The stability of the boundary layers can be assumed to be controlled by the Rayleigh number based on the thermal boundary layer width, $\delta _T$,

(1.7)\begin{equation} Ra_{\delta_T} = Ra \left ( \frac{\delta_T} {L} \right )^3 \sim Ra Nu^{{-}3}. \end{equation}

If it is assumed that the boundary layers will remain stable, $Ra_{\delta _T}$ must stay below a critical threshold and we must have $\gamma \ge 1/3$. On the other hand, as long as they are stable they will grow by diffusion and, therefore, $\delta _T$ will adjust so that $Ra_{\delta _T}$ is just below the threshold, implying that $\gamma = 1/3$. In the same tradition, extensive efforts have been made to rigorously establish upper bounds on the heat transfer by maximising $Nu$ over all possible flow fields. The most impressive result is the upper bound $Nu \le 0.152 Ra^{1/3}$, derived in the infinite $Pr$ limit by Chan (1971). The techniques used in such analyses have been reviewed by Howard (1972) and Busse (1978).

A vast number of experiments have been carried out on the Nusselt–Rayleigh number problem. In most experiments $\gamma$ fall in the range $[2/7, 1/3]$ (e.g. Goldstein & Tokuda 1980; Castaing et al. 1989; Wu & Libchaber 1991; Chavanne et al. 1997; Lui & Xia 1998; Ashkenazi & Steinberg 1999; Glazier et al. 1999). The experiment that has reached the highest $Ra$ was carried out by Niemela et al. (2000) who used helium that was cooled close to the lambda point. A log-log plot extending over more than ten decades shows an almost perfect straight line with slope $\gamma = 0.31$. When the data are plotted as $Nu/Ra^{1/3}$ in a lin-log plot, there is a discernible deviation from a straight line. In the range $Ra = 10^{6}\unicode{x2013}10^{14}$ there is a convex approach to a straight line. In the range $Ra = 10^{14}\unicode{x2013}10^{17}$ there is a concave dip below a straight line. The Prandtl number in the experiment was constant, $Pr \approx 0.7$, up to $Ra \approx 10^{13}$, above which it was continuously increasing to $Pr \approx 30$ at $Ra = 10^{17}$. Niemela et al. (2000) concluded that the Prandtl number effect on the Nusselt number is small. Subsequently, Niemela & Sreenivasan (2006) found a way to introduce a slight correction of the data, which changed the scaling exponent to $\gamma = 0.32$. The plot given by Niemela & Sreenivasan (2006) with the corrected exponent is reproduced in figure 1. As pointed out by Niemela et al. (2000), there is no sign of a transition to the asymptotic high $Ra$ range predicted by Kraichnan (1962). In the authors eyes this plot is absolutely fascinating, not because the corrected exponent is $0.32$ rather than $0.31$ or exactly $1/3$, but because the scaling extends over ten decades of $Ra$. Such an extended scaling range suggests that the underlying physics is ruled by some permanent principle.

Figure 1. Nusselt number versus Rayleigh number plotted using slightly corrected data from the experiment by Niemela et al. (2000). Reproduced from Niemela & Sreenivasan (2006).

In the wake of the experiment by Niemela et al. (2000) the $\gamma = 1/2$ prediction of Kraichnan (1962) has been the object of a great number of investigations as well as a lively debate. Chavanne et al. (2001) claimed that they found experimental evidence of a regime with $\gamma \approx 0.38$ at $Ra > 10^{11}$ and interpreted this result as evidence of a transition to the ‘ultimate’ Kraichnan regime, despite the difference between $0.38$ and $0.5$. Funfschilling, Bodenschatz & Ahlers (2009) carried out experiments up to $Ra \approx 3\times 10^{14}$ and reported results consistent with Niemela et al. (2000), with no evidence of a transition to the ultimate regime. Roche et al. (2010) assessed data from different experiments and claimed that there is clear evidence of a transition to the ultimate regime in the range $Ra = 10^{11}\unicode{x2013}10^{15}$, displaying $\gamma = 0.33\unicode{x2013}0.43$. He et al. (2012) also claimed that they found experimental evidence of a transition at $Ra \approx 10^{13}\unicode{x2013}10^{14}$, where they measured $\gamma \approx 0.38$, which they interpreted as an approach to the ultimate regime. These claims were questioned by Urban, Musilovà & Skrbek (2011) and Urban et al. (2012) who found that $\gamma = 1/3$ is retained if the data are accurately analysed. This was in turn questioned by He et al. (2013) who suggested that the results of Urban et al. (2012) were artifacts of particular parameter choices that led to strong non-Boussinesq effects. Skrbek & Urban (2015) made a systematic investigation of virtually all experimental studies until 2015 and found that the evidence for a transition is weak. The claim of He et al. (2012) was also criticised by Doering (2020b) who scrutinised their curve fitting procedure and found that it was far from satisfactory. In an attempt to reconcile different data sets in which some of them seem to show evidence of a transition Roche (2020) developed a model of a transition from a semi-turbulent boundary layer to a fully turbulent boundary layer in accordance with the prediction of Kraichnan (1962). Iyer et al. (2020) investigated the issue by carrying out direct numerical simulations (DNS) in a cylindrical cell with aspect ratio $1/10$ and found no evidence of a transition. On the contrary, the $\gamma = 1/3$ scaling was observed up to $Ra = 10^{15}$.

In a recent meta study by Ahlers et al. (2022), data from a great number of experiments were compiled and a correction with respect to the aspect ratio of the convection cell was proposed, which improved the collapse of the data. The data were divided into two sets, where $Pr \approx 4.4$ in the first and $Pr \approx 0.8$ in the second. Two lin-log plots show $Nu/Ra^{1/3}$ versus $Ra$ uncorrected and two similar plots show the corrected data. The $Pr \approx 0.8$ plots are reproduced in figure 2. All four curves have a maximum around $Ra = 10^{6}$ and a convex decrease at higher $Ra$, seemingly towards a straight line, which would correspond to a perfect $\gamma = 1/3$ range. For $Pr \approx 4.4$, with $Ra_{max} \approx 10^{12}$, the curves have still a clear negative slope in the high $Ra$ end. For $Pr \approx 0.8$, on the other hand, the curves seem to level off to a straight line in the high $Ra$ end, $Ra \sim 10^{12}\unicode{x2013}10^{14}$, although there is no range of extension where the curve is a perfect straight line. There is no discernible difference between the two sets at different Prandtl number, which is also consistent with recent analyses of DNS data by Pandey & Sreenivasan (2021). On the contrary, it appears that there would be an almost perfect collapse of the two sets if they had been presented in the same plot. In the data from three of the low Prandtl number experiments the curves show upwards bending short tails in the high $Ra$ end with $\gamma \approx 0.38$. In the uncorrected data the three curves start to bend upwards at different $Ra$. Only after the data have been corrected the curves show a common point of transition. Nevertheless, Ahlers et al. (2022) interpret these tails as evidence of a transition to the ultimate Kraichnan regime taking place at $Ra \approx 10^{14}$, in spite of the criticism of Skrbek & Urban (2015) and Doering (2020b) as well as the numerical evidence presented by Iyer et al. (2020).

Figure 2. (a) Lin-log plot of the compensated Nusselt number versus Rayleigh number from different data sets with $Pr \approx 0.8$. (b) The same data plotted using an aspect ratio correcting model. The inset shows an enlargement at the highest $Ra$, illustrating a supposed transition to the ultimate regime. Reproduced from Ahlers et al. (2022).

It is now time to put the debate aside and develop an independent analysis. In the following we will first take a closer look at the hypothesis of Kraichnan (1962) and formulate an alternative hypothesis. Then we will analyse the Nusselt–Rayleigh number problem from a slightly different point of view as compared with previous studies and suggest that the wide scaling range seen in figure 1 is a manifestation of the fact that the boundary layers are of a special type in Rayleigh–Bénard convection. In the end, we will also argue that the $\gamma \approx 1/3$ scaling seen in figure 1 is the ultimate regime.

2. An alternative hypothesis

Kraichnan (1962) makes the hypothesis that ‘big eddies develop boundary layers in strict analogy with the boundary layer of a fully developed shear flow’. In a high-Reynolds-number turbulent flow a standard turbulent shear boundary layer will therefore form in the vicinity of a wall even in the absence of a mean velocity. The role played by the mean velocity will be taken over by the turbulent velocity, $u$. Accordingly, there is a viscous layer whose width scales as $\delta _u \sim \nu /u_{\tau }$, where $u_{\tau } = \sqrt {\tau _w/\rho }$ is the friction velocity based on the wall shear stress $\tau _w$, associated with the turbulence. We use the term ‘viscous layer’ where Kraichnan used ‘viscous sublayer’ to denote the region close to the wall where the dynamics is strongly influenced by viscous forces. Instead of $\delta _ u = 20 \nu /u_\tau$ as used by Kraichnan we set $\delta _u = 50 \nu /u_{\tau }$ for a standard shear boundary layer, so that the viscous layer consists of the whole region below the log layer, including both the two layers that in modern text books (see e.g. Pope 2000) are called the ‘viscous sublayer’ and the ‘buffer layer’. Outside the viscous layer Kraichnan (1962) assumes that there is a log layer where

(2.1)\begin{equation} u = u_{\tau} \left ( \frac{1}{\bar{\kappa}} \ln \left (\frac{u_{\tau} z}{\nu} \right ) + B \right), \end{equation}

and $\bar {\kappa }$ is the Kármán constant. All readers who are familiar with the basics of turbulent boundary layers will appreciate that the Kraichnan (1962) analogy is indeed very ‘strict’. If we let $L$ be the outer length scale of the boundary layer and $u$ the corresponding velocity scale of the largest eddies, we will have

(2.2)\begin{equation} u \sim u_{\tau} \ln(Re_\tau) \Rightarrow Re \sim Re_{\tau} \ln(Re_{\tau}), \end{equation}

where $Re_{\tau } = u_{\tau } L /\nu$ is the friction Reynolds number and $Re = u L/\nu$. The width of the viscous layer thus scales as

(2.3)\begin{equation} \delta_u \sim Re_{\tau}^{{-}1} L \sim \ln(Re_{\tau}) Re ^{{-}1} L. \end{equation}

This strong scaling with $Re$ means that the viscous layer is extremely thin in comparison to virtually all other length scales in fluid mechanics. The width of a laminar boundary layer scales as $Re^{-1/2} L$ and the Kolmogorov scale, $\eta = (\nu ^3/\epsilon )^{1/4}$, as $Re^{-3/4} L$, where $L$ in each case is the relevant macroscopic length scale. The author can only think of one example of a stronger scaling with $Re$ and that is a shock wave whose width scales exactly as $Re^{-1} L$ (Whitham 1970).

Since $\eta$ generally is considered to be the smallest length scale that can develop in an incompressible flow, it seems paradoxical that the viscous layer with respect to scaling with $Re$ more resembles a shock wave than $\eta$. The paradox is resolved if we take into consideration that there is a continuous decrease of $\eta$ through the log layer down to the wall. The reason for this is that there is local balance between turbulence production and dissipation in the log layer, expressed as (see e.g. Pope 2000)

(2.4)\begin{equation} \epsilon ={-} \frac{\partial U}{\partial z} \langle uw \rangle = \frac{u_\tau^3}{\bar{\kappa} z}, \end{equation}

where $U$ is the mean velocity and $\langle uw \rangle$ is the Reynolds stress. As the wall is approached, the Kolmogorov scale decreases as $\sim z^{1/4}$, due to increased turbulence production by mean shear close to the wall. The dependence of $\eta$ on wall distance can be expressed as

(2.5)\begin{equation} \eta = \bar{\kappa}^{1/4} Re_z^{{-}3/4} z, \end{equation}

where $Re_z = u_\tau z/\nu$. The width of the viscous layer, corresponding to $Re_z = 50$, can thus be estimated as

(2.6)\begin{equation} \delta_u \approx 20 \eta |_{z=\delta_u}, \end{equation}

stating that it ends where $z$ is equal to a fixed number of local Kolmogorov scales. The continuous decrease of dissipation with wall distance, in accordance with (2.4), is thus crucial for the existence of the thin viscous wall layer in a shear flow. Without such a decrease, a viscous layer obeying the scaling (2.3) would be thinner than the local Kolmogorov scale, which can be deemed impossible.

In the following it will be argued that dissipation does not necessarily scale as (2.4) above the viscous wall layer in a turbulent flow and that Kraichnan's hypothesis therefore cannot be generally valid. The argument is based on a comparison between a turbulent flow confined by walls and a corresponding flow where the walls are absent. First, we consider a flow in a cubic tank with side $L$ where turbulent kinetic energy with characteristic velocity $u$ and integral length scale ${\sim }L$ is produced by some mechanical device in the central region. Second, we consider a flow that is generated by letting the tank be surrounded by identical tanks with identical flows that in turn are surrounded by identical tanks ad infinitum, after which all walls are removed. In this way, an almost homogeneous and isotropic turbulent field will be generated. In particular, the dissipation will be almost homogeneous. Since both flows are in a statistically stationary state, the total dissipation will be the same since energy injection is the same. It is also reasonable to assume that the characteristic velocity as well as the dissipation in the central region of a cube will be almost equal in the two flows. The dissipation in the central region can be estimated using the fundamental scaling law of three-dimensional turbulence, stating that

(2.7)\begin{equation} \frac{\epsilon}{u^3/L} = {\mbox{Constant}} \quad {\mbox{as}}\quad Re = \frac{uL}{\nu} \rightarrow \infty. \end{equation}

We thus obtain the standard estimate $\eta \sim Re^{-3/4} L$ for the Kolmogorov scale. According to the hypothesis of Kraichnan (1962), viscous layers would form close to the walls in the first flow and the outer length scale of the boundary layers would scale with $L$. The ratio between the width of the viscous layer and the Kolmogorov scale in the central region would be

(2.8)\begin{equation} \delta_u/\eta \sim Re^{{-}1/4} \ln Re_\tau. \end{equation}

In the high-Reynolds-number limit we would thus have $\delta _u \ll \eta$. Outside the viscous layers there would be log layers in which dissipation is continuously decreasing away from the walls in accordance with (2.4). The width of the log layers would scale with the turbulent integral length scale and, thus, be of the order of $L$, as also pointed out by Grossmann & Lohse (2011). This is analogous to channel flow where the log layer ends at a fixed ratio of the channel hight. Since dissipation would decrease as $z^{-1}$ through the bulk of the flow, total dissipation would be much larger in the first flow as compared with the second. However, this is contradicting the fact that total dissipation must be equal. It can be concluded that the mere presence of walls in a turbulent field is not sufficient for the formation of standard shear boundary layers. In order for the effect of the walls on total dissipation in the tank to be negligible they must have a much weaker influence on the flow field as a whole.

An alternative hypothesis can be formulated by assuming that dissipation is uniform above the viscous layers, in which case $\delta _u$ must obey a different scaling than (2.3). In analogy with (2.6) it is reasonable to assume that $\delta _u$ also in this case ends at a fixed number of Kolmogorov scales, with the important difference that $\eta$ now is independent of wall distance above $\delta _u$. The width of the viscous layer, in the following called the boundary layer, will then scale as

(2.9)\begin{equation} \delta_u \sim Re^{{-}3/4} L, \end{equation}

which is the important difference compared with hypothesis (2.3) of Kraichnan (1962). In the boundary layers turbulent kinetic energy is strongly attenuated by increased dissipation. Outside the boundary layers, the main influence of the walls is to induce an increasing degree of anisotropy with decreasing wall distance. As the wall is approached, the wall normal velocity component is damped and the tangential component is enforced by pressure-strain, while total turbulent kinetic energy remains approximately independent of wall distance, just as in the log layer of a classical shear boundary layer (see Pope 2000). This means that the logarithmic scaling (2.1) of the turbulent velocity proposed by Kraichnan (1962) cannot hold.

In a convection cell the production of turbulent kinetic energy by buoyancy is approximately independent of wall distance above the boundary layer. It is therefore reasonable to assume that the flow more closely resembles the flow in a tank rather than, for example, a channel or pipe flow and that the scaling (2.9) therefore holds rather than (2.3). It may be objected that if a mean flow emerges in the cell it may sustain the turbulence in the same way as in (2.4). However, even if a mean flow emerges it cannot sustain the turbulence, since unlike a channel or pipe flow, there is no mean pressure work at the boundaries of a convection cell that can sustain the mean flow. For the thermal boundary layer thickness, a hypothesis corresponding to $\delta _u \sim \eta$ would be that it scales with the dissipation scale of thermal fluctuations. Batchelor (1959) argues convincingly that the dissipation scale of thermal fluctuations is $\eta _B = (\nu \kappa ^{2}/\epsilon )^{1/4}$. The Batchelor scale is formed under the assumption that the dissipation scale of thermal fluctuations can be estimated as $(\kappa /\omega )^{1/2}$, where $\omega$ is a characteristic value of the vorticity that in turn can be estimated as $\omega \sim (\epsilon /\nu )^{1/2}$. Batchelor (1959) only develops the argument in the limit of high $Pr$. In the limit of low $Pr$ he finds that the thermal dissipation scale is $\kappa ^{3/4}/\epsilon ^{1/4}$. However, we find his argument for $\eta _B$ convincing not only in the high $Pr$ limit but as long as the low $Pr$ limit is not approached. It is more justified to assume that the thermal boundary layer width is determined by the Batchelor scale than by the Kolmogorov scale. In the case $Pr > 1$, the thermal field can remain fully turbulent within the outer part of the velocity boundary layer, in spite of the fact that the velocity field is not fully turbulent in the same region. Conversely, in the case $Pr < 1$, the velocity field can remain fully turbulent within the outer part of the thermal boundary layer. Our hypothesis can thus be formulated in the following way. In the high-Reynolds-number limit, the Kolmogorov and Batchelor scales of the turbulence outside the boundary layers are the smallest length scales of the flow in a convection cell, determining both the dissipation length scales outside the boundary layers and the widths of the boundary layers.

3. Nusselt–Rayleigh number scaling

The mean temperature equation is

(3.1)\begin{equation} \frac{\partial \langle \theta \rangle} {\partial t } ={-} \frac{\partial}{\partial z} \langle w \theta \rangle + \kappa \frac{\partial^2 \langle \theta \rangle} {\partial z^2}, \end{equation}

where $w$ is the vertical velocity component. In a stationary state the equation can be integrated to give

(3.2)\begin{equation} \langle w \theta \rangle = \kappa \frac{\partial \langle \theta \rangle}{\partial z} - \kappa \left.\frac{\partial \langle \theta \rangle}{\partial z} \right|_{z=0}. \end{equation}

At a position well above the thermal boundary layer, the first term on the right-hand side can be neglected and we have

(3.3)\begin{equation} \langle w \theta \rangle = \Delta \theta \frac{\kappa}{L} Nu. \end{equation}

The equation for the mean kinetic energy per unit mass, $E = \langle \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {u} \rangle /2$, can be written as

(3.4)\begin{equation} \frac{\partial E}{\partial t} ={-} \frac{1}{2} \frac{\partial }{\partial z} \langle w \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} \rangle - \frac{1}{\rho} \frac{\partial}{\partial z} \langle w p \rangle + g\langle w \theta \rangle - \epsilon + \nu \frac{\partial^2 E}{\partial z^2}, \end{equation}

where $\epsilon = \nu \langle \partial _i u_j \partial _i u_j \rangle$ is the mean kinetic energy dissipation per unit mass and the last term is the viscous transport of mean kinetic energy (see e.g. Pope 2000). The energy equation (3.4) can be integrated to give (Shraiman & Siggia 1990)

(3.5)\begin{equation} \epsilon_V = \nu \kappa^2 (Nu-1)Ra/L^4, \end{equation}

where $\epsilon _V$ is the volume-averaged dissipation. Relation (3.5) is an expression of the global dissipation–production balance. Relation (3.3) shows that production is approximately independent of wall distance in the central region of the convection cell. Assuming that both $\epsilon$ and $g\langle w \theta \rangle$ are of leading order above the boundary layers ($z > \delta _u, \delta _T$), we obtain

(3.6)\begin{equation} \epsilon \sim g\langle w \theta \rangle \sim \frac{\nu \kappa^2}{L^4} Nu Ra. \end{equation}

The relation (3.6) is the correspondence to the local dissipation–production balance (2.4) in the log layer of a classical shear boundary layer. There is both experimental and numerical evidence that the flow reaches a state in which dissipation is independent of wall distance and approximately equal to production outside the boundary layers already at $Ra \sim 10^{7}$. Deardorff & Willis (1966) studied the kinetic energy budget by means of hot-wire anemometry in a high aspect ratio convection cell filled with air and found an increasing degree of local dissipation–production balance with increasing $Ra$, with almost perfect balance outside the boundary layers at $Ra = 10^7$ (see their figures 21–23). Kerr (1996, figure 10) plots the dissipation/production ratio from a spectral DNS at $Ra = 2\times 10^{7}$, $Pr = 0.7$ and aspect ratio $\varGamma = 6$. Outside the boundary layers the ratio is almost perfectly constant and just below unity. Pandey et al. (2022, figure 6b) plots the dissipation as a function of wall distance from high aspect ratio DNS at $Ra = 10^7$ and different $Pr$. Again, outside the boundary layers dissipation is almost perfectly constant. Direct numerical simulations at high $Ra$ (e.g. Iyer et al. 2020, figure 3c) show that turbulent kinetic energy is uniform above $\delta _u$, just as in the log layer of a standard shear boundary layer. This result may seem counterintuitive, since the energy associated with the vertical velocity must decrease continuously as the wall is approached, even outside the boundary layers. However, this decrease is compensated by an increase of the energy associated with the horizontal component, in such a way that the total kinetic energy is independent of wall distance (see Kerr 1996, figure 9). If it is assumed that the transport terms in (3.4) associated with pressure and advection can be accurately modelled by a simple gradient diffusion expression, just as in the $k-\epsilon$ model (see e.g. Pope 2000), the dissipation–production balance can be understood from the fact that there is no kinetic energy gradient. Indeed, there seems to be a close analogy between the log layer and the central region of a convection cell, but of a different kind than envisioned by Kraichnan (1962). Turbulent kinetic energy is approximately independent of wall distance, transport terms are subleading in the turbulent kinetic energy equation and there is an approximate local dissipation–production balance.

Traditionally, the relation between the Nusselt number and the thermal boundary layer width is written as

(3.7)\begin{equation} Nu = \frac{1}{2} \frac{L}{\delta_T}. \end{equation}

The factor one half has no special significance. The basic assumption behind (3.7) is that the total temperature drop in the cell takes place over the boundary layers and that the mean temperature drop close to the wall therefore can be written in the form of a ‘law of the wall’,

(3.8)\begin{equation} \theta_0 - \langle \theta(z) \rangle = \tfrac{1}{2}\Delta \theta f(z/\delta_T). \end{equation}

The Nusselt number is then related to $\delta _T$ as

(3.9)\begin{equation} Nu = \frac{1}{2} f^{\prime}|_{z=0} \frac{L}{\delta_T}, \end{equation}

where $f^{\prime } |_{z=0}$, in the absence of an independent definition of $\delta _T$, can be fixed to unity. Combining (3.6) and (3.7) with the hypothesis $\delta _T \sim \eta _B$, we obtain

(3.10)\begin{equation} Nu \sim \frac{L}{\eta_B} = \frac{\epsilon^{1/4} L}{\nu^{1/4}\kappa^{1/2}} \sim Ra^{1/4}Nu^{1/4} \Rightarrow Nu \sim Ra^{1/3}, \end{equation}

independent of Prandtl number. The boundary layer widths scale as

(3.11)$$\begin{gather} \frac{\delta_T}{L} \sim Ra^{{-}1/3}, \end{gather}$$

(3.12)$$\begin{gather}\frac{\delta_u}{L} \sim Ra^{{-}1/3} Pr^{1/2}. \end{gather}$$

The thickness of the thermal boundary layers can be estimated from the prefactor in $Nu = a Ra^{1/3}$, as $\delta _T \approx 0.5 a^{-3/4} \eta _B$. With $a \approx 0.05$ (Iyer et al. 2020; Ahlers et al. 2022), we obtain $\delta _T \approx 5 \eta _B$, which seems reasonable. Sun, Cheung & Xia (2008) report $\delta _T = 0.58\,{\mbox {mm}}$ and $\eta = 0.4\,{\mbox {mm}}$ from an experiment at $Ra = 2.5 \times 10^{10}$ with water ($Pr = 4.3$). These numbers give us $\delta _T \approx 3 \eta _B$, lending direct experimental support to the hypothesis that the thermal boundary layer width is of the order of the Batchelor scale. They also report $\delta _T/L \sim Ra^{-0.32}$, in close agreement with (3.11). The results for $\delta _u$ are less conclusive. Using two different operational definitions of $\delta _u$, they obtain $\delta _u/L \sim Ra^{-0.27}$ and $\delta _u/L \sim Ra^{-0.37}$ in each case, respectively. Also when DNS data are analysed, the scaling of $\delta _u$ displays a large degree of variation. Using five different definitions of $\delta _u$, Scheel & Schumacher (2014) obtain a whole range of different scalings, ranging from $\delta _u/L \sim Ra^{-0.21}$ to $\delta _u/L \sim Ra^{-0.38}$.

In case the Prandtl number is far from unity – exactly how far is difficult to say – the arguments we have developed are probably no longer valid. If $Pr \gg 1$, the thermal boundary layer will be so separated from the turbulence outside the velocity boundary layer, so its width can probably not be determined by the Batchelor scale of the turbulence outside the boundary layers. If $Pr \ll 1$, the thermal boundary layer will be so thick that a substantial fraction of the total dissipation will take place within the thermal boundary layer and $\delta _T$ will probably be decoupled from $\delta _u$. If we, nevertheless, would assume that the analysis can be extended beyond the $Pr \sim 1$ regime, we would need to use $\delta _T \sim \kappa ^{3/4} \epsilon ^{-1/4}$ in the low-Prandtl-number limit, in line with the analysis of Batchelor. We then obtain $Nu \sim Ra^{1/3} Pr^{1/3}$ in this limit, which is consistent with the $Pr < 0.1$ regime of Kraichnan (1962) at moderate $Ra$.

To define a Reynolds number, we need a velocity that we take as the turbulent velocity, $u$. We also assume that the fundamental scaling law (2.7) of three-dimensional turbulence holds. Under this assumption it is straightforward to derive the relation

(3.13)\begin{equation} Re \sim (Nu Ra Pr^{{-}2})^{1/3} = (F(Ra, Pr) Ra Pr^{{-}2})^{1/3}, \end{equation}

which must hold in the high-Reynolds-number limit, irrespectively of what form of $Nu = F(Ra,Pr)$ we assume to be valid. It is worth pointing out that (3.13) also can be formulated as

(3.14)\begin{equation} \frac{u}{u_f} \sim Nu^{1/3} Ra^{{-}1/6} Pr^{{-}1/6}, \end{equation}

where $u_f = \sqrt {gL\Delta \theta }$ is the free fall velocity. According to (3.14), $Nu \sim Ra^{1/2} Pr^{1/2}$ is the only scaling for which the ratio $u/u_f$ looses its dependence on both $Ra$ and $Pr$. In our case we find that

(3.15)$$\begin{gather} Re = \frac{\epsilon^{1/3} L^{4/3}}{\nu} \sim Ra^{4/9} Pr^{{-}2/3}, \end{gather}$$

(3.16)$$\begin{gather}\frac{u}{u_f} \sim Ra^{{-}1/18} Pr^{{-}1/6}. \end{gather}$$

The relation (3.15) has been derived by a number of other investigators (e.g. Kraichnan 1962; Siggia 1994) under different assumptions and seems to be in quite good agreement with experimental data. Ashkenazi & Steinberg (1999) give the experimentally determined dependence $Re \sim Ra^{0.43} Pr^{-0.75}$. In this case, the $Ra$ dependence is in very good agreement with the data while the $Pr$ dependence is deviating a little bit. Grossmann & Lohse (2000) present a plot of data extracted from Chavanne et al. (1997), in which $Re \sim Ra^{1/2} Pr^{-2/3}$. In this case, on the other hand, the $Pr$ dependence is in exact agreement with the data while the $Ra$ dependence is deviating.

It is worth pointing out that only (3.6) and (3.7) combined with $\delta _T \sim \eta _B$ were used to derive (3.10) and that the fundamental scaling law (2.7) was not used. In the following argumentation we will make use of (2.7), which means that Bolgiano scaling is ruled out. The scaling was originally proposed by Bolgiano (1959) for stratified turbulence and is derived under the assumption that (2.7) does not hold. Several authors (e.g. Procaccia & Zeitak 1989; Yakhot 1992) have proposed that Bolgiano scaling should hold for convective turbulence, but the evidence for this is weak (a review is given by Lohse & Xia 2010). The issue has been investigated theoretically and numerically by Verma, Kumar & Pandey (2017) who found that Bolgiano scaling does not hold for convective turbulence. Moreover, the validity of the fundamental scaling relation (2.7) was numerically investigated by Pandey et al. (2022) who found that it is very well satisfied already at $Ra \approx 10^{6}$, with the constant on the right-hand side being approximately equal to $0.2$.

The assumptions we made to arrive at (3.10)–(3.12) can be expressed as $\delta _T \sim Pr^{-1/2} Re^{-3/4}L$ and $\delta _u \sim Re^{-3/4} L$. Let us fix the Prandtl number to unity and make an ansatz that is more general with respect to the scaling with $Re$. We thus set $\delta _T \sim \delta _u \sim Re^{-\alpha } L$, in order to deduce $\gamma$ in $Nu \sim Ra^{\gamma }$ in terms of $\alpha$. Putting $Ra \sim Nu^{1/\gamma }$ and $\kappa = \nu$ in (3.6), we obtain

(3.17)\begin{equation} \epsilon \sim \frac{\nu^3}{L^4} Nu^{(1+ \gamma)/\gamma} \Rightarrow Nu \sim \left ( \frac{\epsilon}{\nu^3} \right )^{\gamma/(1+\gamma)} L^{4\gamma/(1+\gamma)}. \end{equation}

Using $\delta _u \sim \delta _T$ in (3.7) together with (3.17) and the first equality in (3.15), we find that

(3.18)\begin{align} \delta_u \sim Nu^{{-}1} L \sim \left(\frac{\nu^3}{\epsilon}\right)^{\gamma/(1+\gamma)} L^{{-}4\gamma/(1+\gamma)}L = \left(\frac{\nu}{\epsilon^{1/3} L^{4/3}} \right )^{3\gamma/(1+\gamma)} L = Re^{{-}3\gamma/(\gamma +1)} L. \end{align}

We therefore have

(3.19)\begin{equation} \alpha = \frac{3\gamma}{\gamma +1} \Rightarrow \gamma = \frac{\alpha}{3-\alpha}. \end{equation}

For $\alpha = 1/2$, we have $\gamma = 1/5$ and for $\alpha = 1$, $\gamma = 1/2$. The value $\gamma = 2/7$ predicted by Castaing et al. (1989) and Shraiman & Siggia (1990) corresponds to $\alpha = 2/3$. It is interesting to note that $\alpha = 1/2$ is the case when the boundary layer width is equal to the Taylor microscale of the adjacent turbulence and that $\alpha = 2/3$ implies that the boundary layers are thinner than the Taylor microscale but thicker than the Kolmogorov scale. As a historical curiosity, it is worth mentioning that $\alpha = 2/3$ corresponds to the dissipation scale of the direct interaction approximation theory of Kraichnan (1959), which he subsequently abandoned for the Lagrangian-history approximation theory (Kraichnan 1965), which is consistent with the Kolmogorov theory. In the case $\alpha = 1$, we would have

(3.20)\begin{equation} \delta_u \sim ({\eta}/{L})^{1/3} \eta, \end{equation}

which means that the boundary layer width would be much smaller than the Kolmogorov scale.

Grossmann & Lohse (2000, 2011) argue that a mean flow, $U$, develops in a convection cell and that this can play a crucial role in the development of either classical laminar or turbulent boundary layers. In the case of a convection cell with a very large aspect ratio the mean can neither be defined as a spatial mean over a horizontal plane nor as a time average, since such mean values would be zero. The mean must be defined in some other way, presumably identifying it with coherent structures with a much longer life time than the background turbulent eddies. In the authors opinion it seems unlikely that such a mean flow could determine the scaling properties of the boundary layers. The hypothesis can be tested by working out the scaling of the mean/turbulent ratio $U/u$ under the assumption $\delta _u \sim Re_U ^{-\beta } L$,

(3.21)\begin{equation} \frac{U}{u} = Re_U Re^{{-}1} \sim (\delta_u/L)^{{-}1/\beta} Re^{{-}1} \sim Re^{3\gamma/\beta(\gamma +1)-1} \sim Ra^{4(3\gamma - \beta \gamma -\beta)/9\beta(\gamma +1)}, \end{equation}

where we have used (3.18) and (3.15) and assumed that $Pr=1$. If we take $\gamma = 1/3$ (close to some experimental values), we find that

(3.22)\begin{equation} \left.\begin{gathered} \beta = 1/2 \Rightarrow \frac{U}{u} \sim Ra^{2/9},\\ \beta = 1 \Rightarrow \frac{U}{u} \sim Ra^{{-}1/9}. \end{gathered}\right\} \end{equation}

If we take $\gamma = 2/7$ (close to other experimental values), we find that

(3.23)\begin{equation} \left.\begin{gathered} \beta = 1/2 \Rightarrow \frac{U}{u} \sim Ra^{4/27},\\ \beta = 1 \Rightarrow \frac{U}{u} \sim Ra^{{-}4/27}. \end{gathered}\right\} \end{equation}

With laminar boundary layers ($\beta = 1/2$) the mean flow would thus become much stronger than the turbulence in the limit of large $Ra$, which seems implausible. With standard turbulent boundary layers ($\beta = 1$, neglecting the logarithmic correction) it would become much weaker, which leads us into a contradiction. If $U/u \ll 1$, there is no good reason in the first place to assume that $\delta _u$ scales with $Re_U$.

A boundary layer obeying $\alpha = 3/4$ scaling is laminar/turbulent to the same degree as the viscous layer of a standard turbulent boundary layer. We may call it a semi-turbulent boundary layer. As demonstrated in the DNS study by Shi, Emran & Schumacher (2012) and illustrated by visualisations from DNS data by Iyer et al. (2020), it has a rich three-dimensional internal structure and bears little resemblance with a laminar boundary layer. Assuming that the mean flow does not exist, alternatively that it exists and $U/u \sim 1$, we must have $\gamma = 1/5$ for laminar, $\gamma = 1/3$ for semi-turbulent and $\gamma = 1/2$ (with some logarithmic correction) for standard turbulent boundary layers. The laminar regime is consistent with the second regime of Grossmann & Lohse (2000) whose analysis in this case is based on the assumption that the boundary layers scale as if they were laminar.

The horizontal velocity at the top of the boundary layer scales as $u \sim (L\epsilon )^{1/3}$. A characteristic value of the local boundary layer dissipation can thus be estimated as

(3.24)\begin{equation} \epsilon_{bl} \sim \nu \frac{u^2}{\eta^2} \sim \nu^{{-}1/2} L^{2/3} \epsilon^{7/6} = \left( \frac{L}{\eta} \right )^{2/3} \epsilon \sim Ra^{2/9} \epsilon, \end{equation}

implying that $\epsilon _{bl} \gg \epsilon$, which is intuitively reasonable. On the other hand, the ratio between the integrated boundary layer dissipation and the total dissipation can be estimated as $(\eta /L)^{1/3} \sim Ra^{-1/9} \sim Re^{-1/4}$, implying that only a minor fraction of the total dissipation takes place in the boundary layers. Clearly, it is a consistency requirement of our hypothesis that this ratio must go to zero in the limit of high $Re$. The estimate (3.24) is in quite good agreement with DNS data of Kunnen et al. (2008); see figure 3(a) in Lohse & Xia (2010). A similar estimate gives us the scaling of the friction velocity based on the turbulent wall shear stress

(3.25)\begin{equation} u_{\tau} \sim Re^{{-}1/8} (\epsilon L)^{1/3} \sim Re^{{-}1/8} u. \end{equation}

The wall shear stress is weaker in a flow with a semi-turbulent boundary layer than it is according to the shear boundary layer analogy (2.1) of Kraichnan (1962) and stronger than in a laminar boundary layer where $u_{\tau } \sim Re^{-1/4}U$.

The prediction $Nu \sim Ra^{1/3}$ of the present analysis is consistent with the prediction of the marginal stability theory of Malkus (1954) and Howard (1966). Obviously, the two approaches are very different. The marginal stability theory treats the boundary layer as a laminar structure whose width is determined by the length scale characteristic of disturbances that are marginally stable. In the present approach, the boundary layer is treated as a semi-turbulent structure whose width is determined by the length scale at which fully developed turbulence cannot survive. It is unclear to the author if the two approaches can be reconciled. The present approach does not exclude the possibility that the boundary layers are convectively unstable. In the authors view, it seems more likely that they are just above than below the threshold of convective instability and remain there by feedback from the turbulence. In the case $Ra_{\delta _T}$ increases far above the threshold, they become more unstable, resulting in increased seeding of the turbulence in the middle region. Stronger turbulence will lead to a smaller Batchelor scale, pushing back the boundary layers, so that $Ra_{\delta _T}$ is restored. If, on the other hand, $Ra_{\delta _T}$ goes below the threshold the turbulence in the central region is weakened permitting the boundary layers to grow since the Batchelor scale will increase with weakened turbulence. In this way, $Ra_{\delta _T}$ is restored above the threshold. The semi-turbulent boundary layer may thus be marginally unstable rather than marginally stable, which may be helpful in explaining the origin of thermal plumes.

4. Conclusion

A classical boundary layer is identified as the region close to a wall where a mean flow is attenuated. In the present paper it is suggested that a boundary layer in a convection cell should rather be conceived as the region close to a wall where the turbulence is attenuated. The wide scaling range seen in figure 1 can then be seen as a manifestation of the scaling of such a boundary layer, which is intermediate to the scalings of a laminar and a turbulent shear boundary layer. The arguments leading to this conclusion are based on high-Reynolds-number turbulence phenomenology. It is therefore clear that these arguments also imply that the $Nu \sim Ra^{1/3}$ regime is the ultimate regime. To see that the laboratory measurements are indeed carried out in a high-Reynolds-number regime, it is instructive to estimate the Taylor microscale Reynolds number, $Re_{\lambda } = u^2 \sqrt {15/(\nu \epsilon )}$, from the experimental data used in figure 2. Using (3.6) together with $Nu = a Ra^{1/3}$ and $\epsilon = \beta u^3/L$, we obtain $Re_{\lambda } = \beta ^{-2/3} a^{1/6} \sqrt {15} Pr^{-1/3} Ra^{2/9}$. Using $\beta = 0.45$, as in isotropic turbulence (Sreenivasan 1998) or $\beta = 0.2$, as found in DNS of Rayleigh–Bénard convection (Pandey et al. 2022), $a = 0.05$, $Pr = 0.8$ and $Ra > 10^{12}$, we obtain $Re_\lambda > 2000$. Even in the case we would land on the more conservative estimate $Re_{\lambda } > 1000$, we can conclude that the highest Rayleigh number experiments are carried out at some of the highest $Re_{\lambda }$ ever obtained in a laboratory flow, to be compared with $Re_{\lambda } = 852$ measured at the centreline of a jet in the experiment by Anselmet et al. (1984), or $Re_{\lambda } = 1450$ of the turbulent boundary layer in the experiment carried out by Saddoughi & Veeravalli (1994) in the NASA Ames open wind tunnel. As pointed out in the previous section, the width of the thermal boundary layer can be estimated as $\delta _T \approx 5 \eta _B$. A corresponding estimate for the velocity boundary layer would be $\delta _u \approx 5 \eta$. We are thus faced with a wall bounded turbulent flow at $Re_\lambda > 1000$, with a viscous wall layer whose width is about five Kolmogorov scales. It seems unlikely that such a flow could undergo a laminar–turbulent transition. The flow is already as turbulent as it can be and the boundary layer is already as thin as it possibly can be. In the author's interpretation, the upwards bending short tails in the large $Ra$ end of three of the curves of Ahlers et al. (2022) seen in figure 2 can therefore not be significant. They are most likely caused by experimental uncertainties and should therefore be discarded, which is in line with the criticism put forward by Skrbek & Urban (2015) and Doering (2020b). In the plots presented by Ahlers et al. (2022) there is no discernible difference between the data compiled at $Pr \approx 0.8$ and $Pr \approx 4.4$. This is consistent with the prediction that $Nu$ is independent of $Pr$, when $Pr$ is not too far from unity. It is noteworthy that this prediction is a direct consequence of the assumption that the thermal boundary layer width is determined by the Batchelor scale.

The Nusselt–Rayleigh problem does not seem to be a hundred percent clean problem from an experimental point of view. It is likely that different experiments will give slightly different values of $\gamma$ also in the future, even if they are conducted with great care. Direct numerical simulation is a clean complementary method in which the experimental difficulties are circumvented and a high degree of accuracy is attained. The low aspect ratio DNS of Iyer et al. (2020) indicate that there is no transition occurring in the range $Ra = 10^{14}\unicode{x2013} 10^{15}$. In a not too distant future, it will be possible to confirm this result in higher aspect ratio simulations and carry out low aspect ratio simulations at even higher $Ra$. Another interesting research project would be to investigate whether the hypothesis put forward in the present paper is correct to leading order at smaller $Ra$ and to systematically investigate what the sources are of the finite-Rayleigh-number correction to $Nu \sim Ra^{1/3}$. It is the hope of the author that such research will be undertaken in the near future.