Bounds for Rayleigh-B\'enard convection between free-slip boundaries with an imposed heat flux

We prove the first rigorous bound on the heat transfer for three-dimensional Rayleigh-B\'enard convection of finite-Prandtl-number fluids between free-slip boundaries with an imposed heat flux. Using the auxiliary functional method with a quadratic functional, which is equivalent to the background method, we prove that the Nusselt number $Nu$ is bounded by $Nu \leq 0.5999 R^{1/3}$ uniformly in the Prandtl number, where $R$ is the Rayleigh number based on the imposed heat flux. In terms of the Rayleigh number based on the mean vertical temperature drop, $Ra$, we obtain $Nu \leq 0.4646 Ra^{1/2}$. The scaling with Rayleigh number is the same as that of bounds obtained with no-slip isothermal, free-slip isothermal, and no-slip fixed flux boundaries, and numerical optimization of the bound suggests that it cannot be improved within our bounding framework. Contrary to the two-dimensional case, therefore, the $Ra$-dependence of rigorous upper bounds on the heat transfer obtained with the background method for three-dimensional Rayleigh-B\'enard convection is insensitive to both the thermal and the velocity boundary conditions.


Introduction
Rayleigh-Bénard (RB) convection, the buoyancy-driven motion of a fluid confined between horizontal plates, is a cornerstone of fluid mechanics. Its applications include atmospheric and oceanic physics, astrophysics, and industrial engineering (Lappa 2010, Chapter 3), and due to its rich dynamics it has also become a paradigm to investigate pattern formation and nonlinear phenomena (Ahlers et al. 2009).
One of the fundamental questions in the study of convection is to which extent the flow enhances the transport of heat across the layer. Precisely, one would like to relate the Nusselt number Nu (the nondimensional measure of the heat transfer enhancement) to the parameters of the fluid and the strength of the thermal forcing. These are described, respectively, by the Prandtl and Rayleigh numbers Pr = ν/κ and Ra = αgh 3 ∆/(νκ), where α is the fluid's thermal expansion coefficient, ν is its kinematic viscosity, κ is its thermal diffusivity, h is the dimensional height of the layer, g is the gravitational acceleration, and ∆ is the average temperature drop across the layer. It is generally expected that for large Rayleigh numbers the Nusselt number obeys a simple scaling law of the form Nu ∼ Pr a Ra b . However, different phenomenological arguments predict different scaling exponents in the ranges −1/4 a 1/2 and 2/7 b 1/2 (see Tables  I and II in Ahlers et al. 2009), and the available experimental evidence in the high-Ra regime is controversial (Ahlers et al. 2009).
Discrepancies in the measurements are often attributed to differences in the boundary conditions (BCs) or in the Prandtl number. From the modelling point of view, eight basic configurations of RB convection can be identified depending on the Prandtl number (finite or infinite), the BCs for the fluid's temperature (fixed temperature or fixed flux), and the BCs for its velocity (no-slip or free-slip). Two-dimensional simulations (Johnston & Doering 2009;Goluskin et al. 2014;van der Poel et al. 2014) have shown that changing the thermal BCs for given velocity BCs has no quantitative effect on Nu, while replacing no-slip boundaries with free-slip ones can dramatically reduce the heat transfer through the appearance of zonal flows. However, zonal flows have not been observed in three dimensions (van der Poel et al. 2014) and how different BCs affect the Nu-Ra-Pr relationship in general remains an open problem.
In the absence of extensive numerical result for the high-Ra regime in three dimensions, one way to make progress is through rigorous analysis of the equations that ostensibly describe RB convection. A particularly fruitful approach is to use the background method (Doering & Constantin 1992, 1994, 1996Constantin & Doering 1995) and derive rigorous bounds of the form Nu f (Ra, Pr ) for each of the eight configurations described above.
In contrast, the only bounds available for free-slip velocity BCs are for RB convection between isothermal plates. All identities and estimates used in the no-slip analysis of Doering & Constantin (1996) hold also for free-slip boundaries, so one immediately obtains Nu Ra 1/2 at finite Pr . This result can be tightened to Nu Ra 5/12 in two dimensions and at infinite Pr in three dimensions by explicitly taking advantage of both the stress-free and the isothermal BCs (Whitehead & Doering 2011.
Free-slip conditions pose a challenge for the background method when a constant heat flux κβ, rather than a fixed boundary temperature, is imposed. The reason is that the analysis usually relies on at least one of the temperature and horizontal velocities being fixed at the top and bottom boundaries, which is not the case with free-slip and fixed flux BCs. In this short paper we show that such lack of "boundary control" for the dynamical fields can be overcome with a simple symmetry argument and thereby prove the first rigorous upper bound on Nu for RB convection between free-slip boundaries with imposed heat flux.
The exposition is organised as follows. Section 2 reviews the Boussinesq equations used to model the system. We formulate a bounding principle for Nu in §3, and prove our main result in §4. Finally, §5 offers further discussion and conclusive remarks.

The model
We model the system using the Boussinesq equations and make all variables nondimensional using h, h/κ, and hβ, respectively, as the length, time and temperature scales (Otero et al. 2002). The nondimensional velocity u(x, y, z, t), pressure p(x, y, z, t), and perturbations θ(x, y, z, t) from the conductive temperature profile T c = −z then satisfy (Otero et al. 2002) where e z is the unit vector in the z direction and R = αgβh 4 /(νκ) is the Rayleigh number based on the imposed boundary heat flux. Note that R is related to the Rayleigh number based on the (unknown) mean temperature drop, Ra, by R = Ra Nu (Otero et al. 2002). The domain is periodic in the horizontal (x, y) directions and the vertical BCs are Since the average vertical heat flux across the layer is fixed to 1 in nondimensional units, convection reduces the mean temperature difference between the top and bottom plates and hence the mean conductive heat flux −∂ z T = 1 − ∂ z θ (here and throughout this work overlines denote averages over infinite time, while angle brackets denote volume averages). The Nusselt number-the ratio of the average vertical heat flux and the mean conductive flux-is then given by (see also Otero et al. 2002) (2.3)

Upper bound formulation
When R < 120 conduction is globally asymptotically stable and Nu = 1 (Chapman & Proctor 1980;Goluskin 2016b). For R > 120 convection sets in (Hurle et al. 1967) and we look for a positive lower bound L on 1 − ∂ z θ , implying Nu 1/L. To find L we use the background method (Doering & Constantin 1994, 1996Constantin & Doering 1995) but we formulate it in the language of the auxiliary functional method (Chernyshenko et al. 2014;Chernyshenko 2017) because of its conceptual simplicity: it relies on one simple inequality, rather than a seemingly ad hoc manipulation of the governing equations.
The analysis starts with the observation that any uniformly bounded and differentiable time-dependent functional V(t) = V{θ(·, t), u(·, t)} satisfies dV/dt = 0. Consequently, to prove that 1 − ∂ z θ L it suffices to show that at any instant in time Using the ideas outlined by Chernyshenko (2017), it can be shown that constructing a background temperature field in the "classical" background method analysis is equivalent to finding constants a, b and L and a function ϕ(z) such that (3.1) holds for We assume that u and θ are sufficiently regular in time to ensure differentiability of this functional, while uniform boundedness can be proven using estimates similar to those presented in this paper (we do not give a full proof in this work due to space limitations, but outline the argument in appendix A). The functional S{θ(·, t), u(·, t)} corresponding to (3.2) can be expressed in terms of u and θ using (2.1a)-(2.1c). Integrating the volume average u·(2.1a) by parts using incompressibility and the BCs yields d dt Averaging θ×(2.1c) and ϕ×(2.1c) in a similar way gives d dt Combining expressions (3.3)-(3.5) and rearranging we find To prove that (3.1) holds at all times, we make one key further simplification: we drop the equation of motions and choose a, b, L and ϕ(z) such that S{θ, u} 0 for any time-independent fields θ = θ(x, y, z) and u = u(x, y, z) that satisfy (2.1b) and the BCs. Hereafter, we also assume that a, b > 0 to ensure that S{θ, u} is bounded below.
Incompressibility can be incorporated explicitly in (3.6) upon substitution of the horizontal Fourier expansions is the horizontal position vector and k = (k x , k y ) is the wavevector. The z-dependent Fourier amplitudes u k , θ k satisfy the same vertical BCs as the full fields in (2.2). Using the Fourier-transformed incompressibility constraint one can show that (Doering & Constantin 1996;Otero et al. 2002) In these equations and in the following we write k 2 = k 2 x + k 2 y , · denotes the standard Lebesgue L 2 norm on the interval (0, 1), andw k is the complex conjugate of w k .
The right-hand side of (3.8) is clearly non-negative if S 0 0 and S k 0 for all wavevectors k = (0, 0). (A standard argument based on the consideration of fields θ and u with a single Fourier mode shows that these conditions are also necessary, so enforcing the positivity of each S k individually does not introduce conservativeness. However, necessity is not required to proceed with our argument so we omit the details for brevity.) In particular, given a, b, and ϕ the largest value of L for which S 0 0 is found upon completing the square (in the L 2 norm sense) in (3.9a), so we set . (3.10) We will try to maximise this expression over a, b and ϕ subject to the non-negativity of the functional S k in (3.9b) for all wavevectors k = (0, 0). Note that S k and the righthand side of (3.10) reduce, respectively, to the quadratic form and the bound obtained by Otero et al. (2002) using the "classical" background method analysis if we let a = b−1 and identify [ϕ (z) − 2b + 1]/(2b) with the derivative of the background temperature field. We also remark that our analysis appears more general because the choice a = 1 − b is unjustified at this stage, but its optimality (at least within the context of our proof) will be demonstrated below.

An explicit bound
Let δ 1/2 and consider the piecewise-linear profile ϕ(z) shown in figure 1, whose derivative is (4.1) To show that a, b, and δ can be chosen to make the quadratic form S k {θ k , w k } in (3.9b) positive semidefinite we rewrite θ k and w k as the sum of functions that are symmetric and antisymmetric with respect to z = 1/2. In other words, we decompose (4.3a,b) (The subscripts + and − denote, respectively, the symmetric and antisymmetric parts.) Since ϕ (z) is symmetric with respect to z = 1/2 by construction we obtain i.e., we can split the quadratic form S k {θ k , w k } into its symmetric and antisymmetric components also. Symmetric and antisymmetric Fourier amplitudes θ k and w k -for which one term on the right-hand side of (4.4) vanishes-are also admissible, so S k is nonnegative if and only if it is so for arguments that are either symmetric or antisymmetric with respect to z = 1/2 (and, of course, satisfy the correct BCs). As before, the "only if" statement is not needed to proceed but guarantees that no conservativeness is introduced. The decomposition into symmetric and antisymmetric components is the essential ingredient of our proof. In fact, contrary to the case of no-slip boundaries considered by Otero et al. (2002), the free-slip and fixed flux BCs cannot be used to control the indefinite term in S k {θ k , w k } via the usual elementary functional-analytic estimates. However, θ ± and w ± (or the appropriate derivatives) are known not only at the boundaries, but also on the symmetry plane. In particular, for small δ the indefinite term in (3.9b) can be controlled without recourse to the free-slip, fixed-flux BCs using w ± (0) = w + (1/2) = θ − (1/2) = 0. (4.5) To prove this, recall that for any symmetric or antisymmetric quantity q(z) (4.6) Symmetry, (4.1), and the identity |θ ±w± | = |θ ± w ± | (w ± is the complex conjugate of w ± ) yield with M := (1 + a + b)/b. Since w ± (0) = 0 the product θ ± w ± vanishes at z = 0 and for any z δ 1/2 the fundamental theorem of calculus implies Using the fact that w ± (0) = 0 once again, the fundamental theorem of calculus for ξ 1/2, the Cauchy-Schwarz inequality, and (4.6) we also obtain Furthermore, the conditions in (4.5) imply that the product θ ± w ± vanishes at the symmetry plane, so similar estimates as above yield (4.10) Upon inserting (4.9) and (4.10) into (4.8), applying the Cauchy-Schwarz inequality, and using (4.6) we arrive at (4.11) Substituting this estimate into (4.7) and integrating gives an estimate for the indefinite term in (3.9b), and after dropping the term ak 2 w ± 2 /(bR) we conclude that Recalling the definition of M and that a quadratic form αu 2 + βuv + γv 2 is positive semidefinite if β 2 4αγ, the right-hand side of (4.12) is non-negative if we set Having chosen δ to ensure the non-negativity of S k , all is left to do is optimise the eventual bound Nu L −1 over a and b as a function of R. Substituting (4.1) into (3.10) for our choice of δ yields (4.14) In order to maximise this expression with respect to a, b > 0 we set the partial derivatives ∂L/∂a and ∂L/∂b to zero. After some rearrangement it can be verified that ∂L ∂a = 0 ⇔ Ab 1/4 (7a + b + 1) − 4R 1/4 a 3/4 (1 + a + b) 1/2 = 0, (4.15a) A few lines of simple algebra show that setting to zero the second factor in (4.15b) leads to a solution with negative a or b, so we must choose b = 1 + a where a > 0 satisfies A 4 (1 + 4a) 4 − 64R(1 + a)a 3 = 0. (4.16) (No positive roots exist if R 4A 4 ≈ 43.92, but we are only interested in R 120 because conduction is globally asymptotically stable otherwise. It can also be checked that this stationary point is a maximum; the algebra is straightforward but lengthy and uninteresting, so we do not report it for brevity). In particular, when R tends to infinity (4.16) admits an asymptotic solution of the form a = a 1 R −1/3 + O(R −2/3 ). Substituting this expansion into (4.16) and solving for the leading order terms gives a 1 = A 4/3 /4. We then set b = 1 + a and a = a 1 R −1/3 in (4.14), simplify, and estimate Note that this bound is sharp as R → ∞. Consequently, Recalling that R = Nu Ra (Otero et al. 2002) we can also express this bound in terms of the Rayleigh number Ra based on the average temperature drop across the layer:

Discussion
The bound proven in this work is the first rigorous result for three-dimensional RB convection between free-slip, fixed flux boundaries (but note that our proof holds also in the two-dimensional case). Key to the result is a symmetry argument that overcomes the loss of boundary control for the trial fields when the no-slip velocity conditions are replaced with free-slip ones. Our approach is fully equivalent to the "classical" application of the background method to the temperature field, and the scaling of our bound with Ra is the same as obtained for no-slip BCs (irrespectively of the thermal BCs, see Doering & Constantin 1996;Otero et al. 2002) and for free-slip isothermal BCs (Doering & Constantin 1996). Modulo differences in the prefactor, therefore, rigorous upper bounds on the the heat transfer obtained with the background method for three-dimensional RB convection at finite Pr are insensitive to both the velocity and the thermal BCs.
Whether convective flows observed in reality exhibit the same lack of sensitivity to the BCs, however, remains uncertain. Two-dimensional simulations indicate that the thermal BCs make no quantitative difference for given velocity BCs (Johnston & Doering 2009;Goluskin et al. 2014), while replacing no-slip with free-slip leads to zonal flows with reduced vertical heat transfer (Goluskin et al. 2014;van der Poel et al. 2014). Partial support for such observations comes from the improved bound Nu Ra 5/12 obtained with free-slip isothermal boundaries in two dimensions (Whitehead & Doering 2011). It does not seem unreasonable to expect that a symmetry argument similar to that of this paper will extend the result to the fixed flux case, but we leave a formal confirmation to future work. On the other hand, zonal flows have not been observed in three dimensions (van der Poel et al. 2014). More extensive three-dimensional numerical simulations should be carried out to reveal if and how free-slip conditions affect the Nu-Ra relationship, as well as whether the thermal BCs can have any influence. Should numerical simulations in three dimensions suggest that Nu grows more slowly than Ra 1/2 , the challenge will be to improve the scaling exponents in (4.18)-(4.19). The argument by Whitehead & Doering (2012) may be adapted to study the infinite-Pr limit, but cannot be used at finite Pr . Moreover, at finite Pr it does not seem sufficient to consider a more sophisticated choice of a, b, and ϕ(z) in the functional (3.2). To provide evidence of this fact, we used quinopt (Fantuzzi et al. 2017a,b) to maximise the constant L (and, consequently, minimise the eventual bound Nu L −1 ) over all constants a, b and functions ϕ(z) that make the functionals in (3.9a) and (3.9b) positive semidefinite. We considered domains with period 2π and 10π in both horizontal directions, respectively, and the corresponding optimal bounds on Nu are compared to the analytic bound (4.18) in figure 2(a). For both values of the horizontal period a least-square power-law fit to the numerical results for R 10 6 returns L −1 ≈ 0.325 R 0.33 . Moreover, as illustrated in figure 2(b) for R = 10 5 , the optimal ϕ(z) closely resembles the analytical profile sketched in figure 1: it is approximately linear with slope a + b in the bulk and it decreases near the top and bottom boundaries. This strongly suggests that carefully tuning a, b, and ϕ(z) can only improve the prefactor in (4.18).
Lowering the scaling exponent for three-dimensional RB convection at finite Prandtl number, if at all possible, will therefore demand a different approach. Recently,  have proven that the auxiliary functional method gives arbitrarily sharp bounds on maximal time averages for systems governed by ordinary differential equations. This gives hope that progress may be achieved in the context of RB convection if a more general functional than (3.2) is considered. The resulting bounding problem will inevitably be harder to tackle with purely analytical techniques, but the viability of this approach may be assessed with computer-assisted investigation based on sum-of-squares programming (see, e.g., Goulart & Chernyshenko 2012;Chernyshenko et al. 2014;Fantuzzi et al. 2016;Goluskin 2016a). Another option is to try and lower the bound proven here through the study of optimal "wall-to-wall" transport problems (Hassanzadeh et al. 2014;. Exactly how much these alternative bounding techniques can improve on the background method and advance our ability to derive a rigorous quantitative description of hydrodynamic systems is the subject of ongoing research. We are indebted to D. Goluskin and J. P. Whitehead, who introduced us to the problem studied in this paper. We thank them, C. R. Doering, A. Wynn, and S. I. Chernyshenko for their encouragement and helpful comments. Funding by an EPSRC scholarship (award ref. 1864077) and the support and hospitality of the Geophysical Fluid Dynamics program at Woods Hole Oceanographic Institution are gratefully acknowledged.