2. Boundary layer thickness control

The aim of this section is to derive explicit expressions for the forcing amplitude $f(t)$ in (1.8) under which either of the boundary layer thicknesses $\delta ^\ast (t)\equiv 1$ or $\theta (t) \equiv 1$ can be maintained for solutions to (1.4) which satisfy the periodic boundary conditions (1.5). Figure 1 shows a schematic overview of how such a forcing $f(t)$ is defined. It consists of two components

(2.1) \begin{equation} f = K( e) + F( U , \langle uv \rangle ), \end{equation}

where $K$ and $F$ are nonlinear functionals, both of which will be defined subsequently.

Figure 1. A schematic overview of the proposed boundary layer thickness control scheme. The idea is to force either $\delta ^\ast$ or $\theta$ to converge to a reference value of $1$ , by choice of the nonlinear control laws $K$ and $F$ . The symbol $\varSigma$ denotes the summation of two signals in the feedback loop.

The first term, $K$ , is an error feedback controller that uses the deviation $e=e_\delta =1-\delta ^\ast$ (or $e = e_\theta = 1- \theta$ ) from the desired unity value of boundary layer thickness. The second term, $F$ , which depends only on the mean streamwise velocity $U$ and the stresses $\langle uv\rangle$ , is a feed-forward control term that will be constructed from analysis of the governing equations. Importantly, this will allow us to identify an explicit formula for the forcing $f$ once the flow is in a statistically steady state.

2.1. Displacement thickness control

To determine appropriate choices of the controllers $K$ and $F$ , first let $e_\delta (t) = 1-\delta ^\ast (t)$ . It is shown in Appendix A that

(2.2) \begin{equation} \frac {{\textrm{d}} e_\delta }{{\textrm{d}}t} = f(t)\delta ^\ast (t) - \frac {1}{2} C_{\!f}(t), \end{equation}

where we define the instantaneous streamwise- and spanwise-averaged skin friction by

(2.3) \begin{equation} C_{\!f}(t) := \frac {2}{\textit{Re}} \frac {\partial U}{\partial y}(0,t), \qquad t \geqslant 0. \end{equation}

Given (2.2), it is not difficult to see that asymptotic decay, i.e. $e_\delta (t) \rightarrow 0$ , of the error can be ensured by defining the feedback and feed-forward components of the forcing amplitude $f(t)$ by

(2.4) \begin{equation} K(e) = -\frac {ke}{1-e} = - \frac {ke_\delta }{\delta ^\ast }, \end{equation}

and

(2.5) \begin{equation} F(U) = \frac {\frac { \partial U}{\partial y}(0,t)}{\textit{Re} \int _0^\infty (1-U)\, {\textrm{d}}y} = \frac {C_{\!f}}{2\delta ^\ast }, \end{equation}

respectively, where $k\gt 0$ is a constant gain parameter. This observation gives the following result.

Lemma 1. Let $e_\delta = 1-\delta ^\ast$ . Suppose that the forcing amplitude ( 1.8 ) is chosen such that

(2.6) \begin{equation} f(t) := \frac {-k e_\delta (t) + \frac 12 C_{\!f}(t)}{\delta ^{*}(t)}, \end{equation}

for some $k \gt 0$ . Then, for any solution to ( 1.4 ) with boundary conditions ( 1.5 ) the displacement thickness satisfies $\delta ^{*}(t) \rightarrow 1$ as $t \rightarrow \infty$ .

An important corollary of Lemma 1 is that

(2.7) \begin{equation} \lim _{t\rightarrow \infty } \left |f(t) - \frac 12 C_{\!f}(t) \right | =0, \end{equation}

which implies that once the flow has reached a statistically steady state it must satisfy the PDE

(2.8) \begin{align} \begin{split} \frac {\partial \boldsymbol{u}}{\partial t}+(\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u} + \boldsymbol{\nabla }p &= \frac {1}{\textit{Re}} \left ( \Delta \boldsymbol{u} + y \frac {\partial \boldsymbol{u}}{\partial y}\frac {\partial U}{\partial y}(0,t) \right )\!, \\ \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &= 0, \end{split} \end{align}

subject to the boundary conditions (1.5).

The fact that (2.8) is an explicit form of the governing equations for a periodic boundary layer flow and will be useful later in § 3.2 when discussing the dependence of the resulting flow on the forcing trajectory $f(t)$ . In contrast, the numerical scheme first proposed in Biau (Reference Biau2023) defines the forcing $f(t)$ implicitly, via the solution to an optimisation problem, and this makes a theoretical analysis of the governing equations more challenging. Nonetheless, we will show in § 5 that a numerical implementation of (1.4) with $f(t)$ defined by either (2.6), or by using Biau’s implicit method, give rise to the same flows.

2.2. Momentum thickness control

Momentum thickness can be controlled in an analogous manner to the displacement thickness. Letting $e_\theta (t) = 1-\theta (t)$ , it is shown in Appendix A that

(2.9) \begin{equation} \frac {{\textrm{d}}e_\theta }{{\textrm{d}}t} = f(t) \theta (t) + \frac 12 C_{\!f}(t) - \frac {2}{\textit{Re}} \int _0^\infty \left ( \frac { \partial U}{\partial y}\right )^2 {\textrm{d}}y + 2 \int _0^\infty \langle uv \rangle \frac {\partial U}{\partial y} \, {\textrm{d}}y. \end{equation}

In this case, the error feedback is defined by

(2.10) \begin{equation} K(e_\theta ) = \frac {-ke_\theta }{1-e_\theta } = \frac {-ke_\theta }{\theta }, \end{equation}

and the feed-forward term given by

(2.11) \begin{equation} F(U,\langle uv \rangle ) = \frac { \frac {1}{\textit{Re}} \frac {\partial U}{\partial y}(0) - \frac {2}{\textit{Re}} \int _0^\infty \left ( \frac { \partial U}{\partial y}\right )^2 {\textrm{d}}y + 2 \int _0^\infty \langle uv \rangle \frac {\partial U}{\partial y} \, {\textrm{d}}y}{ \int _0^\infty U(1- U) {\textrm{d}}y}. \end{equation}

Letting $f(t) = K(e) + F(U,\langle uv \rangle )$ and substituting into (2.9) then gives the following result.

Lemma 2. Let $e_\theta (t) = 1-\theta (t)$ . Suppose that the forcing amplitude ( 1.8 ) is chosen such that

(2.12) \begin{equation} f(t) := \frac {-ke_\theta (t) - \frac 12 C_{\!f}(t) + \frac {2}{\textit{Re}} \int _0^\infty \left ( \frac { \partial U}{\partial y}\right )^2 {\mathrm{d}}y - 2 \int _0^\infty \langle uv \rangle \frac {\partial U}{\partial y} \, {\mathrm{d}}y}{\theta (t)}, \end{equation}

for some constant $k\gt 0$ . Then for any solution to ( 1.4 ) satisfying the boundary conditions ( 1.5 ), the momentum thickness satisfies $\theta (t) \rightarrow 1$ as $t \rightarrow \infty$ .

A consequence of Lemma 2 is that an explicit formula can be found for the asymptotic forcing

(2.13) \begin{equation} \lim _{t \rightarrow \infty } \left | f(t) - \left [ \frac 12 C_{\!f}(t) - 2 \int _0^\infty \frac {1}{\textit{Re}} \left ( \frac {\partial U}{\partial y} \right )^2 - \frac {\partial U}{\partial y} \langle uv \rangle \, {\textrm{d}}y \right ] \right | =0. \end{equation}

This result will be supported by numerical evidence presented in § 5, that confirms the asymptotic limit if $f(t)$ is given either by (2.12), or computed implicitly using the method of Biau (Reference Biau2023).

3. Rigorous bounds on skin friction

We now address the question of whether solutions to (1.4), (1.5) on the periodic domain $\varOmega$ have comparable statistical properties to those of a canonical, spatially evolving, boundary layer. In particular, our aim is to prove rigorous upper bounds on the time-averaged value of skin friction

(3.1) \begin{equation} C_{\!f} := \limsup _{t \rightarrow \infty } \frac {1}{T} \int _0^T C_{\!f}(t)\, {\textrm{d}}t, \end{equation}

and to understand how these bounds scale with $(\textit{Re})$ . It is important to note that, even though the results of § 2 ensure that the boundary layer thickness is kept constant, there is no a priori guarantee that the underlying fluid velocity field itself, or its statistics such as $C_{\!f}$ , are well behaved. In other words, it is not clear that controlling a property of the flow’s outer layer will necessarily give inner layer behaviour that resembles a section of a canonical, spatially evolving, boundary layer. We now investigate this question from a theoretical perspective, before presenting numerical evidence in § 5.

In this section, we consider solutions to (1.4), (1.5) with the forcing $f(t)$ given by (2.6). This ensures that the displacement thickness satisfies $\delta ^\ast (t) \rightarrow 1$ as $t \rightarrow \infty$ . Consequently, the Reynolds number should be interpreted as

(3.2) \begin{equation} \textit{Re} = \textit{Re}_{\delta ^\ast } = \frac {U_\infty \delta ^\ast }{\nu } = \frac {U_\infty }{\nu }. \end{equation}

Throughout this section, we also only consider solutions to (1.4) that satisfy both $C_{\!f}(t) \geqslant 0, t \geqslant 0$ , and the decay condition

(3.3) \begin{equation} \lim _{y \rightarrow \infty } y(\boldsymbol{e}_x-\boldsymbol{u}(x,y,z)) = 0, \qquad 0 \leqslant x \leqslant L_x, \, 0 \leqslant z \leqslant L_z, \end{equation}

which places only a weak constraint on the following analysis. A similar assumption was made in Kumar & Garaud (Reference Kumar and Garaud2020) when bounding the drag coefficient for flow past a finite-length flat plate.

We begin by taking the dot product of (1.4) with $\boldsymbol{u}-\boldsymbol{e}_x$ , integrating by parts, using the boundary conditions (1.5), and the assumption (3.3) to obtain the energy equation

(3.4) \begin{equation} \frac 12 \frac {\textrm{d}}{{\textrm{d}}t} \| \boldsymbol{u} - \boldsymbol{e}_x \|^2 + \frac {1}{{\textit{Re}}}\| \boldsymbol{\nabla }\boldsymbol{u} \|^2 = \frac 12 C_{\!f}(t) - \frac 12 f(t) \|\boldsymbol{u}-\boldsymbol{e}_x\|^2, \end{equation}

a proof of which is given in Appendix A. For clarity, in the above equation, we have used the notation

(3.5) \begin{equation} \|\boldsymbol{u}\|^2 = \int _\varOmega \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{u} \, {\textrm{d}}\boldsymbol{x}. \end{equation}

Next, consider a decomposition of the velocity field

(3.6) \begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = (1-\psi (y))\boldsymbol{e}_x + \boldsymbol{v}(\boldsymbol{x},t), \end{equation}

where $\psi :[0,\infty ) \rightarrow {\mathbb{R}}$ is any function satisfying

(3.7) \begin{equation} \psi (0)=1 \quad \text{and} \quad \lim _{y \rightarrow \infty } y\psi (y) = 0. \end{equation}

The perturbation $\boldsymbol{v} = v_x \boldsymbol{e}_x + v_y \boldsymbol{e}_y + v_z \boldsymbol{e}_z$ then inherits the periodic boundary conditions in the $\boldsymbol{e}_x,\boldsymbol{e}_z$ directions, is incompressible, and also satisfies

(3.8) \begin{equation} \boldsymbol{v}(x,0,z,t)=0 \quad \text{and} \quad \lim _{y \rightarrow \infty } y\boldsymbol{v}(x,y,z,t)=0, \end{equation}

for any fixed $x,z$ and $t$ .

Employing a flow decomposition of the form (3.6) is the base step of applying the ‘background method’ (Doering & Constantin Reference Doering and Constantin1992) for rigorous flow analysis. This decomposition is in terms of a ‘background profile’ $(1-\psi (y))\boldsymbol{e}_x$ about which perturbations $\boldsymbol{v}$ are considered. The profile $\psi$ can be arbitrarily chosen, so long as the boundary conditions (3.7) are satisfied, and its selection is a key part of the analysis.

To explain how to select $\psi$ in more detail, we begin with the energy equation

(3.9) \begin{align} \frac 12 \frac {\textrm{d}}{{\textrm{d}}t}\|\boldsymbol{v}\|^2 =- \frac {1}{\textit{Re}} \|\boldsymbol{\nabla }\boldsymbol{v}\|^2 + \int v_x v_y \psi ' {\textrm{d}}\boldsymbol{x} + \frac {1}{\textit{Re}}\int \frac {\partial v_x}{\partial y} \psi ' {\textrm{d}}\boldsymbol{x} + f(t) \int y \boldsymbol{v} \boldsymbol{\cdot }\frac {\partial \boldsymbol{u}}{\partial y} {\textrm{d}}\boldsymbol{x}, \end{align}

which is satisfied by the perturbations $\boldsymbol{v}$ . The meaning of each term in (3.9) can be made more transparent if it is simplified by using two identities. The first

(3.10) \begin{equation} \|\boldsymbol{\nabla }\boldsymbol{u}\|^2 = \|\boldsymbol{\nabla }\boldsymbol{v}\|^2 + \|\psi '\|^2 - 2 \int \frac {\partial v_x}{\partial y} \psi ' {\textrm{d}}\boldsymbol{x}, \end{equation}

follows directly from (3.6); while the second

(3.11) \begin{equation} \int y \boldsymbol{v} \boldsymbol{\cdot }\frac {\partial \boldsymbol{u}}{\partial y} {\textrm{d}}\boldsymbol{x} = -\frac 12 \|\boldsymbol{v}\|^2 - \int y v_x \psi ' {\textrm{d}}\boldsymbol{x}, \end{equation}

can be proven using (3.6) and the fact that $\int y \boldsymbol{v} \boldsymbol{\cdot }({\partial \boldsymbol{v}}/{\partial y}) {\textrm{d}}\boldsymbol{x} = - (1/2) \|\boldsymbol{v}\|^2$ which follows from integration by parts. Upon substituting (3.10) and (3.11) into (3.9), the energy equation for the perturbations $\boldsymbol{v}$ to the background profile $(1-\psi (y))\boldsymbol{e}_x$ is equivalent to

(3.12) \begin{align} &\frac 12 \frac {\textrm{d}}{{\textrm{d}}t}\|\boldsymbol{v}\|^2 + \frac {1}{2\textit{Re}}\big ( \| \boldsymbol{\nabla }\boldsymbol{v} \|^2 + \|\boldsymbol{\nabla }\boldsymbol{u}\|^2 \big ) + \frac 12 f(t) \|\boldsymbol{v}\|^2 \nonumber \\ &\qquad \qquad = \frac {1}{2\textit{Re}} \|\psi '\|^2 +\int v_x v_y \psi ' {\textrm{d}}\boldsymbol{x} - f(t) \left ( \int y v_x \psi ' {\textrm{d}}\boldsymbol{x}\right )\!. \end{align}

The final two terms on the left-hand side of (3.12) are dissipative (since $f(t)$ will be chosen to be positive), while the three terms on the right-hand side provide forcing to the energy equation. The aim of the background method is to select $\psi$ so that the forcing terms can be appropriately dominated by the dissipative terms. If this is possible then, after time averaging (3.12), one can conclude both that the perturbation energy is bounded and, furthermore, obtain quantitative bounds on the time average of the dissipative terms. Importantly, since Lemma 1 implies that $f(t) \rightarrow (1/2) C_{\!f}(t)$ , this opens the door to obtaining bounds on the skin-friction coefficient of the flow.

A key challenge in applying the background method is that the magnitude of the dissipative terms in (3.12) decreases with increasing Reynolds number. Consequently, $\psi$ must be chosen to be $\textit{Re}$ -dependent in such a way as to also make the forcing terms on the right-hand side of (3.12) decrease with increasing Reynolds number. The fact that the perturbations $\boldsymbol{v}$ are constructed to have homogeneous boundary conditions (3.8) is useful to achieve this. Since $\boldsymbol{v}$ is small near the boundaries, if $\psi '$ is constructed to be non-zero only near the boundaries, then the magnitude of the final two forcing terms in (3.12) can be controlled. Ideally, one would like to pick $\psi ' \equiv 0$ to eliminate all the forcing terms. However, this choice cannot be made since $\psi$ must satisfy the boundary conditions (3.7). These state that $\psi$ is unity at the wall and must decrease to zero far from it, which necessitates $\psi '\neq 0$ . Instead, the selection of a good background profile $\psi$ must balance the competing aims of minimising $(2\textit{Re})^{-1}\|\psi '\|^2$ (which creates energy in (3.12)), while simultaneously concentrating $\psi$ close to the wall (which necessarily increases $\psi '$ ). Optimal choices of $\psi$ , as will be constructed subsequently, then exhibit a natural $\textit{Re}$ -dependent boundary layer structure near the wall. It is important to note that the best choice of $\psi$ for proving bounds will not necessarily correspond to any physical flow profile, such as the mean streamwise velocity.

We now proceed with the technical details for constructing a background profile $\psi$ which achieves the aims of the previous heuristic discussion. To begin, we combine the energy equations (3.12) and (3.4) to obtain

(3.13) \begin{align} \frac{\textrm{d}}{{\textrm{d}}t}& \big [ 2\|\boldsymbol{v}\|^2 - \|\boldsymbol{u}-\boldsymbol{e}_x\|^2 \big ] +(C_{\!f}(t) - 2f(t)) \nonumber\\ &\qquad +\, 2f(t) \left (1 - \frac 12\|\boldsymbol{u}-\boldsymbol{e}_x\|^2 + \|\boldsymbol{v}\|^2 + 2 \int y v_x \psi ' {\textrm{d}}y \right ) + 2 Q_\psi (\boldsymbol{v}) = \frac {2}{\textit{Re}} \|\psi '\|^2, \end{align}

where

(3.14) \begin{equation} Q_\psi (\boldsymbol{v}):= \frac {1}{\textit{Re}} \|\boldsymbol{\nabla }\boldsymbol{v}\|^2 - 2\int v_x v_y \psi ' {\textrm{d}}\boldsymbol{x}, \end{equation}

is a balance between dissipation and production of perturbation energy. As discussed above, the plan is to construct $\psi$ so that the final two terms on the left-hand side of (3.13) are dissipative, i.e. such that $Q_\psi \geqslant 0$ and such that the bracketed terms multiplying the positive function $f$ are themselves positive. After making such a choice of $\psi$ , taking a time average of (3.13) will then facilitate a quantitative upper bound on the skin-friction coefficient since, by Lemma 1, we have $\overline {2f(t)} = C_{\!f}$ . To explain precisely how time averaging the identity (3.13) leads to such an upper bound, we now discuss each of the four terms on the left-hand side of (3.13) in turn.

(i) The time average of the derivative term $({\textrm{d}}/{{\textrm{d}}t}) [2\|\boldsymbol{v}\|^2-\|\boldsymbol{u}-\boldsymbol{e}_x\|^2]$ is zero. This follows because it is generally true that

(3.15) \begin{equation} \lim _{T \rightarrow \infty } \frac {1}{T} \int _0^T \frac {{\textrm{d}}g}{{\textrm{d}}t} {\textrm{d}}t = \lim _{T \rightarrow \infty } \frac {g(T)-g(0)}{T}=0, \end{equation}

for any bounded function $g$ and, using the assumption that $C_{\!f}(t) \geqslant 0$ , it follows from (3.4) that both $\|\boldsymbol{u}-\boldsymbol{e}_x\|^2$ and $\|\boldsymbol{v}\|^2$ are uniformly bounded in time.

(ii) For the second term $C_{\!f}(t)-2f(t)$ , we use Lemma 1 to infer that $f(t) \rightarrow (1/2) C_{\!f}(t)$ as $t\rightarrow \infty$ . Hence, taking the time average of the third term in (3.13) gives

(3.16) \begin{equation} \overline { C_{\!f}(t) - 2f(t)} = 0. \end{equation}

(iii) The importance of the third term

(3.17) \begin{equation} 2f(t) \left (1 - \frac 12\|\boldsymbol{u}-\boldsymbol{e}_x\|^2 + \|\boldsymbol{v}\|^2 + 2 \int y v_x \psi ' {\textrm{d}}y \right )\!, \end{equation}

is that its time average will produce a term proportional to the skin-friction coefficient since, again by Lemma 1, we have $2f(t) \rightarrow C_{\!f}(t)$ as $t \rightarrow \infty$ . However, a difficulty in analysing this expression is that it involves products of time-varying quantities. This difficulty can be avoided if the bracketed, flow-dependent, term multiplying $f(t)$ in (3.17) can be appropriately estimated. The aim is to show that $\psi$ can be picked such that this term is always positive, for any velocity field $\boldsymbol{u}$ and its corresponding perturbation $\boldsymbol{v}$ . Once achieved, this will imply that time averaging the term involving $f(t)$ will give a quantity proportional to the skin-friction coefficient.

To achieve this aim, first use the definition (3.6) of the perturbation $\boldsymbol{v}$ to obtain

(3.18) \begin{equation} \|\boldsymbol{v}\|^2 - \frac 12 \|\boldsymbol{u}-\boldsymbol{e}_x\|^2 = \frac 12 \|\boldsymbol{v}\|^2 + \int \psi v_x {\textrm{d}}\boldsymbol{x} - \frac 12 \|\psi '\|^2. \end{equation}

Hence, the bracketed term multiplying $f(t)$ in (3.17) becomes

(3.19) \begin{align} 1 - \frac 12\|\boldsymbol{u}-\boldsymbol{e}_x\|^2 + \|\boldsymbol{v}\|^2 &+ 2 \int y v_x \psi ' \, {\textrm{d}}\boldsymbol{x} \nonumber \\ &= 1 + \frac 12 \|\boldsymbol{v}\|^2 - \frac 12 \|\psi \|^2 + \int v_x(\psi + 2y \psi ')\, {\textrm{d}}\boldsymbol{x}. \end{align}

The final term in (3.19) can be estimated using the Cauchy–Schwarz inequality and the identity $2ab \leqslant a^2+b^2$ by

(3.20) \begin{equation} \left | \int v_x (\phi + 2y\psi ') {\textrm{d}}\boldsymbol{x} \right | \leqslant \| v_x \| \|\psi + 2y \phi '\| \leqslant \frac 12 \|v_x\|^2 + \frac 12 \|\psi + 2y \psi '\|^2. \end{equation}

Using this estimate and (3.19) implies that, for any velocity field $\boldsymbol{u}$ and its corresponding perturbation $\boldsymbol{v}$ ,

(3.21) \begin{align} 1 - \frac 12\|\boldsymbol{u}-\boldsymbol{e}_x\|^2 &+ \|\boldsymbol{v}\|^2 + 2 \int y v_x \psi ' \, {\textrm{d}}\boldsymbol{x} \nonumber \\ &\geqslant 1 + \frac 12\|\boldsymbol{v}\|^2 - \frac 12 \|\psi \|^2 - \left | \int v_x (\phi + 2y\psi ')\, {\textrm{d}}\boldsymbol{x} \right | \nonumber \\ &\geqslant 1 + \frac 12\|\boldsymbol{v}\|^2 - \frac 12 \|\psi \|^2 - \frac 12 \|v_x\|^2 - \frac 12 \|\psi + 2y\psi '\|^2 \nonumber \\ &=1 + \frac 12\big(\|v_y\|^2 + \|v_z\|^2 \big) - \frac 12 \big( \|\psi \|^2 + \| \psi + 2y \psi '\|^2\big) \nonumber \\ & \geqslant 1 - \frac 12 \big(\|\psi \|^2 + \| \psi + 2y \psi '\|^2\big), \end{align}

where the last inequality follows since $\|v_y\|,\|v_z\| \geqslant 0$ . For a final simplification, expand the final term in (3.21) to obtain

(3.22) \begin{equation} \|\psi + 2y\psi '\|^2 = \|\psi \|^2 + 4 \int y \psi \psi ' {\textrm{d}}\boldsymbol{x} + 4 \|y\psi '\|^2. \end{equation}

Now, using integration by parts and the boundary conditions (3.7) on $\psi$

(3.23) \begin{equation} \int _0^\infty y \psi \psi '{\textrm{d}}y = \left [ y \psi (y) \right ]_{y=0}^\infty - \|\psi \|^2 - \int _0^\infty y\psi \psi ' {\textrm{d}}y \Rightarrow \int _0^\infty y \psi \psi ' {\textrm{d}}y = -\frac 12 \|\psi '\|^2. \end{equation}

Hence, using (3.21), (3.22) and (3.23) implies that the lower bound

(3.24) \begin{equation} \left ( 1 - \frac 12\|\boldsymbol{u}-\boldsymbol{e}_x\|^2 + \|\boldsymbol{v}\|^2 + 2 \int y v_x \psi ' \, {\textrm{d}}\boldsymbol{x}\right ) \geqslant 1 - 2 \|y\psi '\|^2, \end{equation}

holds for any velocity field $\boldsymbol{u}$ and its corresponding perturbation $\boldsymbol{v}$ .

The importance of this lower bound is that if the profile $\psi$ is chosen such that $1-2\|y\psi '\|^2 \gt 0$ , we can use the fact that $f(t) \geqslant 0$ to infer that

(3.25) \begin{equation} 2f(t) \left ( 1 - \frac 12\|\boldsymbol{u}-\boldsymbol{e}_x\|^2 + \|\boldsymbol{v}\|^2 + 2 \int y v_x \psi ' \, {\textrm{d}}\boldsymbol{x}\right ) \geqslant 2f(t)(1-2\|y\psi '\|^2), \end{equation}

for any velocity field $\boldsymbol{u}$ and its corresponding perturbation $\boldsymbol{v}$ . Time averaging the above equation, and using the fact that $|f(t)- (1/2) C_{\!f}(t)| \rightarrow 0$ from Lemma 1, then gives

(3.26) \begin{equation} \overline {2f(t) \left ( 1 - \frac 12\|\boldsymbol{u}-\boldsymbol{e}_x\|^2 + \|\boldsymbol{v}\|^2 + 2 \int y v_x \psi ' \, {\textrm{d}}\boldsymbol{x}\right )} \geqslant C_{\!f} (1-2\|y\psi '\|^2), \end{equation}

where $C_{\!f} = \overline {C_{\!f}(t)}$ is the time-averaged skin-friction coefficient.

To summarise the manipulations performed so far, we have shown that time averaging (3.13) removes the first two terms on the left-hand side and, by (3.26), then gives

(3.27) \begin{equation} 2\overline {Q_\psi (\boldsymbol{v})} + C_{\!f}(1-2\|y\psi '\|^2) \leqslant \frac {2}{\textit{Re}}\|\psi '\|^2. \end{equation}

(iv) For the fourth term, suppose that it is possible to choose $\psi$ such that $Q_\psi (\boldsymbol{v})$ is positive, that is,

(3.28) \begin{equation} Q_\psi (\boldsymbol{v}) = \frac {1}{\textit{Re}} \|\boldsymbol{\nabla }\boldsymbol{v}\|^2 - 2 \int v_x v_y \psi ' {\textrm{d}}\boldsymbol{x} \geqslant 0, \end{equation}

holds for any vector field $\boldsymbol{v}$ satisfying the boundary conditions (3.8). Since any solution to the PDE (1.4) gives rise to such a field $\boldsymbol{v}$ , it would then follow from (3.27) that

(3.29) \begin{equation} C_{\!f}(1-2\|y\psi '\|^2) \leqslant \frac {2}{\textit{Re}}\|\psi '\|^2. \end{equation}

The previous four steps analysing the effect of time averaging of (3.13) imply that if $\psi$ can be picked to satisfy the three conditions:

1. (C1) $Q_\psi (\boldsymbol{v}) \geqslant 0$ for any field $\boldsymbol{v}$ satisfying the boundary conditions (3.8);

2. (C2) $1-2\|y\psi '\|^2 \gt 0$ ; and

3. (C3) $\psi$ satisfies the boundary conditions (3.7);

then the skin-friction coefficient must satisfy the upper bound

(3.30) \begin{equation} C_{\!f} \leqslant \frac {2\|\psi '\|^2}{\textit{Re}(1 - 2 \|y\psi '\|^2)}. \end{equation}

Given this observation, the challenge is to determine whether any profile $\psi$ exists which satisfies the three conditions (C1)–(C3) and, if so, to find the profile which implies the tightest possible (i.e. smallest upper) bound on $C_{\!f}$ via (3.30). Heuristically, one can observe from the definition (3.28) that if $\psi '$ is sufficiently small then (C1) should hold. The same observation applies trivially to condition (C2) and, in addition, smaller choices of $\psi '$ give better upper bounds via (3.30). However, as previously discussed, all of these three observations must work against the boundary conditions (C3) that $\psi$ must satisfy.

Theorem1, which is the main result of the paper, shows that it is possible to choose $\psi$ to satisfy the competing constraints (C1)–(C3) and, consequently, that an upper bound on skin friction for solutions to the governing (1.4) can be proven.

Theorem 1. Let $\textit{Re} \gt 1/\sqrt {2}$ . For any solution of ( 1.4 ) satisfying the boundary conditions ( 1.5 ) for which $C_{\!f}(t) \geqslant 0$ and ( 3.3 ) hold, the time-averaged skin friction satisfies

(3.31) \begin{equation} C_{\!f} \leqslant \frac {\textit{Re}}{2(\sqrt {2}\textit{Re} -1)}. \end{equation}

Proof. Let $\psi (y)$ be defined by

(3.32) \begin{equation} \psi (y) = e^{-\frac {\textit{Re}}{\sqrt {2}}y}, \qquad y \geqslant 0. \end{equation}

Then (C3) is satisfied since $\psi (0)=1$ and $y\psi (y)\rightarrow 0$ as $y \rightarrow \infty$ . To verify (C1), note first that for any differentiable function $f:[0,\infty ) \rightarrow {\mathbb{R}}$ satisfying $f(0)=0$ , the Cauchy–Schwarz inequality can be used to show that

(3.33) \begin{equation} | f(y)| = \left | \int _0^y f'(z)\, {\textrm{d}}z \right | \leqslant \sqrt {y} \left ( \int _0^y |f'(y)|^2 {\textrm{d}}y \right )^{\frac 12}. \end{equation}

Applying this pointwise estimate to each of the velocity components $v_x,v_y$ and using the definition of $\|\boldsymbol{\nabla }\boldsymbol{v}\|^2$ then gives

(3.34) \begin{equation} \left | \int v_x v_y \psi ' {\textrm{d}}\boldsymbol{x}\right | \leqslant \frac {1}{2\sqrt {2}} \|\boldsymbol{\nabla }\boldsymbol{v} \|^2 \int y |\psi '(y)|\, {\textrm{d}}y = \frac {1}{2\textit{Re}} \|\boldsymbol{\nabla }\boldsymbol{v}\|^2, \end{equation}

for any $\boldsymbol{v}$ satisfying $\boldsymbol{v}(x,0,z)=0$ . Consequently, condition (C1) is satisfied for this choice of $\psi$ . Finally, It then follows from (3.30) that

(3.35) \begin{equation} C_{\!f} \leqslant \frac {2\|\psi '\|^2}{\textit{Re}(1 - 2 \|y\psi '\|^2)} = \frac {\textit{Re} / (2\sqrt {2})}{1 -1 / (\sqrt {2}\textit{Re})} = \frac {\textit{Re}}{2(\sqrt {2}\textit{Re} -1)}. \end{equation}

An interesting corollary of Theorem1 can be obtained by retaining the perturbation energy term $E_{\boldsymbol{v}}(t):= (1/2) ( \|v_y\|^2 + \|v_z\|^2 )$ in the estimate (3.21). In Theorem1 these terms are estimated from below by zero, but one would expect that their time average satisfies $\overline {E_v} \gt 0$ for any turbulent solution to the governing equations. To make use of this observation, we introduce the notion of the covariance of two time series, $a(t)$ and $b(t)$ , letting

(3.36) \begin{equation} \text{cov}(a,b) := \lim _{T\rightarrow \infty } \frac {1}{T} \int _0^T ( a(t) - \bar {a})(b(t)-\bar {b})\, {\textrm{d}}t. \end{equation}

It follows from this definition that the time average of the product $a(t)b(t)$ can be expressed as

(3.37) \begin{equation} \overline {a b} = \bar {a}\bar {b} + \text{cov}(a,b). \end{equation}

Applying this formula when time averaging (3.13) then gives

(3.38) \begin{equation} \overline {Q_\psi (\boldsymbol{v})} + C_{\!f} \big(1 + \overline {E_{\boldsymbol{v}}} -2\|y\psi '\|^2 \big) \leqslant \frac {2\|\psi '\|^2}{\textit{Re}} - \text{cov}(C_{\!f},E_{\boldsymbol{v}}). \end{equation}

We then have the following corollary, whose proof uses the same function $\psi$ as in Theorem1.

Corollary 1. For any solution of ( 2.8 ) for which $C_{\!f}(t) \geqslant 0$ and ( 3.3 ) hold, the time-averaged skin friction satisfies

(3.39) \begin{equation} C_{\!f} \leqslant \frac {1}{2\sqrt {2}(1 + \overline {E_{\boldsymbol{v}}})} - \frac {{\textrm{cov}}(C_{\!f},E_{\boldsymbol{v}})}{1 + \overline {E_{\boldsymbol{v}}}}, \qquad \textit{Re} \gg 1. \end{equation}

3.1. Discussion of the theoretical results

One should not expect to match the value of the constant in the analytically provable bound $C_{\!f} \leqslant (2\sqrt {2})^{-1}$ with the true value of $C_{\!f}$ observed in numerical simulations, although numerically optimising $\psi$ in Constantin–Doering–Hopf bounding arguments can often improve the absolute value of such constants by an order of magnitude (see, for example, studies by Plasting & Kerswell (Reference Plasting and Kerswell2003) and Fantuzzi & Wynn (Reference Fantuzzi and Wynn2016)). However, one may instead hope to capture the correct asymptotic scaling of $C_{\!f}$ with Reynolds number. From this perspective, the result Theorem1 is more important. That $C_{\!f}$ is provably bounded, independent of the Reynolds number, is within a logarithmic correction of empirically observed scaling expected for spatially evolving boundary layers.

Another view is provided if we define a friction Reynolds number based on the displacement thickness by

(3.40) \begin{equation} \textit{Re}_\tau = \frac {U_\tau \delta ^{*}}{\nu } = \textit{Re}_{\delta ^\ast } \sqrt { \frac {C_{\!f}}{2}}, \end{equation}

where $U_\tau = ( \nu \overline {\partial U/\partial y(0,t)} )^{ {1}/{2}}$ is the friction velocity. Theorem1 then implies that

(3.41) \begin{equation} \textit{Re}_\tau \leqslant 2^{-\frac 54} \textit{Re}_{\delta ^{\ast }}. \end{equation}

This is not dissimilar to the observation of Schlatter & Orlu (Reference Schlatter and Orlu2010) that for canonical, spatially evolving, boundary layers

(3.42) \begin{equation} {\textit{Re}}_\tau \approx c \textit{Re}_{\delta ^\ast }^{0.843}, \end{equation}

where $c = 1.13 (\delta ^{\ast })^{0.157} \theta ^{0.843} \delta ^{-1} = {\mathcal{O}}(1)$ is simply a constant accounting for the various possible definitions of the boundary layer thickness, and $\delta$ is the boundary layer thickness defined in terms of $99\,\%$ of the free-stream velocity.

The exponent $0.843$ in this observation should be treated with caution, being based on an empirical fit to data in a small (in the context of understanding asymptotic scalings) range of Reynolds numbers. Indeed, the high Reynolds number scaling corresponding to the logarithmic friction law would be of the form

(3.43) \begin{equation} \textit{Re}_\tau \sim \frac {\textit{Re}_{\delta ^\ast }}{\log {\textit{Re}_{\delta ^\ast }}}, \qquad \textit{Re}_{\delta ^\ast } \gg 1. \end{equation}

Whether an analytical proof can bridge this ‘logarithmic’ gap is an important open question in theoretical fluid mechanics. It is therefore of interest that Corollary 1 provides partial evidence that the origin of such a correction can now be understood, at least for periodic boundary layers. Indeed, the improved bound

(3.44) \begin{equation} C_{\!f} \leqslant \frac {1}{2\sqrt {2}(1 + \overline {E_{\boldsymbol{v}}})} - \frac {{\textrm{cov}}(C_{\!f},E_{\boldsymbol{v}})}{1 + \overline {E_{\boldsymbol{v}}}}, \end{equation}

would imply such a logarithmic correction, if it could be proven that the

(3.45) \begin{align} \sup _{\textit{Re}_{\delta ^\ast } \gt 0} |{\textrm{cov}}(C_{\!f},E_{\boldsymbol{v}})| \lt \infty \quad \text{and} \quad \lim _{\textit{Re}_{\delta ^\ast } \rightarrow \infty } \left ( \frac {\overline {E_{\boldsymbol{v}}}}{\log {\textit{Re}_{\delta ^\ast }}} \right )\gt 0, \end{align}

that is, if the perturbations $\overline {E_{\boldsymbol{v}}}$ (which, via (3.6), are defined in terms of perturbations from the constructed flow field $(1-\psi )\boldsymbol{e}_x$ ) themselves exhibit a logarithmic growth of energy and are only weakly correlated with the instantaneous skin friction.

3.2. Flow dependence on the forcing amplitude $f(t)$

Given the definitions of the forcing functions (2.6) and (2.12) for displacement and momentum control, respectively, an interesting observation is that there are many forcing amplitude trajectories $(f(t))_{t \geqslant 0}$ that can achieve the desired asymptotic unity value of displacement or momentum thickness. Indeed, both (2.6) and (2.12) are parameterised by a gain $k\gt 0$ which gives, in each case, a whole family of forcing amplitudes which achieve asymptotic control of the boundary layer thickness. This raises the natural question of to what extent the resulting flow is dependent on the particular choice of forcing trajectory.

To begin to answer this question consider, for example, the case of displacement thickness control. It follows from (2.2) that any forcing trajectory $f(t)$ which achieves $\delta ^\ast (t) \rightarrow 1$ must necessarily satisfy

(3.46) \begin{equation} \lim _{t \rightarrow \infty } \left | f(t) - \frac 12 C_{\!f}(t)\right | =\lim _{t \rightarrow \infty } \left | f(t) - \frac {1}{\textit{Re}} \frac {\partial U}{\partial y}(0,t) \right | = 0. \end{equation}

This implies that any forcing trajectory that achieves asymptotic convergence of the displacement thickness must converge asymptotically to the same functional form (and the formula (2.6) merely gives concrete examples of some trajectories achieving this). A more subtle question to ask, therefore, is to what extent the resulting flow is dependent on the different possible transient trajectories of $f(t)$ . In other words, asking whether the nature of the journey can affect the final destination.

This is a challenging question, but we already have some evidence that supports the assertion that the eventual flow is not strongly dependent upon the transient forcing trajectory. First, the theoretical scaling law for $C_{\!f}$ proven in Theorem 1 applies equivalently to any forcing trajectory that achieves $\delta ^\ast (t) \rightarrow 1$ . More generally, suppose that the forcing transient has converged in the sense that both $|f(t) - (1/2) C_{\!f}(t)|$ and $|\delta ^\ast (t)-1|$ are negligible. It then follows from the discussion following Lemma 1 that the flow field must satisfy the PDE (2.8). Importantly, the form of this PDE is independent of the choice of forcing trajectory, since it depends only on the asymptotic functional form $(\textit{Re})^{-1} U_y(0,t)$ . Consequently, the dependence of the solution of (1.4) upon the transient forcing trajectory $f(t)$ is equivalent to considering solutions to forcing-independent PDE (2.8) initialised from different initial conditions satisfying $\delta ^\ast (0)=1$ .

It is in this subtle sense that solutions to the periodic boundary layer equations (1.4) are independent of the forcing function $f(t)$ . Of course, one would hope that (2.8) is ergodic in the sense that solutions initialised from an appropriate subset of initial conditions all have the same time-averaged statistics. However, this is a deep question in pure mathematics, which still remains open for the three-dimensional Navier–Stokes equations, and will require significant future work to resolve.

4. Numerical implementation

To perform DNS and ILES of the periodic boundary layer equations (1.4) we use the Xcompact3d framework, a suite of fluid flow solvers dedicated to the study of turbulent flows on supercomputers (Bartholomew et al. Reference Bartholomew, Deskos, Frantz, Schuch, Lamballais and Laizet2020), which has been extensively validated for wall-bounded turbulent flows (Diaz-Daniel, Laizet & Vassilicos Reference Diaz-Daniel, Laizet and Vassilicos2017; Mahfoze & Laizet Reference Mahfoze and Laizet2021; O’Connor et al. Reference O’Connor, Diessner, Wilson, Whalley, Wynn and Laizet2023). The ILES are based on a strategy that introduces targeted numerical dissipation at the small scales through the discretisation of the second derivatives of the viscous terms (Lamballais, Fortuné & Laizet Reference Lamballais, Fortuné and Laizet2011; Dairay et al. Reference Dairay, Lamballais, Laizet and Vassilicos2017). It was shown in these studies that it is possible to design a high-order finite-difference scheme in order to mimic a subgrid-scale model for ILES based on the concept of spectral vanishing viscosity, at no extra computational cost and with excellent performance for wall-bounded turbulent flows (Mahfoze & Laizet Reference Mahfoze and Laizet2021).

The incompressible flow solver within Xcompact3d is based on sixth-order compact finite-difference schemes (Laizet & Lamballais Reference Laizet and Lamballais2009) for the spatial discretisation and a fractional-step method using a semi-implicit approach that combines the Crank–Nicholson and third-order Adams–Bashforth methods. Within the fractional-step method, the incompressibility condition is dealt with by directly solving a Poisson equation in spectral space using three-dimensional (3-D) fast Fourier transforms (FFTs) and the concept of the modified wavenumber (Lele Reference Lele1992). The velocity-pressure mesh arrangement is half-staggered to avoid spurious pressure oscillations (Laizet & Lamballais Reference Laizet and Lamballais2009).

The simplicity of the mesh allows an easy implementation of a 2-D domain decomposition based on pencils (Laizet & Ning Reference Laizet and Ning2011). The computational domain is split into a number of sub-domains (pencils) which are each assigned to an MPI-process. The derivatives and interpolations in the x-direction (y-direction, z-direction) are performed in X-pencils (Y-pencils, Z-pencils), respectively. The 3-D FFTs required by the Poisson solver are also broken down as a series of 1-D FFTs computed in one direction at a time. Global transpositions to switch from one pencil to another are performed with the MPI command MPI_ALLTOALL(V). The flow solvers within Xcompact3d can scale well with up to hundreds of thousands of MPI-processes for simulations with several billion mesh nodes (Laizet & Ning Reference Laizet and Ning2011; Bartholomew et al. Reference Bartholomew, Deskos, Frantz, Schuch, Lamballais and Laizet2020). All the simulations in this study were carried out on ARCHER2, the UK supercomputing facility. It is equipped with nodes based on dual AMD EPYC $^\textit{TM}$ 7742 processors running at 2.25 GHz, totalling 128 cores per node.

In this numerical study, the following incompressible Navier–Stokes equations are solved

(4.1) \begin{equation} \begin{aligned} \frac {\partial \boldsymbol{u}}{\partial t}+ \frac {1}{2} [\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u} \otimes \boldsymbol{u})+(\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}] &= -\boldsymbol{\nabla }p + (I-\boldsymbol{\varGamma })\mathcal{D}_1+\boldsymbol{\varGamma }\mathcal{D}_2 +\boldsymbol{f}(t), \\ \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &= 0. \end{aligned} \end{equation}

With a slight abuse of notation, $\mathcal{D}_1$ indicates that the diffusion term $\nu \Delta \boldsymbol{u}$ is directly implemented using a conventional sixth-order finite-difference scheme for its second-order spatial derivatives; while $\mathcal{D}_2$ refers to a numerical implementation of $\nu \Delta \boldsymbol{u}$ which adds targeted numerical dissipation at small scales via a customised sixth-order finite-difference scheme. Full details can be found in Lamballais et al. (Reference Lamballais, Fortuné and Laizet2011), Dairay et al. (Reference Dairay, Lamballais, Laizet and Vassilicos2017) and Mahfoze & Laizet (Reference Mahfoze and Laizet2021). The operator $\boldsymbol{\varGamma }$ is used to balance the weight of the two diffusive terms and can depend on the both local geometry and flow velocities via $(\boldsymbol{\varGamma }\boldsymbol{u})(x) = \varGamma (x,\boldsymbol{u})\boldsymbol{u}(x)$ , where $\varGamma (x,\boldsymbol{u}) \in {\mathbb{R}}^{3 \times 3}$ . The choice $\varGamma =0$ corresponds to a DNS implementation, in which (4.1) is mathematically equivalent to (1.4). The case $\varGamma = I$ corresponds to ILES, for which, the unknowns $\mathbf{u}(x,t)$ and $p(x,t)$ should be interpreted as the large-scale component of velocity and pressure. Note finally that the advection terms in (4.1) are written in skew-symmetric form in order to reduce aliasing errors (Kravchenko & Moin Reference Kravchenko and Moin1997).

4.1. Numerical implementation of the forcing term

To allow for a validation with the numerical study of Biau (Reference Biau2023), we will perform simulations in the case of momentum thickness control (see § 2.2) where the forcing amplitude $f(t)$ is given by (2.12). The distinction between the two considered numerical methods (i.e. DNS and ILES) introduces an extra subtlety in terms of how $f(t)$ must be chosen to achieve a desired unity value of the momentum thickness. For DNS, we use the explicit formula (2.12), letting

(4.2) \begin{equation} f(t) = f_{\textit{DNS}}(t) := \frac {-ke_\theta (t) - \frac 12 C_{\!f}(t) + \frac {2}{\textit{Re}} \int _0^\infty \left ( \frac { \partial U}{\partial y}\right )^2 {\textrm{d}}y - 2 \int _0^\infty \langle uv \rangle \frac {\partial U}{\partial y} \, {\textrm{d}}y}{\theta (t)}. \end{equation}

Lemma 2 then guarantees that $\lim _{t \rightarrow \infty } \theta (t) =1$ for any solution to PDE (1.4), meaning that one would expect a DNS of (4.1) to have the same behaviour. We show in § 5 that this is indeed the case.

In the case of ILES it is not appropriate to define the forcing amplitude by (2.12) since this expression is only accurate for solutions of the Navier–Stokes equations (1.4), while ILES only gives an approximate solution via (4.1). This subtlety can be avoided by adding an ‘integral control’ term to the forcing and using

(4.3) \begin{equation} f_{\textit{ILES}}(t) := f_{\textit{DNS}}(t) - k^2 \frac {z(t)}{\theta (t)}, \end{equation}

where $z(t) \in {\mathbb{R}}$ is defined as the integral of the momentum thickness via

(4.4) \begin{equation} \frac {{\textrm{d}}z}{{\textrm{d}}t} = e_\theta (t) = 1- \theta (t). \end{equation}

This accounts, with the addition of only one extra scalar-valued variable, for the small difference between the theoretical value (2.12) and the forcing required to maintain a unity value of momentum thickness in ILES. For simplicity, a gain value of $k=1$ is used throughout this paper.

4.2. Computational domains, data collection and computational cost

The DNS are performed at $\textit{Re}_{\theta }=1000, 2000$ and the ILES at $\textit{Re}_{\theta }=4060, 6500, 8300$ . These Reynolds numbers have been selected in order to allow for a direct comparison with published numerical data of a spatially evolving turbulent boundary layer. The spatial discretisation, temporal discretisation and computational cost of each simulation run are given in table 1. With increasing Reynolds number it is necessary to use both a longer temporal window and a larger domain in the streamwise direction to obtain converged statistics. For example, $L_{x}$ is chosen to be $40/\theta$ for $\textit{Re}_\theta =1000$ and is increased by a factor of three for $\textit{Re}_\theta = 8300$ .

Table 1. Simulation details for the DNS and ILES of (4.1).

The flow statistics are computed by first performing a space/time average to obtain the mean streamwise velocity profile $U(y)$ and, subsequently, the wall shear stress $\tau _{w}=\nu U'(0)$ and friction velocity $U_{\tau }:=\sqrt {\tau _{w} / \rho } = \sqrt {\tau _w}$ . Using $l_\ast =\nu / U_{\tau }$ as a length scale and $t_\ast = l_\ast / U_{\tau }$ as a time scale, all variables can be expressed in wall units, e.g. $u^{+}=U/U_{\tau }$ , $y^{+}=y/l_\ast$ , and $t^{+}=t/t_\ast$ . In addition to the displacement thickness $\delta ^\ast$ and momentum thickness $\theta$ , we will also consider the ‘shape factor’ $H_{12}:=\delta ^\ast /\theta$ , and the wake parameter

(4.5) \begin{equation} C_{w\textit{ake}} = \int _0^\infty \big(U^{+}_{\infty }- U^{+} \big) d(y/\delta ), \end{equation}

introduced by Coles (Reference Coles1956), where $\delta$ is the boundary layer thickness defined in terms of $99\,\%$ of the free-stream velocity.

The mean values of all flow statistics reported in subsequent sections are collected after the flow was observed to have transitioned to a statistically steady state. For the five Reynolds numbers considered in table 1, this was deemed to occur at non-dimensional times $t_0^+ = 1976, 3418, 7534, 12\,208$ and $16\,511$ , ordered in terms of increasing $\textit{Re}_\theta$ . Each flow statistic was then averaged over a window $t^+ \in [t_0^+,T_{f}^+]$ , respectively (see table 1), in order of increasing $\textit{Re}_\theta$ . To describe statistical convergence, suppose that $g(t)$ is a given time-dependent flow property of interest, and $\bar {g}$ is its average over the window $[t_0^+,T_{f}^+]$ . Letting

(4.6) \begin{equation} G(t^+) = \frac {1}{t^+-t_0^+} \int _{t_0^+}^{t^+} g(s) {\textrm{d}}s, \end{equation}

be the moving average of $g$ , we define the maximum percentage deviation of the moving average from the reported mean over the final half of the averaging window by

(4.7) \begin{equation} E = \max \left \{ 100\,\% \boldsymbol{\cdot }\left |\frac {G(t^+) - \bar {g}}{\bar {g}} \right | : \frac { t_0^+ + T_f^+}{2} \leqslant t^+ \leqslant T_f^+ \right \}\!. \end{equation}

The maximum value of $E$ for any of the flow statistics $\delta ^\ast , U_\tau , f, C_{w\textit{ake}}$ that we report subsequently in table 3 is, $1.2\,\%,0.9\,\%,0.9\,\%,1.1\,\%$ and $1.5\,\%$ , respectively, for the five cases considered in order of increasing $\textit{Re}_\theta$ .

Finally, we note that there are significant differences in computational cost, detailed in 1, of the proposed numerical scheme in comparison with that of the highlighted reference studies listed in table 2. The DNS of Sillero et al. (Reference Sillero, Jiménez and Moser2013), achieving a maximum Reynolds number of $\textit{Re}=6500$ , is computationally expensive, and required 45 million core hours while using a large number of cores (32 768). The LES of Eitel-Amor et al. (Reference Eitel-Amor, Orlu and Schlatter2014), achieving a maximum Reynolds number of $\textit{Re}_\theta = 8300$ , used only around 1 million core hours and 4096 cores. In this study, we sought a balance between computational cost and efficiency. At the highest Reynolds number considered in this study, $\textit{Re}_\theta =8300$ , the proposed numerical scheme used 1 60 384 core hours on 8192 cores, a sixfold reduction compared with the LES of Eitel-Amor et al. (Reference Eitel-Amor, Orlu and Schlatter2014).

Table 2. The DNS and LES grid resolutions, and figure formatting conventions.