2. Formulation of the problem

Consider the Navier–Stokes equations in Cartesian coordinates $(x,y,z)$

(2.1) \begin{eqnarray} (\partial _t +\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}=-\boldsymbol{\nabla }\!p+R^{-1}{\nabla} ^2 \boldsymbol{u}+R^{-1}Q\boldsymbol{e}_x,\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \end{eqnarray}

where $t$ denotes time, $\boldsymbol{e}_x$ represents the unit vector in the $x$ -direction,  and $\boldsymbol{\nabla}$ is the gradient operator. The velocity $\boldsymbol{u}=(u,v,w)$ and the pressure $p$ are defined in a channel with $y\in [-1,1]$ . The combination of the boundary conditions $\boldsymbol{u}=(\pm 1,0,0)$ at $y=\pm 1$ and the normalised pressure gradient $Q=0$ corresponds to plane Couette flow, whose base flow is $\boldsymbol{u}=(y,0,0)$ . Imposing $\boldsymbol{u}=\boldsymbol{0}$ at $y=\pm 1$ and $Q=2$ yields plane Poiseuille flow, with base flow $\boldsymbol{u}=(1-y^2,0,0)$ . For both cases, our Reynolds number $R$ is based on the characteristic velocity of the base flow and the channel half width.

We enforce periodicity in $z$ and quasi-periodicity in $x$ and $t$ . More specifically, the flow is assumed to depend on $y$ , $z$ , $\theta _m$ the phase variables, and $m$ denotes their index $\theta _m=\alpha _m x+\omega _mt$ , $m=1,2,\ldots ,M$ , each taken to be $2\pi$ periodic. Here, $\alpha _m$ and $\omega _m$ are real numbers representing the wavenumbers and frequencies, respectively. We use $\beta$ to indicate the $2\pi /\beta$ periodicity of the flow in $z$ . Travelling-wave, periodic and quasi-periodic solutions in a standard periodic box can all be captured within this set-up.

Now, let us apply a ‘Reynolds decomposition’ $\boldsymbol{u}=\overline {\boldsymbol{u}}+\tilde {\boldsymbol{u}}$ and $p=\overline {p}+\tilde {p}$ , defined by averaging over all phase variables. Here and hereafter, the overline denotes this averaging operation (with respect to the variables $\theta _1,\theta _2,\ldots$ ), while the tilde indicates the fluctuation (or wave) component. The streamwise mean velocity $\overline {u}$ represents the streaks (i.e. spanwise modulation of the streamwise velocity), while the cross-stream components $\overline {v}$ and $\overline {w}$ capture the rolls (i.e. streamwise vortices). Substituting the decomposition into (2.1) and applying the averaging yields

(2.2) \begin{eqnarray} \big[\overline {v}\partial _y +\overline {w}\partial _z-R^{-1}(\partial _y^2+\partial _z^2)\big] \left [ \begin{array}{c} \overline {u}\\ \overline {v}\\ \overline {w} \end{array} \right ] + \left [ \begin{array}{c} -R^{-1}Q\\[3pt] \partial _y \overline {p} \\[3pt] \partial _z \overline {p} \end{array} \right ] = \boldsymbol{F},\,\,\, \partial _y\overline {v}+\partial _z\overline {w}=0, \end{eqnarray}

where $\boldsymbol{F}=-\overline {(\tilde {\boldsymbol{u}}\boldsymbol{\cdot }{\boldsymbol{\nabla }})\tilde {\boldsymbol{u}}}$ is the Reynolds stress. Taking the difference between (2.1) and (2.2) gives the fluctuating part of the momentum equations

(2.3) \begin{eqnarray} (\partial _t+\overline {\boldsymbol{u}}\boldsymbol{\cdot }{\boldsymbol{\nabla }})\tilde {\boldsymbol{u}}+ (\tilde {\boldsymbol{u}}\boldsymbol{\cdot }{\boldsymbol{\nabla }})\overline {\boldsymbol{u}} + \big\{(\tilde {\boldsymbol{u}}\boldsymbol{\cdot }{\boldsymbol{\nabla }})\tilde {\boldsymbol{u}}-\overline {(\tilde {\boldsymbol{u}}\boldsymbol{\cdot }{\boldsymbol{\nabla }})\tilde {\boldsymbol{u}}} \big\} =-{\boldsymbol{\nabla }}\! \tilde {p}+R^{-1}{{\nabla} }^2 \tilde {\boldsymbol{u}},\,\,\, \end{eqnarray}

and ${\boldsymbol{\nabla }}\boldsymbol{\cdot }\tilde {\boldsymbol{u}}=0$ . At this stage, no approximations have been made to the equations.

The terms in curly brackets in (2.3) are the wave–wave interactions, which we hereafter refer to as WW. Removing them makes (2.3) linear in the fluctuations. Substituting the normal mode expansions $\tilde {\boldsymbol{u}}=\sum _{m=1}^M\hat {\boldsymbol{u}}_m(y,z)\exp (\text{i}\theta _m)+\text{c.c.}$ and $\tilde {p}=\sum _{m=1}^M \hat {p}_m(y,z)\exp (\text{i}\theta _m)+\text{c.c.}$ , with c.c. indicating the complex conjugate, the equations governing $\hat {\boldsymbol{u}}_m=(\hat {u}_m,\hat {v}_m,\hat {w}_m)$ and $\hat {p}_m$ can be readily obtained as

(2.4a) \begin{align} \mathscr{L}\hat {u}_m+\hat {v}_m\partial _y \overline {u}+\hat {w}_m\partial _z \overline {u}=-\text{i}\alpha _m \hat {p}_m, \\[-28pt] \nonumber \end{align}

(2.4b) \begin{align} \mathscr{L}\hat {v}_m+\hat {v}_m\partial _y \overline {v}+\hat {w}_m\partial _z \overline {v}=-\partial _y \hat {p}_m, \\[-28pt] \nonumber \end{align}

(2.4c) \begin{align} \mathscr{L}\hat {w}_m+\hat {v}_m\partial _y \overline {w}+\hat {w}_m\partial _z \overline {w}=-\partial _z \hat {p}_m, \\[-28pt] \nonumber \end{align}

(2.4d) \begin{align} \text{i}\alpha _m \hat {u}_m+\partial _y\hat {v}_m+\partial _z\hat {w}_m=0, \\[6pt] \nonumber \end{align}

where $\mathscr{L}=\text{i}\alpha _m(\overline {u}-c_m)+\overline {v}\partial _y+\overline {w}\partial _z-R^{-1}(\partial _y^2+\partial _z^2-\alpha _m^2)$ . The phase speeds $c_m=-\omega _m/\alpha _m$ are real, so the waves are neutral. At large $R$ , the behaviour of the above equations is well approximated by the Rayleigh equation, which has a singularity at the critical level where $\overline {u}=c_m$ (see Hall & Sherwin Reference Hall and Sherwin2010).

Among the terms in the wave (2.4), we refer to those involving the roll components $\overline {v},\overline {w}$ as the WR (wave–roll interaction) terms. Also, the streamwise component of the Reynolds stress

(2.5) \begin{eqnarray} \boldsymbol{F}= -\sum _{m=1}^M\left [ \begin{array}{c} \partial _y(\hat {u}_m\hat {v}_m^*)+\partial _z(\hat {u}_m\hat {w}_m^*) \\[5pt] \partial _y(\hat {v}_m\hat {v}_m^*)+\partial _z(\hat {v}_m\hat {w}_m^*)\\[5pt] \partial _y(\hat {w}_m\hat {v}_m^*)+\partial _z(\hat {w}_m\hat {w}_m^*)\end{array} \right ]+\text{c.c.}, \end{eqnarray}

is referred to as the SR (streamwise component of the Reynolds stress) term. The QL-VWI approach solves (2.4) and (2.2) with the right-hand side replaced by (2.5). In Beaume et al. (Reference Beaume, Chini, Julien and Knobloch2015), Deguchi & Hall (Reference Deguchi and Hall2016) and Pausch et al. (Reference Pausch, Yang, Hwang and Eckhardt2019), the same system was solved but with $M=1$ ; the approach of Blackburn et al. (Reference Blackburn, Hall and Sherwin2013) is a further simplification in which the WR and SR terms are neglected. The common feature between the above models and the QL model mentioned in § 1 is that the WW terms are omitted. However, while the former assumes the waves to be neutral and seeks approximations to ECS, the latter allows for the temporal evolution of the waves directly through numerical integration, so that the solution to the fluctuation equations corresponds to neutral Floquet/Lyapunov mode (for periodic/chaotic case). There are numerous variants of the QL type model; for example, Thomas et al. (Reference Thomas, Lieu, Jonanovic, Farrell, Ioannou and Gayme2014) considered multiple waves within a certain wavenumber band, whereas Alizard & Biau (Reference Alizard and Biau2019) used only a single wave component. In the generalised QL approach, some of the WW terms are retained (see Hernandez et al. Reference Hernandez, Yang and Hwang2022).

Asymptotic analysis of the quasi-periodic solutions of the Navier–Stokes equations naturally leads to a multiple-wave extension of VWI. This conclusion remains unchanged even when the governing equations are replaced by QL-VWI, meaning that our model retains all of the dominant terms associated with VWI-type coherent structures at large Reynolds numbers. The asymptotic expansion has leading-order scalings $\overline {u}=O(R^0),{} \overline {v}=O(R^{-1}), \overline {w}=O(R^{-1})$ , $\hat {\boldsymbol{u}}_m=O(\delta ^{1/2}R^{-1})=O(R^{-7/6})$ , except within the critical layers of thickness $\delta =R^{-1/3}$ . These layers emerge around the critical levels to regularise the singularity of the Rayleigh equation. Within its corresponding critical layer, the $m$ th wave is amplified as $\hat {\boldsymbol{u}}_m=O(\delta ^{-1/2}R^{-1})=O(R^{-5/6})$ . The streamwise vorticity is most pronounced within the critical layer and is therefore a suitable quantity for visualising the waves. The Reynolds stress due to the waves is concentrated within the critical layer and influences the outer roll flow through stress jump conditions across the critical level (see Hall & Sherwin Reference Hall and Sherwin2010 for further details). These jump conditions are difficult to impose numerically; therefore, Blackburn et al. (Reference Blackburn, Hall and Sherwin2013) derived their model using viscous regularisation of the VWI system. The interaction of two mirror-symmetric oblique Tollmien–Schlichting waves can also be computed using a QL approach (Dempsey et al. Reference Dempsey, Deguchi, Hall and Walton2016). However, a single Tollmien–Schlichting wave develops a nonlinear critical layer that requires retaining the WW terms at leading order (see Deguchi & Walton Reference Deguchi and Walton2018 and references therein).

It should be emphasised that, after substituting the asymptotic expansion, the leading-order system is determined solely by the assumption of large $R$ . At this stage, we cannot freely choose which terms to keep or drop, unlike in other modelling approaches. This order-of-magnitude comparison confirms that the WW, WR and SR terms are indeed negligible. We also remark that the VWI theory can be extended so that the overlined components evolve on a slow time scale while the tilde components evolve on a fast time scale (see Hall & Sherwin Reference Hall and Sherwin2010). However, the simultaneous temporal evolution of roll–streaks and waves makes the separation between fast and slow scales unclear. The asymptotic analysis is therefore justified only for ECS; yet their properties can indirectly influence the more complex surrounding dynamics, as suggested by dynamical systems theory (see Page et al. Reference Page, Norgaard, Brenner and Kerswell2024; Wang et al. Reference Wang, Ayats, Deguchi, Meseguer and Mellibovsky2025).

The results in the next section are obtained by solving the QL-VWI system using the Newton–Raphson method. The parameters $R$ , $\alpha _m$ and $\beta$ are prescribed, whereas $c_m$ are determined as part of the solution. We use a Chebyshev–Fourier spectral method for discretisation, which allows us to adapt the fully nonlinear Navier–Stokes solver of Deguchi, Hall & Walton (Reference Deguchi, Hall and Walton2013) with only minor modifications. To fully resolve the sharp critical-layer structure, we typically employ up to 300 Chebyshev polynomials in the $y$ -direction and up to 50 Fourier harmonics in the $z$ -direction. At the selected points, we increased these numbers to 400 and 60, respectively, and confirmed that the results are well converged.

3. Numerical results

Section 3.1 shows that the parameter range of plane Couette flow studied in Deguchi & Hall (Reference Deguchi and Hall2016) serves as an ideal test ground of the QL-VWI approach. The key observation is that the gentle periodic orbit (GPO) solution found by Kawahara & Kida (Reference Kawahara and Kida2001) exhibits a standing-wave-like structure, realised through two counter-propagating symmetric waves. Section 3.2 aims to compute the QL-VWI approximation of quasi-periodic states in plane Poiseuille flow. Such quasi-periodic solutions are difficult to obtain directly from the full Navier–Stokes equations, yet they are essential for forming multiscale structures.

3.1. Gentle periodic orbit in plane Couette flow

Figure 1 examines how varying the streamwise period of the periodic box $[0,{2\pi }/{\alpha }]\times [-1,1]\times [0, {2\pi }/{\beta }]$ affects the ECS at $R=3000$ . The black open circles show the drag on the wall $D$ of the steady solution of the Navier–Stokes equations known as NBW (Nagata Reference Nagata1990; Clever & Busse Reference Clever and Busse1992; Waleffe Reference Waleffe1998). When $\alpha$ is increased (i.e. the streamwise period of the computational box is shortened) above that of the minimal box (Hamilton et al. Reference Hamilton, Kim and Waleffe1995), the NBW disappears at a saddle-node bifurcation point and is replaced by the GPO (filled circles), which emanates from a degenerate Hopf bifurcation.

Figure 1.

Comparison of Navier–Stokes results (symbols) and QL-VWI results (lines) in plane Couette flow at $R=3000$ . Our choice of $\beta =1.667$ ensures that, at $\alpha =1.14$ (indicated by the vertical line), the solutions fit the minimal box used in Hamilton et al. (Reference Hamilton, Kim and Waleffe1995). The black plots show the drag on the lower wall $D$ , i.e. the $x$ – $z$ – $t$ average of $\partial _y u$ at $y=-1$ . Dashed lines and open circles represent the Nagata-Busse-Waleffe state (NBW), while solid lines and filled circles represent the GPO. The frequency $\omega$ of the GPO is shown by the blue plots. The left and right colour maps show $\tilde {u}_{\textit{rms}}$ of the NBW ( $\alpha =1$ ) and the GPO ( $\alpha =1.5$ ), respectively, computed using the QL-VWI. The plotted domain spans $y\in [-1,1]$ and $z\in [0,2\pi /\beta ]$ . The green lines denote the critical levels. The grey contours are the level curves of $\overline {u}$ .

As reported by Deguchi & Hall (Reference Deguchi and Hall2016), the NBW result is well approximated by the QL-VWI with $M=1$ , $\alpha _1=\alpha$ (black dashed line). The black solid line in the figure illustrates our novel result: the QL-VWI system with $M=2$ , $\alpha _1=\alpha _2=\alpha$ , provides an excellent approximation of the GPO drag. This solution branch is obtained via bifurcation analysis of the NBW. The two waves in the approximated GPO exhibit a symmetry, satisfying $c_1=-c_2$ . The blue open circles show the frequency $\omega =|c_1\alpha _1|$ of the QL-VWI solution. This agrees well with $\omega =2\pi /T$ computed from the period $T$ of the full Navier–Stokes GPO, providing unequivocal validation of our approximation. We repeated a similar QL-VWI computation at $R=20\,000$ and found that, for $\alpha \lt 1.7$ , the result at $R=3000$ is already asymptotically converged. Computing the full Navier–Stokes GPO at $R=20\,000$ is challenging, but we observed excellent agreement with the QL-VWI solution at several representative values of $\alpha$ (e.g. for $\alpha =1.25$ , the error in the approximation of $D$ and $\omega$ is less than 0.2 %).

The colour maps of figure 1 show the root-mean-square of the streamwise wave component, $\tilde {u}_{\textit{rms}}(y,z)=(2\sum _{m=1}^M|\hat {u}_m|^2)^{1/2}$ . As the left inset shows, the NBW is supported by a steady wave that is concentrated around the critical level where the streak $\overline {u}$ vanishes (green dashed line). On the other hand, as shown in the right inset, the GPO has two critical levels at $\overline {u}=c_1=-0.162$ (green dashed line) and $\overline {u}=c_2=0.162$ (green dot-dashed line). At the bifurcation point, the two critical levels merge into one, and the GPO coincides with the NBW.

Figure 2.

Flow visualisation of the NBW (a) and the GPO (b) at $R=20\,000$ obtained by the QL-VWI. The colour maps denote the streamwise velocity of the streak field and the blue/green isosurfaces denote the 50 % maximum/minimum value of the streamwise vorticity of the fluctuation, $\partial _y\tilde {w}-\partial _z \tilde {v}$ . The dark yellow isosurfaces show the 50 % maximum value of $|\tilde {u}|$ .

As $R$ increases, the wave becomes increasingly concentrated near the critical level, consistent with the VWI theory. Figure 2 shows three-dimensional plots of the QL-VWI solutions obtained at $R=20\,000$ . The dark yellow transparent surfaces represent the magnitude of the streamwise velocity, $|\tilde {u}|$ , while the red and blue surfaces show the streamwise vorticity, $\partial _y\tilde {w}-\partial _z \tilde {v}$ . The Navier–Stokes results are almost indistinguishable from those shown in figure 2 and are therefore omitted. This observation implies that the model captures the self-sustained process of the solution: the streaks induce (single or multiple) waves, which in turn drive the rolls, and the rolls then modify the streaks through the lift-up mechanism.

3.2. Multiscale solution in plane Poiseuille flow

Our starting point consists of two Navier–Stokes travelling-wave solutions. The first is the solution independently discovered by Gibson & Brand (Reference Gibson and Brand2014) and Deguchi (Reference Deguchi2015), which we call TWA in this paper. The second is the solution recently found by Song & Deguchi (Reference Song and Deguchi2025), referred to as MS-A3, which we denote TWB here for simplicity. Throughout this section, we set $R=10^5$ .

Figure 3.

Bifurcation diagram of plane Poiseuille flow solutions. Lines show the QL-VWI results and symbols are the Navier–Stokes results. The colour maps use the same format as figure 1; all panels show the domain $(y,z) \in [-1,0]\times [0,2\pi /8]$ , except for one: (a) TWA with $\alpha /\beta =3/4$ ; (b) TWB with $\alpha /\beta =1/2$ ; (c) $M=2$ QL-VWI computation with $\alpha _1/\beta =3/4=6/8$ and $\alpha _2/\beta =3/2=12/8$ ; (d) $M=3$ QL-VWI solution with $\alpha _1/\beta =6/8$ , $\alpha _2/\beta =12/8$ and $\alpha _3/\beta =7/2=28/8$ . The colour map of the solution indicated by the blue diamond shows the domain $(y,z) \in [-1,0]\times [0,2\pi /6]$ .

Figure 3(a) shows the bifurcation diagram of TWA, obtained by varying $\beta$ while fixing the streamwise–spanwise aspect ratio of the wave, $\alpha /\beta$ . The symbols show the phase speed $c$ of the travelling waves, and the solid lines are the result obtained using the $M=1$ version of the QL-VWI system (write $c_1=c$ , $\alpha _1=\alpha$ hereafter). As observed by Deguchi (Reference Deguchi2015), TWA in the high-wavenumber regime exhibits coherent structures localised near the centre of the channel. The colour map in panel (a) shows the same type of visualisation as in figure 1, restricted to the lower half of the channel. By symmetry of the solution, a pair of streaks appears within one $z$ -period. Figure 3(b) shows analogous computational results for TWB. We find that, in the high-wavenumber range, this solution localises near the wall, as illustrated by the colour map.

Note that the colour maps in panels (a) and (b) are plotted over the same domain, $[-1,0]\times [0,2\pi /8]$ . Here, TWA at the green triangle corresponds to $\beta =8$ and $\alpha =(3/4)\beta =6$ , while TWB at the green circle has $\beta =24$ and $\alpha =(1/2)\beta =12$ . Thus, the colour map for the latter solution includes 3 copies within the domain.

The next step is to construct a new QL-VWI solution by superposing these two travelling-wave solutions. To reduce computational cost, it is desirable that the target solution possesses symmetry about the channel centreline; therefore, we also include a $y$ -symmetric copy of TWB, denoted TWBr, in the superposition. In the QL-VWI with $M=2$ , we set $\beta =8$ , $\alpha _1=6$ and $\alpha _2=12$ , and make the following substitutions:

1. (i) assign the fluctuation part of TWA to $\hat {\boldsymbol{u}}_1, \hat {p}_1$ ;

2. (ii) assign the sum of the fluctuation parts of TWB and TWBr to $\hat {\boldsymbol{u}}_2, \hat {p}_2$ ;

3. (iii) assign the sum of the mean parts of TWA, TWB and TWBr to $\overline {\boldsymbol{u}}, \overline {p}$ ;

4. (iv) assign the phase speeds of TWA and TWB to $c_1$ and $c_2$ , respectively.

After the substitution, the residual of the equations is small but non-zero, as the system is nonlinear. We treat this residual as a forcing term, multiply it by a constant $A$ , and gradually reduce $A$ as a continuation parameter in the Newton–Raphson tracking of the solution branch (i.e. a homotopy approach, see figure 4 a). The resulting solution is the green square in figure 3(c). Note that the substitution of TWB and TWBr actually uses three repeated copies of each solution so that they fit into the $\beta =8$ computational domain. This repetitiveness is, of course, slightly broken after the homotopy.

Figure 4.

Homotopy continuation connecting the superposed state ( $A=1$ ) to the unforced state ( $A=0$ ). Solution symbols are the same as in figure 3.

The solution branch in figure 3(c) is obtained by fixing $\alpha _1/\beta$ and $\alpha _2/\beta$ . After the turning point, the critical layer of the second wave moves toward the channel centre. This creates space near the wall, allowing us to insert the localised solution at the red circle in figure 3(b) ( $(\alpha ,\beta )=(28,56)$ , so its colour map contains 7 repeated copies). The superposition of the red square solution, the red circle solution and its symmetric counterpart can be carried out in QL-VWI with $M=3$ . By choosing $(\alpha _1,\alpha _2,\alpha _3)=(6,12,28)$ and applying the homotopy method (figure 4 b), we obtain the red diamond solution shown in figure 3(d).

Figure 5.

The flow field of the solution at the blue diamond in figure 3(d). (a) Visualisation in the box $[0,2\pi /1.5]\times [-1,0.25] \times [0,2\pi /6]$ . See figure 2 caption for definitions of the surfaces and colour map. (b) Energy of the streak associated with the $l$ th spanwise Fourier mode. (c) The $x$ – $z$ – $t$ mean streamwise velocity $U(y)$ . The green vertical lines indicate the projected locations of the critical layers, whose horizontal positions correspond to the associated phase speeds.

Note that the QL-VWI approximation breaks down at the Kolmogorov scale (see Deguchi Reference Deguchi2015 for details). Indeed, a close inspection of figure 3(b) shows that the red circle lies close to the scale at which the approximation begins to fail. Therefore, for safety, in figure 3(d) we vary $\beta$ from 8 to 6 fixing $\alpha _1/\beta ,\alpha _2/\beta$ and $\alpha _3/\beta$ . Inspecting the flow field of the solution marked by the blue diamond (i.e. $(\alpha _1,\alpha _2,\alpha _3)=(4.5,9,21)$ , the streamwise period becomes $2\pi /1.5$ ), we find that this final adjustment has the side benefit of producing nearly equispaced critical layers. Figure 5(a) shows a three-dimensional visualisation of this solution in the periodic domain, plotted in the same format as figure 2. Its implications for shear-flow turbulence are discussed in the next section.