Universal scaling law in turbulent Rayleigh–Bénard convection with and without roughness

Abstract

Heat-transfer measurements published in the literature seem to be contradictory, some showing a transition for the dependance of the Nusselt number (${\textit{Nu}}$) with the Rayleigh number (${\textit{Ra}}$) behaviour at ${\textit{Ra}} \approx 10^{11}$, some showing a delayed transition at higher ${\textit{Ra}}$, or no transition at all. The physical origin of this discrepancy remains elusive, but is hypothesised to be a signature of the multiple possible flow configurations for a given set of control parameters, as well as the sub-critical nature of the transition to the ultimate regime (Roche 2020 New J. Phys. vol. 22, 073056; Lohse & Shishkina 2023 Phys. Today vol. 76, no. 11, 26–32). New experimental and numerical heat-flux and velocity measurements, both reaching ${\textit{Ra}}$ up to $10^{12}$, are reported for a wide range of operating conditions, with either smooth boundaries, or mixed smooth–rough boundaries. Experiments are run in water at $40\,^\circ \textrm {C}$ (Prandtl number, ${\textit{Pr}}$ is 4.4), or fluorocarbon at $40\,^\circ \textrm {C}$ (${\textit{Pr}}$ is 12), and aspect ratios 1 or 2. Numerical simulations implement the Boussinesq equations in a closed rectangular cavity with a Prandtl number 4.4, close to the experimental set-up, also with smooth boundaries, or mixed smooth–rough boundaries. In the new measurements in the rough part of the cell, the Nusselt number is compatible with a ${\textit{Ra}}^{1/2}$ scaling (with logarithmic corrections), hinting at a purely inertial regime. In contrast to the ${\textit{Nu}}$ vs ${\textit{Ra}}$ relationship, we evidence that these seemingly different regimes can be reconciled: the heat flux, expressed as the flux Rayleigh number, ${\textit{Ra}}\textit{Nu}$, recovers a universal scaling with Reynolds number, which collapses all data, both our own and those in the literature, once a universal critical Reynolds number is exceeded. This universal collapse can be related to the universal dissipation anomaly, observed in many turbulent flows (Dubrulle 2019 J. Fluid Mech. vol. 867, no. P1, 1).

1. Introduction

Turbulent Rayleigh-Bénard Convection is a model system where a layer of fluid is heated from below and cooled from above. It is controlled by three dimensionless parameters: the Rayleigh number ${\textit{Ra}}$

(1.1) \begin{equation} \textit{Ra} = \frac {g\alpha \Delta \textit{TH}^3}{\nu \kappa }, \end{equation}

the Prandtl number ${\textit{Pr}}$

(1.2) \begin{equation} \textit{Pr} = \frac {\nu }{\kappa }, \end{equation}

and the aspect ratio $\varGamma = L/H$ , where $g$ is the acceleration due to gravity, $\alpha$ is the thermal expansion coefficient, $\Delta T$ is the temperature difference across the fluid layer, $\nu$ is the kinematic viscosity, $\kappa$ is the thermal diffusivity, $H$ is the height of the cell and $L$ is the width of the cell. A long-standing endeavour consists of finding universal scaling laws for the heat transfer, expressed as the dimensionless parameter Nusselt number ${\textit{Nu}}$ (Chillà & Schumacher 2012; Stevens et al. 2013; Lohse & Shishkina 2023)

(1.3) \begin{equation} \textit{Nu} = \frac {\dot {q}\textit{H}}{\lambda \Delta T}, \end{equation}

where $\dot {q}$ is the heat flux, and $\lambda$ is the thermal conductivity. This would allow the extrapolation of results from laboratory experiments to natural systems where highly turbulent natural convection plays a role. These include in particular convection in the ocean at the pole or in subglacial lakes (Couston & Siegert 2021), resulting in Rayleigh numbers beyond $10^{14}$ , and in the `ultimate regime’ of convection where ${\textit{Nu}} \sim Ra^{1/2}$ (with logarithmic corrections) (Kraichnan 1962). However, experimental heat-transfer measurements published in the literature for ${\textit{Ra}} \gt 10^{11}$ seem to be contradictory, some showing a transition at ${\textit{Ra}} \approx 10^{11}$ (Chavanne et al. 1997; Niemela & Sreenivasan 2003a ), some showing a delayed transition at higher ${\textit{Ra}}$ (He et al. 2012) or no transition at all (Niemela et al. 2000; Urban et al. 2019). These experiments are within Oberbeck–Boussinesq conditions (Roche et al. 2010) and share similar geometries, but differ in Prandtl numbers and details of the set-up. New recent theoretical models point out that the transition to the ultimate regime may be a subcritical process, and may depend on aspect ratio and Prandtl numbers (Roche 2020; Ahlers et al. 2022; Lohse & Shishkina 2023).

In this work, we evidence a way to recover universality. We report new velocity data obtained from shadowgraph imaging in parallelepiped Rayleigh–Bénard cells with or without roughness and using either deionised water or Fluorocarbon FC-770 as the working fluid, as well as new Direct Numerical Simulations (DNS), carried out in a similar geometric configuration filled with water. The Rayleigh numbers range between $5\,\times 10^{8}$ and $3\,\times 10^{12}$ .

In these conditions, the flow is known to be able to reach various regimes and several possible scaling exponent for ${\textit{Nu}}$ versus ${\textit{Ra}}$ . This is a range of parameters where apparent disagreement is observed between the data from Grenoble, Oregon, Brno and Göttingen cited above. And yet, we show that all the data, both from our new measurements and simulations, and from the literature, can be collapsed using the appropriate scaling of Reynolds and Prandtl numbers.

2. Experimental and numerical set-up

The cell dimensions are $41.5\,\textrm {cm} \times 10.5\,\textrm {cm} \times 41.5\,\textrm {cm}$ for the aspect ratio 1 cell, and $41.5\,\textrm {cm} \times 10.5\,\textrm {cm} \times 20\,\textrm {cm}$ for the aspect ratio 2 cell. The details of the cell are described in Méthivier et al. (2021). The Rayleigh number spans a wide range between $10^{9}$ and $2.5\times 10^{12}$ , the Prandtl number between 4.3 and 7.0 in water and 11 and 14 in FC-770 and the aspect ratio is either 1 or 2. The experimental cells can be set up with roughness on the bottom plate, consisting in square roughness elements of height $h_0 = 2\, \textrm {mm}$ , machined in the plate, as described in Belkadi et al. (2020).

DNS are performed for two cavities at ${\textit{Pr}}=4.4$ (same as water), with depth to height aspect ratios $\varGamma ^\star$ of one half and one quarter, and a width to height aspect ratio $\varGamma = 1$ . Three surface states for the bottom wall have been modelled: smooth, rough with $h_0/H=0.04$ , rough with $h_0/H=0.003$ ). ( $h_0/H$ is only for rough plates). Rayleigh numbers range from $5\,\times 10^{8}$ to $10^{12}$ . The computational grids resolve the Kolmogorov scale in the whole domain, and follow the resolution criteria set in Shishkina et al. (2010). The simulation at ${\textit{Ra}} = 10^{12}$ is performed on a $3328\, \times 896\, \times 5760$ grid. Statistical convergence of the integral heat-flux quantities is less than 2 %. Details of the numerical code, called SUNFLUIDH, can be found in Belkadi et al. (2021).

3. Heat-flux estimates

The experimental smooth cell (later referred to as `SS’) has entirely smooth boundaries, as described in Méthivier et al. (2021). For both experiments and DNS, in the rough case, the cell (later referred to as `RS’), is asymmetrical, with a rough bottom plate and a smooth top plate.

As shown in our previous works (Tisserand et al. 2011; Rusaouen et al. 2018), this allows us to define Rayleigh and Nusselt numbers of the rough and smooth half-cells, provided the bulk temperature $T_{\textit{bulk}}$ is also known. The temperature drop across the boundary layer of the smooth plate, $\Delta T_s = T_{\textit{bulk}} - T_{\textit{top}}$ , and the temperature drop across the boundary layer of the rough plate, $\Delta T_r = T_{\textit{bottom}} - T_{\textit{bulk}}$ , allow us to estimate the total temperature difference across a corresponding virtual symmetric cell with identical smooth plates as $2\Delta T_s$ , or with identical rough plates as $2\Delta T_r$ , from which ${\textit{Ra}}_s, {\textit{Nu}}_s, {\textit{Ra}}_r, {\textit{Nu}}_r$ are derived as

(3.1) \begin{align} {\textit{Ra}}_s &= \frac {2g\alpha (T_{\textit{bulk}}-T_{\textit{top}})H^3}{\nu \kappa },\\[-12pt] \nonumber \end{align}

(3.2) \begin{align} {\textit{Nu}}_s &= \frac {\dot {q}\!H}{2\lambda (T_{\textit{bulk}}-T_{\textit{top}})},\\[-12pt] \nonumber \end{align}

(3.3) \begin{align} {\textit{Ra}}_r &= \frac {2g\alpha (T_{\textit{bottom}}-T_{\textit{bulk}})H^3}{\nu \kappa },\\[-12pt] \nonumber \end{align}

(3.4) \begin{align} {\textit{Nu}}_r &= \frac {\dot {q}\!H}{2\lambda (T_{\textit{bottom}}-T_{\textit{bulk}})}.\\[12pt] \nonumber \end{align}

We have now operated the experimental rough configuration (RS) with Fluorocarbon FC-770, thus significantly extending the previous heat-transfer results in water (Salort et al. 2014). Figure 1 shows the new experimental data, as well as the new DNS data, compared with previous data. The Grossmann–Lohse model is also shown for comparison (with the prefactors from Stevens et al. 2013), as well as the experimental data of Wei et al. (2014) obtained in a symmetrical rough cell (`RR’).

Figure 1. Heat-transfer measurements in experiments and in the numerical simulations. The RS cell in water from Salort et al. (2014). The SS cells in FC-770 from Méthivier et al. (2021). The SS cell in water, RS cell in FC-770 and DNS: new data. Solid black line: ${\textit{Ra}}^{1/2}$ . The effective exponent is closer to ${\textit{Ra}}^{0.39}$ due to the logarithmic corrections. Purple down-pointing triangles: SS cell of Wei et al. (2014) (top and bottom half-cells). Dashed black line: Grossmann–Lohse model (Stevens et al. 2013).

In the fluorocarbon cell, the thermal boundary layer, estimated as $\delta _{\textit{th}} = H / (2 \textit{Nu})$ , ranges between $250$ and $500\,\unicode{x03BC} \textrm {m}$ , much smaller than the roughness height $h_0$ . The kinetic boundary layer, estimated as $\delta _v \approx \delta _{\textit{th}}{\textit{Pr}}^{1/3}$ , ranges between $600$ and $1100\,\unicode{x03BC} \textrm {m}$ , smaller than the roughness height $h_0$ . Therefore, the rough cell lies in regime III as defined by Xie & Xia (2017). By the same definition, DNS configurations can be either in regime II or regime III. At ${\textit{Ra}}=10^{10}$ , the simulation with the smaller roughness size ( $h_0^\star = 0.003$ ) is in regime II (i.e. the thermal boundary layer is smaller than the roughness height $h_0$ , whereas the kinetic one is larger), while the simulation with the larger roughness size ( $h_0^\star = 0.04$ ) is in regime III. As the thickness of the thermal boundary layer decreases with the Rayleigh number, there is a threshold beyond which the roughness becomes active and induces an enhancement of heat transfer.

Even when the rough half-cell has a ${\textit{Nu}} \sim Ra^{1/2}$ scaling (with logarithmic corrections), the Nusselt number in the smooth half-cell remains in good agreement, over more than 3 decades of Rayleigh number, with heat-transfer measurements in the SS cell, as well as with the Grossmann–Lohse model (Stevens et al. 2013). This shows that the top and bottom plates remain independent, even very far from the roughness-triggered enhancement threshold, and this holds even when the Nusselt number of the bottom half is nearly twice as large as that of the top half. This result is in agreement with our previous measurements (Tisserand et al. 2011; Rusaouen et al. 2018), and extends them further from the roughness-enhanced transition. This is a major result that shows how the boundary layers adapt themselves to maintain the heat flux as a constant. It is different from the set-up of Wei et al. (2014) where asymmetry between top and bottom plates is observed, although this is also the case in their SS case.

4. Velocity estimates

In experiments, we computed the velocity fields using Correlation Image Velocimetry (CIV) on shadowgraph recordings in a wide range of ${\textit{Ra}}$ and ${\textit{Pr}}$ , in both for RS and SS cells, yielding consistent mean velocity and Reynolds number estimates for all configurations. The method consists in deriving velocity estimates from the shadowgraph pattern using the same algorithm we use for Particle Image Velocimetry (PIV) images, i.e. the image is divided into smaller boxes that are correlated at $t$ and $t + \Delta t$ , and the displacement that produces the maximum correlation is used as an estimate for the local velocity. In practice, we use a two-pass algorithm from FluidImage (Augier, Mohanan & Bonamy 2019) with initial box size $128\,\textrm {px}\times 128\,\textrm {px}$ and a multipass zoom coefficient of 2, except for the RS FC-770 experiments which use an initial box size $64\,\textrm {px}\times 64\,\textrm {px}$ .

Figure 2. Snapshots of experimental shadowgraph recordings (a,c), and mean velocity fields computed from CIV (b,d), in both SS and RS experimental configurations in water: SS water ${\textit{Ra}}=6.3\,\times 10^{10}$ , $\textit{Re} = 1.6\,\times 10^{4}$ ; RS water ${\textit{Ra}}= 5.6\,\times 10^{10}$ , $\textit{Re} = 1.4\,\times 10^{4}$ .

Figure 3. Snapshots of experimental shadowgraph recordings (a,c), and mean velocity fields computed from CIV (b,d), in both SS and RS experimental configurations in FC-770: SS FC-770 ${\textit{Ra}}= 1.7\,\times 10^{12}$ , $\textit{Re} = 3.3\,\times 10^{4}$ ; RS FC-770 ${\textit{Ra}}= 1.6\,\times 10^{12}$ , $\textit{Re} = 4.0\,\times 10^{4}$ .

Figure 4. Snapshot of the experimental shadowgraph recording in the SS cell in water at ${\textit{Ra}}= 6.3\,\times 10^{10}$ , after histogram equalisation with the CLAHE method. It is the same snapshot as in figure 2 (a,c), but allows us to visualise that plumes fill the full cavity and the pattern is almost homogeneous, once the inhomogeneity of contrast has been removed.

Examples of raw shadowgraph patterns and mean velocity fields are shown in figures 2 and 3, for SS and RS cells. Using the shadowgraph image for CIV is possible at high ${\textit{Ra}}$ , because plumes fill the full cavity, and the correlation algorithm works everywhere. Plumes at the centre of the image are not visible on the raw shadowgraph because the contrast is much lower in the centre of the cell. However, we use a high-dynamics 14 bit camera, which allows us to resolve the plume pattern everywhere, including at the centre. Figure 4 shows a processed image for which the contrast limited adaptive histogram equalisation (CLAHE) method, implemented in the OpenCV Library Bradski (2000), has been applied. It is even more difficult to guess the mean flow pattern, but it clearly shows how homogeneous plume patterns actually are, once spatial variation of contrast is removed. The same holds even more for experiments in FC-770, but the patterns are so much smaller that they are difficult to render on a paper-sized figure.

This velocimetry method, based on the shadowgraph pattern, has previously been validated against standard PIV in experimental conditions where both were possible (Méthivier et al. 2021). One should be aware that the shadowgraph patterns result from the deviation of light across the cell depth. Therefore, the obtained velocity is integrated over the depth, in contrast to PIV which is computed on a plane.

In numerical simulations, slices at mid-plane are saved, allowing us to compute the mean velocity field at mid-depth, see figure 5. In both experiments and numerics, the large-scale Circulation (LSC) may go clockwise or anti-clockwise. To make comparison easier, we applied a mirror to the anti-clockwise fields before plotting, so that all fields appear clockwise. The DNS field at ${\textit{Ra}}=10^{10}$ is in good agreement with the PIV obtained previously (Liot et al. 2017).

The experimental estimates from CIV on the shadowgraph pattern are fairly similar to the mid-depth mean velocity estimates from the DNS, except in the bottom right and top left corners where the corner flows are highly three-dimensional: both rising and falling structures can be observed in the shadowgraph sequence, hence the CIV estimate tends to go to zero.

We find that the mean velocity field has the same structure at ${\textit{Ra}} = 10^{10}$ and ${\textit{Ra}} = 10^{12}$ , and also does not change significantly in rough cells. The raw shadowgraph pictures show, however, stronger density gradients, and smaller structures, filling the cell further from the walls. Beyond ${\textit{Ra}}=10^{12}$ , the shadowgraph pattern is much less inhomogeneous. When roughness is added to the bottom plate, mean velocity is increased between $5\,\%$ and $20\,\%$ in the range of parameters explored in this work, but the overall structure of the mean flow is not significantly changed. However, it does not preclude stronger effects on the velocity fluctuations, as shown by Liot et al. (2017).

Figure 5. Mean velocity fields from the DNS at mid-depth for both depth to height aspect ratios $\varGamma ^{\star }$ (see text, § 2). (a--c) Rayleigh number dependence. (d--f) comparison SS vs RS.

In the following, we focus on the velocity of the upwelling and downwelling jets to allow comparison with the literature where velocity is estimated with pairs of thermometers at mid-height (Wu & Libchaber 1992; Chavanne et al. 2001; Niemela et al. 2001; He et al. 2015; Musilová et al. 2017). To make the Reynolds number values quantitatively comparable, the values from Wei et al. (2014) have been corrected. Indeed, they defined their Reynolds numbers as

(4.1) \begin{equation} \textit{Re} = \frac {V_{{\textit{osc}}}H}{\nu }, \qquad V_{{\textit{osc}}} = \frac {4H}{\tau }, \end{equation}

where $\tau$ is an oscillation period computed from the autocorrelation of local thermistors. However, as shown by Niemela & Sreenivasan (2003b ), the length of the roll is not $4H$ but $\alpha \boldsymbol{\cdot }H$ , where $\alpha$ is between 2 and 4, and strongly depends on ${\textit{Ra}}$ for ${\textit{Ra}}\lt 10^{11}$ . In the range of Rayleigh numbers of Wei et al. (2014), $\alpha \approx 3.4$ , and therefore we apply a corrective prefactor $3.4/4.0$ to their Reynolds numbers. With this correction, the Reynolds numbers of their SS cell collapse fairly well with those of the literature and of the present work. The Reynolds numbers of their rough cells are higher.

In symmetrical cells, we average the velocity of the downwelling and upwelling jets. In RS cells, we associate the downwelling jet with the top plate, and the upwelling jet with the bottom plate. Indeed, we use the following definition for the Reynolds numbers, ${\textit{Re}}_{r,s}$ , of the rough bottom and the smooth top:

(4.2) \begin{equation} {\textit{Re}}_{r,s} = \frac {U_{{\textit{up}},{\textit{do}w\textit{n}}}H}{\nu }, \end{equation}

where $U_{{\textit{up}},{\textit{do}w\textit{n}}}$ is the maximum of the profile along $x$ of vertical velocity in the plane at mid-height.

The Reynolds numbers in the various configurations, as well as those of the literature, show deviations of more than a factor 2, even when differences in the definition have been accounted for, see figure 6(a). These deviations, although relatively small, are significant for quantities that scale with a power of the Reynolds number. For example, the Reynolds number in the RS cell using fluorocarbon is larger than in the SS case, and additionally is asymmetric (velocity is larger near the rough plate). This Reynolds number enhancement triggered by roughness was not visible in our previous measurements using PIV in water at lower Prandtl and Rayleigh numbers (Liot et al. 2017). Indeed, the Rayleigh number of the new data obtained in fluorocarbon is a decade larger and further from the threshold where roughness-enhanced heat transfer is triggered (regime II) than the rough cell using water. This is also observed in the DNS results. The higher Prandtl number also further separates the thermal and viscous boundary layers, which probably plays a role in the boundary layer response to plate roughness.

Figure 6. (a) Reynolds number measurements in FC-770 (green), water (blue) and from DNS (red), in RS (squares) and SS (circles) cells. Data from Wu & Libchaber (1992), Chavanne et al. (2001), Niemela et al. (2001), Brown, Funfschilling & Ahlers (2007), Wei et al. (2014), He et al. (2015) and Musilová et al. (2017) are plotted for comparison (triangles). A corrective factor has been applied to the Reynolds numbers of Wei et al. (2014) (see text). (b) Friction coefficient for all data points (same symbols). The heat flux, ${\textit{Ra}}\textit{Nu}$ , collapses for all data with ${\textit{Ra}}\textit{Nu} \sim 0.2 {\textit{Re}}^3{\textit{Pr}}^2$ .

5. Discussion

One interesting quantity is the friction coefficient

(5.1) \begin{equation} \frac {\textit{RaNu}}{{\textit{Re}}^3{\textit{Pr}}^2}, \end{equation}

which was used by Chavanne et al. as an indicator of the transition to turbulence in the boundary layers (Chavanne et al. 2001), and can alternatively be interpreted as proportional to the ratio, $\mathcal{R}_{\epsilon }$ , between the kinetic energy dissipation rate $\epsilon$ and the contribution of the bulk in kinetic energy dissipation rate $\epsilon _{u, bulk}$ (Grossmann & Lohse 2000)

(5.2) \begin{equation} \epsilon = \frac {\nu ^3}{H^4}(\textit{Nu} - 1){\textit{RaPr}}^{-2} = \epsilon _{u, \textit{BL}} + \epsilon _{u, \textit{bulk}}, \end{equation}

and in the case of a turbulent bulk

(5.3) \begin{equation} \epsilon _{u,\textit{bulk}} \propto \frac {\nu ^3}{H^4}{\textit{Re}}^3. \end{equation}

Hence, in the turbulent regime,

(5.4) \begin{equation} \mathcal{R}_{\epsilon } = \frac {\epsilon }{\epsilon _{u, bulk}} \sim \frac {\textit{NuRa}{\textit{Pr}}^{-2}}{{\textit{Re}}^3}. \end{equation}

In the range of parameters explored in this work, we evidence a transition at a critical Reynolds number, ${\textit{Re}}_c \approx 10^{4}$ . Beyond this transition, the friction coefficient no longer depends on the Reynolds number, or does so only very weakly, see figure 6(b). In the DNS, the ratio $\mathcal{R}_{\epsilon } = \epsilon /\epsilon _{u, bulk}$ can be estimated directly from the three-dimensional velocity fields using the local definition of $\epsilon$

(5.5) \begin{align} \epsilon (\boldsymbol {x}, t) &= \frac {\nu }{2}\sum _{i=1}^3\sum _{j=1}^3{\left ( \partial\! _ju_i + \partial _iu\!_j\right )}^2, \end{align}

(5.6) \begin{align} \epsilon &= \frac {1}{V}\int _V\left \langle \epsilon (\boldsymbol {x}, t)\right \rangle _t \textrm{d}^3\boldsymbol {x}, \end{align}

(5.7) \begin{align} \epsilon _{u, \textit{bulk}} &= \frac {1}{V}\int _{V_{\textit{bulk}}}\left \langle \epsilon (\boldsymbol {x}, t)\right \rangle _t \textrm{d}^3\boldsymbol {x}, \\[9pt] \nonumber \end{align}

where $V = \varGamma\! \varGamma ^{\star }H^3$ is the volume of the convection cell and $V_{\textit{bulk}}$ is the volume of the bulk.

Figure 7. Inverse of the bulk fraction of kinetic dissipation, $R_{\epsilon }$ , in the DNS estimated from the three-dimensional velocity gradients (see text).

The actual value of the ratio $\mathcal{R}_{\epsilon }$ depends on the definition of the bulk region. In the following estimates, we use a practical definition based on the efficiency of mixing in the turbulent bulk: we define the bulk as the region where the mean temperature gradient vanishes (in practice $\partial _z\theta \lt 10^{-3}$ ). Indeed, in the turbulent bulk, mixing makes the mean temperature homogeneous.

The values are shown in figure 7. The plot is similar to that in figure 6(b). The value at the highest $\textit{Re}$ when the plateau is reached, of order $\mathcal{R}_{\epsilon } \sim 1.3$ , suggests that more than $75\,\%$ of the kinetic dissipation occurs outside the boundary layers, and confirms that the total dissipation is dominated by the turbulent bulk.

This transition may correspond to the mixing transition in fully developed turbulent flows described by Dimotakis (2000), which is thought to occur for an outer-scale Reynolds number of order $10^{4}$ (or $R_{\lambda }$ of order $10^{2}$ ). In his paper, he suggested that this mixing transition coincides with the soft turbulence –hard turbulence transition predicted by Castaing et al. (1989) in turbulent thermal convection. In their original paper, Castaing et al. (1989) claimed that they observed the soft–hard transition at a lower Rayleigh number of order $4\,\times 10^{7}$ , based on the observation of a change in the temperature fluctuation distribution (Heslot, Castaing & Libchaber 1987), which was later found to not be a good proxy. We propose that the hard turbulence regime, which corresponds to the mixing transition, is actually reached at higher Reynolds numbers, when the friction coefficient reaches a plateau.

More generally, this plateau is a signature of the dissipation anomaly in turbulence (Dubrulle 2019), which was also observed numerically by Pandey et al. (2022) in turbulent Rayleigh–Bénard convection, and has been observed in various turbulent systems, such as grid flows where it occurs at $R_{\lambda } \sim 50$ (Sreenivasan 1984), or von Kármán flows where estimates based on torque measurements yields a transition at $\textit{Re} \sim 300$ (Saint-Michel et al. 2014). Our results show that the plateau is very robust and consistent with both numerical and experimental observations in turbulent Rayleigh–Bénard convection, as well as previous results in the literature.

The bulk of Rayleigh–Bénard convection in this range of forcing is completely determined by turbulent mixing, regardless of whether or not the Nusselt versus Rayleigh number relationship shows a change in their scaling relationship. Note that, while the critical Reynolds number appears to be universal, as well as the plateau in the turbulent regime, the path towards turbulence does not have to be unique, as noted by Dimotakis (2000). While the data for $\textit{Re} \lt {\textit{Re}}_c$ reveal some non-universality, the important point is that all the data collapse beyond ${\textit{Re}}_c$ .

Therefore, we get a relationship between the dimensionless heat flux, ${\textit{Ra}}\textit{Nu}$ , and the Reynolds and Prandtl numbers, where

(5.8) \begin{equation} \textit{RaNu} = \frac {g \rho ^2\alpha c\!_p H^4}{\mu \lambda ^2} \dot {q}, \end{equation}

where $\rho$ is the density, $c\!_p$ is the specific heat capacity and $\mu = \rho \nu$ is the dynamical viscosity. The phenomenological relationship is

(5.9) \begin{equation} \textit{RaNu} \approx 0.2 {\textit{Re}}^3{\textit{Pr}}^2. \end{equation}

While this scaling has been reported previously for individual sets of experiments (Chavanne et al. 2001; Musilová et al. 2017), the main result of this work is that the prefactor is universal, and does not change when the ${\textit{Nu}}$ vs ${\textit{Ra}}$ scaling law is changed, either by roughness or by the transition to the ultimate regime, hence collapsing data from experiments thought to be in quantitative disagreement.

Indeed, it is remarkable that this relationship, (5.9), holds for the Grenoble data which evidence a transition to the ultimate regime, as well as data which do not evidence such a transition, and also in the case of rough cells in which ${\textit{Nu}}\sim Ra^{1/2}$ , as well as smooth cells with a classical scaling. This shows that all these experiments lie in the turbulent regime, nevertheless, at moderate Prandtl numbers, ${\textit{Nu}}(Ra)$ it is not unique because the Nusselt number depends on $\epsilon _{\theta }$ , which in these regimes depends on the local form of the thermal boundary layer which has certainly a strong interaction with LSC but also depends on boundary details. Although this should not be a surprise, as theories need to set both $\epsilon$ and $\epsilon _{\theta }$ to make a prediction on the Nusselt number (Kraichnan 1962; Castaing et al. 1989; Shraiman & Siggia 1990; Grossmann & Lohse 2000, 2011), the critical shear Reynolds number is often brought up as a control parameter for the onset of the ultimate regime (Stevens et al. 2013; Lohse & Shishkina 2023).

The present result shows that this may not be the discriminating criterion to understand why several datasets in the literature do not evidence a transition, or evidence a delayed one. Indeed, all these high- ${\textit{Ra}}$ seemingly conflicting datasets already lie beyond the transition to turbulence where the kinetic dissipation rate scales like ${\textit{Re}}^3$ . The condition on the Reynolds is therefore fulfilled. It appears, however, that this is not sufficient.

In this work, we have used the outer Reynolds number as a proxy, while it would be more physically appropriate to use $R_{\lambda }$ . The reason is that we do not have access to $R_{\lambda }$ for most of the data discussed in this paper. However, since the outer Reynolds number is determined by the wind, and the Reynolds number is not the discriminating factor to understand the difference between experiments which evidence the ultimate regime and those that do not, it follows that the mean wind is not the limiting factor for the onset of the ultimate regime in these experiments.

This is consistent with previous works: (i) the onset of the transition to the ultimate regime in Grenoble is not modified when the mean wind is changed (Roche et al. 2010), (ii) the velocity profiles have a logarithmic profile even at relatively low Rayleigh numbers (Liot et al. 2016a , b ) in the classical regime and (iii) the influence of the mean wind intensity on the global heat transport is very weak (Ciliberto, Cioni & Laroche 1996; Ahlers, Grossmann & Lohse 2009), even when the wind is highly depleted (Carbonneau et al. 2025), suggesting that a Reynolds number based on fluctuations would be more relevant than one based on the mean wind, in agreement with Roche (2020).