Nonlinear distortions and short-wavelength secondary instability directly induced by distributed roughness in three-dimensional boundary layers

Abstract

Surface roughness of fairly small (micron-sized) height is known to influence significantly three-dimensional boundary-layer transition. In this paper, we investigate this sensitive effect from the viewpoint that roughness alters the base flow thereby inducing new instabilities. We consider distributed roughness in the form of a wavy wall with its height being taken to be of $\mathit{O} (R^{-1/3 } \delta ^{\ast })$, where the Reynolds number $R$ is defined using the local boundary-layer thickness $\delta ^{\ast }$. Despite having a height much smaller than $\delta ^{\ast }$, the roughness is high enough to induce nonlinear responses. The roughness-distorted boundary-layer flow is characterised by a wall layer (WL) – a thin layer adjacent to the surface – the main layer and a critical layer (CL) – the vicinity of a special position at which a singularity of the Rayleigh equation occurs. The widths of both the WL and CL are of $\mathit{O} (R^{-1/3} \delta ^{\ast })$. Surface roughness alters the base flow significantly, leading to $\mathit{O} (1)$ vorticity distortions in these layers. We show for the first time that the nonlinearly distorted flows in these layers support small-scale local instabilities due to the roughness-induced $\mathit{O} (1)$ vorticities. Two types of modes, CL and WL modes, are identified. The CL modes have short wavelengths and high frequencies, with the spatial and temporal instabilities being governed by essentially the same equation. Thus, we focus on the former, which can be formulated as a linear generalised eigenvalue problem. The WL modes have short wavelengths but $\mathit{O} (1)$ frequencies. The temporal WL mode is governed by a linear eigenvalue problem similar to that for the CL modes, while the spatial WL mode is described by a nonlinear eigenvalue problem. The onset of these small-scale fluctuations could form a crucial step in the transition to turbulence.

1. Introduction

It is well acknowledged that surface roughness can significantly influence boundary-layer transition through different mechanisms, depending on its form, location and scale. Experimental observations of surface roughness effects are reviewed first before summarising related numerical and theoretical investigations.

Early experiments focused on identifying the critical height below which surface roughness has no substantial effect on transition (Fage 1943; Tani 1961). Aimed at obtaining empirical formulae, those experiments offered little insight into the mechanism involved. Klebanoff & Tidstrom (1972) performed first detailed experiments to shed light on the effect of two-dimensional isolated surface roughness on boundary-layer transition. The roughness heights were set as $0.7$ – $0.8 \delta ^{\ast }$ with $\delta ^{\ast }$ being the local boundary-layer thickness. They attributed the mechanism to the flow distortion in the wake, which features an inflectional profile and is thus susceptible to strong inviscid instability.

Two-dimensional boundary layers are known to be susceptible to Tollmien–Schlichting (T–S) instability, which is of viscous nature. Corke, Bar-Sever & Morkovin (1986) studied, by using hot-wire measurements and smoke-wire visualisation, the effect of distributed surface roughness with heights approximately $0.5 \delta ^{\ast }$ . The T–S waves grew more rapidly in the boundary layer over the rough wall than on the smooth wall. However, no inflectional point appears in boundary-layer profiles over the rough wall, which possessed a similar shape factor as that for the smooth wall. This suggested that the enhanced growth was not due to inviscid instability, but rather likely caused by the distorted base flow, whereas the excitation of T–S waves by the roughness interacting with free stream turbulence led to increased amplitude. Ma’mun, Asai & Inasawa (2014) investigated effects of both two- and three-dimensional sinusoidal surface corrugations on the growth of T–S waves in a zero-pressure-gradient boundary layer. The wavelengths of surface corrugations are of the same order as that of T–S modes, while the roughness heights are approximately $0.1 \delta ^{\ast }$ . They found the two-dimensional surface corrugation significantly enhanced the growth of T–S waves, while the three-dimensional corrugation had little effect despite its greater height of $0.2 \delta ^{\ast }$ . However, the three-dimensional corrugation induced flow distortions, which interacted with the planar T–S waves, generating pairs of oblique waves. The latter may undergo rapid growth, facilitated by fundamental resonance, leading to transition of Klebanoff-type, which is characterised by aligned peak-valley structures.

Three-dimensional boundary layers may support four types of instabilities including attachment-line instability, centrifugal (Görtler) instability, streamwise (T–S) instability and cross-flow instability (Saric, Reed & White 2003); the last has attracted much attention since it plays a leading role in destabilisation and transition of boundary layers over swept wings. Cross-flow instability comprises two modes: stationary vortices and travelling vortices. Stationary vortices dominate transition in low free stream disturbance environments such as flight conditions, while travelling vortices become important in the presence of high-level free stream turbulence (Deyhle & Bippes 1996). It is well recognised that surface roughness played a crucial role in the excitation of stationary cross-flow vortices. Reibert et al. (1996) studied the impact of various surface roughness forms on a three-dimensional boundary layer over a swept wing. They found that submicron-sized irregular surface roughness located near the attachment line of the wing could excite irregular stationary cross-flow vortices, while spanwise periodic roughness arrays excited regular stationary vortices. The roughness spacings had considerable influence on the excitation and growth of these vortices. Further experimental investigations by Carrillo, Reibert & Saric (1997) and Saric, Carrillo & Reibert (1998) confirmed the dramatic sensitivity of stationary cross-flow vortices to surface roughness. Radeztsky, Reibert & Saric (1999) examined micron-sized roughness in both distributed (painted and polished surfaces) and isolated forms, showing that stationary cross-flow instability is strongly altered by surface roughness while barely affected by acoustic disturbances.

Theoretically, four possible mechanisms have been identified or proposed by which surface roughness affects boundary-layer transition: receptivity, local scattering, modal interactions and alteration of the base flow thereby changing flow stability. Receptivity refers to the process by which instability modes are excited by external disturbances. For T–S waves in two- and three-dimensional boundary layers, the excitation requires a scale conversion/tuning process, which includes roughness-acoustic (Ruban 1984; Goldstein 1985) and roughness-vortical disturbances interactions (Duck, Ruban & Zhikharev 1996; Wu 2001a , 2001b ). For the excitation of travelling cross-flow vortices, it is generally considered that the roughness-induced perturbation also needs to interact with unsteady free-stream disturbances to generate suitable forcing. In contrast, stationary cross-flow vortices can be excited alone by surface roughness with comparable length scales. A mathematical description of this receptivity was formulated by Crouch (1993) and Choudhari (1994) using the finite-Reynolds-number theory based on the Orr–Sommerfeld (O–S) equation. However, non-parallelism cannot be incorporated in this approach, whose applicability is restricted to roughness with sufficiently small heights. For receptivity involving roughness with large height to induce nonlinear flow distortion, a high-Reynolds-number asymptotic approach needs to be adopted as was done by Choudhari & Duck (1996). Using this approach, Butler & Wu (2018) investigated the receptivity of three-dimensional boundary layer to chordwise-localised, spanwise-periodic surface roughness located near the leading edge, where non-parallelism plays a leading-order role. They found that roughness with sharper edges could generate stronger vortices and confirmed ‘superlinearity’ behaviour found in numerical simulations (Kurz & Kloker 2014).

A localised roughness element (hump) that has a length scale comparable to the wavelength of the instability and is positioned in the main unstable regions can influence transition by scattering the oncoming instability waves (Wu & Hogg 2006). The local stability analysis fails in this case because of the strong non-parallelism associated with the roughness-induced mean-flow distortion. The scattering problems of a small-amplitude T–S wave by a shallow and large roughness that induces linear and nonlinear flow distortions, were studied by Wu & Hogg (2006) and Wu & Dong (2016), respectively, using the triple-deck theory (Smith 1979). They introduced, for the first time, the concept of a transmission coefficient, which is defined as the ratio of the transmitted wave amplitude to that of the incident wave. Xu et al. (2016) investigated the effect of scattering by surface hump induced localised distortion using direct numerical simulation (DNS) and confirmed the prediction by the asymptotic theory. A finite-Reynolds-number formulation for the local scattering problem was developed by Huang & Wu (2017). The theoretical predictions agreed well with DNS results and experimental data.

A possible explanation for the extreme sensitivity of three-dimensional boundary-layer transition to micron-sized surface roughness, as reported in Reibert et al. (1996), Carrillo et al. (1997), Saric et al. (1998) and Radeztsky et al. (1999), is the modal interactions between cross-flow instability modes and ‘roughness modes’, which are distortions induced by distributed roughness rather than usual eigenmodes (He, Butler & Wu 2019). Resonant interactions, including a generalised Bragg scattering and a triadic resonance, could take place among an appropriate roughness mode and one or more cross-flow vortices, enhancing the amplification of the latter. This mechanism was demonstrated for the Falkner–Skan–Cooke (FSC) flow (He et al. 2019) and the boundary layer over the NLF(2)-0415 swept wing (Xu & Wu 2022).

In addition to the mechanisms above, surface roughness located in the main unstable region, which is a pervasive situation, can distort the base flow thereby changing stability characteristics, or even inducing new instabilities. Theoretical studies for this mechanism include calculating the roughness-distorted flow and performing stability analysis for the new base state. The impact of a localised roughness on a two-dimensional boundary layer has been examined by Nayfeh, Ragab & Al-Maaitah (1988), Cebeci & Egan (1989) and Gao, Park & Park (2011). In general, the presence of a surface hump was found to promote the growth of T–S waves. Wie & Malik (1998) investigated the effect of surface waviness, which is widely used as a model for distributed surface roughness, on two-dimensional boundary-layer transition. The distorted mean flow was calculated using the interactive boundary layer (IBL) method. Then the method of linear parabolised stability equations (PSE) was adopted to predict T–S wave amplification. The effect is characterised by the $N$ -factor increment $\Delta N$ , which was found to scale as $ h^{\ast 2} / \lambda ^{\ast }$ , where $h^{\ast }$ and $\lambda ^{\ast }$ are the roughness height and wavelength, respectively. Masad (1996) considered the effect of streamwise localised and spanwise uniform surface waviness on three-dimensional boundary-layer transition. The roughness-distorted mean flow was also calculated using the IBL method, and the transition location, predicted using linear stability analysis coupled with the $\textrm{e}^{N}$ -method, was found to be shifted upstream. Floryan (1997) studied the stability of wall-bounded shear layers subjected to distributed wall suctions, which were introduced to model distributed surface roughness with short wavelengths comparable to the shear-layer thickness. A linear approximation was made to calculate the flow modification, and Floquet theory was adopted to study the stability of the modified flow. For all cases considered (the Poiseuille flow, Couette flow and Blasius boundary layer), a new type of instability in the form of streamwise vortices is induced.

Recently, Hall (2021) identified a streamwise vortex instability mechanism in boundary layers caused by a two-dimensional wavy wall with the wavelength comparable to the boundary-layer thickness. The streamwise vortex instability ‘mode’, which has the long length scale of the underlying boundary layer, interacts with the wavy-wall induced flow in the wall layer (WL) to generate a short-scale motion. The latter interacts in turn with the waviness-induced flow to produce a spanwise velocity, which serves as an effective slip condition on the streamwise vortex residing in the main boundary layer. This is essentially a steady streaming effect, and renders the vortex to amplify over the same streamwise length scale as that of the base flow when the waviness height is $\mathit{O} (R^{-2/3} / \sqrt {\ln R})$ , where $R$ is Reynolds number based on the boundary-layer thickness.

Stationary cross-flow vortices in three-dimensional boundary layers resemble nonlinear distortions induced by distributed roughness, and like the latter do not lead to transition directly. Experiments show that under the influence of nonlinearity these vortices saturate well before transition happens (Reibert et al. 1996; Bippes 1999). Due to the presence of finite wall-normal vorticity and/or substantially altered spanwise vorticity, the evolving nonlinearly distorted flow is susceptible to secondary instability, the rapid amplification of which causes transition (Kohama, Saric & Hoos 1991; White & Saric 2005; Serpieri & Kotsonis 2016).

Using the FSC flow as a model, Malik, Li & Chang (1994) performed a secondary instability analysis for saturated stationary cross-flow vortices, which were calculated using nonlinear PSE, and demonstrated the occurrence of high-frequency secondary instability. Högberg & Henningson (1998) conducted DNS to investigate the development of stationary cross-flow vortices, the secondary instability of which was assessed by linear eigenvalue calculations. They identified a high-frequency secondary mode and a low-frequency interaction mode. The latter was first discovered by Fischer & Dallmann (1991) as a ‘secondary instability’, but may more appropriately be treated as an interaction between stationary and travelling vortices. Högberg & Henningson (1998) suggested that the high- and low-frequency secondary modes may cause transition in low- and high-turbulence environments, respectively. Instead of FSC flow, Malik et al. (1999) analysed the secondary instability of stationary cross-flow vortices in a boundary layer over a swept wing, and identified two families of secondary instability modes, named as $y$ - and $z$ -modes, which were found to be associated with the wall-normal and spanwise shears, respectively, and both have high frequencies. Wassermann & Kloker (2002) identified a high-frequency and a low-frequency modes, both being of $z$ -mode, and showed that the secondary instability is convective, as is the primary instability. Travelling cross-flow vortices also support secondary instability modes in different frequency ranges, shown by Wassermann & Kloker (2003). These modes can reach a pronounced saturation, which is ended by rapid secondary amplification of a medium-frequency mode, prior to final breakdown.

In this paper, we examine the effect of a wavy-wall roughness of small height on a three-dimensional boundary layer. We consider the mechanism of roughness altering the base flow thereby inducing new instability. As mentioned above, in order to pursue a mechanism of this kind it is necessary to calculate first the roughness-distorted base flow. Previous studies either made linearisation approximation (which is inappropriate), or employed numerical methods such as the IBL method and Navier–Stokes (N–S) solvers. The streamwise length scale of the roughness considered is much greater than the boundary-layer thickness so that local stability analysis may be performed based on parallel-flow approximation. In the present study, the roughness length scale is comparable to the boundary-layer thickness, and we adopt the self-consistent high-Reynolds-number asymptotic approach to calculate the nonlinear roughness-distorted flow. The main response generated by the roughness is in the WL, but unique to three-dimensional boundary layers, strong response arises also in the critical layer (CL) located in the bulk of the flow. The roughness height is taken to be much smaller than the boundary-layer thickness but just large enough to induce $\mathit{O} (1)$ vorticity in these two layers. We then study secondary instability of this distorted base flow using an eigenvalue and initial-value calculations. The instability modes are different from any of those previously known in that they have short wavelengths comparable to the width of the WL and CL. The instability in the CL also has high frequencies. They are completely trapped, and determined by the local flows, in these layers, unaffected by the overall boundary-layer flow.

The rest of the paper is organised as follows. The roughness-distorted flow is analysed in § 2 based on the high-Reynolds-number assumption, under which the flow field acquires a structure consisting of a WL, the main layer and a CL. Both the WL and CL are nonlinearly distorted. The numerical results for the roughness-distorted flow are presented in § 3. The secondary instability analyses for nonlinearly distorted CL and WL are formulated in § 4, where both the temporal and spatial secondary instabilities are investigated. In § 5, we present the numerical results for the small-scale local secondary instability. A summary is given in § 6, where future investigations are discussed. The supplementary material is available at https://doi.org/10.1017/jfm.2025.10794, wherein details of derivation and computation, and additional numerical results are presented.

2. Problem formulation for the roughness-distorted flow

We consider a three-dimensional boundary-layer flow over a wall with distributed periodic surface roughness on it. The flow is described using the Cartesian coordinates $ (x^{\ast }, y^{\ast }, z^{\ast } )$ with its origin located at the leading edge of the wall, where $x^{\ast }$ , $y^{\ast }$ and $z^{\ast }$ denote the distances in the chordwise, wall-normal and spanwise directions, respectively. The superscript $\ast$ signifies a dimensional quantity. We consider the downstream region where the chordwise distance to the leading edge is $L^{\ast }$ and the local boundary-layer thickness is $\delta ^{\ast }$ . The coordinates and the velocities are non-dimensionalised by $\delta ^{\ast }$ and the chordwise slip velocity $U_{\infty L}^{\ast }$ outside the boundary layer at $L^{\ast }$ , respectively. The Reynolds number is defined as

(2.1) \begin{equation} R = U_{\infty L}^{\ast } \delta ^{\ast } / \nu ^{\ast }, \end{equation}

with $\nu ^{\ast }$ being the kinematic viscosity of the fluid. Since $L^{\ast }$ is chosen to be fixed, $R$ is a constant. We take $R \gg 1$ to be asymptotically large.

The periodically distributed surface roughness is chosen to be of the form

(2.2) \begin{equation} y_w := h R^{-1/3} F(x,z) = h R^{-1/3} \sum _{n=-\infty }^{\infty } F_{n} \textrm{e}^{\textrm{i} n (\alpha _w x + \beta _w z)}, \end{equation}

where $\alpha _w$ and $\beta _w$ are the chordwise and spanwise wavenumbers of surface waviness, respectively, and the parameter $h \sim \mathit{O} (1)$ controls the roughness height. The $\mathit{O} (h R^{-1/3})$ height is much smaller than the boundary-layer thickness for $R \gg 1$ and might correspond to a more-or-less micron-sized height in the dimensional setting, given that typical local boundary-layer thickness $\delta ^{\ast } \sim 1$ mm, Reynolds number $R \sim 10^{3}$ and $h = 0.4$ (at which vigorous small-scale instability occurs according to our calculations whilst its onset may take place at a fraction of this height). We take $F_{0} = 0$ without losing generality since a uniform elevation can be removed by a shift of the $y$ -coordinate.

The flow field in the boundary layer is characterised by an asymptotic structure including a WL, the main layer and a CL as is sketched in figure 1. The immediate response of the boundary layer to surface roughness is in a viscous WL, which is introduced to satisfy the no-slip boundary conditions. The inviscid main layer occupies the bulk of the boundary layer and has a singular point $y_c$ . The CL is introduced to regularise the singularity. We will analyse each layer below. The flow is governed by the incompressible N–S equations,

(2.3a,b) \begin{gather} {\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{u} = 0, \quad \partial _{t} \boldsymbol{u} + (\boldsymbol{u} \boldsymbol{\cdot }{\boldsymbol{\nabla }}) \boldsymbol{u} = - {\boldsymbol{\nabla }} p + R^{-1} {\boldsymbol{\nabla} }^{2} \boldsymbol{u}, \end{gather}

where $\boldsymbol{u} = (u, v, w)$ with $u$ , $v$ and $w$ denoting the chordwise, wall-normal and spanwise velocities, respectively, and $p$ denotes the pressure non-dimensionalised by $\rho ^{\ast } U_{\infty L}^{\ast 2}$ , with $\rho ^{\ast }$ being the reference density.

Figure 1. Asymptotic structure and scaling of the distorted flow field. Surface roughness significantly alters the base flow within the WL (§ 2.1) and CL (§ 2.3), leading to the emergence of $\mathit{O} (1)$ vorticities $\bar {\varOmega }$ and $\tilde {\varOmega }$ , which render the flows in these layers susceptible to small-scale secondary instability. The subscript $B$ denotes the base-flow quantities, while the subscripts $s$ and $m$ denote the streaming and the forced perturbation, respectively, in the main layer (§ 2.2).

We assume that the base flow is spanwise uniform, i.e. independent of the spanwise coordinate $z$ . It varies slowly in the chordwise direction and is described by a spatial variable,

(2.4) \begin{equation} \bar {x} = R^{-1} x. \end{equation}

The base flow scale as

(2.5) \begin{equation} ( \boldsymbol{u}_{B}, p_{B} ) = \bigl (U_{\!B} (\bar {x},y), R^{-1} V_{B} (\bar {x},y), W_{\!B} (\bar {x},y), P_{B} (\bar {x}) \bigr ). \end{equation}

Substituting (2.5) into (2.3), we obtain the steady three-dimensional boundary-layer equations,

(2.6a) \begin{align} U_{\!B, \bar {x}} + V_{B, y} &= 0, \end{align}

(2.6b) \begin{align} U_{\!B} U_{\!B, \bar {x}} + V_{B} U_{\!B, y} &= - P_{B, \bar {x}} + U_{\!B, yy}, \end{align}

(2.6c) \begin{align} U_{\!B} W_{B,\bar {x}} + V_{B} W_{B, y} &= W_{B,yy}, \end{align}

and $P_{B, y} = 0$ , which allows the pressure gradient to be determined as

(2.7) \begin{equation} P_{B, \bar {x}} = - U_{\infty } U_{\infty , \bar {x}}, \end{equation}

with $U_{\infty }$ being the chordwise slip velocity outside the boundary layer.

As in numerous previous studies, the FSC flow is used as a convenient model for three-dimensional boundary layers. The chordwise and spanwise slip velocities are given by

(2.8a,b) \begin{gather} U_{\infty }^{\ast } (x^{\ast }) = U_{\infty L}^{\ast } \bigl ( x^{\ast } / L^{\ast } \bigr )^{m}, \quad W_{\infty }^{\ast } = \textrm{const}, \end{gather}

where $m$ is the acceleration parameter. The local sweep angle at $L^{\ast }$ is defined as $\phi = \arctan (W_{\infty }^{\ast } / U_{\infty L}^{\ast }) = \arctan (W_{\infty })$ . The FSC flow admits a similarity solution, which can be expressed in terms of a similarity variable $\eta$ ,

(2.9) \begin{equation} \eta = \left ( \frac {m+1}{2} \frac {U_{\infty }^{\ast } (x^{\ast })}{\nu ^{\ast } x^{\ast }} \right )^{1/2} y^{\ast } = \left ( \frac {m+1}{2} \frac {R}{L} \right )^{1/2} \bar {x}^{(m-1)/2} y \end{equation}

where the constant $L$ is to be determined below (Högberg & Henningson 1998). According to its definition, the local boundary-layer thickness $\delta ^{\ast }$ is expressed as

(2.10a,b) \begin{gather} \delta ^{\ast } = \left ( \frac {m+1}{2} \frac {U_{\infty L}^{\ast }}{\nu ^{\ast } L^{\ast }} \right )^{-1/2} c, \quad c = \int _{0}^{\infty } \left (1 - f' (\eta ) \right ) \textrm{d} \eta , \end{gather}

where the first implies $L = (m+1) R / (2 c^{2})$ and $f (\eta )$ is a function introduced to express the stream function, whose expression is omitted here. We then rewrite $\eta$ as $\eta = c \bar {x}^{(m-1)/2} y$ . The velocity components of the base flow can be expressed as

(2.11a,b) \begin{align} U_{\!B} = \bar {x}^{m} f' (\eta ), \quad V_{B} &= - \left ( \frac {m+1}{2} \right )^{1/2} \bar {x}^{(m - 1) / 2} \left ( f (\eta ) + \frac {m - 1}{m + 1} \eta f' (\eta ) \right )\!, \end{align}

(2.11c) \begin{align} W_{\!B} &= W_{\infty } g(\eta ). \qquad \qquad \qquad \end{align}

Substituting (2.11) into (2.6) gives the equations governing $f$ and $g$ ,

(2.12a,b) \begin{align} f^{\prime \prime \prime } + f f^{\prime \prime } + \frac {2 m }{m + 1} \left ( 1 - f^{\prime 2} \right ) = 0, \quad g^{\prime \prime } + f g^{\prime } = 0, \end{align}

which are subject to the boundary and far-field matching conditions

(2.13a–e) \begin{align} f = f' = g = 0 \quad \textrm{at} \, \, \, \eta = 0; \quad f' \to 1, \quad g \to 1 \quad \textrm{as} \, \, \, \eta \to \infty . \end{align}

The reader is reminded that although a simple boundary layer with a self-similar solution is chosen, the ensuing mathematical formulations of the nonlinear distortions and the stability of the distorted flows are valid for a general boundary layer admitting no similarity solution.

We decompose the boundary-layer flow field into the base flow and the roughness-induced distortion,

(2.14) \begin{equation} \boldsymbol{u} = \boldsymbol{U}_{\!B} + \epsilon \boldsymbol{u}_{d} + \mathit{O} (\epsilon ^2), \end{equation}

where $\epsilon \ll 1$ is the magnitude of the disturbance. Next, we determine $\epsilon$ in terms of $R$ for the roughness height specified by (2.2). Near the wall, the distortion $\boldsymbol{u}_{d}$ varies much more rapidly in the wall-normal direction than the base flow $\boldsymbol{U}_{B}$ , which allows us to Taylor expand the base flow about $y = 0$ . We take the chordwise component as an example, which satisfies the no-slip boundary condition at the wall

(2.15) \begin{equation} 0 = \big (\big . U_{\!B} \big \rvert _{y=0} + \big . U_{\!B,y} \big \rvert _{y=0} h R^{-1/3} F(x,z) + \mathit{O}\big(h^2 R^{-2/3}\big) \big ) + \epsilon \left . u_{d} \right \rvert _{y = h R^{-1/3} F} + \mathit{O}(\epsilon ^2), \end{equation}

which suggests the asymptotic scale of the induced disturbance,

(2.16) \begin{equation} \epsilon \sim h R^{-1/3}. \end{equation}

When $h \ll 1$ the boundary condition of $u_{d}$ can be imposed at $y = 0$ ,

(2.17) \begin{equation} \left . u_{d} \right \rvert _{y = 0} = - \left . U_{\!B,y} \right \rvert _{y = 0} F(x,z). \end{equation}

2.1. Wall layer

The immediate response of the flow to the surface roughness is in a viscous WL with an $\mathit{O}(R^{-1/3})$ width, which is derived by balancing the inertia and viscous terms in the chordwise momentum equations in (2.3). We thus introduce in this layer a local wall-normal coordinate,

(2.18) \begin{equation} Y = R^{1/3} y = \mathit{O} (1), \end{equation}

in terms of which the wavy wall is specified as $Y_{w} = h F(x,z)$ . Note that with $h = \mathit{O} (1)$ , the chordwise and spanwise velocities of the induced disturbance are comparable to those of the background flow, both being of $\mathit{O} (R^{-1/3})$ . It follows that the total flow field in the WL is scaled as

(2.19) \begin{equation} (u,v,w,p ) = \bigl (R^{-1/3} U_{w}, R^{-2/3} V_{w} , R^{-1/3} W_{w}, R^{-2/3} P_{w} \bigr ). \end{equation}

Substitution of (2.19) into (2.3) gives boundary-layer equations,

(2.20a) \begin{align} U_{w, x} + V_{w, Y} + W_{w, z} = & 0, \end{align}

(2.20b) \begin{align} U_{w} U_{w, x} + V_{w} U_{w, Y} + W_{w} U_{w, z} = & - P_{w, x} + U_{w, Y Y}, \end{align}

(2.20c) \begin{align} 0 = & P_{w, Y}, \end{align}

(2.20d) \begin{align} U_{w} W_{w, x} + V_{w} W_{w, Y} + W_{w} W_{w, z} = & - P_{w, z} + W_{w, Y Y}. \end{align}

The no-slip and impermeability boundary conditions at the roughness surface read

(2.21a–c) \begin{align} U_{w} = V_{w} = W_{w} = 0 \quad \textrm{at} \, \, \, Y = Y_{w}; \end{align}

while the matching conditions read

(2.22a,b) \begin{align} ( U_{w}, W_{w} ) \to ( \lambda _{1}, \lambda _{3} ) Y + h ( \bar {U}_{s}, \bar {W}_{s}) \quad \textrm{as} \, \, \, Y \to \infty , \end{align}

where $\lambda _{1}$ and $\lambda _{3}$ denote the chordwise and spanwise wall shears, respectively; $\bar {U}_{s}$ and $\bar {W}_{s}$ are the chordwise and spanwise streaming caused by nonlinear interactions in the WL and are independent of $Y$ , as is shown in the supplementary material (§ S1).

In order to calculate the flow field in a rectangular domain, we apply the Prandtl transformation,

(2.23a,b) \begin{align} \bar {Y} = Y - Y_{w}, \quad \bar {V}_w = V_w - h \left ( \bar {U}_w F_x + \bar {W}_w F_z \right )\!, \end{align}

where an overbar signifies variables under the transformed coordinates. The governing equations for the transformed variables remain the same form as (2.20), while the boundary conditions are transformed into

(2.24a–c) \begin{align} \bar {U}_w = \bar {V}_w = \bar {W}_w = 0 \quad &\textrm{at} \, \, \, \bar {Y} = 0; \end{align}

(2.24d,e) \begin{align} ( \bar {U}_{w}, \bar {W}_{w} ) \to ( \lambda _{1}, \lambda _{3} ) (\bar {Y} + h F) + h ( \bar {U}_{s}, \bar {W}_{s}) \quad &\textrm{as} \, \, \, \bar {Y} \to \infty . \ \end{align}

We then decompose the total flow field into the base flow and the disturbance,

(2.25) \begin{align} (\bar {U}_w,\bar {V}_w,\bar {W}_w,\bar {P}_w) &= (\lambda _1 (\bar {Y}+h F),- h \lambda _{\textit{eff}} F_{\zeta } (\bar {Y}+h F),\lambda _3 (\bar {Y}+ h F),P_B)\nonumber\\&\quad + h ( \bar {U}^{\dagger }, \bar {V}^{\dagger }, \bar {W}^{\dagger },\bar {P}^{\dagger }), \end{align}

and express the latter as Fourier series,

(2.26) \begin{equation} (\bar {U}^{\dagger },\bar {V}^{\dagger },\bar {W}^{\dagger },\bar {P}^{\dagger }) = \sum _{n=-\infty }^{\infty } (\bar {U}_n(\bar {Y}),\bar {V}_n (\bar {Y}) ,\bar {W}_n (\bar {Y}),\bar {P}_n (\bar {Y})) \textrm{e}^{\textrm{i} n \alpha _w \zeta }, \end{equation}

where we introduce the skewed coordinate $\zeta = x + (\beta _w / \alpha _w) z$ and the effective base-flow wall shear $\lambda _{\textit{eff}} = \lambda _1 + (\beta _w / \alpha _w) \lambda _3$ . To keep the physical quantities real, we require

(2.27) \begin{equation} \bar {q}_{-n} = (\bar {q}_n)_{cc}, \end{equation}

where $\bar {q}_n$ stands for any of the Fourier components on the right-hand side of (2.26). Substituting (2.25)–(2.26) into (2.20) gives the governing equations for the Fourier components,

(2.28a) \begin{align} \textrm{i} n \alpha _w \bar {U}_n + \bar {V}_{n, \bar {Y}} + \textrm{i} n \beta _w \bar {W}_n & = 0, \end{align}

(2.28b) \begin{align} \textrm{i} n \alpha _w \lambda _{\textit{eff}} \bar {Y} \bar {U}_n + \lambda _1 \bar {V}_n + \textrm{i} n \alpha _w \bar {P}_n - \bar {U}_{n, \bar {Y} \bar {Y}} & = - h ( \lambda _{\textit{eff}} \mathcal{R}_n^x + \mathcal{N}_n^x ), \end{align}

(2.28c) \begin{align} \bar {P}_{n, \bar {Y}} & = 0, \end{align}

(2.28d) \begin{align} \textrm{i} n \alpha _w \lambda _{\textit{eff}} \bar {Y} \bar {W}_n + \lambda _3 \bar {V}_n + \textrm{i} n \beta _w \bar {P}_n - \bar {W}_{n, \bar {Y} \bar {Y}} & = - h ( \lambda _{\textit{eff}} \mathcal{R}_n^z + \mathcal{N}_n^z ), \end{align}

where

(2.29a) \begin{align} \mathcal{R}_n^x &= \left \{ \widehat {F \bar {U}_{\zeta }^{\dagger }} \right \}_n + \frac {\lambda _1}{\lambda _{\textit{eff}}} \left ( \left \{ \widehat { F_{\zeta } \bar {U}^{\dagger }} \right \}_n + \frac {\beta _w}{\alpha _w} \left \{ \widehat { F_{\zeta } \bar {W}^{\dagger }} \right \}_n \right ) - \bar {Y} \left \{ \widehat { F_{\zeta } \bar {U}_{\bar {Y}}^{\dagger } } \right \}_n - h \left \{ \widehat {F F_{\zeta } \bar {U}_{\bar {Y}}^{\dagger }} \right \}_n , \end{align}

(2.29b) \begin{align} \mathcal{R}_n^z &= \left \{ \widehat {F \bar {W}_{\zeta }^{\dagger }} \right \}_n + \frac {\lambda _3}{\lambda _{\textit{eff}}} \left ( \left \{ \widehat { F_{\zeta } \bar {U}^{\dagger }} \right \}_n + \frac {\beta _w}{\alpha _w} \left \{ \widehat { F_{\zeta } \bar {W}^{\dagger }} \right \}_n \right ) - \bar {Y} \left \{ \widehat { F_{\zeta } \bar {W}_{\bar {Y}}^{\dagger } } \right \}_n - h \left \{ \widehat {F F_{\zeta } \bar {W}_{\bar {Y}}^{\dagger }} \right \}_n , \end{align}

(2.29c) \begin{align} \mathcal{N}_n^x &= \textrm{i} n \alpha _w \left \{ \widehat {\bar {U}^{\dagger } \bar {U}^{\dagger }} \right \}_n + \left \{ \widehat {\bar {V}^{\dagger } \bar {U}^{\dagger } } \right \}_{n,\bar {Y}} + \textrm{i} n \beta _w \left \{ \widehat { \bar {W}^{\dagger } \bar {U}^{\dagger } } \right \}_n, \end{align}

(2.29d) \begin{align} \mathcal{N}_n^z &= \textrm{i} n \alpha _w \left \{ \widehat {\bar {U}^{\dagger } \bar {W}^{\dagger } } \right \}_n + \left \{ \widehat { \bar {V}^{\dagger } \bar {W}^{\dagger }} \right \}_{n,\bar {Y}} + \textrm{i} n \beta _w \left \{ \widehat {\bar {W}^{\dagger } \bar {W}^{\dagger }} \right \}_n, \end{align}

with $\widehat {\qquad }$ denoting the Fourier transform with respect to $\zeta$ and the subscript $n$ representing the $n$ th harmonic. The boundary conditions are

(2.30a–c) \begin{align} \bar {U}_{n} = - \lambda _1 F_n, \quad \bar {V}_{n} = h \lambda _{\textit{eff}} \left \{ \widehat {F F_{\zeta }} \right \}_n, \quad \bar {W}_{n} = - \lambda _3 F_n \quad \textrm{at} \, \, \, \bar {Y} = 0. \end{align}

Note that $\bar {U}_{0} (0) = \bar {W}_{0} (0) = 0$ since $F_{0} = 0$ , and $\bar {V}_{0} (0) = 0$ since $ \{ \widehat {F F_{\zeta }} \}_{n} = \textrm{i} n \alpha _{w} \{ \widehat {F^{2}} \}_{n} / 2$ . For all $n \ne 0$ Fourier components, we have the matching conditions

(2.31a,b) \begin{align} \bar {U}_{n} \to 0, \quad \bar {W}_{n} \to 0 \quad \textrm{as} \, \, \, \bar {Y} \to \infty . \end{align}

Each Fourier component of the pressure $\bar {P}_{n}$ is a constant according to (2.28c ). Combining (2.28b ) and (2.28d ) and taking the limit $\bar {Y} \to \infty$ , we obtain the expression for a $n \ne 0$ Fourier component of the pressure

(2.32) \begin{equation} \bar {P}_{n} = - \alpha _w^2 k_w^{-2} \lambda _{\textit{eff}} \bigl [ \bar {V}_{n,\infty } / (\textrm{i} n \alpha _w) + h F_n \bar {U}_{s \textit{eff}} \bigr ], \end{equation}

where $\bar {V}_{n,\infty } \equiv \bar {V}_{n} \rvert _{\bar {Y} \to \infty }$ and

(2.33a,b) \begin{align} k_{w} = \sqrt {\alpha _w^2 + \beta _w^2}, \quad \bar {U}_{s \textit{eff}} = \bar {U}_{s} + (\beta _{w} / \alpha _{w}) \bar {W}_{s}. \end{align}

The Fourier component with $n = 0$ , also referred to as the mean-flow distortion, is caused by nonlinear interactions of the roughness-induced disturbance within the WL. It exhibits a different asymptotic behaviour as $\bar {Y} \to \infty$ so that a separate formulation is required. Taking $n = 0$ in (2.28a ) gives $\bar {V}_{0, \bar {Y}} = 0$ , integration of which with respect to $\bar {Y}$ with use of the boundary condition (2.30b ) leads to $\bar {V}_{0} \equiv 0$ . Setting $n = 0$ in (2.28b ) and (2.28d ), we obtain the governing equations for $\bar {U}_{0}$ and $\bar {W}_{0}$ ,

(2.34a,b) \begin{align} \bar {U}_{0, \bar {Y} \bar {Y}} = h \bigl ( \mathcal{R}_0^x + \mathcal{N}_0^x \bigr ), \quad \bar {W}_{0, \bar {Y} \bar {Y}} = h \bigl ( \mathcal{R}_0^z + \mathcal{N}_0^z \bigr ). \end{align}

Integrating (2.34) twice, we obtain the solutions for $\bar {U}_{0}$ and $\bar {W}_{0}$ ,

(2.35a) \begin{align} \bar {U}_{0} = h \int _{0}^{\bar {Y}} ( \bar {Y} - Y_1 ) \bigl ( \mathcal{R}_0^x (Y_1) + \mathcal{N}_0^x(Y_1) \bigr ) \textrm{d} Y_1 + \bar {U}_{0,\bar {Y}} (0) \bar {Y}, \end{align}

(2.35b) \begin{align} \bar {W}_{0} = h \int _{0}^{\bar {Y}} ( \bar {Y} - Y_1 ) \bigl ( \mathcal{R}_0^z (Y_1) + \mathcal{N}_0^z(Y_1) \bigr ) \textrm{d} Y_1 + \bar {W}_{0,\bar {Y}} (0) \bar {Y}, \end{align}

where

(2.36a,b) \begin{align} \bar {U}_{0,\bar {Y}} (0) = - h \frac {\alpha _w^2}{k_w^2} \lambda _{\textit{eff}} \sum _{\substack {k=-\infty \\ k \ne 0}}^{\infty } F_{-k} \bar {V}_{k, \infty } , \ \bar {W}_{0,\bar {Y}} (0) = - h \frac {\alpha _w \beta _w}{k_w^2} \lambda _{\textit{eff}} \sum _{\substack {k=-\infty \\ k \ne 0}}^{\infty } F_{-k} \bar {V}_{k, \infty }, \end{align}

which are fixed by considering the large- $\bar {Y}$ asymptotes of $\bar {U}_{0}$ and $\bar {W}_{0}$ , as discussed in the supplementary material (§ S1). It can be shown that

(2.37a,b) \begin{align} \bar {U}_{0} \to \bar {U}_s, \quad \bar {W}_{0} \to \bar {W}_s \quad \textrm{as} \, \, \, \bar {Y} \to \infty , \end{align}

with both $\bar {U}_{s}$ and $\bar {W}_{s}$ being constants. This nonlinear streaming effect induces slip velocities which act at the bottom of the main layer, driving a mean-flow distortion of $\mathit{O} (R^{-1/3})$ there.

The WL generates a ‘blowing velocity’ of $\mathit{O} (R^{-2/3})$ as $\bar {Y} \to \infty$ , which can be derived by integrating (2.20a ) and substituting into it the matching conditions (2.24d , e ),

(2.38) \begin{equation} \bar {V}_{w} \to - h ( \lambda _{1} F_{x} + \lambda _{3} F_{z} ) \bar {Y} + \bar {V}_{t} (x, z) \quad \textrm{as} \, \, \, \bar {Y} \to \infty , \end{equation}

where $\bar {V}_{t}$ has the expression

(2.39) \begin{equation} \bar {V}_{t} = - h \sum _{n = - \infty }^{\infty } \biggl [ \int _{0}^{\infty } ( \textrm{i} n \alpha _w \bar {U}_n + \textrm{i} n \beta _w \bar {W}_n ) \textrm{d} Y \biggr ] \textrm{e}^{\textrm{i} n \alpha _w \zeta }. \end{equation}

Transforming the above quantities back to the original coordinates, we can calculate the blowing velocity of the WL,

(2.40) \begin{equation} V_w \to \bar {V}_t + h^2 \left (\lambda _{\textit{eff}} F + \bar {U}_{s \textit{eff}} \right ) F_{\zeta } \quad \textrm{as} \, \, \, Y \to \infty , \end{equation}

which acts at the bottom of the main layer to induce a response of the same order there; the first and second terms on right-hand side of (2.40) represent, respectively, the transpiration velocity induced by the viscous motion in the WL and the pure geometric deflection of the streamlines by the roughness topography. Physically, this blowing velocity can be interpreted as a blowing/suction unit on the wall, which generates three-dimensional vorticity in the main layer.

The WL dynamics is strongly nonlinear with $\mathit{O} (R^{-1/3})$ chordwise and spanwise velocities. The associated vorticity is even greater with an $\mathit{O} (1)$ amplitude, making the WL susceptible to small-scale secondary instability as will be shown in § 4.

2.2. Main layer

The main layer occupies the bulk of the boundary layer, where $y = \mathit{O} (1)$ . The $\mathit{O} (R^{-2/3})$ blowing velocity and the $\mathit{O} (R^{-1/3})$ chordwise and spanwise streaming produce corresponding responses in the main layer. The flow field there is expressed as

(2.41) \begin{align} (\boldsymbol{u}, p ) &= \big( U_{\!B}, R^{-1} V_{B}, W_{\!B},P_B \big) \nonumber \\ & \quad + R^{-1/3} \Big( U_{s}^{\langle 1 \rangle }, R^{-1} V_{s}^{\langle 1 \rangle },W_{s}^{\langle 1 \rangle }, P_{s}^{\langle 1 \rangle } \Big) + R^{-2/3} \Big( U_{s}^{\langle 2 \rangle }, R^{-1} V_{s}^{\langle 2 \rangle }, W_{s}^{\langle 2 \rangle }, P_{s}^{\langle 2 \rangle } \Big) \nonumber \\ & \quad + R^{-2/3} \sum _{\substack {n=-\infty \\ n \ne 0}}^{\infty } ( u_{n} \left ( y \right )\!, v_{n}\left ( y \right )\!, w_{n}\left ( y \right )\!, p_{n}\left ( y \right ) ) e^{i n \alpha _w \zeta } + {\cdots} , \end{align}

where the subscript $s$ denotes the main-layer streaming driven by the mean-flow distortion from the WL and the CL, which will be shown later in § 2.3. This streaming, though larger, is not our main interest since it does not affect the dynamics in the CL nor the secondary instability. We focus on the $\mathit{O} (R^{-2/3})$ forced disturbance, which is governed by

(2.42a) \begin{align} \textrm{i} n \alpha _w u_{n} + v_{n,y}+ \textrm{i} n \beta _w w_{n} = 0, \end{align}

(2.42b) \begin{align} \textrm{i} n \alpha _w U_{\!B\textit{eff}} u_{n} + v_{n} U_{\!B,y} + \textrm{i} n \alpha _w p_{n} = 0, \end{align}

(2.42c) \begin{align} \textrm{i} n \alpha _w U_{\!B\textit{eff}} v_{n} + p_{n,y} = 0, \end{align}

(2.42d) \begin{align} \textrm{i} n \alpha _w U_{\!B\textit{eff}} w_{n} + v_{n} W_{B,y} + \textrm{i} n \beta _w p_{n} = 0, \end{align}

where we have introduced the effective base flow

(2.43) \begin{equation} U_{\!B \textit{eff}} = U_B + (\beta _w / \alpha _{w}) W_B. \end{equation}

Equations in (2.42) can be reduced to the steady Rayleigh equation

(2.44) \begin{equation} v_{n, y y} - \left (n^2 k_w^2 + \frac {U_{\!B \textit{eff}, yy}}{U_{\!B \textit{eff}}} \right ) v_{n} = 0 \quad n \ne 0, \end{equation}

which is subject to the inhomogeneous boundary conditions,

(2.45a,b) \begin{align} v_{n} \to V_{w \infty ,n} \quad \textrm{as} \, \, \, y \to 0; \quad v_{n} \to 0 \quad \textrm{as} \, \, \, y \to \infty . \end{align}

Here $V_{w \infty ,n}$ represents the Fourier component on the right-hand side of (2.40).

The Rayleigh (2.44) becomes singular at the critical level $y_c$ , where $U_{\!B \textit{eff}} (y_c) = 0$ while $U_{\!B \textit{eff}, yy} (y_c) \ne 0$ in general. Near $y_c$ , a local solution to (2.44) can be constructed using the Frobenius method, which gives

(2.46) \begin{equation} v_{n} = a_{n}^{\pm } \phi _{an} (y)+b_{n} \Bigl [\phi _{bn} (y) + \frac {U_c^{\prime\prime}}{U_c^{\prime}} \ln {\left |y-y_c\right |} \phi _{an}(y) \Bigr ] \end{equation}

with

(2.47a,b) \begin{align} \phi _{an} = \tilde {y} + \frac {1}{2} \frac {U_c^{\prime \prime }}{U_c^{\prime }} \tilde {y} ^ 2 + {\cdots} , \quad \phi _{bn} = 1 + \frac {1}{2} \left (n^2 k_w^2 + \frac {U_c^{\prime \prime \prime }}{U_c^{\prime }} - 2 \frac {U_c^{\prime \prime 2}}{U_c^{\prime 2}} \right ) \tilde {y}^2 + {\cdots} , \end{align}

and $\tilde {y} = y - y_{c}$ . The subscript $c$ signifies the quantities related to $U_{\!B \textit{eff}} (y_{c})$ , and a prime denotes the derivative with respect to $y$ . There are velocity jumps $(a_n^+ - a_n^-)$ across $y_c$ , which will be determined by analysing the CL dynamics.

In the main layer, we can also express $\{ u_n,w_n,p_n \}$ in terms of $v_n$ and its derivative,

(2.48a) \begin{align} u_{n} &= \frac {\textrm{i} \alpha _w}{n k_w^2} \left (v_{n,y} - \frac {U_{\!B\textit{eff},y}}{U_{\!B\textit{eff}}} v_{n} \right ) + \frac {\textrm{i} U_{\!B,y}}{n \alpha _w U_{\!B\textit{eff}}}v_{n}, \end{align}

(2.48b) \begin{align} w_{n} &= \frac {\textrm{i} \beta _w}{n k_w^2} \left (v_{n,y} - \frac {U_{\!B\textit{eff},y}}{U_{\!B\textit{eff}}} v_{n} \right ) + \frac {\textrm{i} W_{B,y}}{n \alpha _w U_{\!B\textit{eff}}}v_{n}, \end{align}

(2.48c) \begin{align} p_{n} &= - \frac {\textrm{i} \alpha _w}{n k_w^2} \left ( v_{n,y} - \frac {U_{\!B\textit{eff},y}}{U_{\!B \textit{eff}}} v_{n} \right ) U_{\!B \textit{eff}}. \end{align}

2.3. Critical layer

To resolve the singularity of (2.44) at the critical level $y_c$ , we need to consider a CL around $y_c$ with an $\mathit{O} (R^{-1/3})$ width, within which the viscous effect is reinstated. We introduce a new local coordinate in the CL,

(2.49) \begin{equation} \tilde {Y} = R^{1/3} (y - y_c) = \mathit{O} (1). \end{equation}

In order to determine the form of expansion for the disturbance in this layer, the main-layer solutions near the critical level, (2.46) and (2.48), are rewritten in terms of $\tilde {Y}$ ,

(2.50a) \begin{align} u &\to U_{\!B} (y_c) + R^{-\frac {1}{3}} \biggl [ U_{\!B}^{\prime} (y_c) \tilde {Y} + U_{s}^{\langle 1 \rangle } (y_c) + \sum _{\substack {n = -\infty \\ n \ne 0}}^{\infty } \frac {\textrm{1}}{\textrm{i} n \alpha _w} \left (\frac {\alpha _w^2}{k_w^2} - \frac {U_{\!B}^{\prime}(y_\textit{c})}{U_c^{\prime}} \right ) b_n \frac {1}{\tilde {Y}} \textrm{e}^{\textrm{i} n \alpha _w \zeta } \biggr ] \nonumber \\ &\quad + R^{-\frac {2}{3}} \biggl [\frac {U_{\!B}^{\prime\prime} (y_\textit{c})}{2} \tilde {Y}^2 + U_{s}^{\langle 2 \rangle } (y_c) - \sum _{\substack {n = -\infty \\ n \ne 0}}^{\infty } \frac {1}{\textrm{i} n \alpha _w U_c^{\prime}} \biggl ( U_{\!B}^{\prime} (y_c) a_{n}^{\pm } + U_{\!B}^{\prime\prime} (y_c) b_n \nonumber \\ &\quad + \frac {U_c^{\prime\prime}}{2} \left (\frac {\alpha _w^2}{k_w^2} - \frac {U_{\!B}^{\prime}(y_\textit{c})}{U_c^{\prime}} \right ) b_n + \frac {U_{\!B}^{\prime}(y_c) U_c^{\prime\prime}}{U_c^{\prime}} b_n \ln \left | \tilde {Y} \right | \biggr ) \textrm{e}^{\textrm{i} n \alpha _w \zeta } \biggr ] + {\cdots} , \end{align}

(2.50b) \begin{align} v &\to R^{-1} V_{B}(y_c) + R^{-\frac {2}{3}} \sum _{\substack {n = \infty \\ n \ne 0}}^{\infty } b_n \textrm{e}^{\textrm{i} n \alpha _w \zeta } + R^{-1} \sum _{\substack {n = \infty \\ n \ne 0}}^{\infty } \biggl ( a_{n}^{\pm } + \frac {U_c^{\prime\prime}}{U_c^{\prime}} b_n \ln \left | \tilde {Y} \right | \biggr ) \tilde {Y} \textrm{e}^{\textrm{i} n \alpha _w \zeta } + {\cdots} , \end{align}

(2.50c) \begin{align} w &\to \! W_{\!B} (y_c) + R^{-\frac {1}{3}} \biggl [\! W_{\!B}^{\prime} (y_c) \tilde {Y} + W_{s}^{\langle 1 \rangle } (y_c) + \sum _{\substack {n = -\infty \\ n \ne 0}}^{\infty } \frac {\textrm{1}}{\textrm{i} n \alpha _w} \!\left (\!\frac {\alpha _w \beta _w}{k_w^2} - \frac {W_{\!B}^{\prime}(y_\textit{c})}{U_c^{\prime}} \!\right )\! b_n \frac {1}{\tilde {Y}} \textrm{e}^{\textrm{i} n \alpha _w \zeta }\! \biggr ] \nonumber \\ &\quad + R^{-\frac {2}{3}} \biggl [\frac {W_{\!B}^{\prime\prime} (y_\textit{c})}{2} \tilde {Y}^2 + W_{s}^{\langle 2 \rangle } (y_c) - \sum _{\substack {n = -\infty \\ n \ne 0}}^{\infty } \frac {1}{\textrm{i} n \alpha _w U_c^{\prime}} \biggl ( W_{\!B}^{\prime} (y_c) a_{n}^{\pm } + W_{\!B}^{\prime\prime} (y_c) b_n \nonumber \\ &\quad + \frac {U_c^{\prime\prime}}{2} \left (\frac {\alpha _w \beta _w}{k_w^2} - \frac {W_{\!B}^{\prime}(y_\textit{c})}{U_c^{\prime}} \right ) b_n + \frac {W_{\!B}^{\prime}(y_c) U_c^{\prime\prime}}{U_c^{\prime}} b_n \ln \left | \tilde {Y} \right | \biggr ) \textrm{e}^{\textrm{i} n \alpha _w \zeta } \biggr ] + {\cdots} , \end{align}

(2.50d) \begin{align} p &\to P_{B} (y_c) + R^{-\frac {2}{3}} \sum _{\substack {n = -\infty \\ n \ne 0}}^{\infty } \frac {\textrm{i} \alpha _w}{n k_w^2} U_c^{\prime} b_n \textrm{e}^{\textrm{i} n \alpha _w \zeta } + {\cdots} , \end{align}

where the terms associated with $\ln R$ are excluded as they are the natural continuations of the outer solutions and remain regular in the CL. We note that $u$ and $w$ have a simple-pole singularity, while $v$ has a logarithmic singularity. Based upon the asymptotic matching principle, we express the flow field in the CL as

(2.51) \begin{align} (\boldsymbol{u}, p ) &= (U_B (y_{c}),R^{-1} V_B (y_{c}), W_B (y_{c}),P_B ) \nonumber \\ & \quad + R^{-\frac {1}{3}} ( U_{\!B}^{\prime} (y_{c}), 0, W_{\!B}^{\prime} (y_c), 0 ) \tilde {Y} + \frac {1}{2} R^{-\frac {2}{3}} ( U_{\!B}^{\prime\prime} (y_{c}), 0, W_{\!B}^{\prime\prime} (y_c), 0 ) \tilde {Y}^{2} \nonumber \\ & \quad + R^{-1/3} ( \tilde {U}, R^{-1/3} \tilde {V}, \tilde {W}, R^{-1/3} \tilde {P} ), \end{align}

where

(2.52) \begin{align} \left( \tilde {U}, R^{-1/3} \tilde {V}, \tilde {W}, R^{-1/3} \tilde {P} \right) &= \sum _{n = -\infty }^{\infty } \left [ \left(\tilde {u}_{n}^{\langle 1 \rangle } (\tilde {Y}), R^{-1/3} b_{n}, \tilde {w}_{n}^{\langle 1 \rangle } (\tilde {Y}), R^{-1/3} \tilde {p}_{n}^{\langle 1 \rangle } \right)\right . \nonumber \\ &\quad \left . + R^{-1/3} \left(\tilde {u}_{n}^{\langle 2 \rangle } (\tilde {Y}), R^{-1/3} \tilde {v}_{n}^{\langle 2 \rangle } (\tilde {Y}), \tilde {w}_{n}^{\langle 2 \rangle } (\tilde {Y}), R^{-1/3} \tilde {p}_{n}^{\langle 2 \rangle } \right) \right ]\nonumber \\ &\quad \times \textrm{e}^{\textrm{i} n \alpha _w \zeta }. \end{align}

By matching with the main-layer solution (2.50), we derive

(2.53a,b) \begin{align} b_{0} &= 0, \quad \tilde {v}_{0}^{\langle 1 \rangle } = 0; \end{align}

(2.53c–e) \begin{align} \tilde {p}_{0}^{\langle 1 \rangle } &= 0, \quad \ \tilde {p}_{n}^{\langle 1 \rangle } = (\textrm{i} \alpha _w) / (n k_w^2) U_c^{\prime} b_n \ (n \ne 0); \qquad \tilde {p}_{n}^{\langle 2 \rangle } \equiv 0. \end{align}

Substituting (2.51)–(2.52) into the N–S (2.3), we obtain, at leading order, the governing equations for each Fourier component $(n \ne 0)$ of the regularised disturbance,

(2.54a) \begin{align} \alpha _w \tilde {u}_n^{\langle 1 \rangle } + \beta _w \tilde {w}_n^{\langle 1 \rangle } &= 0, \end{align}

(2.54b) \begin{align} \textrm{i} n \alpha _w U_c^{\prime} \tilde {Y} \tilde {u}_n^{\langle 1 \rangle } + b_{n} U_{\!B}^{\prime} (y_c) + \textrm{i} n \alpha _{w} \tilde {p}_{n}^{\langle 1 \rangle } - \tilde {u}_{n,\tilde {Y} \tilde {Y}}^{\langle 1 \rangle } &= -\tilde {\mathcal{N}}_{n}^{ \left \langle 1 x \right \rangle }, \end{align}

(2.54c) \begin{align} \tilde {p}_{n,\tilde {Y}}^{\langle 1 \rangle } &= 0, \end{align}

(2.54d) \begin{align} \textrm{i} n \alpha _w U_c^{\prime} \tilde {Y} \tilde {w}_n^{\langle 1 \rangle } + b_{n} W_{B}^{\prime} (y_c) + \textrm{i} n \beta _{w} \tilde {p}_{n}^{\langle 1 \rangle } - \tilde {w}_{n,\tilde {Y} \tilde {Y}}^{\langle 1 \rangle } &= -\tilde {\mathcal{N}}_{n}^{ \left \langle 1 z \right \rangle }, \end{align}

where

(2.55a) \begin{align} \tilde {\mathcal{N}}_{n}^{ \left \langle 1 x \right \rangle } = & \sum _{\substack {k=-\infty \\ k \ne 0, n }}^{\infty }b_{n-k} \tilde {u}_{k,\tilde {Y}}^{\langle 1 \rangle } + b_n \tilde {u}_{0,\tilde {Y}}^{\langle 1 \rangle } + \left (\tilde {u}_0^{\langle 1 \rangle } + \frac {\beta _w}{\alpha _w} \tilde {w}_0^{\langle 1 \rangle } \right ) \textrm{i} n \alpha _w \tilde {u}_n^{\langle 1 \rangle }, \end{align}

(2.55b) \begin{align} \tilde {\mathcal{N}}_{n}^{ \left \langle 1 z \right \rangle } = & \sum _{\substack {k=-\infty \\ k \ne 0, n }}^{\infty }b_{n-k} \tilde {w}_{k,\tilde {Y}}^{\langle 1 \rangle } + b_n \tilde {w}_{0,\tilde {Y}}^{\langle 1 \rangle } + \left (\tilde {u}_0^{\langle 1 \rangle } + \frac {\beta _w}{\alpha _w} \tilde {w}_0^{\langle 1 \rangle } \right ) \textrm{i} n \alpha _w \tilde {w}_n^{\langle 1 \rangle }. \end{align}

The spanwise momentum (2.54d ) is proportional to the chordwise momentum (2.54b ) according to the continuity (2.54a ). The solution for $\tilde {u}_n^{\langle 1 \rangle } (n \ne 0 )$ is determined by solving (2.54b ) numerically subject to the boundary conditions

(2.56) \begin{equation} \tilde {u}_n^{\langle 1 \rangle } \to \textrm{i} / (n \alpha _w) \left( U_{\!B}^{\prime}(y_c) / U_{c}^{\prime} - \alpha _w^2 k_{w}^{-2} \right) b_n / \tilde {Y} \qquad \textrm{as} \, \, \, \tilde {Y} \to \pm \infty , \end{equation}

which is derived by matching the CL solution with its counterpart in the main layer.

The governing equations for the leading-order mean-flow distortion $\tilde {u}_{0}^{\langle 1 \rangle } (\tilde {Y})$ and $\tilde {w}_{0}^{\langle 1 \rangle } (\tilde {Y})$ in the CL are

(2.57a,b) \begin{align} \tilde {u}_{0, \tilde {Y} \tilde {Y}}^{\langle 1 \rangle } =\sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {u}_{k, \tilde {Y}}^{\langle 1 \rangle }, \quad \tilde {w}_{0, \tilde {Y} \tilde {Y}}^{\langle 1 \rangle } = \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {w}_{k, \tilde {Y}}^{\langle 1 \rangle }. \end{align}

Integrating (2.57a ) from $-\infty$ to $\tilde {Y}$ gives

(2.58) \begin{equation} \tilde {u}_{0, \tilde {Y}}^{\langle 1 \rangle } (\tilde {Y}) = \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {u}_{k}^{\langle 1 \rangle } (\tilde {Y}) + \left . \tilde {u}_{0, \tilde {Y}}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} \to -\infty }, \end{equation}

and then integrating again from $0$ to $\tilde {Y}$ leads to

(2.59) \begin{equation} \tilde {u}_{0}^{\langle 1 \rangle } (\tilde {Y}) = \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \int _{0}^{\tilde {Y}} \tilde {u}_{k}^{\langle 1 \rangle } (Y_1) \textrm{d} Y_1 + \left . \tilde {u}_{0, \tilde {Y}}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} \to -\infty } \tilde {Y} + \left . \tilde {u}_{0}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} = 0}. \end{equation}

The first integral term in (2.59) approaches a constant as $\tilde {Y} \to \infty$ , which can be seen by substituting into it the matching condition (2.56). To prevent the non-physical $\mathit{O} ( \tilde {Y} )$ growth as $\tilde {Y} \to \infty$ in (2.59), we need to set

(2.60) \begin{equation} \left . \tilde {u}_{0, \tilde {Y}}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} \to -\infty } = 0. \end{equation}

It follows from (2.59) that there is a velocity jump in the leading-order mean-flow distortion

(2.61) \begin{equation} \left . \tilde {u}_{0}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} \to \infty } - \left . \tilde {u}_{0}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} \to -\infty } = \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} -\!\!\!\!\!\int _{-\infty }^{\infty } \tilde {u}_{k}^{\langle 1 \rangle } (Y_1) \textrm{d} Y_1, \end{equation}

where the integral must be evaluated as a Cauchy principal value. Through matching the solutions in the CL and main layer, we have the streaming jump across the CL,

(2.62) \begin{equation} U_{s}^{\langle 1 \rangle } (y_c^+) - U_{s}^{\langle 1 \rangle } (y_c^-) = \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} -\!\!\!\!\!\int _{-\infty }^{\infty } \tilde {u}_{k}^{\langle 1 \rangle } (Y_1) \textrm{d} Y_1. \end{equation}

Similarly, we can derive

(2.63) \begin{equation} \tilde {w}_{0, \tilde {Y}}^{\langle 1 \rangle } (\tilde {Y}) = \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {w}_{k}^{\langle 1 \rangle } (\tilde {Y}) \end{equation}

and the spanwise streaming jump. It follows from (2.54a ), (2.58) and (2.63) that

(2.64) \begin{equation} \tilde {u}_{0,\tilde {Y}}^{\langle 1 \rangle } (\tilde {Y}) + (\beta _{w} / \alpha _{w}) \tilde {w}_{0,\tilde {Y}}^{\langle 1 \rangle } (\tilde {Y}) = 0, \end{equation}

indicating that $\tilde {u}_{0}^{\langle 1 \rangle } + (\beta _w / \alpha _w) \tilde {w}_{0}^{\langle 1 \rangle }$ should be a constant, which we denote as

(2.65) \begin{equation} U_{s \textit{cef}}^{\langle 1 \rangle } \equiv \tilde {u}_{0}^{\langle 1 \rangle } + (\beta _{w} / \alpha _{w}) \tilde {w}_{0}^{\langle 1 \rangle }. \end{equation}

With the use of (2.65) in (2.55a ), the latter is expressed as

(2.66) \begin{equation} \tilde {\mathcal{N}}_{n}^{ \left \langle 1 x \right \rangle } = \sum _{\substack {k=-\infty \\ k \ne 0, n }}^{\infty }b_{n-k} \tilde {u}_{k,\tilde {Y}}^{\langle 1 \rangle } + b_n \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {u}_{k}^{\langle 1 \rangle } + U_{s \textit{cef}}^{\langle 1 \rangle } \textrm{i} n \alpha _w \tilde {u}_n^{\langle 1 \rangle }. \end{equation}

The governing equations for each Fourier component $(n \ne 0)$ of the second-order terms in expansion (2.52) are found to be

(2.67a) \begin{align} \textrm{i} n \alpha _w \tilde {u}_n^{\langle 2 \rangle } + \tilde {v}_{n, \tilde {Y}}^{\langle 2 \rangle } + \textrm{i} n \beta _w \tilde {w}_n^{\langle 2 \rangle } = & 0, \end{align}

(2.67b) \begin{align} \textrm{i} n \alpha _w U_c^{\prime} \tilde {Y} \tilde {u}_n^{\langle 2 \rangle } + \frac {1}{2} \textrm{i} n \alpha _w U_{c}^{\prime\prime} \tilde {Y}^2 \tilde {u}_n^{\langle 1 \rangle } + \tilde {v}_n^{\langle 2 \rangle } U_{\!B}^{\prime}(y_c) \qquad \qquad & \nonumber \\ + V_B(y_c) \tilde {u}_{n,\tilde {Y}}^{\langle 1 \rangle } + b_{n} U_{\!B}^{\prime\prime}(y_c) \tilde {Y} - \tilde {u}_{n,\tilde {Y} \tilde {Y}}^{\langle 2 \rangle } = & - \tilde {\mathcal{N}}_{n}^{ \left \langle 2 x \right \rangle }, \end{align}

(2.67c) \begin{align} \textrm{i} n \alpha _w U_c^{\prime} \tilde {Y} \tilde {w}_n^{\langle 2 \rangle } + \frac {1}{2} \textrm{i} n \alpha _w U_{c}^{\prime\prime} \tilde {Y}^2 \tilde {w}_n^{\langle 1 \rangle } + \tilde {v}_n^{\langle 2 \rangle } W_{B}^{\prime}(y_c) \qquad \qquad & \nonumber \\ + V_B(y_c) \tilde {w}_{n,\tilde {Y}}^{\langle 1 \rangle } + b_{n} W_{B}^{\prime\prime}(y_c) \tilde {Y} - \tilde {w}_{n,\tilde {Y} \tilde {Y}}^{\langle 2 \rangle } = & - \tilde {\mathcal{N}}_{n}^{ \left \langle 2 z \right \rangle }, \end{align}

where use has been made of (2.53e ) and

(2.68a) \begin{align} \tilde {\mathcal{N}}_{n}^{ \left \langle 2 x \right \rangle } &= \sum _{\substack {k=-\infty \\ k \ne 0, n }}^{\infty } \left [\left (\tilde {u}_{n-k}^{\langle 2 \rangle } + \frac {\beta _w}{\alpha _w} \tilde {w}_{n-k}^{\langle 2 \rangle } \right ) \textrm{i} k \alpha _w \tilde {u}_{k}^{\langle 1 \rangle } + b_{n-k} \tilde {u}_{k,\tilde {Y}}^{\langle 2 \rangle } + \tilde {v}_{n-k}^{\langle 2 \rangle } \tilde {u}_{k,\tilde {Y}}^{\langle 1 \rangle } \right ] + U_{s \textit{cef}}^{\langle 1 \rangle } \textrm{i} n \alpha _w \tilde {u}_{n}^{ \langle 2 \rangle } \nonumber \\ &\quad + \tilde {u}_{0\textit{eff}}^{\langle 2 \rangle } \textrm{i} n \alpha _w \tilde {u}_{n}^{ \langle 1 \rangle } + b_n \tilde {u}_{0,\tilde {Y}}^{\langle 2 \rangle } + \tilde {v}_n^{\langle 2 \rangle } \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } b_{-k} \tilde {u}_{k}^{\langle 1 \rangle }, \end{align}

(2.68b) \begin{align} \tilde {\mathcal{N}}_{n}^{ \left \langle 2 z \right \rangle } &= \sum _{\substack {k=-\infty \\ k \ne 0, n }}^{\infty } \left [\left (\tilde {u}_{n-k}^{\langle 2 \rangle } + \frac {\beta _w}{\alpha _w} \tilde {w}_{n-k}^{\langle 2 \rangle } \right ) \textrm{i} k \alpha _w \tilde {w}_{k}^{\langle 1 \rangle } + b_{n-k} \tilde {w}_{k,\tilde {Y}}^{\langle 2 \rangle } + \tilde {v}_{n-k}^{\langle 2 \rangle } \tilde {w}_{k,\tilde {Y}}^{\langle 1 \rangle } \right ] + U_{s \textit{cef}}^{\langle 1 \rangle } \textrm{i} n \alpha _w \tilde {w}_{n}^{ \langle 2 \rangle } \nonumber \\ &\quad + \tilde {u}_{0\textit{eff}}^{\langle 2 \rangle } \textrm{i} n \alpha _w \tilde {w}_{n}^{ \langle 1 \rangle } + b_n \tilde {w}_{0,\tilde {Y}}^{\langle 2 \rangle } + \tilde {v}_n^{\langle 2 \rangle } \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } b_{-k} \tilde {w}_{k}^{\langle 1 \rangle } \end{align}

with $\tilde {u}_{0\textit{eff}}^{\langle 2 \rangle } = \tilde {u}_{0}^{\langle 2 \rangle } + (\beta _{w} / \alpha _{w}) \tilde {w}_{0}^{\langle 2 \rangle }$ . Combining (2.67b ) and (2.67c ), and making use of (2.67a ) followed by differentiating the resulting equation with respect to $\tilde {Y}$ , we obtain

(2.69) \begin{equation} Q_{n,\tilde {Y} \tilde {Y}} - \textrm{i} n \alpha _w U_{c}^{\prime} \left ( \tilde {Y} + \frac {U_{s \textit{cef}}^{\langle 1 \rangle }}{U_{c}^{\prime}} \right ) Q_{n} + \textrm{i} n \alpha _w U_c^{\prime\prime} b_{n} = n \sum _{\substack {k = -\infty \\ k \ne 0, n}}^{\infty } \frac {1}{k} b_{n-k} Q_{k, \tilde {Y}} - \textrm{i} n \alpha _w b_n \tilde {u}_{0\textit{eff},\tilde {Y} \tilde {Y}}, \end{equation}

where $Q_{n} = \tilde {v}_{n, \tilde {Y} \tilde {Y}}^{\langle 2 \rangle }$ .

The mean-flow distortions at $\mathit{O} (R^{-2/3})$ are governed by

(2.70a) \begin{align} \tilde {u}_{0, \tilde {Y} \tilde {Y}}^{\langle 2 \rangle } = & V_{B}(y_c) \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {u}_{k}^{\langle 1 \rangle } + \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } \left ( \tilde {v}_{-k}^{\langle 2 \rangle } \tilde {u}_{k}^{\langle 1 \rangle } + b_{- k} \tilde {u}_{k}^{\langle 2 \rangle } \right )_{\tilde {Y}}, \end{align}

(2.70b) \begin{align} \tilde {w}_{0, \tilde {Y} \tilde {Y}}^{\langle 2 \rangle } = & V_{B}(y_c) \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {w}_{k}^{\langle 1 \rangle } + \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } \left ( \tilde {v}_{-k}^{\langle 2 \rangle } \tilde {w}_{k}^{\langle 1 \rangle } + b_{- k} \tilde {w}_{k}^{\langle 2 \rangle } \right )_{\tilde {Y}}, \end{align}

which can be integrated from $- \infty$ to $\infty$ to show that

(2.71a) \begin{align} \left . \tilde {u}_{0, \tilde {Y}}^{\langle 2 \rangle } \right \rvert _{-\infty }^{\infty } = & V_{B}(y_c) \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} -\!\!\!\!\!\int _{-\infty }^{\infty } \tilde {u}_{k}^{\langle 1 \rangle } (Y_1) \textrm{d} Y_1 + \frac {\textrm{i} \alpha _w}{k_w^2} \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } \frac {1}{k} \left ( a_{k}^{+} - a_{k}^{-} \right ) b_{-k}, \end{align}

(2.71b) \begin{align} \left . \tilde {w}_{0, \tilde {Y}}^{\langle 2 \rangle } \right \rvert _{-\infty }^{\infty } = & V_{B}(y_c) \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} -\!\!\!\!\!\int _{-\infty }^{\infty } \tilde {w}_{k}^{\langle 1 \rangle } (Y_1) \textrm{d} Y_1 + \frac {\textrm{i} \beta _w}{k_w^2} \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } \frac {1}{k} \left ( a_{k}^{+} - a_{k}^{-} \right ) b_{-k}, \quad \end{align}

where use has been made of (2.50). Equation (2.71) are the gradient jumps of the mean-flow distortion at $\mathit{O} (R^{-2/3})$ . We combine (2.70a ) and (2.70b ) using the continuity (2.54a ) and (2.67a ) to obtain

(2.72) \begin{equation} \tilde {u}_{0\textit{eff},\tilde {Y} \tilde {Y}} = \sum _{\substack {k = - \infty k \ne 0}}^{\infty } \frac {\textrm{i}}{k \alpha _w} b_{-k} Q_{k}, \end{equation}

which can be inserted into the right-hand side of (2.69), giving

(2.73) \begin{equation} Q_{n,\tilde {Y} \tilde {Y}} - \textrm{i} n \alpha _w U_{c}^{\prime} \left ( \tilde {Y} + \frac {U_{s \textit{eff},c }^{\langle 1 \rangle }}{U_{c}^{\prime}} \right ) Q_{n} + \textrm{i} n \alpha _w U_c^{\prime\prime} b_{n} = n \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } \frac {1}{k} \bigl ( b_{n-k} Q_{k, \tilde {Y}} + b_n b_{-k} Q_{k} \bigr ). \qquad \end{equation}

To eliminate the effect of the steady streaming $U_{s \textit{cef}}^{\langle 1 \rangle }$ from the calculations of the CL, we perform a coordinate shift in the transverse direction

(2.74) \begin{equation} \mathcal{Y} = \tilde {Y} + U_{s \textit{cef}}^{\langle 1 \rangle } / U_{c}^{\prime}, \end{equation}

and rewrite (2.73) as

(2.75) \begin{equation} Q_{n,\mathcal{Y} \mathcal{Y}} - \textrm{i} n \alpha _w U_{c}^{\prime} \mathcal{Y} Q_{n} + \textrm{i} n \alpha _w U_c^{\prime \prime } b_{n} = n \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } \frac {1}{k} \bigl ( b_{n-k} Q_{k, \mathcal{Y}} + b_n b_{-k} Q_{k} \bigr ). \end{equation}

The above system of coupled nonlinear equations is to be solved subject to the far-field boundary conditions

(2.76) \begin{equation} Q_{n} \to \big(U_{c}^{\prime \prime } / U_{c}^{\prime }\big) b_{n} / \mathcal{Y} \quad \textrm{as} \, \, \, \mathcal{Y} \to \pm \infty . \end{equation}

We denote the velocity jumps across the CL as

(2.77) \begin{equation} J_{n} = a_{n}^{+} - a_{n}^{-} \quad (n \ne 0), \end{equation}

and matching $Q_n$ with the solution in the main layer gives

(2.78) \begin{equation} J_{n} = \,\kern-1.5pt\hbox{--}\kern0.5pt\!\!\!\!\!\int _{- \infty }^{\infty } Q_{n} \left ( Y_1 \right ) \textrm{d} Y_1. \end{equation}

The governing equations and boundary conditions for $Q_{n}$ (2.75)–(2.76) and the matching condition (2.78) along with the Rayleigh (2.44) and the Frobenius solution near the critical level (2.46) are used simultaneously to resolve the coupled dynamics between the main layer and CL.

The CL is strongly nonlinear as (2.54b ) and (2.73) indicate. All harmonics of the roughness-induced disturbance are generated at the same order. Importantly, due to three-dimensionality of the disturbance and the algebraic singularity of their chordwise and spanwise velocities, the regularised disturbance velocity in the CL has a larger amplitude of $\mathit{O} (R^{-1/3})$ than its $\mathit{O}(R^{-2/3})$ counterpart in the main layer. Furthermore, the associated vorticity in the CL is even larger with an $\mathit{O} (1)$ amplitude. This can give rise to high-frequency short-wavelength secondary instability, as will be shown in § 4.

To summarise, the WL converts surface displacement into a blowing velocity of $\mathit{O} (R^{-2/3})$ , which drives a forced disturbance in the main layer. The forced disturbance is governed by the steady Rayleigh equation with a singularity at the critical level $y = y_{c}$ , which is regularised by introducing viscous and nonlinear effects in the CL. The disturbance inside the CL has a larger, $\mathit{O} (R^{-1/3})$ , amplitude, and leads to $\mathit{O} (1)$ vorticities. Physically, the effect of such surface roughness is significant near the wall and around the critical level $y_{c}$ , where the flows are nonlinearly distorted.

3. Numerical results for the roughness-distorted flow

3.1. Numerical methods

For the roughness-induced disturbance in the WL, we solve (2.28) using the second-order finite-difference scheme to calculate $n \ne 0$ Fourier components. Only half of the Fourier components $(n \gt 0)$ need to be computed with the other half $(n \lt 0)$ being given by (2.27). The $n = 0$ Fourier component is calculated by integrating (2.35b ) using the trapezoidal formula. The roughness interaction and nonlinear terms $\{ \mathcal{R}_{n}^{x}, \mathcal{R}_{n}^{z},\mathcal{N}_{n}^{x},\mathcal{N}_{n}^{z} \}$ on the right-hand side of (2.28) are evaluated using a pseudospectral method, which transforms the quantities in the spectral space to the physical space to carry out the multiplications and then transforms the products back to the spectral space. The $3/2$ rule (Canuto 1988) is applied to eliminate the aliasing error. An iteration process is needed for these nonlinear calculations. We initially simplify (2.28) by setting the roughness interaction and nonlinear terms on the right-hand side to zero, which enables us to compute numerically the linear solution or to derive the analytical solution. The linear solution then serves as the initial guess for starting the iteration process solving the nonlinear problem for a small roughness height $h$ . The iteration continues until the criterion

(3.1) \begin{equation} \| \boldsymbol{\varPhi } - \boldsymbol{\varPhi }_{\textit{old}} \|_{2} / \| \boldsymbol{\varPhi } \|_{2} \le 10^{-7} \end{equation}

is satisfied, where $\boldsymbol{\varPhi }$ is the solution matrix consisting of all Fourier components of the velocity $\{ \bar {U}^{\dagger },\bar {V}^{\dagger },\bar {W}^{\dagger } \}$ . Subsequently, we employ a parameter continuation method, incrementally increasing the roughness height and using the convergent solution as the initial guess to advance the calculation for larger $h$ . To enhance the convergence, we use the under-relaxation technique (Xu, Zhang & Wu 2017).

For the main-layer dynamics, our primary interest is in the forced disturbance driven by the blowing velocity (2.40) from the WL motion. The forced disturbance is governed by the steady Rayleigh equation (2.44) subject to inhomogeneous boundary conditions (2.45). Apart from the singularity at $y = y_c$ , the Rayleigh equation is also singular at the lower boundary $y = 0$ because the effective base flow satisfies the no-slip boundary condition while $U_{\!B \textit{eff}, yy} (0) \ne 0$ . We need to construct a Frobenius solution to (2.44) at $y = 0$ similar to that near $y = y_c$ , namely

(3.2) \begin{equation} v_{wn} = a_{n}^{\left \langle w \right \rangle } \phi _{an}^{\left \langle w \right \rangle } (y) + b_{n}^{\left \langle w \right \rangle } \Bigl [ \phi _{b n}^{\left \langle w \right \rangle }(y) + \frac {U_{w}^{\prime\prime}}{U_{w}^{\prime}} \ln { \left |y \right | } \phi _{an}^{\left \langle w \right \rangle }(y) \Bigr ], \end{equation}

where

(3.3a,b) \begin{gather} \phi _{an}^{\left \langle w \right \rangle } = y + \frac {1}{2} \frac {U_w^{\prime \prime }}{U_w^{\prime }} y ^ 2 + {\cdots} , \quad \phi _{bn}^{\left \langle w \right \rangle } = 1 + \frac {1}{2} \left (n^2 k_w^2 + \frac {U_w^{\prime \prime \prime }}{U_w^{\prime }} - 2 \frac {U_w^{\prime \prime 2}}{U_w^{\prime \prime 2}} \right )y^2 + {\cdots} , \end{gather}

with the subscript $w$ of $U_{w}$ signifying the quantities related to $U_{\!B \textit{eff}}$ and its derivatives evaluated at $y = 0$ . Imposing the lower boundary condition (2.45a ) on (3.2) gives

(3.4) \begin{equation} b_{n}^{\left \langle w \right \rangle } = V_{w \infty , n}. \end{equation}

We start by calculating the linear solution $v_{n}$ in the main layer. The solution to the steady Rayleigh equation (2.44) has the far-field asymptote, $v_{n} \to C_{n}^{\left \langle \infty \right \rangle } \textrm{exp}(- n k_w) y$ as $y \to \infty$ , where $C_{n}^{\left \langle \infty \right \rangle }$ is an unknown constant. We first use the MATLAB solver ODE15s to integrate (2.44) downwards from a large $y$ position, denoted as $y_{\infty }$ , to $y_{c} + \delta y$ with the initial conditions

(3.5a,b) \begin{gather} \check {v}_{n} = \exp \left (- n k_{w} y_{\infty } \right )\!, \quad \check {v}_{n}' = - n k_{w} \exp \left (- n k_{w} y_{\infty } \right )\!. \end{gather}

The resulting numerical solution at $y_{c} + \delta y$ , denoted as $\{ \check {v}_{n} (y_{c} + \delta y), \check {v}_{n}' (y_{c} + \delta y) \}$ , is equated to, after being multiplied by $C_{n}^{\left \langle \infty \right \rangle }$ , the local Frobenius solution (2.46) at $y = y_{c} + \delta y$ ,

(3.6a) \begin{align} C_{n}^{\left \langle \infty \right \rangle } \check {v}_{n} (y_{c} + \delta y) = a_{n}^{+} \phi _{an} (y_{c} + \delta y) + b_{n} \varphi _{bn} (y_{c} + \delta y), \end{align}

(3.6b) \begin{align} C_{n}^{\left \langle \infty \right \rangle } \check {v}_{n}' (y_{c} + \delta y) = a_{n}^{+} \phi _{an}^{\prime} (y_{c} + \delta y) + b_{n} \varphi _{bn}^{\prime} (y_{c} + \delta y), \end{align}

where $\varphi _{bn} (y) = \phi _{bn} (y) + (U_{c}^{\prime\prime} / U_{c}^{\prime}) \ln |y - y_{c}| \phi _{an} (y)$ . We also integrate (2.44) downwards from $y_{c} - \delta y$ to $\delta y$ with the initial conditions,

(3.7a,b) \begin{gather} \hat {v}_{an} = \phi _{an} (y_c - \delta y), \quad \hat {v}_{an}^{\prime} = \phi _{an}^{\prime} (y_c - \delta y), \end{gather}

and

(3.8a,b) \begin{gather} \hat {v}_{bn} = \varphi _{bn} (y_c - \delta y), \quad \hat {v}_{bn}^{\prime} = \varphi _{bn}^{\prime} (y_c - \delta y), \end{gather}

respectively, and the numerical results at $y = \delta y$ are denoted as $\{ \hat {v}_{an} (\delta y), \hat {v}_{an}^{\prime} (\delta y) \}$ and $\{ \hat {v}_{bn} (\delta y), \hat {v}_{bn}^{\prime} (\delta y) \}$ . The linear combination of the above two, $a_{n}^{-} \hat {v}_{an} + b_{n} \hat {v}_{bn}$ , gives the total solution. Equating it and its derivative to the local Frobenius solution (3.2) at $y = \delta y$ , we obtain

(3.9a) \begin{align} a_{n}^{-} \hat {v}_{an} (\delta y) + b_{n} \hat {v}_{bn} (\delta y) = a_{n}^{\left \langle w \right \rangle } \phi _{an}^{\left \langle w \right \rangle } (\delta y) + b_{n}^{\left \langle w \right \rangle } \varphi _{wn} (\delta y), \end{align}

(3.9b) \begin{align} a_{n}^{-} \hat {v}_{an}^{\prime} (\delta y) + b_{n} \hat {v}_{bn}^{\prime} (\delta y) = a_{n}^{\left \langle w \right \rangle } \phi _{an}^{\left \langle w \right \rangle }{'} (\delta y) + b_{n}^{\left \langle w \right \rangle } \varphi _{wn}^{\prime} (\delta y), \end{align}

where $\varphi _{wn} (y) = \phi _{bn}^{\left \langle w \right \rangle } ( y) + (U_{w}^{\prime\prime} / U_{w}^{\prime}) \ln | y | \phi _{an}^{\left \langle w \right \rangle } ( y)$ . Using (3.6) and (3.9) as well as the linear velocity jump condition across the CL,

(3.10) \begin{equation} a_{n}^{+} - a_{n}^{-} = \mathcal{K}_{c} \pi b_{n}, \end{equation}

where $\mathcal{K}_{c} = \textrm{i} (U_{c}^{\prime\prime} / U_{c}^{\prime}) \textrm{sign} (\alpha _{w} U_{c}^{\prime})$ , we can form a linear system to solve the coefficients $\{ a_n^{\pm },b_n,a_n^{\left \langle w \right \rangle } \}$ in the Frobenius solutions (2.46) and (3.2) as well as the coefficient in the large- $y$ asymptote, $C_{n}^{\left \langle \infty \right \rangle }$ , for different modes $n$ ( $n \ne 0$ ),

(3.11) \begin{equation} \boldsymbol{L}_{n} \boldsymbol{x}_{n} = \boldsymbol{B}_{n}, \end{equation}

where

(3.12a) \begin{align} \boldsymbol{L}_{n} & = \begin{pmatrix} \phi _{an}(y_c + \delta y) & \varphi _{bn} (y_{c} + \delta y) & 0 & 0 & - \check {v}_n(y_c + \delta y) \\ \phi _{an}^{\prime}(y_c + \delta y) & \varphi _{bn}^{\prime} (y_{c} + \delta y) & 0 & 0 & - \check {v}_{n}' (y_c + \delta y) \\ 1 & - \mathcal{K}_{c} \pi & -1 & 0 & 0 \\ 0 & \hat {v}_{b n} (\delta y) & \hat {v}_{an} (\delta y) & -\phi _{an}^{\left \langle w \right \rangle } (\delta y) & 0 \\ 0 & \hat {v}_{b n}' (\delta y) & \hat {v}_{an}^{\prime} (\delta y) & - \phi _{an}^{\left \langle w \right \rangle }{'} (\delta y) & 0 \end{pmatrix}\!, \qquad \end{align}

(3.12b) \begin{align} \boldsymbol{x}_{n} &= {\begin{pmatrix} a_n^+ & b_n & a_n^- & a_{n}^{\left \langle w \right \rangle } & C_{n}^{\left \langle \infty \right \rangle } \end{pmatrix}}^{\textrm{T}}, \end{align}

(3.12c) \begin{align} \boldsymbol{B}_{n} &= \begin{pmatrix} 0 & 0 & 0 & b_{n}^{\left \langle w \right \rangle } \varphi _{wn} (\delta y) & b_{n}^{\left \langle w \right \rangle } \varphi _{wn}^{\prime} (\delta y) \end{pmatrix}^{\textrm{T}}. \end{align}

After obtaining the linear main-layer solution, we can calculate the linear solution for $Q_n$ in the CL by solving (2.75) with the right-hand side terms being set to $0$ . The derivative in (2.75) is approximated using the second-order finite-difference scheme. The linear solution for $Q_{n}$ can also be found analytically, but the detail is neglected here.

The linear results for $\{ b_{n}, Q_{n} \}$ are used as the initial guesses for nonlinear calculations for a small roughness height. Similar to solving the nonlinear equations in the WL, we adopt a parameter continuation method to incrementally increase the roughness height, using the convergent result for the previous slightly smaller height as the initial guesses to start iterations. The nonlinear terms are evaluated using the pseudospectral method. After obtaining $Q_{n}$ , we can calculate the velocity jump $J_{n}$ according to (2.78). With the nonlinear jump $J_n$ , we use a modified linear system compared with (3.11) to update $\boldsymbol{x}_{n}$ ,

(3.13) \begin{equation} \boldsymbol{L}_{n} \boldsymbol{x}_{n} = \boldsymbol{B}_{n}^{\left \langle N \right \rangle }, \end{equation}

where

(3.14) \begin{equation} \boldsymbol{B}_{n}^{\left \langle N \right \rangle } = \begin{pmatrix} 0 & 0 & J_{n} - \mathcal{K}_{c} \pi b_{n,0} & b_{n}^{\left \langle w \right \rangle } \varphi _{wn} (\delta y) & b_{n}^{\left \langle w \right \rangle } \varphi _{wn}^{\prime} (\delta y) \end{pmatrix}^{\textrm{T}}, \end{equation}

and $b_{n,0}$ is the value of $b_{n}$ calculated in the previous iteration. The blowing velocity $V_{w\infty ,n}$ for the current roughness height is applied on the last two entries of (3.14) according to (3.4). We implement the iteration for the above calculations to resolve the coupled CL and main-layer dynamics until the following criteria are satisfied:

(3.15a,b) \begin{gather} \| \boldsymbol{b} - \boldsymbol{b}_{\textit{old}} \|_2 / \| \boldsymbol{b} \|_2 \le 10^{-7} \quad \textrm{and} \quad \| \boldsymbol{Q} - \boldsymbol{Q}_{\textit{old}} \|_2 / \| \boldsymbol{Q} \|_2 \le 10^{-7}, \end{gather}

where $\boldsymbol{b} = \{ b_{n} \}$ and $\boldsymbol{Q} = \{ Q_{n} \}$ are the solution vector and matrix consisting of all Fourier components, respectively. We then calculate $\tilde {u}_{n}^{\langle 1 \rangle }$ in the CL by solving numerically (2.54b ), and further obtain $\tilde {w}_{n}^{\langle 1 \rangle }$ using the continuity (2.54a ).

3.2. Numerical results for disturbances in the WL

We consider the simplest roughness configuration with $F_{1} = F_{-1} = 1$ and $F_{n} = 0 \ (n \ne \pm 1)$ in (2.2) for our numerical calculations. A diagram of the surface roughness is shown in figure 2. Two pairs of roughness wavenumbers, Case I with $(\alpha _{w},\beta _{w}) = (-0.525,0.606)$ and Case II with $(\alpha _{w},\beta _{w}) = (-0.262,0.303)$ , are chosen for numerical calculations. They were chosen based on simple guiding principles: they are more-or-less of $\mathit{O} (1)$ as assumed in the theory, and fall within the wavenumber band where stationary cross-flow vortices arise. The first pair was somewhat close to that of the neutral mode of the stationary cross-flow vortices ( $(\alpha , \beta ) = (-0.6376,0.7621)$ ), while the second took the half-values of the first. We present numerical results for Case I here, while those for Case II are shown in the supplementary material (§ S2). For the base FSC flow, we set the spanwise slip velocity $W_{\infty } = 1$ and the acceleration parameter $m = 0.5$ , as it characterises the portion of typical swept-wing boundary-layer flows where a strong favourable pressure gradient is present.

Figure 2. Diagram of the surface roughness used in calculations.

We start our numerical calculations in the WL from a small roughness height ( $h = 0.1$ ) and gradually increase it. In the wall-normal direction, the domain is $\bar {Y} \in [0,30]$ with the normal grid size being set as $\Delta \bar {Y} = 0.05$ , and 17 Fourier modes $(-8 \le n \le 8)$ are included. These discretisation choices were chosen based on resolution checks, with due consideration of managing computational cost. This adequacy was supported by several observations: the results near the truncated computational domain matched the corresponding asymptotes; the current results agree with those calculated using a halved mesh size; the magnitudes of the highest Fourier mode ( $n = \pm 8$ ) remain below the iteration stopping criteria ( $10^{-7}$ ). Our numerical codes can converge up to $h = 0.825$ for Case I roughness. Figures 3 and 4 show the profiles of several harmonics of the chordwise and wall-normal velocities of the disturbance in the WL calculated for different roughness heights. The profiles of the spanwise velocities are similar to those of the chordwise velocities, and hence are omitted here. For the fundamental mode $n = 1$ , the profiles for different roughness heights, normalised by $h$ (see (2.25)), almost overlap. This is due to the fact that the flow is primarily disturbed by the fundamental mode of the roughness shape function $F$ . For higher harmonics of the velocity components, their profiles also nearly overlap when appropriately normalised, i.e. by dividing $\bar {U}_{2}$ here by $h$ and $\bar {U}_{3}$ by $h^{2}$ . Figure 4 shows that the Fourier components of the wall-normal velocity approach different constants as $\bar {Y} \to \infty$ . These constants (or equivalently, the combination of the chordwise and spanwise velocity components through the continuity (2.28a )) are used to calculate the transpiration velocity $\bar {V}_t$ , and furthermore the blowing velocity $V_w |_{\bar {Y} \to \infty }$ , which drives the $\mathit{O} (R^{- 2/3})$ disturbances in the main layer.

Figure 3. Profiles of the chordwise velocity harmonics $\bar {U}_{n}$ in the WL for different roughness height.

Figure 4. Profiles of the wall-normal velocity harmonics $\bar {V}_{n}$ in the WL for different roughness height.

Figure 5 shows the mean-flow distortions (2.35b ) in the WL. They are caused by nonlinear interactions of the disturbance, and approach constants as $\bar {Y} \to \infty$ . Such streaming in the nonlinear WL drives the $\mathit{O} (R^{-1/3})$ mean-flow distortion in the main layer. Note that the chordwise velocity is approximately four times as large as the spanwise velocity.

Figure 5. Profiles of the mean-flow distortions in the WL for different roughness height.

The $\mathit{O} (R^{-1/3})$ disturbance velocity components, along with the base flow, correspond to $\mathit{O} (1)$ vorticities in the WL, namely

(3.16) \begin{equation} \left(\partial _{x}, R^{1/3} \partial _{\bar {Y}}, \partial _{z}\right) \times \left( R^{-1/3} \bar {U}_{w}, R^{-2/3} \bar {V}_{w}, R^{-1/3} \bar {W}_{w}\right) = \left(\bar {W}_{w, \bar {Y}},0, -\bar {U}_{w, \bar {Y}}\right) + {\cdots} . \end{equation}

To aid the study of the secondary instability of the distorted WL to be conducted in § 4, we rotate the $x$ – $z$ plane of the original coordinate system so that the $x$ -axis aligns with the direction of $(\bar {\alpha },\bar {\beta })$ , where $\{ \bar {\alpha }, \bar {\beta } \}$ , introduced in (4.16), are wavenumbers of the secondary instability mode. The chordwise and spanwise velocity components, $\{ \bar {U}_{w}, \bar {W}_{w} \}$ , are projected onto the new coordinate axes, as $\{ \bar {\mathfrak{U}}_{w}, \bar {\mathfrak{W}}_{w} \}$ ,

(3.17) \begin{equation} \begin{pmatrix} \bar {\mathfrak{U}}_{w} \\ \bar {\mathfrak{W}}_{w} \end{pmatrix} = \begin{pmatrix} \cos \bar {\theta } & \sin \bar {\theta } \\ - \sin \bar {\theta } & \cos \bar {\theta } \end{pmatrix} \begin{pmatrix} \bar {U}_{w} \\ \bar {W}_{w} \end{pmatrix} \end{equation}

with $\bar {\theta } = \arctan (\bar {\beta } / \bar {\alpha })$ . The vorticity along the direction of $(\bar {\alpha },\bar {\beta })$ can be calculated as

(3.18) \begin{equation} \bar {\varOmega } = \bar {\mathfrak{U}}_{w,\bar {Y}} + {\cdots} = - \left ( \bar {\alpha } \Big / \sqrt {\bar {\alpha }^{2} + \bar {\beta }^{2}} \right ) \Bigl [ ( \bar {\beta } / \bar {\alpha } ) \bar {U}_{w} - \bar {W}_{w} \Bigr ]_{\bar {Y}} + {\cdots} . \end{equation}

Figure 6 displays contours of this skewed vorticity $\bar {\varOmega }$ without the contributions from base-flow wall shears for two roughness heights. The overall features are similar, but the vorticity induced by higher roughness is stronger, and acquires a more or less $\mathit{O} (1)$ intensity, as expected. The vorticity in the direction perpendicular to $(\bar {\alpha },\bar {\beta })$ (the direction of $(-\bar {\beta }, \bar {\alpha })$ ) controls the secondary instability, and its contours for the roughness with $h = 0.825$ will be shown later in § 5.2.1, where instability results are presented.

Figure 6. Contours of the skewed WL vorticity $\bar {\varOmega }$ without the contributions from base-flow wall shears in Case I: $(a)$ $h = 0.2$ and $(b)$ $h = 0.825$ .

The blowing velocity from the WL is of importance as it drives the response in the main layer and CL. Figure 7 shows the modulus of the first three Fourier components in the blowing velocity versus the roughness height. The fundamental mode has a much more significant modulus than higher harmonics.

Figure 7. Modulus of Fourier components in the blowing velocity versus the roughness height.

3.3. Numerical results for disturbances in the main and critical layers

For the calculations pertaining to the main layer, we first need to ascertain the position of the critical level $y_{c}$ , which is $y_{c} = 2.22081$ . Then we calculate the linear solution to the Rayleigh equation (2.44). In practice, we solve (2.44) subject to unity blowing velocity, $v_{n} (0) = b_{n}^{\left \langle w \right \rangle } = 1$ . Figure 8 presents the profiles of the first three harmonics of $v_{n} / V_{w \infty , n}$ , the linear solutions normalised by the blowing velocity. The resulting solution, multiplied by corresponding Fourier components of the blowing velocity induced by the roughness with a small height, serves as initial guesses to start the nonlinear calculations for the coupled main-layer and CL dynamics. Numerically, the linear solutions in the main layer are obtained by solving the system (3.11).

Figure 8. The normalised solutions to Rayleigh equation (2.44).

Next, we consider the coupled main-layer and CL dynamics by solving the system (3.13) along with (2.75)–(2.78). The coupling is facilitated through the velocity jumps $J_{n}$ across the CL, as shown in (3.14). Since the mean-flow distortion in the main layer is not our primary concern, and its effect on the CL dynamics amounts to a simple coordinate shift (2.74), we set $U_{s \textit{cef}}^{\langle 1 \rangle } = 0$ in our calculations and revert the wall-normal coordinate from $\mathcal{Y}$ to $\tilde {Y}$ in the CL. The system governing CL dynamics is truncated in the region $-30 \le \tilde {Y} \le 30$ , with the mesh size $\Delta \tilde {Y} = 0.05$ . The Fourier series for $b_{n}$ and $Q_{n}$ consist of 17 harmonics, i.e. $-8 \le N \le 8$ . These numerical settings were determined with both accuracy and computational efficiency in mind. Our numerical codes for this coupled dynamics can converge up to $h = 0.43$ for Case I roughness. As with the WL calculations, the nonlinear numerical codes fail to converge beyond a certain threshold, possibly due to the limitation of the pseudospectrum-based iteration; for higher $h$ , full Newton iteration may be needed. The dynamics within the CL changes from nearly linear to strongly nonlinear as the height of surface roughness is increased.

In our calculations, the Fourier components of the roughness shape function are only non-zero for $n = \pm 1$ , and the wavenumbers of Case I roughness are fairly close to those of the neutral mode of the stationary cross-flow vortices, which results in particularly significant response in the $n = \pm 1$ Fourier components of the disturbance. The nonlinear interactions of the fundamental modes generate the second and higher harmonics. We monitor only the first two harmonics of the velocity jumps $J_{n}$ . To scale out the coefficient $\mathcal{K}_{c}$ and $b_{n}$ in (3.10), we introduce a normalised jump,

(3.19) \begin{equation} \tilde {J}_{n} = J_{n} / (\mathcal{K}_{c} b_{n}), \quad n = 1, 2. \end{equation}

Figure 9 illustrates the variations of $\tilde {J}_1$ and $\tilde {J}_{2}$ versus the roughness height $h$ . When $h$ is small, $\tilde {J}_1 \to \pi$ ; whereas for sufficiently large $h$ , $\textrm{Re} \{ \tilde {J}_1 \}$ approaches zero but $\textrm{Im} \{ \tilde {J}_1 \}$ does not, which indicates a changeover from a linear viscous CL to a nonlinear viscous CL. As shown in the zoomed-in portion of figure 9(a), $\tilde {J}_{2}$ approaches a constant near $\pi$ , but not exactly $\pi$ , when $h \to 0$ due to nonlinearity. The variation of the velocity jump with $h$ can be non-monotonic.

Figure 9. Normalised velocity jumps $\tilde {J}_{1}$ and $\tilde {J}_{2}$ across the CL.

Figure 10 shows the profiles of different harmonics of $Q_{n}$ in the CL. The fundamental and second harmonic have comparable amplitudes, while the third harmonic is appreciably weaker. Higher harmonics appear to be more oscillatory since it is more substantially affected by nonlinearity. After obtaining $Q_{n}$ and $b_{n}$ , we can evaluate the nonlinear terms on the right-hand side of (2.54b ) and calculate $\tilde {u}_n^{\langle 1 \rangle }$ $(n \ne 0)$ , the profiles of which are shown in figure 11. The spanwise velocity $\tilde {w}_n^{\langle 1 \rangle }$ $(n \ne 0)$ is proportional to $\tilde {u}_{n}^{\langle 1 \rangle }$ as indicated by (2.54a ). The fundamental has a much greater amplitude than the second and third harmonics. This is due to the fact that its wavenumbers are close to those of the neutral stationary cross-flow mode, which leads to a near resonance. The direct forcing of the second and third harmonics by the blowing velocity is relatively weak, while the excitation by nonlinear interactions of the fundamental is more substantial. This also explains their more oscillatory profiles.

Figure 10. Profiles of $Q_{n}$ in the CL for different roughness height.

Figure 11. Profiles of the leading-order chordwise velocity harmonics $\tilde {u}_{n}^{\langle 1 \rangle }$ in the CL for different roughness height.

Figure 12 $(a)$ shows the profiles of the leading-order chordwise velocity of the mean-flow distortion (2.59) in the CL. The distribution features a high shear near the critical level $y = y_{c}$ $(\tilde {Y} = 0)$ . Note that although the velocity of the distortion is $\mathit{O} (R^{-1/3})$ , the associated shear, or vorticity, is $\mathit{O} (1)$ . The mean-flow distortion is not confined because it exhibits a jump across the CL whereby driving a mean-flow distortion of the same order-of-magnitude in the main layer. The latter must be considered to determine the constant $ \tilde {u}_{0}^{\langle 1 \rangle } |_{\tilde {Y} \to \infty }$ . Figure 12 $(b)$ displays the variations with roughness height $h$ of the velocity jump (2.61) and shear jump (2.71a ) of the leading- and second-order mean-flow distortions. The magnitudes of both jumps increase with $h$ .

Figure 12. The mean-flow distortion in the CL: $(a)$ the leading-order chordwise velocity profiles for different roughness height and $(b)$ jumps of the mean-flow distortion versus roughness height.

A significant feature of the CL dynamics is that the regularised disturbance velocity, which is of $\mathit{O} (R^{-1/3})$ , significantly exceeds that in the main layer. Similar to the WL, the amplitude of the associated vorticity in the CL can be of $\mathit{O} (1)$ . We also rotate the $x$ – $z$ plane to align the $x$ -axis with the direction of $(\tilde {\alpha },\tilde {\beta })$ , where $\{ \tilde {\alpha },\tilde {\beta } \}$ , introduced in (4.1), are wavenumbers of the secondary instability mode in the distorted CL flow. The chordwise and spanwise velocity components are projected accordingly, as specified in (3.17) with $\bar {\theta }$ being replaced by $\tilde {\theta } = \arctan (\tilde {\beta } / \tilde {\alpha })$ . The vorticity along the direction of $(\tilde {\alpha },\tilde {\beta })$ can be calculated as

(3.20) \begin{equation} \tilde {\varOmega } = - \left ( \tilde {\alpha } \Big / \sqrt {\tilde {\alpha }^{2} + \tilde {\beta }^{2}} \right ) \left [ (\tilde {\beta } / \tilde {\alpha }) \left( U_{\!B}^{\prime} (y_{c}) + \tilde {U}_{\tilde {Y}} \right) - \left( W_{\!B}^{\prime} (y_{c}) + \tilde {W}_{ \tilde {Y}} \right) \right ] + {\cdots} . \end{equation}

Figure 13 displays contours of this skewed CL vorticity without the contributions from base-flow wall shears. For small $h$ , the vorticity is rather weak. However, for $h = 0.43$ (figure 13 b), a significantly distorted flow field is observed. The vorticity in the direction perpendicular to $(\tilde {\alpha },\tilde {\beta })$ controls the local secondary instability in the CL, and its contour plot for the roughness with $h = 0.43$ will be displayed later in § 5.1, where instability results are discussed.

Figure 13. Contours of the skewed CL vorticity $\tilde {\varOmega }$ without the contributions from the base-flow shear in Case I: $(a)$ $h = 0.1$ and $(b)$ $h = 0.43$ .

4. Secondary instability of roughness-distorted three-dimensional boundary layers

We have demonstrated that the amplitudes of the vorticities in both the nonlinearly distorted CL and WL are of $\mathit{O} (1)$ as figure 1 shows. The stability of the distorted flow may be fundamentally altered despite the fact that the velocity distortion is still small, being of $\mathit{O} (R^{-1/3})$ . Next, we show that these layers are susceptible to small-scale secondary instability, the onset of which may form a crucial step in transition to turbulence.

We seek secondary instability in the CL, i.e. CL modes, and in the WL, i.e. WL modes in this section. The CL modes have short-wavelength high-frequency characters and can be formulated as linear generalised eigenvalue problems. The temporal and spatial CL modes are found to be essentially equivalent. The WL modes, on the other hand, have short wavelengths but $\mathit{O} (1)$ frequencies. The temporal WL mode is formulated similarly to the CL modes, while the spatial WL mode is described by a nonlinear eigenvalue problem.

4.1. Secondary instability in the CL

We start with investigating spatial secondary instability in the CL. The wall-normal length scale of the local flow is $\mathit{O} (R^{-1/3})$ . The instability is taken to have characteristic wavelength of this order. Since the local velocity $(U_{\!B} (y_c), W_{\!B} (y_c) )$ is of $\mathit{O} (1)$ , the instability would have $\mathit{O} (R^{1/3})$ frequencies. With such a perception in mind, we perturb the nonlinearly distorted CL flow by a disturbance of $\mathit{O} (\epsilon _{s})$ , and express the locally perturbed flow as

(4.1) \begin{align} (u,v,w,p) &= \bigl (U_B(y_c) + R^{-1/3} U_{\!B}^{\prime}(y_c) \tilde {Y}, R^{-1} V_B(y_c), W_B(y_c) + R^{-1/3} W_{B}^{\prime}(y_c) \tilde {Y}, P_B \bigr ) \nonumber \\ & \quad + \left (R^{-1/3} \tilde {U} (x, \tilde {Y},z), R^{-2/3} \tilde {V} (x, \tilde {Y}, z), R^{-1/3} \tilde {W} (x, \tilde {Y},z), R^{-2/3} \tilde {P} (x, z) \right ) \nonumber \\ & \quad + \epsilon _{s} \left [ \left (\tilde {U}^{\langle s \rangle } (x, \tilde {Y}, z), \tilde {V}^{\langle s \rangle } (x, \tilde {Y}, z), \tilde {W}^{\langle s \rangle } (x, \tilde {Y}, z), R^{-1/3} \tilde {P}^{\langle s \rangle } (x, \tilde {Y}, z) \right )\right . \nonumber \\ & \quad \left . \times \textrm{exp} \left ( \textrm{i} R^{1/3} ( \tilde {\alpha } x + \tilde {\beta } z - \tilde {\omega } t ) \right ) + \textrm{c.c.} \right ]\!, \end{align}

where the last term represents a secondary instability mode and $\epsilon _{s} \ll 1$ . Note that the factor $R^{1/3}$ in the exponent indicates the high-frequency short-wavelength character of the instability with $\tilde {\omega }$ and $(\tilde {\alpha }, \tilde {\beta })$ representing the rescaled frequency and wavenumbers. Furthermore, we require that

(4.2) \begin{equation} \tilde {\omega } = \tilde {\alpha } U_{\!B} (y_c) + \tilde {\beta } W_{\!B} (y_c) \end{equation}

such that the secondary instability is trapped within the CL. For this spatial secondary instability, we fix the high-frequency $\tilde {\omega }$ , or equivalently the wavenumbers $\{ \tilde {\alpha }, \tilde {\beta } \}$ , to be real, and seek $\mathit{O} (1)$ modal growth in the $x$ -direction. Substituting (4.1) into the N–S (2.3), we obtain the governing equations for the spatial CL mode,

(4.3a) \begin{align} \textrm{i} \tilde {\alpha } \tilde {U}^{\langle s \rangle } + \tilde {V}_{\tilde {Y}}^{\langle s \rangle } + \textrm{i} \tilde {\beta } \tilde {W}^{\langle s \rangle } = 0, \end{align}

(4.3b) \begin{align} \tilde {\mathcal{L}} \tilde {U}^{\langle s \rangle } + \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} ) \tilde {U}^{\langle s \rangle } + (U_{\!B}^{\prime}(y_c)+\tilde {U}_{\tilde {Y}}) \tilde {V}^{\langle s \rangle } + \textrm{i} \tilde {\alpha } \tilde {P}^{\langle s \rangle } = 0, \end{align}

(4.3c) \begin{align} \tilde {\mathcal{L}} \tilde {V}^{\langle s \rangle } + \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} ) \tilde {V}^{\langle s \rangle } + \tilde {P}_{\tilde {Y}}^{\langle s \rangle } = 0, \end{align}

(4.3d) \begin{align} \tilde {\mathcal{L}} \tilde {W}^{\langle s \rangle } + \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} ) \tilde {W}^{\langle s \rangle } + (W_{B}^{\prime}(y_c)+\tilde {W}_{\tilde {Y}} ) \tilde {V}^{\langle s \rangle } + \textrm{i} \tilde {\beta } \tilde {P}^{\langle s \rangle } = 0, \end{align}

where we have introduced the operator

(4.4) \begin{equation} \tilde {\mathcal{L}} = U_B(y_c) \frac {\partial }{\partial x} + W_B(y_c) \frac {\partial }{\partial z} + \textrm{i} \bigl [ \tilde {\alpha } U_{\!B}^{\prime}(y_c) + \tilde {\beta } W_{B}^{\prime}(y_c) \bigr ] \tilde {Y}. \end{equation}

By eliminating $\tilde {U}^{\langle s \rangle }$ , $\tilde {W}^{\langle s \rangle }$ and $\tilde {P}^{\langle s \rangle }$ among the equations in (4.3), the system is reduced to a single equation for $\tilde {V}^{\langle s \rangle }$ ,

(4.5) \begin{equation} \left [ \tilde {\mathcal{L}} + \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} ) \right ] \left (\frac {\partial ^2}{\partial \tilde {Y}^2} - \tilde {k}^{2} \right ) \tilde {V}^{\langle s \rangle } - \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} )_{ \tilde {Y} \tilde {Y} } \tilde {V}^{\langle s \rangle } = 0, \end{equation}

where $\tilde {k} = \sqrt {\tilde {\alpha }^{2} + \tilde {\beta }^{2}}$ .

To formulate the eigenvalue problem for the spatial CL mode, we assume a modal dependence of $x$ and express the instability solution in the form,

(4.6) \begin{equation} ( \tilde {U}^{\langle s \rangle }, \tilde {V}^{\langle s \rangle }, \tilde {W}^{\langle s \rangle }, \tilde {P}^{\langle s \rangle } ) = ( \tilde {u}^{\langle s \rangle }, \tilde {v}^{\langle s \rangle }, \tilde {w}^{\langle s \rangle }, \tilde {p}^{\langle s \rangle } ) \textrm{exp} ( \tilde {\sigma } x ), \end{equation}

where $\tilde {\sigma } = \tilde {\sigma }_{r} + \textrm{i} \tilde {\sigma }_{i}$ is to be sought as the eigenvalue, with $\tilde {\sigma }_{r}$ representing the growth rate, which is of $\mathit{O} (1)$ . Substituting (4.6) into (4.5), we obtain the governing equation

(4.7) \begin{equation} \Bigl [ \tilde {\sigma } U_{\!B} (y_{c}) + \tilde {\mathcal{L}} + \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} ) \Bigr ] \left (\frac {\partial ^2}{\partial \tilde {Y}^2} - \tilde {k}^{2} \right ) \tilde {v}^{\langle s \rangle } - \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} )_{ \tilde {Y} \tilde {Y} } \tilde {v}^{\langle s \rangle } = 0. \end{equation}

Since the leading-order roughness-induced disturbance concentrates in the CL, (4.7) collapses to $(\partial _{\tilde {Y} \tilde {Y}}^{2} - \tilde {k}^{2} ) \tilde {v}^{\langle s \rangle } \to 0$ as $\tilde {Y} \to \pm \infty$ , which is further simplified to give the boundary conditions for the CL mode

(4.8) \begin{equation} \left ( \frac {\partial }{\partial \tilde {Y}} \mp \tilde {k} \right ) \tilde {v}^{\langle s \rangle } = 0 \quad \textrm{as} \, \, \, \tilde {Y} \to \mp \infty ; \end{equation}

these serve to exclude the non-physical unbounded (exponentially large) solutions.

The roughness-induced disturbance velocities $\tilde {U}$ and $\tilde {W}$ are expressed as Fourier series, as indicated by (2.52), and we only need to retain the leading-order terms here,

(4.9) \begin{equation} (\tilde {U}, \tilde {W}) = \sum _{n = - \infty }^{\infty } \left ( \tilde {u}_{n}^{\langle 1 \rangle } (\tilde {Y}), \tilde {w}_{n}^{\langle 1 \rangle } (\tilde {Y}) \right ) \textrm{e}^{\textrm{i} n \alpha _{w} \zeta } + {\cdots} . \end{equation}

The coefficients of the (4.7) are periodic functions of $\zeta$ , and thus according to Floquet theory we express $\tilde {v}^{\langle s \rangle } (x, \tilde {Y}, z)$ in the form of Fourier series with respect to $\zeta$ ,

(4.10) \begin{equation} \tilde {v}^{\langle s \rangle } (x, \tilde {Y}, z) = \ \sum _{n = - \infty }^{\infty } \tilde {v}_{n}^{\langle s \rangle } (\tilde {Y}) \textrm{e}^{\textrm{i} n \alpha _{w} \zeta }; \end{equation}

in general the right-hand side may have an exponential factor $\textrm{exp} (\textrm{i} q \alpha _{w} \zeta )$ with $0 \le q \le 1/2$ being a Floquet-exponent-like parameter. It turns out that the resulting extra term $\textrm{i} q \beta _{w} W_{\!B} (y_c)$ can be absorbed into $\tilde {\sigma }$ , and thus we set $q = 0$ without losing generality. Substituting (4.9) and (4.10) into (4.7) with the definition of the critical $y_{c}$ leads to the governing equation for $\tilde {v}_{n}^{\langle s \rangle }$ ,

(4.11) \begin{align} & \textrm{i} ( \tilde {\alpha } U_{\!B}^{\prime}(y_c) + \tilde {\beta } W_{B}^{\prime}(y_c)) \tilde {Y} \left (\!\frac {\textrm{d}^2}{\textrm{d} \tilde {Y}^2} - \tilde {k}^{2} \!\right ) \tilde {v}_{n}^{\langle s \rangle } + \sum _{k = -\infty }^{\infty } \textrm{i} \left [ \!\left(\tilde {\alpha } \tilde {u}_{n-k}^{\langle 1 \rangle } + \tilde {\beta } \tilde {w}_{n - k}^{\langle 1 \rangle } \!\right) \!\left (\!\frac {\textrm{d}^2}{\textrm{d} \tilde {Y}^2} - \tilde {k}^{2} \!\right ) \tilde {v}_{k}^{\langle s \rangle } \right . \nonumber \\ & \quad \left . - ( \tilde {\alpha } \tilde {u}_{n-k}^{\langle 1 \rangle } + \tilde {\beta } \tilde {w}_{n-k}^{\langle 1 \rangle } )_{ \tilde {Y} \tilde {Y} } \tilde {v}_{k}^{\langle s \rangle } \right ] = - \tilde {\sigma } U_{\!B} (y_c) \left (\frac {\textrm{d}^2}{\textrm{d} \tilde {Y}^2} - \tilde {k}^{2} \right ) \tilde {v}_{n}^{\langle s \rangle }. \end{align}

The above equations are subject to boundary conditions (4.8). Discretisation of (4.11) and (4.8) gives rise to a linear generalised eigenvalue problem,

(4.12) \begin{equation} \tilde {\boldsymbol{A}} \tilde {\boldsymbol{V}}^{\langle s \rangle } = \tilde {\sigma } \tilde {\boldsymbol{B}} \tilde {\boldsymbol{V}}^{\langle s \rangle }, \end{equation}

where the solution vector $\tilde {\boldsymbol{V}}^{\langle s \rangle }$ , the matrices $\tilde {\boldsymbol{A}}$ and $\tilde {\boldsymbol{B}}$ are defined in the supplementary material (§ S4), and further details can be found there.

For temporal secondary instability, the perturbation in (4.1) is replaced by the CL mode of the form,

(4.13) \begin{align} & \epsilon _{s} \Bigl [ \bigl (\check {U}^{\langle s \rangle } (x, \tilde {Y}, z), \check {V}^{\langle s \rangle } (x, \tilde {Y}, z), \check {W}^{\langle s \rangle } (x, \tilde {Y}, z), R^{-1/3} \check {P}^{\langle s \rangle } (x, \tilde {Y}, z) \bigr ) \nonumber \\ & \qquad \qquad \qquad \qquad \qquad \qquad \times \textrm{exp} \bigl ( \check {\sigma } t + \textrm{i} R^{1/3} ( \tilde {\alpha } x + \tilde {\beta } z - \tilde {\omega } t ) \bigr ) + \textrm{c.c.} \Bigr ], \end{align}

where the condition (4.2) is required to hold, and $\check {\sigma }$ stands as the eigenvalue, with its real part representing the temporal growth rate. By expressing $\check {V}^{\langle s \rangle } (x, \tilde {Y}, z)$ as a Fourier series with respect to $\zeta$ ,

(4.14) \begin{equation} \check {V}^{\langle s \rangle } (x, \tilde {Y}, z) = \sum _{k = - \infty }^{\infty } \check {V}_{n}^{\langle s \rangle } \textrm{e}^{\textrm{i} n \alpha _{w} \zeta }, \end{equation}

the governing equation for the Fourier component $\check {V}_{n}^{\langle s \rangle }$ can be derived and turns out to be the same as that for spatial instability, i.e. (4.11), with $\check {\sigma }$ playing exactly the role of $\tilde {\sigma } U_{\!B} (y_{c})$ , leading to the equivalence

(4.15) \begin{equation} \check {\sigma } = U_{\!B} (y_c) \tilde {\sigma }. \end{equation}

Therefore, we will not treat the temporal instability separately, and only show numerical results for spatial CL modes. As this small-scale local instability is newly identified, verification through initial-value calculations is necessary; the formulation and numerical results are presented in the supplementary material (§ S3).

4.2. Secondary instability in the WL

We now analyse instability of the roughness-distorted WL flow to small-scale perturbation. The WL has an $\mathit{O} (R^{-1/3})$ width, and we seek WL modes with characteristic wavelength of this order. Unlike the CL, the local velocity is $\mathit{O} (R^{-1/3})$ , and so the characteristic frequency of WL modes is of $\mathit{O} (1)$ .

4.2.1. Temporal modes

We first study the temporal WL mode. The instantaneous WL flow, consisting of the roughness-distorted part and the $\mathit{O} (\epsilon _{s})$ instability mode, is expressed as

(4.16) \begin{align} (u, v, w, p) &= \left(\! R^{-1/3} \bar {U}_{w} (x, \bar {Y}\!, z), R^{-2/3} \bar {V}_{w} (x, \bar {Y}\!, z), R^{-1/3} \bar {W}_{w} (x, \bar {Y}\!, z), R^{-2/3} \bar {P}_{w} (x, z)\!\right) \nonumber \\ & \quad + \epsilon _{s} \bigl [ \bigl ( \bar {U}^{\langle s \rangle } (x, \bar {Y}, z), \bar {V}^{\langle s \rangle } (x, \bar {Y}, z), \bar {W}^{\langle s \rangle } (x, \bar {Y}, z), R^{-1/3} \bar {P}^{\langle s \rangle } (x, \bar {Y}, z) \bigr ) \nonumber \\ & \quad \times \textrm{exp} \bigl ( \bar {\sigma } t + \textrm{i} R^{1/3}(\bar {\alpha } x + \bar {\beta } z ) \bigr ) + \textrm{c.c.} \bigr ], \end{align}

where the last term represents a WL mode with its scaled chordwise and spanwise wavenumbers $\bar {\alpha }$ and $\bar {\beta }$ being taken as the real parameters, and $\bar {\sigma } = \bar {\sigma }_{r} + \textrm{i} \bar {\sigma }_{i}$ is to be sought as the eigenvalue, whose real and imaginary parts represent the temporal growth rate and frequency of the WL mode, respectively. By inserting (4.16) into (2.3), the governing equations for the temporal WL mode are derived as

(4.17a) \begin{align} \textrm{i} \bar {\alpha } \bar {U}^{\langle s \rangle } + \bar {V}_{\bar {Y}}^{\langle s \rangle } + \textrm{i} \bar {\beta } \bar {W}^{\langle s \rangle } = 0, \end{align}

(4.17b) \begin{align} \bar {\sigma } \bar {U}^{\langle s \rangle } + \textrm{i} \left ( \bar {\alpha } \bar {U}_{w} + \bar {\beta } \bar {W}_{w} \right ) \bar {U}^{\langle s \rangle } + \bar {V}^{\langle s \rangle } \bar {U}_{w, \bar {Y}} + \textrm{i} \bar {\alpha } \bar {P}^{\langle s \rangle } = 0, \end{align}

(4.17c) \begin{align} \bar {\sigma } \bar {V}^{\langle s \rangle } + \textrm{i} \left ( \bar {\alpha } \bar {U}_{w} + \bar {\beta } \bar {W}_{w} \right ) \bar {V}^{\langle s \rangle } + \bar {P}_{\bar {Y}}^{\langle s \rangle } = 0, \end{align}

(4.17d) \begin{align} \bar {\sigma } \bar {W}^{\langle s \rangle } + \textrm{i} \left ( \bar {\alpha } \bar {U}_{w} + \bar {\beta } \bar {W}_{w} \right ) \bar {W}^{\langle s \rangle } + \bar {V}^{\langle s \rangle } \bar {W}_{w, \bar {Y}} + \textrm{i} \bar {\beta } \bar {P}^{\langle s \rangle } = 0. \end{align}

Equations in (4.17d ) can be reduced into a single equation for $\bar {V}^{\langle s \rangle }$ ,

(4.18) \begin{equation} \Bigl [ \bar {\sigma } + \textrm{i} ( \bar {\alpha } \bar {U}_{w} + \bar {\beta } \bar {W}_{w} ) \Bigr ] \left ( \frac {\partial ^2}{\partial \bar {Y}^2} - \bar {k}^{2} \right ) \bar {V}^{\langle s \rangle } - \textrm{i} \left ( \bar {\alpha } \bar {U}_{w} + \bar {\beta } \bar {W}_{w} \right )_{\bar {Y} \bar {Y}} \bar {V}^{\langle s \rangle } = 0, \end{equation}

where $\bar {k} = \sqrt { \bar {\alpha }^{2} + \bar {\beta }^{2} }$ . Since the roughness-distorted WL flow, given by (2.25) and calculated in § 3, is expressed in the form of Fourier series of $\zeta$ , the WL mode can be expressed in the same form,

(4.19) \begin{equation} \bar {V}^{\langle s \rangle } \left (x, \bar {Y}, z \right ) = \sum _{n = -\infty }^{\infty } \bar {V}_{n}^{\langle s \rangle } \textrm{e}^{\textrm{i} n \alpha _w \zeta }. \end{equation}

Substitution of (2.25) and (4.19) into (4.18) leads to

(4.20) \begin{align} & \textrm{i} ( \bar {\alpha } \lambda _1 + \bar {\beta } \lambda _3 ) \bar {Y} \left ( \frac {\textrm{d}^2}{\textrm{d} \bar {Y}^2} - \bar {k}^{2} \right ) \bar {V}_{n}^{\langle s \rangle } \nonumber \\ & \quad + h \sum _{k = -\infty }^{\infty } \textrm{i} \biggl \{ \Bigl [ \bar {\alpha } \left ( \lambda _1 F_{n-k} + \bar {U}_{n-k} \right ) + \bar {\beta } \left ( \lambda _3 F_{n-k} + \bar {W}_{n-k} \right ) \Bigr ] \left ( \frac {\textrm{d}^2}{\textrm{d} \bar {Y}^2} - \bar {k}^{2} \right ) \bar {V}_{k}^{\langle s \rangle } \nonumber \\ & \quad - ( \bar {\alpha } \bar {U}_{n-k} + \bar {\beta } \bar {W}_{n-k} )_{\bar {Y} \bar {Y}} \bar {V}_{k}^{\langle s \rangle } \biggr \} = -\bar {\sigma } \left ( \frac {\textrm{d}^2}{\textrm{d} \bar {Y}^2} - \bar {k}^{2} \right ) \bar {V}_{n}^{\langle s \rangle }. \end{align}

The secondary instability mode is trapped within the WL, that is, it decays exponentially as $\bar {Y} \to \infty$ , as can be inferred from (4.18). However, unlike the CL mode, which attenuates exponentially beneath the CL, the WL mode has to satisfy the impermeability condition at the wall surface. Therefore, the boundary conditions for (4.20) are

(4.21a,b) \begin{align} \left ( \frac {\textrm{d}}{\textrm{d} \bar {Y}} + \bar {k} \right ) \bar {V}_{n}^{\langle s \rangle } = 0 \quad \textrm{as} \, \, \, \bar {Y} \to \infty ; \qquad \bar {V}_{n}^{\langle s \rangle } = 0 \quad \textrm{at} \, \, \, \bar {Y} = 0. \end{align}

Again, discretisation of (4.20) and (4.21) leads to a linear generalised eigenvalue problem similar to (4.12),

(4.22) \begin{equation} \bar {\boldsymbol{A}} \bar {\boldsymbol{V}}^{\langle s \rangle } = \bar {\sigma } \bar {\boldsymbol{B}}\bar {\boldsymbol{V}}^{\langle s \rangle }, \end{equation}

where the matrices and the details of the formulation are provided in the supplementary material (§ S4).

4.2.2. Spatial modes

Finally, we consider the spatial WL mode. The temporal mode in (4.16) is replaced by

(4.23) \begin{align} & \epsilon _{s} \bigl [ (\bar {\bar {U}}^{\langle s \rangle } (x, \bar {Y}, z), \bar {\bar {V}}^{\langle s \rangle } (x, \bar {Y}, z), \bar {\bar {W}}^{\langle s \rangle } (x, \bar {Y}, z), R^{- 1 / 3} \bar{\bar{P}}^{\langle s \rangle} (x, \bar {Y}, z)) \nonumber \\ & \quad \times \textrm{exp} \bigl (\textrm{i} R^{1/3} (\bar {\bar {\alpha }} x + \bar {\bar {\beta }} z ) - \textrm{i} \bar {\bar {\omega }} t \bigr ) + \textrm{c.c.} \bigr ], \end{align}

where $\bar {\bar {\beta }}$ and $\bar {\bar {\omega }}$ are taken as real parameters, and $\bar {\bar {\alpha }} = \bar {\bar {\alpha }}_{r} + \textrm{i} \bar {\bar {\alpha }}_{i}$ is the eigenvalue with $\bar {\bar {\alpha }}_{r}$ and $- \bar {\bar {\alpha }}_{i}$ representing the wavenumber and the spatial growth rate, respectively. Note that the spatial growth rate is of $\mathit{O} (R^{1/3})$ . The governing equation for $\bar {\bar {V}}^{\langle s \rangle }$ can be derived as

(4.24) \begin{equation} \Bigl [ - \bar {\bar {\omega }} + \left( \bar {\bar {\alpha }} \bar {U}_{w} + \bar {\bar {\beta }} \bar {W}_{w} \right) \Bigr ] \left ( \frac {\partial ^2}{\partial \bar {Y}^2} - \bar {\bar {k}}^{2} \right ) \bar {\bar {V}}^{\langle s \rangle } - ( \bar {\bar {\alpha }} \bar {U}_{w} + \bar {\bar {\beta }} \bar {W}_{w} )_{\bar {Y} \bar {Y}} \bar {\bar {V}}^{\langle s \rangle } = 0, \end{equation}

where $\bar {\bar {k}} = \sqrt {\bar {\bar {\alpha }}^{2} + \bar {\bar {\beta }}^{2}}$ . We express $\bar {\bar {V}}^{\langle s \rangle }$ as a Fourier series of $\zeta$ and substitute it, along with (2.25), into (4.24). The resulting system and the boundary conditions form a nonlinear eigenvalue problem, which is, after truncation and discretisation, converted into an algebraic system

(4.25) \begin{equation} \bar {\bar {\boldsymbol{M}}} \bar {\bar {\boldsymbol{V}}}^{\langle s \rangle } = 0. \end{equation}

The requirement of a non-trivial solution leads to an implicit dispersion relation for the spatial secondary instability

(4.26) \begin{equation} \mathcal{D} (\bar {\bar {\alpha }}; \bar {\bar {\beta }}, \bar {\bar {\omega }}) \equiv \textrm{det} \left ( \bar {\bar {\boldsymbol{M}}} (\bar {\bar {\alpha }}; \bar {\bar {\beta }}, \bar {\bar {\omega }}) \right ) = 0. \end{equation}

The derivation and the matrix $\bar {\bar {\boldsymbol{M}}}$ are given in the supplementary material (§ S4).

5. Numerical results for the secondary instability

5.1. Secondary instability in the CL

For the secondary instability in the CL, we only calculate the spatial CL mode by solving the linear eigenvalue formulation (4.12) directly using an appropriate MATLAB eigenvalue solver. It is useful to note symmetry properties of the eigenvalues. First, (4.11) indicates that for each wavenumber pair $(\tilde {\alpha }, \tilde {\beta })$ , $\tilde {\sigma }$ obeys the relation

(5.1) \begin{equation} \tilde {\sigma } (-\tilde {\alpha }, -\tilde {\beta } ) = - \tilde {\sigma }(\tilde {\alpha }, \tilde {\beta }). \end{equation}

Second, if (4.7) with boundary conditions (4.8) admits an eigenvalue $\tilde {\sigma }$ with eigenfunction $\tilde {\boldsymbol{V}}^{\langle s \rangle }$ , then $- (\tilde {\sigma })_{cc}$ and $(\tilde {\boldsymbol{V}}^{\langle s \rangle })_{cc}$ also represent an eigenvalue and eigenfunction, which follows from taking complex conjugate of (4.7). For a given pair $(\tilde {\alpha }, \tilde {\beta })$ , we can calculate the most unstable eigenmode whose eigenvalue $\tilde {\sigma }$ has the largest positive real part, while $-(\tilde {\sigma })_{cc}$ is another eigenvalue corresponding to the most damped mode. On the other hand, according to (5.1), we have $-(\tilde {\sigma })_{cc} = [ \tilde {\sigma } (-\tilde {\alpha },-\tilde {\beta }) ]_{cc}$ . The latter is the most unstable eigenmode for $(-\tilde {\alpha }, -\tilde {\beta })$ . Therefore, it suffices to consider only the upper-half of the $\tilde {\alpha }$ – $\tilde {\beta }$ parameter plane, i.e. restrict the calculations to $\tilde {\beta } \gt 0$ .

Since the matrix formulation in (4.12) involves two-dimensional discretisation in the transverse direction $\tilde {Y}$ and Fourier spectral space, the resulting matrices are large, making direct eigenvalue calculation challenging. We begin by calculating the eigenvalue of (4.12) using a coarse grid with the direct method (MATLAB function eig). The results from the coarse grid are then used as initial guesses to initialise iterations in the Arnoldi method (MATLAB function eigs) to obtain more accurate eigenvalues on a refined grid.

We calculate the spatial CL mode supported by the nonlinearly distorted CLs pertinent to Cases I and II roughness, with heights $h = 0.43$ and $h = 1.135$ , respectively (which are the greatest heights for which convergent solutions for roughness-induced flows can be computed). The numerical results corresponding to Case II roughness are presented in the supplementary material (§ S2). Figure 14 presents contours of the growth rates $ \tilde {\sigma }_{r}$ in the upper-half of the $\tilde {\alpha }$ – $\tilde {\beta }$ plane. There exist two families of modes with the maximal growth rate $\tilde {\sigma }_{r} = 0.0280$ located at $(\tilde {\alpha }, \tilde {\beta }) = (-0.21,0.05)$ . It is noteworthy that the leading-order spanwise velocity of the roughness-distorted flow $\tilde {w}_{n}^{\langle 1 \rangle }$ is proportional to its chordwise counterpart $\tilde {u}_{n}^{\langle 1 \rangle }$ by the factor $- \alpha _{w} / \beta _{w}$ , as indicated by (2.54a ). This implies that when $(\tilde {\alpha }, \tilde {\beta })$ aligns with $(\alpha _{w}, \beta _{w})$ , the leading-order roughness-induced flow controlling the secondary instability is cancelled out in (4.11), which becomes degenerated. The present dominant instability no longer exists, and a weaker instability might be possible if the higher-order terms, $\tilde {u}_{n}^{\langle 2 \rangle }$ and $\tilde {w}_{n}^{\langle 2 \rangle }$ , are taken into account.

Figure 14. Contours of the growth rate $ \tilde {\sigma }_{r}$ of the spatial CL mode on the nonlinear CL flow with $h = 0.43$ . The maximal growth rate is indicated by the red cross. The red dashed lines represent the direction parallel to $(\alpha _{w}, \beta _{w})$ , along which the leading-order effect of the roughness-induced disturbance cancels out.

We now present quantitative information about the instability characteristics. In figure 15, the eigenvalues are plotted for three typical values of $\tilde {\beta }$ . For each $\tilde {\beta }$ , there exist two branches of modes, one branch with negative $\tilde {\alpha }$ featuring the dominant peak and a subdominant peak confined to small but negative $\tilde {\alpha }$ , and another branch with $\tilde {\alpha }$ mostly being positive. The two branches appear to be disconnected with the left-hand having $\tilde {\sigma }_{i} \lt 0$ , and the right-hand having $\tilde {\sigma }_{i} \gt 0$ . Furthermore, for the left-hand branch $\tilde {\sigma }_{i}$ varies almost linearly with respect to $\tilde {\alpha }$ , indicating that the mode is nearly non-dispersive. In contrast, the right-hand branch mode is dispersive. For $\tilde {\beta } = \{ 0.03, 0.05,0.07 \}$ chosen in figure 15, the roughness wavenumber-aligned $\tilde {\alpha }$ are $\{ -0.026, -0.043, -0.061 \}$ , respectively, the position of which are marked in figure 15 $(a)$ . The instability disappears when approaching these points. Initial-value calculations were performed to confirm some of eigenvalue results shown in figure 15 (see the supplementary material (§ S3)). When the spatial growth rates shown in figure 15 are converted using (4.15), the corresponding temporal growth rates are in the range of $(0, 0.025)$ , which is approximately twice that of the secondary instability of stationary cross-flow vortices at the same Reynolds number (Malik et al. 1994).

Figure 15. Instability characteristics of the spatial CL mode for different values of $\tilde {\beta }$ . The black-filled symbols on the $\tilde {\alpha }$ -axis mark the values of $\tilde {\alpha }$ aligning with the roughness wavenumber for each $\tilde {\beta }$ .

Figure 16 displays contours of the normalised eigenfunction $\tilde {\boldsymbol{V}}^{\langle s \rangle }$ of the most unstable mode for $(\tilde {\alpha }, \tilde {\beta }) = (-0.21,0.05)$ and $(\tilde {\alpha }, \tilde {\beta }) = (0.03,0.05)$ , respectively. Superposed are the contours of the $\mathit{O} (1)$ CL vorticity, defined as

(5.2) \begin{equation} \tilde {\varOmega }_{\perp } = (\tilde {\alpha } / \tilde {k}) \Bigl \{ U_{\!B}^{\prime} (y_c) + (\tilde {\beta } / \tilde {\alpha }) W_{\!B}^{\prime} (y_c) + \sum _{n = -\infty }^{\infty } \bigl [ \tilde {u}_{n}^{\langle 1 \rangle } + (\tilde {\beta } / \tilde {\alpha }) \tilde {w}_{n}^{\langle 1 \rangle } \bigr ]_{\tilde {Y}} \textrm{e}^{\textrm{i} n \alpha _{w} \zeta } \Bigr \}, \end{equation}

which represents the vorticity component projected to the direction perpendicular to $(\tilde {\alpha },\tilde {\beta })$ . The term in summation in (5.2) corresponds to the roughness-induced vorticity. It is the component that fundamentally alters the stability (despite small $\mathit{O} (R^{-1/3})$ velocity distortions), and indeed controls the secondary instability as (4.11) indicates. The CL mode with negative $\tilde {\alpha }$ appears rather local, residing primarily in a relatively small region where $\tilde {\varOmega }_{\perp }$ is negative (figure 16 a). With its short-scale carrier wave being accounted for, rapid oscillations in time and space are expected to be observed in this region. In contrast, the CL mode with positive $\tilde {\alpha }$ spreads out in the planar directions, thus exhibiting a distributed character (figure 16 b).

Figure 16. Contours of the normalised eigenfunction corresponding to the most unstable spatial CL mode (solid lines), superposed onto the eigenfunction contours are the skewed vorticity $\tilde {\varOmega }_{\perp }$ of the roughness-distorted CL flow (dashed lines and labelled with values): (a) $(\tilde {\alpha },\tilde {\beta }) = (-0.21,0.05)$ ; (b) $(\tilde {\alpha },\tilde {\beta }) = (0.03,0.05)$ .

5.2. Secondary instability in the WL

5.2.1. Temporal WL mode

For the temporal WL mode, the linear generalised eigenvalue problem (4.22) is first solved using a MATLAB eigenvalue solver. As the temporal mode obeys the same symmetries as the CL mode, we restrict our calculations to $\bar {\beta } \gt 0$ in the $\bar {\alpha }$ – $\bar {\beta }$ parameter plane.

We sought the temporal WL mode using the numerical solutions for the nonlinear roughness-distorted WL flow in Case I with $h = 0.825$ and Case II with $h = 1.135$ , which are the largest heights for which converged solutions can be computed. The numerical results for Case II are shown in the supplementary material (§ S2). Figure 17 shows contours of the growth rate $\bar {\sigma }_{r}$ in the upper-half of the $\bar {\alpha }$ – $\bar {\beta }$ plane, revealing a narrow band of modes exhibiting significant growth rates. The maximal growth rate, $\bar {\sigma }_{r} = 0.0071$ , is identified, attained at $(\bar {\alpha },\bar {\beta }) = (-0.092,0.185)$ . Figure 18 displays the instability characteristics for three different values of $\bar {\beta }$ . For each $\bar {\beta }$ , there exists a band of unstable modes, an interesting feature of which is the presence of a dominant peak and a subdominant peak; this is also shown in figure 17 $(b)$ , where there is a gap between two slanted regions. The linear variation of $\bar {\sigma }_{i}$ with $\bar {\alpha }$ indicates that the WL modes have almost the same phase speeds and are thus weakly dispersive, similar to the CL modes with $\tilde {\alpha } \lt 0$ .

Figure 17. Contours of the growth rate $\bar {\sigma }_{r}$ of the temporal WL mode (a) and local contours (b) marked by the rectangle in (a). The maximum is indicated by the red cross. A red diamond in the right panel indicates a subpeak for $\bar {\beta } = 0.18$ shown in figure 18 $(a)$ .

Figure 18. Instability characteristics of the temporal WL mode for different values of $\bar {\beta }$ .

We then plot in figure 19 contours of the normalised eigenfunction $\bar {\boldsymbol{V}}^{\langle s \rangle }$ of an unstable temporal WL mode, $(\bar {\alpha },\bar {\beta }) = (-0.10,0.20)$ . Superposed in the background are the contours of the $\mathit{O} (1)$ vorticity of the roughness-distorted WL flow, defined as

(5.3) \begin{equation} \bar {\varOmega }_{\perp } = (\bar {\alpha } / \bar {k}) \Bigl \{ \lambda _{1} + (\bar {\beta } / \bar {\alpha }) \lambda _{3} + h \sum _{n = - \infty }^{\infty } \bigl [ \bar {U}_{n} + (\bar {\beta } / \bar {\alpha }) \bar {W}_{n} \bigr ]_{\bar {Y}} \textrm{e}^{\textrm{i} n \alpha _{w} \zeta } \Bigr \}, \end{equation}

which represents the vorticity component projected to the direction perpendicular to that of $(\bar {\alpha },\bar {\beta })$ . The summation term in (5.3) corresponds to the roughness-induced vorticity in the WL, and it is this term that alters the WL stability, and indeed governs the temporal WL mode as can be seen in (4.20). Similar to the CL mode with negative $\tilde {\alpha }$ , the WL mode appears also fairly localised in its horizontal extent that is a function of the wavelength of the roughness arrays.

Figure 19. Contours of the normalised eigenfunction corresponding to an unstable temporal WL mode (solid lines), superimposed are contours of the skewed vorticity $\bar {\varOmega }_{\perp }$ of the roughness-distorted flow (dashed lines and labelled with values).

5.2.2. Spatial WL mode

The spatial WL mode amounts to a nonlinear eigenvalue problem (4.26). The eigenvalue is a root of (4.26), which needs to be obtained by an iteration process with appropriate initial guesses. The discretisation in two dimensions leads to the size of matrix $\bar {\bar {\boldsymbol{M}}}$ being too large to evaluate its determinant directly. We circumvent this difficulty by employing instead an indirect method via a rank-revealing QR factorisation of $\bar {\bar {\boldsymbol{M}}}$ (Kublanovskaya 1970),

(5.4) \begin{equation} \bar {\bar {\boldsymbol{M}}} \boldsymbol{P} = \boldsymbol{Q} \boldsymbol{R}, \end{equation}

with $\boldsymbol{P}$ being a permutation matrix which ensures that the diagonal entries in $\boldsymbol{R}$ form a decreasing sequence according to their magnitude so that the last diagonal entry of $\boldsymbol{R}$ , denoted as $r_{br}$ , becomes zero before any other diagonal entry does. It follows that $\textrm{det} (\bar {\bar {\boldsymbol{M}}}) = \textrm{det} (\boldsymbol{R}) = \prod \textrm{diag}(\boldsymbol{R})$ , which implies that the roots of $r_{br} = 0$ are also the roots of $\textrm{det} (\bar {\bar {\boldsymbol{M}}}) =0$ (Güttel & Tisseur 2017). Therefore, we apply Muller’s method (Muller 1956) to

(5.5) \begin{equation} r_{br} (\bar {\bar {\alpha }}) = 0. \end{equation}

Figure 20. Comparisons between the transformed and directly calculated spatial WL mode for $\bar {\bar {\beta }} = 0.20$ .

On the other hand, a spatial instability mode can be related to a temporal mode through certain transformations, such as the well-known Gaster transformation. As the validity and accuracy of the Gaster transformation are somewhat uncertain, we apply the consistent second-order temporal–spatial transformation recently derived (Xu et al. 2023). We introduce $\bar {\omega } = \textrm{i} \bar {\sigma }$ to rewrite the exponent in (4.16) as

(5.6) \begin{equation} \textrm{exp} \bigl ( \bar {\sigma } t + \textrm{i} R^{1/3} (\bar {\alpha } x + \bar {\beta } z ) \bigr ) = \textrm{exp} \bigl ( \textrm{i} R^{1/3} (\bar {\alpha } x + \bar {\beta } z) - \textrm{i} \bar {\omega } t \bigr ). \end{equation}

The consistent temporal–spatial transformation of Xu et al. (2023), which links the temporal and spatial growth rates for modes of the same real frequency, i.e. $\bar {\bar {\omega }} - \bar {\omega } = - \textrm{i} \bar {\omega }_{i}$ , gives the chordwise wavenumber as

(5.7) \begin{equation} \bar {\bar {\alpha }}_{0} = \bar {\alpha } + \biggl [ -1 + \biggl ( 1 - 2 \textrm{i} \bar {\omega }_{i} \frac {\textrm{d}^{2} \bar {\omega } }{\textrm{d} \bar {\alpha }^{2}} \Big / \left ( \frac {\textrm{d} \bar {\omega }}{\textrm{d} \bar {\alpha }} \right )^{2} \biggr ) ^{1/2} \biggr ] \frac {\textrm{d} \bar {\omega }}{\textrm{d} \bar {\alpha }} \Big / \left ( \frac {\textrm{d}^{2} \bar {\omega }}{\textrm{d} \bar {\alpha }^{2}} \right )\!. \end{equation}

The transformed wavenumber $\bar {\bar {\alpha }}_{0}$ may be taken as an approximation of the spatial eigenvalue, or serve as the initial guess for the iteration process, constructed using Muller’s method, to refine the eigenvalue $\bar {\bar {\alpha }}$ .

We apply the consistent transformation to temporal instability results for $\bar {\bar {\beta }} = 0.20$ . Since the temporal–spatial transformation requires information on the derivatives $\bar {\omega }$ , we need to ensure the temporal instability modes stay on a continuous branch without entangling with other branches. The transformed chordwise wavenumbers $\bar {\bar {\alpha }}_{0}$ are shown in figure 20, where comparison is made with the refined results using the iteration. The good agreement suggests that use of the transformation to temporal instability provides spatial instability characteristics of satisfactory accuracy. Spatial mode is also non-dispersive. For a typical Reynolds number $R = 10^{3}$ , the (unrescaled) spatial growth rate of WL mode, $R^{1/3} (- \bar {\alpha }_{i})$ , is approximately $0.028$ (figure 20), comparable to that of the CL mode (figure 15). We note that these are significantly greater than (approximately twice) the typical spatial growth rates of the secondary instability of stationary cross-flow vortices reported in Högberg & Henningson (1998).

Our calculations suggest that surface roughness may induce small-scale secondary instability within local regions. The present study is, to the best of our knowledge, the first to have demonstrated the existence of such small-scale secondary instabilities, although some of previous numerical calculations found modes appear to be local in their character. Specifically, for the secondary instability of stationary cross-flow vortices, Malik et al. (1999) found that the $y$ -mode and $z$ -mode reside in the bulk of the boundary layer (albeit with the former being at somewhat higher positions), and both have rather high frequencies (4500 and 3500 Hz, respectively). Wassermann & Kloker (2002) identified a high-frequency family and a low-frequency family of modes, and the former is located at appreciably higher positions than the latter, but both were deemed to be the $z$ -mode. A further computational study by Bonfigli & Kloker (2007) revealed two high-frequency modes, which occupy the main boundary layer and correspond to $z$ - and $y$ -modes, respectively, and both of these modes appear to be concentrated in a thin region in the wall-normal direction. A low-frequency $z$ -mode developing close to the wall was also found and it exhibits the feature of a wall mode. These modes were detected in experiments (White & Saric 2005; Serpieri & Kotsonis 2016). Despite visual local characters, the precise local nature of the secondary instability remains uncertain. This is because in all secondary instability analyses, the primary cross-flow vortices in the entire boundary layer were taken as the base flow, and the length scale was taken to be comparable to the local boundary-layer thickness, with the boundary conditions being imposed at the wall and at the outer edge of the background boundary layer, respectively, unlike the present formulation, which specifies explicitly the short wavelength and is pertinent to the local flow field. The boundary conditions (4.8) and (4.21) of the CL and WL modes indicate that these modes are trapped within the CL and WL, respectively, unaffected by the overall boundary-layer flow.

The instability identified in this paper would be manifested as amplification of small-scale high-frequency disturbances, possibly leading to premature transition and ‘flow tripping’. In the presence of spatially extended and periodic roughness arrays as assumed, both temporal and spatial instabilities are likely to be relevant.

From the temporal instability perspective, when the roughness height exceeds a critical value, inevitable spatially extended incidental initial perturbations within the boundary layer, which in general consist of modal components, would undergo collective amplification in time, causing an otherwise laminar state turbulence. This of course does not happen when the roughness height is below the critical value for instability. The instability may thus be associated with flow tripping by roughness. We note that the phenomenon of tripping by streamwise isolated roughness has been explained in terms of instability and resulting transition of the distorted flow (Kurz & Kloker 2016). Tripping by spatially distributed roughness array is another possibility, and could be explained in a similar vein.

From the perspective of spatial instability, one may envisage that the predicted small-scale high-frequency modes are excited at a certain upstream location, and their subsequent amplification could lead to onset of turbulence at a downstream position. The process of amplification and transition are likely to be similar to that associated with the secondary instability of stationary cross-flow vortices.

Given that roughness-induced distortions and the directly induced small-scale instability resemble primary cross-flow vortices and their secondary instability, respectively, we think that experimental investigations of the former would, in principle, be similar to those of the latter. In addition to constructing appropriate patterned roughness arrays, a possible challenge is to avoid or minimise primary cross-flow vortices. This might be achieved by deploying arrays of sufficient height in a region relatively close to the leading edge so that the roughness-induced distortions and the directly induced small-scale instability dominate.

Unfortunately, to the best of our knowledge, there has not been any experimental or direct numerical investigation aimed specifically at the local (CL) structures of the distorted boundary-layer flows by periodically distributed roughness arrays, or the resulting small-scale instability. Nevertheless, we note that there have been works on using regularly distributed roughness to suppress primary cross-flow vortices and their secondary instability (Wassermann & Kloker 2002; Saric et al. 2019; Ide, Hirota & Tokugawa 2021; Zoppini et al. 2022; Suzuki et al. 2024). Those experimental and computational techniques could be adapted to investigate/check the scenarios/findings presented here. On the other hand, we anticipate that our mathematical framework as well as theoretical results would prompt and guide future experimental and numerical investigations into this mechanism.

Of close relevance to our work are the recent studies by Hirota and his colleagues. Ide et al. (2021) investigated the transition-delay effect of sinusoidal roughness elements (SREs) using DNS. Unlike discrete roughness elements, these roughness arrays were distributed sinusoidally along a direction inclined at an angle to the leading edge, in which sense rather similar to the form considered in the present paper, but their overall height is modulated by a chordwise Gaussian envelope with its maximum at a distance downstream. Intriguingly, the inclination is chosen to nearly overlap the streamline at the outer edge of the boundary layer. The DNS results suggest that compared with discrete roughness elements, SREs are more viable and effective because with greater heights they do not cause tripping (i.e. direct transition to turbulence) while exciting in their wake the desired cross-flow vortices of appropriate amplitude (the so-called control mode) to inhibit the targeted most dangerous cross-flow vortices and delay the onset of the secondary instability. The effectiveness and advantage of SREs were confirmed by experiments (Suzuki et al. 2024) while detailed mechanisms remain elusive. Hirota, Ide & Hattori (2024) further investigated the role of the inclination direction of roughness arrays, showing that the desired direction must be close to that of the inviscid free stream streamline, beyond which secondary instability takes place in the immediate wake of the roughness arrays causing rapid transition. From the mathematical viewpoint, this direction is special for (at least) two reasons. Firstly, it pertains to a distinguished cross-flow instability regime, in which stationary vortices are ‘free-streamline aligned’ and of viscous wall-mode type (Choudhari 1995). Secondly, the boundary-layer response to roughness arrays with this orientation does not have a CL in the bulk. As a result, the local CL instability as described in our paper does not arise, and one might thus speculate that this may be the reason why such SREs did not cause tripping. The present theory does not apply to this case, and a new formulation is required to characterise the roughness-induced distortion and possible secondary instability. Furthermore, mathematical theories could be developed for chordwise modulated arrays employed in experiments so that the excitation of the control mode, and the interactions among the targeted incoming most dangerous mode, as well as the roughness-induced distortion could be investigated.

We believe that such local instabilities are rather common and play a pivotal role in boundary-layer transition. In particular, we believe that secondary instabilities leading to breakdown of primary stationary and travelling cross-flow vortices (eigenmodes) are likely to be of local nature, and the present asymptotic approach may be adapted to provide a firm answer as well. Moreover, in fully developed shear turbulence, such small-scale instabilities that are local in space may operate to facilitate direct energy transfer from large-scale structures to much smaller-scale fluctuations, which is non-local in spectral space. Such a process may supplement the continuous energy cascade across scales.

6. Summary and conclusions

We have carried out a theoretical study of effects of distributed surface roughness in the form of a wavy wall on the three-dimensional boundary-layer instability. Despite its height being much smaller than the local boundary-layer thickness $\delta ^{\ast }$ , the roughness is capable of inducing strongly nonlinear distortions within localised regions in the wall-normal direction, the CL and WL, where the distorted flows become susceptible to small-scale secondary instability.

The present analysis is formulated on the basis of a high-Reynolds-number asymptotic framework. We considered a wavy-wall surface roughness with the chordwise and spanwise length scales being $\mathit{O} ( \delta ^{\ast } )$ , and the height being much smaller, of $\mathit{O} (R^{-1/3} \delta ^{\ast })$ , which may cover micron-sized roughness in the dimensional setting. The boundary-layer flow is characterised by a WL, the main layer and a CL. The viscous WL, having a width of $\mathit{O} (R^{-1/3} \delta ^{\ast })$ , accommodates the immediate response to the surface roughness. The surface displacement is converted into an $\mathit{O} (R^{-2/3})$ blowing velocity by the WL. Moreover, the disturbance in the WL is strongly nonlinear and generates a streaming, which is characterised by an $\mathit{O} (R^{-1/3})$ mean-flow distortion. In the bulk of the boundary layer, the disturbance is composed of the $\mathit{O} (R^{-1/3})$ streaming and the $\mathit{O} (R^{-2/3})$ forced disturbance. Of great significance is the forced disturbance, which is governed by the steady Rayleigh equation with the blowing velocity from the WL acting as the lower boundary condition. The Rayleigh equation becomes singular at a special position, denoted as the critical level $y_c$ . This singularity is resolved by introducing a CL around $y_{c}$ , where the viscous effect is reintroduced and its width is of $\mathit{O} (R^{-1/3} \delta ^{\ast })$ . The regularised chordwise and spanwise disturbance velocities in the CL are of $\mathit{O} (R^{-1/3})$ , much larger than those in the main layer. The disturbance in the CL is strongly nonlinear. Most crucially, the vorticities of the nonlinearly distorted flows in the CL and WL, are of $\mathit{O} (1)$ , and it is these quantities that render the wall and CLs susceptible to small-scale secondary instability.

We demonstrated for the first time that a roughness-distorted three-dimensional boundary layer supports small-scale local secondary instability, which is inviscid and biglobal in its nature (Theofilis 2011). The CL modes have high frequencies and short wavelengths comparable to the CL width, and spatial and temporal instabilities were found to be governed by essentially the same linear generalised eigenvalue problems. Calculations were performed and results presented for the spatial CL mode. In order to validate the eigenvalue results of the CL mode, direct initial-value calculations were performed using both the eigenfunction and an arbitrary velocity profile as initial conditions. The temporal WL mode, which has short $(\mathit{O} (R^{-1/3} \delta ^{\ast }))$ wavelengths but $\mathit{O} (1)$ frequencies, was analysed following a similar procedure; a linear generalised eigenvalue problem was solved numerically, and validated through initial-value calculations. The computation of the spatial WL mode is more challenging since it is described by a nonlinear eigenvalue problem, which was solved using Muller’s iterative method. The initial guesses for iteration processes were obtained by applying the consistent temporal–spatial transformation (Xu et al. 2023) to the temporal WL mode. It is worth stressing that our work appears to be the first to establish the existence of such spatially local short-wavelength (high-frequency) instability, which probably plays a fundamental role in facilitating spectrally non-local energy transfer from large to much smaller scales in turbulent shear flows.

The roughness configuration adopted in this research is quite simple, which has a single aligned orientation with respect to the chordwise direction. Hence further work could examine the effects of surface roughness in more general forms such as those consisting of multiple components in different orientations corresponding to a multitude of pairs of the roughness wavenumbers $(\alpha _{w},\beta _{w})$ . As a result, a continuum of CLs would emerge instead of a single one. The resulting secondary instability would be different but of interest. The present asymptotic framework with Fourier series expansions is limited to spatially regular (periodic or quasiperiodic) roughness configurations, but it can be extended to randomly distributed roughness, in which case an appropriate statistical characterisation of the roughness is necessary. Moreover, the secondary instability analysis in this paper may be adapted to investigate the detailed structure of the short-wavelength (high-frequency) local secondary instability of nonlinearly evolving stationary and travelling cross-flow vortices.