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. Reference Zampiron, Cameron and Nikora2020). 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 Reference Perkins1970; Speziale Reference Speziale1982; Bottaro et al. Reference Bottaro, Soueid and Galletti2006), 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 Reference Speziale1982; Chen, Lien & Leschziner Reference Chen, Lien and Leschziner1997). Here we utilise the QCR nonlinear model presented in Spalart (Reference Spalart2000), 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 (Reference Spalart2000). 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 Reference del Álamo and Jiménez2006; Pujals et al. Reference Pujals, Garcìa-Villalba, Cossu and Depardon2009; Hwang & Cossu Reference Hwang and Cossu2010; Morra et al. Reference Morra, Semeraro, Henningson and Cossu2019 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 Reference Spalart and Allmaras1994) 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 Reference Herring and Mellor1968; Spalart & Allmaras Reference Spalart and Allmaras1994), 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 Reference Spalart and Allmaras1994), 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. Reference Durbin, Medic, Seo, Eaton and Song2000; Suga, Craft & Iacovides Reference Suga, Craft and Iacovides2006; Aupoix Reference Aupoix2007; Brereton & Yuan Reference Brereton and Yuan2018; Prakash & Laurendeau Reference Prakash and Laurendeau2020, 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 Reference Rotta1962). Relying on the outer-layer similarity hypothesis of Townsend (Reference Townsend1976), 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 Reference Grigson1992), 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 Reference Aupoix2007). 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 Reference Knoll and Keyes2004), using the ‘hookstep’ approach of Viswanath (Reference Viswanath2007) 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 Reference Reynolds and Hussain1972). 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. Reference Wangsawijaya, Baidya, Chung, Marusic and Hutchins2020). 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. Reference Chung, Monty and Hutchins2018) or by a lateral variation of the transversal slip length (Neuhauser et al. Reference Neuhauser, Schäfer, Gatti and Frohnapfel2022). 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. Reference Wangsawijaya, Baidya, Chung, Marusic and Hutchins2020). 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 Reference Brandt2014). 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. (Reference Zampino, Lasagna and Ganapathisubramani2022)), 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 Reference Perkins1970). 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 Reference Chavarin and Luhar2020; Ran et al. Reference Ran, Zare and Jovanović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. Reference Neuhauser, Schäfer, Gatti and Frohnapfel2022), decays more rapidly with the distance from the wall (Meyers et al. Reference Meyers, Ganapathisubramani and Cal2019). 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 Reference Mejia-Alvarez and Christensen2013; Barros & Christensen Reference Barros and Christensen2014) 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 Reference Chernyshenko and Baig2005), 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. Reference Hutchins, Ganapathisubramani, Schultz and Pullin2023; Frohnapfel et al. Reference Frohnapfel, von Deyn, Yang, Neuhauser, Stroh, Órlú and Gatti2024). 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 Reference Barros and Christensen2014), 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 Reference Zampino2023).

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. (Reference Frohnapfel, von Deyn, Yang, Neuhauser, Stroh, Órlú and Gatti2024), 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. Reference Wangsawijaya, Baidya, Chung, Marusic and Hutchins2020) and numerical simulations (Chung et al. Reference Chung, Monty and Hutchins2018; Neuhauser et al. Reference Neuhauser, Schäfer, Gatti and Frohnapfel2022). 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. (Reference Wangsawijaya, Baidya, Chung, Marusic and Hutchins2020) 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. (Reference Chung, Monty and Hutchins2018), 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. (Reference Neuhauser, Schäfer, Gatti and Frohnapfel2022) 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. Reference Frohnapfel, von Deyn, Yang, Neuhauser, Stroh, Órlú and Gatti2024).

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. (Reference Chung, Monty and Hutchins2018), 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. (Reference Chung, Monty and Hutchins2018) 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. Reference Stroh, Hasegawa, Kriegseis and Frohnapfel2016; Castro & Kim Reference Castro and Kim2024), 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. (Reference Chung, Monty and Hutchins2018) 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. (Reference Chung, Monty and Hutchins2018) 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 Reference Mejia-Alvarez and Christensen2013; Barros & Christensen Reference Barros and Christensen2014; Willingham et al. Reference Willingham, Anderson, Christensen and Barros2014), 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 Reference Jiménez2004), 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. Reference Wangsawijaya, Baidya, Chung, Marusic and Hutchins2020; Frohnapfel et al. Reference Frohnapfel, von Deyn, Yang, Neuhauser, Stroh, Órlú and Gatti2024), 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. (Reference Chung, Monty and Hutchins2018) 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 Reference Bou-Zeid, Parlange and Meneveau2007). 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 (Reference Yang and Anderson2018), 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. Reference Neuhauser, Schäfer, Gatti and Frohnapfel2022).

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 Reference Reynolds and Hussain1972; del Álamo & Jiménez Reference del Álamo and Jiménez2006; Hwang & Cossu Reference Hwang and Cossu2010; Morra et al. Reference Morra, Semeraro, Henningson and Cossu2019; Pickering et al. Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021). 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. Reference Neuhauser, Schäfer, Gatti and Frohnapfel2022). 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.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. Reference Zampino, Lasagna and Ganapathisubramani2022) 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. (Reference Zampino, Lasagna and Ganapathisubramani2022). 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 Reference Spalart and Allmaras1994; Aupoix & Spalart Reference Aupoix and Spalart2003). 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 Reference Aupoix and Spalart2003) 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. Reference Zampino, Lasagna and Ganapathisubramani2022). 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. (Reference Zampino, Lasagna and Ganapathisubramani2022), 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. (Reference Chung, Monty and Hutchins2018) 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 Reference Castro and Kim2024) 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 (Reference Wang and Cheng2006) (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. (Reference Frohnapfel, von Deyn, Yang, Neuhauser, Stroh, Órlú and Gatti2024) 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.