Linear models of strip-type roughness

Abstract

Prandtl's secondary flows of the second kind generated by laterally varying roughness are studied using the linearised Reynolds-averaged Navier–Stokes approach proposed by Zampino et al. (J. Fluid Mech., vol. 944, 2022, p. A4). The momentum equations are coupled to the Spalart–Allmaras model while the roughness is captured by adapting established strategies for homogeneous roughness to heterogeneous surfaces. Linearisation of the governing equations yields a framework that enables a rapid exploration of the parameter space associated with heterogeneous surfaces, in the limiting case of small spanwise variations of the roughness properties. Channel flow is considered, with longitudinal high- and low-roughness strips arranged symmetrically. By varying the strip width, it is found that linear mechanisms play a dominant role in determining the size and intensity of secondary flows. In this setting, secondary flows may be interpreted as the time-averaged output response of the turbulent mean flow subjected to a steady forcing produced by the wall heterogeneity. In fact, the linear model predicts that secondary flows are most intense when the strip width is about 0.7 times the half-channel height, in excellent agreement with available data. Furthermore, a unified framework to analyse combinations of heterogeneous roughness properties and laterally varying topographies, common in applications, is discussed. Noting that the framework assumes small spanwise variations of the surface properties, two separate secondary-flow-inducing source mechanisms are identified, i.e. the lateral variation of the virtual origin from which the turbulent structure develops and the lateral variation of the streamwise velocity slip, capturing the acceleration/deceleration perceived by the bulk flow over troughs and crests of non-planar topographies.

1. Introduction

Prandtl's secondary flows of the second kind (Prandtl 1952) emerge when a turbulent flow develops over an heterogeneous surface with a lateral variation of its properties. Two equivalent standpoints, based on the analysis of the Reynolds-averaged equations, explain the formation of such currents. One standpoint considers secondary currents as the product of the imbalance between production and dissipation of turbulent kinetic energy induced by the roughness heterogeneity (Hinze 1973), whereby turbulence-rich fluid is advected towards low-turbulence regions. The second standpoint considers the streamwise vorticity balance (Perkins 1970), whereby cross-stream gradients of the Reynolds stresses arising from the cross-stream velocity components induce a turbulent torque that acts as a source term in the streamwise vorticity equation (Castro & Kim 2024). Overall, such mechanisms produce large-scale counter-rotating longitudinal rolls appearing in the time-averaged wall-bounded flow. The associated upwelling and downwelling motions produced by the rolls induce a lateral distortion of the boundary layer height (Barros & Christensen 2014), together with alternating high- and low-streamwise-momentum regions (Mejia-Alvarez & Christensen 2013; Willingham et al. 2014; Anderson et al. 2015), arranged analogously to the classical roll-streak pattern in shear flows (Brandt 2014).

Secondary flows are commonly observed in many industrial and environmental applications, where surfaces are either characterised by lateral variations of the topography, i.e. the elevation, or of the friction, e.g. by means of varying roughness properties. These two types of heterogeneity have been idealised in the literature as ridge-type and strip-type roughness configurations, respectively. The first type consists of longitudinal ribs located on a smooth, planar surface having rectangular or more complex cross-sections (Goldstein & Tuan 1998; Hwang & Lee 2018; Zampiron, Cameron & Nikora 2020; Castro et al. 2021; Long, Wang & Pan 2023; Zampino, Lasagna & Ganapathisubramani 2023; Zhdanov, Jelly & Busse 2024), or alternatively smooth sinusoidal modulations of the wall (Wang & Cheng 2006; Vidal et al. 2018). The second type, the focus of this work, consists of alternating longitudinal strips of high and low roughness. Secondary motions over such an arrangement have been extensively characterised experimentally (Bai et al. 2018; Wangsawijaya et al. 2020; Wangsawijaya & Hutchins 2022; Frohnapfel et al. 2024), and in numerical simulations (Willingham et al. 2014; Anderson et al. 2015; Chung, Monty & Hutchins 2018; Forooghi, Yang & Abkar 2020; Stroh et al. 2020; Neuhauser et al. 2022; Schäfer et al. 2022).

Despite the burgeoning interest in these flows and the intense examination of their characteristics, there is a number of aspects clearly documented in the literature for which a physics-based, mechanistic model is not available. The first aspect is related to the marked dependence of the size and intensity of secondary flows on one or more spanwise length scales characterising the surface heterogeneity (Medjnoun, Vanderwel & Ganapathisubramani 2018; Yang & Anderson 2018; Wangsawijaya et al. 2020). For strip-type roughness this length scale is usually expressed by the width $S$ of the strips. Consensus has emerged on the existence of three separate regimes as $S$ varies in relation with the average boundary layer thickness $\delta$ (Chung et al. 2018). When the strip width is much smaller than the boundary layer thickness, $S\ll \delta$, secondary flows are confined to the vicinity of the surface and do not strongly influence the outer region. Conversely, when $S\gg \delta$ secondary flows are localised in regions where the surface properties vary more rapidly, and wide areas of local flow homogeneity are observed away from such regions (Wangsawijaya et al. 2020). When $S \approx \delta$, the secondary flows are most intense and can significantly influence the flow structure. Nevertheless, a model that captures the nature of these regimes and identifies boundaries between them is not available at present. In addition, most studies have considered strips of equal width, but the width ratio between high- and low-roughness strips is certainly important, as it is for rectangular ridges where the ratio of recessed and elevated area influences the flow structure (Medjnoun, Vanderwel & Ganapathisubramani 2020; Zampino, Lasagna & Ganapathisubramani 2022).

A second aspect is related to the occurrence of the so-called tertiary flows. These are weaker longitudinal roll structures adjacent to the dominant rolls often associated to a reversal of the vertical flow direction at the centre of the high- (or low-)roughness strip, or at the centre of the ridge (or trough) (Vanderwel & Ganapathisubramani 2015). Tertiary flows are commonly observed over surfaces with longitudinal ridges (e.g. Medjnoun et al. 2020), especially when the width of the troughs or of the ridges is large enough to accommodate multiple streamwise vortices next to each other. For heterogeneous rough surfaces, however, tertiary flows have not been observed. In fact, for wide strips, cross-stream motions have been observed to be mostly confined in a roughly square region around the transition between the strips. This applies to both boundary layer (Wangsawijaya et al. 2020) and channel flows (Chung et al. 2018; Neuhauser et al. 2022). One explanation may be that tertiary flows over roughness strips might be difficult to discern in the mean flows obtained from experiments or simulations, especially when instantaneous structures meander quite significantly in the longitudinal direction (Kevin, Monty & Hutchins 2019; Zampiron, Cameron & Nikora 2021), smearing weak cross-stream motions. Neuhauser et al. (2022) hypothesised that the boundary conditions utilised in numerical simulations to capture the roughness effect may also play a role, although this hypothesis does not appear to explain why tertiary flows are not seen in experiments.

A third aspect that still lacks a robust mechanistic explanation is motivated by features of realistic surfaces in engineering and natural applications, whereby lateral changes of the roughness height are almost invariably accompanied by a lateral change in the elevation (Stroh et al. 2020; Schäfer et al. 2022). Decoupling these two effects may be easier in numerical simulations where the roughness heterogeneity is modelled by suitable spanwise heterogeneous boundary conditions applied to an otherwise planar boundary of the numerical domain (Chung et al. 2018; Neuhauser et al. 2022), but requires care when setting up experiments with, e.g. sandpaper strips or in roughness-resolving numerical simulations (Frohnapfel et al. 2024). One explicit attempt to study the coupling between these two effects was carried out by Stroh et al. (2020) and then later by Schäfer et al. (2022) who performed a series of direct numerical simulations over surfaces characterised by alternating rough and smooth regions. In their paper, Stroh et al. (2020) completely resolved the surface roughness using an immersed boundary method and studied three different configurations: the mean roughness height is (i) lower, (ii) equal to and (iii) higher than the elevation of the smooth surface. The authors observed a change in the flow organisation moving from case (i) to (iii) and vice versa. This behaviour was not reproduced in more recent roughness-resolving simulations (Frohnapfel et al. 2024), which was attributed to the importance of the strip width, relative to the roughness height.

One last aspect for which a model does not seem to be available concerns the relation between naturally occurring very-large-scale motions (VLSMs), populating the log-layer over homogeneous surfaces, and secondary flows (Chung et al. 2018; Lee, Sung & Adrian 2019; Wangsawijaya & Hutchins 2022). It has been speculated that secondary flows may be interpreted as VLSMs locked in place by the surface heterogeneity, given some similarity in their features. This can readily explain why secondary flows are most intense when $S \approx \delta$, because the strip width is commensurate with the spanwise length scale of such motions. Evidence shows that VLSMs and secondary flows do indeed coexist and do interact to a significant extent, since energy from the former appear to leak into the latter (Wangsawijaya et al. 2020; Zampiron et al. 2020). However, the specific mechanism for which large-scale structures residing in the outer layer should be locked in place so effectively by the roughness heterogeneity is not fully clear.

In recent work (Zampino et al. 2022), we developed a predictive framework to understand how far can linear mechanisms go in explaining these aspects, focusing on ridge-type roughness. The framework originates from the long line of work that relies on the Reynolds-averaged Navier–Stokes (RANS) equations, augmented with a turbulent viscosity model and linearised about the turbulent mean, to explain the structure of smooth-wall turbulence (see del Álamo & Jiménez 2006; Pujals et al. 2009; Hwang & Cossu 2010; McKeon & Sharma 2010 and references therein), or as a systematic tool to investigate flow control strategies (Moarref & Jovanović 2012; Luhar, Sharma & McKeon 2015) and patterned surfaces (Chavarin & Luhar 2020; Ran, Zare & Jovanović 2020). Differently to previous efforts (Meyers, Ganapathisubramani & Cal 2019), the framework utilises the Spalart–Allmaras (SA) equation (Spalart & Allmaras 1994) to capture turbulent viscosity transport phenomena in combination with the nonlinear quadratic constitutive relation (QCR; see Spalart 2000) to model the anisotropy of the Reynolds stress tensor, required to produce secondary motions (see Speziale 1982; Speziale, Sarkar & Gatski 1991; Bottaro, Soueid & Galletti 2006). It also assumes that spanwise variations of the surface topography are infinitesimally small. This allows the mean response of the turbulent flow to be obtained using linear equations where different spanwise wavenumber components have decoupled. The relevance of the assumption for surfaces with finite-amplitude topographies remains to be examined, although recent work on rectangular ridges (Castro & Kim 2024) suggests that this approximation may only be acceptable for ridges with moderate height. This, in turn, indicates that the mean response of the turbulent flow to a perturbation of the surface topography may be quite nonlinear. Nevertheless, one advantage of the linear framework is its computationally efficiency. It thus enables the vast parameter space characterising heterogeneous surfaces to be explored rapidly, for instance to unravel the effect of ridge geometry (Zampino et al. 2023). A second key advantage of the framework is that it provides a perspective of secondary motions as being the output response of the turbulent mean subjected to a steady perturbation produced by the surface heterogeneity.

In this paper, we bring the same framework to bear on the problem of strip-type roughness. We assume that the spanwise variation of the roughness is small, so that linear equations governing the response of the flow can be obtained. The effect of the surface roughness is introduced following well-established modelling strategies for rough walls (Aupoix 2007; Prakash & Laurendeau 2020). Briefly, such strategies consist of modifying the virtual origin from which the turbulent flow develops, in order to obtain the desired shift of the logarithmic velocity profile. The overall aim of this paper is to characterise the formation and structure of secondary flows developing above strip-type roughness by means of the proposed linear framework. This will allow fundamental insight into the linear mechanisms that control such flows to be generated. We apply the proposed linear framework to flows in channels, and examine the role that the surface arrangement plays on: (i) the strength of secondary motions as a function of the strip width, identifying the three regimes discussed in Chung et al. (2018); (ii) the occurrence of tertiary flows as the relative width of the high- and low-roughness strips is varied; (iii) the structure of low- and high-momentum pathways; and, finally, (iv) the combination of roughness and surface elevation effects.

The modelling framework, and its extension to rough surfaces, is presented in § 2. Results are then reported in § 3. In § 4, the framework is generalised to more-complex surface heterogeneities, combining the effects of roughness and surface elevation. Finally, conclusions are summarised in § 5.

2. Methodology

2.1. Governing equations

The incompressible flow of a fluid with kinematic viscosity $\nu$ and density $\rho$ is considered in a pressure-driven channel with half-height $h$ and subjected to a streamwise pressure gradient $\varPi$. The friction velocity $u_\tau =\sqrt {\tau _w/\rho }$, with $\tau _w = h \varPi$ the mean friction, yields the friction Reynolds number $Re_\tau =u_\tau h/\nu$. Index notation is used for the Cartesian coordinates $x_i$ and velocities components $u_i$. Quantities are generally normalised by $h$ and $u_\tau$. The superscript $({\cdot })^+$ is omitted in the following to reduce clutter, unless necessary to identify a length scaled by the viscous length. The channel walls are covered by alternating strips of high and low roughness having width $S_h$ and $S_l$, respectively, as shown in figure 1. The strips are aligned streamwise and are placed symmetrically on the two walls. The pattern repeats with spanwise periodicity $\varLambda = S_h + S_l$, the fundamental length scale. We also introduce the duty cycle $DC = S_h/\varLambda$ to characterise the relative width of the strips, and refer to $S$ as the strip width when $S_h = S_l$, i.e. for $DC=0.5$.

Figure 1. Sketch of the bottom half of the channel with high- and low-roughness strips, of width $S_h$ and $S_l$, respectively, placed symmetrically on the bottom and upper (not shown) walls and aligned with the streamwise direction $x_1$. This pattern repeats laterally with period $\varLambda = S_h+S_l$. These dimensions are normalised with the channel half-height $h$. The coordinate system is anchored at the bottom plane, at the centre of one of the high-roughness strips. Owing to the symmetry of the strip arrangement on the two walls, the shaded red area defines the fundamental repeating flow unit in which the flow structure is visualised later in the paper.

The continuity and momentum equations are Reynolds-averaged and made non-dimensional using $h$, $u_\tau$ and $\rho$. Average and fluctuation quantities are denoted by an overbar and a prime. For streamwise-aligned strips, we assume a streamwise-independent time-averaged flow, i.e. $\partial ({\cdot })/\partial x_1 \equiv 0$, which filters out the meandering of secondary currents (Zampiron et al. 2020). As a result, the mean pressure can be eliminated by considering the mean streamwise vorticity equation and introducing the streamfunction $\bar {\psi }$, satisfying $\nabla ^2 \bar {\psi }=\bar {\omega }_1$ with

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

the mean streamwise vorticity. 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}$. The Reynolds-averaged equations for the streamwise momentum and the streamfunction are then

(2.2a)\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.2b)\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( \displaystyle \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}

where $\tau _{ij}= - \overline {u_i'u_j'}$ is the Reynolds stress tensor.

2.2. Turbulence modelling

When the linear Boussinesq hypothesis is used to express the deviatoric component of the Reynolds stresses as a function of the mean velocity gradients, namely

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

where $\nu _t$ is the turbulent eddy viscosity and $S_{ij}$ is the symmetric component of the mean velocity gradient tensor

(2.4)\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}

the Reynolds stresses in (2.2b) do not depend on the streamwise velocity. Then, the streamfunction equation decouples from the streamwise momentum equation and its solution is trivially $\bar {\psi }\equiv 0$, i.e. no secondary flows are generated.

As discussed extensively in the literature (see e.g. Perkins 1970; Speziale 1982; Bottaro et al. 2006), a nonlinear stress model is needed to predict Prandlt's secondary flows of the second kind, produced by spatial gradients of the anisotropy of the Reynolds stresses. Several approaches have been proposed in literature (e.g. Speziale 1982; Chen, Lien & Leschziner 1997). Here we utilise the QCR nonlinear model presented in Spalart (2000), whereby the deviatoric component of the Reynolds stresses becomes

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

where $O_{ij}$ is the normalised rotation tensor defined as

(2.6)\begin{equation} O_{ij} = {2W_{ij}}\left/{\sqrt{{\frac{\partial \bar{u}_m}{\partial x_n} \frac{\partial \bar{u}_m}{\partial x_n}}}}\right., \end{equation}

and $W_{ij}$ is the antisymmetric part of the velocity gradient tensor, with $m$ and $n$ being summation indices. The QCR model depends on a tuning single constant, whose value $c_{r1}=0.3$ was calibrated to match the anisotropy of the outer region of wall-bounded turbulent flows in Spalart (2000). The default value is used throughout the paper.

To close the momentum equations, a model for the eddy viscosity $\nu _t$ is necessary. Previous studies that have utilised the linearised Navier–Stokes equations have adopted analytical eddy viscosity profiles to analyse smooth-wall turbulent flows (see del Álamo & Jiménez 2006; Pujals et al. 2009; Hwang & Cossu 2010; Morra et al. 2019 among others). Here, a complete transport model is preferred over such analytical ansätze, as it is not clear a priori how the eddy viscosity field should change when the mean flow structure is significantly distorted by secondary currents, or when roughness effects are important. For this purpose, the SA turbulence model (Spalart & Allmaras 1994) is employed in this work. The SA model is preferred here over other commonly employed two-equation models because it can be linearised relatively easily. In addition, the SA model was developed for attached shear flows, hence it should provide satisfactory predictions for the present case. The steady SA model defines a transport equation for the modified eddy viscosity $\tilde {\nu }$, normalised with $u_\tau$ and $h$. This quantity is related to the turbulent viscosity by the relation

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

where $f_{v1}={\chi ^3}/({\chi ^3+c_{v1}^3})$ 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. The term $f_{v1}$ ensures the correct decay of the turbulent viscosity in the viscous sublayer (Herring & Mellor 1968; Spalart & Allmaras 1994), although $\tilde {\nu }$ behaves linearly in the log layer down to the wall, which is advantageous for numerical reasons. The transport equation is

(2.8)\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}

where the terms model advection, production, diffusion and destruction, respectively. In the production term, the quantity $\tilde {{\mathcal {S}}}$ is defined as

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

with $\kappa$ the von Kármán constant. The destruction term in (2.8) 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.10)\begin{equation} f_{w} = g \left[ \frac{1+ c_{w3}^6}{ g^6 + c_{w3}^6}\right]^{1/6}, \end{equation}

with

(2.11a,b)\begin{equation} g = r+c_{w2} ( r^6 - r) \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$ and $c_{w3}=2$ are used (Spalart & Allmaras 1994), with $c_{w1}=c_{b1}/\kappa ^2+(1+c_{b2})/\sigma$ to balance production, diffusion and destruction in the log-layer, with $\kappa =0.41$.

2.3. Roughness model for homogeneous surfaces

Many rough-wall modelling strategies for homogeneous roughness for RANS simulations rely on the notion of equivalent sandgrain roughness $k_s^+$ (e.g. Durbin et al. 2000; Suga, Craft & Iacovides 2006; Aupoix 2007; Brereton & Yuan 2018; Prakash & Laurendeau 2020, among others). These strategies, described in this section, link the equivalent sandgrain roughness to suitable non-zero turbulence quantities (the modified eddy viscosity for the SA model) at the smooth, planar boundary of the numerical domain, to capture the increased turbulence activity near the rough surface and obtain the desired shift of the logarithmic velocity profile as the main effect of the surface roughness. No-slip boundary conditions are applied for the velocity. The mean turbulence structure is assumed to develop from a new virtual origin, displaced beneath the numerical boundary by a suitable distance $d_{0}^+$, to be determined (Rotta 1962). Relying on the outer-layer similarity hypothesis of Townsend (1976), far away from the surface the law of the wall is preserved, and the shift of the streamwise velocity profile observed over a rough surface is captured by the empirical relation

(2.12)\begin{equation} \left.\frac{\partial \bar{u}_{1,r} }{\partial x_2^+}\right|_{x_2^+} = \left.\frac{\partial \bar{u}_{1, s}}{\partial x_2^+}\right|_{x_2^+ + d_0^+}, \end{equation}

where the subscripts $({\cdot })_r$ and $({\cdot })_{s}$ denote quantities over the rough and smooth walls. Integrating this relation from the wall with $\bar {u}_{1,r}(x_2^+=0) = 0$ yields

(2.13)\begin{equation} \bar{u}_{1, r}(x_2^+) = \bar{u}_{1, s}(x_2^+ + d_0^+) - \bar{u}_{1, s}(d_0^+), \end{equation}

which evaluated far away from the surface gives the logarithmic shift

(2.14)\begin{equation} \Delta \bar{u}_1 = \bar{u}_{1, r}(x_2^+ \gg d_0^+) - \bar{u}_{1, s}(x_2^+ \gg d_0^+) = \bar{u}_{1, s}(d_0^+), \end{equation}

i.e. the distance $d_0^+$ can be found as the wall-normal coordinate where the velocity over the smooth wall is equal to the desired velocity shift $\Delta \bar {u}_1$.

Then a roughness function that links the equivalent sandgrain roughness to the shift of the log region is required. Among various options available, here we use the Colebrook–Grigson roughness function (Grigson 1992), given by

(2.15)\begin{equation} \Delta \bar{u}_1= \frac{1}{\kappa} \log \left( 1+ \frac{k_s^+}{ \exp(3.25 \kappa)} \right), \end{equation}

with $\kappa$ being the von Kármán constant. Although significant variations can be observed in the transitional regime, little practical difference are found for the fully rough regime in using this and other models, such as Nikuradze's roughness function (Aupoix 2007). Knowing the smooth-wall velocity profile $\bar {u}_{1, s}(x_2^+)$ and equating the relations (2.15) and (2.14) allows the displacement $d_0^+$ to be expressed as a function of the desired sandgrain roughness $k_s^+$.

With such information, a solution consistent with (2.12) can be found when the eddy viscosity satisfies

(2.16)\begin{equation} \nu_{t, r}(x_2^+) = \nu_{t, s}(x_2^+ + d_0^+). \end{equation}

This is achieved in two steps. First, the inhomogeneous boundary condition

(2.17)\begin{equation} \tilde{\nu}_r(x_2^+ = 0) = \tilde{\nu}_s(d_0^+) = d_0^+ \kappa/Re_\tau \end{equation}

is enforced to the modified eddy viscosity in the SA model, where the second equality stems from the fact that, in the SA model, $\tilde {\nu }$ varies linearly as $\kappa x_2$ near the wall by construction. Second, the distance $d$ between any point in the computational domain and the nearest wall appearing in the SA model, as a fundamental field variable that controls the balance between the production and destruction terms, needs to be updated to reflect the location of the new virtual origin, slightly below the numerical domain boundary. For instance, for the lower half of the channel, $d = x_2 + d_0$. Overall, this procedure allows the SA model to produce the desired shifted logarithmic velocity profile, consistent with the updated boundary condition (2.17) on the modified eddy viscosity.

Figure 2 illustrates example results at $Re_{\tau }=1000$, at which most of the results presented in later sections were obtained. Solutions were obtained numerically with an in-house RANS code, based on a Chebyshev-collocation discretisation method. Mesh independence studies, omitted here, showed that 252 collocation points where sufficient to obtain mesh-independent results. The nonlinear system of algebraic equations formed by the streamwise momentum equation and the SA equation was solved using a Jacobian-free Newton–Krylov technique (Knoll & Keyes 2004), using the ‘hookstep’ approach of Viswanath (2007) to improve convergence. Initial guesses for the streamwise velocity were obtained by first solving the momentum equation using Cess's analytical eddy viscosity profile (Reynolds & Hussain 1972). For a desired equivalent sandgrain roughness $k_s^+$ for the homogeneous surface, the Colebrook–Grigson roughness function in figure 2(a) is first used to obtain $\Delta \bar {u}_1$. Using the smooth-wall velocity profile obtained from the SA model (figure 2c), the virtual origin $d_0^+$ is obtained upon applying (2.14) (figure 2b). Repeating this procedure for several equivalent sandgrain roughness yields the curve ‘SA-CG’ in figure 2(b). Clearly, the choice of the smooth-wall velocity profile influences the results. For instance, coupling the Colebrook–Grigson formula to the log-law $\bar {u}_{1, s}(x_2^+) = \log (x_2^+)/\kappa + A$, with $\kappa = 0.41$ for consistency with the standard SA model and $A=5.1$, yields

(2.18)\begin{equation} d_0^+ = \left(1+\frac{k_s^+}{\exp(3.25 \kappa)}\right) \exp(-\kappa A) \simeq 0.0326 k_s^+, \end{equation}

denoted as ‘Log-CG’ in the figure. This virtual origin is then used for the boundary condition (2.17) and in the SA model. Overall, this yields shifted velocity profiles (figure 2c), that produce the desired $\Delta \bar {u}_1$ as a function of $k_s^+$. This is demonstrated in figure 2(a), which compares the Colebrook–Grigson roughness function (solid line) used at the first step with the logarithmic shift obtained at the last step of this procedure for three values of $k_s^+$, denoted by the circles. It is worth pointing out that small absolute variations of the turbulent viscosity distribution in figure 2(d) are sufficient to produce relatively significant alterations of the mean velocity profile. Numerically, this makes the equations relatively stiff to solve.

Figure 2. (a) The analytical Colebrook–Grigson roughness function (solid line) and the logarithmic shift obtained via the procedure described in the text (open circles). (b) Virtual origin $d_0^+$ as a function of the equivalent sandgrain roughness $k_s^+$ obtained by using the smooth-wall SA velocity profile or the standard log-law (black lines). The orange lines denote the derivative of $d_0^+$ with respect to $k_s^+$ (axes on the right-hand side). Numerical solutions of the streamwise momentum equation and the SA turbulence model for channels with smooth and homogeneous rough walls for $Re_\tau = 1000$. (c) Streamwise velocity. (d) Turbulent viscosity. The circles in panel (a) indicate the logarithmic shift associated to the profiles of panel (c), calculated using the procedure in the text. For the logarithmic law, the constants $\kappa = 0.41$ and $A = 5.1$ are used.

2.4. Roughness model for heterogeneous surfaces

The equivalent sandgrain roughness is a dynamic parameter that is non-trivially related to the roughness geometry. For homogeneous roughness, it can be readily estimated from correlations once the shift of the velocity profile is known. However, for heterogeneous roughness, e.g. the present surface with alternating strips, it is not immediately clear how one should assign an equivalent sandgrain roughness to the two strips from velocity measurements, as the flow structure and, thus, the logarithmic shift also depend on the spatial distribution of the roughness properties (Wangsawijaya et al. 2020). This conundrum is fundamentally the same as discussed in numerical simulation studies in which the roughness is not resolved but suitable boundary conditions are applied at the smooth, planar boundary of the numerical domain. In such strategies, the roughness heterogeneity can be modelled directly by a lateral variation of the shear stress (Chung et al. 2018) or by a lateral variation of the transversal slip length (Neuhauser et al. 2022). These strategies require a model that links the boundary conditions to the desired logarithmic shift. Such models are generally derived for homogeneous surfaces and their applicability to heterogeneous surfaces may be questioned.

Here, given the lack of a better strategy, the aforementioned approach is adopted. Specifically, the alternating strips are defined by a spanwise variation of the equivalent sandgrain roughness, following the expression

(2.19)\begin{equation} k_s^+(x_3) = k_s^{(0)}+\epsilon k_s^{(1)}(x_3), \end{equation}

where $k_s^{(0)}$ is the reference, spatially constant roughness height and $k_s^{(1)}(x_3)$ captures the variation of the roughness properties over the two strips. For a unitary $k_s^{(1)}$ amplitude of the roughness pattern, the spanwise variation is defined by the unitary peak-to-peak amplitude, zero-mean function

(2.20)\begin{equation} k_s^{(1)}(x_3) = \begin{cases} 1 - DC & 0 \leq x_3 \leq S_h/2 \ \mathrm{and} \ \varLambda - S_h/2 \leq x_3 \leq \varLambda\\ - DC & S_h/2 \leq x_3 \leq \varLambda - S_h/2, \end{cases} \end{equation}

as demonstrated in figure 3. Because the difference in roughness between the two strips defined by $k_s^{(1)}(x_3)$ is unitary, the parameter $\epsilon$ in (2.19) controls the actual difference in roughness between the two strips, although it is not related to the physical structure of the roughness. This definition is preferred over specifying the roughness of the two strips, or considering roughness strips separated by smooth regions (Wangsawijaya et al. 2020). This choice is motivated by the fact that we consider in the present analysis the asymptotic limit when $\epsilon$ tends to zero, so that linearised equations governing the response of the turbulent shear flow developing over a rough surface to a small spanwise variation of the roughness properties can be obtained.

Figure 3. The hashed diagram on the bottom shows the heterogeneous rough surface with roughness strips having higher or lower roughness than average, covering about two periods of the laterally repeating pattern. The duty cycle is $0.8$. The diagram at the top shows the distribution of the perturbation roughness height $k_s^{(1)}(x_3)$. Note that this term represents the deviation from the reference roughness height $k_s^{(0)}(x_3)$, as defined by (2.19).

2.5. Linearisation of the Reynolds-averaged equations and of the roughness model

The streamwise momentum, streamfunction and SA equations form a coupled system of three nonlinear partial differential equations that can be solved for any desired strip configuration. However, when the difference between the properties of the high- and low-roughness strips is small, in the limit when $\epsilon \ll 1$, the resulting flow structure can be thought of as the sum of the flow in a channel with the homogeneous reference roughness $k_s^{(0)}$ and a small perturbation, produced by the surface roughness heterogeneity $k_s^{(1)}(x_3)$ and capturing the heterogeneous flow structure of secondary flows. This small perturbation obeys a set of linear equations which is much easier to solve, and only captures linear input–output mechanisms. To derive such equations, a generic time-averaged flow quantity $\bar {q}(x_2, x_3)$ is first expanded in series as

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

Higher-order terms in (2.21) are neglected within the current framework. The convergence of the series for finite $\epsilon$ and the validity of the resulting predictions must still be verified. However, this approach is motivated by the goal of calculating the linear response of the turbulent mean flow to a small, nearly infinitesimal, perturbation in surface attributes, to assess how well linear mechanisms can account for the formation of secondary structures over heterogeneous surfaces.

Substituting this expression for all mean quantities in the Reynolds-averaged equations and in the SA equation, and taking terms at order zero in $\epsilon$, leads to the nonlinear equations governing the flow over the homogeneous rough surface. The streamwise vorticity equation is trivially satisfied by $\psi ^{(0)} = 0$. The streamwise momentum equation is

(2.22)\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}

and it is coupled to the SA transport equation via the definition of the Reynolds stress $\tau ^{(0)}_{12}$. These two equations are solved in a coupled fashion using the approach discussed in § 2.3. At first order, the equations governing the perturbation of the streamwise velocity and the streamfunction are

(2.23a)$$\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.23b)$$\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} \left( \tau_{33}^{(1)}-\tau_{22}^{(1)}\right)+\left( \frac{\partial^2}{\partial x_2^2}-\frac{\partial^2}{\partial x_3^2}\right)\! \tau_{23}^{(1)}, \end{gather}$$

where we define $\varGamma =\partial u_1^{(0)}/\partial x_2$, showing that the zero-order solution, through the mean velocity gradient $\varGamma$, needs to be available for the solution of the first-order equations. The first term on the left-hand side of (2.23a), analogous to the off-diagonal coupling operator in the Orr–Sommerfeld–Squire linearised equations, is the only coupling term explicitly appearing in this set of equations. Physically, this term captures the interaction between the mean shear and the perturbation velocity and underpins energy extraction mechanisms in shear flows via the lift-up effect (Brandt 2014). All other terms obtained from the nonlinearity vanish because the streamfunction at order zero is identically zero. Secondary currents also introduce an alteration of the spatial organisation of the turbulent viscosity through an alteration of the balance of the transport terms in the SA equation. The linearised SA equation governing such organisation is coupled to the streamwise momentum and the streamfunction equations and contributes to the perturbation of the Reynolds stresses entering (2.23). Linearisation of the SA model is tedious and leads to complex expressions. More detail on the linearisation procedure is reported in our previous work (see Appendix B of Zampino et al. (2022)), and is omitted here for brevity.

It is worth noting that the streamfunction equation contains the perturbation of the Reynolds stresses originating from the cross-stream velocity components, as is well known (Perkins 1970). Although at order zero these terms exhibit negligible influence, at order one the perturbation of the Reynolds stress tensor becomes pivotal to couple the two equations in the differential system (2.23). Here, the first-order stresses are found by expanding the nonlinear Reynolds stress model (2.5) in a Taylor series in $\epsilon$, leading to

(2.24)\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.24), the individual perturbation Reynolds stresses appearing in (2.23) are

(2.25a)$$\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.25b)$$\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.25c)$$\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.25d)$$\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.25e)$$\begin{gather}\tau_{33}^{(1)}=2 \nu_t^{(0)} \frac{\partial^2 \psi^{(1)}}{\partial x_2 \partial x_3}. \end{gather}$$

Except for $\tau _{33}^{(1)}$, which coincides with its linear Boussinesq's definition, all other stresses contain an additional term specific to the QCR model and proportional to the $c_{r1}$ constant. In particular, the stresses appearing in the streamfunction equation contain spatial gradients of the streamwise velocity, and vice versa. These terms result in a tighter, two-way coupling between the streamfunction and streamwise velocity equations, now able to sustain secondary currents.

To obtain boundary conditions for the field variables, the wall roughness treatment model discussed in § 2.4 needs to be linearised. The key idea is that small spanwise perturbations of the equivalent sandgrain roughness are modelled as small spanwise variations of the virtual origin. More formally, over the heterogeneous surface given by (2.19), the shift of the virtual origin varies according to

(2.26)\begin{equation} d_{0}^+(x_3) = d_{0}^{(0)} + \epsilon d_{0}^{(1)}(x_3), \end{equation}

where $d_0^{(0)}$ is the shift of the virtual origin of the reference homogeneous surface with equivalent sandgrain roughness $k_s^{(0)}$. On the other hand, the first-order term can be found by differentiating numerically the curve reported in figure 2(b) at $k_s^+= k_s^{(0)}$ and using (2.20), since

(2.27)\begin{equation} d_0^{(1)} = \left. \frac{\mathrm{d} d_0^+}{\mathrm{d} k^+_s}\right|_{k_s^{(0)}} k_s^{(1)}. \end{equation}

Asymptotically, for large $k_s^{(0)}$ and considering the log law, the derivative of the curve in figure 2(b) tends to about 0.0326, implying that a peak-to-peak variation of the equivalent sandgrain roughness of $1/0.0326 \approx 30$ is necessary to obtain a peak-to-peak variation of the virtual origin equal to the viscous length scale. Once $d_{0}^{(1)}$ is known, linearising (2.17) yields the wall condition for the modified eddy viscosity

(2.28)\begin{equation} \tilde{\nu}^{(1)}(x_2 = 0) = d_0^{(1)} \kappa/Re_\tau, \end{equation}

showing that the spanwise variation of the roughness properties is modelled as a lateral change in the eddy viscosity at boundary of the numerical domain. Finally, homogeneous boundary conditions are used for the streamwise velocity perturbation and the streamfunction perturbation and its wall-normal derivative.

One last remark is in order. A sensible question is whether the inclusion of the linearised turbulence model to describe the perturbation of the turbulent viscosity is really necessary, given that this is not customary in many previous studies using linearised Navier–Stokes equations. On the one hand, the strength of the mean flow response to a lateral perturbation of the surface attributes may likely depend on the well-known selective amplification properties of the linearised Navier–Stokes operator, which are largest when such perturbation occurs at a specific spanwise length scale. On the other hand, the SA model provides a means to model realistically the effect of the surface heterogeneity, because it provides clear insight into how the perturbation of the effective distance $d$ influences the perturbation of the turbulent viscosity field. As described in § 4, the lateral perturbation of the distance $d$ is the dominant source mechanism that leads to secondary flows in the present framework. In principle, one could first introduce an ansatz on the perturbation of the turbulent viscosity and then only solve the linearised Navier–Stokes equations. However, given the range of transport phenomena modelled by the SA equation, defining the correct ansatz does not appear to be a straightforward task.

2.6. Numerical solution of the linearised equation

The spanwise variation of the equivalent roughness height can be modelled as a square wave approximated by the cosine series

(2.29)\begin{equation} k_s^{(1)}= \sum_{n=1}^{\infty} k_s^n \cos \left (n \frac{2{\rm \pi} }{\varLambda} x_3\right). \end{equation}

The coefficients $k_s^n$ can be calculated analytically for each combination of widths $S_h$ and $S_l$ and the corresponding coefficients $d_0^n$ for the spanwise variation of the virtual origin $d_0^{(1)}$ are found by using (2.27).

Expanding the unknown field variables at first order in series, e.g. for the streamwise velocity

(2.30)\begin{equation} u_1^{(1)} (x_2, x_3) = \sum_{n=1}^{\infty} \hat{u}_1(x_2; n) \cos \left (n \frac{2{\rm \pi} }{\varLambda} x_3 \right), \end{equation}

and substituting these expressions in the linearised equations leads to one set of three linear ordinary differential equations in $x_2$ for each integer wavenumber. As opposed to previous studies considering the linearised Navier–Stokes equations (Chavarin & Luhar 2020; Ran et al. 2020), each set of three ordinary different equations is independent of all other wavenumbers and can be solved in isolation. This would not be the case if higher-order terms had been retained in (2.21), and a larger problem would need to be solved taking into account harmonic interactions. A Chebyshev-collocation method was used for the discretisation. Although the field variable $d$ in the SA model has a sharp cusp at $x_2=0$ and, hence, a spectral technique is not ideal for the solution of this problem, we have observed that the numerical method is robust enough to provide accurate results when a sufficiently fine grid is used. In the following simulations, we used up to 252 collocation points. For the spanwise discretisation, we observed that solutions converge relatively rapidly with the number of Fourier modes retained in the expansion (2.29). This can be motivated by the observation that, far away from the wall, only large-scale perturbations of the surface features can influence the flow structure, while the effect of small-scale perturbations, i.e. sharp gradients of the boundary conditions (Neuhauser et al. 2022), decays more rapidly with the distance from the wall (Meyers et al. 2019). The number of spanwise modes required increases with the fundamental wavelength $\varLambda$. We always checked that results did not change visibly when doubling the number of modes. As a reference, 20 modes were sufficient at $\varLambda \approx 1$ to obtain a converged description of the perturbation velocity field.

The final solution is then found by combining the solutions at each wavenumber, as the superposition principle applies. One important implication of this property is that the flow structure over surfaces with complex topographic/roughness characteristics (Mejia-Alvarez & Christensen 2013; Barros & Christensen 2014) may be rationalised and better understood by decomposing the surface forcing into its constitutive components. The strength of the mean flow response at each spanwise length scale will then depend on the amplitude of such components times a factor that captures the selective amplification properties of the linearised Navier–Stokes operator (Chernyshenko & Baig 2005), augmented by the linearised SA equation.

It is worth noting that the spatially constant component at $n=0$ does not appear in (2.29) because of the assumption that the spanwise variation of the roughness height given by (2.20) is zero mean. It does also not appear in the solution (2.30), because the linearity of the model implies that a perturbation of the surface properties at wavenumber $n$ only produces a distortion of the time-averaged flow at the same wavenumber. A corollary of this property is that the present model does not predict any change in mean friction drag, the subject of several recent studies (Hutchins et al. 2023; Frohnapfel et al. 2024). In fact, the spanwise-constant component $\hat {u}_1(x_2; 0)$ and, thus, the perturbation of the bulk velocity computed from this profile, which would allow calculating the change in friction coefficient at constant friction velocity, is identically zero. The model does indeed capture the spanwise modulation of the streamwise velocity distribution, i.e. high- and low-momentum pathways (Barros & Christensen 2014), but second-order effects in $\epsilon$ that produce interactions between harmonics are necessary to obtain a velocity perturbation at wavenumber $n=0$ from surface perturbations at $n>1$ and, thus, capture the change in friction (Zampino 2023).

3. Structure and strength of secondary currents

The volume-averaged kinetic energy of the cross-sectional velocity components ${\mathcal {K}}$, defined as

(3.1)\begin{equation} {\mathcal{K}} = \frac{1}{4\varLambda} \int_{0}^{2} \int_{0}^{\varLambda} [ u_2^{(1)}(x_2, x_3)^2 + u_3^{(1)}(x_2, x_3)^2] \,\mathrm{d}\kern0.07em x_3 \,\mathrm{d}\kern0.07em x_2, \end{equation}

is used here to characterise the strength of the secondary flows. We also use the streamfunction peak $\max _{x_2,x_3} |\psi ^{(1)}(x_2,x_3)|$ to quantify the cross-stream flow rate associated with the vortices, as in other studies (Vidal et al. 2018). Note that these variables are scaled with $u_\tau$ and $h$. The solution of the linearised equations for a given strip configuration can be obtained quite rapidly, which enables a rapid exploration of the parameter space $(S_h, S_l)$. Results are reported in figure 4(a,b) for $Re_{\tau }=1000$, and using $k_s^{(0)} = 180$. Figure 4(c,d) shows cuts along lines for three duty cycles as a function of the fundamental length scale $\varLambda$.

Figure 4. Contours of the volume-averaged kinetic energy of the cross-stream velocities ${\mathcal {K}}$ (a) and the streamfunction peak value $\max _{x_2,x_3}|\psi ^{(1)}|$ (b) as a function of the width of the high- and low-roughness strips. The Reynolds number is $Re_\tau =1000$ and $k_s^{(0)} = 180$. Cases at constant spatial fundamental length scale $\varLambda$ are identified by the dashed diagonal lines with negative slope. Markers identify cases discussed later in the text. (c,d) The same two quantities for $DC=0.2$, 0.3 and 0.5, as a function of $\varLambda$.

We note that, as in our previous application of these techniques to surfaces with longitudinal ridges (Zampino et al. 2022), the results of the linearised model become asymptotically Reynolds number independent for high Reynolds numbers, somewhat supporting the weak Reynolds number dependence documented in the literature (Wangsawijaya & Hutchins 2022). This ultimately stems from known properties of the SA model (Spalart & Allmaras 1994), which is designed to produce an eddy viscosity distribution consistent with the log law. Hence, the discussion presented here can also be applied to higher-Reynolds number-flows relevant to applications. Note also that flow variables, such as the velocity or streamfunction perturbations, are computed in the present linear modelling framework per unit variation of the equivalent sandgrain roughness height $k_s^{(1)}$ (scaled in inner units) in analogy to what was described for ridge-type roughness in Zampino et al. (2022) where the same quantities are obtained per unitary ridge height (scaled in outer units). Given that experiments on secondary flows over heterogeneous surfaces are often conducted on roughness strips with a considerable difference in roughness properties (Chung et al. 2018; Wangsawijaya et al. 2020), the numerical values reported here will appear quite small. For graphical convenience, quantities are pre-multiplied by a large factor, e.g. $10^8$ in figure 4(a), in the figures.

Regardless of the quantity used for measuring the strength of the secondary currents, a peak is observed for $S_h = S_l \simeq 0.7$, corresponding to a fundamental length scale $\varLambda \approx 1.4$, although the streamfunction peaks slightly later. Further, the quantities in figure 4 are symmetric with respect to the line $DC=0.5$, where the strength peaks. For $\varLambda \gtrsim 2.5$ two peaks are observed, located symmetrically with respect to the line $DC=0.5$. Examination of the flow structure for some of these cases indicates that the secondary flows observed over the high- and low-roughness regions for a generic configuration $(S_h, S_l)$ are identical in strength but opposite in flow direction when the width of the two strips is swapped. The strip width at which secondary currents are most intense reflects previous observations. For instance, Chung et al. (2018), using large-eddy simulations, and Wangsawijaya et al. (2020), using experiments above spanwise-alternating smooth and rough strips, reported a maximum intensity for the secondary flows when the width of the strips is comparable with the boundary layer thickness. In particular, Wangsawijaya et al. (2020) found that the swirl strength was largest among the cases they considered for $S_h = S_l = 0.62$.

The contours of the kinetic energy around the peak region appear elongated along the line $S_h+S_l={\rm const}$. As a result, the three cuts along lines at constant duty cycle all display a peak for $\varLambda \simeq 1.4$, which may interpreted as $\varLambda$ being the relevant length scale. This is partly correct, as the response to sinusoidal perturbations of the sandgrain roughness does indeed peak for this length scale, as for ridge-type roughness (Zampino et al. 2022). However, there is a marked effect of the relative size of the strips away from the peak and the two widths are indeed necessary to correctly characterise the response. The maps of figure 4 show strong similarities with the maps displayed in Zampino et al. (2022, see figure 10) and Zampino et al. (2023) (see figure 2) for the ridge-type roughness, as the peak amplification occurs for similar values of $\varLambda$. This similarity suggests that the selective amplification of secondary flows is an intrinsic property of the mean flow, and perhaps less strongly an effect of the type of forcing, e.g. whether it is produced by elevation (ridges) or roughness (strips) variations. In this regard, there has been recent discussion on the relation and coexistence between secondary currents and VLSMs (e.g. Lee et al. 2019). One speculation is that secondary currents are naturally occurring VLSMs that are phase locked spatially by the heterogeneous surface (Chung et al. 2018; Wangsawijaya et al. 2020) and emerge in the time-averaged flow. The present approach, based on the Reynolds-averaged equations where the concept of VLSMs does not apply immediately, suggests that secondary currents may be interpreted as the time-averaged response of a forcing localised near the wall and produced by gradients of the turbulent stresses. The linearised Navier–Stokes operator with its selective amplification properties then produces more intense time-averaged structures at specific forcings wavelengths, when $\varLambda \simeq 1.4$. This hypothesis leverages the same physics used in transient growth analysis studies (del Álamo & Jiménez 2006; Cossu, Pujals & Depardon 2009; Pujals et al. 2009) to explain the formation of coherent structures in shear flows from properties of the Orr–Sommerfeld–Squire equations, except that here we consider the steady response to a given steady perturbation localised at the wall rather than the transient amplification of an optimal initial perturbation.

3.1. Flow structure

To visualise the cross-stream structure of the secondary currents, fields of the wall-normal and spanwise velocity perturbation are reported in figure 5, for duty cycle $DC=0.5$, i.e. for equal length of the high- and low-roughness strips, for a range of strip widths $S$. These configurations correspond to the star markers in figure 4(a). Contour lines of the streamfunction are also reported. The roughness strips produce two counter-rotating vortices inducing a downwelling over the high-roughness regions and an upwelling over the low-roughness regions. For narrow strips, even narrower than what shown here, the vortices are confined in the near-wall region and the flow appears homogeneous at distances from the wall larger than the strip width. Increasing the strip width, the vortices grow with $S$ until they occupy the full half-height of the channel. For $S \simeq 0.7$, the cross-stream components and, in particular, the wall-normal component is more intense than all other cases considered. This configuration identifies the peak amplification region of the maps of figure 4. Increasing $S$ further, the strength of the wall-normal motions decreases slightly, but the volume-averaged strength of such motions decreases much further as the flow structure converges to an idealised wide-strip asymptotic limit where moderately intense cross-stream motions are only found in the immediate vicinity of the transition between strips, with fluid is at rest in the ‘homogeneous’ regions above the centres of the roughness strips. This explains the trend of the two quantities in figure 4(c,d). While ${\mathcal {K}}$ is a volume-averaged quantity and decreases with $S$ for $S\gg 1$, the streamfunction peak is a local quantity that measures the strength of the individual vortex cores. This quantity shows a first peak for $S \simeq 0.75$, and eventually tends to an asymptotic value (with a larger amplitude) for large $S$, characterising the strength of the ‘isolated vortex’ regime.

Figure 5. Colour maps of the wall-normal velocity perturbation $u_2^{(1)}$ over roughness strips at $DC=0.5$ for several strip widths $S$ (a). The field in the fundamental domain (see figure 1) is repeated four and two times for the two narrowest strip cases. Contour lines of the perturbation of the streamfunction $\psi ^{(1)}$ are also reported to better describe the secondary flows. Negative $\psi ^{(1)}$ are indicated by using dashed lines. Contours of the spanwise velocity perturbation $u_3^{(1)}$ (b). The Reynolds number is $Re_{\tau }=1000$, for $k_s^{(0)} = 180$. The darker/lighter patches denote the high/low roughness strips.

Near the transition between strips, the spanwise component is particularly intense in the near-wall region (figure 5b). The spanwise velocity peak is more intense than the vertical velocity peak in agreement with observations Frohnapfel et al. (2024), and is localised in the near-wall region. Given that the spanwise velocity obeys no-slip condition, this results in very high streamwise vorticity localised at the transition between strips. Analyses not reported here show that the wall-normal location of the spanwise velocity peak scales in inner units when the Reynolds number is increased, while its magnitude becomes $Re_\tau$ independent.

From a qualitative viewpoint, the flow structure predicted by the present linearised model resembles previous experimental observations (Wangsawijaya et al. 2020) and numerical simulations (Chung et al. 2018; Neuhauser et al. 2022). One aspect of discussion concerns the spanwise location of the streamwise-aligned vortices with respect to the alternating pattern of roughness. In the present case, the vortices are symmetrically located above the interface between the strips. This is because the linearity of the governing equations preserves the symmetry that exists across the jump. By contrast, the experiments of Wangsawijaya et al. (2020) show that the centres of the vortices are typically found over the low-roughness region. The same was predicted by the simulations of Chung et al. (2018), using an inhomogeneous shear-stress boundary condition to model the roughness. The puzzling aspect is that one would initially attribute the displacement of the vortices to nonlinear convective effects not captured by the linear model, as if the vortices were transported towards the low-roughness region by the relatively intense spanwise velocities near the wall associated to the streamwise vorticity field. However, such a displacement is not observed in the simulations of Neuhauser et al. (2022) who applied a Navier slip boundary condition for the spanwise velocity to model the roughness, or occurs in the opposite direction in more-recent roughness-resolving simulations of submerged roughness strips (Frohnapfel et al. 2024).

From a quantitative viewpoint, a comparison with published results is slightly less straightforward given the particular set-up considered in this paper. For this purpose, we use the channel-flow simulations of Chung et al. (2018), at $Re_\tau = 590$. The roughness strips are modelled by setting the shear stress to 50 % more and 50 % less than the average shear stress. Assuming that the low- and high-roughness strips correspond to the smooth and rough wall regions, respectively, Chung et al. (2018) estimates that the equivalent sandgrain roughness of the roughness patches is $k_s^+ = 205$. With such settings, they observe maximum wall-normal velocities that peak between $0.3u_\tau$ and $0.4 u_\tau$ (see their figure 9a). To match the roughness properties of these simulations, one needs to recall that the solution produced by the present linear model is defined per unit variation of the equivalent sandgrain roughness between the strips. From the results of figure 5, maximum velocities on the order of $17 \times 10^{-4} u_\tau$ are obtained, for the optimal width $S$. Multiplying this value by $k_s^+ = 205$ we obtain velocities on the order of $0.35 u_\tau$, in very good quantitative agreement with the numerical simulations. Given that the intensity of the cross-stream velocity components characterises somehow the equilibrium between source and sink mechanisms of the streamwise vorticity balance (Stroh et al. 2016; Castro & Kim 2024), the favourable agreement with simulations suggests that such mechanisms are correctly captured by the linearised RANS model. However, the maximum wall-normal velocity in Chung et al. (2018) is obtained for a relatively wide strip, $S = 1.57$, while the present model indicates that the peak occurs at $S\approx 0.7$, and lower velocities are observed for $S=1.57$. It is argued that this is not a Reynolds number effect, but it is due to the vortices in Chung et al. (2018) being, as discussed, closer to each other than the strip width $S$ would suggest, resulting in larger induced velocities. Evidence for this is given by the fact that the wall-normal velocities depend on the duty cycle $DC$, as shown later in figure 9, and peak when the vortices are artificially pushed together by a narrow low-roughness strip, $S_L < S_h$.

3.2. High- and low-momentum pathways

A unique characteristic of flows over heterogeneous surfaces is that the longitudinal secondary currents are flanked by high- and low-momentum regions, produced by the vertical ‘pumping’ of high- and low-momentum fluid, respectively, induced by the vortical motions (Mejia-Alvarez & Christensen 2013; Barros & Christensen 2014; Willingham et al. 2014), which determine a significant spanwise alteration of the boundary layer depth. Although a wall-bounded flow above an homogeneous rough surface displays a positive deficit of the streamwise velocity due to the flow deceleration induced by the surface roughness (Jiménez 2004), the colocation between the high- and low-momentum regions and the roughness strips is counterintuitive as faster flow can be found for certain conditions above the high-roughness strips. In this section, we demonstrate that the present framework clearly captures this phenomenon.

In figure 6, contours of the streamwise velocity perturbation $u_1^{(1)}$ are shown for the same configurations of figure 5. This visualisation differs to what customarily reported in previous work in that it shows the velocity deviation from the streamwise velocity distribution $u_1^{(0)}$ observed over the homogeneous surface. The difference between the maps in the two columns is that in figure 6(a) the secondary currents were artificially ‘turned off’ by setting the constant $c_{r1}$ in the nonlinear stress model to zero, so that spanwise and wall-normal gradients of the mean streamwise velocity do not produce any of the Reynolds stresses in the streamwise vorticity equation necessary to sustain longitudinal vortices. In figure 6(b), solutions for the standard value $c_{r1} = 0.3$ are reported. Without secondary currents, the flow experiences a net deceleration above the high-roughness strips, especially in the near-wall region. However, further away from the wall, e.g. at the centre of the channel, the change in streamwise velocity depends strongly on the strip width, because this parameter controls the depth at which the roughness-induced deceleration ‘diffuses’ in the shear flow from the wall due to the turbulent viscosity field. For narrow strips, the velocity deficit at the channel centre is small, and only becomes significant when the strip width is at least half the half-height of the channel. When secondary currents are ‘turned back on’ (figure 6b), the streamwise velocity over the high-roughness strips is now generally positive because of the downwelling motion in this region, and vice versa for the low-roughness strips. This only applies for $x_2 \gtrsim 0.1$, because nearer to the wall the local roughness properties control whether the flow is faster/slower. These motions produce dispersive stresses, e.g. $u_1^{(0)}u_2^{(0)}$, that alter the equilibrium in the streamwise direction and result in a non-trivial dependence of the streamwise velocity from the wall. This, fundamentally, implies that the logarithmic velocity distribution is significantly altered by the addition of the dispersive stresses. However, as the strips become wider, secondary currents are not intense enough to produce any significant alteration of the streamwise momentum equilibrium and the deceleration effect produced by the high roughness begins to dominate, starting from the region closest to the wall. Overall, this is the same behaviour observed experimentally (Wangsawijaya et al. 2020; Frohnapfel et al. 2024), where a downwards/upwards bulging of the contours of the streamwise velocity are observed at the edge of the boundary layer/near the wall.

Figure 6. Contours of the streamwise velocity perturbation $u_1^{(1)}$ over roughness strips at $DC=0.5$ for several strip widths $S$, without QCR model (a) and with QCR model (b). The field in the fundamental domain (see figure 1) is repeated four and two times for the two narrowest strip cases. Contour lines of the perturbation of the streamfunction $\psi ^{(1)}$ are also reported to visualise the secondary flows. Negative $\psi ^{(1)}$ contours are indicated by dashed lines. The Reynolds number is $Re_{\tau }=1000$, for $k_s^{(0)} = 180$. The darker/lighter patches denote the location of the high-/low-roughness strips.

The influence of the strip width on the perturbation of the streamwise velocity field is summarised in figure 7. The streamwise velocity profile at the centre of the high-roughness strip is extracted from several calculations with $S$ in the range $[0.1, 20]$, with and without the nonlinear Reynolds stress model. These profiles are concatenated together to form the maps in figures 7(a) and 7(b), respectively. For the case with $c_{r1}=0.3$, we also report a similar plot for the wall-normal component, in figure 7(c). Given that a logarithmic velocity shift is not a meaningful quantity to compute, we report in figure 7(d) the velocity deviation in the centre of the channel, above the centre of the high-roughness strips. In figure 7(c), the solid red line denotes the wall-normal location of the perturbation streamfunction peak. The three regimes discussed in Chung et al. (2018) can be clearly identified. For wide strips, the velocity deficit tends to a value controlled by the roughness function at $k_s^{(0)}$ for both cases, as the turbulent structure over each strip tends to its equivalent over an homogeneous surface given that the influence of neighbouring patches and of the secondary currents localised at the transition between strips vanishes. In this regime, the model predicts that the centre of the rolls is located at a distance from the wall of about $0.42$, i.e. the rolls are space-filling in the vertical direction. The streamwise velocity perturbation is also roughly constant as a function of the distance from the wall for $S \gtrsim 8$, producing the expected shift of the logarithmic velocity profile over each strip, regardless of whether the nonlinear Reynolds stress is active or not.

Figure 7. Streamwise velocity extracted on a vertical line at the centre of the high-roughness strip for increasing strip width $S$, for $c_{r1}=0$ (a) and the standard value 0.3 (b). Wall-normal velocity extracted on the same vertical line for $c_{r1}=0.3$ (c). The dashed lines indicate wall-normal distances growing linearly with $S$ for $DC=0.5$. Streamwise velocity perturbation above the centre of the high-roughness strip, at the mid-plane, for the same two cases as a function of the strip width (d). The Reynolds number is $Re_{\tau }=1000$ and $k_s^{(0)} = 180$.

For narrow strips, the height of the channel is much larger than the ‘depth’ at which the effect of the wall inhomogeneity is perceived, analogously to the blending height concept for spatially varying roughness discussed in (Bou-Zeid, Parlange & Meneveau 2007). For the case $c_{r1}=0$, without rolls, this depth (estimated from the contours of the streamwise velocity perturbation) varies linearly with $S$, as one would expect to see in a diffusion-driven problem, given that the turbulent viscosity does not vary significantly except for near the wall. For the case $c_{r1}=0.3$, with rolls, the depth measured in terms of the wall-normal velocity also appears to increase linearly with $S$, at least for the smallest $S$ analysed here. However, the linear scaling (the dashed lines) appears to lose accuracy relatively rapidly, at $S\approx 0.25$, as soon as the rolls occupy about half of the half-channel height. In between these two regimes, in the ‘transitional regime’ of Yang & Anderson (2018), the streamwise velocity perturbation display a complex behaviour, highlighting the significant lack of flow homogeneity. At $S=0.7$ the rolls are most intense and induce the maximum perturbation at a wall-normal location close to the centres of the rolls. However, the velocity perturbation reaches its peak at the mid-plane only at $S=0.95$ and then changes sign from $S \gtrsim 1.75$. Interestingly, the model never predicts flow reversal when the strip width is increased, in agreement with observations in the literature (Neuhauser et al. 2022).

Overall, the present framework appears to capture the three flow regimes documented in the literature correctly. The implication is that linear mechanisms, whereby secondary flows may be interpreted as the output response of the mean shear flow to a steady forcing localised at the wall, may be sufficient to predict the size and strength of the rolls. Based on the streamwise velocity evaluated at the mid-plane, the boundaries may be located at $S \approx 0.4$ and $S \approx 8$, but differences may arise with alternative criteria that consider the bulk of the flow.

3.3. Eddy viscosity perturbation

In many studies that have examined the properties of the linearised Navier–Stokes equations, velocity perturbations resulting from optimal or stochastic forcing are computed by assuming that the turbulent viscosity distribution is not affected by the flow perturbation and follows an analytical or empirical distributions (Reynolds & Hussain 1972; del Álamo & Jiménez 2006; Hwang & Cossu 2010; Morra et al. 2019; Pickering et al. 2021). In the present case, and unlike previous work, the governing equations must include a transport equation for the turbulent viscosity. It is speculated that the peak location predicted by the present linear model (see figure 4) is determined in large part by the selectivity of the linearised Navier–Stokes operator, rather than by the specifics of the turbulence model adopted. However, the inclusion of such a model and the resulting perturbation of the turbulent viscosity $\nu _t^{(1)}$ is key to capture the influence of the heterogeneous surface roughness and the generation of secondary flows, as described in § 2.2. In practice, spanwise gradients of the eddy viscosity produce, through the QCR model, Reynolds stresses that act as source terms for the streamwise vorticity equation, resulting in secondary motions.

Figure 8 shows contours of the perturbation turbulent viscosity for several strip widths $S$. In figure 8(a), the QCR constant $c_{r1}$ has been set to zero, to characterise the eddy viscosity distribution in the absence of cross-stream motions and highlight the interaction with the streamwise velocity fields of figure 6. In figure 8(b), the QCR constant is $c_{r1} = 0.3$. In general, positive eddy viscosity perturbations are observed above the high-roughness strips, and vice versa, reflecting the boundary condition (2.28) and the altered distance $d$ from the wall. For narrow strips, the eddy viscosity perturbation is confined near the wall, given that rapid spatial variations of the eddy viscosity tend to be damped by the diffusion term in the linearised SA equation. As the strip gets wider, more-intense eddy viscosity distributions are observed, reflecting the increased acceleration/deceleration of the flow over the low and high roughness strips, respectively. Small, but likely significant, changes are observed when secondary flows are ‘turned on’, see figure 8(b). It can be observed that the cross-stream motions produce a further distortion of the eddy viscosity distributions. This is the result of two mechanisms that result from the analysis of the SA transport model: (a) the advection of turbulent viscosity of the background flow operated by the vertical and lateral velocities and (b) the altered production of turbulent viscosity due to the alteration of wall-normal and spanwise gradients of the streamwise velocity.

Figure 8. Contours of the eddy viscosity perturbation $\nu _t^{(1)}$ over roughness strips at $DC=0.5$ for several strip widths $S$, without the QCR model (a) and with the QCR model (b). The field in the fundamental domain (see figure 1) is repeated four and two times for the two narrowest strip cases. Contour lines of the perturbation of the streamfunction $\psi ^{(1)}$ are also reported to visualise the secondary flows. Negative $\psi ^{(1)}$ contours are indicated by dashed lines. The Reynolds number is $Re_{\tau }=1000$, for $k_s^{(0)} = 180$. The darker/lighter patches denote the location of the high-/low-roughness strips.

3.4. Tertiary structures and role of the duty cycle

The present model shows that the duty-cycle appears to play an important role in controlling the formation of tertiary structures, which have so far not been observed for strip-type roughness (Neuhauser et al. 2022). This is shown in figure 9, where contours of the wall-normal and streamwise perturbation velocities are shown for increasing duty cycles, at $\varLambda =2.8$. These cases correspond to the triangular marker in figure 4. For each of these solutions, as well as for other solutions in the interval $DC = [0, 1]$, the wall-normal and streamwise velocity profiles are extracted at the centre of the low-roughness strip. All these profiles are then combined together to form the colour map reported in the bottom panels of figure 9. This visualisation suggests that flow reversal over the low-roughness strip may begin to appear in practice as the duty cycle is decreased to about $0.4$ and would be the most intense for $DC \approx 0.23$ (i.e. for $S_h \approx 0.3 S_l$). Incidentally, the occurrence of flow reversal and tertiary structures explains the two ‘ears’ of the contours of the kinetic energy density map of figure 4(a), where the cross-stream motions are relatively intense, compared to other duty cycles. This analysis also suggests that the vertical velocities over the low-roughness strip are most intense for a duty cycle equal to about 0.8 and not for the symmetric case $DC=0.5$ which the kinetic energy density maps would suggest. Arguably, this can be attributed to the constructive interference of the wall-normal velocities induced by two neighbouring vortices, pushed closer to each other by the decreasing width of the low-roughness strip. The streamwise velocity field is also particularly affected by the duty cycle, as faster or slower flow over the low-roughness strip can be found depending on the $DC$. The transport of fast/slow fluid operated by the cross-stream velocities is particularly visible. For instance, for $DC=0.23$ the left vortex, rotating counter-clockwise transports low-momentum fluid from the near-wall region at $x_3 = 0$ to its right flank, producing a negative velocity streak at $x_3 \approx 0.2$, and similarly for the other longitudinal vortex.

Figure 9. Contours of the wall-normal (a–c) and streamwise (d–f) velocity perturbation for increasing duty cycles, indicated in the figures, and $\varLambda =2.8$. Panels (g) and (h) show the wall-normal and streamwise velocity, respectively, extracted on the vertical line (shown in panels a–f) at the centre of the low-roughness strip $x_3/\varLambda = 0.5$ for duty cycles ranging from 0 (narrow high-roughness strip) to 1 (wide high-roughness strips). The lines labelled (a)–( f) in panels (g) and (h) refer to the corresponding panels. The Reynolds number is $Re_{\tau }=1000$ and $k_s^{(0)} = 180$.

4. Generalising the framework to complex surface heterogeneities

It has been demonstrated that when conducting experiments (Wangsawijaya et al. 2020) or roughness-resolving simulations (Stroh et al. 2020) over realistic heterogeneous rough surfaces it is pivotal to ensure that the shear-increasing effects of roughness are decoupled from the inevitable variation of the mean surface height. Both roughness and elevation heterogeneity produce secondary currents and, therefore, the combination of such effects can have a significant influence on the strength and potentially the direction of the resulting secondary motions (Schäfer et al. 2022; Frohnapfel et al. 2024). In this section, we analyse this aspect through the lens of the linearised model, to initiate the formulation of a unifying framework for flows over complex heterogeneous surfaces.

4.1. Secondary-flow-inducing source mechanisms

In the present linearised framework, in the limit case where the spanwise variation of the roughness or the elevation is small, ridge-type and strip-type roughness are modelled with the same approach. In both cases, the flow–surface interaction develops through three separate source mechanisms corresponding to three different inhomogeneous terms acting as forcing in the linearised equations. To illustrate these mechanisms, it is instructive to examine ridge-type roughness considered in our previous work (Zampino et al. 2022) where the mechanisms are all active. In such case, the lateral variation of the elevation was defined by a unitary peak-to-peak, zero-mean function $f(x_3)$ so that, e.g. the bottom wall of the channel is located at $x_2 = \epsilon f(x_3)$ and the small parameter $\epsilon$ controls the actual amplitude of the topography. The first source mechanism, denoted as $A$ in what follows, is mediated by the linearised boundary condition on the streamwise velocity, e.g. on the lower wall,

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

derived in Zampino et al. (2022). Physically, this condition produces a velocity slip that captures the acceleration and deceleration perceived by the bulk flow above the troughs and the crests of the non-planar topography, respectively. What leads to the formation of secondary flows is the resulting spanwise gradient of the streamwise velocity in the near-wall region. This gradient induce, via the QCR model, spanwise gradients of anisotropic Reynolds stresses (see (2.25)) that then create secondary flows. This mechanism does not seem to have been discussed previously in the literature of secondary flows. More generally, we are not aware of studies that consider the somewhat artificial surface arrangement consisting of longitudinal, flush-mounted belts moving upstream/downstream alternated with regions of solid wall. This set-up would capture source mechanism $A$ directly, and we predict that it could generate relatively intense secondary currents with a downwelling over the upstream-moving belts. When modelling strip-type roughness, this source mechanisms is not active since no-slip boundary conditions for the streamwise velocity are used. This is equivalent to stating that the mean height of the two roughness strips is the same and the boundary of the numerical domain is at some suitable location where the strip-averaged streamwise velocity goes to zero.

The second source mechanism, denoted as $B$, is active for both types of heterogeneity. It is mediated by the destruction term of the SA transport equation, where the inverse of the squared distance $d$ between a point in the numerical domain and the nearest ‘wall’ models the blocking effect of the wall (Spalart & Allmaras 1994; Aupoix & Spalart 2003). With such term, the SA model predicts an accurate log-layer and, thus, lateral perturbations of the distance $d$ induced by the elevation (via the function $f(x_3)$) or by the displacement of the virtual origin (via the term $d_0^{(1)}(x_3)$), produce a perturbation of the log-layer and spanwise gradients of the turbulent viscosity field. Crucially, the sign of this source mechanism is opposite for strip-type and ridge-type roughness: while locally increasing elevations correspond to locally reducing distances $d$, increasing roughness produces increasing distances, as the virtual origin is further displaced downwards beneath the boundary of the numerical domain. In this regard, the framework suggests that the mean roughness height is not the important factor. Rather, the displacement of the virtual origin, a dynamic parameter that ultimately depends on the drag of the surface, is what controls the intensity of the forcing and of the resulting secondary flows. As a side note, while source mechanism $A$ acts as a boundary condition, source mechanism $B$ acts at all wall distances, as it captures the perturbed development of the wall-bounded flow from a different origin.

The third source mechanism, denoted as $B^\prime$, is again active for both types of heterogeneities. It is mediated by the inhomogeneous boundary condition on the modified turbulent viscosity. For ridge-type roughness, this is

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

For strip-type roughness, the lateral variation of the virtual origin $d_0^{(1)}(x_3)$ plays the same role of the function $f(x_3)$ describing the ridge topography, but with an opposite effect as dictated by the boundary condition (2.28). This condition must be applied consistently with source mechanism $B$, so that the shifted eddy viscosity profile produced by such mechanism is consistent with the boundary condition (4.2). Studies on modelling roughness in RANS simulations (Aupoix & Spalart 2003) have shown (and we confirm it later) that this third mechanism is quite weak, because capturing the overall development of the turbulent structure from a different virtual origin is more important than applying non-zero boundary conditions for the turbulent quantities.

4.2. Combining elevation and roughness variations

To examine the relative strength and the combination of these three source mechanisms, we first perform linearised calculations for $Re_\tau = 1000$ for a smooth sinusoidal wall where $f(x_3) = \cos (2{\rm \pi} /\varLambda x_3)$, with period $\varLambda = 1$, by activating only one source mechanism each time. Results are reported in figure 10. The wall-normal velocity profiles, taken at $x_3=0$ over the crest of the topography, show that source mechanism $A$ produces a downwelling flow over the region where the slip velocity is negative. Conversely, the decrease of the distance from the wall over the crest, source mechanism $B$, produces an upwelling of slightly greater magnitude, while source mechanism $B^\prime$ is much weaker than the first two, as discussed. Because the proposed model is linear, the superposition principle applies and the effect of varying simultaneously the roughness properties and the elevation can be obtained easily by combining appropriately solutions obtained in the two cases. For ridge-type roughness, source mechanisms $A$, $B$ and $B^\prime$ are all active, while for strip-type roughness only source mechanisms $B$ and $B^\prime$ should be retained, after inverting the sign of the induced flow given the different orientations of $f(x_3)$ and $d_0^{(1)}$ (see figure 12).

Figure 10. Wall-normal velocity profiles over the crest of a sinusoidal topography, for $\varLambda =1$ and $Re_\tau = 1000$, obtained by activating each of the three secondary-flow-inducing source mechanisms in turn.

To better demonstrate the relative importance of these mechanisms, we then consider a configuration where strips and rectangular ridges are arranged in phase, and the high-roughness regions are placed over the ridges (see figure 11). The width $S_l$ coincides with the gap between the ridges while $S_h$ coincides with the ridge width (denoted as $G$ and $W$ in Zampino et al. 2022). To characterise the relative strength of the two effects, we introduce the parameter $\beta$ as the ratio between the displacement of the virtual origin produced by the strips and the topography, so that

(4.3)\begin{equation} d_0^{(1)}(x_3) = \beta f(x_3), \end{equation}

with the caveat that positive displacements are in different directions depending on the type of heterogeneity. Case $\beta =0$ corresponds to smooth ridges, leading to secondary flows produced by lateral variations of the elevation. As discussed in Zampino et al. (2022), the linearised RANS model predicts an upwelling over the high-elevation regions. On the other hand, case $\beta =1$ corresponds to the combination of the two roughness heterogeneities where the downward displacement of the virtual origin produced by the shear-increasing roughness is, in theory, fully compensated by a increased elevation.

Figure 11. Superposition of high- and low-roughness strips with smooth rectangular ridges modelling the effect of lateral variation of elevation over the bottom wall of the channel. The lateral variation of the virtual origin of the rough surface is captured by $d_0^{(1)}(x_3)$ (positive if the virtual origin is below the mean height at $x_2=0$), while the lateral variation of the elevation is captured by the function $f(x_3)$ (positive if it is above $x_2=0$).

Figure 12 shows results of this analysis, where we conduct calculations over surfaces with combined roughness and elevation as a function of $\beta$. The kinetic energy density obtained at each composite surface is normalised with the reference value at $\beta =0$ and is shown in figure 12(g). Results are shown for $Re_\tau = 1000$ and $k_s^{(0)} = 180$, for strip (and ridge) widths configurations characterised by constant $S_h = S_l = 0.7$. The kinetic energy shows a minimum for $\beta \approx 0.3$, where the cross-stream velocity components vanish and no secondary flows are predicted. The location of the minimum does not seem to depend greatly on the strip configuration (not shown for brevity). For $\beta > 0.3$ the effect of the lateral variation of the roughness becomes dominant and the associated kinetic energy density can be several times higher than the reference value. The resulting flow structure as $\beta$ is increased, i.e. as the effect of roughness is increased, is shown in figures 12(a–f). Maps of the perturbation streamwise and wall-normal velocity components are shown, for $\beta =0$, in the region where the ridge-type roughness is dominant, $\beta =0.45$ close to the minimum of ${\mathcal {K}}(\beta )/{\mathcal {K}}(0)$, and 1, where the strip-type roughness is dominant. For $\beta =0$ (figure 12a,d), the flow topology shows an upwelling over the ridges, as the upward displacement of the wall produced by the ridge-type roughness in this region is dominant. For $\beta =0.45$ (figure 12b,e), secondary currents have changed in direction but their strength is relatively weak. For $\beta =1$ (figure 12c, f), the flow topology shows a downwelling over the high-roughness strip, as the effect of the downward displacement of the virtual origin produced by the roughness prevails in this region over the influence of the increased elevation.

Figure 12. Contours of the perturbation streamwise (a–c) and wall-normal (d–f) velocity components for combinations of ridge-type and strip-type roughness for $S_h = S_l = 0.7$. The parameter $\beta$ is 0.0, 0.45 and 1, as we move from left to right with configurations where the effect of the elevation heterogeneity dominates to cases where the roughness heterogeneity prevails. Normalised kinetic energy density of the combination of both ridge-type and strip-type roughness as a function of the parameter $\beta$ ( f). The Reynolds number is fixed at 1000 and $k_s^{(0)} = 180$.

The key result of this analysis is that lateral variations of the roughness properties and of the elevation do not have the same effect on the strength of the resulting secondary motions, because of the ‘damping’ effect of the spanwise variation of the slip velocity produced by source mechanism $A$, active for heterogeneous elevation but not for heterogeneous roughness. The location of the minimum of ${\mathcal {K}}(\beta )$ suggests that the secondary-flow-inducing effect of roughness is about three times stronger than that of ridges, for the same displacement of the virtual origin in absolute terms. However, this prediction is clearly no better than the prediction of the strength of secondary flows for the two heterogeneities. The comparison with the heterogeneous roughness simulations of Chung et al. (2018) reported earlier suggests that the linear model captures quite well the strength of secondary motions over such surfaces. However, for ridges, recent work (Castro & Kim 2024) has suggested that the linear model over-predicts the intensity of secondary motion for tall ridges, owing to the importance of nonlinear effects near the corners of the ridges, while predictions can be more accurate for short ridges that do not protrude excessively in the wall-bounded flow. Overall, this indicates that the response of the wall-bounded flow to a perturbation of the surface elevation is far from being linear. To the best of the authors’ knowledge, evidence for this claim was perhaps first given in Wang & Cheng (2006) (see their figure 20), who showed how the vertical velocity produced by secondary motions saturates rather quickly as the height of the topography is increased. This last piece of evidence suggests that the secondary-flow-inducing effect of the roughness might be even stronger than what the linear model suggest here, although a precise quantification might require dedicated experimental work. Overall, this could also justify the recent results of Frohnapfel et al. (2024) who considered the same surface arrangement considered in this section, with roughness strips located in phase with the ridges. These authors increased the ridge height but did not observe a reversal of the flow direction above the ridges, which was dominated by the downwelling caused by the roughness. Using their data, we calculate the height of the ridges to be 5.16 % the height of the channel and the downward displacement of the virtual origin to be 1.08 %, from their homogeneous roughness data, for a ratio $\beta \simeq 0.21$. This ratio is clearly to the left of the minimum in figure 12, reinforcing the idea that the response of the wall-bounded flow to finite lateral variations of the elevation are less intense than predicted by a linear model for infinitesimal perturbations.

5. Discussion and concluding remarks

In this paper, we proposed a linearised RANS-based framework to predict the structure of Prandlt's secondary flows of the second kind developing over laterally heterogeneous rough surfaces. The work extends our previous efforts on modelling smooth non-planar surfaces (Zampino et al. 2022), e.g. surfaces with longitudinal ridges. The model couples the linearised RANS equations with the SA transport equation to capture the altered turbulent structure. Rough surfaces with alternating streamwise-aligned strips of high and low roughness are modelled using established RANS modelling strategies available in the literature for homogeneous rough surfaces (Aupoix 2007). Briefly, these strategies adopt a virtual origin framework, whereby the shift of the logarithmic profile is obtained by displacing beneath the boundary of the numerical domain the origin from which turbulent quantities develop (Rotta 1962). This results in altered boundary conditions as well as a domain forcing term when the distance from the wall appears in the turbulence model's transport equations. The framework also employs a nonlinear Reynolds stress model, i.e. the QCR (see Spalart 2000), so that secondary currents induced by the inhomogeneity of anisotropic turbulent stresses can be predicted (Speziale 1982).

There are several aspects that the model predicts remarkably well in agreement with previous observations. One first aspect is the presence of three separate flow regimes as the strip width is increased (Chung et al. 2018). For narrow strips, the linear model supports a flow structure consisting of rolls localised in the vicinity of the wall and having wall-normal size scaling linearly with the strip width. For wider strips, the flow structure tends to an ‘isolated-vortex’ regime, where streamwise vorticity is concentrated in a roughly square region localised around the transition between strips, while the bulk flow in the centre of the roughness strips tends to its homogeneous rough-wall flow dictated by the local roughness properties. In the intermediate regime, secondary currents are most intense when the high- and low-roughness strips have the same widths, and are equal to about $0.7$ of the half-channel height. The model also provides adequate quantitative predictions of the intensity of the cross-stream velocity component, compared with, e.g. the numerical simulations of (Chung et al. 2018). A second aspect concerns tertiary structures and flow reversal, not observed in previous studies on roughness strips that have most often examined high- and low-roughness strips of equal width. The linear model predicts that these phenomena only appear when the strips have different width. For instance, for $\varLambda = 2.8$, flow reversal is strongest on the low-roughness strip when this strip is about four times wider than the high-roughness strip. It would be interesting to confirm this prediction through experiments or simulations. A third aspect concerns the occurrence of low- and high-momentum pathways flanking the longitudinal rolls, where high-speed flow may be found on the high-roughness strip (and vice versa) in regions dominated by the vertical velocities induced by the rolls. Away from the rolls, or for wide strips, the expected relationship between surface roughness and streamwise velocity defect is recovered.

From a practical standpoint, the advantage of the present approach is, undoubtedly, its computational efficiency. However, the ability to rapidly probe the parameter space has enabled progress to be made on a more fundamental standpoint. Specifically, previous work has suggested that secondary currents may be the time-averaged picture of naturally occurring large-scale motions locked in place by the heterogeneity. The robustness of the above-mentioned similarities between the flow structure predicted by the present framework and previous observations suggests an alternative input–output perspective whereby these currents are the output response of the Navier–Stokes operator linearised about the turbulent mean and subjected of a steady, streamwise-independent forcing localised at the wall, associated to the lateral perturbation of the surface characteristics. This perspective complements the well-accepted viewpoint that instantaneous coherent structures in wall-bounded turbulence may be described to a satisfactory degree by the output properties of the linearised operator subjected to a random forcing (Hwang & Cossu 2010; McKeon & Sharma 2010). Admittedly, this viewpoint does not explain the observed interplay between large-scale motions and secondary structures (Wangsawijaya & Hutchins 2022), whereby energy of the former leeches into the latter. One possible explanation that is worth exploring further is that the mean flow distortion produced by the secondary currents may locally alter the selective amplification properties of the linearised operator, producing a spanwise modulation of the nature of large-scale motions that may be interpreted as an energy interaction between secondary currents and large-scale motions.

A second key output of this study is that it offers a unified perspective to examine both ridge-type and strip-type roughness. Examination of these two cases within the present framework shows that secondary motions produced by complex surface heterogeneities, e.g. arbitrary combinations of elevation and roughness properties, may be seen as originating from two separate source mechanisms. The first is a lateral variation of the virtual origin from which the mean turbulence structure develops. The sign of this variation is opposite for ridge-type and strip-type roughness: the virtual origin is shifted upwards by ridges, and downwards by higher roughness. The second source mechanism is mediated by the lateral variation of the streamwise slip velocity in the vicinity of the wall, and the associated spanwise gradients of the streamwise velocity. This source mechanism captures the acceleration/deceleration perceived by the bulk flow above the troughs and crests of a non-planar topography, respectively, or when the mean roughness height varies laterally. In ridge-type roughness, we have shown that this source mechanism damps the first so that the resulting secondary motions are weaker compared with those that would be predicted from the first mechanism. In other words, for the same lateral variation of the virtual origin strip-type roughness produces more intense secondary flows than ridge-type roughness. The caveat is that this perspective applies, in the limit where the lateral variation of surface attributes is ‘small’, in the region of validity of the linear model. For finite-amplitude perturbations of the surface attributes, these predictions would need to be further verified. In any case, the present modelling framework suggests that the mean roughness height, a geometric quantity used in previous studies that considered the coupling of roughness and elevation is not the important quantity to be monitored when investigating combinations of elevation and roughness. Instead, we suggest that the notion of the virtual origin, a dynamical parameter associated to the downward shift of the logarithmic distribution, should be considered.

Like all other analyses of the linearised Navier–Stokes operator, the present approach yields useful insight but has its limitations. The ability of a one-equation turbulence model equipped with a nonlinear Reynolds stress model to capture the unsteady motion of secondary structures may be questioned. A more extensive assessment of alternative turbulence modelling strategies and a comparison with high-fidelity simulations is certainly warranted in future work, although it must be pointed out that the SA-QCR model used here does seem to capture fairly well the anisotropic nature of the Reynolds stress tensor in square-duct flow (Modesti 2020). In this regard, it is speculated that such alternatives may not necessarily lead to significant qualitative changes in the predictions obtained here. In fact, we argue that the selectivity of the linearised Navier–Stokes operator may be dominant, with minimal effects of the turbulence modelling strategy adopted. This speculation is supported by the extensive evidence that useful predictions have been made using linearised Navier–Stokes equation approaches using much simpler turbulence modelling techniques than used here. Such a speculation may be verified simply by examining the response obtained from other turbulence modelling strategies. Given that linearisation of complex models is a tedious task, nonlinear calculations using existing solvers may be performed for sufficiently small roughness variations that the linear approximation is reasonably valid. A further limitation is that it is unclear how sensible is the streamwise independence assumption used here when secondary motions display a strong meandering behaviour (Kevin et al. 2019).

Further, the importance of nonlinear effects neglected in the present linear framework in altering the structure of secondary flows is not clearly understood. Convergence of the expansion (2.21) for finite-amplitude surface perturbations should be examined. Analyses of the streamwise vorticity budget for ridge-type roughness (Castro & Kim 2024) shows that nonlinear effects appear significant near sharp geometric features, and similar conclusions might apply to regions where roughness properties vary sharply. Nonlinear effects must also be introduced in order to predict the change in drag (Zampino 2023), since the one further limitation of the present linear approach is that it offers no insight into the dependence of drag on the surface properties. This would be highly useful for predicting the drag of real-world surfaces, which is currently a topic of active research. Clearly, the question is whether classical turbulence models and established rough-wall treatment strategies may adequately capture the drag characteristics of heterogeneous surfaces. A possible avenue forward would be to bypass these modelling strategies and the small-amplitude assumption and instead adopt the framework reviewed in Zare, Georgiou & Jovanović (2020) to model the second-order statistics of the velocity fluctuations from the response of the linearised Navier–Stokes subjected to a structured stochastic forcing. However, the extension of this framework to rough surfaces would need to be considered first.