Optimal wall shapes and flows for steady planar convection

Abstract

We compute steady planar incompressible flows and wall shapes that maximize the rate of heat transfer ($Nu$) between hot and cold walls, for a given rate of viscous dissipation by the flow ($Pe^2$), with no-slip boundary conditions at the walls. In the case of no flow, we show theoretically that the optimal walls are flat and horizontal, at the minimum separation distance. We use a decoupled approximation to show that flat walls remain optimal up to a critical non-zero flow magnitude. Beyond this value, our computed optimal flows and wall shapes converge to a set of forms that are invariant except for a $Pe^{-1/3}$ scaling of horizontal lengths. The corresponding rate of heat transfer $Nu \sim Pe^{2/3}$. We show that these scalings result from flows at the interface between the diffusion-dominated and convection-dominated regimes. We also show that the separation distance of the walls remains at its minimum value at large $Pe$.

1. Introduction

The transfer of heat from solid boundaries to adjacent fluid flows is fundamental to many natural and technological processes. Important examples include the formation and evolution of stars and planets, and energy consumption in buildings, computer heat sinks and organisms (Blundell & Blundell 2010; Bergman et al. 2011). For flows with large spatial scales and/or large flow speeds, the heat transfer is often controlled by thin viscous and thermal layers at the solid boundaries. Consequently, the rate of heat transfer is sensitive to the features of the solid boundaries, and in particular to their shape or geometry (Webb & Kim 2005; Lienhard 2013; Tobasco 2022; Wen, Goluskin & Doering 2022b; Wen et al. 2022a; Song, Fantuzzi & Tobasco 2023).

The effect of wall shape modulations or roughness on thermal transport has been studied many times for both natural and forced convection. Toppaladoddi, Succi & Wettlaufer (2017) considered natural convection in a fluid layer between hot and cold sinusoidal surfaces. They fixed the roughness amplitude at one-tenth the fluid layer height, and found that heat transfer was maximized when the roughness amplitude and wavelength were similar. This is one example of a large body of numerical and experimental work on the effect of wall roughness on the rate of heat transfer in turbulent natural convection with various roughness geometries – sinusoidal (Toppaladoddi, Succi & Wettlaufer 2015; Toppaladoddi et al. 2017; Zhu et al. 2017), triangular (Zhang et al. 2018), cubic (Rusaouën et al. 2018), ratchet-shaped (Jiang et al. 2018), ring-shaped (Emran & Shishkina 2020) and fractal (Toppaladoddi et al. 2021) and other multiscale roughness profiles (Zhu et al. 2019; Sharma, Chand & De 2022). A main focus of these works is the asymptotic scaling of the rate of heat transfer with the Rayleigh number and its dependence on the geometric form of the roughness and its length scales including wavelengths and heights of the roughness profiles (Yang et al. 2021). An experimental work that considered pyramid-shaped roughness elements of varying aspect ratio related the heat transfer enhancement to the dynamics of thermal plumes near the roughness elements (Xie & Xia 2017).

Variations of these problems include the effect of tilting the rough walls, for triangular (Chand, Sharma & De 2022) and ratchet-shaped (Jiang et al. 2023) surfaces, and the heat transfer enhancement that can be obtained by moving the rough plates (Jin et al. 2022). Zhang et al. (2018) showed computationally that in some cases with small roughness heights, roughness can actually decrease the rate of heat transfer. On the theoretical side, Goluskin & Doering (2016) used the ‘background method’ to obtain upper bounds on the rate of heat transfer in fluid layers between upper and lower rough walls whose profiles correspond to single-valued functions of the horizontal coordinate. Bleitner & Nobili (2022) used a similar approach to derive upper bounds on heat transfer for Navier-slip rough boundaries.

A related body of work has relaxed the requirement that the flow solve the Boussinesq equations for natural convection, and instead optimized the rate of heat transfer over the larger class of all incompressible flows between hot and cold boundaries (Hassanzadeh, Chini & Doering 2014; Souza 2016; Tobasco & Doering 2017; Marcotte et al. 2018; Motoki, Kawahara & Shimizu 2018a,b; Doering & Tobasco 2019; Souza, Tobasco & Doering 2020; Kumar 2022; Alben 2023). These works consider a simple geometry – inspired by Rayleigh–Bénard convection – consisting of a layer of fluid between flat horizontal walls at different fixed temperatures. Optimal flows have also been calculated for domains obtained by conformal mappings (Alben 2017b) and in flows through channels (Alben 2017a). All incompressible flows are solutions to the incompressible Navier–Stokes equations with a certain distribution of force per unit volume applied over the flow domain (Alben 2017a). Although the forcing distribution may be difficult to create exactly in an experiment, simple approximate distributions may be sufficient to obtain large improvements from previous flows (Alben 2017a). The flows can be used to identify typical features of efficient flows, such as branching structures (Zimparov, Da Silva & Bejan 2006; Tobasco & Doering 2017; Motoki et al. 2018a; Kumar 2022; Alben 2023).

A large body of engineering work has studied the effect of wall roughness in simple forced convection scenarios. The geometry of walls separating two fluids in a heat exchanger is often manipulated to increase the rate of heat transfer for a given amount of power needed to drive the flow (Webb & Kim 2005; Bergles & Manglik 2013). A typical strategy is the insertion of helical coils in tubes with axial flow (Gee & Webb 1980). Also well studied is the enhancement of heat transfer due to sinusoidal wall modulations in channel flow, and the fluid-dynamical mechanisms underlying the enhancement (Stone, Stroock & Ajdari 2004; Kanaris, Mouza & Paras 2006; Guzman et al. 2009; Castelloes, Quaresma & Cotta 2010; Muthuraj & Srinivas 2010; Sui et al. 2010; Gong et al. 2011; Ramgadia & Saha 2013; Grant Mills et al. 2014).

In this work we again consider optimal incompressible flows for heat transfer between hot and cold walls (Hassanzadeh et al. 2014; Souza 2016; Tobasco & Doering 2017; Motoki et al. 2018a; Souza et al. 2020; Kumar 2022; Alben 2023), but now in the presence of wall roughness. In particular, we extend the computational approach of Alben (2023) to search for optimal flows together with optimal boundary wall shapes, in order to maximize the rate of heat transfer ($Nu$) given a certain rate of power consumption by the flow ($Pe^2$). At zero power consumption, flat walls are shown to be optimal theoretically in § 3. Below a certain power consumption rate, the computed optimal flows are rectangular convection rolls with flat walls, as in previous studies (Souza et al. 2020; Alben 2023). In § 4 this is shown theoretically by showing that the leading-order effects of the flows and the roughness are decoupled when both are small.

In § 5 we describe the numerical methods for computing temperature fields and optimal flows. Using these methods, in § 5.1 we show the solution features – sharp temperature gradients – that limit the accuracy of the computations, and how the accuracy depends on key wall shape parameters. We also use the computations to test the accuracy of the decoupled approximation at small $Pe$ and wall deflection amplitudes, in § 5.2.

We present computed optima at moderate and large $Pe$ in § 6. Above a critical $Pe$, the computed optima with wavy walls outperform those with flat walls. At large power consumption rates, the optimal flow streamlines and wall shapes fall within a large but well-defined set of typical configurations. The configurations are invariant at large $Pe$, except for a power-law scaling of their horizontal period together with an $O(1)$ vertical roughness of the walls. Thus the convection rolls are very elongated vertically in the limit of large $Pe$. The scaling of the horizontal period, $L_x \sim Pe^{-1/3}$, is close to that seen for the flat-wall optima in previous work, and is shown here to result in flows at the interface between diffusion-dominated and convection-dominated regimes. The corresponding heat transfer rate $Nu$ scales as $Pe^{2/3}$, which was shown to be an upper bound in the flat-wall case (Souza 2016) in two-dimensional (2-D) and three-dimensional (3-D) flows. The scaling was observed numerically (Motoki et al. 2018a) and shown analytically (Kumar 2022) for 3-D flows with flat walls but only up to logarithmic corrections in 2-D flows with flat walls (Tobasco & Doering 2017). The $Pe^{2/3}$ scaling corresponds to an upper bound for flows with rough walls shown theoretically by Goluskin & Doering (2016).

In § 7 we explain why the optima have the minimum possible separation between the hot and cold walls. In § 8 we present results from alternative optimization problems, and § 9 presents the conclusions and additional context for the results.

2. Model

A 2-D layer of fluid is contained in the region $-\infty < x < \infty$ and $y_{bot}(x) < y < y_{top}(x)$ (see figure 1). The curvilinear walls at $y = y_{bot}(x)$ and $y = y_{top}(x)$ have temperatures 1 and 0, respectively. The two walls are periodic in $x$ with a common period $L_x$ that is determined during the optimization. The walls (black lines in figure 1) are constrained so $y_{bot} \leq 0$ and $y_{top} \geq 1$, respectively; i.e. the walls may not cross the blue dashed lines in figure 1. Without this constraint, the distance between the walls could approach zero and then the temperature gradient would diverge, generally along with the rate of heat transfer, regardless of the fluid flow. The constraint allows us to ask which aspects of the walls’ shapes other than their separation distance are useful for heat transfer by the intervening fluid. One can also interpret the constraint as choosing a characteristic length scale for the problem, which is $H_{min}$, the minimum allowed separation distance between the walls. The separation distance $H \equiv \min y_{top}(x) - \max y_{bot}(x)$ is constrained to be ${\geq H_{min}}$, or 1 when we non-dimensionalize all distances by $H_{min}$, going forward. Although $H > 1$ is possible, we find that $H = 1$ for all computed optima. Intuitively, minimizing the separation distance between the walls seems advantageous for heat transfer, and in § 7 we use a scaling argument to explain why this choice might be optimal. Also, the choice of wall temperatures (1 and 0) is equivalent to setting the characteristic temperature scale to be the difference between the wall temperatures. In addition to the constraint $H \geq 1$, the walls’ vertical positions must be single-valued functions of $x$ but are otherwise essentially arbitrary.

Figure 1. Example of a domain that is periodic in $x$ with period $L_x$. The bottom boundary is $y_{bot}(x) \leq 0$ and the top boundary is $y_{top}(x) \geq 1$. The boundaries do not cross the blue dashed lines. The distance between the walls’ inward extrema is $H \geq H_{min}$. The grey region is indicated for the discussion in § 3.

Between the walls, the temperature field is determined by the steady advection–diffusion equation

(2.1)\begin{equation} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} T - \nabla^2T = 0, \end{equation}

driven by an incompressible flow field $(u,v) = \boldsymbol {u} = (\partial _y\psi,-\partial _x\psi )$, with $\psi (x,y)$ the stream function. We have set the prefactor of $\nabla ^2T$ (the thermal diffusivity) to unity, corresponding to non-dimensionalizing velocities by $\kappa /H_{min}$, with $\kappa$ the dimensional thermal diffusivity.

The set-up is similar to Rayleigh–Bénard convection, but instead of solving the Boussinesq equations for the flow field, we consider all incompressible flows and find those that maximize the Nusselt number, the heat transfer from the hot lower boundary (per unit horizontal distance):

(2.2)\begin{equation} Nu = \frac{1}{L_x}\int_{C = \{(x, y_{bot}(x))\}} \partial_n T \,{\rm d}s . \end{equation}

Integrating $\nabla ^2 T$ over the whole region and using periodicity at the side walls, we see that $Nu$ is also the heat transfer into the upper (cold) boundary.

We maximize $Nu$ over flows with a given spatially averaged rate of viscous dissipation $Pe^2$:

(2.3)\begin{align} Pe^2 &\equiv \frac{1}{L_x}\int_0^{L_x}\int_{y_{bot}(x)}^{y_{top}(x)} 2\boldsymbol{\mathsf{E}}:\boldsymbol{\mathsf{E}} \, {{\rm d}\kern0.05em y} \, {{\rm d}\kern 0.06em x} \nonumber\\ &=\frac{1}{L_x}\int_0^{L_x}\int_{y_{bot}(x)}^{y_{top}(x)} (\nabla^2 \psi)^2 \, {{\rm d}\kern0.05em y} \, {{\rm d}\kern 0.06em x},\quad \boldsymbol{\mathsf{E}} \equiv \frac{1}{2}(\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^T). \end{align}

Here, $Pe^2$ is the average rate of viscous dissipation per unit horizontal width (Batchelor 1967), non-dimensionalized by $\mu \kappa ^2/H_{min}^3$, with $\mu$ the viscosity. Tensor $\boldsymbol{\mathsf{E}}$ is the rate-of-strain tensor and $\boldsymbol{\mathsf{E}}:\boldsymbol{\mathsf{E}}$ is the ‘matrix scalar product’, i.e. the sum of the squares of the entries of $\boldsymbol{\mathsf{E}}$. The expression for $Pe^2$ in terms of $\psi$ in (2.3) holds for incompressible flows given the no-slip boundary conditions at the top and bottom walls and periodicity in $x$ we have here (Lamb 1932; Milne-Thomson 1968).

An incompressible flow field that maximizes $Nu$ for a given $Pe$ can also be considered a solution to the Boussinesq or Navier–Stokes equations with a suitable forcing function – whatever forcing is needed to balance the remaining terms in the equations. The optimization problem here is essentially the same as in Souza et al. (2020) and Alben (2023) except that the walls are no longer flat and horizontal but must be determined together with the flow field and the horizontal period.

Like the walls, the flow $\psi$ is periodic in $x$ with period $L_x$. In terms of $\psi$, the no-slip boundary conditions are $\psi = \partial _n\psi = 0$ at $y = y_{bot}(x)$ and $\psi = \psi _{top}$ and $\partial _n\psi = 0$ at $y = y_{top}(x)$. The constant $\psi _{top}$ is the net horizontal fluid flux through the channel. In most of our computed optima including those with the highest $Nu$ at a given $Pe$, $\psi _{top}$ is very close to 0.

3. No convection ($Pe = 0$)

First we consider the limiting case $Pe = 0$, the pure-conduction problem. We show that the maximum $Nu$ is achieved with flat walls with the minimum separation, $y_{bot}(x) \equiv 0$ and $y_{top}(x) \equiv 1$. Using $\nabla ^2 T = 0$, we show that the net vertical conductive heat flux across all horizontal cross-sections is constant over the portion of $y$ values that the walls do not cross, $0 \leq y \leq 1$:

(3.1)\begin{equation} \frac{{\rm d}}{{\rm d}\kern0.05em y}\int_0^{L_x} \partial_y T(x,y)\, {{\rm d}\kern 0.06em x} = \int_0^{L_x} \partial_{yy} T(x,y) \,{{\rm d}\kern 0.06em x} = \int_0^{L_x} -\partial_{xx} T(x,y) \,{{\rm d}\kern 0.06em x} = 0, \quad 0 \leq y \leq 1, \end{equation}

using the periodicity in $x$ for the last equality. Thus

(3.2)\begin{align} \int_0^{L_x} \partial_y T {{\rm d}\kern 0.06em x} &= \{ \text{constant in } y: 0 \leq y \leq 1 \} = \int_0^1 \int_0^{L_x} \partial_y T \,{{\rm d}\kern 0.06em x} \, {{\rm d}\kern0.05em y} \end{align}

(3.3)\begin{align} & = \int_0^{L_x} \int_0^1 \partial_y T \,{{\rm d}\kern0.05em y} \, {{\rm d}\kern 0.06em x} = \int_0^{L_x} T(x,1) - T(x,0)\, {{\rm d}\kern 0.06em x}. \end{align}

If, at some point $x_0$, $y_{bot}(x_0) < 0$ or $y_{top}(x_0) > 1$, then for some open $x$ interval containing $x_0$ we have $T(x,0) < 1$ or $T(x,1) > 0$, respectively, by the maximum principle for harmonic functions – the maximum and minimum of $T$ are taken on the boundary. In that case $T(x,0) - T(x,1)$, which is $\leq$1 for any wall shape by the maximum principle, is in fact strictly less than 1 on some $x$ interval. Therefore its average (over $x$) is strictly less than 1:

(3.4)\begin{equation} 1 > \frac{1}{L_x}\int_0^{L_x} T(x,0) - T(x,1) \, {{\rm d}\kern 0.06em x} = \frac{1}{L_x}\int_0^{L_x} -\partial_y T|_{y = 0}\, {{\rm d}\kern 0.06em x}. \end{equation}

The last equality in (3.4) follows by (3.2)–(3.3). We can show that the last quantity in (3.4) is in fact $Nu$, by again integrating $\nabla ^2 T$, this time over the grey region in figure 1. We have

(3.5)\begin{gather} 0 = \frac{1}{L_x}\iint \nabla^2 T \,{\rm d}A = \frac{1}{L_x} \oint \partial_n T \,{\rm d}s, \end{gather}

(3.6)\begin{gather}0 = \frac{1}{L_x}\int_{C = \{(x, y_{bot}(x))\}} \partial_n T \,{\rm d}s + \frac{1}{L_x} \int_0^{L_x} \partial_y T|_{y = 0} \,{{\rm d}\kern 0.06em x}. \end{gather}

The contributions to the closed contour integral in (3.5) from the vertical sides of the grey region cancel by $x$ periodicity. Combining (2.2), (3.4) and (3.6), we have $Nu < 1$ for non-flat walls, whereas $Nu = 1$ for flat walls (in which case $-\partial _y T \equiv 1$), the optimal solution with no flow ($Pe = 0$).

4. Small convection ($0 < Pe \ll 1$)

Now we consider the case $0 < Pe \ll 1$. We denote the deviations of the bottom and top walls from their flat states by $H_1$ and $H_2$:

(4.1a,b)\begin{equation} H_1(x) \equiv{-}y_{bot} \geq 0,\quad H_2(x) \equiv y_{top}(x) - 1 \geq 0. \end{equation}

For this asymptotic analysis and for the subsequent computations, we define $(p,q)$ coordinates, in which the flow domain is fixed as a unit square, for all wall shapes and all $L_x$:

(4.2a,b)\begin{equation} p \equiv \frac{x}{L_x} ,\quad q \equiv \frac{y - y_{bot}(x)}{y_{top}(x) - y_{bot}(x)} = \frac{y + H_1(x)}{1+H_1(x) + H_2(x)}. \end{equation}

In Appendix A we give further details on how the objective and constraint equations, (2.1)–(2.3), are written in $(p,q)$ coordinates.

In the limit $Pe \to 0$, the optimal combination of flow and wall shapes must have $H_1, H_2 = o(1)$, i.e. tend to case of flat walls, by the argument in the previous section. Otherwise $Nu$ would be strictly less than 1 in the limit, whereas the flat-wall case has $Nu \geq 1$ for any incompressible flow field (Tobasco & Doering 2017).

When we write the advection–diffusion equation (2.1) in $(p,q)$ coordinates, the small perturbation to flat walls corresponds to adding small space-varying coefficients in front of the differential operators. At small $Pe$ the flow also takes the form of a small space-varying coefficient in (2.1). We write

(4.3a–c)\begin{equation} \boldsymbol{u} = \boldsymbol{u}_1 + O(Pe^2),\quad \nabla_{x,y} = \nabla_{p,q} + \nabla_1 + \cdots,\quad \nabla_{x,y}^2 = \nabla_{p,q}^2 + \nabla^2_1 + \cdots, \end{equation}

where $\boldsymbol {u}_1 = O(Pe)$ and $\nabla _1$ and $\nabla ^2_1$ are first- and second-order differential operators in $p$ and $q$ with coefficients that are linear in $H_1$ and $H_2$, and therefore $o(1)$ as $Pe \to 0$. Terms with quadratic and higher-order coefficients are subdominant and not written explicitly. Given the flow and wall shape, the optimal temperature field has an expansion

(4.4)\begin{equation} T = T_0 + T_1 + \cdots, \end{equation}

with $T_0 = 1-q$ and $T_1 = o(1)$. We insert the expansions of $\boldsymbol {u}$, $T$ and the differential operators into the advection–diffusion equation (2.1). At $O$($Pe^0$),

(4.5)\begin{equation} -\!\nabla_{p,q}^2 T_0 = 0, \end{equation}

so $T_0 = 1-q$, the flat-wall steady-conduction solution. We then take the next-order terms, linear in $Pe$ or the wall-perturbation amplitude:

(4.6)\begin{equation} -\!\nabla_{p,q}^2 T_1 = \nabla_1^2 T_0 - \boldsymbol{u}_1 \boldsymbol{\cdot} \nabla_{p,q} T_0. \end{equation}

Thus $T_1$ satisfies Poisson's equation forced by a superposition of the wall perturbation and the small flow. We decompose $T_1$ as the sum of the perturbations from each forcing term separately:

(4.7)\begin{gather} T_1 = T_{1A} + T_{1B}, \end{gather}

(4.8)\begin{gather}-\nabla_1^2 T_{1A} = \nabla_1^2 T_0 , \end{gather}

(4.9)\begin{gather}-\nabla_1^2 T_{1B} ={-} \boldsymbol{u}_1 \boldsymbol{\cdot} \nabla_{p,q} T_0. \end{gather}

Here $T_{1A}$ is the leading-order change in the temperature field with the wall perturbation and zero flow, while $T_{1B}$ is the leading-order change in the temperature field with the $O(Pe)$ flow and flat walls.

The Nusselt number

(4.10)\begin{equation} Nu = \frac{1}{L_x}\int_0^{L_x}\hat{\boldsymbol{n}}\boldsymbol{\cdot} \boldsymbol{\nabla} T \left.\frac{{\rm d}s}{{\rm d}\kern 0.06em x} \right|_{y = y_{bot}(x)} \,{{\rm d}\kern 0.06em x}= \int_0^{1} \hat{\boldsymbol{n}}\boldsymbol{\cdot} \boldsymbol{\nabla} T \left.\frac{{\rm d}s}{{\rm d}\kern 0.06em x} \right|_{q=0} \,{\rm d}p \end{equation}

is expanded similarly, using the expansion for $T$. The wall perturbation results in an expansion:

(4.11)\begin{equation} \frac{{\rm d}s}{{\rm d}\kern 0.06em x}\hat{\boldsymbol{n}}\boldsymbol{\cdot} \boldsymbol{\nabla} |_{q=0} =\partial_q + \left(\frac{{\rm d}s}{{\rm d}\kern 0.06em x}\hat{\boldsymbol{n}}\boldsymbol{\cdot} \boldsymbol{\nabla}\right)_1 + \cdots. \end{equation}

The subscript 1 again denotes the leading-order perturbation from the flat-wall term, $\partial _q$ in this case since ${\rm d}s/{{\rm d}\kern 0.06em x} = 1$ for the flat wall. The leading-order terms in the expansion of $Nu$ are

(4.12)\begin{align} Nu &= \int_0^{1} \partial_q T_0 |_{q=0} \,{\rm d}p + \left.\int_0^{1} \left(\frac{{\rm d}s}{{\rm d}\kern 0.06em x}\hat{\boldsymbol{n}}\boldsymbol{\cdot} \boldsymbol{\nabla}\right)_1 T_0 \right|_{q=0} \,{\rm d}p\nonumber\\ &\quad + \int_0^{1} \partial_q T_{1A} \,{\rm d}p + \int_0^{1} \partial_q T_{1B} \,{\rm d}p + \cdots \end{align}

(4.13)\begin{align} &= Nu_0 + Nu_{1A} + Nu_{1B} + \cdots. \end{align}

The first integral in (4.12) is $Nu_0$, the purely conductive heat transfer with flat walls, and equals 1. The sum of the second and third integrals in (4.12) is $Nu_{1A}$, the leading-order correction due to non-flat walls, and the fourth integral in (4.12) is $Nu_{1B}$, the leading-order correction due to the flow. For small wall perturbations and small $Pe$, the leading-order correction due to non-flat walls, $Nu_{1A}$, is independent of the flow. It is the change in the pure conduction heat transfer due to non-flat walls. In the previous section ($Pe = 0$) it was shown that this change can only decrease $Nu$, i.e. $Nu_{1A} < 0$.

Furthermore, the leading-order correction due to the flow, $Nu_{1B}$, is independent of the wall perturbation. It is the leading-order change in $Nu$ due to the optimal flow with flat walls at small $Pe$, which is a periodic array of almost square convection rolls (Hassanzadeh et al. 2014; Souza 2016; Alben 2023). It turns out that even though $T_{1B} = O(Pe)$, $Nu_{1B}$ is smaller, $O(Pe^2)$ (Hassanzadeh et al. 2014; Souza et al. 2020; Alben 2023). However, any $o(1)$ form of $Nu_{1B}$ is consistent with the expansion (4.13), which gives the decoupled effect of the walls and flow at leading order. The other main ingredient for (4.13) was that the optimal wall perturbation amplitude is $o(1)$ as $Pe \to 0$, as shown previously.

The main consequence of (4.13) that we highlight is that non-flat walls can only be optimal when $Pe$ is large enough for the remainder terms in (4.13) to be comparable in magnitude to $Nu_{1B}$, so they outweigh its negative effect on $Nu$. Therefore, the transition to optima with non-flat walls occurs only above some critical $Pe > 0$ and not at arbitrarily small $Pe$. The same expansion can be used to show that a flat wall shape maximizes $Nu$ for sufficiently weak natural convection, i.e. at a Rayleigh number $Ra$ sufficiently close to $Ra_c$, the critical value for the Rayleigh–Bénard instability (Drazin & Reid 2004). The argument requires that $\boldsymbol {u} = o(1)$ as $Ra \to Ra_c$.

In the next section, we describe our computational methods for solving for temperature fields, testing the accuracy of the decoupled approximation (4.13), and computing optimal flows and wall shapes across the transition to non-flat walls.

5. Computational methods

In order to solve (2.1) with a variety of wall shapes and horizontal periods, we use $(p,q)$ coordinates (4.2) so the computational domain is square. We write each of (2.1)–(2.3) in $(p,q)$ coordinates, and the differential operators and integrals in $(x,y)$ are replaced with differential operators and integrals in $(p,q)$ multiplied by functions of $L_x$, $y_{bot}(p)$, $y_{top}(p)$ and $q$ and their first and second derivatives with respect to $p$ and $q$.

We define the wall deformations to be single-signed and bounded using auxiliary functions $h_1$ and $h_2$:

(5.1)\begin{equation} \left.\begin{array}{c} \displaystyle y_{bot}(p) ={-}A \left(\dfrac{1+\sin{A_1}}{2}\right)\dfrac{h_1(p)^2}{\|h_1(p)^2\|_6},\\ \displaystyle y_{top}(p) = A \left(\dfrac{1+\sin{A_2}}{2}\right)\dfrac{h_2(p)^2}{\|h_2(p)^2\|_6}, \end{array}\right\} \end{equation}

where $A$ is approximately the maximum amplitude of $y_{bot}$ and $y_{top}$, approximate because $\|\cdot \|_6$ instead of $\|\cdot \|_\infty$ appears in the denominators. This is done to allow for a smooth dependence of $y_{bot}$ and $y_{top}$ on $h_1$ and $h_2$, which is needed to calculate first derivatives in our quasi-Newton optimization method. An approximate maximum amplitude is sufficient because the optima become insensitive to $A$ above a certain threshold, and we run the simulations for a range of large $A$. With $\|\cdot \|_2$ and $\|\cdot \|_{10}$ in place of $\|\cdot \|_6$ in (5.1), we obtain optimal flow and wall configurations and $Nu$ values that are very similar to those obtained with $\|\cdot \|_6$. The $Nu$ values with $\|\cdot \|_6$ and $\|\cdot \|_{10}$ were somewhat better than with $\|\cdot \|_2$, presumably because constraining the maximum wall deflection (approximated by $\|\cdot \|_6$ and $\|\cdot \|_{10}$) is more important for numerical accuracy than constraining the root mean square ($\|\cdot \|_{2}$). The parameters $A_1$ and $A_2$ allow the approximate amplitudes to take any value between 0 and the maximum, $A$. Terms $h_1(p)$ and $h_2(p)$ are defined as Fourier series in $p$ with wavenumbers ranging from 0 to $M_1$, which results in $2M_1 + 1$ modes (sines and cosines) for each function. The Fourier coefficient vectors are denoted $\boldsymbol {B}_1$ and $\boldsymbol {B}_2$, respectively.

For $L_x$ we use the expression

(5.2)\begin{equation} L_x = \min(5.4 Pe^{{-}0.37},0.5)+1+\sin(L_0), \end{equation}

which, for arbitrary real $L_0$, constrains $L_x$ to be non-negative and bounded below by one of the terms in the min function. The first of these is chosen to be about half the value for the flat-wall case at large $Pe$ from Alben (2023), and the second is used at small $Pe$. The non-zero lower bound is used to avoid spurious optima with very inaccurate $\partial _n T$ that appear at very small $L_x$.

We express the flow $\psi$ using the same modes as in Alben (2023) but in $(p,q)$ coordinates instead of $(p,y)$ coordinates. We also include modes that give a net flow through the channel for greater generality, though their coefficients turn out to be essentially zero in the optima, as was found in preliminary computations in Alben (2023). For $\psi$, we first form $\tilde {\psi }$, a linear combination of functions with p-period 1 that have zero first derivatives in $p$ and $q$ at the walls. The functions are products of Fourier modes in $p$ and linear combinations of Chebyshev polynomials in $q$:

(5.3)\begin{gather} \tilde{\psi}(p,q) = \tilde{A}(3q^2-2q^3) + \sum_{j = 0}^{M}\sum_{k = 1}^{N-3} A_{jk}Q_k(q)\cos(2{\rm \pi} j p) + B_{jk}Q_k(q)\sin(2{\rm \pi} j p), \end{gather}

(5.4)\begin{gather}Q_k(q) \in \langle T_0(2q-1), \ldots, T_{k+4}(2q-1) \rangle. \end{gather}

The first term, with coefficient $\tilde {A}$, gives a net flux through the channel, and the remaining terms modify the flow distribution without changing the net flux. Term $Q_k$ is a linear combination of Chebyshev polynomials of the first kind up to degree $k+4$ that have zero values and first derivatives at $q = 0$ and 1. Its computation is described in appendix B of Alben (2023). The functions $\tilde {\psi }$, $y_{bot}$ and $y_{top}$ can be shifted by the same arbitrary amount in $p$ without changing the solution. To remove this degree of freedom we set $B_{11} = 0$. Since $\partial _p \tilde {\psi } = \partial _q \tilde {\psi } = 0$ at the walls, the first derivatives in the tangential and normal directions ($\partial _s \tilde {\psi }$, $\partial _n \tilde {\psi }$) are zero there as well, so no-slip conditions are obeyed.

We define a grid that is uniform in $p$ with $m$ points, $\{0, 1/m,\ldots, 1-1/m\}$. Typically $m = 256$. We concentrate points near the boundaries in $q$, in case sharp boundary layers appear in the optimal flows as in Alben (2023). This is done by starting with a uniform grid for $\eta \in [0, 1]$ with spacing $1/n$, and mapping to the $q$-grid by

(5.5)\begin{equation} q = \eta -\frac{q_f}{2{\rm \pi} }\sin(2{\rm \pi} \eta), \end{equation}

with $q_f$ a scalar. The $q$-spacing is maximum, ${\approx } (1+q_f)/n$, near $q = 1/2$, and minimum, ${\approx } (1-q_f)/n$, near $q = 0$ and 1. We take $q_f = 0.997$ and $n = 256- 1024$, giving a grid spacing $\Delta q \approx 3 \times 10^{-6}$–$1.2 \times 10^{-5}$ at the boundaries. The derivative operators are discretized with second-order finite differences on these grids.

We obtain $\psi$ from $\tilde {\psi }$ by normalizing it so the flow has the power dissipation rate $Pe^2$ in (2.3). We do this in discrete form, discretizing the derivatives in (2.3) with second-order finite differences and the integrals with the trapezoidal rule in $(p,q)$.

The discretized $\psi$ is written $\mathbf {\varPsi }$, a vector of values at the $m(n-1)$ interior mesh points for $(p,q) \in [0, 1) \times (0,1)$. To form $\mathbf {\varPsi }$, we arrange the modes in (5.3) as columns of an $m(n-1) \times (2M+1)(N-3)-1$ matrix $\boldsymbol{\mathsf{V}}$. Here we take $M = 5m/32$ and $N = 5n/32$, so we limit the modes to those whose oscillations can be resolved by the mesh. We set $\tilde {\mathbf {\varPsi }} = \boldsymbol{\mathsf{V}}\boldsymbol {c}$, a linear combination of the discretized modes with coefficients $\boldsymbol {c}$. To form a $\mathbf {\varPsi }$ with power $Pe^2$, we discretize the integral in (2.3) as a quadratic form $\boldsymbol {a}^T\boldsymbol{\mathsf{M}}\boldsymbol {a}$, where the vector $\boldsymbol {a}$ stands for a discretized $\psi$ in (2.3), and the matrix $\boldsymbol{\mathsf{M}}$ gives the effect of the discretized derivatives and integrals in (2.3). Then

(5.6)\begin{equation} \mathbf{\varPsi} = \frac{Pe\, \boldsymbol{\mathsf{V}}\boldsymbol{c}}{\sqrt{(\boldsymbol{\mathsf{V}}\boldsymbol{c})^T\boldsymbol{\mathsf{M}}\boldsymbol{\mathsf{V}}\boldsymbol{c}}} \end{equation}

has power $Pe^2$, as can be seen by evaluating $\mathbf {\varPsi }^T\boldsymbol{\mathsf{M}}\mathbf {\varPsi }$. Such a $\mathbf {\varPsi }$ automatically gives an incompressible flow by the stream function definition, and automatically satisfies the power constraint.

The various constraints on the walls’ shapes and the flows have been enforced implicitly in the chosen forms of the functions that describe them. We are left with an unconstrained optimization problem to find the values of the parameters or coefficients that appear in these functions, and which may take any real values. Thus we maximize $Nu$ ((2.2), discretized at second order) over $\boldsymbol {c} \in \mathbb {R}^{(2M+1)(N-3)-1}$, $\{\boldsymbol {B_1},\boldsymbol {B_2}\} \in \mathbb {R}^{2M_1+1}$ and $\{A_1, A_2, L_0\} \in \mathbb {R}$. We compute optima over a range of $Pe$, and use various combinations for $M_1$ and $A$ (the maximum wall deformation amplitude, approximately). We find empirically that for a given $M_1$, if $A$ is too large, the algorithm converges to spurious optima that are underresolved. We study this phenomenon using model problems in § 5.1.

We solve the optimization problem using the Broyden–Fletcher–Goldfarb–Shanno method (Martins & Ning 2021), a quasi-Newton method that requires evaluations of the objective function $Nu$ and its gradient with respect to the design parameters, $\{\boldsymbol {c}, \boldsymbol {B_1},\boldsymbol {B_2}, A_1, A_2, L_0\}$. Parameter $Nu$ is computed by forming $\mathbf {\varPsi }$ from the design parameters, then computing $\boldsymbol {u}$ and solving (2.1) for the discretized temperature field $\boldsymbol {T}$ using second-order finite differences. We use a second-order rather than spectral discretization to obtain a sparse matrix in the advection–diffusion equation, allowing for relatively fast solutions.

The gradient can be computed efficiently using the adjoint method (Martins & Ning 2021). The procedure is the same as that described in Alben (2023), but now including the wall shape parameters. In Appendix B we present formulae for the gradient of $Nu$ with respect to the parameters using the adjoint variable.

We initialize with about 100 random sets of $\{\boldsymbol {c}, \boldsymbol {B_1},\boldsymbol {B_2}, A_1, A_2, L_0\}$ with $M_1 = 1- 4$ at various $A$. We run the optimization until the 2-norm of the gradient of $Nu$ is less than $10^{-10}$ or the number of iterations reaches 10 000. The latter case is more common at larger $Pe$, where the gradient is larger in the initial random state, and where convergence to very small gradients is more difficult. Here the gradient norm is often ${\approx }10^{-2}$ after 10 000 iterations. At large $Pe$, the initial gradient norm is typically between $10^3$ and $10^4$, so a gradient norm of $10^{-2}$ corresponds to $10^{-5}$–$10^{-6}$ relative to the initial norm. The discretized advection–diffusion equation becomes increasingly ill-conditioned as $Pe$ increases, causing a decrease of accuracy in $Nu$ and its gradient, and slowing convergence (Kelley 1999). Nonetheless, we are able to obtain accurate $Nu$ values up to $Pe = 10^{7}$, the largest value used with non-flat walls in this study. It is possible to obtain accurate optima at larger $Pe$, but the large-$Pe$ behaviour is already clear at $Pe = 10^{7}$.

5.1. Boundary heat flux resolution

It turns out to be computationally challenging to compute solutions accurately both with large wall-perturbation amplitudes $A$ and with many modes (large $M_1$), or with small $L_x$. This challenge occurs at small $L_x$ even though we use a horizontal coordinate $p$ that is $x$ divided by $L_x$, so the number of grid points per period is fixed. The challenge is illustrated numerically by considering the simplest case of no flow, i.e. Laplace's equation with wavy walls. Figure 2(a–g) shows temperature and heat flux distributions in this and two other simple cases. Figure 2(h–m) shows measures of the numerical errors in these cases with 256 grid points uniformly spaced in the $p$ coordinate.

Figure 2. Examples of the temperature fields and heat flux near wavy walls, and corresponding relative errors. For sinusoidal wavy walls with $A = 0.4$ and no flow, (a) the temperature field, (b) the norm of its gradient and (c) the heat flux density along the bottom wall. (d) Streamlines for a steady flow through the wavy channel with $Pe = 10^4$. (e) Heat flux density along the bottom wall corresponding to the flow in (d). (f) Streamlines for steady convection rolls with $Pe = 10^4$. (g) Heat flux density along the bottom wall corresponding to the flow in (f). (h–m) Relative errors in $\partial T/\partial n$ and $Nu$, computed on a 256-by-257 mesh relative to a 512-by-257 mesh, for various choices of wall perturbation amplitude $A$ and horizontal period $L_x$. Specifically, panels (h,i) plot the maximum relative errors in $\partial T/\partial n$ along the bottom wall and the relative error in its mean, for the case of no flow (corresponding to a–c). We plot the same error quantities for the flow through the channel (panel d) in (j,k), and for the convection rolls (panel f) in panels (l,m). The colourbar at the bottom left shows the error values.

To begin, figure 2(a) shows the pure diffusion temperature field for $A =L_x = 0.4$. The temperature gradient is nearly uniform in the central region, $0 \leq y \leq 1$, and nearly zero above and below this region. Figure 2(b) shows this clearly by plotting values of the temperature gradient norm, which transitions from 1 in the central region to 0 in the remainder of the domain. However, there are small regions at the wall inward extrema (near $y = 0$ and 1) where the gradient norm rises sharply to slightly above 4. Figure 2(c) plots the wall heat flux along the bottom wall, which has sharp peaks at the wall's inward extrema. The wall shape is sinusoidal with two periods per domain period $L_x$, and is easy to resolve with 256 points per $L_x$. But the resulting heat flux peaks are much sharper and more difficult to resolve. For Laplace's equation there are many alternative methods, such as boundary integral methods, that would give better resolution of the boundary heat flux with all the grid points concentrated along the boundary. However, in general we have the steady advection–diffusion equation with arbitrary flow fields in the advection-dominated limit, which requires a fine grid in the interior, not just on the boundary. Figure 2(d,e) shows the streamlines and boundary heat fluxes for a second case, a Poiseuille-like flow that conforms to the boundary, $\psi = c(3q^2-2q^3)$, with $c$ chosen so $Pe = 10^4$. Figure 2(e) shows that the peaks in $\partial _n T$ decrease and become asymmetric, but are still quite sharp. Figure 2(f,g) shows a third case, $\psi = c(3q^2-2q^3)\cos {2{\rm \pi} p}$, with $c$ again chosen so $Pe = 10^4$. The streamlines (figure 2f) show that we have chosen convection rolls aligned with the walls’ extrema, and that between the convection rolls, vertical jets of fluid run from one wall extremum to the other. The configuration is similar to some of the optimal flows and wall shapes that are presented later. The heat flux from the bottom wall (figure 2g) is greatly increased from the previous cases at its maximum, where the downward jet impinges on the lower boundary, as well as at the smaller peak where the upward jet emanates from the lower boundary. All three cases are qualitatively similar in having sharp peaks in $\partial _n T$ that require a fine grid to resolve.

Figure 2(h–m) plots estimates of the errors in $\partial _n T$ in the three cases. For each case, two measures of error are used: the max norm ($\|\cdot \|_\infty$) and the relative error in $Nu$. The error is estimated by comparing $\partial _n T$ at the lower wall, computing $T$ with a 512-by-257 ($\,p,q$) grid and the 256-by-257 grid that omits every other point in the $p$ direction. Figures 2(h) and 2(i) show

(5.7)\begin{equation} \left.\begin{array}{c} \displaystyle \varDelta_{rel}\|\partial_n T\|_{\infty} \equiv \dfrac{\max_p{|\partial_n T_{256}(p,0) - \partial_n T_{512}(p,0)|}}{\max_p|\partial_n T_{512}(p,0)|} ,\\ \displaystyle \varDelta_{rel} Nu \equiv \left|\dfrac{Nu_{256} - Nu_{512}}{Nu_{512}}\right|, \end{array}\right\} \end{equation}

respectively, for the case of no flow (figure 2a–c). Figures 2(j,k) and 2(l,m) show the same quantities for the second and third cases (figures 2d,e and 2f,g, respectively). In each panel in the bottom row, the errors have a similar distribution in $L_x$–$A$ space. Errors are largest at small $L_x$ and large $A$, and are roughly constant along curves $A = L_x^r$ for $r$ slightly larger than 1. The $\varDelta _{rel}\|\partial _n T\|_{\infty }$ errors (figure 2h,j,l) are significantly larger than the corresponding $\varDelta _{rel} Nu$ errors (figure 2i,k,m), particularly in the second and third cases (figures 2j versus 2k and figures 2l versus 2m). This is perhaps not surprising since $\varDelta _{rel}\|\partial _n T\|_{\infty }$ is a measure of local error and is more sensitive to how well the peak of $\partial _n T$ is resolved.

Our optimization algorithm computes $Nu$ over a much wider range of flows than the three cases here, but when we check the accuracy of $Nu$ for our optimal flows, the same general trend holds: for a given grid size and a given level of error, there is a limit as to how large $A$ can be and how small $L_x$ can be. We also find large errors if we exceed relatively small integer values of $M_1$, the number of modes that describe the wall shape. Increasing $M_1$ allows for finer length scales in the wall shape, similarly to decreasing $L_x$, except that the grid spacing is proportional to $L_x$ but does not change with $M_1$. In this study we perform the optimization with $M_1$ ranging from 1 to 4, with a smaller maximum $A$ required at larger $M_1$. The limit on $A$ can be (very roughly) approximated by a power law of the form $M_1^R$ with $R \approx -1.7$. Fortunately, we find in our optimal flow computations that the optimal flows give significantly lower relative errors at a given ($A$, $L_x$) pair than the values shown in figure 2(h–m). The reason is not entirely clear, but the lower errors are seen across a wide range of different optimal flow and wall configurations. Specifically, at $Pe = 10^2, 10^3, \ldots, 10^7$, we take the 2–4 flow and wall shape optima with largest $Nu$ at each $M_1 = 1- 4$, yielding 10–12 optima in total at each $Pe$. The optima are computed with $m = 256$ and $n$ ranging from 256 to 1024, with larger $n$ at larger $Pe$. For each case, we compute the error quantities (5.7) by doubling $m$ and keeping $n$ fixed. We find that $\varDelta _{rel} Nu \leq 0.012$ across all 78 optima, and $\leq$0.005 in 74 out of 78 cases. Also, $\varDelta _{rel}\|\partial _n T\|_{\infty } \leq 0.13$ in all 78 cases with a median value of 0.044. The maximum error occurs at a sharp peak of $\partial _n T$ and is much smaller elsewhere, which is why the error in $Nu$ is much smaller. Since $Nu$ is our main focus of interest, $m = 256$ is a reasonable value to use for the optimization, since much larger $m$ would greatly slow the computations.

We also study the analogous error quantities when $n$ is doubled and $m$ is fixed. The maximum of $\varDelta _{rel} Nu$ is somewhat larger than previously, 0.033 over the 78 cases, though it is below 0.0076 in 74 out of 78 cases. By contrast, the maximum $\varDelta _{rel}\|\partial _n T\|_{\infty }$ is much smaller than previously, 0.041 over the 78 cases, with a median value of 0.0021.

Having described the computational framework, we now use it to validate the decoupled approximation (4.13) that was used in § 4 to explain the transition to optima with non-flat walls at a critical $Pe > 0$.

5.2. Testing the accuracy of the decoupled approximation

In § 4 we derived an expansion for $Nu$ for small flow and wall perturbations in which the leading-order effects of the perturbations are decoupled. In figure 3 we show that the expansion (4.13) is accurate for a large ensemble of flows and wall shapes over ranges of small $Pe$ and $A$. The expansion should hold for any flows and wall shapes – not just optima – as long as $Pe$ and $A$ are sufficiently small.

Figure 3. Validation of (4.13) which shows the decoupled effect of small wall and flow perturbations on Nusselt number at leading order. (a) Example of a wall shape and flow streamlines using randomly generated coefficients. (b) Streamlines of the optimal flow for $Pe = 10^{1.5}$ (with flat walls), stretched vertically to correspond to a sinusoidal wall perturbation. (c,d) Contour maps of the relative errors (log base 10) in the linear decoupled approximation to the Nusselt number for the streamlines and wall shapes in (a,b), but with a range of flow and wall perturbation amplitude $Pe$ and $A$, respectively. The relative error is $\log _{10} (|Nu-(Nu_A+Nu_B-Nu_0)|/|Nu-Nu_0|)$, with the quantities as defined in the main text. (e) Contour maps of the maximum of the relative error over the values in (c,d), and an ensemble of 50 other cases described in the main text.

To test this hypothesis, we take $m = n = 256$, set $M_1 = 4$ giving 9 modes for each wall shape and use 15 flow modes (products of the first three modes in the $p$ direction with the first five modes in the $q$ direction, corresponding to $M = 1$ and $N = 8$). We generate 50 random sets of coefficients for all the modes to obtain an ensemble of 50 sets of random wall shapes and flows (one example is shown in figure 3a, with $A = 10^{-0.5}$), and add two more non-random choices. These are the optimal flows at $Pe = 10^{1.5}$ and $10^2$ with flat walls, i.e. convection rolls $\psi (p,y)$ on a rectangular domain $0 \leq p, y \leq 1$. We set $y_{bot}(p)$ and $y_{top}(p)$ to sinusoidal walls of varying $A$, and substitute $q$ for $y$ in $\psi$, essentially stretching the rectangular convection rolls to fit the wavy walls. Figure 3(b) shows an example – the optimal flow at $Pe = 10^{1.5}$ with $A = 10^{-0.5}$.

For the resulting set of 52 flow and wall-shape combinations, we vary $Pe$ from $10^{-1.5}$ to $10^1$, vary $A$ from $10^{-4}$ to $10^0$ and compute $Nu$ in each case. In order to validate (4.13) we need to compute $Nu_{1A}$ and $Nu_{1B}$ and compare the sum with $Nu - 1$. We do this in a way that highlights the decoupled effect of the flow and wall perturbations. Parameter $Nu$ with both the flow and wall perturbations should be approximately $Nu_0 + Nu_{1A} + Nu_{1B}$; $Nu$ with the same wall perturbation but no flow – call it $Nu_A$ – should be approximately $Nu_0 + Nu_{1A}$; and $Nu$ with the same flow perturbation but no wall perturbation – call it $Nu_B$ – should be approximately $Nu_0 + Nu_{1B}$. Summing the last two expressions, $Nu_A+Nu_B-Nu_0$ should be approximately $Nu$. In figure 3(c–e) we plot the absolute value of the ratio of the difference of $Nu$ and $Nu_A+Nu_B-Nu_0$ – which should be the omitted higher-order terms in (4.13) – relative to $Nu-Nu_0$, which approximates the first-order terms. If the ratio is small, the first-order terms that arise from decoupling the wall and flow perturbations are indeed dominant over the remaining terms. Figure 3(c) gives the relative error (log base 10) for the flow and wall configuration in figure 3(a), figure 3(d) gives the relative error for figure 3(b) and figure 3(e) gives the maximum of the relative errors in figure 3(c,d) and the 50 other flow and wall shape combinations in the ensemble of 52 described earlier. Figure 3(c–e) shows that the error in the decoupled approximation is indeed small when $A$ and $Pe$ are small. Figure 3(e) shows that the relative errors become significant in some cases when $Pe$ reaches $10^{0.3}$. The transition to optima with wavy walls occurs at higher $Pe$, between $10^{1}$ and $10^{2}$.

It is possible that the largest relative errors in figure 3(e) are due to flows that are unlike the optimal flows and wall shapes near the transition. We will see in figure 4 that near the transition the optima are relatively smooth, as in figure 3(b), and compared with the less smooth case in figure 3(a), we can see smaller relative errors in the decoupled approximation for the smoother case (contours of a given value in figure 3d are shifted upwards and rightwards relative to those in figure 3c). Therefore, for the optimal flows, perhaps the decoupled approximation remains accurate up to higher $Pe$ than for some of the random flows studied here.

Figure 4. Optimal flows and wall shapes, temperature fields and heat flux densities at $Pe = 10^2$. Eight cases are shown, four in the top three rows and four in the bottom three rows. Each case is labelled A–H at the top, followed by the Nusselt number value. Below the label are three rows with a contour plot of the streamlines (top row), the temperature field colour plot (middle row, with colourbar at left) and a plot of the heat flux distribution versus the horizontal coordinate along the lower wall (bottom row). The numbers of modes describing the walls’ shapes ($M_1$) are (A, B) 1, (C, D) 2, (E–G) 3; H shows the flat-wall optimum for comparison.

6. Moderate-to-large convection ($Pe \gtrsim 1$)

The transition to optima with non-flat walls at a critical non-zero $Pe$ value, predicted theoretically in § 4, is consistent with our optimization calculation. We find only optima with flat walls for $Pe \leq 10^1$, but at $Pe \geq 10^2$, we find optima with non-flat walls that outperform those with flat walls – by an increasing margin as $Pe$ increases, as we will show. Using plots of $Nu$ versus $Pe$ (like figure 6, shown later) for the best local optima with flat and non-flat walls, we find the two curves meet at $Pe\approx 10^{1.6}$, which approximates the transition $Pe$ value where non-flat walls become optimal. We find non-flat optima at $Pe = 10^{1.6}$ but not at $10^{1.2}$ and $10^{1.4}$, so the transition value seems to lie between $10^{1.4}$ and $10^{1.6}$.

For $Pe = 10^2$, the top-performing optima are shown in figure 4 (cases A–G), and the flat-wall case, also a local optimum, is shown as case H. For each optimum, three panels are shown. The top panels, labelled A–H and followed by the $Nu$ values, show the streamlines and wall shapes. The middle panels show the temperature fields and the bottom panels show the heat flux per unit horizontal width along the bottom walls. These optima have different numbers of modes ($M_1$) for the walls’ shapes – (A, B) 1, (C, D) 2, (E–G) 3 – and the $Nu$ values, shown at the top, do not correlate strongly with $M_1$.

The convection rolls with non-flat walls (A–G) roughly resemble those with flat walls (H), but are stretched by various amounts vertically, either upward or downward, or both. The non-flat cases have slightly larger $L_x$ values, 1.9–2.2, versus 1.7 for the flat case.

The temperature fields (middle panels) are generally similar in all cases, with one upward hot plume and one downward cold plume. Cases A–C and E have large wall deformations but near the peaks of the deformations there is almost no flow, and the temperature is almost constant. These regions could be extended farther downward (for the bottom-wall peaks) or upward (for the top-wall peaks) without significantly altering the flow and temperature elsewhere, the net power dissipation or the net heat transfer. Therefore, there is some insensitivity to wall perturbations beyond a certain magnitude.

The distributions of heat flux per horizontal width from the bottom wall (bottom panels) clearly show the effect of waviness in the bottom wall (A, B, D–G) compared with the cases with flat bottom walls (C, H). In the flat cases, there is a peak in heat flux where the cold plume reaches the bottom wall, a trough where the hot plume emerges from the wall and a monotonic change between them. For the non-flat walls, there is one large peak where the downward jet impinges on an upward peak in the bottom wall. There is usually a smaller ‘dimpled’ peak (in A, B, D and F) where the upward jet emanates from another upward peak on the bottom wall. Between these peaks, the rest of the heat flux distribution is essentially zero, in the troughs along the bottom wall. These are in the aforementioned ‘dead zones’ of no flow and nearly constant temperature corresponding to nearly zero heat flux density. In each of these cases, the upward peaks of the bottom wall usually align with the upward and downward jets, i.e. the streamlines that connect stagnation points on the top and bottom walls. Although the heat flux distributions are quite different from those of case H, the net increases in $Nu$ are modest, about 1 %–5 % at this $Pe$. Cases E and G show similar but less symmetric changes in the extrema in the heat flux distributions, because the upward peaks of the bottom wall are asymmetric, sloped on one side (where heat transfer density is larger) and flattened on the other (where heat transfer density is smaller).

Next, we compare the best optima with non-flat and with flat walls across a range of $Pe$ extending up to five orders of magnitude larger than that of figure 4. Each column of figure 5 compares these two optima at $Pe = 10^3$, $10^4$, $10^5$, $10^6$ and $10^7$ (in order from left to right). In each column, the top six panels arrayed vertically show the streamlines (first and second rows, $Nu$ values at the top), temperature fields (third and fourth rows, colourbar at left) and flow speed distributions (fifth and sixth rows, colourbar at left), with the non-flat optima above the flat optima in each case. The seventh row shows the bottom-wall heat flux distributions (flat-wall case in red, non-flat-wall case in blue).

Figure 5. Comparisons of the best optima with non-flat walls and flat walls. Each column compares the two cases for $Pe = 10^3$, $10^4$, $10^5$, $10^6$ and $10^7$ (in order from left to right). In each column, the top two rows show streamline contour plots for the optima with non-flat and flat walls, respectively, each labelled at the top by the corresponding $Nu$ value. The third and fourth rows show the temperature fields (colourbar at left). The fifth and sixth rows show the flow speed distributions. The seventh row compares the heat flux distributions along the bottom wall for the optima with non-flat walls (blue) and flat walls (red).

Figure 6. Values of $Nu$ and $L_x$ for optima with non-flat walls (black plus signs and circles, respectively), with flat walls (red crosses and circles, respectively) and for two special flow and wall configurations A and B (green and blue symbols, respectively), for $Pe$ from $10^2$ to $10^7$.

First, we note from the $Nu$ values at the top of the streamline plots (top two rows) that $Nu$ is larger for the non-flat walls, 37 % larger at $Pe = 10^3$ and a factor of 3.8 larger at $Pe = 10^7$. The streamline distributions show a series of convection rolls for both flat and non-flat walls. Unlike in figure 4, the $x$ axes are dilated relative to the $y$ axes. This makes the flow structures visible at large $Pe$ and in all cases with non-flat walls, whose deflections are very large. The streamlines are somewhat concentrated away from the walls in the non-flat cases, leaving ‘dead zones’ near the walls. The non-flat cases’ streamlines show a few different types of configurations. The first and third columns have convection rolls positioned horizontally between the walls’ inward extrema, though the two rolls are much more symmetric in the third column. In the second column, the convection rolls touch the inward extrema on one side only, and lie below the outward extrema on the other side. The fourth and fifth columns’ rolls are similar but skewed, so they meet the inward extrema at an oblique angle. We will show later that similar versions of these different optima and many more seem to occur in all $Pe$ with wavy walls. That is, the same types of optimal flow and wall configurations occur at all $Pe \gtrsim 10^3$. The main difference is simply that the horizontal period $L_x$ decreases with $Pe$ as a power law (to be shown later). Both walls are wavy at these higher $Pe$, unlike for most cases at $10^2$ in figure 4.

The corresponding flat-wall optima, shown in the second row, have predominantly rectangular convection rolls with $L_x$ that also decreases as a power law. The $L_x$ values in the flat cases are about a factor of 4 smaller than in the non-flat cases at the largest $Pe$, but the difference is closer to a factor of 2 when one considers that there are two pairs of convection rolls in the top row of the fourth and fifth columns. There is a strong skewness in the flat walls’ convection rolls in the last two columns, but this is exaggerated by the $x$ dilation, and does not seem to influence heat transfer much (Alben 2023). Actually, there is a key difference in the flat-wall optima that is difficult to see in the streamline plots. They have a thin boundary layer whose thickness scales as $Pe$ to a power, close to that of $L_x$. The boundary layer switches from a horizontal flow along the wall to branching flows in the last two columns, too thin to be seen clearly. No similar boundary layer is seen in the wavy-wall optima.

The presence or absence of boundary layers is somewhat apparent in the temperature fields, shown in the third and fourth rows. In the third row, as in figure 4, there are hot and cold temperature plumes aligned with the downward and upward jets that run between the walls. The full range of temperatures from 0 to 1 can be seen clearly in all the columns of the third row, across $Pe$. The fourth row also has hot and cold temperature plumes, but moving leftward toward larger $Pe$, the temperature values are more concentrated in the middle range, near 0.5. Temperatures close to 0 and 1, the wall temperatures, are confined to ever thinner boundary layers. The flow speed fields, in the fifth and sixth rows, show the boundary layer phenomenon as well. In the fifth row, with wavy walls, the full range of flow speeds occurs throughout the domain and is equally visible at small and large $Pe$. In the sixth row, by contrast, the largest flow speeds (yellow) are less visible at large $Pe$ because they are more confined to the boundary layer. Evidence for boundary layers is not particularly clear in the seventh row, the bottom-wall heat flux distributions, but clear differences are seen between the wavy-wall and flat-wall distributions (blue and red, respectively). Logarithmic scales are used instead of the linear scales in figure 4, but otherwise the flat-wall distributions are similar, with intervals of monotonic change, more numerous for the branching flows in the fourth and fifth columns. The wavy-wall distributions show peaks near the convection rolls and decay sharply in the ‘dead zones’ where temperatures are nearly constant.

In figure 6 we plot $Nu$ and $L_x$ (red crosses and circles) for the best flat-wall optima at $Pe = 10^2, 10^3, \ldots, 10^6, 10^{6.5}, 10^7$, and at each $Pe$, $Nu$ and $L_x$ values (black plus signs and circles) for an ensemble of 8–12 non-flat optima that include the 2–4 best optima at each $M_1$ from 1 to 4, including all the cases in figures 4 and 5. The fit lines showing $Nu\sim Pe^{0.575}$ and $L_x \sim Pe^{-0.36}$ are shown in red next to the corresponding flat-wall data. There are subtle but noticeable changes in the positions of the data relative to the fit lines when $Pe$ exceeds the branching transition value for the flat-wall optima, between $Pe = 10^5$ and $10^{5.25}$ (Alben 2023). At lower $Pe$, the flat-wall $Nu$ are better fit by $Pe^{0.54}$, while the $L_x$ values (red circles) lie slightly below the red fit line versus slightly above at larger $Pe$. The $Nu$ values for the non-flat optima (black plus signs) seem to follow a single trend, $Pe^{2/3}$, more consistently. This is the scaling that was found for computations of 3-D optimal flows with flat walls by Motoki et al. (2018a). This is also an a priori upper bound that was proved using the background method for the 2-D flat-wall case (Souza 2016), and for the same class of non-flat walls considered here by Goluskin & Doering (2016). In Appendix C we compare the $Nu$ values in figure 6 with those from recent simulations of Rayleigh–Bénard convection with wavy and flat walls. Perhaps not surprisingly, the optimal flows with flat and wavy walls give significantly larger $Nu$ than both steady and unsteady flows from natural convection at the same $Pe$.

The $L_x$ values for the non-flat optima are shown by black circles in figure 6, and they show more scatter than the $Nu$ values (the black plus signs). One reason is that for these optima the horizontal periods may contain one or two pairs of convection rolls. Differences in flow and wall configurations like in the top row of figure 5 also account for a significant amount of the scatter in $L_x$ values, as can be seen subsequently in figure 7, by comparing $L_x$ for flows at the same $Pe$. The $L_x$ values in the flat cases scale approximately as $Pe^{-0.36}$, while in the non-flat cases, the bottom envelope of the black circles scales approximately as $Pe^{-0.35}$. These scalings are much closer than those for the $Nu$ values.

Figure 7. Typical flow and wall configurations of optima. Each purple box (labelled A–H at the top) shows a pair of streamline contour plots for optima with similar flow and wall configurations, one at high $Pe$ (left) and one at low $Pe$ (right). At the top of each contour plot is a label showing ($Pe$ value): ($Nu$ value). Following the eight pairs are eight individual examples of other optimal flows and wall shapes (the two rightmost plots in the third row and all six of the plots in the fourth row).

By extrapolating the plots of $Nu$ versus $Pe$ for flat and non-flat walls to smaller $Pe$, we find the two curves meet at $Pe \approx 10^{1.6}$, which approximates the transition $Pe$ value where non-flat walls become optimal.

In figure 6 we also plot $Nu$ and $L_x$ values (green and blue, labelled A and B in the legend) for two specific flow configurations to be discussed later.

When discussing the best optima with wavy walls that we found, shown in figure 5, we suggested that all the optima are different members of a single class that occurs at all $Pe \gtrsim 10^3$, and are simply scaled by different $L_x$ at each $Pe$. In figure 7 we provide evidence for this claim. Starting from the top, we have eight purple boxes, labelled A–H, showing a pair of optima. The first occurs at a large $Pe$ ($10^7$, labelled at the top) and the second at a much smaller $Pe$, $10^3$–$10^5$. The two are chosen from among the four top optima at each $M_1$ ranging from 1 to 4 (except for the second member of pair H which is eighth best at $Pe = 10^5$ and $M_1 = 3$). In A–D, $M_1 = 1$; in E and F, $M_1 = 2$; in G and H, $M_1 = 3$ except for the first member of pair H, which has $M_1 = 4$. The remaining flows, following the eight pairs, are chosen simply for their differences from the preceding flows and from each other, to show the diversity of flows and wall configurations that occurs among the optima.

Pair A are optima that show essentially the same pattern of streamlines and wall shapes (but with reversed flow directions), even though $Pe$ differs by a factor of $10^3$, $Nu$ by a factor of about $10^2$ and $L_x$ by about a factor of 10. Pair B have streamlines that are more confined vertically than pair A, and with more asymmetric walls, but very similar $Nu$ values, showing that very different flows can achieve heat transfer close to the top values we have found. Pair C have a more complicated wall shape and less symmetric vortices than pairs A and B. The second member of pair C is approximately obtained by reflecting the first member about a horizontal line and reversing the flow direction, and increasing $L_x$ by a factor of about 20. Pair D is an almost identical pair of wall shapes and flows, except for the difference in $L_x$. In E and F, $M_1$ is increased to 2, allowing for somewhat more complicated wall shapes and vortices. Pairs G and H have wall shapes that commonly occur at $M_1 = 3$ and 4, with peaks resembling indented hats. Pairs A–H are just a small sample of the similar flows and wall shapes that occur across $10^3 \leq Pe \leq 10^7$. Examination of a larger set of 100 optima at each $Pe = 10^2, 10^3, \ldots, 10^6, 10^{6.5}$ and $10^7$ shows that most flow and wall configurations recur throughout this $Pe$ range with little change apart from $L_x$.

The last eight flows and wall shapes in figure 7 are chosen to illustrate optimal configurations that are clearly different from those that have already been shown. In order from left to right and top to bottom, these have $M_1 = 1$, 1, 2, 2, 2, 3, 2 and 2, respectively. The convection rolls appear somewhat oblique and asymmetric in some cases, but when the horizontal dilation of the figures is taken into account – i.e. that $L_x \ll 1$ in most cases – the convection rolls are seen to be very thin horizontally and elongated vertically. The walls have various positions but in all cases the sides of the convection rolls run alongside the walls vertically for $O(1)$ distances at all $Pe$, allowing for heat transfer over a large surface length.

Having presented a variety of optimal flow and wall configurations, we now attempt to understand the basic scalings of the optima with $Pe$, i.e. $Nu \sim Pe^{2/3}$ (the black plus signs in figure 6) and $L_x \sim Pe^{-0.35}$ (the bottom envelope of the black circles in figure 6). We consider a simple model flow and wall configuration, the same one we considered in figure 3(b), but with $A$ increased to 1. This configuration, shown in figure 8(a), is a model for some of the optima we have seen, e.g. at the top of the first and third columns in figure 5, and pairs B, G and H in figure 7. The ‘dead zones’ of little flow and nearly constant temperature seen in some of those cases can essentially be treated as though the wall were at the inward boundary of the dead zone. There is little viscous dissipation in the dead zone and the temperature at its inward boundary is almost the same as at the wall. To the right of panel (a) in figure 8, a four-by-four array of temperature fields is shown. These result from the configuration in figure 8(a) but with the stream function scaled to have four different $Pe$. At each $Pe$, the flow and wall shape are scaled horizontally to have four different $L_x$. Thus, from left to right, the four columns correspond to $Pe = 10^5$, $10^6$, $10^7$ and $10^8$, respectively. Within each column, four $L_x$ values are used, increasing from top to bottom (listed at the right-hand end of the $x$ axis in each case). The $Nu$ values are shown at the top of each temperature field, and the maximum values occur in the third row.

Figure 8. For the flow and wall configuration in (a), variations of temperature fields and heat flux distributions with $L_x$ are shown at four large $Pe$ values. The same streamline pattern, shown in the upper left panel, is used in all cases. In the four-by-four grid of temperature fields, $Pe$ increases from left to right with the values $10^5$, $10^6$, $10^7$ and $10^8$. Within each column, $Pe$ is fixed and $L_x$ increases from top to bottom, with the values given at the right-hand ends of the $x$ axes. The column at the far right gives the heat flux distribution along the bottom wall for the rightmost temperature distributions, $Pe = 10^8$. (b) An alternative flow configuration, described in the main text.

Each column shows a transition between two types of temperature fields as $L_x$ increases. First, we note that all the temperature fields have a basic structure similar to those already described in figures 4 and 5. In each case there is an upward plume of hot fluid at the left and a downward plume of cold fluid at the right, running between the inward extrema of the two walls. In the top row, the temperature fields resemble those with diffusion only (e.g. figure 2a) in two ways. First, the temperature is almost constant in the regions near the walls’ outward extrema. Second, at each $x$ location the temperature field decreases almost monotonically from the bottom wall to the top wall. The temperature gradient is close to linear in the middle part of the flow, $0 \leq y \leq 1$, particularly at the largest $Pe$ (the upper right corner, $Nu = 3750$). For the rightmost temperature distributions ($Pe = 10^8$), the column at the far right shows the distributions of heat flux at the bottom wall. Along most of the wall, the temperature gradient is nearly zero, but is large close to the walls’ inward extrema, so $Nu$ is much greater than the value with diffusion only, 1. However, with further decreases in $L_x$, the temperature field tends to the diffusion-only case and $Nu \to 1$.

In the bottom row at the largest $L_x$, the temperature is more well-mixed in the middle part of the flow, with temperatures close to 1/2. This well-mixed fluid penetrates close to the walls even at their outward extrema, so there is significant heat flux along most of the bottom wall, as shown by the plot adjacent to the bottom-right temperature field. Whereas the top row resembles the diffusion-dominated limit, the bottom row is a more convection-dominated case. As the flow speed increases further, the temperature field in most of the domain becomes closer and closer to 1/2, because the flow on streamlines adjacent to the walls exchanges only a small amount of heat with the hot and cold walls during the short time it passes along either one. The hot upward and cold downward plumes shrink to very thin jets concentrated along the streamlines that emanate from the stagnation points on the upper and lower walls. This extreme is suboptimal for efficiency in two ways. First, viscous power is expended for the flow in the middle of the domain, but it transfers little heat because the temperature is almost uniform there, so the temperature gradient is small. Second, the upward and downward plumes that do transport heat are more widely separated in the bottom row than in the rows above. Therefore there is less heat exchange per unit horizontal width. By comparison, the top row has more upward and downward plumes per unit horizontal width, but the small values of $L_x$ there mean that the flow needs to be much slower in the top row than in the bottom row to have the same viscous dissipation rate. Hence a more viscous-dominated temperature field occurs. Optimal heat transfer, in the third row, seems to occur at the interface of the viscous-dominated and convection-dominated regimes. The value of $L_x$ is large enough, and thus the flow speed is large enough, to allow cold fluid emanating from the top wall to remain relatively cold as it approaches and passes close to the bottom wall, giving a large temperature gradient there (likewise for the hot fluid emanating from the bottom wall, when it nears the top wall). But $L_x$ is also small enough to have many plumes per unit width.

We consider again the heat flux distribution plots in the rightmost column. The largest peak values occur in the second row, about twice those of the smallest peak values, in the fourth row. However, both plots have about the same average, which is $Nu$ ($1.38\times 10^{4}$ versus $1.36\times 10^{4}$, shown above the temperature fields to the left). The larger peaks in the second row have much smaller widths than in the fourth row. The third row has almost the peak values of the second row combined with almost the widths of the fourth row, so it is better than both. The ratio in peak values between the third and fourth row is roughly the ratio of $\textrm {d}s/{\textrm {d}\kern 0.06em x}$ values (the inverse of the ratio of $L_x$ values), so the greater heat transfer in the third row is due to having larger $\textrm {d}s/{\textrm {d}\kern 0.06em x}$ with about the same $\textrm {d}T/\textrm {d}n$, in other words, more surface area of the bottom wall per unit horizontal width.

For each $Pe = 10^3, 10^4,\ldots, 10^6$, $10^{6.5}$ and $10^7$ we find the $L_x$ that maximizes $Nu$ for the flow in figure 8(a). The $L_x$ values are plotted as green asterisks in figure 6, and the $Nu$ values as green triangles. We do the same for the flow in figure 8(b), which has the convection rolls shifted by a quarter period horizontally. This configuration, with the jets running between the farthest points on the two walls, is more different than flow A from the optimal flows we have seen so far. The $Nu$ values, shown as blue squares in figure 6, are well below those of flow A, but both have the same scaling, ${\sim }Pe^{2/3}$ (the slopes are within 0.5 % of 2/3 for each of the line segments connecting the data points at the four largest $Pe$). The $L_x$ data fit a $Pe^{-1/3}$ scaling nearly as well for both flow configurations. Although flows A and B and the more general optima all have $Nu \sim Pe^{2/3}$, there is a significant difference in the prefactors. The highest $Nu$ for the more general optima is about 30 % higher than that of flow A, which is in turn about 30 % higher than that of flow B. Such differences are perhaps not surprising given the differences in the streamlines of the different flows. Furthermore, all these flows approximate $L_x \sim$ $Pe^{-1/3}$.

To explain these scalings, we note that flows A and B have the form

(6.1)\begin{equation} \psi = c\tilde{\psi}(x/L_x, y) = c\tilde{\psi}(p, y). \end{equation}

That is, they each have a single pattern of streamlines given by $\tilde {\psi }$ that is scaled horizontally by $L_x$, and with $c$ chosen to give the power dissipation rate $Pe^2$. The optima from the general optimization problem also seem to have this form, as seen in pairs A–H of figure 7 for example. For a flow of form (6.1) the power dissipation constraint (2.3) is

(6.2)\begin{equation} c^2 \int_0^1 \int_{y_{bot}(p)}^{y_{top}(p)} \left(\frac{1}{L_x^2} \partial_{pp} \tilde{\psi} + \partial_{yy} \tilde{\psi} \right)^2 {\rm d}p\, {{\rm d}\kern0.05em y} = Pe^2. \end{equation}

For flows A and B, $y_{bot}$, $y_{top}$ and $\tilde {\psi }$ are fixed with respect to $Pe$ and $L_x$, and the same holds approximately for pairs A–H of figure 7. The size of $L_x$ determines whether one term in parentheses in (6.2) is dominant, and we use this to determine the scaling of $c$ in terms of $Pe$ and $L_x$:

(6.3a,b)\begin{equation} L_x \ll 1 \rightarrow c \sim Pe\, L_x^2,\quad L_x \gtrsim 1 \rightarrow c \sim Pe. \end{equation}

We can estimate $Nu$, the rate of heat transfer per horizontal width, in the case of a thin thermal boundary layer on the walls, e.g. in the last row on temperature fields in figure 8. In those cases, the walls are nearly vertical except at their extrema. Near the wall, the flow velocity is approximately vertical, and the vertical component is

(6.4)\begin{align} v ={-}\partial_x \psi ={-}\frac{c}{L_x} \partial_{p} \tilde{\psi} &={-}\frac{c}{L_x} (\partial_{p} \tilde{\psi}|_{p = p_{wall}} + (\partial_{pp} \tilde{\psi}|_{p = p_{wall}})p +O(p^2)) \end{align}

(6.5)\begin{align} &={-}\frac{c}{L_x} \left( (\partial_{pp} \tilde{\psi}|_{p = p_{wall}})\frac{x}{L_x} +O\left(\frac{x}{L_x}\right)^2\right). \end{align}

In (6.4) we have Taylor-expanded $\partial _{p} \tilde {\psi }$ about a point on the wall. By the no-slip condition, $\partial _{p} \tilde {\psi }|_{p = p_{wall}} = 0$, so removing this term and writing $p = x/L_x$ gives (6.5). In the advection–diffusion equation (2.1) we can follow the classical boundary layer analysis (Lévêque 1928; Bergman et al. 2011; Shah & London 2014) and neglect $\partial _{yy}T$ relative to $\partial _{xx}T$ in the thin thermal boundary layer near a vertical wall, because $\partial _{xx}T \sim 1/\delta ^2$, where $\delta$ is the thickness of the boundary layer over which the temperature changes by $O(1)$ (from 0 or 1 at the wall to about 1/2 in this problem), while $\partial _{yy}T \sim O(1)$, since the temperature changes by $O(1)$ over an $O(1)$ distance along the wall. Since $u \ll v$ near the wall, $u\partial _x T$ and $v\partial _y T$ are of the same order (Lévêque 1928; Bergman et al. 2011; Shah & London 2014). Using the Taylor expansion (6.5) and $\partial _{y}T \sim O(1)$,

(6.6)\begin{equation} \frac{c}{L_x^2}x \sim v\partial_y T \sim \partial_{xx} T \sim \frac{1}{\delta^2}. \end{equation}

In (6.6) we have $x \sim \delta$ because we are in the thermal boundary layer, so $\delta \sim c^{-1/3} L_x^{2/3}$. The two estimates for $c$ in (6.3) then give two estimates for $\delta$:

(6.7a,b)\begin{equation} L_x \ll 1 \rightarrow \delta \sim Pe^{{-}1/3},\quad L_x \gtrsim 1 \rightarrow \delta \sim Pe^{{-}1/3}\, L_x^{2/3}. \end{equation}

One could perhaps obtain more detailed formulae including prefactors using the similarity solution of Lévêque (1928), but we focus only on the scaling with $Pe$ here. Parameter $Nu$ is given by (2.2), and we estimate it assuming the total length of approximately vertical walls per horizontal length $L_x$ is $O(1)$, so

(6.8)\begin{gather} Nu \sim \frac{1}{L_x} \partial_x T \sim \frac{1}{L_x} \frac{1}{\delta}, \end{gather}

(6.9a,b)\begin{gather}L_x \ll 1 \rightarrow Nu \sim Pe^{1/3}\,L_x^{{-}1},\quad L_x \gtrsim 1 \rightarrow Nu \sim Pe^{1/3}\, L_x^{{-}5/3}. \end{gather}

In both estimates in (6.9), smaller $L_x$ gives larger $Nu$. Therefore we do not consider the possibility that $L_x$ is sufficiently large that the walls are no longer approximately vertical because this case is clearly suboptimal. The assumption of a thin thermal boundary layer constrains how small $L_x$ can be. If $L_x \ll \delta$, the thermal boundary layer is much larger than the horizontal period, i.e. the spacing between approximately vertical segments on the same wall. This is the case in the top row of figure 8, with the smallest $L_x$. In the bottom row, the temperature varies from 0 or 1 to 1/2 over the thermal boundary layer, but in the top row, the boundary layer is so large relative to the spacing between adjacent vertical wall segments that the temperature does not change much from 0 or 1 in the region between the two wall segments. Then there is little heat transfer from the vertical wall segments; it is confined to small regions near the inward wall extrema, also seen in figure 2(c). This is the diffusion-dominated limit, in which $Nu$ is much smaller, $O(1)$ at large $Pe$. The upper boundary of this regime, $L_x \sim \delta$, is the smallest that $L_x$ can be. Thus by (6.7) and (6.9) we have

(6.10a,b)\begin{equation} L_x \sim \delta \sim Pe^{{-}1/3},\quad Nu \sim Pe^{2/3}. \end{equation}

Before concluding, we discuss the observation that the hot and cold wall separation distance is always close to the minimum, 1, in the computed optima.

7. Optimality of minimal wall separation

We mentioned in § 2 that the separation distance between the hot and cold walls is constrained to be $\geq$1, but it turns out to be 1 (to within about $10^{-4}$) for all the computed optima. We now use a scaling argument to explain why this is so. Let $\{\psi _1(x,y); L_{x,1}; y_{bot,1}(x); y_{top, 1}(x)\}$ maximize

(7.1)\begin{equation} Nu = \frac{1}{L_x}\int_{C = \{(x, y_{bot}(x))\}} \partial_n T \,{\rm d}s, \end{equation}

with power

(7.2)\begin{equation} Pe^2 \equiv \frac{1}{L_x}\int_0^{L_x} \int_{y_{bot}(x)}^{y_{top}(x)} (\nabla^2 \psi)^2\, {{\rm d}\kern0.05em y} \, {{\rm d}\kern 0.06em x}, \end{equation}

such that

(7.3)\begin{equation} -\!\nabla^\perp \psi \boldsymbol{\cdot} \boldsymbol{\nabla} T - \nabla^2 T = 0, \end{equation}

with temperatures 1 and 0 on the bottom and top walls, respectively, and with $\max y_{bot,1}(x) = 0$ and $\min y_{top,1}(x) = 1$. That is, this is an optimum over all periodic wall shapes with separation distance exactly 1. Write the resulting values of $Nu$ and the horizontal period as $Nu_1(Pe)$, $L_{x,1}(Pe)$.

Now change the separation distance from 1 to $Y$, any positive value. Let $\{\psi _Y(x,y); L_{x,Y}; y_{bot,Y}(x); y_{top, Y}(x)\}$ maximize $Nu$ with the same power as before, $Pe^2$, and with $\max y_{bot,Y}(x) = 0$, but $\min y_{top,Y}(x) = Y$, so this is the optimum over walls whose separation distance is exactly $Y$ now. If we divide all distances by $Y$, we see that the resulting $Nu$ can be written

(7.4)\begin{equation} Nu_Y = \frac{1}{Y}\frac{Y}{L_{x,Y}}\int_{C = \{x/Y, y_{bot}(x/Y)/Y\}} \partial_n T \,{\rm d}s, \end{equation}

and the power can be written

(7.5)\begin{equation} Pe^2 \equiv \frac{1}{Y^3}\frac{Y}{L_{x,Y}}\int_0^{L_{x,Y}/Y} \int_{y_{bot}(x/Y)/Y}^{y_{top}(x/Y)/Y} (\nabla_{(x/Y),(y/Y)}^2 )^2 \,{\rm d}(y/Y) \, {\rm d}(x/Y), \end{equation}

and (7.3) still holds in the rescaled coordinates. Since $Y$ is fixed, an optimum of problem (7.4)–(7.5) is also an optimum of problem (7.1)–(7.2) but with power $Pe^2Y^3$ instead of $Pe^2$. Comparing (7.4)–(7.5) with (7.1)–(7.2) we have

(7.6)\begin{gather} Y Nu_Y = Nu_1(Pe \, Y^{3/2}), \end{gather}

(7.7)\begin{gather}L_{x,Y}/Y = L_{x,1}(Pe \, Y^{3/2}). \end{gather}

If we assume power-law scalings

(7.8a,b)\begin{equation} Nu_1(Pe) = a\,Pe^\alpha,\quad L_{x,1}(Pe) = b\, Pe^\beta, \end{equation}

then

(7.9a,b)\begin{equation} Nu_Y(Pe) = Y^{{-}1+3\alpha/2}Nu_1(Pe), \quad L_{x,Y}(Pe) = Y^{1+3\beta/2} L_{x,1}(Pe). \end{equation}

For the problem of optimizing the flow only but keeping the walls horizontal in Alben (2023), we have $\alpha \lesssim 0.575$, at least up to $Pe = 10^7$. Thus $-1+3\alpha /2 < 0$, so $Nu_Y$ decreases as $Y$ increases. In this paper we optimize the flow and the wall shapes, and find asymptotically $\alpha = 2/3$. Thus $-1+3\alpha /2 = 0$, so $Nu_Y$ is independent of $Y$. The power-law behaviour is only approximate, however. For flows A and B in figure 8, with $Nu$ and the corresponding optimal $L_x$ shown in figure 6, $Nu$ increases slightly more slowly than $Pe^{2/3}$ but approaches this behaviour at large $Pe$. Therefore, $Nu_Y$ decreases slightly as $Y$ increases, and this is confirmed by comparing $Nu_Y$ for flows A and B with $Y = 0.75$, 1 and 1.5. For the peak $Nu$ among the computed wavy-wall optima, i.e. the highest black plus signs at each $Pe$ in figure 6, the slopes of the polygonal curve connecting adjacent points fluctuate more strongly about 2/3, from 0.64 to 0.7, than for the $Nu$ data for flows A and B. This may be because the highest black plus signs correspond to different flows with different wall shapes and streamline distributions as $Pe$ varies, and thus the $Nu$ values are more strongly affected by the finiteness and randomness of the ensemble of optima. For example, different slopes (but in a similar range, $[0.64, 0.69]$) are obtained if instead the second-highest black plus signs at each $Pe$ are used. However, if the same streamline and wall-shape configuration were used at each $Pe$, as with flows A and B, we might again obtain a slight decrease in $Nu_Y$ with $Y$, consistent with the wall separation distance always assuming the smallest allowed value, 1.

As a check, we also compute $L_{x,Y}$ for flows A and B with $Y = 0.75$, 1 and 1.5, and find very good agreement with $L_{x,Y} \sim Y^{0.5}$, consistent with (7.9), with $\beta = -1/3$ from (6.10).

8. Alternative optimization problems

We now describe a few alternative optimization problems that give additional insights. In the first alternative, we computed solutions that optimize the average rate of heat transfer per unit arc length of the boundary instead of per unit horizontal length. That is, we replaced $1/L_x$ in (2.2) with $1/s_{total}(L_x)$, where

(8.1)\begin{equation} s_{total}(L_x) = \int_0^{L_x} \sqrt{1+\frac{{\rm d}y_{bot}}{{\rm d}\kern 0.06em x}^2} + \sqrt{1+\frac{{\rm d}y_{top}}{{\rm d}\kern 0.06em x}^2}\, {{\rm d}\kern 0.06em x}. \end{equation}

We found that all the optima had nearly flat walls. It is not surprising that with this performance measure, the flat-wall optima would outperform the wavy-wall optima with $Nu$ plotted in figure 6, at least in the limit of large $Pe$. The ratio $L_x/s_{total}(L_x) \sim Pe^{-1/3}$, so the $Pe^{2/3}$ scaling in figure 6 would be reduced to $Pe^{1/3}$, while the flat-wall cases would remain the same, ${\approx } Pe^{0.575}$.

A second alternative that we considered (at the suggestion of a referee) is to optimize the shape of one wall only, keeping the other wall flat. Here we find that the optima combine the features of the flat-wall and wavy-wall optima. For $Pe \geq 10^5$ they are branched at the flat wall, but remain unbranched near the wavy wall. Examples of optima and further details are given in Appendix D.

A third alternative that we considered (also at the suggestion of a referee) is to study uniformly differentiable wall shapes, which more strongly limits the non-flatness of the walls. We modify (5.1) to

(8.2)\begin{equation} \left.\begin{array}{c} \displaystyle y_{bot}(p) ={-}L_x A \left(\dfrac{1+\sin{A_1}}{2}\right)\dfrac{h_1(p)^2}{\|2 h_1(p) h'_1(p)\|_\beta} ,\\ \displaystyle y_{top}(p) = L_x A \left(\dfrac{1+\sin{A_2}}{2}\right)\dfrac{h_2(p)^2}{\|2 h_2(p) h'_2(p)\|_\beta}, \end{array}\right\} \end{equation}

which corresponds to bounding $\|\textrm {d}y_{bot}/{\textrm {d}\kern 0.06em x}\|_\beta$ and $\|\textrm {d}y_{top}/{\textrm {d}\kern 0.06em x}\|_\beta$ by $A$. We use $\beta = 2$ and 6 to test bounds on the root-mean-square slope and an approximate maximum slope. For each $\beta$ we take $A = 1$ with $M_1 = 1$, 2 and 3 and $A = 3$ with $M_1 = 1$.

Examples of optima at $Pe = 10^5$ and $10^6$ are shown in figures 9(a–d) and 9(e–h), respectively, and within each row the optima have different choices of norms (i.e. $\beta$), numbers of modes ($M_1$) and bounds on the norms ($A$). At the larger $A$ value of 3 (figure 9a,c,e,g), $Nu$ is larger (shown above each optimum) and there is less branching, particularly in the top row. Thus the optima seem to interpolate between those with flat walls, with more branching, and the optima with the less restrictive constraint (5.1), which did not have branching. The latter constraint (5.1) bounded the deflection but allowed for much larger slopes (and deflections) than (8.2). Without branching (i.e. figures 7 and 9a), the convection rolls are often aligned with the peaks in wall deflection. In the branched flows of figure 9(b–h), the branches often meet the walls at local maxima of wall deflection, but not in every case.

Figure 9. Optimal wall shapes with bounds on the wall slope. (a–d) Optima with $Pe = 10^5$ and bounds on the 2-norm (a,b) or 6-norm (c,d) of the wall slope of $A = 3$ (a,c) with $M_1 = 1$ or $A = 1$ (b,d) with $M_1 = 3$. (e–h) Optima at the same parameters as (a–d) but with $Pe$ increased to $10^6$.

A fourth alternative that we considered (again at the suggestion of a referee) is $\alpha$-Hölder types of bounds (Toppaladoddi et al. 2021) on the wall roughness. For this case we add another version of (5.1) with

(8.3)\begin{equation} \left.\begin{array}{c} \displaystyle y_{bot}(p) ={-}A L_x^\alpha \left(\dfrac{1+\sin{A_1}}{2}\right)\dfrac{h_1(p)^2}{\|h_1(p)^2\|_6},\\ \displaystyle y_{top}(p) = A L_x^\alpha \left(\dfrac{1+\sin{A_2}}{2}\right)\dfrac{h_2(p)^2}{\|h_2(p)^2\|_6}. \end{array}\right\}\end{equation}

In the large-$Pe$ regime, $L_x \to 0$. Since $h_1$ and $h_2$ are superpositions of just a few Fourier modes, their slopes are uniformly bounded as $L_x \to 0$. Thus, given two points $p_1$ and $p_2$, $|y_{top/bot}(p_1)-y_{top/bot}(p_2)| \sim L_x^\alpha$ while $|x(p_1)-x(p_2)|^\alpha \sim L_x^\alpha$, since the horizontal domain size is $L_x$. Thus (8.3) corresponds to Hölder exponent $\alpha$ (Toppaladoddi et al. 2021). The original bound (5.1) corresponds $\alpha = 0$. We now consider $\alpha = 0.25$, 0.5 and 1. For each $\alpha$ we take $A = 1$ with $M_1 = 3$ and $A = 3$ with $M_1 = 1$. With $\alpha = 0.25$ and 0.5, the optima are generally not branched, and resemble those with the original constraint (5.1), but have smaller wall deflection amplitudes and smaller $Nu$. More distinct optima are found at $\alpha = 1$, particularly at large $Pe$, and examples are shown in figure 10 at $Pe = 10^6$ and $10^7$. At each $Pe$, the largest $Nu$ values occur with the largest wall deflections, $A = 3$ (figure 10a,e), and the flows are not branched. In the remaining cases, which have $A = 1$, the flows are branched, and the branches are much more strongly aligned with the peaks in wall deflection than in the previous case (figure 9). Taking $\alpha = 1$ is approximately the same as bounding the wall slope, but the different forms of (8.2) and (8.3) lead to distinct types of optima.

Figure 10. Optima with $\alpha$-Hölder constraints, for $\alpha = 1$, and $Pe = 10^6$ (a–d) and $10^7$ (e–g). Panels (a,e) have a relatively larger amplitude and smaller number of modes ($A = 3$, $M_1 = 1$) than (b–d) and (f,g), respectively. These panels have the parameter values reversed ($A = 1$, $M_1 = 3$).

In figure 11 we plot the $Nu$ values for optima with the alternative forms of the wall shapes (8.2) and (8.3) together with the values for flat walls (red crosses) and the original constraint on the wall deflection (5.1) (black plus signs). Here $Nu$ is divided by $Pe^{2/3}$, so the black plus signs follow a nearly horizontal trend. With the slope constraint (8.2), i.e. the examples in figure 9, the green squares and plus signs correspond to $A = 1$, and are not much higher than the flat-wall optima (red crosses). Increasing $A$ to 3, the optima are less constrained and thus have larger $Nu$, shown by the blue squares and plus signs. With the $\alpha$-Hölder constraint (8.3), the $Nu$ values (yellow, cyan and magenta circles and triangles) are higher as a group. The allowed wall deflection magnitudes increase as $\alpha$ decreases from 1 to 0.25, and thus $Nu$ for the yellow symbols are generally below $Nu$ for the cyan symbols, while the magenta symbols have the largest $Nu$ within this class, at least at the highest $Nu$. There is considerable overlap among these three sets because the allowed wall deflections depend strongly on $A$, which is 1 or 3. Overlap also results from the randomness of the initialization in the optimization routine. Other local optima were found with $Nu$ values below the red crosses, but these are omitted. Figure 11 shows that with these alternative constraints the optima interpolate between the flat-wall optima and those with the less stringent deflection constraint (5.1). The power laws corresponding to these alternative constraints probably depend on the various parameters ($A$, $\alpha$, $\beta$, $M_1$) in complicated ways, and we do not attempt to quantify the dependences.

Figure 11. Values of $Nu$ divided by $Pe^{2/3}$ for optimal flows with flat walls (red crosses) and with optimal wall shapes of various functional forms. The green and blue symbols correspond to the constraint on the wall slope (figure 9) with $A = 1$ and 3, respectively. The yellow, cyan and magenta symbols correspond to walls with Hölder exponents 1, 0.5 and 0.25 (figure 10). The black plus signs correspond to the original (mild) constraint on the wall deflection (5.1) with $A \leq 6$.

9. Discussion and conclusions

We have studied the problem of finding incompressible flows and wall shapes that maximize heat transfer between two walls, one hot and one cold, with a fixed rate of viscous power dissipation for the flow. In order to avoid the singularity when the walls touch, we defined a minimum allowed separation distance between the walls and used it as the characteristic length scale. We showed that with no flow, flat walls maximize $Nu$. With small flow strength ($Pe$), we showed that the effect of the walls and the flow approximately decouple, so flat walls are also optimal over some interval of $Pe$ values starting from zero. We then described our computational optimization method, and explained how limits in resolving the temperature gradient relate to limits in the walls’ horizontal periods, maximum amplitudes and the wavenumbers of the components defining the walls’ shapes. We also verified the accuracy of the decoupled approximation.

At moderate-to-large $Pe$, we described the transition to optima with wavy walls, and then compared the best optima with wavy and flat walls. The best wavy-wall optima are essentially the same as $Pe$ increases, except for a $Pe^{-1/3}$ scaling of the horizontal period $L_x$. By contrast, the flat-wall optima have an additional boundary-layer structure near the walls, that transitions to branching flows above $Pe = 10^5$. For the optima with wavy walls, $Nu$ scales as $Pe^{2/3}$. We explained the $L_x$ and $Pe$ scalings by studying model flows and showing that the optimal $L_x$ correspond to flows at the interface between the diffusion-dominated and advection-dominated regimes.

We now briefly explain that the $Pe^{2/3}$ scaling corresponds to an upper bound for all incompressible flows with wavy walls, shown by Goluskin & Doering (2016) using the background method. First, Souza (2016) showed that integrating the energy balance equation for Rayleigh–Bénard convection at a given Rayleigh number $Ra$ over time and space shows that

(9.1)\begin{equation} Ra(Nu -1) = Pe^2. \end{equation}

Goluskin & Doering (2016) used the background method to show that

(9.2)\begin{equation} Nu - 1 \leq C\, Ra^{1/2} \end{equation}

for Rayleigh–Bénard convection with rough walls (with $y$ a single-valued function of the horizontal coordinate, as assumed here). Combining (9.1) and (9.2),

(9.3)\begin{equation} Nu -1 \leq C\,Pe(Nu -1)^{{-}1/2}, \end{equation}

so $Nu-1 \leq C\,Pe^{2/3}$, and thus the optimal flow and wall configurations in this work achieve the maximum scaling with $Pe$ for rough walls.

We also used a scaling argument to show that the optima should have the minimum allowed separation between the hot and cold walls, at least for the model flows. These have fixed streamline configurations that are scaled horizontally by the optimal $L_x$ at each $Pe$, like the computed optimal flows and wall shapes.

We considered several alternative optimization problems. In the first, we optimized the average rate of heat transfer per unit arc length of the wall, and found that all the optima had nearly flat walls. In the second alternative, we went back to optimizing the average rate of heat transfer per unit horizontal length, but we optimized the shape of one wall only, keeping the other flat. The optima included branching at the flat wall but not at the wavy wall. In the third and fourth alternatives, we posed the wall shapes to have constrained slopes and Hölder exponents, respectively. As the constraints were relaxed, the optima transitioned from branching flows with flat walls to non-branching flows and larger $Nu$ for highly corrugated walls, as in the original optimization problem considered here.

Previous work on natural convection studied combinations of wall roughness amplitude and wavelength that maximized the scaling exponent of $Nu$ with $Ra$ (Toppaladoddi et al. 2017; Zhu et al. 2017). It was found that for a given roughness amplitude, the optimal roughness wavelength approximately equals the roughness amplitude. The optimal scaling $Nu \sim Ra^{1/2}$ occurs over an intermediate range of $Ra$ in which the thermal boundary layer thickness is comparable to the roughness amplitude. At larger and smaller $Ra$, the scaling is the same as that without roughness. Although we do not consider the Rayleigh–Bénard system here, the present results are connected in the sense that in the $Nu$-maximizing configurations, the wall roughness has a horizontal length scale of the same order as the thickness of the thermal boundary layers running along the almost vertically sloped sections of the walls.