2. The instability equations

Consider a shear flow of an incompressible viscous fluid above a wavy wall of wavelength $b$ having crests inclined at an angle ${{\rm \pi} }/{2}-\theta$ to the main flow direction; see figure 1. The amplitude of the undulations is taken to be $2\epsilon b$ where $\epsilon$ is small, a more precise size for $\epsilon$ will be given later. We suppose that when $\epsilon =0$ the flow is predominantly in the $x^*$ direction with respect to a Cartesian coordinate system $(x^*,y^*,z^*)$. We use the wall wavelength as a typical length scale and so define dimensionless variables $(x,y,z)=(x^*,y^*,z^*)b^{-1}$.

Figure 1. The basic configuration with a uniform shear flow $u=y$ in the $x$ direction over a rigid obliquely orientated wavy wall.

We define the friction velocity $u_s=\sqrt {\nu {u_b^*}' (0) }$, where $\nu$ is the kinematic viscosity and $u_b^*(y^*)$ is the flow in the $x^*$ direction when $\epsilon =0$. We take $U^*={u_b^*}' (0)b$ as the unperturbed velocity a distance $b$ from the mean wall position as a typical velocity and define the Reynolds number

(2.1)\begin{equation} R = \frac {{U^*} b }{\nu}= {\frac{u^2_s {b^2}}{\nu^2}} =R_f^2 , \end{equation}

where $R_f$ is the friction Reynolds number. The wavelength $b$ of the wall is assumed to be small compared with a typical length scale of the flow in the wall-normal direction so that the basic state with $\epsilon =0$ can be taken to be a uniform shear flow in the $x$ direction. We scale time on ${1}/{{u_b^*}'(0) }$ so that the dimensionless equations of motion are

(2.2)$$\begin{gather} \frac{{\rm D}{\boldsymbol{ {u}}}}{{\rm D}t} ={-}\boldsymbol{\nabla} p+\frac{1}{R}\nabla^{2}{\boldsymbol{ {u}}}, \end{gather}$$

(2.3)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{{u}} =0. \end{gather}$$

Here, ${p}$ is the pressure and the Reynolds number is assumed large. In dimensionless variables the wall is defined by

(2.4)\begin{equation} y=2\epsilon \cos[ {\alpha x +\beta z } ] , \end{equation}

and since we have used the wall wavelength as a typical length scale we can take $\alpha ^2+\beta ^2=1.$ At the wall the fluid velocity must vanish and far from the wall the flow at leading order is a uniform shear flow in the $x$ direction. Thus, for $y\gg 1$, we require that ${\boldsymbol {u}} \sim (y,0,0),$ see figure 1 for a sketch of the configuration. Note that the problem formulated above applies to any predominantly uni-directional shear flow over a wavy wall with the global properties of the unperturbed state entering only through the friction velocity. The relationship of the problem formulated above to the finite problem where an upper wall exists will be commented on later.

2.1. Vortex wavelengths comparable to wall wavelength

Let us summarise the key relevant results from H2 where streamwise vortex disturbances of wavelength comparable to the wall wavelength are considered. In that case the instability is associated with a VWI with the wave forced by the wavy boundary. The disturbance field then splits into a two layer structure with a thin viscous layer and an $O(1)$ region adjacent to the wall. Away from the wall there is a streamwise vortex with wavelength ${2{\rm \pi} }/{K}$ driven by a VWI process in the wall layer. If we define the amplitude parameter

(2.5)\begin{equation} \varGamma_0 =\epsilon^2R^{{4}/{3}}\log R, \end{equation}

then the neutral value of $\varGamma _0$ for a disturbance with wavenumber $K$ is given by (3.5) of H2. If the vortex is not neutrally stable so that its streamwise velocity component $\sim u_0(y)\exp ({{\rm i}Kz+\sigma _0 t})$ then (3.5) of H2 becomes

(2.6)\begin{equation} \varGamma_0=\frac{2([1+K^2]^2-4\beta^2K^2)\sqrt{\sigma_0+K^2}(\sqrt{\sigma_0 +K^2} +K)}{ Ai'^2(0)\alpha^{4/3}K^2[(1+K^2)(10 +9\beta^2K^2 )-38\beta^2K^2 ]}+\cdots,\end{equation}

where $\alpha ^2+\beta ^2=1$ and $Ai$ is the Airy function. For given values of $K, \varGamma _0$ (2.6) is a quadratic equation for $\sqrt {\sigma _0+K^2}$ which can be solved to give the growth rate $\sigma _0$.

Figure 2(a) shows neutral values of $\varGamma _0$ as a function of $K$ for a range of values of $\alpha$. Note that $\alpha =1$ corresponds to the aligned case where the wave crests are perpendicular to the basic flow. We observe that for $K\gg 1$ the right-hand side of (2.6) with $\sigma _0=0$ tends to a constant so that, for any non-zero value of $\beta$, instability is possible for $0< K<\infty$ for a large enough value of $\varGamma _0$. In fact, the approach to that constant is from above or below depending on whether $\beta$ is greater or less than $\beta _c={(47-\sqrt {769})}/{72}\sim 0.53.$ Thus, for the unbounded problem, instability first occurs as $\epsilon ^2R^{4/3}\log R$ is increased at a finite value of $K$. Note that for the finite problem instability can occur first as $K\rightarrow \infty$; see H2.

Figure 2. (a) Values of $\varGamma _0$ as a function of $K$ for different values of the wavenumber $\alpha.$ Note that $\alpha =1, \beta =0$ corresponds to the aligned case. (b) The growth rate $\sigma _0$ as a function of $K$ for different values of $\alpha$ when $\varGamma _0=70$.

Figure 2(a) shows the growth rate $\sigma _0$ as a function of $K$ for the same range of values of $\alpha$ with $\varGamma _0=70.$ We see that the aligned case $\alpha =1$ gives instability over only a finite range of $K$ but $\alpha = 0.4,0.6,0.8$ are unstable for all $K$. The case $\alpha =0.2$ is unstable apart from a small interval near $K=1$. We see that for small $K$ the growth rate $\sigma _0\sim K$ whilst for $K\gg 1$ we find $\sigma _0\sim K^2$ so that the largest growth rates always correspond to the short-wave limit. The primary aim of the present paper is to determine the fate of the latter short-wave instability as $K$ increases. In doing so we will uncover a new instability mechanism which enables disturbances originating in the $K=O(1)$ problem to persist to indefinitely large $K$ for any non-zero $\beta$.

2.2. The modified basic state induced by the wavy wall

We shall in the first instance describe how, at large Reynolds numbers, the flow adjusts to the presence of the wavy wall in a viscous layer of thickness $R^{-{1/3}}$ next to the wall. The adjustment depends crucially on the size of this viscous layer compared with the wave amplitude. Here, we assume the amplitude is small compared with the viscous layer depth. Thus, we require $\epsilon \ll R^{-{1/3}}$; for a discussion of the flow at larger amplitudes see, for example, Smith (Reference Smith1982). It is a common misconception in the wavy-wall instability literature that flow reversal is needed for instability, but at the amplitudes discussed here there is no flow reversal when instability occurs. In terms of the wall layer variable $Y=R^{1/3}y$, the perturbed flow for $Y=O(1)$ is

(2.7)$$\begin{gather} {\boldsymbol{u}}={\boldsymbol{u_b}}=(R^{-{1/3}}Y,0,0)+[\epsilon E(u_1(Y), R^{-{1/3}} v_1(Y), w_1(Y)) +{\rm c.c.}]+\cdots, \end{gather}$$

(2.8)$$\begin{gather}p=p_b+R^{-{1/3}}[\epsilon Ep_1+{\rm c.c.} ]+\cdots. \end{gather}$$

Here, $E=\exp ({{\rm i}[\alpha x+ \beta z]})$, c.c. denotes the complex conjugate and $p_1$ is a constant. The order-$\epsilon$ correction above must satisfy $u_1(0)=-1, v_1(0)=w_1(0)=0$, and, because the flow in the outer $y=O(1)$ region satisfies Rayleigh's equation, which has no non-trivial solution with normal velocity vanishing at the wall, we require that $v_1'(\infty )=0$. Substituting the above into the equations of motion written in terms of the wall layer variable $Y$, following the steps given in H2 it is readily shown that

(2.9a,b)$$\begin{gather} v_1 ={-}3[{\rm i} \alpha]^{2/3}\int_0^{({\rm i} \alpha)^{1/3}Y}\,{\rm d}\theta\int_\infty^\theta Ai (\psi)\,{\rm d}\psi,\quad w_1 =\frac{3\alpha \beta Ai'(0)}{\alpha^2+\beta^2}{\mathcal{L}}([{\rm i} \alpha]^{1/3}Y) , \end{gather}$$

(2.10a,b)$$\begin{gather}p_1={-}\frac{3[{\rm i} \alpha]^{5/3}Ai'(0)}{\alpha^2+\beta^2},\quad {\rm i}\alpha u_1 ={-} \frac{{\rm d}v_1}{{\rm d}Y}-{\rm i}\beta w_1. \end{gather}$$

Here, $Ai$ is the Airy function and ${\mathcal {L}}(\psi )$ is the Scorer function which satisfies

(2.11)$$\begin{gather} {\mathcal{L}}''-\psi {\mathcal{L}} =1, \end{gather}$$

(2.12)$$\begin{gather}{\mathcal{L}}(0)={\mathcal{L}}(\infty) = 0. \end{gather}$$

Before deriving the equations governing the instability of the flow modified by the wall undulation we note that the wall wavelength has been used as the length scale in the problem and so the discussion above is valid for $O(1)$ wavelengths. If an upper boundary is present then constraints on the layer depth and wavelength must be applied; see H2.

2.3. The equations for near-wall disturbances

The instability of the wave-modified flow described above to streamwise vortex disturbances of wavelength comparable to the wall wavelength is discussed in H1 for $\beta =0$ and more generally in H2. In the oblique case $\alpha, \beta \ne 0$ we noted above that a short-wave instability occurs with the temporal growth rate increasing like the square of the vortex wavenumber at high wavenumbers. In H2 is is shown that, for flows with an upper wall, the short-wave instability can be the first to occur as the Reynolds number is increased at a given wall amplitude. Thus, in for example Couette flow where we must now satisfy the no-slip condition at $y=1$, the short-wave instability is the most dangerous one for oblique walls of wavelength greater than 5.83 channel depths. But in both the bounded and unbounded cases the largest growth rates are associated with the short-wave limit.

We first note that when $K$ is increased in (2.6) the associated disturbance field becomes localised in a layer of depth $K^{-1}$ near the wall. But the vortex in the $y=O(1)$ region is driven by a VWI process in the viscous wall layer of thickness $O(R^{-{1/3}})$ and so a new distinguished limit emerges when $K\sim R^{1/3}.$ In other words the short-wave limit of H2 must take on a new structure when the vortex wavelength scales with the wall layer thickness, we investigate that problem here. By subsequently taking the long-wavelength limit of the new distinguished limit we show how that situation connects with the short-wave limit of H2. We also discuss the short-wavelength limit of the new distinguished limit and find that the disturbance then takes on a self-similar form of increasing complexity.

The size of $\epsilon$ at which instability occurs in the new regime here can be inferred from H2. However, it can also be trivially obtained from well-known results for three-dimensional boundary layer instabilities. We will see that the instability can be interpreted as a modulated cross-flow vortex instability; see for example Gregory, Stuart & Walker (Reference Gregory, Stuart and Walker1955) or Hall (Reference Hall1986) in the context of the rotating disc boundary layer.

Cross-flow vortices occurring in three-dimensional boundary layers are inviscid Rayleigh modes of the composite velocity profile ${\hat {\alpha }} \bar {u} +{\hat {\beta }}{\bar {w}}.$ Here, $\hat {\alpha }, \hat {\beta }$ are the cross-flow vortex wavenumbers in the $x,z$ directions and $\bar {u}, \bar {w}$ are the boundary layer velocities in those directions. On swept wings, or the flow over a rotating disc, the latter velocities are of comparable size and instability is readily obtained by choosing the wavenumbers so that the composite velocity has an inflection point. However, as discussed by Bassom & Hall (Reference Bassom and Hall1991), cross-flow vortices persist when $\bar {w} \ll \bar {u}$. Within the wall layer (2.7) shows that the streamwise velocity is of size $R^{ -{1/3}}$ whilst the spanwise velocity is of size $\epsilon$. Therefore, if we perturb the flow in the wall layer to a disturbance with wavelength ${2{\rm \pi} }/{\alpha }$ in the $x$ direction and wavelength scaled on $R^{ -{1/3}}$ the wall layer thickness in the $z$ direction, then the composite velocity associated with the wall layer flow will have terms of comparable size if $\epsilon =O(R^{-{2/3}})$. The same scaling can be deduced from H2 by taking the long-wavelength limit of the problem considered there. Thus we anticipate a cross-flow vortex instability if $\epsilon =O(R^{-{2/3}})$, but the situation is slightly more complex than the usual scenario. Since we are seeking a disturbance in the wall layer with wavelength of size $R^{-{1/3}}$ in the $z$ direction, and the base flow depends on $Y$ and $\alpha x +\beta z$, the disturbance depends on the variables $\alpha x+\beta z, Y=R^{ {1/3}}y$ and $R^{ {1/3}}z$.

The wall layer is viscous and so, unlike the classical cross-flow vortex scenario, we cannot ignore viscosity. But the composite basic velocity profile here is a function of $\alpha x+ \beta z$ and so, unlike the usual cross-flow vortex case, the base flow and instability vary on the same streamwise length scale. Therefore a parallel flow approximation cannot be made but, since the base flow is periodic in $\alpha x+\beta z$, that is taken care of by appealing to Floquet theory. Here, unlike the VWI case, the leading-order stability problem involves all the streamwise harmonics. The upshot then is that the instability we describe is best interpreted as a viscous cross-flow instability of a spatially periodic flow. But we will see that when the vortex wavenumber $k$ scaled on the reciprocal of the wall layer thickness tends to zero, the instability reverts to the VWI type. However, for $k\rightarrow \infty$ a structure apparently not related to known instabilities develops. But we see that at all disturbance wavenumbers neither of the two often-assumed views that wavy-wall-induced instabilities are centrifugal in nature, or instabilities associated with reversed flow, are correct.

It is convenient to now introduce $\varPhi =\alpha x+\beta z$ and define $Z=R^{1/3}z.$ Thus, we seek a structure with the disturbance having a $z$ dependence through the slow variable $\varPhi$ and the fast variable $Z$. The disturbance varies only slowly in the $x$ direction so that the only $x$ dependence is through the slow variable $\varPhi$. Within the wall layer we perturb the basic flow given by (2.7)–(2.12) by writing

(2.13)$$\begin{gather} u=u_b+\Delta U^+(\varPhi,Y,Z)+\cdots, \end{gather}$$

(2.14)$$\begin{gather}v=v_b+\Delta R^{-{1/3}} V^+(\varPhi,Y,Z)+\cdots, \end{gather}$$

(2.15)$$\begin{gather}w=w_b+\Delta R^{-{1/3}} W^+(\varPhi, Y,Z)+\cdots, \end{gather}$$

(2.16)$$\begin{gather}p=p_b +\Delta R^{{-}1} P^+(\varPhi, Y,Z)+\cdots, \end{gather}$$

where the disturbance amplitude $\varDelta \ll 1$ and the basic state ${\boldsymbol {u}}_b, p_b$ is as defined by (2.7)–(2.12). If we take $\epsilon R^{2/3}=O(1)$ and substitute the above expansions into the equations of motion expressed in terms of wall layer variables, linearising in $\varDelta$, then the leading order approximation for $R\gg 1$ is

(2.17)$$\begin{gather} \left[\boldsymbol{\nabla}^2_2- \alpha Y \frac{\partial }{\partial \varPhi}-\epsilon R^{2/3} {\widehat{w_1}}\frac{\partial }{\partial Z}\right]U^+ = V^+ , \end{gather}$$

(2.18)$$\begin{gather}\left[\boldsymbol{\nabla}^2_2- \alpha Y \frac{\partial }{\partial \varPhi}-\epsilon R^{2/3} {\widehat{w_1}}\frac{\partial }{\partial Z}\right]V^+ = \frac{\partial P^+}{\partial Y} , \end{gather}$$

(2.19)$$\begin{gather}\left[\boldsymbol{\nabla}^2_2 - \alpha Y \frac{\partial }{\partial \varPhi}- \epsilon R^{2/3}{\widehat{w_1}}\frac{\partial }{\partial Z}\right]W^+ = \frac{\partial P^+}{\partial Z} +\epsilon R^{2/3} \left[V^+\frac{\partial {\widehat{w_1}}}{\partial Y}+\alpha U^+ \frac{\partial {\widehat{w_1}}}{\partial \varPhi}\right], \end{gather}$$

(2.20)$$\begin{gather}\alpha \frac{\partial U^+}{\partial \varPhi}+\frac{\partial V^+}{\partial Y} +\frac{\partial W^+}{\partial Z}=0, \end{gather}$$

where $\boldsymbol {\nabla }^2_2=[{\partial ^2}/{\partial Y^2 }+{\partial ^2}/{\partial Z^2}]$, ${\widehat {w_1}}=w_1E+{\rm c.c.}$ and $w_1$ is given by (2.9a,b). The streamwise wavenumber $\alpha$ can be scaled out of the above system by the transformation $(Y,Z,U^+,V^+,W^+,P^+)\rightarrow ({Y}/{\alpha ^{ {1/3}}} , {Z}/{\alpha ^{ {1/3}}}, U^+, \alpha ^{{2/3}}V^+, \alpha ^{ {2/3}}W^+,\alpha P^+)$ and defining the wall amplitude parameter $\varGamma$ by

(2.21)\begin{equation} \varGamma={\varGamma}(\alpha)= 3 \alpha^{2/3}[1-\alpha^2]^{1/2}\epsilon R^{2/3}Ai'(0)= 3\epsilon R^{2/3}Ai'(0)\sin \theta \cos^{2/3} \theta ,\end{equation}

where, as defined earlier, ${{\rm \pi} }/{2}-\theta$ is the angle between wave crests and the unperturbed flow direction. Note that if there exists a neutral solution for a given $\varGamma$ then there exists a corresponding solution for $-\varGamma$ since switching the sign of $\epsilon$ corresponds to a shift of ${\rm \pi}$ in $\varPhi$. We see that $\tan \theta =\sqrt {{3}/{2}}$ minimises the value of $\epsilon R^{2/3}$ needed for a given value of $\varGamma$ so, for a given shear flow and wall, the flow is most unstable when the wave crests are inclined at an angle of approximately $39^\circ$ to the flow direction. Thus, with the orientation of the wall built into the definition of $\varGamma$, (2.17)–(2.20) and the associated boundary conditions can be replaced by

(2.22)$$\begin{gather} \left[\boldsymbol{\nabla}^2_2- Y \frac{\partial }{\partial \varPhi}-\varGamma{\bar{w}_1}\frac{\partial }{\partial Z}\right] U^+ = V^+ , \end{gather}$$

(2.23)$$\begin{gather}\left[\boldsymbol{\nabla}^2_2-Y \frac{\partial }{\partial \varPhi}-\varGamma{\bar{w}_1} \frac{\partial }{\partial Z}\right] V^+ = \frac{\partial P^+}{\partial Y} , \end{gather}$$

(2.24)$$\begin{gather}\left[\boldsymbol{\nabla}^2_2- Y \frac{\partial }{\partial \varPhi}-\varGamma{\bar{w}_1}\frac{\partial }{\partial Z}\right]W^+ = \frac{\partial P^+}{\partial Z}+\varGamma\left[V^+\frac{\partial {\bar{w}_1}}{\partial Y}+ U^+\frac{\partial {\bar{w}_1}}{\partial \varPhi}\right], \end{gather}$$

(2.25)$$\begin{gather}\frac{\partial U^+}{\partial \varPhi}+\frac{\partial V^+}{\partial Y} +\frac{\partial W^+}{\partial Z}=0, \end{gather}$$

(2.26a,b)$$\begin{gather}U^+ = V^+ = W^+ = 0, \quad Y=0, \infty, \end{gather}$$

where $\bar {w}_1={\mathcal {L}}({\rm i}^{1/3}Y)E+{\rm c.c.}$ The above system is linear and the coefficients are independent of $Z$ so we can seek solutions proportional to ${\rm e}^{\pm {\rm i} k Z}$ in which case the system simplifies to a partial differential system in the variables $\varPhi, Y.$ The $\varPhi$ dependence is periodic and so the above system is similar to that governing the centrifugal instability of a Stokes layer with the Scorer function playing the role of the Stokes layer velocity profile; see Seminara & Hall (Reference Seminara and Hall1976). We now take the disturbance to be periodic in $Z$ and write

(2.27)\begin{equation} (U^+,V^+)={\boldsymbol{U}}^+(\varPhi, Y, Z)={\boldsymbol{U}}(\varPhi, Y)\,{\rm e}^{{\rm i} k Z}+{\rm c.c.}, \end{equation}

and eliminate the disturbance pressure and spanwise velocity component from the above equations to give the following system to determine ${\boldsymbol {U}}$:

(2.28)$$\begin{gather} \left[\frac{\partial^2}{\partial Y^2}-k^2 - Y \frac{\partial }{\partial \varPhi}-{\rm i}k\varGamma{\bar{w}_1} \right]U=V , \end{gather}$$

(2.29)$$\begin{gather}\left[\frac{\partial^2}{\partial Y^2}-k^2 - Y \frac{\partial }{\partial \varPhi}-{\rm i}k\varGamma{\bar{w}_1}\right] \left[\frac{\partial^2}{\partial Y^2}-k^2\right] V={-}2{\rm i}k \varGamma \frac{\partial}{\partial Y}\left[\frac{\partial {\bar{w}}_1}{\partial \varPhi} U \right] - {\rm i}k\varGamma \frac{\partial^2 {\bar{w}}_1}{\partial Y^2}V, \end{gather}$$

(2.30a,b)$$\begin{gather}U=V =\frac{\partial V}{\partial Y}=0,\quad Y= 0, \infty. \end{gather}$$

Notice that, if $\bar {w}_1$ is independent of $\varPhi$, we can solve for $V$ independently of $U$, which is then found from (2.28). That situation corresponds to the viscous Orr–Sommerfeld case of cross-flow vortex instabilities for the composite profile $\alpha ^+Y+{\rm i}k \varGamma w_1$ for a disturbance periodic in $\varPhi$ with wavenumber $\alpha ^+.$ In the present problem, instability arises from the first term on the right-hand side of (2.29) so the instability owes its origin to the periodicity in $\varPhi$ of the basic state.

The form of the $\varPhi$ dependence of ${\bar {w}}_1$, and the fact that wherever it appears above it is multiplied by ${\rm i}k$ means that (2.28)–(2.30a,b) admit Floquet solutions of the form

(2.31)$$\begin{gather} U= {\rm e}^{\mu \varPhi}\sum^\infty_{n={-}\infty}[U_{2n}(Y)E^{2n} \cos kZ+U_{2n-1}(Y)E^{2n-1} \sin kZ], \end{gather}$$

(2.32)$$\begin{gather}V= {\rm e}^{\mu \varPhi} \sum^\infty_{n={-}\infty}[V_{2n}(Y)E^{2n} \cos kZ+V_{2n-1}(Y)E^{2n-1} \sin kZ] , \end{gather}$$

where $\mu$ is the Floquet exponent. The corresponding form of $W$ is similar but with $\cos kZ$ and $\sin kZ$ interchanged. We can substitute the above Fourier expansions into (2.28)–(2.30a,b) to obtain an infinite system of coupled ordinary differential equations of similar form to those found by Seminara & Hall (Reference Seminara and Hall1976) for the Stokes layer problem. However, it can be seen from (2.29) for large $Y$ that, far from the wall, $V_n$ has solutions ${\rm e}^{ \pm {2/3}{\rm i}^{1/2}n^{1/2}y^{3/2}}$ and ${\rm e}^{\pm kY}$ so that the infinite system is stiff and not readily solved by shooting. In fact it was found that shooting works for small $k$ where the solution is dominated by the $n=0$ terms.

Instead of shooting we found solutions by marching (2.28)–(2.30a,b) in $\varPhi$ using a finite difference scheme until the solution settled down to the mode with the largest growth rate $\mu _r$. In the following section we present results obtained using that procedure, but we note here that all solutions found have $\mu$ purely real and have period $2{\rm \pi}$ in $\varPhi$.

3. Numerical results for $O(1)$ values of $k$ and $\varGamma$

The linear disturbance equations (2.28)–(2.30a,b) have coefficients periodic in $\varPhi$ and therefore have Floquet solutions of the form

(3.1)\begin{equation} (U,V)={\rm e}^{ \mu \varPhi }({{\bar{U}}}(Y,\varPhi), {{\bar{V}}}(Y,\varPhi)). \end{equation}

Here, ${{\bar {U}}}(Y,\varPhi )$ and ${{\bar {V}}}(Y,\varPhi )$ are periodic functions of $\varPhi$ and the Floquet exponent $\mu$ is a function of the spanwise wavenumber $k$ and $\varGamma$. For each $(k, \varGamma )$ we expect there will be an infinite discrete spectrum of eigenvalues $\mu$. If we march (2.28)–(2.29) forward in $\varPhi$ over a sufficiently long range the solution will be dominated by the Floquet solution with largest positive real part $\mu _r$. The system (2.28)–(2.30a,b) is parabolic in $\varPhi$ and is solved using a finite difference scheme; see for example Hall (Reference Hall1983) for an appropriate scheme first used for the Görtler vortex problem. For a given $(k, \varGamma )$ an initial disturbance is imposed at $\varPhi =0$ and then marched forward in $\varPhi$. The initial disturbance is chosen so that the solution does not have a singularity at the wall for small $\varPhi$.

Numerically, it was confirmed that, after integrating over a sufficiently long interval in $\varPhi$, the solution is dominated by the least stable Floquet mode. Once the solution has converged to the dominant Floquet solution the Floquet exponent is calculated by comparing the solution over an interval of $2{\rm \pi}$ in $\varPhi.$ Because most of the energy of the disturbance is in the streamwise velocity component the calculations reported on here used

(3.2)\begin{equation} I=\int_0^\infty\vert{{U}}(Y,\varPhi)\vert^2\,{\rm d}Y, \end{equation}

as the disturbance property used to find $\mu.$ We found that the dominant Floquet mode was always real and synchronous with the forcing in $\varPhi$ so that $\mu$ is real and given by

(3.3)\begin{equation} \mu=\frac{1}{4{\rm \pi}}[\log{I(\varPhi_0+2{\rm \pi})}-\log I({\varPhi_0})] \end{equation}

for some suitably large $\varPhi _0$. Typically, we found it necessary to integrate over some 20–30 periods in $\varPhi$ in order to obtain $\mu$ correct to two decimal places. Alternative disturbance properties such as the maximum disturbance amplitude and the wall shear were also used to monitor disturbance growth and gave consistent results. However, for reasons which will become apparent later when we give results, the wall shear is not an appropriate disturbance measure at higher wavenumbers.

In the calculations it is found that the value of $Y=Y_\infty$ needed to approximate $\infty$ depends on the spanwise wavenumber $k$. For small $k$, the disturbance develops a two layer structure in $Y$ with an $O(1)$ layer at the wall and an outer layer of depth $O(k^ {-1})$; see § 4 for details of this structure. Thus in the computations it is necessary to take $kY_\infty \gg 1$ whilst keeping the step length in $Y$ small enough to resolve the $O(1)$ layer at the wall. Typically, it was found necessary to use $2000$ points in $Y$ with $Y_\infty =100$. For the smaller values of $k$ reported below it was necessary to increase the number of points to $10^4$ with infinity at $500$. Calculations at higher values of $k$ showed that the disturbance in that limit localises in a layer which moves progressively further away from the wall and develops an increasingly rapid variation with respect to $\varPhi$. The disturbance structure for $k\gg 1, \varGamma \gg 1$ will be discussed in § 5 but we note here that to capture it the numerical scheme needs to be modified. The analysis given in § 5 shows that the neutral disturbance for large $\varGamma$ is exponentially small apart from a layer of depth $\varGamma ^{1/12}$ located a distance $\varGamma ^{1/3}$ from the wall. Therefore, the numerical scheme used for large values of $\varGamma$ employed a grid taking account of that structure.

Figure 3(a) shows the Floquet exponent $\mu$ as a function of $k$ for different values of $\varGamma.$ We see that, as expected, viscous effects stop the monotonic increase with wavenumber of the growth rate found in H2 for disturbances with wavenumber comparable to the wavelength of the undulating wall. For a given value of $k$ we vary $\varGamma$ until the Floquet exponent vanishes; by then varying $k$ a neutral curve in the $k\unicode{x2013}\varGamma$ plane is constructed.

Figure 3. (a) The growth rate as a function of $k$ for different values of $\varGamma$. (b) The neutral curve where $\mu =0$ in the $k\unicode{x2013}\varGamma$ plane and the small $k$ approximation to the neutral curve described in § 4.

That neutral curve is shown in figure 3(b). Also shown in the figure is a small $k$ asymptotic prediction of the neutral curve to be derived in § 4. We see there is rather good agreement between the small $k$ and full numerical solutions when $k\rightarrow 0.$ For the larger values of $k$ shown $\varGamma$ increases; later we will see that for $k\gg 1$ the neutral curve has $k\sim \varGamma ^{1/6}.$ That limit will be discussed in § 5 and further results given there.

Figure 3(a) shows that the most amplified disturbance moves to smaller $k$ as $\varGamma$ increases. That trend continues at even higher values of $\varGamma$ but we do not pursue it further here.

Now let us discuss how the structure of the disturbance changes along the neutral curve. In order to see the changes we plot contours of constant values of the disturbance functions $(\bar {U},\bar {V})$ in the $\varPhi \unicode{x2013}Y$ plane. We also plot isosurfaces of the streamwise vorticity

(3.4)\begin{equation} \omega_x=\left[\frac{\partial^2}{\partial Y^2}-k^2\right]{\bar{V}}+\frac {\partial^2{\bar{U}}}{\partial \varPhi \partial Y}. \end{equation}

Figure 4 shows contours of constant values of $\bar {U}, \bar {V},$ in the $\varPhi \unicode{x2013}Y$ plane and isosurfaces of $\omega _x$ for the neutral case with $\varGamma =1.38.$ For this case the velocity field is predominantly real and independent of $\varPhi$ with a small imaginary part varying with $\varPhi.$ That is consistent with the small $k$ solution described in the next section or the large wavenumber limit of H2. The vorticity isosurfaces where the vorticity is $\pm \%60$ of its maximum value are illustrated in the lower part of figure 4. We note that $\varPhi$ is defined by scaling $x^*,z^*$ with the wall wavelength $b$ whereas $Y$ is scaled on $b R^ {-{1/3}}$ so that in dimensional coordinates the humps in the figure will be elongated in the streamwise direction. Once again the isosurfaces are consistent with the small $k$ structure to be described in § 4. Figure 5 gives the corresponding results further to the right along the neutral curve where $\varGamma =2.55.$ Here, the isosurfaces of the vorticity have not changed significantly but the real and imaginary parts of the velocity field have become comparable to a significant variation in $\varPhi$ having developed.

Figure 4. Results for $\varGamma =1.38$ (a–d) shown across the real and imaginary parts of $\bar {U}, \bar {V}$ in the first and second rows, respectively. (e) The $\pm 60\,\%$ isosurfaces of the streamwise vorticity for $\varGamma =1.38$.

Figure 5. Results for $\varGamma =2.55$ (a–d) shown across the real and imaginary parts of $\bar {U}, \bar {V}$ in the first and second rows, respectively. (e) The $\pm 60\,\%$ isosurfaces of the streamwise vorticity for $\varGamma =2.55$.

Figures 6 and 7 show results for the neutral cases with $\varGamma =10,100$ and here we see a distinct new structure emerging. The velocity contours in figure 6 have begun to move away from the wall and figure 7 clearly shows the velocity field now localising away from the wall. The form of the isosurfaces seen in figures 4 and 5 is shown in figure 6 to be changing significantly before now taking on a new distinctive structure in figure 7. Calculations at higher values of $\varGamma$ reveal that the layer in which the disturbance localises moves further away from the wall and thickens slightly as $\varGamma$ increases. In addition a short-scale dependence in $\varPhi$ develops, we will discuss that structure further in § 5 where we take the limit $\varGamma \rightarrow \infty$. Surprisingly, we will find an analytical form for the leading-order asymptotic approximation to the increasingly complicated short-scale structure of the disturbance.

Figure 6. Results for $\varGamma =10.$ (a–d) Results across the real and imaginary parts of $\bar {U}, \bar {V}$ in the first and second rows, respectively. (e) The $\pm 60\,\%$ isosurfaces of the streamwise vorticity for $\varGamma =10$.

Figure 7. Results for $\varGamma =100.$ (a–d) Results across the real and imaginary parts of $\bar {U}, \bar {V}$ in the first and second rows, respectively. (e) The $\pm 60\,\%$ isosurfaces of the streamwise vorticity for $\varGamma =100$.

Let us now consider how the wall shear associated with the disturbance varies along the neutral curve. We find that

(3.5)\begin{equation} \mu_0= {\frac{\partial U^+ } { \partial Y}}_{\vert_Y=0}=\left[\frac{\partial}{\partial Y}\left(\bar{U}(\varPhi=\alpha x+\beta z,Y)\right) {\rm e}^{{\rm i} R^{1/3}z} \right] _{\vert_Y=0}. \end{equation}

Because of the two length scales in $z$ appearing in the expression for $\mu _0$ we must specify a value for the Reynolds number to evaluate $\mu _0$ as a function of $x$ and $z$. Before doing so, we give results for $\mu _0$ along $z=0$, these results will illustrate why $\mu _0$ as a function of $x,z$ becomes somewhat complicated for large $\varGamma.$

Figure 8 gives $\mu _0$ along $z=0$ for $\varGamma =10$ at different values of $k$. We observe the increasing variation with $\varPhi$ of the shear as $k$ increases. Note that, here and in the following, shear plots the disturbance have been normalised such that $\int ^\infty _0\vert \bar {U}\vert ^2 \,{\rm d}Y=1.$ Figure 9 shows the shear at the wall along $z=0$ for an increasing set of relatively large values of $\varGamma.$ We see that, as $\varGamma$ increases, the shear at the wall decreases rapidly and takes on a characteristic shape. In fact we will see in § 5 that, for large $\varGamma$, the disturbance decays exponentially away from the location in $Y$ where it localises, thus the disturbance is exponentially small near the wall; this is the reason why the shear decreases rapidly with $\varGamma$.

Figure 8. The real and imaginary parts of the shear $\mu _0=\bar {U}_Y(0,\varPhi )$ along $z=0$ as a function of $\varPhi$ for $\varGamma =10$ and, from top to bottom, $k=0.1,0.5,0.98$.

Figure 9. The real and imaginary parts of the shear $\mu _0=\bar {U}_Y(0,\varPhi )$ along $z=0$ as a function of $\varPhi$ on the neutral curve for, from top to bottom, $\varGamma =10,100, 200$.

Figure 10 shows the wall shear in the $x\unicode{x2013}z$ plane for different values of $\varGamma$ along the neutral curve. The values of $\alpha, \beta$ have been chosen to correspond to the most dangerous orientation with the crests at an angle of approximately $39^\circ$ to the basis state. Moving along the wave crests corresponds to fixing $\varPhi$ and varying $R^{1/3}z$ so we observe the fast $z$ dependence arising from the dependence of $U^+$ on $Z$.

Figure 10. Contours of constant values of the streamwise wall shear in the $x\unicode{x2013}z$ plane associated with the disturbance at values of $k, \varGamma.$ Note the horizontal axis is $x$ and the vertical axis is $z$. (a,b) Results for $\varGamma =2.1, 4.5$ and (c,d) results for $\varGamma =15, 30$. The value of $R$ used is $10^3$ in all cases.

Figure 10(a,b) corresponds to $\varGamma =2.1, 4.5.$ For the smallest value of $k$ we see a slight waviness between cell boundaries. However, in panel (b) we see that the slight cell boundary waviness seen in the first plot has now magnified significantly. The waviness is similar to the wavy vortex pattern which appears in the Taylor vortex problem; see for example Davey, DiPrima & Stuart (Reference Davey, DiPrima and Stuart1968). The waviness arises for small $k$ because, as can be seen in figure 4 and will be demonstrated asymptotically in the following section, the functions $\bar {U}, \bar {V}$ have dominant parts independent of $\varPhi$ with small imaginary parts which are wavy in $\varPhi.$ In panel (e) the values of $\varGamma =15,30$ and so the disturbance has begun concentrating away from the wall so that the wall shear falls significantly. Clearly, a distinct large $\varGamma$ pattern has already emerged for these relatively moderate values of $\varGamma.$ However, although that pattern persists to much larger values of $\varGamma$, its magnitude falls exponentially to zero and so the presence of the instability at large $\varGamma$ perhaps cannot be detected by measurements near the wall.

Note that the streamwise vorticity isosurfaces shown earlier in the $\varPhi -Y-Z$ space can in principle be similarly used to produce isosurfaces for a given $R$ in the $x-Y-z$ space. However, doing so produces patterns as shown already in the direction perpendicular to the wave crests but to extend the isosurfaces of the previous figures, each of which needed $500^3$ points to resolve the figure, is not possible over many wave crests.

4. The small wavenumber–small roughness limit

Now let us investigate the nature of the solutions of (2.28)–(2.30a,b) in the limit of long spanwise wavelengths. For small $k$ the wall layer develops a double layer structure with regions of depth $O(1)$ and $O(k^{-1})$. The key feature of the solution in this regime is that the Fourier expansions (2.31)–(2.32) now become power series in $k$ with the leading-order term independent of $\varPhi$, the second proportional to ${\rm e}^{\pm {\rm i} \varPhi }$ and so on. In the outer region the operators on the left-hand sides of (2.28)–(2.29) are dominated by the inviscid terms proportional to $Y$ and so the $\varPhi$-dependent part of the disturbance is inviscid in nature whilst the $\varPhi$ independent part, i.e. the roll-streak flow, is viscous.

4.1. The outer layer solution

In the outer region we define $\psi =kY$ and seek a solution of the form

(4.1)$$\begin{gather} U= [ {\hat{U}}_0 {} + \cdots]\cos kz+\{[k^3{\tilde{U}}_0 {} +\cdots]E+{\rm c.c.}\} \sin kz+\cdots, \end{gather}$$

(4.2)$$\begin{gather}V=[{k^2\hat{V}}_0 {} + \cdots]\cos kz+\{[k^2{\tilde{V}}_0 {} +\cdots]E+{\rm c.c.}\}\sin kz+\cdots, \end{gather}$$

(4.3)$$\begin{gather}W=[k^2{\hat{W}}_0 {} + \cdots]\sin kz+\{[k^2{\tilde{W}}_0 {}+ \cdots]E+{\rm c.c.}\} \cos kz+\cdots, \end{gather}$$

(4.4)$$\begin{gather}P=[{k^3\hat{P}}_0 {} + \cdots]\cos kz+\{[{\tilde{P}}_0 {} +\cdots]E+{\rm c.c.}\}\sin kz +\cdots, \end{gather}$$

where once again $E={\rm e}^{{\rm i}\varPhi }$ and ${\hat {U}}_0 {} , {\hat {U}}_1,$ etc. are functions of only the outer variable $\psi =kY.$ Note that the disturbance in the outer layer at leading order is a streamwise vortex with no modulation in $\varPhi$. The relative size of the terms proportional to $E$ above is fixed by the interaction in the wall layer. Substituting the above into (2.28)–(2.29) written in terms of $\psi$ we find that the leading-order problem satisfies

(4.5a,b)$$\begin{gather} \left[\frac{{\rm d}^2}{{\rm d} \psi^2}-1\right]{\hat{U}}_0= {\hat{V}}_0,\quad \left[\frac{{\rm d}^2}{{\rm d} \psi^2}-1\right]{\hat{V}}_0= \frac{{\rm d}{\hat{P}}_0}{{\rm d}\psi}, \end{gather}$$

(4.6a,b)$$\begin{gather}\left[\frac{{\rm d}^2}{{\rm d}\psi^2}-1\right]{\hat{W}}_0 ={-} {\hat{P}}_0, \quad \frac{ {\rm d}{\hat{V}}_0}{{\rm d}\psi}+{\hat{W}}_0=0. \end{gather}$$

The above linear equations are solved subject to

(4.7a–d)\begin{equation} {\hat{U}}_0={\hat{V}}_0=0,\quad \psi=0,\quad {{\hat{\boldsymbol{U}}_0}}\rightarrow 0, \quad \psi \rightarrow \infty. \end{equation}

Here, we are anticipating that, as discussed in H2, the spanwise vortex of the vortex will be driven by an interaction in the wall layer so we do not impose the condition ${\hat {W}}_0(0)= 0$. Normalising the solution such that ${\hat {W}}_0(0)=-1$ we find that

(4.8a–c)\begin{equation} {\hat{U}}_0={-}\frac{[\psi^2+\psi]}{4}{\rm e}^{-\psi},\quad {\hat{V}}_0=\psi {\rm e}^{-\psi}, \quad {\hat{W}}_0=(\psi-1){\rm e}^{-\psi}. \end{equation}

The next-order approximation to (2.28)–(2.29) determines the leading-order pressure and velocity fields proportional to $E, E^{-1}$ in the outer layer. As mentioned earlier, viscosity is not important for the latter field and so ${\tilde {V}}_0$ satisfies Rayleigh's equation with wavenumber unity and a uniform shear flow in $\psi$. It follows that ${\tilde {V}}_0 =Be^{-\psi }$ there with $B$ a constant fixed by matching with the flow in the wall layer. Note here that an unbounded uniform shear flow is inviscidly stable and so Rayleigh's equation has no solutions with the perturbation vanishing at the wall. Now we consider the flow in the wall layer and construct a solution consistent with (4.8a–c) in the limit $\psi \rightarrow 0$; that will determine the zeroth-order approximation to the neutral curve for $k$ small.

4.2. The inner layer solution for $Y=O(1)$

Here, we seek a solution of (2.28)–(2.29) taking the form

(4.9)$$\begin{gather} U=[k\lambda_0 Y {} + \cdots]\cos kz+\{[k^2{\tilde{u}}_0 {} +\cdots]E+{\rm c.c.}\} \sin kz+\cdots, \end{gather}$$

(4.10)$$\begin{gather}V=[{k^3\hat{v}}_0 {} + \cdots]\cos kz+\{[k^2{\tilde{v}}_0 {} +\cdots]E+{\rm c.c.}\}\sin kz+\cdots, \end{gather}$$

(4.11)$$\begin{gather}W=[k^2{\hat{w}}_0 {} + \cdots]\sin kz+\{[k {\tilde{w}}_0 {} +\cdots]E+{\rm c.c.}\} \cos kz+\cdots, \end{gather}$$

(4.12)$$\begin{gather}P=[{k^3\hat{p}}_0 {} + \cdots]\cos kz+\{[ {\tilde{p}}_0 {} +\cdots]E+{\rm c.c.}\}\sin kz +\cdots. \end{gather}$$

Here, $\lambda _0={\hat {U}}'_0(0)$ and the sizes of the first terms on the right-hand sides of the above expansions are fixed by the limiting form of the outer expansion as the wall is approached. In particular, the first term in the expansion of $U$ is the first term in the Taylor series expansion near the wall of the corresponding term in the outer expansion. Also, the first term on the right-hand size of the expansion of $W$ must match onto the corresponding term in the outer flow which, as noted above, cannot be reduced to zero as the wall is approached. The leading-order terms proportional to $E^{\pm 1}$ are driven by the terms proportional to $\varGamma$ in (2.28)–(2.29) and that fixes the scaling used in the expansions. Substituting the above expansions (2.28)–(2.29) we find that the leading-order system for terms proportional to $E$ is

(4.13)$$\begin{gather} \left[\frac{{\rm d}^2}{{\rm d} Y ^2}-{\rm i}Y\right]{\tilde{u}}_0 = {\tilde{v}}_0-\lambda_0Y\varGamma {\mathcal{L}}({\rm i}^{1/3}Y), \end{gather}$$

(4.14)$$\begin{gather}\frac{{\rm d} {\tilde{p}}_0}{{\rm d}Y}=0, \end{gather}$$

(4.15)$$\begin{gather}\left[\frac{{\rm d}^2}{{\rm d}Y^2}-Y\right]{\tilde{w}}_0 = {\tilde{p}}_0+{\rm i}\lambda_0 Y\varGamma {\mathcal{L}}({\rm i}^{1/3}Y), \end{gather}$$

(4.16)$$\begin{gather}{\rm i}{\tilde{u}}_0+ \frac{{\rm d}{\tilde{v}}_0}{{\rm d}Y}- {\tilde{w}}_0 =0. \end{gather}$$

The above equations must be solved subject to ${\tilde {u}}_0= {\tilde {v}}_0= {\tilde {w}}_0=0$ at $Y=0$ and ${\tilde {u}}_0={\partial {{\tilde {v}}_0}}/{\partial Y}= {\tilde {w}}_0=0$ at $Y=\infty$. We obtain

(4.17a,b)\begin{align} {\tilde{{\boldsymbol{u}}}}_0={-}\lambda_0\varGamma \left({\rm i}\left[\frac{{\mathcal{L}}}{2}-\frac{{\mathcal{L}}''}6\right], -\frac{1}{2{\rm i}^{1/3}}[{\mathcal{L}}''-1],-\left[\frac{{\mathcal{L}}}{2}+\frac{{\mathcal{L}}'''}{3}\right]\right),\quad {\tilde{p}}_0={-}\frac{3\lambda_0 \varGamma {\rm i}^{2/3}}{2}, \end{align}

where ${\mathcal {L}}({\rm i}^{1/3}Y)$ is the Scorer function. Now let us consider the leading-order terms independent of $E$ when (4.9)–(4.12) are substituted into (2.28)–(2.29). In fact, it is sufficient here to consider only the spanwise momentum equation, which yields

(4.18)\begin{equation} \frac{{\rm d}^2 {\hat{w}}_0}{{\rm d}Y^2}={-}\varGamma \left[{2 {\mathcal{L}}^* {\tilde{w}}_0}-\frac{{\rm d}}{{\rm d}Y}({\mathcal{L}}^*{\tilde{v}}_0)\right]+{\rm c.c.}, \end{equation}

where $*$ denotes the complex conjugate. For large $Y$ we can use the above expression for ${\tilde {{\boldsymbol {u}}}}_0$ to show that, in that limit, the equation becomes

(4.19)\begin{equation} \frac{{\rm d}^2 {\hat{w}}_0}{{\rm d}Y^2}\simeq{-}\frac{ 3\varGamma^2\lambda_0 }{ Y^2}. \end{equation}

Replacing $\lambda _0$ by ${\hat {U}}'_0(0)$ it follows for large $Y$ that

(4.20)\begin{equation} {\hat{w}}_0 \simeq 3\varGamma^2 {\hat{U}}'_0(0) \log Y . \end{equation}

If we now replace $\log Y$ by $(\log \psi -\log k)$ it follows that the $\varPhi$ independent spanwise velocity in the wall layer is consistent with the limiting form of (4.8a–c) for small $\psi$ if

(4.21)\begin{equation} \varGamma^2=\frac{4 }{3 \vert \log k \vert } .\end{equation}

The small wavenumber prediction given above is shown in figure 2(b) and we see good agreement with the numerical solution. Equation (2.21) then gives

(4.22)\begin{equation} \epsilon^2 R^{4/3} = \frac{4 } {27\alpha^{4/3}\beta^2 \vert [Ai '(0)]^2 \log k \vert }.\end{equation}

The above result applies to the long-wave limit of the neutral curve for streamwise vortices of wavelength scaled on the viscous wall layer. Now compare the above result with the large $K$ form of (2.6) derived in H2 for disturbances of wavenumber $K$ where the vortex wavelength is scaled on the wall wavelength. It follows that $k=R^{1/3}K$ so expressing (4.22) in terms of $K$ and taking $\log K\ll \log R$ we obtain

(4.23)\begin{equation} \epsilon^2 R^{4/3}\log R = \frac{4 \left(1+3\dfrac{\log K}{\log R}+\cdots\right)}{9\alpha^{4/3}\beta^2 \vert [Ai '(0)]^2 \vert,}\end{equation}

and we find that the leading-order term of (2.6) expanded for $K\gg 1$ coincides with the leading-order term in the above expansion. Moreover, the order-$K^{-2}$ correction term in (2.6) expanded for large $K$ comes from the structure within the viscous wall layer whilst the correction term in (4.23) comes from the outer layer. Thus the two expansions can be combined to give the result

(4.24)\begin{equation} \epsilon^2 R^{4/3}\log R = \frac{4 \left(1+3\dfrac{\log K}{\log R}+\dfrac {1}{9\beta^2K^2}\left\{47-36\beta^2-\dfrac{10}{\beta^2}\right\}+\cdots\right)} {9\alpha^{4/3}\beta^2 \vert [Ai '(0)]^2 \vert }. \end{equation}

The term in the curly brackets is negative if $\beta <\sim 0.52$ and positive otherwise, so for a disturbance of wavelength small compared with the wall wavelength the amplitude needed for a give large Reynolds number has a local minimum at $K=K_c$ where

(4.25)\begin{equation} K_c^2 \simeq \frac{2\log{R} \left\{47-36\beta^2-\dfrac{10}{\beta^2}\right\}}{27\beta^2\alpha^{4/3}}, \end{equation}

if $\beta \gtrsim 0.52$ whilst for $\beta \lesssim 0.52$ the required amplitude increases monotonically. We comment further on this result in the conclusion.

5. The large roughness limit $\varGamma \gg 1$

We saw in § 3 that, for large numerical values of $\varGamma$, a distinctive disturbance structure concentrated at increasingly large distances from the wall emerges. Of course, $\varGamma \rightarrow \infty$ is a second limit having first taken $\epsilon R^{2/3}=O(1)$, as shown in (2.21). More precisely, the second limit must be taken such that $\epsilon R^{1/3}$ remains small so that the interactive region is not established. Here, we explain the emergence of the characteristic scales of the disturbance in that limit and give an asymptotic description of the disturbance. Calculations showed that, for large $\varGamma$, the disturbance extended a large distance from the wall but varied on a shorter length scale in the $Y$ direction. Therefore, we anticipate that the disturbance should be described by the so-called WKBJ procedure due to Wentzel, Kramers, Brillouin and Jeffreys. If the disturbance extends a distance $\varGamma ^\mu$ from the wall we assume a WKBJ expansion which has $U, V$ in (2.28)–(2.29) of the form

(5.1)\begin{equation} U,V\sim \exp({\varGamma^m \varTheta(\eta,\varPhi)}), \end{equation}

where $\eta =\varGamma ^{-\mu }Y$ is a slow variable in the $Y$ direction and $\varTheta (\eta,\varPhi )$ is a phase function. At this stage, the constants $\mu, m$ and the size of $k$ are unknown. For large $\varGamma$ inviscid effects dominate and so the dominant balance in (2.28)–(2.29) will be between the operators $Y({\partial }/{\partial \varPhi })$ and ${\rm i}k\varGamma \bar {w}_1.$ Noting that, for large $Y$ the function $\bar {w}_1\sim Y^{-1}$, it follows that

(5.2)\begin{equation} \varGamma^{\mu+m} \sim k\varGamma^{1-\mu}. \end{equation}

Next, we note that if diffusions in the $Y,Z$ directions are comparable then

(5.3)\begin{equation} k\sim\varGamma^{m-\mu}. \end{equation}

We deduce from the two above equations that $\mu =\frac {1}{3}$, however, $m$ remains unknown but once $m$ is known $k$ follows using $k\sim \varGamma ^{m-{1/3}}$. Now, consider the dominant balances in (2.28)–(2.29) when $Y({\partial }/{\partial \varPhi })+\varGamma \bar {w}_1({\partial }/{\partial Z})$ vanishes at leading order. In (2.28) this gives $k^2U\sim V$, whilst (2.29) requires $k^4V\sim k^2 \varGamma ^{2/3}$ so that $k\sim \varGamma ^{1/6},$ and then $m=\frac {1}{2}$. Thus we define

(5.4)\begin{equation} k= {\varGamma^{1/6}}k^{+} + {\varGamma^{-{1/3}}} k^{+{+}}+\cdots \end{equation}

and we are led to the following WKBJ expansion of $(U,V)$:

(5.5)\begin{align} \left(U,\frac{V}{ \varGamma^{1/3}}\right)&=[(U_0(\varPhi, \eta), V_0(\varPhi, \eta)+\varGamma^{-{1/2}}(U_1(\varPhi, \eta), V_1(\varPhi, \eta)+\cdots]\nonumber\\ &\quad \times \exp\left({\varGamma^{1/2}\varTheta(\varPhi, \eta)}\right). \end{align}

The leading-order approximations to (2.28) and (2.29) both require that

(5.6)\begin{equation} \eta \varTheta_{\varPhi}-2{\rm i} k^+ \eta^{{-}1}\cos\left[\varPhi- \frac{\rm \pi}{6}\right]=0. \end{equation}

It follows that

(5.7)\begin{equation} \varTheta = \frac {2{\rm i} k^+}{\eta^2}\sin\left[\varPhi-\frac{ {\rm \pi}}{6}\right]+\theta(\eta), \end{equation}

where $\theta$ is to be determined. At next order we obtain the following pair of equations for $(U_0,V_0)$:

(5.8)$$\begin{gather} \left[ \eta\frac{\partial}{\partial \varPhi}+{k^+}^2-S^2(\varPhi,\eta)\right]U_0 ={-}V_0, \end{gather}$$

(5.9)$$\begin{gather}\left[\eta\frac{\partial}{\partial \varPhi}+{k^+}^2-S^2(\varPhi,\eta)(S^2-{k^+}^2)V_0\right]= \frac{4{{\rm i}k^+}}{\eta}\sin\left[\varPhi- \frac{ {\rm \pi}}{6}\right]SU_0. \end{gather}$$

Here, $S$ is defined by

(5.10)\begin{equation} S=\left(\theta'-\frac{4{\rm i}k^+}{\eta^3}\sin \left[\varPhi- \frac{ {\rm \pi}}{6}\right]\right) . \end{equation}

The above equations can be combined to give a second-order ordinary differential equation with coefficients periodic in $\varPhi$ and we require solutions of that equation which are periodic in $\varPhi$. The quantity $\theta '$ plays the role of a Floquet exponent in that equation and for given values of $\eta$ and $k^+$ a complex value of $\theta '$ must be chosen so that a periodic solution can be found. Of course, for any $(k^+,\eta )$ there will be a number of possible values of $\theta '$ and so (5.5) should be summed over the different values in order to satisfy the boundary conditions at $0, \infty$. The system (5.8)–(5.9) has the property that if $\theta '$ is a real eigenvalue then $-\theta '$ is also an eigenvalue and for our purpose here we concentrate on solutions having $\theta '$ real, complex eigenvalues exist but they are not relevant in the large $\varGamma$ limit.

Figure 11(a) shows some real eigenvalues $\theta '>0$ for different values of $k^+$ over a range of values of $\eta$. We see that when $\eta$ decreases the eigenvalue decreases towards a minimum before becoming large at small vales of $\eta$. When $k^+$ is decreased then the minimum approaches zero and when $k^+=\bar {k} = 0.797407,\eta =\bar {\eta } = 1.73885$ we find that $\theta '$ vanishes as the eigenvalues $\pm \theta '$ coalesce. Figure 11(b) illustrates the behaviour of the eigenvalues just before coalescence. Thus the WKBJ solution has a second-order turning point at $\eta =\bar {\eta }$ for $\hat {k}=\bar {k}$ and therefore the solution will localise there. A similar structure is observed in the small wavelength Görtler vortex problem; see Hall (Reference Hall1982).

Figure 11. (a) The function $\theta '$ as a function of $\eta$ for different values of the wavenumber $k^+.$ (b) An enlargement of panel (a) and the corresponding negative root near the point of coalescence. (Black, blue, green, red correspond to $k^+=0.7977, 0.7976, 0.7975 , 0.797407$, respectively.)

The localised solution is found by constructing a transition layer connecting the branch with $\theta '\sim -(\eta -\bar {\eta })$ to the right of $\eta =\bar {\eta }$ with the same branch to the left of $\bar {\eta }$. That connection produces a localised solution since exponential decay is achieved when $\eta \rightarrow \pm \infty$; see Hall (Reference Hall1982) for a detailed discussion of the related Görtler problem.

It can be shown that the transition layer is of thickness $\varGamma ^{1/12}$ so the appropriate transition layer variable is $\zeta =\varGamma ^{-{1/12}}(Y-\bar {Y})$, where $\bar {Y}=\varGamma ^{1/3}\bar {\eta }.$ Within the transition layer the disturbance amplitude satisfies the parabolic cylinder equation and, written in terms of $Y$, the leading eigenfunction is simply ${\rm e}^{{-c\varGamma ^{-{1/6}}} (Y-\bar {Y})^2}$ where $c={\theta '(\bar {\eta } )}/{2}$, which can be estimated from figure 11(b). It follows that, as usual, at the second-order turning point the solution in terms of the transition layer variable is exactly the asymptotic form predicted by the WKB solution to the left and right of the turning point.

Figure 12(a) shows the functions $U_0, V_0$ as functions of $\varPhi$ and panel (b) shows a comparison of the neutral value of $k$ computed from the full equations with the first term of the large $\varGamma$ result

(5.11)\begin{equation} k=0.797407\varGamma^{1/6}+ {\varGamma^{-{1/3}}} k^{+{+}}+\cdots.\end{equation}

We see there is excellent agreement between the asymptotic and exact results. The second term in (5.11) is fixed by requiring exponential decay for large $\vert \zeta \vert$ of the solution in the transition layer. More precisely, it is found that, within the transition layer, the eigenfunction will satisfy the equation

(5.12)\begin{equation} \frac{{\rm d}^2 f}{{\rm d}\zeta^2}-\left[a-\frac{c\zeta^2}{4}\right]f=0, \end{equation}

where $a$ is proportional to $k^{++}.$ The equation for $f$ must be solved solved with $f\rightarrow 0, \zeta \rightarrow \pm \infty$ so that there is an infinite sequence of eigenvalues $a=[-n-\frac {1}{2}]c^{1/2}$ for $n=0,1,2,\ldots$.

Figure 12. (a) The eigenfunctions $U_0, V_0$ as functions of $\varPhi$ for the turning point values of $k$ and $\eta.$ (b) A comparison of the large $\varGamma$ asymptotic prediction with the finite $\varGamma$ calculations.

The least stable mode corresponds to $n=0$ and that is the mode which connects with our finite $\varGamma$ calculations. However, our analysis of the limit $k\ll 1$ showed that there is a single neutral value of $\varGamma \ll 1$ in that limit which connects with the single neutral mode found in H2. It follows that the higher-order eigenvalues of the large $\varGamma$ solution must connect with neutral curves which do not connect with solutions of the full problem with vortex wavenumbers comparable to the wall wavelength.

We recall that, for finite $k, \varGamma$ we march (2.28)–(2.29) in $\varPhi$ until the solution is dominated by the least stable mode. That procedure produces a neutral curve which, from figure 12, is seen to connect with (5.11). All the calculations carried out in the unstable regime suggest that higher modes present for large $k$ are always less unstable than the mode corresponding to the neutral curve of figure 2(b). Thus, the fate of the higher modes found at large $k$ requires the calculation of all eigenvalues of (2.28)–(2.29) by a direct method. In fact, as mentioned earlier, before solving (2.28)–(2.29) by marching we used a shooting method to solve the eigenvalue problem. That approach had limited success because the system is stiff, but certainly at moderate values of $k$ there was evidence that multiple neutral modes exist. But for the reasons discussed above the higher modes cannot have $\varGamma \rightarrow 0$ when $k\rightarrow 0$, so it is likely that the higher modes all have $\varGamma \rightarrow \infty$ when $k \rightarrow 0.$

Within the transition layer the leading-order approximation to $(U,V)$ and the streamwise vorticity written in terms of $Y$ are found to be

(5.13)$$\begin{gather} (U,V)\sim(U_0,\varGamma^{1/3}V_0) e^{Q(\varPhi,Y)}+\cdots, \end{gather}$$

(5.14)$$\begin{gather}\omega_x\sim \left( 4\frac{ \varGamma^{5/3}}{Y^5}\sin\left(2\left[\varPhi-\dfrac{\rm \pi}{6}\right]\right)U_0 -\left(1+\dfrac{16\varGamma^2\sin^2\left[\varPhi-\dfrac{\rm \pi}{6}\right]}{Y^6}\right)V_0\right) {\rm e}^{Q(\varPhi,Y)}+\cdots, \end{gather}$$

where

(5.15)\begin{equation} Q(\varPhi, Y)= { \frac {2{\rm i}\varGamma {k }}{Y^2}\sin \left[\varPhi-\frac{ {\rm \pi}}{6}\right]{-c\varGamma^{-{1/6}}} (Y-\bar{Y})^2 }. \end{equation}

Figure 13 shows contours of $U,V$ in the $\varPhi -Y$ plane and isosurfaces of the streamwise vorticity as calculated from the asymptotic forms (5.13)–(5.14) for the point on the neutral curve with $\varGamma =100.$ The results compare well with the results shown in figure 7 as calculated by marching from the full equations. Equations (5.13)–(5.15) explain the origin of the increasingly rapid variations of the velocity and vorticity fields with $\varPhi$ at large $\varGamma$ as $Y$ decreases. The variation is due to the first term on the right-hand side of (5.15) which, since $Y$ is of order $\varGamma ^{1/3}$, is of size $\varGamma ^{1/2}$ and increasing in size as $Y$ decreases. This explains why the variations in $\varPhi$ increase in the lower part of the layer where the disturbance is trapped.

Figure 13. Results for $\varGamma =100$ as predicted by the large $\varGamma$ solution. (a–d) The real and imaginary parts of $U,V$ in the first and second rows, respectively. (e) The $\pm 60\,\%$ isosurfaces of the streamwise vorticity for $\varGamma =100$.

In order to demonstrate this effect even more clearly we have in figure 14(a) shown the isosurfaces of the streamwise vorticity as calculated from the full equations on the neutral curve with $\varGamma =1000$. Figure 14(b) shows the corresponding vorticity isosurfaces as calculated from the asymptotic solution. We see excellent agreement between the two results. In figure 14(c) we have plotted the isosurfaces calculated from the asymptotic solution at $\varGamma =10\,000$. We observe the increasingly large number of foliations due to the first term on the right-hand side of (5.15).

Figure 14. Isosurfaces of 60 % of the minimum and maximum vorticity for $\varGamma = 1000.$ (a) The results from the full eigenvalue calculation and (b) the corresponding high $\varGamma$ prediction. (c) The asymptotic prediction of the same isosurfaces for $\varGamma =10\,000.$