Wall-to-wall optimal transport in two dimensions

TL;DR

This paper develops gradient ascent methods to find two-dimensional incompressible flows that maximize heat transfer between walls, revealing a transition in scaling laws and near-separable optimal flow structures up to very high Péclet numbers.

Contribution

It introduces a novel gradient ascent approach to compute optimal flows for heat transfer, demonstrating a change in scaling behavior and near-separable flow structures at high Péclet numbers.

Findings

Transport scaling transitions from Pe^2 to Pe^0.54

Optimal flows are approximately separable functions

Derived an upper bound consistent with computational results

Abstract

Gradient ascent methods are developed to compute incompressible flows that maximize heat transport between two isothermal no-slip parallel walls. Parameterizing the magnitude of velocity fields by a P\'eclet number $\text{Pe}$ proportional to their root-mean-square rate-of-strain, the schemes are applied to compute two-dimensional flows optimizing convective enhancement of diffusive heat transfer, i.e., the Nusselt number $\text{Nu}$ up to $\text{Pe} \approx 10^5$. The resulting transport exhibits a change of scaling from $\text{Nu}-1 \sim \text{Pe}^{2}$ for $\text{Pe} < 10$ in the linear regime to $\text{Nu} \sim \text{Pe}^{0.54}$ for $\text{Pe} > 10^3$. Optimal fields are observed to be approximately separable, i.e., products of functions of the wall-parallel and wall-normal coordinates. Analysis employing a separable ansatz yields a conditional upper bound $\lesssim \text{Pe}^{6/11} = \text{Pe}^{0.\overline{54}}$ as $\text{Pe} \rightarrow \infty$ similar to the computationally achieved scaling. Implications for heat transfer in buoyancy-driven Rayleigh-B\'enard convection are discussed.

Equations

$${\partial }_{t}T+\nabla \cdot \left(\mathbf{u}T-\kappa \nabla T\right)$$

$${\partial }_{t}T+\nabla \cdot \left(\mathbf{u}T-\nabla T\right)=0,$$

$$\text{Nu}\{\mathbf{u}\}$$

$$⟨f⟩$$

$${\text{Pe}}^2$$

$$\text{Maximise}⟨u_{3}T⟩\text{subject to }$$

$${\partial }_{t}T+\mathbf{u}\cdot \nabla T=\Delta T,$$

$$\nabla \cdot \mathbf{u}=0,⟨|\nabla \mathbf{u}{|}^2⟩={\text{Pe}}^2,$$

$$\text{and the boundary conditions}.$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}+\nabla p=\text{Pr}\Delta \mathbf{u}+\text{Pr}\text{Ra}T\hat{z},$$

$$\nabla \cdot \mathbf{u}=0,\text{ and}$$

$${\partial }_{t}T+\mathbf{u}\cdot \nabla T=\Delta T.$$

$$Pr=\nu /\kappa \text{ and }\text{Ra}=\frac{\alpha g(T_1-T_0)(L_z{)}^3}{\kappa \nu }$$

$${\partial }_{\tau }\xi$$

$$0$$

$${\partial }_{\tau }\mathbf{u}$$

$$0$$

$${\partial }_t\theta +\mathbf{u}\cdot \nabla \theta =\Delta \theta +u_3.$$

$$\mathcal{F}$$

$$\mathcal{S}$$

$$\Delta \eta$$

$${\mathcal{S}}_{\eta }$$

$$\Delta \xi$$

$${\mathcal{S}}_{\xi }$$

$${\mathcal{S}}_{\mathbf{u}}$$

$$\mu \Delta \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$${\mathcal{S}}_{\mathbf{u}}$$

$$\Delta \xi$$

$$\Delta \eta$$

$$\Delta \xi$$

Full text

1 Introduction

Heat transfer via fluid advection is a critical component of atmospheric, oceanographic, geophysical, and astrophysical dynamics, as well as being the basis of cooling systems in engineering applications. Numerous studies on how to design systems that achieve enhanced heat transfer by either manipulation of domain geometry or through the discovery of suitable flow structures have recently appeared in the literature Toppaladoddi et al. (2017); Alben (2017); Marcotte et al. (2018); Motoki et al. (2018a, b). A particularly fruitful approach to discovering flow structures was first introduced by Hassanzadeh et al. (2014) where it was formulated via an optimal control approach.

The original motivation for Hassanzadeh et al. (2014) was to develop a new tool for obtaining upper bounds on thermal transport by buoyancy-driven flows, e.g., for Rayleigh-Bénard convection. The analysis and derivation of upper bounds on transport properties plays a prominent role in expanding our knowledge of fundamental fluid dynamics and serves a complimentary role to other methods of inquiry, i.e., direct numerical simulations of the underlying equations of motion, invoking closure models to determine statistical properties, or postulating phenomenological models of turbulent transport. The first proof of upper bounds on heat transfer by Rayleigh-Bénard convection was achieved in Howard (1963). The complementary “background method” was subsequently introduced in Doering & Constantin (1996). Both approaches leverage certain bulk integral constraints derived from the equations of motion and yield bounds which apply to a strictly larger class of flows. It remains unknown whether the resulting bounds are realizable by buoyancy driven flows.

Unlike those previous approaches, the wall-to-wall transport problem introduced in Hassanzadeh et al. (2014) fully enforces the advection-diffusion equation for the temperature field pointwise in space and time. Admissible incompressible advecting flow fields do not (necessarily) satisfy an equation of motion, but are instead subject to a bulk integral intensity constraint, i.e., fixed finite magnitude of energy or enstrophy, and suitable boundary conditions. This allows consideration of the following question: Amongst all possible incompressible flow fields subject to a fixed intensity constraint and relevant boundary conditions, which ones maximize thermal transport?

As usual, we model heat transport with the advection-diffusion equation

[TABLE]

where the coefficient $\kappa$ is the thermal diffusivity. The two-dimensional spatial domain is $\Omega=[0,L_{x})\times[0,L_{z}]$ and the temperature field $T$ is periodic in the horizontal $x$ direction, “hot” on the bottom boundary where $T(z=0)=T_{0}$ and “cool” on the top boundary where $T(z=L_{z})=T_{1}$ with $T_{0}>T_{1}$. The advecting flow field $\mathbf{u}=u_{1}\hat{x}+u_{3}\hat{y}$ is divergence-free ($\nabla\cdot\mathbf{u}=0$) with no penetration through the boundaries, i.e., $u_{3}(z=0)=0=u_{3}(z=L_{z})$. In this paper the velocity is restricted to satisfy no-slip boundary conditions on the top and bottom boundaries, $u_{1}(z=0)=0=u_{1}(z=L_{z})$. Both components are periodic in the horizontal $x$ direction. Other boundary conditions can be handled using similar methods. Initial ($t=0$) data for the temperature field is provided to formally pose the evolution problem for $t>0$. Using units $L_{z}$ and $L_{z}^{2}/\kappa$ for space and time and changing the temperature $T\rightarrow\frac{T-T_{0}}{T_{1}-T_{0}}$ transforms the system to

[TABLE]

with $x\in[0,\Gamma)$ where $\Gamma=L_{x}/L_{z}$, $z\in[0,1]$, and with $T(z=0)=1$ and $T(z=1)=0$. We consider (2) henceforth.

Given a flow field $\mathbf{u}$ the non-dimensional measure of thermal transport is the space- and long time-average of the convective heat flux in the vertical direction, i.e., the Nusselt number

[TABLE]

where $\langle\cdot\rangle$ denotes the space-time average

[TABLE]

The boundary conditions imply unit mean conductive heat flux, i.e., $\langle-\partial_{z}T\rangle=1$.

In this work, the goal is to determine the largest possible value of Nu as a function of $\mathbf{u}$, which we parameterize by a Péclet number defined as the root-mean-square vorticity (the square root of the mean enstrophy density) $\text{Pe}\equiv\langle|\nabla\times\mathbf{u}|^{2}\rangle^{1/2}$. Incompressibility and the boundary conditions imply that Pe is simply related to the norm of $\nabla\mathbf{u}$ and the mean-square rate of strain:

[TABLE]

The optimal wall-to-wall transport problem is then to maximize Nu as a function of Pe. Explicitly, we take on the task to

[TABLE]

Note in the wall-to-wall problem the velocity field $\mathbf{u}$ is not required to satisfy conservation of momentum. Nevertheless, in Hassanzadeh et al. (2014) its solution was shown to inform bounds on buoyancy-driven driven transport; indeed, the original motivation for introducing the wall-to-wall problem was to find new upper bounds on the Nusselt number in Rayleigh-Bénard convection modeled by the Boussinesq equations

[TABLE]

Here, the Prandtl number $Pr$ and the Rayleigh number Ra are given by

[TABLE]

where $\nu$ is the kinematic viscosity, $g$ is the acceleration due to gravity, and $\alpha$ is the thermal expansion coefficient. The challenge for Rayleigh-Bénard is to derive upper bounds on the convective heat transport of the form $\text{Nu}-1\leq f(\text{Pr},\text{Ra})$ that hold for all solutions of the Boussinesq system. The connection between the Péclet number Pe and the Rayleigh number Ra is obtained by dotting $\mathbf{u}$ into (9) and averaging over space and time (utilising integration by parts with the given boundary conditions and (10)) to obtain the identity $\text{Pe}^{2}=\text{Ra}\left(\text{Nu}-1\right)$. As a result, any upper bound for the wall-to-wall problem implies an upper bound for convective heat transport among solutions of the Boussinesq equations. For example, an upper bound of the form $\text{Nu}-1\leq c\text{Pe}^{\beta}$ implies an upper bound of the form $\text{Nu}-1\leq c^{\frac{2}{2-\beta}}\text{Ra}^{\frac{\beta}{2-\beta}}$.

The optimal wall-to-wall transport problem is evidently non-convex — there may be many local maxima and global extrema — and only by evaluating the global maximum are we assured of an upper bound for Rayleigh-Bénard. Nevertheless, to make progress, we seek local maxima numerically in this paper and discuss the resulting flows. These flows are of interest in their own right as mechanisms to significantly enhance heat transport, e.g., as targets for control.

The rest of this paper is organized as follows. We first introduce Lagrange multipliers to implement the constraints of the wall-to-wall optimal transport problem and examine the structure of the resulting functional in §2. From insights gained on manipulations of the functional, we develop time-stepping algorithms in §3 to solve the Euler-Lagrange equations. Solutions to the saddle point conditions and resulting transport scalings are presented for time-independent two-dimensional flow fields with no-slip boundaries in §4. Upon investigation of the numerical solutions we see that the fields are to a very high degree separable, i.e., the computed stream functions satisfy $\psi(x,z)\approx\phi(x)\Psi(z)$ and similarly for the other fields. This numerical observation motivates an analytic examination of upper bounds on the wall-to-wall problem with an additional separable ansatz — which is apparently almost satisfied by solutions of the wall-to-wall Euler-Lagrange equations — in §5 leading to conditional upper bounds on heat transport of the form $\text{Nu}\lesssim\text{Pe}^{6/11}$ or, in terms of Rayleigh-Bénard, $\text{Nu}\lesssim\text{Ra}^{3/8}$.

Along the way, we discuss the relationship of the wall-to-wall optimal transport problem to the background method of Doering & Constantin (1996) in §2.4, and to the Howard-Busse-Malkus problem Malkus (1954); Howard (1963); Busse (1969) in §2.6. The perspective developed in those sections inspires the design of a time-stepping algorithm for computing optimal flows, similar to that in Wen et al. (2015) for computing optimal background fields. Our methods of temporal and spatial discretization are described, respectively, in §3.2 and Appendix A. In particular we find in §3.1 that evolving the equations

[TABLE]

forward in pseudo-time $\tau$ (for a fixed constant $\mu$) subject to homogeneous boundary conditions for $\eta$ and $\xi$ and no-slip boundary conditions for $\mathbf{u}$ yields local maxima of the wall-to-wall problem. One of the many benefits of the algorithms described here is that optimal flow fields may be computed for other geometries, e.g., cylinders, given suitable Poisson and Stokes’ equations solvers.

2 Theory

In this section we utilise Lagrange mutlipliers to rewrite the wall-to-wall optimisation problem as one of finding saddle points of a certain unconstrained functional $\mathcal{F}$. We then describe various manipulations that can be performed on $\mathcal{F}$ involving its maximization or minimization or both in §2.1-§2.3. This leads us to a direct comparison between the background method and the wall-to-wall problem in §2.4. In particular we conjecture that there exists a duality gap between the two problems in §2.5 and provide a simple polynomial example to illustrate why one should expect the problems to be distinct. Additionally we note the connection to the Howard-Busse-Malkus problem in §2.6. All of these considerations aid us in the development of numerical algorithms for producing candidate optimizers in §3, and in the proof of our conditional upper bounds in §5.

We begin by introducing a new temperature variable $\theta=T-(1-z)$ and rewrite the advection diffusion equation as

[TABLE]

Next we introduce Lagrange multipliers $\mu$ (a positive real number), $p(x,z,t)$ and $\varphi(x,z,t)$ and the unconstrained functional

[TABLE]

the saddle points of which are candidates for the maximization problem (6). The variables $\varphi$ and $p$ come equipped with their natural boundary conditions. Namely, we impose periodic and homogeneous boundary conditions in the $x$ and $z$ directions respectively for $\varphi$, and the usual implicit boundary conditions for $p$.

Given initial data one could search for time-dependent optimal flow fields of the wall-to-wall problem, but we restrict ourselves to time-independent flow fields for a number of reasons. First, steady fields are far easier to compute numerically and time dependence greatly expands the scope and range of our current endeavor. Second, in the context of simpler models such as the Lorenz equations (with a “heat transport” functional analogous to the Nusselt number), time-dependence was found to never increase transport Souza & Doering (2015a, b). Third, preliminary attempts at computing time-dependent flow fields for the wall-to-wall problem have yielded essentially time-independent results, suggesting that time-dependence may not play a role in significantly enhancing heat transport. More precisely, taking the initial condition of the temperature field to be the conductive state $1-z$ we found that the result of the time-dependent optimization was to move the conductive state into a (locally) optimal steady state, and to hold it there. Finally, for time-independent flows we are guaranteed that maximisers exist (and that the functional $\mathcal{F}$ is differentiable) while as of now there is no such assurance for time-dependent flows.

With these considerations in mind, from this point on we focus on the time-independent functional

[TABLE]

where now the brackets $\langle\cdot\rangle$ are understood to give the spatial average only. We seek the saddle points of (19) and for this it will be useful to consider alternative coordinate systems or “constraint manifolds” that pass through these. Many of the manipulations introduced below extend naturally to time-dependent and/or stress-free flow fields.

The first manipulation we perform is to change variables by $\theta=\xi+\eta$ and $\varphi=\xi-\eta$, following the “symmetrization method” described in Tobasco & Doering (2017) and Doering & Tobasco (2019). Various integrations by parts yield the functional

[TABLE]

The advantage of the new $(\xi,\eta)$ coordinates lies in exposing the underlying geometry of the wall-to-wall problem. Indeed, the functional $\mathcal{S}$ is convex with respect to $\eta$ and concave with respect to $\mathbf{u}$ or $\xi$ when the other is held fixed, whereas the functional $\mathcal{F}$ fails to have such properties. It is a bit like choosing to study the polynomial expression $x^{2}-y^{2}$ instead of $(x+y)(x-y)$.

As a warm up to the manipulations that follow, let us consider this simple example a bit more and remark on how one might search for the saddle points of $s(x,y)=x^{2}-y^{2}$. Gradient ascent/descent procedures are problematic on their own, but can be successfully combined with constraints picking out certain curves. For example, the curve $y=0$ passes through the saddle point $(0,0)$ and the resulting function $s(x,0)=x^{2}$ can be minimized by gradient descent. Thinking procedurally, this particular constraint curve is found by taking the derivative of $s$ with respect to $y$ and setting the result equal to zero, i.e., $\frac{\partial s}{\partial y}=2y=0$.

Returning to the functional $\mathcal{S}$, we proceed in §2.1 and §2.2 to derive various constraint manifolds that pass through its saddle points. We do so by setting the variations of $\mathcal{S}$ with respect to $\eta$, $\xi$, or $\mathbf{u}$ equal to zero and relating these to optimisations of $\mathcal{S}$. Then, in §2.3 and §2.4 we show how the wall-to-wall optimal transport problem and the background method arise from two such optimisation procedures, thereby producing insights into the relationship between the two.

2.1 Variations with respect to $\eta$

We start by taking the variation of $\mathcal{S}$ with respect to $\eta$ and setting it equal to zero. This results in the Euler-Lagrange equations

[TABLE]

for $\eta$. Substituting this back into $\mathcal{S}$ results in the constrained functional

[TABLE]

Stated differently, constraining the variable $\eta$ to satisfy (21) preserves the saddles of (20) and yields

[TABLE]

2.2 Variations with respect to $\xi$ and $\mathbf{u}$

Next, we take the variation of $\mathcal{S}$ with respect to the variable $\xi$ and set it equal to zero. This yields

[TABLE]

and after substituting it back into $\mathcal{S}$ we find

[TABLE]

As before, the same functional is obtained by maximizing $\mathcal{S}$ in $\xi$, i.e.,

[TABLE]

Finally, we take variations over all incompressible flow fields $\mathbf{u}$. This produces the constrained functional

[TABLE]

where $S^{-1}\left(\xi\hat{z}-\xi\nabla\eta\right)$ denotes the unique flow field $\mathbf{u}$ solving

[TABLE]

Note that

[TABLE]

These observations serve as a starting point in §5 for establishing upper bounds on transport and motivate our choice of numerical methods in §3. But first, let us see how the wall-to-wall problem comes out of these manipulations.

2.3 Finding the wall-to-wall problem

Consider the structure of the $S_{\eta}=\min_{\eta}\mathcal{S}$ and $S_{\xi}=\max_{\xi}\mathcal{S}$ functionals for a fixed incompressible velocity field $\mathbf{u}$ with enstrophy $\langle|\nabla\mathbf{u}|^{2}\rangle=\text{Pe}^{2}$. The functional $S_{\eta}$ is concave in $\xi$; likewise $S_{\xi}$ is convex in $\eta$. Thus finding the maximum of the former with respect to $\xi$, or the minimum of the latter with respect to $\eta$, is equivalent to enforcing their Euler-Lagrange equations.

The Euler-Lagrange equation for $S_{\eta}$ in $\xi$ is

[TABLE]

Similarly for $S_{\xi}$ in $\eta$ we have

[TABLE]

These can be written more concisely as the system

[TABLE]

gotten by setting

[TABLE]

into (29) and (30). Hence,

[TABLE]

for each fixed velocity field $\mathbf{u}$.

The formula (35), which first appeared in Tobasco & Doering (2017) and was further analyzed in Doering & Tobasco (2019), is an exact variational characterization of the Nusselt number. It allows the optimal wall-to-wall transport problem to be stated succinctly as

[TABLE]

where $\mathbf{u}$ satisfies the given boundary conditions and intensity constraints. In particular, we note that gradient ascent may be applied with impunity to the constrained functional $\mathcal{S}_{\eta}$ to compute local maxima. This functional can also be used to prove lower bounds on the Nusselt number without having to solve the advection-diffusion equation. Indeed, plugging in any $\xi$ admissible in the previous manipulations into $\mathcal{S}_{\eta}$ yields the lower bound

[TABLE]

In fact, reinterpreting for the moment angle brackets as space and time averages, we note that the bound (37) holds for time-independent flows as well, an observation that was exploited in Tobasco & Doering (2017) and Doering & Tobasco (2019) with “branching” trial functions to prove the scaling $\max\,Nu\sim Pe^{2/3}$ up to logarithmic corrections as $Pe\to\infty$. We return to discuss this asymptotic result in the context of our numerical results much further below. Next, we consider the relationship between the wall-to-wall problem and the background method.

2.4 Finding the background method

Consider now the background method which guarantees the absolute upper bound

[TABLE]

where

[TABLE]

The reader may not immediately recognise this as the familiar background method as it has been applied to Rayleigh-Bénard convection Doering & Constantin (1996). Nevertheless, (38) does follow from applying the usual argument to the wall-to-wall problem. (The resulting bounds carry over to the time-dependent case.)

Let us recall the argument now. Starting with the advection-diffusion equation

[TABLE]

and decomposing $T$ as $T=\xi+1-z+\eta$ yields

[TABLE]

Multiplying through by $\xi$ and integrating by parts yields the balance relation

[TABLE]

Now, utilising

[TABLE]

we subtract twice (42) from (43) to conclude that

[TABLE]

Introducing a Lagrange multiplier $\mu/2$ for the enstrophy constraint $\langle|\nabla\mathbf{u}|^{2}\rangle=Pe^{2}$ yields

[TABLE]

and upon performing the operations $\min_{\eta}\max_{\mathbf{u},\xi}$ we deduce (38).

2.5 A possible duality gap

We can now discuss the relationship between the background method and the wall-to-wall optimal transport problem. Combining the definition of $\mathcal{S}$ from (20) and the identity (45) we see that the background method bound (38) can be alternatively written as

[TABLE]

On the lefthand side appears the wall-to-wall optimal transport problem, while on the righthand side appears the background method. Note this inequality is consistent with the results of §2.3 since in any case

[TABLE]

regardless of the definition of $\mathcal{S}$. Now if $\mathcal{S}$ were convex in $\eta$ and jointly concave in $(\mathbf{u},\xi)$ one would be lead on general grounds via convex duality to conjecture that equality should hold between the lefthand and righthand sides above, in which case the wall-to-wall problem and the background method would turn out to be equivalent. We instead propose that the opposite situation is true and that

[TABLE]

Were such an equality true, it would not preclude the possibility that these quantities achieve the same asymptotic scaling as $\text{Pe}\to\infty$, a situation suggested for 3D wall-to-wall optimal transport by the recent numerical scaling $\max\,\text{Nu}\sim\text{Pe}^{2/3}$ reported for a finite range of Pe in Motoki et al. (2018a). It would also be consistent with the dimension-independent logarithmic lower bound $\max\,\text{Nu}\geq C^{\prime}Pe^{2/3}/(\log\text{Pe})^{4/3}$ proved for all large enough Pe in Tobasco & Doering (2017) and Doering & Tobasco (2019).

Let us illustrate the possibility that (49) holds by considering how the previous manipulations operate on the polynomial

[TABLE]

which we see as analogous to (20). The variable $\tau_{1}$ is analogous to the zero’th Fourier mode of $\eta$, while $\tau_{2}$ and $\tau_{3}$ are to the non-zero Fourier modes.111This polynomial was not derived as a modal truncation of $\mathcal{S}$. That would produce a more complicated example. The variables $u$ and $v$ are analogous to $\xi$ and $\mathbf{u}$. The fact is that

[TABLE]

and that the optimizers for the “background method” $\min$ $\max$ problem are not saddle points of $p$. In particular, the critical point equations

[TABLE]

fail to be satisfied by solutions of the $\min$ $\max$ problem.

To see why this is the case, consider the background method procedure wherein the maximum occurs first. If any of the eigenvalues of (51) are positive then the maximum over $u$ and $v$ yields infinity; thus we must calculate the eigenvalues of (51) to see when this occurs. For fixed $(\tau_{1},\tau_{2},\tau_{3})$ the eigenvalues of the matrix in (51) are $\lambda=2(1-\tau_{1})\pm\sqrt{4(\tau_{2})^{2}+(\tau_{3})^{2}}$. The only way these eigenvalues are nonpositive is if $2(1-\tau_{1})+\sqrt{4(\tau_{2})^{2}+(\tau_{3})^{2}}\leq 0$, or equivalently $2\tau_{1}\geq 2+\sqrt{4(\tau_{2})^{2}+(\tau_{3})^{2}}$. Hence,

[TABLE]

It follows immediately that $\min\max p=1$, and that the minimizer satisfies $\tau_{1}=1$ and $\tau_{2}=\tau_{3}=0$. However, such $\tau_{1}$, $\tau_{2}$ and $\tau_{3}$ cannot be a saddle point of $p$: if $\tau_{2}=\tau_{3}=0$ then from (53) we see that $u^{2}=v^{2}$ and $uv=0$ so that $u=v=0$, but then the $\tau_{1}$ equation cannot be satisfied since $\tau_{1}=1\neq u^{2}+v^{2}$.

Proceeding in the reverse order we find that

[TABLE]

The minimizing $\tau$ satisfy

[TABLE]

The maximum over $u,v$ is given by $u=\pm\sqrt{\frac{2}{5}}$ and $v=\pm\sqrt{\frac{2}{5}}$. Hence, $\max\min p=4/5$.

While in this example the background method procedure (maximum followed by minimum) yields results that are incompatible with the saddle points of (50), the wall-to-wall procedure (minimum followed by maximum) does produce saddle points. Returning to the actual wall-to-wall optimal transport problem, we note that the optimal flow fields reported in §4 exhibit non-trivial non-zero Fourier modes for the variable $\eta$, whereas in the background method optimizers must satisfy $\eta=\eta(z)$. Indeed, if $\eta$ satisfies the spectral stability constraint $\langle\mathcal{Q}\rangle\leq 0$ then so does its periodic average $\overline{\eta}$ in $x$, while by Jensen’s inequality $\langle|\frac{d}{dz}\overline{\eta}|^{2}\rangle\leq\langle|\nabla\eta|^{2}\rangle$ with equality if and only if $\eta=\overline{\eta}$. These observations strongly suggest that the conjectured gap (49) between the wall-to-wall problem and the background method should hold, and in particular that the spectral stability constraint should fail to be satisfied by the true saddle points of $\mathcal{S}$.

2.6 Comparison with the Howard-Busse-Malkus problem

We would be remiss if we did not additionally state the connection of the previous discussions on the wall-to-wall and background method approach with the classic Howard-Busse-Malkus approach put forth in Howard (1963). To see the connection between these, start with $\mathcal{S}$ and restrict attention to incompressible flows $\nabla\cdot\mathbf{u}=0$ with enstrophy $\langle|\nabla\mathbf{u}|^{2}\rangle=\text{Pe}^{2}$ and functions $\eta(z)$. At this point it is useful to introduce notation for the horizontal average of a function,

[TABLE]

Computing the minimum of $\mathcal{S}$ with respect to $\eta(z)$ yields the optimality condition

[TABLE]

whose solution is

[TABLE]

Employing this relation in $\mathcal{S}$ yields

[TABLE]

Making the change of variables $\xi=\alpha\sigma$ for a soon to be determined scalar $\alpha$ transforms this to

[TABLE]

a quadratic function of $\alpha$. Maximizing with respect to $\alpha$ determines the optimal choice

[TABLE]

utilising $\alpha^{*}$ in the above and taking $\mathbf{u}=\text{Pe}\frac{\mathbf{v}}{\sqrt{\langle|\nabla\mathbf{v}|^{2}\rangle}}$ results in the multiplicative form of the functional

[TABLE]

exactly as in Howard (1963) under the appropriate retranscriptions. The functional in (65) is homogeneous with respect to $\sigma$ and $\mathbf{v}$, i.e., $\mathcal{M}[\lambda\sigma,\lambda\mathbf{v},\text{Pe}]=\mathcal{M}[\sigma,\mathbf{v},\text{Pe}]$. It has in the past served as a starting point for the analysis of maximal heat transport, in particular leading in Howard (1963) to bounds on transport under the author’s assumptions of homogeneity and statistical similarity. The connection of this functional to the background method has been pointed out before Kerswell (1998).

However, we point out now that the resulting bounds on transport can be improved beyond those obtained using (65) due to the variational representation (35) of the wall-to-wall problem. After a similar series of manipulations utilising all possible $\eta(x,z)$ instead of functions of $z$ alone, we deduce the improved formula

[TABLE]

where

[TABLE]

The crucial difference between (67) and (65) lies in the fact that

[TABLE]

since $\langle\left(\overline{v_{3}\sigma}-\langle v_{3}\sigma\rangle\right)^{2}\rangle$ is just the zero’th Fourier mode contribution to $\langle\left|\nabla\Delta^{-1}\left(\mathbf{v}\cdot\nabla\sigma\right)\right|^{2}\rangle$. Hence,

[TABLE]

for all $(\sigma,\mathbf{v})$ and the inequality is strict in most cases. Bounds on transport obtained using the functional $\tilde{\mathcal{M}}$ are, therefore, at least as tight as those that have been obtained using Howard’s functional $\mathcal{M}$. We take this as additional evidence of the conjectured duality gap (49) between the wall-to-wall problem and the background method/Howard-Busse-Malkus approach.

3 Gradient flow

The theoretical developments of the previous sections give insight into numerical methods for computing the saddle points of the functionals $\mathcal{F}$ in (19) and $\mathcal{S}$ in (20). In this section we exploit their structure to derive time-stepping methods for solving the Euler-Lagrange equations. An advantage of the approach adopted here is that it is only necessary to have a Poisson or Stoke’s solver to compute candidate optimizers to the wall-to-wall problem.

The Euler-Lagrange equations for the wall-to-wall problem are of the form

[TABLE]

To solve this numerically, we introduce a time-derivative on the lefthand side,

[TABLE]

transforming (70) into a dynamical system where every local maximum is an attracting fixed point. When applied to the functional $\mathcal{S}_{\eta}=\min_{\eta}\mathcal{S}$ from (22) this yields its local maximizers, thereby producing saddle points for $\mathcal{S}$.

Composing (70) with a locally invertible function such that $P(\mathbf{0})=\mathbf{0}$ and then introducing a time derivative defines an alternate system $\dot{\mathbf{x}}=P(\mathbf{f}(\mathbf{x}))$. The danger and boon of choosing such a preconditioner is that the stability of a fixed point may change: we could be computing local maxima, minima, or saddles of our original function $\mathbf{f}$ in (70). For example, Newton’s method may be viewed as choosing the negative inverse Jabobian $J^{-1}(\mathbf{x})$ of $\mathbf{f}$, $P(\mathbf{f}(\mathbf{x}))=-J^{-1}(\mathbf{x})\mathbf{f}(\mathbf{x})$, along with a choice of optimal step size ($\Delta t=1$) upon temporal discretisation. With regards to the functional $\mathcal{F}$ we will take the preconditioner approach implemented in a way that is similar to the algorithm in Wen et al. (2015). With regards to the Euler-Lagrange equations of (18) and (20) we ultimately solve

[TABLE]

or equivalently

[TABLE]

3.1 Gradient ascent methods

Here we outline various methods for solving the Euler-Lagrange equations of the wall-to-wall problem described above. One method is to evolve

[TABLE]

forward in time. This strictly enforces the advection-diffusion equation and its adjoint at every time-step, Utilising $\frac{\delta\mathcal{F}}{\delta\mathbf{u}}$ to compute corrections to the flow field. This approach was taken in Motoki et al. (2018b), and the rate limiting step is the solution of the advection-diffusion equation and its adjoint. In the present work, we take a different approach and compute numerical solutions to (72) and (77) Utilising two different algorithms.

The first of our algorithms involves a time-stepping procedure of the form

[TABLE]

This procedure may be understood as follows. For fixed $\mathbf{u}$, the equations for $\partial_{\tau}\theta$ and $\partial_{\tau}\varphi$ evolve towards the steady state solutions $\frac{\delta\mathcal{F}}{\delta\varphi}=0$ and $\frac{\delta\mathcal{F}}{\delta\theta}=0$, hence these evolutions are a relaxation of fully solving the Euler-Lagrange equations for $\theta$ and $\varphi$. The equation for $\mu$ guarantees that we flow towards a flow field with the desired enstrophy Pe. The last condition (for a fixed $\theta$ and $\varphi$) evolves to a solution of the optimality condition. We enforce incompressibility at every time-step and evolve all fields at once.

There is an alternative description of this algorithm in terms of $\mathcal{S}$. Focusing on the $\theta$ and $\varphi$ equations, we see that evolving

[TABLE]

is equivalent to evolving

[TABLE]

after taking sums and differences (making use of $\theta=\xi+\eta$ and $\varphi=\xi-\eta$) and rescaling time. Thus our time-stepping procedure for $\theta$ and $\varphi$ is equivalent to simultaneously applying gradient ascent for the concave variable $\xi$ and gradient descent for the convex variable $\eta$ in the $\mathcal{S}$ functional. This is similar to the philosophy of Wen et al. (2015). In fact, the only modification required to compute two-dimensional no-slip background fields would be to project the $\eta$ variable to the zero’th Fourier mode at each time-step, by taking the $x$ average of the righthand side of the $\eta$ equation.

The second time-stepping method involves computing the saddles of $\mathcal{S}$ via

[TABLE]

while fixing $\mu$. The resulting enstrophy depends implicitly on $\mu$. Enforcing $\frac{\delta\mathcal{F}}{\delta\eta}=0$ at each time-step yields gradient ascent for the local optima of (22), an unconstrained variational problem. The reason why this algorithm is efficient is due to the many existing algorithms for quick inversion of the Laplacian and the Helmholtz operator.

We found numerically that both (83) and (86) yield the same results. We also implemented various other ascent procedures with different preconditioners. For example, we evolved equations of the form

[TABLE]

forward in time, as well as

[TABLE]

evolving different components on different time-scales. Amongst all of our results, the ones presented in §4 maximize Nu for a given Pe.

3.2 Temporal discretisation

Each of the gradient ascent procedures described in §3 are of the form

[TABLE]

where $x$ is the state vector, $\mathcal{L}$ is an “easily” invertible linear operator, $\mathcal{N}$ is a nonlinear operator, and $f$ is a forcing function. We follow Viswanath & Tobasco (2013) and consider linear multi-step schemes as follows:

[TABLE]

where $s$ is the order of the time-stepping scheme, $a_{i}$ and $b_{i}$ are parameters, and $\Delta t$ is the time-step size. The parameters values for orders $s=1,2$ and $3$ are

[TABLE]

For example with $s=1$ and the advection-diffusion equation,

[TABLE]

we use

[TABLE]

Thus for each time-step we must solve a modified Poisson’s equation of the form

[TABLE]

where $c\geq 0$ and we have made the transcriptions

[TABLE]

Analogous discretisations are used for $\varphi,\eta,$ and $\xi$. For updating the optimality condition with $s=1$ one option is to use

[TABLE]

Each time-step involves solving a modified Stokes equation

[TABLE]

where $c\geq 0$ and we have made the transcriptions

[TABLE]

We solve these boundary value problems by Utilising a pseudo-spectral method, where the wall-bounded direction is represented by Chebyshev polynomials and the periodic directions are represented by Fourier Series; however our use of Chebyshev polynomials utilizes spectral integration in the same way as Viswanath (2015). The Stokes equation (100) was solved Utilising the Kleiser-Schuman algorithm Kleiser & Schumann (1980). See Appendix A for details regarding implementation. A theoretical discussion of the approach is in Viswanath (2015). Let us highlight some of the benefits of spectral integration here:

Memory efficiency: The total memory occupied is linear in the number of points required to describe the state variable $\mathbf{u}$. 2. 2.

Speed: Solving the linear equations involve only tridiagonal matrices which are put into LU form at the beginning of the gradient ascent procedure. 3. 3.

Discretisation Accuracy: Everything is represented using Fourier and Chebyshev modes, allowing for spectral accuracy. 4. 4.

Machine Accuracy: Utilising spectral integration in factored form (see Appendix A) allows one to avoid taking wall-bounded derivatives and only take derivatives in periodic direction, without passing to the nodal domain.

At this point all the pieces are place to compute flow fields maximizing heat transfer. By Utilising gradient ascent or other time-marching schemes of §3.1 it is possible to adapt existing Rayleigh-Bénard codes to find steady maximizing flow fields and compare their thermal transfer to that of natural flows.

4 Two-dimensional computations

For every decade of Péclet we computed approximately 20 logarithmically spaced solutions to the Euler-Lagrange equations (72) and (77). We then performed numerical differentiation of $\log{\text{Nu}}$ and $\log{\Gamma}$ with respect to $\log{\text{Pe}}$ to examine (local) scaling relations. The largest enstrophy satisfied $\text{Pe}\approx 10^{5}$, with an $x$ Fourier, $z$ Chebyshev grid size of $512\times 1025$. The time-stepping code slowed substantially at larger Pe.

In addition to numerical continuation from small Péclet ($\text{Pe}\approx 0$), we started at different points in function space in an attempt to find more optimal solutions. We utilized solutions to other related variational problems as “educated guesses” for new flows at which we started our gradient ascent procedure. We computed solutions to the Euler-Lagrange equations for fixed aspect ratio $\Gamma$ as well as optimizing over all $\Gamma$. The solutions presented here achieved the largest transport.

Figure 1 shows visualizations of the temperature field and the stream function with no-slip boundary conditions. In the low Péclet regime the solutions are convection cells while at higher Pe “recirculation zones” develop on the bottom and top boundaries. Additionally the optimal aspect ratio of a cell size shrinks with ever increasing Pe. For a fixed aspect ratio optimal solutions contain a multiplicity of convection cells. One of the benefits of our numerical approach is the automatic computation of the optimal domain size with little additional overhead.

In Figure 2 we report the Nu-Pe and $\Gamma$-Pe relations and local scaling relations for the best known optimizers (for no-slip boundary conditions). After leaving the “linear” regime where $\text{Nu}\sim\text{Pe}^{2}$ we enter a “fully nonlinear” regime where $\text{Nu}\sim\text{Pe}^{0.54}$. In view of the results of Tobasco & Doering (2017) and Doering & Tobasco (2019) this scaling cannot persist for the globally maximal transport at sufficiently large Pe. Nevertheless, it does appear to be optimal for the computed range of Péclet in two spatial dimensions.222The $\text{Nu}\sim\text{Pe}^{0.54}$ optimal transport scaling in this regime was first reported in Souza (2016) and thereafter confirmed by independent computation by Motoki et al. (2018b) and Ding & Kerswell (2019).

In the linear regime there is little change in the optimal aspect ratio, i.e., the optimal $\Gamma\approx\text{constant}$ while in the fully nonlinear regime a nontrivial $\Gamma\sim\text{Pe}^{-0.37}$ scaling emerges. Interestingly, the local exponent exhibits a somewhat oscillatory relaxation rather than a perhaps more expected monotonic convergence.

4.1 Singular value decomposition

Given the structure of the solutions depicted in Figure 1 one may wonder if the corresponding fields are to first approximation separable. With regards to the stream function $\psi$ defined by $(-\partial_{z}\psi,\partial_{x}\psi)=(u_{1},u_{3})$ this would mean

[TABLE]

We investigate this possibility by minimizing the functional

[TABLE]

subject to appropriate constraints on $\Psi(z)$ and $\phi(x)$, and similarly for other fields such as $\eta$ or $\xi$.333E.g., $\phi(x)$ should be periodic and $\Psi(z)$ should satisfy no-slip boundary conditions. If the minimal value is sufficiently small then we may say that the flow fields are nearly separable.

Upon discretization,

[TABLE]

where $x_{i}=\Gamma i/n$ for $i=0,...,n-1$ and $z_{j}=\frac{1}{2}\left(1+\cos(\pi j/m)\right)$ for $j=0,...,m$ are the collocation points of the Fourier and Chebyshev discretisations, and $w_{ij}$ are weights for approximating the integral by a discrete sum. Spectral accuracy for the approximation (105) was obtained by setting

[TABLE]

where $a_{j}=\pi j/m$ for $j=0,...,m$, see Boyd (2001). Since we are dealing with the average value of the integral there is no factor of $\Gamma$ in the $\Delta x_{i}$.

In the discrete problem corresponding to (105) we minimised the weighted Frobenius norm of the matrix $\psi_{ij}$ with respect to an outer product decomposition, hence for $w_{ij}=1$ the solution to the discrete problem is to take $\Psi$ and $\phi$ to be equal to the largest singular vectors of the matrix $\psi_{ij}$ multiplied by the largest singular value. For $w_{ij}\neq 1$ we rescaled the problem taking $\psi_{ij}^{\prime}=\psi_{ij}\sqrt{w_{ij}}$ and found the largest singular vectors of $\psi_{ij}^{\prime}$. The difference between uniform and nonuniform $w_{ij}$ may be interpreted as minimising with respect to different weighted energy norms and, in our numerical experiments, we found little qualitative difference between utilising one over the other. All singular value decompositions reported here were computed using $w_{ij}=1$.

Figure 3 shows the first three modes in the singular value decomposition of the stream function $\psi$. Figure 4 shows the same for the $\xi$ field and the vertical velocity field $u_{3}$. The modes with smaller singular values may be viewed as a preliminary manifestation of a branching-like pattern, perhaps similar to the ones constructed in Tobasco & Doering (2017) and Doering & Tobasco (2019). The structures corresponding to the largest singular value are the direct analogue for the wall-to-wall problem of the single-wavenumber solutions for the Howard-Busse-Malkus problem produced in Howard (1963).

As it turns out, these dominant structures carry $\approx 99\%$ of the overall heat transport in the computed range of Pe. More precisely, writing

[TABLE]

our solutions achieved $(N_{1}-N_{2})/N_{1}\leq 0.01$ uniformly over all computed Péclet. (This error was even less for smaller Péclet.) The fact that (110) is equal to $\text{Nu}-1$ for solutions to the wall-to-wall problem follows from the Euler-Lagrange equation (77).

Figure 5 depicts the first singular vectors for $\xi$ and $u_{3}$. While in the Howard-Busse-Malkus problem the corresponding $x$ dependence is perfectly sinusoidal, for the wall-to-wall problem the $x$ dependence of the first singular vectors resemble Jacobi elliptic functions. The similarity between $\phi_{1}$ and $\phi_{2}$ motivates the “separable ansatz”

[TABLE]

which we consider further in §5.1 below. There we derive upper bounds extending the single wavenumber analysis from Howard (1963) to the wall-to-wall problem. We will see that the functional from §2.6 can be bounded from above by $\text{Pe}^{6/11}=\text{Pe}^{0.\overline{54}}$ amongst all separable ansatzes in accord with the numerically computed $\text{Nu}\sim\text{Pe}^{0.54}$ scaling.

5 Upper bound

As discussed in the previous section, the flow and temperature fields found to maximize heat transport are nearly separable and achieve the scaling $Nu\sim Pe^{0.54}$ within the range of computed Pe. Furthermore, most of the heat transport scaling was attained via a separable approximation. We now make use of the theoretical developments from §2 to argue for the upper bound $Nu\leq CPe^{0.\overline{54}}$ under the assumption of a perfectly separable ansatz. We do this by rederiving the Howard functional from Howard (1963) (as in §2.6) and assuming a separable ansatz at the outset. It should be emphasized that this is not a rigorous proof for the observed numerical scaling from the previous section since it must still be proven that all solutions in the range $\text{Pe}\lesssim 10^{5}$ are indeed separable, or what is more likely the case, that their non-separable part contributes negligibly to transport.

5.1 Upper bounds within a separable ansatz

To obtain upper bounds on the Nusselt number within a separable ansatz, we start by recalling that the Howard-Busse-Malkus problem bounds the wall-to-wall problem from above. Indeed, as was shown in §2.3 and §2.6,

[TABLE]

since the minimum is only made smaller upon enlarging the class of admissible functions. Thus we are led to consider the functional

[TABLE]

from §2.6, while also restricting ourselves to the separable ansatz

[TABLE]

where ′ denotes an ordinary derivative. Note, while the argument that follows is similar to those of Howard (1963) and Doering & Constantin (1996), it is more general in scope since here we allow $\phi(x)$ to be any periodic function rather than some well-chosen Fourier mode. Our choice to study this more general case is motivated by the previous discussion of numerical results, which showed how in the wall-to-wall problem non-sinusoidal $\phi(x)$ arise. Our goal now is to bound the functional $\mathcal{M}$ from above under the assumptions of separability (117)-(119). Again, we note that this does not provide the unconditional upper bound on $\mathcal{M}$ required to deduce bounds on max Nu according to (115). What we lack is a proof that (117)-(119) are indeed valid assumptions over the given range of Pe. Nevertheless, we proceed.

Note that the velocity field is incompressible, and normalize $\phi$ by requiring

[TABLE]

Given the homogeneity of the functional we may impose that

[TABLE]

Given these, we seek to bound the denominator of (116) from below so as to produce the desired upper bound.

Note the identities

[TABLE]

where for convenience we have abbreviated $\int_{0}^{1}f(z)dz$ as $\int f$. Using these we obtain the lower bound

[TABLE]

by an elementary Young’s inequality. The interpolation inequality

[TABLE]

follows by combining Cauchy-Schwarz, integration by parts, and the normalization (120). Similarly, there holds

[TABLE]

by Cauchy-Schwarz and the net flux constraint.

Proceeding, we deduce that

[TABLE]

by way of Young’s inequality in the form $a+b\geq 3^{1/3}\left(\frac{2}{3}\right)^{2/3}a^{1/3}b^{2/3}$ and (128). Utilising Howard’s lemma (the estimate referred to as such in Doering & Constantin (1996)), we can bound the remaining term by

[TABLE]

Taking $\delta=\left(\int(\Psi^{\prime\prime})^{2}\int(\Xi^{\prime})^{2}\right)^{-1/4}$ we see that the denominator of the functional (116) — in the separable ansatz — is bounded below by

[TABLE]

after minimizing over $\delta$.

Hence, the functional $\mathcal{M}$ is bounded from above by $\text{Pe}^{6/11}=\text{Pe}^{0.\overline{54}}$ in the separable ansatz (117)-(119), in accord with the numerically computed solution scaling of $\text{Pe}^{0.54}$. While this is suggestive of the bound $\text{Nu}\leq C\text{Pe}^{6/11}$, it remains to be shown that the separable ansatz holds up to terms negligible in their transport for the given range of Pe.

6 Summary and Conclusions

We developed and applied time-stepping methods for solving the wall-to-wall problem for steady flows. Along the way we developed theoretical tools illuminating connections between the wall-to-wall problem and both the background method of Doering & Constantin (1996) and Howard’s classic formulation of heat transport bounds for Rayleigh-Bénard convection from Howard (1963). We computed optimal steady two-dimensional flow fields with no-slip boundary conditions whose heat transport scaling is $\text{Nu}\sim\text{Pe}^{0.54}$ for Péclet numbers between $10^{3}$ and $10^{5}$. Upon examining these two-dimensional flows we computed their singular value decomposition, revealing concentration onto one dominant mode. This motivated considering consequences of the additional assumption that the flow fields are separable. Within such an ansatz, we proved a conditional upper bound scaling as $\text{Pe}^{6/11}$ corresponding to the $\text{Ra}^{3/8}$ “single wavenumber” upper bound scaling for Rayleigh-Bénard with no-slip boundary conditions.

In light of Tobasco & Doering (2017) and Doering & Tobasco (2019), however, these scalings for the optimal wall-to-wall problem cannot persist for global optimisers as $\text{Pe}\rightarrow\infty$. There it was proved that wall-to-wall optimal transport in two- and three-dimensions, obeying no-slip or stress-free boundaries conditions, satisfies $\max\,\text{Nu}\geq C^{\prime}\text{Pe}^{2/3}/(\log\text{Pe})^{4/3}$ within logarithms of the a priori upper bound $Nu\leq C\,\text{Pe}^{2/3}$. Curiously, this lower bound coincides to a degree with the numerical results in the range of Péclet numbers explored: logarithmic differentiation yields

[TABLE]

not far from the numerically computed $\text{Nu}\sim\text{Pe}^{0.54}$ result. Although at first we found no direct numerical evidence for the branching patterns constructed in Tobasco & Doering (2017) and Doering & Tobasco (2019), a singular value decomposition revealed branching in its inception. We wonder exactly how large Pe must be for the higher modes to play a significant role in optimizing heat transport in a two-dimensional fluid layer.

The three-dimensional computations reported in Motoki et al. (2018a) show a stark difference in flow structures and transport scaling as compared to the two-dimensional flows for both stress-free boundaries from Hassanzadeh et al. (2014) and no-slip boundaries from §4 of the present paper (see also Motoki et al. (2018b)). Three-dimensional versions of branching appear to achieve the optimal scaling $\text{Nu}\sim\text{Pe}^{2/3}$ already for $\text{Pe}\in[10^{3},10^{4}]$ with a prefactor within $10\%$ of the background upper bound computations from Plasting & Kerswell (2003). We wonder if this $10\%$ difference is a manifestation of our proposed duality gap between the wall-to-wall problem and the background method bounds. In any case, the evidence is that optimal flows take advantage of the presence of a third spatial dimension to maximize thermal transport.

We conclude with a list of five fundamental questions remaining for the optimal wall-to-wall transport problem:

Are steady flows optimal? 2. 2.

Is it possible to prove the a priori upper bound $\text{Nu}\leq CPe^{2/3}/\log(\text{Pe})^{4/3}$ for two-dimensional flows as $\text{Pe}\rightarrow\infty$, and does the upper bound $\text{Nu}\leq CPe^{6/11}$ hold for moderate $Pe$ instead? 3. 3.

Do there exist three-dimensional flow fields achieving $\text{Nu}\sim\text{Pe}^{2/3}$ as $\text{Pe}\rightarrow\infty$? 4. 4.

What do optimally transporting flows look like in other geometries, such as cylinders or domains with holes? 5. 5.

Do structural properties of optimal flows resemble those of buoyancy-driven steady and/or statistically stationary turbulent Rayleigh-Bénard convection?

7 Acknowledgments

We thank Rich Kerswell and Divakar Viswanath for many illuminating conversations. This work was supported in part by NSF Awards DMS-1344199 (ANS), DGE-0813964 and DMS-1812831 (IT), DMS-1515161 and DMS-1813003 (CRD), a Van Loo Postdoctoral Fellowship (IT) and a Guggenheim Foundation Fellowship (CRD).

Appendix A Spatial Discretisation

In what follows we provide details of implementing spectral integration with regards to solving the Poisson and Stokes equation with respect to our boundary conditions and domain.

A.1 Spectral Methods

Taking the Fourier transform of (96) and (100) in the horizontal directions leads to the following set of ODE’s to solve

[TABLE]

and

[TABLE]

for $n,\ell\in\mathbb{Z}$, where $D$ the derivative in the vertical direction, and

[TABLE]

The square root denotes the principle branch. Although we could discretise the spatial coordinates using Chebyshev matrices, we instead use spectral integration. This method of solving boundary value problems of the form

[TABLE]

subject to boundary conditions has numerous advantages over the differentiation matrix approach as in Trefethen (2001). With spectral integration the operators that must be inverted have bounded condition numbers and are banded matrices as opposed to dense matrices with unbounded condition numbers.

Instead of solving for functions on $z\in[0,1]$ it is more convenient to take $z\in[-1,1]$ and then convert results back to the original domain. Solutions in the $z\in[0,1]$ domain–denoted by subscripted $1$’s as in $\theta_{1},\mathbf{u}_{1}$–are related to solutions in the $z\in[-1,1]$ domain–denoted by subscripted $2$’s as in $\theta_{2},\mathbf{u}_{2}$–via the following relations

[TABLE]

When performing calculations we use the $[-1,1]$ domain but report results in terms of the original $z\in[0,1]$ domain.

Computing averages (such as the Nusselt number) can be achieved with spectral accuracy. As mentioned in Boyd (2001) the trapezoidal rule is spectrally accurate for periodic functions and for bounded domains there are quadrature weight formulas both in terms of the Chebyshev nodal and modal values. Furthermore all products are computed by taking the inverse transforms, multiplying in real space, and converting back to spectral space, as is common for pseudo-spectral methods.

A.2 Spectral Integration

We have seen that gradient ascent for the wall-to-wall problem reduces to solving differential equations of the form

[TABLE]

where $k\in\mathbf{R}$ and $z\in[-1,1]$. The differential equation has the solution

[TABLE]

where $C$ enforces boundary or integral constraints. This reduces the problem to quadrature and indeed the method that we adopt implicitly constructs the solution in this manner as noted in Viswanath (2015).

To solve (146) we use a modern form of spectral integration developed by Viswanath Viswanath (2015). The general principle is remarkably simple. First compute the homogeneous solution

[TABLE]

and then the particular solution

[TABLE]

so that the general solution is then a linear combination of the particular and homogeneous solution

[TABLE]

where $C$ is a constant that enforces boundary conditions or integral constraints. This basic decomposition of the general solution of an ordinary differential equation serves as the primary building block in the construction of solutions to (96) and (100).

We now discuss how to construct $y^{h}$ and $y^{p}$ in the domain $z\in[-1,1]$. First write $y$ as Chebyshev expansion of the form

[TABLE]

where $P_{n}(x)$ $n=0,1,2,...$ are the Chebyshev polynomials defined by

[TABLE]

The factor of $1/2$ in front of the $y_{0}$ term is standard and convenient for formulas later on.

Let $\mathcal{T}_{n}(y)$ denote the $n$’th Chebyshev coefficent of $y$, e.g.

[TABLE]

where $\int y$ here denotes a particular anti-derivative of $y$. The coefficients may be computed by the linear operator

[TABLE]

Instead of working with (146) directly we work with the equation in integral form

[TABLE]

where $C$ is a constant of integration. An important observation is that if $f=Dg$ for some function $g$, then we do not need to differentiate $g$ to write it in the form (156) but rather we can directly use

[TABLE]

In other words we can avoid numerically differentiating $g$.

Use the relation (154) to construct (156) as a system of equations for the Chebyshev coefficents, i.e.

[TABLE]

The $f_{n}$ are the Chebyshev coefficents of the forcing function $f$. Note that any choice of $C$ yields a particular solution to the problem, but does not enforce the proper boundary conditions. This is an infinite dimensional tridiagonal system of equations for the Chebyshev coefficients of $y$.

In order to be amenable to computation any such system of equations must be truncated. Thus we assume that the solution $y$ and forcing function $f$ is well represented by a finite truncation

[TABLE]

For example, with $N=6$ we would have the following system of equations to solve,

[TABLE]

The use of this finite representation has additional benefits. The Chebyshev spectral coefficients of $y$ are related to the nodal values of $y$ evaluated at $\cos(\pi j/(N-1))$ for $j=0,1,...,N-1$ via the fast-cosine transform. This allows us to quickly convert spectral coefficients into real space, and herein lies one of the many advantages of using a Chebyshev series.

We are now in a position to show how to construct numerical homogeneous and particular solutions $y^{h}$ and $y^{p}$. To construct $y^{h}$ we solve the following problem

[TABLE]

subject to $\mathcal{T}_{0}(v)=0$. The general solution to (164) is

[TABLE]

where $v^{p}=-1/2$ is the particular solution. The condition $\mathcal{T}_{0}(v)=0$ guarantees that $C\neq 0$ since

[TABLE]

and $\mathcal{T}_{0}(v^{h})\neq 0$. Thus the homogeneous solution is $y^{h}=v+1/2$. Denoting the Chebyshev series of $v$ by $v_{n}$ for $n=0,...,N-1$, we find $v$ by solving the system of equations

[TABLE]

We set the $N-1$’st Chebyshev coefficent of $v$ and $f$ as zero for convenience. For example, for $N=6$ we solve the system of equations,

[TABLE]

for the Chebyshev coefficients of $v$ and then add $1$ to $v_{0}$ to construct $y^{h}$.

For the construction of the particular solution $y^{p}$ we solve (149) subject to $\mathcal{T}_{0}(y^{p})=0$. Thus we solve the system of equations

[TABLE]

For $N=6$ implies solving the tridiagonal system

[TABLE]

Now that we have constructed the homogeneous and particular solutions we can enforce boundary conditions. Given $y(a)=b$ for $a\in[-1,1]$ the value of the constant $C$ is determined:

[TABLE]

Typically we enforce the boundary conditions at the endpoint $z=\pm 1$ for which we have readily available formulas to compute the values of $y^{h}$ and $y^{p}$ in terms of their Chebyshev coefficents. For $f=f_{0}/2+\sum_{n=1}^{N-2}f_{n}P_{n}+P_{N-1}f_{N-1}/2$ we compute

[TABLE]

We now show how to construct solutions to the second order equation

[TABLE]

For this problem we have two homogeneous solutions and one particular solution so that the general solution to this differential equation is of the form

[TABLE]

We solve this as a system of two different equations

[TABLE]

First we find the particular and homogeneous solution to

[TABLE]

as described previously so that we have

[TABLE]

Then we find $y^{h_{2}}$ and $y^{p}$ by solving

[TABLE]

subject to $\mathcal{T}_{0}(y^{h_{2}})=0$ and $\mathcal{T}_{0}(y^{p})=0$ via spectral integration for the particular solution. Note that $y^{h_{2}}$ constructed in this manner is linearly independent from $y^{h_{1}}$. To enforce boundary conditions we must now invert a matrix for the coefficients $C_{1}$ and $C_{2}$. For example for boundary condition enforced at the endpoints $z=\pm$ we have

[TABLE]

for the coefficients $a$ and $b$.

Once we have our solution $y$ we find derivatives without recourse to differentiation matrices or Fourier transform methods. Indeed, in the first-order case, i.e. for $y$ that satisfy

[TABLE]

we find the derivative by simply rearranging the equation

[TABLE]

Hence differentiating is obtained by merely summing the the forcing function and the solution multiplied by $k$. For the second order case, i.e. for $y$ that satisfy

[TABLE]

we compute

[TABLE]

It seems that computing derivatives in this manner lead to an order of magnitude improvement of the relative error over other methods (either differentiation matrices or Fourier methods). Using these formulas we compute the derivatives of $y$ at the endpoints in an alternative manner. Instead of using (178) and (179) to the solution $y$, we apply (176) and (176) to (199). From numerical experimentation it does not seem to matter which way the derivatives were evaluated at the endpoint.

We now have all the pieces to solve the modified Poisson’s equation (96). As was the case with the Chebyshev series we must truncate the number of horizontal Fourier wave-modes in hopes that we achieve a good representation of our solution. This problem is linear thus we solve

[TABLE]

subject to $\theta_{n\ell}(z=\pm 1)=0$, mode by mode. Here

[TABLE]

and $N,L$ are even. Thus we solve $N\times L$ second order boundary value problems. By using symmetry or realness of the variables we reduce computation. The $n=\ell=0$ mode for $c=0$ must be handled separately depending on boundary conditions. For example, if $\partial_{z}\theta_{00}(z=\pm 1)=0$ there is the solvability requirement that $0=\int_{z=-1}^{1}f_{00}(z)dz$ otherwise no solution exists. Changing the value of $\theta_{00}$ by any constant still produces a solution, thus instead we replace one of the boundary conditions with a normalization condition such as $\int_{-1}^{1}\theta_{00}dz=0$.

Solving the modified Stokes equation (100) is more complicated and is the subject of the next section.

A.3 Kleiser-Schumann Algorithm

The Kleiser-Schumann algorithm is a method for solving the modified Stokes problem (100), see Kleiser & Schumann (1980); Viswanath & Tobasco (2013). Since the modified Stokes problem is linear we solve wave-number by wave-number equations of the form

[TABLE]

where

[TABLE]

Here $\imath=\sqrt{-1}$ where the square-root is interpreted as the principle branch. From now on we drop the subscript $n\ell$ with the understanding that the modified Stokes problem must be solved for all values of $n$ and $\ell$ individually. For now we consider the $k\neq 0$ case and discuss how to handle this mode separately at the end.

The single wave-number problem is to solve

[TABLE]

Taking the divergence of the first equation yields the following equation for $p$

[TABLE]

As we have seen before we may write the general solution to this problem as the sum of two homogeneous terms and a particular solution

[TABLE]

Hence we may write the equation for $\mathbf{u}$ as

[TABLE]

We split the problem of finding solutions to $\mathbf{u}$ into three parts. Let $\mathbf{u}\hskip 2.84526pt^{i}$ for $i=1,2,3$ be solutions to

[TABLE]

where each $\mathbf{u}\hskip 2.84526pt^{i}$ for $i=1,2,3$ satisfies the boundary conditions for $\mathbf{u}$. With this we write $\mathbf{u}$ as

[TABLE]

To find $C_{1}$ and $C_{2}$ we use auxiliary conditions derived by enforcing incompressibility on the boundary.

For no-slip boundary conditions we have $\partial_{z}w(z=\pm 1)=0$ and for stress-free boundary conditions $\partial_{zz}w(z=\pm 1)=0$. That is to say (225) applies to each component hence to the vertical velocity $w$

[TABLE]

For example with no-slip boundary conditions we apply the vertical derivative $D$ to both sides and solve the following system of equations for $C_{1}$ and $C_{2}$

[TABLE]

For fixed time-step sizes the matrix for are precomputed and factorized only once. Furthermore $\mathbf{u}\hskip 2.84526pt^{1}$ and $\mathbf{u}\hskip 2.84526pt^{2}$ is precomputed. Hence at each time step we only need to solve for $\mathbf{u}\hskip 2.84526pt^{3}$ and the coefficients $C_{1}$ and $C_{2}$.

For this problem the $k=0$ case must be handled separately as well. This is due to the fact that the homogeneous solution for the pressure is of the form $p^{h_{1}}=1/2$ and $p^{h_{2}}=z/2$. Letting $\mathbf{u}\hskip 2.84526pt^{1}=(u^{1},v^{1},w^{1})$ the $k=0$ case implies that $\mathbf{u}\hskip 2.84526pt^{1}=0$ for no-slip boundary conditions and stress-free boundary conditions. For $\beta=0$ one has $(u^{1},v^{1},w^{1})=(D_{1},D_{2},0)$ in the stress-free case where $D_{1}$ and $D_{2}$ are arbitrary constants. By specifying that $(u^{1},v^{1})$ are mean zero we may set these arbitrary constants to zero. The constant $C_{1}$ becomes a free parameter that doesn’t affect the physical flow field $\mathbf{u}$. This corresponds to the fact that the pressure may be changed by an arbitrary constant. We may choose this value such that the average of $p$ is zero, but this is by no means necessary. Now we must find $C_{2}$. From incompressibility and from the boundary conditions we see that $Dw=0\Rightarrow w=0$ which implies444Note that for $\beta=0$ we have $w^{2}=(z^{2}-1)/4$ and for $\beta\neq 0$ we have $w^{2}=-\frac{1}{2\beta^{2}}+\frac{\sinh(\beta)\cosh(\beta z)}{\beta^{2}\sinh(2\beta)}$. $C_{2}w^{2}=w^{3}$. Taking the derivative we find that $Dw^{3}(z\pm 1)/Dw^{2}(z\pm 1)=C_{2}$. This may appear overconstrained, but since we are guaranteed that that the functions $w^{2}$ and $w^{3}$ are proportional to one another we could take either conditions to evaluate the constant $C_{2}$.

Appendix B Optimality Condition for Domain Size

Prior work (Hassanzadeh et al., 2014) indicates that there seems to be an optimal domain size in the periodic direction $x\in[0,\Gamma]$. The derivative of the functional $\mathcal{F}$ with respect to domain size $\Gamma$, is easiest to compute by rewriting $\mathcal{F}=\langle\mathcal{L}\rangle$ with

[TABLE]

Differentiating this functional with respect to $\Gamma$ evaluated on a solution to the Euler-Lagrange equations yields

[TABLE]

The latter part is the spatial Hamiltonian density of $\mathcal{L}$. It may seem that $\frac{\delta\mathcal{F}}{\delta\Gamma}$ is a function of $x$, but the Lagrangian $\mathcal{L}$ has no explicit spatial dependence hence the quantity $\int_{0}^{1}\mathcal{H}^{x}dz$ is independent of the horizontal variable $x$. In other words $\int_{0}^{1}\mathcal{H}^{x}dz=\langle\mathcal{H}^{x}\rangle$. This allows us to simplify the formula and compute555If this calculation is unfamiliar to the reader we refer to Souza (2016) but the essence of the calculation comes from the first exercise of Feynman and Hibbs Feynman & Hibbs (1965).

[TABLE]

The optimal domain size condition for $\mathcal{S}$ is similar but not exactly the same due to a redefinition of the pressure term.

We include an aspect ratio flow of the form $\partial_{\tau}\Gamma=\frac{\delta\mathcal{F}}{\delta\Gamma}$ in addition to the previous time-stepping procedures. There are subtleties associated with implementing this flow to make the procedure successful. Empirically it is found to be useful to withhold evolution of the aspect ratio until one was “close” to a solution of the Euler-Lagrange equations. Furthermore taking a time-step in tandem with the other evolution seemed prohibitively slow, thus the evolution is only taken every $n$’th time-step with respect to the other evolutions, where $n=5$ works well for our purposes. The overhead of changing the domain slows down computations by a negligible amount once these modifications are implemented.