The time evolution equation for advective heat transport as a constraint for optimal bounds in Rayleigh-B\'enard convection

TL;DR

This paper introduces a new constraint based on the time evolution of advective heat transport to improve bounds on heat transfer in Rayleigh-Bénard convection, refining previous theoretical limits.

Contribution

It presents a novel constraint from advective heat transport evolution, enhancing the accuracy of bounds on heat transfer in convection models.

Findings

Improved bounds on heat transport in Rayleigh-Bénard convection.

Enhanced understanding of the role of advective heat transport constraints.

Refinement of theoretical limits for convective heat transfer.

Abstract

Upper bounds on the heat transport and other quantities of interest in Rayleigh-B\'enard convection are derived in previous work from constraints resulting from the equations of time evolution for kinetic energy, the root mean square of temperature, and the temperature averaged over horizontal planes. Here, we investigate the effect of a new constraint derived from the time evolution equation for the advective heat transport. This additional constraint leads to improved bounds on the toroidal dissipation.

Equations

$$\frac{1}{\mathrm{Pr}}\left({\partial }_t\mathbf{v}+\mathbf{v}\cdot \nabla \mathbf{v}\right)$$

$${\partial }_t\theta +\mathbf{v}\cdot \nabla \theta -v_z$$

$$\nabla \cdot \mathbf{v}$$

$$\mathrm{Pr}=\frac{\nu }{\kappa },\mathrm{Ra}=\frac{g\alpha \Delta T{h}^3}{\kappa \nu }$$

$$\frac{1}{\mathrm{Pr}}\left({\partial }_t{\nabla }^2{\Delta }_2\varphi +\hat{\mathbf{z}}\cdot \nabla \times \nabla \times \left[(\nabla \times \mathbf{v})\times \mathbf{v}\right]\right)$$

$$\frac{1}{\mathrm{Pr}}\left({\partial }_t{\Delta }_2\psi -\hat{\mathbf{z}}\cdot \nabla \times \left[(\nabla \times \mathbf{v})\times \mathbf{v}\right]\right)$$

$${\partial }_t⟨\theta {⟩}_A=⟨{\partial }_z(\theta {\Delta }_2\varphi ){⟩}_A+⟨{\partial }_z^2\theta {⟩}_A.$$

$${\partial }_t⟨\frac{1}{2}{\theta }^2⟩=-⟨\theta {\Delta }_2\varphi ⟩-⟨|\nabla \theta {|}^2⟩.$$

$${\partial }_t⟨\frac{1}{2}{\mathbf{v}}^2⟩=-\mathrm{Pr}\mathrm{Ra}⟨\theta {\Delta }_2\varphi ⟩-\mathrm{Pr}⟨|(\hat{\mathbf{z}}\times \nabla ){\nabla }^2\varphi {|}^2+|\nabla {\partial }_x\psi {|}^2+|\nabla {\partial }_y\psi {|}^2⟩.$$

$${\partial }_t⟨v_z\theta ⟩=⟨v_z^2⟩+\mathrm{Pr}\mathrm{Ra}⟨{\theta }^2⟩-(1+\mathrm{Pr})⟨\nabla \theta \cdot \nabla v_z⟩-⟨\theta {\partial }_{z}p⟩.$$

$${\nabla }^2{p}_1$$

$${\nabla }^2{p}_2$$

$$\frac{1}{\mathrm{Pr}\mathrm{Ra}}\int \theta {\partial }_z{p}_{2}dV$$

$$⟨\theta {\partial }_z{p}_2⟩=\mathrm{Pr}\mathrm{Ra}⟨|{\partial }_z\nabla {\Delta }^{-1}\theta {|}^2⟩$$

$${\partial }_t⟨{\phi }_n\theta ⟩=⟨{\phi }_n{\partial }_z(\theta {\Delta }_2\varphi )⟩+⟨{\phi }_n{\nabla }^2\theta ⟩.$$

$$F({\lambda }_1...{\lambda }_N,{\lambda }_R,{\lambda }_E,{\lambda }_M,\theta ,\varphi ,\psi )=\sum _{n=1}^N{\lambda }_n{\partial }_t⟨{\phi }_n\theta ⟩-{\lambda }_R{\partial }_t⟨\frac{1}{2}{\theta }^2⟩-\frac{{\lambda }_E}{\mathrm{Pr}}{\partial }_t⟨\frac{1}{2}{\mathbf{v}}^2⟩+{\lambda }_M{\partial }_t⟨v_z\theta ⟩$$

$$-Z+{\lambda }_0+F({\lambda }_1...{\lambda }_N,{\lambda }_R,{\lambda }_E,{\lambda }_M,\theta ,\varphi ,\psi )\ge 0$$

$${\lambda }_0$$

$$-Z+{\lambda }_0+\sum _{n=1}^N{\lambda }_n\left[⟨{\phi }_n{\partial }_z(\theta {\Delta }_2\varphi )⟩+⟨{\phi }_n{\nabla }^2\theta ⟩\right]+{\lambda }_R\left[⟨\theta {\Delta }_2\varphi ⟩+⟨|\nabla \theta {|}^2⟩\right]$$

$$+{\lambda }_E\left[\mathrm{Ra}⟨\theta {\Delta }_2\varphi ⟩+⟨|(\hat{\mathbf{z}}\times \nabla ){\nabla }^2\varphi {|}^2⟩+d^2\right]$$

$$+{\lambda }_M\left[⟨v_z^2⟩+\mathrm{Pr}\mathrm{Ra}\left(⟨{\theta }^2⟩-⟨|{\partial }_z\nabla {\Delta }^{-1}\theta {|}^2⟩\right)-(1+\mathrm{Pr})⟨\nabla \theta \cdot \nabla v_z⟩-⟨\theta {\partial }_z{p}_1⟩\right]\ge 0$$

$$\int (z-\frac{1}{2}){\partial }_z{p}_{1}dV=\int dz(z-\frac{1}{2}){\partial }_z\int \int dxdy{p}_1$$

$$\int \int dxdy{\nabla }^2{p}_1$$

$$\int (z-\frac{1}{2}){\partial }_z{p}_{1}dV=-\int dz(z-\frac{1}{2}){\partial }_z\int \int dxdy{v}_z^2=\int v_z^{2}dV$$

$${\lambda }_M\int (T-\frac{1}{2}){\partial }_z{p}_{1}dV\le |{\lambda }_M|\int |T-\frac{1}{2}|\cdot |{\partial }_z{p}_1|dV\le \frac{1}{2}|{\lambda }_M|\int |{\partial }_z{p}_1|dV$$

$$⟨|{\partial }_z{p}_1|⟩\le \frac{1}{2}⟨\sum _{i,j}({\partial }_i{v}_j{)}^2⟩.$$

$${\lambda }_0$$

$$-Z+{\lambda }_0+\sum _{n=1}^N{\lambda }_n\left[⟨{\phi }_n{\partial }_z(\theta {\Delta }_2\varphi )⟩+⟨{\phi }_n{\nabla }^2\theta ⟩\right]+{\lambda }_R\left[⟨\theta {\Delta }_2\varphi ⟩+⟨|\nabla \theta {|}^2⟩\right]$$

$$+{\lambda }_E\left[\mathrm{Ra}⟨\theta {\Delta }_2\varphi ⟩+⟨|(\hat{\mathbf{z}}\times \nabla ){\nabla }^2\varphi {|}^2⟩+d^2\right]$$

$$+{\lambda }_M\left[\mathrm{Pr}\mathrm{Ra}\left(⟨{\theta }^2⟩-⟨|{\partial }_z\nabla {\Delta }^{-1}\theta {|}^2⟩\right)-(1+\mathrm{Pr})⟨\nabla \theta \cdot \nabla v_z⟩\right]\ge \frac{1}{4}{\lambda }_{\mathrm{abs}}⟨\sum _{i,j}({\partial }_i{v}_j{)}^2⟩$$

$${\lambda }_{\mathrm{abs}}\ge {\lambda }_M\ge -{\lambda }_{\mathrm{abs}}$$

$$⟨\sum _{i,j}({\partial }_i{v}_j{)}^2⟩=⟨|(\hat{\mathbf{z}}\times \nabla ){\nabla }^2\varphi {|}^2⟩+d^2$$

$${\lambda }_0$$

$$-Z+{\lambda }_0+\sum _{n=1}^N{\lambda }_n\left[⟨{\phi }_n{\partial }_z(\theta {\Delta }_2\varphi )⟩+⟨{\phi }_n{\nabla }^2\theta ⟩\right]+{\lambda }_R\left[⟨\theta {\Delta }_2\varphi ⟩+⟨|\nabla \theta {|}^2⟩\right]$$

$$+{\lambda }_E\left[\mathrm{Ra}⟨\theta {\Delta }_2\varphi ⟩+⟨|(\hat{\mathbf{z}}\times \nabla ){\nabla }^2\varphi {|}^2⟩+d^2\right]$$

$$+{\lambda }_M^{\prime }\left[\mathrm{Ra}\left(⟨{\theta }^2⟩-⟨|{\partial }_z\nabla {\Delta }^{-1}\theta {|}^2⟩\right)-⟨\nabla \theta \cdot \nabla v_z⟩\right]\ge 0$$

$$\theta =\sum _{n=1}^N\sum _{k_x}\sum _{k_y}{\hat{\theta }}_{n,k_x,k_y}{T}_n(2z-1)e^{i(k_{x}x+k_{y}y)}.$$

Full text

††preprint: APS/123-QED

I Introduction

Thermal convection is among the best studied problems of fluid mechanics. Its most basic idealization is Rayleigh-Bénard convection in which a fluid fills a plane layer infinitely extended in the horizontal directions, heated from below and cooled from above. The applied temperature difference, or the amplitude of the driving of the flow, is commonly expressed in terms of the Rayleigh number. Numerical simulations and experiments can investigate convection only up to a certain Rayleigh number. It is also possible to derive rigorous bounds on, for example, the heat transport across the convecting layer (Howard, 1963; Busse, 1969; Doering and Constantin, 1996; Plasting and Ierley, 2005; Seis, 2015). The bounds derived in these references are valid at all Rayleigh numbers, but they tend to overestimate the actual heat transport. If power laws are fitted to both the bounds and numerically computed time averages, the fit to the bounds have a larger prefactor and usually a larger exponent. Furthermore, the bounds do not show any Prandtl number dependence.

These deficiencies arise because the derivation of the bounds does not exploit the full equations of evolution, but only a few integrals deduced from them: the energy budget, the relation between advective heat transport and dissipation, and the temperature equation integrated over horizontal planes. Seis Seis (2015) also makes use of the maximum principle for temperature. An improvement of the bounds necessarily requires to take advantage of additional constraints. This process is already well understood for systems of ordinary differential equations such as the Lorenz model Tobasco, Goluskin, and Doering (2018). In the context of convection, the implementation of further constraints is more cumbersome. Vitanov and Busse Vitanov and Busse (2000) split the energy budget into two equations, one for the poloidal and one for the toroidal energy. The resulting optimization problem preserves a dependence on the Prandtl number (and also on rotation if one is interested in rotating convection). But there is a price to pay in this approach. The optimization problem is not convex, and its Euler-Lagrange equations have to be solved numerically. Because of extensive coupling between different modes in a spectral decomposition of these equations, a numerical solution can only be obtained at low resolution, and hence at low Rayleigh numbers.

The approach pursued in this paper is to solve a semidefinite program (SDP). Previously obtained bounds on Rayleigh-Bénard convection can be reproduced by this method Tilgner (2017a). It has already been demonstrated how one can obtain increasingly sharp bounds on solutions of systems of ordinary differential equations by including more and more constraints Tobasco, Goluskin, and Doering (2018); Goluskin (2018). While the same systematic procedure is in principle possible for the Navier-Stokes equation Chernyshenko et al. (2014), it becomes unpractical for Rayleigh-Bénard convection even at modest Rayleigh numbers because all energy unstable modes need to be retained which leads to a large SDP. Fantuzzi et al. Fantuzzi, Pershin, and Wynn (2018) study bounds for Bénard-Marangoni convection and in this context discuss in depth limitations of the SDP approach and possible future lines of investigation. At present, the most straightforward approach remains to formulate an optimization problem in the form of an SDP with many decoupled linear matrix inequalities rather than a problem with a few, but large linear matrix inequalities. The goal of the present paper is to derive a Prandtl number dependent optimization problem for Rayleigh-Bénard convection which is convex and whose solution does not imply a significantly larger computational burden than optimization problems solved previously to find bounds on various observables of interest in Rayleigh-Bénard convection.

II The optimization problem

Let us consider the problem of Rayleigh-Bénard convection within the Boussinesq approximation for stress free boundaries. This choice of boundary conditions is not essential for obtaining the new constraint but it simplifies the calculation. A Cartesian coordinate system $(x,y,z)$ is chosen such that a plane layer is infinitely extended in the $x-$ and $y-$ directions and the boundaries of the layer are separated by the distance $h$ along the $z-$direction. Gravitational acceleration $g$ is acting along the negative $z-$direction. The layer is filled with fluid of density $\rho$, kinematic viscosity $\nu$, thermal diffusivity $\kappa$, and thermal expansion coefficient $\alpha$. Top and bottom boundaries are held at the fixed temperatures $T_{\mathrm{top}}$ and $T_{\mathrm{top}}+\Delta T$, respectively. We will consider the equations of evolution immediately in nondimensional form, choosing for units of length, time, and temperature deviation from $T_{\mathrm{top}}$ the quantities $h$, $h^{2}/\kappa$ and $\Delta T$. With this choice, the equations within the Boussinesq approximation for the fields of velocity $\bm{v}(\bm{r},t)$, temperature $T(\bm{r},t)$ and pressure $p(\bm{r},t)$ become:

[TABLE]

In these equations, $T=\theta+1-z$, so that $\theta$ represents the deviation from the conduction profile. The Prandtl number $\mathrm{Pr}$ and the Rayleigh number $\mathrm{Ra}$ are given by

[TABLE]

and $\bm{\hat{z}}$ denotes the unit vector in $z-$direction. The boundary conditions on temperature require that $\theta=0$ at $z=0$ and 1 and the stress free conditions chosen here lead to $\partial_{z}v_{x}=\partial_{z}v_{y}=v_{z}=0$ at the boundaries.

It will be helpful to reduce the number of dependent variables by introducing poloidal and toroidal scalars $\phi$ and $\psi$ such that $\bm{v}=\nabla\times\nabla\times(\phi\bm{\hat{z}})+\nabla\times(\psi\bm{\hat{z}})$ and $\nabla\cdot\bm{v}=0$ is satisfied by construction. The $z-$component of the curl and the $z-$component of the curl of the curl of eq. (1) yield the equations of evolution for $\phi$ and $\psi$,

[TABLE]

with $\Delta_{2}=\partial^{2}_{x}+\partial^{2}_{y}$. For brevity, $\bm{v}$ is not replaced by its expression in terms of $\phi$ and $\psi$ in these equations. The stress free boundary conditions translate into $\phi=\partial^{2}_{z}\phi=\partial_{z}\psi=0$.

The currently known methods for obtaining upper bounds on fluid dynamic quantities cannot exploit the equations of evolution in full detail but extract averages from the original equations. Several types of averages will become important: the average over the entire volume, denoted by angular brackets without subscript, the average over an arbitrary plane $z=\mathrm{const.}$, denoted by $\langle...\rangle_{A}$, and the average over time, which will be signaled by an overline.

Three useful relations between averages can now directly be obtained. The average of eq. (2) over planes $z=const.$ leads to

[TABLE]

Multiplication of eq. (2) with $\theta$ and a subsequent volume average yields

[TABLE]

Finally, the dot product of $\bm{v}$ with eq. (1), followed by a volume average, leads to

[TABLE]

These three equations form the basis of the previous work on bounds on, for instance, the Nusselt number $\mathrm{Nu}$ given by $\mathrm{Nu}=1-\langle\overline{\theta\Delta_{2}\phi}\rangle=1+\langle\overline{v_{z}\theta}\rangle$.

An additional relation will be derived below by computing $\partial_{t}\langle v_{z}\theta\rangle$ from the time evolution equation of the advective heat transport $v_{z}\theta$. There are at least two reasons why such a relation looks promising. Bounds have been computed for double diffusive convection. In this problem, salinity $S$ drives convection together with temperature, and salinity obeys the same advection-diffusion equation as temperature except for a different diffusion constant. This problem looks superficially identical to ordinary convection and one may think that adding the equation for $\partial_{t}\langle\frac{1}{2}S^{2}\rangle$ to the equations (7-9) provides us with enough information to derive bounds on double diffusive convection. In fact, it does not. It is necessary to include the equation for $\partial_{t}\langle\theta S\rangle$ to derive bounds Balmforth et al. (2006). This demonstrates the usefulness of considering cross products of different physical quantities and motivates us to also look at $\partial_{t}\langle v_{z}\theta\rangle$. Another type of flows to which bounding methods were applied in the past are flows in periodic volumes driven by a body force $\bm{f}$. Doering and Foias (2002); Childress, Kerswell, and Gilbert (2001); Rollin, Dubief, and Doering (2011); Tilgner (2017b) In these problems, it is essential to include the equation for $\partial_{t}\langle\bm{f}\cdot\bm{v}\rangle$ obtained by forming the dot product of the momentum equation with $\bm{f}$. In Rayleigh-Bénard convection, the buoyancy force plays the role of the body force in the momentum equation. From this analogy, we have another reason to look at $\partial_{t}\langle v_{z}\theta\rangle$.

If we multiply the $z-$component of eq. (1) by $\theta$, eq. (2) by $v_{z}$, add the two equations and average the sum over the volume, we obtain the time evolution equation for the advective heat transport:

[TABLE]

Thanks to the stress free boundary conditions and the boundary condition on temperature, one deduces from eq. (1) that $\partial_{z}p=0$ at the boundaries. The divergence of eq. (1) yields $\nabla^{2}p=-\nabla\cdot[(\bm{v}\cdot\nabla)\bm{v}]+\mathrm{Pr}\mathrm{Ra}\partial_{z}\theta$. It is therefore possible to split the pressure $p$ into two terms $p_{1}$ and $p_{2}$ such that $p=p_{1}+p_{2}$ and

[TABLE]

with the boundary conditions that $\partial_{z}p_{1}=\partial_{z}p_{2}=0$ at $z=0$ and 1.

It will now be shown that the expression $\langle\theta\partial_{z}p_{2}\rangle$ is quadratic in $\theta$ and positive. To this end, we first consider a layer with periodic boundary conditions and finite periodicity length in the lateral directions. At the end of the calculation, we will send the periodicity length to infinity in order to obtain the desired result for the infinitely extended layer. It helps to introduce two functions $f$ and $g$ defined by $f=\partial_{z}p_{2}/(\mathrm{Pr}\mathrm{Ra})$ (so that $\nabla^{2}f=\partial_{z}^{2}\theta$) and $\nabla^{2}g=\theta$ with the boundary conditions $g=0$ at $z=0$ and

1. These definitions imply that $f$ and $g$ are periodic in $x$ and $y$ and that $f=\partial_{z}^{2}g=g=0$ at $z=0$ and 1. It is possible to compute $\int\theta\partial_{z}p_{2}dV$, where the integration extends over the volume $V$ of a periodicity cell in the layer, by a sequence of integration by parts in which the boundary integrals all vanish:

[TABLE]

We can now divide this equation by the volume $V$ and take the limit $V\rightarrow\infty$ to conclude that

[TABLE]

where $\Delta^{-1}$ denotes the inverse Laplacian for homogeneous Dirichlet boundary conditions.

Before proceeding further, let us look at the form of the optimization problem which yields optimal bounds from eqs. (7-10). There are various ways to formulate such an optimization problem. The formulation presented here is most closely related to the method of auxiliary functions Chernyshenko et al. (2014); Chernyshenko (2017); Tobasco, Goluskin, and Doering (2018). To facilitate the exposition of the numerical implementation later, the formulation presented here is identical to ref. Tilgner, 2017a. Suppose we want to bound some objective function $Z$ which is defined in terms of the velocity and temperature fields. We chose test functions $\varphi_{n}(z)$, $n=1..N$, which depend only on $z$ and on which we project eq. (2):

[TABLE]

We now construct the functional $F(\lambda_{1}...\lambda_{N},\lambda_{R},\lambda_{E},\lambda_{M},\theta,\phi,\psi)$ as

[TABLE]

and replace all time derivatives by their expressions given in eqs. (8), (9), (15), and (10). The task now is to find a set of coefficients $\lambda_{0},\lambda_{1}...\lambda_{N},\lambda_{R},\lambda_{E},\lambda_{M}$ such that the inequality

[TABLE]

holds for all fields $\theta(\bm{r})$, $\phi(\bm{r})$ and $\psi(\bm{r})$ which obey the free slip boundary conditions and $\theta=0$ at $z=0$ and 1. If such a set of $\lambda$’s is found, the above inequality holds in particular for the fields taken from an actual time evolution. Taking the time average of relation (17) thus leads to $\overline{Z}\leq\lambda_{0}-\overline{F}=\lambda_{0}$. The time average of $F$ is zero because it is the linear combination of time derivatives. Replacing the time derivatives by their expressions in eqs. (8), (9), (15), and (10), the best upper bound for $\overline{Z}$ is given by the $\lambda_{0}$ which solves the optimization problem

[TABLE]

with $d^{2}=\langle|\nabla\partial_{x}\psi|^{2}+|\nabla\partial_{y}\psi|^{2}\rangle$. The minimization occurs over all $\lambda$’s and the inequality needs to hold for all $d$ and all eligible $\theta$ and $\phi$.

This is nearly in a form which leads to an SDP after a suitable discretization. The problematic term is $\langle\theta\partial_{z}p_{1}\rangle$ which involves a triple product in terms of $\theta$, $\phi$ and $\psi$. However, this term is readily reduced to a quadratic term by invoking the maximum principle which guarantees that in a statistically stationary state, the temperature takes on values in the interval bounded by the temperatures at the top and bottom planes of the layer, which means that $0\leq T=\theta+1-z\leq 1$, or equivalently, $|T-\frac{1}{2}|\leq\frac{1}{2}$.

To make best use of the maximum principle, we want to compute $\langle\theta\partial_{z}p_{1}\rangle$ in the form $\langle(T-\frac{1}{2})\partial_{z}p_{1}\rangle+\langle(z-\frac{1}{2})\partial_{z}p_{1}\rangle$. This expression can be simplified with the same artifice as before: we first integrate over a volume $V$ with periodic boundary conditions in the horizontal before taking the limit $V\rightarrow\infty$. The second term of the above expression requires us to compute

[TABLE]

in which only the horizontal average of $p_{1}$ appears which obeys

[TABLE]

We can thus solve $\partial_{z}^{2}\int\int dx~{}dy~{}\partial_{z}p_{1}=-\partial_{z}^{2}\int\int dx~{}dy~{}v_{z}^{2}$ with the boundary conditions that $\int\int dx~{}dy~{}\partial_{z}p_{1}=0$ at $z=0$ and 1 to find $\int\int dx~{}dy~{}\partial_{z}p_{1}=-\int\int dx~{}dy~{}v_{z}^{2}$ since $v_{z}=0$ at $z=0$ and 1. Inserting this into eq. (19) yields

[TABLE]

after another integration by parts.

The only triple product left in problem (18) is $\lambda_{M}\langle\theta\partial_{z}p_{1}\rangle$. This can be reduced to a quadratic term only at the price of an inequality:

[TABLE]

The last expression is quadratic in $\bm{v}$ and could be combined with (18) to obtain a convex optimization problem which can be represented as an SDP for numerical purposes. However, the solution of this problem would be very expensive. The success of the previous applications of semidefinite programming to Rayleigh-Bénard convection Tilgner (2017a) relied on the fact that if the dependence in $x$ and $y$ is represented as Fourier series, the equations decouple in the wavenumber of the Fourier modes and only a small number of amplitudes of Fourier modes needs to be taken into account in a numerical computation. The term $\int|\partial_{z}p_{1}|dV$ unfortunately destroys that decoupling.

The decoupling can be restored at the expense of a further inequality. It is shown in the appendix that

[TABLE]

Inserting this into eq. (22), and (22) together with (21) into (18), leads to the following optimization problem for the optimal bound $\lambda_{0}$ of the objective function $\overline{Z}$, where the minimum is sought over all $\lambda$’s:

[TABLE]

The last line implies $|\lambda_{M}|\leq\lambda_{\mathrm{abs}}$. Since the right hand side of the first inequality constraint in (24) is positive, $\lambda_{0}$ is smallest when $\lambda_{\mathrm{abs}}$ is chosen as small as possible, which implies that $\lambda_{\mathrm{abs}}=|\lambda_{M}|$ at optimum. The velocity field is left as a variable in (24) so that it is easier to trace the various terms to the preceding development, but eventually, $v_{z}$ is replaced by $v_{z}=-(\partial_{x}^{2}+\partial_{y}^{2})\phi$ and

[TABLE]

.

Problem (24) contains eq. (10) as a constraint which keeps $\mathrm{Pr}$ as a parameter in the optimization. The effect of $\mathrm{Pr}$ is expected to be strongest for $\mathrm{Pr}\rightarrow\infty$. It will be of interest to also study a reduced optimization problem which results from (24) by introducing $\lambda_{M}^{\prime}=\lambda_{M}\mathrm{Pr}$, $\lambda_{\mathrm{abs}}^{\prime}=\lambda_{\mathrm{abs}}\mathrm{Pr}$. The dissipation is known to be bounded uniformly in $\mathrm{Pr}$, Howard (1963); Busse (1969); Doering and Constantin (1996) so that in the limit $\mathrm{Pr}\rightarrow\infty$, the right hand side of the inequality constraint disappears and we are left with a simpler problem:

[TABLE]

The virtue of this formulation is that it is independent of any sloppiness that may have eased the derivation of eq. (23). The optimization (25) would not be changed by a sharper inequality than (23).

The numerical solution of problems (24) and (25) proceeds in exactly the same way as in ref. Tilgner, 2017a so that a brief summary will suffice here. The variables $\phi$ and $\theta$ are decomposed into $N$ Chebychev polynomials $T_{n}$ for the $z-$direction and into plane waves in $x$ and $y$, as for example in

[TABLE]

The test functions $\varphi_{n}$ are chosen as delta functions centered at the collocation points $z_{n}$ defined by

[TABLE]

Inserting all this into the optimization problems transforms the inequality constraints into the condition that some symmetric matrix be positive semidefinite. This is the standard form of an SDP. The constraints decouple in $k^{2}=k_{x}^{2}+k_{y}^{2}$. Only a small number of wavenumbers actually constrain the solution. This set of wavenumbers is determined automatically Tilgner (2017a). The matrix occurring in the SDP would be much larger if the constraints did not decouple in $k$, hence the effort put into obtaining eq. (23). Some improvements of the basic method accelerate the computation but are not essential. These include a partial integration of eq. (15) and the exploitation of symmetry in $z$. Tilgner (2017a) The resulting SDP was solved with the package cvxopt.

III Results

The optimization problem (24) distinguishes itself from previous problems by the terms multiplied by $\lambda_{M}$ and $\lambda_{\mathrm{abs}}$, which include all the terms containing $\mathrm{Pr}$. In order to test the power of the new constraint, it seems best to choose an objective function which varies dramatically as a function of $\mathrm{Pr}$. The dissipation of the non poloidal components fulfills this criterion, because the flow is purely poloidal at infinite $\mathrm{Pr}$, whereas it contains both poloidal and toroidal components at finite $\mathrm{Pr}$. The choice $Z=d^{2}$ thus promises to be a gratifying objective function.

Call $D$ the optimal $\lambda_{0}$ of problem (24) for $Z=d^{2}$. Fig. 1 shows $D$ as a function of $\mathrm{Pr}$ for different $\mathrm{Ra}$. It is seen that $D$ is constant for low $\mathrm{Pr}$ until it monotonically decreases as a function of $\mathrm{Pr}$ and asymptotes towards a value different from zero at high $\mathrm{Pr}$. The optimization problem (24) thus is not powerful enough to show that the toroidal dissipation disappears at infinite $\mathrm{Pr}$. Let us denote with $D_{0}$ the value of $D$ at small $\mathrm{Pr}$ and with $D_{\infty}$ the limiting value of $D$ as $\mathrm{Pr}$ tends to infinity. $D_{0}$ and $D_{\infty}$ are functions of $\mathrm{Pr}$.

$D_{0}$ is equal to the upper bound one obtains without the constraint derived from $\partial_{t}\langle v_{z}\theta\rangle$, or with $\lambda_{M}=\lambda_{\mathrm{abs}}=0$ in (24). This bound is the same as the one computed in ref. Tilgner, 2017a and obeys approximately $D_{0}=0.021\times\mathrm{Ra}^{3/2}$ at high $\mathrm{Ra}$. The new constraint introduced in this paper becomes active and reduces the upper bound for $\mathrm{Pr}>\mathrm{Pr}_{a}$. From plots like fig. 1 one finds $\mathrm{Pr}_{a}$ as a function of $\mathrm{Ra}$, which is shown in fig. 2. $\mathrm{Pr}$ needs to be larger than $\mathrm{Pr}_{a}=1.18\times\mathrm{Ra}^{1/2}$ to obtain improved bounds.

Just as $D_{0}$ can be computed from a simplified optimization problem, so can $D_{\infty}$ be computed from the reduced problem (25) which has the advantage that it does not depend on $\mathrm{Pr}$ any more, and it is independent of how sharp the estimate in (23) is. $D_{\infty}(\mathrm{Ra})$ is the lowest bound found for any $\mathrm{Pr}$ at a given $\mathrm{Ra}$, so that $(D_{0}-D_{\infty})/D_{0}$ is a measure of the maximum fractional improvement of the bound due to the additional constraint (10) at the given $\mathrm{Ra}$. This fraction is shown in fig. 3. As can be seen from this figure, the constraint (10) can improve the bound by more than a factor of 2 at small $\mathrm{Ra}$, but dishearteningly, the improvement vanishes for large $\mathrm{Ra}$. Bounding methods should be useful at high $\mathrm{Ra}$ when ordinary time integrations become too expensive. But it is precisely in this limit that the newly added constraint does not improve the previously known bound.

The objective function most often considered in the context of optimum theory is $Z=\langle v_{z}\theta\rangle$, which is the Nusselt number minus one. With this choice, the fractional improvement of previous known bounds is found to be less than $2\times 10^{-3}$, and since this is less than the tolerances and errors in the numerical procedure, the improvement could be exactly zero. The potential improvement is at any rate so small that it was not considered worthwhile to determine it accurately, or to show that it vanishes.

IV Conclusion

Previous work has derived bounds on the Nusselt number or other flow quantities in convection from the time evolution equations for $\langle\theta\rangle_{A}$, $\langle\frac{1}{2}\theta^{2}\rangle$, and $\langle\frac{1}{2}\bm{v}^{2}\rangle$. The present paper adds $\langle v_{z}\theta\rangle$ to the list. This additional constraint does not improve the bounds on the Nusselt number by an unambiguously detectable amount. The bounds on the toroidal dissipation on the other hand can be improved by more than $50\%$, but there is no improvement at large $\mathrm{Ra}$. This fact is remarkable because it contradicts the behavior one may intuit from other results in the optimum theory of turbulence. In order to derive bounds for flows in 3D periodic boxes driven by a body force $\bm{f}$, it is necessary to take into account the time evolution equation for $\langle\bm{f}\cdot\bm{v}\rangle$. For convection within the Boussinesq approximation, the momentum equation contains a driving term proportional to $\theta\bm{\hat{z}}$. The term analogous to $\langle\bm{f}\cdot\bm{v}\rangle$ of the 3D periodic box is therefore $\langle v_{z}\theta\rangle$. For an immediate analogy with previous work, the forcing should be solenoidal, so that we should first consider $\langle\bm{v}\cdot(\theta\bm{\hat{z}}-\nabla{\tilde{p}})\rangle$, where ${\tilde{p}}$ is chosen such that $\nabla\cdot(\theta\bm{\hat{z}}-\nabla{\tilde{p}})=0$. However, for a solenoidal velocity field with zero normal component at the boundaries, the expression $\langle\bm{v}\cdot(\theta\bm{\hat{z}}-\nabla{\tilde{p}})\rangle$ reduces to $\langle v_{z}\theta\rangle$. It is therefore surprising that there is not a more striking improvement of bounds for convective flows if the time evolution equation of $\langle v_{z}\theta\rangle$ is included in the optimization problem. At finite $\mathrm{Pr}$, the bounding problem relies on inequality (23) which can possibly be improved. In the limit of large $\mathrm{Pr}$, however, the results are independent of this inequality and no improvement of the results presented here can be achieved without adding yet another constraint.

Appendix A

The goal of this appendix is to prove eq. (23). To this end, we have to compute $\frac{1}{V}\int|\partial_{z}p_{1}|dV$, where $V$ is an arbitrarily large volume in the plane layer, and $p_{1}$ is given by

[TABLE]

with $\partial_{z}p_{1}=0$ at $z=0$ and 1. The variable $q$ is used as an abbreviation for $q=-\nabla\cdot[(\bm{v}\cdot\nabla)\bm{v}]$. We will proceed by finding the Green’s function $G(\bm{r},\bm{r}^{\prime})$ such that $\nabla^{2}G(\bm{r},\bm{r}^{\prime})=\delta(\bm{r}-\bm{r}^{\prime})$ so that $p_{1}(\bm{r})=\int G(\bm{r},\bm{r}^{\prime})q(\bm{r}^{\prime})d^{3}\bm{r}^{\prime}$ where the integral extends over the horizontally infinitely extended layer. Once the Green’s function is known, we can estimate the desired integral from

[TABLE]

and taking the limit $V\rightarrow\infty$:

[TABLE]

Because of the translational invariance, the last factor in eq. (31) simplifies to

[TABLE]

The Green’s function for eq. (29) with Neumann boundary conditions describes for example the potential flow out of a point source in a plane layer. The maximization in expression (32) asks for the distance $z^{\prime}$ from the lower boundary of the layer at which one has to place the point source so that the absolute value of the $z-$component of the velocity integrated over the entire layer is maximum. One intuitively expects that the maximum is reached when the point source is located on one of the boundaries of the layer. It seems very likely that the Green’s function for exactly this problem is given somewhere in the existing literature, but I was not able to find a suitable reference. However, the Green’s function for eq. (29) with the Dirichlet boundary conditions $p_{1}=0$ at $z=0$ and 1 describes the electrostatic potential of a point charge between two parallel infinitely extended metallic plates kept at zero potential, and the Green’s function for this electrostatic problem can be found in the textbook by Jackson Jackson (1999). Adapting the result in this textbook from Dirichtlet to Neumann boundary conditions leads to

[TABLE]

where the position $\bm{r}$ is given in cylindrical coordinates $(s,\varphi,z)$ and $C$ is a constant left unspecified by the Neumann boundary conditions and which is irrelevant because we are only interested in $\partial_{z}G$. $K_{0}$ is a modified Bessel function of the second kind in the usual notation.

We will abstain from formal proofs for two properties which we will deduce from numerical evaluation of eq. (33) and which match the intuition about $G$ one may have from its interpretation as a potential flow out of a point source. The first observation concerns the $z^{\prime}$ at which the maximum in (32) is realized. Fig. 4 shows $\int d^{3}\bm{r}|\partial_{z}G|$ as a function of $z^{\prime}$ computed numerically from the first 40 terms of eq. (33). As expected, this integral is maximal for $z^{\prime}=0$ (or $z^{\prime}=1$ by symmetry). The value of the maximum is numerically close to $1/2$.

This second point can be made more precise with the help of a weaker observation. Let us introduce for brevity $g(\bm{r})$ as the Green’s function corresponding to a point source located at the origin, $g(\bm{r})=G(\bm{r},\bm{0})$. A contour plot of $g(\bm{r})$ in cylindrical coordinates is shown in fig. 5. The second property to extract from numerical calculation and which is again expected from the interpretation of $g(\bm{r})$ as potential flow, is that $\partial_{z}g$ has the same sign everywhere. This second property is important because it implies that we can avail ourselves of the absolute values in expression (32). Using the formulas $\int_{0}^{\infty}sK_{0}(s)ds=1$ and $\sum_{n=1}^{\infty}\frac{1}{n}\sin(n\alpha)=\frac{\pi-\alpha}{2}$, one finds

[TABLE]

so that eq. (31) leads to

[TABLE]

For solenoidal fields $\bm{v}$, this can be brought into a form which is more pleasant for the SDP:

[TABLE]

which in the limit $V\rightarrow\infty$ leads to eq. (23).