Flow organization and heat transfer in turbulent wall sheared thermal convection

We perform direct numerical simulations of wall sheared Rayleigh-B\'enard (RB) convection for Rayleigh numbers up to $Ra=10^8$, Prandtl number unity, and wall Reynolds numbers up to $Re_w=10000$. Using the Monin-Obukhov length $L$ we identify three different flow states, a buoyancy dominated regime ($L \lesssim \lambda_{\theta}$; with $\lambda_{\theta}$ the thermal boundary layer thickness), a transitional regime ($0.5H \gtrsim L \gtrsim \lambda_{\theta}$; with $H$ the height of the domain), and a shear dominated regime ($L\gtrsim 0.5H$). In the buoyancy dominated regime the flow dynamics are similar to that of turbulent thermal convection. The transitional regime is characterized by rolls that are increasingly elongated with increasing shear. The flow in the shear dominated regime consists of very large-scale meandering rolls, similar to the ones found in conventional Couette flow. As a consequence of these different flow regimes, for fixed $Ra$ and with increasing shear, the heat transfer first decreases, due to the breakup of the thermal rolls, and then increases at the beginning of the shear dominated regime. For $L \gtrsim 0.5H$ the Nusselt number $Nu$ effectively scales as $Nu \sim Ra^{\alpha}$, with $\alpha \ll 1/3$ while we find $\alpha \simeq 0.31$ in the buoyancy dominated regime. In the transitional regime the effective scaling exponent is $\alpha>1/3$, but the temperature and velocity profiles in this regime are not logarithmic yet, thus indicating transient dynamics and not the ultimate regime of thermal convection.


Introduction
Rayleigh-Bénard (RB) convection, i.e. the flow in a box heated from below and cooled from above, is one of the paradigmatic fluid dynamical systems (Ahlers et al. 2009;Lohse & Xia 2010;Chilla & Schumacher 2012;Xia 2013). The dynamics of RB convection are controlled by the Rayleigh number Ra, which is the non-dimensional temperature difference between the horizontal plates, and the Prandtl number P r, defined as the ratio of momentum and thermal diffusivities. For a definition of these and other parameters we refer the reader to section 2.
For strong enough thermal driving, i.e. high enough Ra, the flow in the bulk region becomes fully turbulent. For even stronger thermal driving, beyond some critical Ra number Ra c , also the boundary layers become turbulent and the system reaches the regime of so-called ultimate convection (Kraichnan 1962;Grossmann & Lohse 2000, 2001. This ultimate regime sets in when the shear Reynolds number at the boundary layers is sufficiently high so that the boundary layer becomes turbulent, leading to a strong increase in the heat transport, which in dimensionless form is quantified by the Nusselt number N u. Ahlers et al. (2012) found that the transition to the ultimate regime sets in around Ra c ∼ O(10 14 ). While in the classical regime one generally finds N u ∼ Ra 0.31 , in the ultimate regime N u ∼ Ra 0.38 , in agreement with theoretical predictions (Grossmann & Lohse 2011).
The transition to the ultimate regime has also been observed in direct numerical simulations (DNS) of two-dimensional RB convection (Zhu et al. 2018a). In Taylor-Couette flow, which is a very analogous system, experiments and DNS have observed the ultimate regime as well (Grossmann et al. 2016). However, so far the ultimate regime has not yet been achieved in DNS of three-dimensional RB flows (Stevens et al. 2010(Stevens et al. , 2011 as the required computational time to achieve this is still out of reach. In an attempt to trigger the transition to the ultimate regime, here we add a Couette type shearing to the RB system to increase the shear Reynolds number of the boundary layers. In Couette flow the top and bottom walls move in opposite directions (Thurlow & Klewicki 2000;Barkley & Tuckerman 2005;Tuckerman & Barkley 2011) and just as in other examples of wall bounded turbulence (Jiménez 2018;Smits et al. 2011;Smits & Marusic 2013) the flow is dominated by elongated streaks, which have been observed in experiments (Kitoh & Umeki 2008) and DNS (Lee & Kim 1991;Tsukahara et al. 2006), even at relatively low Reynolds numbers (Chantry et al. 2017). Pirozzoli et al. (2011Pirozzoli et al. ( , 2014 and Orlandi et al. (2015) showed that these streaks in Couette flow have much longer characteristic length scales than in Poiseuille flow, where the flow is forced by a uniform pressure gradient rather than by wall shear. Rawat et al. (2015) showed that these large-scale flow structures even survive when the small-scale structures are artificially surpressed. Recently, Lee & Moser (2018) found that the streak length increases with increasing Reynolds number and that some correlation in the streamwise direction remains visible up to a length of almost 80 times the distance between the plates.
Investigating the interaction between buoyancy and shear effects is also very important to better understand oceanic and atmospheric flows (Deardorff 1972;Moeng 1984;Khanna & Brasseur 1998). For example early experiments on sheared thermal convection by Ingersoll (1966) and Solomon & Gollub (1990) showed the appearance of large-scale structures. Fukui & Nakajima (1985) showed that in channel flow unstable stratification increases the longitudinal velocity fluctuations close to the wall, while in the bulk region the temperature fluctuations are drastically lowered. Channel flow with stable stratification was further investigated by Armenio & Sarkar (2002) and García-Villalba & delÁlamo (2011) who observed that for high Reynolds numbers or weak stratification the buoyancy in the bulk was dominant while only for strong stratification large laminar patches emerged in the near-wall region. Furthermore, recent experiments by Shevkar et al. (2019) investigated the plume spacing in sheared convection and found a scaling law which connects the mean spacing of the plumes with Re w , Ra, and P r of the flow.
Early simulations of sheared convection were performed by Hathaway & Somerville (1986) and Domaradzki & Metcalfe (1988) for Ra O(10 5 ). Domaradzki & Metcalfe (1988) found that in Couette-RB flow the addition of shear at low Ra initially increases the heat transport. However, for Ra 150.000 the heat transport decreases as the added shear breaks up the large-scale structures. More recently, Scagliarini et al. (2014Scagliarini et al. ( , 2015 showed that also in Poiseuille-RB the heat transfer first decreases when the applied pressure gradient is increased. The reason is that for intermediate forcing the longitudinal wind disturbs the thermal plumes which therefore lose their coherence. Only with a strong enough pressure gradient, a heat transfer enhancement is found. The Richardson number, to be defined below, quantifies the ratio between the buoyancy  and shear forces in Couette-RB and Poiseuille-RB based on the applied temperature difference and wall shear Reynolds number. Another way to quantify the buoyancy and shear forces is to determine the Monin-Obukhov length (Monin & Obukhov 1954;Obukhov 1971), which indicates up to which distance from the wall the flow is dominated by shear, based on the observed flow properties. Pirozzoli et al. (2017) found that the Monin-Obukhov length scales as L ≈ 0.15H/Ri 0.85 for channel flow with unstable stratification.
In this study we investigate the effect of an additional Couette type shearing on the heat transfer in RB convection in an attempt to trigger the boundary layers to become fully turbulent and hence observe the transition to the ultimate regime. Figure 1 shows a flow visualization of the temperature field obtained from one of our simulations, which reveals large-scale meandering streaks that are formed near the hot plate. We performed simulations over a wide parameter range, spanning 10 6 Ra 10 8 and 0 Re w 10 4 , while P r = 1 has been used in all cases, see figure 2a. In spite of the very strong forcing for the largest Ra and Re w , we did not yet achieve ultimate turbulence in this study. We were limited by our own requirement of using large domain sizes as recommended by Pirozzoli et al. (2017) to ensure convergence of the main flow properties.
The remainder of this manuscript is organized as follows. In §2 we present the simulation method. We discuss the heat transfer and skin friction measurements in §3.1 and §3.2, respectively. A discussion of the identified flow regimes is given in §4. The concluding remarks follow in §5.

Simulation details
We numerically solve the three-dimensional incompressible Navier-Stokes equations within the Boussinesq approximation, which in non-dimensional form read: with u and θ the non-dimensionalized velocity and temperature. Since for this problem many parameters have to be considered, it is worthwile to distinguish between input and output of the flow. Input parameters are the Rayleigh number Ra (2.3), with H the distance between the plates, β the thermal expansion coefficient of the fluid, g the gravitational acceleration, ∆ the temperature difference between top and bottom plate, and κ and ν the thermal and kinetic diffusivities, respectively, the Prandtl number P r (2.4), the wall shear Reynolds number Re w (2.5), with u w the velocity of the wall, and the bulk Richardson number Ri (2.6). Output parameters are the Nusselt number N u (2.7), where Q = q/(ρc p ) is the total vertical heat flux, with q the turbulent heat flux and c p the constant pressure specific heat, the Taylor Reynolds number Re τ (2.8), where u τ = τ w /ρ is the friction velocity, with τ w the mean wall shear stress and ρ the density of the fluid, the Monin-Obukhov length L (2.9), and the skin friction coefficient C f (2.10).
To solve (2.1) -(2.2) we employ the second-order finite difference code AFiD (van der Poel et al. 2015), which has been validated many times against other numerical and experimental results (Verzicco & Orlandi 1996;Verzicco & Camussi 1997, 2003Stevens et al. 2010Stevens et al. , 2011Ostilla-Mónico et al. 2014;Kooij et al. 2018). The code uses periodic boundary conditions with uniform mesh spacing in the horizontal directions and supports a non-uniform grid distribution in the wall-normal direction. Here we adopted the vertical grid distribution used by Pirozzoli et al. (2014) for Couette flow. For this study we used the GPU version of the code (Zhu et al. 2018b) to allow efficient execution of many largescale simulations. The Couette flow forcing is realized by moving both walls in opposite directions with speed u w and the results for the pure Couette flow case match excellently with the results by Pirozzoli et al. (2014). For example, figure 2b shows that for Couette flow Re τ ∼ Re 0.85 w , which agrees very well with the Couette data of Pirozzoli et al. (2014) and Avsarkisov et al. (2014).
All simulations in this study were performed in a large 9πH × 4πH × H box, in streamwise, spanwise and wall-normal direction (Tsukahara et al. 2006;Pirozzoli et al. 2014), which is required to capture the large-scale structures formed in the Couette (Pirozzoli et al. 2014;Avsarkisov et al. 2014;Lee & Moser 2018). The grid resolutions are based on those used by Pirozzoli et al. (2014). To verify that the adopted grid spacing is adequate for our code, we performed a grid refinement check for the Ra = 10 8 and Re w = 10000 case, i.e. the most challenging simulation of this study. To keep this resolution test manageable it is performed in a smaller 2πH × πH × H domain, see table 1. Figure 3 confirms that the simulations are fully resolved for the chosen resolution. As a further validation, we evaluate the Couette data from Pirozzoli et al. (2014) in §3.2, which collapses very well with our data. Table 2 shows the simulation parameters for the main cases presented in this study.

Global flow characteristics
3.1. Effective scaling of the Nusselt number Figure 4 shows that the heat transfer increases with increasing Ra and Re w and that for given Ra number a minimum heat transfer is obtained at some intermediate Re w . Figure  5a shows the corresponding cross sections for constant Ra which clearly reveal that the location of the minimum heat transfer at constant Ra shifts towards higher Re w with increasing Ra.
The results indicate that this mechanism is influenced by the ratio of the buoyancy and shear forces. Therefore the bulk Richardson number Ri or the Monin-Obukhov length L, which take the ratio of these forces into account, are natural parameters to identify the different flow regimes. As L can be compared to other important length   scales in the flow, we decide to use it for this purpose. From the data in table 2, we find L ≈ 0.16/Ri 0.91 (see appendix C). In figure 6a the Monin-Obukhov length is compared to the thermal boundary layer thickness λ θ and the arbitrary threshold 0.5H is reported for later discussion. Since L/H is the fraction of the domain in which the thermal forcing is dominant in the flow, L 0.5H is the threshold from which on the flow is completely shear dominated since the wall generated shear affects at least half of the fluid layer thickness. This allows us to define three different flow regimes, namely a buoyancy dominated regime (L λ θ ), a transitional regime (0.5H L λ θ ), and a shear dominated regime (L 0.5H). A similar behavior has also been observed in convective boundary layers, where Salesky et al. (2017) find a cell dominated regime for z i /L > 20, where z i is the convective boundary layer thickness, a cell and roll dominated regime as transitional state, and a roll dominated regime for z i /L < 5. Figure 6b shows that the heat transfer in the buoyancy dominated regime scales as For the shear dominated regime we find that the effective scaling exponent α in N u ∼ Ra α is α 1/3 and in the transitional regime we find α > 1/3. An effective scaling exponent larger than 1/3 is one of the characteristics of the ultimate regime. It should occur when the boundary layers have transitioned to the turbulent state, which is indicated by their logarithmic profiles. Our analysis in §4 will show that this is not yet the case in this transitional regime. Instead, for intermediate shear, the heat transfer is decreased with respect to the RB case. The locally larger effective scaling exponent simply is a consequence of the fact that with increasing Ra the heat transfer must again converge to the RB case.

Skin Friction
In figure 7 we compare the measured skin friction coefficient for different Re w and Ra with Prandtl's turbulent friction law (Schlichting & Gersten 2000): Following Pirozzoli et al. (2014) we use a von Kármán constant K = 0.41 and C = 5. The figure shows that the skin friction increases with Ra and decreases with Re w . At fixed Ra the relative strength of the thermal forcing decreases for high Re w and therefore the obtained friction coefficient converges to the Prandtl law. In figure 7b we focus on the data for small Re w . The skin friction in pure Couette flow follows the expected laminar result C f = 4/Re w (Pope 2000) until a transition to the turbulent state occurs around Re w = 650−700. Cerbus et al. (2018) discuss that in pipe flow this jump is caused by the formation of puffs and slugs. Brethouwer et al. (2012) attribute this discontinuous jump in C f to the lack of restoring forces in plane Couette flow (similar to pipe, channel and boundary layer flows). For the Couette-RB case we do not observe such a discontinuous jump. Instead this sheared RB case is another example, next to the application of Coriolis,  buoyancy, and Lorentz forces discussed by Brethouwer et al. (2012), which shows that restoring forces can prevent a discontinuous jump in C f (Re w ). Chantry et al. (2017), on the other hand, claim that all transitions to turbulence should be continuous if the used box size is large enough.

Organization of turbulent structures
To further investigate the dynamics of the different regimes we show visualizations of the temperature field at midheight for all simulations in figure 8 and Appendix A. In the thermally dominated regime the primary flow structure resembles the large-scale flow found in RB convection ). In the transitional regime (L 0.5H), the thermal forcing dominates part of the bulk where large elongated thermal plumes transform into thin straight elongated streaks when L approaches 0.5H. In figure 8 and in the appendix this manifests itself as a very visual line diagonally through the diagram, splitting the more thermal-and the more shear dominated cases. In the shear dominated regime (L 0.5H) we find large-scale meandering structures, similar to the ones found in pressure-driven channel flow with unstable stratification (Pirozzoli et al. 2017). This significant change in flow structure can be linked to the minimum in N u in figure 5. The reason for the minimum is that at intermediate shear the thermal convection rolls are broken up, while the shear is not yet strong enough to increase the heat transfer directly. This observation is in agreement with earlier works described above (Domaradzki & Metcalfe 1988;Scagliarini et al. 2014Scagliarini et al. , 2015Pirozzoli et al. 2017). Figure 9 combines these findings with the above observation that in the shear dominated regime the effective scaling exponent α in N u ∼ Ra α is much smaller than 1/3, in the transitional regime α > 1/3, and in the thermally dominated regime α 0.31. When we compare the regime transitions with the results in figure 5, it becomes clear that the lowest heat transfer for a given Ra occurs at the end of the transitional regime just before the emergence of the thin straight elongated streaks. Due to the large computational time that is required for each simulation the number of considered cases is limited, which makes it difficult to pinpoint exactly when the heat transfer is minimal and what the  flow structure looks like in that case. However, we note that the onset of the shear dominant regime corresponds to the point where the heat transfer starts to increase as the additional shear can then more effectively enhance the overall heat transport.
To get more insight into the boundary layer dynamics in the different regimes, we show the temperature and streamwise velocity at boundary layer height for Ra = 4.6 × 10 6 in figure 10. At this Rayleigh the flow is in the transitional regime for Re w = 2000 and Re w = 3000, and in the shear dominant regime for Re w 4000. For all cases we observe a clear imprint of the large-scale structures observed at midheight, see figure  8 and Appendix A. This indicates that the large-scale dynamics have a pronounced influence on the flow structures in the boundary layers . The figure also reveals that in the transitional and shear dominated regime the lowest temperatures at boundary layer height are observed in the high speed streak regions, which indicates that the regions with the highest shear contribute most to the overall heat flux.

Flow statistics
We now present the streamwise temperature variance spectra E θ (k) in figure 11 to analyze the size of the large-scale structures as a function of the Monin-Obukhov length. The position of the peak in the temperature spectrum indicates the wavelength of the most prominent thermal structure. In panel (b) and (c) we plot the evolution of the wavelength of these structures in relation to the absolute size of the flow field, where L x = 2π/k x peak and L y = 2π/k y peak . If the spectrum does not show a clear peak, but keeps growing for small k, the structure size is set to the limits of the simulation box, which is Γ x = 9π in streamwise (figure 11b) and Γ y = 4π in spanwise direction (figure 11c) in this manuscript. For L → 0 we expect L x , L y ≈ 2πH ) and here we observe that for L → ∞ L x = Γ x H and L y ≈ 0.8πH ≈ 0.2Γ y . In panel (b) we observe that the large-scale structures mainly extend over the whole boxlength Γ x H. When we compare this to figure 2 and figure 14, it becomes very clear that for small L, i.e. in the transitional regime as the RB case (buoyancy dominated regime) is not shown due to the logarithmic axis, the large-scale structures are elongated over the whole streamwise length. For pure RB convection, where L = 0, L x is much lower in agreement with the above expectation . Only once the meandering behavior starts, the structure size drops to about half the box length. If L is further increased, the flow undergoes a transition towards pure Couette flow and L x → Γ x H. For L y in panel (c) we can see a different behavior. For small L, we see a transition from L y = Γ y H to much smaller structures. Here we can observe very thin, elongated streaks in the flow field. When L is large enough to reach the shear dominated regime, the meandering structures can be identified by the local peak in L y . For L → ∞, L y converges towards 0.2Γ y H. Due to the very limited number of datapoints, it is not possible to fully assess the behavior of L x and L y vs L for all Ra and Re. Nevertheless the minimum in L x and peak in L y in the shear dominated regime are very distinct.
To further quantify the cases shown in figure 10, we study their flow statistics in figure  12. It becomes clear that both the temperature and streamwise velocity profiles are not logarithmic in the transitional regime. This indicates that the boundary layers are not turbulent in this state. Hence, the higher N u scaling in the transitional regime does not seem to be caused by triggering the ultimate regime. In the shear dominated regime the streamwise velocity and temperature profiles seem to converge to a logarithmic profile with increasing Re w .
In figure 13 we show the same statistical quantities as in figure 12, but now for fixed Re w = 6000. For Ra 10 7 the flow is in the transitional regime and for Ra 10 7 the flow undergoes a transition into the shear dominated regime. Just as in figure 12 we observe that the temperature and streamwise velocity profiles are not logarithmic in the transitional regime. As the Richardson number decreases with decreasing Ra, we see that the profiles converge towards a logarithmic behavior. From a comparison with table 2 we find that Ri 0.2 seems to be required to achieve logarithmic temperature and velocity profiles. For the parameter regime under investigation the effective scaling exponent α in this regime is well below 1/3. In both figures we can detect a non-monotonic behavior of both u + and T + for low Re w and high Ra. This is connected to an effect of flow layering in the transitional regime where the large-scale thermal plumes get distorted by the shear, but not enough for the meandering structures to evolve and break up this effect. Several further statistical quantities for both constant Ra and constant Re have been calculated and can be found in appendix B.

Concluding remarks
We performed direct numerical simulations of turbulent thermal convection with Couette type flow shearing. We presented cases in a range 10 6 Ra 10 8 and 0 Re w 10 4 , Figure 13. (a) Mean streamwise velocity and (b) temperature profiles, where u + = u/uτ and T + = T /Tτ with Tτ = Q/uτ for Rew = 6000. T + Ra=0 was determined through a passive-scalar temperature field.
achieving up to Re τ ≈ 370. For fixed Rayleigh number we obtain a non-monotonic progression of N u similarly to what was previously observed in unstable stratification with a pressure gradient (Scagliarini et al. 2014). The addition of imposed shear to thermal convection first leads to a reduction of the heat transport by disrupting the turbulent system before the shear becomes strong enough to create meandering streaks that efficiently transport the heat away from the wall. As the impact of the thermal plumes on the flow decreases with increasing shear, the skin friction coefficient at constant Ra drops with increasing Re w .
We find that three flow regimes can be identified in Couette-RB using the Monin-Obukhov length L and the thermal boundary layer thickness λ θ . In the buoyancy dominated regime (L λ θ ) the flow is dominated by large thermal plumes. With decreasing Richardson number we first find a transitional regime (0.5H L λ θ ), before the shear dominated flow regime with large-scale meandering streaks is obtained. For a given Ra the minimum heat transport is found just before the onset of this shear dominated regime when thin straight elongated streaks dominate the flow. We find that in the transitional regime the effective scaling exponent α in N u ∼ Ra α is larger than 1/3. An analysis of the flow characteristics shows that the temperature and streamwise velocity profiles are not logarithmic in this transitional regime, which one would expect when this high scaling exponent would indicate the onset of the ultimate regime. Since it is possible to recover logarithmic profiles for low Richardson number flows we want to investigate in future studies whether it is possible to increase the thermal and sheared forcing far enough to trigger ultimate convection in Couette-RB.