1. Introduction
One of the most fascinating features in turbulent flows is the emergence of complicated chaotic structures involving a wide range of length scales behind which ‘typical’ flow patterns are still recognizable. The state of motion is too complex to allow a detailed description of the fluid velocity and experimental or numerical measurements of certain system quantities appear disorganized and unpredictable. Yet, some statistical properties are reproducible (Frisch Reference Frisch1995). One of the challenges in theoretical fluid dynamics is thus the derivation of quantitative statements on turbulent flows. Many of the approaches consist of various approximation procedures, the imposition of physically motivated but ad hoc assumptions (e.g. ‘closure’) or the introduction of scaling hypotheses. Rigorous results beginning with the equations of motion directly are therefore indispensable for checking the validity of the imposed simplifications and for justifying secondary models and theories.
The idea of extracting information about driven turbulent flows via bounds on physical quantities through mathematically justifiable operations and without imposing ad hoc assumptions goes back to the pioneering works of Malkus (Reference Malkus1954), Howard (Reference Howard1972) and Busse (Reference Busse1970, Reference Busse1979). These authors applied variational approaches for the derivation of bounds on the energy dissipation rate in models for shear flow and heat convection. In the 1990s, Constantin and Doering introduced a practical framework for estimating physical quantities rigorously and directly from the equations of motion, which they called the ‘background flow method’ (Doering & Constantin Reference Doering and Constantin1994, Reference Doering and Constantin1996; Constantin & Doering Reference Constantin and Doering1995).
The background flow method is an extremely robust method for constructing bounds in fluid dynamics and builds on techniques developed by Hopf (Reference Hopf1941) to generalize Leray solutions of the Navier–Stokes equations to finite geometries with physical boundary conditions. In this method one manipulates the equations of motion relative to a steady trial background state, which satisfies the forcing conditions. On decomposing the quantity of interest, e.g. the energy dissipation, into a background and a fluctuating component, the background part yields an upper bound if the fluctuation term satisfies a certain non-negativity condition, which is often referred to as the spectral constraint. In a certain sense, finding the least upper bound using Constantin and Doering’s method resembles a variational saddle point problem. The equations imposed on the fluctuations necessarily include the equations of motion. To simplify the derivation of the Euler–Lagrange equations corresponding to the variational problem for the fluctuation term, however, fluctuations are often chosen in a much larger class of functions. In other words, the spectral constraint in the background method is required to hold for an infinite-dimensional set of vector fields, that strictly contains the solutions of the equations of motion. In that case, enforcing the spectral constraint may yield an overestimation of the quantity of primal interest.
After its introduction, the background flow method was the upper bound method with applications ranging from various problems in turbulent heat convection and boundary-force-driven and body-force-driven turbulence to idealized models in magnetohydrodynamics. The method was soon improved by Nicodemus, Grossmann & Holthaus (Reference Nicodemus, Grossmann and Holthaus1997) who introduced an additional balance parameter, and Kerswell (Reference Kerswell1998) showed that this improved method is actually ‘equivalent’ to the approaches of Busse and Howard.
Apart from its practical performance, for many years the background flow method was considered as a rigorous manifestation of Malkus’s marginally stable boundary layer theory (Malkus Reference Malkus1954). The latter is based on the assumption that turbulent boundary-driven flows organize themselves into marginally stable configurations. If the well-mixed core is bounded by thin laminar boundary layers, the thickness of these layers is determined by the condition of marginal stability. The association of the background flow method with Malkus’s theory relies on the surprising observation that the spectral condition in the background flow method resembles a nonlinear stability condition on the background flow. A recent work of C. Nobili & F. Otto (Personal observations), however, proves the failure of this association – at least in the context of infinite Prandtl number Rayleigh–Bénard convection: the authors compute the least upper bound on the Nusselt number (the quantity of interest in Rayleigh–Bénard convection) within the framework of the background flow method. This bound, however, exceeds the bound derived by Otto & Seis (Reference Otto and Seis2011) using completely different methods. In the context of Rayleigh–Bénard convection, it thus seems that this physical interpretation of the background flow method is misleading. Whether the background flow method indeed gives physically relevant information (apart from scaling bounds) in different problems of fluid dynamics can only be speculated.
In this paper, we focus on the energy dissipation rate as an example of one specific physical quantity and present a new method for its rigorous estimation directly from the equations of motion. (In fact, the method has already been introduced in Otto & Seis (Reference Otto and Seis2011), but its universality was not seen at that time.) To allow a straight comparison with the background flow method, the method is tested on the problems considered by Constantin and Doering in Doering & Constantin (Reference Doering and Constantin1994), Constantin & Doering (Reference Constantin and Doering1995), Doering & Constantin (Reference Doering and Constantin1996, Reference Doering and Constantin1998). More precisely, we study the following classical fluid dynamics problems: shear flow, channel flow, Rayleigh–Bénard convection and porous medium convection. All of these problems can be considered as model problems for boundary-force-driven or body-force-driven flows or for thermal convection. In fact, we will recover the same results as Constantin and Doering in the above mentioned papers.
Our new approach presented in this paper is entirely different from the background flow method in that it is based on (local) conservation laws. More precisely, at the heart of our method conservation laws for certain flux components are established. The conservation laws are local in the direction of the symmetry axis and the conserved quantity can be explicitly expressed in terms of the energy dissipation rate. Averaging over small boundary layers and applying elementary estimates yields a bound on the energy dissipation rate. How the new method applies to fluids with different boundary conditions, e.g. no-stress, remains to be seen.
The author believes that the strength of the new method relies on the fact that it uses the equations of motion and some secondary derived physical laws directly instead of working with a rigid ‘upper bound construction’ which in some cases is too restrictive to yield the optimal result, cf. Otto & Seis (Reference Otto and Seis2011). The mathematical operations involved are elementary.
In order to advertise the new method as an alternative to the most widely used background flow method, we choose quick applications at the expense of not optimizing the numerical constants in our bounds. To give a flavour of how competitive results can be obtained by the method, we compute the numerical prefactor only in the first example (shear flow). We moreover caution the reader that our results are only ‘formally’ true in the sense that many of the mathematical operations performed apply only to sufficiently smooth solutions of the Navier–Stokes equations. Of course, the results can be made rigorous if the analysis is performed on suitable weak solutions, e.g. Leray solutions, in the appropriate framework.
The article is organized as follows: after fixing the notation in § 2, we will derive upper bounds on the energy dissipation for shear flows (§ 3), channel flows (§ 4), Rayleigh–Bénard convection (§ 5) and porous medium convection (§ 6).
2. Notation
In what follows, we will try to be as consistent as possible with regard to the notation even though different physical problems will be considered.
Our models are non-dimensionalized.
Our system is a layer of fluid in the box
$[0,L)^{d-1}\times [0,1]$
where
$L$
is an arbitrary positive number that will not enter our analysis. Throughout the article, we will refer to the first
$d-1$
coordinates as the horizontal ones and the last coordinate as the vertical one. We write
$x=(y,z)\in \boldsymbol{R}^{d-1}\times \boldsymbol{R}$
accordingly, and denote by
$\{e_{1},\dots ,e_{d-1},e_{d}\}$
the canonical basis of
$\boldsymbol{R}^{d}=\boldsymbol{R}^{d-1}\times \boldsymbol{R}$
.
We assume periodic boundary conditions in the horizontal directions for all variables involved. The horizontal boundaries are rigid and the imposed conditions will depend on the particular physical problem under consideration.
The fluid velocity is denoted by
$u$
, and we write
$u=(v,w)\in \boldsymbol{R}^{d-1}\times \boldsymbol{R}$
to distinguish the horizontal velocity vector from the vertical component. The hydrodynamic pressure is
$p$
and
$T$
is the temperature field.
We consider the rate of energy dissipation
$$\begin{eqnarray}{\it\varepsilon}=\int _{0}^{1}\langle |\boldsymbol{{\rm\nabla}}u|^{2}\rangle \,\text{d}z,\end{eqnarray}$$
where
$\langle \cdot \rangle$
denotes the horizontal and time average, that is,
$$\begin{eqnarray}\langle \,f\rangle =\lim _{{\it\tau}\uparrow \infty }\frac{1}{{\it\tau}}\int _{0}^{{\it\tau}}\frac{1}{L^{d-1}}\int _{[0,L)^{d-1}}f(t,y)\,\text{d}y\,\text{d}t.\end{eqnarray}$$
In general, the long-time average need not exist, even if finite-time averages are bounded, and we could be more careful at this point by choosing the
$\limsup$
instead of the
$\lim$
. However, for the sake of a clearer statement and to simplify the subsequent analysis, we will be quite formal in most of our computations.
3. Shear flow
As the simplest example of a boundary-driven flow, we consider a fluid which is confined between two parallel plates that are moving at a constant speed relative to each other. The equations of motion are the Navier–Stokes equations in a box
$[0,L)^{d-1}\times [0,1]$
,
$$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}u+u\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}u-{\rm\Delta}u+\boldsymbol{{\rm\nabla}}p=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }u=0, & \displaystyle\end{eqnarray}$$
$u$
and the pressure
$p$
are both periodic in the horizontal variables. At the horizontal plates, we assume no-slip boundary conditions for the velocity field. If the upper boundary plate is moving with constant speed
$\mathit{Re}$
in direction
$e_{1}$
while the lower plate is at rest, we must have 
$$\begin{eqnarray}u=\left\{\begin{array}{@{}ll@{}}0 & \text{for}~z=0\\ \mathit{Re }\,e_{1} & \text{for}~z=1.\end{array}\right.\end{eqnarray}$$
Notice that all quantities are non-dimensionalized and
$\mathit{Re}$
is the Reynolds number.
The scaling of the energy dissipation rate
${\it\varepsilon}$
as a function of the Reynolds number is of fundamental importance in many engineering applications since, for steady states, the energy dissipation rate measures the rate at which work must be done by an agent to keep the upper plate moving.
The interest in understanding the energy dissipation rate in boundary-driven fluid flows as a function of the Reynolds number dates back to Stokes. It is well-known (Serrin Reference Serrin and Truesdell1959; Keller, Rubenfeld & Molyneux Reference Keller, Rubenfeld and Molyneux1967) that solutions to the Stokes equation minimize the dissipation rate among all divergence-free vector fields with a fixed velocity at the boundary. In many situations of physical interest, this solution is laminar and thus continues to exist as a solution of the Navier–Stokes equation. In our case, the laminar solution is the so-called Couette flow
$u_{C}(z)=\mathit{Re}ze_{1}$
which dissipates energy
${\it\varepsilon}_{C}=\mathit{Re}^{2}$
. For small Reynolds numbers,
$u_{C}$
is a stable solution of the Navier–Stokes equation in the sense that any sufficiently small perturbation will decrease in time. For large Reynolds numbers, however, perturbations of the laminar steady state are unstable and the flows can become chaotic or turbulent. The energy dissipation rate is a monotone function of the Reynolds number and acts as a measure of how turbulent the flow is. With regard to the presence of steady regular solutions in the high Reynolds number regime also, absolute lower bounds on the energy dissipation will be dictated by non-turbulent flows:
${\it\varepsilon}\geqslant {\it\varepsilon}_{C}=\mathit{Re}^{2}$
. A rigorous scaling theory for the dissipation rate can hence only be an upper bound theory. An upper bound, however, sets limits on the possible turbulent structures of the flow and is thus an indispensable piece of information in the study of turbulence.
Developing a conventional statistical turbulence theory for high Reynolds numbers, Doering & Constantin (Reference Doering and Constantin1994, appendix A) predict the ‘logarithmic friction law’
$$\begin{eqnarray}{\it\varepsilon}\sim \frac{\mathit{Re}^{3}}{(\log \mathit{Re})^{2}}\quad \text{for}~\mathit{Re}\gg 1,\end{eqnarray}$$
which is in accordance with the experimental data derived e.g. by Lathrop, Fineberg & Swinney (Reference Lathrop, Fineberg and Swinney1992). In the same paper, the authors derive a first upper bound which proves the conjectured rate up to the logarithmic factor, that is, they prove that
${\it\varepsilon}\lesssim \mathit{Re}^{3}$
using the background flow method. In the following, we reproduce Constantin and Doering’s bound using our new method. We will first derive a quick scaling bound and then compute numerical constants in a second step.
As a starting point of our analysis, we recall that the energy dissipation can equally be expressed, for instance, as the trace of the vertical derivative of
$v_{1}=u\cdot e_{1}$
at the top plate, that is,
$$\begin{eqnarray}{\it\varepsilon}=\mathit{Re}\langle \partial _{z}|_{z=1}v_{1}\rangle .\end{eqnarray}$$
This is the averaged energy balance. The identity follows from testing (3.1) with
$u$
, integrating by parts and using (3.2) and (3.3). Moreover, multiplying (3.1) by
$e_{1}$
, taking the horizontal and time average, and using (3.2) and the periodicity, we see that
$$\begin{eqnarray}\langle wv_{1}-\partial _{z}v_{1}\rangle =\text{const.}\quad \text{for all}~z\in [0,1],\end{eqnarray}$$
and thus, in view of (3.3) and (3.5), we obtain the formula
$$\begin{eqnarray}{\it\varepsilon}=\mathit{Re}\langle \partial _{z}v_{1}-wv_{1}\rangle \quad \text{for all}~z\in [0,1].\end{eqnarray}$$
In particular, averaging this identity over some boundary layer
$[0,\ell ]$
, where
$0\leqslant \ell \leqslant 1$
has to be determined later, we have
$$\begin{eqnarray}{\it\varepsilon}=\frac{\mathit{Re}}{\ell }\int _{0}^{\ell }\langle \partial _{z}v_{1}\rangle \,\text{d}z-\frac{\mathit{Re}}{\ell }\int _{0}^{\ell }\langle wv_{1}\rangle \,\text{d}z.\end{eqnarray}$$
On the one hand, an application of Jensen’s inequality yields
$$\begin{eqnarray}\frac{1}{\ell }\int _{0}^{\ell }\langle \partial _{z}v_{1}\rangle \,\text{d}z\leqslant \left(\frac{1}{\ell }\int _{0}^{\ell }\langle (\partial _{z}v_{1})^{2}\rangle \,\text{d}z\right)^{1/2}.\end{eqnarray}$$
On the other hand, with the help of the Cauchy–Schwarz and Poincaré inequalities, we estimate
$$\begin{eqnarray}\displaystyle \frac{1}{\ell }\int _{0}^{\ell }\langle wv_{1}\rangle \,\text{d}z & {\leqslant} & \displaystyle \frac{1}{\ell }\left(\int _{0}^{\ell }\langle w^{2}\rangle \,\text{d}z\int _{0}^{\ell }\langle v_{1}^{2}\rangle \,\text{d}z\right)^{1/2}\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \frac{4\ell }{{\rm\pi}^{2}}\left(\int _{0}^{\ell }\langle (\partial _{z}w)^{2}\rangle \,\text{d}z\int _{0}^{\ell }\langle (\partial _{z}v_{1})^{2}\rangle \,\text{d}z\right)^{1/2}.\end{eqnarray}$$
The optimal Poincaré constant
$4\ell ^{2}/{\rm\pi}^{2}$
is computed in the appendix A. Notice, moreover, that the divergence-free condition (3.2) reveals that
$\langle (\partial _{z}w)^{2}\rangle \leqslant \langle |\boldsymbol{{\rm\nabla}}_{y}v|^{2}\rangle$
. The above estimates can thus be combined into
$$\begin{eqnarray}\frac{{\it\varepsilon}}{\mathit{Re}}\leqslant \left(\frac{1}{\ell ^{1/2}}+\frac{4\ell }{\sqrt{2}{\rm\pi}^{2}}\left(\int _{0}^{1/2}\langle (\partial _{z}w)^{2}+|\boldsymbol{{\rm\nabla}}_{y}v|^{2}\rangle \,\text{d}z\right)^{1/2}\right)\left(\int _{0}^{1/2}\langle (\partial _{z}v_{1})^{2}\rangle \,\text{d}z\right)^{1/2}.\end{eqnarray}$$
To obtain a quick bound, we estimate the integral expressions by
${\it\varepsilon}$
and obtain
$$\begin{eqnarray}{\it\varepsilon}\lesssim \frac{\mathit{Re}}{\ell ^{1/2}}{\it\varepsilon}^{1/2}+\mathit{Re}\ell {\it\varepsilon}.\end{eqnarray}$$
Optimizing the last expression in
$\ell$
yields
$\ell \sim {\it\varepsilon}^{-1/3}\ll 1$
, and thus
${\it\varepsilon}\lesssim \mathit{Re}{\it\varepsilon}^{2/3}$
, which entails
$$\begin{eqnarray}{\it\varepsilon}\lesssim \mathit{Re}^{3}.\end{eqnarray}$$
This estimate agrees with the bound derived by Doering & Constantin (Reference Doering and Constantin1994).
To compute a numerical constant, we optimize (3.11) in
$\ell$
and obtain
$$\begin{eqnarray}\frac{{\it\varepsilon}}{\mathit{Re}}\leqslant \frac{3}{2^{1/6}{\rm\pi}^{2/3}}\left(\int _{0}^{1/2}\langle (\partial _{z}w)^{2}+|\boldsymbol{{\rm\nabla}}_{y}v|^{2}\rangle \,\text{d}z\right)^{1/6}\left(\int _{0}^{1/2}\langle (\partial _{z}v_{1})^{2}\rangle \,\text{d}z\right)^{1/2}.\end{eqnarray}$$
With the help of Young’s inequality
$ab\leqslant (1/p)a^{p}+(1/q)b^{q}$
for any
$p,q\in (1,\infty )$
with
$1/p+1/q=1$
, we further have
$$\begin{eqnarray}\displaystyle \frac{{\it\varepsilon}}{\mathit{Re}} & {\leqslant} & \displaystyle \frac{9{\it\delta}}{2^{4/3}{\rm\pi}^{4/3}}\left(\int _{0}^{1/2}\langle (\partial _{z}w)^{2}+|\boldsymbol{{\rm\nabla}}_{y}v|^{2}\rangle \,\text{d}z\right)^{1/3}+\frac{1}{2{\it\delta}}\int _{0}^{1/2}\langle (\partial _{z}v_{1})^{2}\rangle \,\text{d}z\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \sqrt{\frac{27}{2}}\frac{{\it\delta}^{2}}{{\rm\pi}^{2}}+\frac{1}{2{\it\delta}}\int _{0}^{1/2}\langle (\partial _{z}w)^{2}+|\boldsymbol{{\rm\nabla}}_{y}v|^{2}+(\partial _{z}v_{1})^{2}\rangle \,\text{d}z\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \sqrt{\frac{27}{2}}\frac{{\it\delta}^{2}}{{\rm\pi}^{2}}+\frac{1}{2{\it\delta}}\int _{0}^{1/2}\langle |\boldsymbol{{\rm\nabla}}u|^{2}\rangle \,\text{d}z,\end{eqnarray}$$
where
${\it\delta}>0$
is an arbitrary constant. For symmetry reasons, an analogous estimate can be derived in the upper half of the box
$[0,L)^{d-1}\times [0,1]$
. Altogether we have
$$\begin{eqnarray}\frac{{\it\varepsilon}}{\mathit{Re}}\leqslant \sqrt{\frac{27}{2}}\frac{1}{{\rm\pi}^{2}}{\it\delta}^{2}+\frac{1}{4{\it\delta}}{\it\varepsilon}.\end{eqnarray}$$
Optimizing in
${\it\delta}$
finally yields
$$\begin{eqnarray}\frac{{\it\varepsilon}}{\mathit{Re}^{3}}\leqslant \frac{81\sqrt{6}}{128{\rm\pi}^{2}}\approx 0.157.\end{eqnarray}$$
This bound is a factor of two larger than the bound derived by Doering & Constantin (Reference Doering and Constantin1994). To keep the present paper concise, we do not make any attempt at improving the constant.
4. Channel flow
As the simplest example of a body-force-driven flow, we consider a fluid in a rectangular domain, which is driven by a pressure gradient in the direction of one of the horizontal boundary plates. The problem is modelled by the forced Navier–Stokes equation
$$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}u+u\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}u-{\rm\Delta}u+\boldsymbol{{\rm\nabla}}p=\mathit{Gr}e_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }u=0, & \displaystyle\end{eqnarray}$$
$[0,L)^{d-1}\times [0,1]$
, supplemented with no-slip boundary conditions 
$$\begin{eqnarray}u=0\quad \text{for}~z\in \{0,1\}.\end{eqnarray}$$
Here,
$\mathit{Gr}$
denotes the non-dimensional Grashof number. We remark that when the force is specified a priori as in the present situation, the Reynolds number is an emergent quantity.
The channel flow problem has the laminar solution
$u_{P}(z)=(\mathit{Gr}(z^{2}-z)e_{1})/2$
, the so-called Poiseuille flow, for which the energy dissipation rate is
${\it\varepsilon}_{P}=\mathit{Gr}^{2}/12$
. The Poiseuille flow plays a similar role in the channel flow problem as the Couette flow in the shear flow problem:
$u_{P}$
solves the corresponding Stokes equation and is an unstable solution of the Navier–Stokes equation in the high Grashof number regime. When bounds are expressed in terms of the Grashof number instead of the Reynolds number, the Stokes limit represents an upper bound on the energy dissipation rate:
${\it\varepsilon}\leqslant {\it\varepsilon}_{P}=\mathit{Gr}^{2}/12$
.
The scaling of the energy dissipation rate in the high Grashof number regime is expected to obey the modified logarithmic friction law
$$\begin{eqnarray}{\it\varepsilon}\sim \frac{\mathit{Gr}^{3/2}}{(\log \mathit{Gr})^{2}},\end{eqnarray}$$
similar to the shear flow. The best we can hope for is a lower bound on the energy dissipation rate since the Poiseuille flow always provides an upper bound. A first rigorous lower bound on the energy dissipation rate was established by Constantin & Doering (Reference Constantin and Doering1995),
${\it\varepsilon}\gtrsim \mathit{Gr}^{3/2}$
, which is optimal up to the logarithm. For further improvements we refer to Petrov, Lu & Doering (Reference Petrov, Lu and Doering2005) and references therein. In the following, we will recover this scaling bound by applying our new method.
We see from testing (4.1) with
$u$
, integrating by parts and using (4.2) and (4.3) that the energy dissipation rate is proportional to the average flow velocity
$$\begin{eqnarray}{\it\varepsilon}=\mathit{Gr}\int _{0}^{1}\langle v_{1}\rangle \,\text{d}z.\end{eqnarray}$$
Moreover, multiplying (4.1) by
$e_{1}$
, taking the horizontal and time average, and using (4.2) and periodicity, we obtain
$$\begin{eqnarray}\partial _{z}\langle wv_{1}-\partial _{z}v_{1}\rangle =\mathit{Gr}\quad \text{for all}~z\in [0,1].\end{eqnarray}$$
From the no-slip boundary conditions (4.3) we thus infer that
$$\begin{eqnarray}\langle wv_{1}-\partial _{z}v_{1}\rangle =\mathit{Gr}z-\langle \partial _{z}|_{z=0}v_{1}\rangle \quad \text{for all}~z\in [0,1].\end{eqnarray}$$
On the one hand, averaging in
$z$
over the strip
$[1-\ell ,1]$
of width
$\ell \ll 1$
yields
$$\begin{eqnarray}\mathit{Gr}\sim \frac{1}{\ell }\int _{1-\ell }^{1}\langle v_{1}w-\partial _{z}v_{1}\rangle \,\text{d}z+\langle \partial _{z}|_{z=0}v_{1}\rangle .\end{eqnarray}$$
On the other hand, averaging (4.7) over
$[0,\ell ]$
yields
$$\begin{eqnarray}\langle \partial _{z}|_{z=0}v_{1}\rangle \sim \frac{1}{\ell }\int _{0}^{\ell }\langle \partial _{z}v_{1}-v_{1}w\rangle \,\text{d}z+\ell \mathit{Gr}.\end{eqnarray}$$
Since
$\ell \ll 1$
, a combination of the previous two estimates gives
$$\begin{eqnarray}\mathit{Gr}\lesssim \frac{1}{\ell }\int _{0}^{\ell }\langle \partial _{z}v_{1}-v_{1}w\rangle \,\text{d}z+\frac{1}{\ell }\int _{1-\ell }^{1}\langle v_{1}w-\partial _{z}v_{1}\rangle \,\text{d}z.\end{eqnarray}$$
We can apply the same arguments as in the shear flow case considered in the previous section to deduce
$$\begin{eqnarray}\mathit{Gr}\lesssim \frac{{\it\varepsilon}^{1/2}}{\ell ^{1/2}}+\ell {\it\varepsilon}.\end{eqnarray}$$
The optimal
$\ell$
is of size
$\ell \sim {\it\varepsilon}^{-1/3}$
, which implies that
$$\begin{eqnarray}\mathit{Gr}^{3/2}\lesssim {\it\varepsilon}.\end{eqnarray}$$
In view of (4.5), this bound is equivalent to
$$\begin{eqnarray}\mathit{Gr}^{1/2}\lesssim \int _{0}^{1}\langle v_{1}\rangle \,\text{d}z,\end{eqnarray}$$
which agrees with the bound derived by Constantin & Doering (Reference Constantin and Doering1995) and appears to be sharp to within logarithms.
5. Rayleigh–Bénard convection
Rayleigh–Bénard convection is the transport of heat by thermal convection in a fluid layer that is heated from below and cooled from above. The problem is modelled by the equations of the Boussinesq approximation
$$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}T+u\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}T-{\rm\Delta}T=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }u=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{\mathit{Pr}}(\partial _{t}u+u\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}u)-{\rm\Delta}u+\boldsymbol{{\rm\nabla}}p=\mathit{Ra}\,Te_{d}. & \displaystyle\end{eqnarray}$$
$\mathit{Ra}$
and the Prandtl number
$\mathit{Pr}$
. The equations are complemented by the boundary conditions 
$$\begin{eqnarray}T=1\quad \text{on}~z=0,\quad \text{and}\quad T=0\quad \text{on}~z=1,\end{eqnarray}$$
representing heating at the bottom and cooling at the top, and no-slip boundary conditions for the velocity field,
$$\begin{eqnarray}u=0\quad \text{on}~z\in \{0,1\}.\end{eqnarray}$$
The quantity of interest in this model is the so-called Nusselt number, a measure for the average upward heat flux. It is defined by
$$\begin{eqnarray}\mathit{Nu}=\int _{0}^{1}\langle (uT-\boldsymbol{{\rm\nabla}}T)\boldsymbol{\cdot }e_{d}\rangle \,\text{d}z=\int _{0}^{1}\langle wT-\partial _{z}T\rangle \,\text{d}z.\end{eqnarray}$$
The scaling of the Nusselt number in terms of
$\mathit{Ra}$
and
$\mathit{Pr}$
has been the subject of enormous experimental, numerical and theoretical research for over fifty years. For a recent review, we refer to Ahlers, Grossmann & Lohse (Reference Ahlers, Grossmann and Lohse2009) and references therein. For all values of
$\mathit{Ra}$
and
$\mathit{Pr}$
, a laminar solution is given by
$u_{c}=0$
and
$T_{c}=1-z$
, which corresponds to pure conduction. The corresponding Nusselt number is
$\mathit{Nu}_{c}=1$
, and the laminar solution is unstable for large Rayleigh numbers.
In the high Rayleigh number regime,
$\mathit{Ra}\gg 1$
, the scaling of the Nusselt number is proportional to the scaling of the energy dissipation rate. Indeed, testing (5.3) with
$u$
, using the incompressibility assumption (5.2) and invoking the boundary conditions (5.4) and (5.5) for
$T$
and
$u$
yields
$$\begin{eqnarray}{\it\varepsilon}=\mathit{Ra}(\mathit{Nu}-1)\sim \mathit{Ra}\,\mathit{Nu}.\end{eqnarray}$$
The bound on the Nusselt number in the ultimate turbulent regime is expected to be
$$\begin{eqnarray}\mathit{Nu}\sim \mathit{Ra}^{1/2},\end{eqnarray}$$
if
$\mathit{Ra}\gg 1$
, uniformly in
$\mathit{Pr}$
. In the following, we derive an upper bound on this scaling with the help of different representations of the Nusselt number, similar to the approach in the previous two sections. Averaging the heat equation (5.1) and using periodicity and (5.2), it follows that the heat flux is constant on every horizontal slice, that is
$$\begin{eqnarray}\mathit{Nu}=\langle Tw-\partial _{z}T\rangle \quad \text{for all}~z\in [0,1].\end{eqnarray}$$
Now, averaging this identity over a boundary layer of thickness
$\ell \in [0,1]$
, and using the maximum principle on the temperature,
$\max |T|\leqslant 1$
, which is enforced by the boundary conditions (if not initially, then exponentially fast in time), we obtain that
$$\begin{eqnarray}\mathit{Nu}\leqslant \frac{1}{\ell }\int _{0}^{\ell }\langle wT\rangle \,\text{d}z+\frac{1}{\ell }\leqslant \frac{1}{\ell }\int _{0}^{\ell }\langle |w|\rangle \,\text{d}z+\frac{1}{\ell }.\end{eqnarray}$$
We use Poincaré’s and Hölder’s inequalities and (5.5) to bound the integral over the vertical velocity component by the energy dissipation rate
${\it\varepsilon}$
, that is
$$\begin{eqnarray}\int _{0}^{\ell }\langle |w|\rangle \,\text{d}z\leqslant \ell \int _{0}^{\ell }\langle |\partial _{z}w|\rangle \,\text{d}z\leqslant \ell ^{3/2}\left(\int _{0}^{\ell }\langle (\partial _{z}w)^{2}\rangle \,\text{d}z\right)^{1/2}\leqslant \ell ^{3/2}(\mathit{Ra}\,\mathit{Nu})^{1/2}\end{eqnarray}$$
by (5.7), so that
$$\begin{eqnarray}\mathit{Nu}\leqslant \ell ^{1/2}(\mathit{Ra}\,\mathit{Nu})^{1/2}+\frac{1}{\ell }.\end{eqnarray}$$
Optimizing in
$\ell$
yields that
$\ell \sim (\mathit{Ra}\,\mathit{Nu})^{-1/3}$
, which entails that
$$\begin{eqnarray}\mathit{Nu}\lesssim \mathit{Ra}^{1/2}.\end{eqnarray}$$
This is precisely the same scaling law as derived by Doering & Constantin (Reference Doering and Constantin1996).
Notice that Whitehead & Doering (Reference Whitehead and Doering2011) prove the bound
$\mathit{Nu}\lesssim \mathit{Ra}^{5/12}$
for two-dimensional Rayleigh–Bénard convection with free-stress boundary conditions. The proof relies heavily on the two-dimensional structure (no vortex stretching) and the free-slip boundary conditions (homogeneous vorticity boundary conditions). Whether such a bound can be extended to our problem at hand is not clear to the author. Recent numerical simulations at least indicate that the
$5/12$
scaling should be expected for any solution of (5.1)–(5.3) with finite energy dissipation rate (Hassanzadeh, Chini & Doering Reference Hassanzadeh, Chini and Doering2014).
Applying the same method combined with sophisticated maximal regularity arguments, Chuffrut, Nobili & Otto (Reference Chuffrut, Nobili and Otto2014) recently obtained new bounds on
$\mathit{Nu}$
which improve this bound in certain
$\mathit{Ra}$
–
$\mathit{Pr}$
regimes. The results in particular apply to the large Prandtl number regime. The Nusselt number bound can be interpreted as a bound on the average temperature gradient, cf. (6.9) below. Developing techniques similar to those presented in this paper, the author has derived bounds on higher-order derivatives of the temperature field in infinite Prandtl number convection and estimated deviations of the average vertical temperature profile from linearity (Seis Reference Seis2013).
6. Porous medium convection
We finally consider thermal convection in a porous medium. In this case, Darcy’s law approximates the Navier–Stokes equations, and the Rayleigh–Bénard system (5.1)–(5.3) reduces to
$$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}T+u\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}T-{\rm\Delta}T=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }u=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u+\boldsymbol{{\rm\nabla}}p=(\mathit{Ra})Te_{d}. & \displaystyle\end{eqnarray}$$
$\mathit{Ra}$
is the Rayleigh number. The boundary condition satisfied by the fluid velocity is 
$$\begin{eqnarray}w=0\quad \text{on}~z\in \{0,1\},\end{eqnarray}$$
and we suppose that the container is cooled from above and heated from below, modelled by
$$\begin{eqnarray}T=\left\{\begin{array}{@{}ll@{}}1 & \text{on}~z=0,\\ 0 & \text{on}~z=1.\end{array}\right.\end{eqnarray}$$
Again, we assume periodic boundary conditions in all horizontal directions for all quantities involved.
As in the case of classical Rayleigh–Bénard convection, the quantity of interest here is the Nusselt number
$$\begin{eqnarray}\mathit{Nu}=\int _{0}^{1}\langle wT-\partial _{z}T\rangle \,\text{d}z.\end{eqnarray}$$
Before estimating the energy dissipation rate in this example, we start with a study of the Nusselt number. Experiments and numerics suggest that
$$\begin{eqnarray}\mathit{Nu}\sim \mathit{Ra}\quad \text{for}~\mathit{Ra}\gg 1,\end{eqnarray}$$
cf. Hewitt, Neufeld & Lister (Reference Hewitt, Neufeld and Lister2012). Because of the existence of laminar solutions, we can only expect to prove the upper bound
$\mathit{Nu}\lesssim \mathit{Ra}$
, previously established by Doering & Constantin (Reference Doering and Constantin1998). We first need to establish some alternative identities for
$\mathit{Nu}$
. We first recall that on averaging (6.1), we obtain
$\langle wT-\partial _{z}T\rangle =\text{const.}$
, and thus
$$\begin{eqnarray}\mathit{Nu}=\langle wT-\partial _{z}T\rangle \quad \text{for all}~z\in [0,1].\end{eqnarray}$$
In particular,
$\mathit{Nu}=-\langle \partial _{z}|_{z=0}T\rangle$
thanks to (6.4). Now testing (6.1) with
$T$
, integrating by parts and using (6.2), (6.4) and (6.5), we see that
$$\begin{eqnarray}\mathit{Nu}=\int _{0}^{1}\langle |\boldsymbol{{\rm\nabla}}T|^{2}\rangle \,\text{d}z.\end{eqnarray}$$
On the other hand, by the definition of
$\mathit{Nu}$
and (6.5), testing (6.3) with
$u$
and using (6.2) and (6.4), we obtain
$$\begin{eqnarray}\mathit{Nu}=\frac{1}{\mathit{Ra}}\int _{0}^{1}\langle |u|^{2}\rangle \,\text{d}z+1.\end{eqnarray}$$
We are now able to estimate the Nusselt number. Letting
$\ell \in (0,1)$
be an arbitrary number, we estimate
$$\begin{eqnarray}\mathit{Nu}\leqslant \frac{1}{\ell }\int _{0}^{\ell }\langle wT\rangle \,\text{d}z+\frac{1}{\ell }\end{eqnarray}$$
as in the previous section. Because
$|T|\leqslant 1$
by the maximum principle for the temperature, we obtain via Jensen’s inequality
$$\begin{eqnarray}\mathit{Nu}\leqslant \left(\frac{1}{\ell }\int _{0}^{\ell }\langle w^{2}\rangle \,\text{d}z\right)^{1/2}+\frac{1}{\ell }.\end{eqnarray}$$
By the Nusselt number representation (6.10), this yields
$$\begin{eqnarray}\mathit{Nu}\lesssim \frac{(\mathit{Ra}\,\mathit{Nu})^{1/2}}{\ell ^{1/2}}+\frac{1}{\ell }.\end{eqnarray}$$
This bound is optimized by
$\ell \sim 1$
, so that
$$\begin{eqnarray}\mathit{Nu}\lesssim \mathit{Ra}\end{eqnarray}$$
because
$\mathit{Ra}\gg 1$
.
In a final step, we would like to relate the current bound on the Nusselt number to a bound on the viscous dissipation rate. In fact, we will show that
${\it\varepsilon}\lesssim (\mathit{Ra})^{2}\mathit{Nu}$
, so that the above statement turns into
$$\begin{eqnarray}{\it\varepsilon}\lesssim (\mathit{Ra})^{3}.\end{eqnarray}$$
We start noting that Darcy’s law (6.3) together with the boundary conditions (6.4) and (6.5) provides us with Neumann boundary conditions for the pressure:
$\partial _{z}|_{z=0}p=\mathit{Ra}$
and
$\partial _{z}|_{z=1}p=0$
. In particular, differentiating the horizontal velocity components with respect to
$z$
,
$\partial _{z}v=-\boldsymbol{{\rm\nabla}}_{y}\partial _{z}p$
, multiplying by
$v$
, averaging and integrating by parts yields
$$\begin{eqnarray}\langle v\cdot \partial _{z}v\rangle =-\langle v\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{y}\partial _{z}p\rangle =\langle (\boldsymbol{{\rm\nabla}}_{y}\boldsymbol{\cdot }v)\partial _{z}p\rangle .\end{eqnarray}$$
In particular, the above values for
$\partial _{z}p$
and the horizontal periodicity imply that
$$\begin{eqnarray}\langle v\cdot \partial _{z}v\rangle |_{z=0,1}=0.\end{eqnarray}$$
It thus follows via integration by parts that
$$\begin{eqnarray}\int _{0}^{1}\langle |\boldsymbol{{\rm\nabla}}v|^{2}\rangle \,\text{d}z=\langle v\cdot \partial _{z}v\rangle |_{z=0}^{z=1}-\int _{0}^{1}\langle v\cdot {\rm\Delta}v\rangle \,\text{d}z\stackrel{(6.17)}{=}-\int _{0}^{1}\langle v\boldsymbol{\cdot }{\rm\Delta}v\rangle \,\text{d}z.\end{eqnarray}$$
Now notice that
$-{\rm\Delta}v=(\mathit{Ra})\boldsymbol{{\rm\nabla}}_{y}\partial _{z}T$
. Indeed, taking the divergence of (6.3) yields
${\rm\Delta}p=(\mathit{Ra})\partial _{z}T$
, and thus
$-{\rm\Delta}v=\boldsymbol{{\rm\nabla}}_{y}{\rm\Delta}p=(\mathit{Ra})\boldsymbol{{\rm\nabla}}_{y}\partial _{z}T$
. Therefore, (6.18) becomes
$$\begin{eqnarray}\int _{0}^{1}\langle |\boldsymbol{{\rm\nabla}}v|^{2}\rangle \,\text{d}z=\mathit{Ra}\int _{0}^{1}\langle v\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{y}\partial _{z}T\rangle \,\text{d}z=-\mathit{Ra}\int _{0}^{1}\langle (\boldsymbol{{\rm\nabla}}_{y}\boldsymbol{\cdot }v)\partial _{z}T\rangle \,\text{d}z.\end{eqnarray}$$
Since
$|\boldsymbol{{\rm\nabla}}_{y}\boldsymbol{\cdot }v|\lesssim |\boldsymbol{{\rm\nabla}}_{y}v|\leqslant |\boldsymbol{{\rm\nabla}}v|$
, we can use the Cauchy–Schwarz inequality to obtain
$$\begin{eqnarray}\int _{0}^{1}\langle |\boldsymbol{{\rm\nabla}}v|^{2}\rangle \,\text{d}z\lesssim \mathit{Ra}\left(\int _{0}^{1}\langle |\boldsymbol{{\rm\nabla}}v|^{2}\rangle \,\text{d}z\int _{0}^{1}\langle |\boldsymbol{{\rm\nabla}}T|^{2}\rangle \,\text{d}z\right)^{1/2},\end{eqnarray}$$
and thus, via (6.9),
$$\begin{eqnarray}\int _{0}^{1}\langle |\boldsymbol{{\rm\nabla}}v|^{2}\rangle \,\text{d}z\lesssim (\mathit{Ra})^{2}\mathit{Nu}.\end{eqnarray}$$
The estimate of the vertical velocity component is easier because of (6.4): testing
${\rm\Delta}w=(\mathit{Ra}){\rm\Delta}T-\partial _{z}{\rm\Delta}p=(\mathit{Ra}){\rm\Delta}_{y}T$
with
$w$
, integrating by parts and using the Cauchy–Schwarz inequality yields
$$\begin{eqnarray}\int _{0}^{1}\langle |\boldsymbol{{\rm\nabla}}w|^{2}\rangle \,\text{d}z\lesssim (\mathit{Ra})^{2}\mathit{Nu}.\end{eqnarray}$$
Combining the last two estimates finally yields
${\it\varepsilon}\lesssim (\mathit{Ra})^{2}\mathit{Nu}$
as desired.
Acknowledgements
Part of this work was performed when the author was visiting the MPI MIS in Leipzig. The author is grateful for the hospitality. He further acknowledges stimulating discussions with C. Nobili, C. Doering, J. Whitehead and R. Wittenberg.
Appendix A
In this appendix, we derive the Poincaré inequality
$$\begin{eqnarray}\int _{0}^{\ell }f^{2}\,\text{d}z\leqslant \frac{4\ell ^{2}}{{\rm\pi}^{2}}\int _{0}^{\ell }(f^{\prime })^{2}\,\text{d}z,\end{eqnarray}$$
for all functions
$f$
such that
$f(0)=0$
. This estimate is sharp as can be seen by choosing
$f(z)=\sin ({\rm\pi}/2\ell )z$
. Notice that (A 1) is equivalent to the Poincaré inequality
$$\begin{eqnarray}\int _{-\ell }^{\ell }g^{2}\,\text{d}z\leqslant \frac{4\ell ^{2}}{{\rm\pi}^{2}}\int _{-\ell }^{\ell }(g^{\prime })^{2}\,\text{d}z,\end{eqnarray}$$
for all functions
$g$
satisfying
$\int _{-\ell }^{\ell }g\,\text{d}z=0$
in the sense that the Poincaré constants agree. Indeed, if
$f$
satisfies (A 1) and
$g$
is obtained from
$f$
by odd reflection at
$z=0$
, then
$g$
has mean zero on
$(-\ell ,\ell )$
and satisfies (A 2). On the other hand, if
$g$
is an odd function for which (A 2) holds and
$f$
is the restriction of
$g$
to
$(0,\ell )$
then
$f$
satisfies (A 1).
The statement in (A 2) is equivalent to the variational problem
$$\begin{eqnarray}\frac{{\rm\pi}^{2}}{4\ell }=\min \left\{\int _{-\ell }^{\ell }(g^{\prime })^{2}\,\text{d}z:\int _{-\ell }^{\ell }g\,\text{d}z=0,\int _{-\ell }^{\ell }g^{2}\,\text{d}z=\ell \right\},\end{eqnarray}$$
and the expression on the right is minimized by
$g(z)=\sin ({\rm\pi}/2\ell )z$
.
