1. Introduction
Flows between two parallel plates are canonical configurations involving abundant flow stability and turbulence phenomena. Stability analyses of plane Poiseuille (PP) flow driven by an external pressure gradient, and plane Couette (PC) flow induced by bounding, moving walls, are classical examples. For PP flow, asymptotic solutions of the Orr–Sommerfeld equation (Lin Reference Lin1945) provided a critical Reynolds number 
$Re_{cr} \approx 5300$ for modal flow instability, with a more accurate 
$Re_{cr} = 5772.22$ given later by Orszag (Reference Orszag1971). Evidence for linear stability of PC flow was shown by Gallagher & Mercer (Reference Gallagher and Mercer1962) and Deardorff (Reference Deardorff1963). Later, Romanov (Reference Romanov1973) proved that, for modal perturbations, if the perturbation wavenumber is 
$\alpha >0$ and Reynolds number is 
$Re>0$, then the flow is linearly stable.
Subsequent research has focused on extensions such as plane Couette–Poiseuille (PCP) flow, where the plate relative velocity vector and the externally applied pressure gradient are aligned. Stability of PCP was studied by Potter (Reference Potter1966) and Hains (Reference Hains1967), offering numerical evidence for the consensus view that a Couette component stabilizes the flow. Later, Cowley & Smith (Reference Cowley and Smith1985) studied the asymptotic behaviour of neutral curves for PCP flow at large Reynolds numbers, modifying the earlier results of Smith (Reference Smith1979) for PP flow. More recently, Kirthy & Diwan (Reference Kirthy and Diwan2021) conducted an energy budget analysis, pointing out that upper plate movement causes a negative production term, thus leading to stabilization. Guha & Frigaard (Reference Guha and Frigaard2010) and Shankar & Shivakumara (Reference Shankar and Shivakumara2021) investigated PCP stability, emphasizing certain volume flow aspects that correspond to a mechanical alternative for stabilization. Stability of PCP flow for other fluid constitutive relationships has been studied, including non-Newtonian (Nouar et al. Reference Nouar, Kabouya, Dusek and Mamou2007; Sengupta & De Reference Sengupta and De2020; Pal & Samanta Reference Pal and Samanta2023) and non-ideal (Ren, Fu & Pecnik Reference Ren, Fu and Pecnik2019) fluids. Related work focusing on Couette–Poiseuille flow between concentric cylinders was reported by Sadeghi & Higgins (Reference Sadeghi and Higgins1991) and Gittler (Reference Gittler1993) for both inner and outer stationary cylinder motions, respectively. Non-modal operator based methods have been applied to PP, PC and PCP flows (Böberg & Brösa Reference Böberg and Brösa1988; Butler & Farrell Reference Butler and Farrell1992; Reddy & Henningson Reference Reddy and Henningson1993; Bergström Reference Bergström2005; Schmid Reference Schmid2007), providing a novel perspective on transient growth in the early stages of turbulence transition.
A generalized Couette–Poiseuille (GCP) flow was defined by Cheng et al. (Reference Cheng, Pullin, Samtaney and Luo2023) as a channel-like, wall-bounded flow with non-parallel superposition of a PP and a PC flow, thus forming a three-dimensional (3-D) base flow velocity profile as a further extension of PCP flow. Squire's theorem (Squire Reference Squire1933) for PCP flow allows determination of a critical Reynolds number for general perturbations to the plane-parallel base flow by restriction to planar or two-dimensional disturbances. For GCP flow, Squire's theorem cannot be applied directly. In § 2, the GCP flow configuration is defined in terms of superposed base flows, and the relevant stability equation is obtained. Numerical determination of neutral surfaces is described in § 3, and global and local critical points are defined. A mapping from GCP stability to an equivalent PCP stability problem is also identified that allows reduction in the dimensionality of local critical point dependence on the flow parameters. Global GCP stability properties are explored in § 4, which reveals complex behaviour in the dependence of the critical Reynolds number and corresponding wavenumbers as functions of 
$(\theta, \phi )$. Conclusions and future perspectives are discussed in § 5.
2. Problem description: GCP flow
In Cartesian coordinates 
$(\hat {x},\hat {y},\hat {z})$ (‘hat’ denoting dimensional) and velocity components 
$(\hat {U},\hat {V},\hat {W})$, the top and bottom walls at 
$\hat {y}=\pm {h}$ move in the 
$\pm \hat {x}$-directions with 
$\hat {x}$ velocity components 
$\hat {U}_C$ and 
$-\hat {U}_C$, respectively. The base flow is an exact solution of the steady Navier–Stokes equations for an incompressible fluid of constant density and viscosity 
$\rho$, 
$\mu$, respectively, with 
$(\hat {U},\hat {V},\hat {W}) = (\hat {U}(\kern0.07em\hat {y} ),0,\hat {W}(\kern0.07em\hat {y} ))$ in the presence of constant, given, mean pressure gradients in both the 
$\hat {x}$- and 
$\hat {z}$-directions, and with boundary conditions 
$(\hat {U},\hat {V},\hat {W}) = (-\hat {U}_C,0,0)$ at 
$\hat {y} = -h$, and 
$(\hat {U},\hat {V},\hat {W}) = (\hat {U}_C,0,0)$ at 
$\hat {y} =h$. Instead of pressure gradients, we work with bulk mass flow speeds 
$\hat {U}_{M,\hat {x}}=\int _{-h}^{h} \hat {U}(\kern0.07em\hat {y})\,{\rm d}\hat {y} /2h$ and 
$\hat {U}_{M,\hat {z}}= \int _{-h}^{h} \hat {W}(\kern0.07em\hat {y})\,{\rm d}\hat {y} /2h$. The base flow can then be expressed as
\begin{equation} (\hat{U},\hat{V},\hat{W}) =\left( -\frac{3}{2}\,\hat{U}_{M,\hat{x}}\,\left(\left(\frac{\hat{y}}{h}\right)^2-1\right) + \hat{U}_C\left(\frac{\hat{y}}{h} \right) ,0, -\frac{3}{2}\,\hat{U}_{M,\hat{z}}\left(\left(\frac{\hat{y}}{h}\right)^2-1\right)\right). \end{equation}
This is a linear combination of Couette flow in the 
$\hat {x}\unicode{x2013}\hat {y}$ plane and PP flows in both the 
$\hat {x}$- and 
$\hat {z}$-directions. For the present linear base flow, the vector formed by the two bulk flow speeds as components is parallel to the pressure gradient vector. Both form an angle 
$\phi$ to the 
$\hat {x}$-axis, as shown in figure 1.
Configuration of GCP flow with distance 
$2h$ between two walls moving with 
$\hat {x}$ velocities 
${\pm }\hat {U}_C$. The angle 
$\phi$ describes the inclination of the bulk velocity vector to the 
$\hat {x}$-axis.
The parameters 
$(\hat {U}_C,\hat {U}_{M,\hat {x}},\hat {U}_{M,\hat {z}})$ are independent speed variables of the GCP flow. Three independent Reynolds numbers based on these speeds form the all-positive octant of a 3-D Cartesian space. It is convenient instead to work in a spherical coordinate system in this space, with variables 
$(Re,\theta,\phi )$ defined as
\begin{equation} Re=\frac{ \rho
\hat{U}_0 h }{\mu}, \enspace
\theta=\arccos\left( \frac{\hat{U}_C}{\hat{U}_{0}}\right),
\theta \in \left[ 0^\circ, 90^\circ \right], \enspace \phi=\arctan\left(
\frac{\hat{U}_{M,\hat{z}}}{\hat{U}_{M,\hat{x}}}\right),
\phi \in \left[ 0^\circ, 90^\circ
\right],
\end{equation}
where 
$\hat {U}_0= ({\hat {U}_C^2}+ {\hat {U}_{M,\hat {x}}}^2+{\hat {U}_{M,\hat {z}}}^2)^{1/2}$ is a global velocity scale, 
$\mu$ is the dynamic viscosity, 
$\cos \theta$ is a measure of the relative weighting of Couette-type flow and composite base flow, and 
$\phi$ is the angle between the moving plates velocity vector and the bulk velocity vector. The angle 
$\phi$ can also be interpreted as a measure of the bulk velocity ratio in two wall-parallel directions. Dimensionless base flow velocities scaled on 
$\hat {U}_0$ are given by
\begin{align} U=\frac{\hat{U}}{\hat{U}_0}={-}\frac{3}{2}\cos\phi\sin\theta\,(y^2-1)+y\cos\theta,\quad W = \frac{\hat{W}}{\hat{U}_0}={-}\frac{3}{2}\sin\phi\sin\theta \,(y^2-1), \end{align}
with 
$y=\hat {y}/h$. This comprises the non-parallel, 3-D GCP base flow 
$\boldsymbol {U}(y) = (U(y),0,W(y))$, as illustrated in figure 1 in dimensional coordinates.
When 
$\theta \to 0^\circ$, the flow approaches pure PC flow for all 
$\phi$. When 
$\theta \to 90^\circ$, it approaches pure PP flow at an angle 
$\phi$ to the 
$x$-axis. When 
$\phi \to 0^\circ$, the directions of Couette flow and Poiseuille flow coincide. This is PCP flow. When 
$\phi \to 90^\circ$, Couette and Poiseuille flow are orthogonal. The base flow velocity vector always lies in an 
$x\unicode{x2013}z$ plane. For PCP flow, the velocity vector also lies in the 
$x\unicode{x2013}y$ plane, which we refer to as a ‘parallel’ base flow. For 
$\phi >0^\circ$, the velocity vector is generally at a finite angle (which is a function of 
$y$) to the 
$x$-axis, except at 
$y = 1,-1$, where the velocity vector makes angles 
$0^\circ$ and 
$180^\circ$, respectively, with the positive 
$x$-axis. We call this ‘non-parallel’ flow.
2.1. Stability equations
An infinitesimal linear perturbation 
$(\epsilon \tilde {\boldsymbol {u}},\epsilon \tilde {p})$ with 
$\epsilon \ll 1$ is imposed on the base flow 
$(\boldsymbol {U}, P)$, with 
$P$ the base pressure distribution. The perturbed flow is given by
\begin{equation} (\boldsymbol{U}+\epsilon\tilde{\boldsymbol{u}},P+\epsilon\tilde{p}) = ( U+\epsilon\tilde{u}, \epsilon\tilde{v}, W+\epsilon\tilde{w}, P+\epsilon\tilde{p} ). \end{equation}
The disturbance quantities 
$(\tilde {\boldsymbol {u}},\tilde {p})$ are spatially periodic and have the form
\begin{equation} (\tilde{u},\tilde{v},\tilde{w},\tilde{p})= (u(y),v(y),w(y),p(y))\exp({{\rm i}(\alpha{x}+\beta{z}-\varOmega{t})}), \end{equation}
where 
$\varOmega \in {\mathbb {C}}$ is the complex frequency, 
$\alpha,\beta \in \mathbb {R}$ are the 
$x$- and 
$z$-direction wavenumbers, respectively, and 
$u(y)$, 
$v(y)$, 
$w(y)$ and 
$p(y)$ represent perturbation amplitudes. Putting 
$\alpha = k\cos \chi$, 
$\beta = k\sin \chi$, the perturbation can be interpreted as a wave form in the 
$x\unicode{x2013}z$ plane with wavenumber vector oriented at an angle 
$\chi$ to the 
$x$-axis.
The linear modal stability equations are obtained by using the base flow perturbed by (2.5) in the incompressible flow continuity and Navier–Stokes equations, expanding in 
$\epsilon$, and truncating at first order, 
$O(\epsilon )$. This gives
\begin{equation} \left.\begin{gathered} {\rm i}(\alpha{u}+\beta{w})+{D}v = 0, \\ -{\rm i}\varOmega{u}+({\rm i}\alpha{U}+{\rm i}\beta{W})u+v\times {D}U+{\rm i}\alpha{p}=\frac{1}{Re}\,({D}^2-k^2)u, \\ -{\rm i}\varOmega{v}+({\rm i}\alpha{U}+{\rm i}\beta{W})v+{D}{p}=\frac{1}{Re}\,({D}^2-k^2)v, \\ -{\rm i}\varOmega{w}+({\rm i}\alpha{U}+{\rm i}\beta{W})w+v\times {D}W+{\rm i}\beta{p}=\frac{1}{Re}\,({D}^2-k^2)w, \end{gathered}\right\}\end{equation}
where 
$k = (\alpha ^2 + \beta ^2)^{1/2}$ is a resultant wavenumber, and 
$D \equiv \textrm {d}/{\textrm {d}y}$. Since there are three independent spatial variables 
$(x,y,z)$, a stream function for the perturbed flow, which is typically used in PCP stability analysis, does not appear to exist. After some algebra, a modified Orr–Sommerfeld (OS) stability equation is given as
\begin{equation} \left[(U_{P}-{c})+\frac{\alpha}{\alpha_{1}}y\cos\theta\right] ({D}^2-k^{2})v -{v}\times {D}^2{U_{P}}= \frac{1}{{\rm i}{\alpha}_{1}\,Re}\,({D}^2-k^{2})^{2}v, \end{equation}with
\begin{equation} U_{P}(y,\theta)={-}\tfrac{3}{2}(y^2-1)\sin\theta, \quad \alpha_{1} = \alpha\cos\phi+\beta\sin\phi. \end{equation}
Here, 
$U_{P}$ is a characteristic velocity profile denoting the Poiseuille component driven by the bulk mass flow, which gives 
$U=U_{P}\cos \phi +y\cos \theta$ and 
$W=U_{P}\sin \phi$. An orientational wavenumber 
$\alpha _{1}$ is introduced, with geometric angle 
$\phi$ involved. This also suggests perturbation rearrangement by defining the complex phase velocity 
$c= \varOmega /\alpha _{1}$ (
$c={c}_{r}+\textrm {i}{c_{i}}$). We note that 
$\alpha _{1}$ differs from the resultant wavenumber vector 
$k$.
For completeness, a corresponding modified Squire equation is also derived. With the definition 
$\eta = \partial {u}/ \partial {z} - \partial {w}/\partial {x}$, the linearized 
$y$-direction vorticity equation gives
\begin{equation} {-}{\rm i}\varOmega{\eta}+({\rm i}\alpha{U}+{\rm i}\beta{W})\eta= {\rm i}(\alpha\times{D}W-\beta\times{D}U)v+\frac{1}{Re}\,({D}^2-k^2)\eta. \end{equation}Referring to (2.8a,b), a modified Squire equation is derived as
\begin{align} &\left[(U_{P}-c)+\frac{\alpha}{\alpha_{1}}\,y\cos\theta\right]\eta -\frac{1}{{\rm i}\alpha_{1}\,Re}\,({D}^2-k^{2})\eta\nonumber\\ &\quad= \left[\frac{1}{\alpha_{1}}\, (\alpha\sin\phi-\beta\cos\phi)\,{D}U_{P}-\frac{\beta}{\alpha_{1}}\cos\theta\right]v. \end{align}
Boundary conditions for the modified OS and Squire equations are 
$v(\pm {1})=\eta (\pm {1})=0$ and 
${D}v(\pm {1})=0$. In stability analysis, we focus on the sign of the imaginary part of the least stable mode, 
$c_{i,max}$. For Squire modes, it can be shown that
\begin{equation} c_{i,{Squire}} \int\nolimits_{{-}1}^{1} |{\eta}|^{2}\,{{\rm d}y}={-}\frac{1}{\alpha_{1}\,Re} \int_{{-}1}^{1} |D{\eta}|^{2}+k^{2}|{\eta}|^{2} \,{{\rm d}y} <0. \end{equation}The Squire modes are thus always damped, leading the disturbances to attenuate consistently.
The OS equation (2.7) with boundary conditions produces a generalized eigenvalue problem with five parameters: 
$\theta$, 
$\phi$, 
$Re$, 
$\alpha$ and 
$\beta$. Then for fixed 
$\theta$ and 
$\phi$, an isosurface 
$c_i=0$ can be determined in 
$(Re, \alpha,\beta )$ space. This is referred to as the neutral surface.
2.2. Numerical method
A spectral collocation method is adopted where the coefficient matrices are constructed in a Chebyshev domain using 
$v(y)=\sum _{n=0}^N a_{n}\,T_{n}(y)$ and 
$y_{j}=\cos (\kern0.07em{j{\rm \pi} }/{N})$, where 
$a_{n}$ is the coefficient of the 
$n$th Chebyshev polynomial 
$T_{n}(y)$, truncated at order 
$N$. Thus the generalized eigenvalue problem can be stated as 
$\boldsymbol{\mathsf{A}}\, \boldsymbol {a}=c \boldsymbol{\mathsf{B}}\, \boldsymbol {a}$, with 
$\boldsymbol{\mathsf{A}}$ and 
$\boldsymbol{\mathsf{B}}$ the coefficient matrices, and 
$\boldsymbol {a}$ the eigenvectors. The calculation is done using the QZ algorithm. The numbers of collocation points are chosen to guarantee that the maximum relative error of 
$c_i$, the imaginary part of the eigenvalues, for the five least stable modes is below 
$10^{-5}$, specifically at high Reynolds number and small wavenumbers.
Numerical verification is implemented with literature data on PCP flow stability with three parameters 
$(\theta ; Re, \alpha )$. Neutral curves in 
$(\theta ; Re, \alpha )$ space are shown in figure 2(a). Three cuts at constant 
$\theta =90^{\circ }$, 
$81.5^{\circ }$ and 
$73.3^{\circ }$ are also shown. The locus of minimum or critical Reynolds number 
$Re_{cr}(\theta )$ is shown using red triangles. The corresponding critical wavenumber is 
$\alpha _{cr}(\theta )$. Three neutral curves at fixed 
$\theta$ are plotted in the 
$Re\unicode{x2013} \alpha$ plane in figure 2(b). The present neutral curves for 
$\theta =90^{\circ }$, 
$81.5^{\circ }$ and 
$73.3^{\circ }$ agree well with data by Potter (Reference Potter1966), which are marked as symbols.
The PCP (
$\phi =0^\circ$) stability. (a) Neutral curves in 
$(\theta ;Re,\alpha )$ space. Coloured lines are neutral curves defined as surface cuts by planes of constant 
$\theta$. Red triangles indicate the trajectory of local critical points with 
$({Re}_{cr},\alpha _{cr})$ at each 
$\theta$. (b) Comparison of three extracted neutral curves at different 
$\theta$: 
$90^{\circ }$ (red dot-dashed line), 
$81.5^{\circ }$ (blue dashed line) and 
$73.3^{\circ }$ (black solid line). Lines indicate present results, and symbols indicate data from Potter (Reference Potter1966).
3. Global and local critical states
With 
$(\theta,\phi )$ fixed, we calculate 
$c_i(Re,\alpha,\beta )$, from which the neutral surface 
$f(\theta,\phi ;Re,\alpha,\beta )=0$ corresponding to the isosurface 
$c_i=0$ can be determined. The scheme for this is discussed below. Global critical parameters are the minimum 
$Re_{cr}=Re$ on this surface with corresponding 
$\alpha _{cr},\beta _{cr}$, and with stability for 
$Re< Re_{cr}$. An example is shown in figure 3(a) with 
$\phi =45^\circ$, 
$\theta =80^\circ$ where the critical state is marked with a red dot. If we fix 
$\chi >0^\circ$, then this defines a plane in 
$(Re,\alpha,\beta )$ space, 
$\alpha = \beta \cot \chi$. We can then define a local critical state as the minimal 
$Re_{cr,l}= Re$, together with corresponding wavenumbers 
$(\beta _{cr,l} \cot \chi,\beta _{cr,l})$ on the planar curve defined by the intersection of this plane with the neutral surface as defined above. A local critical state thus represents neutral stability of a GCP flow for a perturbation with specified direction 
$\chi$. One example of a local critical point with 
$\chi = {\rm \pi}/4$, marked as a square, is shown in black in figure 3(a).
The GCP stability illustration. (a) Example of a neutral surface in 
$(\alpha,\beta,Re)$ space at fixed 
$\phi =45^\circ$ and 
$\theta =80^\circ$. The red dot indicates the global critical point with 
$(\alpha _{cr},\beta _{cr},Re_{cr})=(0.041,1.016,5414.88)$. The intersection between a plane surface 
$\alpha =\beta \cot \chi$ with 
$\chi ={\rm \pi} /4$ and the neutral surface generates a local neutral curve (black) and a corresponding local critical point (black square) 
$( \alpha _{cr,l},{Re}_{cr,l})=(0.489,11674.89)$. (b) Mapping local critical Reynolds number 
$Re^*_{cr,l}$ in plot of 
$\theta ^*$ versus 
$Re_{cr,l}^*\alpha _{cr,l}/k_{cr,l}$. Values of 
$\cot \chi$: red, 0.1; black, 1; blue, 10. Values of 
$\phi$: square, 
$30^\circ$; circle, 
$45^\circ$; up triangle, 
$60^\circ$; down triangle, 
$90^\circ$. The dotted line indicates PCP results.
3.1. Mapping from GCP to PCP flow stability
Equation (2.7) can be reduced to a 3-D perturbed PCP OS equation (
$\phi =0^\circ$), given by
\begin{equation} \left[(U_{P}-{c})+y\cos\theta\right] ({{D}}^2-k^{2})v -{v}\times {{D}}^2{U_{P}}= \frac{1}{{\rm i}{\alpha}\,{Re}}\,({D}^2-k^{2})^{2}v, \end{equation}
where we still have 
$k = (\alpha ^2 + \beta ^2)^{1/2}$, and now 
$\alpha _1=\alpha$. This represents a PCP stability equation, but Squire's theorem, where one puts 
$\beta =0$, 
$k=\alpha$, has not been utilized.
To draw an analogy to the PCP equation (3.1), multiplying the modified OS equation (2.7) by 
$\sin \theta ^*/\sin \theta$, we can obtain
\begin{equation} \left[({U}^*-{c}^*)+ y\cos{\theta}^*\right] ({{D}}^2-k^{2})v -{v}\times {{D}}^2{{U}^*}= \frac{1}{{\rm i}\alpha \,{{Re}^*} }\, ({D}^2-k^{2})^{2}v, \end{equation}
where 
$\ast$ denotes mapping parameters defined as
\begin{align} \tan\theta^* =\frac{\alpha_{1}}{\alpha}\tan\theta,\quad {c}^* =c\,\frac{\sin{\theta}^*}{\sin\theta},\quad {U}^*={-}\frac{3}{2}\,(y^2-1)\sin{\theta}^* ,\quad {Re}^*= Re\,\frac{\alpha_1}{\alpha}\,\frac{ \sin \theta}{\sin{\theta}^*}. \end{align}
Comparing (3.1) to (3.2) shows that the latter can be interpreted as the OS equation for an equivalent PCP flow with 
${\theta } \to \theta ^*$ and 
${Re} \to Re^*$. Hence the full GCP stability equation can be mapped to that for a defined PCP stability.
For (3.2), for given 
${\theta }^*$, a neutral surface 
${f}^*({\theta }^*;{Re}^*,\alpha,\beta ) = 0$ can be determined. The equivalent PCP critical state 
$({Re}^*_{cr},\alpha _{cr})$ will be at the 
${Re}^*$ minimum of the curve formed by the intersection of this surface with the plane 
$\beta =0$. This is irrelevant to the GCP stability problem because Squire's theorem does not apply. So for given 
$(\theta,\phi )$, the global GCP critical point 
$(Re_{cr},\alpha _{cr},\beta _{cr})$ will map onto a point on the neutral surface 
$f^* =0$, but this will not generally coincide with the equivalent PCP critical state.
The above mapping admits definition of a local critical Reynolds number 
$Re^*_{cr,l}$, which is the mapped value of 
$Re_{cr,l}$. In other words, 
$Re^*_{cr,l}$ is the minimum 
$Re^*$ on the intersection of the plane 
$\alpha = \beta \cot \chi$ with the neutral surface 
$f^*(\theta ^*;Re^*,\alpha,\beta )=0$. So for given 
$\theta ^*$, it is the minimum 
$Re^*$ on the curve in the 
$Re^*\unicode{x2013} \beta$ plane 
$f^*(\theta ^*;Re^*,\beta \cot \chi,\beta ) =0$. This states that 
$Re^*_{cr,l} = G(\theta ^*,\chi )$, where 
$G$ is some function to be determined numerically. This can be simplified by writing 
$1/(\textrm {i}\alpha \,{{Re}^*})$ on the right-hand side of (3.2) as 
$1/(\textrm {i}k\,\widetilde {Re})$, where 
$\widetilde {Re} = Re^*\cos \chi$ (recall that 
$\cos \chi =\alpha /k$). If 
$\chi$ is fixed, corresponding to a local critical state, then (3.2) contains parameters 
$(\theta ^*;\widetilde {Re},k)$ only. It follows that the local critical Reynolds number must be of the form 
$\widetilde {Re}_{cr,l} = H(\theta ^*)$, or equivalently 
$Re^*_{cr,l}\alpha _{cr,l}/k_{cr,l} = H(\theta ^*)$.
In figure 3(b), calculated values of 
$Re^*_{cr,l} \alpha _{cr,l}/k_{cr,l}$ collapse well versus 
$\theta ^*$. The dotted line represents data obtained directly for PCP flow using (3.2). The symbols represent determination of local critical properties 
$(Re^*_{cr,l},\alpha _{cr,l},\beta _{cr,l})$ with 
$\cot \chi =0.1,1,10$, obtained using (2.7) for 
$\phi =30^\circ$, 
$45^\circ$, 
$60^\circ$ and 
$90^\circ$. These are then mapped to 
$^*$ parameters defined by (3.3a–d), where they are seen to collapse onto the PCP curve in figure 3(b). This collapse is expected for all 
$\phi$ and 
$\chi$, which has been further tested numerically for other 
$\chi$, but not shown here. Since the GCP global critical parameters 
$Re_{cr}(\theta,\phi ),\alpha _{cr}(\theta,\phi ),\beta _{cr}(\theta,\phi )$ always correspond to some particular set of local critical parameters, their mapped versions also must lie somewhere on the curve of figure 3(b). It follows that all global critical parameters for GCP stability must lie on this curve.
3.2. Calculation of critical parameters
The above mapping further indicates a method for determination of global critical states for GCP flow. Starting with (3.2), a neutral surface 
$f^*(\theta ^*;Re^*,\alpha,\beta )=0$ can be determined. In 
$(Re^*,\alpha,\beta )$ space, there is a different neutral surface for each 
$\theta ^*$. For given 
$(\theta,\phi ;\alpha,\beta )$, 
$\theta ^*$ can be calculated from (3.3a–d). Since 
$\alpha,\beta$ are specified, 
$Re^*$ can be calculated. With known 
$Re^*$, 
$Re$ can be calculated from the last of (3.3a–d). An algorithm can be summarized as follows.
(i) Specify
$(\theta,\phi )$.(ii) For each
$(\alpha,\beta )$ over some range:
(iii) Over the range of
$(\alpha,\beta )$ chosen, find the minimum $Re$, which is $Re_{cr}$. This will also give the corresponding $(\alpha _{cr},\beta _{cr})$. That is the red point in figure 3(a).(iv) Over the ranges of
$(\theta,\phi )$, this will give $Re_{cr}$, $\alpha _{cr}$ and $\beta _{cr}$ as functions of $(\theta,\phi )$.
4. Global critical parameters for GCP flow instability
Here we discuss global GCP critical state properties. Figure 4 shows 3-D plots of 
$Re_{cr}(\theta,\phi )$, 
$\alpha _{cr}(\theta,\phi )$ and 
$\beta _{cr}(\theta,\phi )$. Coloured isosurfaces are constructed with data at intervals of 
$0.1^\circ$ in both 
$\theta$ and 
$\phi$, with symbols at intervals of 
$10\,^\circ$ overlayed for clarity. As 
$\theta \to 0^\circ$, for PC flow, 
$Re_{cr}\to \infty$. As both 
$\theta$ and 
$\phi$ increase, 
$Re_{cr}$ decreases monotonically, displaying 
$Re_{cr} \to 3848$ at 
$\theta \to 90^\circ$ for pure PP flow. For PCP flow with 
$\phi =0^\circ$, there is singular behaviour with 
$Re_{cr}\to \infty$ for 
$\theta < 62.17^\circ$, with finite 
$Re_{cr}$ for 
$\theta >62.17^\circ$. For all 
$\phi >0^\circ$, 
$Re_{cr}$ is finite except in the PC flow limit 
$\theta =0^\circ$.
Global critical properties in a 3-D plot: (a) global critical Reynolds number 
${Re}_{cr}$; (b) global critical wavenumbers 
$\alpha _{cr}$ and 
$\beta _{cr}$. Surface composed of data in the 
$\phi \unicode{x2013} \theta$ plane with intervals of 
$0.1^\circ$. Dots for data with intervals of 
$10^\circ$ show the general behaviour.
Isosurface plots of corresponding 
$\alpha _{cr}(\theta,\phi )$ and 
$\beta _{cr}(\theta,\phi )$ are shown in figure 4(b). At 
$\phi =0^\circ$ for PCP flow, 
$\beta _{cr}=0$ for 
$\theta >62.17^\circ$ where 
$Re_{cr}$ is finite, indicating zero 
$z$-direction perturbation. At fixed 
$\phi$ when 
$\theta$ decreases from 
$\theta = 90^\circ$, 
$\beta _{cr}$ increases almost discontinuously to near its maximum with a corresponding decrease of 
$\alpha _{cr}$. For the region where 
$\theta \lesssim 60^\circ$, 
$\alpha _{cr}$ is generally small, and finite 
$\beta _{cr}$ dominates critical perturbations. Both maximum values of 
$\alpha _{cr}$ and 
$\beta _{cr}$ are approximately 
$1.02$, the value for PP flow instability (Orszag Reference Orszag1971). Figure 4(b) clearly shows a region of the 
$\theta \unicode{x2013} \phi$ plane where both 
$\alpha _{cr}$ and 
$\beta _{cr}$ display rapidly varying, non-monotonic and even discontinuous features. This is now discussed.
4.1. Scaling of $Re_{cr}$
In figure 5(a), we show the global 
$Re_{cr}$ versus 
$\theta$. Nine cases with different 
$\phi$ are shown, which are divided into three groups. The first set includes 
$\phi =0^\circ, 5^\circ,15^\circ$, relevant to discussion of the deviation behaviour from PCP flow stability. The second set considers 
$\theta =19^\circ, 20^\circ, 20.5^\circ$, with focus on rapid changes in this medium 
$\phi$ region. The third set comprises 
$\theta =30^\circ, 60^\circ, 90^\circ$, useful for discussing general trends and possible scaling. The PCP curve (
$\phi =0^\circ$), as noted previously, shows rapidly increasing 
$Re_{cr}$ when 
$\theta \to \theta _c \approx 62.17^\circ$. For GCP flow, 
$Re_{cr}$ deviates substantially from the PCP curve, indicating non-differentiable variation at 
$\theta \approx 60^\circ$ for 
$\phi \lesssim 20^\circ$. The subregion where both GCP and PCP results collapse shrinks with increasing 
$\phi$. For 
$\phi \gtrsim 20^\circ$, 
$Re_{cr}$ for GCP stability generally diverges from the PCP data trend.
(a) Plot of 
$Re_{cr}$ versus 
$\theta$ for different 
$\phi$. (b) Scaling using 
$Re_{cr}^*(\alpha _{cr}/k_{cr})$ for different 
$\phi$ with 
${Re}^*= Re\,({\alpha _1}/{\alpha })({ \sin \theta }/{\sin {\theta }^*})$. Dashed lines indicate discontinuities with no data points. For 
$\theta <55^{\circ }$, data collapse on 
$3848$. The dotted line indicates PCP results. Values of 
$\phi$: square, 
$0^\circ$; circle, 
$5^\circ$; up triangle, 
$15^\circ$; down triangle, 
$19^\circ$; diamond, 
$20^\circ$; left triangle, 
$20.5^\circ$; right triangle, 
$30^\circ$; 
$\times$, 
$60^\circ$; 
$*$, 
$90^\circ$.
Following our previous hypothesis, the mapped 
$Re_{cr}^{*}$ in (3.3a–d) is employed using critical wavenumbers 
$\alpha _{cr}$ and 
$\beta _{cr}$. In order to avoid sparsely scattered values of 
$Re_{cr}^{*}$ in the range 
$10^3$ to 
$10^{20}$, we instead plot 
$Re_{cr}^* \alpha _{cr}/k_{cr}$ versus 
$\theta$ as shown in figure 5(b), where the region of 
$\theta \lesssim 55^\circ$ is not shown owing to the reasonable collapse with 
$Re_{cr}^* \approx 3848$, corresponding to literature data at 
$Re_{cr}=5772$ via the mapping transformation for PCP flow. Two common, interesting tendencies are apparent in figure 5(b). One is the sharp decrease of 
$Re_{cr}$ that takes place over a relatively small 
$\theta$ range; as 
$\phi$ increases, there is a monotonic decrease in the value of 
$\theta$ at which this sharp decrease in 
$Re_{cr}^*$ occurs. The second tendency is the gradual increase of 
$Re_{cr}^*$ for decreasing 
$\theta$ beginning at 
$\theta =90^\circ$. This increase is strongest for PCP flow at 
$\phi =0^\circ$, and then weakens with increasing 
$\phi$, with the apparent local maximum occurring at increasing 
$\theta$. At 
$\phi \approx 19^\circ$ in the figure, several discontinuities of 
$Re_{cr}$ are apparent, which can be ascribed to competition between 
$x$- and 
$z$-direction wavenumbers 
$\alpha$ and 
$\beta$, respectively.
With further scaling of 
$\theta ^*$, all 
$Re_{cr}^* \alpha _{cr}/k_{cr}\unicode{x2013} \theta ^*$ will collapse onto the PCP results, which is guaranteed for local critical Reynolds numbers with fixed 
$\beta /\alpha$, as demonstrated in figure 3. This collapse of local properties is not useful, however, for identifying behaviour of the global 
$Re_{cr}\unicode{x2013} \theta$ relationship since the GCP–PCP transformation is heavily distorted in parameter space, with multi-valued subdomains, and thus does not provide a clear one-to-one mapping for global critical parameters.
4.2. Critical wavenumbers $\alpha _{cr}$ and $\beta _{cr}$
The complex variation of 
$Re_{cr}^*$ implies an intrinsic competition between 
$\alpha _{cr}$ and 
$\beta _{cr}$. In figure 6, we focus on the flow with 
$\phi =20^\circ$, discussing the competing behaviour of two critical wavenumbers, with neutral surface illustrations at several 
$\theta$. In figure 6(a), we show the 
$\alpha _{cr}$, 
$\beta _{cr}$ and 
$Re_{cr}$ variation with 
$\theta$. It is apparent that 
$Re_{cr}$ is always continuous but with discontinuous derivatives at certain 
$\theta$ where both 
$\alpha _{cr}$ and 
$\beta _{cr}$ show discontinuities. With decreasing 
$\theta$, 
$\alpha _{cr}$ decreases gradually with corresponding increase of 
$\beta _{cr}$. At 
$\theta$ slightly smaller than 
$85^\circ$, 
$\alpha _{cr}$ drops discontinuously to a smaller value, while 
$\beta _{cr}$ also jumps discontinuously. With further decrease of 
$\theta$ to less than 
$77^\circ$, a rebound of 
$\alpha _{cr}$ and a corresponding fall of 
$\beta _{cr}$ can be seen. Further decrease in 
$\theta$ gives a gradual decrease of 
$\alpha _{cr}$. At 
$\theta \lesssim 70^\circ$, another discontinuity appears, with 
$\alpha _{cr}$ recovering to a small value, while 
$\beta _{cr}$ returns to almost its maximum.
General behaviour of global critical parameters for 
$\phi =20^\circ$: (a) 
$\alpha _{cr}$ (dot-dashed line), 
$\beta _{cr}$ (dashed line) and 
$Re_{cr}$ (solid line). Cases at four 
$\theta$ are represented by dotted lines with corresponding data shown in: (b), point ‘b’, 
$\theta =85^\circ$; (c), point ‘c’, 
$\theta =80^\circ$; (d), point ‘d’, 
$\theta =70^\circ$; (e), point ‘e’, 
$\theta =65^\circ$. Blue for the neutral surface, grey plane for 
$\alpha /\beta =\alpha _{cr}/ \beta _{cr}$, intersection line for neutral curve, red dot for corresponding global critical state.
To depict these three discontinuities, we provide four neutral surface plots. The first is at 
$\theta =85^\circ$ in figure 6(b), where the 3-D neutral surface is semi-transparent blue and the grey plane corresponds to the 
$\alpha _{cr}\unicode{x2013} \beta _{cr}$ plane, the global critical state. The intersection of the two surfaces gives the global Reynolds number 
$Re_{cr}$. On the neutral surface, the local Reynolds number varies slowly with increasing 
$\theta$. At the global critical point, 
$(\alpha _{cr},\beta _{cr})$ is at a comparable scale, which indicates that the critical stability is determined by contributions from both PCP flow in the 
$x$-direction and PP flow in the 
$z$-direction.
In the second case, 
$\theta =80^\circ$ in figure 6(c), the middle region of the neutral surface shifts to higher Reynolds number, and two small 
$Re$ regions are shown, while 
$Re_{cr}$ locates on the smaller 
$\alpha _{cr}$ valley (in the sense of a local minimum 
$Re$) or lobe. This is consistent with the first discontinuity noted above. Similarly, in all cases of 
$\theta \in [77.5^\circ,84^\circ ]$, the critical stability is dominated by PP flow in the 
$z$-direction. With 
$\theta$ decreasing to 
$70\,^\circ$ in figure 6(d), both lobes of the neutral surface contract, while 
$Re_{cr}$ jumps onto the relatively larger 
$\alpha$ lobe, producing the second discontinuity in critical wavenumbers. For this case, the dominant part of the flow is PCP flow in the 
$x$-direction. The flow with 
$\theta =65^\circ$, in figure 6(e), shows that 
$Re_{cr}$ jumps back to the small 
$\alpha$ lobe, representing the third and final discontinuity with varying 
$\theta$. Through these rapid changes of flow patterns, the wavenumber 
$\alpha _{cr}$ remains positive but small. The critical stability is again dominated by 
$z$-direction perturbed Poiseuille flow, and this pattern remains as 
$\theta$ further decreases to zero. This interpretation for 
$\phi =20^\circ$ clarifies that the discontinuous behaviour of critical wavenumbers is compatible with the idea of wavenumber jumping/competition.
Critical wavenumber behaviour at different 
$\phi$ can vary from that discussed above. Variation of 
$\beta _{cr}\unicode{x2013} \theta$ is shown in figure 7(a), where nine cases are again adopted. Except for 
$\phi =0^\circ$, all other flows show non-zero values. For 
$\phi =5^\circ$, 
$\beta _{cr}$ shows a mild peak around 
$\theta =86^\circ$, then maintains a relatively small-value range until a sharp increase to approximately unity at 
$\theta \approx 60^\circ$. With increasing 
$\phi$, the increase discontinuity in 
$\theta$ moves towards 
$\theta =90^\circ$, and the small peak migrates towards lower 
$\theta$. The case 
$\phi =15^\circ$ shows behaviour similar to that for 
$\phi =5^\circ$. For 
$\phi =19^\circ$, a second peak region is seen, which adds two additional discontinuities. With further increase in 
$\phi$, the first discontinuity at 
$\theta \approx 84^\circ$ tends to be small, while the other two discontinuities become closer in 
$\theta$. The two discontinuities disappear quickly with further increase in 
$\phi$. With 
$\phi =30^\circ, 60^\circ, 90^\circ$, complex discontinuous behaviour in 
$\beta _{cr}(\theta,\phi )$ is not observed, and variation with 
$\theta$ is smooth and monotonic.
Global critical wavenumbers: (a) 
$\beta _{cr} (\phi,\theta )$, (b) 
$\alpha _{cr} (\phi,\theta )$. Lines with symbols represent different 
$\phi$, referring to figure 5. Vertical dashed lines represent discontinuous jumps that accompany critical state jumps between lobes of the neutral surface.
Figure 7(b) shows 
$\alpha _{cr}$ with corresponding discontinuities. The case 
$\phi =5^\circ$ is quite close to 
$\phi =0^\circ$ except for 
$\theta \lesssim 65^\circ$. This shows that for small 
$\phi$, although 
$\beta _{cr}$ is not zero for any 
$\theta _{cr}$, stability is still 
$\alpha _{cr}$-dominant, gradually deviating from PCP stability.
5. Conclusion
Analysis of the neutral surfaces/curves and critical state properties of the linear, modal stability of generalized Couette–Poiseuille (GCP) flow has been explored. A GCP flow is generated by the shear provided by the velocity-difference vector of two relatively moving parallel plates defining an 
$x$-direction shear flow and a volume flow produced by an externally applied pressure gradient vector inclined at an angle 
$\phi$ to the 
$x$-axis. The base flow consists of a superposition of plane Couette (PC) and plane Poiseuille flows that gives a 3-D velocity field where the mean velocity is always parallel to the channel walls but rotates such that both its magnitude and flow angle are functions of the wall-normal coordinate between plates. A parameter set 
$(Re, \theta, \phi )$ in spherical coordinates, with radius a Reynolds number 
$Re$, polar angle 
$\theta$ and azimuthal angle 
$\phi$, defines a particular GCP flow, with 
$\tan \theta$ describing the relative proportion of Couette-type and Poiseuille-type contributions. Both 
$x$-direction (wavenumber 
$\alpha$) and 
$z$-direction (wavenumber 
$\beta$) perturbations are discussed. The GCP class of flows is unusual in that Squire's theorem does not apply. For the special case of plane Couette–Poiseuille (PCP) flow (
$\phi =0^\circ$), it defines a three parameter 
$(\theta ;Re,\alpha )$ stability problem. The critical coefficients extracted from the present 3-D PCP results agree well with corresponding data of Potter (Reference Potter1966).
Properties of the neutral surface 
$f(\theta,\phi ;Re,\alpha,\beta )=0$ for GCP flow, defined in 
$(Re,\alpha,\beta )$ space for 
$(\theta,\phi )$ fixed, are discussed. A mapping from GCP to an equivalent PCP stability is found, from which local critical scaling relations 
$\theta ^*$ versus 
$Re^*_{cr,l}\,\alpha _{cr,l}/k_{cr,l}$, where 
${Re}^*= Re\,({\alpha _1}/{\alpha })\,({\sin \theta }/{\sin {\theta }^*})$ and 
$\tan {\theta }^*=({\alpha _{1}}/\alpha )\tan {\theta }$ with 
$\alpha _1=\alpha \cos \phi +\beta \sin \phi$, collapse all data for a fixed perturbation wavenumber direction, which is equivalent to a linear relation between local orthogonal scalar wavenumbers. This demonstrates a general failure of Squire's theorem. The global critical 
$Re_{cr}$ shows continuous variation in 
$(\theta,\phi )$ space with certain singular limits that includes PC flow, while 
$\alpha _{cr}$ and 
$\beta _{cr}$ have complex behaviour with discontinuities produced by jumps across different lobes of the neutral surface. In particular, the example illustrated with 
$\phi =20^\circ$ displays three discontinuities in 
$\theta$ of a type not seen for PCP flow stability. This discontinuous behaviour appears for all 
$\phi \lesssim 21^\circ$, while for larger values, both 
$\alpha _{cr}(\theta )$ and 
$\beta _{cr}(\theta )$ are generally smooth and monotonic.
The present study of GCP flow provides a useful reference for stability analysis of two-dimensional perturbations imposed on a base velocity profile with wall-normal twist. An interesting and open question is whether the GCP 
$\to$ PCP mapping described here is applicable to wider classes of 3-D mean velocity profiles, perhaps suggesting further novel, discontinuous neutral surface dynamics with corresponding failure of Squire's theorem and complex critical-property behaviour.
Funding
This work was supported by the National Natural Science Foundation of China (nos 12172352, 12388101), the Strategic Priority Research Program of the Chinese Academy of Sciences (no. XDB0500300) and the National Key R&D Program of China (no. 2022YFA1005200).
Declaration of interests
The authors report no conflict of interest.
