Stability of obliquely driven cavity flow

Abstract The linear stability of the incompressible flow in an infinitely extended cavity with rectangular cross-section is investigated numerically. The basic flow is driven by a lid which moves tangentially, but at yaw with respect to the edges of the cavity. As a result, the basic flow is a superposition of the classical recirculating two-dimensional lid-driven cavity flow orthogonal to a wall-bounded Couette flow. Critical Reynolds numbers computed by linear stability analysis are found to be significantly smaller than data previously reported in the literature. This finding is confirmed by independent nonlinear three-dimensional simulations. The critical Reynolds number as a function of the yaw angle is discussed for representative aspect ratios. Different instability modes are found. Independent of the yaw angle, the dominant instability mechanism is based on the local lift-up process, i.e. by the amplification of streamwise perturbations by advection of basic flow momentum perpendicular to the sheared basic flow. For small yaw angles, the instability is centrifugal, similar as for the classical lid-driven cavity. As the spanwise component of the lid velocity becomes dominant, the vortex structures of the critical mode become elongated in the direction of the bounded Couette flow with the lift-up process becoming even more important. In this case the instability is made possible by the residual recirculating part of the basic flow providing a feedback mechanism between the streamwise vortices and the streamwise velocity perturbations (streaks) they promote. In the limit when the basic flow approaches bounded Couette flow the critical Reynolds number increases very strongly.

shallow (Γ = 0.5) and a deep cavity (Γ = 2) are presented in § 4 and the connection of the critical curves along α is demonstrated. Flow structures and instability mechanisms are investigated by considering the local and global production rates of kinetic perturbation energy. Finally, limiting cases of the yaw angle and common properties of the instabilities found are discussed in § 5.

Formulation of the problem
We consider the incompressible flow of a Newtonian fluid with density ρ and kinematic viscosity ν in a rectangular cavity (figure 1). The depth d in the y-direction and the width h in the x-direction define the aspect ratio Γ = d/h. In the z-direction the cavity is assumed to be infinitely extended. The origin of the coordinate system is placed in the centre of the (x, y) cross-section. The flow is driven by the steady tangential motion of a lid at the top y = d/2 of the cavity. The lid velocity vector U = U(cos αe x + sin αe z ), where e x and e z are the unit vectors in the x and z directions, respectively, is inclined with respect to the x axis with inclination (or yaw) angle α.
(2.2c) Furthermore, we consider the case of a vanishing pressure gradient in the z-direction, ∂p/∂z = 0. The problem is thus defined by three parameters: the aspect ratio Γ , the inclination angle α and the Reynolds number Due to the translation invariance of the problem in z and t the governing equations allow for a steady two-dimensional basic flow u 0 (x, y) which only depends on x and y.
We are interested in the linear stability boundary, expressed by Re c (Γ, α), at which the two-dimensional basic flow becomes unstable to three-dimensional perturbations. Linearising (2.1) with respect to small perturbations of the basic flow yields the linear perturbation equations ∂u ∂t + u 0 · ∇u + u · ∇u 0 = −∇p + ∇ 2 u, (2.4a) where now u and p denote the deviations from the basic state.

Methods of solution
All differential equations are discretized with triangular elements on a rectangular domain (x, y) using the finite element library FEniCS (Alnaes et al. 2015). To properly resolve the flow fields near the boundaries the mesh is refined towards all walls by subsequently doubling the number of grid points within 5 %, 1 % and 0.5 % of the width and the depth of the cavity. Taylor-Hood elements are employed which implement a quadratic interpolation for the velocity fields and a linear interpolation for the pressure. and wave number k c as functions of the grid resolution for α = 0 • . The column labelled 'Grid' refers to the initial grid size before refinement, n T denotes the number of triangles, while n DOF is the number of degrees of freedom.
3.1. Basic state The steady two-dimensional flow (u 0 , p 0 ) is computed using Newton-Raphson iteration already implemented in the FEniCS framework, which only requires the variational formulation and the boundary conditions. Absolute and relative convergence criteria based on the L 2 norm of the residuum are set to 10 −10 and 10 −8 , respectively. During tracking of the stability boundary, the basic state calculation is typically terminated due to the absolute convergence criterion. The converged basic flow field enters the linear stability analysis parametrically.

Linear stability analysis
Once the basic state (u 0 , p 0 ) is computed the linear stability equations (2.6) are solved on the finite element mesh for given α, Re and k using an implicitly restarted Arnoldi method implemented in ARPACK (Lehoucq & Salinger 2001) and called within the subroutine eigs of SciPy. In order to ensure that the method captures all the eigenvalues of interest the dimension of the Krylov space is set to 300, while the number of converged eigenvalues required to assume convergence is set to 50. We noticed that when lowering both these numbers some eigenvalues might not be captured correctly.
To find the neutral curves a sensitivity-based algorithm has been developed which is detailed in Appendix A. The required first-order sensitivity of the eigenvalues with respect to wavelength variations is derived in Appendix B.

Code verification
In a first step a grid-convergence study for the critical Reynolds and wave numbers is carried out. Table 1 shows (Re c , k c ) as a function of the grid resolution for three aspect ratios Γ and α = 0 • . Grid convergence is clearly obtained and the converged results compare very well, i.e. within 1 %, with the reference results of Albensoeder et al. (2001). The data suggest that a basic mesh of 40 × 40 provides already very accurate results for Γ = 1. Note the grid specified represents the grid with equidistant spacing, the actual grid used is refined towards the walls as specified above such that the formal 40 × 40 resolution practically is made of 13 976 elements or 95 328 degrees of freedom. With similar arguments, the initial grids 40 × 80 and 80 × 40 for Γ = 2 and Γ = 1/2, respectively, are used for the stability analysis for α > 0 • . To verify the growth rate σ and the oscillation frequency ω as functions of the wave number k, we consider Γ = 1 and α = 0 • for which reference data are available in the literature. To that end, the most dangerous mode has been computed for Re = 200 and Re = 1000. Figure 2 shows the growth rates and oscillation frequencies of the fastest growing mode for Re = 200 (dashed lines) and Re = 1000 (full lines) in comparison with the results of Ding & Kawahara (1999) ( ) and Albensoeder et al. (2001) (•). For Re = 200, an excellent agreement is found for all k considered, using the basic grid resolution of 40 × 40. The numerical results for Re = 1000 also show a good agreement with the reference data for the frequency ω. Agreement of the growth rate σ obtained for the current resolution with the reference data is acceptable. Typically, our results are in-between the two reference data sets and tend to compare slightly better with those of Ding & Kawahara (1999) than with those of Albensoeder et al. (2001).

Energy analysis
In order to understand the fundamental instability mechanisms it is helpful to evaluate the budget of the kinetic energy of the critical mode u. Taking into account the perturbation flow vanishes on the boundaries and the nonlinear term u · ∇u is energy-preserving, the Reynolds-Orr equation can be written as (Albensoeder et al. 2001) where E kin = V (u 2 /2) dV and V is the volume occupied by the fluid over one period in z. The dissipation D * = V (∇u) 2 dV 0 is positive and always contributes to a reduction of the kinetic energy. If the total energy production rate I = − V u · (u · ∇u 0 ) overcomes the dissipation rate, the perturbation kinetic energy grows and the basic flow is unstable.
It is useful to decompose the perturbation flow u into components parallel and perpendicular to the basic flow (Albensoeder et al. 2001) and define Inserting this decomposition in (3.1) and normalising all terms with D * yields the following local dissipation and production terms: In this formulation the Reynolds-Orr equation reads as where I n = V i n dV. The local and the total energy production rates are i = n i n and I = n I n , respectively. The four local energy production terms i n describe the rate of change of the kinetic energy density due to the transport of basic state momentum u 0 by the perturbation flow u either perpendicular (e ⊥ · ∇) or parallel (e · ∇) to the direction of the basic flow, feeding to the (perpendicular or parallel) perturbation flow itself. These advective transport mechanisms build on the local shear (i 1 , i 2 ) or the local deceleration of the basic flow (i 3 , i 4 ).
Symmetries restrict the energy production terms. For instance, a local acceleration of the basic flow with e · (e · ∇u 0 ) > 0 cannot locally produce kinetic perturbation energy, because this condition renders i 4 < 0. On the other hand, a flow deceleration (e · (e · ∇u 0 ) < 0) can locally increase the kinetic energy by the process i 4 . Furthermore, it is easy to see that i 1 = 0 for unidirectional basic flows, and for parallel plane shear flows also, i 3 = i 4 = 0. Therefore, the energy production term i 2 is the dominant energy production term in most shear-dominated systems.
The term i 2 plays a major role for the centrifugal instabilities of the basic vortex flow in the lid-driven cavity (Albensoeder et al. 2001). In this system the streamlines of the basic flow are locally curved and the momentum of the basic flow decreases radially outward from the vortex core due to the stationary rigid walls. The process i 2 is also important in plane shear flows in which streamwise perturbation vortices (u ⊥ ) can extract considerable energy from the basic flow u 0 and feed this energy to the streamwise velocity perturbation (u ) in form of streaks leading to a considerable transient energy growth (Butler & Farrell 1992). Today the process i 2 is called the lift-up mechanism. The terminology originates from the observation that streamwise vortices seem to lift-up slow-speed streaks away from the wall just before a burst event occurs, initiating the transition to turbulence in a boundary layer flow (Kline et al. 1967;Landahl 1975;Brandt 2014). The energy production term i 2 also includes the Orr mechanism (Farrell & Ioannou 1993;Jiao, Hwang & Chernyshenko 2021), or shear stress mechanism (Butler & Farrell 1992), which essentially describes the shearing of spanwise vortices by the basic flow.
The physical transport mechanisms associated with i 2 are independent of the particular flow system. Here we shall address i 2 as the lift-up term, because it turns out that the critical modes typically arise as streamwise vortices and the contribution of the Orr mechanism is of secondary importance. Perturbations in the form of streamwise vortices extended along the basic flow direction can be identified by u ⊥ . They can efficiently extract energy from the basic shear flow and transfer momentum, hence energy, to the streamwise perturbation flow u via the process i 2 . An example are the curved Görtler vortices in the lid-driven square cavity (Albensoeder et al. 2001). Typically, the energy transfer occurs in the region between counter-rotating streamwise vortices. Therefore, regions (isosurfaces) of large i 2 arise as elongated structures just as the streamwise vortices. In a situation in which the lift-up process i 2 dominates the energy budget of the perturbations the isosurfaces of i 2 are very similar to the isosurfaces of u . In such case, typical for the present investigation, isosurfaces of i 2 can safely be identified as curved streaks.

λ 2 -criterion
In order to detect and visualise vortices in the perturbation flow associated with the critical eigenmodes we use the λ 2 -criterion, introduced by Jeong & Hussain (1995). To that end, the perturbation velocity gradient is decomposed into a symmetric and an anti-symmetric part, respectively, The vortex core is then defined as the connected region in which two of the real eigenvalues of S 2 + 2 are negative. If the eigenvalues (λ 1 , λ 2 , λ 3 )(x) are ordered by size, λ 2 < 0 should be negative in the vortex core. A vortex is then identified as a connected region within which λ 2 < 0 and a vortex core can be visualised by displaying isosurfaces of constant λ 2 < 0.
3.6. Nonlinear numerical simulation For the purpose of an additional verification and for a clarification of the bifurcation character being sub-or supercritical, we also carried out full numerical simulations of the time-dependent three-dimensional flow. To that end, the problem (2.1) was solved by employing the spectral element code NEK5000 (Fischer, Lottes & Kerkemeier 2008).
For these calculations, the flow was assumed periodic in z with a wavelength corresponding to 2π/k c . Using a regular tensor mesh composed of N x × N y × N z = 20 × 20 × 10 elements of polynomial order p = 6 for the velocity and p = 4 for the pressure, simulations are carried out for Γ = 1 with α = 22.5 • . Temporal integration was performed using a third-order Adams-Bashforth scheme with third-order extrapolation of the convective terms.

Basic flow
The basic flow u 0 = u 2-D 0 + u C 0 is the two-dimensional steady solution of (2.1) and (2.2). It can be decomposed into a recirculating two-dimensional cavity flow u 2-D 0 (x, y) driven by the effective Reynolds number Re 2-D = Re cos α, and the parallel bounded Couette flow u C 0 (x, y) = w C 0 (x, y)e z driven by the effective Reynolds number Re C = Re sin α. The recirculating part u 2-D 0 of the flow field is independent of the Couette part of the flow, because ∇ · u 0 = ∇ · u 2-D 0 = 0 and the nonlinear coupling terms u C 0 · ∇u 2-D 0 = 0 vanishes. On the other hand, the parallel Couette part of the flow u C 0 depends on the recirculating part of the flow and results from a linear equation balancing viscous diffusion and advection by u 2-D 0 in the (x, y) plane. The strength of both parts of the flow are related to each other via the Reynolds number and the inclination angle.
In the combined basic flow u 0 fluid elements have helical trajectories. This flow structure also arises in the context of air motion in street canyons driven by oblique wind directions (see e.g. Soulhac, Perkins & Salizzoni 2008;Zajic et al. 2011). The projections of the fluid trajectories onto the (x, y) plane correspond to the closed streamlines of the recirculating part u 2-D 0 of the flow. The pitch of the fluid trajectories is determined by the spanwise component u C 0 . Owing to the strong maximum principle for linear elliptic partial differential equations (see e.g. Evans 2010), the spanwise velocity u C 0 of a fluid element (and also its mean) is always less than the span component Re sin α of the lid velocity. Therefore, the spanwise velocity is considerably stronger near the moving lid than in the bulk of the cavity, and fluid elements are transported in the z-direction mainly in the upper part of the cavity.
In the limit α → 0 the classical lid-driven cavity flow is recovered with u C 0 = 0. The stability boundary has been investigated by several authors with Albensoeder et al. (2001) perhaps providing the most comprehensive stability results quasi-continuously covering the range of aspect ratios Γ ∈ [0.2, 4]. In the other limit of α → π/2, the recirculating part u 2-D 0 = 0 vanishes and the basic flow arises as a pure bounded Couette flow in a rectangular channel which can be written in form of an infinite series (4.1) The stability of this basic flow has been considered by Theofilis et al. (2004). No unstable modes have been found by these authors, even for Reynolds numbers as large as Re = 5000. Our linear stability analysis also indicates the basic flow is linearly stable, at least up to Re = 3000. Since the critical Reynolds numbers for lid-driven cavity flows for α = 0 • and Γ 0.5 satisfy Re c < 10 3 (Albensoeder et al. 2001), a strong stabilisation of the basic flow is expected as α → π/2. The stability boundary Re c (α, Γ ) for intermediate parameter values depends on the inclination angle α and on the cross-sectional aspect ratio Γ . Therefore, calculations have been carried out for selected aspect ratios, varying α quasi-continuously, and for representative yaw angles, varying Γ .     Table 2. Critical data and integral production rates of kinetic energy as functions of the aspect ratio Γ and the inclination angle α of the lid velocity vector.
to qualitatively different critical modes, depending on α. As α approaches π/2, the critical wave number becomes very small (figure 3b), indicating the critical mode becomes nearly two-dimensional. Numerical data for the critical parameters are listed in table 2 for several representative yaw angles α.

Modes I and II
At α = 0 • the classical Taylor-Görtler mode (mode I, Albensoeder et al. 2001) with relatively high wave number is recovered. As the inclination angle increases from zero, the Taylor-Görtler mode I with a small wavelength evolves continuously and changes only slightly due to the Couette part of the basic flow. While the Taylor-Görtler mode I is stationary for α = 0 • , the Görtler vortices drift in the positive z-direction with a phase velocity which increases as α increases. From figure 3(b,c), it can be observed that the phase velocity of the Görtler mode increases nearly linearly with α. When α is increased, the basic flow is slightly stabilised until, at α ≈ 4.3 • , the critical mode I (blue) changes to mode II (grey) which has a similar wave number. Mode II very much resembles mode I and the corresponding neutral stability boundary extends down to α = 0 • (not shown). At α = 0 • mode II is only the second most dangerous mode and, to the best of our knowledge, it has not yet been reported in the literature.
The neutral mode II is illustrated in figure 4 for α = 6 • . Shown is the perturbation velocity field u in the plane y = 0 (figure 4a) and in a plane z = const. with respect to the streamlines of u 2-D 0 . Finally, figure 4(c) shows a three-dimensional view over two wavelengths of the isosurfaces of the energy production rate i at 10 % of its maximum value i max . The region in which i > 0.1 × i max is also indicated by the blue areas in figure 4(b). Comparing figure 4(c) with figure 4(a) it is clear that the banana-shaped regions of high energy transfer are reflecting the perturbation vortex structures on which the energy transfer relies. The perturbation vortices are located just in-between neighbouring isosurfaces of i shown in figure 4(c). Since the energy budget is dominated by I 2 , the isosurfaces shown in figure 4(c) very well approximate isosurfaces of alternating streaks, i.e. of u , produced by the Taylor-Görtler vortices (u ⊥ ). This correspondence is demonstrated further below in figure 6(c).
Mode II very much resembles the classical Taylor-Görtler mode I for α = 0 • (Albensoeder et al. 2001) with strong vortices on the wall at x = −1/2, upstream of the moving lid, where the energy production peaks. On the downstream wall at x = 1/2 the vortices are much weaker. The weak vortices on the downstream wall are slightly offset in the positive z-direction as compared with the vortices on the upstream wall (figure 4a), as a result of the Couette part of the basic flow. Hence, the vortices tend to be slightly spiral with a small pitch. For α < 5.9 • , this mode propagates in the negative spanwise direction as can be seen from figure 3(c). This means that for the small range of α for which this mode is the critical mode, it propagates against the z-component of the lid velocity. Progressively, the phase velocity diminishes and changes sign such that, for α > 5.9, the wave propagates in the same z-direction as the lid. Again, the propagation speed scales approximately linearly with the yaw angle α.

Mode III
Near α = 5.8 • the critical mode II (grey) changes to mode III (orange in figure 3) which has approximately half the wave number as the low-α Taylor-Görtler modes. It can be seen from figure 3 that the neutral mode III originates from α = 0 • and has already been reported (Ding & Kawahara 1999;Albensoeder et al. 2001). Unlike mode I, which is stationary at α = 0 • , mode III is oscillatory at α = 0 • and arises as a pair of waves travelling in the ±z directions. As α increases from zero, the degeneracy of the neutral Reynolds number is removed and the basic flow is strongly destabilised with respect to the mode which propagates in the negative z-direction (opposite to the direction of the Couette part of the basic flow), while the basic state is stabilised with respect to the complex conjugate mode which travels in the positive z-direction. This behaviour can be inferred from the slope ∂Re III n /∂α| α=0 • / = 0 in figure 3(a). After the neutral mode III has become critical for α > 5.8 • the propagation direction of mode III turns in the positive z-direction at α = 9.6 • (figure 3c). For larger α, the magnitude of the oscillation frequency increases monotonically with α, indicating an increasing phase velocity. For α 15 • , the increase of |ω n | is approximately linear in α.
Mode III is illustrated in figure 5 for α = 0 • , 22.5 • and 35 • , showing the same quantities as in figure 4. Similar to mode II the critical mode III arises as vortices which are the strongest near the upstream wall at x = −1/2. The vortices, best seen in figure 5(b) for α = 22.5 • , have a similar extension in the wall-normal direction in both cases. However, different from mode II at α = 6 • , mode III has a much larger wavelength throughout the range of yaw angles over which it is critical (cf. figures 3 and 5). Near the downstream wall we do not find the same vortices. Rather, in the plane y = 0 larger-scale vortex structures occupying the full width of the cavity can be identified. Furthermore, the pitch of the vortices of mode III is larger than that for mode II which can be seen by correlating the vortex structures (e.g. the isolines of v in the plane y = 0) near the two walls at x = ±1/2.
Common to modes I, II and III, they all extract most of their energy from the basic state in the curved boundary layer of u 2-D 0 on the upstream wall (see also Albensoeder et al. 2001). In addition, we also find minor contributions to the energy production near the bottom of the cavity. Correspondingly, the energy budgets of all three modes are very similar (table 2, figure 3d). All modes destabilise the basic flow primarily through the process described by i 2 (3.3c). Therefore, the modes I, II and III may be called spiral Taylor-Görtler vortices.
For constant Reynolds number Re, the effective Reynolds numbers for the recirculating part of the basic flow and for the Couette part of the basic flow scale like Re 2-D ∼ cos(α) and Re C ∼ sin(α), respectively. Therefore, as α increases for Re = const., the velocity magnitude and shear rate of the recirculating part of the basic flow decrease monotonically, because Re 2-D decreases. Since the shear in the two-dimensional curved boundary layer drives the Taylor-Görtler instability, an increase of the yaw angle should result in a certain stabilization of the basic flow with respect to modes building on the Taylor-Görtler mechanism. However, a stabilisation of the basic flow with increasing α is not generally found. For instance, the critical Reynolds number for mode III decreases and reaches a  minimum near α ≈ 22.5 • . This behaviour can be explained by the Couette part of the basic flow w C 0 which exhibits a plateau in the centre of the cavity and significant gradients near the boundaries from which kinetic energy can be extracted. Furthermore, the vortex structures of mode III at α = 22.5 • are larger than at α = 0 • (figure 5a,b). Associated with the structural changes of the neutral mode, also the region of energy production within which i 2 (lift-up) is dominant changes and, for increasing α, extends over the bottom of the cavity up to the wall downstream of the moving lid (figure 5d-f ). Interestingly, the extended isosurfaces of perturbation-energy production of the oscillatory mode III for α = 0 • make an angle of ≈ 25 • with respect to the direction of motion of the lid (figure 5g). As the yaw angle α is increased, the orientation of the production isosurfaces, and, thus, the perturbation vortices, turns into the direction of the lid velocity such that, near the minimum of the critical Reynolds number at α ≈ 22.5 • , the perturbation-energy production surfaces are approximately aligned parallel to the (x, y) plane (figure 5h).

Mode IV
) mode III is less efficient in extracting energy from the basic flow, and for α > 37.4 • , mode IV (red in figure 3) becomes critical with a critical wave numbers k c ≈ 8, slightly larger than the one of mode III. Furthermore, the magnitude of the phase velocity |−ω c /k c | is about a factor of two smaller than the one of mode III. The lift-up mechanism I 2 becomes even more preponderant in the integral perturbation-energy budget, while the contribution due to I 3 (anti-lift-up) decreases to 14 % for α = 40 • (table 2, figure 3d).
Mode IV is visualised in figure 6 for α = 40 • . The critical mode in the (x, y) cross-section arises as a series of counter-rotating vortices aligned, one after the other, along the outer streamlines of the two-dimensional part u 2-D 0 of the basic flow (figure 6b). As the wave propagates, the perturbation vortices travel, in the (x, y) plane, in the direction of u 2-D 0 (clockwise in figure 6b). This property of the perturbation flow in the (x, y) plane is similar to the pure two-dimensional critical mode in a square cavity in which a double row of vortices (vortex street) circulates about the vortex core (Cazemier, Verstappen & Veldman 1998;Auteri, Parolini & Quartapelle 2002a).
In the third dimension the perturbation vortices extend in a spiral fashion wrapping about the basic vortex core of u 2-D 0 . This is illustrated in figure 6(d) by isosurfaces of λ 2 = −3 which are approximately aligned with the three-dimensional basic flow u 0 . The regions of high local energy production i (blue in figure 6b,c) are approximately centred between two neighbouring spiral perturbation vortices. Figure 6(c) reveals that the structures of i, i 2 and |u | superimpose almost perfectly. The close agreement of the isosurfaces of i and i 2 is not surprising, because I 2 represents 82 % of I (table 2). But we also observe that the threads of i are nearly congruent with the isosurfaces of |u | and, hence, with streaks of the perturbation flow. Close to the lid the production isosurfaces show a double structure ( figure 6b,c), an effect which is due to the strong gradients of the basic flow close to the lid. The spiral perturbation vortices wrapping about the swirling basic vortex core are conceptionally similar to the azimuthally periodic spiral vortices arising in spiral Couette and spiral Poiseuille flow (Ludwieg 1964;Meseguer & Marques 2000. The number of spiral perturbation vortices in the present cavity is difficult to quantify exactly, because their diameters in the (x, y) plane vary strongly. However, seven of the vortices are clearly visible in figure 6(b,c). After passing the downstream singular corner and being transported by the basic flow, the perturbation vortices grow stronger and attain their maximum strength when they reach the upstream wall ( figure 6a,b), where also the local production i reaches its maximum. During the acceleration of the basic flow along the moving lid the perturbation vortices are strongly damped, which is confirmed by the large dissipation of perturbation energy d * near the moving lid (not shown). An animation of the vortex motion in the (x, y) plane with a zoom into the downstream corner is provided as a supplementary material available at https://doi.org/10.1017/jfm.2021.804. Unlike modes I, II and III, the axes of the perturbation vortices of mode IV are strongly deflected from the (x, y) plane. In this respect, the critical mode IV is a member of another family of modes whose vorticity is primarily aligned in the z-direction.  Consistent with the trend observed for modes IV and V, the critical wavelength λ c increases with α. Finally, mode VII undergoes a dramatic stabilisation for α 80 • , rapidly reaching critical Reynolds numbers which are beyond those for which our numerical solver has been designed. We anticipate that other modes similar to mode VII become critical as the yaw angle is further increased.

Modes V to VII
Common to the spiral vortices of modes IV to VII they grow in size and strength as they travel in the (x, y) plane with the basic flow u 2-D 0 downstream from the singular corner at (x, y) = (1/2, 1/2). The perturbation vortices reach their maximum size and strength near the wall upstream of the moving lid. As the perturbation vortices travel along the moving lid, they shrink and cannot be unambiguously identified anymore when they pass the downstream singular corner. The damping and size reduction of the perturbation vortices seem to be strongly related to the acceleration of the basic flow along the moving lid during which the length scale of the shear layer of u 2-D 0 with negative vorticity shrinks, whereas the perturbation vortices grow in the widening boundary layer of the decelerating flow u 2-D 0 along the downstream wall which has positive vorticity. Therefore, the helical nature of the perturbation vortices is disrupted at the downstream singular corner and the total number of vortices present at any time in a cross-section at constant z cannot be precisely specified, even though from an analogy with spiral Couette flow an even number of vortices is expected. 4.2.5. General properties of the critical modes for Γ = 1 All critical modes found for Γ = 1 are destabilised by the lift-up mechanism represented by i 2 in (3.3c), whose integral contribution ranges from 68 % for α = 0 • to 94 % for α = 75 • , as can be seen in table 2 or figure 3. The energy production is most pronounced near the upstream wall of the cavity. Therefore, we conclude that the mechanism is essentially a modification of the Taylor-Görtler instability mechanism which is well established for α = 0 • . This interpretation is supported by the expectation that the pure Couette part u C 0 of the basic flow (4.1) is linearly stable and no linear instability mechanism can be derived from this parallel shear flow alone. Since the Taylor-Görtler-like vortices are aligned with the direction of the total basic flow u 0 , the pitch of the vortices of the critical helical modes increases with α. For small α, the diameter of the vortices is small, as they scale with the boundary layer thickness of u 2-D 0 . Therefore, as α increases, more and more helical vortices penetrate the unit cell defined by one wavelength λ c of the perturbation flow (see e.g. figure 6c,d). As the basic flow turns predominantly into the span (z) direction and the critical wavelength increases, the Taylor-Görtler-like vortices become longer and can grow to larger diameters (see figures 6b to 9b), because the characteristic length scale becomes the depth of the cavity. Therefore, the trend is reversed and a lesser number of vortices penetrate the unit cell.
The phase velocity c n = −ω n /k n of each neutral mode is shown in figure 10(a). For each mode, it scales almost linearly with α. Scaling the phase velocity with the spanwise component of the lid velocity Re sin α in figure 10(b), the scaled propagation speed is almost independent of α for the second family of modes IV, V, VI, VII as well as for mode I, which is stationary at α = 0 • . For all these modes, the scaled phase velocity is always less than 0.5.
On the other hand, the scaled propagation speeds (figure 10b) of neutral modes which are travelling at α = 0 • (modes II and III) necessarily diverge as α → 0. However, the rescaled phase velocities nearly saturate for α > 20 • . Moreover, as discussed earlier, the direction of propagation of modes II and III is initially opposing the direction of the spanwise lid motion for small yaw angles α.  Theofilis et al. (2004) who also performed a linear stability analysis of the same flow for α = 22.5 • , 45 • and 67.5 • , but did not find any unstable eigenmode at Re = 800 in the range k ∈ [0, 25]. Since our results are at variance with these previously published data, we carried out an independent nonlinear numerical simulation using NEK5000 for Γ = 1, α = 22.5 • and for a slightly supercritical Reynolds number with 800 > Re = 650 > Re c = 619.9 and periodic boundary conditions in z with period λ = 2π/k c . Impulsively starting the lid motion from a state of rest at t = 0, we find the basic flow to be established near t ≈ 0.5. At about the same time small amplitude oscillations of w become visible and start growing exponentially in an oscillatory fashion (figure 11a). Fitting the signal w(t) within the grey region shown in figure 11(a) by w F (t) = a + be σ F t sin(ω F t + c), we find the growth rate σ F = 5.66 > 0 and the angular frequency |ω F | = 478.5 for spectral element polynomial order p = 6. The fit is not shown, because it cannot be visually distinguished from the simulation data on the scale of figure 11(a). The growth rate σ F compares reasonably well with the real part σ = 4.63 of the eigenvalue of the linear stability problem at the same Reynolds and wave number, and the oscillation frequency |ω F | is in excellent agreement with the imaginary part of the eigenvalue |ω| = 475.7. While relative deviations among growth rates with σ ≈ 0 obtained by different methods may seem large, the deviations between corresponding values of Re c = Re(σ = 0) are not. This is demonstrated in figure 11(b) which shows growth rates obtained by different methods and resolutions. From the growth rates shown we obtain the critical Reynolds numbers Re c = 613.9, 616.9 and 619.9, respectively, by quadratic interpolation of the data obtained with NEK5000 and p = 6 (green), and those of the finite element linear stability analysis with N = 80 (orange) and N = 40 (blue). Thus, the critical Reynolds number varies by less than 1 % upon changing the methods and resolutions of the discretisation. Moreover, the growth rates at Re = 650 obtained using the higher polynomial orders p = 8 and 10 deviate from the growth rate obtained for p = 6 only at the fourth significant digit. Therefore, the three-dimensional simulation can be considered to be fully converged. Our independent nonlinear simulation of the cavity flow thus confirms the instability of the basic flow with a critical Reynolds number Re c < 620, consistent with the present linear stability analysis.

4.3.
Linear stability for Γ = 0.5 An overview on the linear stability analysis for a shallow cavity with Γ = 0.5 is shown in figure 12. For the classical case with α = 0 • , we find that Re c = 706.7, k c = 10.64 and ω c = 818.9. This is in very good agreement (differences less than 1 %) with the result of 819.9 ± 4. Except for k c our result is also in good agreement with the data of Theofilis et al. (2004) (mode T2 from their table 7: Re c = 720.18, k c = 11.40, ω c = 838).

Mode I
The critical mode I at α = 0 • is oscillatory and arises as a pair of waves with relatively short wavelengths which propagate in the positive or negative z-direction. As α increases, the degeneracy of the two waves (and the associated critical parameters) is removed. While the wave propagating in the positive z-direction is stabilised, the wave propagating in the negative z-direction, opposite to the z-component of the lid velocity vector, is destabilised and becomes the critical mode. Increasing α, the phase velocity of the critical mode slows down and the mode becomes stationary at α = 10.1 • . For larger yaw angles, the mode starts propagating again, but now in the direction of the spanwise lid motion. The same qualitative dependence on α of the propagation direction was found before for mode III Γ =1 .
For shallow cavities and elevated Reynolds numbers, the basic flow at α = 0 • arises as a spanwise vortex near the downstream end of the cavity. Similarly as described by Albensoeder et al. (2001) for Γ = 0.25, the present basic flow at α = 0 • and Γ = 0.5 becomes unstable due to a centrifugal instability in the region where the basic vortex flow separates from the bottom wall. As α increases, the Couette part of the basic flow Mode I is the critical mode over a wide range of α ∈ [0 • , 45.7 • ] with a minimum critical Reynolds number of Re c = 599.5 at α = 13.6 • . For α = 0 • , the critical mode I is primarily destabilised by the lift-up process described by i 2 . However, as α increases to α = 15 • , the integral contribution I 2 is reduced and, thus, cannot explain the destabilisation by about Re ≈ 100 compared with α = 0 • . The reason is I 1 gains importance and overcompensates the reduction of I 2 (see table 2 or figure 3d). This trend continues and, at α = 40 • , I 1 and I 2 have become comparable in magnitude with a share of 42 % and 54 %, respectively, of the total energy budget. Different from i 2 , i 1 vanishes in parallel flow, because the energy-transfer process of i 1 requires the direction of the basic flow to change perpendicular to itself. Therefore, i 1 cannot build on gradients of the Couette part of the flow u C 0 . Moreover, as the swirling part of the basic flow u 2-D 0 becomes weaker as α increases, the destabilisation from α = 0 • to α = 15 • cannot be explained by u 2-D 0 alone. Therefore, it is the change of the modal structure accompanied with the increase of α which must be responsible for the ability of the critical mode to extract more energy from u 2-D 0 via i 1 , despite u 2-D 0 becoming weaker with α.
The change of the critical mode is demonstrated by figure 14 for α = 40 • . Compared with α = 15 • (figure 13) the vortex structures of the critical mode for α = 40 • have become stronger near the downstream half of the cavity and weaker in the upstream half. The critical mode now arises mainly as vortices nearly perpendicular (but slightly tilted) to the moving wall and near the downstream end of the cavity ( figure 14a,b). This indicates that the critical mode I for α = 40 • mainly receives its kinetic energy from gradients of the basic flow in the downstream half of the cavity, where the two regions in which i 1 and i 2 dominate are interwoven in a complicated fashion (figure 14c).

Modes II and III
As the Couette part of the basic flow becomes dominant upon an increase of α, the critical mode changes to mode II at α = 45.7 • (green in figure 12). Mode II has a similar critical wave number as mode I and it is illustrated in figure 15 for α = 50 • . The structure of mode II is quite different from that of mode I and resembles the spiral critical modes for Γ = 1 discussed in § 4.2. From table 2, the lift-up mechanism I 2 dominates the energy budget of the critical mode II, similar as for Γ = 1. From figure 15(b) we can see different patches (blue) of localised energy production which are arranged around the periphery of the basic vortex u 2-D 0 . These local production regions extend in a spiral fashion in three dimensions as shown in figure 15(c). The threads of energy production feed energy to the helical-type of perturbation vortices which are visualised in figure 15(d) by isosurfaces of λ 2 . One can identify six spiral vortices in the bulk of the flow, not all of which are visible in an arbitrary cross-section, because, similar as for modes IV Γ =1 to VII Γ =1 for Γ = 1, the spiral perturbation vortices are strongly suppressed when passing the downstream corner. From figure 15(d) we can also recognise two weaker vortices per period of the flow reaching into the upper half of the cavity near the upstream wall (x < −0.2, y > 0.05). Since the basic flow is weak in this region, these perturbation vortex 'appendices' have little effect on the instability and contribute less than ≈5 % to the energy budget.
As the inclination angle is further increased, mode III (pink in figure 12) becomes critical at α = 55.9 • . As an example for mode III, we consider α = 65 • in figure 16. For mode III, the energy production I 2 by the lift-up process is even more important than for Γ = 1 with I 1 , I 3 and I 4 altogether contributing less than 16 % to the total energy transfer (table 2). While mode III is similar to mode II, the wavelength of mode III is about twice as long as that for mode II, and the perturbation vortices are more aligned with the z-direction. The critical mode III is characterised by helical vortices which wind about the recirculating basic vortex u 2-D 0 . Three of the vortices are clearly visible in the cross-section shown in figure 16(b). At the phase depicted, the major perturbation vortex extends somewhat into the more quiescent region (with respect to u 2-D 0 ) near the upstream 928 A25-26 corner of the cavity. In the plane shown the vortices are fed by four patches (blue) of energy production (a double patch near the moving lid). Similar as for modes VI Γ =1 to VII Γ =1 for Γ = 1, the spiral vortices are suppressed near the downstream singular corner and the pitch of the perturbation vortices near the moving lid differs from the one near the bottom of the cavity. Increasing α from α = 55.9 • the critical curve reaches a local minimum at α = 63.4 • . On a further increase of α the basic flow is rapidly stabilised and we did not follow the critical curve beyond α = 75 • . Up to this inclination angle, mode III remains critical and retains the same characteristics as for α = 65 • . 4.4. Linear stability for Γ = 2 An overview on the neutral Reynolds numbers for a deep cavity with Γ = 2 is shown in figure 17. Again, the critical data for the stationary mode at α = 0 • , Re c (α = 0 • ) = 444.90 with a high wave number denoted III Γ =1 is also shown. Its relevance and character will be discussed in § 4.5.2 below. As α is increased from zero, the stationary mode starts drifting in the positive z-direction. The character of the critical mode does not change very much even at α = 50 • , for which the critical mode is shown in figure 18. It is very similar to the stationary mode for α = 0 • reported in figures 20 and 21 of Albensoeder et al. (2001). The most important region of energy production (again I 2 is dominant) is located in the curved boundary layer of u 2-D 0 just before the basic vortex flow separates from the downstream wall at x = 1/2 (figure 18b). In the (x, y) plane the perturbation flow is a vortex slightly offset from the basic state vortex in the direction towards the cavity centre. The perturbation vortex changes its sense of rotation periodically in z, which leads to a modulation of the total finite-amplitude vortex flow as has been observed experimentally for α = 0 • and Γ = 1.6 by Siegmann-Hegerfeld, Albensoeder & Kuhlmann (2013). Associated with the perturbation flow are periodic up-and down-flow (in the y-direction) regions at the midplane y = 0 shown in figure 18(a,b) which arise just at the edge of the basic state vortex.
As the inclination angle is increased beyond α ≈ 45 • , the wavelength of the critical mode I increases and the critical mode changes to mode II at α = 78.9 • . Due to the modal change the wavelength suffers a step reduction, but it grows again with α and reaches λ c = 2π/k c = 9.2 at α = 85 • . Mode II is shown in figure 19 for α = 85 • . Due to the long wavelength the structure of the perturbation flow is stretched in the z-direction. This is a consequence of the wall-bounded Couette part u C 0 of the basic flow. Yet, the region near the separation line of the basic flow from the wall at x = 1/2 remains of crucial importance for the transfer of kinetic energy to the perturbation (figure 19b), now being nearly exclusively due to I 2 . For α = 85 • , the critical mode has a significant spanwise velocity component w. The ratios of the magnitude of the maxima of the perturbation velocity components u and v compared with max(w) for α = 85 • are max(u)/max(w) = 0.1950 and max(v)/max(w) = 0.1456, whereas the corresponding ratios for α = 0 • are max(u)/max(w) = 2.2195 and max(v)/max(w) = 1.8370. This indicates the predominance of streaks in the perturbation flow. The streaks arise as elongated structures of u illustrated by isosurfaces of |u | = 0.5 max |u | in figure 19(c) (the isosurfaces of i look very similar). The isosurfaces are coloured according to the z-component of the perturbation velocity parallel to the basic flow w = e z · u . In addition, u ⊥ is shown by arrows. The velocity components of u ⊥ , representing the streamwise perturbation vortices, are almost confined to the (x, y) plane. Furthermore, the z-component e z · u ⊥ is about ten times smaller in magnitude than the maximum streak velocity max |u |. The framing of the streaks by pairs of streamwise vortices is a visual confirmation for the lift-up effect being the main instability mechanism.
Similar as for Γ = 1 and Γ = 0.5, the basic flow is strongly stabilised with respect to linear perturbations as α → 90 • (figure 17). Finally, common to all aspect ratios is an increase of the critical wavelengths with α as well as an increase of the size of the perturbation flow structures in the cross-sectional (x, y) plane.

4.5.
Common properties of the critical modes and dependence on the aspect ratio 4.5.1. Comparison of results for Γ = 0.5 and Γ = 1 The instability scenario upon a variation of α for Γ = 0.5 is similar to that for Γ = 1. Mode I Γ =0.5 of the former case corresponds to mode III Γ =1 of the latter: the propagation direction for small yaw angles α is opposing the spanwise direction of the lid motion, but progressively aligns with it as α increases. In both cases the critical Reynolds number decreases with α, before increasing to values Re c (α) > Re c (α = 0 • ). Furthermore, modes II Γ =0.5 and III Γ =0.5 seem to correspond to the modes V Γ =1 and VI Γ =1 , respectively. This is also suggested by the wave numbers of modes III Γ =1 and V Γ =1 being very close, a behaviour which is also found for I Γ =0.5 and II Γ =0.5 .
This interpretation is also corroborated by Re n (Γ ) and σ (k) for α = 0 • provided by Albensoeder et al. (2001) in their figures 6 and 16, suggesting that the mode of Ding & Kawahara (1999) (mode III Γ =1 ) is identical, i.e. smoothly connected in Γ , with mode I Γ =0.5 . Besides, modes (VI, VII) Γ =1 and III Γ =0.5 are characterised by common structures with perturbation vortices winding about the core of the basic vortex. The perturbation vortices of these modes are damped so strongly near the moving lid such that they are practically disconnected from the vortices generated downstream of the singular corner.

Critical curves as functions of the aspect ratio
In order to verify whether the modes observed for different aspect ratios are connected mutually, neutral curves are computed, now varying the aspect ratio Γ , for two representative yaw angles α = 25 • and α = 50 • . The mesh used for the computation is initially 80 × 80 and refined in the same way as described in § 3. The computational domain is rescaled upon variation of the aspect ratio, and the number of unknowns remains the same for all Γ . Inferring from table 2, the critical parameters should vary from the previously provided results only from the third or fourth significant digit depending on the aspect ratio. Figure 20 shows the variations of the neutral Reynolds numbers, wave numbers and frequencies of the neutral modes for α = 25 • . It reveals that the critical modes at Γ = 0.5 and Γ = 1 are indeed the same, continuously transforming into each other as the aspect ratio is varied. However, the critical mode observed at Γ = 2 is only critical for Γ > 1.78 and stabilises rapidly as the aspect ratio is decreased. In the range Γ ∈ [1.11, 1.78] mode III Γ =1 (cyan in figure 20) is critical. This mode is only slightly more stable than mode I Γ =2 at Γ = 2 (see also figure 17). The structure of the velocity field and the local energy production rate of mode III Γ =1 are extremely similar to those of mode III Γ =1 . Indeed, the neutral modes III Γ =1 and III Γ =1 are the same. Only in a certain narrow range of aspect ratios near Γ ≈ 1.1 the neutral curve Re n (k) develops two local minima within a small distance in k near k ≈ 5.7. At Γ ≈ 1.11, at which both minima are the same, the critical Reynolds number, belonging to the absolute minimum of Re n (k), switches from one minimum to the other (from k c = 6.38 to k c = 5.06). It is interesting to note that, for α = 0 • , a similar jump of the neutral mode III Γ =1 (alias I Γ =0.5 , alias the mode of Ding & Kawahara 1999) arises at Γ = 0.94 and Re n = 933 with k n switching from k n = 6.1 to k n = 7.6 (all data extracted graphically from figure 6 of Albensoeder et al. 2001). Unlike for α = 0 • where mode III Γ =1 is only a neutral mode, it is the critical mode for α = 25 • over a wide range of Γ . The destabilisation trend with increasing α of this mode can be recognised for all aspect ratios investigated (figures 3, 12 and 17). Thus, the sequence of critical modes upon a variation of Γ for α = 25 • is very similar to the one observed in the classical cavity for α = 0 • (Albensoeder et al. 2001), except for the absence of the high-wave-number Taylor-Görtler mode (mode I Γ =1 ) when α = 25 • .
The dependence of Re c on Γ for α = 50 • is provided in figure 21. As anticipated, the critical modes at Γ = 0.5 and Γ = 1 turn out to be the same, while the critical mode at Γ = 2 is a different branch, not connected to the critical modes for Γ = 1 and Γ = 0.5. Interestingly, the neutral Reynolds numbers, wave numbers and frequencies of the modes IV Γ =1 and V Γ =1 do not change much for Γ 1. Unlike for α = 25 • , mode I Γ =2 is critical over a wider range of aspect ratios, i.e. for Γ > 1.04.

Discussion and conclusion
The linear stability of the steady flow in a rectangular cavity driven by the oblique motion of a lid has been investigated with respect to spatially periodic perturbations. The parameter space for this problem is made of the Reynolds number Re, the inclination angle of the lid α and the cross-sectional aspect ratio Γ . Three representative cavities have been investigated: a cavity with a square cross-section (Γ = 1), a shallow cavity (Γ = 0.5) and a deep cavity (Γ = 2).
The basic flow in the obliquely driven cavity is a swirling flow made of a superposition of two types of motion. One is the well-known two-dimensional cavity flow u 2-D 0 (x, y), driven by the x-component Re cos(α) of the normalised lid velocity which is reduced compared with the absolute normalised lid velocity Re. The other part of the flow field is made by the parallel Couette-type of flow u C 0 (x, y) in the spanwise direction. It is driven by the spanwise component Re sin(α) of the normalised lid velocity. While the recirculating part of the motion is independent of the spanwise motion, the Couette part u C 0 (x, y) = w 0 (x, y)e z of the basic flow is affected by u 2-D 0 which advects the spanwise momentum w 0 .
Critical Reynolds numbers as a function of the yaw angle α have been computed for all three aspect ratios. For α = 0 • , the accurate stability boundaries provided by Albensoeder et al. (2001) are recovered. The slope ∂Re c /∂α| α=0 • = 0 of the critical curve at α = 0 • vanishes for critical modes which are stationary (Γ = 1, Γ = 2), because the isolated real eigenvalue must evolve continuously and symmetrically with respect to α = 0 • . Therefore, the critical Reynolds number increases quadratically as α is increased from zero, and the critical modes start drifting in the direction of the spanwise lid motion (here in the positive z-direction). On the other hand, the degeneracy of the critical Reynolds number for oscillatory eigenmodes at α = 0 • is removed for α / = 0 • and ∂Re c /∂α| α=0 • = ±a / = 0, where a = const., such that the critical Reynolds number is always reduced and the critical mode for α > 0 evolves from one of the travelling waves at α = 0 • . For those latter cases, we find that the critical/neutral mode which destabilises the basic state for small α (modes III Γ =1 and I Γ =0.5 ) is propagating in the spanwise direction opposite (ω > 0) to the spanwise component of the lid motion. As α increases, the critical mode becomes stationary near α ≈ 10 • and turns propagating parallel to the z-component of the lid motion for larger α.
When α is small, the basic flow is dominated by the recirculating part of the flow u 2-D 0 . In this situation all critical modes arise in the curved boundary layer of u 2-D 0 and receive their kinetic energy mainly due to the lift-up process described by i 2 = −D −1 * u · (u ⊥ · ∇u 0 ). The similarity of the modal structures and of the instability mechanism with those of the classical lid-driven cavity at α = 0 • (Koseff & Street 1984a;Albensoeder et al. 2001) suggests calling the modes for α / = 0 • spiral Taylor-Görtler vortices. The conceptual relationship between the Taylor-Görtler vortices in lid-driven cavities and Taylor vortices between concentric rotating cylinders calls for a qualitative comparison between the present results and the flow between rotating cylinders with axial through flow. For intermediate yaw angles, the critical modes in the oblique cavity are spiral waves. Between concentric cylinders azimuthally periodic spiral Taylor-Görtler vortices are known to arise when (a) the inner cylinder translates axially with zero axial pressure gradient (spiral Couette flow) (Ludwieg 1964;Wedemeyer 1967;Meseguer & Marques 2000), (b) the axial flow is purely pressure-driven (spiral Poiseuille flow) (Meseguer & Marques 2002), or (c) the flow is driven both by an axial cylinder motion and a pressure gradient (spiral Couette-Poiseuille flow) (Ali & Weidman 1993;Meseguer & Marques 2000). The main distinction from the obliquely driven cavity is the strong perturbation of the azimuthal symmetry by the rectangular geometry and the discontinuous boundary conditions. Nevertheless, for sufficiently strong spanwise/axial flow, the basic flow is destabilised by spiral vortices whose number in the (x, y) plane typically increases for intermediate values of α (Ali & Weidman 1993;Meseguer & Marques 2002). While the spiral vortices are strictly periodic between concentric cylinders, they grow downstream from the downstream singular corner in the obliquely driven cavity flow. This growth is in parallel with the growth of the thickness of the curved boundary layer downstream from the singular edge. Another similarity, for small yaw angles, is the spanwise wave propagation opposite to the spanwise direction of the lid motion of certain modes without a pronounced spiral character. In the spiral Couette-Poiseuille flow Ali & Weidman (1993) found a similar retrograde drift with respect to the axial motion of the inner cylinder of the toroidal modes for dominant swirl if the radius ratio is sufficiently large.
For large inclination angles α → 90 • , the recirculating part u 2-D 0 (x, y) of the basic flow diminishes and the basic flow tends to the confined Couette flow (4.1). As the basic flow becomes more parallel, the most dangerous modes become elongated in the spanwise, i.e. streamwise, direction. In the limit the energy production terms i n for n = 1, 3 and 4 vanish and only the lift-up term i 2 remains. This trend is also reflected by the integral contributions to the kinetic energy budget of the perturbation flow listed in table 2. As long as even a weak recirculating part of the basic flow can provide a feedback from the streaks to the nearly streamwise vortices which create the streaks, a linear instability is possible. With the recirculating basic flow getting weaker the feedback becomes weaker and the stability boundary Re c (α) increases strongly as α → 90 • . This interpretation is consistent with the previous investigation of Theofilis et al. (2004). There does not seem to be any linear process which could destabilise the wall-bounded Couette flow at α = 90 • .
Although the general interpretation of Theofilis et al. (2004) regarding the flow stabilisation at very large yaw angles is corroborated by the present work, we found that the critical Reynolds numbers for intermediate yaw angles are much lower than the estimates provided by them. For Γ = 1 and all three yaw angles α = 22.5 • , 45 • and 67.5 • , they bracketed the critical Reynolds number to be in the range Re c ∈ [800, 900]. On the contrary, for Γ = 1, we find that Re c (α = 22.5 • ) = 619.9 (confirmed by independent nonlinear simulation), Re c (α = 45 • ) = 753.3 and Re c (α = 67.5 • ) = 643.3. Therefore, the lower bound Re = 800 on Re c specified by Theofilis et al. (2004) seems too large by 29 %, 6 % and 24 %, respectively.
Periodic instabilities of the flow in a cavity infinitely extended in the spanwise (z) direction have been analysed. A natural extension of this problem would be an investigation and the yaw angle α, similar as for the sensitivity with respect to changes in the boundary conditions (Meliga, Sipp & Chomaz 2010).
To find (Re c , k c ) as a function of α using the above iteration, a good initial guess is required. This is obtained by extrapolating the converged critical data obtained for the previous values of α. Based on the N + 1 known data (α i , Re c,i , k c,i ) for i ∈ 0, . . . , N the distance function from the first point (α 0 , Re c,0 , k c,0 ), is evaluated for i = 1, . . . , N, where the Nth data set is the data set found in the last converged iteration. The coefficient a = 0.1 has been selected in order to improve the condition number of the fit, because the Reynolds number is typically two to three orders of magnitude larger than α and k, which also applies to their variations. Based on the parametrisation d i each of the three quantities (α, Re c , k c ) is fitted by a polynomial P f (d i ) of maximum order three, where f ∈ [α, Re c , k c ]. The coefficients are obtained by least-squares minimisation of the functional N i=1 w i [ f i − P f (d i )] 2 , where the weights w i are selected to give preference to the last point by setting w i = i. The three polynomials obtained are then evaluated for d > d N to arrive at the new initial guess for (α N+1 , Re c,N+1 , k c,N+1 ). Using a low polynomial order (P ∈ P 3 ) renders the method stable.
Appendix B. Sensitivity of the eigenvalues with respect to a variation of the wave number Based on the eigenvalue equation (2.6) of the form [γ i M + A(k)]q i = 0, whereq i = (û i ,p i ) T is the vector of the mode variables, we use the classical approach of Marquet, Sipp & Jacquin (2008) who considered the sensitivity of the eigenvalue γ i and eigenvectorq i with respect to a variation of the basic flow. As in Marquet et al. (2008), we use an optimal control framework and defineq i and γ i as state variables and k as a control parameter. The Lagrangian is defined as whereq † i = (û † i ,p † i ) T defines the adjoint variables. Formally the first term on the right-hand side of (B1) is the cost function of the problem, while the second term is enforcing the constraints through Lagrangian multipliers. Cancelling the derivative with respect to the Lagrangian multiplierq † i is equivalent to enforcing the state equation, i.e. the eigenvalue problem. Cancelling the derivative with respect to the state variablesq i and γ i is equivalent to solving the adjoint problem and enforcing a normalisation condition on q i andq † i . Evaluating the derivative of the Lagrangian with respect to the parameter k will give us the gradient of the cost function and, by construction, the sensitivity of the eigenvalue with respect to variations of k.
B.1. Differentiating L with respect to γ i Evaluating the differential of the Lagrangian functional (B1) with respect to δγ i and requiring δ γ i L = 0 yields This leads to the normalisation condition B.2. Differentiating L with respect to any eigenvectorq i Setting the differential of the Lagrangian in the direction δq i to zero, δq i L = 0, we obtain which leads to the adjoint eigenvalue problem B.3. Differentiating L with respect to k Considering δ k L = 0 we obtain In the present formulation A is the right-hand side of (2.6). Taking the derivative of the vectorial form of these equations with respect to k we obtain The first term on the right-hand side derives from the viscous diffusion. The following terms derive from the transport of perturbation momentum in the spanwise direction, the pressure gradient and the continuity equation. Finally, we obtain the sensitivity of the eigenvalue γ i with respect to wave number changes where we used the normalisation û i ,û † i = 1 which derives from (B3).