2.1. Problem set-up and equations of motion

The incompressible flow of a fluid with kinematic viscosity $\nu$ in a pressure-driven channel with fixed streamwise pressure gradient $\varPi$ is examined. The streamwise, wall-normal and spanwise directions, normalised by the channel mean half-height $h$, are identified by the Cartesian coordinates $(x_1, x_2, x_3)$, with the origin of the wall-normal coordinate located at the channel midplane. The friction velocity $u_\tau =\sqrt {\tau _w/\rho }$, with $\tau _w = h \varPi$ the mean wall friction, is used to normalise the velocity components $(u_1, u_2, u_3)$ along the three directions. The reference pressure is $p_{ref}=\rho u_\tau ^2$ and this leads to a non-dimensional pressure gradient $\partial \bar {p}/\partial x_i = \delta _{i1}$, with $\delta _{ij}$ being the Kronecker delta. Reynolds-averaging produces the mean velocity $\bar {u}_i$ and the fluctuation $u^\prime _i$. The superscript $(\cdot )^+$, generally used for inner scaled quantities, is omitted in the following to reduce clutter, unless necessary. With these definitions, the friction Reynolds number is $Re_\tau =u_\tau h/\nu$. We consider channels with streamwise-independent modulations of the wall topography, namely, sinusoidal modulations and rectangular ridges, as illustrated in figure 1.

Figure 1. (a) Sinusoidal and (b) ridge-type topographies considered in this paper. The coordinate system $(x_1,x_2,x_3)$, with origin on the symmetry plane, is shown. The streamwise direction $x_1$ is oriented into the page. When scaled by $h$, the mean channel height is equal to $2$. Symmetric configurations obtained by mirroring the lower wall geometries shown in the diagrams about the midplane $x_2 = 0$ are considered. For sinusoidal topographies, the period of the modulation is denoted by $\lambda _3$. For ridge-type topographies, the spacing between elements (the period) is denoted by $S$, while $W$ and $G$ are used to indicate the ridge width and the gap between elements, respectively.

The time-averaged flow structure in the channel is governed by the non-dimensional Reynolds-averaged continuity and momentum equations

(2.1a)$$\begin{gather} \frac{\partial \bar{u}_i}{\partial x_i} = 0, \end{gather}$$

(2.1b)$$\begin{gather}\bar{u}_j\frac{\partial \bar{u}_i}{\partial x_j}={-} { \delta_{i1}}+\frac{1}{Re_\tau}\frac{\partial^2 \bar{u}_i}{\partial x_j^2}-\frac{\partial \overline{u_i'u_j'}}{\partial x_j}, \end{gather}$$

with no-slip boundary conditions on the two walls. As common, the trace of the Reynolds stress tensor is absorbed in the pressure term and we thus introduce the traceless stress tensor $\tau _{ij} = -\overline {u_i'u_j'} + \frac {1}{3}\overline {u_i'u_j'}\delta _{ij}$. Assuming that a streamwise-independent mean flow (i.e. $\partial (\cdot )/\partial x_1 \equiv 0$) develops over streamwise-independent modulations, the mean pressure can be eliminated by employing a streamwise velocity/stream function formulation, where the stream function $\bar {\psi }$ satisfies $\nabla ^2 \bar {\psi }=\bar {\omega }_1$ with

(2.2)\begin{equation} \bar{\omega}_1= {\frac{\partial \bar{u}_3}{\partial x_2}- \frac{\partial \bar{u}_2}{\partial x_3}},\end{equation}

the streamwise vorticity. With these definitions, the cross-stream velocity components are $\bar {u}_2=-{\partial \bar {\psi }}/{\partial x_3}$ and $\bar {u}_3={\partial \bar {\psi }}/{\partial x_2}$, satisfying automatically the continuity equation reduced to the cross-plane section. The Reynolds-averaged streamwise momentum and stream function equations then become

(2.3a)\begin{align} &\frac{\partial \bar{\psi}}{ \partial x_2} \frac{\partial \bar{u}_1}{ \partial x_3} - \frac{\partial \bar{\psi}}{ \partial x_3} \frac{\partial \bar{u}_1}{ \partial x_2} = 1 + \frac{1}{Re_\tau} \left( \frac{\partial^2 \bar{u}_1}{\partial x_2^2} + \frac{\partial^2 \bar{u}_1}{\partial x_3^2} \right) + \frac{\partial \tau_{12}}{\partial x_2} + \frac{\partial \tau_{13}}{\partial x_3}, \end{align}

(2.3b)\begin{align} &\frac{\partial^2}{\partial x_2 \partial x_3} \left[ \left(\frac{\partial \bar{\psi}}{\partial x_2}\right)^2 -\left( \frac{\partial \bar{\psi}}{\partial x_3}\right)^2 \right] + \left( \frac{\partial^2}{\partial x_3^2} - \frac{\partial^2}{\partial x_2^2}\right) \frac{\partial \bar{\psi}}{\partial x_2} \frac{\partial \bar{\psi}}{\partial x_3} \nonumber\\ & \quad =\frac{1}{Re_\tau} \left( \frac{\partial^2}{\partial x_2^2} +\frac{\partial^2}{\partial x_3^2} \right)^2 \bar{\psi} +\frac{\partial^2}{\partial x_2 \partial x_3} ( \tau_{33} - \tau_{22}) + \left( \frac{\partial^2}{\partial x_2^2} - \frac{\partial^2}{\partial x_3^2} \right) \tau_{23}. \end{align}

2.2. Linearised response model

Without loss of generality, we assume the wall modulation to be spanwise periodic, with fundamental period $\lambda _3$. We only consider zero-mean modulations of the wall geometry since perturbations of the mean channel height are trivially explained as a change in the Reynolds number, or as a wall-normal shift of the flow characteristics in boundary layers. Hence, an arbitrary modulation can be expressed by a function $f(x_3)$, with cosine series

(2.4)\begin{equation} f(x_3) = \sum_{n=1}^{\infty} f^n \cos(n k_3 x_3), \end{equation}

with $k_3=2 \pi /\lambda _3$ the fundamental wavenumber and $f^n$ the amplitude of the $n$th wavenumber mode. Expressions for $f(x_3)$ for the two surfaces considered in the present work are given in (3.1) and (4.1), respectively. Following Russo & Luchini (Reference Russo and Luchini2016), we then assume that the amplitude of the modulation is smaller than any other relevant geometric or flow length scale and we introduce a small parameter $\epsilon \ll 1$. The lower channel wall is then located at $x_2 = - 1 + \epsilon f(x_3)$, while several configurations are possible for the upper wall. In one of the alternatives, the upper wall is located at $x_2 = 1 - \epsilon f(x_3)$, defining symmetric channels where secondary currents occupy at most half the channel height. In antisymmetric channels, where the upper wall is located at $x_2 = 1 + \epsilon f(x_3)$, as in Vidal et al. (Reference Vidal, Nagib, Schlatter and Vinuesa2018), the secondary currents can occupy the entire channel height and interact with the two shear layers developing over the top and bottom walls. In order to model more closely secondary currents in boundary layers (Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015; Hwang & Lee Reference Hwang and Lee2018; Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2020) or open channel flows (Zampiron et al. Reference Zampiron, Cameron and Nikora2020), where the secondary currents develop across the shear flow, symmetric channels are considered in this paper.

In a small-modulation scenario, a generic time-averaged quantity ${q}(x_2, x_3)$ in the channel with modulated walls (dropping the overbar to reduce clutter) can be expanded in a Taylor series in $\epsilon$ as

(2.5)\begin{equation} q(x_2, x_3)=q^{(0)}(x_2)+\epsilon {q}^{(1)}(x_2, x_3)+ O(\epsilon^2),\end{equation}

where ${q}^{(0)}$ denotes the plane channel solution. This expansion implies that the strength of secondary flows produced by a shallow modulation varies linearly with the amplitude $\epsilon$ and the perturbation quantity ${q}^{(1)}$ can be thus interpreted as the flow response (i.e. secondary currents) for a unitary change of the wall geometry given by (2.4).

Substituting the Taylor expansion (2.5) for all flow variables in the Reynolds-averaged equations (2.3b) and considering terms at order zero in $\epsilon$, the time-averaged streamwise momentum equation is

(2.6)\begin{equation} 0=1 + \frac{1}{Re_\tau} \frac{\partial^2 u_1^{(0)}}{\partial x_2^2} + \frac{\partial \tau_{12}^{(0)}}{\partial x_2},\end{equation}

while the stream function equation is trivially satisfied, since $u^{(0)}_2 = u^{(0)}_3 = 0$ in a plane channel. Retaining terms at order one in $\epsilon$, we obtain the set of equations

(2.7a)$$\begin{gather} -\frac{\partial \psi^{(1)}}{\partial x_3} \varGamma = \frac{1}{Re_\tau} \left(\frac{\partial^2 }{\partial x_2^2} + \frac{\partial^2 }{\partial x_3^2}\right) u_1^{(1)} + \frac{\partial \tau_{12}^{(1)} }{\partial x_2} + \frac{\partial \tau_{13}^{(1)} }{\partial x_3}, \end{gather}$$

(2.7b)$$\begin{gather}0 = \frac{1}{Re_{\tau}} \left( \frac{\partial^2 }{\partial x_2^2} + \frac{\partial^2}{\partial x_3^2}\right)^2 \psi^{(1)} + \frac{\partial^2}{\partial x_2 \partial x_3} ( \tau_{33}^{(1)} - \tau_{22}^{(1)}) + \left( \frac{\partial^2}{\partial x_2^2} - \frac{\partial^2}{\partial x_3^2}\right) \tau_{23}^{(1)}, \end{gather}$$

where $\varGamma =\partial u_1^{(0)}/\partial x_2$. These equations describe the new equilibrium between the perturbation of mean flow quantities ($u_1^{(1)}$, $\psi ^{(1)}$) and the perturbation of the turbulent stress tensor $\tau _{ij}^{(1)}$. It is worth pointing out that the term ${\partial \psi ^{(1)}}/ {\partial x_3} \varGamma$, analogous to the off-diagonal coupling operator in the Orr–Sommerfeld–Squire linearised equations (Schmid & Henningson Reference Schmid and Henningson2000), is the only coupling term explicitly appearing in this set of equations. Physically, this terms produces a spanwise modulation of the streamwise velocity as a result of secondary motions in the cross-stream plane.

The key property of these equations is linearity, since second-order perturbation–perturbation terms arising from the convective nonlinearity are neglected at order one. As pointed out in Meyers et al. (Reference Meyers, Ganapathisubramani and Cal2019), neglecting these terms is justified by the fact that the cross-stream velocity components are generally quite weak, i.e. less than 5 % the external velocity scale (Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015; Hwang & Lee Reference Hwang and Lee2018; Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2020), especially at large distances from the wall. The key advantage is that the flow response induced by an arbitrary, spectrally complex modulation $f(x_3)$ can be obtained by appropriately combining solutions of linear equations obtained at each spanwise wavenumber characterising the modulation in the expansion (2.4).

2.3. Nonlinear Reynolds stress model

To close the mean equations at order zero and one, it is now necessary to express the Reynolds stress tensor as a function of other mean quantities. One option is to introduce a linear Boussinesq hypothesis, using the turbulent eddy viscosity $\nu _t$ to derive the linear constitutive relation

(2.8)\begin{equation} \tau_{ij}^L= 2 \nu_t S_{ij}, \end{equation}

with $S_{ij}$ the mean velocity gradient tensor

(2.9)\begin{equation} S_{ij}=\frac{1}{2} \left( \frac{\partial \bar{u}_i}{ \partial x_j}+ \frac{\partial \bar{u}_j}{\partial x_i}\right). \end{equation}

Expanding the turbulent stresses in a Taylor series as in (2.5), the leading terms at order zero and one are

(2.10)$$\begin{gather} \tau_{ij}^{L(0)}=2 \nu_t^{(0)} S_{ij}^{(0)}, \end{gather}$$

(2.11)$$\begin{gather}\tau_{ij}^{L(1)}=2 \nu_t^{(0)} S_{ij}^{(1)}+2 \nu_t^{(1)} S_{ij}^{(0)}, \end{gather}$$

where $\nu _t^{(1)}$ is the unknown perturbation of the eddy-viscosity profile induced by the wall modulation. When a linear relation is used, however, no secondary flows are predicted (Perkins Reference Perkins1970; Speziale Reference Speziale1982; Bottaro et al. Reference Bottaro, Soueid and Galletti2006). In fact, the stresses appearing in (2.7b) would not depend on the streamwise velocity since the stress tensor is isotropic and the stream function equation (2.7b) decouples from the streamwise momentum equation (2.7a). Transient energy amplification from inhomogeneous initial conditions can be observed (del Álamo & Jiménez Reference del Álamo and Jiménez2006; Pujals et al. Reference Pujals, Garcìa-Villalba, Cossu and Depardon2009) but the steady response to an exogenous forcing, for example, from the wall modulation, is trivial, $\psi ^{(1)} \equiv 0$. Hence, a nonlinear Reynolds stress model is necessary. Several approaches have been described in the literature (e.g. Speziale Reference Speziale1991; Speziale, Sarkar & Gatski Reference Speziale, Sarkar and Gatski1991; Chen, Lien & Leschziner Reference Chen, Lien and Leschziner1997). Here we use the quadratic constitutive relation (QCR) nonlinear model introduced by Spalart (Reference Spalart2000), which contains simple terms proportional to the product of the rotation and the strain tensors. This model was recently utilised by Spalart, Garbaruk & Stabnikov (Reference Spalart, Garbaruk and Stabnikov2018) to predict the high-Reynolds number asymptotic properties of secondary flows in square and elliptical ducts, providing a good approximation of the secondary vortical flow topology and of the wall friction coefficient. Compared with other approaches, the QCR model is straightforward to manipulate analytically, and it is thus chosen here to remain in the original spirit of developing a simple predictive model of secondary flows over heterogeneous surfaces.

In the QCR model, the Reynolds stresses become

(2.12)\begin{equation} \tau_{ij}=\tau_{ij}^{L}-C_{r1}[ O_{ik}\tau_{jk}^{L} +O_{jk}\tau_{ik}^{L}], \end{equation}

where the tuning constant $C_{r1}$ controls the anisotropy of the Reynolds stress tensor. Spalart (Reference Spalart2000) suggests using $C_{r1}=0.3$ to match the anisotropy in the outer region of wall-bounded turbulent flows and we follow this indication in this paper. In (2.12), $O_{ij}$ is the normalised rotation tensor

(2.13)\begin{equation} O_{ij} = \frac{2W_{ij}}{\sqrt{ {\dfrac{\partial \bar{u}_m}{\partial x_n} \dfrac{\partial \bar{u}_m}{\partial x_n}}}}, \quad \text{with } W_{ij} = \frac{1}{2}\left( \frac{\partial \bar{u}_i}{\partial x_j}-\frac{\partial \bar{u}_j}{\partial x_i}\right). \end{equation}

At order zero, the nonlinear stress tensor is equal to the expression obtained from the linear constitutive relation. At first order, the Reynolds stress tensor is

(2.14)\begin{equation} \tau_{ij}^{(1)} = \tau_{ij}^{L(1)} - C_{r1} [O_{ik}^{(1)}\tau_{jk}^{L(0)} + O_{ik}^{(0)} \tau_{jk}^{L(1)} + O_{jk}^{(1)}\tau_{ik}^{L(0)} + O_{jk}^{(0)}\tau_{ik}^{L(1)}], \end{equation}

where $O_{ij}^{(1)}$ is the normalised rotation tensor induced by the first-order velocity components (see Appendix A). Developing (2.14), the individual perturbation Reynolds stresses appearing in (2.7) are

(2.15a)$$\begin{gather} \tau_{12}^{(1)}=\nu_t^{(0)} \frac{\partial u_1^{(1)}}{\partial x_2}+ \nu_t^{(1)}\varGamma+2C_{r1} \mathrm{sign}(\varGamma)\nu_t^{(0)} \frac{\partial^2 \psi^{(1)}}{\partial x_2 \partial x_3}, \end{gather}$$

(2.15b)$$\begin{gather}\tau_{13}^{(1)}=\nu_t^{(0)} \frac{\partial u_1^{(1)}}{\partial x_3} -2C_{r1} \mathrm{sign}(\varGamma)\nu_t^{(0)} \frac{\partial^2 \psi^{(1)}}{\partial x_2^2}, \end{gather}$$

(2.15c)$$\begin{gather}\tau_{23}^{(1)}=\nu_t^{(0)}\left( \frac{\partial^2}{\partial x_2^2}- \frac{\partial^2}{\partial x_3^2}\right) \psi^{(1)}+2 C_{r1} \mathrm{sign}(\varGamma) \nu_t^{(0)} \frac{\partial u_{1}^{(1)}}{\partial x_3}, \end{gather}$$

(2.15d)$$\begin{gather}\tau_{22}^{(1)}={-}2 \nu_t^{(0)} \frac{\partial^2 \psi^{(1)}}{\partial x_2 \partial x_3}+2 C_{r1} \left[ \mathrm{sign}(\varGamma)\nu_t^{(0)} \frac{\partial u_{1}^{(1)}}{\partial x_2}+\mathrm{sign}(\varGamma) \nu_t^{(1)}\varGamma\right], \end{gather}$$

(2.15e)$$\begin{gather}\tau_{33}^{(1)}=2 \nu_t^{(0)} \frac{\partial^2 \psi^{(1)}}{\partial x_2 \partial x_3}, \end{gather}$$

where ‘$\mathrm {sign}$’ is the sign function. Except for $\tau _{33}^{(1)}$, which coincides with its linear Boussinesq definition, all other stresses contain an additional term specific to the QCR model, which results in a tighter, two-way coupling between the stream function and streamwise velocity equations, able to sustain secondary currents.

2.4. Eddy-viscosity transport model

The perturbation of the turbulent stresses (2.15) still contains the unknown perturbation eddy viscosity $\nu _t^{(1)}$. Past studies that have utilised linearised RANS equations to examine transient energy amplification in plane turbulent channels (del Álamo & Jiménez Reference del Álamo and Jiménez2006; Pujals et al. Reference Pujals, Garcìa-Villalba, Cossu and Depardon2009) have often used analytical eddy-viscosity profiles (Cess Reference Cess1958; Reynolds & Hussain Reference Reynolds and Hussain1972). In these works, the eddy viscosity was assumed to be constant and not influenced by the growth of the optimal structures. This assumption, however, has little physical justification for a modulated geometry. To provide a better description of the eddy-viscosity distribution in the modulated geometry and capture transport effects, we use in the present paper the one-equation SA turbulence transport model (Spalart & Allmaras Reference Spalart and Allmaras1994), initially developed for attached shear flows. Using the channel half-height and the friction velocity for normalisation, the SA model introduces one transport equation for the transformed eddy viscosity $\tilde {\nu }$ related to the turbulent viscosity by the relation

(2.16)\begin{equation} \nu_t=\tilde{\nu} f_{v1}, \end{equation}

where

(2.17)\begin{equation} f_{v1}=\frac{\chi^3}{\chi^3+c_{v1}^3}, \end{equation}

with $\chi = Re_\tau \tilde {\nu }$ and $c_{v1}$ a tuning constant. The modified eddy viscosity coincides with the turbulent viscosity away from the wall. Additionally, the term (2.17) ensures the correct decay of the turbulent viscosity in the viscous sublayer (Spalart & Allmaras Reference Spalart and Allmaras1994; Herring & Mellor Reference Herring and Mellor1968) when $\tilde {\nu }$ behaves linearly in the log-layer down to the surface, which is advantageous for numerical reasons. The steady transport equation for $\tilde {\nu }$,

(2.18)\begin{equation} \bar{u}_i \frac{\partial \tilde{\nu}}{\partial x_i} =c_{b1} \tilde{\mathcal{S}} \tilde{\nu}+\frac{1}{\sigma}\left\{ \frac{\partial }{\partial x_j} \left[ \left(\frac{1}{Re_\tau}+\tilde{\nu}\right) \frac{\partial \tilde{\nu} }{\partial x_j }\right] + c_{b2} \frac{\partial \tilde{\nu}}{ \partial x_j}\frac{\partial \tilde{\nu}}{ \partial x_j}\right\}-c_{w1} f_w \left( \frac{\tilde{\nu}}{d}\right)^2, \end{equation}

is composed of convection, production, diffusion and destruction terms. In the production term, the quantity $\tilde {\mathcal {S}}$ is defined as

(2.19)\begin{equation} \tilde{\mathcal{S}}=\sqrt{2W_{ij} W_{ij}}+ \frac{\tilde{\nu}}{k^2 d^2} f_{v2} \quad \text{with } f_{v2}= 1- \frac{\chi}{1+ \chi f_{v1}}, \end{equation}

with $k$ the von Kármán constant. The destruction term in (2.18) captures the blocking effect of the wall on turbulent fluctuations and is a function of the distance to the nearest surface $d$. With this term, the model produces an accurate log-layer in wall-bounded flows. It includes a non-dimensional function $f_{w}$ that increases the decay of the destruction term in the outer region. This term reads as

(2.20)\begin{equation} f_{w} = g \left[ \frac{1+ c_{w3}^6}{ g^6 + c_{w3}^6}\right]^{1/6}, \end{equation}

with

(2.21a,b)\begin{equation} g = r+c_{w2} \left( r^6 - r\right) \quad \mathrm{and} \quad r = \frac{\tilde{\nu}}{ \tilde{\mathcal{S}} k^2 d^2}. \end{equation}

Standard values for the calibration constants $c_{v1}=7.1$, $c_{b1}=0.1355$, $\sigma =2/3$, $c_{b2}=0.622$, $c_{w2}=0.3$, $c_{w3}=2$ are used (Spalart & Allmaras Reference Spalart and Allmaras1994), with $c_{w1}=c_{b1}/k^2+(1+c_{b2})/\sigma$ to balance production, diffusion and destruction in the log-layer and with $k=0.41$.

Expanding all flow variables in a Taylor series, the transport equation for the modified eddy viscosity at order zero and one can be obtained. At order zero, the equation is trivially obtained from (2.8) and it is omitted here. At first order, the eddy viscosity $\nu _t^{(1)}$ appearing in the stresses (2.15) can be readily obtained as

(2.22)\begin{equation} \nu_t^{(1)}=\tilde{\nu}^{(1)} f_{v1}^{(0)}+\tilde{\nu}^{(0)} f_{v1}^{(1)},\end{equation}

where $f_{v1}^{(1)}$ and other additional terms appearing at first order are reported in Appendix B. In the linearisation process, it is key to observe that the topographic modulation can be thought of as a perturbation of the distance from the solid wall. This is a key physical parameter in the SA turbulence model as it controls the formation of a log-layer through the balance of production and destruction, where it appears directly. In particular, the distance is expanded as

(2.23)\begin{equation} d(x_2, x_3) = d^{(0)}(x_2) + \epsilon d^{(1)}(x_2, x_3), \end{equation}

with $d^{(0)}$ the original distance in the plane channel and

(2.24)\begin{equation} d^{(1)}(x_2, x_3) = \mathrm{sign}(x_2) f(x_3), \end{equation}

where the sign function in (2.24) captures the symmetric modulation of the walls and models the fact that the distance from the nearest physical wall decreases/increases for points above the crests/troughs of the topography in the lower channel half, as illustrated in figure 2.

Figure 2. Illustration of the effect of topographic modulations on the distance $d$ appearing in the production and destruction terms of the SA transport model. For a point $(x_2, x_3)$ above the trough in the lower channel half, the (positive) distance to the nearest wall increases from $d^{(0)}$, the original distance from the flat lower wall, by an amount $d^{(1)} = -f(x_3)$. Opposite effects are produced on the crests of the topography or in the upper half of the channel.

After algebraic operations, the transport equation for the perturbation of the modified eddy viscosity $\tilde {\nu }^{(1)}$ reads as

(2.25)\begin{align} -\frac{\partial \psi^{(1)}}{\partial x_3} \frac{\partial \tilde{\nu}^{(0)}}{\partial x_2} &= \frac{1}{\sigma}\left( \frac{1}{R_{e_\tau}}+\tilde{\nu}^{(0)}\right) \left( \frac{\partial^2}{\partial x_2^2}+ \frac{\partial^2}{\partial x_3^2}\right) \tilde{\nu}^{(1)}+\frac{1}{\sigma}\frac{\partial^2 \tilde{\nu}^{(0)}}{\partial x_2^2}\tilde{\nu}^{(1)} \nonumber\\ &\quad +\frac{1}{\sigma}(2+2 c_{b2}) \frac{\partial \tilde{\nu}^{(0)}}{\partial x_2} \frac{\partial \tilde{\nu}^{(1)}}{\partial x_2} +c_{b1}\tilde{\nu}^{(0)} \tilde{\mathcal{S}}^{(1)}+c_{b1} \tilde{\nu}^{(1)} \tilde{\mathcal{S}}^{(0)} \nonumber\\ &\quad -2 \tilde{\nu}^{(0)} c_{w1} f_w^{(0)} \frac{\tilde{\nu}^{(1)} d^{(0)}- \tilde{\nu}^{(0)} d^{(1)}}{d^{(0)3}} - c_{w1} f_w^{(1)} \left( \frac{\tilde{\nu}^{(0)}}{d^{(0)}}\right)^2. \end{align}

This equation is coupled to the stream function equation by the convective transport term on the left-hand side, modelling the wall-normal transport of the background turbulent fluctuations by the secondary motions. An additional coupling term with the streamwise momentum equation appears in the production term $\tilde {\mathcal {S}}^{(1)}$, which models the change in the production of turbulent kinetic energy as a result of the distortion of the streamwise velocity profile.

2.5. Linearised boundary conditions

Boundary conditions for the linearised transport equations are now derived using established methods (Busse & Sandham Reference Busse and Sandham2012; Luchini Reference Luchini2013). Assuming that the topographic perturbation is small, we retain the original rectangular geometry of the domain but we introduce inhomogeneous boundary conditions on the perturbation quantities derived by imposing the original conditions on the displaced surface.

Considering the lower wall, expanding the velocity near the surface in a Taylor series and enforcing the no-slip condition we obtain

(2.26)\begin{equation} \bar{u}_i({-}1+\epsilon f(x_3), x_3) = \bar{u}_i|_{x_2={-}1}+ \epsilon f(x_3) \left.\frac{\partial \bar{u}_i}{ \partial x_2} \right|_{x_2={-}1} = 0.\end{equation}

Substituting the expansion (2.5) for the velocity in (2.26), noting that $u_i^{(0)}=0$ at $x_2=-1$, and retaining terms at order one in $\epsilon$ provides

(2.27)\begin{equation} u_{i}^{(1)} |_{x_2={-}1} + f(x_3) \left.\frac{\partial u_i^{(0)}}{\partial x_2}\right|_{x_2={-} 1} = 0, \end{equation}

i.e. the perturbation velocity at the boundary of the numerical domain is proportional to the wall-normal gradient of the velocity in the plane channel to preserve the no-slip condition on the modulated topography. The boundary condition on the streamwise velocity perturbation then becomes

(2.28)\begin{equation} u_{1}^{(1)}(x_2 ={-}1) ={-} f(x_3) \left.\frac{\partial u^{(0)}}{\partial x_2}\right|_{x_2={-} 1} ={-} f(x_3) Re_{\tau}, \end{equation}

while $u_{3}^{(1)}(x_2=-1) = 0$ and $u_{2}^{(1)}(x_2=-1) = 0$. The boundary conditions for the perturbation stream function are

(2.29)\begin{equation} \frac{\partial \psi^{(1)}}{\partial x_2}(x_2 ={-}1) = \psi^{(1)}(x_2 ={-}1) = 0.\end{equation}

Using a similar strategy, and noting that the modified eddy viscosity satisfies homogeneous boundary conditions at the wall (Spalart & Allmaras Reference Spalart and Allmaras1994), the inhomogeneous boundary condition

(2.30)\begin{equation} \tilde{\nu}^{(1)} (x_2={-} 1)={-} f(x_3) \left. \frac{\partial \tilde{\nu}^{(0)}}{\partial x_2}\right|_{x_2={-}1}={-}f(x_3)k,\end{equation}

can be derived for the perturbation of the transformed eddy-viscosity at the lower numerical boundary. The last equality holds since the modified eddy viscosity obeys the linear relation $\tilde {\nu } = k x_2$ near the wall (Spalart & Allmaras Reference Spalart and Allmaras1994). No conditions are required for the eddy viscosity $\nu _t$, since this is not directly associated with a transport equation in the SA model. With a similar procedure, boundary conditions on the upper numerical boundary can be obtained. Equivalently, symmetric boundary conditions can also be applied at the channel centreline when symmetric channels are studied, to reduce computational costs. However, in the present work, the full channel domain with linearised boundary conditions on both upper and bottom surfaces was considered, as it was easily modelled using available Chebyshev discretisation tools.

2.6. Fourier spectral expansion of the solution

When using linearised equations, any arbitrary topography can be analysed by examining each fundamental spanwise length scale separately from the others. The solution of the linearised equations can be first expressed by the Fourier series

(2.31a)$$\begin{gather} u_1^{(1)} (x_2, x_3) =\sum_{n=1}^{\infty} \hat{u}_1(x_2; n) \cos{(n k_3 x_3)}, \end{gather}$$

(2.31b)$$\begin{gather}\psi^{(1)} (x_2, x_3) =\sum_{n=1}^{\infty} \hat{\psi}(x_2; n) \sin{(n k_3 x_3)}, \end{gather}$$

(2.31c)$$\begin{gather}\tilde{\nu}^{(1)} (x_2, x_3) =\sum_{n=1}^{\infty} \hat{\nu}(x_2; n) \cos{(n k_3 x_3)}, \end{gather}$$

where $\hat {u}_1(x_2;n)$, $\hat {\psi }(x_2;n)$ and $\hat {\nu }(x_2;n)$ are the real-valued, wall-normal profiles of the perturbation streamwise velocity, stream function and modified eddy viscosity at each integer spanwise wavenumber $n$. Then, components at different spanwise wavenumbers decouple, forming the set of three ordinary differential equations

(2.32a)$$\begin{gather} -n k_3 \hat{\psi} \varGamma = \frac{1}{Re_\tau} \left(\frac{\mathrm{d}^2 }{\mathrm{d} x_2^2} -n^2 k_3^2 \right) \hat{u}_1+n k_3 \hat{\tau}_{13}+\frac{\mathrm{d} \hat{\tau}_{12}}{\mathrm{d} x_2}, \end{gather}$$

(2.32b) $$\begin{gather}0=\frac{1}{Re_\tau} \left( \frac{\mathrm{d}^2 }{\mathrm{d} x_2^2} -n^2 k_3^2 \right)^2 \hat{\psi}-k_3 \frac{\mathrm{d}}{\mathrm{d} x_2} ( \hat{\tau}_{33} - \hat{\tau}_{22}) + \left( \frac{\mathrm{d}^2}{\mathrm{d} x_2^2} +n^2 k_3^2 \right) \hat{\tau}_{23}, \end{gather}$$

(2.32c)\begin{align} -n k_3 \hat{\psi} \frac{\mathrm{d} \tilde{\nu}^{(0)}}{\mathrm{d}x_2} &= \frac{1}{\sigma} \left(\frac{1}{Re_\tau}+ \tilde{\nu}^{(0)}\right)\left(\frac{\mathrm{d}^2}{\mathrm{d} x_2^2} - n^2 k_3^2\right)\hat{\nu} + \frac{1}{\sigma} \frac{\mathrm{d}^2 \tilde{\nu}^{(0)}}{\mathrm{d}x_2^2}\hat{\nu} \nonumber\\ &\quad+ \frac{1}{\sigma}(2+2 c_{b2})\frac{\mathrm{d} \tilde{\nu}^{(0)}}{\mathrm{d} x_2}\frac{\mathrm{d} \hat{\nu}}{\mathrm{d} x_2}+c_{b1} \tilde{\nu}^{(0)} \widehat{\tilde{\mathcal{S}}} +c_{b1}\tilde{\mathcal{S}}^{(0)}\hat{\nu} \nonumber\\ &\quad- 2 \tilde{\nu}^{(0)} c_{w1} f_{w}^{(0)} \frac{\hat{\nu} d^{(0)}+\tilde{\nu}^{(0)} f(x_3)}{d^{(0)^3}} - c_{w1}f_{w}^{(1)}\left(\frac{\tilde{\nu}^{(0)}}{d^{(0)}}\right), \end{align}

along the wall-normal direction at each integer wavenumber $n = 1, 2, \ldots\,$. In these equations, the wall-normal profiles $\hat {\tau }_{ij}(x_2; n)$ are the components of the Reynolds stress tensor $\tau _{ij}^{(1)}$ obtained by substituting the expansion (2.31) into the definitions of the perturbations (2.15). This leads ultimately to a set of equations that only contains the quantities $\hat {u}_1(x_2; n)$, $\hat {\psi }(x_2; n)$ and $\hat {\nu }(x_2; n)$. Using the boundary conditions ((2.28)–(2.30)), these variables must satisfy

(2.33a)$$\begin{gather} \hat{u}_{1}(x_2={\pm} 1) ={-} f^n Re_\tau, \end{gather}$$

(2.33b)$$\begin{gather}\hat{\psi} (x_2={\pm} 1) = \mathrm{d}\hat{\psi}/\mathrm{d} x_2 (x_2={\pm} 1) = 0, \end{gather}$$

(2.33c)$$\begin{gather}\hat{\nu} (x_2={\pm} 1) ={-} f^n k. \end{gather}$$

Inspection of these boundary conditions and the governing equation shows that the wall topography affects the formation of secondary flows with three separate forcing terms. The first mechanism is mediated by the distance perturbation $d^{(1)}=-f(x_3)$. This term appears directly in the linearised transport equation of the eddy viscosity as a source term, suggesting that the topography modulation is felt throughout the domain as an alteration of the wall-normal development of the turbulent stresses. Crucially, spanwise heterogeneity of the topography produces a spanwise modulation of the eddy viscosity, i.e. of the Reynolds stress, which is known to be a source term in the transport equation of the turbulent kinetic energy (Barros & Christensen Reference Barros and Christensen2014; Hwang & Lee Reference Hwang and Lee2018). The second and third mechanisms are localised at the wall and are controlled by the inhomogeneous boundary conditions on the streamwise velocity and the perturbation eddy viscosity, respectively. The former produces a positive/negative velocity slip on the trough/crests of the modulation and generates a streaky motion with the associated streamwise velocity spanwise gradients. All these forcing terms are proportional to the strength of the coefficient $f^n$ in the series (2.4) characterising the surface geometry, showing the importance of fully characterising the spectral content of the wall topography.

The numerical solution of the system (2.32c) with the boundary conditions (2.33) is obtained by discretising the equations over $x_2 \in [-1, 1]$ using a Chebyshev collocation method. A spectral technique is technically not ideal for this problem, because $d^{(0)}$ has a sharp cusp at $x_2=0$. Nevertheless, we have observed that the spectral technique is robust in practice and provides accurate results when a sufficiently fine collocation grid is utilised. In the following calculations, we used no less than 202 collocation points, progressively increasing the resolution at the higher Reynolds numbers considered. The numerical code was also validated on sinusoidal channels using a nonlinear SA–QCR custom implementation in OpenFoam, with good agreement.

2.7. Reynolds-averaged solution in plane channels

The profiles of the mean streamwise velocity and the eddy viscosity of the plane channel appear in the first-order equations (2.32c) and are thus shown in this section. Profiles of these quantities were obtained by solving the SA equation (2.18) coupled with the streamwise momentum equation (2.1b) on a one-dimensional domain extending in the wall-normal direction using an in-house code. A linear Boussinesq approach is used, as this is sufficient in plane channels. The numerical code is based on a Chebyshev collocation discretisation and uses a Jacobian-free Newton–Krylov technique to solve the nonlinear coupled system of algebraic equations (Knoll & Keyes Reference Knoll and Keyes2004).

Mean streamwise velocity profiles obtained from the RANS solver at $Re_{\tau }=550$ and $Re_{\tau }=5200$ are shown in figures 3(a) and 3(b), respectively, as a function of the wall normal distance $x_2^+$ scaled by the viscous length (dashed red lines). These profiles extend to the channel midplane and are compared with the direct numerical simulation results of Lee & Moser (Reference Lee and Moser2015) (solid blue lines). The SA solution agrees well with the DNS data, especially in the log-layer, although higher velocities are observed in the buffer layer region. Profiles of the turbulent eddy viscosity $\nu _t^{(0)}$ are shown in figures 3(c) and 3(d), for the same Reynolds numbers. The eddy viscosity is extrapolated from the DNS simulation data by dividing the turbulent stress $- \overline {u_1'u_2'}$ with the wall-normal gradient of the streamwise velocity $\varGamma$. Good agreement with the DNS data is observed, although larger deviations are observed for $|x_2| \gtrsim 0.4$.

Figure 3. (a,b) Profiles of streamwise velocity (b,d) and of the turbulent eddy viscosity in plane channel from the SA model ($----$) and from the direct numerical simulation (DNS) (—-) of Lee & Moser (Reference Lee and Moser2015). Data is shown for $Re_{\tau }=550$ in panels (a,c) and $Re_{\tau }=5200$ in panels (b,d).

3.1. Organisation of secondary currents

The flow topology predicted by the linearised model is visualised in figure 4 for $\lambda _3 = 0.2, 0.5, 1, 2$ and $4$, in figure 4(a)–(d), respectively. Contours of the perturbation stream function (dashed contours for negative values) are reported. The colour map shows the wall-normal component $u_2^{(1)}$. Data at a large Reynolds number, $Re_\tau = 5200$, is reported as an illustrative example. Reynolds number effects are discussed later. A sketch of the harmonic topography is also reported below figure 4(e) for $\lambda _3 = 4$. For the symmetric configuration considered here, only data in the lower half of the channel is shown. The predicted secondary structure displays two counter-rotating vortices per period in the lower half of the channel. A similar flow organisation was recently observed by Vidal et al. (Reference Vidal, Nagib, Schlatter and Vinuesa2018) using DNS on wavy channels with antisymmetric walls. However, the present results refer to a symmetric channel where the vortices are confined to the half-channel height. On the contrary, in the simulations carried out by Vidal et al. (Reference Vidal, Nagib, Schlatter and Vinuesa2018) the vortices extend from the bottom to the upper wall. In addition, finite modulation amplitudes, with non-negligible convective effects, are considered Vidal et al. (Reference Vidal, Nagib, Schlatter and Vinuesa2018), unlike in the present case, where the modulation is infinitesimal. Nevertheless, in both cases the vortices flank the crest of the modulation and produce an upwelling motion above the crests. Conservation of mass through the channel then implies that a downwash is observed in the troughs of the topography. The height of the region affected by the secondary motion increases with $\lambda _3$ and, eventually, the vortices occupy the full half-height of the channel for $\lambda _3 \approx 1$. This topology persists from low periods up to $\lambda _3 \approx 6$, beyond which a large-scale flow reversal, where secondary currents rotate in the opposite direction and produce a downwash over the crest, is observed. This phenomenon is not related to the appearance of tertiary flows (Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015; Hwang & Lee Reference Hwang and Lee2018; Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2020) and might be a product of the turbulence model utilised in this paper that would not be observed in DNS or experiments. However, data to validate or disprove this behaviour for modulations with such large period does not seem to be available in the literature and further investigation is warranted.

Figure 4. Contours of the perturbation stream function $\psi ^{(1)}$ in the cross-plane $(x_2, x_3)$ at $Re_\tau =5200$ and varying wavelength: panel (a) $\lambda _3=0.2$; panel (b) $\lambda _3=0.5$; panel (c) $\lambda _3=1$; panel (d) $\lambda _3=2$; panel (e) $\lambda _3=4$. The stream function perturbation is limited to $[-1, 1]$ in panels (a) and (b), to $[-2, 2]$ in panels (c) and (d), to $[-0.2, 0.2]$ in panel (e), for a better representation of the flow structures. Dashed lines are used for negative values. The colour map of the wall-normal velocity perturbation (in units of the friction velocity and per unit of modulation amplitude) is also reported. For wavelengths smaller than $\lambda _3=4$, the flow topology for a single period is repeated to better display the evolution in size and strength of secondary flows. The topography from crest-to-crest is illustrated for the sake of clarity for $\lambda _3=4$ below panel (e).

3.2. Velocity profiles

Wall-normal profiles of the streamwise velocity component for $\lambda _3=0.2, 0.5, 1, 2, 4$ and $10$ are reported in figure 5 at $Re_\tau =550$ in figure 5(a) and at $Re_\tau =5200$ in figure 5(b). These profiles are localised at $x_3=0$, on the crest of the modulation. Velocity profiles at any other spanwise location, for example, over the trough, can be obtained by utilising the expansion (2.31) restricted to a single spanwise wavenumber mode. Given that the amplitude of the wall modulation (3.1) is proportional to $\lambda _3$, the velocity is first scaled by the wavelength and it should thus be interpreted as the flow response per unit amplitude of modulation expressed in terms of $h$. The velocity perturbation at the lower domain boundary is equal to $-Re_\tau$ due to the boundary conditions (2.28). We observe that the streamwise velocity is always negative, for all periods considered, corresponding to a LMP over the crest, as previously observed by other authors, for example Vanderwel & Ganapathisubramani (Reference Vanderwel and Ganapathisubramani2015) among others. The velocity perturbation decreases in magnitude when moving towards the channel centre. The depth of the disturbance increases with $\lambda _3$, as the secondary structures grow in size. The effect of the Reynolds number is moderate and consists of a slight increase in the momentum deficit when comparing corresponding profiles in figure 5(a) and figure 5(b).

Figure 5. Streamwise velocity perturbation per unit of wall-modulation amplitude $u_1^{(1)}(x_2, 0)/\lambda _3$ at $Re_\tau =550$ (a) and $Re_\tau =5200$ (b) from $\lambda _3=0.2$ to $\lambda _3=10$. The velocity profiles are extracted above the modulation crest. The velocity axis is restricted to $[-80, 0]$ for clarity, since the velocity perturbation at the lower domain boundary is $-Re_\tau$.

To better elucidate how the wall modulation alters the spatial structure of the streamwise velocity component, the quantity $\varGamma x_2$ is subtracted from the profiles of figure 5. This quantity attempts to capture the velocity perturbation produced by the shift in the mean velocity profile when the wall is displaced, particularly strong in the near-wall region but formally zero at the midplane. Results are reported in figure 6(a,b). It can be observed that the streamwise velocity perturbation is more pronounced in the near-wall region and relatively less in the channel centre. For short periods, this perturbation is positive, indicating that the near-wall flow over the crests moves faster than it would do over a flat wall. By contrast, for larger periods, the streamwise velocity perturbation is negative, initially in the vicinity of the wall and then gradually across the full channel half-width.

Figure 6. Profiles of the modified streamwise velocity perturbation $u_1^{(1)}(x_2, 0)/\lambda _3 - \varGamma x_2$ at $Re_\tau =550$, panel (a), and $Re_\tau =5200$, panel (b), at different spanwise wavelengths. Profiles are located over the modulation crest. In the figure: —$\circ$— $\lambda _3=0.2$; —$\square$— $\lambda _3=0.5$; —$\triangle$— $\lambda _3=1$; $---$ $\lambda _3=2$; —$\Diamond$— $\lambda _3=4$; —$\text{X}$— $\lambda _3=10$. In panels (c,d), the effect of turning on/off the QCR strain-stress model is shown for the same Reynolds numbers. Symbols are the same as in panel (a) but filled symbols are used for solutions at $C_{r1}=0$.

The change of sign with $\lambda _3$ suggests that two competing mechanisms are at play. The first mechanism is originated from the vertical ‘protrusion’ of the crests towards the midplane, causing higher velocity over the crests ‘exposed’ to the bulk of the flow. The second mechanism is the upwelling/downwelling motion introduced by the secondary structures. As shown in figure 4, these structures transport low momentum fluid from the near-wall region over the crest upwards towards the channel core, causing a local reduction of the flow velocity and vice versa over the troughs. When $\lambda _3$ is sufficiently large so that secondary currents are strong enough and they span a sufficiently large fraction of the channel, this second effect prevails and a low speed streak forms over the crests between the streamwise rolls, similarly to the optimal roll/streak configuration found in shear flows (del Álamo & Jiménez Reference del Álamo and Jiménez2006; Pujals et al. Reference Pujals, Garcìa-Villalba, Cossu and Depardon2009).

To better quantify the strength of these two competing mechanisms, we report in figure 6(c,d) the streamwise velocity profiles obtained from calculations where the QCR constant $C_{r1}$ is set to zero (filled symbols), corresponding to using a linear Boussinesq stress/strain relation. From a practical perspective, this is equivalent to ‘turning off’ secondary motions, so that only the first mechanism is active. The profiles are compared with the reference case at $C_{r1}=0.3$ (open symbols) and data for $\lambda _3=0.2, 1$ at the same Reynolds numbers of figure 6(a,b) is shown. When $C_{r1} = 0$, the velocity perturbation is always positive due to the protrusion of the crests into the bulk of the flow, as just mentioned, but when the QCR model is activated, negative velocities can be observed.

A further remark is that the profiles of the streamwise velocity show that the local wall shear stress perturbation can be significant. However, the perturbation of the spanwise-averaged shear stress predicted by the present linearised model is identically zero. In fact, from the expansions (2.31), it is easy to show that the wall shear stress is simply a harmonic function, with zero mean. A linear method cannot predict changes in spatially averaged quantities for flows obeying translational symmetries as in the present case, and second-order effects (i.e. large perturbations) must be taken into account to uncover, for example, how drag is affected by topography changes. However, having zero-mean velocities does not imply that spatial averaging is trivial for all other quantities. For instance, the spanwise averaged dispersive stresses often reported to characterise the secondary currents (Smith & McLean Reference Smith and McLean1977; Raupach & Shaw Reference Raupach and Shaw1982; Nikora et al. Reference Nikora, Goring, McEwan and Griffiths2001) can be non-zero. It is worth noting that the effective streamwise velocity profile resulting from the wall modulation is given, in our framework, by the sum of the flat channel profile and a small perturbation produced by the wall modulation. The streamwise velocity perturbation reported in figure 6 does not show any logarithmic behaviour, as it is the product of the two competing mechanisms discussed above. This might explain the strong distortion of the log-layer behaviour often observed in experiments or simulations of flows over heterogeneous surfaces (Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2018).

Profiles of the wall-normal and spanwise velocity components at $x_3=0$, on the crest of the modulation, and $x_3=\lambda _3/4$, respectively, are reported in figure 7 for the same Reynolds numbers and wavelengths considered in figure 6. As anticipated, in the lower half of the domain, the linearised RANS model predicts positive wall-normal velocities, indicating an upwash on the crest of the modulation and a downwash in the trough produced by secondary currents induced by the topography. For short periods, these effects are confined near the wall but the depth of the region influenced by this upwelling motion increases with the spanwise period up until $\lambda _3 \approx 1$, where the wall-normal motion involves the entire channel half-height. When the spanwise length scale is further increased, the wall-normal velocity decreases, as the vortical structures do not have additional space to grow.

Figure 7. Comparison of the profiles of the velocity components for different $\lambda _3$, at $Re_\tau =550$ in panels (a,c) and $Re_\tau =5200$ in panels (b,d). In panels (a,b) the wall-normal velocity $u_2^{(1)}(x_2,0)$ is plotted, in (c,d) the spanwise velocity $u_3^{(1)}(x_2,\lambda _3/4)$ is plotted. In the figure: —$\circ$— $\lambda _3=0.2$; —$\square$— $\lambda _3=0.5$; —$\triangle$— $\lambda _3=1$; $---$ $\lambda _3=2$; —$\Diamond$— $\lambda _3=4$; —$\text{X}$— $\lambda _3=10$.

For $\lambda _3=10$, the direction of rotation of the secondary vortices, occupying the entire channel half-height, is opposite to what is observed at the lower periods shown in figure 4 and downwash is now observed over the crest. Interestingly, the streamwise velocity profile in figure 5 is still negative, i.e. the flow displays a LMP associated with a downwash. Although this might appear to contradict established knowledge, the downwash velocity is relatively small in magnitude. It is argued that the negative streamwise velocity is purely a result of the contraction of the channel height over the crest (symmetric channels are considered), since the wall-normal velocity is too weak in this case to justify the observed change in the streamwise velocity. This effect would not be observed in an open flow like a boundary layer, where the wall modulation does not produce a contraction of the available cross-sectional area.

For the spanwise velocity, strong negative values are observed near the wall on the right-hand flank of the harmonic topography, producing a lateral jet-like motion towards the modulation crest. Generally, the negative velocity peak is larger than the peak of positive velocity, due to the confinement of the vortices near the wall (see figure 7c,d). The peak location varies only modestly with $\lambda _3$, but it gets closer to the wall and more intense at larger Reynolds numbers. For $\lambda _3=10$, the spanwise velocity profile shows a positive peak in the near-wall region due to the change of direction of rotation discussed previously.

3.3. Effect of wavelength and Reynolds number on the intensity of secondary flows

We now turn to investigating in more depth the effect of the wavelength $\lambda _3$ and of the Reynolds number on the strength of the secondary flows. For this purpose, we utilise the volume-averaged kinetic energy of the cross-flow components

(3.2)\begin{equation} \mathcal{K} = \frac{1}{4\lambda_3} \int_{{-}1}^{1} \int_{0}^{\lambda_3} [ u_2^{(1)}(x_2, x_3)^2 + u_3^{(1)}(x_2, x_3)^2] \,\mathrm{d} x_3 \,\mathrm{d} x_2, \end{equation}

to characterise the global amplitude of secondary flows. We also utilise the peak value of the perturbation stream function $\max _{x_2, x_3} |\psi ^{(1)}(x_2, x_3)|$, following Vidal et al. (Reference Vidal, Nagib, Schlatter and Vinuesa2018), to characterise the flow rate associated with the vortical flow and the peak wall-normal velocity $\max _{x_2, x_3} |u^{(1)}_2(x_2, x_3)|$. Results are reported in figure 8. In the left-hand panels, the dimensional spanwise period is scaled with the viscous length, i.e. $\lambda ^+_3 = \lambda _3 Re_\tau$, while in the right-hand panels the dimensional spanwise period is scaled with the outer scale $h$. Data for several Reynolds numbers, spanning the range $Re_\tau = 550$ to 5200 are reported. The vertical red lines denote regions where the predicted qualitative behaviour changes and are discussed later on. The key result is that the linearised model predicts two amplification peaks, indicating that the response of the turbulent wall-bounded flow to a harmonic topography modulation is stronger at preferential spanwise length scales. The location of these peaks is weakly dependent on the metric employed. In particular, the location of the first peak collapses when the wavelength is expressed in outer units to a value of $\lambda _3 \approx 1.54$. This peak is associated with large-scale vortical structures that occupy the entire half-height of the channel and produce a significant wall-normal transport through intense upwash/downwash regions on the crests/troughs of the modulation, as described in figure 4. On the other hand, the location of the second peak collapses when scaled in inner units, at $\lambda _3^+ \approx 45$. We have tested that the constant $C_{r1}$ of the QCR model does not affect the location of these peaks, but only their amplitude.

Figure 8. Intensity of secondary flows as a function of the spanwise wavelength. Different intensity metrics are compared. The panels (a,b) display the kinetic energy density $\mathcal {K}$, panels (c,d) the maximum stream function $\max _{x_2,x_3} | \psi ^{(1)}|$ and panels (e,f) the maximum of the wall-normal velocity $\max _{x_2,x_3} | u_2^{(1)}|$. The wavelength is scaled in inner units in (a,c,e) and outer units in (b,d,f). In the figure: —— $Re_\tau =550$; $-----$ $Re_\tau =1000$; —$\bullet$— $Re_\tau =3000$; —$\Diamond$— $Re_\tau =5200$. The vertical lines denote particular spanwise length scales where a change in the flow structure (flow reversal) is predicted.

This behaviour mirrors the predictions of transient growth analysis reported by del Álamo & Jiménez (Reference del Álamo and Jiménez2006) and Pujals et al. (Reference Pujals, Garcìa-Villalba, Cossu and Depardon2009) for plane channels. However, the location of the inner peak predicted in the present case is approximately half of the value found from the transient analysis, i.e. $\lambda _3^+ \approx 100$, which is predictive of the spanwise spacing of near-wall velocity streaks. It is also lower than what proposed by the conceptual model of Vidal et al. (Reference Vidal, Nagib, Schlatter and Vinuesa2018) who suggested an inner peak at $\lambda _3^+\approx 130$.

To further characterise these amplification peaks, profiles of the wall-normal and spanwise velocity components are reported in figure 9 for the outer peak (left-hand panels) and inner peak (right-hand panels). The Reynolds number varies from $550$ to $5200$. There are two major observations. Firstly, the present model predicts that the flow response to the surface modulation becomes, asymptotically, independent of the Reynolds number when scaled with the friction velocity, for both the inner and outer peaks. This is a major difference from transient growth analysis, where the energy gain increases with the Reynolds number. More importantly, this result is also in contrast with the findings of Vidal et al. (Reference Vidal, Nagib, Schlatter and Vinuesa2018) (and references therein) who performed DNS in wavy channels and showed that secondary flow velocities scaled by the bulk velocity are not sensitive to the Reynolds number if the Reynolds number is large enough to prevent marginally turbulent flow effects. The predictions of the present model can be attributed to fundamental properties of the SA model used in this study, as already indicated by Spalart et al. (Reference Spalart, Garbaruk and Stabnikov2018). In fact, the SA model is built in order to obtain a collapse of the eddy viscosity profile in the log-layer, where the transport equation (2.18) has solution $\tilde {\nu } = k x_2$, as well as in the outer layer. This implies that the eddy-viscosity profile, and thus the Reynolds stresses driving the formation of secondary flows of (2.15) are also, asymptotically, Reynolds number independent when scaled with the friction velocity.

Figure 9. Wall-normal (a,b) and spanwise (c,d) velocity profiles for the outer peak at $\lambda _3=1.54$ in (a,c) and inner peak at $\lambda _3^+=46.5$ in (b,d) for increasing Reynolds number. In the figure: —— $Re_\tau =550$; $-\cdot -\cdot -$ $Re_\tau =1000$; —$\bullet$— $Re_\tau =3000$; —$\Diamond$— $Re_\tau =5200$.

The second major observation is that the flow topology predicted by our model for the inner peaks is characterised by a downwash over the crest of the modulation, confined in the near-wall region ($x_2^+ < 30$). In fact, all quantities shown in figure 8 display two low amplification regions: one at $\lambda _3^+\approx 10^{2}$ and one at $\lambda _3\approx 6$, as denoted by the vertical lines in figure 8. At these spanwise length scales, a structural change in the topology predicted by the present model is observed, where a downwash is predicted over the crests of the modulation for either very large or very small wavelengths. While data for very large wavelengths, $\lambda _3 > 6$, does not appear to be presently available in the literature to compare our model with, the flow past surface corrugations at $\lambda _3^+ \approx 50$ is well known (e.g. Choi, Moin & Kim Reference Choi, Moin and Kim1993; Chu & Karniadakis Reference Chu and Karniadakis1993; Goldstein & Tuan Reference Goldstein and Tuan1998) and an upwash is typically observed over the crests of the corrugations. The origin of this discrepancy and of the difference in the location of the inner peak compared with what is found from transient growth analysis, can be attributed to the fact that the present RANS-based model is likely not able to capture correctly the nature of the interaction between the surface modulations and near-wall turbulent structures when these have commensurate lengths.

3.4. Large-scale motions in wall-bounded shear flows and secondary currents

The secondary currents predicted by the present model (figure 4) and their amplification as a function of the spanwise length scale (figure 8) are reminiscent of the optimal structures found with transient growth analysis in flat-wall turbulent channels by various authors (del Álamo & Jiménez Reference del Álamo and Jiménez2006; Pujals et al. Reference Pujals, Garcìa-Villalba, Cossu and Depardon2009). These smooth-wall analyses have demonstrated that the Navier–Stokes operator linearised around the turbulent mean profile and augmented with an eddy-viscosity term can support transient energy amplification at two specific spanwise length scales, scaling in inner and outer units, respectively. Specifically, streamwise-elongated roll-like motions introduced as initial conditions of the initial value problem develop into longitudinal streamwise streaks. These analyses have provided a formal description of the ubiquitous presence of near-wall streaky motions and large-scale structures in the outer layer of turbulent shear flows. The underlying mechanism is well known, i.e. the constructive interaction of nearly parallel stable eigenfunctions of the Orr–Sommerfeld–Squire equations (Butler & Farrell Reference Butler and Farrell1993). It was recently proposed by Chung et al. (Reference Chung, Monty and Hutchins2018) that a lateral variation of surface attributes may act a ‘phase lock’ to hold naturally occurring large-scale, outer layer structures around a fixed spatial location. Our linear operator analysis clarifies the relation between these structures and secondary currents, following the view expressed in, for example, Adrian & Marusic (Reference Adrian and Marusic2012). Spanwise heterogeneity of surface attributes may be interpreted as a steady forcing on the linearised operator, which then produces strong secondary motions with non-zero time average, locked on to the surface modulation. On the contrary, large-scale, outer-layer structures may be interpreted as the manifestation of the amplified response to unsteady turbulent velocity fluctuations, as proposed by del Álamo & Jiménez (Reference del Álamo and Jiménez2006) and Pujals et al. (Reference Pujals, Garcìa-Villalba, Cossu and Depardon2009), without spatial locking and hence with zero time-average. Hence, both types of structures may be interpreted as the independent manifestation of the amplification properties of the same linear operator, subjected to either steady or unsteady perturbations.

4.1. Effect of geometrical parameters

To elucidate the role of relevant parameters, we use the kinetic energy density of the cross-stream components, defined in (3.2), to characterise the global response in the cross-stream plane, and the peak stream function value $\max _{x_2,x_3}| \psi ^{(1)}(x_2, x_3)|$ to characterise the flow rate associated with the cross-stream motions (Vidal et al. Reference Vidal, Nagib, Schlatter and Vinuesa2018). These two quantities are reported in figure 10 as a function of the width $W$ and the gap $G$. Configurations at constant spacing $S = G + W = 1, 2, 3, \ldots$ lie on the white lines with slope $-1$. Note that configurations at constant duty cycle $DC = W/S$, considered as a relevant parameter in, for example, Castro et al. (Reference Castro, Kim, Stroh and Lim2021), lie on straight lines passing through the origin with slope $1/DC-1$. Results for $Re_\tau =5200$ are reported, since, as discussed in § 3.3, the SA-based RANS model predictions are asymptotically Reynolds number independent, and no qualitative changes to the following discussion arise when the response at other Reynolds numbers is examined.

Figure 10. Contours of the volume averaged kinetic energy of the cross-stream plane velocities $\mathcal {K}$ (a) and stream function peak value $\max _{x_2,x_3}|\psi ^{(1)}|$ (b) as a function of the gap $G$ and ridge width $W$. The Reynolds number is $Re_\tau =5200$. In panel (a), cases at constant duty cycle $DC= 0.25$ and $0.5$ are identified by the red lines. Cases at constant spacing $S=1, 2, 3, \ldots$ are identified by the white lines. Dashed lines identify cases at constant gap or width, with markers for configurations discussed later in the text.

Regardless of the metric used, secondary motions are weak for $S < 1$ and their strength peaks at $S \approx 1.34$, close to that obtained for sinusoidal walls and in agreement with predictions obtained in experiments on rectangular ridges (Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2020) but also for secondary flows developing over roughness strips (Chung et al. Reference Chung, Monty and Hutchins2018; Wangsawijaya et al. Reference Wangsawijaya, Baidya, Chung, Marusic and Hutchins2020) and streamwise arrays of roughness elements (Yang & Anderson Reference Yang and Anderson2018). The contours of the response have a preferential orientation whereby weaker changes in the response are observed when the spacing $S$ is held constant at the optimal value and $W$ and $G$ are varied. This occurs because such surfaces have a strong periodic component at the optimal length scale $S\approx 1.34$. This explains why many studies have identified this length scale as producing the largest response, despite significant differences in the ridge width/gap utilised. Nevertheless, our model predicts that the strongest response occurs when gap and width are equal, at $(W, G) \approx (0.67, 0.67)$, i.e. for relatively wide ridges.

For constant $G$ or $W$ equal to 0.67, significant amplification is observed when varying $W$ or $G$, respectively, along the two orthogonal red dashed lines in figure 10. Along these directions, one additional local peak is clearly visible at spacing $S\approx 2.8$, but several other (weaker) peaks occur at higher gaps or widths, at integer multiples of the optimal width $W \approx 0.67$. It is anticipated that these further peaks correspond to configurations with strong tertiary, quaternary of high-order structures (Hwang & Lee Reference Hwang and Lee2018) above/within the ridge/trough, confirming the conceptual model of Medjnoun et al. (Reference Medjnoun, Vanderwel and Ganapathisubramani2020). Nonetheless, further increasing the width (respectively, the gap) at constant gap (respectively, $W$) does not produce major changes in the strength of the response. These configurations tend asymptotically to the isolated ridge (respectively, gap) state, where the interaction between flow structures generated by adjacent ridges (respectively, gaps) can be neglected and the response is constant, regardless of the measure utilised.

A further important observation is that the response shows a symmetry with respect to the line $G = W$. The symmetry arises from the linear nature of the present analysis. For any surface configuration $(W, G)$, the flow topology in the trough is identical but with opposite flow direction to that on the ridge when $G$ and $W$ are swapped. The symmetry of the problem implies that the conceptual model developed by Medjnoun et al. (Reference Medjnoun, Vanderwel and Ganapathisubramani2020) speculating on the formation of tertiary structures over wide ridges can also be employed to describe the formation of tertiary structures in wide troughs induced by ‘virtual roughness element’ as proposed by Vanderwel & Ganapathisubramani (Reference Vanderwel and Ganapathisubramani2015).

Finally, the implication of the response maps of figure 10 is that, despite the spacing $S$ is a relevant length scale to characterise secondary flows, two surface parameters are required to characterise in a complete manner the strength of secondary currents. While many choices are possible, for example, $S$ and $W$ as in Hwang & Lee (Reference Hwang and Lee2018), Vanderwel & Ganapathisubramani (Reference Vanderwel and Ganapathisubramani2015) and Medjnoun et al. (Reference Medjnoun, Vanderwel and Ganapathisubramani2020), or $S/W$ and $S$ as in Castro et al. (Reference Castro, Kim, Stroh and Lim2021), using $G$ and $W$ is particularly convenient as (i) the response has a symmetry with respect to the line $G=W$ and (ii) these two parameters have similar roles when the flow organisation is considered, as we discuss in the next section.

A comparison with harmonic wall modulations is now performed to highlight similarities and differences between the two types of heterogeneity. To enable a direct comparison, we use the period $\lambda _3$ for harmonic walls and the spanwise spacing $S$ for rectangular ridges. Hence, for the latter case, where two geometrical parameters are strictly necessary, we restrict the analysis to specific configurations at duty cycle $DC=0.25,0.5$ and for a fixed width $W=0.67$, corresponding to the special identified by the additional lines in figure 10. Due to the symmetry of response discussed above, cases at duty cycles $DC^\prime$ greater than 0.5 have the same kinetic energy density and stream function peak of a cases at $DC=1 - DC^\prime$. For both types of heterogeneity, the same infinitesimal peak-to-peak modulation amplitude is used. The major similarity between the two types of heterogeneity is that the secondary structures have maximum intensity at a similar spanwise period, regardless of the quantity considered. For sinusoidal walls, the peak occurs at $\lambda _3 \approx 1.54$ while for rectangular ridges the peak is observed at $S \approx 1.34$. The effect of the duty cycle is only moderate, due to the orientation of the contours of the kinetic energy density in figure 10. The small difference between the two types of heterogeneity can be attributed to the fact that the first harmonic mode of the expansion (2.4) for the rectangular ridge geometry given by (4.1) is the largest and hence provides the dominant contribution to the overall response. However, the peak strength depends on the duty cycle. For instance, at $DC=0.25$, the strength of secondary flows is 75 % less intense than the vortices generated for the same spacing at $DC=0.5$. As discussed later on in § 4.2, secondary vortices do not have enough lateral space to develop at a sufficient strength when the ridges are too narrow, i.e. for small duty cycles. Interestingly, the peak value observed in figure 11 for rectangular ridges can be larger than what is observed for sinusoidal walls, provided the duty cycle is selected appropriately to intersect the high amplification region in the map of figure 10. Within the present linear framework, this effect can be explained to be the result of the constructive interference of the individual flow responses at all wavenumber modes defining the rectangular ridge geometry. A second major difference is that the response for rectangular ridges can exhibit a second peak at a larger spanwise period depending on the duty cycle. In particular, a peak at $S\approx 2.5$ can be observed at $DC=0.25$ and for $W=0.67$, while for $DC=0.5$, the second peak is shifted to a higher spacing ($S\approx 4$), associated with configurations at higher $W$ and $G$. It is anticipated that these secondary peaks correspond to ridge configurations where tertiary flows, developing at the centre of the ridge (or troughs), have the maximum strength.

Figure 11. Comparison of the kinetic energy density $\mathcal {K}$ (a) and the maximum of the stream function $\max _{x_2,x_3} |\psi ^{(1)}|$ (b) at $Re_\tau =5200$ as a function of the periodicity $S$. For the rectangular ridges, the quantities are obtained for constant $W=0.67$ and for $DC=0.25$ and $0.5$.

4.2. Topology of secondary flows

Based on the symmetry highlighted from the response maps, we now show how the parameters $W$ and $G$ affect the organisation of secondary flows. We consider flows at $Re_\tau =5200$, at which the response has saturated to its high-$Re$ asymptotic state. In figure 12, contours of the perturbation stream function are reported together with colour maps of the wall-normal velocity perturbation for configurations at constant gap $G=0.67$ and at varying $W = 0.3, 0.67, 1.5$ and $2$ (see triangles in figure 10). The black lines at $x_2 = -1$ define the locations of the ridges. Note that the fields are spanwise periodic and only half of the ridge is shown, as the ridges are centred at $x_3 = 0$.

Figure 12. Flow organisation for $G=0.67$ and width $W=0.3$ (a), 0.67 (b), 1.5 (c) and 2 (d). Results for $Re_{\tau }=5200$ are shown. Contours of the perturbation stream function $\psi ^{(1)}$ between $-2$ and 2 are shown. The dashed lines indicate negative stream function values. The colour map of the perturbation wall-normal velocity component $u_2^{(1)}$ is also reported in the lower half of the channel. The ridges are sketched on the bottom line using bold lines. Note that the ridges are centred at $x_3=0$ and 1 and the fields are spanwise periodic.

Starting from $W=0.3$, the linear model predicts counter-rotating vortical structures elongated in the wall-normal direction and occupying the entire half-width of the channel. These structures are locked in proximity of the ridge edges where the surface discontinuity acts as a strong source term. A downwash inside the troughs and an upwash above the ridges in proximity of the edge is observed. The maximum intensity of these vertical motions at $W=0.67$ is approximately 15$u_\tau$, per unit of ridge height. This means that for a peak-to-peak ridge height of $0.09$ (in units of the boundary layer thickness) as in case HS6 (Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2020) for $Re_{\tau }=3239$, the predicted peak vertical velocity is $3\,\%$ of the bulk velocity, which agrees with the experimental data ($2\,\%$). However, the comparison between the case HS6 of Medjnoun et al. (Reference Medjnoun, Vanderwel and Ganapathisubramani2020) and the present model can only be qualitative because of structural differences between channel flows and boundary layers. In particular, wall-normal secondary motions in a boundary layer are not confined or blocked by the upper wall (or symmetry plane), but can develop freely. This effect is likely more important when the secondary structures have size comparable to the shear layer thickness, and secondary motions can reach the outer edge of the turbulent shear flow. For a short $W=0.3$ (figure 12a), the vortical structures compete for the available space over the ridge, push each other towards the gap centre and are highly elongated in the vertical direction. For $W = 0.67$ (figure 12b), the vortices can now fully extend towards the ridge centre. For $W = 1.5$ (figure 12c), there is sufficient space over the ridge for tertiary flows to emerge in the region immediately above the ridge. In this condition, downwash is observed over the ridge centre. This, however, is associated with a reduction in the strength of the upwash in the vicinity of the edges, strongest at $W=0.67$. Tertiary vortical structures are initially weak but gain strength at $W\approx 2.1$, where they can fully extend to the channel midplane. The strength of the downwash velocity at the ridge centre for $W=2.1$ relative to the downwash velocity over the gap is significant. This is likely exacerbated by confinement effects in the channel, in which the spanwise-averaged vertical mass transport operated by secondary currents is necessarily zero. In boundary layers, no such constraint would exist. Although not shown here, for $W > 3.5$ a further reorganisation is observed, where weak quaternary vortical structures emerge near the ridge centre ($x_3=0$), producing a weak upwash motion.

One important remark is that the present linearised model does not capture correctly flow features observed in the immediate vicinity of the ridge such as, for instance, recirculation regions induced by strong spanwise motions over the ridge top, frequently observed in direct numerical simulations (Hwang & Lee Reference Hwang and Lee2018; Castro et al. Reference Castro, Kim, Stroh and Lim2021). The wall-normal extent of these regions is (i) strongly influenced by the ridge geometry (rectangular, circular, etc.) and (ii) likely scaling with the ridge height, which is always finite in experiments and simulations. In the present linear model, the ridge height is infinitesimal and only large-scale flow features developing far away from the surface are likely to be captured correctly. Localised near-wall effects produced by a finite ridge height and contributing less prominently to the alteration of vertical transport phenomena are unlikely to be accounted for. Nevertheless, the model predicts structures with similar characteristics to those observed in many other studies, where the ridges protrude into the log region. It could be argued that, for a shallow surface modulation, all the mean flow quantities (e.g. the Reynolds stresses) develop in the wall-normal direction according to the same law of the wall, as if the wall was flat. The lateral variation of the origin of these profiles produced by the modulation then produces at any distance from the wall spanwise gradients of the Reynolds stresses, i.e. the required source terms in the streamwise vorticity equation (2.7). This mechanism might be at play regardless of the height of the modulation, although for large protrusions other mechanism became relevant, for example, the wall-normal deflection of spanwise velocity fluctuations (see e.g. Hwang & Lee Reference Hwang and Lee2018).

For completeness, the evolution of the flow organisation for a constant $W=0.67$ as the gap $G$ increases is shown in figure 13. These configurations correspond to the squares in figure 10, and parallel the configurations shown in figure 12. For $G=0.3$ (figure 13a), the vortical structures compete for the available space over the gap and push each other away towards the ridge. As the gap is further increased to $G=1.5$ and then 2, tertiary structures form in the centre of the trough producing vertical velocities weaker than the velocity induced by the secondary structures over the ridge. As anticipated, this behaviour was described by Vanderwel & Ganapathisubramani (Reference Vanderwel and Ganapathisubramani2015), who observed that, when the spacing is large enough, an additional upwelling motion is generated at the centre of the trough as if a ‘virtual’ ridge element was placed between physical ridges.

Figure 13. Flow organisation for $W=0.67$ and gap $G=0.3$ (a), 0.67 (b), 1.5 (c) and 2 (d). Results for $Re_{\tau }=5200$ are shown. Contours of the perturbation stream function $\psi ^{(1)}$ between $-2$ and 2 are shown. The dashed lines indicate negative stream function values. The colour map of the perturbation wall-normal velocity component $u_2^{(1)}$ is also reported is also reported in the lower half of the channel. The ridges are sketched on the bottom line using bold lines. Note that the ridges are centred at $x_3=0$ and 1 and the fields are spanwise periodic.

In conclusion, the secondary structures shown in figures 12 and 13 and predicted by the present model are similar in size and organisation, especially in the region closest to the wall, to the secondary currents observed over strip-type roughness both numerically (Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015; Chung et al. Reference Chung, Monty and Hutchins2018) and experimentally (Hinze Reference Hinze1967, Reference Hinze1973; Mejia-Alvarez & Christensen Reference Mejia-Alvarez and Christensen2010, Reference Mejia-Alvarez and Christensen2013; Barros & Christensen Reference Barros and Christensen2014; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015; Bai, Kevin & Monty Reference Bai, Kevin and Monty2018; Forooghi et al. Reference Forooghi, Yang and Abkar2020; Wangsawijaya et al. Reference Wangsawijaya, Baidya, Chung, Marusic and Hutchins2020). The similarity may be due to the fact that for strip-type heterogeneity, the wall-normal extent of the surface perturbation is localised, loosely speaking, within the roughness height. When scaled in outer units, this height is often much smaller than the typical (and finite) ridge height used in previous work. However, the key difference is that the direction of rotation of secondary structures over strip-type roughness is opposite to what is observed for ridge-type heterogeneity and a downwash motion is predicted over the high roughness region (Stroh et al. Reference Stroh, Schäfer, Frohnapfel and Forooghi2020b; Schäfer et al. Reference Schäfer, Stroh, Forooghi and Frohnapfel2022). The strong similarity suggests that the present linearised framework may also be adapted to study secondary currents produced by other types of complex surfaces, where the lateral variation of the surface attributes may be modelled accurately by suitable boundary conditions that capture the physics of the flow-surface interaction. In addition to the strip-type roughness case just discussed, examples of such surfaces include superhydrophobic surfaces (Busse & Sandham Reference Busse and Sandham2012; Türk et al. Reference Türk, Daschiel, Stroh, Hasegawa and Frohnapfel2014; Stroh et al. Reference Stroh, Hasegawa, Kriegseis and Frohnapfel2016), porous surfaces (Abderrahaman-Elena & Garcìa-Mayoral Reference Abderrahaman-Elena and Garcìa-Mayoral2017; Efstathiou & Luhar Reference Efstathiou and Luhar2018; Rosti, Brandt & Pinelli Reference Rosti, Brandt and Pinelli2018; Bottaro Reference Bottaro2019) or surface lateral variations of the heat flux (Salesky, Calaf & Anderson Reference Salesky, Calaf and Anderson2022).

4.3. Velocity profiles over rectangular ridges

For a better characterisation of the flow structures developing over the ridges, the wall-normal velocity profile at five different locations between the left-hand edge of the ridge and its centre are reported in figure 14. The sketch the sketch beneath figure 14(c,d) illustrates the location where profiles are extracted. The velocity profiles are obtained at $Re_\tau =5200$ and for $DC=0.5$. The ridge width $W$ varies from $0.3$ in figure 14(a) to $2$ in figure 14(d). For the two smaller widths, the velocity is always positive corresponding to the upwash region of figure 12(a,b). For the optimal configuration $W=0.67$, where the vortical structures fit the available space without significant lateral distortion, the velocity above the ridge edge is small. This contrasts with experimental/numerical observations (e.g. Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2020), where intense upwards motions are often observed at the ridge edge. This might be the result of the finite ridge height in these cases, which produces a wall-normal deflection of the spanwise velocity fluctuations. The present linear model, with infinitesimal ridges, does not capture this deflection although it does predict an intense spanwise motion at the ridge edge. Tertiary flows occur for $W = 1.5$ and $2$, where the $u_2^{(1)}$ profile at the centre of the ridges (light red) displays a negative value. However, the negative peak is approximately 50 % less intense than the positive one, although the size of the region interested by the wall normal motion is similar (figure 14c,d). The intensity of the tertiary flows slightly increases with $W$ while the secondary flows appear to be unaffected. This suggests a saturation of the secondary flows, while higher-order structures occur at the ridge centre.

Figure 14. Wall-normal velocity profiles at different locations along the ridge. The Reynolds number is $Re_\tau =5200$ and duty cycle is $DC=0.5$. The width W varies from $0.3$ (a) to 2 (d). The profile locations are also reported in the sketch beneath panels (c,d) using a different colour gradation.

To better characterise the intensity and direction of the secondary flows as a function of the spanwise location, the quantity

(4.2)\begin{equation} \mathcal{I}_2^p(x_3)=\int_{{-}1}^{0} u_2^{(1)^p}(x_2,x_3)\, \mathrm{d} x_2, \end{equation}

is now introduced. We first discuss the case for $p=2$. Results are reported in figure 15 for $DC=0.5$ and $W$ varying from 0.67 to 3. To align all profiles, we use the spanwise coordinate $x_3+W/2$, so that the ridge edge is always located at $0$ while the ridge centre corresponds to $W/2$. For the smaller widths, the quantity $\mathcal {I}_2^2(x_3)$ is quite intense and a single peak produced by the secondary structures locked on the ridge edges is observed, approximately at the ridge centre. Increasing the ridge width to $W=1.5$, the quantity $\mathcal {I}_2^2(x_3)$ shows two peaks: the first in proximity to the ridge edge (with smaller magnitude than at the optimal width $W=0.67$) and the second at the ridge centre, characterising the strength of tertiary flows. When the width is further increased, secondary flows develop fully and only moderate effects on their strength near the ridge edge is observed. Major differences are still observed for tertiary flows at the ridge centre, although the expectation is that such differences would eventually vanish as the ridge width is increased further.

Figure 15. The quantity $\mathcal {I}_2^2(x_3)$ for $Re_\tau =5200$ and duty cycle $DC=0.5$. The ridge width $W$ varies from 0.67 to 3. In particular: —$\circ$— $W=0.67$; —$\square$— $W=0.8$; $---$ $W=1$; —$\Diamond$— $W=1.5$; —$X$— $W=2$; —— $W=2.5$; $--- \triangle ---$ $W=3$.

4.4. Analysis of the wall-normal velocity direction over the ridge centre

To better visualise the region of parameter space where the present linearised model predicts a large-scale change in the flow direction above the centre of the ridges, we use the quantity $\mathcal {I}_2^1(0)$ to quantify the average, or ‘bulk’, wall-normal flow direction at $x_3 = 0$, as a function of the gap $G$ and the width $W$. Results are reported in figure 16. The linearised model indicates that the bulk wall-normal velocity becomes negative for $W \gtrsim 1.2$, with a moderate effect of the gap. The maximum average velocity occurs for $W\approx 0.5, G\approx 0.75$, indicating that optimising the intensity of secondary currents using the strength of the average wall-normal velocity yields narrower ridges than what suggested by using the integral perturbation energy or the stream function peak. The bulk velocity turns positive again for $W \gtrsim 2.8$ when the ridge is wide enough to support the formation of quaternary structures. The quantity $\mathcal {I}_2^{1}(0)$ alone, however, might not be sufficient to capture the change in flow direction that is often observed in the proximity of the obstacle (Castro et al. Reference Castro, Kim, Stroh and Lim2021). The onset of this change is thus also indicated in the figure, by tracing the set of points (solid black line) in parameter space where the wall-normal velocity at the centre of the ridge first changes sign. Due to the aforementioned symmetries, large-scale or incipient change in flow direction in the troughs, observed, for example, by Vanderwel & Ganapathisubramani (Reference Vanderwel and Ganapathisubramani2015) can be characterised by swapping the role of $G$ and $W$ and inverting the sign of $\mathcal {I}_2^{1}(0)$ (computed at $x_3 = S/2$, in the centre of the trough). The region where a change in flow direction is predicted in the troughs by the present model is shown as a dashed black line. The model predicts that the difference between the average and incipient change of flow direction is minimal. However, this difference might be more pronounced for finite height ridges, where the flow topology near the ridge is more complicated than what can be captured by the present linear model.

Figure 16. Colour map of the quantity $\mathcal {I}_2^{1}(0)$ as a function of the gap $G$ and width $W$, for $Re_\tau = 5200$. Configurations studied in the recent literature are denoted by symbols (VG2015 for Vanderwel & Ganapathisubramani (Reference Vanderwel and Ganapathisubramani2015), MVG2020 for Medjnoun et al. (Reference Medjnoun, Vanderwel and Ganapathisubramani2020), CKS2021 for Castro et al. (Reference Castro, Kim, Stroh and Lim2021) and HL2018 for Hwang & Lee (Reference Hwang and Lee2018)). Closed symbols denote configurations where downwash over the ridge has been observed. The black lines delimit the regions where the linear model predicts incipient change of flow direction at the midpoint over the ridge (solid line) and at the centre of the trough (dashed line).

Data from recent numerical and experimental investigations that have considered streamwise rectangular ridges are also reported in figure 16. As a note of caution, most of these cases (with the exception of Castro et al. (Reference Castro, Kim, Stroh and Lim2021) who considered channel flows) are extracted from studies of secondary flows in boundary layers. As discussed previously, secondary flow structures originating in different flow types might display significant topological differences and the following analysis should be regarded as qualitative. Interestingly, a large fraction of experiments and numerical simulations available in the literature is focused on the region of narrow width, relatively far from the optimal configuration predicted by the present model. In the figure, closed symbols denote configurations where a large-scale change in the flow direction (and not simply in the neighbourhood of the ridge) was observed above the ridge or in the trough. These are the case HS6 from Medjnoun et al. (Reference Medjnoun, Vanderwel and Ganapathisubramani2018) and P24S12 from Hwang & Lee (Reference Hwang and Lee2018), where a downwelling motion is observed above the ridge at large distance from the wall, and $S/\delta =1.76$ from Vanderwel & Ganapathisubramani (Reference Vanderwel and Ganapathisubramani2015), where upwelling is measured in the trough. For case HS6, the present model predicts a positive net wall-normal velocity, in contrast to experimental evidence. Inspection of the velocity field for this case in Medjnoun et al. (Reference Medjnoun, Vanderwel and Ganapathisubramani2018) shows that the Reynolds-averaged vortical structures are smaller in size (in both directions) and less coherent than what is predicted by the present model. In turn, this would increase the space available for fluid to reverse its direction. This difference might be due to the different flow type (boundary layer in Medjnoun et al. (Reference Medjnoun, Vanderwel and Ganapathisubramani2018) and channel flow in the present work) or to the finite ridge height in experiments.