# Transition to the ultimate regime in two-dimensional Rayleigh-Bénard convection

###### Abstract

The possible transition to the so-called ultimate regime, wherein both the bulk and the boundary layers are turbulent, has been an outstanding issue in thermal convection, since the seminal work by Kraichnan [Phys. Fluids 5, 1374 (1962)]. Yet, when this transition takes place and how the local flow induces it is not fully understood. Here, by performing two-dimensional simulations of Rayleigh-Bénard turbulence covering six decades in Rayleigh number Ra up to for Prandtl number Pr , for the first time in numerical simulations we find the transition to the ultimate regime, namely at . We reveal how the emission of thermal plumes enhances the global heat transport, leading to a steeper increase of the Nusselt number than the classical Malkus scaling [Proc. R. Soc. London A 225, 196 (1954)]. Beyond the transition, the mean velocity profiles are logarithmic throughout, indicating turbulent boundary layers. In contrast, the temperature profiles are only locally logarithmic, namely within the regions where plumes are emitted, and where the local Nusselt number has an effective scaling , corresponding to the effective scaling in the ultimate regime.

Rayleigh-Bénard (RB) flow, in which the fluid is heated from below and cooled from above, is a paradigmatic representation of thermal convection, with many features that are of interest in natural and engineering applications. Examples range from astrophysical and geophysical flows to industrial thermal flows Ahlers et al. (2009); Lohse and Xia (2010); Chilla and Schumacher (2012). When the temperature difference between the two plates (expressed in dimensionless form as the Rayleigh number Ra) is large enough, the system is expected to undergo a transition from the so-called “classical regime” of turbulence, where the boundary layers (BLs) are of laminar typeSun et al. (2008); Zhou and Xia (2010); Zhou et al. (2010); du Puits and Willert (2016), to the so-called “ultimate regime”, where the BLs are of turbulent type, as first predicted by Kraichnan Kraichnan (1962) and later by others Spiegel (1971); Grossmann and Lohse (2000, 2001, 2011, 2012). In the classical regime, the Nusselt number Nu (dimensionless heat transfer) is known to effectively scale as Ra, with the effective scaling exponent Grossmann and Lohse (2000, 2001); Stevens et al. (2013); Malkus (1954); Priestley (1954). Beyond the transition to the ultimate regime, the heat transport is expected to increase substantially, reflected in an effective scaling exponent Kraichnan (1962); Ahlers et al. (2009); Grossmann and Lohse (2011).

Hitherto, the evidence for the transition to the ultimate regime has only come from experimental measurements. In fact, the community is debating at what Ra the transition starts and even whether there is a transition at all. For example, Niemela & Sreenivasan Niemela and Sreenivasan (2010) observed that first increases above 1/3 around and then decreases back to 1/3 again for . Subsequently, Urban et al. Urban et al. (2014) also reported for . However, Chavanne et al. Chavanne et al. (1997, 2001) found that the effective scaling exponent increases to 0.38 for . In the experiments mentioned above, low temperature helium was used as the working fluid and Prandtl number Pr changes with increasing Ra. In contrast to helium, SF has roughly pressure independent Pr. This allows He et al. He et al. (2012a, b) to achieve the ultimate regime more conclusively. They observed a similar exponent 0.38, but this exponent was found only to start at a much higher Ra (the transition starts at Ra ). This observation is compatible with the theoretical prediction Grossmann and Lohse (2000, 2001) for the onset the ultimate regime. It is also consistent with the theoretical prediction of Refs. Kraichnan (1962); Grossmann and Lohse (2011), according to which logarithmic temperature and velocity BLs are necessary to obtain an effective scaling exponent for that Ra.

The apparent discrepancies among various high Ra RB experiments have been attributed to many factors. The change of Pr, the non-Boussinesq effect, the use of constant temperature or constant heat flux condition, the finite conductivity of the plates, and the sidewall effect can all play different roles Ahlers et al. (2009); Stevens et al. (2011). Direct numerical simulations (DNS), which do not have these unavoidable artefacts as occurring in experiments, can ideally help to understand the transition to the ultimate regime, with the strict accordance to the intended theoretic RB formulations. Unfortunately, high Ra simulations in three dimensions (3D) are prohibitively expensive Shishkina et al. (2010); Stevens et al. (2010). The highest Rayleigh number achieved in 3D RB simulations is Stevens et al. (2011), which is one order of magnitude short of the expected transitional Ra. Two-dimensional (2D) RB simulations, though different from 3D ones in terms of integral quantities for small Pr Schmalzl et al. (2004); van der Poel et al. (2013), still capture the many essential features of 3D RB van der Poel et al. (2013). Consequently, in recent years, 2D DNS has been widely used to test theories, not only for normal RB Huang and Zhou (2013); Zhang et al. (2017), but also for RB in porous media Hewitt et al. (2012). Although also expensive at high Ra, now we have the chance to push forward to using 2D simulations as we will show in this manuscript.

Another advantage of DNS as compared to experiment is that velocity and temperature profiles can be easily measured, to check whether they are logarithmic in the ultimate regime, as expected from theory. Specifically, for the temperature, only a few local experimental measurements were available in the near-sidewall regions of RB cells, which showed logarithmic profiles Ahlers et al. (2012, 2014). Even worse, for velocity, there is almost no evidence for the existence of a logarithmic BL, due to the experimental challenges. For instance, in cylindrical cells with aspect ratio , the mean velocity profile cannot be easily quantified because of the absence of a stable mean roll structure [24]. In situations where stable rolls do exist (e.g. narrow rectangular cells), the highest Ra available are still far below the critical Ra at which logarithmic velocity BLs can manifest themselves Sun et al. (2008); du Puits and Willert (2016).

As numerical simulations provide us with every detail of the flow field which might be unavailable in experiments, they also enable us to reveal the links between the global heat transport and the local flow structures. A few attempts (both 2D and 3D) have been made in the classical regime, in which logarithmic temperature BLs were detected, by selectively sampling the regions where the plumes are ejected to the bulk Ahlers et al. (2012); van der Poel et al. (2015a). However, it is still unclear how these local logarithmic BLs contribute to the attainment of the global heat transport enhancement during the transition to the ultimate regime.

In this work, for the first time in DNS we do observe a transition to the ultimate regime in 2D, namely at , similar as in the 3D RB experiments of Ref. He et al. (2012a). DNS also provides first evidence that the mean velocity profiles follow the log-law of the wall, in analogy to other paradigmatic turbulent flows, e.g. pipe, channel, and boundary flows Perry and Chong (1982); Marusic et al. (2010); Smits et al. (2011). Further, we explore the link between the local and global quantities to reveal the mechanism leading to the increased scaling exponent beyond the transition.

The simulations have been carried out using a well validated second-order finite-difference code Verzicco and Orlandi (1996); van der Poel et al. (2015b). The two control parameters are and , with being the thermal expansion coefficient, the gravitational acceleration, the temperature difference across a fluid layer of depth , the kinematic viscosity, and the thermal diffusivity. In the simulations, Pr were fixed at 1 and aspect ratio were fixed at 2, where is the width of the domain. With this , it has been found that the heat flux approximates the heat flux at infinite aspect ratio Johnston and Doering (2009). The boundary conditions were non-slip for velocity, constant temperature for the bottom and top plates, and periodic horizontally. Nu was calculated from the relation , with being the vertical velocity, the temperature, and the average over a horizontal plane and time. All the cases were well resolved. At the highest , we used a grid with mesh points. For details of the simulations, we refer to Ref. See Supplemental Material at [URL will be inserted by publisher] for numerical details (2017).

We begin by looking at the heat transport Nu as a function of Ra. In Fig. 1, we show Nu(Ra) compensated with Ra, for the range Ra=[]. Up to Ra = 10 (blue symbol), the effective scaling is essentially the same () as has been already observed Johnston and Doering (2009); van der Poel et al. (2013, 2014) in the classical regime where the BLs are laminar Zhou and Xia (2010); Zhou et al. (2010). This trend continues up to the transitional Rayleigh number Ra (green symbol). Beyond this, we witness the start of the transition to the ultimate regime, with a notably larger effective scaling exponent , as evident from the plateau in the compensated plot.

Next, to appreciate how the flow structures are different before and beyond the transition (Ra), we show the respective instantaneous temperature fields (see Fig. 2). The top panel presents a relatively low Ra (below Ra), while the middle panel shows a high Ra (beyond Ra). At low Ra, intense large scale rolls (LSR) are clearly visible. In comparison, at high Ra, the LSR, although still evident, contains much weaker and smaller structures. Interestingly, even at the highest Ra, the temperature field still has both plume ejecting and impacting regions. Additionally, these observations indicate that the spatial extent of plume ejecting regions do not grow in spite of the increase in Ra.

We now focus on the mean temperature and velocity fields at the transitional Ra. Remarkably, even after 500 dimensionless time units, the flow domain still shows a stable mean roll structure, i.e. the rolls are pinned with clearly demarcated plume ejecting and impacting regions (see Fig. 2(c)). The mean temperature and velocity fields display horizontal symmetry, which enables us to average them over a single LSR instead of the whole domain (as the velocity averaged horizontally for the whole domain will be zero).

Figure 3(a) shows the temporally and spatially averaged velocity profiles, performed on one single LSR. We plot the profiles in dimensionless wall units, in terms of and , where and . Here is the friction velocity Pope (2000). Similar to channel, pipe, and boundary layer flows, we can identify two distinct layers: a viscous sub-layer where , followed by a logarithmic region, where the velocity profile follows Pope (2000). The inverse slope gives , which is remarkably close to the Kármán constant in various 3D canonical wall-bounded turbulent flows Marusic et al. (2010); Smits et al. (2011). However, the parameter varies with Ra. With increasing Ra, the logarithmic range grows in spatial extent, until at Ra, it spans one decade in . Similarly, we express the temperature profile in wall units , where is the bottom plate temperature, and a characteristic temperature scale analogous to for the velocity Yaglom (1979). The mean temperature profile shows a similar viscous sub-layer , followed by a rather flat region, without a clear logarithmic dependence. Since the ultimate regime is associated with logarithmic profiles, the key question remains, as to why the mean temperature profile is not logarithmic despite the global scaling relations suggesting a transition in Fig. 1.

To find out, we look back more closely into the flow field of Fig. 2(c) where the mean flow was separated into (a) a plume ejection region, and (b) a plume impacting region. As noted earlier, the spatial extent of these regions does not grow with increasing Ra, and the mean flow field is horizontally symmetric. Therefore, the domain can be divided into plume ejection and impacting regions, enabling us to perform a conditional analysis for the temperature profiles specific to the respective regions. In Fig. 3(c), we plot these profiles separately, for different Ra. All the profiles collapse into a single curve in the viscous sub-layer. Beyond the viscous sub-layer, the impacting and ejecting regions show very different behavior. For the impacting regions, the temperature profile is flat (dotted curves), and remains so for all Ra. However, for the plume ejecting regions, we observe a clear log-layer (solid curves) with a profile , where is the equivalent Kármán constant for temperature, and varies with Ra. Similar to the velocity profiles, the extent of the log-layer increases with Ra. At the transitional Ra, it spans one decade in .

Temperature profiles that are locally logarithmic (in plume ejecting regions) have been observed before for both the classical and the ultimate regimes Ahlers et al. (2012, 2014); van der Poel et al. (2015a). Based on this, one hypothesis regarding how the system undergoes the transition to the ultimate regime is that the fraction of plume emitting regions (or hot spots) will gradually grow with increasing Ra van der Poel et al. (2015a). As speculated, the trend would continue until the entire BL becomes a hot spot, thus leading to a mean logarithmic temperature profile. Our findings indicate that here this is not the case, as even at plume impacting regions do not show a logarithmic temperature profile. The presence of these impacting regions makes the mean temperature profile also non-logarithmic (see Fig. 3(b)).

We now explain how the global heat transport scaling can still undergo a transition to the ultimate regime, though only the local temperature profile is logarithmic, not the globally averaged one. We recall that by definition on the plate surface, Nu . Following the observations from Fig. 2(c), we compute the local Nu on the plate surface from ejecting (Nu) and impacting (Nu) regions separately. These are shown in Fig. 4, compensated by Ra. Up to Ra, both Nu and Nu follow a similar trend, with their respective local scaling exponents and . However, beyond Ra, Nu and Nu diverge. The ejecting regions show an increased heat transport, with , which is precisely the ultimate scaling exponent predicted for with logarithmic BLs. In contrast, the impacting regions have a much lower scaling exponent . This means that the flow is partially in the ultimate regime and partially still in the classical regime. Based on these, we express the global Nusselt number, Nu = Nu + Nu, in analogy to the Grossmann-Lohse approach Grossmann and Lohse (2000, 2011) wherein the dissipation rate was separated into bulk and BL contributions. We write , where is expected to become even larger with increasing Ra Grossmann and Lohse (2011). The above expression asymptotically approaches the ultimate regime scaling when the plume ejecting regions become more and more dominant in transporting the heat with increasing Ra. Thus, with only the local temperature profile being logarithmic (in plume ejecting regions), the system can still undergo a gradual transition to the ultimate regime.

Finally, it is worthwhile to clarify the effect of the imposed two-dimensionality on the heat transfer. As mentioned in the beginning, 2D RB is different from 3D RB. However, the effective scaling exponents observed are identical in 2D and 3D for a wide range of Ra in the classical regime van der Poel et al. (2013), and here we found both in 2D and 3D the transition starts at Ra. Furthermore, the logarithmic BLs are theoretically expected for both 2D and 3D, as the theoretical argument Grossmann and Lohse (2012) is built on the Prandtl equations which are 2D. Also in other 2D canonical flows, logarithmic BLs have been observed, e.g. in channel flow Samanta et al. (2014); L’vov et al. (2009). Therefore, the physical insights gained from this work are useful for understanding the transition to ultimate turbulence in both 2D and 3D flows.

In conclusion, we have used two-dimensional simulations of Rayleigh-Bénard convection to investigate the transition to the ultimate regime of thermal convection. We followed the approach of using the local flow structures to explain the globally observed heat transfer enhancement. A transitional Rayleigh number Ra was found for the 2D RB with Pr , beyond which the mean velocity profile has a log-layer spanning one decade. However, the temperature profile is logarithmic only within the regions where plumes are ejected. The local effective Nusselt scaling exponent increases to 0.38 in the plume ejecting regions, corresponding to the ultimate regime. The transition to the ultimate regime can be understood as the gradual takeover of the global heat transport by the contribution from the regions of plume ejection. In future work we will extend these 2D DNS to smaller (and larger) Pr, to check the predicted Pr-dependence Grossmann and Lohse (2000, 2001) of the transition to the ultimate regime.

Many open questions remain, for example whether wall roughness can trigger a transition to an asymptotic ultimate regime, in which NuRa, i.e., the logarithmic corrections vanish. A previous study that reached Ra = has shown that this was not yet the case Zhu et al. (2017). However, in rough wall Taylor-Couette (TC) simulations (reaching a Taylor number of Ta ) and experiments (reaching Ta ) we did reach the corresponding asymptotic ultimate regime for the angular momentum transport in TC flow thanks to the effect of pressure drag Zhu et al. (2018). As the analog to pressure drag is absent in the heat flux balance for RB flow, such an asymptotic ultimate regime may not exist in RB flow Owen and Thomson (1963).

###### Acknowledgements.

We thank Daniel Chung for discussions and for pointing us to Ref. Owen and Thomson (1963). The work was financially supported by NWO-I, NWO-TTW, the Netherlands Center for Multiscale Catalytic Energy Conversion (MCEC), all sponsored by the Netherlands Organization for Scientific Research (NWO), and the COST Action MP1305. Part of the simulations were carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. We also acknowledge PRACE for awarding us access to Marconi based in Italy at CINECA under PRACE Project No. 2016143351 and the DECI resource Archer based in the United Kingdom at Edinburgh with support from the PRACE aisbl under Project No. 13DECI0246.## References

- Ahlers et al. (2009) G. Ahlers, S. Grossmann, and D. Lohse, Rev. Mod. Phys. 81, 503 (2009).
- Lohse and Xia (2010) D. Lohse and K.-Q. Xia, Ann. Rev. Fluid Mech. 42, 335 (2010).
- Chilla and Schumacher (2012) F. Chilla and J. Schumacher, Eur. Phys. J. E 35, 58 (2012).
- Sun et al. (2008) C. Sun, Y. H. Cheung, and K.-Q. Xia, J. Fluid Mech. 605, 79 (2008).
- Zhou and Xia (2010) Q. Zhou and K.-Q. Xia, Phys. Rev. Lett. 104, 104301 (2010).
- Zhou et al. (2010) Q. Zhou, R. J. A. M. Stevens, K. Sugiyama, S. Grossmann, D. Lohse, and K.-Q. Xia, J. Fluid Mech. 664, 297 (2010).
- du Puits and Willert (2016) R. du Puits and C. Willert, Phys. Fluids 28, 044108 (2016).
- Kraichnan (1962) R. H. Kraichnan, Phys. Fluids 5, 1374 (1962).
- Spiegel (1971) E. A. Spiegel, Ann. Rev. Astron. Astrophys. 9, 323 (1971).
- Grossmann and Lohse (2000) S. Grossmann and D. Lohse, J. Fluid. Mech. 407, 27 (2000).
- Grossmann and Lohse (2001) S. Grossmann and D. Lohse, Phys. Rev. Lett. 86, 3316 (2001).
- Grossmann and Lohse (2011) S. Grossmann and D. Lohse, Phys. Fluids 23, 045108 (2011).
- Grossmann and Lohse (2012) S. Grossmann and D. Lohse, Phys. Fluids 24, 125103 (2012).
- Stevens et al. (2013) R. J. A. M. Stevens, E. P. van der Poel, S. Grossmann, and D. Lohse, J. Fluid Mech. 730, 295 (2013).
- Malkus (1954) M. V. R. Malkus, Proc. R. Soc. London A 225, 196 (1954).
- Priestley (1954) C. H. B. Priestley, Aust. J. Phys. 7, 176 (1954).
- Niemela and Sreenivasan (2010) J. Niemela and K. R. Sreenivasan, New J. Phys. 12, 115002 (2010).
- Urban et al. (2014) P. Urban, P. Hanzelka, V. Musilová, T. Králík, M. La Mantia, A. Srnka, and L. Skrbek, New J. Phys. 16, 053042 (2014).
- Chavanne et al. (1997) X. Chavanne, F. Chilla, B. Castaing, B. Hebral, B. Chabaud, and J. Chaussy, Phys. Rev. Lett. 79, 3648 (1997).
- Chavanne et al. (2001) X. Chavanne, F. Chilla, B. Chabaud, B. Castaing, and B. Hebral, Phys. Fluids 13, 1300 (2001).
- He et al. (2012a) X. He, D. Funfschilling, H. Nobach, E. Bodenschatz, and G. Ahlers, Phys. Rev. Lett. 108, 024502 (2012a).
- He et al. (2012b) X. He, D. Funfschilling, E. Bodenschatz, and G. Ahlers, New J. Phys.. 14, 063030 (2012b).
- Stevens et al. (2011) R. J. A. M. Stevens, D. Lohse, and R. Verzicco, J. Fluid Mech. 688, 31 (2011).
- Shishkina et al. (2010) O. Shishkina, R. J. A. M. Stevens, S. Grossmann, and D. Lohse, New J. Phys. 12, 075022 (2010).
- Stevens et al. (2010) R. J. A. M. Stevens, R. Verzicco, and D. Lohse, J. Fluid Mech. 643, 495 (2010).
- Schmalzl et al. (2004) J. Schmalzl, M. Breuer, S. Wessling, and U. Hansen, Europhys. Lett. 67, 390 (2004).
- van der Poel et al. (2013) E. P. van der Poel, R. J. A. M. Stevens, and D. Lohse, J. Fluid Mech. 736, 177 (2013).
- Huang and Zhou (2013) Y.-X. Huang and Q. Zhou, J. Fluid Mech. 737 (2013).
- Zhang et al. (2017) Y. Zhang, Y.-X. Huang, N. Jiang, Y.-L. Liu, Z.-M. Lu, X. Qiu, and Q. Zhou, Phys. Rev. E 96, 023105 (2017).
- Hewitt et al. (2012) D. R. Hewitt, J. A. Neufeld, and J. R. Lister, Phys. Rev. Lett. 108, 224503 (2012).
- Ahlers et al. (2012) G. Ahlers, E. Bodenschatz, D. Funfschilling, S. Grossmann, X. He, D. Lohse, R. J. A. M. Stevens, and R. Verzicco, Phys. Rev. Lett. 109, 114501 (2012).
- Ahlers et al. (2014) G. Ahlers, E. Bodenschatz, and X. He, J. Fluid Mech. 758, 436 (2014).
- van der Poel et al. (2015a) E. P. van der Poel, R. Ostilla-Mónico, R. Verzicco, S. Grossmann, and D. Lohse, Phys. Rev. Lett. 115, 154501 (2015a).
- Perry and Chong (1982) A. E. Perry and M. S. Chong, J. Fluid Mech. 119, 106 (1982).
- Marusic et al. (2010) I. Marusic, B. J. McKeon, P. A. Monkewitz, H. M. Nagib, A. J. Smits, and K. R. Sreenivasan, Phys. Fluids 22, 065103 (2010).
- Smits et al. (2011) A. J. Smits, B. J. McKeon, and I. Marusic, Ann. Rev. Fluid Mech. 43, 353 (2011).
- Johnston and Doering (2009) H. Johnston and C. R. Doering, Phys. Rev. Lett. 102, 064501 (2009).
- See Supplemental Material at [URL will be inserted by publisher] for numerical details (2017) See Supplemental Material at [URL will be inserted by publisher] for numerical details, (2017).
- Verzicco and Orlandi (1996) R. Verzicco and P. Orlandi, J. Comput. Phys. 123, 402 (1996).
- van der Poel et al. (2015b) E. P. van der Poel, R. Ostilla-Mónico, J. Donners, and R. Verzicco, Computers & Fluids 116, 10 (2015b).
- van der Poel et al. (2014) E. P. van der Poel, R. Ostilla-Mónico, R. Verzicco, and D. Lohse, Phys. Rev. E 90, 013017 (2014).
- Pope (2000) S. B. Pope, Turbulent Flow (Cambridge University Press, Cambridge, 2000).
- Yaglom (1979) A. M. Yaglom, Annu. Rev. Fluid Mech. 11, 505 (1979).
- Samanta et al. (2014) D. Samanta, F. Ingremeau, R. Cerbus, T. Tran, W. I. Goldburg, P. Chakraborty, and H. Kellay, Phys. Rev. Lett. 113, 024504 (2014).
- L’vov et al. (2009) V. S. L’vov, I. Procaccia, and O. Rudenko, Phys. Rev. E 79, 045304 (2009).
- Zhu et al. (2017) X. Zhu, R. A. J. M. Stevens, R. Verzicco, and D. Lohse, Phys. Rev. Lett. 119, 154501 (2017).
- Zhu et al. (2018) X. Zhu, R. A. Verschoof, D. Bakhuis, S. G. Huisman, R. Verzicco, C. Sun, and D. Lohse, Nat. phys. 14, 417 (2018).
- Owen and Thomson (1963) P. R. Owen and W. R. Thomson, J. Fluid Mech. 15, 321 (1963).