Stability Analysis of Sheared Thermal Boundary Layers and its Implication for Modelling Turbulent Rayleigh-Bénard Convection

Predicting the heat flux through a horizontal layer of fluid confined between a hot bottom plate and a cold top one has always spurred theoretical, numerical and experimental work on Rayleigh–B´enard convection. Customarily, the Nusselt number (the heat flux in non-dimensional form) has been modelled in the form of one or several power-laws of three parameters, the Rayleigh, Prandtl and Reynolds numbers. Quantifying the large-scale flow that spontaneously develops in a turbulent Rayleigh–B´enard cell, the Reynolds number, unlike the Rayleigh and Prandtl numbers, is not a control parameter strictly speaking and, depending on the model, is sought as another power-law or introduced as an external input. Whereas balancing the di ff erent transport mechanisms can predict the exponents in these power laws, experimental and numerical results are required to adjust the various prefactors. The early and simple model of Malkus [1] and Howard [2] assumed that the value of the Nusselt number could be directly deduced from the marginal stability of the two sheared thermal boundary layers along the upper and lower plates, interacting via the large-scale flow. Maintaining this simplicity, this work shows that in the classical regime of turbulent convection, considering the linear critical conditions of absolute (as opposed to convective) thermo-convective instabilities alleviates the flaws of the original model. Revisiting available Direct Numerical Simulations from which a Reynolds number can be unambiguously extracted, the present approach then yields the Nusselt number as a function of the Rayleigh and Prandtl numbers agreeing well with the numerical results.


Introduction
Rayleigh-Bénard (RB) convection has attracted the interest of the scientific community for more than a century as, on the one hand, its simple configuration (a fluid confined between two, hot and cold, horizontal plates) facilitates its study and, on the other hand, it encompasses the basic physics occurring in numerous natural or industrial flows.Rayleigh-Bénard convection has therefore played a crucial role in the study of hydro-dynamic instabilities, from the early concepts to more involved ones as spatio-temporal chaos [3] and constitutes a convenient set-up to experimentally and numerically investigate turbulent processes.The driving force of thermal convection, the buoyancy, is quantified by the Rayleigh number Ra = gα∆h 3 /(νκ) based on the thermal expansion coefficient α, kinematic viscosity ν and thermal diffusivity κ of the fluid, and the acceleration due to gravity g, the vertical distance and the temper-ature difference between the hot, bottom and top, cold plates, h and ∆, respectively.The main outcome and metric of RB convection, the vertical heat flux across the cell, is quantified by the Nusselt number Nu = Φh/(λ∆) based on the horizontally averaged heat flux Φ and λ the thermal conductivity of the fluid.Moreover, the fluidspecific competition between heat diffusion and viscosity is quantified by the Prandtl number Pr = ν/κ.For Rayleigh numbers between 10 5 and 10 11 − 10 12 , experimental and numerical results all converge towards a universal functional dependence of Nu with Ra and Pr [4].In addition to Nu, a secondary response parameter is the Reynolds number Re = Vh/ν, characterizing the velocity field observed in the RB cell.A Reynolds number is more difficult to define, investigate and measure than Nu.Indeed, for fixed control parameters and in steady state, Φ is constant and equal to the heat flux through the two horizontal plates, whereas the velocity field results from the complex combination of turbulent motions and coherent structures, as plumes, boundary layers and large-scale circulation in the bulk of the cell, also known as "wind" [5].Consequently, several definitions are possible for the characteristic velocity V, such as the root-mean-square velocity averaged over the whole fluid, V r.m.s.= ⟨V 2 x + V 2 y + V 2 z ⟩, or some magnitude of the wind flow (V w ).Nevertheless, a functional dependence of Re with Ra and Pr also seems to prevail [6,7].These converging results are noteworthy for such a system, in which boundary layers (BLs) and turbulent flows in the bulk interact in a complex fashion.
On the other hand, in configurations where Ra is larger than 10 11 − 10 12 , a controversy has been going on for more than twenty years on the existence and conditions for the appearance of a convection regime called "ultimate regime" [8,9].For these large Rayleigh numbers, experimental results seem to call into question a functional dependence of Nu and Re with Ra and Pr.The present model covers the classical convection regime (Ra ⪅ 10 10 ) and there is no evidence that it might be applicable to the ultimate regime.Indeed, the ultimate regime is associated with boundary layers having transitioned from laminar to turbulent behaviour [8,9], which contradicts the assumptions used in the present theory.
Numerous theories have been developed to predict the functional dependence Nu(Ra, Pr) in the classical regime [see 5, 10, 3, for detailed and thorough accounts].The widely accepted Grossmann and Lohse (GL) model [11,12] is based on splitting the mean kinetic energy and thermal dissipation rates into two contributions each, one from the bulk and one from the boundary layers.It assumes that RB convection is a mixture of eight convection regimes, Nu and Re w = V w h/ν being described by power-laws of Ra and Pr for each of these regimes.The Nusselt and Reynolds numbers are then obtained as functions of Ra and Pr, using the two equations: The function f captures the cross-over of a thermal BL of thickness δ T nested in the kinetic one of thickness δ V to the reciprocal situation.For very large Pr's, the function g describes the saturation of the kinetic BL thickness as Re w decreases below a critical Reynolds number.This critical Reynolds number and the four prefactors c 1 -c 4 have been computed in [13] by fitting (1) to experimental and numerical results.It was recently proposed [7] to use functional forms for the prefactors c 1c 4 to improve the predictions of (1).Note that whether from experiments or simulations, obtaining the velocity of the wind V w remains a major challenge.It is generally assumed that V w = ξ V r.m.s., with ξ = 1, the value of ξ only impacting the model constants [11,12].
Before these fine-tuned quantitative models, Malkus [1] and Howard [2] had suggested that in turbulent convection (Ra ⪆ 10 5 ), the lower and upper thermal boundary layers would self-adjust to reach the critical conditions, above which thermal convection rolls develop.These critical conditions are accessible by linear stability analysis.They deduced that Nu ∝ Ra 1/3 (Eq.5), in agreement with the dimensional analysis proposed by Priestley [14].This exponent 1/3 is indeed the only scaling for Nu ensuring that the heat flux is independent of the height of the cell h.To capture the departure from the 1/3 exponent observed in the experimental results, Castaing et al. [5] proposed a mixing zone model, distinguishing three regions in the RB cell.In addition to the two thermal boundary layers and the centre of the cell, they considered an intermediate region where sheets of fluid are ejected from the BLs.The scalings of the velocities of the ejected fluid with Ra then led to predict that Nu ∝ Ra 2/7 .This mixing zone model was then further developed to take into account the Pr-dependence [15,16], predicting Nu ∝ Ra 2/7 Pr 2/7 for small Pr and Nu ∝ Ra 2/7 Pr −1/7 for large enough Pr.
In order to agree with experimental and numerical results, these models have been increasingly complexified, departing from the seminal simplicity of Malkus input that has to be inferred from numerical simulations.

Revisiting the Malkus and Howard model
Malkus and Howard assumed that, for turbulent convection, the temperature averaged over time is almost uniform in the bulk flow and the mean temperature increases by ∆/2 when crossing the thickness δ T of each thermal boundary layer (see Fig. 1a).They also assumed that the bulk flow only serves to transmit the constant heat flux between the two boundary layers and has no impact on their stability.More precisely, we postulate that the bulk could be seen as a 'conveyor belt' advecting the perturbations of temperature and velocity from one boundary layer to the other.This crude picture advantageously discards the need of a mixing zone and, as underlined in [5], the Reynolds number The correction factor F comes from the fact that the wind velocity at the edge of each thermal BL is not systematically equal to V w (see Fig. 1 in [11]).For Pr ≤ 1, as δ V ≤ δ T , F should be equal to 1.For larger Pr, the relevant velocity should be less than V w , namely about Even for large Pr numbers, though, Direct Numerical Simulations (DNS) in [7] have shown that δ V remains very close to δ T , leading to F very close to unity.We will show later that variations in F have little impact on the results of the model.
For Ra ⪅ 10 10 , previous experimental and numerical works have confirmed the validity of (4) (see Fig. 2a of [7] for example).In contrast, for larger Ra numbers or for the ultimate regime, as a logarithmic mean temperature profile is expected and reported in experiments [18,19], (4) certainly no longer holds.Equations ( 2) and ( 4) then result in Irrespective of the velocity profile and boundary conditions, when thermo-convective instabilities of infinite extension along the direction of the shear are considered, rolls the axes of which are aligned with this direction are the first to become linearly unstable.Their critical Rayleigh number is found to be independent of the shear flow, the well-known value Ra 0 c ≈ 1708 being retrieved, for Dirichlet boundary conditions.Assuming Ra bl = Ra 0 c , equation ( 5) then shows that Nu must vary as Ra 1/3 .This approach, however, has three caveats.
First, the linear stability analysis of these longitudinal rolls show that Ra 0 c is independent of Re bl , so the stability of the BLs is obviously not affected by the wind flow.Then, as Ra 0 c is also found to be independent of the Prandtl number, so is the scaling (5).Finally, experiments usually observe that Ra bl > Ra 0 c and conclude that the BLs should be unstable [20], calling the initial assumption of this model into question.
We show in the forthcoming that it is possible to fix these caveats by considering the linear impulse response of the sheared thermal BLs instead of their stability to infinitely extended perturbations.The Green that the dynamics of turbulent Rayleigh-Bénard convection is driven by the marginal stability of the thermal boundary layers with respect to absolute modes and independent of the bulk.This amounts to check that in existing numerical and experimental results, Nu and Re w are related by equations ( 7)-( 8).

Bénard convection
To check the validity of our model, the critical threshold of absolute thermo-convective instabilities developing in the merged boundary layers Ra abs c (Re bl , Pr) must be computed beforehand.Although this analysis should be restricted to configurations that are homogeneous in the direction of the mean flow, it is assumed here that changes in the BLs along this direction remain sufficiently weak to be neglected in the analysis.
The configuration the stability of which is computed then consists of a fluid confined between two horizontal plates at rest, the lower one at temperature ∆/2 while the upper one, at distance 2δ T above, is set at temperature −∆/2, and swept by mean shear flow, as shown in Fig. 1(b).For the sake of simplicity, this shear flow is chosen in the shape of a piece-wise linear profile and it has been further checked that the results of the stability analysis were only minutely impacted by changing of this profile to a Poiseuille flow.The Rayleigh and Reynolds numbers are defined as above, using respectively (2) and (3), with F = 1.Using the streamwise, spanwise and wall-normal coordinates (x, y, z) and related basis, the velocity, pressure and temperature fields are decomposed into the steady base state V x,b , 0, 0, P b , Θ b , consisting of the laminar shear flow and conduction solution, and temporally evolving perturbations V x,p , V y,p , V z,p , P p , Θ p .
Following the procedure and non-dimensionalization scheme detailed in [22,23] for Rayleigh-Bénard-Poiseuille convection, the double curl of the Navier-Stokes equation and the heat equation, both linearized about the base state, together with the continuity equation in the Boussinesq approximation, yield a system of partial differential equations satisfied by the wall-normal velocity and temperature fields of the perturbation.Seeking these perturbations in the form V z,p , Θ p = v z (z) , θ(z) exp −iωt + ik x x + ik y y recasts these PDE's into the following generalized eigenvalue problem and boundary conditions the operators reading with d z the z-derivative, k 2 = k 2 x + k 2 y and ∆ = d 2 z − k 2 .The eigenvalues of ( 9) are the complex frequencies ω, the imaginary parts of which are the growth rates of the instabilities.
For a given set of parameters and wavenumber Re bl , Pr, Ra bl , k x , k y , problem ( 9) is solved by a taucollocation spectral method using Chebychev polynomials on 32 Gauss-Lobatto collocation points.For fixed values of Re bl and Pr, the critical conditions for the absolute instability are sought after as the set (Ra abs c , k x , k y ), with complex k x and k y , ensuring that 6 These critical conditions are computed using Newton-Raphson algorithms and all the linear algebra involved in those computations is accomplished using NAG (Numerical Algorithms Group, Oxford, UK) routines.
Whereas convective instabilities are known to take the form of longitudinal rolls the axis of which is aligned with the direction of the wind, i.e. k x = 0, absolute instabilities are always found in the form of transverse rolls, i.e. k y = 0. Figure 2(a) shows the results of this parametric stability analysis, namely the function Ra abs c (Re bl , Pr) to be used in equations ( 7)- (8).Note that besides the Reynolds number Re bl , this function would only marginally be affected by changes in the velocity profile V x,b (z) of the stability problem ( 9)- (10).

Comparison between model predictions and DNS results.
To put our interpretation on firmer grounds, we now proceed to test the relation (8) between Nu and Re w in existing numerical results.As for the GL model, it is assumed hereafter that V w = ξ V r.m.s., or equivalently Re w = ξ Re r.m.s., with ξ independent of Ra and Pr numbers.The stability analysis in §1 assumes a twodimensional base flow but, spanning both wavenumbers k x and k y , is fully three-dimensional.Selecting k y = 0, the outcome of this analysis is a two-dimensional flow in the form of transverse rolls.Thus, our model should be able to retrieve both 2D and 3D DNS results of Rayleigh-Bénard convection.However, perhaps counter-intuitively, the bulk flow is more complex in 2D than in 3D because of recirculation loops in the corners [26,27] and the coexistence of multiple statistically stable states when the ratio between horizontal and vertical extension of the 2D cell is greater than one [26,28].For these reasons, in what follows, we will compare our model with the results of 3D DNS for
which we have both: (i) the Reynolds number based on the root-mean-square velocity (V r.m.s. ), and (ii) a large enough aspect ratio and three-dimensional flow, so that a one-to-one dependence of Nu and Re r.m.s. with Ra and Pr prevails [29].
with b = 0.15.Using (11), ξ = 0.22 and F = 1, a nice agreement is observed for Nu between the DNS data and the predictions of Eqs. ( 7)-( 8), depicted by the thick solid blue curve in Fig. 3(b).This agreement was also obtained and extended to Pr 1 using a fit of the numerically computed values for Re r.m.s.instead of (11).
Figure 3(b) also shows that variations of ξ impact this result, particularly as Ra increases.There are no DNS results available so far to unambiguously calculate ξ.It seems reasonable, however, to find a typical wind flow velocity V w lower than the the root-mean-square veloc-ity.The velocity V w could actually be closer to the characteristic horizontal velocity, V h,r.m.s.= ⟨V 2 x + V 2 y ⟩/2.Using the results of [30,31], it is found V h,r.m.s./V r.m.s.≈ 0.5, a value closer to ξ = 0.22.fitted from DNS results as plotted in Fig. 4(a).For Ra Pr ≤ 10 9 , as δ V ≤ δ T , F (δ T /δ V ) = 1 [7] and changes in the ratio δ T /δ V have no effect on the theoretical results.For Ra Pr ≥ 10 9 , F (δ T /δ V ) = δ T /δ V can be calculated using the numerical results of [7]. Figure 4(b) shows that the effect of the ratio δ T /δ V is significant for Ra ≥ 10 9 only and Pr ≥ 10 (magenta solid line: F = 1, and dash-dotted line: F = δ T /δ V ).Even in this case, the effect of the ratio δ T /δ V is limited and numerical results are inconclusive to discriminate the theoretical results.
In agreement with experimental, numerical and previous theoretical approaches, this model shows that Nu does not strictly behave as a power-law of Ra and for curves: GL theory [12] with prefactors from [13].Solid curves in (b): present model ( 7)-( 8) with F = 1, ξ = 0.22 and Re r.m.s.fitted from DNS results depicted by the solid curves in (a).Dash-dotted magenta curve in (b): ( 7)-( 8) with ξ = 0.22, Re r.m.s.represented by magenta curve in (a) and the ratio F = δ T /δ V as given in [7].

Conclusions
It has been found in this study that the linear threshold of absolute instability of the sheared thermal boundary layers could quantitatively relates the heat flux through a RB cell and the magnitude of the large scale wind in this cell.On the basis of computing the critical Rayleigh number of absolute instability, Ra abs c (Re bl , Pr) (these numerical results are available upon request) of the merged sheared thermal boundary layers, equations ( 7) and ( 8) predict the variations of Nu as a function of the three parameters Ra, Pr and Re w .Whereas our approach supports the pivotal role of absolute thermoconvective instabilities in the dynamics of turbulent Rayleigh-Bénard convection, it does not yield, nonetheless, a quantitative model predicting the Nusselt number as a function of the Prandtl and Rayleigh numbers, as the wind Reynolds number Re w remains an external input.Moreover, whereas our results support the idea that the boundary conditions between the bulk and the thermal boundary layers might be unimportant to model the heat flux, they do not exclude that these boundary conditions could be pivotal for the wind, the velocity of which V w remains an external input in our model.Unfortunately, the wind characteristic velocity V w remains ambiguous to define and is seldom reported in the literature.On the contrary, V r.m.s. is univocally defined and frequently addressed, but its exact relation to V w remains an open question.
While two-and three-dimensional DNS and experiments with aspect ratios close to unity exhibit a single large scale circulation in the bulk, several cells separated by plumes have been reported in large aspect ratio two-dimensional DNS [19] and a complex network of plumes circumscribing regions of obvious sheared ther-mal boundary layers is observed experimentally [34].
These plumes obviously carry the heat through the bulk, but it remains unclear whether they contribute or not to the heat flux close to the plates and this latter could still be mostly imposed by the heat flux through the thermal boundary layers.Besides, two-dimensional DNS have also shown that this large-scale organization in cells and plumes is non-unique [28], though the threedimensional generalization of this result remains an open question.The relation between the wind and V r.m.s. is very likely non trivial and V w should probably be described in a statistical fashion [35].Our approach could nonetheless remain valid and useful as an "elementary brick" upon which to build a more complex model.Albeit crude, assuming a linear relation V w = ξV r.m.s to compute Re w from values of Re r.m.s. reported in the literature, our model still retrieves the corresponding Nu.The coefficient ξ is so far the only parameter the adjustment of which is required to validate the model against experimental or DNS results.In the ranges 10 5 ≤ Ra ≤ 10 9 and 0.02 ≤ Pr ≤ 200, setting ξ = 0.22 was found to be good enough to capture Nu as a function of Pr and Ra with a reasonable accuracy.This agreement could suggest that thermal plumes may have a limited impact on the total thermal flux.This ad hoc value ξ = 0.22 seems sensible since we find a typical wind flow velocity V w close to the mean horizontal velocities V h,r.m.s., as reported in DNS [30,31].Though always close to unity, the value of ξ nevertheless depends on the actual velocity profile used in the stability analysis of the thermal boundary layers.
Despite the hefty amount of existing experimental and numerical results, the characterization of the wind, its typical velocity, its fluctuations, its relation to the averaged velocity still require further work.Experimental measurements of the velocity and temperature in the thermal boundary layers remain particularly delicate and three-dimensional DNS in cells with large aspect ratios are still limited by their numerical cost.

and 9 and 10 − 2 <
Howard's interpretation.Furthermore, these improvements do not fix the original flaw of this interpretation: whatever the shear profile and boundary conditions, the linear threshold of thermo-convective instability of the thermal BL is unaffected by the shear and independent of the Prandtl number.This work establish that it is possible to eliminate this flaw by considering the absolute instabilities of the sheared thermal BLs instead of the convective ones, i.e. the impulse response instead of infinitely extended perturbations.After revisiting the early Malkus-Howard model and its shortcomings in §1, this work alleviate these latter by specifically considering the marginal stability of the thermal BLs in the framework of absolutely (as opposed to convectively) unstable sheared Rayleigh-Bénard (sRB) convection in §1.Section 1 shows that the prediction for Nu as a function of Ra, Pr and Re w of this model compares very favourably with existing results from the three-dimensional Direct Numerical Simulations from which both the Nusselt and the Reynolds numbers are available over the ranges 10 5 < Ra < 2•10 Pr < 10 3 .Section 1 finally concludes and further discusses the pivotal roles of the large-scale flow and Reynolds number, which, in our model, remains an

Figure 1 :
Figure 1: (a) Schematic view of the flow in a RB cell for Ra ⪆ 10 5 .Also shown: the mean temperature profile and the two thermal boundary layers of thickness δ T close to the bottom and top plates.∆ is the difference of temperature between the two horizontal plates.In the forthcoming analysis, the top and bottom boundary layers will be modelled by sheared Rayleigh-Bénard flows, the schematic view of which is shown in (b).

3 =/ 3 . ( 8 )Equation ( 8 )
function of the linearized dynamic equations is the natural framework to address the stability of open flows, as these BLs swept by the wind are.It leads to distinguish among the instabilities the convective ones that grow in space and time but are eventually swept out by the wind, and the absolute ones, the growth of which is vigorous enough to overcome this wind and propagate both upand downstream [21].Absolute instabilities are doubly relevant here.First, whereas convective instabilities are usually extrinsic and observed as the result of an exter-nal forcing, absolute instabilities are driven by the intrinsic dynamics of the flow and are more robust.Then, whereas convective instabilities in sheared Rayleigh-Bénard convection retrieve the critical Rayleigh number Ra 0 c ≈ 1708, unaffected by the shear flow and independent of the Prandtl number, it is established hereinafter and shown in Fig. 2 that the critical Rayleigh number of the absolute instabilities (Ra abs c ) substantially varies with Pr and Re bl .So, following the idea of Malkus and Howard, we still suppose that the thermal BLs thickness (δ T ) adjusts itself so as to set Ra bl at the threshold of instability, but we further claim that the relevant threshold pertains to the emergence of absolute instabilities: Ra bl = Ra abs c (Re bl , Pr). (6) Combining (2) and (3) fixes the thickness δ T , and the condition of absolute instability (6) becomes Ra abs c (Re bl , Pr) Re bl Ra (Re w F ) 3 .(7) Once Ra abs c as a function of Re bl and Pr has been computed, as elaborated in §1, solving equation (7) yields Re bl as a function of the three parameters Ra, Pr and Re w , and equations (5) and (6) yield the Nusselt number: Nu = Ra Ra abs c Re bl (Ra, Pr, Re w ) , Pr 1expresses Nu as a function of the 3 parameters Ra, Pr and Re w .Whereas Ra and Pr are external control parameters, Re w results from the dynamics of the flow.Computing Nu from this model thus requires to plug in a value for Re w .Strictly speaking, this model cannot predict the Nusselt number out of the external parameters.It is nonetheless possible to check the validity of the basic assumption of our interpretation:

Figure 3 :
Figure 3: Reynolds (a) and Nusselt (b) numbers in compensated form, a functions of Ra, for Pr = 1.Symbols: DNS from respective references.

3
(a) are consistent with Re r.m.s.behaving as the squareroot of Ra: Re r.m.s.= b Ra 1/2 ,

Figure. 4
Figure.4(b) then shows that the present model also retrieves with a good accuracy the variations of Nu with both Ra and Pr, using again F = 1, ξ = 0.22 and Re r.m.s.