Bounds for Rayleigh–Bénard convection between free-slip boundaries with an imposed heat flux

Abstract

We prove the first rigorous bound on the heat transfer for three-dimensional Rayleigh–Bénard 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 $\mathit{Nu}$ is bounded by $\mathit{Nu}\leqslant 0.5999\mathit{R}^{1/3}$ uniformly in the Prandtl number, where $\mathit{R}$ is the Rayleigh number based on the imposed heat flux. In terms of the Rayleigh number based on the mean vertical temperature drop, $\mathit{Ra}$, we obtain $\mathit{Nu}\leqslant 0.4646\mathit{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 optimisation of the bound suggests that it cannot be improved within our bounding framework. Contrary to the two-dimensional case, therefore, the $\mathit{Ra}$-dependence of rigorous upper bounds on the heat transfer obtained with the background method for three-dimensional Rayleigh–Bénard convection is insensitive to both the thermal and the velocity boundary conditions.

1 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, chap. 3), and due to its rich dynamics it has also become a paradigm to investigate pattern formation and nonlinear phenomena (Ahlers, Grossmann & Lohse 2009).

One of the fundamental questions in the study of convection is to what extent the flow enhances the transport of heat across the layer. Precisely, one would like to relate the Nusselt number $\mathit{Nu}$ (the non-dimensional 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 $\mathit{Pr}=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$ and $\mathit{Ra}=\unicode[STIX]{x1D6FC}gh^{3}\unicode[STIX]{x1D6E5}/(\unicode[STIX]{x1D708}\unicode[STIX]{x1D705})$ , where $\unicode[STIX]{x1D6FC}$ is the fluid’s thermal expansion coefficient, $\unicode[STIX]{x1D708}$ is its kinematic viscosity, $\unicode[STIX]{x1D705}$ is its thermal diffusivity, $h$ is the dimensional height of the layer, $g$ is the gravitational acceleration, and $\unicode[STIX]{x1D6E5}$ 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 $\mathit{Nu}\sim \mathit{Pr}^{a}\mathit{Ra}^{b}$ . However, different phenomenological arguments predict different scaling exponents in the ranges $-1/4\leqslant a\leqslant 1/2$ and $2/7\leqslant b\leqslant 1/2$ (see Tables I and II in Ahlers et al. 2009), and the available experimental evidence in the high- $\mathit{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 $\mathit{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 $\mathit{Nu}$ – $\mathit{Ra}$ – $\mathit{Pr}$ relationship in general remains an open problem.

In the absence of extensive numerical results for the high- $\mathit{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, 1996; Constantin & Doering 1995) and derive rigorous bounds of the form $\mathit{Nu}\leqslant f(\mathit{Ra},\mathit{Pr})$ for each of the eight configurations described above.

The no-slip case has been studied extensively. For fluids with finite Prandtl number the bound $\mathit{Nu}\lesssim \mathit{Ra}^{1/2}$ holds uniformly in $\mathit{Pr}$ irrespective of the thermal BCs (Doering & Constantin 1996; Otero et al. 2002; Wittenberg 2010; Wittenberg & Gao 2010). When $\mathit{Pr}=\infty$ (and $\mathit{Pr}\gtrsim \mathit{Ra}^{1/3}$ with isothermal boundaries), instead, one has $\mathit{Nu}\lesssim \ell (\mathit{Ra})\mathit{Ra}^{1/3}$ , where $\ell (\mathit{Ra})$ is a logarithmic correction whose exact form depends on the thermal BCs (Otto & Seis 2011; Whitehead & Wittenberg 2014; Choffrut, Nobili & Otto 2016).

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 $\mathit{Nu}\leqslant \mathit{Ra}^{1/2}$ at finite $\mathit{Pr}$ . This result can be tightened to $\mathit{Nu}\leqslant \mathit{Ra}^{5/12}$ in two dimensions and at infinite $\mathit{Pr}$ in three dimensions by explicitly taking advantage of both the stress-free and the isothermal BCs (Whitehead & Doering 2011, 2012).

Free-slip conditions pose a challenge for the background method when a constant heat flux $\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FD}$ , 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 $\mathit{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 $\mathit{Nu}$ in § 3, and prove our main result in § 4. Finally, § 5 offers further discussion and conclusive remarks.

2 The model

We model the system using the Boussinesq equations and make all variables non-dimensional using $h$ , $h/\unicode[STIX]{x1D705}$ , and $h\unicode[STIX]{x1D6FD}$ , respectively, as the length, time and temperature scales (Otero et al. 2002). The non-dimensional velocity $\boldsymbol{u}(x,y,z,t)$ , pressure $p(x,y,z,t)$ , and perturbations $\unicode[STIX]{x1D703}(x,y,z,t)$ from the conductive temperature profile $T_{c}=-z$ then satisfy (Otero et al. 2002)

(2.1a ) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{t}\boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}+\unicode[STIX]{x1D735}p=\mathit{Pr}\,\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\mathit{Pr}\,\mathit{R}\,(\unicode[STIX]{x1D703}-z)\boldsymbol{e}_{z}, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.1c ) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D703}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}+w, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{e}_{z}$ is the unit vector in the $z$ direction and $R=\unicode[STIX]{x1D6FC}g\unicode[STIX]{x1D6FD}h^{4}/(\unicode[STIX]{x1D708}\unicode[STIX]{x1D705})$ is the Rayleigh number based on the imposed boundary heat flux. Note that $\mathit{R}$ is related to the Rayleigh number based on the (unknown) mean temperature drop, $\mathit{Ra}$ , by $\mathit{R}=\mathit{Ra}\,\mathit{Nu}$ (Otero et al. 2002). The domain is periodic in the horizontal ( $x$ , $y$ ) directions and the vertical BCs are

(2.2a,b ) $$\begin{eqnarray}\unicode[STIX]{x2202}_{z}u=\unicode[STIX]{x2202}_{z}v=w=0,\quad \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}=0\quad \text{at}~z=0~\text{and}~z=1.\end{eqnarray}$$

Since the average vertical heat flux across the layer is fixed to 1 in non-dimensional units, convection reduces the mean temperature difference between the top and bottom plates and hence the mean conductive heat flux $\overline{\langle -\unicode[STIX]{x2202}_{z}T\rangle }=\overline{1-\langle \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}\rangle }$ (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) $$\begin{eqnarray}\mathit{Nu}=\left(\overline{1-\langle \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}\rangle }\right)^{-1}.\end{eqnarray}$$

3 Upper bound formulation

When $\mathit{R}<120$ , conduction is globally asymptotically stable and $\mathit{Nu}=1$ (Chapman & Proctor 1980; Goluskin 2016). For $\mathit{R}>120$ , convection sets in (Hurle, Jakeman & Pike 1967) and we look for a positive lower bound $L$ on $\overline{1-\langle \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}\rangle }$ , implying $\mathit{Nu}\leqslant 1/L$ . To find $L$ we use the background method (Doering & Constantin 1994; Constantin & Doering 1995; Doering & Constantin 1996) 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 ${\mathcal{V}}(t)={\mathcal{V}}\{\unicode[STIX]{x1D703}(\cdot ,t),\boldsymbol{u}(\boldsymbol{\cdot },t)\}$ satisfies $\overline{\text{d}{\mathcal{V}}/\text{d}t}=0$ . Consequently, to prove that $\overline{1-\langle \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}\rangle }\geqslant L$ it suffices to show that at any instant in time

(3.1) $$\begin{eqnarray}{\mathcal{S}}\{\unicode[STIX]{x1D703}(\cdot ,t),\boldsymbol{u}(\boldsymbol{\cdot },t)\}:=\frac{\text{d}{\mathcal{V}}}{\text{d}t}+1-\langle \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}\rangle -L\geqslant 0.\end{eqnarray}$$

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 $\unicode[STIX]{x1D711}(z)$ such that (3.1) holds for

(3.2) $$\begin{eqnarray}{\mathcal{V}}\{\unicode[STIX]{x1D703}(\cdot ,t),\boldsymbol{u}(\boldsymbol{\cdot },t)\}:=-\frac{a}{2\mathit{Pr}\,\mathit{R}}\langle \left|\boldsymbol{u}\right|^{2}\rangle -\frac{b}{2}\langle \unicode[STIX]{x1D703}^{2}\rangle +\langle \unicode[STIX]{x1D711}\unicode[STIX]{x1D703}\rangle .\end{eqnarray}$$

We assume that $\boldsymbol{u}$ and $\unicode[STIX]{x1D703}$ are sufficiently regular in time to ensure differentiability of this functional, while uniform boundedness can be proved 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 ${\mathcal{S}}\{\unicode[STIX]{x1D703}(\cdot ,t),\boldsymbol{u}(\boldsymbol{\cdot },t)\}$ corresponding to (3.2) can be expressed in terms of $\boldsymbol{u}$ and $\unicode[STIX]{x1D703}$ using (2.1a )–(2.1c ). Integrating the volume average $\langle \boldsymbol{u}\boldsymbol{\cdot }(2.1a)\rangle$ by parts using incompressibility and the BCs yields

(3.3) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\frac{\langle \left|\boldsymbol{u}\right|^{2}\rangle }{2\mathit{Pr}\,\mathit{R}}=-\frac{\langle \left|\unicode[STIX]{x1D735}\boldsymbol{u}\right|^{2}\rangle }{\mathit{R}}+\langle w\unicode[STIX]{x1D703}\rangle .\end{eqnarray}$$

Averaging $\unicode[STIX]{x1D703}\times$  (2.1c ) and $\unicode[STIX]{x1D711}\times$  (2.1c ) in a similar way gives

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\frac{\langle \unicode[STIX]{x1D703}^{2}\rangle }{2}=-\langle \left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\right|^{2}\rangle +\langle w\unicode[STIX]{x1D703}\rangle , & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \unicode[STIX]{x1D711}\,\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D703}\rangle =\langle \unicode[STIX]{x1D711}^{\prime }w\unicode[STIX]{x1D703}\rangle -\langle \unicode[STIX]{x1D711}^{\prime }\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}\rangle . & \displaystyle\end{eqnarray}$$

Combining expressions (3.3)–(3.5) and rearranging we find

(3.6) $$\begin{eqnarray}{\mathcal{S}}\{\unicode[STIX]{x1D703}(\cdot ,t),\boldsymbol{u}(\boldsymbol{\cdot },t)\}=1-L-\langle (\unicode[STIX]{x1D711}^{\prime }+1)\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}\rangle +\left\langle \frac{a}{\mathit{R}}\left|\unicode[STIX]{x1D735}\boldsymbol{u}\right|^{2}+b\left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\right|^{2}+(\unicode[STIX]{x1D711}^{\prime }-a-b)w\unicode[STIX]{x1D703}\right\rangle .\end{eqnarray}$$

To prove that (3.1) holds at all times, we make one key further simplification: we drop the equations of motion and choose $a$ , $b$ , $L$ and $\unicode[STIX]{x1D711}(z)$ such that ${\mathcal{S}}\{\unicode[STIX]{x1D703},\boldsymbol{u}\}\geqslant 0$ for any time-independent fields $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}(x,y,z)$ and $\boldsymbol{u}=\boldsymbol{u}(x,y,z)$ that satisfy (2.1b ) and the BCs. Hereafter, we also assume that $a,b>0$ to ensure that ${\mathcal{S}}\{\unicode[STIX]{x1D703},\boldsymbol{u}\}$ is bounded below.

Incompressibility can be incorporated explicitly in (3.6) upon substitution of the horizontal Fourier expansions

(3.7a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D703}=\mathop{\sum }_{\boldsymbol{k}}\unicode[STIX]{x1D703}_{\boldsymbol{k}}(z)\text{e}^{\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}},\quad \boldsymbol{u}=\mathop{\sum }_{\boldsymbol{k}}\boldsymbol{u}_{\boldsymbol{k}}(z)\text{e}^{\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}},\end{eqnarray}$$

where $\boldsymbol{x}=(x,y)$ is the horizontal position vector and $\boldsymbol{k}=(k_{x},k_{y})$ is the wavevector. The $z$ -dependent Fourier amplitudes $\boldsymbol{u}_{\boldsymbol{k}}$ , $\unicode[STIX]{x1D703}_{\boldsymbol{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)

(3.8) $$\begin{eqnarray}{\mathcal{S}}\{\unicode[STIX]{x1D703},\boldsymbol{u}\}\geqslant {\mathcal{S}}_{0}\{\unicode[STIX]{x1D703}_{0}\}+b\mathop{\sum }_{\boldsymbol{k}\neq (0,0)}{\mathcal{S}}_{\boldsymbol{k}}\{\unicode[STIX]{x1D703}_{\boldsymbol{k}},w_{\boldsymbol{k}}\},\end{eqnarray}$$

with

(3.9a ) $$\begin{eqnarray}{\mathcal{S}}_{0}\{\unicode[STIX]{x1D703}_{0}\}:=b\Vert \unicode[STIX]{x1D703}_{0}^{\prime }\Vert ^{2}-\int _{0}^{1}(\unicode[STIX]{x1D711}^{\prime }+1)\unicode[STIX]{x1D703}_{0}^{\prime }\,\text{d}z+1-L,\end{eqnarray}$$

(3.9b ) $$\begin{eqnarray}\displaystyle {\mathcal{S}}_{\boldsymbol{k}}\{\unicode[STIX]{x1D703}_{\boldsymbol{k}},w_{\boldsymbol{k}}\} & := & \displaystyle \Vert \unicode[STIX]{x1D703}_{\boldsymbol{k}}^{\prime }\Vert ^{2}+k^{2}\Vert \unicode[STIX]{x1D703}_{\boldsymbol{k}}\Vert ^{2}+\frac{a}{b\mathit{R}}\left(\frac{\Vert w_{\boldsymbol{k}}^{\prime \prime }\Vert ^{2}}{k^{2}}+2\Vert w_{\boldsymbol{k}}^{\prime }\Vert ^{2}+k^{2}\Vert w_{\boldsymbol{k}}\Vert ^{2}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\int _{0}^{1}\frac{\unicode[STIX]{x1D711}^{\prime }-a-b}{b}\,\text{Re}(\unicode[STIX]{x1D703}_{\boldsymbol{k}}\tilde{w}_{\boldsymbol{k}})\,\text{d}z.\end{eqnarray}$$

In these equations and in the following we write $k^{2}=k_{x}^{2}+k_{y}^{2}$ , $\Vert \boldsymbol{\cdot }\Vert$ denotes the standard Lebesgue ${\mathcal{L}}^{2}$ norm on the interval $(0,1)$ , and $\tilde{w}_{\boldsymbol{k}}$ is the complex conjugate of $w_{\boldsymbol{k}}$ .

The right-hand side of (3.8) is clearly non-negative if ${\mathcal{S}}_{0}\geqslant 0$ and ${\mathcal{S}}_{\boldsymbol{k}}\geqslant 0$ for all wavevectors $\boldsymbol{k}\neq (0,0)$ . (A standard argument based on the consideration of fields $\unicode[STIX]{x1D703}$ and $\boldsymbol{u}$ with a single Fourier mode shows that these conditions are also necessary, so enforcing the positivity of each ${\mathcal{S}}_{\boldsymbol{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 $\unicode[STIX]{x1D711}$ the largest value of $L$ for which ${\mathcal{S}}_{0}\geqslant 0$ is found upon completing the square (in the ${\mathcal{L}}^{2}$ norm sense) in (3.9a ), so we set

(3.10) $$\begin{eqnarray}L=1-\frac{\Vert \unicode[STIX]{x1D711}^{\prime }+1\Vert ^{2}}{4b}.\end{eqnarray}$$

We will try to maximise this expression over $a$ , $b$ and $\unicode[STIX]{x1D711}$ subject to the non-negativity of the functional ${\mathcal{S}}_{\boldsymbol{k}}$ in (3.9b ) for all wavevectors $\boldsymbol{k}\neq (0,0)$ . Note that ${\mathcal{S}}_{\boldsymbol{k}}$ and the right-hand 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 $[\unicode[STIX]{x1D711}^{\prime }(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.

4 An explicit bound

Let $\unicode[STIX]{x1D6FF}\leqslant 1/2$ and consider the piecewise-linear profile $\unicode[STIX]{x1D711}(z)$ shown in figure 1, whose derivative is

(4.1) $$\begin{eqnarray}\unicode[STIX]{x1D711}^{\prime }(z)=\left\{\begin{array}{@{}ll@{}}\phantom{a}-1,\quad & z\in [0,\unicode[STIX]{x1D6FF}]\cup [1-\unicode[STIX]{x1D6FF},1],\\ a+b,\quad & z\in (\unicode[STIX]{x1D6FF},1-\unicode[STIX]{x1D6FF}).\end{array}\right.\end{eqnarray}$$

To show that $a$ , $b$ and $\unicode[STIX]{x1D6FF}$ can be chosen to make the quadratic form ${\mathcal{S}}_{\boldsymbol{k}}\{\unicode[STIX]{x1D703}_{\boldsymbol{k}},w_{\boldsymbol{k}}\}$ in (3.9b ) positive semidefinite we rewrite $\unicode[STIX]{x1D703}_{\boldsymbol{k}}$ and $w_{\boldsymbol{k}}$ as the sum of functions that are symmetric and antisymmetric with respect to $z=1/2$ . In other words, we decompose

(4.2a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{\boldsymbol{k}}(z)=\unicode[STIX]{x1D703}_{+}(z)+\unicode[STIX]{x1D703}_{-}(z),\quad w_{\boldsymbol{k}}(z)=w_{+}(z)+w_{-}(z),\end{eqnarray}$$

with

(4.3a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{\pm }(z)=\frac{\unicode[STIX]{x1D703}_{\boldsymbol{k}}(z)\pm \unicode[STIX]{x1D703}_{\boldsymbol{k}}(1-z)}{2},\quad w_{\pm }(z)=\frac{w_{\boldsymbol{k}}(z)\pm w_{\boldsymbol{k}}(1-z)}{2}.\end{eqnarray}$$

(The subscripts $+$ and $-$ denote, respectively, the symmetric and antisymmetric parts.) Since $\unicode[STIX]{x1D711}^{\prime }(z)$ is symmetric with respect to $z=1/2$ by construction we obtain

(4.4) $$\begin{eqnarray}{\mathcal{S}}_{\boldsymbol{k}}\{\unicode[STIX]{x1D703}_{\boldsymbol{k}},w_{\boldsymbol{k}}\}={\mathcal{S}}_{\boldsymbol{k}}\{\unicode[STIX]{x1D703}_{+},w_{+}\}+{\mathcal{S}}_{\boldsymbol{k}}\{\unicode[STIX]{x1D703}_{-},w_{-}\},\end{eqnarray}$$

i.e. we can split the quadratic form ${\mathcal{S}}_{\boldsymbol{k}}\{\unicode[STIX]{x1D703}_{\boldsymbol{k}},w_{\boldsymbol{k}}\}$ into its symmetric and antisymmetric components also. Symmetric and antisymmetric Fourier amplitudes $\unicode[STIX]{x1D703}_{\boldsymbol{k}}$ and $w_{\boldsymbol{k}}$ – for which one term on the right-hand side of (4.4) vanishes – are also admissible, so ${\mathcal{S}}_{\boldsymbol{k}}$ is non-negative 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.

Figure 1. Sketch of the piecewise-linear function $\unicode[STIX]{x1D711}(z)$ .

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 ${\mathcal{S}}_{\boldsymbol{k}}\{\unicode[STIX]{x1D703}_{\boldsymbol{k}},w_{\boldsymbol{k}}\}$ via the usual elementary functional-analytic estimates. However, $\unicode[STIX]{x1D703}_{\pm }$ and $w_{\pm }$ (or the appropriate derivatives) are known not only at the boundaries, but also on the symmetry plane. In particular, for small $\unicode[STIX]{x1D6FF}$ the indefinite term in (3.9b ) can be controlled without recourse to the free-slip, fixed-flux BCs using

(4.5) $$\begin{eqnarray}w_{\pm }(0)=w_{+}^{\prime }(1/2)=\unicode[STIX]{x1D703}_{-}(1/2)=0.\end{eqnarray}$$

To prove this, recall that for any symmetric or antisymmetric quantity $q(z)$

(4.6) $$\begin{eqnarray}\int _{0}^{1/2}\left|q(z)\right|^{2}\,\text{d}z=\frac{\Vert q\Vert ^{2}}{2}.\end{eqnarray}$$

Symmetry, equation (4.1), and the identity $\left|\unicode[STIX]{x1D703}_{\pm }\tilde{w}_{\pm }\right|=\left|\unicode[STIX]{x1D703}_{\pm }w_{\pm }\right|$ ( $\tilde{w}_{\pm }$ is the complex conjugate of $w_{\pm }$ ) yield

(4.7) $$\begin{eqnarray}\left|\int _{0}^{1}\frac{\unicode[STIX]{x1D711}^{\prime }-a-b}{b}\text{Re}(\unicode[STIX]{x1D703}_{\pm }\tilde{w}_{\pm })\,\text{d}z\right|\leqslant 2M\int _{0}^{\unicode[STIX]{x1D6FF}}\left|\unicode[STIX]{x1D703}_{\pm }w_{\pm }\right|\,\text{d}z,\end{eqnarray}$$

with $M:=(1+a+b)/b$ . Since $w_{\pm }(0)=0$ the product $\unicode[STIX]{x1D703}_{\pm }w_{\pm }$ vanishes at $z=0$ and for any $z\leqslant \unicode[STIX]{x1D6FF}\leqslant 1/2$ the fundamental theorem of calculus implies

(4.8) $$\begin{eqnarray}\left|\unicode[STIX]{x1D703}_{\pm }(z)w_{\pm }(z)\right|\leqslant \int _{0}^{z}\left|\unicode[STIX]{x1D703}_{\pm }(\unicode[STIX]{x1D709})\right||w_{\pm }^{\prime }(\unicode[STIX]{x1D709})|\,\text{d}\unicode[STIX]{x1D709}+\int _{0}^{z}|\unicode[STIX]{x1D703}_{\pm }^{\prime }(\unicode[STIX]{x1D709})|\left|w_{\pm }(\unicode[STIX]{x1D709})\right|\,\text{d}\unicode[STIX]{x1D709}.\end{eqnarray}$$

Using the fact that $w_{\pm }(0)=0$ once again, the fundamental theorem of calculus for $\unicode[STIX]{x1D709}\leqslant 1/2$ , the Cauchy–Schwarz inequality, and (4.6) we also obtain

(4.9) $$\begin{eqnarray}\left|w_{\pm }(\unicode[STIX]{x1D709})\right|=\left|\int _{0}^{\unicode[STIX]{x1D709}}w_{\pm }^{\prime }(\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702}\right|\leqslant \sqrt{\frac{\unicode[STIX]{x1D709}}{2}}\Vert w_{\pm }^{\prime }\Vert .\end{eqnarray}$$

Furthermore, the conditions in (4.5) imply that the product $\unicode[STIX]{x1D703}_{\pm }w_{\pm }^{\prime }$ vanishes at the symmetry plane, so similar estimates as above yield

(4.10) $$\begin{eqnarray}\displaystyle \left|\unicode[STIX]{x1D703}_{\pm }(\unicode[STIX]{x1D709})w_{\pm }^{\prime }(\unicode[STIX]{x1D709})\right| & = & \displaystyle \left|\int _{\unicode[STIX]{x1D709}}^{1/2}[\unicode[STIX]{x1D703}_{\pm }(\unicode[STIX]{x1D702})w_{\pm }^{\prime \prime }(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D703}_{\pm }^{\prime }(\unicode[STIX]{x1D702})w_{\pm }^{\prime }(\unicode[STIX]{x1D702})]\,\text{d}\unicode[STIX]{x1D702}\right|\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \frac{1}{2}\Vert \unicode[STIX]{x1D703}_{\pm }\Vert \Vert w_{\pm }^{\prime \prime }\Vert +\frac{1}{2}\Vert \unicode[STIX]{x1D703}_{\pm }^{\prime }\Vert \Vert w_{\pm }^{\prime }\Vert .\end{eqnarray}$$

Upon inserting (4.9) and (4.10) into (4.8), applying the Cauchy–Schwarz inequality, and using (4.6) we arrive at

(4.11) $$\begin{eqnarray}\left|\unicode[STIX]{x1D703}_{\pm }(z)w_{\pm }(z)\right|\leqslant \frac{z}{2}\left(\Vert \unicode[STIX]{x1D703}_{\pm }\Vert \Vert w_{\pm }^{\prime \prime }\Vert +\frac{1+\sqrt{2}}{\sqrt{2}}\Vert \unicode[STIX]{x1D703}_{\pm }^{\prime }\Vert \Vert w_{\pm }^{\prime }\Vert \right).\end{eqnarray}$$

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}\Vert w_{\pm }\Vert ^{2}/(b\mathit{R})$ we conclude that

(4.12) $$\begin{eqnarray}\displaystyle {\mathcal{S}}_{\boldsymbol{k}}\{\unicode[STIX]{x1D703}_{\pm },w_{\pm }\} & {\geqslant} & \displaystyle \frac{2a}{b\mathit{R}}\Vert w_{\pm }^{\prime }\Vert ^{2}-\frac{(1+\sqrt{2})M\unicode[STIX]{x1D6FF}^{2}}{2\sqrt{2}}\Vert w_{\pm }^{\prime }\Vert \Vert \unicode[STIX]{x1D703}_{\pm }^{\prime }\Vert +\Vert \unicode[STIX]{x1D703}_{\pm }^{\prime }\Vert ^{2}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{a}{b\mathit{R}k^{2}}\Vert w_{\pm }^{\prime \prime }\Vert ^{2}-\frac{M\unicode[STIX]{x1D6FF}^{2}}{2}\Vert w_{\pm }^{\prime \prime }\Vert \Vert \unicode[STIX]{x1D703}_{\pm }\Vert +k^{2}\Vert \unicode[STIX]{x1D703}_{\pm }\Vert ^{2}.\end{eqnarray}$$

Recalling the definition of $M$ and that a quadratic form $\unicode[STIX]{x1D6FC}u^{2}+\unicode[STIX]{x1D6FD}uv+\unicode[STIX]{x1D6FE}v^{2}$ is positive semidefinite if $\unicode[STIX]{x1D6FD}^{2}\leqslant 4\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FE}$ , the right-hand side of (4.12) is non-negative if we set

(4.13a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}=A\left[\frac{ab}{(1+a+b)^{2}\mathit{R}}\right]^{1/4},\quad A:=\left(\frac{8}{1+\sqrt{2}}\right)^{1/2}.\end{eqnarray}$$

Having chosen $\unicode[STIX]{x1D6FF}$ to ensure the non-negativity of ${\mathcal{S}}_{\boldsymbol{k}}$ , all that is left to do is optimise the eventual bound $\mathit{Nu}\leqslant L^{-1}$ over $a$ and $b$ as a function of $\mathit{R}$ . Substituting (4.1) into (3.10) for our choice of $\unicode[STIX]{x1D6FF}$ yields

(4.14) $$\begin{eqnarray}L=1-\frac{(1+a+b)^{2}}{4b}+\frac{A}{2}\frac{(1+a+b)^{3/2}a^{1/4}}{b^{3/4}R^{1/4}}.\end{eqnarray}$$

In order to maximise this expression with respect to $a,b>0$ we set the partial derivatives $\unicode[STIX]{x2202}L/\unicode[STIX]{x2202}a$ and $\unicode[STIX]{x2202}L/\unicode[STIX]{x2202}b$ to zero. After some rearrangement it can be verified that

(4.15a ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}L}{\unicode[STIX]{x2202}a} & = & \displaystyle 0\quad \Leftrightarrow \quad Ab^{1/4}(7a+b+1)-4R^{1/4}a^{3/4}(1+a+b)^{1/2}=0,\end{eqnarray}$$

(4.15b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}L}{\unicode[STIX]{x2202}b} & = & \displaystyle 0\quad \Leftrightarrow \quad (1+a-b)[2R^{1/4}(1+a+b)^{1/2}-3Aa^{1/4}b^{1/4}]=0.\end{eqnarray}$$

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

(4.16) $$\begin{eqnarray}A^{4}(1+4a)^{4}-64\mathit{R}(1+a)a^{3}=0.\end{eqnarray}$$

(No positive roots exist if $R\leqslant 4A^{4}\approx 43.92$ , but we are only interested in $\mathit{R}\geqslant 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 $\mathit{R}$ tends to infinity (4.16) admits an asymptotic solution of the form $a=a_{1}\mathit{R}^{-1/3}+O(\mathit{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}\mathit{R}^{-1/3}$ in (4.14), simplify, and estimate

(4.17) $$\begin{eqnarray}L=\frac{A^{4/3}}{4\mathit{R}^{1/3}}\left[\sqrt{2}\left(4+\frac{A^{4/3}}{\mathit{R}^{1/3}}\right)^{3/4}-1\right]\geqslant \frac{3A^{4/3}}{4\mathit{R}^{1/3}}.\end{eqnarray}$$

Note that this bound is sharp as $\mathit{R}\rightarrow \infty$ . Consequently,

(4.18) $$\begin{eqnarray}\mathit{Nu}\leqslant \frac{1}{L}\leqslant \frac{4\mathit{R}^{1/3}}{3A^{4/3}}\approx 0.5999\mathit{R}^{1/3}.\end{eqnarray}$$

Recalling that $\mathit{R}=\mathit{Nu}\,\mathit{Ra}$ (Otero et al. 2002) we can also express this bound in terms of the Rayleigh number $\mathit{Ra}$ based on the average temperature drop across the layer:

(4.19) $$\begin{eqnarray}\mathit{Nu}\leqslant \frac{8\mathit{Ra}^{1/2}}{3\sqrt{3}A^{2}}\approx 0.4646\,\mathit{Ra}^{1/2}.\end{eqnarray}$$

5 Discussion

The bound proved 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 $\mathit{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 heat transfer obtained with the background method for three-dimensional RB convection at finite $\mathit{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 $\mathit{Nu}\lesssim \mathit{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 $\mathit{Nu}$ – $\mathit{Ra}$ relationship, as well as whether the thermal BCs can have any influence.

Figure 2. (a) The analytic bound (4.18) (dashed line) versus the numerically optimal bounds for a $2\unicode[STIX]{x03C0}$ -periodic layer (solid line) and a $10\unicode[STIX]{x03C0}$ -periodic layer (dot-dashed line). The vertical dotted line at $R=120$ marks the global stability boundary for pure conduction. (b) The optimal $\unicode[STIX]{x1D711}(z)$ at $\mathit{R}=10^{5}$ . A dotted line with slope $a+b$ for the optimal values of these parameters at $\mathit{R}=10^{5}$ is shown for comparison with the analytical profile sketched in figure 1.

Should numerical simulations in three dimensions suggest that $\mathit{Nu}$ grows more slowly than $\mathit{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- $\mathit{Pr}$ limit, but cannot be used at finite $\mathit{Pr}$ . Moreover, at finite $\mathit{Pr}$ it does not seem sufficient to consider a more sophisticated choice of $a$ , $b$ and $\unicode[STIX]{x1D711}(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 $\mathit{Nu}\leqslant L^{-1}$ ) over all constants $a$ , $b$ and functions $\unicode[STIX]{x1D711}(z)$ that make the functionals in (3.9a ) and (3.9b ) positive semidefinite. We considered domains with period $2\unicode[STIX]{x03C0}$ and $10\unicode[STIX]{x03C0}$ in both horizontal directions, respectively, and the corresponding optimal bounds on $\mathit{Nu}$ are compared to the analytic bound (4.18) in figure 2(a). For both values of the horizontal period a least-squares power-law fit to the numerical results for $\mathit{R}\geqslant 10^{6}$ returns $L^{-1}\approx 0.325\,\mathit{R}^{0.33}$ . Moreover, as illustrated in figure 2(b) for $\mathit{R}=10^{5}$ , the optimal $\unicode[STIX]{x1D711}(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 $\unicode[STIX]{x1D711}(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, Tobasco, Goluskin & Doering (2017) have proved 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 investigations based on sum-of-squares programming (see, e.g. Goulart & Chernyshenko 2012; Chernyshenko et al. 2014; Fantuzzi et al. 2016; Goluskin 2017). Another option is to try and lower the bound proved here through the study of optimal ‘wall-to-wall’ transport problems (Hassanzadeh, Chini & Doering 2014; Tobasco & Doering 2017). 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.