# [

###### Abstract

Results from direct numerical simulation for three-dimensional Rayleigh-Bénard convection in samples of aspect ratio and up to Rayleigh number are presented. The broad range of Prandtl numbers is considered. In contrast to some experiments, we do not see any increase in , neither due to number effects, nor due to a constant heat flux boundary condition at the bottom plate instead of constant temperature boundary conditions. Even at these very high , both the thermal and kinetic boundary layer thicknesses obey Prandtl-Blasius scaling.

and number dependence of heat transport in high number RB] Prandtl and Rayleigh number dependence of heat transport in high Rayleigh number thermal convection Richard J.A.M. Stevens, Detlef Lohse, and Roberto Verzicco] Richard J.A.M. Stevens, Detlef Lohse, and Roberto Verzicco

## 1 Introduction

In Rayleigh-Bénard (RB) convection fluid in a box is heated from below and cooled from above (Ahlers et al. (2009c)). This system is paradigmatic for turbulent heat transfer, with many applications in atmospheric and environmental physics, astrophysics, and process technology. Its dynamics is characterized by the Rayleigh number and the Prandtl number . Here, is the height of the sample, is the thermal expansion coefficient, the gravitational acceleration, the temperature difference between the bottom and the top of the sample, and and the kinematic viscosity and the thermal diffusivity, respectively. Almost all experimental and numerical results on the heat transfer, indicated by the Nusselt number , agree up to (see the review of Ahlers et al. (2009c) for detailed references) and are in agreement with the description of the Grossmann-Lohse (GL) theory (Grossmann & Lohse (2000, 2001, 2002, 2004)). However, for higher the situation is less clear.

Most experiments for are performed in samples with aspect ratios and , where and are the diameter and height of the sample, respectively. The majority of these experiments are performed with liquid helium near its critical point (Chavanne et al. (2001); Niemela et al. (2000, 2001); Niemela & Sreenivasan (2006); Roche et al. (2001, 2002, 2010)). While Niemela et al. (2000, 2001) and Niemela & Sreenivasan (2006) found a increase with , the experiments by Chavanne et al. (2001) and Roche et al. (2001, 2002, 2010) gave a steep increase with . In these helium experiments the number increases with increasing . Funfschilling et al. (2009) and Ahlers et al. (2009a, b) performed measurements around room temperature with high pressurized gases with nearly constant and do not find such a steep increase. Niemela & Sreenivasan (2010) found two branches in a sample. The high number branch is higher than the low number branch. By necessity, increases more steeply in the transition region. The scaling in the transition region happens to be around . There are thus considerable differences in the heat transfer obtained in these different experiments in the high number regime. Very recently, Ahlers and coworkers (see Ahlers (2010) and addendum to Ahlers et al. (2009b)) even found two different branches in one experiments with the steepest branch going as .

There is no clear explanation for this disagreement although it has been conjectured that variations of the number, the use of constant temperature or a constant heat flux condition at the bottom plate, the finite conductivity of the horizontal plates and side wall, non Oberbeck-Boussinesq effects, i.e. the dependence of the fluid properties on the temperature, the existence of multiple turbulent states (Grossmann & Lohse (2011)), and even wall roughness and temperature conditions outside the sample might play a role. Since the above differences among experiments might be induced by unavoidable technicalities in the laboratory set–ups, within this context, direct numerical simulations are the only possibility to obtain neat reference data that strictly adhere to the intended theoretical problem and that could be used as guidelines to interpret the experiments: This is the main motivation for the present study.

## 2 Numerical procedure

We start this paper with a description of the numerical procedure that is used to investigate the influence on the heat transfer for two of the issues mentioned above. First, we discuss the effect of the number on the heat transport in the high number regime. Subsequently, we will discuss the difference between simulations performed with a constant temperature at the bottom plate with simulations with a constant heat flux at the bottom plate. We take the constant heat flux condition only at the bottom plate, because in real setups the bottom plate is in contact with a heater while the top plate is connected to a thermostatic bath. Thus, the condition of constant heat flux applies at most only to the bottom plate (Niemela et al. (2000, 2001) and Niemela & Sreenivasan (2006)). At the top plate constant temperature boundary conditions are assumed to strictly hold, i.e. perfect heat transfer to the recirculating cooling liquid. We will conclude the paper with a brief summary, discussion, and outlook to future simulations.

The flow is solved by numerically integrating the three-dimensional Navier-Stokes equations within the Boussinesq approximation. The numerical method is already described in Verzicco & Camussi (1997, 2003) and Verzicco & Sreenivasan (2008) and here it is sufficient to mention that the equations in cylindrical coordinates are solved by central second-order accurate finite-difference approximations in space and time. We performed simulations with constant temperature conditions at the bottom plate for and in an aspect ratio sample. We also present results for a simulation at at in a sample. In addition, we performed simulation with a constant heat flux at the bottom plate and a constant temperature at the top plate, see Verzicco & Sreenivasan (2008) for details, for , , and . Because in all simulations the temperature boundary conditions are precisely assigned, the surfaces are infinitely smooth, and the Boussinesq approximation is unconditionally valid, the simulations provide a clear reference case for present and future experiments.

In Stevens et al. (2010b) we investigated the resolution criteria that should be satisfied in a fully resolved DNS simulation and Shishkina et al. (2010) determined the minimal number of nodes that should be placed inside the boundary layers. The resolutions used here are based on this experience and we stress that in this study we used even better spatial resolution than we used in Stevens et al. (2010b) to be sure that the flow is fully resolved.

To give the reader some idea of the scale of this study we mention the resolutions that were used in the most demanding simulations, i.e., simulations that take at least 100.000 DEISA CPU hours each. For this are the simulations at , which are performed on either a (azimuthal, radial, and axial number of nodes) grid for and or on a grid for and . The simulations at are performed on a grid. The simulations for were run for at least dimensionless time units (defined as ), while these simulations at cover about time units. The simulation with a constant heat flux condition at the bottom plate and constant temperature condition at the top plate in a sample with at is performed on a grid. The simulation at with in the sample has been performed on a grid, which makes it the largest fully bounded turbulent flow simulation ever. This simulation takes about vectorial CPU hours on HLRS (equivalent to DEISA CPU hours). To store one snapshot of the field (, ,, because follows from continuity) costs 160 GB in binary format. A snapshot of this flow is shown in figure 1. Movies of this simulation are included in the supplementary material.

## 3 Numerical results on

In figure 2a and figure 3 the DNS results for in the and samples are compared with experimental data. The result for in the sample agrees well with the experimental data of Niemela et al. (2000, 2001), Niemela & Sreenivasan (2006), Funfschilling et al. (2009), and Ahlers et al. (2009a, b), while there is a visible difference with the results of Chavanne et al. (2001). A comparison of the results for with the experimental data of Roche et al. (2010) shows that there is a good agreement for higher numbers, while for lower we obtain slightly larger than in those experiments. We again stress that we performed resolution checks for this case (up to ), and in addition considering the good agreement with the results for , we exclude that our DNS results overestimate .

Figure 2b and 3 show that in some experiments the number increases with increasing . This difference in is often mentioned as one of the possible causes for the observed differences in the heat transfer between the experiments. Figure 4 shows the number as function of for different . This figure shows that the effect of the number on the heat transfer decreases with increasing . This means that the differences in the heat transport that are observed between the experiments for , see figure 2a and 3a, are not a number effect. This is in agreement with the theoretical prediction of the GL-model for , which is shown in figure 4.

## 4 Scaling of thermal and kinetic boundary layers

We determined the thermal and kinetic BL thickness for the simulations in the sample. The horizontally averaged thermal BL thickness () is determined from , where is the intersection point between the linear extrapolation of the temperature gradient at the plate with the behavior found in the bulk (Stevens et al. (2010a)). In figure 5a it is shown that the scaling of the thermal BL thickness is consistent with the number measurements when the horizontal average is taken over the entire plate.

The horizontally averaged kinetic BL thickness () is determined from , where is based on the position where the quantity reaches its maximum. We use this quantity as Stevens et al. (2010b) and Shishkina et al. (2010) showed that it represents the kinetic BL thickness better than other available methods. Stevens et al. (2010a, b) also explained that this quantity cannot be used close to the sidewall as here it misrepresents the kinetic BL thickness. For numerical reasons, it can neither be calculated accurately in the center region. To be on the safe side we horizontally average the kinetic BL between , where indicates the diameter of the cell and all given values refer to that used.

Figure 5b reveals that for the kinetic and thermal BL thickness have the same scaling and thickness over a wide range of . In figure 6 the ratio is compared with results of the Prandtl-Blasius BL theory. We find a constant difference of about between the numerical results and the theoretical PB type prediction, see e.g. Shishkina et al. (2010). We emphasize that the deviation of the prefactor of only 15% is remarkably small, given that the PB boundary layer theory has been developed for parallel flow over an infinite flat plate, whereas here in the aspect ratio cell one can hardly find such regions of parallel flows at the top and bottom plates. Nonetheless, the scaling and even the ratio of the kinetic and thermal boundary layer thicknesses for these large numbers is well described by Prandtl-Blasius BL theory. This result agrees with the experimental results of Qiu & Xia (1998) and Sun et al. (2008). Indeed, Qiu & Xia (1998) showed that the kinetic BL near the sidewall obeys the scaling law of the Prandtl-Blasius laminar BL and Sun et al. (2008) showed the same for the boundary layers near the bottom plate. Recently, Zhou & Xia (2010); Zhou et al. (2010) have developed a method of expressing velocity profiles in the time-dependent BL frame and found that not only the scaling obeys the PB expectation, but even the rescaled velocity and temperature profiles From all this we can exclude that at the BL is turbulent.

## 5 Constant temperature versus constant heat flux condition at the bottom plate

It has also been argued that the different boundary conditions at the bottom plate, i.e. that some experiments are closer to a constant temperature boundary condition, and some are closer to a constant heat flux boundary condition, might explain the differences in the heat transport that are observed in the high number regime. Figure 7a compares the number in the simulations with constant temperature and constant heat flux at the bottom plate. The figure shows that the difference between these both cases is small and even decreases with increasing . For large no difference at all is seen within the (statistical) error bars, which however increase due to the shorter averaging time (in terms of large eddy turnovers) at the very large .

Figure 7b shows the time-averaged temperature of the bottom plate in the simulations with constant heat flux at the bottom plate for different . The radial dependence at the lower numbers can be understood from the flow structure in the sample: Due to the large scale circulation the fluid velocities are largest in the middle of the sample. Thus in the middle more heat can be extracted from the plate than close to the sidewall where the fluid velocities are smaller. It is the lack of any wind in the corners of the sample that causes the relative high time-averaged plate temperature there. The figure also shows that this effect decreases with increasing . The reason for this is that the turbulence becomes stronger at higher and this leads to smaller flow structures. Therefore the region close to the sidewall with relative small fluid velocities decreases with increasing and this leads to a more uniform plate temperature at higher . This effect explains that the simulations with constant temperature and constant heat flux condition at the bottom plate become more similar with increasing .

The small differences between the simulations with constant temperature and constant heat flux at the bottom plate shows that the differences between the experiments in the high number regime can not be explained by different plate conductivity properties. This finding is in agreement with the results of Johnston & Doering (2009). In their periodic two-dimensional RB simulations the heat transfer for simulations with constant temperature and constant heat flux (both at the bottom and the top plate) becomes equal at . For the three-dimensional simulations the heat transfer for both cases also becomes equal, but at higher . This is due to the geometrical effect discussed before, see figure 7b, that cannot occur in periodic two-dimensional simulations (Johnston & Doering (2009)).

## 6 Conclusions

In summary, we presented results from three-dimensional DNS simulations for RB convection in cylindrical samples of aspect ratios and up to and a broad range of numbers. The simulation at with in an aspect ratio sample is in good agreement with the experimental results of Niemela et al. (2000, 2001), Niemela & Sreenivasan (2006), Funfschilling et al. (2009), and Ahlers et al. (2009a, b), while there is a visible difference with the results of Chavanne et al. (2001). In addition, we showed that the differences in the heat transfer observed between experiments for can neither be explained by number effects, nor by the assumption of constant heat flux conditions at the bottom plate instead of constant temperature conditions. Furthermore, we demonstrated that the scaling of the kinetic and thermal BL thicknesses in this high number regime is well described by the Prandtl-Blasius theory.

Several questions remain: Which effect is responsible for the observed difference in vs scaling in the various experiments? Are there perhaps different turbulent states in the highly turbulent regime as has been suggested for RB flow by Grossmann & Lohse (2011), but also for other turbulent flows in closed systems by Cortet et al. (2010)? At what number do the BLs become turbulent? As in DNSs both the velocity and temperature fields are known in the whole domain (including in the boundary layers where the transition between the states is suggested to take place), they will play a leading role in answering these questions.

Acknowledgement: We thank J. Niemela, K.R. Sreenivasan, G. Ahlers, and P. Roche for providing the experimental data. The simulations were performed on Huygens (DEISA project), CASPUR, and HLRS. We gratefully acknowledge the support of Wim Rijks (SARA) and we thank the DEISA Consortium (www.deisa.eu), co-funded through the EU FP7 project RI-222919, for support within the DEISA Extreme Computing Initiative. RJAMS was financially supported by the Foundation for Fundamental Research on Matter (FOM). We thank the Kavli institute, where part of the work has been done, for its hospitality. This research wasÊ supported in part by the National Science Foundation under Grant No.ÊPHY05-51164, via the Kavli Institute of Theoretical Physics.

- Ahlers (2010) Ahlers, G. January 2010, lecture at the euromech colloquium in les houches, see www.hirac4.cnrs.fr/hirac4-talksfiles/ahlers.pdf.
- Ahlers et al. (2009a) Ahlers, G., Bodenschatz, E., Funfschilling, D. & Hogg, J. 2009a Turbulent Rayleigh-Bénard convection for a Prandtl number of 0.67. J. Fluid. Mech. 641, 157–167.
- Ahlers et al. (2009b) Ahlers, G., Funfschilling, D. & Bodenschatz, E. 2009b Transitions in heat transport by turbulent convection at Rayleigh numbers up to . New J. Phys. 11, 123001.
- Ahlers et al. (2009c) Ahlers, G., Grossmann, S. & Lohse, D. 2009c Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection. Rev. Mod. Phys. 81, 503.
- Ahlers & Xu (2001) Ahlers, G. & Xu, X. 2001 Prandtl-number dependence of heat transport in turbulent Rayleigh-Bénard convection. Phys. Rev. Lett. 86, 3320–3323.
- Chavanne et al. (2001) Chavanne, X., Chilla, F., Chabaud, B., Castaing, B. & Hebral, B. 2001 Turbulent Rayleigh-Bénard convection in gaseous and liquid he. Phys. Fluids 13, 1300–1320.
- Cortet et al. (2010) Cortet, P., Chiffaudel, A., Daviaud, F. & Dubrulle, B. 2010 Experimental evidence of a phase transition in a closed turbulent flow. Phys. Rev. Lett. 105, 214501.
- Funfschilling et al. (2009) Funfschilling, D., Bodenschatz, E. & Ahlers, G. 2009 Search for the ”ultimate state” in turbulent Rayleigh-Bénard convection. Phys. Rev. Lett. 103, 014503.
- Grossmann & Lohse (2000) Grossmann, S. & Lohse, D. 2000 Scaling in thermal convection: A unifying view. J. Fluid. Mech. 407, 27–56.
- Grossmann & Lohse (2001) Grossmann, S. & Lohse, D. 2001 Thermal convection for large Prandtl number. Phys. Rev. Lett. 86, 3316–3319.
- Grossmann & Lohse (2002) Grossmann, S. & Lohse, D. 2002 Prandtl and Rayleigh number dependence of the Reynolds number in turbulent thermal convection. Phys. Rev. E 66, 016305.
- Grossmann & Lohse (2004) Grossmann, S. & Lohse, D. 2004 Fluctuations in turbulent Rayleigh-Bénard convection: The role of plumes. Phys. Fluids 16, 4462–4472.
- Grossmann & Lohse (2011) Grossmann, S. & Lohse, D. 2011 Multiple scaling in the ultimate regime of thermal convection. Phys. Fluids, in press .
- Johnston & Doering (2009) Johnston, H. & Doering, C. R. 2009 Comparison of turbulent thermal convection between conditions of constant temperature and constant flux. Phys. Rev. Lett. 102, 064501.
- Niemela et al. (2000) Niemela, J., Skrbek, L., Sreenivasan, K. R. & Donnelly, R. 2000 Turbulent convection at very high Rayleigh numbers. Nature 404, 837–840.
- Niemela et al. (2001) Niemela, J., Skrbek, L., Sreenivasan, K. R. & Donnelly, R. J. 2001 The wind in confined thermal turbulence. J. Fluid Mech. 449, 169–178.
- Niemela & Sreenivasan (2010) Niemela, J.J. & Sreenivasan, K.R. 2010 Does confined turbulent convection ever attain the ’asymptotic scaling’ with power? New J. Phys. 12, 115002.
- Niemela & Sreenivasan (2006) Niemela, J. & Sreenivasan, K. R. 2006 Turbulent convection at high Rayleigh numbers and aspect ratio 4. J. Fluid Mech. 557, 411 – 422.
- Qiu & Xia (1998) Qiu, X. L. & Xia, K.-Q. 1998 Viscous boundary layers at the sidewall of a convection cell. Phys. Rev. E 58, 486–491.
- Roche et al. (2001) Roche, P. E., Castaing, B., Chabaud, B. & Hebral, B. 2001 Observation of the 1/2 power law in Rayleigh-Bénard convection. Phys. Rev. E 63, 045303.
- Roche et al. (2002) Roche, P. E., Castaing, B., Chabaud, B. & Hebral, B. 2002 Prandtl and Rayleigh numbers dependences in Rayleigh-Bénard convection. Europhys. Lett. 58, 693–698.
- Roche et al. (2010) Roche, P.-E., Gauthier, F., Kaiser, R. & Salort, J. 2010 On the triggering of the ultimate regime of convection. New J. Phys. 12, 085014.
- Shishkina et al. (2010) Shishkina, O., Stevens, R. J. A. M., Grossmann, S. & Lohse, D. 2010 Boundary layer structure in turbulent thermal convection and its consequences for the required numerical resolution. New J. Phys. 12, 075022.
- Stevens et al. (2010a) Stevens, R. J. A. M., Clercx, H. J. H. & Lohse, D. 2010a Boundary layers in rotating weakly turbulent Rayleigh-Bénard convection. Phys. Fluids 22, 085103.
- Stevens et al. (2010b) Stevens, R. J. A. M., Verzicco, R. & Lohse, D. 2010b Radial boundary layer structure and Nusselt number in Rayleigh-Bénard convection. J. Fluid. Mech. 643, 495–507.
- Sun et al. (2008) Sun, C., Cheung, Y. H. & Xia, K. Q. 2008 Experimental studies of the viscous boundary layer properties in turbulent Rayleigh-Bénard convection. J. Fluid Mech. 605, 79 – 113.
- Verzicco & Camussi (1997) Verzicco, R. & Camussi, R. 1997 Transitional regimes of low-prandtl thermal convection in a cylindrical cell. Phys. Fluids 9, 1287–1295.
- Verzicco & Camussi (2003) Verzicco, R. & Camussi, R. 2003 Numerical experiments on strongly turbulent thermal convection in a slender cylindrical cell. J. Fluid Mech. 477, 19–49.
- Verzicco & Sreenivasan (2008) Verzicco, R. & Sreenivasan, K. R. 2008 A comparison of turbulent thermal convection between conditions of constant temperature and constant heat flux. J. Fluid Mech. 595, 203–219.
- Xia et al. (2002) Xia, K.-Q., Lam, S. & Zhou, S. Q. 2002 Heat-flux measurement in high-Prandtl-number turbulent Rayleigh-Bénard convection. Phys. Rev. Lett. 88, 064501.
- Zhou et al. (2010) Zhou, Q., Stevens, R. J. A. M., Sugiyama, K., Grossmann, S., Lohse, D. & Xia, K.-Q. 2010 Prandtl-Blasius temperature and velocity boundary layer profiles in turbulent Rayleigh-Bénard convection. J. Fluid. Mech. 664, 297Ð312.
- Zhou & Xia (2010) Zhou, Q. & Xia, K.-Q. 2010 Measured instantaneous viscous boundary layer in turbulent Rayleigh-Bénard convection. Phys. Rev. Lett. 104, 104301.