2.1. Nonlinear formulation for the roughness-induced steady flow

We formulate the equations for these two decks, starting with the viscous deck, which is a thin layer with $O(\epsilon ^4)$ width above a curved wall $y^*=h^*F$ . The derivation of the governing equations is facilitated by introducing appropriately rescaled coordinates:

(2.1a–c) \begin{align} x=x^{*}=O(1), \quad Y=\epsilon ^{-4}(y^{*}-h^{*}F)=O(1), \quad Z=\epsilon ^{-3}z^{*}=O(1), \end{align}

where Prandtl’s transformation (Prandtl Reference Prandtl1938) has been implemented to convert the curved boundary to $Y=0$ . The velocities and pressure are scaled as

(2.2) \begin{align} (u^{*},v^{*}-h^{*}u^{*}F_{x^{*}}-h^{*}w^{*}F_{z^{*}},w^{*},p^{*})=(U,\epsilon ^{4} V,\epsilon ^{3} W,\epsilon ^{6} P)(x,Y,Z)+\mbox{h.o.t.}, \end{align}

where h.o.t. stands for high-order terms in the expansion that have not been written out explicitly; such a convention is followed throughout the paper. Substitution of (2.1) and (2.2) into the NS equations leads to the leading-order governing equations:

(2.3a) \begin{align} \frac {\partial U}{\partial x}+\frac {\partial V}{\partial Y}+\frac {\partial W}{\partial Z}&=0, \end{align}

(2.3b) \begin{align} U\frac {\partial U}{\partial x}+V\frac {\partial U}{\partial Y}+W\frac {\partial U}{\partial Z}&=\big(1+h^{2}F^{2}_{Z}\big)\frac {\partial ^{2} U}{\partial Y^{2}}, \end{align}

(2.3c) \begin{align} U\frac {\partial W}{\partial x}+V\frac {\partial W}{\partial Y}+W\frac {\partial W}{\partial Z}&=-\frac {\partial P}{\partial Z}+\big(1+h^{2}F^{2}_{Z}\big)\frac {\partial ^{2} W}{\partial Y^{2}}, \end{align}

plus $\partial P/\partial Y=0$ , indicating that $P=P(x,Z)$ . Here, the subscript $Z$ represents the spanwise derivative, and $h=O(1)$ is the rescaled height, introduced via

(2.4) \begin{equation} h^{*}=\epsilon ^{3}h. \end{equation}

It is worth adding that in the wall-normal momentum equation, the pressure gradient $\partial p/\partial y^* =O(\epsilon ^2)$ in view of (2.1b ) and (2.2), while the inertia is dominated by the centrifugal force due to the streamwise and spanwise curvatures of the roughness arrays, represented by $-u^{*2} h^* F_{x^*x^*}-w^{*2}h^* F_{z^*z^*}$ , which is of $O(\epsilon ^3)$ according to (2.1)–(2.2) and (2.4). It follows that $\partial P/\partial Y =0$ to leading-order accuracy. The no-slip boundary condition and the far-field matching condition are imposed as

(2.5a) \begin{align} &U=V=W=0 \quad \mathrm{at} \quad Y=0; \end{align}

(2.5b,c) \begin{align} &U\rightarrow 1,\quad W\rightarrow W_{\infty } \quad \mathrm{as} \quad Y\rightarrow \infty . \end{align}

The system (2.3) with (2.5) can be solved when the pressure gradient $\partial P/\partial Z$ in (2.3c ) and the velocity $W_{\infty }$ in (2.5c ) are determined by matching with the inviscid-deck solution.

It is worth noting that in (2.3b , c ), the coefficients of the viscous terms are modified to $1+h^{2}F^{2}_{Z}$ , with the extra terms proportional to $h^2F_Z^2$ arising from applying Prandtl’s transformation and being contributed by the spanwise diffusion terms in the momentum equations. This is a geometric effect and may be interpreted as the spanwise modulation of roughness arrays causing enhanced equivalent transverse diffusion relative to the flat-plate case. This is an interesting effect that does not appear to have been noted previously. Although a similar enhanced (modified) transverse diffusion arises due to the streamwise modulation of a 2-D surface roughness (Hanevy et al. Reference Hanevy, Ferguson, Trevelyan and Griffiths2024). In the cited work, for the modification to be a leading-order effect, the roughness height must be very large, of $O(1)$ , whereas in the present problem, the roughness height and spanwise spacing are both of $O(\epsilon ^3)$ (so that the spanwise slope is of $O(1)$ ). Since the extra terms in (2.3b , c ) originate from the spanwise diffusion, the system (2.3) has gone beyond the classical boundary-layer equations. In the streamwise momentum equation (2.3b ), the pressure gradient is absent because it is of $O(\epsilon ^{6})$ , significantly smaller than the $O(1)$ inertial and viscous terms. Consequently, the system (2.3) is streamwise parabolic, and so can be solved numerically using an efficient streamwise marching method. In addition, meshing a domain over a curved boundary is avoided after applying Prandtl’s transformation. Due to these two features, the streaky flow can be obtained at a significantly reduced computation cost than performing DNS.

Prandtl’s transformation in (2.2) indicates that the displacement due to the roughness drives a wall-normal velocity at the outer edge of the boundary layer:

(2.6) \begin{align} v^{*}_{\infty }=\epsilon ^{3} h(F_{x}+W_{\infty } F_{Z}) \quad \mathrm{as} \quad Y\rightarrow \infty . \end{align}

This so-called transpiration velocity forces $O(\epsilon ^{3})$ perturbations outside the boundary layer. The 3-D components reside primarily in the inviscid deck with $O(\epsilon ^3)$ width and attenuate there. However, if the roughness shape contains a 2-D (spanwise-averaged) component, denoted as $F_{x,2D}$ , the induced 2-D disturbance, also of $O(\epsilon ^{3})$ , attenuates more slowly in the extended outer region occupying a vertical thickness of $O(1)$ . We expand and scale the velocities and pressure as

(2.7) \begin{align} (u^{*},v^{*},w^{*},p^{*})=\big (1+\epsilon ^{3}\bar {u}(x,y^{*}),\epsilon ^{3}\bar {v}(x,y^{*}),0,P_{B}+\epsilon ^{3}\bar {p}(x,y^{*})\big )+\mbox{h.o.t}. \end{align}

Substitution of (2.7) into the NS equations leads to the leading-order equations:

(2.8a–c) \begin{align} \frac {\partial \bar {u}}{\partial x}+\frac {\partial \bar {v}}{\partial y^{*}}=0, \quad \frac {\partial \bar {u}}{\partial x}=-\frac {\partial \bar {p}}{\partial x}, \quad \frac {\partial \bar {v}}{\partial x}=-\frac {\partial \bar {p}}{\partial y^{*}}. \end{align}

The matching condition between the outer region and the boundary layer, and the far-field condition are

(2.9a,b) \begin{align} \bar {v}\rightarrow hF_{x,2D}\quad \mathrm{as} \quad y^{*}\rightarrow 0; \quad (\bar u, \bar v, \bar p)\rightarrow 0 \quad \mathrm{as} \quad y^{*}\rightarrow \infty . \end{align}

The inviscid deck outside the boundary layer occupies the region with a width of $O(\epsilon ^3)$ , which is determined by the spanwise spacing of the roughness arrays. The wall-normal velocity in this layer is also of $O(\epsilon ^{3})$ , determined by the transpiration velocity (2.6). Moreover, an order-of-magnitude analysis of the continuity and spanwise momentum equations reveals that the spanwise velocity and pressure are of $O(\epsilon ^{3})$ and $O(\epsilon ^{6})$ , respectively. Thus, we introduce the rescaled wall-normal coordinate $y=\epsilon ^{-3}y^{*}$ , and scale the velocities and pressure as

(2.10) \begin{align} (u^{*},v^{*},w^{*},p^{*})=(1+\epsilon ^{3} u,\epsilon ^{3}v,\epsilon ^{3}w,\epsilon ^{6}p)(x,y,Z)+\mbox{h.o.t}. \end{align}

Substituting (2.10) into the NS equations written in terms of $(x,y,Z)$ results in the simplified equations

(2.11a) \begin{align} \frac {\partial v}{\partial y}+\frac {\partial w}{\partial Z}&=0, \end{align}

(2.11b) \begin{align} \frac {\partial v}{\partial x}+v\frac {\partial v}{\partial y}+w\frac {\partial v}{\partial Z}&=-\frac {\partial p}{\partial y}, \end{align}

(2.11c) \begin{align} \frac {\partial w}{\partial x}+v\frac {\partial w}{\partial y}+w\frac {\partial w}{\partial Z}&=-\frac {\partial p}{\partial Z}, \end{align}

where the streamwise momentum equation is not needed. These equations are strongly nonlinear unlike their counterparts in most multi-layered (triple- or double-deck) structures. As the transverse extent of the present inviscid layer is comparable to the roughness height, the ‘boundary surface’ between the viscous and inviscid decks is located at $y=hF$ . The matching conditions between the two are

(2.12a,b) \begin{align} v= h(F_{x}+W_{\infty }F_{Z}), \quad w= W_{\infty } \quad \mathrm{at} \quad y=hF. \end{align}

The far-field condition reads

(2.13) \begin{align} v\rightarrow hF_{x,2D}, \quad w\rightarrow 0 \quad \mathrm{as} \quad y\rightarrow \infty , \end{align}

where the first indicates the matching with the roughness-induced 2-D component attenuating slowly in the outer region.

Equations (2.3) and (2.11) along with conditions (2.5) and (2.12)–(2.13) form a closed nonlinear system describing the roughness-induced streaky flow. As scaling (2.4) indicates, the system is applicable in principle to roughness arrays with heights greater than the boundary-layer thickness, and to the best of our knowledge has not been derived before. The formulation shows that the elongated roughness arrays distort streamlines thereby generating a transpiration velocity $v$ at the outer edge of the boundary layer; see (2.12a ). The blowing induces the spanwise pressure gradient $p_{Z}$ and velocity $w$ in the inviscid deck, which in turn act on the viscous layer; see (2.3c ) and (2.5c ).

The streaky flow supports viscous instability, whose characteristic streamwise and spanwise wavelengths are of $O(\epsilon ^3)$ . This stability problem is formulated in § 6. Next we consider the numerical methods for solving the nonlinear system governing the streaky flow.

3.1. Inviscid deck

The continuity (2.11a ) allows for introduction of a streamfunction $\psi$ , defined as

(3.1a,b) \begin{align} v=\psi _{Z}, \quad w=-\psi _{y}. \end{align}

Subtracting the derivative of (2.11c ) with respect to $y$ from the derivative of (2.11b ) with respect to $Z$ , followed by substituting (3.1) into the resulting equation, we obtain

(3.2) \begin{align} \frac {\partial {\nabla} ^{2}\psi }{\partial x}-\frac {\partial \psi }{\partial y}\frac {\partial {\nabla} ^{2}\psi }{\partial Z}+\frac {\partial \psi }{\partial Z}\frac {\partial {\nabla} ^{2}\psi }{\partial y}=0. \end{align}

Equation (3.2) indicates that the material derivative $\mathrm{d} \varOmega /\mathrm{d} x=0$ , where the vorticity $\varOmega ={\nabla} ^{2}\psi$ . As $x \rightarrow 0$ , $\varOmega \rightarrow 0$ . It follows that $\varOmega =0$ for all $x\gt 0$ in the region above the surface $y=hF$ , i.e. the inviscid flow field is irrotational and so $\psi$ satisfies the Laplace equation

(3.3) \begin{align} {\nabla} ^{2}\psi =0. \end{align}

In view of the matching condition (2.12a ) or the far-field condition (2.13), which is contingent on the existence of a 2-D component in the roughness arrays, (3.3) admits the solution

(3.4) \begin{align} \psi (y,Z)=a_{0}Z+\sum _{k=-\infty , k \ne 0}^{\infty }a_{k}\mathrm{e}^{-|k|\beta y}\mathrm{e}^{\mathrm{i}k\beta Z}, \end{align}

where the first term represents the 2-D component and $|k|$ ensures attenuation of 3-D components in the far field. The non-periodic term $a_{0}Z$ amounts to a wall-normal velocity $a_0$ , which is a constant in the present inviscid deck, consistent with the fact that this 2-D velocity attenuates slowly over much longer vertical distances. Note that the streamfunction $\psi$ is non-periodic, but the physical quantities, the velocity and pressure, are periodic. The solution (3.4) can be recast in terms of the transformed wall-normal coordinate $\bar {y}=y-hF$ :

(3.5) \begin{align} \psi (\bar {y},Z)&=a_{0}Z+\sum _{k=-\infty , k \ne 0}^{\infty }a_{k}\mathrm{e}^{\mathrm{i}k\beta Z-|k|\beta hF}\mathrm{e}^{-|k|\beta \bar {y}}\nonumber \\ &=a_{0}Z+\sum _{n=-\infty }^{\infty }\sum _{k=-\infty , k \ne 0}^{\infty }a_{k}d_{kn}\mathrm{e}^{-|k|\beta \bar {y}}\mathrm{e}^{\mathrm{i}n \beta Z}, \end{align}

where $d_{kn}$ are obtained by expanding $\mathrm{e}^{\mathrm{i}k\beta Z-|k|\beta hF}$ into the Fourier series

(3.6) \begin{align} \mathrm{e}^{\mathrm{i}k\beta Z-|k|\beta hF}=\sum _{n=-\infty }^{\infty }d_{kn}\mathrm{e}^{\mathrm{i}n \beta Z}. \end{align}

In terms of the shifted coordinate $(\bar {y}$ , $Z)$ , the boundary condition (2.12) becomes

(3.7) \begin{align} \frac {\partial \psi }{\partial Z}=hF_{x} \quad \mathrm{at} \quad \bar {y}=0, \end{align}

which is directly determined by $h$ and $F$ , independent of the hitherto unknown $W_{\infty }$ . This leads to decoupling of the numerical procedures for the inviscid and viscous decks. Substituting (3.5) into (3.7) and expanding $F_{x}(x$ , $Z)$ into a Fourier series, $F_{x}(x,Z)=\sum f_{x,n}\mathrm{e}^{\mathrm{i}n\beta Z}$ , we find that $a_0=hf_{x,0}=hF_{x,2D}$ , and for $n\neq 0$

(3.8) \begin{align} \mathrm{i}n\beta \sum _{k=-\infty , k \ne 0}^{\infty }a_{k}d_{kn}=hf_{x,n}. \end{align}

The coefficients $a_{k}$ ( $k\neq 0$ ) are obtained by solving the truncated system (3.8) with $-K\leqslant k \leqslant K$ , where $K\geqslant 16$ is found to ensure adequate accuracy of the solution. The streamfunction $\psi$ is then reconstructed using (3.5). The velocities $v$ and $w$ can be computed from (3.1), and in particular $W_{\infty }=w|_{y=hF}$ is found. The pressure gradients, $\partial p/\partial y$ and $\partial p/\partial Z$ , are determined using (2.11b , c ). As mentioned earlier, the values of $w$ and $\partial p/\partial Z$ at the ‘boundary surface’ are used to solve the viscous-deck equations (2.3). It is worth noting that the inviscid flow is controlled by $f_x$ , the streamwise slope of the roughness, as (3.8) indicates. It is thus expected that there would be a reversal of the signs of $w$ and $\partial p/\partial Z$ past the peak of the roughness.

3.2. Viscous streaky boundary layer

The system of boundary-layer equations, (2.3), does not admit a self-similar solution due to the spanwise modulation caused by roughness arrays. Nevertheless, it is convenient to rearrange the equations in (2.3) using the similarity variable and the streamfunctions:

(3.9a–c) \begin{align} \eta =x^{-1/2}Y,\quad \frac {\partial f}{\partial \eta }=U,\quad \frac {\partial g}{\partial \eta }=W. \end{align}

Substitution of (3.9) into (2.3a ) followed by integration with respect to $\eta$ yields

(3.10) \begin{align} x\frac {\partial f}{\partial x}-\frac {1}{2}\eta U+\frac {1}{2}f +x^{1/2}V+x\frac {\partial g}{\partial Z}=0, \end{align}

where use has been made of the no-slip condition and the fact that $f=g=0$ at $\eta =0$ . Substituting (3.9) into (2.3b , c ) and using (3.10) to eliminate $V$ , we have

(3.11a) \begin{align} &xU\frac {\partial U}{\partial x}-\left (x\frac {\partial f}{\partial x}+x\frac {\partial g}{\partial Z}+\frac {1}{2}f \right )\frac {\partial U}{\partial \eta }+xW\frac {\partial U}{\partial Z}=\big(1+h^{2}F^{2}_{Z}\big)\frac {\partial ^{2} U}{\partial \eta ^{2}}, \end{align}

(3.11b) \begin{align} &xU\frac {\partial W}{\partial x}-\left (x\frac {\partial f}{\partial x}+x\frac {\partial g}{\partial Z}+\frac {1}{2}f \right )\frac {\partial W}{\partial \eta }+xW\frac {\partial W}{\partial Z}=-xP_{Z}+\big(1+h^{2}F^{2}_{Z}\big)\frac {\partial ^{2} W}{\partial \eta ^{2}}. \end{align}

Furthermore, the roughness-induced streaky flow is decomposed as

(3.12) \begin{align} (f,g,U,W)=(f_{B},0,U_{B},0)+(f_{1},g_{1},U_{1},W_{1}), \end{align}

where the first part represents the Blasius solution, which satisfies the equations

(3.13a,b) \begin{align} -\frac {1}{2}f_{B}\frac {\partial U_{B}}{\partial \eta }=\frac {\partial ^{2} U_{B}}{\partial \eta ^{2}}, \quad \frac {\partial f_{B}}{\partial \eta }=U_{B}, \end{align}

with the boundary conditions $f_{B}(0)=U_{B}(0)=0$ and $U_{B}(\infty )=1$ . The second part in (3.12) denotes the roughness-induced flow with $U_{1}=\partial f_{1}/\partial \eta$ and $W_{1}=\partial g_{1}/\partial \eta$ . Substitution of (3.12) and (3.13) into (3.11) leads to the coupled equations

(3.14a) \begin{align} &xU_{B}\frac {\partial U_{1}}{\partial x}-\frac {1}{2}f_{B}\frac {\partial U_{1}}{\partial \eta }-\frac {\partial ^{2} U_{1}}{\partial \eta ^{2}}=N(U_{1})+L(U_{B}), \end{align}

(3.14b) \begin{align} &xU_{B}\frac {\partial W_{1}}{\partial x}-\frac {1}{2}f_{B}\frac {\partial W_{1}}{\partial \eta }-\frac {\partial ^{2} W_{1}}{\partial \eta ^{2}}=-xP_{Z}+N(W_{1}), \end{align}

where

(3.15a) \begin{align} &N(\bullet )=L(\bullet ) -xU_{1}\frac {\partial \bullet }{\partial x}-xW_{1}\frac {\partial \bullet }{\partial Z}, \end{align}

(3.15b) \begin{align} & L(\bullet )=h^{2}F^{2}_{Z}\frac {\partial ^{2} \bullet }{\partial \eta ^{2}}+\left (\frac {1}{2}f_{1}+x\frac {\partial f_{1}}{\partial x}+x\frac {\partial g_{1}}{\partial Z}\right )\frac {\partial \bullet }{\partial \eta }. \end{align}

The boundary and matching conditions in terms of $\eta$ are specified as

(3.16a,b) \begin{align} U_{1}=W_{1}=0 \quad \mathrm{at} \quad \eta =0;\quad (U_{1},W_{1})\rightarrow (0,\left . w \right |_{\bar {y}=0}) \quad \mathrm{as\ } \quad \eta \rightarrow \infty . \end{align}

After rearranging and decomposing, the resulting equations (3.14) have better numerical stability and are convenient for implementing a streamwise marching scheme.

The parabolic system (3.14) is to be solved with an appropriate upstream condition. For this purpose, we consider roughness arrays whose rescaled height $F(x, Z) \rightarrow 0$ as $x\rightarrow 0$ . For simplicity, the wall shape is taken to be of the variable-separation form:

(3.17) \begin{align} F(x,Z)=F_{\!f}(x) F_g(Z). \end{align}

Roughness arrays of the above form with suitable $F_{\!f}$ and $F_g$ may represent those in experiments (e.g. Downs & Fransson Reference Downs and Fransson2014). The initial conditions for two forms of $F_{\!f}$ and their constructions are provided in Appendix A. It is worth pointing out the mathematical formulation is valid, and numerical solutions can be obtained, for a general roughness shape $F(x,Z)$ , which always assumes the form (3.17) in the upstream limit, as can be shown by Taylor expansion of $F(x,Z)$ for $x\ll 1$ .

Given that the roughness arrays are spanwise-periodic, so is the induced streaky flow. Its solution $(f_{1}$ , $g_{1}$ , $U_{1}$ , $W_{1})$ is thus expanded into Fourier series:

(3.18) \begin{align} (f_{1},g_{1},U_{1},W_{1})=\sum _{n=-\infty }^{\infty }(\hat {f}_{n}, \hat {g}_{n}, \hat {U}_{n}, \hat {W}_{n})\exp (\mathrm{i}n\beta Z), \end{align}

where $\beta$ represents the spanwise wavenumber of the roughness arrays so that the spanwise wavelength $\varLambda =2\pi /\beta$ . Substituting (3.18) into (3.14) results in

(3.19a) \begin{align} &xU_{B}\frac {\partial \hat {U}_{n}}{\partial x}-\frac {1}{2}f_{B}\frac {\partial \hat {U}_{n}}{\partial \eta }-\frac {\partial ^{2} \hat {U}_{n}}{\partial \eta ^{2}}=\left [\mathcal{F}\!\left (N(U_{1})+L(U_{B})\right )\right ]_{n}, \end{align}

(3.19b) \begin{align} &xU_{B}\frac {\partial \hat {W}_{n}}{\partial x}-\frac {1}{2}f_{B}\frac {\partial \hat {W}_{n}}{\partial \eta }-\frac {\partial ^{2} \hat {W}_{n}}{\partial \eta ^{2}}=\left [\mathcal{F}\!\left (-xP_{Z}+N(W_{1})\right )\right ]_{n}, \end{align}

(3.19c,d) \begin{align} &\hat {U}_{n}=\frac {\partial \hat {f}_{n}}{\partial \eta },\quad \hat {W}_{n}=\frac {\partial \hat {g}_{n}}{\partial \eta }, \end{align}

where $\mathcal{F}$ denotes a discrete-Fourier-transformed quantity. Equations (3.19) are truncated by retaining $2K$ terms in the Fourier series, and discretised in the streamwise and wall-normal directions using the Crank–Nicolson scheme.

The solution procedure for solving (3.19) is as follows. First, the right-hand terms (except the pressure gradient) in (3.19a , b ) at the second streamwise position are estimated using the velocity field at the initial streamwise position, given by the similarity solution (A9). Equations (3.19b ) and (3.19a ) are solved sequentially to obtain $W_{1}$ and $U_{1}$ on each $\eta {-}Z$ slice. Then, (3.19c , d ) are used to solve for $f_{1}$ and $g_{1}$ . After the velocities and streamfunctions are predicted, the right-hand terms except the pressure gradient are evaluated again in physical space and subsequently transformed back to spectral space, where an enhanced version of the so-called $2/3$ rule (Orszag Reference Orszag1971) is used to eliminate aliasing errors. Next, due to nonlinearity, the velocities are updated iteratively until the difference between consecutive estimates falls below a predefined tolerance of $10^{-8}$ . This iterative process yields the solution at each streamwise position. The computational domain is $[0.01$ , $2]\times [0$ , $40]$ with step sizes of $\Delta x=0.001$ and $\Delta \eta =0.05$ . In the spanwise direction, use of $2K=64$ Fourier modes proves to suffice for capturing the nonlinear effect at most locations. However, close to the boundary-layer separation, the number of Fourier modes needs to be increased to $2K=256$ . Calculations on different grids are performed to ensure that the numerical results are mesh-independent. The running time for each case is approximately six hours (on an Intel i7-10700 8-core processor).

5.1. Inviscid-deck solution

Figure 2.

Contours of ( $a$ ) the spanwise slip velocity $w$ and ( $b$ ) the spanwise surface pressure gradient $p_{Z}$ for $h=0.11$ , $d=0.2$ and $\beta =0.455\pi$ .

As was shown earlier, the solution in the inviscid deck provides the values of $w$ and $p_{Z}$ at the ‘boundary surface’, which represent the spanwise slip velocity and surface pressure gradient, respectively, and they drive the streaky boundary-layer flow, as is indicated by matching conditions (3.16b ) and spanwise momentum (2.3c ). Contours of these quantities on the $(x, Z)$ plane are displayed in figure 2, indicating that $w$ and $p_{Z}$ on the surface are antisymmetric about $Z/\varLambda =0$ . Downstream of the peak $(x\gt 1)$ , the positive and negative signs of the spanwise velocity and pressure gradient for $Z/\varLambda \lt 0$ and $Z/\varLambda \gt 0$ , respectively, mean that the fluid is brought towards the plane $Z/\varLambda =0$ , and an ‘adverse pressure gradient’ is present on either side of $Z/\varLambda =0$ . As a result, a collision type of separation arises as is shown later. Interestingly, the spanwise surface pressure gradient and slip velocity also exhibit approximate symmetry and antisymmetry about $x = 1$ , respectively, which is due to the fact that the roughness element shape is approximately symmetric about $x = 1$ .

5.2. Streaky boundary-layer flow

The spanwise-periodic roughness arrays cause the initially 2-D background (Blasius) flow to become a 3-D streaky flow, the formation and streamwise evolution of which are illustrated by contours of the roughness-induced streamwise velocity displayed in figure 3. Upstream of the peak ( $x\lt 1$ ), low-speed streaks emerge around $Z/\varLambda =\pm 1/2$ , while a high-speed streak forms along $Z/\varLambda =0$ near the wall. From $x=1$ to $1.1$ , the high-speed region ‘elongates’ in the spanwise direction and splits into two, which outflanks a faint low-speed streak that forms in between. From $x=1.1$ to $1.2$ , the low- and high-speed streaks intensify. Downstream of $x=1.2$ , the central low-speed streak concentrates in the spanwise direction and is lifted up in the wall-normal direction, so that a cusp starts to form in the contours. This feature suggests that a collision type of boundary-layer separation is imminent downstream. With the streamwise pressure gradient absent, the separation is caused by the ‘adverse’ spanwise pressure gradient on either side of $Z=0$ , probably similar to that described in Stewartson & Simpson (Reference Stewartson and Simpson1982) and Goldstein & Leib (Reference Goldstein and Leib1993). The cusp signals shortening of the spanwise length scale and formation of a singularity of focusing kind. The present classical boundary-layer equations cannot describe the flow in the separation region. This is the limitation of the present theory. A detailed analysis of the separation region is required to identify first the singularity structure as well as the physical mechanism that may remove the singularity, and then develop the theory capable of describing the separated flow. This task is beyond the scope of the present study, but we consider it an important topic for future investigation. Our view is that the singularity is likely to be resolved by a local asymptotic theory that accounts for full spanwise diffusion or ellipticity, and it follows that the flow upstream of the separation would not deviate significantly from the results obtained by marching. This viewpoint is supported by the fact that solutions to the boundary-region equations remain regular despite these equations being parabolic in the streamwise direction (Goldstein et al. Reference Goldstein, Sescu, Duck and Choudhari2016; Xu, Zhang & Wu Reference Xu, Zhang and Wu2017). Figure 4 shows the velocity vectors ( $V,W$ ) in the cross-flow plane at the streamwise positions corresponding to those in figure 3. The high-speed streaks along $Z/\varLambda =0$ upstream ( $x\lt 1$ ) are associated with the downward motion on both sides. This switches to upwash in the region downstream $(x\gt 1)$ , leading to formation of the low-speed streak. This reversal of the cross-flow velocity is a direct consequence of $hF_{x}$ , which controls the spanwise slip velocity and pressure gradient as was noted in § 3.1. The change of the sign of $f_x$ at the peak position $(x=1)$ leads to the switchover from downwash to upwash.

Figure 3.

Contours of the streaky flow at (a–f) different streamwise positions for $h=0.11$ , $d=0.2$ and $\beta =0.455\pi$ . Curves: solid, $U_{1}\gt 0$ ; dashed, $U_{1}\lt 0$ ; dotted, $U=U_{B}+U_{1}$ .

Figure 4.

Velocity vectors in the cross-flow plane at ( $a$ ) $x=0.7$ , ( $b$ ) $x=1$ , ( $c$ ) $x=1.1$ and ( $d$ ) $x=1.4$ . The colours represent the modulus of the velocity vector.

The streak amplitude and the boundary-layer thickness are two important physical quantities of the streaky flow. The former is defined using the integral streak amplitude as (Shahinfar et al. Reference Shahinfar, Fransson, Sattarzadeh and Talamelli2013)

(5.1) \begin{align} A_{ST}(x)=\frac {1}{\varLambda }\int _{-\varLambda /2}^{\varLambda /2}\int _{0}^{\infty }\big |U(x,\eta ,Z)-U^{Z}(x,\eta )\big |\,\mathrm{d}\eta \,\mathrm{d}Z, \end{align}

where the superscript ‘ $Z$ ’ denotes a spanwise-averaged quantity, while the latter is taken to be the dimensionless displacement thickness $\delta$ , defined as

(5.2) \begin{align} \delta (x,Z)=\int _{h^{*}F}^{\infty }\left [1-U(x^{*},y^{*},z^{*})\right ]\,\mathrm{d}y^{*}=\textit{Re}^{-1/2}x^{1/2}\int _{0}^{\infty }\left [1-U(x,\eta ,Z)\right ]\,\mathrm{d}\eta . \end{align}

Figure 5( $a$ ) shows the streak amplitude versus $x$ . Interestingly, $A_{ST}$ does not follow the simple trend of monotonic increase–decrease of the roughness height. Instead, two local maxima form at $x=0.92$ and $1.27$ , which are upstream and downstream of the roughness-array peak $(x=1)$ , respectively, while close to the peak height, a local minimum occurs at $x=1.07$ . Further downstream of the roughness elements, the spanwise pressure gradient gradually diminishes but the streak amplitude changes slowly, implying that streaks may persist over a distance much longer than do the roughness arrays. Figure 5( $b$ ) shows the streamwise evolution of the displacement thickness along $Z/\varLambda =0$ , $1/4$ and $1/2$ , together with the result of the Blasius solution. Upstream of $x=1.07$ , the largest displacement thickness is found along $Z/\varLambda =1/2$ , where the distortion concentrates, while the displacement thickness hardly alters along $Z/\varLambda =0$ . In contrast, downstream of $x=1.07$ , the largest displacement thickness is observed along $Z/\varLambda =0$ . This is intrinsically associated with collision-type separation, signalled by the formation of the cusp, shown in figure 3.

Figure 5.

The streamwise development of ( $a$ ) the streak amplitude $A_{ST}$ and ( $b$ ) the non-dimensional displacement thickness $\delta$ for $h=0.11$ , $d=0.2$ , $\beta =0.455\pi$ and $\epsilon =0.2036$ .

Figure 6.

Contours of ( $a$ ) the total wall shear $\lambda _{u}(x,Z)=\partial U/\partial Y(x,0,Z)$ and ( $b$ ) the streak-associated wall shear $\Delta \lambda _{u}(x,Z)=\partial U_{1}/\partial Y(x,0,Z)$ , on the $(x,Z)$ plane for $h=0.11$ , $d=0.2$ and $\beta =0.455\pi$ .

Similar to previous studies (Kátai & Wu Reference Kátai and Wu2020; Zhu & Wu Reference Zhu and Wu2022; Zheng & Wu Reference Zheng and Wu2024), the wall shear controls the lower-branch viscous instability at leading order, but the latter is now also influenced by the shape of the roughness arrays. Figure 6 shows contours of the total wall shear $\lambda _{u}(x,Z)=\partial U/\partial Y(x,0,Z)$ and the deviation from the Blasius flow, $\Delta \lambda _{u}(x,Z)=\partial U_{1}/\partial Y(x,0,Z)$ . Upstream of the peak, a maximum is observed in the central region, with a minimum on either side. Downstream of the peak, the minimum occurs along the centreline ( $Z/\varLambda =0$ ), while the maximum appears around $Z/\varLambda =\pm 1/4$ .

The characteristics of the wall shear are further analysed in the supplementary material available at https://doi.org/10.1017/jfm.2026.11174, where $\lambda _u$ is found to undergo spanwise concentration towards $Z/\varLambda =0$ . The development of the Fourier spectral components of $\lambda _u$ indicates that even though only the fundamental component is present in the roughness shape function, harmonics are generated downstream in $\lambda _u$ due to nonlinear effects, and higher harmonics eventually gain greater amplitudes causing spanwise focusing and eventual separation of collision type.

6.1. Lower deck

The lower deck of the viscous instability is a very thin layer adjacent to the curved surface with an $O(\epsilon ^{5})$ width. Thus the local wall-normal coordinate is introduced as

(6.2) \begin{align} \bar {Y}=\epsilon ^{-5}(y^{*}-h^{*}F)=O(1). \end{align}

The flow is decomposed into the streaky base flow and a small-amplitude disturbance, namely

(6.3) \begin{align} (u^{*},v^{*},w^{*},p^{*})=&\,(\bar {U}_{b},\epsilon ^{4}\bar {V}_{b},\epsilon ^{3} \bar {W}_{b},\epsilon ^{6} \bar {P}_{b})\nonumber \\ &+\epsilon _{d}(\epsilon \bar {U}_{d},\epsilon ^{3}\bar {V}_{d}+\epsilon hF_{Z}\bar {W}_{d},\epsilon \bar {W}_{d},\epsilon ^{2} \bar {P}_{d})E+\mathrm{c.c.}, \end{align}

where the subscripts ‘ $b$ ’ and ‘ $d$ ’ denote the streaky base flow and the disturbance, respectively, and ‘c.c.’ stands for complex conjugate. The disturbance has a small amplitude $\epsilon _{d}\ll 1$ and is of modal form proportional to $E=\exp [\mathrm{i}\alpha X -\mathrm{i}\omega t]$ , where $\alpha =O(1)$ and $\omega =O(1)$ are the normalised wavenumber and frequency, respectively.

The lower deck for the instability mode with an $O(\epsilon ^{5})$ thickness is much thinner than the streaky boundary layer, whose thickness is of $O(\epsilon ^{4})$ . Therefore, we approximate the base-flow velocities by their Taylor series expansions about the wall:

(6.4a–c) \begin{align} \bar {U}_{b}=\epsilon \lambda _{u}\bar {Y}+\cdots\!, \quad \bar {V}_{b}=\epsilon ^{2} \lambda _{v}\bar {Y}^{2}+\cdots\!, \quad \bar {W}_{b}=\epsilon \lambda _{w}\bar {Y}+\cdots\!, \end{align}

where $\lambda _{u}$ and $\lambda _{w}$ are the streamwise and spanwise wall shears, respectively, and $\lambda _{v}=-(\partial \lambda _{u}/\partial x+\partial \lambda _{w}/\partial Z)/2$ . Substituting (6.1), (6.3) and (6.4) along with (2.1c ) into the NS equations and linearising, we obtain the lower-deck disturbance equations:

(6.5a) \begin{align} &\mathrm{i}\alpha \bar {U}_{d}+\frac {\partial \bar {V}_{d}}{\partial \bar {Y}}+\frac {\partial \bar {W}_{d}}{\partial Z}=0, \end{align}

(6.5b) \begin{align} &(\mathrm{i}\alpha \lambda _{u}\bar {Y}-\mathrm{i}\omega ) \bar {U}_{d}+\lambda _{u}\bar {V}_{d}+\lambda _{uZ}\bar {Y}\bar {W}_{d}=-\mathrm{i}\alpha \bar {P}_{d}+(1+h^{2}F_{Z}^{2})\frac {\partial ^{2} \bar {U}_{d}}{\partial \bar {Y}^{2}}, \end{align}

(6.5c) \begin{align} &(\mathrm{i}\alpha \lambda _{u}\bar {Y}-\mathrm{i}\omega )\bar {W}_{d}=-\frac {\partial \bar {P}_{d}}{\partial Z}+(1+h^{2}F_{Z}^{2})\frac {\partial ^{2} \bar {W}_{d}}{\partial \bar {Y}^{2}}, \end{align}

along with $\partial \bar {P}_{d}/\partial \bar {Y}=0$ . The non-parallelism associated with the slow variable, $\partial /\partial X$ , is of $O(\epsilon ^{3})$ and thus negligible in the disturbance momentum equations at leading order for the elongated roughness elements. Given that the wall-normal length scale in the lower deck is much smaller than the streamwise and spanwise length scales, the wall-normal diffusion is a leading-order effect, as usual, but as with the streaky base flow the spanwise diffusion plays a leading role in the streamwise and spanwise momentum equations through the coordinate shift (6.2); $\partial ^2(\bar U_d, \bar W_d)/\partial z^{*2}$ together with (6.2) contributes to the extra terms $\partial ^2(\bar U_d, \bar W_d)/\partial \bar Y ^{2}$ in (6.5b ) and (6.5c ). As a result, the coefficients of the viscous terms are modified to $1+h^{2}F^{2}_{Z}$ . This geometric effect arises due to the roughness arrays having an $O(1)$ spanwise slope, and its presence stands as a main difference from the previous work (Kátai & Wu Reference Kátai and Wu2020; Zheng & Wu Reference Zheng and Wu2024).

6.2. Main-deck solution

In the main deck, where the shifted coordinate $Y=O(1)$ , the velocities and pressure expand as

(6.6) \begin{align} (u^{*},v^{*},w^{*},p^{*})=&\,(U_{b},\epsilon ^{4}V_{b},\epsilon ^{3} W_{b},\epsilon ^{6} P_{b})\nonumber \\ &+\epsilon _{d}\Big (\epsilon U_{d},\epsilon ^{2}(V_{d}+hF_{Z}W_{d}),\epsilon ^{2} W_{d},\epsilon ^{2} P_{d}\Big )E+\mathrm{c.c.} \end{align}

Note that Prandtl’s transformation is applied as indicated by the expression for the wall-normal velocity of the disturbance. Substitution of (6.1), (6.6) and (2.1b , c ) into the NS equations leads to the leading-order equations for the disturbance in the main deck:

(6.7a–c) \begin{align} \mathrm{i}\alpha U_{d}+\frac {\partial V_{d}}{\partial Y}=0, \quad \mathrm{i}\alpha U_{b}U_{d}+V_{d}\frac {\partial U_{b}}{\partial Y}=0,\quad \mathrm{i}\alpha U_{b} W_{d}+\frac {\partial P_{d}}{\partial Z}=0, \end{align}

plus $\partial P_{d}/\partial Y=0$ . The solution is found as

(6.8a–d) \begin{align} U_{d}=A_{1}\frac {\partial U_{b}}{\partial Y},\quad V_{d}=-\mathrm{i}\alpha A_{1}U_{b},\quad W_{d}=-\frac {1}{\mathrm{i}\alpha U_{b}}\frac {\partial P_{d}}{\partial Z}, \quad P_{d}=\bar {P}_{d}(Z), \end{align}

where $A_{1}(X,Z)$ is the displacement function. Note that in terms of the shifted coordinate (2.1b ) and transformed $V_d$ , the solution remains of the same form as that in the flat-plate case.

6.3. Upper-deck analysis

In the upper deck, the velocities and pressure expand as

(6.9) \begin{align} (u^{*},v^{*},w^{*},p^{*})=(1,\epsilon ^{3}v_{b},\epsilon ^{3}w_{b},\epsilon ^{6}p_{b})+\epsilon _{d}\epsilon ^{2}( u_{d}, v_{d},w_{d}, p_{d})E+\mathrm{c.c.} \end{align}

Substituting (6.1), (6.9) and $y=\epsilon ^{-3}y^{*}$ into the NS equations, we obtain

(6.10a–d) \begin{align} \mathrm{i}\alpha u_{d}+\frac {\partial v_{d}}{\partial y}+\frac {\partial w_{d}}{\partial Z}=0,\quad u_{d}=-p_{d},\quad \mathrm{i}\alpha v_{d}=-\frac {\partial p_{d}}{\partial y},\quad \mathrm{i}\alpha w_{d}=-\frac {\partial p_{d}}{\partial Z}, \end{align}

which can be rearranged into a single equation for $p_{d}$ :

(6.11) \begin{align} -\alpha ^{2} p_{d}+\frac {\partial ^{2} p_{d}}{\partial y^{2}}+\frac {\partial ^{2} p_{d}}{\partial Z^{2}}=0. \end{align}

The matching conditions read

(6.12a,b) \begin{align} \left .\frac {\partial p_{d}}{\partial y}\right |_{y\rightarrow hF}\rightarrow -\alpha ^{2}A_{1}+hF_{Z}\frac {\partial P_{d}}{\partial Z}, \quad \left .p_{d}\right |_{y\rightarrow hF}\rightarrow P_{d}(Z), \end{align}

where the first condition follows from (6.10c ) and matching the wall-normal and spanwise velocities in the upper deck, $v_{d}$ and $w_{d}$ , with their main-deck counterparts given in (6.6) and (6.8b , c ).

Somewhat surprisingly, the roughness arrays affect the streaky flow and its instability in different manners. Roughness arrays displace streamlines to induce a pressure and a spanwise slip velocity, which act on the boundary layer, driving a spanwise motion. The latter then generates streaks. The relation between the steady perturbations in the main and outer layers is ‘hierarchical’. In contrast, for the instability, the relation is an interactive one: the wall-normal and spanwise velocities of the mode in the main layer affect the pressure generation in the upper layer. The pressure in turn drives the mode in the lower deck.

6.4. The lower-deck solutions

The lower-deck disturbance equations are solved following similar steps in Smith (Reference Smith1979b ), Walton (Reference Walton1996) and Walton & Patel (Reference Walton and Patel1998). We introduce the rescaled variable

(6.13) \begin{align} \zeta =(1+h^{2}F_{Z}^{2})^{-1/3}(\mathrm{i}\lambda _{u}\alpha )^{1/3}\bar {Y}+\zeta _{0}, \end{align}

where $\zeta _{0}=-\mathrm{i}\omega (1+h^{2}F_{Z}^{2})^{-1/3} (\mathrm{i}\lambda _{u}\alpha )^{-2/3}$ .

Rewritten in terms of $\zeta$ , the spanwise momentum (6.5c ) is simplified to

(6.14) \begin{align} \frac {\partial ^{2} \bar {W}_{d}}{\partial \zeta ^{2}}-\zeta \bar {W}_{d}=(1+h^{2}F_{Z}^{2})^{-1/3}(\mathrm{i}\lambda _{u}\alpha )^{-2/3}\frac {\partial \bar {P}_{d}}{\partial Z}, \end{align}

which satisfies the no-slip and matching conditions:

(6.15a,b) \begin{align} \hspace {-0.5em} \bar {W}_{d}=0 \quad \mathrm{at} \hspace {0.5em} \zeta =\zeta _{0}; \hspace {0.5em} \bar {W}_{d}\rightarrow -(1+h^{2}F_{Z}^{2})^{-1/3}(\mathrm{i}\lambda _{u}\alpha )^{-2/3}\zeta ^{-1}\frac {\partial \bar {P}_{d}}{\partial Z}\quad \mathrm{as} \hspace {0.5em} \zeta \rightarrow \infty . \end{align}

The solution to (6.14)–(6.15) is given by

(6.16) \begin{align} \bar {W}_{d}=\pi \big(1+h^{2}F_{Z}^{2}\big)^{-1/3}(\mathrm{i}\lambda _{u}\alpha )^{-2/3}\frac {\partial \bar {P}_{d}}{\partial Z}\left [\frac {\mathrm{Gi}(\zeta _{0})}{\mathrm{Ai}(\zeta _{0})}\mathrm{Ai}(\zeta )-\mathrm{Gi}(\zeta )\right ]\!. \end{align}

Differentiating (6.5b ) with respect to $\bar {Y}$ to eliminate the pressure and using (6.5a ) to eliminate $\bar {V}_{d}$ in the resulting equation, we obtain the equation governing $\bar {U}_{d}$ :

(6.17) \begin{align} (\mathrm{i}\alpha \lambda _{u}\bar {Y}-\mathrm{i}\omega )\frac {\partial \bar {U}_{d}}{\partial \bar {Y}}-\lambda _{u}\frac {\partial \bar {W}_{d}}{\partial Z}+\lambda _{uZ}\bar {W}_{d}+\lambda _{uZ}\bar {Y}\frac {\partial \bar {W}_{d}}{\partial \bar {Y}}=\big(1+h^{2}F_{Z}^{2}\big)\frac {\partial ^{3} \bar {U}_{d}}{\partial \bar {Y}^{3}}. \end{align}

In order to solve for $\bar {U}_{d}$ , (6.17) is rewritten in terms of $\zeta$ as

(6.18) \begin{align} \frac {\partial ^{3} \bar {U}_{d}}{\partial \zeta ^{3}}-\zeta \frac {\partial \bar {U}_{d}}{\partial \zeta }=-\frac {1}{\mathrm{i}\alpha }\frac {\partial \bar {W}_{d} }{\partial Z} +\frac {\lambda _{uZ}}{\mathrm{i}\alpha \lambda _{u}}\bar {W}_{d}+(\zeta -\zeta _{0})\frac {\lambda _{uZ}}{\mathrm{i}\alpha \lambda _{u}}\frac {\partial \bar {W}_{d}}{\partial \zeta }, \end{align}

on which the boundary and matching conditions are specified as

(6.19a,b) \begin{align} &\bar {U}_{d}=0,\quad \frac {\partial ^{2} \bar {U}_{d}}{\partial \zeta ^{2}} =\big(1+h^{2}F_{Z}^{2}\big)^{-1/3}(\mathrm{i}\lambda _{u}\alpha )^{-2/3}\mathrm{i}\alpha \bar {P}_{d}\quad \mathrm{at} \hspace {0.5em} \zeta =\zeta _{0}; \end{align}

(6.19c) \begin{align} &\bar {U}_{d}\rightarrow \lambda _{u}A_{1} \quad \mathrm{as} \hspace {0.5em}\zeta \rightarrow \infty . \end{align}

The solution to (6.18)–(6.19) is found as

(6.20) \begin{align} \bar {U}_{d}(\zeta )=&D_{0}\int ^{\zeta }_{\zeta _{0}}\mathrm{Ai}(\tau )\mathrm{d}\tau +\Big (D_{1}-\frac {3}{4}D_{3}\Big )\left [\mathrm{Ai}(\zeta )-\mathrm{Ai}(\zeta _{0})\right ] +\frac {D_{3}}{4}\left [\zeta \mathrm{Ai}'(\zeta )-\zeta _{0} \mathrm{Ai}'(\zeta _{0})\right ]\nonumber \\ &+\Big (D_{2}-\frac {3}{4}D_{4}\Big )\left [\mathrm{Gi}(\zeta )-\mathrm{Gi}(\zeta _{0})\right ]+\frac {D_{4}}{4}\left [\zeta \mathrm{Gi}'(\zeta )-\zeta _{0} \mathrm{Gi}'(\zeta _{0})\right ], \end{align}

where the expressions for $D_{\!j}\ (j=0,1,2,3,4)$ are given in Appendix B.

Imposing (6.19c ) on (6.20) leads to the pressure equation:

(6.21) \begin{align} \frac {\partial ^{2} \bar {P}_{d}}{\partial Z^{2}}-K_{1}\frac {\partial \bar {P}_{d} }{\partial Z}-\alpha ^{2}\bar {P}_{d}=K_{2}A_{1}, \end{align}

where

(6.22a) \begin{align} K_{1}=&\left [\mathrm{Ai}'(\zeta _{0})+\zeta _{0}\kappa (\zeta _{0})\right ]\frac {\zeta _{0}}{2\mathrm{Ai}(\zeta _{0})}\left [\frac { \lambda _{uZ}}{\lambda _{u}}+\frac {h^{2}F_{Z}F_{ZZ}}{1+h^{2}F_{Z}^{2}}\right ] +\frac {3}{2}\frac { \lambda _{uZ}}{\lambda _{u} }+\frac {h^{2}F_{Z}F_{ZZ}}{2\big(1+h^{2}F_{Z}^{2}\big)}, \end{align}

(6.22b,c) \begin{align} K_{2}=&\big(1+h^{2}F_{Z}^{2}\big)^{1/3}(\mathrm{i}\lambda _{u}\alpha )^{5/3}\mathrm{Ai}'(\zeta _{0})\kappa (\zeta _{0})^{-1},\quad \kappa (\zeta _{0})=\int ^{\infty }_{\zeta _{0}}\mathrm{Ai}(\tau )\mathrm{d}\tau , \end{align}

with $\mathrm{Ai}(\zeta )$ and $\mathrm{Gi}(\zeta )$ denoting the Airy and Scorer functions, respectively. A closed system describing the viscous lower-branch instability is formed by (6.11), (6.12) and (6.21). In contrast to bi-global stability analysis (Piot et al. Reference Piot, Casalis and Rist2008; Siconolfi et al. Reference Siconolfi, Camarri and Fransson2015; Szabó et al. Reference Szabó, Nagy, De Baets, Vanierschot and Paál2024; Bi et al. Reference Bi, Wang, Xu and Liu2023), the asymptotically reduced system is a one-dimensional eigenvalue problem in the spanwise direction. The coefficients in (6.21) depend solely on $\lambda _{u}$ , $hF$ and their spanwise derivatives, which provides a priori insight that the viscous lower-branch instability is controlled at leading order by the spanwise-dependent wall shear and the shape function of the roughness elements. The latter is an added new feature compared with Kátai & Wu (Reference Kátai and Wu2020), arising due to longer and higher roughness elements being considered.

6.5. The eigenvalue problem

Given that $F_{Z}(x,Z)$ and the streaky base flow $U(x,Y,Z)$ are spanwise-periodic, so are the coefficients $K_{1}$ and $K_{2}$ in (6.21). Thus, we adopt Floquet theory to solve the one-dimensional eigenvalue problem, expressing $p_{d}$ , $\bar {P}_{d}$ and $A_{1}$ as

(6.23) \begin{align} (p_{d},\bar {P}_{d},A_{1})=\exp (\mathrm{i}q\beta Z) \sum ^{\infty }_{n=-\infty }(\bar {p}_{n}(y),\bar {P}_{n},\bar {A}_{n})\exp (\mathrm{i}n\beta Z), \end{align}

where $q$ is the Floquet exponent, restricted to the interval [ $0$ , $1/2$ ] as explained by Kátai & Wu (Reference Kátai and Wu2020). Substitution of (6.23) into (6.11) yields

(6.24) \begin{align} \frac {\partial ^{2} \bar {p}_{n}}{\partial y^{2}} -\big [(n+q)^{2}\beta ^{2}+\alpha ^{2}\big ]\bar {p}_{n}=0. \end{align}

The solution satisfying the attenuating far-field condition for $\bar {p}_{n}$ is

(6.25) \begin{align} \bar {p}_{n}(y) =C_{n}\exp \Big [-\sqrt {(n+q)^{2}\beta ^{2}+\alpha ^{2}}y\Big ], \end{align}

where $C_{n}$ is a coefficient. Substituting (6.23) with (6.25) into (6.12) and writing the resulting equation in terms of $\bar {y}$ lead to

(6.26a) \begin{align} &\bar {P}_{n}=\sum ^{\infty }_{k=-\infty }a_{nk}C_{k}, \end{align}

(6.26b) \begin{align} & \bar {A}_{n}=\sum ^{\infty }_{k=-\infty }\alpha ^{-2}\Big [\sqrt {(k+q)^{2}\beta ^{2}+\alpha ^{2}}(a_{nk}+c_{nk})-\mathrm{i}(k+q)\beta b_{nk}\Big ]C_{k}, \end{align}

where we have performed Fourier series expansions:

(6.27) \begin{align} &\big(1,hF_{Z},h^{2}F_{Z}^{2}\big)\exp \Big [-\sqrt {(k+q)^{2}\beta ^{2}+\alpha ^{2}}hF\Big]\exp (\mathrm{i}k\beta Z)=\sum _{n=-\infty }^{\infty }(a_{nk},b_{nk},c_{nk})\mathrm{e}^{\mathrm{i}n \beta Z}. \end{align}

Similarly, we rewrite $K_{1}$ and $K_{2}$ into Fourier series:

(6.28a,b) \begin{align} K_{1}=\sum _{l=-\infty }^{\infty }R_{l}\mathrm{e}^{\mathrm{i}\beta l Z}, \quad K_{2}=\sum _{l=-\infty }^{\infty }Q_{l}\mathrm{e}^{\mathrm{i}\beta lZ}. \end{align}

Substituting (6.26) with (6.28) into (6.21) and equating the coefficients of the respective Fourier components, we arrive at the equations

(6.29) \begin{align} \big [\alpha ^{2}+\beta ^{2}(n+q)^{2}\big ]\bar {P}_{n} =-\sum _{j=-\infty }^{\infty }\mathrm{i}\beta (j+q)R_{n-j}\bar {P}_{\!j}-\sum _{j=-\infty }^{\infty }Q_{n-j}\bar {A}_{\!j}. \end{align}

Truncating the system (6.29) for $-K\leqslant n\leqslant K$ and eliminating $\bar {A}_{\!j}$ using (6.26), we write (6.29) in matrix form:

(6.30) \begin{align} \boldsymbol{M}\!\bar{\boldsymbol{P}}=0, \end{align}

where $\bar {\boldsymbol{P}}=(\bar {P}_{-K}$ , $\bar {P}_{-K+1}$ , $\dots$ , $\bar {P}_{0}$ , $\dots$ , $\bar {P}_{K})^{T}$ . The implicit dispersion relation is given by

(6.31) \begin{align} \mathcal{D}(\alpha )\equiv \det \big (\boldsymbol{M}(\alpha ,\omega ,\lambda _{u},q,hF)\big )=0. \end{align}

In the absence of roughness arrays, the dispersion relation (6.31) degenerates into

(6.32) \begin{align} D_{B}(\alpha )\equiv \alpha ^{2}\big [\alpha ^{2}+\beta ^{2}(n+q)^{2}\big ]^{1/2} +(\mathrm{i}\alpha \lambda _{B}(x))^{5/3}\mathrm{Ai}'(\zeta _{0})\kappa (\zeta _{0})^{-1}=0, \end{align}

where $\beta _{\textit{TS}}=\beta (n+q)$ is the spanwise wavenumber of TS modes and $\lambda _{B}(x)$ represents the wall shear of the Blasius flow.

We focus on the spatial instability for which the frequency $\omega$ is real, and $\alpha =\alpha _{r}+\mathrm{i} \alpha _{i}$ is complex with $\alpha _{r}$ being the streamwise wavenumber and $-\alpha _{i}$ the growth rate. The implicit dispersion relation (6.31) is solved using Muller’s iteration method (Muller Reference Muller1956). The numerical program used for these calculations is an in-house code that has been validated in previous studies (Kátai & Wu Reference Kátai and Wu2020; Zhu & Wu Reference Zhu and Wu2022; Zheng & Wu Reference Zheng and Wu2024). Furthermore, computations with different values of $K$ ensure that numerical results are independent of $K$ .

Given that sufficiently upstream of the roughness arrays the base flow approaches the state of the flat-plate case, the solutions for (6.31) and (6.32) are the same there. The latter may thus be used as the an initial guess to initiate the iteration for finding the former. A mode in the primary region influenced by roughness can be considered as developing continuously from a usual TS mode, and is computed as such. Detailed descriptions of the solution procedure and physical interpretation can be found in Kátai & Wu (Reference Kátai and Wu2020).

7.1. Local growth rates

In order to investigate the influence of the roughness elements on the TS modes, we first determine the most relevant ones on which to focus. Transition is related to the accumulated growth of disturbances, which is, in the linear regime, measured by the $N$ -factor, defined as

(7.1) \begin{align} N(\omega ,\beta _{\textit{TS}})=-\epsilon ^{-3}\int ^{x}_{x_{n}}\alpha _{i}(x)\,\mathrm{d}x, \end{align}

where $x_{n}$ denotes the neutral position, i.e. $\alpha _{i}(x_n)=0$ . In the incompressible regime, the most unstable and most amplified modes are 2-D. To select appropriate upstream modes to trace, we calculate the local growth rate and $N$ -factors of 2-D TS modes ( $\beta _{\textit{TS}}=0$ ) in the case of a flat plate. Figure 7 shows the streamwise variation of the $N$ -factor and the corresponding frequency of the most amplified mode, i.e. the mode with the largest $N$ at location $x$ . The $N$ -factor and the corresponding frequency $\omega$ increase and decrease with $x$ , respectively. The results indicate that modes with frequency $\omega \in [0.5, 3]$ acquire large $N$ -factors (while modes outside this frequency band would either acquire large $N$ -factors at too distant locations downstream or have considerably smaller $N$ -factors). Based on this observation, we investigate effects of roughness elements on the TS modes with frequencies ranging from 0.5 to 3.

Figure 7.

Streamwise variation of the $N$ -factor (solid line) in the flat-plate case for $\beta _{\textit{TS}}=0$ and the corresponding frequency (dashed line) of the most amplified mode.

Figure 8.

Local growth rates $(-\alpha _{i})$ versus the frequency $\omega$ of planar ( $\beta _{\textit{TS}}=0$ ) and oblique modes at various streamwise positions. (a–c) Symmetric (even) modes; (d–f) antisymmetric (odd) modes; (g–i) $\beta _{\textit{TS}}=0.1\pi$ ; (j–l) $\beta _{\textit{TS}}=0.2\pi$ . Curves: blue solid, Blasius; black solid, $h=0.11$ ; black dashed, $h=0.095$ . For both heights, the arrays have spanwise wavenumber $\beta =0.455\pi$ .

Figure 8 displays the local growth rates for roughness arrays with the same spanwise wavenumber $\beta =0.455\pi$ but different heights, $h=0.095$ and $0.11$ . The first, second and third columns are pertinent to streamwise locations $x = 0.9$ , 1.2 and 1.4, respectively. The TS modes with spanwise wavenumbers $\beta _{\textit{TS}}=0$ , $0.1\pi$ and $0.2\pi$ are computed. The non-zero values of $\beta _{\textit{TS}}$ cover a substantial range of obliqueness angles; e.g. at $x=0.9$ , $\beta _{\textit{TS}}=0.2\pi$ covers obliqueness angles between $35^\circ$ (minimum) and $74^\circ$ (maximum), which are attained by modes with $\omega =0.5$ and $3$ , respectively. The first and second rows in figure 8 show the results for the 2-D modes ( $\beta _{\textit{TS}}=0$ ), which are unique in that their ensuing development may be either symmetric or antisymmetric; they will, for brevity, be referred to as ‘even’ and ‘odd’ modes, respectively. Clearly, low-frequency modes are stabilised by the roughness elements, while higher-frequency modes are destabilised, consistent with the finding in Kátai & Wu (Reference Kátai and Wu2020). Even and odd modes follow similar trends across the frequency range with odd modes having slightly higher growth rates (cf. figures 8 $a$ – $c$ and 8 $d$ – $f$ ). The ‘critical’ frequency for the switching from stabilising to destabilising effects depends on the streamwise location and weakly on the value of $h$ . The above conclusions hold for weakly oblique TS modes, as is indicated in the third row of figure 8 for $\beta _{\textit{TS}}=0.1\pi$ . For stronger 3-D modes, the stabilising effect on low-frequency modes becomes weaker, whereas the destabilising effect extends to lower frequencies and becomes more significant as the fourth row in figure 8 shows. A comparison of the results represented by the black solid and dashed lines indicates that roughness arrays of greater height have stronger effects on linear stability.

Figure 9 shows the local growth rates for roughness elements with the same height $h=0.095$ but different spanwise wavenumbers, $\beta =0.455\pi$ , $0.341\pi$ and $0.227\pi$ . For 2-D modes, roughness elements with larger spanwise wavenumbers have a stronger stabilising/destabilising effect on the linear stability in the lower/higher frequency band. The switchover between the stabilising and destabilising effects exhibits again a weak dependence on $\beta$ . This trend remains for weakly oblique modes with $\beta _{\textit{TS}}=0.1\pi$ , as shown in figure 9(g–i). However, the arrays with $\beta =0.227\pi$ appear somewhat out of line as the switchover occurs at a lower frequency than that for $\beta =0.341\pi$ and $0.455\pi$ . For stronger oblique modes with $\beta _{\textit{TS}}=0.2\pi$ , roughness elements slightly inhibit lower-frequency modes, but appreciably destabilise higher-frequency ones, and this effect is stronger for larger $\beta$ . Overall, roughness arrays having the same height but different spanwise wavenumbers play a qualitatively similar role, but the quantitative extent of their influence depends on $\beta _{\textit{TS}}$ and $\beta$ , the spanwise wavenumbers of the modes and the roughness arrays.

Figure 9.

Local growth rates $(-\alpha _{i})$ versus the frequency $\omega$ of planar ( $\beta _{\textit{TS}}=0$ ) and oblique modes at various streamwise positions. (a–c) Symmetric (even) modes; (d–f) antisymmetric (odd) modes; (g–i) $\beta _{\textit{TS}}=0.1\pi$ ; (j–l) $\beta _{\textit{TS}}=0.2\pi$ . Curves: blue solid, Blasius; black solid, $\beta =0.455\pi$ ; black dashed, $\beta =0.341\pi$ ; black dotted, $\beta =0.227\pi$ . The roughness height $h=0.095$ .

The viscous instability is controlled by the wall shear and the geometry of the roughness elements, which are determined by two parameters, height $h$ and spanwise wavenumber $\beta$ . The studies mentioned above monitored effects of $h$ for a fixed $\beta$ (figure 8) and of $\beta$ for a fixed $h$ (figure 9). A close inspection of the mathematical description suggests that the roughness-induced spanwise pressure gradient and slip velocity driving the spanwise motion in the boundary layer are proportional to $h$ (if $h$ is assumed small), and so the streaks, forced through the continuity equation, must be proportional to $h\beta$ . It follows that the wall shear and the geometric effect (represented by $hF_{Z}$ ) both scale as $h\beta$ . It is thus tempting to consider the product $h\beta$ as the relevant parameter for characterising the effect of roughness arrays on the instability. However, it should be noted that this proposition is underpinned by taking $h\ll 1$ and the observation that $\beta$ appears on its own in the stability equations. In what follows, we address the question: to what degree is the proposition valid?

We performed calculations for roughness arrays with different height $h$ and spanwise wavenumber $\beta$ : $(h,\beta )=(0.055,0.909\pi )$ , $(0.11,0.455\pi )$ and $(0.22,0.227\pi )$ , for all of which the product $h\beta =0.05\pi$ . The local growth rates are shown in figure 10. For planar symmetric modes with low frequencies, the growth rates for the three cases collapse perfectly into a single curve (figure 10 a–c). However, when the frequency exceeds a critical value, the results differ, and roughness arrays with smaller $\beta$ (i.e. large spanwise spacing) have stronger destabilising effects. Surprisingly, the growth rates of odd modes completely overlap for all frequencies (figure 10 d–f), indicating that $h\beta$ is the sole relevant parameter in this case. Weakly oblique modes ( $\beta _{\textit{TS}}=0.1\pi$ ) exhibit a similar dependence on $h\beta$ (figure 10 g–i) to the planar even modes. For stronger oblique modes with $\beta _{\textit{TS}}=0.2\pi$ (figure 10 j–l), the growth rate is primarily dependent on the product $h\beta$ , exhibiting appreciable but insignificant difference from case to case only at the high-frequency end.

Figure 10.

Local growth rates $(-\alpha _{i})$ versus the frequency $\omega$ of planar $(\beta _{\textit{TS}}=0)$ and oblique modes at various streamwise positions. (a–c) Symmetric (even) modes; (d–f) antisymmetric (odd) modes; (g–i) $\beta _{\textit{TS}}=0.1\pi$ ; (j–l) $\beta _{\textit{TS}}=0.2\pi$ . Curves: blue solid, Blasius; black solid line with triangle, $(h,\beta )=(0.055,0.909\pi )$ ; black dashed line with square, $(h,\beta )=(0.11,0.455\pi )$ ; black dotted line with circle, $(h,\beta )=(0.22,0.227\pi )$ .

7.2. Spatial development and the change of N-factor

In order to characterise further the effect of roughness elements on instability, we now trace the development of upstream TS modes into the streaky region. Six modes are selected with frequencies $\omega =0.5$ , $1$ , $1.5$ , $2$ , $2.5$ and $3$ . Their growth rates as a function of $x$ are displayed in figure 11. The results show once again, but from a different perspective, that 2-D symmetric modes ( $\omega =0.5$ and $1$ , shown in figure 11 a,b) are stabilised on entering the streaky region. The stabilising effect operates in the region before or just slightly after the mode acquires its maximum growth rate. However, for $\omega =1.5$ (figure 11 c), the mode is first inhibited and then enhanced by the roughness elements, which play a destabilising role overall. For larger frequencies ( $\omega =2$ , $2.5$ and $3$ , shown in figure 11 d–f), the modes are hardly affected by the roughness arrays before reaching their respective peak growth rates, downstream of which they are destabilised. Similar trends are observed for the 2-D antisymmetric and 3-D modes, the results for which are shown in the supplementary material available at https://doi.org/10.1017/jfm.2026.11174.

Figure 11.

(a–f) Local growth rates $(-\alpha _{i})$ versus $x$ for symmetric planar modes ( $\beta _{\textit{TS}}=0$ ). Curves: blue solid, Blasius; black solid, $h=0.11$ ; black dashed, $h=0.095$ . The spanwise wavenumber of roughness arrays $\beta =0.455\pi$ .

Figure 12.

The $\Delta N$ -factor versus $\omega$ for representative instability modes in the case of the Gaussian-like and sinusoidal roughness. (a–c) $\beta _{\textit{TS}}=0$ ; (d–f) $\beta _{\textit{TS}}=0.1\pi$ ; (g–i) $\beta _{\textit{TS}}=0.2\pi$ . (a,d,g) $h=0.11$ and $0.095$ , $\beta =0.445\pi$ ; (b,e,h) $\beta =0.445\pi$ , $0.341\pi$ and $0.227\pi$ , $h=0.095$ ; (c, f,i) $(h,\beta )=(0.055, 0.909\pi )$ , $(0.11, 0.455\pi )$ and $(0.22, 0.227\pi )$ with $h\beta =0.157$ . The parameter: $\epsilon =0.2036$ .

To investigate the overall effects of roughness elements on the modes, we monitor the change of $N$ -factor, defined as

(7.2) \begin{equation} \Delta N(\omega ,\beta _{\textit{TS}})=\epsilon ^{-3}\int ^{x_{1}}_{x_{0}}[-\alpha _{i}-(-\alpha _{i,0})]\mathrm{d}x, \end{equation}

where $x_{0}$ and $x_{1}$ represent upstream and downstream locations sufficiently far from the roughness elements, respectively. The local growth rates in the streaky and flat-plate cases are denoted by $-\alpha _{i}$ and $-\alpha _{i,0}$ , respectively.

For the Gaussian-like and sinusoidal arrays (Roughness 1), we calculate the variation of $\Delta N$ with $\omega$ for representative modes with $\beta _{\textit{TS}}=0$ , $0.1\pi$ and $0.2\pi$ ; here, the roughness-array parameters correspond to those in figures 8–10. The results are displayed in figure 12. The first column is pertinent to arrays with the same spanwise wavenumber but different heights. The integration range is set to be $(x_{0}, x_{1})=(0.01, 1.48)$ , where $x_{1}=1.48$ is very close to the separation point. The second column presents the $\Delta N$ -factor for arrays having the same height but different spanwise wavenumbers, with the integration range being $(x_{0}, x_{1})=(0.01, 1.94)$ . The third column displays the $\Delta N$ -factor for roughness elements with different $\beta$ and $h$ , but the same product $h\beta =0.05\pi$ ; the integration range is set to be $(x_{0}, x_{1})=(0.01, 1.48)$ . The results shown in the three figures reaffirm that the roughness elements exert an overall stabilising effect on 2-D and 3-D low-frequency modes, while destabilising most of the other modes. Higher roughness arrays produce stronger effects. However, the effect of the spanwise wavenumber $\beta$ does not follow a monotonic trend.

7.3. Comparison with the experimental results of Downs & Fransson (Reference Downs and Fransson2014)

7.3.1. Comparison of $\Delta N$ values

Downs & Fransson (Reference Downs and Fransson2014) showed that the leading and trailing edges of the roughness elements induce an enhanced amplification, the reason for which is not entirely clear. However, over the central part of the roughness elements, 420 mm $\lt x^{+}\lt$ 540 mm and 470 mm $\lt x^{+}\lt$ 720 mm for the two configurations (where $x^{+} = xL$ represents the dimensional streamwise coordinate), the disturbances grow at reduced rates compared with the flat-plate case.

Figures 13–15 present a comparison of experimentally measured and theoretically predicted $\Delta N$ -factors in the central region. Computations are performed for three types of roughness configurations: Gaussian-like and sinusoidal (Roughness 1), step-like and hemispherical-cap (Roughness 2) and step-like and sinusoidal (Roughness 3), of which the second models the experiments closely, and is our main focus. For each configuration, the streamwise length and spanwise spacing may take two values. The experimental data are extracted from figures 10, 13 and 16 in Downs & Fransson (Reference Downs and Fransson2014). Overall, the theoretical predictions agree with the measurements qualitatively in terms of general trends, while appreciable quantitative discrepancies are noted. Theoretical results in figure 13 indicate that the $\Delta N$ -factor decreases with the spanwise spacing of the roughness array. Furthermore, figure 14 demonstrates that for longer roughness arrays, both theoretical and experimental results feature an extended region of negative $\Delta N$ -factor, indicating a more substantial stabilising effect on TS waves. Figure 15 presents the results for higher-frequency modes, which are similar to those in figure 13( $a$ ).

Figure 13.

The $\Delta N$ -factor versus the dimensional coordinate $x^{+}$ for planar symmetric modes ( $\beta _{\textit{TS}}=0$ ). The spanwise spacing of the arrays: ( $a$ ) $\varLambda ^{+}=15$ mm; ( $b$ ) $\varLambda ^{+}=30$ mm. Curves: black solid, step-like and hemispherical-cap roughness ( $l_{x}=0.28$ and $d_{x}=0.1$ ); black dashed, step-like and sinusoidal roughness; black dotted, Gaussian-like and sinusoidal roughness ( $d=0.2$ ); circles, experimental results ( $l_{x}^{+}=250$ mm and $h^{+}=1$ mm) of Downs & Fransson (Reference Downs and Fransson2014). Parameters used in our calculations: $h=0.12$ ( $h^{+}=0.43$ mm), $\varLambda =4.4/8.8$ and $\epsilon =0.2036$ .

Figure 14.

The $\Delta N$ -factor versus $x^{+}$ for planar symmetric modes ( $\beta _{\textit{TS}}=0$ ). Curves: black solid, step-like and hemispherical-cap roughness ( $l_{x}=0.4$ and $d_{x}=0.1$ ); black dashed, step-like and sinusoidal roughness; black dotted, Gaussian-like and sinusoidal roughness ( $d=0.3$ ); circles, experimental results ( $l_{x}^{+}=480$ mm and $h^{+}=1$ mm) of Downs & Fransson (Reference Downs and Fransson2014). Parameters used in our calculations: $h=0.103$ ( $h^{+}=0.43$ mm), $\varLambda =3.6$ and $\epsilon =0.1976$ .

Figure 15.

The $\Delta N$ -factor versus $x^{+}$ for planar symmetric modes ( $\beta _{\textit{TS}}=0$ ) for different frequencies: ( $a$ ) $F_{r}=110$ ; ( $b$ ) $F_{r}=120$ . Curves: black solid, step-like and hemispherical-cap roughness ( $l_{x}=0.28$ and $d_{x}=0.1$ ); black dashed, step-like and sinusoidal roughness; black dotted, Gaussian-like and sinusoidal roughness ( $d=0.2$ ); circles, experimental results ( $l_{x}^{+}=250$ mm and $h^{+}=1$ mm) of Downs & Fransson (Reference Downs and Fransson2014). Parameters used in our calculations: $h=0.12$ ( $h^{+}=0.43$ mm), $\varLambda =4.4$ and $\epsilon =0.2036$ .

The comparison between the three configurations underscores the influence of both the streamwise and spanwise patterns of roughness arrays on the $\Delta N$ -factor. The second configuration turns out to be consistently most effective. Comparing the step-like and sinusoidal arrays with the Gaussian-like and sinusoidal arrays (figures 13 $a$ and 15), we note that the difference of the $\Delta N$ -factor becomes larger with increasing frequency. This behaviour is attributed to the fact that the maximum growth rate of higher-frequency modes is located closer to the leading edge of the flat plate, where the step-like distribution has a greater influence on the base flow than the Gaussian-type distribution.

Note that the heights used in our calculations are a little less than half of the experimental ones, but the theoretically calculated $\Delta N$ values for some cases are higher than the experimental measurements in some cases. There are several possible reasons for the quantitative discrepancy. The first may be attributed to the presence of a weak but appreciable adverse pressure gradient in the experiments. The effect of the latter may be non-negligible since boundary-layer instability is known to be sensitive to it. The second reason may be the inherent non-parallelism of the streaky boundary layer, and different ways of defining the $\Delta N$ -factor in the theory and experiment. In our formulation, the non-parallelism of the streaky flow is neglected and $\Delta N$ is calculated using the local modal growth rate, which is independent of the wall-normal and spanwise positions. In the experiments, the $N$ -factor was first computed using the local amplitude of the disturbance. The latter may be measured differently, and in Downs & Fransson (Reference Downs and Fransson2014) it was taken to be the root mean square of the fluctuation integrated across the boundary layer and one spanwise period. The difference between the $N$ -factors with and without roughness arrays was taken to be $\Delta N$ . Clearly, the $N$ - and $\Delta N$ -factors were affected by the non-parallelism. Furthermore, the fluctuations measured likely consist of non-modal components. We note that the asymptotic theory can be extended to account for effects of both pressure gradient and non-parallelism. Finally, the discrepancy could be due to the rather moderate Reynolds number considered, and/or the limitation of the lower-branch scaling regime adopted. Figure 16( $c$ ) of Downs & Fransson (Reference Downs and Fransson2014) shows that the amplitude is already undergoing attenuation near $x^{+}\approx 550$ mm, indicating the proximity to the upper branch of the neutral curve. The instability may have to be described by a five-deck asymptotic structure (Smith & Bodonyi Reference Smith and Bodonyi1982; Bodonyi, Smith & Gajjar Reference Bodonyi, Smith and Gajjar1983). This prospect is left for future investigation.

Figure 16.

Experimental results of Downs & Fransson (Reference Downs and Fransson2014): (a) contours of the TS-wave amplitude, $A_{T}$ , at the streamwise locations $x=1.17$ ( $x^{+}=500$ mm), where $U_{e}$ is the slip velocity. The reader is referred to Downs & Fransson (Reference Downs and Fransson2014) for the definitions of these quantities. The spanwise periods of the first, second and third rows in each panel are $\varLambda ^{+}=15$ , 20 and 30 mm, respectively. The spacing $\varLambda =$ 4.2, 5.6 and 8.4, respectively. The curves in (b) show the profile of the TS-wave amplitude along $Z/\varLambda =0$ (thin solid line), $Z/\varLambda =1/4$ (dot-dashed line) and $Z/\varLambda =1/2$ (dashed line ). The thick solid line represents the eigenfunction of the undistorted TS mode in the absence of roughness, which is included as a reference. The TS-wave frequency $\omega \approx 1.4$ .

7.3.2. Comparison of eigenfunctions

Figure 16 shows the key measurements of Downs & Fransson (Reference Downs and Fransson2014): contours of TS-wave amplitude (figure 16 a) and TS-wave profiles (figure 16 b) for a mode with $\omega \approx 1.4$ at $x=1.17$ ( $x^{+}=500$ mm). The spanwise spacings of the roughness arrays are $\varLambda ^{+}=15$ , 20 and 30 mm. The relevant characteristics are commented on when we present our theoretical results. To facilitate comparisons with the experiment, the leading-order lower- and main-deck solutions for the streamwise velocity are combined to construct a composite solution, valid in both the main and lower decks:

(7.3) \begin{equation} U_{c}=\bar {U}_{d}(\bar {Y})+U_{d}(Y)-\lambda _{u}A_{1}, \end{equation}

where $\bar {U}_{d}$ and $U_{d}$ are given by (6.8 $a$ ) and (6.20), respectively, while the third term is the common part of $\bar {U}_{d}$ and $U_{d}$ . The solution is written in terms of $\hat {\eta }$ with $\bar Y$ and $Y$ being converted to $\hat {\eta }$ , where $\hat {\eta }=(Y+\epsilon ^{-1}hF)/\sqrt {x}$ . The value of $\epsilon$ used in the conversion is taken to be $0.2036$ , which corresponds to $\textit{Re}=3.4\times 10^{5}$ .

Figure 17.

Eigenfunction of the symmetric mode with $\omega =1.4$ and $\beta _{\textit{TS}}=0$ for the step-like and hemispherical-cap roughness ( $l_{x}=0.28$ , $d_{x}=0.1$ ). (a,c,e,g) Contours of $|U_{c}|$ in the $\hat {\eta }{-}Z$ plane; (b,d, f,h) profiles of $|U_{c}|$ at $Z/\varLambda =0$ (solid line), $1/4$ (dotted line) and $1/2$ (dashed line). The blue solid line represents the eigenfunction of the undistorted TS mode in the absence of roughness, which is included as a reference. Parameters used in our calculations: $h=0.12$ ( $h^{+}=0.43$ mm) and $\varLambda =4.4$ .

Figure 18.

Eigenfunction of the symmetric mode with $\omega =1.4$ and $\beta _{\textit{TS}}=0$ for the step-like and hemispherical-cap roughness ( $l_{x}=0.28$ , $d_{x}=0.1$ ). (a,c,e,g) Contours of $|U_{c}|$ in the $\hat {\eta }{-}Z$ plane; (b,d, f,h) profiles of $|U_{c}|$ at $Z/\varLambda =0$ (solid line), $1/4$ (dotted line) and $1/2$ (dashed line). The blue solid line represents the eigenfunction of the undistorted TS mode in the absence of roughness, which is included as a reference. Parameters used in our calculations: $h=0.12$ ( $h^{+}=0.43$ mm) and $\varLambda =6$ .

Figure 19.

Eigenfunction of the symmetric mode with $\omega =1.4$ and $\beta _{\textit{TS}}=0$ for the step-like and hemispherical-cap roughness ( $l_{x}=0.28$ , $d_{x}=0.1$ ). (a,c,e,g) Contours of $|U_{c}|$ in the $\hat {\eta }{-}Z$ plane; (b,d, f,h) profiles of $|U_{c}|$ at $Z/\varLambda =0$ (solid line), $1/4$ (dotted line) and $1/2$ (dashed line). The blue solid line represents the eigenfunction of the undistorted TS mode in the absence of roughness, which is included as a reference. Parameters used in our calculations: $h=0.12$ ( $h^{+}=0.43$ mm) and $\varLambda =8.8$ .

Figures 17–19 show the distribution of the eigenfunction $U_{c}$ of a symmetric mode with $\omega =1.4$ and $\beta _{\textit{TS}}=0$ in the case of step-like and hemispherical-cap roughness ( $l_{x}=0.28$ , $d_{x}=0.1$ and $h=0.12$ ) with spanwise spacing $\varLambda =4.4$ (figure 17), 6 (figure 18) and 8.8 (figure 19). The results include contours in the $\hat {\eta }{-}Z$ plane (left-hand column) and typical wall-normal profiles (right-hand column) at $x=0.8$ , 1, 1.1 and 1.17.

Upstream of the roughness peak, there is a single maximum per spanwise wavelength with the eigenfunction concentrating at $Z/\varLambda =0$ (panel ( $a$ ) of figures 17–19). The modal shape undergoes significant distortion downstream, the process of which depends on the form of the spanwise spacing and height of the roughness arrays. Distortion occurs from $x=1$ to 1.17. The central maximum splits into two off-centre maxima, followed by the emergence of a new central maximum as shown in panels ( $c$ – $g$ ) of figures 17–19. The overall feature, the presence of three maxima, resembles the experimental measurements of Downs & Fransson (Reference Downs and Fransson2014) (figure 16). However, differences in details can be noted. The disturbance predicted exhibits high spanwise gradients between the adjacent maxima, and becomes highly localised at $Z/\varLambda =0$ , where eigenfunction contours feature a cusp (see panel ( $g$ ) of figures 17 and 18), which is due to the like focusing behaviour of the streaks as shown in figure 3. Comparison of figures 17–19 suggests that the sharpness of the cusp decreases with the spanwise spacing, and in particular a cusp does not appear for the widest spacing (figure 19 $g$ ).

Now we turn to profiles of modal disturbances. The predicted profiles at various spanwise positions bear broadly similar characteristics to those observed in the experiments (cf. profiles shown in figures 17 h and 18 h and figure 16). The profiles differ appreciably from that of the usual undistorted TS mode. The profile at $Z/\varLambda =0$ stands out as it has a greater magnitude, and features a prominent peak at larger distances to the wall ( $\hat {\eta }\approx 3$ ) compared with those at other spanwise locations.

While the salient overall features of the experimental measurements of Downs & Fransson (Reference Downs and Fransson2014) are captured qualitatively by the theory, differences in details can be noted. The disturbance predicted exhibits significant spanwise focusing, leading to high spanwise gradients between the adjacent maxima, which is not present in the experimental data. The difference may be due to the limitation of the present asymptotic approach of boundary-layer type, which does not fully account for the spanwise diffusion (ellipticity). This approximation is no longer valid when the focusing-induced spanwise gradients become strong and a new asymptotic theory that accounts for the full effect of the shorter spanwise scale must be developed. This issue is left as a topic for future investigation.