Linear stability analysis of oblique Couette–Poiseuille flows

Abstract

We perform a detailed numerical study of modal and non-modal stability in oblique Couette–Poiseuille profiles, which are among the simplest examples of three-dimensional boundary layers. Through a comparison with the Orr–Sommerfeld operator for the aligned case, we show how an effective wall speed succinctly characterizes modal stability. Large-scale parameter sweeps reveal that the misalignment between the pressure gradient and wall motion is, in general, destabilizing. For flows that are sufficiently oblique, the instability is found to depend exclusively on the direction of wall motion and not on its speed, a conclusion supported, in part, by the perturbation energy budget and the evolution of the critical layers. Closed forms for the critical parameters in this regime are derived using a simple analysis. From a non-modal perspective, pseudoresonance is examined through the resolvent and the $\epsilon$-pseudospectra. An analysis of the unforced initial value problem shows that the maximum energy gain is highly dependent on both the magnitude and direction of the wall velocity. However, the strongest amplification is always achieved for configurations that are only weakly skewed. Finally, the optimal perturbations appear to develop via a lift-up effect enhanced by an Orr-like mechanism, the latter driven by cross-flow shear.

1. Introduction

Shear flows demonstrate rich dynamical behaviour and underpin a variety of technological applications, ranging from micro-fluidics and turbo-machinery to large-scale aerodynamics and meteorology. Quantifying the transition to turbulence in these flows is a multi-faceted problem, and despite much concentrated effort in the past few decades, a unified description is yet to be agreed upon. The state-of-the-art on the topic can be found in references such as Kerswell (2005), Manneville (2015), Barkley (2016), Eckhardt (2018) and Avila, Barkley & Hof (2023). Unfortunately, the insight afforded by prevailing theories is rather limited since a vast majority, primarily for the sake of simplification, focus on flows that are statistically two-dimensional, with a mean flow direction independent of the wall-normal coordinate. In contrast, most practical flow scenarios suffer from non-equilibrium effects that introduce anisotropy and symmetry-breaking, inducing a three-dimensional boundary layer described by a skewing of the mean velocity vectors and a flow direction that varies as a non-constant function along the wall-normal axis (Johnston & Flack 1996). In this work, we systematically examine the stability of a relatively under-explored class of three-dimensional internal flows that are both simple in their description and physically representative in their context.

Usually, the investigation of instability in a fluid system derives from the Navier–Stokes equations, appropriately linearized around some time-invariant state of interest. The classical (modal) approach focuses on the eigenvalues of the associated linear operator, and the flow is considered unstable with respect to some dimensionless quantity – usually a Reynolds number $Re$ – if there exists a mode with a positive growth rate. These disturbances are capable of achieving exponential-in-time amplification, a type of so-called primary instability, before saturating due to nonlinear phenomena. The result is either transition or the development of a steady, modified, base flow susceptible to secondary instabilities. To varying degrees of success, this framework has been applied to numerous laminar profiles, such as rectilinear pressure-driven flow (Thomas 1953; Gage & Reid 1968; Orszag 1971; Zhang et al. 2013), plane Couette flow (Drazin & Reid 2004; Zou et al. 2023), classic Rayleigh–Bénard convection (Rayleigh 1916; Chandrasekhar 1961), Hagen–Poiseuille flow (Salwen, Cotton & Grosch 1980; Schmid & Henningson 2001) and the asymptotic suction boundary layer (Hughes & Reid 1965; Fransson & Alfredsson 2003). Contextualizing these calculations against experimental observations, however, is not always straightforward, particularly because the most unstable disturbance, as predicted by modal theory, can only be realized in asymptotic time. On the other hand, significant non-modal energy growth can occur on much shorter time scales and is, therefore, not captured (Trefethen et al. 1993; Trefethen 1997; Schmid & Henningson 2001; Schmid 2007). A potential model for this behaviour lies in the non-normality of the linearized Navier–Stokes operator and its, in general, non-orthogonal eigenfunctions. In particular, within the basis expansion of an arbitrary initial disturbance, the contribution of these modes can grow or decay at different rates, allowing for a transient amplification of energy. In many canonical flows and their variants, this non-modal growth has been shown to be substantial, even in linearly stable parameter regimes (Farrell 1988; Reddy, Schmid & Henningson 1993; Schmid & Henningson 1994; Hristova et al. 2002; Meseguer 2002; Meseguer & Trefethen 2003; Liu & Liu 2012).

A flow that has arguably received limited attention in the general fluids literature is the plane Couette–Poiseuille (PCP) flow, generated by the interaction between a pressure gradient and the prescribed motion of one or both boundaries. The PCP configurations are standard in geophysical fluid mechanics (for example, in modelling asthenospheric counterflows; see Turcotte & Schubert 2002), flow in ducts (Owolabi, Dennis & Poole 2019) and tribology. Comprehensive stability analyses for PCP flows are somewhat sparse in number, and the first few treatments can be found in Potter (1966), Reynolds & Potter (1967) and Hains (1967). Described by a wall speed $\xi$ often made non-dimensional with the Poiseuille maximum, the superposition of a (parallel) Couette component with an otherwise pressure-driven flow is, in general, stabilizing, at least in terms of a critical Reynolds number $Re_c$ below which modal instability is absent. Furthermore, beyond a threshold value for $\xi$, the base flow is capable of achieving complete modal stability against infinitesimal perturbations, $Re_c\to \infty$. Cowley & Smith (1985), using a weakly nonlinear analysis, determined this ‘cutoff’ velocity to be $\xi \approx 0.7$. From a non-modal perspective, Bergström (2004) showed that the peak in transient energy amplification can depend heavily on the relative influence of the Poiseuille and Couette components. Modifications to the geometry and rheology of PCP flows have also been considered and their linear response analysed, for example, in Nouar & Frigaard (2009), Guha & Frigaard (2010), Chokshi et al. (2015) and Samanta (2020). More recently, using a zero-mean advection velocity set-up, Klotz et al. (2017) experimentally probed the subcritical transition in PCP flows. Uniform turbulence originating from a natural transition was observed beyond $Re\approx 780$ (with the Reynolds number based on the wall velocity), which stands in reasonable agreement with the results of Tsanis & Leutheusser (1988).

However, despite their individual merits, the previous literature on the transitional regime in PCP flows has remained somewhat restrictive in its applicability. For convenience in modelling or due to experimental limitations, the pressure gradient and wall velocity vectors are almost always assumed to be perfectly coincident. Although this unidirectional assumption enables a fairly tractable analysis, it breaks down in more practical scenarios, e.g. wind–ocean interactions, where the direction of the pressure gradient in the bulk flow need not be aligned with that of the wind shear. In these cases, a cross-flow must be taken into account, and the flow angle cannot be assumed to be constant, potentially affecting, among other aspects, the onset of instability. A linear modal and non-modal analysis of such three-dimensional PCP flows is yet to be performed, at least to the best of our knowledge. Therefore, the primary goal here is to contribute to this gap.

We note that linear instability in three-dimensional boundary layers has already been the subject of extensive prior investigation, with excellent reviews published in Reed & Saric (1989) and Saric, Reed & White (2003). In most flow situations, the introduction of mean three-dimensionality allows for multiple competing mechanisms for transition. For example, in flows over swept wings, local pressure gradients induce a near-wall cross-flow component that is inflectional and, therefore, unstable to the inviscid amplification of the so-called cross-flow modes (Gray 1952; Gregory, Stuart & Walker 1955; Bippes 1999). These modes are co-rotating and, following nonlinear saturation, incite breakdown to turbulence via high-frequency secondary instabilities (White & Saric 2005). Such cross-flow instabilities are also common in, for example, the boundary layers that form on rotating disks (Fedorov et al. 1976; Kobayashi, Kohama & Takamadate 1980; Malik 1986a). In particular, using von Kármán's similarity solution for a swirling flow over an infinitely long rotating disk, Malik, Wilkinson & Orszag (1981) determined the critical Reynolds number associated with these cross-flow modes to be $Re_c\approx 170$. Taking into account the streamline curvature and Coriolis effects, the authors refined this estimate in the same work to $Re_c\approx 290$, which was in better agreement with their experimental results. Later, Malik (1986b) and Balakumar & Malik (1990) found a second minimum in the neutral stability curve at $Re\approx 450$, which they associated with a viscous-type instability. Other classic studies on the transition of three-dimensional boundary layers include those of swept cylinders (Poll 1985; Kohama, Ukaku & Ohta 1988; Itoh 1996), rotating cones (Kobayashi 1981; Kobayashi & Izumi 1983; Kobayashi, Kohama & Kurosawa 1983) and the Ekman layer (Lilly 1966; Melander 1983; Spall & Wood 1984). On the topic of non-modal disturbances, Corbett & Bottaro (2001) found that swept boundary layers subjected to adverse or favourable pressure gradients were prone to larger transient growth than unswept cases. The authors also determined that, similar to two-dimensional flows, streamwise-elongated streaks comprised the most dangerous initial perturbation. More recently, Hack & Zaki (2014) probed the effects of time-harmonic spanwise wall motion on transitional boundary layers using direct numerical simulation. They observed that the wall motion could either suppress or accelerate transition. Using a frozen-in-phase as well as a Floquet approach, Hack & Zaki (2015) attributed this to the competition between the non-modal amplification of boundary layer streaks and the inviscid growth of inflectional modes introduced by the spanwise Stokes profile.

Interestingly, studies on skewed turbulent Couette–Poiseuille flows seem not to be uncommon, though almost every work so far has focused exclusively on wall motion that is purely orthogonal to the pressure gradient. In this community, such cases fall into the category of ‘shear-driven’ or ‘viscous-induced’ three-dimensional boundary layers. Coleman, Kim & Le (1996), Howard & Sandham (1997) and Le, Coleman & Kim (2000), for example, explored the variation in turbulent statistics of a two-dimensional channel flow after the sudden imposition of a spanwise wall motion. Kannepalli & Piomelli (2000) displaced only a finite section of the wall, focusing on the contrast between the initial response to the perturbation and the subsequent relaxation to a two-dimensional equilibrium turbulence. More recently, Holstad, Andersson & Pettersen (2010) investigated near-wall coherent structures in a turbulent Couette flow skewed by a spanwise pressure gradient. A common topic of emphasis within these studies seems to be the counterintuitive structural changes that occur between two-dimensional and three-dimensional (both equilibrium and non-equilibrium) boundary layers. For example, although the addition of mean shear in the two-dimensional case is known to enhance turbulence, it usually leads to a reduction in turbulent stresses in the three-dimensional setting. Furthermore, Townsend's structure parameter is also often quoted to decrease, indicating a decline in the efficiency of turbulent kinetic energy production by the mean flow. For relevant reviews on the topic, we direct the reader to Olcmen & Simpson (1993), Eaton (1995) and Johnston & Flack (1996). Given the intricate nature of three-dimensional boundary layers and our limited understanding of their physics, it is hoped that this work will supplement and invigorate ongoing investigations in this area.

We structure the remainder of the paper as follows. Section 2 introduces our base flow and develops our analysis frameworks. Sections 3 and 4 explore, respectively, modal and non-modal perturbations. Section 5 offers conclusions and a discussion of relevant future work.

2. Problem formulation

2.1. Governing equations and base profiles

We use the standard equations of motion for an incompressible Newtonian fluid. In the absence of body forces, these can be expressed in dimensional format as follows:

(2.1)\begin{equation} \rho\left[\frac{\partial\tilde{\boldsymbol{u}}}{\partial\tilde{t}} +(\tilde{\boldsymbol{u}}\boldsymbol{\cdot}\tilde{\boldsymbol{\nabla}})\tilde{\boldsymbol{u}}\right] ={-}\tilde{\boldsymbol{\nabla}}\tilde{p} + \mu\tilde{\nabla}^2\tilde{\boldsymbol{u}}, \end{equation}

where $\tilde {\boldsymbol {u}} = (\tilde {u}\ \tilde {v} \ \tilde {w})^\intercal$ is the Eulerian velocity field, $\tilde {p}$ the hydrodynamic pressure, $\rho$ the fluid density and $\mu$ the dynamic viscosity. The flow of interest in this study is illustrated in the schematic presented in figure 1. Two rigid surfaces, infinite in the wall-parallel directions and located at $\tilde {y} = \pm h$, confine an incompressible fluid subject to a fixed streamwise pressure gradient $\mathcal {G} < 0$. A cross-flow is established by additionally translating the top wall with a constant velocity $U_w$ at an angle $\theta$ with respect to the positive $\tilde {x}$ axis. The steady laminar profile $\tilde {\boldsymbol {U}}$ satisfies

(2.2a–c)\begin{equation} {-}\tilde{\boldsymbol{\nabla}} \tilde{P} + \mu\tilde{\nabla}^2\tilde{\boldsymbol{U}} = 0, \quad \tilde{\boldsymbol{U}} = \begin{pmatrix} \tilde{U}(\tilde{y}) & 0 & \tilde{W}(\tilde{y}) \end{pmatrix}^\intercal, \quad\tilde{\boldsymbol{\nabla}}\tilde{P} = \begin{pmatrix} \mathcal{G} & 0 & 0 \end{pmatrix}^\intercal, \end{equation}

subject to the boundary conditions

(2.3a–c)\begin{align} \tilde{U}(\tilde{y}=h) = U_w\cos\theta,\quad \tilde{W}(\tilde{y}=h) = U_w\sin\theta, \quad\tilde{U}(\tilde{y}={-}h) = \tilde{W}(\tilde{y}={-}h) = 0. \end{align}

In particular, we can find

(2.4a,b)\begin{align} \tilde{U}(\tilde{y}) ={-}\frac{h^2\mathcal{G}}{2\mu}\left(1 - \frac{\tilde{y}^2}{h^2}\right) + \frac{U_w}{2}\left(1+\frac{\tilde{y}}{h}\right)\cos\theta,\quad \tilde{W}(\tilde{y}) = \frac{U_w}{2}\left(1+\frac{\tilde{y}}{h}\right)\sin\theta. \end{align}

The resulting system is, therefore, a viscous-induced three-dimensional boundary layer, for which the flow angle, defined as

(2.5)\begin{equation} \phi(\tilde{y}) = \tan^{{-}1}\left(\frac{\tilde{W}(\tilde{y})}{\tilde{U} (\tilde{y})}\right), \end{equation}

varies with the wall-normal direction. These configurations are herein referred to as oblique Couette–Poiseuille flows (OCPfs) and, to our knowledge, have not received prior treatment in the stability literature, despite being among the simplest three-dimensional flows capable of retaining homogeneity in the streamwise and spanwise directions. Respectively, $\tilde {U}$ and $\tilde {W}$ are Couette–Poiseuille and Couette profiles, their relative strengths modulated by the direction of wall movement. In the limit $U_w\to 0$, standard Poiseuille flow is recovered. On the other hand, for $\theta \to 0$ and $U_w\neq 0$, the cross-flow vanishes and the system reduces to the well-known aligned Couette–Poiseuille flow (ACPf), in which the pressure gradient and wall motion coincide exactly.

Figure 1. A sketch of the three-dimensional flow geometry for oblique Couette–Poiseuille flows; here, $\mathrm {d}p/\mathrm {d}{\kern0.8pt}x = \mathcal {G} < 0$ is the constant streamwise pressure gradient. The wall at $\tilde {y} = h$ translates with velocity $U_w$ at an angle $\theta \neq 0$ to the streamwise direction, inducing a three-dimensional shear flow.

The parameter space characterizing OCPfs is rather complex, and, as is the case for ACPf, there exist multiple routes to rendering the governing equations non-dimensional. An obvious candidate is $U_p$, the so-called Poiseuille velocity scale, which is the streamwise maximum computed in the absence of wall motion. The other option is $U_{max}$, the ‘actual’ streamwise maximum, and is preferred if non-equilibrium effects are expected to significantly distort the streamwise profile away from $U_p$. However, in all possible realizations of OCPf, the boundedness of $\cos \theta$ and $\sin \theta$ ensures that $\tilde {U}$ is $O(U_p)$. Therefore, to facilitate comparison with the previous literature, we choose to scale with $U_p$. More specifically, the following non-dimensionalization scheme is adopted:

(2.6a–d)\begin{equation} \boldsymbol{x} = \frac{\tilde{\boldsymbol{x}}}{h},\quad \boldsymbol{u} = \frac{\tilde{\boldsymbol{u}}}{U_p},\quad t = \frac{\tilde{t}}{h/U_p},\quad p = \frac{\tilde{p}}{\rho U_p^2}, \end{equation}

which yields the dimensionless form of the momentum equations

(2.7)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u} ={-}\boldsymbol{\nabla} p + \frac{1}{Re}\nabla^2\boldsymbol{u}, \end{gather}$$

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

Here, (2.8) represents the incompressibility constraint, and $Re = \rho U_p h/\mu = U_ph/\nu$ is a Reynolds number, with $\nu$ being the kinematic viscosity. The base velocity profiles become

(2.9a,b)\begin{equation} U(y) = 1-y^2 + \frac{\xi}{2}(1+y)\cos\theta, \quad W(y) = \frac{\xi}{2}(1+y)\sin\theta, \end{equation}

where by defining $Re_w = U_wh/\nu$, we can interpret $\xi =U_w / U_p = Re_w/Re$ as the non-dimensional wall speed. In this setting, the influence of the shear angle on the base profiles becomes more apparent. Suppose that $\xi$ is fixed and $\theta$ is varied; while $W$ maintains its Couette nature, $U$ evolves continuously as a one-parameter homotopy between ACPf and the plane Poiseuille flow ($\theta = n{\rm \pi} /2$ for odd $n$). Therefore, it is reasonable to limit attention to pairs $(\xi, \theta )\in [0, 1]\times [0, 2{\rm \pi} ]$, the former due to its physical relevance and the latter due to the periodicity of the base profiles that can be expected to permeate the forthcoming calculations. For select values of the flow parameters, the associated non-dimensional profiles are offered in figure 2.

Figure 2. From left to right, plots of the streamwise and spanwise velocities as well as the flow direction $\phi$ (normalized by ${\rm \pi}$) against the wall-normal coordinate $y$: (a–c) $\theta = {\rm \pi}/4$ and $\xi \in \{0.2, 0.4, 0.6, 0.8, 1\}$; (d–f) $\xi =0.5$ and $\theta \in \{{\rm \pi} /8, {\rm \pi}/4, 3{\rm \pi} /8, {\rm \pi}/2\}$. Formally, $\phi$ is singular near the lower wall, where $U$ and $W$ both vanish due to the no-slip condition. However, from l'Hopital's rule, the limit can be computed as $\phi (y\to -1) = \tan ^{-1}(\xi \sin \theta /(4+\xi \cos \theta ))$, evidently the angle between the wall shear stresses.

2.2. The linearized system

This section follows standard monographs on hydrodynamic stability, and we refer the reader to the works of Schmid & Henningson (2001) or Drazin & Reid (2004), for example. In operator format, the Navier–Stokes equations can be rewritten as

(2.10)\begin{equation} \frac{\partial \boldsymbol{u}^*}{\partial t} = \mathcal{N}(\boldsymbol{u}^*), \end{equation}

where $\mathcal {N}$ is a nonlinear function of the state vector $\boldsymbol {u}^* = (\boldsymbol {u} \ p)^\intercal$. We decompose $\boldsymbol {u}^*$ as $\boldsymbol {u}^* = \boldsymbol {U}^* + {\boldsymbol {u}^\prime }^*$, where $\boldsymbol {U}^*$ is a time-independent base state superposed by a set of infinitesimal fluctuations ${\boldsymbol {u}^\prime }^* = (\boldsymbol {u}^\prime \ p^\prime )^\intercal$. In particular, we have

(2.11a–c)\begin{equation} \boldsymbol{U}^* = \begin{pmatrix} \boldsymbol{U} & P \end{pmatrix}^\intercal, \quad \boldsymbol{U} = \begin{pmatrix} U(y) & 0 & W(y) \end{pmatrix}^\intercal, \quad \boldsymbol{\nabla} P = \begin{pmatrix} -2/Re & 0 & 0 \end{pmatrix}^\intercal. \end{equation}

By Taylor expanding $\mathcal {N}$ around $\boldsymbol {U}^*$ and neglecting terms that are $O(\lVert {\boldsymbol {u}^\prime }^* \rVert ^2)$, we obtain a linearized system of evolution equations for the perturbation variables. To reduce computational complexity and the size of the matrices dealt with, the usual procedure here is to eliminate the pressure. This yields a rephrased system based only on fluctuations in the wall-normal velocity/vorticity $\boldsymbol {q} = (v^\prime \ \eta ^\prime )$:

(2.12)$$\begin{gather} \left[\left(\frac{\partial}{\partial t} + U\frac{\partial}{\partial x} + W\frac{\partial}{\partial z}\right)\nabla^2 - \frac{\mathrm{d}^2 U}{\mathrm{d} y^2}\frac{\partial}{\partial x}- \frac{\mathrm{d}^2 W}{\mathrm{d} y^2}\frac{\partial}{\partial z}- \frac{1}{Re}\nabla^4\right]v' = 0, \end{gather}$$

(2.13)$$\begin{gather}\left[\frac{\partial}{\partial t} + U\frac{\partial}{\partial x} + W\frac{\partial}{\partial z} - \frac{1} {Re}\nabla^2\right]\eta'-\frac{\mathrm{d}W}{\mathrm{d} y}\frac{\partial v'}{\partial x}+\frac{\mathrm{d}U}{\mathrm{d} y}\frac{\partial v'}{\partial z}= 0, \end{gather}$$

where $\nabla ^2$ is the usual Laplacian in a Cartesian coordinate system and $\nabla ^4\langle \cdot \rangle \equiv \nabla ^2(\nabla ^2\langle \cdot \rangle )$ is the bi-harmonic operator. Here on, for notational brevity, we drop the prime notation. Note that, contrary to the case of a purely streamwise base flow for which $W = 0$, the so-called Squire equation, (2.13), is now forced by mean shear from both the streamwise and spanwise profiles, which are, in general, non-zero. The spatial homogeneity can be exploited via a Fourier transform,

(2.14)\begin{equation} \bar{\boldsymbol{q}}(y,t; \alpha, \beta) = \iint_{-\infty}^\infty\boldsymbol{q}(x,y,z,t)\, {\rm e}^{-{\rm i}(\alpha x+\beta z)}\,\mathrm{d}{\kern0.8pt}x\,\mathrm{d}z, \end{equation}

to obtain the canonical form of the Orr–Sommerfeld–Squire (OSS) system. Here, $\alpha,\beta \in \mathbb {R}$ are the real-valued wavenumbers in the $x$ and $z$ directions and $\bar {\boldsymbol {q}} = (\bar {v} \ \bar {\eta })$ is a block vector of Fourier coefficients. The transformed equations can be compactly written as

(2.15)\begin{equation} {{\boldsymbol{\mathsf{L}}}} \bar{\boldsymbol{q}} ={-}\frac{\partial}{\partial t}{{\boldsymbol{\mathsf{M}}}} \bar{\boldsymbol{q}}, \end{equation}

where, by denoting $\mathcal {D}\equiv \mathrm {d}/\mathrm {d} y$ and $k^2 = \alpha ^2+\beta ^2$, we have defined

(2.16a,b)\begin{equation} {{\boldsymbol{\mathsf{L}}}} = \begin{pmatrix} \mathcal{L}_{OS} & 0 \\ {\rm i} \beta \mathcal{D}U - {\rm i} \alpha\mathcal{D}W & \mathcal{L}_{SQ} \end{pmatrix},\quad {{\boldsymbol{\mathsf{M}}}} = \begin{pmatrix} \mathcal{D}^2 - k^2 & 0 \\ 0 & 1 \end{pmatrix}. \end{equation}

The Orr–Sommerfeld (OS) and Squire operators, $\mathcal {L}_{OS}$ and $\mathcal {L}_{SQ}$ respectively, are given by

(2.17)$$\begin{gather} \mathcal{L}_{OS} = ({\rm i} \alpha U+{\rm i} \beta W)(\mathcal{D}^2-k^2) - {\rm i} \alpha \mathcal{D}^2 U-{\rm i} \beta \mathcal{D}^2 W-\frac{1}{Re}(\mathcal{D}^2-k^2)^2, \end{gather}$$

(2.18)$$\begin{gather}\mathcal{L}_{SQ} = {\rm i} \alpha U +{\rm i} \beta W -\frac{1}{ Re}(\mathcal{D}^2-k^2). \end{gather}$$

Equation (2.15) forms an initial-value problem for the Fourier-transformed state vector $\bar {\boldsymbol {q}}$ in wavenumber space, where the associated boundary conditions can be obtained by applying no-slip/impermeability at both walls. Whenever necessary, the velocity–vorticity formulation of the OSS problem can be recast into one for the primitive fluctuations using the transformation

(2.19)\begin{equation} \begin{pmatrix} \bar{u} \\ \bar{v} \\ \bar{w} \end{pmatrix} = \frac{1}{k^2}\begin{pmatrix} {\rm i} \alpha\mathcal{D} & -{\rm i} \beta\\ k^2 & 0\\ {\rm i} \beta\mathcal{D} & {\rm i} \alpha \end{pmatrix}\begin{pmatrix} \bar{v} \\ \bar{\eta} \end{pmatrix}. \end{equation}

Details of the numerical discretization of (2.15) can be found in Appendix A.

2.3. Modal analysis

For a modal or eigenvalue analysis, an additional Fourier transform is conducted in time:

(2.20)\begin{equation} \hat{\boldsymbol{q}}(y; \alpha, \beta, \omega) = \int_{-\infty}^\infty \bar{\boldsymbol{q}}(y, t;\alpha, \beta)\, {\rm e}^{{\rm i} \omega t}\,\mathrm{d}t, \end{equation}

where $\omega = \omega _r + {\rm i} \omega _i \in \mathbb {C}$ is the complex wave frequency. Equation (2.15) then reduces to a generalized eigenvalue problem described by the linear operator pencil $({{\boldsymbol{\mathsf{L}}}},{{\boldsymbol{\mathsf{M}}}})$:

(2.21)\begin{equation} {{\boldsymbol{\mathsf{L}}}}\hat{\boldsymbol{q}} = {\rm i} \omega{{\boldsymbol{\mathsf{M}}}}\hat{\boldsymbol{q}}, \end{equation}

with eigenvalues corresponding to ${\rm i} \omega = {\rm i} \omega _r - \omega _i$. Note that this is equivalent to solving for the eigensystem of ${{\boldsymbol{\mathsf{S}}}}^\prime = {{\boldsymbol{\mathsf{M}}}}^{-1}{{\boldsymbol{\mathsf{L}}}}$. In general, the spectrum is a function of $\{\alpha, \beta, Re, \xi, \theta \}$, and exponential amplification occurs over time if $\omega _i > 0$. Consequently, we seek the manifold of marginal stability, designated by

(2.22)\begin{equation} \omega_i(\alpha, \beta, Re, \xi, \theta) = 0. \end{equation}

We note that the presence of a non-zero spanwise velocity in OCPfs prevents an application of Squire's theorem in its usual form. Although a two-dimensional problem may well be constructed (see e.g. Mack 1984; Schmid & Henningson 2001), the ‘effective’ base velocity depends on both spatial wavenumbers and there is no a priori indication of the appropriate search space. Therefore, for a given configuration $(\xi, \theta )$, since a full stability portrait requires a sweep through the $(\alpha, \beta, Re)$-space, a numerical approach will inevitably be marred by a lack of resolution. While this is a valid criticism, we point out that most canonical shear flows only become linearly unstable at modest wavenumbers, if at all. Furthermore, in § 3, we demonstrate that from the perspective of modal stability, OCPfs are essentially continuations of the aligned variant. Therefore, the results of a sufficiently broad numerical search, as conducted here, are likely global.

Before proceeding, we make some key observations. First, as is true for strictly streamwise base flows, the Squire modes remain damped. The proof proceeds in the usual way by converting to a formulation involving the $x$-phase speed, $c=\omega /\alpha$, multiplying the homogeneous Squire equation by the complex conjugate of the fluctuating normal vorticity and integrating over $y$. Therefore, for a modal analysis, it suffices to consider only the OS operator, (2.17). Furthermore, since neither component of the base velocity is inflectional, OCPfs do not admit an inviscid cross-flow-like instability as observed, for example, over swept wings or rotating disks. In particular, in the inviscid limit, Rayleigh's criterion can be modified to require the following expression to hold at some wall-normal location:

(2.23)\begin{equation} \mathcal{D}^2U + \gamma\mathcal{D}^2W = 0, \end{equation}

where $\gamma = \beta /\alpha$. Although (2.23) will, for general flows, vary in wavenumber space, the linearity of $W$ implies that $\mathcal {D}^2W = 0$ for OCPfs. Thus, since $\mathcal {D}^2U = -2$, the instability must be viscous in nature.

2.4. Non-modal analysis

For most shear flows, a spectral analysis of the linearized Jacobian as in § 2.3 rarely agrees with experiment. In fact, quite often, the transition to turbulence is observed at subcritical $Re$, that is, below the threshold predicted by modal theory (Trefethen et al. 1993). This behaviour is now well understood to be a consequence of the highly non-normal nature of the OSS operator ${{\boldsymbol{\mathsf{S}}}}^\prime$, which, in part, arises from the off-diagonal term $({\rm i} \beta \mathcal {D} U - {\rm i} \alpha \mathcal {D}W)$ driving the Squire equation; see (2.13) and (2.16a,b). In general OCPfs, this forcing can evidently comprise both the streamwise and spanwise mean shear.

A non-normal operator such as ${{\boldsymbol{\mathsf{S}}}}^\prime$ admits eigenfunctions that are non-orthogonal in the underlying Hilbert space. When arbitrary initial states are transformed into the basis of these eigenfunctions, they can suffer from large cross-terms in the induced norm (Schmid 2007). An immediate consequence is that while a modal analysis might suggest asymptotic decay, energy amplification can still occur over finite time horizons. In shear flows, the transition to turbulent regimes has often been attributed to these transient phenomena, providing a potential explanation for the so-called bypass transition (Butler & Farrell 1992). Furthermore, there is no guarantee that the long-time eigenmode is even realized, in spite of the most careful calibration, since sufficiently strong transient amplification will likely excite nonlinear mechanisms in the flow and violate the linear assumption (Waleffe 1995b; Trefethen 1997)

To explore the implications of non-normality in OCPfs, we first solve the initial-value problem in (2.15) exactly to yield

(2.24)\begin{equation} \bar{\boldsymbol{q}}(t) = \varPhi(t, 0)\bar{\boldsymbol{q}}_0, \end{equation}

where $\varPhi (t,0)\equiv {\rm e}^{{\rm i} {{\boldsymbol{\mathsf{S}}}} t}$ is the state-transition operator, ${{\boldsymbol{\mathsf{S}}}} = i {{\boldsymbol{\mathsf{S}}}}^\prime$, and $\bar {\boldsymbol {q}}_0$ is the state of the system at the initial time $t=0$. Under appropriate norms in the input and output spaces, the gain can be defined as

(2.25)\begin{equation} G(\alpha, \beta, Re, \xi, \theta, t) = \sup_{\bar{\boldsymbol{q}}_0\neq 0}\frac{\lVert \bar{\boldsymbol{q}}\rVert_{{out}}^2}{\lVert \bar{\boldsymbol{q}}_0\rVert_{in}^2}, \end{equation}

where, due to its physical significance, we let $\lVert \cdot \rVert _{out}=\lVert \cdot \rVert _{in}=\lVert \cdot \rVert _E$ be an energy norm,

(2.26)\begin{equation} \lVert \bar{\boldsymbol{q}}\rVert_E^2 = \int_{{-}1}^1\bar{v}^{{{\dagger}}}\bar{v} + \frac{1}{k^2}\left(\bar{\eta}^{{{\dagger}}}\bar{\eta} + \frac{\partial\bar{v}^{{{\dagger}}}}{\partial y}\frac{\partial\bar{v}}{\partial y}\right)\,\mathrm{d} y\simeq \bar{\boldsymbol{q}}^{\dagger}{{\boldsymbol{\mathsf{E}}}}\bar{\boldsymbol{q}}, \end{equation}

over the volume $V$ defined by the Cartesian product $(x,y,z)\in [0,2{\rm \pi} /\alpha ]\times [-1,1]\times [0,2{\rm \pi} /\beta ]$. In this way, the energy of one full wavelength of a disturbance can be captured (see Butler & Farrell 1992). Here, $\langle \cdot \rangle ^{\dagger}$ denotes a conjugate transpose operation, and the operator ${{\boldsymbol{\mathsf{E}}}}$ is positive-definite and incorporates the Clenshaw–Curtis quadrature weights (Trefethen 2000). With a Cholesky decomposition, we may write ${{\boldsymbol{\mathsf{E}}}} = {{\boldsymbol{\mathsf{F}}}}^{\dagger} {{\boldsymbol{\mathsf{F}}}}$ so that

(2.27)\begin{equation} \lVert \bar{\boldsymbol{q}}\rVert_E^2 \simeq \bar{\boldsymbol{q}}^{\dagger}{{\boldsymbol{\mathsf{F}}}}^{\dagger}{{\boldsymbol{\mathsf{F}}}}\bar{\boldsymbol{q}} = \lVert {{\boldsymbol{\mathsf{F}}}} \bar{\boldsymbol{q}}\rVert_2^2. \end{equation}

It immediately follows that

(2.28)\begin{equation} G = \sup_{\bar{\boldsymbol{q}}_0\neq 0}\frac{\lVert{{\boldsymbol{\mathsf{F}}}}\varPhi(t,0) \bar{\boldsymbol{q}}_0\rVert_2^2}{\lVert{{\boldsymbol{\mathsf{F}}}} \bar{\boldsymbol{q}}_0\rVert_2^2} = \sup_{\bar{\boldsymbol{q}}_0\neq 0}\frac{\lVert{{\boldsymbol{\mathsf{F}}}}\varPhi(t,0) {{\boldsymbol{\mathsf{F}}}}^{{-}1}{{\boldsymbol{\mathsf{F}}}} \bar{\boldsymbol{q}}_0\rVert_2^2}{\lVert{{\boldsymbol{\mathsf{F}}}} \bar{\boldsymbol{q}}_0\rVert_2^2} =\lVert{{\boldsymbol{\mathsf{F}}}}\varPhi(t,0){{\boldsymbol{\mathsf{F}}}}^{{-}1}\rVert_2^2, \end{equation}

which can be computed trivially via the singular value decomposition (note, in fact, that $G = \lVert \varPhi (t,0)\rVert _E^2$). The associated right and left singular functions represent, respectively, the initial condition and response pair for which the gain at time $t$ is realized.

Intuitively, no energy growth is expected if $G\leq 1$. An equivalent condition can be expressed in terms of the resolvent of ${{\boldsymbol{\mathsf{S}}}}$. Consider an exogenous harmonic forcing profile $\boldsymbol {H}(x,y,z,t) = \boldsymbol {h}(x,y,z)\, {\rm e}^{-{\rm i}\zeta t}$ with frequency $\zeta \in \mathbb {C}$ to the linearized system, appropriately transformed into wavenumber space:

(2.29)\begin{equation} \bar{\boldsymbol{h}}(y; \alpha, \beta) = \iint_{-\infty}^\infty\boldsymbol{h}(x,y,z)\, {\rm e}^{-{\rm i}(\alpha x+\beta z)}\,\mathrm{d}{\kern0.8pt}x\,\mathrm{d}z. \end{equation}

The long-time response, assuming asymptotic stability, can easily be verified to be

(2.30)\begin{equation} \bar{\boldsymbol{q}} = {\rm i}{\rm e}^{-{\rm i}\zeta t}(\zeta {{\boldsymbol{\mathsf{I}}}} - {{\boldsymbol{\mathsf{S}}}})^{{-}1}\bar{\boldsymbol{h}}, \end{equation}

where the operator ${{\boldsymbol{\mathsf{R}}}} \equiv (\zeta {{\boldsymbol{\mathsf{I}}}}-{{\boldsymbol{\mathsf{S}}}})^{-1}$ is known as the resolvent. From an input–output perspective, ${{\boldsymbol{\mathsf{R}}}}$ serves as a transfer function between the excitation and its response. The quantity $\mathcal {R} = \lVert {{\boldsymbol{\mathsf{R}}}}\rVert _E$ is, therefore, of particular interest here, since for a non-normal system, it can be large even if the forcing is pseudoresonant, that is, $\zeta \notin \varLambda ({{\boldsymbol{\mathsf{S}}}})$, the spectrum of ${{\boldsymbol{\mathsf{S}}}}$ (Trefethen & Embree 2005). Such a paradigm is especially informative for the receptivity of the flow to external disturbances (Brandt 2014), and if $\zeta$ is restricted to real values, a physical interpretation of the resolvent is the perturbed operator that can result, for example, from external vibrations or planar imperfections (Trefethen et al. 1993). By further generalizing to the complex plane, one recovers the $\epsilon$-pseudospectra, the set of values defined as

(2.31)\begin{equation} \varLambda_\epsilon({{\boldsymbol{\mathsf{S}}}}) = \{\zeta\in\mathbb{C}\colon\mathcal{R}\geq \epsilon^{{-}1}\}. \end{equation}

For non-normal operators, $\varLambda _\epsilon$ can protrude deep into the upper-half of the complex plane, and the more pronounced this effect, the greater the potential for transient growth irrespective of the presence of linear instability. More rigorously, the Hille–Yosida theorem states that $G\leq 1$ if and only if the $\epsilon$-pseudospectra lie sufficiently close to the lower half-plane (Reddy et al. 1993). For further details, we refer the reader to that paper, the citations within and the text of Trefethen & Embree (2005).

2.5. Energy budget analysis

An investigation of the perturbation energy budget can reveal the mechanism of instability in OCPfs. Throughout this section, the Einstein convention is implied via repeated indices. We define the perturbation energy density $\mathcal {E}$ as

(2.32)\begin{equation} \mathcal{E} = \tfrac{1}{2}\boldsymbol{u}^{\dagger}\boldsymbol{u} = \tfrac{1}{2}u_i^{\dagger} u_i = \tfrac{1}{2}(|u|^2 + |v|^2 + |w|^2). \end{equation}

By multiplying equation (2.7) throughout by $\boldsymbol {u}^{\dagger}$ and integrating over $V$, evolution equations for the total energy are recovered:

(2.33)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\int_V\mathcal{E}\,\mathrm{d}V =\int_V \frac{\mathrm{d}\mathcal{E}}{\mathrm{d}t}\,\mathrm{d}V={-}\int_V \frac{1}{2}(u_i^{\dagger} u_j+u_i u_j^{\dagger})\frac{\partial U_i}{\partial x_j}\,\mathrm{d}V -\frac{1}{Re}\int_V\frac{\partial u_i^{\dagger}}{\partial x_j}\frac{\partial u_i}{\partial x_j}\,\mathrm{d}V, \end{equation}

where we have assumed spatial periodicity of the disturbance field in $x$ and $z$. Under the normal mode ansatz, (2.20), the above expression reduces to

(2.34)\begin{equation} 2\omega_i\int_{{-}1}^1\frac{1}{2}\hat{u}_i^{\dagger}\hat{u}_i\,\mathrm{d} y = \int_{{-}1}^1\mathcal{P}\,\mathrm{d} y-\int_{{-}1}^1\varepsilon\,\mathrm{d} y, \end{equation}

where we have defined

(2.35)$$\begin{gather} \mathcal{P} =\underbrace{-\overbrace{\frac{1}{2}(\hat{u}^{\dagger}\hat{v} + \hat{u}\hat{v}^{\dagger})}^{\tau_u}\frac{\partial U}{\partial y}}_{\mathcal{P}_{u}}\underbrace{-\overbrace{\frac{1}{2}(\hat{w}^{\dagger}\hat{v} + \hat{w}\hat{v}^{\dagger})}^{\tau_w}\frac{\partial W}{\partial y}}_{\mathcal{P}_{w}}, \end{gather}$$

(2.36)$$\begin{gather}\varepsilon =\frac{1}{Re}[(\mathcal{D}\hat{u}_i)^{\dagger}\mathcal{D}\hat{u}_i + k^2 \hat{u}_i^{\dagger}\hat{u}_i]. \end{gather}$$

Two contributions to the disturbance kinetic energy can be identified: $\mathcal {P}$, the production against the background shear(s), and $\varepsilon$, the viscous dissipation. The former can be further separated into terms representing the transfer of energy from the base streamwise and spanwise flows, respectively, to the perturbation field through the action of the associated Reynolds stresses, $\tau _{u}$ and $\tau _{w}$. These have been denoted by $\mathcal {P}_{u}$ and $\mathcal {P}_{w}$. In general, (positive) production destabilizes, whereas dissipation stabilizes the disturbance field.

3. Modal analysis

3.1. Characteristics of the eigenspectra

We begin by investigating the dynamics of the eigenspectra in OCPfs. For a sample wavenumber combination, figure 3 illustrates the loci of the first $\approx 50$ least stable modes as the non-dimensional wall speed $\xi$ is varied at $\theta = {\rm \pi}/6$. The results have been presented in terms of $c = \omega /\alpha$. A familiar Y-shaped distribution can be observed, with three distinct branches reminiscent of the spectrum for plane Poiseuille flow (pPf). As the wall speed increases, this structure collectively translates further into the right half-plane, and the most unstable mode monotonically stabilizes. In doing so, the shape of the $S$-branch, comprising the so-called mean modes related to the mean velocity, remains relatively undistorted. On the other hand, a sharper change occurs in the $A$-branch – the wall modes – which separate into two distinct subsets associated, respectively, with each wall. In a somewhat similar manner, starting from its bottom half, the $P$-branch of centre modes also begins to split into two noticeable sub-branches. Together, these observations are indicative of the increased Couette contribution to the base flow, since the spectra for various flavours of Couette flow are usually scattered symmetrically within two $A$-branches (e.g. Duck, Erlebacher & Hussaini 1994; Schmid & Henningson 2001; Zou et al. 2023).

Figure 3. The locus of the eigenspectrum for $(\alpha,\beta ) = (1, 0.5)$ at $Re = 5700$ and $\theta = {\rm \pi}/6$ for (a–f) $\xi \in \{0, 0.2, 0.4, 0.6, 0.8, 1\}$. The $A$, $P$ and $S$ branches have been appropriately labelled. On each plot, a grey dashed line denotes the stability boundary, $c_i = 0$.

In OCPfs, the distribution of this Couette component between the base velocities is directly controlled by the shear angle $\theta$. However, its impact at the level of the OS equation is rather subtle. Since $W$ is linear, $\mathcal {D}^2W = 0$, and the OS operator, simplified from (2.17), becomes

(3.1)\begin{equation} \mathcal{L}_{OS} = \overbrace{({\rm i} \alpha U+{\rm i} \beta W)(\mathcal{D}^2-k^2)}^{\mathcal{O}_1} - \overbrace{{\rm i} \alpha \mathcal{D}^2 U}^{\mathcal{O}_2}-\overbrace{\frac{1}{Re}(\mathcal{D}^2 - k^2)^2}^{\mathcal{O}_3}, \end{equation}

where $U$ and $W$ retain their definitions from (2.9a,b), instantiated with some wall speed $\xi$. We immediately observe, despite the three-dimensionality of the flow, that the spanwise velocity appears only in a single term, $\mathcal {O}_1$, in (3.1). In particular, for a spanwise-independent mode, $\beta =0$, the effects of obliqueness in the base flow are, in a sense, ‘shut off’, since the corresponding OS operator

(3.2)\begin{equation} \mathcal{L}_{OS} = {\rm i} \alpha U(\mathcal{D}^2-\alpha^2) - {\rm i} \alpha \mathcal{D}^2 U-\frac{1}{Re}(\mathcal{D}^2 - \alpha^2)^2 \end{equation}

reduces precisely to that for ACPf under the umbrella of Squire's theorem (excluding, of course, the factor of $\cos \theta$ in $U$, which can essentially be lumped into the wall speed). To extend this analogy to more general disturbances, a modification must first be made. Consider the generic three-dimensional (that is, prior to an application of Squire's result) OS operator for ACPf

(3.3)\begin{equation} \mathcal{L}_{OS}^{ACPf} = \overbrace{{\rm i} \alpha U_{{ACPf}}(\mathcal{D}^2-k^2)}^{\mathcal{A}_1} - \overbrace{{\rm i} \alpha \mathcal{D}^2 U_{{ACPf}}}^{\mathcal{A}_2}-\overbrace{\frac{1}{Re}(\mathcal{D}^2 - k^2)^2}^{\mathcal{A}_3}, \end{equation}

where

(3.4)\begin{equation} U_{ACPf} = 1-y^2 + \frac{\xi_{ACPf}}{2}(1+y). \end{equation}

Comparing the two operators in (3.1) and (3.3) allows us to identify crucial similarities in structure. Specifically, for a constant wave triplet $(\alpha,\beta, Re)$, while $\mathcal {O}_3\equiv \mathcal {A}_3$ is immediate, $\mathcal {O}_2\equiv \mathcal {A}_2$ follows from the fact that $\mathcal {D}^2U=-2=\mathcal {D}^2 U_{{ACPf}}$. Therefore, $\mathcal {L}_{OS}$ and $\mathcal {L}_{OS}^{ACPf}$ differ exclusively in their terms $\mathcal {O}_1$ and $\mathcal {A}_1$, respectively. However, since $k^2$ has also been fixed by our choice of wavenumbers, $\mathcal {O}_1 \equiv \mathcal {A}_1$ can be made possible by requiring

(3.5)\begin{equation} {\rm i} \alpha U_{ACPf} = {\rm i} \alpha U + {\rm i} \beta W\implies \xi_{ACPf} = \xi(\cos\theta + \gamma\sin\theta), \end{equation}

where $\gamma = \beta /\alpha$. Therefore, for an arbitrary OCPf, the OS problem at any wavenumber pair can be exactly mapped to one for ACPf via the ‘effective’ wall speed $\xi _{{eff}}$:

(3.6)\begin{equation} \xi_{{eff}} \equiv \xi(\cos\theta + \gamma\sin\theta). \end{equation}

With the corollary

(3.7)\begin{equation} \omega(\alpha,\beta,Re,\xi,\theta\neq0) = \omega(\alpha,\beta,Re,\xi_{{eff}},\theta=0), \end{equation}

we conclude that the stability of any OCPf can be prescribed entirely by comparison with the appropriate ACPf configuration(s). A stronger result, and one perhaps in the same spirit as Squire's theorem, is as follows: if $\mathfrak {O}$ denotes the set of all possible OS operators for OCPf and $\mathfrak {A}$ the equivalent set for ACPf, then $\mathfrak {O}\subseteq \mathfrak {A}$.

The influence of the shear angle on modal behaviour can now be made precise. We start by noting that $\xi _{eff}$ is $2{\rm \pi}$-periodic and

(3.8)\begin{equation} {-}\xi k/\alpha \leq \xi_{eff} \leq \xi k/\alpha, \end{equation}

so that it varies strongly even throughout wavenumber space. Mathematically, at a fixed triplet $(\alpha, \beta, Re)$, its action in $\theta$ seems to be to accentuate or mask the strength of the wall speed. As an example, figure 4 shows how changes in $\xi _{eff}$ with $\theta$ affect the real and imaginary components of the most unstable eigenvalue for an arbitrarily chosen wavenumber pair. It is evident that the periodicity of $\xi _{eff}$ directly translates to that of the spectrum, which itself becomes, at a minimum, $2{\rm \pi}$-periodic. Furthermore, we found (not shown here; refer to figure 3) that variations in $\xi _{eff}$ modified the distribution of the eigenmodes in the complex plane in much the same fashion as variations in $\xi$ for fixed $\theta$, e.g. increasing $\xi _{eff}$ increased $c_r$, and vice versa. When taken together, these observations, combined with (3.7) and the interpretation of $\xi _{eff}$, suggest that the manifold of marginal stability for OCPfs is contained wholly within that for ACPf. Accounting for a non-trivial directionality in the flow affects perhaps only the subset of the latter that is ultimately accessed.

Figure 4. At $Re=10\ 000$ and $\xi =0.2$, the variation with $\xi _{eff}$ of the least stable eigenmode for $(\alpha, \beta ) = (1,0.25)$. The only dashed grey line marks the boundary $c_i=0$. Both components change in tandem with $\xi _{eff}$, and when juxtaposed with the information in figure 3, lend weight to $\xi _{eff}$ serving as an effective wall speed. Note that $c_i$ is, in fact, ${\rm \pi}$-periodic, the underlying mechanism being precisely that which allows for symmetric growth rates around $\xi = 0$ for ACPf (see § 3.2).

3.2. Exploring criticality in OCPfs

In this section, we present the findings of a comprehensive investigation into the modal stability of OCPfs. We introduce the critical Reynolds number, denoted $Re_c$, which represents the minimum Reynolds number below which the flow remains linearly stable. At this value, at least one disturbance, characterized by the critical wavenumbers $(\alpha _c,\beta _c)$, must achieve neutral stability. When analysing two-dimensional flows, Squire's theorem (Squire 1933) allows us to focus solely on disturbances that are independent of the spanwise direction, that is, $\beta _c = 0$. However, for general three-dimensional profiles, an accurate assessment of stability necessitates the consideration of modes with non-zero $\beta$. Consequently, in the case of an OCPf $(\xi,\theta )$, a thorough exploration of the entire three-dimensional $(\alpha, \beta, Re)$-space is required.

To reduce the degree of computation, we now consider some important simplifications. First, we note that in the stability literature for ACPf, the analysis for $\xi <0$ is typically neglected, since the modal growth rates are symmetric around $\xi = 0$ (although the corresponding real parts might not be). Potter (1966) rationalized this by adopting the coordinate transformation $y\to -y$. Since $\xi _{eff}(\theta + {\rm \pi}) = -\xi _{eff}(\theta )$, a similar argument allows us to restrict our attention to $\theta \in [0, {\rm \pi}]$. However, a second reduction is also possible and can be achieved by noting that, at a constant $\xi$, if $(\alpha _c, \beta _c, Re_c)$ is the critical tuple for $\theta =\theta ^\prime$, then $(\alpha _c, -\beta _c, Re_c)$ is necessarily the critical tuple for $\theta = {\rm \pi}- \theta ^\prime$. This result is immediate from the definition of $\xi _{eff}$ in (3.6), since

(3.9)\begin{equation} \xi_{{eff}}(\alpha, \beta, \xi, \theta^\prime) ={-}\xi_{{eff}}(\alpha, -\beta, \xi, {\rm \pi}-\theta^\prime), \end{equation}

where we have assumed $\alpha > 0$. Thus, it suffices to explore the range $\theta \in [0, {\rm \pi}/2]$. For the Fourier wavenumbers, we focused on small to intermediate values, in particular, $(\alpha, \beta ) \in [-3, 3]\times [-3, 3]$. This is generally the subspace of the wavenumber plane within which linear instability is first encountered in most canonical shear flows, and particularly for ACPfs (Potter 1966). In total, $O(10^{10})$ different parameter combinations were investigated, and our numerical procedure, including our method for traversing such an unwieldy space, is outlined in Appendix A. In what follows, all results are presented for values of $\theta$ in degrees rather than in radians.

In general, the introduction of skewness in Couette–Poiseuille flows was found to be destabilizing, at least relative to ACPf. However, two qualitative regimes could still be identified in $\theta$. The first, denoted $\varTheta _1$, comprises $0^\circ <\theta \lessapprox 20^\circ$ and is arguably the most interesting of the two, as it exhibits drastic changes in stability throughout its extent. Since OCPfs reduce to the standard aligned case as $\theta \to 0$, it is natural to expect the stability characteristics of ACPf to continue at least to modest $\theta$. Figure 5 supports this intuition. For all $\theta \in \varTheta _1$, a short range of stabilization in the $Re_c$ curves is followed by an inflection point between $0.2\lessapprox \xi \lessapprox 0.4$ and then further growth, a trend that is precisely reminiscent of ACPf (Potter 1966).

Figure 5. The critical Reynolds number $Re_c$ against $\xi$ for $\varTheta _1\equiv (0, 20^\circ ]$. Throughout this figure, a dashed line indicates the equivalent plot for ACPf. The insets magnify regions of particular interest that have been discussed in the text. A circle in inset (b) denotes the crossing point $\mathcal {I}$. In this range of shear angles, a typical $Re_c$ curve mimics that for ACPf when $\xi \lessapprox \xi _A$, but appears to have been ‘dragged’ down from infinity when $\xi > \xi _A$, yielding a finite $Re_c$ even beyond this threshold wall speed.

Perhaps the most striking feature is the fact that linear instability seems to persist throughout the entire range of wall speeds considered here. This behaviour was observed even for ‘small’ angles, such as $\theta \in \{1^\circ, 1.5^\circ \}$ (and even down to $\theta \in \{0.5^\circ, 0.75^\circ \}$, not shown here), where the wall motion is approximately parallel to the pressure gradient. This is in stark contrast to ACPf, which achieves unconditional linear stability, $Re_c\to \infty$, against infinitesimal disturbances beyond the so-called cutoff wall speed $\xi _A\approx 0.7$ (Potter 1966; Cowley & Smith 1985). A crude explanation for this is that the inclusion of a spanwise velocity makes $\beta$ a relevant stability parameter, providing, in light of the effective wall speed $\xi _{{eff}}$ and the analysis outlined at the end of § 3.1, an additional buffer for an OCPf to return to a region of linear instability for ACPfs. The absence of an equivalent cutoff wall speed for OCPfs seems to manifest itself in terms of the appearance of a limiting regime in the critical parameters. Here, linear instability seems to become entirely independent of $\xi$, as evidenced, in part, by the flattening of the $Re_c$ curves in figure 5.

We note that before achieving the respective asymptotes in their $Re_c$ curves, the destabilization experienced by OCPfs in $\varTheta _1$ at any wall speed is not necessarily monotonic with $\theta$. Beginning with inset (a) in figure 5, we see that increasing the shear angle is conclusively destabilizing up to $\xi \approx 0.2$. However, between approximately $0.2<\xi \lessapprox 0.275$, some of the more modest angles, say $\theta \gtrapprox 10^\circ$, tend to stabilize, while the values of $Re_c$ for even larger angles, $\theta \gtrapprox 18.5^\circ$, remain below those of ACPf (see figure 5b). Around $\xi \approx 0.3$ lies the crossing point $\mathcal {I}$, which initiates a region of monotonic stabilization for $\theta \gtrapprox 10^\circ$, a pattern that persists until around $\xi \approx 0.375$. Here, as seen in figure 5(c), all $Re_c$ curves experience a turning point, which occurs at increasing wall speeds with $\theta$, and begin to ascend towards their eventual plateaus. Beyond this location, the stabilization becomes monotonic throughout $\theta >0$, as can be verified in figure 5(d).

Figure 6 presents the spatial wavenumbers at criticality for $\theta \in \varTheta _1$. Consistent with the trends observed for ACPf, the critical streamwise wavenumber $\alpha _c$ generally displays an initial monotonic decline towards $\alpha _c = 0$, the latter limit corresponding to the complete loss of linear instability in ACPf beyond $\xi = \xi _A$. However, even for the smallest shear angles treated here, we see that the critical streamwise wavenumber for OCPfs is only close to, but never exactly, zero. Furthermore, in conjunction with the critical Reynolds number, the $\alpha _c$ curves also level out at sufficiently high wall speeds, reinforcing the presence of a limiting regime in modal stability. Meanwhile, the critical spanwise wavenumber $\beta _c$ is generally non-zero, a reminder of the three-dimensional nature of the base flow and the consequent inapplicability of Squire's theorem. An especially interesting behaviour is observed in figure 6(b) for a short range around $\xi \approx 0.3$, where $\beta _c\approx 0$. Here, as shown in figure 5, the $Re_c$ curves for OCPfs in $\varTheta _1$ also experience an inflection point. The latter is a key stability feature in ACPf, and, in this range, Potter (1966) had noted that $\xi \approx c_{r,c}$, the real part of the $x$-phase speed at criticality, hinting at some sort of link between this equality and the accompanying destabilization. We were able to verify this relation for $\varTheta _1$ as well, suggesting that its secondary effect here is a preference for spanwise-independent modes (note that this is implicit for ACPf). Finally, just before this inflectional region, for $\theta =20^\circ$, we can resolve a rather dramatic trough in the $\alpha _c$-curve, which seems to temporarily terminate at the corresponding asymptotic ($\xi$-independent value) before recovering to its original trajectory. We interpret this as the first indication of an imminent departure from the modal characteristics of ACPf, which naturally leads to a discussion of $\varTheta _2$, the second stability regime defined by $20^\circ < \theta \leq 90^\circ$.

Figure 6. Curves of the critical wavenumbers (a) $\alpha _c$ and (b) $\beta _c$ versus $\xi$ for $\varTheta _1\equiv (0, 20^\circ ]$. As usual, a dashed line represents the equivalent plots for ACPf (note that for the latter, Squire's theorem implies $\beta _c = 0$ for all linearly unstable wall speeds). Asymptotic behaviour similar to the curves for $Re_c$ is observed for all $\theta$. Furthermore, at wall speeds beyond the cutoff value $\xi _A$ for ACPf, the $\alpha _c$ curves appear to once again have been pulled away from $\alpha _c = 0$ as $\theta$ increases.

In particular, figure 7 highlights for $\varTheta _2$ a noticeable shift in the stability characteristics of OCPfs. The $Re_c$ curves lose their inflectional nature as in $\varTheta _1$, and while increasing the shear angle still induces destabilization, the decrease in $Re_c$ at any $\xi$ is uniform in $\theta$. Furthermore, all critical parameters reach their asymptotic values at smaller wall speeds, doing so by following relatively smoother trajectories (compare, for example, with the $\beta _c$ curves in figure 6). From the perspective of the critical Reynolds number, the most unstable OCPf configurations occur as $\theta \to 90^\circ \in \varTheta _2$, when the wall motion is perfectly orthogonal to the direction of the pressure gradient. In this limit, the streamwise and spanwise velocities reduce to

(3.10a,b)\begin{equation} U = 1-y^2, \quad W = \frac{\xi}{2}(1+y), \end{equation}

which are, respectively, pPf and Couette profiles. At this angle, we found that the critical parameters remained invariant for all the wall speeds studied here, approaching, in fact, the equivalent tuple for pPf. Specifically, we had

(3.11)\begin{equation} (\alpha_c, \beta_c, Re_c)_{\theta=90^\circ} = (1.02, 0, 5773.22), \end{equation}

effectively indicating for this $\theta$ a superposition of the stability of the individual velocity components (note that the Couette flow is always linearly stable; see Romanov 1973). To rationalize this, we recall that the influence of the spanwise cross-flow on the OS operator is partially modulated by the wavenumbers, particularly through the effective wall speed $\xi _{{eff}}$. Since the corresponding eigenvalue problem can always be mapped to an equivalent one for ACPf, one would wish to somehow negate the Couette contribution, which is known to be stabilizing, in order to ‘maximize’ instability. This is achieved most optimally in disturbances with $\beta = 0$, immediately reducing the OS operator to that for pPf under Squire's theorem and yielding the critical tuple in (3.11).

Figure 7. The critical Reynolds numbers and Fourier wavenumbers plotted against $\xi$ for some choices of $\theta \in \varTheta _2\equiv (20^\circ, 90^\circ ]$. The black arrow depicts the direction of increasing $\theta$ in increments of $10^\circ$ from $\theta = 30^\circ$ to $\theta = 90^\circ$ (perfect orthogonality). At the latter angle, the critical triplet is constant in $\xi$ and equal to that obtained from an analysis of the two-dimensional OS equation for Poiseuille flow.

For full contour plots of the critical parameters in the $(\xi, \theta )$ plane as well as a discussion of the critical phase speeds, we refer the reader to Appendix B.

3.3. The limiting regime of modal stability

In the previous section, it was observed that when the wall speed is sufficiently high, the stability of OCPfs becomes independent of $\xi$. The values of $\xi _f$, which represents the approximate wall speed that initiates this limiting regime, are shown in figure 8. It can be seen that $\xi _f$ generally decreases with $\theta$, following a roughly linear relationship within the range $\varTheta _2$, where OCPfs demonstrate the strongest deviations from the stability characteristics of ACPf. In this section, by adopting a simple juxtaposition with known results on the linear stability of pPf and ACPf, we aim to derive analytical formulae for the asymptotic values of the critical parameters. A small part of the following argument was briefly mentioned earlier when discussing criticality for $\theta =90^\circ$, but will now be elaborated upon.

Figure 8. The variation with $\theta$ of $\xi _f$, the wall speed at which the critical parameters asymptote. A dashed line indicates the linear law $\xi _f\sim \theta$, which holds well in $\varTheta _2$.

For purely streamwise velocity profiles, a classic argument due to Squire (1933) is that transversal ($\beta = 0$) modes must become unstable at a lower $Re$ than both longitudinal ($\alpha =0$) and oblique ($\alpha,\beta \neq 0$) modes. The proof proceeds by defining a two-dimensional ($\beta\,{=}\,0$) OS problem and arguing that, at criticality, the corresponding Reynolds number $Re_{c, 2D}$ cannot be larger than that of the three-dimensional ($\beta \neq 0$) case $Re_{c,3D}$. Now, consider ACPf, $\theta =0$, and recall that in the limit $\xi \to 0$, pPf can be recovered. In particular, an application of Squire's theorem to both these flows yields

(3.12a,b)\begin{equation} Re_{c, 3D}^{pPf} > Re_{c, 2D}^{pPf},\quad Re_{c, 3D}^{ACPf} > Re_{c, 2D}^{ACPf}. \end{equation}

However, from the work of Potter (1966), we know that $Re_{c, 2D}^{ACPf} \geq Re_{c, 2D}^{pPf}$, which immediately allows us to conclude that

(3.13)\begin{equation} Re_{c, 3D}^{ACPf} > Re_{c, 2D}^{pPf}, \end{equation}

as well. For OCPfs, due to the mean three-dimensionality, Squire's theorem is no longer valid. However, (3.7) implies that there exists a one-to-one mapping of the OS eigenproblem for an arbitrary OCPf, initialized with any combination of the wavenumbers, to an equivalent (in general) three-dimensional one for ACPf. Therefore, we have that

(3.14)\begin{equation} Re_c (\theta\neq 0, \xi\neq 0; \alpha,\beta) = Re_c (\theta = 0, \xi = \xi_{eff}; \alpha,\beta)\geq Re_{c, 2D}^{pPf}. \end{equation}

Here, the merits of the effective wall speed $\xi _{eff}$ are once again apparent, as it can be made as large or as small as possible due to its dependence on $\theta$ and, more importantly, on the wavenumbers themselves. Therefore, to optimize the instability for an OCPf, which is equivalent to ‘$\geq$’ in (3.14) approaching equality, the OS operator should degenerate precisely into the two-dimensional analogue for pPf. This can happen if and only if

(3.15a–c)\begin{equation} \xi_{eff} = \xi(\cos\theta+\gamma \sin\theta) = 0, \quad \alpha^2+ \beta^2 = \alpha_{c,{pPf}}^2, \quad \alpha Re = \alpha_{c,{pPf}}Re_{c,{pPf}}, \end{equation}

where $\gamma = \beta /\alpha$ and $(\alpha _{c,{pPf}}, Re_{c,{pPf}})\approx (1.02, 5773.22)$ (see Orszag 1971). Assuming $\xi \neq 0$, the first of these three equations yields

(3.16)\begin{equation} \gamma ={-}\cot\theta, \end{equation}

a result that can be combined with the remaining constraints in (3.15a–c) to obtain the following closed solutions:

(3.17a–c)\begin{equation} \alpha = \alpha_{c,{pPf}} |\sin \theta|, \quad \beta ={-}\alpha_{c,{pPf}} \,\mathrm{sgn}(\sin \theta)\cos \theta ,\quad Re = Re_{c,{pPf}}|\csc\theta|, \end{equation}

where $\mathrm {sgn}$ represents the signum function and we have selected the positive solution for $\alpha$. Figure 9 presents the asymptotic values of the critical parameters, denoted by the subscript $\langle \cdot \rangle _f$, obtained from our numerical results overlaid with the analytical solution in (3.17a–c). An almost exact match is observed up to the resolution error of the domain sweep.

Figure 9. The asymptotic values of ($a$,$b$) the critical streamwise and spanwise wavenumbers and (c) the critical Reynolds number versus $\theta$. The solid line denotes the theoretical estimate provided in (3.17a–c).

Some crucial remarks can now be made. As the shear angle approaches zero, (3.17a–c) claims that the asymptotic streamwise wavenumber vanishes, while the asymptotic spanwise wavenumber experiences a discontinuity. Although this may seem erroneous at first glance, we note that the critical parameters in ACPf continuously vary with wall speed, showing no limiting behaviour, so (3.17a–c) has no meaning in this limit anyway. On the other hand, figure 6(a) seems to suggest that a flattening of the critical parameters for ACPf, if it occurs, should do so for $\alpha \to 0$. Meanwhile, we note that for $\theta \to {\rm \pi}/2$, an asymptotic spanwise wavenumber $\beta _f=0$ is predicted, which validates our findings in § 3.2. Finally, an additional interesting consequence arises when considering the wavenumber vector $\boldsymbol {k} = (\alpha \ \beta )^\intercal$. In wave theory, the wavenumber vector encodes the direction of wave motion, and a wave with wavenumber vector $\boldsymbol {k}$ propagates at an angle $\psi$ to the positive streamwise direction, where $\tan \psi = \beta /\alpha = \gamma$. From (3.16), we can then conclude that the asymptotic eigenmode propagates at an angle $\psi =\theta - {\rm \pi}/2$ to the pressure gradient, that is, exactly perpendicular to the wall motion.

3.4. Eigenmodes and linear energetics

Figure 10 illustrates for $\xi =0.35$ the spatial distributions of $u$ and $w$, which are, respectively, the streamwise and spanwise velocity perturbations associated with the most unstable eigenmode at criticality. For $\theta =10^\circ \in \varTheta _1$, it is observed that both disturbance components propagate approximately parallel to the streamwise direction. This can be attributed, in part, to the critical spanwise wavenumber being close to zero, which, as shown in figure 6(b), is typically the case for intermediate wall speeds within this range of angles. The streamwise fluctuations bear some resemblance to a Tollmien–Schlichting wave, suggesting a possible similarity between the mechanisms of modal transition in ACPf and weakly skewed OCPfs. However, the spanwise fluctuations are non-zero and consist of weakly parallel flattened structures that are localized near the lower wall. The exact role of these structures in the transition process is not immediately clear and requires further numerical investigation, which is beyond the scope of this paper. On the other hand, for $\theta = 45^\circ \in \varTheta _2$ and $\theta = 135^\circ$, $\xi =0.35\approx \xi _f$, indicating that the stability characteristics of both flows are close to the asymptotic regime. Therefore, the most unstable wavenumber pair satisfies (3.16) and is oblique. Consistent with the argument presented in § 3.3, we observe that the associated eigenmode propagates exactly perpendicular to the direction of motion of the wall. The streamwise and spanwise fluctuations are qualitatively similar, both consisting of vortices tilted slightly away from each end of the channel. The cross-sections of these vortices for $u$ are somewhat distorted, while for $w$ they have a more regular, elliptical shape. Furthermore, we found that between $\theta$ and $180^\circ - \theta$, the support of these structures moved from the lower to the upper wall, although this effect was not very noticeable.

Figure 10. Iso-surfaces of the (a,c,e,g) streamwise $u$ and (b,d, f,h) spanwise $w$ velocity fluctuations for the most unstable eigenmode at $\xi =0.35$ for different $\theta$. For each case, the blue and red contours represent 25 % of the (signed) minimum and maximum values of the perturbations, respectively.

The most interesting behaviour is observed for $\theta = 90^\circ$, when the pressure gradient and the wall velocity vectors are exactly orthogonal. First, we recall from (3.11) that despite the value of $\xi$ for this configuration, $\beta _c = 0$, and the OS operator can always be identified with that for the pPf. Thus, the wall-normal and streamwise components (see (2.19)) of the disturbance are identically invariant with the wall speed, precisely reducing to the Tollmien–Schlichting instability found in pPf. However, unlike the latter flow, for which $\beta \to 0$ induces a vanishing normal vorticity in the (resulting) homogeneous Squire equation, the spanwise fluctuations are no longer zero due to the cross-flow shear in the off-diagonal forcing term, (2.16a,b). They comprise bands of transverse arch-like structures that, to the best of our knowledge, have not been previously recorded in the linearized analysis of any canonical shear flow. Moreover, varying $\xi$ had little effect on the shape of these modes.

Figure 11 shows the linear energy budget at criticality for the angle $\theta = 30^\circ$. As the wall speed increases from $\xi =0.1$ to $\xi =0.2$, the streamwise production $\mathcal {P}_u$ decreases somewhat dramatically in the lower half of the channel, consistent with the initial rise in $Re_c$ seen in figure 7(a). However, for $\xi \geq \xi _f$, it eventually stagnates near the stationary wall, whereas a continuous, although noticeably slow, decline is observed near the moving wall. On the other hand, while $\mathcal {P}_w$ seems to always suffer from a region of negative production near the upper wall, it operates at least one order of magnitude lower than $\mathcal {P}_u$. As a result, the total production $\mathcal {P}\approx \mathcal {P}_u$ remains positive throughout most of the channel. Viscous dissipation decreases with increasing $\xi$ and, similar to $\mathcal {P}_u$ (and, therefore, $\mathcal {P}$), it converges to some extent for high wall speeds. Therefore, even from the standpoint of the linear energy budget, there is a clear indication of modal stability in OCPfs approaching asymptotic regimes, where the potential for exponential amplification at criticality seems to become entirely agnostic to the wall speed.

Figure 11. At criticality, the spatial distribution of the perturbation energy budget terms for $\theta =30^\circ$ and wall speeds between $\xi =0.1$ and $\xi =0.8$ in increments of 0.1. (a) The streamwise $\mathcal {P}_u$ and (b) the spanwise $\mathcal {P}_w$ production, (c) the viscous dissipation $\varepsilon$, (d) the total production ($\mathcal {P}=\mathcal {P}_u+\mathcal {P}_w$) and the Reynolds stresses (e) $\tau _u$ and ( f) $\tau _w$.

The spatial variation of the Reynolds stresses is also shown in figure 11(e, f). Regarding $\tau _u$, we see that while it is primarily negative near the stationary wall, it is always positive near the moving wall. This is in stark contrast to ACPf and many of its variants (e.g. Sadeghi & Higgins 1991; Nouar & Frigaard 2009; Guha & Frigaard 2010), for which an increasing wall speed also generates a region of negative stress in the upper half of the channel. Since the base streamwise shear in ACPfs is (typically) also negative throughout this region, the latter phenomenon decreases the overall energy production, stabilizing the flow (see Appendix A). A related commentary can be made on the changes and eventual disappearance of the critical layers at each wall, which, for strictly streamwise flows, are the wall-normal locations where the streamwise velocity matches the real part $c_r$ of the $x$-phase speed. In ACPf, the critical layer near the moving wall is known to vanish as $\xi$ increases, which is often identified with stabilization. For OCPfs, due to the structure of the OS equation, an analogous argument can be constructed using the effective velocity profile:

(3.18)\begin{equation} U_{eff} = U + \gamma W = 1-y^2 + \frac{\xi_{eff}}{2}(1+y), \end{equation}

which essentially represents a projection of the base velocity in the direction of the wavenumber vector (Schmid & Henningson 2001). Here, we are interested in the wall-normal location(s) $y_c$ such that $U_{eff}(y_c) = c_r$ (note that since we are focusing on criticality, $c = c_r$). In general, since the effective velocity is quadratic in $y$, two such points can exist, associated with each wall, and are explicitly given in closed form as follows:

(3.19)\begin{equation} y_c = \tfrac{1}{4}\left(\xi_{eff} \pm \sqrt{(4+\xi_{eff})^2 - 16c_r}\right). \end{equation}

The distance of each layer from the wall can then be expressed as $\delta _c = 1 - |y_c|$ (see e.g. Nouar & Frigaard 2009), and is shown in figure 12 for criticality at $\theta =10^\circ \in \varTheta _1$ and $\theta =60^\circ \in \varTheta _2$. For the former case, similar to ACPf (see Appendix A), both layers initially approach the channel boundaries. Eventually, past the point of inflection in the associated $Re_c$ curve, the critical layer near the moving wall vanishes (i.e. $\delta _c\to 0$), whereas the one near the lower wall begins to move further towards the centreline, stabilizing the flow. Mathematically, this can be attributed to changes in the asymmetry of $U_{eff}$, which, in turn, are influenced by variations in the effective wall speed $\xi _{eff}$. However, a particularly intriguing behaviour is observed when $\xi \geq \xi _f$. In this case, $\delta _c$ for the moving wall experiences a sudden increase from zero to a roughly constant value approximately equal to that of the lower wall. This can be explained by the constraints on the asymptotic wavenumber pair, (3.15a–c), which enforce $\xi _{eff}\approx 0$ in the asymptotic regime and reduce $U_{eff}$ to a symmetric parabolic profile. As a result, the absolute values of the roots in (3.19) coalesce and remain $|y_c| \leq 1$, which appears to be consistent with the plateauing observed in the critical parameters (figures 5 and 6). In a similar vein, for $\theta =60^\circ$, because the asymptotic regime is accessed earlier (that is, at smaller $\xi$), the effective velocity profile is almost always perfectly symmetric, so that the distance $\delta _c$ remains roughly identical and non-zero for both walls throughout most of the $\xi$ range explored here. On a separate note, for $\theta =60^\circ$, figure 12(b) highlights an initial increase in $\delta _c$ near the stationary wall, which seems to support the monotonicity of the corresponding $Re_c$ curve in this range (figure 7a). Furthermore, considering the behaviour of the critical layers when the wall speed approaches the point of inflection in the $Re_c$ curves for ACPf and OCPfs in $\varTheta _1$, it is likely that the absence of such a feature for $\varTheta _2$ is related to the fact that both critical layers remain intact at these angles.

Figure 12. For representative $\theta$, the development of $\delta _c$ for each critical layer. A solid versus dashed line is used to distinguish the lower, stationary, wall from the upper, moving, wall. In each case, an inset illustrates the effective mean velocity profiles $U_{eff}$ for values of the wall speed between $\xi = 0.2$ and $\xi =1$ in increments of 0.1. (a) For $\theta =10^\circ$, an arrow depicts the direction of increasing $\xi$ (note that $\xi \to \xi _f$ implies $\xi _{eff} \to 0$). (b) For $\theta =60^\circ$, while not immediately apparent, the effective velocity profiles for the wall speeds chosen here coincide almost exactly.

4. Non-modal analysis

We now proceed to investigate the potential for non-modal growth in OCPfs, focusing first on the resolvent ${{\boldsymbol{\mathsf{R}}}}$ and its associated energy norm $\mathcal {R}\equiv \lVert {{\boldsymbol{\mathsf{R}}}}\rVert _E$. For some representative pairs of wavenumbers, figure 13 visualizes $\mathcal {R}_{max}$, defined as

(4.1)\begin{equation} \mathcal{R}_{max} = \max_{\zeta} \mathcal{R}(\alpha,\beta, Re, \xi, \theta, \zeta), \end{equation}

for $Re = 1500$, which is subcritical for all possible configurations. As is the case with many of the quantities explored in § 3, the periodicity of the base flow in $\theta$ once again embeds itself in $\mathcal {R}_{max}$, which appears to be ${\rm \pi}$-periodic. On the other hand, the forcing frequency $\zeta$ that gives rise to $\mathcal {R}_{max}$ was found to be $2{\rm \pi}$-periodic, though it is not shown here. Comparing figure 13(a,b), we immediately see that the response of spanwise-independent disturbances is relatively damped compared with disturbances with $\alpha =0$, a difference that spans approximately an order of magnitude. In the former case, allowing for a weak $z$ dependence by steadily varying $\beta$ produced only a negligible increase. However, for disturbances of the second kind, changes in the streamwise wavenumber appeared to have a more diverse effect. In particular, while $\mathcal {R}_{max}$ decreased for, say, $\theta \lessapprox 90^\circ$, it actually increased for $90^\circ \lessapprox \theta \lessapprox 165^\circ$. This is significantly different from many classical two-dimensional flows, where it is usually disturbances completely independent of $x$ that elicit the most vigorous response (Trefethen et al. 1993; Schmid & Henningson 2001).

Figure 13. Plots of $\mathcal {R}_{max}$, the maximum of the resolvent energy norm across all forcing frequencies $\zeta$ for $\xi =0.35$ and $Re = 1500$: (a) $\alpha =1$ and $\beta$ varied in increments of $0.1$; (b) $\beta =1$ and $\alpha$ varied in increments of $0.05$. In (b), the inset zooms in on the region where $\mathcal {R}_{max}$ appears to increase in conjunction with $\alpha$. Note the periodicity in $\theta$.

For more general forcing frequencies, $\zeta \in \mathbb {C}$, figure 14 illustrates for $\theta = 30^\circ$ the contours of $\lVert {{\boldsymbol{\mathsf{R}}}}^{-1}\rVert _{-E}$ in the complex plane, where $\lVert \cdot \rVert _{-E}$ is the ‘inverse’ energy norm:

(4.2)\begin{equation} \lVert {{\boldsymbol{\mathsf{R}}}}^{{-}1}\rVert_{{-}E} = \sigma_{min}({{\boldsymbol{\mathsf{F}}}}{{\boldsymbol{\mathsf{R}}}}^{{-}1}{{\boldsymbol{\mathsf{F}}}}^{{-}1}), \end{equation}

in which $\sigma _{min}$ denotes the smallest singular value of the operator ${{\boldsymbol{\mathsf{F}}}}{{\boldsymbol{\mathsf{R}}}}^{-1}{{\boldsymbol{\mathsf{F}}}}^{-1}$. Here, we have chosen to focus on wall speeds $\xi \geq \xi _f$ and the associated asymptotic parameters $(\alpha _f, \beta _f, Re_f)$ (§ 3.3). Within the paradigm of the $\epsilon$-pseudospectra, the interpretation of these plots is as follows. Since the properties of the $2$-norm imply

(4.3)\begin{equation} \mathcal{R} = \sigma_{max}({{\boldsymbol{\mathsf{F}}}}{{\boldsymbol{\mathsf{R}}}}{{\boldsymbol{\mathsf{F}}}}^{{-}1}) = \sigma_{min}^{{-}1}({{\boldsymbol{\mathsf{F}}}}{{\boldsymbol{\mathsf{R}}}}^{{-}1}{{\boldsymbol{\mathsf{F}}}}^{{-}1}) = \lVert {{\boldsymbol{\mathsf{R}}}}^{{-}1}\rVert_{{-}E}^{{-}1}, \end{equation}

the set $\varLambda _\epsilon$ in (2.31) admits the alternative definition

(4.4)\begin{equation} \varLambda_\epsilon = \{\zeta\in\mathbb{C}\colon\lVert {{\boldsymbol{\mathsf{R}}}}^{{-}1}\rVert_{{-}E} \leq \epsilon\}. \end{equation}

Thus, within the level curve $\lVert {{\boldsymbol{\mathsf{R}}}}^{-1}\rVert _{-E}=\epsilon$, $O(\epsilon ^{-1})$ amplification can be realized. Additionally, the extent to which these contours protrude into the upper half-plane can be connected to the potential for transient energy amplification. In particular, Reddy et al. (1993) showed that such growth cannot occur (that is, $G\leq 1$) if and only if $\beta _\epsilon \leq \epsilon$ for all $\epsilon \geq 0$, where

(4.5)\begin{equation} \beta_\epsilon = \sup_{\zeta\in\varLambda_\epsilon({{\boldsymbol{\mathsf{S}}}})}\Im(\zeta) \end{equation}

is the pseudospectral radius. The significance of the restriction on $\beta _\epsilon$ lies in noting that for a normal operator, the 2-norm $\epsilon$-pseudospectra comprise closed balls of radius $\epsilon$ centred around the eigenvalues (Trefethen & Embree 2005). Therefore, for a linearly (not necessarily asymptotically) stable normal operator, $\beta _\epsilon$ reaches its maximum at $\beta _\epsilon = \epsilon$, specifically for a marginally stable mode. More concretely, the pseudospectral radius can also provide a reasonably sharp lower bound for unforced algebraic growth (Trefethen et al. 1993). Returning to figure 14, we observe for all wall speeds pseudoresonance down to $\epsilon \approx 10^{-8}$. Larger $\xi$ only expand and widen the pseudospectral contours, particularly toward the unstable half of the complex plane, which indicates an increasing susceptibility to non-modal mechanisms despite the convergence of modal features. Note also that in figure 14, the associated spectra appear to be invariant with $\xi$. This is, of course, a consequence of the asymptotic triplet enforcing a reduction to the ‘same’ eigenvalue problem (recall § 3.3 and see also Appendix C). An immediate consequence is that all variations in non-modal growth within this regime must stem exclusively from the off-diagonal forcing term, since any $\xi$ dependence in $\mathcal {L}_{OS}$ and $\mathcal {L}_{SQ}$ necessarily drops out.

Figure 14. For $\theta =30^\circ$ and (a–d) various $\xi \geq \xi _f$, the logarithmic level curves at $(\alpha _f,\beta _f,Re_f)$ for the $\epsilon$-pseudospectra from $\log \epsilon = -1$ to $\log \epsilon = -8$ (outer to inner) in decrements of $-1$. The OS and Squire modes in the $x$-phase speed formulation are depicted via circles. A dashed line indicates the stability boundary $c_i = 0$.

We now turn our attention to the dynamics of the unforced initial value problem, (2.15). In particular, we are interested in $G_{max}$, defined for an OCPf configuration $(\xi, \theta )$ as the maximal amplification in time and wavenumber space:

(4.6)\begin{equation} G_{max}(Re, \xi, \theta) = \max_{\alpha, \beta, t} G(\alpha, \beta, Re, \xi, \theta, t). \end{equation}

Figure 15 outlines the findings of a large parameter sweep for $G_{max}$ at $Re=1000$, slightly above the range for transition in ACPf as quoted, for example, by Tsanis & Leutheusser (1988) and Klotz et al. (2017) (note, however, that their Reynolds numbers are based on the wall velocity). As in § 3, all results are presented relative to those of ACPf, $\theta = 0$, which experiences a monotonic increase in $G_{max}$ with $\xi$. The introduction of a weak misalignment maintains this trend, but allows for greater amplification throughout the full range of wall speeds explored here. This effect was determined to be most pronounced at $\theta \approx 4.5^\circ$. At even larger shear angles, two different regimes can be identified in $\xi$. In particular, while $G_{max}$ continues to grow with $\theta$ (albeit slowly) for $0<\xi \lessapprox 0.15$, it tends to decrease quite rapidly for $0.2\lessapprox \xi \leq 1$. Furthermore, at least for wall speeds in this range, no asymptotic behaviour was resolved for $G_{max}$, which is in sharp contrast to our modal calculations.

Figure 15. For various $\theta$, curves of $G_{max}$, the largest possible energy gain exhibited by OCPfs across time and wavenumber space. The black dashed line indicates the equivalent plot for the aligned case, $\theta = 0$. In general, the largest amplification is realized for small but non-zero angles, peaking at $\theta =4.5^\circ$ for most $\xi$. A greater degree of skewness in the flow tends to suppress the amplification, particularly for modest to large wall speeds.

Interestingly, within the paradigm of transient growth, it is apparent that larger values of $\theta$ are typically the most ‘stable’, with $\theta = 90^\circ$ providing the strongest reduction in $G_{max}$ for a wide range of wall speeds. The latter observation stands, of course, in strong opposition to the results presented in the previous section, particularly (3.11), which claims that a perfectly orthogonal OCPf configuration is, in fact, capable of minimizing $Re_c$ in the $(\xi,\theta )$ plane. An antagonistic effect, therefore, appears to be at play here, since, individually, both ACPf and the standard Couette flow support strong transient responses, yet for sufficiently skewed OCPfs, $G_{max}$ can drop to as low as $46\,\%$ of the equivalent value for pPf at this Reynolds number ($G_{max}\approx 196$). From a mathematical perspective, one can attribute this to the nonlinearity of the operator norm or to the fact that, contrary to modal analysis, we are now investigating the full OSS system, for which $\xi _{{eff}}$ no longer constitutes an informative parameter. Physically, however, an intriguing analogy can be drawn to fully turbulent three-dimensional boundary layers, for which increasing skewness, in the mean sense, is known to dampen the generation of Reynolds stresses and, therefore, the production of turbulent kinetic energy relative to the two-dimensional case (Eaton 1995; Coleman et al. 1996; Johnston & Flack 1996; Lozano-Durán et al. 2020). Although the physical mechanisms involved are not yet well understood, it is often believed that the addition of a mean spanwise strain detracts large momentum-carrying eddies from their optimal alignment (Van Den Berg et al. 1975; Bradshaw & Pontikos 1985). In the context of laminar OCPfs as treated here, one can partially quantify the existence of such an ideal configuration by considering $\theta _{max}$:

(4.7)\begin{equation} \theta_{max}(Re, \xi) = \{\theta^\prime\mid G_{max}(Re, \xi, \theta^\prime) = \max_{\theta}G_{max}(Re, \xi, \theta)\}, \end{equation}

which, at a given $\xi$, represents the angle of wall motion that achieves the most vigorous algebraic amplification. Figure 16(a) highlights that $\theta _{max}$ decays primarily as a power law. More importantly, as shown in figure 16(b,c), the cross-flow component $W$ associated with $\theta _{max}$ is quite weak, allowing the flow direction $\phi$ to collapse throughout most of the channel and experience rapid variation only near the upper wall. To further visualize this, an average skewness $\langle \phi \rangle$, defined as

(4.8)\begin{equation} \langle \phi\rangle = \tfrac{1}{2}\int_{{-}1}^1 \phi(y)\,\mathrm{d} y, \end{equation}

is also plotted in the inset of figure 16(b) and remains small ($\lessapprox 6^\circ$ at best) for all $\xi$. Thus, the optimal configuration for energy growth in OCPfs appears to be an approximately collateral boundary layer, with a flow direction roughly constant in $y$ and almost aligned with the streamwise axis.

Figure 16. (a) The variation in $\xi$ of $\theta _{max}$, the shear angle that optimizes $G_{max}$. Some combinations of $(\xi,\theta _{max})$ have been selected and the associated flow directions and cross-flow profiles highlighted with the appropriate colour in (b,c), respectively. Our conclusions are robust to the choice of these pairs. The inset in (b) shows the $y$-averaged deviation $\langle \phi \rangle$ of the optimal net base flow from the streamwise direction, (4.8).

In figure 17, we present the contours of the wavenumbers and the time $t$ that achieve the maximum amplification $G_{max}$. For purely streamwise flows, this optimal gain is generally observed for longitudinal modes, $\alpha _{max} = 0$ (Trefethen et al. 1993; Schmid & Henningson 2001). However, this may or may not be the case for three-dimensional flows. Indeed, for OCPfs, we found the maximum amplification to occur for small but often non-zero streamwise wavenumbers, with little overall variation in the $(\xi, \theta )$-space. In some cases, streamwise-invariant disturbances remained optimal, but these were the exception rather than the rule. On the contrary, the optimal spanwise wavenumbers fluctuated more strongly, varying from $\beta _{max} \approx 2.05$ to $\beta _{max}\approx 2.25$ and even down to $\beta _{max}\approx 1.3$. A rudimentary explanation for these trends can be obtained by recalling that the non-normality of the OSS operator ${{\boldsymbol{\mathsf{S}}}}$ is partially tied to the mean shear(s) coupling the velocity perturbations in the Squire equation through the operator ${\rm i} \beta \mathcal {D} U - {\rm i} \alpha \mathcal {D}W$ (see (2.13) and (2.16a,b)). Thus, in the same spirit as the effective wall speed $\xi _{{eff}}$, one can attempt to split this operator into its ‘Poiseuille’ and ‘Couette’ constituents:

(4.9)\begin{equation} {\rm i} \beta\mathcal{D} U - {\rm i} \alpha \mathcal{D}W = {\rm i}({-}2\beta y + \xi_{off}), \end{equation}

where we have defined

(4.10)\begin{equation} \xi_{{off}} = \frac{\xi}{2}(\beta\cos\theta - \alpha\sin\theta). \end{equation}

Since the Poiseuille contribution in (4.9) is agnostic to the wall motion and does not, at any rate, favour either half of the channel, maximizing the effective mean shear is equivalent to maximizing $|\xi _{{off}}|$. Then, assuming $\beta >0$, smaller shear angles would be biased towards larger spanwise wavenumbers because $\sin \theta \approx 0$. On the contrary, for larger $\theta$, $\cos \theta$ is small and $\alpha$ should ideally become increasingly negative, though figure 17(a) shows that this preference seems to emerge only at higher wall speeds. In figure 17(d), we see that $\xi _{{off}}$ somewhat emulates the changes in $G_{max}$ at this Reynolds number. For example, while its maximum in the $(\xi, \theta )$ plane occurs in the purely parallel, high-$\xi$ limit, its minimum is realized as $\theta \to 90^\circ$. Of course, however, this correlation is bound to be imperfect, particularly because $\xi _{{off}}$, in this case, cannot be a definitive statistic. The nature of the remaining blocks in ${{\boldsymbol{\mathsf{S}}}}$ is equally important and does not admit a simple interpretation, calling for a more complex analysis that is outside the scope of this article. Finally, the time $t_{max}$ taken to achieve the maximum growth appears to be the longest when an OCPf is weakly oblique and generally decreases in the direction of increasing $\xi$ and $\theta$. Note that, for the latter combination of flow parameters, $G_{max}$ also tends to a minimum, indicating that transient phenomena for these configurations operate on shorter time scales and are likely suppressed by the viscosity before reaching sufficient amplitudes to trigger further instability. We remark in passing that while a search was conducted for $\theta \in (90^\circ, 180^\circ ]$, $G_{max}$ and all associated optimal parameters were found to be symmetric around $\theta =90^\circ$. An exception to this was $\alpha _{max}$, which was determined to be anti-symmetric, a result that is not suggested by any operator-level symmetries.

Figure 17. The variation in $(\xi, \theta )$-space of (a) the streamwise wavenumber $\alpha _{max}$, (b) the spanwise wavenumber $\beta _{max}$ and (c) the time $t_{max}$ at which the maximum energy amplification $G_{max}$ is attained. In the case of $\alpha _{max}$, a dashed line (not obeying the colour-bar) indicates the level curve $\alpha _{max} = 0$. (d) Contours of $\xi _{{off}}$, which seem to be positively correlated to $G_{max}$.

Figure 18 concludes this section by visualizing for $\theta ={\rm \pi} /4$ and $\theta = {\rm \pi}/2$ the initial condition and response pair associated with $G_{max}$. Respectively, these are the first right and left singular functions of the state transition operator $\varPhi (t,0)$. In two-dimensional flows, the optimal initial field, typically observed for $\alpha = 0$, is characterized by weak counter-rotating streamwise vortices that evolve through a redistribution of horizontal mean momentum by the normal velocity fluctuations to form high-energy streaks at $t = t_{max}$. Originally proposed in the works of Ellingsen & Palm (1975) and Landahl (1975, 1980), this process is commonly referred to as the lift-up effect, in which a linear amplification in time proportional to the streamwise shear can be achieved for a streamwise-independent disturbance, at least in the inviscid limit. In the viscous alternative, this growth would continue to persist, but only to the leading order before decaying due to viscosity (Brandt 2014). Although the three-dimensionality of our flow introduces additional nuance, Ellingsen & Palm (1975) had suggested that the lift-up process could remain viable even in skewed boundary layers, arguing, however, that streak growth would substantially decrease, particularly so in the case of OCPfs because the streamwise shear itself decreases as $\theta \to {\rm \pi}/2$. Indeed, Corbett & Bottaro (2001), for example, found that streamwise streaks developing via lift-up remained the optimal disturbance in Falkner–Skan–Cooke boundary layers.

Figure 18. The optimal initial condition (a,d) and response (b,e) pair for $\theta = {\rm \pi}/4$ (a,b) and $\theta = {\rm \pi}/2$ (d,e) at $\xi = 0.25$. Black-to-white shades denote the streamwise vorticity, blue-to-red the streamwise velocity $u$ and arrows the cross-stream components $v$ and $w$. All quantities represent disturbance variables. (c,f) The time evolution of $E_w$, the domain-averaged cross-flow (i.e. spanwise) perturbation energy; here, the dashed line represents the energy-optimal initial condition for ACPf.

Therefore, it is not surprising that the optimal initial conditions in figure 18 comprise weak streamwise vortices whose amplification at the optimal time decreases in response to an increase in flow obliqueness. However, what is more interesting is that these vortices initially oppose and are eventually rotated in the direction of the mean spanwise shear. This, as noted in Blesbois et al. (2013) and Hack & Zaki (2015), allows for an additional non-modal energy gain in the cross-flow perturbation component through the inviscid down-gradient Reynolds stress mechanism proposed by Orr (Orr 1907; Butler & Farrell 1992). One particular physical characterization of this process invokes conservation of circulation around material circuits (e.g. Waleffe 1995a), and is visualized in figure 18 through the disturbance $x$-vorticity (note that $\alpha _{max}\approx 0$ here, but this discussion can be made robust by adopting a coordinate system aligned with the streak direction – see Blesbois et al. (2013)). Sheets of perturbation streamwise vorticity lean against the background cross-flow shear at initial time (figure 18a,d) and are subsequently realigned to a down-shear orientation (figure 18b,e), driving a transient amplification in the cross-flow disturbance energy (figure 18c, f). This growth is absent in the energy-optimal initial conditions for two-dimensional base flows, such as ACPf, which develop exclusively via lift-up and involve no feedback on the cross-stream components, which then must monotonically decay due to viscous dissipation. In short, it appears that at least up to the stage of primary instability, the route to transition for OCPfs is dominated by a lift-up process interspersed by an Orr-type mechanism. Furthermore, since the Orr mechanism is enhanced by the presence of increasing (spanwise) shear and because the ultimate contribution of the cross-flow perturbation energy to $G_{max}$ is somewhat negligible (see figure 15), it is likely that the trends observed in the maximal amplification can be attributed to a decrease in the effectiveness of the lift-up process.

5. Conclusion

We performed a comprehensive modal and non-modal stability analysis in OCPfs, which are described by a wall motion at an angle $\theta$ to the pressure gradient. These are generalizations of the traditional aligned case and, to the best of our knowledge, have not received prior attention in the stability literature.

We derive the corresponding OSS system, identifying by a simple analogy an effective wall speed that completely characterizes modal solutions. A large-scale numerical sweep reveals that, in general, a misalignment between the pressure gradient and the wall velocity is destabilizing, at least relative to the aligned case. Considerations of symmetry and periodicity allow for a restriction of the parameter space to $\theta \in [0^\circ, 90^\circ ]$ and, in this range, two regimes are identified. For shear angles $0^\circ <\theta \lessapprox 20^\circ$, almost all stability features in OCPfs continue from the aligned case. On the other hand, the range $20^\circ \lessapprox \theta \leq 90^\circ$ demonstrates sharper differences, including, in particular, the lack of a trademark inflection point in the $Re_c$ curves as observed for ACPf. The linear energy budget and the movement of the critical layers generated by the effective velocity profile seem to confirm these trends. Modal instability is optimized by the perfectly orthogonal configuration, $\theta = 90^\circ$, which exhibits a constant critical tuple for all $\xi$. For all $\theta \neq 0$, we find that unstable modes persist throughout $\xi \in [0, 1]$, notably distinct from ACPf, where $\xi \approx 0.7$ marks the transition to a regime of unconditional linear stability. This behaviour is accompanied by a convergence of the critical parameters starting at the threshold $\xi = \xi _f$, which appears to decrease with $\theta$. A simple theoretical analysis explains the latter phenomenon and derives the exact asymptotic values of the critical parameters in the limiting regime.

On the topic of non-modal disturbances, OCPfs, through the non-normality endowed to the OSS operator by the streamwise and cross-flow shears, exhibit trends that conflict with those commonly quoted for two-dimensional flows. For example, the resolvent norm is not necessarily maximized for disturbances with $\alpha = 0$. Meanwhile, the $\epsilon$-pseudospectra reveal that even if modal stability converges beyond $\xi = \xi _f$, non-modal mechanisms might continue to be amplified by changes in the wall speed. Finally, considering the unforced initial-value problem, the maximum energy amplification $G_{max}$ appears to decrease strongly with the skewness of the base profiles, implying that the imposition of three-dimensionality is generally detrimental to energy growth. This is reminiscent of fully turbulent three-dimensional boundary layers, where increased skewness is known to suppress turbulent energy production. Note that while Hack & Zaki (2015) observed similarly declining energy gains for their Blasius–Stokes flow, Corbett & Bottaro (2001) calculated a stronger transient growth relative to the two-dimensional case in their study on swept boundary layers, suggesting a strong dependence on the particular mechanism enforcing skewness in the base flow. At all wall speeds, however, the configuration that optimizes energy amplification appears to be a weakly three-dimensional collateral boundary layer. Finally, the most energetic initial perturbations seem to develop via a lift-up process complemented by an Orr-like mechanism driven by the spanwise shear, the latter being absent in streak amplification for two-dimensional flows.

A natural extension of this article seems to be through an investigation of OCPfs in the turbulent regime, which is perhaps where a majority of practical applications reside. Although there have been some previous studies, they have primarily focused on the perfectly aligned case (Kim & Lee 2018; Kim et al. 2020; Cheng et al. 2023) or the perfectly orthogonal case (Coleman et al. 1996; Howard & Sandham 1997; Kannepalli & Piomelli 2000; Le et al. 2000), with little or no attention devoted to intermediate $\theta$. Some initial work on the latter configurations has been conducted (e.g. Zhang et al. 2023), but additional effort is needed and will likely contribute well to our overall understanding of the physics in three-dimensional boundary layers.