Linearised Reynolds-averaged predictions of secondary currents in turbulent channels with topographic heterogeneity

Abstract

A rapid predictive tool based on the linearised Reynolds-averaged Navier–Stokes equations is proposed in this work to investigate secondary currents generated by streamwise-independent surface topography modulations in turbulent channel flow. The tool is derived by coupling the Reynolds-averaged momentum equation to the Spalart–Allmaras transport equation for the turbulent eddy viscosity, using a nonlinear constitutive relation for the Reynolds stresses to capture correctly secondary motions. Linearised equations, describing the steady flow response to arbitrary surface modulations, are derived by assuming that surface modulations are shallow. Since the equations are linear, the superposition principle holds and the flow response induced by an arbitrary modulation can be obtained by combining appropriately the elementary responses obtained over sinusoidal modulations at multiple spanwise length scales. The tool permits a rapid exploration of large parameter spaces characterising structured surface topographies previously examined in the literature. Here, channels with sinusoidal walls and with longitudinal rectangular ridges are considered. For sinusoidal walls, a large response is observed at two spanwise wavelengths scaling in inner and outer units respectively, mirroring the amplification mechanisms in turbulent shear flows observed from transient growth analysis. For longitudinal rectangular ridges, the model suggests that the analysis of the response and the interpretation of the topology of secondary structures is facilitated when the ridge width and the gap between ridges are used instead of other combinations proposed in the literature.

1. Introduction

When a wall-bounded turbulent flow develops over a surface with heterogeneous attributes, for example, with lateral variations of the topography or of the roughness properties, secondary currents emerge in the form of coherent streamwise-aligned vortices. These flows, named by Prandtl as secondary flows of the second kind (Prandtl 1952), have attracted significant interest since the first experiments in rectangular ducts with heterogeneous rough surfaces conducted by Hinze (1967, 1973). In fact, these flows are highly relevant in many industrial and environmental applications, where aerodynamic surfaces are rarely smooth and homogeneous. Despite being relatively weak, with velocities of a few per cent of the external velocity scale, these currents can alter natural wall-normal transport properties of wall-bounded turbulent flows (Volino, Schultz & Flack 2011; Mejia-Alvarez & Christensen 2013; Vanderwel & Ganapathisubramani 2015; Hwang & Lee 2018; Medjnoun, Vanderwel & Ganapathisubramani 2020; Zampiron, Cameron & Nikora 2020) and can thus increase friction and heat transfer (Stroh et al. 2020a) and modify the performance of aerodynamic surfaces (Mejia-Alvarez & Christensen 2013; Barros & Christensen 2014).

Broadly speaking, the heterogeneity can be distinguished between topographical variations, i.e. alternating regions of high/low relative elevation (Hwang & Lee 2018; Medjnoun, Vanderwel & Ganapathisubramani 2018; Medjnoun et al. 2020; Castro et al. 2021) and skin-friction variations, where the local wall shear stress varies as a consequence of changes in the surface attributes, such as the roughness properties (Barros & Christensen 2014; Chung, Monty & Hutchins 2018; Forooghi, Yang & Abkar 2020; Stroh et al. 2020b) or over superhydrophobic surfaces (Türk et al. 2014; Stroh et al. 2016). Combinations of these two have also been considered, (e.g. Vanderwel & Ganapathisubramani 2015; Yang & Anderson 2018; Stroh et al. 2020b). However, in all cases, the flow topology observed above such surfaces is characterised by alternating high-momentum pathways (HMP), corresponding to a downwash motion, and low-momentum pathways (LMPs), correlated to an upwash motion, as observed by Mejia-Alvarez & Christensen (2013) and Willingham et al. (2014). This alternance of HMP and LMPs is observed both experimentally (Barros & Christensen 2014; Anderson et al. 2015; Vanderwel & Ganapathisubramani 2015) and numerically (Stroh et al. 2016; Chung et al. 2018). Even though the instantaneous field is highly complex (Vanderwel et al. 2019), these motions are associated, in a Reynolds-averaged sense, with large-scale streamwise vortical structures, driven by a turbulent torque produced by lateral variations of the (anisotropic) Reynolds stress tensor (Perkins 1970; Bottaro, Soueid & Galletti 2006).

The lateral organisation and intensity of HMP and LMPs and of the associated vortical structures is often discussed in relation to a characteristic spanwise length scale of the heterogeneity, such as the spacing between longitudinal ridges or the width of roughness strips or patches of superhydrophobic surface. Many authors have performed parametric studies and have demonstrated that secondary motions are most intense when this characteristic length scale is of the order of the thickness of the turbulent shear layer (Vanderwel & Ganapathisubramani 2015; Chung et al. 2018; Yang & Anderson 2018). However, significant changes in the flow topology, for example, the appearance of tertiary flows, have also observed when other surface parameters are varied, such as the width of the ridges or the ridge geometry. In an effort to quantify these aspects, Medjnoun et al. (2020) introduced the ratio between the cross-sectional areas above and below the mean surface height as the key surface parameter that distinguishes different topographies and the observed flow structure. They showed that the circulation of the time-averaged vortical structures is proportional to this ratio. However, a complete description of how surface characteristics influence the structure and intensity of secondary motions is still lacking. In fact, this endeavour has been hindered by the high-dimensional nature of the parameter space that characterises heterogeneous surfaces, which is costly to fully explore using experiments or scale-resolving simulations.

The overarching aim of this work is to develop a rapid predictive tool to aid the exploration of such spaces. In this paper, we restrict our attention to surfaces with lateral variations of the topography, but extensions to other types of heterogeneity are possible. The proposed tool is based on the steady linearised Reynolds-averaged Navier–Stokes (RANS) equations, augmented by a turbulent eddy-viscosity term. These equations have been used in past work to clarify key mechanisms of wall bounded turbulence. For instance, the characteristic spanwise length of near-wall streaks and large-scale motions in turbulent shear flows is well captured by the energy amplification properties of the Orr–Sommerfeld–Squire equations (del Álamo & Jiménez 2006; Pujals et al. 2009; Hwang & Cossu 2010). Luchini & Charru (2010) and Russo & Luchini (2016) used linearised RANS equations to model flows over undulated bottoms or to examine the response to volume forcing. Meyers, Ganapathisubramani & Cal (2019) utilised the linearised RANS equations to predict the decay rate of dispersive stresses associated with secondary motions in the outer-layer region. Unlike in some of the previous literature, where simple analytical profiles for the eddy viscosity have been used, here the Reynolds-averaged momentum equations are coupled with the Spalart–Allmaras (SA) transport equation for the turbulent eddy viscosity (Spalart & Allmaras 1994), to capture more faithfully the variable topography. Linearised equations are then derived by assuming that the topography is shallow when compared with any inner or outer length scale. For shallow modulations, the nonlinear convective terms are negligible and arbitrary surface topographies can be modelled using inhomogeneous linearised boundary conditions (Luchini 2013). Using these equations, the response of the shear flow to an arbitrary, spectrally complex surface topography can be obtained by applying the superposition principle, i.e. by appropriately combining the elementary responses obtained for all the harmonic components defining the given surface. Channels with sinusoidal walls (Vidal et al. 2018) and with longitudinal rectangular ridges are considered in this paper as two paradigmatic configurations that have received significant attention in the recent literature.

The modelling technique and the linearisation of the governing equations is discussed in § 2. The approach is first applied to sinusoidal modulations in § 3, to clarify the fundamental role of the spanwise length scale on the strength and structure of secondary motions. With this insight, channels with rectangular ridges are considered in § 4. Finally, conclusions are reported in § 5.

2. Methodology

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 (2016), 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. (2018), 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 2015; Hwang & Lee 2018; Medjnoun et al. 2020) or open channel flows (Zampiron et al. 2020), 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 2000), 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. (2019), 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. 2015; Hwang & Lee 2018; Medjnoun et al. 2020), 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 1970; Speziale 1982; Bottaro et al. 2006). 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 2006; Pujals et al. 2009) 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 1991; Speziale, Sarkar & Gatski 1991; Chen, Lien & Leschziner 1997). Here we use the quadratic constitutive relation (QCR) nonlinear model introduced by Spalart (2000), which contains simple terms proportional to the product of the rotation and the strain tensors. This model was recently utilised by Spalart, Garbaruk & Stabnikov (2018) 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 (2000) 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 2006; Pujals et al. 2009) have often used analytical eddy-viscosity profiles (Cess 1958; Reynolds & Hussain 1972). 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 1994), 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 1994; Herring & Mellor 1968) 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 1994), 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 2012; Luchini 2013). 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 1994), 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 1994). 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 2014; Hwang & Lee 2018). 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 2004).

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 (2015) (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 (2015). Data is shown for $Re_{\tau }=550$ in panels (a,c) and $Re_{\tau }=5200$ in panels (b,d).

3. Secondary flows in sinusoidal channels

Secondary flows in symmetric channels with sinusoidal walls (see figure 1a) are now considered to elucidate the fundamental role of the spanwise length scale on the generation of secondary flows. This insight can then be used to analyse surfaces with complex spatial characteristics (Barros & Christensen 2014; Anderson et al. 2015). We consider modulations expressed by the cosine law

(3.1)\begin{equation} f(x_3) = \lambda_3 \cos(k_3 x_3). \end{equation}

Scaling the amplitude with the period $\lambda _3$ ensures that the aspect ratio of the modulation (peak-to-peak amplitude to spanwise length scale) remains constant, i.e. we follow the shallow-roughness limit introduced in Luchini (2013).

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. (2018) 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. (2018) 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. (2018), 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 2015; Hwang & Lee 2018; Medjnoun et al. 2020) 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 (2015) 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 2006; Pujals et al. 2009).

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 1977; Raupach & Shaw 1982; Nikora et al. 2001) 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. 2018).

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. (2018), 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 (2006) and Pujals et al. (2009) 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. (2018) 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. (2018) (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. (2018). 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 1993; Chu & Karniadakis 1993; Goldstein & Tuan 1998) 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 2006; Pujals et al. 2009). 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 1993). It was recently proposed by Chung et al. (2018) 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 (2012). 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 (2006) and Pujals et al. (2009), 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. Secondary flows above rectangular ridges

Secondary flows above rectangular ridges are now considered. As shown in figure 1(b), the geometrical parameters considered are the spacing between the ridges $S$ and the ridge width $W$. The gap between the elements is $G=S-W$. Linearised flow solutions in this geometry are obtained wavenumber-by-wavenumber as discussed in § 2.6. Except for very near the wall, the solution is smooth and the Fourier expansion (2.31) converges rapidly. To improve the convergence of our spectral code in the near-wall region, the discontinuous wall geometry is approximated by the smooth function

(4.1)\begin{equation} f(x_3)=\frac{1}{\arctan(\alpha)} \arctan\left( \alpha \left[ \cos(k_3 x)-\cos \left (k_3 \frac{W}{2} \right) \right]\right), \end{equation}

where $\alpha$ is used to round the corners of the ridges and to increase the roll-off of the coefficients $f^n$ of its cosine series (2.4). Here, $\alpha$ is chosen so that $\mathrm {d}f / \mathrm {d} x_3(W/2) = 2 \times 10^4$. The surface geometry is then discretised with at least $150$ cosine waves, ensuring that the ratio $|\,f^1/f^{150}|$ is no less than 300. We have repeated some calculations at finer resolutions, and no appreciable change in the structure of large-scale motions developing over this geometry has been observed.

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. 2018). 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. (2021), 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. 2020) but also for secondary flows developing over roughness strips (Chung et al. 2018; Wangsawijaya et al. 2020) and streamwise arrays of roughness elements (Yang & Anderson 2018). 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 2018) above/within the ridge/trough, confirming the conceptual model of Medjnoun et al. (2020). 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. (2020) 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 (2015).

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 (2018), Vanderwel & Ganapathisubramani (2015) and Medjnoun et al. (2020), or $S/W$ and $S$ as in Castro et al. (2021), 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. 2020) 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. (2020) 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 2018; Castro et al. 2021). 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 2018).

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 (2015), 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. 2015; Chung et al. 2018) and experimentally (Hinze 1967, 1973; Mejia-Alvarez & Christensen 2010, 2013; Barros & Christensen 2014; Anderson et al. 2015; Bai, Kevin & Monty 2018; Forooghi et al. 2020; Wangsawijaya et al. 2020). 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. 2020b; Schäfer et al. 2022). 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 2012; Türk et al. 2014; Stroh et al. 2016), porous surfaces (Abderrahaman-Elena & Garcìa-Mayoral 2017; Efstathiou & Luhar 2018; Rosti, Brandt & Pinelli 2018; Bottaro 2019) or surface lateral variations of the heat flux (Salesky, Calaf & Anderson 2022).

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. 2020), 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. 2021). 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 (2015) 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 (2015), MVG2020 for Medjnoun et al. (2020), CKS2021 for Castro et al. (2021) and HL2018 for Hwang & Lee (2018)). 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. (2021) 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. (2018) and P24S12 from Hwang & Lee (2018), where a downwelling motion is observed above the ridge at large distance from the wall, and $S/\delta =1.76$ from Vanderwel & Ganapathisubramani (2015), 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. (2018) 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. (2018) and channel flow in the present work) or to the finite ridge height in experiments.

5. Conclusions

A rapid tool for the prediction of secondary currents developing in turbulent channels with streamwise-independent surface modulations has been presented. The approach is based on the linearisation of the steady RANS equations, coupled to the SA equation for the transport of the turbulent eddy viscosity. The linearisation of these equations is based on the assumption that the surface modulation is small when compared with any relevant geometric or physical length scale. The influence of the surface modulation is then modelled using inhomogeneous boundary conditions for the streamwise velocity component and the turbulent eddy viscosity. Because of the linearity, the superposition principle applies and the flow response induced by an arbitrary surface with spectrally complex topography can be obtained by appropriately combining the elementary responses to harmonic modulations at each spanwise wavelength.

The computational efficiency of the tool allows large parameter spaces characterising complex surfaces to be explored at little cost. In this paper, two canonical surface configurations are studied, namely, harmonic modulations and rectangular ridges. For harmonic modulations, characterised by a single spanwise length scale, the wavelength $\lambda _3$, the turbulent shear flow is found to have the largest response at two spanwise wavelengths, scaling in inner and outer units, respectively. The outer peak is found at $\lambda _3 \approx 1.54$, in units of the half-channel width, and corresponds to large-scale secondary vortices that occupy the entire half-channel width. These produce an upwelling motion over the crests and a downwelling motion over the troughs, with no tertiary vorticity observed. The inner peak, of much lower intensity, is found at $\lambda _3^+ \approx 45$ and corresponds to small scale vortices extending by approximately 30 viscous units in the wall normal direction. The presence of two peaks mirrors the results of transient growth analysis in turbulent channels by del Álamo & Jiménez (2006) and Pujals et al. (2009) and suggests that surface topography modulation of the right spanwise length scale can excite a strong, steady response by leveraging amplification mechanisms intrinsic to the turbulent shear flow. However, a major difference with the optimal structures found by these works is that the strength of the steady response to surface modulations predicted by the present tool becomes asymptotically Reynolds number independent when the cross-plane velocities are scaled with the friction velocity. Fundamentally, this is due to the SA transport model utilised in this work, designed to produce the law of the wall and in which the turbulent eddy viscosity (and the Reynolds stresses driving secondary currents) become, asymptotically, Reynolds number independent.

For rectangular ridges, the present model suggests that (i)  both the ridge width $W$ and the gap between ridges $G$ are key parameters to quantify the response and that (ii) the analysis is more revealing when these two parameters are used and not other combinations previously used in the literature. More importantly, the largest response is found at a symmetric configuration where $W=G=0.67$, i.e. a rather large ridge width for a spanwise spacing of $S=G+W \approx 1.34$. For other ridge configurations, the secondary vortices compete for the available space with structures developing on the same ridge or over neighbouring ridges or are weakened by tertiary structures appear at large gaps or widths.

It is important to mention that the proposed approach has its limitations and it should not be seen as a replacement for DNS or experiments to obtain precise quantitative predictions. Secondary currents have been shown to display highly unsteady behaviour (Vanderwel et al. 2019) and meander in the spanwise direction (Hutchins & Marusic 2007; Kevin, Monty & Hutchins 2019; Zampiron et al. 2020). The present model assumes steady, streamwise-independent perturbations, and cannot fully capture any of these phenomena. Secondly, it is likely that surfaces characterised by prominent ridges, or surfaces with rapid lateral variations of the geometry (i.e. surfaces with sharp corners or with large lateral slope) cannot be satisfactorily modelled using linearised boundary conditions. For example, several authors (e.g. Hwang & Lee 2018) have observed that the wall-normal deflection of spanwise velocity fluctuations operated by the vertical sides of rectangular ridges is an important mechanism for the generation of secondary currents. This effect may be important, especially when the ridges protrude significantly in the log-layer and convective effects result in complex flow topologies, as in Castro et al. (2021). This mechanism is clearly not accounted for in the present model, where the ridge height is infinitesimal. Nevertheless, the model does produce secondary structures that resemble those observed in DNS or experiments. Most importantly, it correctly predicts the spanwise ridge spacing at which peak strength is achieved, in agreement with previous observations. This suggests that while different generation mechanisms might be at play, the linearised Navier–Stokes operator still provides an adequate description of the response to an external forcing.

With appropriate modelling assumptions, the present approach would also enable a rapid exploration of the vast parameter space characterising other surface heterogeneities that have been recently considered in the literature, for example, strip-type roughness (Willingham et al. 2014; Anderson et al. 2015; Chung et al. 2018), superhydrophobic surfaces (e.g. Türk et al. 2014; Stroh et al. 2016) or combinations of topography and roughness, as in, for example, Stroh et al. (2020b) and Schäfer et al. (2022). However, modelling the physics of the flow–surface interaction within a linear framework is less straightforward than modelling a modulation of the surface height, as in the present case. These configurations are currently being considered and results will be reported in future work.