## Abstract

High–Rayleigh number convective turbulence is ubiquitous in many natural phenomena and in industries, such as atmospheric circulations, oceanic flows, flows in the fluid core of planets, and energy generations. In this work, we present a novel approach to boost the Rayleigh number in thermal convection by exploiting centrifugal acceleration and rapidly rotating a cylindrical annulus to reach an effective gravity of 60 times Earth’s gravity. We show that in the regime where the Coriolis effect is strong, the scaling exponent of Nusselt number versus Rayleigh number exceeds one-third once the Rayleigh number is large enough. The convective rolls revolve in prograde direction, signifying the emergence of zonal flow. The present findings open a new avenue on the exploration of high–Rayleigh number turbulent thermal convection and will improve the understanding of the flow dynamics and heat transfer processes in geophysical and astrophysical flows and other strongly rotating systems.

## INTRODUCTION

Thermally driven turbulent flows widely occur in geophysical flows and industrial processes. Examples are thermal convection in the atmospheric and mantle convection (*1*, *2*), in the ocean (*3*), and in many industrial processes (*4*). Rayleigh-Bénard convection (RBC), a fluid layer heated from below and cooled from above, is an ideal model for the study of thermally driven turbulent flows (*5*–*7*). The main challenge for thermal turbulence studies is to explore the flow dynamics and heat transfer of the system in a wide range of the control parameters. The main attention has been on how the heat transfer depends on the Rayleigh number, which is the dimensionless temperature difference and measures the intensity of the thermal driving of the system. The Rayleigh number is defined as *g* is the gravitational acceleration, β is the isobaric thermal expansion coefficient, ν is the kinematic viscosity, κ is the thermal diffusivity, Δ is the temperature difference between the hot plate and the cold plate, and *L* is the thickness of the fluid layer between the aforementioned plates.

To achieve high Rayleigh numbers, numerous strategies have been proposed in the past decades. Looking at the definition of *Ra*, one can see that large *Ra* can be reached in several ways. One approach is to use a working fluid for which the parameter *8*–*11*) or mercury (*12*). Recently, it was found that pure gas, particularly gasses with high molecular weight and under high pressure (up to 19 bar), is also an effective way to push *Ra* to high values at a constant *Pr* number (*13*, *14*). Another commonly used approach in the field of turbulent thermal convection is to increase the system thickness *L* (*10*, *11*, *13*–*18*).

The fact that hitherto not much attention was on the changing the gravitational acceleration provides the motivation of the present work. In this study, we propose a novel system, annular centrifugal RBC (ACRBC), which is a cylindrical annulus with cooled inner and heated outer walls under a solid body rotation. In this system, the buoyancy force can be efficiently enhanced by replacing the gravitational acceleration (*g*) by the centrifugal acceleration, and consequently, *Ra* can be increased by increasing the rotation rate for a given fluid and temperature difference. Another advantage of this system is that the thermal convection in rapidly rotating cylindrical annulus has been recognized as a good model for the study of the flows in planetary cores and stellar interiors (*19*–*24*) and also in engineering applications under strong rotation, i.e., flow in gas turbines (*25*–*27*). The flow dynamics in this system can give insights into flows in geophysical context and flows in industrial processes. In this work, we aim to study the heat transfer and the turbulent structures in ACRBC.

## GOVERNING EQUATIONS, SIMULATIONS, AND EXPERIMENTS

### Governing equations and numerical simulations

Using the parameter definitions as shown in Fig. 1B, the governing dimensionless Boussinesq equations in a rotating reference frame (*28*) can be expressed as**u** is the velocity vector normalized by the free fall velocity *t* is the dimensionless time normalized by *T _{h}* −

*T*. Here, ω,

_{c}*R*,

_{i}*R*, and Δ are the rotational speed, the thickness of the fluid layer between the two cylinders, the outer radius of the inner cylinder, the inner radius of the outer cylinder, and the temperature difference between the hot (

_{o}*T*

_{h}) and cold (

*T*

_{c}) cylinders, respectively.

From the above governing equations, the relevant control parameters in ACRBC are the Rayleigh number (characterizing the thermal driving strength) *Pr* = ν/κ, and the radius and aspect ratio (geometric properties) η = *R _{i}*/

*R*and Γ =

_{o}*H*/

*L*. The key response parameter of the system is the Nusselt number (measuring the ratio of the total heat transport over the conductive one)

*Nu*= −

*Q*ln (η)/(αΔ2π

*H*). Here,

*Q*is the measured heat input through the outer cylinder into the system per unit of time; α = κρ

*c*is the thermal conductivity of the working fluid, with ρ and

_{p}*c*being the density and the specific heat capacity of the fluid, respectively; and

_{p}*H*is the height of the gap between two cylinders. It should be noted that the definition of

*Nu*in ACRBC is slightly different from that in classical RBC because of the cylindrical geometry with a heat flux in the radial direction, and the detailed derivations of the

*Nu*for ACRBC are documented in the Supplementary Materials.

Direct numerical simulations (DNS) are performed using an energy-conserving second-order finite-difference code (*29*–*31*). The code has been extensively validated and used in previous studies in both the Cartesian and cylindrical coordinates for convective systems (*30*–*33*). In all numerical simulations and experiments, the radius ratio is fixed at η = 0.5. No-slip boundary conditions were used for the velocity and constant temperature boundary conditions for the inner and outer cylinders. Periodic boundary conditions were used for the top and bottom surfaces. For three-dimensional simulations, the aspect ratio was chosen the same as in the experiments Γ = *H*/*L* = 1 (the only two exceptions are for the cases at *Ra* = 1.16 × 10^{9} and *Ra* = 2.2 × 10^{9}, where the height of the domain is *L*/4 and *L*/8, respectively, as the flows here are quasi–two-dimensional). For the cases of which the flows are quasi–two-dimensional at large *Ra*, we use two-dimensional simulations. *Pr* was fixed at 4.3 for all the simulations. The adequate resolution was ensured for all cases, for example, at *Ra* = 4.7 × 10^{8}, 4608 × 384 × 384 grid points were used and at *Ra* = 4.7 × 10^{10} with grid points 18,432 × 1536. To obtain sufficient statistics, each simulation was run at least 200 free fall time units. As reported in (*34*, *35*), the boundary layers in rotation system are expected to be thinner compared with classical RBC. Hence, the verification of boundary layer grid resolution is necessary. A posteriori check of grid resolution shows that at *Ra* = 4.7 × 10^{10}, there are 48 grid points inside thermal boundary layers and 64 grid points inside viscous boundary layers, which guarantees to resolve the boundary layer adequately. Besides, we have conducted a set of grid independence studies and checked the spatial resolution in bulk region to resolve all relevant scales (Kolmogorov scale and Batchelor scale). The parameter space that was explored can be found in Fig. 1C. More details about the DNS are documented in the Supplementary Materials.

### Experiments

Experiments are performed in a cylindrical annulus with solid-body rotation as sketched in Fig. 1A. The two concentric cylinders are machined from a solid piece of copper to control the symmetry of the system, and their surfaces are electroplated with a thin layer of nickel to prevent oxidation. The inner cylinder has an outer radius of *R _{i}* = 120 mm, and the outer cylinder has an inner radius of

*R*= 240 mm, resulting in a gap of

_{o}*L*= 120 mm and a radius ratio of η =

*R*/

_{i}*R*= 0.5. The gap with a height

_{o}*H*= 120 mm is sandwiched by a top plexiglass plate and a bottom teflon plate, resulting in an aspect ratio of Γ =

*H*/

*L*= 1. Plexiglass with an excellent transparency is used as the material of the top lid, allowing us to visualize the flow field, whereas teflon with an excellent corrosion resistance and a good strength is machined as the bottom base. These end plates and cylinders are fixed together and leveled on a rotating aluminum frame with a rotation rate up to 546 rpm. Four silicone rubber film heaters are attached to the outside of the outer cylinder and are supplied by a DC power supply (Ametek, XG 300-5) with 0.005% long-term stability. We use water as our working fluid, which has

*Pr*≈ 4.3.

Inside the cold inner cylinder, 16 channels are machined for the coolant fluid to pass through, and the inner cylinder is regulated at constant temperature using a temperature-controlled circulating bath (PolyScience, AP45R-20-A12E). The coolant pipes and electric cables go through the hollow stainless steel shaft and connect with a slip ring and a rotary union (Moflon, MEPH200), which has 2 channels for liquids, 6 channels for powers (2220 W), and 48 channels for electrical signals. The shaft is driven by a toothed belt, and the pulley has a gear ratio of 60 : 32 ≈ 1.88 : 1. The whole system is powered by a servomotor (Yaskawa, SGM7G-1AA), which has a rated power of 11 kW.

For high-precision temperature and heat flux measurements, it is essential to minimize the heat leakage from the experimental apparatus to the surroundings. Various thermal shields that are regulated at appropriate temperatures are installed in the system to prevent heat losses. Furthermore, the whole system is placed in a big toughened glass box where the temperature is controlled by a proportional-integral-derivative (PID) controller to match the mean temperature of the bulk fluid. For more detailed descriptions on the experimental setup, we refer to the Supplementary Materials. The measurement techniques are documented in Methods.

## RESULTS

### Effects of the Coriolis force on the heat transport and flow structures

A series of numerical simulations are performed to explore the influence of the Coriolis force on the heat transfer and flow structures. In numerical experiments, *Ro*^{−1} varies from 10^{−2} to 4 × 10^{2} with several fixed Rayleigh numbers *Ra* = 10^{6}, 2.2 × 10^{6}, 4.7 × 10^{6}, 10^{7}, 2.2 × 10^{7}, 4.7 × 10^{7}, 10^{8}, 2.2 × 10^{8}, and 4.7 × 10^{8}. As shown in Fig. 2A, for these *Ra* studied, the influence of Coriolis force can be divided into three regimes. In regime I (*Nu* does not depend on *Ro*^{−1}. In regime III (*36*, *37*), resulting in a quasi–two-dimensional flow state and thus a reduction of *Nu* compared with regime I. For regime II (*Ra* and *Ro*^{−1} with rich flow states. In this regime, the general trend is that *Nu* decreases with *Ro*^{−1}. The instantaneous temperature fields for several *Ro*^{−1} at *Ra* = 10^{8} (see Fig. 2, B, C, E, and F) from DNS support the explanation of this *Ro* dependence on the heat transport. As illustrated in Fig. 2B, the flow is three-dimensional at *Ro*^{−1} = 0.1 because the Coriolis force is too small to influence the flow effectively. However, with Ro^{−1} increasing, the enhanced Coriolis force tends to suppress vertical variation of the convection flow, and the flow gradually becomes two-dimensional, which is a manifestation of the Taylor-Proudman theorem. The two-dimensionalization of the flow field should be responsible for the reduction of heat transport at *38*). As is evident in Fig. 2F, the flow is nearly two-dimensional at *Ro*^{−1} saturates the influence of the Coriolis force, which explains the nearly constant *Nu* when *u _{z}*)

_{rms}〉

_{r, φ, z}to measure the influence of the Taylor-Proudman theorem, as shown in Fig. 2D. It shows that 〈(

*u*)

_{z}_{rms}〉

_{r, φ, z}has nearly a constant value when

*Ro*

^{−1}, and finally approaches to around zero after

*u*)

_{z}_{rms}〉

_{r, φ, z}versus Ro

^{−1}is consistent with the dependence of

*Nu*on

*Ro*

^{−1}in general.

We have noticed that the critical *Ra*. As illustrated in Fig. 2D, for *Ra* = 10^{7}, we have *Ra* = 10^{8}, we have *u _{z}*)

_{rms}〉

_{r, φ, z}, and

*Ra*cases are available in the Supplementary Materials.)

### Global turbulent heat transport

As shown in the section above, *Nu* does not depend much on *Ro*^{−1} for high *Ro*^{−1} regime, in which we will study how heat transfer depends on *Ra* in this section. We have performed 48 experiments and more than 130 numerical simulations. Figure 1C shows the parameter space of the experimental and numerical studies. In the experimental studies, we note that the existence of Earth’s gravity and lids is unavoidable. Several numerical simulations have been performed to study the influences of Earth’s gravity and lids, which show that their effects on *Nu* are small. In addition, we note that the centrifugal force increases linearly with radial location. To study the influence of nonuniform driving force, we have performed two sets of simulations with the radial-dependent centrifugal acceleration ω^{2}*r* and with a constant acceleration ω^{2}(*R _{o}* +

*R*)/2. The results show that the effect of the radial-dependent gravity has a similar role as the non–Oberbeck-Boussinesq effect (

_{i}*39*–

*41*), which does not have much influence on the heat transport and flow structures in the current parameter regime (see the Supplementary Materials for details).

Figure 3A shows the measured *Nu* as a function of *Ra* for different *Ro*. Each measurement lasts at least 4 hours after the system has reached a thermally steady situation to ensure a statistically stationary state (for the detailed measurement procedure, see the Supplementary Materials). As being evident in Fig. 3A, the experiments and numerical simulations are in excellent agreement. Combining experiments and simulations, this study covers more than four decades of *Ra*, i.e., from 10^{6} to 6.5 × 10^{10}. The range of *Ro*^{−1} is from 18 to 58 in experiments, and the data series show a consistent dependence of *Nu* on *Ra*. For comparison, we also plot the data calculated from Grossmann-Lohse (G-L) theory (*39*) in the corresponding *Ra* range. To better reveal the local exponent, we plot the compensated data in Fig. 3B. The experimental and numerical data have a lower amplitude as compared to the classical RBC (G-L line) due to the different flow geometry. However, the scaling dependence of *Nu* versus *Ra* shows a good agreement between the data and G-L theory at *Ra* ⪅ 10^{10} with a scaling exponent γ = 0.27 ± 0.01, which is close to the typical value found in two-dimensional RBC.

It is very unexpected that the local effective exponent γ of *Nu* ∝ *Ra*^{γ} even exceeds one-third at *Ra* ≳ 4 × 10^{10} in experiments. Is this scaling regime connected to the appearance of ultimate regime of Rayleigh-Bénard turbulence (*42*–*44*)? If it is not in the ultimate regime, what causes this steep scaling exponent? As our current *Ra* is much lower than the transition *Ra* for ultimate turbulence observed in the classical RBC, we cannot provide a concrete answer on this question. The possible answers are that (i) because of the different way of the driving and the geometry, the transition *Ra* to the ultimate regime at the current system is lower than that in the classical RBC. (ii) Another possibility is that the enhanced local slope is due to the emergence of new flow states. As demonstrated in Fig. 1C, with *Ra* increasing, the location of (*Ra*, *Ro*^{−1}) where the enhanced local effective scaling appears seems to move from regime III to regime II. Thus, the flow may change from a two-dimensional flow state to a three-dimensional flow state due to the strong thermal driving, resulting in a higher *Nu* in this *Ra* range. The two-dimensional simulations for high *Ra* coincide well with the experiments, indicating that within the *Ra* range explored in two-dimensional simulations, the flow is still quasi–two-dimensional. Unfortunately, we cannot push the simulations to the *Ra* regime with the effective exponent larger than one-third, in which the flow may not be in the quasi–two-dimensional state. We incline to anticipate that the changing from a two-dimensional flow state to a three-dimensional flow state may be responsible for the enhancement of local effective exponent at *Ra* ≳ 4 × 10^{10}. We notice that the scaling range with γ > 1/3 is very narrow in the current work; therefore, further work at even higher *Ra* is needed to clarify this issue.

### Zonal flow

We now study the dynamics of zonal flow experimentally and numerically. The “zonal flow” phenomenon has been investigated in experiments of geophysical and astrophysical flows (*45*, *46*). Figure 4B shows the time series of local temperature fluctuations in water at *Ra* = 6.6 × 10^{9}, *Pr* = 4.3, and *Ro*^{−1} = 18, which unexpectedly shows a noticeable periodicity. To prove that this novel phenomenon is connected to the azimuthal movement of the coherent flow, we perform flow visualization with aqueous glycerol solution as mentioned in Methods to demonstrate the flow pattern.

Figure 4A shows some typical streak images from flow visualization at three different times (for video, see movie S2), where there are four pairs of rolls parallel to the rotating axis. For reference, we highlight one of the convective rolls using a yellow ellipse. From Fig. 4A (I to III), unexpectedly, the convection rolls move clockwise around the center with a faster rotation rate than the background rotation of the experimental system, which we name zonal flow. The flow visualization further justifies that the revolution of the cold plume arms triggers the periodic temperature signals measured by the thermistors.

Next to the experiments, we also find this phenomenon in the simulations (see Fig. 4C and movie S3). We use a black ellipse to highlight a selected cold plume. As shown in Fig. 4C, on the reference of the (clockwise) rotating frame, the plume arms still evolve in a clockwise direction, indicating a net rotation of the coherent flow (zonal flow). This zonal flow is further quantified through the averaged (axially and azimuthally in space and over time) azimuthal velocity profile versus the radial distance from the inner cylinder wall *r*, as shown in Fig. 4D.

What causes this net rotation of the convection rolls? The dashed lines in Fig. 4C (I) mark the direction of hot/cold plumes without Coriolis force. But because of the effects of Coriolis force, the plumes deflect to their left when the system rotates clockwise. As is evident in Fig. 4C (I) where the plumes are just formed in the beginning stage, the deflected angle of hot plumes (∠*a*) is approximately equal to the deflected angle of cold plumes (∠*b*). However, because of the different curvature of the cylinders, the similar deflected angles of the hot and cold plumes induce the different effects. The hot plumes directly affect the position (*A*) where the cold plumes are ejected, resulting in the clockwise rotation of the cold plumes and flow, whereas the impact of the cold plumes does not directly affect the motion of the hot plumes due to a relatively large distance between the impact position of the cold plumes (*B*) and the ejecting position of the hot plumes. Thus, the hot plumes win and push the overall flow to move in the same direction of the background rotation. We also perform the experiments and simulations with the opposite direction of the background rotation, and the results are consistent. In addition, several numerical simulations have been performed to test the dependence of zonal flow on the radius ratio η, which shows that the zonal flows become weaker and weaker with η increasing from 0.4 to 0.9. The difference in the curvature of the two cylinders plays the key role for the net rotation of convection rolls. One may expect to observe different types of zonal flows at different η, such as two large-scale winds with the opposite directions near the plates (*46*). This remarkable net rotation of the coherent flow structures along the same direction of the background rotation could be connected to many relevant flow phenomena in nature, which deserves systematical studies in the future.

## SUMMARY

In summary, we have experimentally and numerically studied the global heat transport and turbulent flow structures in a rotating annulus with hot outer cylinder and cold inner cylinder, i.e., ACRBC. In experiments, the mean centrifugal acceleration covers the range from 5 to 60 times gravitational acceleration. We show that the effective scaling exponent of global heat transport transitions to γ = 1/3 at *Ra* ≈ 10^{10} and finally exceeds γ > 1/3 for *Ra* ⪆ 4 × 10^{10}. Unexpectedly, we observed the faster azimuthal motion of the coherent overall flow structure faster than the background rotation of the system, which we call zonal flow, and provided a physical understanding on this remarkable net rotation motion of the coherent structures. This novel experimental approach sheds new lights on how to efficiently extend the parameter regimes for the study of buoyancy-driven turbulent flows. Furthermore, the findings in the current system driven by the centrifugal acceleration can help us understand phenomena in geophysical and astrophysical flows.

## METHODS

### Measurement techniques

In our experiments, all the temperature measurements are based on negative temperature coefficient thermistors. We use *47*). Two types of thermistors are used: One with a head diameter of 2.5 mm (Omega, 44131) is used to measure the temperature of the inner and outer cylinders, and the other with a head diameter of 300 μm (Measurement Specialties, GAG22K7MCD419) is inserted into the convection cell to resolve fast temperature fluctuations, which is located 30 mm away from the cold cylinder and at mid-height in the vertical direction.

To visualize the flow field, we use nylon fibers with length *l* = 3 mm and diameter *d* = 0.5 mm. Nylon has a density of ρ* _{p}* ≈ 1.14 g/cm

^{3}that is a little heavier than water. In the flow visualization experiments, we use aqueous glycerol solution with 54% glycerol by volume to match the density of nylon. The properties of the aqueous glycerol solution are listed in the Supplementary Materials for reference. Four light-emitting diodes are used as light source, and a charge-coupled device camera (Ximea, MD028MU-SY) is put on the top of the cell to record images. As the system rotates rapidly, it is difficult to make the camera rotate synchronously. So we keep the camera fixed and then set the frame rate of the camera the same as the rotational speed. Processing the images taken by camera through a MATLAB code, we can get streak images to visualize the flow field. More processing details can be found in the Supplementary Materials.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/40/eabb8676/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank Y. Yang, Y. Liu, Q. Zhou, and S. Liu for insightful discussions and G.-W. Bruggert for the technical assistance with the setup. We acknowledge the National Supercomputing Center in Guangzhou (NSCC-GZ) for the allocation of computing time.

**Funding:**This work was financially supported by the Natural Science Foundation of China under grant nos. 11988102, 91852202, 11861131005, and 11672156; the Tsinghua University Initiative Scientific Research Program (20193080058); and the Netherlands Organisation for Scientific Research (NWO) through the Multiscale Catalytic Energy Conversion (MCEC) research center. X.Z. acknowledges the financial support from the Center for Mathematical Sciences and Applications at the Harvard University.

**Author contributions:**C.S. conceived the project. H.J. and D.W. conducted the experiments. X.Z. and H.J. performed the simulations. H.J., S.G.H., and C.S. designed the setup. All authors discussed the physics and contributed to the writing of the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).