3.1. The modified Orr–Sommerfeld equation

The linear stability of this system is governed by a modified version of the Orr–Sommerfeld equation. To derive this equation, we perturb the basic flow with the expansions

(3.1a–c)\begin{equation} u=U(\kern0.7pt y)+\varepsilon\hat{u}(\kern0.7pt y){\rm e}^{\mathrm{i}\xi},\quad v=\frac{\eta}{\textit{Re}}+\varepsilon\hat{v}(\kern0.7pt y){\rm e}^{\mathrm{i}\xi},\quad p=\varepsilon\hat{p}(\kern0.7pt y){\rm e}^{\mathrm{i}\xi}, \end{equation}

where $\xi =\alpha (x-ct)$, $\alpha \in \mathbb {R}$ is the wavenumber, $c\in \mathbb {C}$ is the corresponding wavespeed and $\varepsilon$ is a small parameter. To obtain the governing equation, (3.1a–c) is substituted into (2.3), terms of $O(\varepsilon ^2)$ are neglected and $\hat {u},\hat {p}$ are both eliminated giving

(3.2a)\begin{equation} \left(U(\kern0.7pt y)-c+\frac{\eta}{\mathrm{i}\alpha\textit{Re}}\frac{\mathrm{d}}{\mathrm{d}y}\right) \left(\frac{\mathrm{d}^2}{\mathrm{d}y^2}-\alpha^2\right)\hat{v} -\frac{\mathrm{d}^2U}{\mathrm{d}y^2}\hat{v} =\frac{1}{\mathrm{i}\alpha\textit{Re}}\left(\frac{\mathrm{d}^2}{\mathrm{d}y^2}-\alpha^2\right)^2\hat{v}, \end{equation}

along with the boundary conditions

(3.2b)\begin{equation} \hat{v}(\pm1)=\left.\frac{\mathrm{d}\hat{v}}{\mathrm{d}y}\right|_{y=\pm1}=0. \end{equation}

To solve (3.2), a Chebyshev collocation approach pioneered by Orszag (Reference Orszag1971) was used. For this particular base flow, an analytic solution in terms of hypergeometric functions is known due to Baldwin (Reference Baldwin1970), but this is somewhat cumbersome to evaluate and it is preferable for the nonlinear extension in § 4 to have a purely numeric solution. Equation (3.2) has been considered in previous studies, applied to different base flows (see Sheppard Reference Sheppard1972; Fransson & Alfredsson Reference Fransson and Alfredsson2003; Sun et al. Reference Sun, Yalcin and Oberlack2024). The results obtained in this paper have been compared with these studies (using the appropriate base flow) to validate the numerical procedure.

Figure 2(a) shows linear neutral curves found for different values of the throughflow number. Recall that when $\eta =0$, we recover Couette flow which is linearly stable for all Reynolds numbers. The relationship between $\alpha$ and $c$ on this neutral curve is shown in figure 2(b) and will be explored in more detail below. The smallest Reynolds number at which linear instability occurs, the linear critical Reynolds number, is denoted $\textit {Re}^{\ell }_c$ and the corresponding wavenumber is denoted $\alpha ^{\ell }_c$. Both are shown in figure 3 as $\eta$ increases. From this, we can see that there exists a critical throughflow number, $\eta ^{\ell }_c$, such that for $\eta >\eta ^{\ell }_c$, linear neutral stability curves exist. From figure 3, we find $\eta ^{\ell }_c\approx 3.35$. This estimate will be refined shortly. The flow is most unstable ($\textit {Re}^{\ell }_c$ is minimised) when $\eta \approx 4.90$. Both of these throughflow numbers agree with Sun et al. (Reference Sun, Yalcin and Oberlack2024) after appropriate rescaling.

Figure 2.

Linear neutral curves for different values of the throughflow number $\eta$. (a) Relationship between $\textit {Re}$ and $\alpha$, and (b) relationship between $c$ and $\alpha$ on the linear neutral curve.

Figure 3.

(a) Linear critical Reynolds number and (b) corresponding wavenumber as $\eta$ is varied. The dashed line is the asymptote as $\eta \to \eta ^{\ell }_c$ from above.

To improve our estimate of $\eta ^{\ell }_c$ and better understand the structure of this solution, we proceed as follows. Consider the limit of (3.2) as $\alpha \to 0$ and $\textit {Re}\to \infty$ such that $\gamma =\alpha \textit {Re}$ remains $\mathit {O}(1)$. Note this is a formal limit and $\gamma$ may still be numerically large or small provided it is finite. This simplifies (3.2a) to

(3.3)\begin{equation} \frac{1}{\mathrm{i}\gamma} \left(\frac{\mathrm{d}^4\hat{v}}{\mathrm{d}y^4} -\eta\frac{\mathrm{d}^3\hat{v}}{\mathrm{d}y^3}\right) -(U-c)\frac{\mathrm{d}^2\hat{v}}{\mathrm{d}y^2} +\frac{\mathrm{d}^2U}{\mathrm{d}y^2}\hat{v}=0. \end{equation}

We can solve (3.3) with the boundary conditions (3.2b), again using Chebyshev collocation, to obtain the linear neutral curve in $(\eta,\gamma )$ space. The resulting neutral curve is shown in figure 4.

Figure 4.

Linear neutral curve in the limit $\alpha \to 0$, $\textit {Re}\to \infty$, $\gamma =\alpha \textit {Re}=\mathit {O}(1)$.

We can use this system to efficiently find $\eta ^{\ell }_c$ to a high accuracy. By the above definition, $\eta ^{\ell }_c$ is the minimum $\eta$ on the curve in figure 4, which is found to be $\eta ^{\ell }_c\approx 3.3511$. The fact that our chosen limit still produces a neutral curve demonstrates that the curves in figure 2 behave such that $\alpha \propto \textit {Re}^{-1}$ as $\textit {Re}\to \infty$. As a result, we do not expect a multi-deck upper and lower branch structure to form that we can analyse analytically in this limit. However, figure 4 also shows that as $\eta$ increases, $\gamma \to \infty$. This implies that in the $\eta \gg 1$ limit, new upper and lower branch structures will emerge, and we study these below.

3.2. Lower and upper branch structures

To analyse the lower and upper branch structures found when $1\ll \eta \ll \mathit {O}(\textit {Re})$, a similar approach to Smith (Reference Smith1979) and Bodonyi & Smith (Reference Bodonyi and Smith1981) is taken for the lower and upper branch, respectively. The coefficients found in the expansions of $\alpha$ and $c$ are similar to those of the asymptotic suction profile in Hughes & Reid (Reference Hughes and Reid1965, p. 13) and Dempsey & Walton (Reference Dempsey and Walton2017, p. 9) ensuring confidence in the results. Due to the similarity to previous studies, only an outline of the derivation will be given, with the aim of illustrating why these structures only emerge in the $\eta \gg 1$ limit.

To perform this analysis, the small parameter $\varepsilon$ is introduced which we will find in terms of the Reynolds number for both structures separately. For the lower branch, we expect a triple deck structure to form as shown in figure 5(a). For the upper branch, six layers are anticipated, as shown in figure 5(b). Note that the coordinate system of figure 5 is flipped when compared with figure 1, due to the transformation (2.6). Since the analysis is completed entirely using the inner coordinates $(\tilde {x},\tilde {y})$, this flipped perspective avoids unnecessary negative signs and results in a closer correspondence to the asymptotic suction profile's upper/lower branch structures. In terms of $\textit {Re}$, the small parameter $\varepsilon$ is not the same for both structures, so a label has been added for clarity in figure 5. When $\varepsilon$ is used without a label, the statement will apply to both structures. Figure 3 shows linear growth of both $\textit {Re}$ and $\alpha$ in $\eta$, so we expect to see this in our distinguished limit. The wavelike variable $\xi$ should also remain $\mathit {O}(1)$, which is achieved by rescaling $x$ and $t$. Hence, we define

(3.4a–d)\begin{equation} \widetilde{\textit{Re}}=\frac{\textit{Re}}{\eta},\quad \tilde{\alpha}=\frac{\alpha}{\eta},\quad \tilde{x}=\eta x,\quad\tilde{t}=\eta t, \end{equation}

to be our variables in this distinguished limit. Recall that when both $\eta,\textit {Re}\gg 1$, we assume $\eta \ll \mathit {O}(\textit {Re})$, in which case we still have $\widetilde {\textit {Re}}\gg 1$. The wavenumber and wavespeed are expanded in the small parameter $\varepsilon$ as

(3.5a)\begin{gather} \tilde{\alpha}=\varepsilon\alpha_1+\varepsilon^2\alpha_2 +\varepsilon^3\alpha_3+\cdots, \end{gather}

(3.5b)\begin{gather}c=c_0-\varepsilon c_1-\varepsilon^2c_2 -\varepsilon^3c_3+\cdots. \end{gather}

We only expect this structure to exist in the $\eta \gg 1$ limit, where the base flow splits into an inner and outer region as stated in (2.5). Substituting the scalings (3.4a–d) for the outer region into (2.3), we find that the perturbation in this outer region is exponentially small. Hence, we restrict our analysis to take place entirely within the inner region with $\mathit {O}(1)$ variable $\tilde {y}$ given in (2.6). In doing so, our solution must decay to zero exponentially as $\tilde {y}\to \infty$ to match the outer region. The expansions

(3.6a–c)\begin{equation} u=U(\kern0.7pt \tilde{y})+\Delta\hat{u}(\kern0.7pt \tilde{y}){\rm e}^{\mathrm{i}\xi},\quad v=\frac{\eta}{\textit{Re}}-\Delta\hat{v}(\kern0.7pt \tilde{y}){\rm e}^{\mathrm{i}\xi},\quad p=\Delta\hat{p}(\kern0.7pt \tilde{y}){\rm e}^{\mathrm{i}\xi}, \end{equation}

are used where $\xi =\tilde {\alpha }(\tilde {x}-c\tilde {t})=\alpha (x-ct)$ and $\varDelta \ll 1$. Substituting (3.6a–c) into (2.3) and neglecting terms of $\mathit {O}(\varDelta ^2)$, we derive the governing equations of a linear perturbation in the inner region as

(3.7a)\begin{gather} \mathrm{i}\tilde{\alpha}\hat{u}+\frac{\mathrm{d}\hat{v}}{\mathrm{d}\tilde{y}}=0, \end{gather}

(3.7b)\begin{gather}\mathrm{i}\tilde{\alpha}(U-c)\hat{u}+\frac{\mathrm{d}U}{\mathrm{d}\tilde{y}}\hat{v} ={-}\mathrm{i}\tilde{\alpha}\hat{p} +\frac{1}{\widetilde{\textit{Re}}}\left( \frac{\mathrm{d}^2\hat{u}}{\mathrm{d}\tilde{y}^2} +\frac{\mathrm{d}\hat{u}}{\mathrm{d}\tilde{y}} -\tilde{\alpha}^2\hat{u}\right), \end{gather}

(3.7c)\begin{gather}\mathrm{i}\tilde{\alpha}(U-c)\hat{v} ={-}\frac{\mathrm{d}\hat{p}}{\mathrm{d}\tilde{y}} +\frac{1}{\widetilde{\textit{Re}}}\left( \frac{\mathrm{d}^2\hat{v}}{\mathrm{d}\tilde{y}^2} +\frac{\mathrm{d}\hat{v}}{\mathrm{d}\tilde{y}} -\tilde{\alpha}^2\hat{v}\right). \end{gather}

We aim to find $c_0,c_1,\alpha _1$ and $\varepsilon$, which is achieved by finding leading order solutions across all of the layers. These solutions are subject to no slip and no penetration at $\tilde {y}=0$, exponential decay as $\tilde {y}\to \infty$ and matching between adjacent layers.

Figure 5.

Layer structures for the (a) lower branch and (b) upper branch solutions. The scalings are given in terms of the inner variables $\tilde {x}$ and $\tilde {y}$, resulting in a flipped perspective compared with figure 1.

The inviscid cores are perhaps most important, as the leading order vertical velocity, $\hat {v}(\kern0.7pt \tilde {y})=\varepsilon G_0(\kern0.7pt \tilde {y})+\mathit {O}(\varepsilon ^2)$, can be found to be given by

(3.8)\begin{equation} G_0(\kern0.7pt \tilde{y})=\mathrm{i}\alpha_1A_0|c_0-{\rm e}^{-\tilde{y}}|, \end{equation}

where $A_0$ is a constant of integration. This is significant as both structures in figure 5 require a thinner layer below the inviscid core. For this to be achieved, it can be verified that we require $G_0(\kern0.7pt \tilde{y})\to 0$ as $\tilde{y}\to 0$, and hence that $c_0=1$. It transpires that a solution with unit wavespeed is only possible in the $\eta \gg 1$ limit, explaining why these upper and lower branch structures only exist in this limit. It can be seen in figure 2(b) that the largest value of $c$ is approaching unity as $\eta$ increases but does not attain this value. With $\eta =10$, $c$ only obtains a maximum of $0.956$, resulting in the upper/lower branch structures breaking down quickly. A large value of $\eta$ is needed to maintain the structures as $\widetilde {\textit {Re}}$ increases.

The scalings given in figure 5 are found by matching between layers as appropriate. Both $\varepsilon _l$ and $\varepsilon _u$ can be found in terms of $\widetilde {\textit {Re}}$ by considering the balance between viscous and inertial terms in appropriate layers. Doing so gives

(3.9a,b)\begin{equation} \varepsilon_l=\widetilde{\textit{Re}}^{{-}1/4},\quad \varepsilon_u=\widetilde{\textit{Re}}^{{-}1/6}. \end{equation}

Completing the analysis and matching between each layer allows us to obtain $\alpha _1$ and $c_1$. It should be noted that for the upper branch solution, the existence of the critical layer results in a phase jump of $+{\rm \pi} {\rm i}$ from above the critical layer to below it. Linear critical layers are well studied and the $+$ sign comes from the fact $U'(\kern0.7pt y_c)<0$, where $y_c$ is the location of the critical layer. For a more detailed explanation of this jump condition, see Ruban et al. (Reference Ruban, Gajjar and Walton2023).

Overall, it is found that in the lower branch structure,

(3.10a)\begin{equation} \alpha^l_1\approx1.0006525,\quad c^l_1\approx2.2967186, \end{equation}

and for the upper branch structure,

(3.10b)\begin{equation} \alpha^u_1=c^u_1=(2{\rm \pi}^2)^{{-}1/6}. \end{equation}

Combining these constants with the appropriate $\varepsilon$, we can approximate the neutral curves in figure 2 with

(3.11a)\begin{gather} \tilde{\alpha}=\alpha^l_1\widetilde{\textit{Re}}^{{-}1/4} \quad\iff\quad \alpha=\alpha^l_1\eta^{5/4}\textit{Re}^{{-}1/4} \end{gather}

(3.11b)\begin{gather} \tilde{\alpha}=(2{\rm \pi}^2)^{{-}1/6}\widetilde{\textit{Re}}^{{-}1/6} \quad\iff\quad \alpha=(2{\rm \pi}^2)^{{-}1/6}\eta^{7/6}\textit{Re}^{{-}1/6}. \end{gather}

This is a leading order description based on $\widetilde {\textit {Re}}$ being asymptotically large, meaning we only expect this approximation to be accurate for very large values of $\widetilde {\textit {Re}}$.

3.3. Comparison of asymptotic theory and numerical solutions

Now that we have obtained leading order upper and lower branch asymptotic approximations to the curves in figure 2(a), we wish to analyse their accuracy. To do so, we rescale the solutions in figure 2(a) to the $(\widetilde {\textit {Re}},\tilde {\alpha })$ coordinates. We also include a neutral curve with $\eta =100$, as for any finite $\eta$, our upper and lower branch structures will eventually break down and the curve will tend to the appropriate limit from figure 4. Increasing $\eta$ results in the multi-deck structure holding for larger values of $\widetilde {\textit {Re}}$, enabling (3.11) to be a better approximation. The results of this rescaling, as well as the asymptotes from (3.11), are shown in figure 6.

Figure 6.

Linear neutral curves, rescaled into the parameters $\widetilde {\textit {Re}}=\textit {Re}/\eta$ and $\tilde {\alpha }=\alpha /\eta$. The black curves represent the upper and lower branch asymptotic approximations.

As $\widetilde {\textit {Re}}$ increases, the curves diverge greatly. For $\eta \leqslant 10$, the asymptotic structures quickly disappear and the behaviour shown in figure 4 dominates. The value $\eta =100$ is sufficiently large that our upper and lower branch structures do not break down within the domain plotted, allowing us to see a convergence to the appropriate asymptotes.

As an aside, we see that as $\eta$ increases, the curves coincide near $\widetilde {\textit {Re}}^{\ell }_c=\textit {Re}^{\ell }_c/\eta$. This implies that in the large $\eta$ limit, the curves in figure 6 are accurate approximations of the linear neutral curve for asymptotic suction flow near the linear critical Reynolds number and we can use this to estimate the critical Reynolds number for this flow. Using the $\eta =100$ curve in figure 6, it is estimated that the linear critical Reynolds number of the asymptotic suction profile is 54 378.60263 which is in very close agreement with the value of 54 378.62032 from Yalcin, Turkac & Oberlack (Reference Yalcin, Turkac and Oberlack2021), lending further credence to the reliability of these results.

4.1. Nonlinear numerical method

Given our previous linear analysis, we now aim to perform a nonlinear analysis of our stability problem. This will be performed at finite Reynolds number and will involve looking for travelling wave solutions to (2.3). The numerical method used in this section is based on that employed by Wong & Walton (Reference Wong and Walton2012). We seek solutions of the form

(4.1) \begin{equation} (u,v,p)=\left(U(\kern0.7pt y)+\hat{u}_0(\kern0.7pt y),\frac{\eta}{\textit{Re}},\hat{p}_0(\kern0.7pt y)\right) +\sum_{n=1}^\infty(\hat{u}_n(\kern0.7pt y),\hat{v}_n(\kern0.7pt y),\hat{p}_n(\kern0.7pt y)) {\rm e}^{\mathrm{i}n\hat{\xi}} +\text{c.c.}, \end{equation}

where $U(\kern0.7pt y)$ is given by (2.2a,b) and $\hat {\xi }=\alpha (x-ct)$. The terms $\hat {u}_0(\kern0.7pt y)$ and $\hat {p}_0(\kern0.7pt y)$ are the mean flow distortion and mean pressure distortion, respectively. For simplicity of notation, we define $\hat {v}_{-n}(\kern0.7pt y)=\overline {\hat {v}_n(\kern0.7pt y)}$. Substituting (4.1) into (2.3) gives our governing equations. To ensure that the problem is computationally feasible, only the first $N$ modes are retained with the value of $N$ taken sufficiently large to achieve satisfactory resolution of the solutions. The equation involving the mean pressure distortion uncouples from the rest of the system and need not be considered further. Taking the resulting system of equations, every instance of $\hat {u}_n$ and $\hat {p}_n$ for $n\geqslant 1$ can be eliminated to give the following system:

(4.2a)\begin{gather} \hat{u}_0''-\eta\hat{u}_0' =\sum_{k=1}^N\frac{\mathrm{i}\textit{Re}}{\alpha k} (\hat{v}_k''\hat{v}_{{-}k} -\hat{v}_k\hat{v}_{{-}k}''),\end{gather}

(4.2b)\begin{gather} {-}\frac{1}{\alpha^2n^2\textit{Re}} \mathrm{D}_n^2\hat{v}_n +\frac{\mathrm{i}}{\alpha n} \left(U+\hat{u}_0-c +\frac{\eta}{\mathrm{i}\alpha n\textit{Re}}\frac{\mathrm{d}}{\mathrm{d}y}\right) \mathrm{D}_n\hat{v}_n -\frac{\mathrm{i}}{\alpha n} (U''+\hat{u}_0'')\hat{v}_n\nonumber\\ =\sum_{k=1}^{n-1}S_{k,n}^{(1)}(\kern0.7pt y) +\sum_{k=n+1}^NS_{k,n}^{(2)}(\kern0.7pt y), \end{gather}

(4.2c)\begin{gather} \hat{u}_0(\pm1)=\hat{v}_n(\pm1)=\hat{v}_n'(\pm1)=0, \end{gather}

where

(4.2d)\begin{gather} S_{k,n}^{(1)}(\kern0.7pt y) =\left(1-\frac{k}{n-k}\right)\hat{v}_k\hat{v}_{n-k}' +\frac{1}{\alpha^2n(n-k)} (\hat{v}_k''\hat{v}_{n-k}' +\hat{v}_k'\hat{v}_{n-k}'')\nonumber\\ -\frac{1}{\alpha^2kn} (\hat{v}_k'''\hat{v}_{n-k} +\hat{v}_k''\hat{v}_{n-k}'), \end{gather}

(4.2e)\begin{gather} S_{k,n}^{(2)}(\kern0.7pt y) =\left(1-\frac{k}{n-k}\right)\hat{v}_k\hat{v}_{n-k}' +\left(1-\frac{n-k}{k}\right)\hat{v}_k'\hat{v}_{n-k} +\frac{1}{\alpha^2k(n-k)} (\hat{v}_k''\hat{v}_{n-k}' +\hat{v}_k'\hat{v}_{n-k}'')\nonumber\\ -\frac{1}{\alpha^2nk} (\hat{v}_k'''\hat{v}_{n-k} +\hat{v}_k''\hat{v}_{n-k}') -\frac{1}{\alpha^2n(n-k)} (\hat{v}_k\hat{v}_{n-k}''' +\hat{v}_k'\hat{v}_{n-k}''), \end{gather}

(4.2f)\begin{equation} \mathrm{D}_n\equiv\frac{\mathrm{d}^2}{\mathrm{d}y^2}-\alpha^2n^2, \end{equation}

and $'$ denotes differentiation with respect to $y$. It can be shown, given a solution $\hat {v}_n(\kern0.7pt y)$ to (4.2a–c), that $\hat {v}_n\mapsto \hat {v}_n{\rm e}^{\mathrm {i}n\theta }$ is also a valid solution to (4.2) for any $\theta \in \mathbb {R}$; this means that we need to impose a condition to fix the value of $\theta$, which is achieved with the equation

(4.2g)\begin{equation} \sum_{n=1}^N {\rm Re}(\hat{v}_n(0))=0. \end{equation}

It can be shown that $\exists \theta \in [0,2{\rm \pi} ]$ such that (4.2g) can always be satisfied. For given values of $\textit {Re}$ and $\eta$, we expect our numerical approach to yield a curve of neutral solutions in the $(\alpha,c)$ plane. To identify different solutions along this curve, an additional input to our system is required which can be perturbed to find different solutions at fixed values of $\textit {Re}$ and $\eta$. This new input is given in the form of a nonlinear amplitude

(4.2h)\begin{equation} \hat{A}^2=\sum_{n=1}^N {\rm Re}(\hat{v}_n(0))^2+{\rm Im}(\hat{v}_n(0))^2. \end{equation}

It should be emphasised that the only approximation made in the derivation of (4.2) is the truncation of the Fourier series in $\hat {\xi }$. This means a fully resolved solution to (4.2) with $\hat {A}>0$ corresponds to a nonlinear travelling wave solution to the full Navier–Stokes equations (2.3).

To solve (4.2) numerically, we again use Chebyshev collocation. The solution will be evaluated using $M+1$ Chebyshev polynomials evaluated at the Gauss points $y_j=\cos ({\rm \pi} (j-1)/M)$ for $j=1,2,\ldots,M+1$. We express $\hat {u}_0(\kern0.7pt y)$ and $\hat {v}_n(\kern0.7pt y)$ as

(4.3a,b)\begin{equation} \hat{u}_0(\kern0.7pt y)=\sum_{m=1}^{M+1}\gamma_mT_{m-1}(\kern0.7pt y),\quad \hat{v}_n(\kern0.7pt y)=\sum_{m=1}^{M+1}(a_m^{(n)}+\mathrm{i}b_m^{(n)})T_{m-1}(\kern0.7pt y), \end{equation}

where $T_m(\kern0.7pt y)$ is the $m$th Chebyshev polynomial of the first kind and $\gamma _m,a_m^{(n)},b_m^{(n)}$ are real coefficients of the Chebyshev series. This gives a total of $(2N+1)(M+1)+2$ parameters for our system and (4.2) corresponds to $(2N+1)(M+1)+2$ equations after substitution of (4.3a,b) and splitting into real/imaginary components. This allows us to express (4.2) as a vector equation of the inputs $\alpha,c,a_m^{(n)},b_m^{(n)},\gamma _m$ and Newton iteration can be used to find a solution.

While this is a well-posed numerical problem, in practice, we require a very close initial guess to a correct solution for the procedure to converge. This means that we will not be able to find any solutions without a systematic way to find a starting guess. Observe that if $n=1$ and $\hat {u}_0\equiv 0$, the left-hand side of (4.2b) is equivalent to (3.2a) after moving all the terms in (3.2a) to the left-hand side and multiplying by $\mathrm {i}/\alpha$. Suppose $\hat {v}(\kern0.7pt y),\alpha _{or},c_{or}$ is a solution to (3.2) for a given $\textit {Re},\eta$ normalised so that $\hat {v}(0)=1$. Then, by this observation,

(4.4)\begin{equation} \hat{v}_1(\kern0.7pt y)=\mathrm{i}\hat{A}\hat{v}(\kern0.7pt y),\quad \alpha=\alpha_{or},\quad c=c_{or},\quad \hat{u}_0(\kern0.7pt y)=\hat{v}_n(\kern0.7pt y)\equiv0\quad\text{for }n\geqslant2, \end{equation}

is a solution to (4.2b), (4.2c), (4.2g) and (4.2h). Importantly, (4.4) is a solution for arbitrary $\hat {A}\geqslant 0$. While (4.4) does not satisfy (4.2a) unless $\hat {A}=0$, setting $\hat {A}$ close to zero makes (4.4) arbitrarily close to a non-trivial solution. This observation means we can take a neutral solution to the modified Orr–Sommerfeld equation and bifurcate from this curve into non-zero amplitude space.

This will allow us to find solutions to (4.2), but due to the complexity of the system, convergence would require taking very small steps in $\hat {A}$. To remedy this and make the algorithm more robust, the Jacobian of (4.2) has been calculated analytically. This requires taking derivatives of each equation in (4.2) with respect to each parameter $\alpha,c,a_m^{(n)},b_m^{(n)},\gamma _m$. These derivatives, while important, are algebraically cumbersome and so have not been included for brevity.

This method requires a linear neutral curve from which to bifurcate to generate the initial solutions. This means that this method can only be used for $\eta \geqslant \eta ^{\ell }_c$. For values outside of this range, or to investigate neutral solutions with Reynolds numbers smaller than $\textit {Re}^{\ell }_c$, a different approach is used. Results can be obtained by stepping through values of $\eta$ or $\textit {Re}$ in $\hat {A}>0$ space. This allows an investigation into a solution space that was not previously accessible. Indeed, it is theoretically possible to decrease $\eta$ to zero, allowing a nonlinear investigation of Couette flow. Although this procedure was used to obtain some of the results, solutions in these regimes require far more Fourier modes to achieve accurate convergence. Using more than $N=40$ modes resulted in unreasonable computation times and so a detailed investigation of these highly nonlinear solutions has not been included.

4.2. Weakly nonlinear analysis

To validate the results described subsequently in § 4.3, a weakly nonlinear analysis is used for comparison. The result of this analysis, the Landau–Stuart equation, was first conceived by Landau (Reference Landau1944) and subsequently formally derived by Stuart (Reference Stuart1960) and Watson (Reference Watson1960). Weakly nonlinear theory has been successfully used to classify the bifurcation from the linear neutral curve in a variety of flows including circular Couette flow (Davey Reference Davey1962) and Poiseuille flow (Reynolds & Potter Reference Reynolds and Potter1967). For a full account of the theory, the reader is referred to Ruban et al. (Reference Ruban, Gajjar and Walton2023). Our aim is to see how the nonlinear amplitude $\hat {A}$ varies with $\alpha$ at fixed values of $\textit {Re}$ and $\eta$. This approximation should correctly identify whether there is a subcritical or supercritical bifurcation. Suppose $\alpha _0,c_0$ is a neutral solution to (3.2) at some $\textit {Re}$ and $\eta$. We wish to analyse the system at $\alpha =\alpha _0+\delta a$, where $a=\mathit {O}(1)$ and $\delta \ll 1$.

To undertake this analysis, we expand in the small parameter $\varepsilon$, which will be chosen as a function of $\delta$, using the streamfunction

(4.5)\begin{equation} \psi(x,y,t) =\varPsi(x,y) +\varepsilon\psi_1(\bar{\xi},y,\bar{t}) +\varepsilon^2\psi_2(\bar{\xi},y,\bar{t})+\cdots, \end{equation}

where

(4.6a,b)\begin{equation} \varPsi(x,y)=\int_0^y U(z)\,\mathrm{d}z-\frac{\eta}{\textit{Re}}x,\quad \psi_1(\bar{\xi},y,\bar{t})=\bar{A}(\bar{t}){\rm e}^{\mathrm{i}\bar{\xi}}\phi_1(\kern0.7pt y)+\text{c.c.}, \end{equation}

and our variables are defined as $\bar {\xi }=\alpha _0(x-c_0t)$ and $\bar {t}=\delta t$. From (4.6a,b), we see there are two time scales, $t$ corresponding to oscillations of the perturbations and $\bar {t}$ to describe the slow growth of amplitude. We want the simplest non-trivial solution for the amplitude and it transpires this will occur when $\delta =\varepsilon ^2$. The analysis involves the first three orders of equations from substituting (4.5) into (2.3). At leading order, we find that $\phi _1(\kern0.7pt y)$ is governed by (3.2), making any solution arbitrary up to a multiplicative constant. To normalise the problem, it was chosen to enforce $\phi _1(0)=1$. From this analysis, we obtain the Landau–Stuart equation,

(4.7)\begin{equation} \frac{\mathrm{d}|\bar{A}|^2}{\mathrm{d}\bar{t}} =2a\kappa_r|\bar{A}|^2-\ell_r|\bar{A}|^4, \end{equation}

where $\kappa _r$ and $\ell _r$ are real constants found as part of the analysis. The expressions of $\kappa _r,\ell _r$ are similar to those of Ruban et al. (Reference Ruban, Gajjar and Walton2023), after propagating the non-trivial effects of $\eta \neq 0$ throughout and varying $\alpha$ instead of $\textit {Re}$. These expressions have been omitted explicitly for brevity.

Our aim is to compare the growth of $\bar {A}$ (the weakly nonlinear amplitude) and $\hat {A}$ (the strongly nonlinear amplitude). We expect a neutral solution in our weak formulation if $|\bar {A}|$ is constant, i.e.

(4.8)\begin{equation} |\bar{A}|=\sqrt{\frac{2a\kappa_r}{\ell_r}}, \end{equation}

from (4.7). To compare the weak and strong formulations, note that the leading orders in (4.1) and (4.5) are analogous,

(4.9)\begin{equation} \hat{v}_1(\kern0.7pt y){\rm e}^{\mathrm{i}\hat{\xi}}+\text{c.c.} \approx{-}\mathrm{i}\varepsilon\alpha_0\bar{A}{\rm e}^{\mathrm{i}\bar{\xi}}\phi_1(\kern0.7pt y)+\text{c.c.} \end{equation}

Since we expect a weakly nonlinear solution to only match the strongly nonlinear solution near the linear neutral curve, we approximate $\hat {\xi }\approx \bar {\xi }$. Recall from (4.4) that close to the linear neutral curve, we have $\hat {v}_1(0)\approx \mathrm {i}\hat {A}$. Finally, as part of the weakly nonlinear analysis, it was chosen that $\phi _1(0)=1$, so evaluating (4.9) at $y=0$ gives us

(4.10)\begin{equation} \varepsilon ={-}\frac{\hat{A}}{\alpha_0\bar{A}} =\frac{\hat{A}}{\alpha_0}\sqrt{\frac{\ell_r}{2a\kappa_r}}, \end{equation}

after using (4.8) to find $\bar {A}$ so that $\varepsilon \in \mathbb {R}$ and $\varepsilon >0$. Using (4.10) and the fact $\varepsilon ^2=\delta =(\alpha -\alpha _0)/a$, we find the relationship

(4.11)\begin{equation} \hat{A}^2\approx\frac{2\kappa_r}{\ell_r}\alpha_0^2(\alpha-\alpha_0). \end{equation}

This is our approximation of the nonlinear amplitude and will allow a partial check of the numerical results in § 4.3. This is a linear approximation for $\hat {A}^2$, so when plotting $\hat {A}^2$ against $\alpha$, the approximation will be, at best, tangent to the curve at $\hat {A}=0$.

4.3. Numerical results

In this section, we examine results generated from the numerical approach discussed in § 4.1. The neutral curves in figure 2 exist at relatively high Reynolds numbers, which generally means that many Fourier modes and Chebyshev polynomials are required to accurately generate nonlinear continuations of these curves. To reduce this computational burden, most of the simulations are performed with Reynolds numbers as small as possible, most commonly of magnitude $10^5$.

To begin, the nonlinear amplitude $\hat {A}$ is compared with the approximation from the weakly nonlinear analysis given in (4.11). The amplitudes have been normalised relative to the maximum value attained by $\hat {A}$, which is denoted $\hat {A}_{{max}}$. That is to say, we have defined $\hat {A}_n=\hat {A}/\hat {A}_{{max}}$. Equation (4.11) is rescaled to approximate $\hat {A}_n$ instead of $\hat {A}$. To aid the comparison, $\hat {A}_n^2$ has been plotted as this makes the weak approximation linear in $\alpha$. In figure 7, we see the neutral curves at a Reynolds number of $\textit {Re}={360\,000}$ and a throughflow number of $\eta =6$. It appears that the linear approximation is tangent to the nonlinear amplitude at $\hat {A}=0$, as expected. Both curves in figure 7 correspond to a subcritical bifurcation. Here, $\hat {A}_n$ is of the magnitude $10^{-3}$, which is very small compared with its maximum value of 1. This limited accuracy of the linear approximation can also be observed in other studies (see for example Deguchi & Walton Reference Deguchi and Walton2018, figure 2). Figure 8 shows how, at a fixed Reynolds number of 500 000, increasing the throughflow strength decreases the slope $\kappa _r/\ell _r$ from (4.11). Due to the linear neutral curves only coming into existence at $\eta ^{\ell }_c\approx 3.35$, the growth has only been tracked from $\eta =4$, but it is clear that increasing throughflow makes $\alpha$ grow more rapidly when varying amplitude.

Figure 7.

Normalised amplitude growth predicted by the weakly nonlinear approximation and the strongly nonlinear simulation when bifurcating from (a) the upper branch and (b) the lower branch of the linear neutral curve. A throughflow number of $\eta =6$ was used and a Reynolds number of $\textit {Re}={360\,000}$.

Figure 8.

Slope of the initial amplitude growth, as predicted by the weakly nonlinear theory, of the upper branch at a Reynolds number of 500 000 as a function of the throughflow number.

A selection of nonlinear neutral curves with $\eta =6$ are shown in figure 9. Between $N=10$ and $N=20$ Fourier modes and $M=150$ Chebyshev polynomials were used to generate the various curves in figure 9. Note that the lower branch is only shown for a limited range of amplitudes beyond which it proved impossible to obtain a converged solution. The inability to computationally continue the lower branch into the strongly nonlinear regime has been observed before, e.g. for boundary-layer flow (Hall Reference Hall1995, figure 2) and plane Poiseuille flow (Deguchi & Walton Reference Deguchi and Walton2018, figure 2).

Figure 9.

Nonlinear neutral curves for $\eta =6$. Both (a,c) a two-dimensional (2-D) representation and (b,d) a three-dimensional (3-D) representation are shown with the black curves being the appropriate linear neutral curve from figure 2. The relationship between the amplitude and wavenumber are shown in panels (a,b) and the relationship between the wavenumber and wavespeed in panels (c,d).

Both the relationship between $\hat {A}$ and $\alpha$ (figure 9a,b) and the relationship between $\alpha$ and $c$ (figure 9c,d) are shown. The lower branch can be observed leaving the linear neutral curves in figures 9(b) and 9(d) for $\textit {Re}={360\,000}, {500\,000}$, but these solutions have only been calculated for small $\hat {A}$. While the choice of amplitude measure is somewhat arbitrary, the relationship between $\alpha$ and $c$ is objective and should be used as the primary view of these neutral curves. All graphs in figure 9 show a rapid growth in wavenumber for small perturbations in $\hat {A}$ or $c$ near the linear neutral curve.

The curves in figure 9 form part of the nonlinear neutral surface that contains neutral solutions far below the linear limit of $\textit {Re}^{\ell }_c\approx {357\,769}$. If $\textit {Re}_c$ denotes the minimum $\textit {Re}$ on the nonlinear neutral surface, when $\eta =6$, it was found that $\textit {Re}_c\approx {19\,200}$. The nonlinear critical Reynolds number is shown for a range of values of $\eta$ in figure 10. There is a small amount of curvature; however, the curve appears to be an approximately linear relationship. Crucially, figure 10 implies that faster throughflow leads to a larger range of stable modes. So while $\textit {Re}_c$ is an order of magnitude smaller than its linear counterpart $\textit {Re}^{\ell }_c$, both show that sufficiently strong throughflow stabilises Couette flow.

Figure 10.

Nonlinear critical Reynolds number as the throughflow number varies.

Individual neutral solutions can also be investigated, some of which are shown in figure 11. The functions shown are the mean flow distortion $u_0(\kern0.7pt y)$ and the travelling wave components $v_n(\kern0.7pt y){\rm e}^{\mathrm {i}n\hat {\xi }}+\text {c.c.}$. Since the $v_n(\kern0.7pt y)$ are complex-valued, we plot twice their imaginary part. The reason for this choice is discussed in § 6 when we compare these functions to their analytic counterparts obtained from an asymptotic theory. From each figure, we observe an ordering where $u_0>v_1>v_2>v_3$ with the disturbance concentrated over a region located at the top of the domain. This suggests the existence of a structure that we will look for in § 5. As the amplitude increases, this ordering and structure breaks down. This is most noticeable in figure 11(d), where $v_2$ and $v_3$ are of a similar order of magnitude and $u_0$ does not exponentially decay to zero as $y\to -1$. This behaviour is due to figure 11(d) using an amplitude which is near the maximum for the corresponding neutral curve in figure 9(a).

Figure 11.

The first few modes of the nonlinear neutral solutions. Each graph uses $\eta =6$ and different values for $\textit {Re}$ and $\hat {A}$, as indicated in the bottom right. For the complex-valued modes, $2{\rm Im}(f(\kern0.7pt y))$ is plotted.

The individual components shown in figure 11 can be appropriately combined to visualise the total perturbation to the base flow, as shown in figure 12. The flow visualised is for the parameters shown in figure 11(c). Over one period, we can see that the majority of the activity happens towards the top of the domain. This is due to the mean flow distortion (the dominant flow term) being significant only for $y>0$.

Figure 12.

A nonlinear neutrally stable perturbation to the base flow with $\eta =6$, $\textit {Re}={500\,000}$, $\hat {A}=0.005$. The individual components are shown in figure 11(c).

These neutral solutions allow a much more thorough and realistic investigation into the behaviour of this flow than their linear counterparts. It is clear by comparing the values of $\textit {Re}^{\ell }_c$ and $\textit {Re}_c$ in figures 3(a) and 10, respectively, that these nonlinear results differ greatly from the linear theory. However, obtaining these solutions is computationally expensive, especially at large Reynolds numbers. Furthermore, trying to track $\textit {Re}_c$ becomes computationally infeasible as $\eta$ decreases, preventing us from numerically finding the critical value of $\eta$ beyond which these nonlinear solutions exist. The results in figure 11 seem to indicate an asymptotic structure which would give a pathway to analytically investigate this nonlinear problem.

5.1. Inviscid core

In the inviscid core, we expand our fluid parameters as

(5.2a)\begin{gather} u(x,y,t) =U(\kern0.7pt y)+\varepsilon U_M(\kern0.7pt y)+\varepsilon^2AF(\kern0.7pt y)\cos\xi+\cdots, \end{gather}

(5.2b)\begin{gather}v(x,y,t) =\frac{\eta}{\textit{Re}}+\varepsilon^2AG(\kern0.7pt y)\sin\xi+\cdots, \end{gather}

(5.2c)\begin{gather}p(x,y,t) =\varepsilon^2AP(\kern0.7pt y)\cos\xi+\cdots, \end{gather}

where $\xi =\alpha (x-ct)$. Substituting (5.2) into (2.3) gives us the governing equations

(5.3a)\begin{gather} \alpha F(\kern0.7pt y) =G'(\kern0.7pt y), \end{gather}

(5.3b)\begin{gather}U''(\kern0.7pt y)G(\kern0.7pt y) =(U(\kern0.7pt y)-c)(G''(\kern0.7pt y)-\alpha^2G(\kern0.7pt y)), \end{gather}

(5.3c)\begin{gather}P(\kern0.7pt y) =\alpha^{{-}1}(U'(\kern0.7pt y)G(\kern0.7pt y)-(U(\kern0.7pt y)-c)G'(\kern0.7pt y)), \end{gather}

(5.3d)\begin{gather}U_M''(\kern0.7pt y) =\eta U_M'(\kern0.7pt y), \end{gather}

where $'$ denotes a derivative with respect to $y$. These equations are subject to the inviscid conditions $G(\pm 1)=0$, as well as jump conditions across the critical layer. As discussed in § 2, the $\eta <0$ case can be transformed into the $\eta >0$ case by appropriate choice of coordinate system and frame of reference, so we can restrict ourselves to $\eta >0$ without loss of generality. This means that from (5.3d), we have a piecewise exponential mean flow distortion. We expect continuity of $U_M(\kern0.7pt y)$ over the critical layer; however, the derivative of $U_M(\kern0.7pt y)$ is necessarily discontinuous as we shall see below. Solving (5.3d) subject to these boundary conditions gives

(5.4) \begin{align} U_M(\kern0.7pt y)=\begin{cases} A_1(1-\exp({\eta(\kern0.7pt y+1)})) & \text{if }-1\leqslant y< y_c,\\ A_1\left(\dfrac{1-\exp({\eta(\kern0.7pt y_c+1)})}{1-\exp({\eta(\kern0.7pt y_c-1)})}\right)(1-\exp({\eta(\kern0.7pt y-1)})) & \text{if }y_c< y\leqslant1. \end{cases} \end{align}

It is clear that $U_M(\kern0.7pt y)$ cannot have a continuous derivative and be non-trivial, and hence the constant $A_1$ is determined by the jump condition of $U_M'(\kern0.7pt y)$ over the critical layer. The size of this jump is given by

(5.5)\begin{align} J_{CL} =U_M'(\kern0.7pt y_c^+)-U_M'(\kern0.7pt y_c^-) =\eta A_1{\rm e}^{\eta y_c} \bigg[ {\rm e}^\eta -\bigg(\frac{1-{\rm e}^{\eta(\kern0.7pt y_c+1)}}{1-{\rm e}^{\eta(\kern0.7pt y_c-1)}}\bigg) {\rm e}^{-\eta} \bigg]. \end{align}

Next, the Rayleigh equation (5.3b) is solved to find $G(\kern0.7pt y)$. We also have the boundary conditions $G(\pm 1)=0$ and a jump condition over the critical layer. It is expected that when solving for $y>y_c$, there will be a term of the form $(\kern0.7pt y-y_c)\ln (\kern0.7pt y-y_c)$ in the Frobenius solution for $G(\kern0.7pt y)$ that leads to a logarithmic singularity in $F(\kern0.7pt y)$. Since we are assuming a strongly nonlinear critical layer, the appropriate jump condition to continue the solution to $y< y_c$ will be to replace this logarithmic term with

(5.6)\begin{equation} (\kern0.7pt y-y_c)\ln|y-y_c|, \end{equation}

as done in other studies (see for example Bodonyi et al. Reference Bodonyi, Smith and Gajjar1983; Walton Reference Walton2003; Walton & Barnes Reference Walton and Barnes2023). A full explanation of this nonlinear phase shift condition can be found from Ruban et al. (Reference Ruban, Gajjar and Walton2023). While it is common to solve the Rayleigh equation numerically, the specific form of the base flow (2.2a,b) allows us to take a semi-analytic approach. This approach is described in Appendix A and culminates in (A4) defining the relationship between $\alpha$ and $c$. The normalising condition used was $G(\kern0.7pt y_c)=1$.

The relationship between $\alpha$ and $c$ that results from (A4) is shown in figure 14. Perhaps the most important property to discuss is the existence of solutions. None of the curves in figure 14 span the full domain $0\leqslant c\leqslant 1$. Furthermore, it is known that solutions do not exist at $\eta =0$ (pure plane Couette flow), so it is important to establish the parameter space for which solutions can be found. We will now show that solutions exist given sufficiently strong throughflow, so we define $\eta _c$ such that solutions in figure 14 exist for all $\eta \geqslant \eta _c$.

Figure 14.

Relation found between $\alpha$ and $c$ from the solvability condition (A4) for the Rayleigh equation (5.3b).

To identify $\eta _c$, we consider the limit $\alpha \ll 1$. From figure 14, it seems that to leading order, the wavespeed remains $\mathit {O}(1)$ and so we set $c=c_0$ to be determined. Substituting this limit into (5.3b) leads to the simplified equation

(5.7)\begin{equation} G(\kern0.7pt y)U''(\kern0.7pt y) =(U(\kern0.7pt y)-c_0)G''(\kern0.7pt y), \end{equation}

which can be solved analytically. Solving (5.7) subject to $G(\pm 1)=0$ and the jump condition (5.6) for a non-zero solution results in the solvability condition

(5.8a,b)\begin{equation} c_0(1-c_0) \ln\frac{\kappa_0(1-c_0)}{(1+\kappa_0)c_0} +\kappa_0+c_0 =0,\quad \kappa_0=\frac{{\rm e}^{-\eta}}{2\sinh\eta}. \end{equation}

For a given $\eta$, (5.8a,b) can be solved for $c_0$ and two branches of solutions are shown in figure 15. These two branches correspond to the upper and lower values of $c$, between which we found solutions in figure 14. The behaviour of these branches for large $\eta$ can be derived as

(5.9a,b)\begin{equation} c_0 =1+\frac{1}{W_{{-}1}(-{\rm e}^{{-}2\eta})},\quad c_0 =\zeta {\rm e}^{{-}2\eta}, \end{equation}

on the upper and lower branches, respectively, where $\zeta =1/W_0(1/e)\approx 3.5911$ is the solution to the equation $\zeta (\ln \zeta -1)=1$ and $W_k(x)$ is the $k$th branch of the Lambert $W$-function. These approximations are also shown in figure 15 and visually agree with the true solution for $\eta \gtrsim 4$. The rates at which the upper and lower branches tend to 1 and 0, respectively, explain why as $\eta$ increases, we see $c\to 0$ more rapidly than $c\to 1$ in figure 14.

Figure 15.

Leading order wavespeed in the long-wavelength limit for varying throughflow strength. The dotted lines are the asymptotic approximations given in (5.9a,b).

It can be shown that when varying $c_0$, the minimum of the left-hand side of (5.8a,b) always occurs at $c_0=0.5$. It follows that solutions to (A4) only exist when the left-hand side of (5.8a,b) with $c_0=0.5$ is $\leqslant 0$. Setting $c_0=0.5$ and $\kappa _c={\rm e}^{-\eta _c}/2\sinh \eta _c$ gives the equation

(5.10)\begin{equation} \ln\frac{\kappa_c}{1+\kappa_c}+4\kappa_c+2=0, \end{equation}

which can be solved numerically for $\kappa _c$ and hence $\eta _c$. This allows us to efficiently find $\eta _c$ to arbitrary precision and it is found that $\eta _c\approx 1.1997$. This is significantly smaller than the critical throughflow number for the linear stability analysis at finite Reynolds number $\eta ^{\ell }_c\approx 3.3511$. As well as the relation in figure 14, solving (5.3b) has given us an expression for $G(\kern0.7pt y)$ and its derivative which can be used in other layers as appropriate. The values of $G'(\pm 1)$ are of great importance and have been shown in figure 16. We always have $G'(1)<0$ and $G'(-1)>0$, but the relative sizes of both values varies with $c$ and $\eta$.

Figure 16.

Values of (a) $G'(1)$ and (b) $G'(-1)$ for each solution to the Rayleigh equation shown in figure 14.

To complete our leading order description of this system, we must find $A$, $J_{CL}$ as defined in (5.5) and the expansion parameter $\varepsilon$ in terms of $\textit {Re}$. Doing so requires a leading order description of the wall layers and critical layer; we start with the former.

5.2. Wall layers

In the lower wall layer, the order one coordinate used is $Y_-=\textit {Re}^{1/2}(\kern0.7pt y+1)$ and our expansions are

(5.11a)\begin{gather} u(x,Y_-,t) =U(\kern0.7pt y)+\varepsilon^2\tilde{u}^{(-)}(\xi,Y_-)+\cdots+\varepsilon\textit{Re}^{{-}1/2}u_M^{(-)}(Y_-)+\cdots, \end{gather}

(5.11b)\begin{gather}v(x,Y_-,t)=\frac{\eta}{\textit{Re}}+\varepsilon^2\textit{Re}^{{-}1/2}\tilde{v}^{(-)}(\xi,Y_-)+\cdots, \end{gather}

(5.11c)\begin{gather}p(x,Y_-,t)=\varepsilon^2\tilde{p}^{(-)}(\xi,Y_-)+\cdots. \end{gather}

This system is subject to no penetration and no slip at $Y_-=0$ and matching to the inviscid core as $Y_-\to \infty$. It can be shown that the leading order velocities take the form

(5.12a)\begin{gather} \tilde{u}^{(-)}(\xi,Y_-) =\tfrac{1}{2}\tilde{U}^{(-)}(Y_-){\rm e}^{\mathrm{i}\xi}+\text{c.c.}, \end{gather}

(5.12b)\begin{gather}\tilde{v}^{(-)}(\xi,Y_-) =\tfrac{1}{2}\tilde{V}^{(-)}(Y_-){\rm e}^{\mathrm{i}\xi} +\text{c.c.}, \end{gather}

where $U^{(-)}(Y_-)$ and $V^{(-)}(Y_-)$ can be found analytically. We also require the mean-flow distortion solution $u_M(Y_-)$, which requires inspection of the next order of equations. Substituting (5.11) into (2.3b), the terms of order $\varepsilon ^4$ are investigated. The inertial contribution is

(5.13) \begin{align} \alpha\tilde{u}^{(-)}\frac{\partial\tilde{u}^{(-)}}{\partial\xi} +\tilde{v}^{(-)}\frac{\partial\tilde{u}^{(-)}}{\partial Y_-} =\frac{1}{4}( (\mathrm{i}\alpha\tilde{U}^{(-)2} +\tilde{V}^{(-)}\tilde{U}^{(-)\prime}) {\rm e}^{2\mathrm{i}\xi} +\tilde{U}^{(-)\prime}\overline{\tilde{V}^{(-)}} +\text{c.c.} ), \end{align}

which must be balanced by viscosity and pressure gradient effects. Crucially, (5.13) contains a term independent of $\xi$, which can only be balanced by the mean flow distortion. For this to occur, we need a balance between the $\varepsilon ^4$ terms and the leading order mean flow distortion of magnitude $\varepsilon \textit {Re}^{-1/2}$ from (5.11a) giving the scaling

(5.14)\begin{equation} \varepsilon=\textit{Re}^{{-}1/6}. \end{equation}

This scaling is typical of a nonlinear critical layer study (see for example Smith & Bodonyi Reference Smith and Bodonyi1982b; Walton Reference Walton2003; Walton & Barnes Reference Walton and Barnes2023). The lower wall analysis can then be completed by finding $u_M^{(-)}(Y_-)$. The upper wall layer uses the order one coordinate $Y_+=\textit {Re}^{1/2}(1-y)$ and the same method to find the leading order description.

5.3. Critical layer and amplitude equation

This leaves only the critical layer in which to determine the leading order solution. To complete our analysis, the first five orders must be considered which leads to a lengthy calculation. This allows us to quantify the jump in the mean flow distortion over the critical layer, $J_{CL}$. Then, by considering a net vorticity jump across all layers, an equation for the amplitude $A$ can be obtained. The details can be found in other studies (Kumar & Walton Reference Kumar and Walton2019; Walton & Barnes Reference Walton and Barnes2023) with adjustments being made due to the addition of throughflow. It is found that the jump in the derivative of the mean flow distortion over the critical layer is given by

(5.15)\begin{equation} J_{CL} =\frac{2^{5/2}{\rm \pi} q_2(\kern0.7pt y_c)\varTheta} {(\alpha q_1(\kern0.7pt y_c))^{1/2}} A^{1/2}G(\kern0.7pt y_c)^{1/2}, \end{equation}

and the amplitude is given by

(5.16) \begin{equation} A =8\alpha^{1/3}\biggl( -{\rm \pi}\varTheta q_2(\kern0.7pt y_c) \sqrt{\frac{G(\kern0.7pt y_c)}{q_1(\kern0.7pt y_c)}} \biggr)^{2/3} \left( \frac{G'(1)^2}{\mu^{(+)}}-\frac{G'({-}1)^2}{\mu^{(-)}} \right)^{{-}2/3}, \end{equation}

where

(5.17a) \begin{equation} q_n(\kern0.7pt y) =\frac{1}{n!}\frac{\mathrm{d}^nU}{\mathrm{d}y^n},\quad \mu^{(+)} =\sqrt{\frac{\alpha(1-c)}{2}},\quad \mu^{(-)} =\sqrt{\frac{\alpha c}{2}} \end{equation}

and

(5.17b) \begin{equation} \varTheta =\int_1^\infty\biggl(\int_0^{2{\rm \pi}} (s-\cos\xi)^{{-}1}\,\mathrm{d}\xi \biggr)^{{-}1} -\frac{1}{2{\rm \pi} s^{1/2}}\,\mathrm{d}s -\frac{1}{\rm \pi} \approx{-}0.3103. \end{equation}

Using the solution found to (5.3b), all terms in (5.16) are known and $A$ can be calculated. The results are shown in figure 17. Observe that the amplitude is only real-valued for a subset of the values for which we found solutions in figure 14, as this requires the final bracketed term in (5.16) to be positive. It can be seen from figure 16 that this will not always be true.

Figure 17.

Amplitude $A$ plotted for varying values of $c$ and $\eta$.

With this additional information, the relationship between $\alpha$ and $c$ is shown again in figure 18, but now it is indicated where a real amplitude exists with a solid line. As $\eta$ is increased, a solution exists for a greater proportion of each curve. This tells us that faster throughflow results in our asymptotic structure existing for a greater range of wavespeeds.

Figure 18.

Relationship found between $\alpha$ and $c$ from the Rayleigh equation (5.3b). The lines are solid where a non-negative real amplitude exists and dotted where no such amplitude exists.

By considering the $1\ll \eta \ll \mathit {O}(\textit {Re})$ limit of this asymptotic structure, it can be demonstrated that a real amplitude exists for the entire domain $0\leqslant c\leqslant 1$. Thus, asymptotically large throughflow can enable our asymptotic structure to exist for the entire domain of wavespeeds available. It is not until $\eta \sim \mathit {O}(\textit {Re})$ that the structure breaks down.