Enhancing wall-to-wall heat transport with unsteady flow perturbations

Abstract

We determine unsteady time-periodic flow perturbations that are optimal for enhancing the time-averaged rate of heat transfer between hot and cold walls (i.e. the Nusselt number Nu), under the constraint of fixed flow power (Pe$^2$, where Pe is the Péclet number). The unsteady flows are perturbations of previously computed optimal steady flows and are given by eigenmodes of the Hessian matrix of Nu, the matrix of second derivatives with respect to amplitudes of flow mode coefficients. Positive eigenvalues of the Hessian correspond to increases in Nu by unsteady flows, and occur at $Pe\geqslant 10^{3.5}$ and within a band of flow periods $\tau \sim Pe^{-1}$. For $\tau {\textit{Pe}}\leqslant 10^{0.5}$, the optimal flows are chains of vortices that move along the walls or along eddies enclosed by flow branches near the walls. At larger $\tau {\textit{Pe}}$, the vorticity distributions are often more complex and extend farther from the walls. The heat flux is enhanced at locations on the walls near the unsteady vorticity. We construct an iterative time-spectral solver for the unsteady temperature field, and find increases in Nu of up to 7 % at moderate-to-large perturbation amplitudes.

1. Introduction

One of the most widely studied applications of fluid dynamics is convective heat transfer. This field is closely connected to fundamental areas of fluid dynamics such as boundary layers, turbulence and multiphase flows (Layton et al. 1988; Kaviany 1994). Many studies have considered how to manipulate boundary layers and turbulence to enhance convective heat transfer from solid surfaces (Webb & Kim 2005). ‘Passive’ techniques include roughening the solid surfaces and inserting vortex generators or other obstacles to enhance fluid and thermal mixing within the boundary layer (Webb & Kim 2005; Glezer, Mittal & Alben 2016). ‘Active’ techniques involve applying local forcing to the fluid (or the boundaries) using external electric fields, jet impingement, fluid injection, stirring, or surface vibrations (Yabe, Mori & Hijikata 1996; Webb & Kim 2005; Laohalertdecha, Naphon & Wongwises 2007; Fang & Khan 2013; Léal et al. 2013; Shank & Tiari 2023; Wang & Alben 2025).

Hassanzadeh, Chini & Doering (2014) formulated an optimisation problem in order to study upper bounds on the rate of convective heat transfer in a simple geometry. They considered two-dimensional (2-D) steady incompressible fluid flows in a layer between hot and cold horizontal walls. They solved the Euler–Lagrange equations numerically, and in certain limits analytically, to determine the flows that maximise the rate of heat transfer (the Nusselt number Nu) subject to a given rate of viscous power dissipation (Pe $^2$ , with Pe the Péclet number). The problem set-up was motivated by the long-standing search for upper bounds on heat transfer in natural (Rayleigh–Bénard) convection (Sondak, Smith & Waleffe 2015; Wen et al. 2020, 2022; Lohse & Shishkina 2024). The optimal flows in Hassanzadeh et al. (2014) are constrained only by incompressibility and fixed power dissipation, not by the buoyancy-driven Navier–Stokes equations (the Oberbeck–Boussinesq equations) of natural convection. In fact, the optimal flows are solutions to the Navier–Stokes equations with a suitable forcing term that may be taken as a model for the forcing in active heat transfer enhancement. More generally, the structures of these optimal flows improve our understanding of the fluid mechanics of efficient heat transfer in both natural and forced convection.

The ‘wall-to-wall’ optimisation problem of Hassanzadeh et al. (2014) was studied and extended by several subsequent works. Souza and others extended the work to no-slip instead of free-slip boundary conditions, and higher-power dissipation (i.e. Pe) (Souza 2016; Souza, Tobasco & Doering 2020). They also studied optimal transport in the Lorenz equations, a low-mode truncation of the wall-to-wall problem (Souza & Doering 2015a , b ). Tobasco and Doering constructed analytical branching flows that saturate an upper bound on the rate of heat transfer up to a logarithmic correction (Tobasco & Doering 2017; Doering & Tobasco 2019). Alben computed 2-D optimal flows using an adjoint method up to Pe = 10 $^{7.5}$ and found branched structures for Pe $\gtrsim$ 10 $^{4.5}$ (Alben 2023). Previously, Motoki and others had solved the Euler–Lagrange equations for optimal steady three-dimensional (3-D) flows computationally using a continuation method (Motoki et al. 2018a , b ). For the wall-to-wall problem, they found a bifurcation from 2-D convection rolls to 3-D optima at Pe $\approx$ 80 (Motoki et al. 2018a ). The flows show a tree-like pattern with branching near the walls for Pe $\gtrsim 500$ , much lower than in two dimensions. Kumar (2024) constructed 3-D branching ‘pipe’ flows analytically that attain the upper bound without the logarithmic correction. The upper bound was also found for 2-D walls with elongated convection rolls by allowing for wavy wall shapes (Alben 2024). The same optimisation framework was used to study optimal flows through channels (Alben 2017a ) and other geometries (Alben 2017b ). A body of work has also found improved heat transfer in natural convection from rough walls (Toppaladoddi et al. 2017, 2021). Other recent work has considered optimal cooling of the interior of domains (Marcotte et al. 2018; Song, Fantuzzi & Tobasco 2023).

In this work, we extend the wall-to-wall optimisation problem to the unsteady case, in 2-D space. Solving the discretised steady advection–diffusion equation in the advection-dominated regime requires a fine mesh in order to resolve very thin temperature boundary layers and thin plumes that move far from the boundary, into the main flow. In the adjoint-based optimisation method of Alben (2023), the problem was formulated as an unconstrained optimisation problem in which the Nusselt number was maximised using BFGS, a Hessian-free iterative method. At each iteration, a backtracking line search was performed in which the steady advection–diffusion equation was solved multiple times, to find a flow that yields a larger Nu. About $10^3{-}10^4$ iterates were computed for each initial guess for the optimal flow. Multiple initial guesses were used at each Pe to search for multiple local optima. A wide range of Pe values was used, $[0,10^{7.5}]$ . Altogether, there were many solves of the advection–diffusion equation, and the total computational cost was substantial.

In the unsteady case, this optimisation approach would be orders of magnitude more expensive. Adding the time dimension to the two spatial dimensions makes solving the advection–diffusion equation much more expensive. An explicit time-stepping scheme converges very slowly to the time-periodic solution in the advection-dominated regime. Instead, we adopt an implicit time-spectral method here, which results in a sparse linear system with $O(10^6{-}10^7)$ unknowns, and requires a iterative solver. The condition number of the linear system is large, again because we work in the advection-dominated regime, so many iterations are required in some cases. In addition, the search occurs in a much higher dimensional space than the steady case, with the additional coefficients needed to describe the time dependence of the flow modes. In order to mitigate the computational cost, we restrict the optimisation problem to the limit of small unsteady flow perturbations, using the approach of Alben, Prabala & Godek (2024). However, we evaluate the performance of the optimal perturbations with amplitudes ranging from small to large.

The unsteady 2-D wall-to-wall heat transfer optimisation problem was considered previously by Souza (2016). A gradient ascent method was proposed, but results were limited to the steady case due to computational challenges. For the unsteady problem, a low-mode approximation was considered: the Lorenz equations, with one flow mode and two temperature modes. The ‘background’ method and an optimal control method were used to show that steady flows are optimal in this approximation. The same methods were used to show that steady flows are also optimal for the ‘double Lorenz’ equations, with four flow modes and four temperature modes. In this paper, we investigate whether steady flows are optimal for the unsteady advection–diffusion equation.

We compute optimal flow modes in the limit of small unsteady perturbations of the best currently known steady 2-D flows from Alben (2023). We show that unsteady flows outperform these steady flows above a critical Pe value. The unsteady flows are defined by a finite set of eigenmodes of a Hessian matrix. We test the unsteady perturbations with amplitudes ranging from small to large, and find modest but noticeable increases in the rate of heat transfer (Nu).

In § 2, we introduce the problem. Section 3 explains the perturbative approach for computing optima. The numerical discretisation is given in § 3.1, followed by the modal expansions of the perturbations in § 4. Section 5 shows the perturbations that can increase Nu, indicated by positive eigenvalues of the Hessian in the quadratic term in the expansion of Nu. The corresponding eigenmodes’ vorticity and temperature fields are described. Section 6 shows the temperature fields and Nu values when the perturbations are employed with a range of amplitudes from small to large. Section 7 gives the conclusions.

2. Problem set-up

The ‘wall-to-wall’ problem has been described in a series of works (Hassanzadeh et al. 2014; Souza 2016; Tobasco & Doering 2017; Motoki et al. 2018a ; Doering & Tobasco 2019; Souza et al. 2020; Alben 2023), so we give a brief description of the problem – its geometry, equations and boundary conditions. A layer of fluid is contained between horizontal walls at $y = 0$ and 1, with temperatures 1 and 0, respectively (an example is shown in figure 1 a). We have chosen to non-dimensionalise all quantities with dimensions of length (including $y$ ) by the layer height $H$ , and quantities with units of temperature by the temperature difference between the walls $\Delta T$ . Given an incompressible flow field $\boldsymbol{u} = (u(x,y,t),v(x,y,t))$ , we solve the unsteady advection–diffusion equation for the temperature $T$ :

(2.1) \begin{align} \partial _tT + u\,\partial _x T + v\,\partial _y T - {\nabla} ^2T = 0. \end{align}

In the dimensional version of (2.1), the prefactor of ${\nabla} ^2T$ is the thermal diffusivity $\kappa$ ; it is set to unity in the dimensionless (2.1), corresponding to non-dimensionalising velocities by $\kappa /H$ . Figure 1(a) shows the temperature field for an example flow in figure 1(b) that is steady and horizontally periodic with period 1.72. The flow has counter-rotating convection rolls (anticlockwise on the left, clockwise on the right) that result in an upward plume of hot fluid in the middle of the $x$ range in figure 1(a), and a downward cold plume that is centred at the periodic boundary at $x$ = 0 and $x$ = 1.72.

Figure 1. Temperature field (a) corresponding to a horizontally periodic flow with streamlines shown in (b). The flow has two counter-rotating convection rolls per unit cell, and Pe = $10^2$ .

We search for a 2-D incompressible flow field that maximises the Nusselt number Nu, the rate of heat transfer out of the bottom wall, averaged over time and horizontal period length:

(2.2) \begin{align} {\textit{Nu}} &= \frac {1}{\tau }\int _0^\tau \frac {1}{L_x}\int _0^{L_x} -\partial _y T \big |_{y = 0}\, {\rm d}x \, {\rm d}t. \end{align}

Here, we assume a time-periodic temperature field with period $\tau$ (i.e. due to a flow with the same period). Using (2.1) and the divergence theorem, one can show that Nu is also the average rate of heat transfer into the top wall. We maximise Nu over flows that satisfy certain constraints. First, we require the flow to be incompressible, automatically imposed by setting $(u,v) = \boldsymbol{u} = (\partial _y\psi ,-\partial _x\psi)$ , for a given stream function $\psi (x,y,t)$ . We require zero flow velocity (the no-slip condition) at the walls: $\psi = \partial _y\psi = 0$ at $y = 0$ , and $\psi = \psi _{top}$ and $\partial _y\psi = 0$ at $y = 1$ . The constant $\psi _{top}$ is the net horizontal fluid flux through the channel. We will take $\psi _{top} = 0$ to be consistent with the previous steady optimal flows, discussed further below. One of the main constraints is a fixed average rate of viscous energy dissipation, i.e. the time- and space-averaged power needed to maintain the flow, in dimensionless form:

(2.3) \begin{align} \mbox{power} = \frac {1}{\tau } \int _0^\tau \int _0^1 \frac {1}{L_x} \int _0^{L_x} 2e_{\textit{ij}}^2\, {\rm d}x\, {\rm d}y\, {\rm d}t, \end{align}

with the spatial integral taken over a period cell. Here, $e_{\textit{ij}} = (\partial _{x_i}u_{\!j} + \partial _{x_{\!j}}u_i)/2$ is the symmetric part of the velocity gradient tensor, and $e_{\textit{ij}}^2$ is the sum of the squares of its entries. Thus

(2.4) \begin{align} 2 e_{\textit{ij}}^2 &= 2\, \partial _x u^2 + (\partial _y u + \partial _x v)^2 + 2\,\partial _y v^2. \end{align}

We have non-dimensionalised the space- and time-averaged power in (2.3) by $\mu \kappa ^2/H^4$ , where $\mu$ is the fluid viscosity. We now write $u$ and $v$ in terms of $\psi$ to obtain

(2.5) \begin{align} \mbox{power} &= \frac {1}{\tau }\int _0^\tau \int _0^1 \frac {1}{L_x}\int _0^{L_x} \left (\partial _{xx}\psi - \partial _{yy}\psi \right)^2 + 4\,\partial _{xy}\psi ^2 \,{\rm d}x\,{\rm d}y\,{\rm d}t \nonumber\\ & = \frac {1}{\tau }\int _0^\tau \int _0^1 \frac {1}{L_x} \int _0^{L_x} \left ({\nabla} ^2\psi \right)^2 + 4\left (\partial _{xy}\psi ^2-\partial _{xx}\psi\, \partial _{yy}\psi \right){\rm d}x\,{\rm d}y\,{\rm d}t . \end{align}

The second term in the integrand of (2.5) can be written as the divergence of $2(-u\,\partial _y v+v\,\partial _y u,-v\,\partial _x u+u\,\partial _x v)$ , and after applying the divergence theorem, this part integrates to zero when the flow is periodic or has zero tangential component on the sides of a rectangular domain, as is the case here. (A similar derivation appears in Lamb 1932, art. 329). The basic problem is to maximise Nu over a set of incompressible flows with power Pe $^2$ , a constant:

(2.6) \begin{align} \max _\psi {\textit{Nu}}(\psi) \quad \mbox{subject to} \quad \begin{cases} \partial _tT + \partial _y\psi\, \partial _x T -\partial _x\psi\, \partial _y T - {\nabla} ^2T = 0, \\ \frac {1}{\tau }\int _0^\tau \int _0^1 \frac {1}{L_x} \int _0^{L_x} \left ({\nabla} ^2\psi \right)^2 {\rm d}x\,{\rm d}y\,{\rm d}t = {\textit{Pe}}^2, \end{cases} \end{align}

where Pe stands for the Péclet number, a measure of the dimensionless flow velocity. The boundary conditions for $T$ in problem (2.6) are given above (2.1). The class of $\psi$ is left unspecified because different choices are possible. In Alben (2023), we computed optimal steady flows $\psi$ , with no-slip boundary conditions at the top and bottom walls, and horizontal periodicity.

In this paper, we compute optimal unsteady flows, assuming that they are small unsteady perturbations of the optimal steady flows. Extending the algorithm in Alben (2023) directly to unsteady flows is too expensive, as we will discuss later, but the perturbative assumption makes the computations feasible. However, we will evaluate the resulting unsteady flow perturbations at both small and large amplitudes. The optimal steady flows in Alben (2023) were computed for Pe $\in [0, 10^{7.5}]$ , and the horizontal period $L_x$ was optimised together with $\psi$ . We found $\psi _{top} = 0$ to be optimal in all cases. Since our unsteady flows will be perturbations of the steady optima, for simplicity we will assume $\psi _{top} = 0$ and that the flow is periodic in $x$ with the same period as the optimal steady flow at a given Pe, for $L_x$ given in Alben (2023).

In summary, the problem has the geometry and boundary conditions of Rayleigh–Bénard convection, but instead of solving the Oberbeck–Boussinesq equations for the flow field, we determine the flow field that maximises the Nusselt number, the average rate of heat transfer from the boundaries. The optimal flows are solutions of the Navier–Stokes equations with a certain force per unit volume distributed over the flow domain, in other words, to a forced convection flow. The forcing function is not computed explicitly here, but a method for doing so is discussed in Appendix A of Alben (2017a ). However, the average rate of work done by the forcing function corresponds to the average rate of viscous dissipation Pe $^2$ .

3. Unsteady perturbations of steady optima

The expansion of the unsteady flows as perturbations of the steady flows, the corresponding Nu expansion, and the computation of optimal unsteady flow modes, follows the same procedure as in Alben et al. (2024), with a few adjustments. The steady flows in Alben et al. (2024) were different – horizontal flows through a channel with cold fluid flowing inwards at the left boundary, heated walls (with $T$ = 1) at the top and bottom boundaries, and outflow at the right boundary, instead of the horizontally periodic problem here in which heat moves from the bottom to the top.

Figure 2. Streamlines of optimal steady flows computed in Alben (2023) at Pevalues (a) 10 $^2$ , (b) 10 $^{2.5}$ , (c) 10 $^3$ , (d) 10 $^{3.5}$ , (e) 10 $^4$ , (f) 10 $^{4.5}$ , (g) 10 $^5$ , (h) 10 $^{5.5}$ , (i) 10 $^6$ , (j) 10 $^{6.5}$ , (k) 10 $^7$ . For (h–k), inset panels to the right of the main flows show the streamlines near the top and bottom walls.

In figure 2, we show the base states of the perturbation expansion – the optimal steady flows from Alben (2023) at Pe ranging from 10 $^2$ to 10 $^7$ . These are the best flows found from 30–50 random initialisations using a BFGS iterative method. The flows consist of rectangular convection rolls with horizontal period $L_x \sim$ Pe $^{-0.36}$ . Small closed eddies occur near the walls at Pe = $10^4$ (figure 2 e) and above, and branching flows are optimal at Pe $\gt 10^5$ (figures 2 h–k). Up to and including Pe = $10^4$ (figure 2 a–e), most initialisations converged to the global optimum, with just a few resulting in a suboptimal flow that is slightly perturbed, usually near the walls. Between Pe = 10 $^{4.5}$ and Pe = 10 $^6$ (figure 2 f–i), the global optimum was found multiple times, but with decreasing frequency as Pe increased. Here, a wider range of suboptimal flows was found, with different branching patterns near the walls. At Pe = 10 $^{6.5}$ and 10 $^7$ (figure 2 j,k), the branching patterns are more complex with multiple levels, and several flows with small to moderate variations in branching are within 0.1 % of the maximum Nu found. Only the flows with the maximum Nu are shown in figure 2, and these are used as the base flows for the perturbation expansions.

Figure 3. Temperature fields corresponding to the flows in figure 2.

Figure 3 shows the temperature fields corresponding to the optimal steady flows. Upward and downward plumes are located at the streamlines that emanate from the walls, separating neighbouring convection rolls or enclosed eddies. Computing Nu accurately requires resolving very thin temperature boundary layers with fine meshes as discussed below in § 3.1.

Figure 4. Nusselt number versus Péclet number (black crosses) for steady optimal flows and temperature fields, including those shown in figures 2 and 3. The green and red lines show Nu $\sim$ Pe $^{0.575}$ and Nu $\sim$ Pe $^{2/3}$ (log Pe) $^{-4/3}$ , respectively.

Figure 4 shows the values of Nu for the optimal flows over a range of Pe, including those in figure 2. The green and red lines show Nu scalings proportional to Pe $^{0.575}$ and Pe $^{2/3}$ (log Pe) $^{-4/3}$ , respectively. The latter form was given for analytically constructed flows with a self-similar branched structure in Tobasco & Doering (2017) and Doering & Tobasco (2019). One of the main questions that we address is whether we can find unsteady flows with larger Nu than those of the steady optima in figure 4.

At a given Pe, we denote the optimal steady flow as $\psi _0(x,y)$ and consider an unsteady flow of the form

(3.1) \begin{align} \psi (x,y,t) &= \frac {1}{\sqrt {1+c_1\epsilon ^2}}\psi _0(x,y) + \frac {\epsilon }{\sqrt {1+c_1\epsilon ^2}}f_A(x,y)\cos (2\pi t/\tau)\notag\\&\quad + \frac {\epsilon }{\sqrt {1+c_1\epsilon ^2}}f_B(x,y)\sin (2\pi t/\tau). \end{align}

Here, $0 \lt \epsilon \ll 1$ and the constant $c_1$ will be chosen so that $\psi$ has the same time-averaged power as $\psi _0(x,y)$ , i.e. Pe $^2$ . We will now compute Nu for $\psi$ by expanding both $\psi$ and $T$ in powers of $\epsilon$ ,

(3.2) \begin{align} \psi (x,y,t) &= \psi _0(x,y) + \epsilon\, \psi _1(x,y,t) + \epsilon ^2\,\psi _2(x,y,t) + \cdots, \end{align}

(3.3) \begin{align} T(x,y,t) &= T_0(x,y) + \epsilon\, T_1(x,y,t) + \epsilon ^2\, (T_{2s}(x,y) + T_{2u} (x,y,t)) + \cdots , \end{align}

plugging both series into (2.1), and solving the equation at each power of $\epsilon$ . Then we plug the $T$ expansion (3.3) into (2.2) to obtain a similar expansion for Nu:

(3.4) \begin{align} {\textit{Nu}} &= {\textit{Nu}}_0 + \epsilon\, {\textit{Nu}}_1 + \epsilon ^2\, {\textit{Nu}}_2 + \cdots . \end{align}

In (3.3), we have separated $T_2$ into its time average $T_{2s}$ (the steady part) and the remainder $T_{2u}$ (the unsteady part); only the steady part will be needed for Nu, which is a time average. At $O(\epsilon ^0)$ , (2.1) is

(3.5) \begin{align} \partial _y\psi _0\, \partial _xT_0 - \partial _x\psi _0\,\partial _yT_0 - {\nabla} ^2T_0 = 0, \end{align}

thus Nu $_0$ is the Nusselt number for the optimal steady flow $\psi _0$ . At $O(\epsilon ^1)$ , (2.1) is

(3.6) \begin{align} \partial _t T_1 + \partial _y\psi _0\, \partial _xT_1 - \partial _x\psi _0\,\partial _yT_1 - {\nabla} ^2T_1 = - \partial _y\psi _1\, \partial _xT_0 + \partial _x\psi _1\,\partial _yT_0. \end{align}

Since the right-hand side has components proportional to $\cos (2\pi t/\tau)$ and $\sin (2\pi t/\tau)$ , the solution can be written as $T_1 = T_{1A}(x,y)\cos (2\pi t/\tau) + T_{1B}(x,y)\sin (2\pi t/\tau)$ . Inserting this expression into (3.6) gives a coupled system for $T_{1A}$ and $T_{1B}$ :

(3.7) \begin{align} \frac {2\pi }{\tau }T_{1B} + \partial _y\psi _0\, \partial _xT_{1A} - \partial _x\psi _0\,\partial _yT_{1A} - {\nabla} ^2T_{1A} &= -\partial _y f_A\, \partial _xT_0 + \partial _x f_A\,\partial _yT_0, \end{align}

(3.8) \begin{align} -\frac {2\pi }{\tau } T_{1A} + \partial _y\psi _0\, \partial _xT_{1B} - \partial _x\psi _0\,\partial _yT_{1B} - {\nabla} ^2T_{1B} &= -\partial _y f_B\, \partial _xT_0 + \partial _x f_B\,\partial _yT_0. \end{align}

Since $T_1$ has zero time average, Nu $_1 = 0$ . We need to compute higher-order terms in Nu to determine whether the unsteady flow perturbation can increase Nu. At $O(\epsilon ^2)$ , (2.1) is

(3.9) \begin{align} \partial _t T_2 + \partial _y\psi _0\, \partial _xT_2 - \partial _x\psi _0\,\partial _yT_2 - {\nabla} ^2T_2 &= - \partial _y\psi _1\, \partial _xT_1 + \partial _x\psi _1\,\partial _yT_1 \notag\\ &\quad - \partial _y\psi _2\, \partial _xT_0 + \partial _x\psi _2\,\partial _yT_0. \end{align}

Equation (3.1) gives $\psi _2 = -c_1\psi _0/2$ . Since $\psi _1$ and $T_1$ are proportional to $\cos (2\pi t/\tau)$ and $\sin (2\pi t/\tau)$ , the right-hand side of (3.9) has steady terms and terms proportional to $\cos (4\pi t/\tau)$ and $\sin (4\pi t/\tau)$ . The unsteady terms have zero time average, so only the steady terms make a non-zero contribution to Nu $_2$ . We compute it by first solving for $T_{2s}$ as

(3.10) \begin{align} \partial _y\psi _0\, \partial _xT_{2s} - \partial _x\psi _0\,\partial _yT_{2s} - {\nabla} ^2T_{2s} = - \frac {1}{2}\partial _yf_A\, \partial _xT_{1A} + \frac {1}{2}\partial _xf_A\,\partial _yT_{1A} \nonumber \\ - \frac {1}{2}\partial _yf_B\, \partial _xT_{1B} + \frac {1}{2}\partial _xf_B\,\partial _yT_{1B} + \frac {c_1}{2}\partial _y\psi _0\, \partial _xT_0 - \frac {c_1}{2}\partial _x\psi _0\,\partial _yT_0, \end{align}

and then computing

(3.11) \begin{align} {\textit{Nu}}_2 &= \frac {1}{L_x}\int _0^{L_x} -\partial _y T_{2s} \big |_{y = 0}\, {\rm d}x. \end{align}

The time average is dropped in (3.11) because the integrand is steady. To summarise, the leading-order change in Nu due to the unsteady flow perturbation is $\epsilon ^2\,Nu_2$ , and to calculate Nu $_2$ , we need to calculate $T_{1A}$ , $T_{1B}$ and $T_{2s}$ . This requires solving three elliptic partial differential equations – (3.7), (3.8) and (3.10) – two of which are coupled. In § 3.1, we discuss the numerical discretisation of the equations, and in § 4 we explain how the calculation of Nu $_2$ is used to determine the optimal unsteady perturbation as a superposition of flow modes.

3.1. Discretisation

We solve (3.7), (3.8) and (3.10) by representing each component of $\psi$ and $T$ in (3.2) and (3.3) on tensor product grids in $(x,y)$ space that are concentrated at the boundaries, the same as in Alben (2023). We start with unstretched coordinates $p$ and $Y$ . We define uniform grids on $[0, 1]$ in $p$ and $Y$ with spacings $1/m$ and $1/n$ , respectively, where $m$ = 256 and $n$ = 256, 512 or 1024. We then map these points to $(x,y)$ space:

(3.12) \begin{align} x = L_x p, \quad y = Y - \eta \frac {1}{2\pi }\sin (2\pi Y), \end{align}

with $\eta$ a stretching factor that is set to 0.997. We obtain a stretched grid with points concentrated near the boundary in $y$ , to resolve the boundary layers. The grid spacing in $y$ ranges from $(1-\eta)/n \approx 3\times 10^{-6}$ – $10^{-5}$ at the boundaries, to $(1+\eta)/n \approx 2\times 10^{-3}$ – $8\times 10^{-3}$ at the domain centre. We discretise all spatial derivatives with second-order finite differences using three or four grid points, one-sided at the boundaries. Integrals such as those in (2.2) and (2.5) are discretised with the trapezoidal rule.

4. Flow modes

Our overall goal is to maximise Nu over a large space of unsteady flow perturbations. The particular form of the perturbation in (3.1) is determined by the functions $f_A$ and $f_B$ , and the scalar $\tau$ . To represent a wide range of flows, we expand $f_A$ and $f_B$ using a large number of products of Fourier modes in $x$ and Chebyshev polynomials in $y$ , the same as for the steady optima in Alben (2023). The $2M+1$ Fourier modes are arranged as alternating cosines and sines, from small to large wavenumbers:

(4.1) \begin{align} X_{2j+1}(x) &= \cos (2\pi j x/L_x),\ j = 0,\ldots , M,\nonumber\\ X_{2j}(x) &= \sin (2\pi j x/L_x),\ j = 1,\ldots , M. \end{align}

We take $M = 5m/64$ , much smaller than $m$ , so the modes are resolved well by the horizontal mesh. The vertical modes $Y_k(y)$ , $k = 1, \ldots , N-3$ , are linear combinations of Chebyshev polynomials that satisfy no-slip at the walls. They are generated using an orthogonal projection procedure described in Alben (2023, app. B). We take $N = n/8$ , so the most oscillatory mode ( $k = N-3$ ) is resolved well by the vertical mesh.

We form all the products $X_{\!j}(x)\,Y_k(y)$ on the spatial grid, and assemble them as $(2M+1)(N-3)$ columns of a matrix $\boldsymbol{Z}$ . The problem is further simplified by making the flow modes orthogonal with respect to the power inner product (2.5), in discretised form. For a discretised flow given as a column vector $\boldsymbol{Z}_i$ multiplied by either $\cos (2\pi t/\tau)$ or $\sin (2\pi t/\tau)$ , we can write the discretised power integral (2.5) as a weighted inner product $(\sqrt {\boldsymbol{W}}\boldsymbol{L}\boldsymbol{Z}_i)^{\rm T}\sqrt {\boldsymbol{W}}\boldsymbol{L}\boldsymbol{Z}_i$ , where $\boldsymbol{L}$ is the discrete Laplacian, and $\boldsymbol{W}$ is the diagonal matrix of weights corresponding to the compound trapezoidal rule for the spatial integrals in (2.5), multiplied by $1/2$ for the time integral. We perform a QR factorisation

(4.2) \begin{align} \frac {1}{{Pe}}\boldsymbol{\sqrt {W}L}\boldsymbol{Z} = \boldsymbol{Q}\hat {\boldsymbol{R}}. \end{align}

Then $\boldsymbol{V} \equiv \boldsymbol{Z}\hat {\boldsymbol{R}}^{-1}$ has columns $\boldsymbol{V}_{\!j}$ that obey $(\sqrt {\boldsymbol{W}}\boldsymbol{L}\boldsymbol{V}_i)^{\rm T}\sqrt {\boldsymbol{W}}\boldsymbol{L}\boldsymbol{V}_{\!j}=Pe^2\,\delta_{\textit{ij}}$ . Thus the $\boldsymbol{V}_{\!j}$ are orthogonal with respect to the power inner product, and have power Pe $^2$ . These are the modes that are used as a basis for $f_A$ and $f_B$ in the unsteady part of the flow given by (3.1).

We use the modes of $\boldsymbol{V}$ to write a discrete version of the flow in (3.1):

(4.3) \begin{align} \boldsymbol{\varPsi } = \frac {1}{\sqrt {1+\|\boldsymbol{a}\|^2}} \boldsymbol{\varPsi }_0 + \sum _{j = 1}^{N_m} \frac {{a}_{j}}{\sqrt {1+\|\boldsymbol{a}\|^2}} \boldsymbol{V}_{\!j}\cos \left (\frac {2\pi t}{\tau }\right) + \sum _{j = 1}^{N_m} \frac {{a}_{j+N_m}}{\sqrt {1+\|\boldsymbol{a}\|^2}} \boldsymbol{V}_{\!j}\sin \left (\frac {2\pi t}{\tau }\right). \end{align}

Here, $\boldsymbol{V}_{\!j}$ are columns of $\boldsymbol{V}$ representing $N_m$ different spatial modes, and $\boldsymbol{a}$ is a vector of coefficients for the modes. The terms in (4.3) are mutually orthogonal with respect to the power inner product. Those with different temporal functions are orthogonal due to the time integral in (2.5), and those with different $\boldsymbol{V}_{\!j}$ are orthogonal by the definition of $\boldsymbol{V}$ . Hence $\boldsymbol{\varPsi }$ has power Pe $^2$ for any real vector $\boldsymbol{a}$ . When $\|\boldsymbol{a}\| = O(\epsilon)$ , the flow (4.3) is a discrete version of the flow (3.1). In this case, we can write a Taylor expansion for Nu like (3.4) but in powers of $\boldsymbol{a}$ :

(4.4) \begin{align} {\textit{Nu}}(\boldsymbol{a}) = {\textit{Nu}}(\boldsymbol{0}) + \boldsymbol{a}^{\rm T} D{\textit{Nu}}\big |_{\boldsymbol{a} = \boldsymbol{0}}+\frac {1}{2}\boldsymbol{a}^{\rm T} D^2{\textit{Nu}}\bigg |_{\boldsymbol{a} = \boldsymbol{0}} \boldsymbol{a} + \cdots, \end{align}

where ${\textit{Nu}}(\boldsymbol{0})$ = Nu $_0$ is Nu for the steady flow, and $D{\textit{Nu}}$ and $D^2{\textit{Nu}}$ are the gradient vector and Hessian matrix of first and second derivatives of Nu with respect to components of $\boldsymbol{a}$ , respectively. We compute the gradient and Hessian as the limits of finite differences. The gradient entries are

(4.5) \begin{align} D{\textit{Nu}}_i \big |_{\boldsymbol{a} = \boldsymbol{0}} = \frac {\partial {\textit{Nu}}}{\partial {a}_i}\bigg |_{\boldsymbol{a} = \boldsymbol{0}} = \lim _{\epsilon \to 0} \frac {{\textit{Nu}}(\epsilon \boldsymbol{e}_i)-{\textit{Nu}}(\boldsymbol{0})}{\epsilon } = {\textit{Nu}}_1(\epsilon \boldsymbol{e}_i)-{\textit{Nu}}_1(\boldsymbol{0}) = 0. \end{align}

The limit of the finite difference quotient in (4.5) is computed by inserting the perturbation expansion of Nu from (3.4) for the flows corresponding to $\boldsymbol{a} =\epsilon \boldsymbol{e}_i$ and $\boldsymbol{a} = \boldsymbol{0}$ . As explained below (3.8), Nu $_1$ = 0 for $\boldsymbol{a} =\epsilon \boldsymbol{e}_i$ for all $i$ , since the first-order temperature field has zero time average. We also have Nu $_1$ = 0 for $\boldsymbol{a} =\boldsymbol{0}$ since Nu $(\boldsymbol{0}) \equiv$ Nu $_0$ , so the remaining terms in (3.4) vanish. The entries of the Hessian matrix are also computed using (3.4):

(4.6) \begin{align} D^2{\textit{Nu}}_{\textit{ij}} \big |_{\boldsymbol{a} = \boldsymbol{0}} &= \frac {\partial ^2 {\textit{Nu}}}{\partial {a}_i\,\partial {a}_{\!j}}\bigg |_{\boldsymbol{a} = \boldsymbol{0}} = \lim _{\epsilon \to 0} \frac {{\textit{Nu}}(\epsilon \boldsymbol{e}_i\!+\!\epsilon \boldsymbol{e}_{\!j})-{\textit{Nu}}(\epsilon \boldsymbol{e}_i)-{\textit{Nu}}(\epsilon \boldsymbol{e}_{\!j})+{\textit{Nu}}(\boldsymbol{0})}{\epsilon ^2}\nonumber\\ &= {\textit{Nu}}_2(\epsilon \boldsymbol{e}_i\!+\!\epsilon \boldsymbol{e}_{\!j})-{\textit{Nu}}_2(\epsilon \boldsymbol{e}_i)-{\textit{Nu}}_2(\epsilon \boldsymbol{e}_{\!j}). \end{align}

These entries are non-zero in general, so the quadratic term in (3.4) is the leading-order change to Nu due to the unsteady flow perturbation. Thus for a given (small) magnitude of $\|\boldsymbol{a}\|_2$ , $\eta$ say, we maximise the quadratic term

(4.7) \begin{align} \max _{\|\boldsymbol{a}\|_2 = \eta } \frac {1}{2}\boldsymbol{a}^{\rm T} D^2{\textit{Nu}}\bigg |_{\boldsymbol{a} = \boldsymbol{0}} \boldsymbol{a}. \end{align}

The optimal flow perturbation $\boldsymbol{a}$ is the eigenvector of the Hessian with the largest eigenvalue, which is real since the Hessian is symmetric. Therefore, we compute the Hessian for a range of $\tau$ and Pe, compute the largest eigenvalues in each case, and examine the corresponding eigenvectors (i.e flows). We will examine not only the flow corresponding to the largest eigenvalue, but also the flows corresponding to other positive eigenvalues, to consider a wider range of beneficial flows.

For flows of the form (4.3), the Hessian is $2N_m\times 2N_m$ , with $N_m = (2M+1)(N-3)$ . The number of entries in each Hessian ranges from $ \approx 6 \times 10^6$ for $n = 256$ to $\approx 10^8$ for $n = 1024$ , reduced by a factor of 4 by symmetry. We compute the Hessian at a large number of $\tau$ and Pe values, so an efficient method is used to compute a single Hessian matrix, i.e. the large number of entries given by (4.7) for all $i$ and $j$ . The computing time ranges from 1 h to 2 h for $n$ = 256, to approximately 8 days for $n$ = 1024 on a single processor.

The sequence of computations for the Hessian is described in detail in Alben et al. (2024), so here we give a summary. First, given $\psi _0$ , we solve the zeroth-order problem (3.5) for $T_0$ . Next, we compute Nu $_2(\epsilon \boldsymbol{e}_i)$ for $i = 1,\ldots , N_m$ by setting the flow mode coefficient vector $\boldsymbol{a} = \epsilon \boldsymbol{e}_i$ for $i = 1,\ldots , N_m$ . Therefore, in (3.1), we set $f_A$ to column $i$ of $\boldsymbol{V}$ , $f_B$ to $\boldsymbol{0}$ , and $c_1$ to 1. Then we follow the remaining steps in that section to get Nu $_2$ : we solve the system (3.7)–(3.8) for $(T_{1A},T_{1B})$ , then solve (3.10) for $T_{2s}$ and compute Nu $_2$ from (3.11). To compute Nu $_2(\epsilon \boldsymbol{e}_i)$ for $i = N_m+1,\ldots , 2N_m$ , the steps are the same, but we now have a single sine term instead of a cosine term in (4.3). So we set $f_A$ to $\boldsymbol{0}$ , $f_B$ to column $i$ of $\boldsymbol{V}$ , and $c_1$ to 1. Having computed Nu $_2(\epsilon \boldsymbol{e}_i)$ in (3.11) for all $i$ , the same set of values gives Nu $_2(\epsilon \boldsymbol{e}_{\!j})$ for all $j$ .

For Nu $_2(\epsilon \boldsymbol{e}_i\!+\!\epsilon \boldsymbol{e}_{\!j})$ , we have two unsteady terms in (4.3) that may appear in either the cosine or sine summand. We put this flow in the form (3.1) by setting $c_1 = 2$ , and setting $f_A$ to the sum of the columns of $\boldsymbol{V}$ (if any) that appear in the cosine summand, and $f_B$ to the sum of the columns of $\boldsymbol{V}$ (if any) that appear in the sine summand. Then $(T_{1A},T_{1B})$ is simply the sum of those for Nu $_2(\epsilon \boldsymbol{e}_i)$ and Nu $_2(\epsilon \boldsymbol{e}_{\!j})$ , since the right-hand sides of (3.7)–(3.8) are linear in $f_A$ and $f_B$ . We then form the right-hand side of (3.10), solve it for $T_{2s}$ , and compute Nu $_2$ from (3.11).

Actually, we can obtain Nu $_2$ for each right-hand side of (3.10) without having to solve the matrix equation. In discretised form, (3.10)–(3.11) form a linear mapping from the right-hand-side vector of (3.10) to the scalar Nu $_2$ , and this mapping is given by an inner product of a fixed vector with the right-hand side of (3.10), as explained in Alben et al. (2024). So the dominant cost is simply forming the $\sim 10^7$ – $10^8$ right-hand-side vectors, which may still take several days on a single CPU due to the large number of them.

We mentioned that only about one-quarter of the Hessian entries need to be computed by symmetry. We show in Appendix A of Alben et al. (2024) that the Hessian is block skew-symmetric,

(4.8) \begin{align} D^2{\textit{Nu}} =\left [\begin{array}{cc} A &B \\ -B&A \\ \end{array}\right ], \end{align}

and since the Hessian is also symmetric, the blocks $A$ and $B$ are symmetric and skew-symmetric, respectively. Therefore, to form the Hessian, we only need to compute about half the entries of $A$ and $B$ , or about one-quarter of the entries of the Hessian.

We have described how the Hessian is computed for flows of the form (4.3), with just a single period. If the flow has unsteady perturbation modes with multiple periods, then we show in Appendix B of Alben et al. (2024) that the Hessian is block diagonal, with the blocks consisting of the Hessians for each period alone. Therefore, the problem is decoupled, with the eigenvalues and eigenvectors given by those for each period alone. Next, we compute the single-period blocks of the Hessian for various values of the period, and the corresponding eigenvectors and eigenvalues. We focus on cases with positive eigenvalues, so the eigenvectors give flow perturbations that increase Nu.

5. Eigenvalue distributions and eigenvectors (flow modes)

In figure 5, we plot the distributions of eigenvalues for a range of Pe values, from 10 $^2$ to 10 $^7$ (one value per panel). Within each panel, we plot the eigenvalue distributions across ranges of $\tau {\textit{Pe}}$ from 10 $^{-1}$ to 10 $^{2.5}$ , where the largest eigenvalues occur. We use $\tau {\textit{Pe}}$ instead of $\tau$ because we find that the peak eigenvalues occur at $\tau {\textit{Pe}}\in [10^0, 10^{1.5}]$ , so the optimal $\tau$ seems to scale approximately as Pe $^{-1}$ , as was found for channel flows also (Alben et al. 2024). We will discuss this scaling further shortly. In each panel, across $\tau {\textit{Pe}}$ , we plot selected eigenvalues ranging from the largest ( $\lambda _1$ ) to the smallest ( $\lambda _{N_m}$ ). The eigenvalues occur in identical pairs due to the block skew-symmetry of the Hessian (4.8), as explained by Alben et al. (2024). Physically, the pairs correspond to a phase shift by $1/4$ period in time, which takes the $\cos (2\pi t/\tau)$ , $\sin (2\pi t/\tau)$ components of a flow to $\sin (2\pi t/\tau)$ , $-\cos (2\pi t/\tau)$ components, but results in the same Nu.

Figure 5. Distributions of Hessian eigenvalues at Pe values (a) 10 $^2$ , (b) 10 $^3$ , (c) 10 $^{3.5}$ , (d) 10 $^4$ , (e) 10 $^{4.5}$ , (f) 10 $^5$ , (g) 10 $^{5.5}$ , (h) 10 $^6$ , (i) 10 $^{6.5}$ , (j) 10 $^7$ .

In some cases, the pairs approximately coalesce into a group or four (or more) nearly identical values. In figure 5, the leading three quartets of eigenvalues are shown with solid blue, red and green lines, respectively. In some cases (the right-hand sides of figures 5 a–f), two distinct blue lines can be seen, corresponding to two distinct pairs of values. In other cases (the left-hand sides of figure 5 a–f), the two blue lines overlap, sometimes together with the two red lines. These coalescences into groups of four or more correspond to symmetric groups of flows for the associated eigenvectors. In other words, the eigenvectors’ flow fields map into each other under reflections about the horizontal and/or vertical midlines of the domain, corresponding to the symmetries of the steady flows in figures 2(a–h). We will give concrete examples of such flows after discussing this figure.

After the leading three quartets, lines are given in figure 5 for every 20th eigenvalue up to $\lambda _{100}$ , then $\lambda _{200}$ , and the last three lines give a sense of the eigenvalue distribution near the bottom of the spectrum. These 15 lines (including two for each of the three leading quartets) give a sense of the overall distribution of the $N_m \approx$ 2400–10 000 eigenvalues. At Pe = 10 $^2$ (figure 5 a) and Pe = 10 $^3$ (figure 5 b), all eigenvalues are negative. At Pe = 10 $^{3.5}$ (figure 5 c), a band of positive eigenvalues appears. The dashed line in this and subsequent panels marks a range of values near zero that is omitted from the plot. This omission allows us to give one logarithmic scale for positive eigenvalues (above the dashed line) and another for negative eigenvalues (below the dashed line). The range of omitted values is chosen to be small enough to exclude all the values on the 15 lines. The peak eigenvalue in figure 5(c) occurs at $\tau {\textit{Pe}}$ = 10 $^{0.5}$ . A second peak occurs in figures 5(d–f) at $\tau {\textit{Pe}}$ = 10 $^{1.5}$ , followed by a single peak in figures 5(g–j), at the Pe with branched steady flows. As Pe increases from figures 5(c–j), both the number of positive eigenvalues at each $\tau {\textit{Pe}}$ and their maxima increase. These last two properties are also shared by the eigenvalues in the channel flow problem (Alben et al. 2024). There, the eigenvalue distributions are simpler and smoother, with a single peak that is more symmetrical about its maximum. This may be related to the simpler unidirectional form of the base steady channel flows.

Most of the eigenvalues lie between $\lambda _{200}$ and $\lambda _{N_m-600}$ , given by the black dotted lines and cyan solid lines. These are relatively close together on the negative part of the vertical axis, so there are no more than 200 positive eigenvalues in all these cases, out of $N_m \approx 5000$ in figures 5(c–f) and 10 000 in figures 5(g–j). At the bottom of the range, the eigenvalues decrease sharply from $\lambda _{N_m-600}$ to $\lambda _{N_m}$ , at $\tau {\textit{Pe}}\gtrsim 10^0$ .

Figure 6. Vorticity fields for a selection of leading eigenmodes at Pe = $10^{3.5}$ – $10^5$ , when the steady base flows are unbranched. From left to right, the four columns correspond to Pe = 10 $^{3.5}$ (a–c, $\tau {\textit{Pe}}$ = 10 $^{0}$ , 10 $^{0.5}$ , 10 $^{1}$ ), Pe = 10 $^{4}$ (d–g, $\tau {\textit{Pe}}$ = 10 $^{0}$ , 10 $^{0.5}$ , 10 $^{1}$ ), Pe = 10 $^{4.5}$ (h–k, $\tau {\textit{Pe}}$ = 10 $^{0}$ , 10 $^{0.5}$ , 10 $^{1}$ , 10 $^{1.5}$ ), and Pe = 10 $^{5}$ (l–q, $\tau {\textit{Pe}}$ = 10 $^{-0.5}$ , 10 $^{0}$ , 10 $^{0.5}$ , 10 $^{1}$ , 10 $^{1.5}$ , 10 $^{2}$ ). Thus the letter at the start of each panel label identifies Pe and $\tau {\textit{Pe}}$ . Following the letter is the mode number (in order from largest to smallest eigenvalue). In a few cases, the number is followed by ‘s’, denoting the $\sin (2\pi t/\tau)$ component. All other cases show the $\cos (2\pi t/\tau)$ component. In each column, the black dashed lines separate $\tau {\textit{Pe}}\leqslant 10^{0.5}$ and $\tau {\textit{Pe}}\geqslant 10^{1}$ .

We now show the corresponding eigenvectors, and only a small selection of those with positive eigenvalues. Even in this small set, we find a range of patterns. Figure 6 shows the vorticity fields for eigenvectors with Pe = $10^{3.5}$ – $10^5$ , corresponding to the unbranched steady flows (figure 2 d–g). The four columns from left to right show the top eigenmodes for Pe = 10 $^{3.5}$ , 10 $^{4}$ , 10 $^{4.5}$ , 10 $^{5}$ , respectively. Within each column, $\tau {\textit{Pe}}$ increases from top to bottom. Each vorticity field is labelled with a letter identifying the Pe and $\tau {\textit{Pe}}$ values, followed by the mode number (1 denoting the mode with largest eigenvalue), and then in a few cases ‘s’ for the $\sin (2\pi t/\tau)$ component. All cases without ‘s’ show the $\cos (2\pi t/\tau)$ component.

In the first column, letters (a–c) label modes at Pe = $10^{3.5}$ with $\tau {\textit{Pe}}$ = 10 $^{0}$ , 10 $^{0.5}$ , 10 $^{1}$ , respectively; Pe = $10^{3.5}$ is the lowest Pe in figure 5 with positive eigenvalues. The flows consist of chains of vortices of alternating sign near the walls, very similar to those for the channel flows in Alben et al. (2024) at small to moderate $\tau {\textit{Pe}}$ . In that work, the top modes consisted of a single vortex chain extending along the channel wall. Figure 6(a1) here shows the $\cos (2\pi t/\tau)$ component of the top mode, which has two vortex chains, corresponding to the two convection rolls of the steady flow, with streamlines shown by black dotted lines. The steady flow is rightwards and leftwards along the left- and right-hand sides of the bottom wall, respectively. Figure 6(a1s) shows the $\sin (2\pi t/\tau)$ component of the top mode, which also has two vortex chains, but shifted spatially in phase from those of the cosine mode by one-quarter of the spatial period of the vortex chain. When the cosine and sine components are superposed, the spatial and temporal shifts result in two counterpropagating travelling waves of vorticity moving at the speed of the steady flow near the vortices, as explained in Alben et al. (2024) for the case of a single travelling wave. Movie 1 in the supplementary movies are available at https://doi.org/10.1017/jfm.2025.10931 shows the same phenomenon at a higher Pe, 10 $^{4.5}$ . Figure 6(b1,b1s) show the corresponding modes with $\tau {\textit{Pe}}$ increased from 10 $^{0}$ to 10 $^{0.5}$ . The vorticity field is similar to figure 6(a1,a1s), but the vortices are fewer and larger. The enlargement of vortices with increasing $\tau$ Pe was also seen in Alben et al. (2024). Figure 6(c1), with $\tau$ Pe increased to 10 $^1$ , shows vorticity regions enlarged to half the horizontal domain size. This could be seen both as a continuation of the trend of larger discrete vortices and as part of a common transition to different vorticity patterns (usually more global) at $\tau$ Pe = 10 $^1$ and larger. This transition will be seen in the next groups of flows at higher Pe. We use black dashed lines in each column of figure 6 and the next three figures to indicate this transition with respect to $\tau$ Pe.

The vortex chains are at the bottom walls in figure 6(a–c). These eigenvectors actually occur in pairs (not shown), with the same vortex chains at either the top or bottom walls, and the same eigenvalues. Considering also the symmetry of the cosine and sine components of the eigenvectors mentioned earlier, we actually have quartets of identical eigenvalues.

Moving on to the second column, with the next larger Pe, 10 $^4$ , we now have positive eigenvalues at four $\tau$ Pevalues, with modes shown in figure 6(d–g). In figure 6(d,e), we show both the first and ninth modes, the latter representing the third quartet of modes, those with the smallest positive eigenvalues here. In figure 6(d9), the vortex chains are broken into two clusters each on the left- and right-hand halves, similarly to the eigenmodes below the top mode for the channel flows of Alben et al. (2024). In figure 6(f1), at $\tau$ Pe = 10 $^1$ , we have the aforementioned transition to different vorticity patterns – here, a pointlike vortex at the middle of the domain. At $\tau$ Pe = 10 $^{1.5}$ , vortices fill the domain both horizontally and vertically in the cosine component (figure 6 g1), with a rather different and weaker global vorticity pattern in the sine component (figure 6 g1s). The larger vortices of figure 6(g1,g1s) mix the boundary layer fluid more strongly with fluid in the bulk.

The sequence is generally similar at Pe = 10 $^{4.5}$ , shown in the third column, figure 6(h–k). Vortex chains occur in figure 6(h–i), at smaller $\tau$ Pe, and then a pointlike vortex in figure 6(j1) is followed by domain-filling vortices in figure 6(k1,k1s). This pattern is seen again at Pe = 10 $^{5}$ , the fourth column, figure 6(l–q). Vortex chains are seen in figure 6(l–n), with $\tau$ Pe $\leqslant$ 10 $^{0.5}$ . In figure 6(o–q), with $\tau$ Pe $\geqslant$ 10 $^{1}$ , we have vertically elongated vorticity patterns together with pointlike vortices in figure 6(o5). In figure 6(o–q), the vertically elongated vorticity is shifted from the steady convection rolls by one-quarter of the horizontal period. These unsteady vortices thus transport fluid from one steady convection roll to the other. We note that the change of eigenmodes from vortex chains to other vorticity distributions at $\tau$ Pe = 10 $^{1}$ for all Pe in figure 6 coincides with the local minimum of the top two quartets of eigenvalues in figure 5(d–f).

Figure 7. Perturbations in temperature fields ( $T_{2s}$ ) and heat flux distributions ( $-\partial _yT_{2s}$ ) at leading (quadratic) order for the flow modes in figure 6. The black lines in each panel are graphs of $-\partial _yT_{2s}$ versus $x$ for the wall shown (top or bottom). These graphs are scaled vertically to fit within the panels, with the value 0 at the vertical midpoint of each panel.

In the next section, we will consider the effect of the unsteady flow perturbations on heat transfer, for perturbation amplitudes ranging from small to large. First, we examine the leading-order (quadratic) effect on the time-averaged temperature and on Nu, in the limit of small perturbations. In figure 7, the colour fields show $T_{2s}$ for the flow perturbations in figure 6, in the same regions of space. In each panel, the black line is a graph of the corresponding heat flux out of the bottom wall or into the top wall, with the value 0 at the vertical midpoint of each panel. For both top and bottom walls, enhanced heat transfer occurs where the black line lies above the midpoint, i.e. $-\partial _yT_{2s} \gt 0$ .

In figures 7(a1–b1s), the black graphs show two main peaks in the perturbation’s average heat flux from the bottom wall $-\partial _yT_{2s}$ . The peaks lie above dark blue regions in $T_{2s}$ near the lower (hot) wall, showing increased cooling there. These are approximately where the vortex chains lie in figures 6(a1–b1s). The enhancement is partly cancelled by the black curves’ troughs near the periodic boundary. The black curves have similar peaks near the vortex chains in figure 7(d1,e1,h1,i1,l1,m1,n1), though the curves differ in some respects, e.g. the troughs’ shapes. In general, the troughs correspond to regions that the vortex chains move away from, as seen in supplementary movies 1–4. The heat flux perturbations for pointlike vortices – figure 7(f1,j1,o5), all with $\tau$ Pe = 10 $^1$ – have a common form, different from those of the vortex chains. The heat flux peaks near the pointlike vortices and is nearly zero elsewhere. Figure 7(c1) also fits this pattern, though the vortices are less pointlike. For the global vorticity distributions that occur at $\tau$ Pe $\gt$ 10 $^1$ (figure 7 k1,p1,q1), the heat flux distributions typically have cusp-like peaks that lie between convection rolls of the perturbation flows.

Figure 8. Vorticity fields for the cosine components of the top eigenmodes at Pe = $10^{5.5}$ – $10^7$ , when the steady base flows are branched. From left to right, the four columns correspond to Pevalues (a–g) 10 $^{5.5}$ , (h–n) 10 $^{6}$ , (o–u) 10 $^{6.5}$ , and (v–z,A,B) 10 $^{7}$ . The seven panels in each column correspond to $\tau$ Pe = 10 $^{-0.5}$ , 10 $^{0}$ , …, 10 $^{2.5}$ , from top to bottom. In each column, the black dashed lines separate $\tau$ Pe $\leqslant 10^{0.5}$ and $\tau$ Pe $\geqslant 10^{1}$ .

We now consider the top eigenmodes at Pe $\geqslant 10^{5.5}$ , beyond the branching transition of the steady flows. The eigenvalue distributions are shown in figure 5(g–j). Now there are many positive eigenvalues at the seven $\tau$ Pe from 10 $^{-0.5}$ to 10 $^{2.5}$ . In figure 8, we present just the eigenmode of the largest eigenvalue, and just its cosine component, which is often similar to the sine component. The four columns from left to right show the top eigenmodes for Pe = 10 $^{5.5}$ , 10 $^{6}$ , 10 $^{6.5}$ , 10 $^{7}$ . The seven rows from top to bottom show the top eigenmodes for $\tau$ Pe = 10 $^{-0.5}$ , 10 $^{0}$ , …, 10 $^{2.5}$ . At $\tau$ Pe = 10 $^{-0.5}$ (top row), the modes consist of vortex chains along the top or bottom wall. The vortices become ever smaller as Pe increases, but the main difference with the previous lower-Pe figure is that the vortex chains are shorter horizontally. In particular, the vortex chains avoid the regions on the wall adjacent to the enclosed eddies between the branches. In figure 8(a), these regions occur at the centre and the two sides of the panel (seen more clearly by the black dotted streamlines in figure 8(b), which are the same in figure 8 a). In figure 8(h), the vortices are clustered on either side of the enclosed eddy on the right-hand side (seen more clearly in figure 8 i), and in figures 8(o,v), the eddies are concentrated on single flow branches next to enclosed eddies.

The second row, with $\tau$ Pe = 10 $^{0}$ , shows a feature unique to the branched flows – the vortex chains follow the streamlines off the wall, along enclosed eddies near the wall. The sine components again have vortex chains that are shifted by one-quarter spatial period, so the vortex chains move as travelling waves around the eddies (see supplementary movie 5). The vortices in the vortex chains are larger at $\tau$ Pe = 10 $^{0.5}$ (third row), and then in the fourth row, the individual vortices lose their circular shape and become spread out along the eddy, showing again that there is a transition between $\tau$ Pe = 10 $^{0.5}$ and 10 $^{1}$ , marked with dashed lines. The fifth row is similar to the fourth. In the sixth and seventh rows, the vortices are more elongated and curve around the eddies, sometimes appearing in elongated pairs of opposite-signed vorticity. These features persist in the seventh row, but strong vorticity is also concentrated in pointlike regions near the wall that are difficult to see clearly. We have shown only the top eigenmodes, but briefly comment on the qualitative features of the second to 16th eigenmodes. At smaller $\tau$ Pe (the first four rows), the subsequent modes have vortex chains like the top eigenmodes, but they are placed at different eddies or branches of the steady flows. At larger $\tau$ Pe (the last three rows), the subsequent modes have similar elongated vortices superposed on large regions of non-zero vorticity that extend far from the walls.

Figure 9. Perturbations in temperature fields ( $T_{2s}$ ) and heat flux distributions ( $-\partial _yT_{2s}$ ) at leading (quadratic) order for the flow modes in figure 8. The black lines in each panel are graphs of $-\partial _yT_{2s}$ versus $x$ for the wall shown (top or bottom). These graphs are scaled vertically to fit within the panels, with the value 0 at the vertical midpoint of each panel.

Figure 9 shows the leading perturbations in the temperature fields ( $T_{2s}$ ) and heat flux ( $-\partial _yT_{2s}$ , black curves) corresponding to the flows in figure 8. The top row shows peaks in the heat flux perturbation near the vortex chains along the wall, as before. The second and third rows also have peaks in heat flux perturbations near the vortex chains, even though the vortex chains now extend away from the walls while the heat flux occurs at the wall. The fourth and fifth rows show sharp heat flux peaks even though the corresponding vorticity in figure 8 is more spread out across the wall. However, there are faint opposite-signed vortices visible further from the wall – slightly yellow regions in figures 9(d,r,s,y), and slightly blue regions in figure 9(k,l,z) – and the heat flux peaks seem to approximately correspond to these regions, in both number and horizontal location. In the sixth and seventh rows, the heat flux peaks generally correspond to regions of strong vorticity or interfaces between opposite-signed vorticity, though the relationship is less clear, particularly in figure 9(u,B).

Now that we have described the eigenvalue distributions and some of the top eigenmodes and the resulting heat flux patterns, we have some understanding of the optimal flow perturbations and how they increase heat flux. The flows share some features of the optimal perturbations of channel flows, such as vortex chains, but with additional complexity due to the structure of the steady optimal flows here, i.e. thin rectangular convection rolls with branching. The resulting heat flux distributions have peaks and troughs that generally correspond to features of the vorticity patterns.

Next, we show the effect of the unsteady flow perturbations on the Nusselt number as we vary the perturbation magnitude $\epsilon$ from small to large. This requires a more capable unsteady advection–diffusion solver than the one used for channel flows in Alben et al. (2024), as we explain next.

6. Unsteady simulations from small to large amplitude

We now examine the performance of the optimal eigenmodes in the small-to-large amplitude regime. We solve the unsteady advection–diffusion (2.1) for flows of the form (3.1), with $\psi _0$ set to one of the steady optima in figure 2, ${f}_{\! A}$ and ${f}_{\! B}$ corresponding to one of the top eigenmodes, $c_1 = 1$ , and the amplitude parameter $\epsilon$ ranging from $10^{-3}$ to $10^{-0.25}$ . At the upper end of this range, the unsteady perturbation is large, with energy close to that of the steady flow. We have cross-checked the simulations and eigenvalue computations, verifying that for small $\epsilon$ , Nu $\approx$ Nu $_0$ + $\lambda _i \epsilon ^2$ /2, in agreement with the expansion (4.4), when $\boldsymbol{a}$ is the Hessian eigenvector with norm $\epsilon$ corresponding to eigenvalue $\lambda _i$ . We primarily use the simulations to determine the range of $\epsilon$ for which the perturbations are effective, and characterise the resulting temperature fields.

6.1. Numerical method

On the finest grid ( $256\times1024$ ), the steady version of the discretised advection–diffusion equation can be solved in a few seconds using a direct solver (i.e. ‘backslash’ or mldivide in Matlab). In Alben et al. (2024), we solved the unsteady advection–diffusion equation with Crank–Nicolson time stepping for 20 periods of the unsteady flow, by which time the unsteady flows had usually reached a periodic steady state. With 200 time steps per period, 4000 time steps were required, with a cost equivalent to approximately 4000 steady solves. Unlike the channel flows, the wall-to-wall flows are still quite far from the periodic state after 20 periods, unfortunately, particularly at Pe = 10 $^6$ and 10 $^7$ , but also at smaller Pe such as 10 $^4$ . The difference may be due to underlying steady flows. In Alben et al. (2024), the steady flows are unidirectional flows through a channel, unlike the thin convection rolls of the wall-to-wall problem that have fine eddies and branched structures. For the channel flows, fluid parcels with initial deviations from the periodic steady-state temperature solution are rapidly advected out of the channel, while in the wall-to-wall flows, they circulate within the flow domain and may cause more extensive and long-lasting temperature deviations throughout the domain.

Because of the very slow convergence of the time-stepping approach, we instead solve the larger linear system for the temperature field that is fully coupled in time over the flow period. A space–time multigrid method has been developed for the time-periodic diffusion equation (Hackbusch 1981) and used to solve an advection–diffusion problem for channel flow at Pe = O(10 $^2$ –10 $^4$ ) (Wang 2014). The latter problem had a smoother flow field, which allowed for a uniform spatial grid. Such methods can display poor convergence for strongly advection-dominated problems (De Sterck et al. 2025), and are complicated to implement with non-uniform grids, so we adopt a simpler approach here.

First, we consider the Crank–Nicolson second-order finite-difference time discretisation of the advection–diffusion equation (2.1):

(6.1) \begin{align} &\left (1 + \frac {\Delta t}{2}\left (u\,\partial _x+v\,\partial _y-{\nabla} ^2\right)\right)T^{j+1} \notag\\ &\qquad = \left (1 - \frac {\Delta t}{2}\left (u\,\partial _x+v\,\partial _y-{\nabla} ^2\right)\right)T^{j}, \quad j = 1, \ldots , n_t, \end{align}

where $T^j(x,y)$ is the temperature field at time step $j$ , which ranges over the $n_t$ time steps of a period. Periodicity implies $T^{n_t+1}\equiv T^{1}$ , resulting in a closed system of equations for $\{T^{1},\ldots , T^{n_t}\}$ . With the unknowns ordered first with increasing $x$ at a given $y$ and $t$ , then by increasing $y$ at a given $t$ , then by $t$ , the sparsity pattern of the linear system is shown in figure 10(a). Each block row and block column corresponds to the equations and unknowns at certain time steps, respectively.

Figure 10. Matrix sparsity patterns for (a) Crank–Nicolson discretisation and (b) the time-spectral method.

The second approach that we consider is the time-spectral method. We expand the temperature field as

(6.2) \begin{align} T(x,y,t) = 1-y+\sum _{k = 0}^{N_t} A_k(x,y)\cos \left (\frac {2\pi k t}{\tau }\right) + \sum _{k = 1}^{N_t} B_k(x,y)\sin \left (\frac {2\pi k t}{\tau }\right). \end{align}

The velocity field corresponding to (3.1) is of the form

(6.3) \begin{align} u(x,y,t) &= u_{St}(x,y) + u_A(x,y) \cos \left (\frac {2\pi t}{\tau }\right) + u_B(x,y) \sin \left (\frac {2\pi t}{\tau }\right), \nonumber \\ v(x,y,t) &= v_{St}(x,y) + v_A(x,y) \cos \left (\frac {2\pi t}{\tau }\right) + v_B(x,y) \sin \left (\frac {2\pi t}{\tau }\right). \end{align}

We insert (6.2) and (6.3) into (2.1), and set to zero the inner products of the residual with cosines and sines with frequencies up to $2\pi N_t/\tau$ (i.e. the Galerkin method). When the unknowns are ordered $\{A_0, A_1, B_1, \ldots , A_{N_t}, B_{N_t}\}$ , with the corresponding ordering of the equations for the modes, the linear system is banded, with the sparsity pattern shown in figure 10(b). The products of sines and cosines in $u\,\partial _xT$ and $v\,\partial _yT$ introduce a coupling of $\{A_k,B_k\}$ with $\{A_{k-1},B_{k-1}\}$ and $\{A_{k+1},B_{k+1}\}$ . The couplings between $A_k$ and $B_{k+1}$ and $B_k$ and $A_{k-1}$ result in the blocks farthest from the main diagonal, on the third block diagonals above and below it respectively, giving 7 diagonals of blocks.

In table 1(a) and table 1(b) we show the peak memory (RAM) used and the total run time to solve the linear systems in figures 10(a) and 10(b), respectively, using Matlab’s backslash. We use $m = n = 256$ , the coarsest spatial mesh, which is used at the smallest Pe, $\leqslant 10^4$ . The cost is much higher when $n$ is increased to 512 and 1024. In the last line of each table, the peak RAM is starred because the solver was close to the maximum RAM available, therefore virtual memory, a much slower form of memory, was used. From this point, the run times grow much more rapidly.

Table 1. Peak RAM usage and total run time to solve (a) the Crank–Nicolson linear system (figure 10 a) for various $n_t$ , and (b) the time-spectral linear system (figure 10 b) for various $N_t$ . The spatial grid parameters are $m = n = 256$ .

Table 1 shows that both methods for the direct solution of the unsteady system require large memory, even for a small number of time steps per period $n_t$ (for Crank–Nicolson) and a small number of temporal frequencies $N_t$ (for the time-spectral method). This problem is well suited to the time-spectral method because the advection–diffusion equation (2.1) has a very smooth time dependence, through the flow velocities only, which have just the single frequency $2\pi /\tau$ . In such cases, exponential decay of $A_k$ and $B_k$ with $k$ in (6.2) is typical (Boyd 2001). Furthermore, in the small- $\epsilon$ expansion for $T$ (3.3), the frequency $4\pi /\tau$ does not appear until the $\epsilon ^2T_{2u}$ term, and higher frequencies appear at still higher powers of $\epsilon$ . Hence one may hope for good accuracy with a small value of $N_t$ .

We solve the time-spectral system with $N_t = 10$ , i.e. 21 temporal modes, which is still too large for a direct solution, considering that we take $n = 1024$ for Pe $\geqslant 10^{5.5}$ and we solve at $O(10^3)$ parameter combinations, varying Pe, $\tau$ Pe, the eigenmode number, and $\epsilon$ simultaneously. Therefore, we use the GMRES iterative method, which requires a good preconditioner since the system is ill-conditioned. For $m = 256$ , $n = 256$ –1024, and various Pe and $\tau$ Pe used here, the condition numbers of just the upper left block (which is the steady advection–diffusion matrix) and the upper left $3\times3$ array of blocks are $10^8$ – $10^{10}$ .

We use a block-Gauss–Seidel preconditioner that approximates the matrix by its lower triangular and diagonal blocks, and then solves for the modes at each frequency sequentially from low to high, i.e. we solve the first block row of equations for $A_0$ assuming $A_1 = B_1 = 0$ , then solve the next two block rows of equations simultaneously for $A_1$ and $B_1$ using the just-computed $A_0$ , assuming $A_2 = B_2 = 0$ , and so on. This approximation is particularly good at small $\epsilon$ , where continuing the asymptotic calculation in § 3 shows that the frequency- $j$ terms are $O(\epsilon ^j)$ .

6.2. Numerical results

We run the iterative solver for $T$ at eight values of Pe (10 $^{3.5}$ , 10 $^{4}$ , …, 10 $^{7}$ ), eight values of $\tau$ Pe (10 $^{-1}$ , 10 $^{-0.5}$ , …, 10 $^{2.5}$ ), four mode numbers (1, 5, 9, 13), and six values of $\epsilon$ (10 $^{-3}$ , 10 $^{-2}$ , 10 $^{-1}$ , 10 $^{-0.75}$ , 10 $^{-0.5}$ , 10 $^{-0.25}$ ). The mode numbers are spaced apart by four to avoid repetition of flows, as the top eigenvalues and eigenmodes often occur in nearly identical quartets.

Figure 11. (a–c) Number of iterations needed for the unsteady solver to converge to a relative residual of 10 $^{-10}$ and (d–f) norm of highest wavenumber modes $\|(A_{10},B_{10})\|_2$ relative to the norm of the time-constant mode $\|A_0\|_2$ , at three $\epsilon$ values: (a,d) $10^{-2}$ , (b,e) $10^{-1}$ , and (c,f) 10 $^{-0.5}$ . The four squares clustered at each (Pe, $\tau$ Pe) pair give the values for the four modes 1, 5, 9, 13, at the upper left, upper right, lower left, and lower right squares in each cluster, respectively.

In some cases, particularly at higher Pe and higher $\epsilon$ , the solver fails to converge within a few days of runtime, corresponding to approximately 700 iterations of preconditioned GMRES. The criterion for convergence is that the norm of the residual vector relative to the norm of the right-hand-side vector for the preconditioned system reaches 10 $^{-10}$ . The effect of parameters on convergence is shown in figure 11(a–c), which show the number of iterations needed to converge. The perturbation amplitudes are $\epsilon = 10^{-2}$ , 10 $^{-1}$ and 10 $^{-0.5}$ . The coloured squares are missing where the iteration did not converge within a few days, most often at large Pe, where the condition number is larger, and at large $\epsilon$ , where the unsteady components of the solution are more significant and the Fourier components $A_k$ and $B_k$ decay more slowly with $k$ . This is shown figure 11(d–f), which plots the norm of the highest-frequency components of the temperature, $(A_{10}, B_{10})$ arranged as a vector, relative to the norm of the steady component $A_0$ . Throughout figure 11(d–f), this relative norm is generally small, $\leqslant 10^{-2}$ , but tends to be larger near the parameters where the iterations fail to converge. There is also a noticeable effect of $\tau$ Pe on convergence. The largest number of failures occur at $\tau$ Pe = $10^1$ and adjacent values, which happens to be where the flows transition from vortex chains to other vorticity patterns. We omit Pe = 10 $^{3.5}$ and 10 $^{4}$ on the horizontal axis because the pattern was similar to Pe = 10 $^{4.5}$ , i.e. relatively rapid convergence in most cases.

Next, we discuss the Nu values for the cases that converged. For modes with positive eigenvalues, Nu increases above the steady value for some sequence of $\epsilon$ starting from 10 $^{-3}$ , then reaches a maximum, then decreases. In some cases, the computations fail to converge before a decrease is seen. For modes with negative eigenvalues, we find that Nu decreases below the steady value as expected.

Figure 12. Properties of the unsteady perturbations that achieve the largest improvement in Nu over the steady optima, across ranges of Pe and $\tau$ Pe: (a) ${\textit{Nu}}_{\textit{rel}}$ , the maximum factor of increase in Nu relative to the steady optimum at the same Pe (see (6.4)); (b) the value of $\epsilon$ at which the maximum increase in Nu occurs; (c) which of the four tested modes achieves the maximum increase in Nu.

Figure 12 shows values of Nu for the unsteady flows relative to those of the corresponding steady optima, i.e.

(6.4) \begin{align} {\textit{Nu}}_{\textit{rel}} (\tau , \textit{Pe},\epsilon , \text{mode no.})\equiv {\textit{Nu}}(\tau , \textit{Pe},\epsilon , \text{mode no.})/{\textit{Nu}}_0(\textit{Pe}). \end{align}

Figure 12(a) shows the maximum value of Nu $_{\textit{rel}}$ with respect to $\epsilon$ and mode number, at each Pe and $\tau$ Pe, among the cases that computed successfully. The increases in Nu $_{\textit{rel}}$ above 1 (the steady value) were very small at Pe = 10 $^{3.5}$ , 10 $^{4}$ , so we restrict to larger Pe on the horizontal axis. The colour of each box gives the maximum Nu $_{\textit{rel}}$ over the cases that computed successfully out of the four mode numbers and six $\epsilon$ values at that (Pe, $\tau$ Pe) pair. The maximum Nu $_{\textit{rel}}$ of 1.07 occurs at Pe = 10 $^5$ and $\tau$ Pe = 10 $^0$ , and is significantly smaller than the maximum value 1.56 seen for the channel flows (Alben et al. 2024). Since the steady channel flows were constrained to be unidirectional, this difference would be reduced by allowing for fully 2-D steady channel flows. Since the geometries and boundary conditions are different, some difference is expected in any case. A Nu $_{\textit{rel}}$ value of 1.03 or higher mainly occurs in the band of small-to-moderate $\tau$ Pe = 10 $^{-0.5}$ –10 $^{1}$ , where the flows are a mixture of vortex chains and global vorticity patterns. A selection of cases with the largest Nu $_{\textit{rel}}$ are numbered 1–16, and the vorticity and temperature fields for these cases are shown in the supplementary movies with the same number. In some coloured boxes, mainly at larger Pe, a black square outline is shown, indicating that Nu $_{\textit{rel}}$ occurred at the largest $\epsilon$ that computed successfully, so somewhat larger Nu $_{\textit{rel}}$ might occur at larger $\epsilon$ . For channel flows, by contrast, the largest Nu $_{\textit{rel}}$ increased more strongly with Pe, and lack of convergence at large Pe was less significant.

Figure 12(b) shows $\epsilon _{{max}}$ , the $\epsilon$ value at which the maximum Nu $_{\textit{rel}}$ occurred. At $\tau$ Pe = 10 $^{-1}$ , all eigenvalues are negative, so the smallest $\epsilon$ is best. Just above this $\tau$ Pe, a fairly large $\epsilon$ of 10 $^{-0.5}$ is often best, though Nu $_{\textit{rel}}$ is not much above 1. At $\tau$ Pe $\geqslant 10^1$ , smaller $\epsilon$ values, 10 $^{-2}$ and 10 $^{-1}$ , are often best, though many cases at larger $\epsilon$ did not converge. Figure 12(c ) shows the mode numbers that achieved the largest Nu $_{\textit{rel}}$ . At smaller Pe, the first mode is best. At larger Pe, the fifth mode is sometimes best at larger $\tau$ Pe, while the ninth mode is sometimes best at smaller $\tau$ Pe. As occurred for the channel flows, the perturbation with the largest eigenvalue is not always best at moderately large $\epsilon$ .

Figure 13. Snapshots of temperature fields corresponding to the flows with the largest Nu values at six (Pe, $\tau$ Pe) combinations. Panels a–b, c–d, e–f, g–h, i–j, k–l correspond to supplementary movies 3, 6, 11, 12, 15, 16, respectively, and to the cases in figure 12(a) labelled with the same numbers. The black dotted lines show the instantaneous streamlines.

To illustrate the effect of the flow perturbations at finite amplitude, in figure 13 we show snapshots of the temperature fields for six of the cases in figure 12 – those labelled 3, 6, 11, 12, 15, 16, shown in panels a–b, c–d, e–f, g–h, i–j, k–l, respectively. Each pair of panels shows the temperature field at integer and half-integer multiples of the period $\tau$ , i.e. at opposite phases. The superposed black dotted lines show the instantaneous streamlines. The corresponding supplementary movies (with the same numbers as in figure 12) give a clearer sense of the motion, but the panels illustrate the basic features.

Figure 13(a,b) show the temperature field for the top mode at Pe = 10 $^5$ and $\tau$ Pe = 10 $^0$ , which consists of two arrays of vortices moving along the bottom wall (see supplementary movie 3). Each array moves towards the hot plume that emanates on the right-hand side near $x = 0.12$ , one from the left, and the other from the right (across the periodic boundary). The vortices create an array of small hot plumes along the bottom wall that grow and curl as they move towards the larger plume. The temperature field for the corresponding steady flow is shown in figure 3(g). Figures 13(c,d) show the temperature field for the top mode at Pe = 10 $^{5.5}$ and $\tau$ Pe = 10 $^{0.5}$ , for which the steady flow is branched. The unsteady perturbation consists of chains of vortices that move around the eddies enclosed by the branches on the top and bottom walls (see supplementary movie 6). As a result, the three steady temperature plumes on the top and bottom walls in figure 3(h), two small and one large, oscillate from side to side in figure 13(c,d). A similar side-to-side oscillation occurs in figures 13(e,f), which correspond to the top mode at Pe = 10 $^{6}$ and $\tau$ Pe = 10 $^{1}$ . However, the flow mode is quite different from the previous one. Rather than a chain of vortices, we have one or two patches of vorticity that fill the entire horizontal width of the domain near the walls (see supplementary movie 11). These patches are zones of leftward and rightward shear flows along the walls.

Figures 13(g,h) are for the top mode at the next higher Pe, 10 $^{6.5}$ , but a smaller $\tau$ Pe, 10 $^{-0.5}$ . Here, two chains of tiny vortices move along portions of the bottom wall, approaching two of the three larger hot plumes (shown for the steady case in figure 3 j) from the left. The resulting temperature fields resemble figure 13(a,b). The last two pairs of panels, figure 13(i,j) and 13(k,l), correspond to the largest Pe, 10 $^{7}$ , and two moderate $\tau$ Pe, 10 $^{0.5}$ and 10 $^{1}$ . In the steady case, figure 3(k), there are several plumes of different sizes on the top and bottom walls. The flow perturbations for figure 13(i,j) and 13(k,l), shown in supplementary movies 15 and 16, respectively, are qualitatively similar. Vortices are predominantly at one wall – the top for figure 13(i,j), and the bottom for figure 13(k,l) – and alternate between small patches that move along the steady streamline branches, and larger regions of vorticity that cover the horizontal range near the wall. As previously, the most visible effect is the oscillation and in some cases pinch-off of the plumes. The net increase in Nu of 3–7 % relative to the steady case corresponds to an average thinning of the temperature boundary layers along the wall by the same amount. This relatively small thinning generally corresponds to the vortices moving colder fluid from the outer region towards the wall.

Some features of the temperature fields in figure 13 are surprisingly similar to those that have been observed for natural convection in porous media (Hewitt, Neufeld & Lister 2012; Wen, Corson & Chini 2015). At sufficiently large Rayleigh number, the unsteady temperature fields consist of alternating tall thin columns of hot and cold fluid, similar to the plumes in figure 3 that extend from one wall to the other. Closer to the walls, small fluctuating plumes are seen, like those in figure 13 and the associated supplementary movies. In Hewitt et al. (2012), e.g. the supplementary movie, these small plumes are advected along the wall towards the centres of the large plumes, as in supplementary movies 3 and 12 of the present work. We do not know how close the present computed optima are to the global optima at each Pe, but it would be interesting if they had the same structure as those that arise naturally due to buoyancy in porous media. The present optimisation framework does not include the momentum equation with or without a porous medium, but the computed optima are less similar to large-Rayleigh-number flows for natural convection in fluid without a porous medium, which do not possess the same columnar structure (van der Poel et al. 2012).

7. Conclusions

We have calculated optimal unsteady flow perturbations for increasing the rate of heat transfer Nu from the values for the best currently known steady flows for the wall-to-wall problem, across a range of Pe (square root of flow power) from $10^2$ to $10^7$ . The optimal perturbations are the leading eigenvectors of the Hessian matrix of Nu, the matrix of second derivatives with respect to amplitudes of perturbation modes. The resulting changes in Nu are given by the corresponding eigenvalues, which become positive starting at a Pe value between $10^{3.5}$ and $10^{4}$ . In this regime, unsteady perturbations can outperform the steady flows.

The same perturbation method was used previously in the case of optimal unidirectional flows through a heated channel. The eigenvalue distributions and eigenvector flow patterns have some similarities in the two cases. In both cases, eigenvalues first become positive above a critical Pe near $10^3$ – $10^4$ . The number of positive eigenvalues increases from 1–2 to tens and hundreds as Pe increases to $10^5$ – $10^7$ . Positive eigenvalues occur in bands of flow periods $\tau$ that scale as Pe $^{-1}$ , i.e. $\tau$ Pe $\sim O(1)$ . For the channel flows, the eigenvalues have a single maximum with respect to $\tau$ Pe, at 10 $^{0}$ –10 $^{0.5}$ . For the wall-to-wall flows, there are two maxima, one at $\tau$ Pe = 10 $^{0}$ –10 $^{0.5}$ , and the other at $\tau$ Pe = 10 $^{1.5}$ , which coalesce into a single maximum at Pe $\geqslant$ 10 $^{5.5}$ .

For the channel flows, the corresponding eigenmodes are chains of vortices of alternating sign that move with the steady flow as travelling waves. At Pe $\leqslant 10^5$ , where the steady wall-to-wall flows are rectangular convection rolls, and at $\tau$ Pe $\leqslant$ 10 $^{0.5}$ , the wall-to-wall eigenmodes are also vortex chains near a wall, one per convection roll. At $\tau$ Pe = 10 $^{1}$ , the wall-to-wall eigenmodes are instead pointlike vortices. At larger $\tau$ Pe, they become larger and more complicated vorticity distributions that spread vertically into the bulk flow.

At Pe $\geqslant 10^{5.5}$ , where the steady wall-to-wall convection rolls undergo branching at the walls, the top eigenmodes are vortex chains along the walls at $\tau$ Pe = 10 $^{-0.5}$ , and then along the eddies enclosed by the branches near the walls at larger $\tau$ Pe. They again become more complex at $\tau$ Pe $\geqslant$ 10 $^{1}$ , with a vorticity pattern that has some correlation with the branching pattern of the underlying steady flow.

At leading (quadratic) order, the time-averaged wall heat flux distributions show peaks near the vortex chains and other vorticity structures near the walls.

We also solved for the temperature and Nu when the unsteady flow perturbations have moderate-to-large amplitudes. Unlike for the channel flows, a time-stepping method is inefficient for attaining the time-periodic solution to the unsteady advection–diffusion equation in the wall-to-wall case. Therefore, we solved the equation in time-spectral form using the GMRES method with a block-Gauss–Seidel preconditioner. Convergence was rapid at smaller perturbation amplitudes ( $\epsilon$ ) and smaller Pe, but often did not occur at larger values. Nonetheless, for large ranges of Pe and $\tau$ Pe, we were able to find an interior maximum of Nu within the discrete set of $\epsilon$ used. The peak Nu ranged up to 7 % higher than that of the steady optimal flows. Due to computational time constraints, the search was limited to a somewhat sparse set of parameters. Given the sensitivities of the flows and temperature fields to small changes in parameters, it is possible that significantly larger increases in Nu could be found at nearby parameters. The 7 % increase occurred at Pe = 10 $^5$ with a relatively large $\epsilon$ of $10^{-0.5}$ , while at larger Pe the peak enhancement often occurs at a much smaller $\epsilon$ of $10^{-2}$ , and drops off sharply when $\epsilon$ is increased to $10^{-1}$ . Therefore, the enhancement is often quite sensitive to $\epsilon$ . A finer grid of $\epsilon$ values, with multiple values per order of magnitude, seems likely to turn up some larger increases.

The present unsteady optima were computed as perturbations of the best currently known steady optima of Alben (2023), which extended the computations of Souza et al. (2020) to Pe beyond the branching transition. The branching pattern had been found previously in three dimensions by Motoki et al. (2018a ) and analytically in two dimensions by Tobasco & Doering (2017). It would be interesting if unsteadiness could also be used to improve Nu in analytical branching flows like those of Tobasco & Doering (2017), with flow patterns either similar to or different from those in this paper.