Exhausting the background approach for bounding the heat transport in Rayleigh-B\'enard convection

TL;DR

This paper explores the limits of bounding heat transport in 2D Rayleigh-Bénard convection using a background field approach, showing that the bounds cannot be improved beyond a certain scaling despite complexifying the models.

Contribution

It exhaustively analyzes the background approach for bounding heat transport, demonstrating its limitations and the inability to lower the established bounds with more complex models.

Findings

Optimal bounds approach $Ra^{1/2}$ scaling.

Imposing full heat and velocity constraints does not lower the bound.

Adding a velocity background field does not improve the bound.

Abstract

We revisit the optimal heat transport problem for Rayleigh-B\'enard convection in which a rigorous upper bound on the Nusselt number, $Nu$, is sought as a function of the Rayleigh number $Ra$. Concentrating on the 2-dimensional problem with stress-free boundary conditions, we impose the full heat equation as a constraint for the bound using a novel 2-dimensional background approach thereby complementing the `wall-to-wall' approach of Hassanzadeh \etal \,(\emph{J. Fluid Mech.} \textbf{751}, 627-662, 2014). Imposing the same symmetry on the problem, we find correspondence with their result for $Ra \leq Ra_c:=4468.8$ but, beyond that, the optimal fields complexify to produce a higher bound. This bound approaches that by a 1-dimensional background field as the length of computational domain $L\rightarrow\infty$. On lifting the imposed symmetry, the optimal 2-dimensional temperature background field reverts back to being 1-dimensional giving the best bound $Nu\le 0.055Ra^{1/2}$ compared to $Nu \le 0.026Ra^{1/2}$ in the non-slip case. % We then show via an inductive bifurcation analysis that imposing the full time-averaged Boussinesq equations as constraints (by introducing 2-dimensional temperature {\em and} velocity background fields) is also unable to lower this bound. This then exhausts the background approach for the 2-dimensional (and by extension 3-dimensional) Rayleigh-Benard problem with the bound remaining stubbornly $Ra^{1/2}$ while data seems more to scale like $Ra^{1/3}$ for large $Ra$. % Finally, we show that adding a velocity background field to the formulation of Wen \etal\, (\emph{Phys. Rev. E.} \textbf{92}, 043012, 2015), which is able to use an extra vorticity constraint due to the stress-free condition to lower the bound to $ Nu \le O(Ra^{5/12})$, also fails to improve the bound.

Equations

$$\displaystyle(\mbox{\boldmath${\cal N}$})$$

$$(\mathcal{H})$$

$$\sigma :=\nu /\kappa \&Ra:=g\beta \delta T{d}^3/\nu \kappa$$

$${\mathscr{L}}:=\overline{\langle|\mbox{\boldmath$\nabla$}T|^{2}\rangle}^{t}-\overline{\langle a\mbox{\boldmath$\nu$}\cdot(\mbox{\boldmath${\cal N}$})\rangle}^{t}-\overline{\langle b\theta({\cal H})\rangle}^{t}$$

$${\overline{()}}^t:=\underset{\mathcal{T}\to \infty }{\lim }\sup \frac{1}{\mathcal{T}}{\int }_0^{\mathcal{T}}()dt,⟨⟩:=\frac{1}{V}\int dV.$$

$$\mbox{\boldmath$\phi$}(\mathbf{x}):=\mathbf{u}(\mathbf{x},t)-\mbox{\boldmath$\nu$}(\mathbf{x},t),\quad\tau(\mathbf{x}):=T(\mathbf{x},t)-\theta(\mathbf{x},t)$$

$$\mathcal{L}$$

$${\mathscr{L}}_{s}:=\langle|\mbox{\boldmath$\nabla$}T|^{2}\rangle-a\langle\mathbf{u}\cdot(\mbox{\boldmath${\cal N}$})_{s}\rangle+a\langle\mbox{\boldmath$\phi$}\cdot(\mbox{\boldmath${\cal N}$})_{s}\rangle-b\langle T({\cal H})_{s}\rangle+b\langle\tau({\cal H})_{s}\rangle$$

$${\mathscr{L}}_{s}:=\langle|\mbox{\boldmath$\nabla$}T|^{2}\rangle-a\langle\mathbf{u}\cdot(\mbox{\boldmath${\cal N}$})_{s}\rangle-b\langle T({\cal H})_{s}\rangle+b\int^{d}_{0}\tau(z)\overline{({\cal H})_{s}}^{x,y}\,dz$$

$${\overline{()}}^{x,y}:=\underset{L\to \infty }{\lim }\frac{1}{L^2}{\int }_{-L/2}^{L/2}{\int }_{-L/2}^{L/2}()dxdy$$

$$\frac{\partial u}{\partial z}=w=0,$$

$$T=\tau (x,z)+\theta (x,z,t).$$

$$Nu:=\lim_{T\rightarrow\infty}\frac{1}{T}\int^{T}_{0}\frac{1}{L}\int^{L/2}_{-L/2}\,\frac{\partial T}{\partial z}\biggl{|}_{z=1}\biggr{.}\,dx\,dt=\langle|\boldsymbol{\nabla}T|^{2}\rangle=1+\langle wT\rangle$$

$$⟨(\dots )⟩:=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T{\int }_0^1\frac{1}{L}{\int }_{-L/2}^{L/2}(\dots )dxdzdt$$

$$\mathcal{L}$$

$$\mathcal{L}=\frac{1}{1-a}\left[⟨{|\nabla \tau |}^2⟩-a\right]-\frac{1}{1-a}\mathcal{G}$$

$$\mathcal{G}:=⟨\frac{a}{Ra}|\nabla \mathbf{u}{|}^2+|\nabla \theta {|}^2+2\theta \mathbf{u}\cdot \nabla \tau ⟩$$

$$Nu\le \frac{1}{1-a}\left[⟨{|\nabla \tau |}^2⟩-a\right].$$

$$(\tau ,a)\in \Omega :=\{(\tau ,a)|\mathcal{G}(\mathbf{u},\theta ;\tau ,a)\ge 0\forall (\mathbf{u},\theta )\in \Pi \}.$$

$$\mathcal{G}=⟨\frac{a}{Ra}|{\nabla }^2\psi {|}^2+|\nabla \theta {|}^2+2\theta J(\tau ,\psi )⟩\ge 0$$

$$J(A,B):=\frac{\partial A}{\partial x}\frac{\partial B}{\partial z}-\frac{\partial A}{\partial z}\frac{\partial B}{\partial x}$$

$$\lambda \theta$$

$$\lambda {\nabla }^2\psi$$

$$0$$

$$⟨{|\nabla \tau |}^2⟩-1$$

$$f(\tau ,a):=\frac{1}{1-a}\left[⟨{|\nabla \tau |}^2⟩-a\right]$$

$$\frac{1}{(1-a{)}^2}⟨\nabla |(1-a)\delta \tau +(\tau +z-1)\delta a{|}^2⟩,$$

$$\psi =\sum _{m=0}^M{\psi }_m(z)\sin (m\alpha x),\theta =\sum _{m=0}^M{\theta }_m(z)\cos (m\alpha x),\tau =\sum _{m=0}^M{\tau }_m(z)\cos (m\alpha x).$$

$$[{\psi }_m(z),{\theta }_m(z),{\tau }_m(z)]=\sum _{n=0}^N[{\psi }_{mn},{\theta }_{mn},{\tau }_{mn}]T_n(z)$$

$$0$$

Full text

1 Introduction

In this paper we consider the fundamental problem of assessing how the heat flux per unit area behaves as a function of the Rayleigh number, $Ra$, in Rayleigh-Benard convection where a layer of fluid is heated from below and cooled from above. This situation is ubiquitous in Nature and consequently the focus of a huge body of ongoing research work (e.g. Ahlers et al. 2009). The particular focus here is the use of variational methods which seek an upper bound on the heat flux in the hope that this bound will capture the correct high-$Ra$ scaling for turbulent convection. This approach involves constructing an optimisation problem constrained by information gleaned from the governing equations. Inevitably, the constraints actually imposed form a strict subset of those implied by the governing equations so that any maximum which emerges is an upper bound on what can actually be realised. This approach has its roots in the work of Malkus (1954) who hypothesized that the fluid selects the flow state from all those possible states which maximises the heat transport. The subsequent mathematical formulation by Howard (1963) and Busse (1969) was as a maximization problem (see the early reviews by Howard (1972) and Busse (1978) ). In the 1990s, an alternative complementary approach - the background method - was introduced by Doering & Constantin (1992,1994,1995,1996) which takes the form of a minimization problem. This has the considerable advantage that even a trial solution can yield an upper bound which, experience seems to indicate, yields the same scaling as the proper optimal (e.g. in shear flow and convection see Doering & Constantin 1992,1996 respectively compared to Plasting & Kerswell 2003, hereafter PK03).

In both approaches, however, the outstanding challenge has been to add further dynamical information to improve (lower) the scaling law (e.g. see Ierley & Worthing (2001) for efforts in the Howard-Busse maximization problem). The best current bound on the Nusselt number - the ratio of actual heat flux to the conductive value - is $Nu\leq 0.02634Ra^{1/2}$ as $Ra\rightarrow\infty$ (PK03) whereas most of the current experimental data suggests $Nu\sim Ra^{0.31}$ (see the discussion in Waleffe et al. 2015) and so is more consistent with the simple theoretical prediction of $Nu\sim Ra^{1/3}$ (Malkus 1954, Priestley 1954) with some dependence on the Prandtl number also possible (Grossmann & Lohse 2000). A natural way of incorporating further information exists in the background method through simply extending the definitions of the background fields. To see this, recall that the Malkus-Howard-Busse (maximization) approach and the Doering-Constantin (minimization) approach are dual problems seeking to find an appropriate saddle point of a functional of the velocity and temperature fields (Kerswell 1998, 2001). To explain further we introduce the problem to be considered.

Let a Newtonian fluid be confined between two infinite isothermal plates at $z=0$ and $z=d$ with the lower plate maintained at a constant temperature $\delta T$ hotter than that of the upper plate (gravity is $-g\mathbf{\hat{z}}$ where $g\approx 9.8ms^{-2}$). Using the gap width $d$, $d^{2}/\kappa$ ($\kappa$ is the thermal diffusivity) and $\delta T$ as units of length, time and temperature together with adopting the Boussinesq approximation, the governing equations are

[TABLE]

with \mbox{\boldmath\nabla}\cdot\mathbf{u}=0 where

[TABLE]

are the Prandtl and Rayleigh numbers respectively ($\nu$ is the kinematic viscosity and $\beta$ is the thermal expansion coefficient). The background method starts by writing down the functional

[TABLE]

where the first term on the right is the long-time-averaged Nusselt number $Nu$, \mbox{\boldmath\nu}(\mathbf{x},t) and $\theta(\mathbf{x},t)$ are Lagrange multipliers imposing the momentum and heat equations as constraints respectively (the seemingly redundant extra scalars $a$ and $b$ play a key role later) and the time and volume averages are defined as follows

[TABLE]

The crucial next step is to define steady ‘background’ fields

[TABLE]

which connect the Lagrange multipliers with the physical fields and such that they carry any inhomogeneous boundary conditions (so here just those on the temperature field). Changing variables from (\mathbf{u},T,\mbox{\boldmath\nu},\theta) to (\mathbf{u},T,\mbox{\boldmath\phi},\tau),

[TABLE]

makes it clear that choosing the largest stationary value of ${\mathscr{L}}$ finds the largest long-time-averaged Nusselt number subject to the long-time-averaged power and entropy balances (Lagrange multipliers $a$ and $b$ respectively) and projected information from momentum and heat flux balances (Lagrange multipliers $\phi$ and $\tau$ respectively). Since it can be shown that all the time derivative terms in these constraints vanish under long time averaging, the variational problem can be couched in terms of steady fields only. In particular, the goal is to evaluate the largest stationary value of the functional

[TABLE]

where the subscript $s$ indicates the steady version of the unsubscripted quantity. So far only the minimal choice (\tau,\mbox{\boldmath\phi})=(\tau(z),{\bf 0}) has been explored (Doering & Constantin 1996) which leads to the simplified expression

[TABLE]

(where

[TABLE]

is a horizontal average). This choice turns out to give the dual problem to the Howard-Busse approach (Howard 1963, Busse 1969) and so produces the same Nusselt number bound (Kerswell 2001, PK03). However, here, beyond the total power and entropy balances and insisting that the fluid is incompressible and the boundary conditions are satisfied, only the horizontally-averaged heat equation is imposed as a constraint. It seems reasonable to suppose that imposing further constraints from the governing equations by extending the definitions of the background fields should lower this current best bound for Boussinesq convection. Probing this hypothesis is the motivation for this paper.

This issue is quite general applying to bounds developed in other canonical flows such as plane Couette flow (Doering & Constantin 1992), channel flow (Doering & Constantin 1994) and pipe flow (Plasting & Kerswell 2005). However, the most bounding work has been performed in convection partly because of its wide application and partly because the current best bound appears to have the wrong exponent and therefore calls for the most improvement. Of particular interest in recent efforts to lower the bound has been the introduction of the ‘wall-to-wall’ approach by Hassanzadeh et al. (2014) (see also Souza 2016 and Souza et al. 2019). Here the full heat equation has been imposed as a constraint with some incompressible boundary-compliant flow field which, apart from an overall amplitude, is otherwise unconstrained and a maximization problem is solved. This appears to give a much improved upper bound of $Nu\sim Ra^{5/12}$ for stress-free boundary conditions in 2D convection with Souza (2016) finding a yet stronger (lower) bound of $Nu\sim Ra^{0.371}$ for non slip boundary conditions. Later work by Tobasco & Doering (2017), however, has demonstrated through designing a sophisticated trial function that the upper bound must be at least $Nu\sim Ra^{1/2}$ up to logarithms for non-slip boundary conditions. This directly contradicts the conclusions of Souza (2016) and indirectly those of Hassanzadeh et al. (2014) (the heat flux for stress-free walls should be higher than for non-slip walls). Resolving this paradox by tackling the complementary background formulation - (\tau,\mbox{\boldmath\phi})=(\tau(x,z),{\bf 0}) which imposes the full heat equation in 2D convection and builds a minimization problem - is a good starting point for our study.

Concurrent work by Souza et al. (2019) has considered how the background method is connected to the wall-to-wall approach and speculated that there could be a ‘duality gap’ between them. Coming from a different perspective (the specific details of solving the variational equations), we share this speculation and confirm it here beyond a certain Rayleigh number. Motoki et al. (2018) have also built upon Hassanzadeh et al’s work by extending the maximization search to 3 dimensions. Interestingly they find a 3D optimal solution which scales like $Ra^{1/2}$ with a numerical coefficient just 7.2% below the bound of PK03 (see their figure 2). This 3D result (using stress-free boundary conditions) and the 2D work of Tobasco & Doering (2017) (using non-slip boundary conditions) clearly beg the question whether further information from the momentum equation can be used to rule out the $Ra^{1/2}$ scaling which clearly persists despite imposing the full heat equation. This also will be addressed here.

A further motivation for exploring the addition of further dynamical constraints is the hope that ultimately, the full steady governing equations can be imposed and then finally a direct connection forged between the optimal solution of an upper bounding variational problem and an actual solution of the (steady) governing equations. On the one hand, it seems clear that the best bound possible using the background approach can’t be lower than the highest heat flux attained by any of the many steady solutions of the governing equations whereas on the other, as argued above, the optimal solution should ultimately be a solution of the steady equations. This possibility has been brought into focus by recent computations tracking the (simple) 2D convection roll solution which initially bifurcates from the conductive state up to very high Rayleigh numbers of $O(10^{9})$ (Waleffe et al. 2015, Sondak et al. 2015). Provided the aspect ratio of the rolls is optimized over, a heat flux relationship of $Nu\sim Ra^{0.31}$ is found which is intriguingly close to 3D turbulent convection measurements and $Nu\sim Ra^{1/3}$ which is what many believe might be the ultimate scaling law although not all (e.g. Zhu et al. 2018).

A synopsis of the paper is as follows. Section 2 describes the set-up of 2D Boussinesq convection (§2.1), explains how a bound can be found using the background approach (§2.2) and then discusses the convexity of the optimization problem for a general temperature background field which ensures a unique optimal (§2.3). Section 2.4 explains how the numerical computations are performed with a choice having to be made between a branch continuation approach (PK03) and a time stepping method (Wen et al. 2013, 2015). Section 3 describes the results of tackling the upper bounding problem with the full heat equation imposed in the presence of the same symmetry as used in Hassanzadeh et al. (2014). The appearance of a second fluctuation mode becoming ‘spectrally unstable’ at $Ra=Ra_{c}:=4468.8$ means: a) that Hassanzadeh et al’s result is no longer a bound for $Ra>Ra_{c}$ (i.e. their result becomes only a local maximum and the duality gap suggested by Souza et al. 2019 is realised); and b) a new formulation for how the optimal is tracked needs to be introduced compared to previous work (e.g. PK03). Section 4 discusses this new formulation which is significant because the various background and fluctuation optimal fields can no longer be used to define a set of physical temperature and velocity fields. In particular, the optimal fields do not satisfy the steady heat equation even though this is explicitly imposed as a constraint. Using this reformulation, section 5 shows how the optimal bound behaves for $Ra>Ra_{c}$. The size of the computational domain becomes important in the 2D background problem and it is found that the highest bound is only achieved in the infinite domain limit when the background field becomes increasingly 1D. Removing the symmetry used by Hassanzadeh et al. restores the translational invariance of the problem in which case the optimal has to be 1D and a bound of $Ra\leq 0.055Ra^{1/2}$ is found compared to the well-known result of $0.026Ra^{1/2}$ for non-slip walls (PK03). Having found that imposing the full heat equation does not improve the bound, we then consider adding extra information from the momentum equation by introducing a background velocity field \mbox{\boldmath\phi}(x,z). Now the optimization problem is no longer convex and so we are unable to invoke uniqueness to dismiss non-vanishing $\phi$. Instead we use an inductive bifurcation analysis to show that if \mbox{\boldmath\phi}={\bf 0} before a bifurcation then it remains ${\bf 0}$ after it too meaning that the continuous branch of optimals found by branch tracking out of the energy stability point always has \mbox{\boldmath\phi}={\bf 0}. Noting the one caveat that it’s not impossible that there is an unconnected branch of optimals with \mbox{\boldmath\phi}\neq{\bf 0}, this strongly suggests the surprising result that imposing the full Boussinesq equations does not improve the bound over that obtained using the horizontally-averaged Boussinesq equations. Finally, section 7 observes that adding a velocity background temperature field to the formulation of Wen et al. (2015), which has an additional vorticity constraint, also fails to improve matters. A discussion follows in section 8.

2 Mathematical formulation

2.1 Set-up

We consider the 2 dimensional version of the Boussinesq equations (1) & (2) where $\boldsymbol{u}=u\boldsymbol{\hat{x}}+w\boldsymbol{\hat{z}}$ over a box $(x,z)\in[-{\textstyle{1\over 2}}L,{\textstyle{1\over 2}}L]\times[0,1]$ together with the following stress-free boundary conditions

[TABLE]

following Hassanzadeh et al. (2014). Applying the background method, we decompose the temperature field as

[TABLE]

where the (steady) background temperature $\tau$ carries the boundary conditions of $T$ (i.e. $\tau|_{z=0}=1$ and $\tau|_{z=1}=0$) so that the perturbation field $\theta$ vanishes at $z=0,1$. The time-averaged heat transport is characterized by the time-averaged Nusselt number $Nu$

[TABLE]

in which henceforth

[TABLE]

is the spatial-temporal average. To find the maximum heat transport possible over all solutions to the Boussinesq equations, we construct the Lagrangian

[TABLE]

where $a/\sigma Ra$ is the Langrange multiplier imposing the global constraint \langle\boldsymbol{u}\cdot\mbox{\boldmath{\cal N}}\rangle=0, $b$ is a Lagrange multiplier imposing the global constraint $\langle T\,{\cal H}\rangle=0$ and $b\tau(x,z)$ is the Lagrange multiplier field imposing the time-averaged heat equation pointwise in the domain. The inclusion of $b$ is actually redundant given the constraint imposed by $\tau$ implies $\langle T\,{\cal H}\rangle=0$ so the value of $b$ is chosen for convenience. Expression (16) can be rewritten using integration by parts and the fact that $\langle wT\rangle={\cal L}-1$ (see (14)) for solutions of the Boussinesq equations as

[TABLE]

where setting $b=2$ makes

[TABLE]

a purely quadratic form in $\theta$ and $\boldsymbol{u}$.

2.2 Bounds on $Nu$

The key realisation is that if $\mathscr{G}\geq 0$ for all $(\boldsymbol{u},\theta)\in\Pi$ (the set of incompressible velocity and temperature fields which satisfy homogeneous versions of the boundary conditions (11) and (12) ), which is a spectral constraint on $\tau$ and $a\in(0,1)$, we then have the bound

[TABLE]

The challenge is then to find the lowest such bound by minimizing over all $(\tau,a)$ which satisfy this spectral constraint, i.e.

[TABLE]

After introducing a streamfunction, $(u,w)=(\partial\psi/\partial z,-\partial\psi/\partial x$), the constraint that

[TABLE]

where

[TABLE]

is equivalent to requiring that all of the eigenvalues $\lambda$ of the linear problem

[TABLE]

(with boundary conditions $\psi=d^{2}\psi/dz^{2}=\theta=0$ for $z=\{0,1\}$) are negative semi-definite.

2.3 Convexity & Uniqueness

The Euler-Lagrange equations for stationarizing the Lagrangian ${\mathscr{L}}$ in (17) are

[TABLE]

and, as a nonlinear set of equations, can have many solutions. However, only solutions with $(\tau,a)\in\Omega$ yield a bound through the value of ${\mathscr{L}}$ generated. Due to the convexity of $\Omega$ (i.e. if $(\tau_{1},a_{1})$ and $(\tau_{2},a_{2})$ are in $\Omega$ then so is $\lambda(\tau_{1},a_{1})+(1-\lambda)(\tau_{2},a_{2})$ for $\lambda\in(0,1)$), and the fact that the objective functional

[TABLE]

to be minimized is a strictly convex functional (the terms second order in $\delta\tau$ and $\delta a$ in the difference $f(\tau+\delta\tau,a+\delta a)-f(\tau,a)$, specifically

[TABLE]

are positive definite), there is in fact at most one solution which satisfies the spectral constraint. This solution, hereafter referred to as the optimal solution, is what is sought.

2.4 Numerical approach

Recently, [Wen et al. (2015)] have proved that when $\tau$ is 1-dimensional, i.e. $\tau=\tau(z)$, appropriately augmenting the (steady) Euler-Lagrange equations with time derivatives leads to a system where the optimal solution is a unique attracting steady state. This proof carries over to 2-dimensional background fields $\tau=\tau(x,z)$ in the 3-dimensional Rayleigh-Benard problem but not in the 2-dimensional problem (see appendix A for details) where the dimensionality of the background field then matches that of the physical fields. This means that any steady attractor which emerges from time-stepping using $\tau(x,z)$ in the 2D problem is not guaranteed to be the required optimal solution. The time stepping approach can still be used if it is married with a spectral constraint check but then there is always the prospect of rerunning with different initial conditions until the optimal solution is found. Given this, we chose instead to use the branch continuation approach - Newton’s method with parametric continuation - starting from the energy stability bifurcation point as performed in PK03. While very robust, this has the general disadvantage of only being able to continuously trace optimal solutions from the energy stability bifurcation as $Ra$ varies meaning that any new unconnected optimals cannot be found at a given $Ra$. This is not a problem here as the aforementioned uniqueness of the optimal solution means that no other optimal solution branches exist.

We consider periodic boundary conditions in $x$ and, exactly as in Hassanzadeh et al. (2014), assume that the streamfunction $\psi$ is odd (or antisymmetric), while $\theta$ and $\tau$ are even (or symmetric) about $x=0$ by seeking the solution of (25)-(28) in the following form:

[TABLE]

We will find that this choice prevents a 1D background optimal even though this is allowed by the boundary conditions and imposed symmetry. Here $\alpha:=2\pi/L$ and $\psi_{m}$, $\theta_{m}$, $\tau_{m}$ are expanded in Chebyshev polynomials, $T_{n}$,

[TABLE]

where $T_{n}(z):=\cos(n\cos^{-1}(2z-1))$. Resolution varies from $(N,M)=(30,30)$ to $(80,80)$ to ensure numerical accuracy as $Ra$ increases and $L$ changes .

3 Connecting to Hassanzadeh et al.(2014)

The conductive temperature profile $\tau=1-z$ is a spectrally-stable background field until $Ra=27\pi^{4}/4$ in a domain of size $L=2\sqrt{2}$ when energy instability first starts. Pinning the marginal fluctuation fields $(\theta,\psi)$ - hereafter a mode - as was done in PK03, the optimal solution was then tracked up to $Ra=10^{7}$ with the domain size $L=L^{*}(Ra)$ simultaneously optimized to yield the highest heat flux at a given $Ra$: see figure 1. The calculated $Nu$ values correspond exactly with those found by Hassanzadeh et al. (2014) (as do flow fields computed at $Ra=10^{5}$ and $10^{6}$; see the inset of 1). This indicates that Hassanzadeh et al.’s (2014) wall-to-wall transport approach is equivalent to the background method when a single mode is considered.

In their wall-to-wall optimal control approach, however, Hassanzadeh et al. (2014) had no way of identifying whether their local optimal was in fact the global optimal. It should be sufficiently close to the energy stability point but experience in other related problems (e.g. PK03) suggests that further modes in the spectral constraint eventually become marginal as $Ra$ increases. The optimal solution should subsequently modify itself to keep these new modes marginal with concomitant adjustments in the $Nu$-scaling. Fortunately, in the background approach, the spectral constraint provides a check on whether a given Euler-Lagrange solution is the optimal solution. Solving the eigenvalue problem (23)-(24) for disturbances which are also periodic over $[0,L^{*}(Ra)]$ demonstrates that the eigenvalue ($\lambda_{1}$ in figure 2) of the first mode is pinned at 0, while a second mode becomes marginal at $Ra=4468.8$ for an aspect ratio $L^{*}=2.234$. This means that [Hassanzadeh et al.(2014)]’s result is not a bound for $Ra>4468.8$. This apparently is not because any further 2D bifurcation has been missed (Chini, private communication) but more a reflection of the ‘duality gap’ suggested by Souza et al. (2019) being realised. Figure 3(a,b) shows the first mode $(\psi_{1},\theta_{1})$ with wavenumber $\alpha_{1}:=2\pi/L^{*}$ so the flow field contains one pair of convection cells. The second mode ($\psi_{2},\theta_{2}$) with $\alpha_{2}=2\alpha_{1}$ illustrated in figure 3(c,d) has two pairs of convection cells. The optimal background field at $Ra=2000$ is shown in figure 4 and for the now non-optimal 1-mode solution at $Ra=20,000$. In both cases the field is weakly 2-dimensional indicating that the first mode consists of non-monochromatic (i.e. non-single $\alpha$) velocity and temperature fields. The emergence of the second mode at $Ra=4468.8$ indicates that the background profile is now degenerate in a way which has important implications for solving the Euler-Lagrange equations for higher $Ra$ while respecting the spectral constraint. We discuss this in the following section.

4 Multi-modal optimals

When a new mode becomes marginal in the spectral constraint as the background field $\tau$ evolves with $Ra$, a further ‘pinning’ constraint needs to be added to keep the new mode marginal in the spectral constraint as $Ra$ increases further. This procedure is thoroughly discussed in [Doering & Constantin (1996)] and implemented in PK03 for a background field of lower dimensionality than the fluctuation field. In this situation, an example of which is using a 1D case $\tau=\tau(z)$ in the 2D Rayleigh-Benard problem, the fluctuation field can be Fourier-transformed over the spatial dimension(s) across which $\tau$ is invariant and then considered parameterized by the Fourier wavenumber $k$. Different spectrally-marginal fluctuation fields have different $k$ and are then naturally orthogonal under averaging over this spatial dimension. This means that the Euler Lagrange equations (25)-(28),

[TABLE]

(the overbar represents averaging over $x$) can simply be extended to include the new marginal mode

[TABLE]

when it appears. Equivalently, the Lagrangian is just

[TABLE]

where $\mathscr{G}$ naturally partitions into the contributions from the various marginal modes as follows

[TABLE]

where

[TABLE]

The appearance of a new spectrally-marginal mode merely extends the set of wavenumbers contributing to the definition of the fluctuation field by one,

[TABLE]

Importantly, this means it is possible to talk about the unique optimal solution of the variational problem which satisfies the imposed physical constraints as being

[TABLE]

i.e. the spectral constraint is satisfied at a saddle point of ${\mathscr{L}}$.

This pleasing situation in which the marginal fluctuation fields have a physical interpretation changes, however, when the dimensionality of the background field equals the dimensionality of the problem (the case here), or, pathologically, there is more than one marginal mode for a given wavenumber (see chapter 3 of Fantuzzi 2018). In these scenarios, the natural orthogonality property of different marginal fluctuation fields disappears with the result that the physical meaning of the fluctuation fields is lost. To see this, the key is to realise that pinning the marginal fluctuation fields is done (Doering & Constantin 1996) as before by writing the Lagrangian as

[TABLE]

The constraint that each $\mathscr{G}_{j}$ vanishes pins the $j^{th}$ mode to be marginal (the Lagrange multiplier imposing this is absorbed into the amplitude of the $j^{th}$ marginal fluctuation field) while $\mathscr{G}>0$ for all other fluctuation fields. However, since the modes $(\psi_{j},\theta_{j})$ are not now orthogonal,

[TABLE]

for $N\geqslant 2$ where

[TABLE]

is taken as the total optimal fluctuation field. In fact, $\mathscr{G}>0$ and so this total optimal field is not a solution of the Heat equation. The clear implication is that the spectral constraint is not satisfied for $N\geqslant 2$ at any saddle point of the Lagrangian (17) where the steady heat equation is imposed. Consequently, the optimization procedure is forced to find an optimal away from the saddle points of the Lagrangian (17) where the spectral constraint is satisfied to deliver a bound.

From a different perspective, Souza et al. (2019) have also recently argued that this should happen when exploring the connection between the wall-to-wall approach (a max-min problem) with the associated background method (a min-max problem). A duality gap means that

[TABLE]

(making the connection $\eta=\tau-(1-z)$ and $\zeta=\theta$ with the variables used by Souza et al 2019) where the optimal solution to the wall-to-wall problem is achieved at a stationary point of $\mathscr{L}$ thereby implying that to the background method is not. They also supply a nice simple quadratic polynomial in 5 variables to illustrate the phenomenon. The calculations described in the next section confirm that this gap starts to exist as soon as $N=2$.

5 Extending Hassanzedah et al. with a symmetric 2D background field $\tau(x,z)$

To explore multi-modal bounding solutions, a first series of computations were done in the fixed domain $L=2\sqrt{2}$. In this geometry, the first mode appears at $Ra=27\pi^{4}/4$ (the energy stability threshold), the second mode at $Ra=3,075$ and the third mode at $Ra=24,650$. The 1-mode and 2-mode optimal solution branches could be easily continued up to $Ra=10^{5}$ whereas the 3-mode solution branch proved difficult to continue much beyond $Ra>40,000$ due to numerical issues: see figure 5(a). The 3-mode solution, which provides an upper bound in this geometry over at least the range $24,650\leq Ra\leq 40,000$, presents only a modest correction to the 2-mode optimal solution which is no longer a bound for these $Ra$.

A second series of computations were then carried out to investigate the dependence of the $Nu$-bound on the aspect ratio $L$. Three different $Ra$ values were chosen to explore the dependence of the bound on $L$: $Ra=5000$ and $10,000$ where the bound is given by a 2-mode solution, and $Ra=25,000$ where the bound is given by a 3-mode solution. In all three cases, the largest bound is achieved as the aspect ratio $L\rightarrow\infty$: see figure 5(b). This is very different from the optimal control results of [Hassanzadeh et al.(2014)] where the optimal aspect ratio scales like $Ra^{-1/4}$ and so vanishes as $Ra\rightarrow\infty$.

Figure 6 shows the structure of the two modes at $Ra=10^{4}$. The fluctuation fields $\psi_{i}$ and $\theta_{i}$ for both $i=1$ and $2$ have a convection roll structure and increasing $L$ just means that more of the rolls fit into the domain. On closer inspection it is clear that the rolls are slightly different near to $x=0$ and $x=\pm\frac{1}{2}L$) where they are forced to have a certain symmetry (symmetry around $x=0$ and periodicity over a length $L$ force symmetry about $x=\pm\frac{1}{2}L$ as well). When the domain is short, e.g. $L=\pi$, the background field is clearly two-dimensional as seen in figure 7. However, as $L$ increases to $L=8\pi$, the background field become predominantly one-dimensional ($1$-$D$) away from the imposed lines of symmetry at $x=0$ and $x=\pm\frac{1}{2}L$ (the ends of the domain shown). Plotting the streamfunctions $\psi_{1}$ and $\psi_{2}$ over this long domain - see figure 8 - confirms that the convection cells are similar away from the symmetry lines (‘zone 1’ in figure 8) where $\tau$ is predominantly $1$-$D$ but are quite different close to the symmetry lines (‘zone 2’) where $\tau$ is clearly $2$-$D$.

The structure of the optimal fields (both background and fluctuation) and the fact that the bound is maximised as $L\rightarrow\infty$ indicate that the optimal solution is trying to minimise the effect of the imposed symmetry requirements at $x=0$ and $x=\pm\frac{1}{2}L$. Without this imposed symmetry, the problem becomes translationally invariant and the optimal solution must be $1$-$D$ by the convexity result in §2.3. There is another simple way to see this. Since the bounding functional $f(\tau,a)$ (see (29) ) is strictly convex in both $\tau(x,z)$ and $a$, any $2$-$D$ solution $(\tau_{2D}(x,z),a)\in\Omega$

[TABLE]

by Jensen’s inequality. Taking the limit of $N\rightarrow\infty$ in the left hand side and using translational invariance of the problem in the right hand side leads to

[TABLE]

so that a $1$-$D$ background field always produces a better bound than a $2$-$D$ field. The results in figure 7 indicate that this is what the optimal solution is trying to achieve.

5.1 Lifting the symmetry: 1D background field

Calculation of the optimal solution assuming from the onset that the background field is $1$-$D$ simplifies the computation since the fluctuation fields can then be parametrised by their single wavenumber in $x$ (as in (37)). In this case, rather than setting a domain periodicity and insisting the fluctuation wavenumbers be consistent with this, the wavenumbers themselves can be optimised over as real continuous variables meaning, in effect, that $L$ is infinite. For example, the Euler-Lagrange equation corresponding to the $m^{th}$ wavenumber $k_{m}$ is

[TABLE]

With this formulation, Newton’s method with branch continuation proved much faster than the time-stepping approach. It took around $4$ hours cputime on a 2.6Ghz laptop using Newton’s method to obtain the optimal solution from $Ra=27\pi^{4}/4$ up to $Ra=5\times 10^{8}$ while the time-stepping approach took at least a day to generate a single point at $5\times Ra=10^{8}$. However, when the domain is fixed, Newton’s method becomes very inefficient as the critical wavenumbers $k_{m}$ are discrete and cannot be tracked using a (continuous) continuation method: in this case, time-stepping is the better choice. The numerical solution of the one-dimensional background problem gives the upper bound of $Nu\leq 0.055Ra^{1/2}$ as shown in figure 9(a) with 5 critical modes present by $Ra=10^{9}$ ( see figure 9(b) ). This result has the same scaling exponent as the non-slip result $Nu\leq 0.026Ra^{1/2}$ of PK03 but with a larger numerical coefficent as should be expected for stress-free boundary conditions. The prior work of [Wen et al. (2015)] indicates that adding a further enstrophy constraint (possible only in stress-free 2D convection) significantly improves the bound obtained here down to $Nu\leq 0.106Ra^{5/12}$.

It is worth briefly discussing how the critical wavenumbers which appear in the 1-D and 2-D background field calculations are related. Figure 9(b) indicates that at $Ra=10^{4}$ there is only one critical wavenumber $k_{1}=3.284$ for the 1-D background problem. However, for the 2-D (symmetric) background problem, there are two critical modes as seen in figure 6 and figure 8 with both having an approximate wavenumber $\approx 3.3$. Both these modes are forced to be antisymmetric (and so in phase) about $x=0$ and $x=\pm L/2$ (zone 2 in figure 8) but away from these points endeavour to be approximately $\pi/2$ out of phase (zone 1 in figure 8). With this phase difference together with matching amplitudes so

[TABLE]

the nonlinear term in Eq.(27)

[TABLE]

is $x$ independent and therefore can only drive a 1-D background field. From another perspective, the 1-D background field problem really has two modes with $k=3.284$ but only one needs to be tracked as the nonlinear term is horizontally averaged (e.g. see (35) ) ensuring that the background field stays 1-D.

At $Ra=25,000$ there is even a third mode in the 2-D background problem compared to still only 1 mode in the 1-D background problem (the second wavenumber $k_{2}$ appears at $Ra\approx 26,450$). Figure 10 shows that in fact $\psi_{3}$ is only significant in zone 2 where the imposed symmetries dominate. In zone 1 where the background field is essentially 1-D profile, $\psi_{3}$ vanishes.

The conclusion of the computations so far is that imposing the full heat equation in the bounding calculation does not improve (lower) the bound over that obtained using the horizontally-averaged heat equation. The next obvious question is whether this is also true when imposing the full momentum equation. The next section addresses this.

6 Imposing the full momentum equation: \mbox{\boldmath\phi}\neq{\bf 0}

In this section, we attempt to improve the bound by using a background temperature field and a background velocity field of the same dimension as the physical problem which means that the full momentum equation and heat equation are imposed as constraints. Importantly, the optimization problem is no longer convex and so we are unable to invoke uniqueness to dismiss $\phi$. Instead we use an inductive bifurcation analysis to show that if \mbox{\boldmath\phi}={\bf 0} before a bifurcation then it remains ${\bf 0}$ after it too meaning that the continuous branch of optimals found by branch tracking out of the energy stability point always has \mbox{\boldmath\phi}={\bf 0}.

The analysis begins by constructing the following Lagrangian:

[TABLE]

which, after introducing the extended background decomposition

[TABLE]

can be rewritten as

[TABLE]

where

[TABLE]

(note $\mathscr{G}$ depends parametrically on $\tau$, $\boldsymbol{\phi}$, $a$, $\sigma$ and $Ra$ but this is suppressed for clarity). If $\inf_{\boldsymbol{v},\theta}\mathscr{G}$ exists (and necessarily $0<a<1$), a bound is then given by

[TABLE]

The difficulty here is that the objective functional is no longer convex and so it’s unclear how to establish a priori that the optimal solution takes the form (\tau,\mbox{\boldmath\phi})=(\tau(z),{\bf 0}). Hence we consider what has to happen at bifurcation points. Minimization of $\mathscr{G}$ with respect to incompressible $\boldsymbol{v}$ and $\theta$ requires

[TABLE]

[TABLE]

[TABLE]

the solution of which is denoted as $(\boldsymbol{v}_{0},\theta_{0})$. The Lagrangian can then be written as

[TABLE]

where $\boldsymbol{v}_{0}$, $\boldsymbol{v}_{i}$ ($i=1\ldots N$) and $\boldsymbol{\phi}$ are incompressible fields and

[TABLE]

is a purely quadratic functional of $(\boldsymbol{v},\theta)$ which must be positive semi-definite - the spectral constraint - for $\inf\mathscr{G}$ to exist. The fields $(\boldsymbol{v}_{i},\theta_{i})$ are marginal in that $\mathscr{H}(\boldsymbol{v}_{i},\theta_{i})=0$ and their number $N$ increases with $Ra$. The aim is to minimize the upper bound over $\tau$, $\boldsymbol{\phi}$ and $a$ at fixed $\sigma$ and $Ra$ subject to this spectral constraint. The Euler-Lagrange equations are: the spectral constraint equations for $\boldsymbol{v}_{i}$ and $\theta_{i}$

[TABLE]

[TABLE]

the forced field equations for $\boldsymbol{v}_{0}$ and $\theta_{0}$

[TABLE]

[TABLE]

the background field equations

[TABLE]

[TABLE]

and finally the balance parameter equation

[TABLE]

The pressure-like quantities have been rescaled $Ra\,p_{i}\rightarrow p_{i}$, $Ra\,p_{0}\rightarrow p_{0}$ and $(1-a)Ra\,q\rightarrow q$) and incompressibility conditions on $\boldsymbol{v}_{0}$, $\boldsymbol{v}_{i}$ and $\boldsymbol{\phi}$ are left implicit. A key point here is that the forced field pair $(\boldsymbol{v}_{0},\theta_{0})$ is not marginal in the spectral constraint i.e. $\mathscr{H}(\boldsymbol{v}_{0},\theta_{0})>0$.

6.1 The first bifurcation point

The solution at the first critical point $Ra_{c}=27\pi^{4}/4$ is $\tau=1-z$ and $\boldsymbol{\phi}=0$, $(\boldsymbol{v}_{0},\theta_{0})=({\bf 0},0)$ and $a=1$. At $Ra=Ra_{c}$, the spectral constraint becomes marginal for the first time, i.e. there is a non-trivial solution to the spectral problem

[TABLE]

[TABLE]

[TABLE]

There are two different modes (using symmetries):

[TABLE]

[TABLE]

Since $d\tau/dz=-1$, the structure in $z$ is simple: $U:=\pi\cos(\pi z)$, $W:=-k\sin(\pi z)$ and $\Theta:=-\frac{\sqrt{2}}{3\pi}\sin(\pi z)$ where $k=\pi/\sqrt{2}$. Slightly away from the critical point, $Ra=Ra_{c}+\varepsilon$, the fields need to be expanded as follows

[TABLE]

where $\tau_{0}:=1-z$ and $a_{0}:=1$.

6.1.1 Leading Order

To leading order, the spectral constraint is satisfied by $(\boldsymbol{v}^{0}_{i},\theta^{0}_{i})$ defined in (73-74) and these fields force the other leading order equations for the background fields

[TABLE]

which are coupled with the forced field equations

[TABLE]

[TABLE]

The forcing term in (75) is

[TABLE]

and in (76)

[TABLE]

(recall $k=\pi/\sqrt{2}$) so that simply

[TABLE]

[TABLE]

where $c_{1},c_{2}$ and $c_{3}$ are specific constants. Finally, the leading order balance (which is at $O(\varepsilon)$) in the balance parameter equation (69) is

[TABLE]

which relates $A_{1}^{2}+A_{2}^{2}$, $a_{1}$ and the higher order unknown $\boldsymbol{\phi}_{2}$.

6.1.2 Next Order

A further piece of information to identify the leading order fields comes from a solvability condition on the spectral constraint equations at next order ($(O(\varepsilon^{3/2})$) which is

[TABLE]

[TABLE]

Formally, this has two solvability conditions:

[TABLE]

but they are equivalent since $\tau_{1}$ is 1-dimensional (i.e. solely a function of $z$) with the resulting condition linking $a_{1}$ and $A_{1}^{2}+A_{2}^{2}$. Only after $\boldsymbol{\phi}_{2}$ is found can $a_{1}$ and $A_{1}^{2}+A_{2}^{2}$ be fully determined from (83) and (86).

The fields $(\boldsymbol{v}_{i}^{1},\theta_{i}^{1})$ are linearly dependent on $\boldsymbol{v}_{i}^{0}$, $\theta_{i}^{0}$ and so can be written as

[TABLE]

[TABLE]

where $\mathcal{U}\neq U$, $\mathcal{W}\neq W$. These fields along with $(\boldsymbol{v}_{i}^{0},\theta_{i}^{0})$ drive the higher order equations governing further corrections to the background fields and the forced fields. These are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The apparent driving term $-a_{1}/\sigma\sum_{i=1}^{2}\boldsymbol{v}_{i}^{0}\cdot\boldsymbol{\nabla}\boldsymbol{v}_{i}^{0}$ can be balanced by the pressure term (see (80) ) as can $([a_{1}Ra_{c}+1]\tau_{1}+a_{1}\tau_{0})\boldsymbol{e}_{z}$. Also importantly for what follows, $Ra_{c}a_{2}\tau_{0}\boldsymbol{e}_{z}$ in (90) can also be absorbed into the pressure term which means that $\tau_{2}$ and $\boldsymbol{\phi}_{2}$ do not depend on $a_{2}$ (this is crucial for the argument surrounding (102) below). This leaves the driving term for the 2D background temperature field

[TABLE]

and the driving term in Eq.(90) for $\boldsymbol{v}_{0}^{2}$

[TABLE]

Given the form of these driving terms, $\tau_{2}$ can be split into two parts: a 1D part which depends only on $z$ proportional to $A_{1}^{2}+A_{2}^{2}$, and a 2D part proportional to $A_{1}^{2}-A_{2}^{2}$ which has both $x$ and $z$ dependence whereas the remaining corrections $\boldsymbol{\phi}_{2}$, $\boldsymbol{v}_{0}^{2}$ and $\theta_{0}^{2}$ only have a 2D part proportional to $A_{1}^{2}-A_{2}^{2}$, so

[TABLE]

(the expressions for $\boldsymbol{v}_{0}^{2}$ and $\theta_{0}^{2}$ are not needed in what follows and hence suppressed). At this point $\boldsymbol{\phi}_{2}$ is now known as a function of $A_{1}^{2}-A_{2}^{2}$ and even with the previously derived relations (83) and (86), it is still not possible to identify $A_{1}$, $A_{2}$ and $a_{1}$ without further information about $A_{1}$ and $A_{2}$. This comes from solvability conditions at the next order of the spectral constraint ($O(\varepsilon^{5/2})$).

Before pursuing this, we remark that the next order ($O(\varepsilon^{2})$) of the balance equation (69) involves the higher order unknown $\boldsymbol{\phi}_{3}$ and so at this order $a_{2}$ is unspecified. In fact $a_{2}$ is also set by solvability conditions at $O(\varepsilon^{5/2})$ of the spectral constraint to which we now turn.

6.1.3 Solvability at $O(\varepsilon^{5/2})$

The spectral constraint equations at $O(\varepsilon^{5/2})$ are

[TABLE]

[TABLE]

Identifying the operator on the left hand side of (97)$-$(98) as $\mathcal{L}$, then since it is self adjoint and annihilates the leading order fields $(\boldsymbol{v}_{j}^{0},\theta_{j}^{0})$ $j=1,2$, there are solvability conditions for $(\boldsymbol{v}_{0}^{2},\theta_{0}^{2})$

[TABLE]

Taking $i=j$ (the $i\neq j$ conditions vanish trivially), this can be rearranged to

[TABLE]

where

[TABLE]

Crucially Term${}_{1}(1)$=Term${}_{1}(2)$ whereas Term${}_{2}(1)$=$-$Term${}_{2}(2)$ so that (100) implies that

[TABLE]

The unspecified coefficient $a_{2}$ only enters Term1 and so is set by the condition this vanishes. In contrast, there are no free constants in Term2 which is non-zero in our computations (although we have been unable to prove this is always the case). In this situation, $A_{1}^{2}-A_{2}^{2}=0$ is forced instead which eliminates at a stroke all 2-dimensional fields in the bifurcation analysis. Consequently, the background flow field remains zero and the background temperature field stays 1D after the first bifurcation.

6.2 Subsequent bifurcations

Now we consider subsequent bifurcations to establish that if $\tau=\tau(z),\boldsymbol{\phi}=0$ exists before then that situation persists after the bifurcation. The approach is inductive: assume $\tau=\tau(z),\boldsymbol{\phi}=0$ after $m$ bifurcations and consider the $(m+1)^{th}$ bifurcation at $Ra=Ra_{c}^{(m+1)}$ where two new neutral modes appear so that there are now $2(m+1)$ critical modes in the spectral constraint. Defining $\varepsilon:=Ra-Ra_{c}^{(m+1)}$ we expand:

[TABLE]

[TABLE]

where the leading fields $\tau_{0}(z)$, $(\boldsymbol{v}_{i}^{0},\theta_{i}^{0})$ ($i=1,\ldots 2m+2$) and $a_{0}$ are all known. In particular, the $i^{th}$ wavenumber $k_{i}$ ($i=1,2,\dots,m$), is associated with two modes which, to leading order, are

[TABLE]

and

[TABLE]

where $A_{i}^{2}=B_{i}^{2}$ for $i=1,2,...,m$. The two new modes emerging at the $(m+1)^{th}$ bifurcation point can be assumed to have the following general 3D form:

[TABLE]

and

[TABLE]

where $\boldsymbol{k}_{m+1}=(k_{x},k_{y},0)$. The objective in what follows is to show that $A_{m+1}^{2}=B_{m+1}^{2}$ after the bifurcation so that the optimization problem remains 1D.

At leading order ($O(\epsilon)$) in the forced field and background field equations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where again it is implicit that $\boldsymbol{v}_{0}$, $\boldsymbol{v}_{i}$ and $\boldsymbol{\phi}$ are incompressible fields. The spectral constraint for modes $i=1,2,...,2m$ at $O(\varepsilon)$ and for modes $i=2m+1,2m+2$ at $O(\varepsilon^{3/2})$ is

[TABLE]

[TABLE]

The system of equations (119)-(124) is linear in $\boldsymbol{v}_{0}^{0}$, $\theta_{0}^{0}$, $\boldsymbol{\phi}_{0}$, $\boldsymbol{v}_{i}^{1}$ and $\theta_{i}^{1}$ ($i=1,2,...,2m+2$). The emergent critical modes at $Ra_{c}^{m+1}$ give rise to the new driving terms in (121) & (122)

[TABLE]

which have a 1D part proportional to $A_{m+1}^{2}+B_{m+1}^{2}$ and a non-1D part proportional to $A_{m+1}^{2}-B_{m+1}^{2}$. If $A^{2}_{m+1}=B^{2}_{m+1}$ then $\tau_{1}=\tau_{1}(z)$ and $\boldsymbol{\phi}={\bf 0}$ using the arguments of section 6.1. Hence for $A^{2}_{m+1}\neq B^{2}_{m+1}$ and using equation (122), we can assume a solution structure of the form

[TABLE]

where $\tau_{1}^{*}(x,z)$ and $\boldsymbol{\Phi}^{*}(x,z)$ collect all the other wavenumber dependence on $x$ in $\tau_{1}$ and $\boldsymbol{\phi}_{0}$ respectively (this is unimportant in what follows). The key is then examining the solvability conditions

[TABLE]

($i,j,=1,\ldots m+1$) on the corrections $\boldsymbol{v}_{i}^{1}$ to all the critical modes $\boldsymbol{v}_{i}^{0}$ of the spectral constraint. These set the amplitudes $A_{i}^{2}\,(=B_{m}^{2})$ ($i=1,\ldots m$) and $(A_{m+1},B_{m+1})$ ($a_{1}$ is determined by the balance parameter equation (69) at $O(\varepsilon)$).

To establish that the background fields stay 1D, it is sufficient to focus on the solvability conditions for the new critical modes ($i=2m+1$ and $2m+2$). Here, the solvability condition is explicitly

[TABLE]

where

[TABLE]

with $c:=A_{m+1}^{2}$ for $i=2m+1$ or $B_{m+1}^{2}$ for $i=2m+2$. Crucially, it is straightforward to show that

[TABLE]

so that, as in subsection (6.1), we must have

[TABLE]

Our computations indicate ${\rm Term}(2m+1)\neq 0$ (although, as in subsection (6.1), we have been unable to prove this in general) forcing $A_{m+1}^{2}=B_{m+1}^{2}$. This forces $\tau_{1}=\tau_{1}(z)$ and $\boldsymbol{\phi}_{0}={\bf 0}$ so that the optimal solution stays 1D after the $(m+1)^{th}$ ($m\geq 0$) bifurcation if it is of this form before. Taken together with the first bifurcation analysis in subsection (6.1), this means that the optimal solution is 1D for all $Ra$ and so, surprisingly, there is no benefit of imposing the full momentum and heat balances in the upper bound problem.

7 Adding a background velocity field to [Wen et al. (2015)]

Following the success of adding an enstrophy constraint in 2D stress-free convection [Wen et al. (2015)], an interesting question is whether adding a 1-D background velocity field by using the decomposition

[TABLE]

would improve the bound further since this imposes additional information from the Navier-Stokes equations. To maximize the heat flux, the Lagrangian

[TABLE]

is constructed where $\boldsymbol{\omega}=\omega(x,z)\boldsymbol{e}_{y}:=(\boldsymbol{\nabla}\times\boldsymbol{v})$. After some integration by parts, judicious use of boundary conditions, and building in the fact that $Nu=1+\langle wT\rangle$, this leads to the expression

[TABLE]

The two linear terms in the second line of this expression - $v_{1}\phi^{\prime\prime}$ and $\omega\phi^{\prime\prime\prime}$ - mean that optimization over the fluctuation fields $\boldsymbol{v}$ and $\theta$ will give rise to a non-zero contribution to be added to $\langle|\tau^{\prime}|^{2}\rangle$. This complication can be avoided (or rather made more explicit) by defining a shifted variable

[TABLE]

which is possible if $\boldsymbol{v}$ and $\phi$ are both assumed to satisfy (natural) homogeneous boundary conditions and allows both linear terms to be absorbed into perfect squares. As a result of this, the expression becomes

[TABLE]

where

[TABLE]

is a purely quadratic functional of $\hat{\boldsymbol{v}}$, $\hat{\omega}:=\boldsymbol{e}_{y}\cdot\boldsymbol{\nabla}\times\hat{\boldsymbol{v}}=\omega+{\textstyle{1\over 2}}\phi^{\prime}$ and $\theta$. Provided the background fields $\tau$ and $\phi$ are chosen such that $\mathscr{G}\geq 0$ for all permissable $\hat{\boldsymbol{v}}$, $\hat{\omega}$ and $\theta$, then a bound follows on $Nu$. Now it is clear that: 1) the objective functional is convex in the background fields and 2) the set of allowable background fields is convex (if $(\tau_{1},\phi_{1})$ and $(\tau_{2},\phi_{2})$ ensure $\mathscr{G}\geq 0$ so does $(\tau,\phi)=\lambda(\tau_{1},\phi_{1})+(1-\lambda)(\tau_{2},\phi_{2})$ with $0\leq\lambda\leq 1$). This implies that the optimizer is unique and is attained for $(\tau,\phi)=(\tau,0)$ i.e. the background velocity field vanishes indicating that the extra information this folds into the optimization is, in fact, unimportant. Physically, a bifurcation analysis shows that the fluctuation fields are always such as to produce zero Reynolds stress so that no background field is generated.

A numerical solution shown in figure 11 using Newton’s method on the Euler-Lagrange equations in an infinitely long domain confirms that $\phi=0$ as does a bifurcation analysis developed in the previous section. The bound compares well with the earlier results of [Wen et al. (2015)] who considered a fixed domain of $L=2\sqrt{2}$ indicating further that the bound is not that sensitive to the domain size. The fashion in which the necessarily discretized critical modes found by [Wen et al. (2015)] cluster around the (continuous) optimal wave numbers in our study confirms this conclusion.

8 Discussion

This paper has revisited the optimal heat transport problem in two-dimensional Rayleigh-Bénard convection with stress-free boundary conditions using an extended background method. The key novelty has been to consider background temperature and velocity fields whose dimensional dependence matches that of the physical problem (so 2D here). This situation needs a reformulation in the way the variational equations are solved which has the significant consequence of breaking any link between the optimal fields which emerge and a single physical temperature and velocity field. In particular, this means that the optimal fields do not obviously satisfy the heat equation even though that is explicitly imposed when the background temperature field is allowed to be fully 2D. This is due to the spectral constraint (that ensures a bound) which means the optimal bound found does not correspond with the highest stationary point of the Lagrangian (i.e. the Euler-Lagrange equations are not all satisfied) but is strictly above it. In other words, there is a gap between the highest heat flux attained by a steady solution of the governing equations imposed and the best (lowest) bound because of the additional spectral constraint. Unfortunately and importantly, this means that there is no direct connection between the optimal solution in the background method built around the steady governing equations and a steady solution of the governing equations (here the Boussinesq equations but clearly more generally true). This realisation removes the possibility, for example, that the simple 2D roll solution computed by Waleffe et al. (2015) could actually be the optimal solution to the background bounding problem. It now seems clear that it would be spectrally unstable.

In revisiting the exact 2D Rayleigh-Benard problem treated by Hassanzadeh et al. (2014), we have shown that their maximal heat flux result is only a global maximum up to $Ra\leq Ra_{c}:=4468.8$. Beyond this, the optimal solution complexifies over a given spatial domain. If this domain is extended, the optimal becomes increasingly 1D. Removing the symmetry imposed by Hassanzadeh et al. (2014) and reinstating translational invariance in the horizontal direction by making the domain unbounded, the optimal solution is then provably just 1D and the classic scaling result of $Nu\sim Ra^{1/2}$ is recovered albeit with the larger numerical coefficient of 0.055 as opposed to the already known 0.026 (PK03) for non-slip boundary conditions. The conclusion is then that imposing the full heat equation in the bounding calculation does not improve (lower) the bound over that obtained using the horizontally-averaged heat equation. We then considered adding extra information from the momentum equation to the upper bounding problem by introducing a background velocity field \mbox{\boldmath\phi}(x,z). Now the optimization problem is no longer convex and we use an inductive bifurcation analysis to show that if \mbox{\boldmath\phi}={\bf 0} before a bifurcation then it remains ${\bf 0}$ after it too. This means that the continuous branch of optimals found by branch tracking out of the energy stability point always has \mbox{\boldmath\phi}={\bf 0}. Noting the caveats that a) it’s not impossible that there is an unconnected branch of optimals with \mbox{\boldmath\phi}\neq{\bf 0} and b) $Term(2m+1)$ in (139) could serendipitiously vanish at a subsequent bifurcation beyond our calculations, this strongly suggests the surprising result that imposing the full Boussinesq equations does not improve the bound over that obtained using the horizontally-averaged Boussinesq equations.

The ‘take-home’ message from this study is that the background method of seeking an upper bound on heat flux in Rayleigh-Benard convection has been exhausted with disappointingly no improvement possible over the minimal choice of a 1D background temperature field originally made in 1996 by Doering and Constantin. It’s hard not to imagine this realisation also generalising to the analogous background formulations for shear flows too, e.g. plane Couette flow (Doering & Constantin, 1992), channel flow (Constantin and Doering, 1995) and pipe flow (Plasting & Kerswell, 2005). Simply extending the definition of the background fields ostensibly folds in more information from the governing equations but not in a fruitful way. However, it seems generating extra information by differentiating the governing equations can help. Whitehead & Doering (2011) (see also Wen et al. 2015) used an extra vorticity constraint to significantly lower the bound from $Nu\sim Ra^{1/2}$ to $Nu\sim Ra^{5/12}$ but only in the 2D situation with stress-free boundary conditions. Interestingly, this approach can be inverse-engineered into the form of background method by loosening the connection between the Lagrange multiplier \mbox{\boldmath\nu}(\mathbf{x},t) and the velocity field $\mathbf{u}(\mathbf{x},t)$ from \mathbf{u}(\mathbf{x},t)-\mbox{\boldmath\nu}(\mathbf{x},t)=\phi(z)\hat{\mathbf{x}} to

[TABLE]

where $c$ is a new scalar Lagrange multiplier imposing the global vorticity constraint

[TABLE]

This clearly extracts something more from the governing equations than just taking projections. Maybe there is some mileage in exploring this but a shortage of boundary conditions is the usual impediment to this approach. Looking ahead, an emergent Sum-of-Squares approach to bounding (e.g. Fantuzzi et al. 2016, Goluskin & Fantuzzi 2019) offers far greater potential for progress since it extends the quadratic constraints used here to more general polynomials albeit at the expense of a fully numerical approach.

Acknowledgements: The authors are very grateful to Andre Souza and Charlie Doering for helpful discussions and sharing their recent preprint (Souza et al. 2019). The authors also acknowledge the support of EPSRC under grant EP/P001130/1.

Appendix A Time stepping for a 2D background temperature field in 2D

Here we show that the time-marching method of Wen et al. (2015) is not guaranteed to have the optimal solution as the unique steady attractor when $\tau=\tau(x,z)$ has the same spatial dimensionality as the physical temperature field $T(x,z,t)$. Time-stepping would work, however, for a 2-dimensional $\tau(x,z)$ in a 3-dimensional problem. To explain this, we revisit the proof of Wen et al. (2015). The time-stepping approach consists of adding time derivatives for $\theta$, $\nabla^{2}\psi$ and $\tau$ to the left hand sides of (25)-(27) respectively. Small disturbances $(\theta^{\prime},\boldsymbol{u}^{\prime},\tau^{\prime},p^{\prime})$ on top of a solution to the Euler-Lagrange equations, $(\theta,\boldsymbol{u},\tau,p)$, then evolve according to the following equations

[TABLE]

at fixed balance parameter $a$. Then $\langle\,\theta^{\prime}(\ref{T_1})-\psi^{\prime}(\ref{T_2})+\tau^{\prime}(\ref{T_3})\,\rangle$ gives

[TABLE]

In the 1-dimensional background field case, $\tau=\tau(z)$, (147) becomes

[TABLE]

(where the overbar represents averaging over $x$) and the possible fluctuation fields can, after a Fourier transform, be assumed to have a specific wavenumber in $x$. There are then two types of fluctuation fields: 1) those with wavenumbers which don’t overlap with those in the optimal solution ($\lambda<0$ in the spectral constraint) and therefore do not generate any concomitant disturbance $\tau^{\prime}$, and 2) those which do have a non-vanishing $\tau^{\prime}$ but necessarily have $\lambda=0$ (the optimal solution is unique for any given balance parameter $a\in(0,1)$ by the same arguments presented in the main text and, by construction, includes any fluctuation fields $(\theta,\psi)$ which are neutral in the spectral constraint). In both cases, the fluctuation fields have to decay, in the former case because $\lambda<0$ and in the latter through the $\tau^{\prime}$ component generated in (149). The unique solution is therefore an attractor but the key step is proving that it is the only such. This follows by realising that if a solution to the Euler-Lagrange equations does not satisfy the spectral constraint, then there is an unstable eigenfunction of the linear time-stepping operator defined in (145)-(147) which consists of the fluctuation field $(\theta^{\prime},\psi^{\prime})$ which makes $\mathscr{G}<0$. This is because a fluctuation field with $\lambda\neq 0$ does not overlap under $x$-averaging with the underlying state and so does not generate a $\tau^{\prime}$ component via (149). This argument can clearly be extended to 2-dimensional $\tau(x,z)$ in 3-dimensional Rayleigh-Benard convection since orthogonality in $x$ is replaced by orthogonality of $y$ but breaks down for 2-dimensional Rayleigh-Benard convection. In the latter situation, fluctuation fields which violate the spectral constraint will generate a $\tau^{\prime}$ component via (147) and may not then represent a growing eigenfunction for the time stepping procedure. The implication of this is that some saddles of ${\mathscr{L}}$ may also be local attractors so if the time-stepping procedure leads to a steady state it is not guaranteed to be the optimal solution. Preliminary numerical tests demonstrated this multistability with the final steady state depending on the initial condition used.