Hostname: page-component-84b7d79bbc-g5fl4 Total loading time: 0 Render date: 2024-07-31T02:25:08.781Z Has data issue: false hasContentIssue false

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

Published online by Cambridge University Press:  05 January 2018

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.

Type
JFM Rapids
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© 2018 Cambridge University Press

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

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. Reference Ahlers, Grossmann and Lohse2009), and the available experimental evidence in the high- $\mathit{Ra}$ regime is controversial (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009).

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 Reference Johnston and Doering2009; Goluskin et al. Reference Goluskin, Johnston, Flierl and Spiegel2014; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Verzicco and Lohse2014) 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. Reference van der Poel, Ostilla-Mónico, Verzicco and Lohse2014) 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 Reference Doering and Constantin1992, Reference Doering and Constantin1994, Reference Doering and Constantin1996; Constantin & Doering Reference Constantin and Doering1995) 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 Reference Doering and Constantin1996; Otero et al. Reference Otero, Wittenberg, Worthing and Doering2002; Wittenberg Reference Wittenberg2010; Wittenberg & Gao Reference Wittenberg and Gao2010). 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 Reference Otto and Seis2011; Whitehead & Wittenberg Reference Whitehead and Wittenberg2014; Choffrut, Nobili & Otto Reference Choffrut, Nobili and Otto2016).

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 (Reference Doering and Constantin1996) 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 Reference Whitehead and Doering2011, Reference Whitehead and Doering2012).

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. Reference Otero, Wittenberg, Worthing and Doering2002). 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. Reference Otero, Wittenberg, Worthing and Doering2002)

(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. Reference Otero, Wittenberg, Worthing and Doering2002). 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. Reference Otero, Wittenberg, Worthing and Doering2002)

(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 Reference Chapman and Proctor1980; Goluskin Reference Goluskin2016). For $\mathit{R}>120$ , convection sets in (Hurle, Jakeman & Pike Reference Hurle, Jakeman and Pike1967) 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 Reference Doering and Constantin1994; Constantin & Doering Reference Constantin and Doering1995; Doering & Constantin Reference Doering and Constantin1996) but we formulate it in the language of the auxiliary functional method (Chernyshenko et al. Reference Chernyshenko, Goulart, Huang and Papachristodoulou2014; Chernyshenko Reference Chernyshenko2017) 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 (Reference Chernyshenko2017), 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 Reference Doering and Constantin1996; Otero et al. Reference Otero, Wittenberg, Worthing and Doering2002)

(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. (Reference Otero, Wittenberg, Worthing and Doering2002) 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. (Reference Otero, Wittenberg, Worthing and Doering2002), 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. Reference Otero, Wittenberg, Worthing and Doering2002) 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 Reference Doering and Constantin1996; Otero et al. Reference Otero, Wittenberg, Worthing and Doering2002) and for free-slip isothermal BCs (Doering & Constantin Reference Doering and Constantin1996). 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 Reference Johnston and Doering2009; Goluskin et al. Reference Goluskin, Johnston, Flierl and Spiegel2014), while replacing no-slip with free-slip leads to zonal flows with reduced vertical heat transfer (Goluskin et al. Reference Goluskin, Johnston, Flierl and Spiegel2014; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Verzicco and Lohse2014). 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 Reference Whitehead and Doering2011). 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. Reference van der Poel, Ostilla-Mónico, Verzicco and Lohse2014). 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 (Reference Whitehead and Doering2012) 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. Reference Fantuzzi, Wynn, Goulart and Papachristodoulou2017a ,Reference Fantuzzi, Wynn, Goulart and Papachristodoulou 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 (Reference Tobasco, Goluskin and Doering2017) 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 Reference Goulart and Chernyshenko2012; Chernyshenko et al. Reference Chernyshenko, Goulart, Huang and Papachristodoulou2014; Fantuzzi et al. Reference Fantuzzi, Goluskin, Huang and Chernyshenko2016; Goluskin Reference Goluskin2017). Another option is to try and lower the bound proved here through the study of optimal ‘wall-to-wall’ transport problems (Hassanzadeh, Chini & Doering Reference Hassanzadeh, Chini and Doering2014; Tobasco & Doering Reference Tobasco and Doering2017). 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.

Acknowledgements

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.

Appendix A. Boundedness of ${\mathcal{V}}$

The Cauchy–Schwarz inequality and the estimate $\langle \left|\unicode[STIX]{x1D703}\right|^{2}\rangle =\langle \left|T+z\right|^{2}\rangle \leqslant 2\langle \left|T\right|^{2}\rangle +2/3$ imply that the functional in (3.2) is bounded if $\langle \left|\boldsymbol{u}\right|^{2}\rangle ,\langle \left|T\right|^{2}\rangle <\infty$ . Following ideas by Doering & Constantin (Reference Doering and Constantin1992) and Hagstrom & Doering (Reference Hagstrom and Doering2014), this holds if velocity and temperature perturbations $\hat{\boldsymbol{u}}:=\boldsymbol{u}-\unicode[STIX]{x1D74D}$ and $\unicode[STIX]{x1D717}:=T-\unicode[STIX]{x1D70F}$ from steady background fields $\unicode[STIX]{x1D74D}$ and $\unicode[STIX]{x1D70F}$ satisfy $\langle |\hat{\boldsymbol{u}}|^{2}\rangle ,\langle |\unicode[STIX]{x1D717}|^{2}\rangle <\infty$ . Below we briefly outline how to find suitable $\unicode[STIX]{x1D74D}$ and $\unicode[STIX]{x1D70F}$ .

Let $T(\cdot ,0)$ and $\boldsymbol{u}(\cdot ,0)$ be given initial conditions. The volume-averaged temperature and horizontal velocities ( $\langle T\rangle$ , $\langle u\rangle$ and $\langle v\rangle$ ) are conserved, e.g. $\langle T(\cdot ,t)\rangle =\langle T(\cdot ,0)\rangle$ . This follows after taking the volume average of the Boussinesq equations using the divergence theorem, incompressibility, and the BCs. Then, let $\unicode[STIX]{x1D74D}:=\langle u(\cdot ,0)\rangle \boldsymbol{e}_{x}+\langle v(\cdot ,0)\rangle \boldsymbol{e}_{y}$ and set $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}(z)$ with

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

for some $\unicode[STIX]{x1D6FF}>0$ to be determined and the constant of integration chosen such that $\int _{0}^{1}\unicode[STIX]{x1D70F}(z)\,\text{d}z=\langle T(\cdot ,0)\rangle$ . It follows that $\hat{\boldsymbol{u}}$ satisfies the same BCs as the full velocity field, $\unicode[STIX]{x1D717}$ satisfies $\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D717}|_{z=0}=0=\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D717}|_{z=1}$ , and $\langle \unicode[STIX]{x1D717}\rangle =\langle \hat{u} \rangle =\langle \hat{v}\rangle =0$ at all times.

Since $\unicode[STIX]{x1D74D}$ and $\unicode[STIX]{x1D70F}$ are independent of time and $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D74D}=0$ , for any constant $C>0$ we can use incompressibility, the BCs, and the Boussinesq equations to write

(A 2) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\left\langle \frac{\left|\unicode[STIX]{x1D717}\right|^{2}}{2}+\frac{\left|\hat{\boldsymbol{u}}\right|^{2}}{2\mathit{Pr}\,\mathit{R}}\right\rangle =-\left\langle \left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\right|^{2}+\frac{\left|\unicode[STIX]{x1D735}\hat{\boldsymbol{u}}\right|^{2}}{\mathit{R}}+(\unicode[STIX]{x1D70F}^{\prime }-1){\hat{w}}\unicode[STIX]{x1D717}+(\unicode[STIX]{x1D70F}^{\prime }+1)\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D717}+C\right\rangle +C.\end{eqnarray}$$

The task is then to find $\unicode[STIX]{x1D6FF}$ in (A 1), $C>0$ , and a constant $\unicode[STIX]{x1D6FE}>0$ such that

(A 3) $$\begin{eqnarray}\left\langle \left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\right|^{2}+\frac{\left|\unicode[STIX]{x1D735}\hat{\boldsymbol{u}}\right|^{2}}{\mathit{R}}+(\unicode[STIX]{x1D70F}^{\prime }-1){\hat{w}}\unicode[STIX]{x1D717}+(\unicode[STIX]{x1D70F}^{\prime }+1)\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D717}+C\right\rangle -\unicode[STIX]{x1D6FE}\left\langle \frac{\left|\unicode[STIX]{x1D717}\right|^{2}}{2}+\frac{\left|\hat{\boldsymbol{u}}\right|^{2}}{2\mathit{Pr}\,\mathit{R}}\right\rangle \geqslant 0\end{eqnarray}$$

for all time-independent trial fields $\hat{\boldsymbol{u}}$ and $\unicode[STIX]{x1D717}$ with $\langle \unicode[STIX]{x1D717}\rangle =\langle \hat{u} \rangle =\langle \hat{v}\rangle =0$ and $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}}=0$ that satisfy the BCs. In fact, combining (A 2) and (A 3) shows that $\langle |\unicode[STIX]{x1D717}|^{2}/2+|\hat{\boldsymbol{u}}|^{2}/(2\mathit{Pr}\,\mathit{R})\rangle$ decays when it is large, remaining bounded. Hence, $\langle |\hat{\boldsymbol{u}}|^{2}\rangle$ and $\langle \left|\unicode[STIX]{x1D717}\right|^{2}\rangle$ are also bounded.

Inequality (A 3) can be proved wavenumber by wavenumber upon considering horizontal Fourier expansions for $\unicode[STIX]{x1D717}$ and $\hat{\boldsymbol{u}}$ provided that (i) $\unicode[STIX]{x1D6FE}<\min \{4,\,4\mathit{Pr},\,2k_{m}^{2},\,2\mathit{Pr}k_{m}^{2}\}$ with $k_{m}^{2}:=\min _{\boldsymbol{k}\neq (0,0)}k^{2}$ (here $k^{2}$ is the magnitude of the horizontal wavevector, cf. § 3; the minimum is strictly positive because we work in a finite periodic domain), (ii)  $C$ is sufficiently large, and (iii) $\unicode[STIX]{x1D6FF}$ is sufficiently small. Non-zero wavevectors can be analysed using estimates similar to those of § 4, while the case $\boldsymbol{k}=(0,0)$ is handled using Poincaré-type inequalities deduced using the zero-average conditions $\langle \unicode[STIX]{x1D717}\rangle =\langle \hat{u} \rangle =\langle \hat{v}\rangle =0$ .

References

Ahlers, G., Grossmann, S. & Lohse, D. 2009 Heat transfer and large scale dynamics in turbulent Rayleigh–Bénard convection. Rev. Mod. Phys. 81 (2), 503537.Google Scholar
Chapman, C. J. & Proctor, M. R. E. 1980 Nonlinear Rayleigh–Bénard convection between poorly conducting boundaries. J. Fluid Mech. 101 (4), 759782.Google Scholar
Chernyshenko, S. I.2017 Relationship between the methods of bounding time averages. arXiv:1704.02475 [physics.flu-dyn].Google Scholar
Chernyshenko, S. I., Goulart, P. J., Huang, D. & Papachristodoulou, A. 2014 Polynomial sum of squares in fluid dynamics: a review with a look ahead. Phil. Trans. R. Soc. Lond. A 372 (2020), 20130350.Google Scholar
Choffrut, A., Nobili, C. & Otto, F. 2016 Upper bounds on Nusselt number at finite Prandtl number. J. Differ. Equ. 260 (4), 38603880.Google Scholar
Constantin, P. & Doering, C. R. 1995 Variational bounds on energy dissipation in incompressible flows. II. Channel flow. Phys. Rev. E 51 (4), 31923198.Google Scholar
Doering, C. R. & Constantin, P. 1992 Energy dissipation in shear driven turbulence. Phys. Rev. Lett. 69 (11), 16481651.Google Scholar
Doering, C. R. & Constantin, P. 1994 Variational bounds on energy dissipation in incompressible flows: shear flow. Phys. Rev. E 49 (5), 40874099.Google Scholar
Doering, C. R. & Constantin, P. 1996 Variational bounds on energy dissipation in incompressible flows. III. Convection. Phys. Rev. E 53 (6), 59575981.Google Scholar
Fantuzzi, G., Goluskin, D., Huang, D. & Chernyshenko, S. I. 2016 Bounds for deterministic and stochastic dynamical systems using sum-of-squares optimization. SIAM J. Appl. Dyn. Syst. 15 (4), 19621988.Google Scholar
Fantuzzi, G., Wynn, A., Goulart, P. J. & Papachristodoulou, A. 2017a Optimization with affine homogeneous quadratic integral inequality constraints. IEEE Trans. Autom. Control 62 (12), 62216236.Google Scholar
Fantuzzi, G., Wynn, A., Goulart, P. J. & Papachristodoulou, A.2017b quinopt, v2.1. Available from: https://github.com/aeroimperial-optimization/QUINOPT.Google Scholar
Goluskin, D. 2016 Internally Heated Convection and Rayleigh–Bénard Convection. Springer.Google Scholar
Goluskin, D. 2017 Bounding averages rigorously using semidefinite programming: mean moments of the Lorenz system. J. Nonlinear Sci. (in press); doi:10.1007/s00332-017-9421-2.Google Scholar
Goluskin, D., Johnston, H., Flierl, G. R. & Spiegel, E. A. 2014 Convectively driven shear and decreased heat flux. J. Fluid Mech. 759, 360385.Google Scholar
Goulart, P. J. & Chernyshenko, S. I. 2012 Global stability analysis of fluid flows using sum-of-squares. Phys. D 241 (6), 692704.Google Scholar
Hagstrom, G. I. & Doering, C. R. 2014 Bounds on surface stress-driven shear flow. J. Nonlinear Sci. 24 (1), 185199.Google Scholar
Hassanzadeh, P., Chini, G. P. & Doering, C. R. 2014 Wall to wall optimal transport. J. Fluid Mech. 751, 627662.Google Scholar
Hurle, D. T. J., Jakeman, E. & Pike, E. R. 1967 On the solution of the Bénard problem with boundaries of finite conductivity. Proc. R. Soc. Lond. A 296 (1447), 469475.Google Scholar
Johnston, H. & Doering, C. R. 2009 Comparison of turbulent thermal convection between conditions of constant temperature and constant flux. Phys. Rev. Lett. 102 (6), 064501.Google Scholar
Lappa, M. 2010 Thermal Convection: Patterns, Evolution and Stability. Wiley.Google Scholar
Otero, J., Wittenberg, R. W., Worthing, R. A. & Doering, C. R. 2002 Bounds on Rayleigh–Bénard convection with an imposed heat flux. J. Fluid Mech. 473, 191199.Google Scholar
Otto, F. & Seis, C. 2011 Rayleigh–Bénard convection: Improved bounds on the Nusselt number. J. Math. Phys. 52 (8), 083702.Google Scholar
Tobasco, I. & Doering, C. R. 2017 Optimal wall-to-wall transport by incompressible flows. Phys. Rev. Lett. 118 (26), 264502.Google Scholar
Tobasco, I., Goluskin, D. & Doering, C. R.2017 Optimal bounds and extremal trajectories for time averages in nonlinear dynamical systems. Phys. Lett. A (in press); doi:10.1016/j.physleta.2017.12.023.Google Scholar
van der Poel, E. P., Ostilla-Mónico, R., Verzicco, R. & Lohse, D. 2014 Effect of velocity boundary conditions on the heat transfer and flow topology in two-dimensional Rayleigh–Bénard convection. Phys. Rev. E 90 (1), 013017.Google Scholar
Whitehead, J. P. & Doering, C. R. 2011 Ultimate state of two-dimensional Rayleigh–Bénard convection between free-slip fixed-temperature boundaries. Phys. Rev. Lett. 106 (24), 244501.Google Scholar
Whitehead, J. P. & Doering, C. R. 2012 Rigid bounds on heat transport by a fluid between slippery boundaries. J. Fluid Mech. 707, 241259.Google Scholar
Whitehead, J. P. & Wittenberg, R. W. 2014 A rigorous bound on the vertical transport of heat in Rayleigh–Bénard convection at infinite Prandtl number with mixed thermal boundary conditions. J. Math. Phys. 55 (9), 093104.Google Scholar
Wittenberg, R. W. 2010 Bounds on Rayleigh–Bénard convection with imperfectly conducting plates. J. Fluid Mech. 665, 158198.Google Scholar
Wittenberg, R. W. & Gao, J. 2010 Conservative bounds on Rayleigh–Bénard convection with mixed thermal boundary conditions. Eur. Phys. J. B 76 (4), 565580.Google Scholar
Figure 0

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

Figure 1

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.