2 Derivation of the governing equations

2.1 The vortex/TS wave interaction equations

Our starting point is the incompressible Navier–Stokes equations

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{t}+(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}p+\mathit{Re}^{-1}{\rm\nabla}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

non-dimensionalized on the channel half-width and centreline velocity of the unperturbed state, where $(x,y,z,t)$ are the streamwise, normal, spanwise and temporal coordinates. In plane Poiseuille flow an imposed constant pressure gradient creates a parabolic basic flow $\boldsymbol{u}=(1-y^{2},0,0)$ . Here we consider the large-Reynolds-number framework where a three-dimensional TS wave of the lower branch type nonlinearly interacts with a roll/streak flow. The interaction is initially generated by a pair of neutral oblique TS waves which produces the symmetries

(2.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle [u,v,w,p](x,y,z)=[u,v,-w,p](x,y,-z), & \displaystyle\end{eqnarray}$$

(2.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle [u,v,w,p](x,y,z)=[u,-v,w,p](x+L_{x}/2,-y,z), & \displaystyle\end{eqnarray}$$

(2.2c ) $$\begin{eqnarray}\displaystyle & \displaystyle [u,v,w,p](x,y,z)=[u,v,w,p](x+L_{x}/2,y,z+L_{z}/2), & \displaystyle\end{eqnarray}$$

where $L_{x}$ and $L_{z}$ are the streamwise and spanwise wavelengths of the wave, respectively. When calculating numerical solutions to both the VWI and Navier–Stokes solutions we will assume these symmetries, although the majority of the formulation of the VWI equations which follows is independent of this assumption. At high Reynolds number the TS wave disturbance of lower branch type is characterized by streamwise and spanwise wavenumbers of $O(\mathit{Re}^{-1/7})$ with an associated wavespeed of $O(\mathit{Re}^{-2/7})$ : see Lin (Reference Lin1945). We therefore suppose that, for example, the normal component of the wave perturbation has the form

(2.3) $$\begin{eqnarray}{\it\delta}\,\widehat{v}(y,Z)E+\text{c.c}.;\quad E\equiv \exp [\text{i}{\it\alpha}\,(X-c\,T)],\end{eqnarray}$$

where $\text{c.c}.$ denotes complex conjugate, ${\it\delta}$ is a small parameter to be determined in terms of the Reynolds number, and $(X,Z)=\mathit{Re}^{-1/7}(x,z)$ , $T=\mathit{Re}^{-3/7}t$ are the spatial and temporal scales of the wave. We will concentrate on nonlinear travelling wave disturbances for which the scaled streamwise wavenumber ${\it\alpha}$ and phasespeed $c$ are real $O(1)$ amplitude-dependent quantities to be determined as part of the solution. We restrict our attention to solutions which have period $2{\rm\pi}/{\it\beta}$ in $Z$ , with scaled spanwise wavenumber ${\it\beta}$ to be specified. Bennett et al. (Reference Bennett, Hall and Smith1991) derived sets of interaction equations for both $O(1)$ and $O(\mathit{Re}^{1/7})$ spanwise length scales. The analysis in this section concerns the latter case and thus the spanwise wavenumber has the same scaling as the streamwise wavenumber.

We will now identify the critical size of ${\it\delta}$ in terms of $\mathit{Re}$ such that the nonlinear self-interaction of the wave induces an $O(1)$ correction to the streamwise component of velocity. In VWI parlance the components independent of $x$ constitute the vortex part of the flow, while in self-sustaining process terminology the cross-stream contribution is referred to as the roll, while the streamwise component represents a streak. Here we consider the wave and roll–streak interacting in the core, which forms the bulk of the fluid motion away from the channel walls, and where the normal variable $y$ is of $O(1)$ . We shall separately consider the wall layer where the wave is generated later. First we note that, in view of the largeness of the Reynolds number, it is possible for a tiny roll field to act as a convective mechanism for both the roll itself and an $O(1)$ streak, while still ensuring that viscous effects are present at leading order. Thus, from an inertial–viscous balance of terms in the $x$ -averaged streamwise momentum equation, it follows that the $y$ component of the roll velocity ${\sim}\,\mathit{Re}^{-1}$ in the core. Next, the continuity balance implies a $z$ component of the roll velocity ${\sim}\,\mathit{Re}^{-6/7}$ , in view of the spanwise scaling assumed above. A pressure gradient–viscous balance in the $x$ -averaged spanwise momentum equation then yields a core pressure scaling ${\sim}\,\mathit{Re}^{-12/7}$ . We now suppose that this roll–streak flow is forced by the nonlinear self-interaction of the wave. Consideration of the $x$ -averaged $y$ -momentum equation leads to the core pressure ${\sim}\,{\it\delta}^{2}$ , in view of (2.3). Comparing this with the earlier estimate for the core pressure results in

(2.4) $$\begin{eqnarray}{\it\delta}=\mathit{Re}^{-6/7}\end{eqnarray}$$

as the critical size of the wave for this strongly nonlinear interaction.

The full core flow expansion can now be expressed as

(2.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle u=U(y,Z)+{\it\epsilon}^{5}[\hat{u} (y,Z)E+\text{c.c}.]+\cdots \,, & \displaystyle\end{eqnarray}$$

(2.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle v={\it\epsilon}^{6}[\hat{v}(y,Z)E+\text{c.c}.]+{\it\epsilon}^{7}V(y,Z)+\cdots \,, & \displaystyle\end{eqnarray}$$

(2.5c ) $$\begin{eqnarray}\displaystyle & \displaystyle w={\it\epsilon}^{6}W(y,Z)+{\it\epsilon}^{7}[{\hat{w}}(y,Z)E+\text{c.c}.]+\cdots \,, & \displaystyle\end{eqnarray}$$

(2.5d ) $$\begin{eqnarray}\displaystyle & \displaystyle p=-{\it\epsilon}^{7}2x+{\it\epsilon}^{7}[\hat{p}(y,Z)E+\text{c.c}.]+{\it\epsilon}^{12}P(y,Z)+\cdots \,, & \displaystyle\end{eqnarray}$$

where we have defined ${\it\epsilon}=\mathit{Re}^{-1/7}$ . Here $U$ , $(V,W)$ and $\hat{\boldsymbol{u}}$ are the streak, roll and wave components, respectively. These expansions are then substituted into the Navier–Stokes equations (2.1). In expansions (2.5) the precise scaling of the wave component is chosen so that it is driven by the streak and satisfies the linear inviscid balances

(2.6a−d ) $$\begin{eqnarray}\text{i}{\it\alpha}\hat{u} +\hat{v}_{y}=0,\quad \text{i}{\it\alpha}U\hat{u} +\hat{v}U_{y}=0,\quad \text{i}{\it\alpha}U\hat{v}=-\hat{p}_{y},\quad \text{i}{\it\alpha}U{\hat{w}}=-\hat{p}_{Z}.\end{eqnarray}$$

We note here how, under our scaling, the smaller spanwise wave component uncouples from the $(\hat{u} ,\hat{v})$ flow. This is an expected result because, as we shall see shortly, the size of ${\hat{w}}$ increases upon approach to the wall in order to satisfy the full three-dimensional continuity equation in the near wall sublayer. Solving (2.6a−d ) up to an amplitude function $A(Z)$ we obtain

(2.6e−h ) $$\begin{eqnarray}\hat{u} =-AU_{y},\quad \hat{v}=\text{i}{\it\alpha}AU,\quad {\hat{w}}=-(\text{i}{\it\alpha}U)^{-1}\hat{p}_{Z},\quad \hat{p}_{y}={\it\alpha}^{2}AU^{2}.\end{eqnarray}$$

This leads to a pressure-amplitude (displacement) relation

(2.7) $$\begin{eqnarray}\hat{p}(1,Z)-\hat{p}(-1,Z)={\it\alpha}^{2}A\int _{-1}^{1}[U(s,Z)]^{2}\,\text{d}s.\end{eqnarray}$$

Next we turn to the equations governing $(U,V,W)$ . Owing to the scale difference between the wall normal and spanwise directions, viscous effects are negligible in the $y$ -momentum equation which, as a consequence, reduces to a simple balance between the normal pressure gradient and the nonlinear wave forcing. This can be used to eliminate the roll pressure gradient in the spanwise momentum equation and we are left with the balances

(2.8a ) $$\begin{eqnarray}\displaystyle & \displaystyle V_{y}+W_{Z}=0, & \displaystyle\end{eqnarray}$$

(2.8b ) $$\begin{eqnarray}\displaystyle & \displaystyle VU_{y}+WU_{Z}=2+U_{yy}, & \displaystyle\end{eqnarray}$$

(2.8c ) $$\begin{eqnarray}\displaystyle & \displaystyle VW_{yy}+WW_{yZ}=2{\it\alpha}^{2}(|A|^{2})_{Z}UU_{y}+W_{yyy}, & \displaystyle\end{eqnarray}$$

where (2.8c ) is arrived at by using (2.6e−h ) and differentiating with respect to $y$ to eliminate an unknown function of $Z$ . This roll/streak combination is to be solved subject to

(2.8d ) $$\begin{eqnarray}U=V=W=0\quad \text{on}~y=\pm 1.\end{eqnarray}$$

To close the system (2.8) we need to determine the amplitude function $A(Z)$ . This is found by considering the behaviour of the wave in the viscous wall layers. From the symmetry (2.2b ), it suffices for our purposes to consider only the lower wall layer at $y=-1$ . First, it is important to note from (2.8) that, when $y+1\sim h\ll 1$ , the core roll–streak components have magnitudes $U\sim O(h),V\sim O(h^{2}),W\sim O(h)$ , $P\sim O(1)$ . Thus, the behaviour of the streak and (2.6e−h ) imply the core wave sizes $\hat{u} \sim O(1),\hat{v}\sim O(h),{\hat{w}}\sim O(h^{-1})$ , $\hat{p}\sim O(1)$ in this limit. In the wall layer the viscous wave is convected by a streak of $O(h)$ . The thickness of this layer in terms of the Reynolds number can then be established to be $h\sim O({\it\epsilon}^{2})$ , by balancing the inertial operator $u\partial /\partial x\sim O({\it\epsilon}h)$ against the viscous contribution $\mathit{Re}^{-1}\partial ^{2}/\partial y^{2}\sim O({\it\epsilon}^{7}h^{-2})$ . This motivates the introduction of the lower wall layer variable, $Y$ , defined to be

(2.9) $$\begin{eqnarray}Y={\it\epsilon}^{-2}(y+1).\end{eqnarray}$$

Then, as suggested by (2.5) and the near-wall core flow behaviour mentioned above, it is found that the flow variables in the lower wall layer expand as

(2.10a ) $$\begin{eqnarray}\displaystyle & \displaystyle u={\it\epsilon}^{2}{\it\lambda}(Z)Y-{\it\epsilon}^{4}Y^{2}+{\it\epsilon}^{5}[\bar{u}(Y,Z)E+\text{c.c}.]+\cdots \,, & \displaystyle\end{eqnarray}$$

(2.10b ) $$\begin{eqnarray}\displaystyle & \displaystyle v={\it\epsilon}^{8}[\bar{v}(Y,Z)E+\text{c.c}.]+{\it\epsilon}^{11}\bar{V}(Y,Z)+\cdots \,, & \displaystyle\end{eqnarray}$$

(2.10c ) $$\begin{eqnarray}\displaystyle & \displaystyle w={\it\epsilon}^{5}[\bar{w}(Y,Z)E+\text{c.c}.]+{\it\epsilon}^{8}\bar{W}(Y,Z)+\cdots \,, & \displaystyle\end{eqnarray}$$

(2.10d ) $$\begin{eqnarray}\displaystyle & \displaystyle p=-{\it\epsilon}^{7}2x+{\it\epsilon}^{7}[\bar{p}(Z)E+\text{c.c}.]+{\it\epsilon}^{12}\bar{P}(Z)+\cdots \,, & \displaystyle\end{eqnarray}$$

where ${\it\lambda}(Z)$ is the streak shear stress on the lower wall, defined by

(2.11) $$\begin{eqnarray}{\it\lambda}(Z)=\partial U/\partial y\quad \text{on }y=-1,\end{eqnarray}$$

and the second term in (2.10a ) balances the constant pressure gradient. The wave components satisfy the following linearized unsteady boundary-layer-type equations:

(2.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}{\it\alpha}\bar{u}+\bar{v}_{Y}+\bar{w}_{Z}=0, & \displaystyle\end{eqnarray}$$

(2.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}{\it\alpha}({\it\lambda}Y-c)\bar{u}+{\it\lambda}\bar{v}+{\it\lambda}_{Z}Y\bar{w}=-\text{i}{\it\alpha}\bar{p}+\bar{u}_{YY}, & \displaystyle\end{eqnarray}$$

(2.12c ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}{\it\alpha}({\it\lambda}Y-c)\bar{w}=-\bar{p}_{Z}+\bar{w}_{YY}, & \displaystyle\end{eqnarray}$$

where $\bar{p}=\bar{p}(Z)=\hat{p}(-1,Z)$ ; we note that the wave in the wall layer satisfies the full continuity equation in contrast to the core equations (2.6a−d ). Note also that the $\mathit{Re}^{-1/7}$ wavenumber scaling of the lower branch ensures the balancing of the pressure gradient term in the streamwise equation (2.12b ). The solution of equations (2.12a–c ) must satisfy the no slip boundary conditions

(2.12d ) $$\begin{eqnarray}\bar{u}=\bar{v}=\bar{w}=0\quad \text{on }Y=0,\end{eqnarray}$$

and match to the core flow solution (2.6e−h ) via the outer conditions

(2.12e,f ) $$\begin{eqnarray}\bar{u}\rightarrow -{\it\lambda}A,\quad \bar{w}\sim -\frac{\bar{p}_{Z}}{\text{i}{\it\alpha}{\it\lambda}Y}\quad \text{as }Y\rightarrow \infty .\end{eqnarray}$$

Here we can see that the quantity $A(Z)$ represents an effective wall displacement for the streamwise wave component due to the presence of the viscous layers. The corresponding upper wall layer solutions can be found by using symmetries (2.2b ) and in particular we find $\hat{p}(1,Z)=-\hat{p}(-1,Z)$ . Equations (2.7) and (2.12) can be manipulated into one equation for $A(Z)$ (see Smith Reference Smith1979), namely

(2.13) $$\begin{eqnarray}(JA)^{\prime \prime }-({\it\lambda}^{\prime }/{\it\lambda})\mathscr{F}({\it\xi})(JA)^{\prime }-{\it\alpha}^{2}(JA)=2({\it\alpha}{\it\lambda})^{5/3}\mathscr{G}({\it\xi})A/{\it\alpha}^{2},\end{eqnarray}$$

where in the above we have defined

(2.14) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\it\xi}(Z)=-\frac{(\text{i}{\it\alpha}{\it\lambda})^{1/3}}{{\it\lambda}}c,\quad \mathscr{F}({\it\xi})=\frac{3}{2}+\frac{{\it\xi}({\it\xi}{\it\kappa}({\it\xi})+\text{Ai}^{\prime }({\it\xi}))}{2\text{Ai}({\it\xi})},\quad \mathscr{G}({\it\xi})=\text{i}^{5/3}\frac{\text{Ai}^{\prime }({\it\xi})}{{\it\kappa}({\it\xi})},\\ \displaystyle {\it\kappa}({\it\xi})=\int _{{\it\xi}}^{\infty }\text{Ai}(s)\,\text{d}s,\quad J(Z)=\int _{-1}^{1}[U(y,Z)]^{2}\,\text{d}y,\end{array}\right\}\end{eqnarray}$$

with ${\it\lambda}(Z)$ given in (2.11) and $\text{Ai}$ representing the usual Airy function. Equation (2.13) is a complex nonlinear eigenvalue problem that determines the wavenumber and phasespeed of the nonlinear wave part of the flow-field. A full derivation of (2.13), (2.14) from (2.7), (2.12) is given in the appendix A.

The interaction problem has now been reduced to solving the roll–streak equations (2.8), coupled with the wave equation (2.13), with the two systems of equations being linked by the amplitude function $A(Z)$ and the streak shear ${\it\lambda}(Z)$ . This coupling allows us to describe the self-sustaining mechanism behind the VWI states. First, the streak $U$ drives the viscous wall wave displacement $A$ in equation (2.13) through the wall shear stress ${\it\lambda}(Z)$ . Second, the displacement creates an inviscid core wave of amplitude $A$ which drives the roll $(V,W)$ in (2.8a ) and (2.8c ) through the nonlinear self-interaction of the wave. Finally, closing the loop, the roll drives the streak through the advection terms in (2.8b ). This provides a framework in which equilibrium states can be described.

We must also prescribe the amplitude $\mathscr{A}$ of the wave, which we define as

(2.15) $$\begin{eqnarray}\mathscr{A}^{2}=\int _{0}^{2{\rm\pi}/{\it\beta}}|A(Z)|^{2}\,\text{d}Z.\end{eqnarray}$$

The behaviour of the solution for small $\mathscr{A}$ is investigated in the next section. It should be noted that an alternative approach would be to prescribe the axial wavenumber ${\it\alpha}$ and determine the corresponding amplitude: the choice made here is more convenient with regard to the computation of our asymptotic problem. To give a physical sense to the size of the states we also introduce the average flux deficit of the perturbation:

(2.16) $$\begin{eqnarray}\mathscr{Q}=-\frac{{\it\beta}}{4{\rm\pi}}\int _{0}^{2{\rm\pi}/{\it\beta}}\int _{-1}^{1}\widetilde{U}(y,Z)\,\text{d}y\,\text{d}Z,\end{eqnarray}$$

where $\widetilde{U}$ is the perturbation to the basic state streak, and is defined by

(2.17) $$\begin{eqnarray}\widetilde{U}(y,Z)=U(y,Z)-(1-y^{2}).\end{eqnarray}$$

2.2 The near-wall roll–streak equations

The previous subsection fully describes the interaction to leading order, determining the wavenumber and phasespeed for given $\mathscr{A}$ , and the flow in the core. However, since the thickness of the wall layer is predicted to be of size $O({\it\epsilon}^{2})=O(\mathit{Re}^{-2/7})$ , this remains a significant part of the flow-field, even at quite large Reynolds numbers. As we will see in states computed from the full Navier–Stokes equations, the largest wall layer structures are found in the roll components. From substitution of (2.10) into the Navier–Stokes equations (2.1), this part of the flow is governed by the equations

(2.18a ) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{V}_{Y}+\bar{W}_{Z}=0, & \displaystyle\end{eqnarray}$$

(2.18b ) $$\begin{eqnarray}\displaystyle & \displaystyle [\text{i}{\it\alpha}\bar{u}^{\ast }\bar{w}+\bar{v}\bar{w}_{Y}^{\ast }+\bar{w}\bar{w}_{Z}^{\ast }+\text{c.c}.]=\bar{W}_{YY}, & \displaystyle\end{eqnarray}$$

where $\ast$ represents complex conjugate. These equations are to be solved subject to no slip on the wall and far field matching conditions. A consideration of the roll part of the flow in the core of the channel reveals that the appropriate matching condition is

(2.18c ) $$\begin{eqnarray}\bar{W}\sim W_{y}(-1,Z)Y\quad \text{as }Y\rightarrow \infty .\end{eqnarray}$$

The near-wall roll flow can therefore be computed once the core flow is found from the interaction equations (2.8) and (2.13). We will return to these equations in § 5.2, when we compare VWI states with Navier–Stokes computations across the entire channel width.

3 Linear and weakly nonlinear analysis

In general the nonlinear nature of the interaction equations makes any analytical progress impossible. However, when the amplitude of the wave is small, the interaction is in a weakly nonlinear state, and we can take advantage of this to construct solutions where the perturbation to the basic flow is small, and the streamwise wavenumber and phasespeed are close to the linear neutral values for the undisturbed flow. This allows us to find analytic expressions which describe how the solutions bifurcate from the linear neutral point into finite amplitude states and also provides a means of validating numerical results as we will see in § 4.

We consider small-amplitude solutions, $|\mathscr{A}|\sim {\it\Delta}\ll 1$ , where ${\it\alpha},c$ are close to their linear neutral values. The forcing in the core is of $O({\it\Delta}^{2})$ from (2.8c ), and so it follows that the perturbations away from the basic flow are also of this order. Therefore, we expand

(3.1) $$\begin{eqnarray}A(Z)={\it\Delta}\cos {\it\beta}Z+{\it\Delta}^{3}(A_{3}^{(1)}\cos {\it\beta}Z+A_{3}^{(3)}\cos 3{\it\beta}Z)+\cdots \,,\end{eqnarray}$$

and, upon substitution into (2.8), we find the resulting perturbed roll–streak flow. For example the streak expands in the form

(3.2) $$\begin{eqnarray}U(y,Z)=1-y^{2}+{\it\Delta}^{2}U_{2}(y)\cos 2{\it\beta}Z+\cdots \,,\end{eqnarray}$$

where $U_{2}$ is a polynomial in $y$ that can be determined analytically. We find similarly that ${\it\alpha},c$ expand as

(3.3a,b ) $$\begin{eqnarray}{\it\alpha}={\it\alpha}_{0}+{\it\Delta}^{2}{\it\alpha}_{2}+\cdots \,,\quad c=c_{0}+{\it\Delta}^{2}c_{2}+\cdots \,,\end{eqnarray}$$

where $({\it\alpha}_{0},c_{0})$ corresponds to the linear neutral point, which satisfies the equation

(3.4) $$\begin{eqnarray}15(2{\it\alpha}_{0})^{5/3}\mathscr{G}({\it\xi}_{0})+8{\it\alpha}_{0}^{2}({\it\alpha}_{0}^{2}+{\it\beta}^{2})=0,\end{eqnarray}$$

from substitution of (3.1), (3.2) into (2.13), with the function $\mathscr{G}$ given in (2.14). Here the quantity ${\it\xi}_{0}=-(2\text{i}{\it\alpha}_{0})^{1/3}c_{0}/2$ is the leading order value of the quantity ${\it\xi}$ defined in (2.14). Taking the imaginary part of (3.4), we find that the imaginary part of $\mathscr{G}({\it\xi}_{0})$ is zero and thus ${\it\xi}_{0}$ and $\mathscr{G}({\it\xi}_{0})$ are determined independently of the spanwise wavenumber. Letting ${\it\xi}_{0}=-\text{i}^{1/3}s_{0}$ , i.e.  $s_{0}=(2{\it\alpha}_{0})^{1/3}c_{0}/2$ , we find

(3.5a,b ) $$\begin{eqnarray}s_{0}\simeq 2.297,\quad \text{and}\quad \mathscr{G}({\it\xi}_{0})\simeq -1.001,\end{eqnarray}$$

with each given to three decimal places, and we note that $\mathscr{G}({\it\xi}_{0})$ is purely real as required. The dispersion relation (3.4) defines the lower branch linear neutral point derived (for the two-dimensional case) by Lin (Reference Lin1945) and given in his equation (13.2). For example, when ${\it\beta}=1$ , we find

(3.6) $$\begin{eqnarray}({\it\alpha}_{0},c_{0})\simeq (1.9426,2.9225).\end{eqnarray}$$

The linear stability problem at finite Reynolds numbers can be reduced to solving the Orr–Sommerfeld equation, which was performed here using the Chebyshev collocation method: see Orszag (Reference Orszag1971). In figure 1, for a scaled spanwise wavenumber ${\it\beta}=1$ , we compute the neutral value of the scaled streamwise wavenumber on the lower branch of the linear neutral curve from the solution of the Orr–Sommerfeld equation. The dotted curve shows the Orr–Sommerfeld quantity as a function of Reynolds number, while the solid horizontal line is the asymptotic value ${\it\alpha}={\it\alpha}_{0}$ established in (3.6). Convergence to the large Reynolds number result is observed, although good agreement is not obtained until $\mathit{Re}\gtrsim 10^{7}$ . This suggests that our nonlinear theory, relying as it does on this asymptotic framework, will only begin to exhibit quantitative agreement with Navier–Stokes solutions when $\mathit{Re}$ becomes similarly large. We will see in § 5.2 that this is indeed the case.

Figure 1. The lower branch of the linear neutral curve for scaled spanwise wavenumber ${\it\beta}=1$ . The solid line represents the asymptotic result (3.6), while the dashed line shows the results of Orr–Sommerfeld computations.

Figure 2. The perturbation to the scaled streamwise wavenumber as a function of scaled spanwise wavenumber ${\it\beta}$ from the numerical solution of the weakly nonlinear equation (3.7).

One of the most important consequences of this analysis is that it describes the nature of the bifurcation from the linear neutral point. The leading order analysis at $O({\it\Delta})$ gives the dispersion relation (3.4), then at $O({\it\Delta}^{3})$ we find, by equating terms proportional to $\cos {\it\beta}Z$ , the solvability condition which relates $({\it\alpha}_{2},c_{2})$ , namely

(3.7) $$\begin{eqnarray}{\it\nu}_{1}{\it\alpha}_{2}+{\it\nu}_{2}c_{2}=1,\end{eqnarray}$$

where ${\it\nu}_{1}$ and ${\it\nu}_{2}$ are complex constants which depend on ${\it\beta}$ . For example, by solving (3.7) we find that for ${\it\beta}=1$ , the bifurcation for small amplitudes is described (to $O({\it\Delta}^{2})$ ) by

(3.8a,b ) $$\begin{eqnarray}{\it\alpha}\simeq 1.9426-0.0090{\it\Delta}^{2},\quad c\simeq 2.9225+0.0057{\it\Delta}^{2}.\end{eqnarray}$$

Numerical computations of (3.7) (figure 2) indicate that in (3.3), ${\it\alpha}_{2}<0$ for all ${\it\beta}$ , i.e. the bifurcation appears to be subcritical for all spanwise wavenumbers.

4 Numerical method for solving the VWI equations

The equations to be solved are (2.8), (2.13) subject to some prescribed wave amplitude, (2.15). We will seek solutions which are periodic in $Z$ with period $2{\rm\pi}/{\it\beta}$ where ${\it\beta}$ is the spanwise wavenumber. The equations are coupled through the wall shear stress ${\it\lambda}(Z)$ , defined in (2.11), and wave forcing terms. Equation (2.8a ) motivates the introduction of a stream function ${\it\Psi}$ , where

(4.1) $$\begin{eqnarray}(V,W)=(\partial _{Z},-\partial _{y}){\it\Psi},\end{eqnarray}$$

and (2.8a )–(2.8c ) now become

(4.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\Psi}_{Z}U_{y}-{\it\Psi}_{y}U_{Z}=2+U_{yy}, & \displaystyle\end{eqnarray}$$

(4.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\Psi}_{y}{\it\Psi}_{yyZ}-{\it\Psi}_{Z}{\it\Psi}_{yyy}=FUU_{y}-{\it\Psi}_{yyyy}, & \displaystyle\end{eqnarray}$$

(4.2c ) $$\begin{eqnarray}F(Z)=2{\it\alpha}^{2}\left(|A|^{2}\right)_{Z}.\end{eqnarray}$$

It is convenient to split the streak into the basic flow and a perturbation denoted by $\widetilde{U}$ , defined in (2.17). The boundary conditions for the perturbed variables are the no-slip conditions

(4.3) $$\begin{eqnarray}\widetilde{U}={\it\Psi}={\it\Psi}_{y}=0\quad \text{on }y=\pm 1,\end{eqnarray}$$

with the use of symmetry (2.2a ) allowing us to impose that ${\it\Psi}$ vanishes on both walls.

The most natural choice, given the flow geometry, is to use spectral methods. We use a Chebyshev/Fourier basis in the wall normal/spanwise direction with $T_{n}$ representing the $n$ th Chebyshev polynomial. The flow variables are thus represented by the series

(4.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \widetilde{U}=\mathop{\sum }_{n=1}^{N}\mathop{\sum }_{k=1}^{K}\widetilde{U}_{n,k}T_{n-1}(y)\cos [2(k-1){\it\beta}Z],\\ \displaystyle {\it\Psi}=\mathop{\sum }_{n=1}^{N}\mathop{\sum }_{k=1}^{K}{\it\Psi}_{n,k}T_{n-1}(y)\sin [2k{\it\beta}Z],\quad A=\mathop{\sum }_{k=1}^{K}A_{k}\cos [(2k-1){\it\beta}Z],\end{array}\right\}\end{eqnarray}$$

with the values of $N$ and $K$ determined by requiring sufficient decay in the spectral coefficients. In order to discretize the system we apply the collocation method in the wall-normal direction with collocation points at $y_{j}=\cos [j{\rm\pi}/N]$ for $j=0,\ldots ,N$ , using $j=1,\ldots ,N-1$ for (4.2a ), and $j=2,\ldots ,N-2$ for (4.2b ), with the remaining points reserved for the application of the no slip boundary conditions (4.3). We apply the Galerkin method in the spanwise direction. A pseudo-spectral approach is used to deal with the nonlinear terms, and de-aliasing is performed by zero-padding.

In the current formulation, the equations are invariant under an arbitrary streamwise phase shift. To account for this degree of freedom we set the imaginary part of $A_{1}$ to be zero. We now have $NK$ unknowns for both $\tilde{U}$ and ${\it\Psi}$ , $2K$ for the real and imaginary parts of $A$ , and $2$ for ${\it\alpha},c$ , resulting in a total of $2NK+2K+2$ unknowns. Then for the core equations (4.2) (including the boundary conditions) we have $2NK$ equations, $2K$ equations for (2.13), and 2 equations from (2.15) and the phase lock condition, resulting in a total of $2NK+2K+2$ nonlinear coupled equations. This system is solved using a quasi-Newton method, where the Jacobian matrix required as part of the scheme is approximated with finite differences. The Newton iteration is continued until the L1-norm of the residual of the discretized equations becomes less than $10^{-8}$ .

Resolution tests were performed in which the amplitude $\mathscr{A}$ was fixed and the number of spectral modes used in the wall-normal and spanwise directions were varied in order to verify that the numerical results were independent of the number of modes used. An example of such a resolution test is described here. This was performed at unit spanwise wavenumber and amplitude $\mathscr{A}^{2}=20$ , which is sufficiently large that the interactions are nonlinear. The streamwise wavenumber was chosen as a representative quantity to monitor as the number of modes was varied. Changing the number of Fourier modes, $K$ , was observed to have the largest effect, but once $K$ exceeded 50 the change in ${\it\alpha}$ was in the 7th decimal place, and by the time $K$ exceeded 110 the change was in the 12th decimal place. The structure of the solution in the wall-normal direction is moderate and a value of $N$ in the range $20$ – $25$ was found to be sufficient for all of the results in this work.

Although some key quantities such as ${\it\alpha},c,\mathscr{A}$ and $\mathscr{Q}$ converge excellently, it was found that an unphysical oscillation contaminates the flow field because of the development of a sharp structure in $Z$ , as we shall see in the next section. The behaviour of the Fourier coefficients shown in figure 3 for larger values of $K$ clearly shows that the slightly upward ‘tail’, where the coefficients near $k=K$ assume relatively large values, is responsible for the unphysical oscillation. We shall see in § 6 that the tail is caused by the loss of regularity in the VWI equations but that, provided the wave amplitude is below a critical value, this irregularity does not prevent the VWI equations from providing the correct leading-order behaviour. The analysis presented will confirm that for the values of $\mathscr{A}$ under consideration, the apparently converged Fourier coefficients ahead of the tail provide a good approximation to the solution of the VWI equations, thereby justifying the neglect of those coefficients associated with the tail.

Figure 3. The spanwise spectral convergence of the numerical solution for $\mathscr{A}^{2}=20$ . The curves with different line styles represent numerical solutions of the VWI equations with different numbers of Fourier modes $K$ , with the legend showing the value of $K$ .

Figure 4. Amplitude $\mathscr{A}$ of the nonlinear state as a function of streamwise wavenumber ${\it\alpha}$ . The solid curve represents the fully nonlinear numerical results from the VWI equations for ${\it\beta}=1$ . The dashed curve is the corresponding weakly nonlinear analysis, which is valid for small $\mathscr{A}$ .

Figure 5. (a,c,e) Represent the core perturbation streak $\widetilde{U}$ of the VWI solutions and (b,d,f) show the accompanying wall shear stress ${\it\lambda}(Z)$ (solid curve) and the wave displacement function $|A|^{2}$ (dotted line): (a,b)  $\mathscr{A}^{2}=1~(\mathscr{Q}\approx 5.7\times 10^{-6})$ ;(c,d)  $\mathscr{A}^{2}=20~(\mathscr{Q}\approx 1.8\times 10^{-3})$ ; (e,f)  $\mathscr{A}^{2}=30~(\mathscr{Q}\approx 3.8\times 10^{-3})$ .

5 Numerical results and comparison

5.1 Computational results for the strongly nonlinear vortex/TS wave interaction

First we will focus on a scaled spanwise wavenumber ${\it\beta}=1$ , although the behaviour for other spanwise wavenumbers, which will be discussed shortly, is qualitatively similar. For small amplitudes, solutions are found to have the structure predicted by the weakly nonlinear analysis as shown in figure 4. Figure 5(a,b) shows the solution for a representative small amplitude $\mathscr{A}^{2}=1$ , where the interaction is still in the weakly nonlinear regime with colour denoting the size of the perturbation to the streak in the core. The structure is mild, with a small spanwise modulation in the roll–streak giving rise to a fairly uniform vortex structure. Increasing the amplitude further results in stronger wave forcing, inducing larger perturbations to the flow and fully nonlinear interactions arise. The numerical scheme described in § 4 allows us to continue the solution branch to much larger amplitudes compared with the typical maximum amplitude $\mathscr{A}\approx 1$ that proved attainable when employing a classical fixed-point iteration technique in an earlier version of our code. The solutions shown in figure 5(c,d) for $\mathscr{A}^{2}=20$ become localized in the spanwise direction, with the wall shear stress developing steep gradients at particular spanwise locations. The flow undergoes further localization for $\mathscr{A}^{2}=30$ as shown in figure 5(e,f), and an ever larger number of Fourier modes is required to achieve satisfactory resolution. For larger amplitudes the resolution required for the spectral convergence becomes beyond our computational resources. Here we remark that we encounter this difficulty before the vortex structure has become fully isolated. In this sense our use of the word ‘localization’ is different from that used in recent numerical investigations of fully isolated vortex structures, e.g. Schneider, Gibson & Burke (Reference Schneider, Gibson and Burke2010), Deguchi et al. (Reference Deguchi, Hall and Walton2013) and Gibson & Brand (Reference Gibson and Brand2014).

Figure 6. The perturbation streak $\widetilde{U}$ from the numerical solution of the VWI equations for (a) ${\it\beta}=2$ , (b) ${\it\beta}=3$ , (c) ${\it\beta}=5$ and (d) the large ${\it\beta}$ limit (5.2). A scaled amplitude $\widehat{\mathscr{A}}\equiv {\it\beta}^{-5}\mathscr{A}=0.1$ is chosen for these plots.

Similar computations are now presented for larger values of ${\it\beta}$ . The solutions at all of these spanwise wavenumbers are found to be qualitatively similar to the ${\it\beta}=1$ case, exhibiting localization in the spanwise coordinate as the amplitude is increased. In figure 6(a–c) we show the perturbation to the streak in the core for ${\it\beta}=2,3$ and  $5$ . For these plots the amplitude is selected as $\mathscr{A}=0.1{\it\beta}^{5}$ . It can be seen that with this choice of amplitude there is very little difference between the solutions for ${\it\beta}=3$ and ${\it\beta}=5$ . Indeed beyond ${\it\beta}=5$ the solution becomes practically independent of the choice of ${\it\beta}$ in terms of this scaling. The reason for this behaviour at large ${\it\beta}$ is that under the rescaling

(5.1a−f ) $$\begin{eqnarray}\widehat{Z}={\it\beta}Z,\quad \widehat{{\it\Psi}}={\it\beta}{\it\Psi},\quad \widehat{{\it\alpha}}={\it\beta}^{6}{\it\alpha},\quad \widehat{c}={\it\beta}^{-2}c,\quad \widehat{A}={\it\beta}^{-5}A,\quad \widehat{\mathscr{A}}={\it\beta}^{-5}\mathscr{A},\end{eqnarray}$$

the roll/streak equations (4.2) remain unchanged in terms of the newly scaled variables, while in (2.13) the final term on the left-hand side becomes ${\it\beta}^{-14}\widehat{{\it\alpha}}^{2}J\widehat{A}$ . This rescaling suggests that in the limit ${\it\beta}\rightarrow \infty$ the VWI remains governed essentially by (4.2), (2.13) with ${\it\beta}=1$ but with the final term on the left-hand side of (2.13) absent, i.e. we have the modified system

(5.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{{\it\Psi}}_{\widehat{Z}}U_{y}-\widehat{{\it\Psi}}_{y}U_{\widehat{Z}}=2+U_{yy}, & \displaystyle\end{eqnarray}$$

(5.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{{\it\Psi}}_{y}\widehat{{\it\Psi}}_{yy\widehat{Z}}-\widehat{{\it\Psi}}_{\widehat{Z}}\widehat{{\it\Psi}}_{yyy}=\widehat{F}UU_{y}-\widehat{{\it\Psi}}_{yyyy}, & \displaystyle\end{eqnarray}$$

(5.2c ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{F}(\widehat{Z})=2\widehat{{\it\alpha}}^{2}|\widehat{A}|_{\widehat{Z}}^{2}, & \displaystyle\end{eqnarray}$$

(5.2d ) $$\begin{eqnarray}\displaystyle & \displaystyle (J\widehat{A})_{\widehat{Z}\widehat{Z}}-({\it\lambda}_{\widehat{Z}}/{\it\lambda})\mathscr{F}(\widehat{{\it\xi}})(J\widehat{A})_{\widehat{Z}}=2(\widehat{{\it\alpha}}{\it\lambda})^{5/3}\mathscr{G}(\widehat{{\it\xi}})\widehat{A}/\widehat{{\it\alpha}}^{2}, & \displaystyle\end{eqnarray}$$

where we now seek solutions with period $2{\rm\pi}$ in $\widehat{{\it\beta}}$ . It is also worth remarking here that in this large spanwise wavenumber limit, the quantity $\widehat{{\it\xi}}=-(\text{i}\widehat{{\it\alpha}}{\it\lambda})^{1/3}\widehat{c}/{\it\lambda}$ remains $O(1)$ . Equations (5.2) were solved computationally, by the same numerical method described in § 4, but now prescribing $\widehat{\mathscr{A}}$ (defined in (5.1f )) and tracking the state which bifurcates from the linear neutral point to finite values of $\widehat{\mathscr{A}}$ . In figure 6(d) we present the streak arising from the solution of (5.2) for $\widehat{\mathscr{A}}=0.1$ . The flow field can be seen to be almost indistinguishable from the ${\it\beta}=5$ case shown in figure 6(c). In figure 7 we see a plot of $\widehat{\mathscr{A}}$ as a function of $\widehat{{\it\alpha}}$ for ${\it\beta}=2,2.5,3$ and $5$ , along with the corresponding result from the solution of (5.2). Convergence to this solution with increasing ${\it\beta}$ can clearly be seen in this figure, with comparisons for other flow quantities exhibiting a similar convergence.

Figure 7. The scaled amplitude $\widehat{\mathscr{A}}$ of the VWI solutions as a function of the scaled wavenumber $\widehat{{\it\alpha}}$ for various finite values of spanwise wavenumber ${\it\beta}$ and the large- ${\it\beta}$ limit.

Equations (5.2) break down when a distinguished limit ${\it\beta}=O(\mathit{Re}^{1/7})$ is reached, which corresponds to $O(1)$ spanwise variations in $z$ . In this limit the large- ${\it\beta}$ scaling (5.1) implies that ${\it\alpha}\rightarrow O(\mathit{Re}^{-1})$ . Simultaneously, the streamwise wave component grows to be as large as the streak, while in the spanwise and wall-normal directions, both the wave and roll become $O(\mathit{Re}^{-1})$ in size. The interaction is therefore governed by the so-called boundary-region equations in which the roll, streak and wave now evolve over a single $O(\mathit{Re})$ streamwise lengthscale. Such an asymptotic problem was recently solved for plane Couette flow by Deguchi et al. (Reference Deguchi, Hall and Walton2013) where the solution was found to strongly localize in the spanwise direction while also exhibiting streamwise localization as the amplitude of the motion was increased. Although we have not investigated the solution of the corresponding boundary-region problem for plane Poiseuille flow, we would expect to observe similar localization for this flow configuration.

5.2 Detailed visualization and comparison with Navier–Stokes computations

In this section we present detailed flow structures computed from the VWI equations and compare the results with finite Reynolds number calculations of the Navier–Stokes equations, restricting our attention to a spanwise wavenumber ${\it\beta}=1$ . The travelling wave solutions of the full Navier–Stokes equations are directly solved by applying Newton’s method to the discretized system. The discretization is performed using a spectral method utilizing a Chebyshev expansion in the wall-normal direction and a Fourier series representation in the streamwise and spanwise directions. More details on the scheme can be found in Deguchi et al. (Reference Deguchi, Hall and Walton2013). The primal seed for the Newton method is a pair of neutral oblique TS waves numerically obtained from the linear stability analysis. A similar starting condition was used by Ehrenstein & Koch (Reference Ehrenstein and Koch1991): however, in our case typically $300$ Fourier modes are taken in the spanwise direction, compared with the $2$ used by the aforementioned authors. In addition, $60$ Chebyshev modes are used in the normal direction to fully resolve the solutions. It is widely regarded that the number of modes in Ehrenstein & Koch (Reference Ehrenstein and Koch1991) is insufficient to obtain reliable solutions even at the moderate Reynolds numbers they considered. Here our computations confirm that well-resolved solutions can be continued to high Reynolds number and finite amplitude. By way of contrast, four or fewer streamwise Fourier modes are typically sufficient to preserve the accuracy of the solutions, consistent with the fact that in VWI only one streamwise Fourier mode is present at leading order. All Navier–Stokes quantities are subjected to the VWI scalings in order to compare with the results given in the last subsection.

Figure 8. Comparison of the solution branch for the VWI with finite Reynolds numbers results from the full Navier–Stokes equations at ${\it\beta}=1$ : (a) phasespeed $c$ versus streamwise wavenumber ${\it\alpha}$ ; (b) the average flux deficit $\mathscr{Q}$ versus ${\it\alpha}$ . The solid circle in (a) is the linear critical point (3.6).

Figure 8(a) shows a comparison of the streamwise wavenumber–phasespeed curves for the VWI (solid curve) and the full Navier–Stokes computations for various Reynolds numbers (dashed curves). In the figure we can see the full Navier–Stokes results eventually converge to the VWI results with increasing Reynolds number, thereby confirming that vortex/TS wave interaction is indeed operational in high-Reynolds-number shear flows. Figure 8(b) shows a similar comparison involving the flux $\mathscr{Q}$ where the convergence to the VWI state can also be seen as the Reynolds number increases.

Figure 9. Comparison of a Navier–Stokes computation and the corresponding VWI state for scaled spanwise wavenumber ${\it\beta}=1$ , showing $U$ in (a,d), $V$ in (b,e), and $W$ in (c,f): (a–c) Navier–Stokes result at $\mathit{Re}=10^{8}$ ; (d–f) composite solution for the VWI state at $\mathit{Re}=10^{8}$ . Both results were computed with $\mathscr{Q}\approx 3\times 10^{-4}$ , which corresponds to $\mathscr{A}^{2}=8$ .

Figure 10. The core and wall layer parts of the VWI solution for $\mathscr{A}^{2}=8$ , ${\it\beta}=1$ : (a–d) the normal and spanwise components of the roll flow; (a,b) show the core solutions $V$ and $W$ ; (c,d) depict the wall layer solutions $\bar{V}$ and $\bar{W}$ .

In figure 9, at $\mathscr{Q}\approx 3\times 10^{-4}$ (which corresponds to $\mathscr{A}^{2}=8$ ), we compare the roll–streak fields for $Re=10^{8}$ . In figure 9(a–c) we show the appropriate components computed from the full Navier–Stokes equations. In the flow field we see near-wall structures which are particularly large for the roll, and especially for the spanwise component $W$ . In the VWI framework these structures are not captured in the core flow to leading order, and in fact are confined to the wall layer as seen in figure 10, where the VWI wall layer solution has been obtained from the numerical solution of (2.18). In these figures we can see how the roll flow matches between the two regions. Adding the core and rescaled wall layer solutions and then subtracting the common part of them, we can construct a composite solution for a given value of $\mathit{Re}$ . Returning now to figure 9(d–f), we show the VWI state throughout the whole channel, created by the composite solution technique. Excellent agreement can be found by comparing the composite solution with the corresponding Navier–Stokes result shown in figure 9(a–c).

Figure 11. Isosurface of streamwise vorticity near the bottom wall for $Re=10^{8}$ , ${\it\beta}=1$ and $\mathscr{Q}\approx 3\times 10^{-4}$ . The red (light grey)/blue (dark grey) surface represents 40 % of the maximum/minimum value: (a) Navier–Stokes solution; (b) VWI solution.

Figure 12. The perturbed streak arising from the Navier–Stokes solution for $Re=10^{8}$ and ${\it\beta}=1$ . The value of flux is chosen as $\mathscr{Q}\approx 3.8\times 10^{-3}$ to compare with the VWI version in figure 5(e).

The wave part of the solution dominates the streamwise vorticity as we can see from the asymptotic expansion (2.9) in the wall layers. In figure 11(a,b) we compare the streamwise vorticity arising from the Navier–Stokes solution with that from the VWI computation. Here we only show a small interval in $y$ near the bottom wall since the major vortex activity is concentrated in the near-wall regions. In the bottom VWI plot we only use the leading order term proportional to $\bar{w}_{Y}$ and its conjugate. The figure displays an excellent agreement between the Navier–Stokes and the asymptotic results, showing a pair of flat counter-rotating vortices emerging from the wall to sustain the core structure.

With the comparison made here, we are satisfied that vortex/TS-wave interaction is successfully describing the large $\mathit{Re}$ asymptotic limit of nonlinear travelling wave states bifurcating from the linear neutral curve in plane Poiseuille flow. For other values of amplitude similar agreement can be observed. In particular, the spanwise localization seen in the previous section is also operational in the full Navier–Stokes computations: see figure 12 in which the Navier–Stokes streak field is shown at $\mathscr{Q}\approx 3.8\times 10^{-3}$ . This figure should be compared with figure 5(e) where the corresponding VWI solution is shown. In the next section we will uncover the details of the localization process and discover that the VWI equations gradually lose regularity at certain spanwise locations as the wave amplitude increases, with a singularity developing at a critical value of $\mathscr{A}$ .

6 The details of the localization process

As we have seen in the previous section, the flow field begins to localize as the wave amplitude increases, with steep gradients developing around particular spanwise locations. This is evident, for example, in figure 5(b,d,f) where there are two such locations: $Z=Z_{1}={\rm\pi}/2{\it\beta}$ and $Z=Z_{2}=3{\rm\pi}/2{\it\beta}$ . At both of these points we have $A=0$ . In this section we shall provide a comprehensive theoretical explanation for the origin of the localization and demonstrate how the VWI solution gradually loses regularity at these local points as the flow becomes increasingly nonlinear, eventually producing a singularity at a critical value of $\mathscr{A}$ which we will determine. The key to the development of the singularity is that, although the full VWI system is elliptic in $Z$ , once the wavenumber and phasespeed are determined from the wave equation (2.13) subject to periodicity requirements, the calculation of the roll–streak field for a given wave forcing arises from the solution of a system (4.2) which is only first order in $Z$ . It is therefore then sufficient to provide the conditions

(6.1a−c ) $$\begin{eqnarray}\displaystyle U|_{Z=Z_{1}}=U|_{Z=Z_{2}},\quad {\it\Psi}_{Z}|_{Z=Z_{1}}={\it\Psi}_{Z}|_{Z=Z_{2}},\quad A|_{Z=Z_{1}}=A|_{Z=Z_{2}}=0, & & \displaystyle\end{eqnarray}$$

in order to solve the interaction problem in $Z\in [Z_{1},Z_{2}]$ . As a consequence, continuity of the higher-order derivatives is not guaranteed in general and, as we shall see, the degree of regularity of the solution can indeed be predicted for a given wave amplitude by analysing the behaviour near the local point in detail.

Without loss of generality we can set $Z=0$ at the local point by applying an appropriate spanwise shift. By the symmetries (2.2), if there is no singularity we must have a regular expansion at $Z=0$ of the form

(6.2a,b ) $$\begin{eqnarray}\displaystyle U=U_{0}(y)+U_{2}(y)Z^{2}+\cdots \,,\quad {\it\Psi}=Z\{{\it\Psi}_{0}(y)+{\it\Psi}_{2}(y)Z^{2}+\cdots \}, & & \displaystyle\end{eqnarray}$$

(6.2c,d ) $$\begin{eqnarray}\displaystyle F=Z\{F_{0}+F_{2}Z^{2}+\cdots \,\},\quad A=Z\{A_{0}\text{e}^{\text{i}{\it\theta}_{0}}+A_{2}\text{e}^{\text{i}{\it\theta}_{2}}Z^{2}+\cdots \}, & & \displaystyle\end{eqnarray}$$

where $U_{m}(y)$ and ${\it\Psi}_{m}(y)$ are even and odd functions, respectively.

Let us suppose now that only the first $2M$ spanwise derivatives of $U$ are continuous at $Z=0$ . We therefore pose the expansions

(6.3a,b ) $$\begin{eqnarray}U=\mathop{\sum }_{m=0}^{M}U_{2m}Z^{2m}+U_{{\it\sigma}}|Z|^{{\it\sigma}}+\cdots \,,\quad {\it\Psi}=Z\left\{\mathop{\sum }_{m=0}^{M}{\it\Psi}_{2m}Z^{2m}+{\it\Psi}_{{\it\sigma}}|Z|^{{\it\sigma}}+\cdots \right\},\end{eqnarray}$$

(6.3c,d ) $$\begin{eqnarray}F=Z\left\{\mathop{\sum }_{m=0}^{M}F_{2m}Z^{2m}+F_{{\it\sigma}}|Z|^{{\it\sigma}}+\cdots \right\},\quad A=Z\left\{\mathop{\sum }_{m=0}^{M}A_{2m}\text{e}^{\text{i}{\it\theta}_{2m}}Z^{2m}+A_{{\it\sigma}}\text{e}^{\text{i}{\it\theta}_{{\it\sigma}}}|Z|^{{\it\sigma}}+\cdots \right\},\vphantom{\mathop{\sum }_{\sum }}\end{eqnarray}$$

${\it\sigma}$

$M$

$2M<{\it\sigma}$

(6.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\Psi}_{0}U_{0}^{\prime }=2+U_{0}^{\prime \prime }, & \displaystyle\end{eqnarray}$$

(6.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\Psi}_{0}^{\prime }{\it\Psi}_{0}^{\prime \prime }-{\it\Psi}_{0}{\it\Psi}_{0}^{\prime \prime \prime }=F_{0}U_{0}U_{0}^{\prime }-{\it\Psi}_{0}^{\prime \prime \prime \prime }, & \displaystyle\end{eqnarray}$$

(6.4c ) $$\begin{eqnarray}U_{0}={\it\Psi}_{0}={\it\Psi}_{0}^{\prime }=0\quad \text{on }y=\pm 1,\end{eqnarray}$$

while the first irregular contributions satisfy the balances

(6.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle U_{{\it\sigma}}^{\prime \prime }+{\it\sigma}{\it\Psi}_{0}^{\prime }U_{{\it\sigma}}-{\it\Psi}_{0}U_{{\it\sigma}}^{\prime }-({\it\sigma}+1)U_{0}^{\prime }{\it\Psi}_{{\it\sigma}}=0, & \displaystyle\end{eqnarray}$$

(6.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\Psi}_{{\it\sigma}}^{\prime \prime \prime \prime }+\{{\it\Psi}_{0}^{\prime \prime }{\it\Psi}_{{\it\sigma}}+({\it\sigma}+1){\it\Psi}_{0}^{\prime }{\it\Psi}_{{\it\sigma}}^{\prime \prime }\}-\{{\it\Psi}_{0}{\it\Psi}_{{\it\sigma}}^{\prime \prime }+({\it\sigma}+1){\it\Psi}_{0}^{\prime \prime }{\it\Psi}_{{\it\sigma}}\} & \displaystyle \nonumber\\ \displaystyle & \displaystyle \quad =F_{{\it\sigma}}U_{0}U_{0}^{\prime }+F_{0}(U_{0}^{\prime }U_{{\it\sigma}}+U_{0}U_{{\it\sigma}}^{\prime }), & \displaystyle\end{eqnarray}$$

(6.5c ) $$\begin{eqnarray}U_{{\it\sigma}}={\it\Psi}_{{\it\sigma}}={\it\Psi}_{{\it\sigma}}^{\prime }=0\quad \text{on }y=\pm 1.\end{eqnarray}$$

Similarly, by substitution of (6.3) into the wave equation (2.13) we find

(6.6) $$\begin{eqnarray}({\it\sigma}+1)\{J_{{\it\sigma}}A_{0}^{2}+J_{0}A_{{\it\sigma}}A_{0}\text{e}^{\text{i}{\it\theta}_{{\it\sigma}}}\}-\frac{{\it\lambda}_{{\it\sigma}}}{{\it\lambda}_{0}}\{\mathscr{F}_{0}^{(\text{r})}+\text{i}\mathscr{F}_{0}^{(\text{i})}\}J_{0}A_{0}^{2}=0,\end{eqnarray}$$

where

(6.7) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\it\lambda}_{0}=U_{0}^{\prime }|_{y=-1},\quad {\it\lambda}_{{\it\sigma}}=U_{{\it\sigma}}^{\prime }|_{y=-1},\quad \mathscr{F}_{0}^{\text{(r)}}+\text{i}\mathscr{F}_{0}^{\text{(i)}}=\mathscr{F}|_{Z=0},\\ \displaystyle J_{0}=\int _{-1}^{1}U_{0}^{2}\,\text{d}y,\quad J_{{\it\sigma}}=\int _{-1}^{1}2U_{0}U_{{\it\sigma}}\,\text{d}y,\\ \displaystyle F_{0}=4{\it\alpha}^{2}A_{0}^{2},\quad F_{{\it\sigma}}=4({\it\sigma}+2){\it\alpha}^{2}A_{0}A_{{\it\sigma}}\cos ({\it\theta}_{{\it\sigma}}-{\it\theta}_{0}).\end{array}\right\}\end{eqnarray}$$

Using the expressions for $F_{0}$ and $F_{{\it\sigma}}$ given above, we can rewrite (6.6) as

(6.8) $$\begin{eqnarray}F_{{\it\sigma}}=\frac{({\it\sigma}+2)}{({\it\sigma}+1)}\left\{\frac{{\it\lambda}_{{\it\sigma}}}{{\it\lambda}_{0}}\mathscr{F}_{0}^{(\text{r})}-({\it\sigma}+1)\frac{J_{{\it\sigma}}}{J_{0}}\right\}F_{0}.\end{eqnarray}$$

Substituting this expression for $F_{{\it\sigma}}$ into equations (6.5), we find that the equation governing the irregular contribution can be represented in the form

(6.9) $$\begin{eqnarray}L\left[\begin{array}{@{}c@{}}U_{{\it\sigma}}\\ {\it\Psi}_{{\it\sigma}}\end{array}\right]=0,\end{eqnarray}$$

where $L$ is a linear operator which depends on $({\it\sigma},F_{0},{\it\alpha},c)$ . Thus, $L$ must be singular to have non-trivial $U_{{\it\sigma}},{\it\Phi}_{{\it\sigma}}$ . Standard techniques from linear algebra can then be used to find the smallest value of ${\it\sigma}$ that produces a zero eigenvalue of $L$ for a given $(F_{0},{\it\alpha},c)$ , and hence to determine the degree of regularity of the solution for a given wave amplitude.

Figure 13. The solid curve denotes the zeroth order local analysis solution. The crosses are full VWI results for ${\it\beta}=1$ , with selected values of $\mathscr{A}^{2}$ indicated by the arrows.

Figure 14. The value of ${\it\sigma}$ associated with the singularity at the local point.

First, as shown in figure 13, we compare the leading order local solution (computed from (6.4) for a given $F_{0}$ using a Runge–Kutta/shooting approach) with the full VWI results computed in the last section. In the full VWI computation $F_{Z}$ and ${\it\lambda}$ evaluated at the local point correspond to the quantities $F_{0}$ and ${\it\lambda}_{0}$ , respectively. Excellent agreement can be seen in the figure, thereby confirming the validity of both analyses. Note that the parameter $F_{0}$ is effectively a measure of the nonlinearity of the system.

Next we examine the regularity of the solution. In figure 14 we plot the smallest value of ${\it\sigma}$ which gives a zero eigenvalue of $L$ for a given $F_{0}$ . The values of $({\it\alpha},c)$ used in the calculation are extracted from the full VWI computation and are dependent on the particular value of $F_{0}$ under consideration. In this way, global effects enter into our local analysis. For $F_{0}\rightarrow 0$ , i.e. the zero amplitude limit, the value of ${\it\sigma}$ tends to infinity and the solution is regular. As $F_{0}$ is increased, the value of ${\it\sigma}$ decreases and as a consequence the VWI solution gradually loses regularity at $Z=0$ . For example, provided $F_{0}\gtrsim 90$ we have ${\it\sigma}<2$ and ${\it\lambda}_{ZZ}$ becomes singular, while if $F_{0}\gtrsim 220$ we find that ${\it\sigma}<1$ and now in addition ${\it\lambda}_{Z}$ develops a singularity at $Z=0$ . Finally when $F_{0}=F_{c}$ with $F_{c}$ in the range 750–800, the value of ${\it\sigma}$ becomes zero and the size of the irregular terms in the expansions (6.3) become comparable with the leading-order terms. This value of $F_{c}$ corresponds to a wave amplitude $\mathscr{A}^{2}$ in the range 30–35 and is consistent with the fact that beyond this value of $\mathscr{A}^{2}$ there is poor agreement between the full VWI solution and the leading-order local analysis even when as many as $2000$ Fourier modes are used in the VWI computation.

Provided we are in the regime where $F_{0}<F_{c}$ , and therefore the solution is still regular at leading order, it is possible to remove the singularity occurring in the higher-order terms by inserting a thin layer around $Z=0$ in which spanwise diffusion of the roll flow is now as important as diffusion in the normal direction. The balancing of these two effects in the Navier–Stokes equations requires $\partial _{yy}\sim {\it\epsilon}^{2}\partial _{ZZ}$ and hence implies $Z\sim O({\it\epsilon})$ and therefore an $O(1)$ spanwise thickness in terms of $z$ . A simple order of magnitude analysis reveals that the appropriate balances in the spanwise diffusion layer are the $O(1)$ spanwise wavelength equations given in Bennett et al. (Reference Bennett, Hall and Smith1991) and referred to as a ‘square vortex’ type interaction in the subsequent literature. However here we should employ matching conditions to the outer flow in the spanwise direction rather than the periodic conditions used in previous studies. The appropriate expansions in the diffusion layer are

(6.10) $$\begin{eqnarray}\displaystyle \left[\begin{array}{@{}c@{}}U\\ {\it\Psi}\\ \hat{p}_{y}\end{array}\right]=\left[\begin{array}{@{}c@{}}U_{0}(y)\\ {\it\epsilon}z{\it\Psi}_{0}(y)\\ {\it\epsilon}{\it\alpha}^{2}U_{0}^{2}zA_{0}\text{e}^{\text{i}{\it\theta}_{0}}\end{array}\right]+\cdots +\left[\begin{array}{@{}c@{}}{\it\epsilon}^{{\it\sigma}}{\check{u}}(y,z)\\ {\it\epsilon}^{{\it\sigma}+1}\check{{\it\psi}}(y,z)\\ {\it\epsilon}^{{\it\sigma}+1}{\it\alpha}^{2}{\check{a}}(y,z)\end{array}\right]+\cdots \,, & & \displaystyle\end{eqnarray}$$

where the leading-order terms $U_{0},{\it\Psi}_{0}$ can be determined from the local problem (6.4), which remains valid in the spanwise diffusion layer. At $O({\it\epsilon}^{{\it\sigma}})$ we find that the components associated with the singular behaviour of the outer flow simply satisfy the square vortex interaction equations linearized about the leading order terms in the expansion (6.10). In principle we can solve such linear equations subject to the usual no-slip conditions in $y$ and the following matching conditions for large $|z|$ :

(6.11a−c ) $$\begin{eqnarray}{\check{u}}\sim |z|^{{\it\sigma}}U_{{\it\sigma}},\quad \check{{\it\psi}}\sim z|z|^{{\it\sigma}}{\it\Psi}_{{\it\sigma}},\quad {\check{a}}\sim z|z|^{{\it\sigma}}(A_{{\it\sigma}}U_{0}^{2}\text{e}^{\text{i}{\it\theta}_{{\it\sigma}}}+2A_{0}U_{0}U_{{\it\sigma}}\text{e}^{\text{i}{\it\theta}_{0}})\quad \text{as }z\rightarrow \pm \infty ,\end{eqnarray}$$

where $U_{{\it\sigma}}(y)$ , ${\it\Psi}_{{\it\sigma}}(y)$ , $A_{{\it\sigma}}$ , ${\it\theta}_{{\it\sigma}}$ are fixed by the outer solution.

An important consequence of the spanwise diffusion layer analysis is that when $F_{0}<F_{c}$ , the solution ( ${\check{u}},\check{{\it\psi}},{\check{a}}$ ) satisfies a linear problem and does not influence the determination of the leading-order outer problem. The local structure is therefore completely passive, in the sense that we can solve the outer solution independently as we have done in the previous section. In the numerical computation of the full VWI problem, a Fourier expansion is employed in $Z$ , implicitly assuming the regularity of the solution. Thus, in the computation the singularity is smoothed because a finite number of bases are used and this is responsible for the creation of the ‘tail’ in the Fourier spectrum due to Gibbs oscillations. Except for this unphysical oscillation, our scheme accurately captures the solution of the outer problem outside of the numerically thin layer whose thickness is comparable with the frequency of the finest Fourier mode; see figure 3 for the convergence of the Fourier spectrum aside from the unphysical ‘tail’.

When $F_{0}$ reaches the critical value $F_{c}$ (and ${\it\sigma}$ approaches zero), the size of the singular part of the local expansion (6.10) becomes comparable with that of the leading-order term and thus the scenario elucidated here breaks down. In this limit the spanwise diffusion layer equations will become nonlinear and may play an active role in the entire interaction problem. At such large amplitudes the present VWI formulation is no longer valid and a new interactive structure must come into play.