Classical 1/3 scaling of convection holds up to Ra = 1015.

The global transport of heat and momentum in turbulent convection is constrained by thin thermal and viscous boundary layers at the heated and cooled boundaries of the system. This bottleneck is thought to be lifted once the boundary layers themselves become fully turbulent at very high values of the Rayleigh number [Formula: see text]-the dimensionless parameter that describes the vigor of convective turbulence. Laboratory experiments in cylindrical cells for [Formula: see text] have reported different outcomes on the putative heat transport law. Here we show, by direct numerical simulations of three-dimensional turbulent Rayleigh-Bénard convection flows in a slender cylindrical cell of aspect ratio [Formula: see text], that the Nusselt number-the dimensionless measure of heat transport-follows the classical power law of [Formula: see text] up to [Formula: see text] Intermittent fluctuations in the wall stress, a blueprint of turbulence in the vicinity of the boundaries, manifest at all [Formula: see text] considered here, increasing with increasing [Formula: see text], and suggest that an abrupt transition of the boundary layer to turbulence does not take place.

Turbulent Rayleigh–Bénard convection (RBC), triggered in a fluid layer which is uniformly heated from below and cooled from above, is considered the paradigm for buoyancy-driven turbulence reaching from stellar convection (1) via atmospheric (2) and oceanic (3) turbulence to cooling blankets in nuclear engineering (4) or the storage of renewable energy in liquid metal batteries (5). A central question in RBC is the fundamental law governing the global turbulent transport of heat and momentum across the convection layer as a function of dimensionless parameters characterizing the flow (6–11). These parameters are the Rayleigh number , which quantifies the relative magnitude of the thermal driving to the viscous and diffusive forces of the fluid motion and is directly proportional to the temperature difference , maintained between the top and bottom plates, and the Prandtl number , which is the ratio of momentum and thermal diffusivities denoted by and , respectively. Here, , , and are the acceleration due to gravity, the isobaric thermal expansion coefficient, and the fluid layer height , respectively. A third parameter is the aspect ratio for closed (cylindrical) containers of diameter , the typical setup in which turbulent convection is studied in the laboratory. For Rayleigh numbers sufficiently above the onset value of , the fluid flow develops thin thermal and viscous boundary layers at the top and bottom walls that sandwich a turbulent bulk layer with a well-mixed temperature of . The heat which is supplied at the bottom has to be transmitted through the boundary layers that act as bottlenecks, which are thought to be removed when the thin boundary layers become turbulent themselves (12, 13), leading, potentially, to the so-called “ultimate” state of turbulent convection. It is not known precisely at what Rayleigh numbers this transition is presumed to occur, but high-Rayleigh-number laboratory experiments in cylindrical cells (14–20) have reported different outcomes on the turbulent transport law for .

In this paper, we assess the laws of turbulent heat and momentum transport in RBC by direct numerical simulations (DNS) of the Boussinesq equations that couple the (turbulent) velocity field with the temperature field , with . We cover Rayleigh numbers in the range at the Prandtl number held fixed at unity, such that the momentum and thermal boundary layers have similar thicknesses and scale the same way with . We used massively parallel computations that apply the spectral element method (21, 22) and experimented with increasing grid resolutions and finer time steps until the veracity of the data was established (see and Table 1 for more details). Despite the heftiness of the computations, our goal of pushing toward the highest possible Rayleigh numbers meant that we had to compromise on the aspect ratio , which we held at . Even for this low aspect ratio, the boundary layer thickness at the highest is about the diameter of the cylindrical cell, and the heat transport results agree closely, where they overlap, with previous results for aspect ratio unity, obtained using the same numerical methods (23). A comparison with the experimental data of Niemela and Sreenivasan (15) at also shows that the two are quite close.

Table 1.

Summary of spectral element simulation parameters

Ra	Nu	Re	t	M	N	N_e
10^8	29.94±0.04	601±1	245	46	3	19,200
10^9	58±10	1,920±14	396	89	5	19,200
10^10	107±11	5,570±50	419	104	5	19,200
10^11	229±13	17,000±200	136	58	9	19,200
10^12	503±25	52,500±2,000	58	99	11	537,600
10^13	1,075±44	133,000±2,000	24	83	13	925,200
10^14	2,228 ± 100	410,000 ± 10,000	14	83	11	17,145,600
10^15	4,845 ± 200	1,094,000 ± 10,000	7	82	11	17,145,600

The table lists the Rayleigh number Ra, the Nusselt number Nu obtained from heat flux through the top/bottom plates, the Reynolds number Re, the total simulation time t in free-fall time units t , the number of three-dimensional simulation snapshots M, the polynomial order N of the Lagrangian interpolation polynomials for each of the three space directions, and the total number of spectral elements Ne. This results in a total number of mesh cells of N × N which results in more than 2.2 × 1010 mesh cells for the biggest simulation run. The Prandtl number is unity for all cases. Simulation time t is the time in a statistically stationary state at the highest spectral resolution in each run. It is this time that is used for the statistical analysis of the main text.

Global Heat and Momentum Transfer

Our principal results are shown in Fig. 1. In Fig. 1 , we plot the Nusselt number against the Rayleigh number . The Nusselt number is obtained for each integration time step as an average of the vertical diffusive heat current taken at the bottom and top walls over the cell cross-section . The data at bottom and top walls are subsequently combined in a time average to give and the associated error bars. In , we compare this method of determination with a second, which is based on the combined volume–time average of the sum of the convective and diffusive heat currents, and find very good consistency. The earlier simulation results at lower Rayleigh number for and (23) are also shown; together, the data cover almost 10 orders of magnitude in with an overlap of the two records for . Fig. 1 shows that a good power-law scaling exists for , with the exponent given by . For lower Ra, the exponent has a smaller value of , as is well known from past considerations (11) and which is close to 2/7 (8). These data are in good agreement with the experimentally measured Nusselt numbers in cryogenic helium (15), as shown in . Recent experiments by de Wit et al. (24) in a water column at for report .

Fig. 1.

Global scaling laws of turbulent transport. (A) Turbulent heat transport law for two datasets on log–log coordinates. The present data (open circles) are for and for . Open squares representing DNS at and for are taken from ref. 23. Also added is the power-law fit of which is obtained for . (B) Linear–log plot of the data of A, displayed in the form of compensated classical power law of slope, with error bars computed from the standard deviation of the time series , obtained at (bottom) and (top). The plot shows that the power-law exponent is smaller () for and is in the high-Ra range, with no tendency to a larger slope as would be required of the possible progression to the ultimate state. The dashed line corresponds to prefactor 0.0525, in close agreement with the classical prediction of 0.073 (6, 7). (C) Compensated power-law plot of the turbulent momentum transport law . A fit over the same Rayleigh number range as in B gives . Symbols mean the same in A–C.

Fig. 1 replots the Nusselt number compensated by the scaling, . The classical scaling is satisfied over five decades of starting from . We stress that the prefactor of is also close to the theoretical prediction of 0.073 (6, 7) based on the assumption that marginally stable boundary layers maximize the heat transport. The present Nusselt numbers match, within error bars, with the earlier data of ref. 23 for the cell aspect ratio of unity, suggesting that the present aspect ratio, although small, is not too small for purposes of heat transport (see discussion of Fig. 2).

Fig. 2.

Convection flow inside the boundary layers. Snapshots of the temperature and velocity fields are displayed for three different Rayleigh numbers. (A–C) Temperature field at a Rayleigh number of (A) in a horizontal plane at , (B) for , and (C) for . The heights of these planes correspond to half of the thermal boundary layer thickness, . The temperature range is given by the color bar below each panel. (D–F) Streamlines of the velocity field which is projected onto the cutting plane, that is, shown at (D) , (E) , and (F) . The distance from the bottom wall is the same as in the corresponding temperature plots. (G–I) Contour plot of the wall shear stress magnitude for at (G) , (H) , and (I) .

The turbulent momentum transport is given by the Reynolds number (see for details). Throughout this work, tildes denote dimensionless physical quantities expressed in characteristic units composed of height , free-fall velocity , and/or the wall-to-wall temperature difference . Times are thus given in units of the free-fall time . Fig. 1 shows the compensated plot for the two sets of data. The Reynolds number scaling in the slender cell is given by for , with no apparent tendency toward a qualitative change toward the highest values of . For , this same power law has an exponent of (23); in the overlapping region of , the magnitudes are larger by a factor of 3 to 5. We shall discuss this result toward the end. Note that a scaling with an exponent of does not necessarily herald the “ultimate” convection state.

Boundary Layer Structure

Fig. 2 illustrates the complex structure of the turbulent fields inside the thermal boundary layers for three values of the Rayleigh number. The mean thickness of the thermal boundary layer is given by . Fig. 2 shows horizontal cuts through the temperature field and the velocity field at . Fig. 2 G–I shows contours of the magnitude of the wall shear stress field at the bottom plate, defined in RBC for in characteristic units asFor all three datasets, the fields develop ever-finer structures and filaments as the Rayleigh number increases from to . The wall shear stress magnitude (Fig. 2, G–I) and velocity field streamlines (Fig. 2, D–F) illustrate the complexity of interactions of the flow near the bottom plate with the fine ridges of thermal plumes (Fig. 2, A–C).

Fig. 3 shows the average behaviors of the thermal and the mean-square velocity fields. Fig. 3 shows the mean temperature profiles at various , Fig. 3 shows the corresponding boundary layer behavior, and Fig. 3 shows the root-mean-square of the velocity fluctuations. They are very similar in behavior to those in higher aspect ratios.

Fig. 3.

Statistical analysis within the boundary layer. (A) Mean dimensionless temperature profiles versus cell height . The profiles are obtained as averages over the circular cross-section at each and time. (B) Zoom-in of the near-wall region. Data are the same as in A, but replotted now as with . An additional arithmetic average over both halves of the mean profiles is performed to improve statistics. This is possible due to the up–down symmetry in the Boussinesq approximation with respect to the midplane. The range displayed here corresponds to , with the mean thermal boundary layer thickness taken at . (C) Same zoom of the mean velocity fluctuation profiles in units of the free-fall velocity . The root-mean-square value is determined with respect to all three velocity components, and profiles are obtained in the same way as in B. The legend is the same for B and C. (D) Probability density function (PDF) of the dimensionless wall shear stress magnitude given by Eq. . In each dataset, the argument is rescaled by the corresponding stress value of the maximum of the PDF, . All PDFs are also normalized to unity. (E) PDF of the single partial derivative at given in units of the corresponding root-mean-square value. The legend is the same as D. (F) Root-mean-square value versus Rayleigh number . In D–F, the analysis is limited to an inner section of the plate with .

Our conclusion, then, is that the behavior of the temperature field and the fluctuating velocity field in the present simulations is very similar to those at higher aspect ratios, essentially because the boundary layers are very thin compared to the cell diameter. This is illustrated further in Fig. 4 by a direct comparison of boundary layer snapshots at aspect ratios 1 and at . All of the heat supplied at the bottom plate has to pass through thin boundary layers; as already stated, this bottleneck is thought to be removed when the boundary layers become fully turbulent themselves, leading to the “ultimate” state (12). However, the boundary layer dynamics in RBC is distinct from that of the canonical isothermal boundary layers in channels or over flat plates with a unidirectional mean flow (25). The boundary layer motion in RBC lumps together hotter (colder) segments of the thermal boundary layer which detach from the bottom (top) plate and rise (fall) into the well-mixed bulk, thus disrupting the near-wall fluid motion by updrafts at all Rayleigh numbers; it is also subject to plume impacts and couples to the large-scale motion in the bulk. The perturbations in the convection boundary layer are large even at moderate Rayleigh numbers, and a distinct event, such as the classical boundary layer transition to turbulence, seems unlikely in convection. This expectation is supported in Fig. 3, which shows that the wall shear stress fluctuations possess no sudden changes with , and that their distributions at different can be collapsed very well between and with a linear tail for the smallest amplitudes. The vertical derivatives for a single velocity component at the plates show similar broadband behavior (Fig. 3), and their relatively large root-mean-square values grow continuously (Fig. 3), suggesting that the boundary layers at these fluctuation levels are already turbulent for Rayleigh numbers . This intrinsic reason undercuts the notion that the heat transport would change to a different behavior at higher Rayleigh numbers. We speculate that the argument may hold for larger aspect ratios as well.

Fig. 4.

Comparison of the boundary layer structure at different aspect ratios but same Rayleigh number. A, D, and G show temperature slices, and B, E, and H show the projected streamlines at . C, F, and I show the wall shear stress magnitude at similar to Fig. 2. The data in G–I correspond to , , and . A–C and zooms in D–F are for , , but (23). The zoom-up pictures with a box side length of show that they are very similar to the case of the lower aspect ratio.

Large-Scale Flow Organization

We now return to the Reynolds number results (see Fig. 1 again). The difference in the momentum transport between the results for (or larger aspect ratios) and the present results can be attributed to the particular large structure in the slender cylinder. Fig. 5 displays isosurfaces of the vertical velocity component. They indicate a helical large-scale structure (similar to a barber pole) for all three Rayleigh numbers shown; indeed, this structure persists throughout the whole range of starting from and serves some of the same purposes as the large-scale circulation in high aspect ratio convection, in linking the bottom and top walls directly. This helical flow structure collapses sometimes, and the winding number alters, but the feature is generally robust. This can be seen clearly at , for which two views are shown. As a result of this flow, the mean temperature profile obeys a nearly linear form in the midsection, just as for higher aspect ratio cases (Fig. 3). While this large-scale flow has some effect in transporting momentum, as seen by the disparity in Reynolds numbers in the overlap range between the present data and those of ref. 23, it seems to have negligible impact on the turbulent heat transport, as discussed earlier. Fig. 5 provides horizontal cuts through the cylinder in the central region, which confirm the appearance of ever-finer vertical velocity filaments that leave the large-scale flow essentially unaffected. As far as the small-scale motion is concerned, the aspect ratio of seems to have no obvious effect. Thus, the different large-scale flow structures in different aspect ratio cylindrical cells seem to matter for turbulent momentum transport but apparently not for the small structure and the turbulent heat transport.

Fig. 5.

Large-scale flow organization in slender cell and the increasingly small structure with increasing . (A) Front view of an isosurface plot showing a snapshot of the vertical velocity field component at the levels (blue, negative; red, positive) at . The symbol to the right indicates the specific plan view. (B) Replot of the data with the front view rotated by 180° with respect to A in order to highlight the helical nature of the large scale. (C and D) Isosurface plots for (C) and (D) for . (E and F) Contour plot of the vertical velocity field at (E) for , and (F) at for the same snapshot. (G and H) Contour plot of the vertical velocity field at for (G) and (H) . The vertical positions of the expanded horizontal planes in E–G are highlighted by a yellow ring in each of the front views (in A–D). The three snapshots also correspond to those displayed in Fig. 2.

In conclusion, our DNS show that, at least in the low- case, there is no tendency toward the “ultimate” state up to . We cannot exclude that the present large-scale flow shifts its appearance at higher , or that some new heat transport law might set in. However, the fact that the boundary layers at these Rayleigh numbers are already strongly fluctuating suggests that no new state of turbulence is likely to appear, thus making it unlikely that the heat transport law will depart from the classical power; this is consistent with the conclusion of ref. 15 from experiments in a cell of .

Data Availability.

Data are available upon request from the corresponding author.