1. Introduction
Flow separation and recirculation are of great interest as they play an important role in the phenomena involved in transport and mixing processes. The flow past open cavities is a prototype of geometrical configurations characterized by a finite region of separated flow. The identification of the flow characteristics related to instability mechanisms (e.g. coherent structures) is also of practical importance since these may lead to resonances, acoustic noise or structural vibrations. Rockwell & Naudascher (Reference Rockwell and Naudascher1978) classified the unstable behaviour of this kind of flow into fluid-dynamic, fluid-resonant, and fluid-elastic.
Acoustic resonance (i.e. fluid-resonant behaviour) has received a remarkable amount of attention in the past due its relation with the noise generation process (Rossiter Reference Rossiter1964; Yamouni, Sipp & Jacquin Reference Yamouni, Sipp and Jacquin2013). In this case, there exists a feedback mechanism between the unstable shear layer susceptible to Kelvin–Helmholtz instability and the pressure waves (Rowley, Colonius & Basu Reference Rowley, Colonius and Basu2002). Gharib & Roshko (Reference Gharib and Roshko1987) observed experimentally that the increase of the cavity length to depth ratio
$(L/D)$
led to a different kind of instability, the so-called ‘wake mode’. This global instability relies on a purely hydrodynamic mechanism (the oscillation Strouhal number is weakly dependent on the Mach number) and is characterized by a large-scale vortex shedding.
1.1. Experimental and numerical investigations of three-dimensional instabilities
Three-dimensional structures in cavities with large span-to-chord ratios were first reported by Maull & East (Reference Maull and East1963) using oil flow and static-pressure measurements. These authors observed stable cell formations inside open rectangular cavities of different aspect ratios
$(L/D)$
. Rockwell & Knisely (Reference Rockwell and Knisely1980) observed a strong coupling between the growth of the primary vortices and the increase of the vorticity of the longitudinal vortices past the downstream edge of the cavity. They concluded that, because of this interaction, the cavity core vortex presents an ordered spanwise modulation. These coherent structures were analysed in detail in the laminar regime by Faure et al. (Reference Faure, Adrianos, Lusseyran and Pastur2007, Reference Faure, Pastur, Lusseyran, Fraigneau and Bisch2009) in cavities with an aspect ratio varying between 0.5 and 2. Using smoke visualizations, Faure et al. (Reference Faure, Adrianos, Lusseyran and Pastur2007) report mushroom-like counter-rotating vortical structures and suggest that the resulting inherent flow is the consequence of a centrifugal instability related to the cavity core vortex.
Brés & Colonius (Reference Brés and Colonius2007b ) performed numerical simulations of the linearized compressible Navier–Stokes equations to investigate the instability features in open cavity flows. These authors identify the onset of the first bifurcation over a wide range of Mach numbers and cavity aspect ratios and discuss the variations of the spanwise modulation with the cavity depth. Brés & Colonius (Reference Brés and Colonius2007a ) showed, furthermore, that the full nonlinear results, computed by direct numerical simulation (DNS), agree very well with the features of the three-dimensional global mode provided by linear stability analysis. Brés & Colonius (Reference Brés and Colonius2008) accurately discuss the properties, the structure and the nature of such instability. Using the generalized Rayleigh criterion (Bayly Reference Bayly1988), they show that the instability is a centrifugal instability associated with the closed streamlines inside the cavity.
Faure et al. (Reference Faure, Pastur, Lusseyran, Fraigneau and Bisch2009) performed an experimental investigation aimed at understanding the three-dimensional flow topology inside cavities of different shapes and for several Reynolds numbers. These authors also performed a secondary instabilities analysis and identified the relevant shear-layer and inner-cavity flow scales. Zhang & Naguib (Reference Zhang and Naguib2006, Reference Zhang and Naguib2008, Reference Zhang and Naguib2011) carried out a systematic study of the effect of the sidewalls on the unsteady open cavity flow at low Mach number. The presence of the sidewalls is shown to lead to strong amplifications of the pressure fluctuations inside the cavity. Lasagna et al. (Reference Lasagna, Donelli, Gregorio and Iuso2011) investigated the effects of a trapped vortex cell (TVC) on the aerodynamic performance of a wing model and found a three-dimensional organization of the flow inside the cell according to the value of the angle of attack and the Reynolds number. di Cicca et al. (Reference di Cicca, Martinez, Haigermoser and Onorato2013) reported time-resolved tomographic particle image velocimetry (PIV) measurements in rectangular cavities having the length-to-depth ratio equal to 2, 3 and 4.
Very recently, de Vicente et al. (Reference de Vicente, Basley, Meseguer-Garrido, Soria and Theofilis2014) examined both experimentally and numerically the instabilities over a rectangular open cavity of aspect ratio
$L/D=2$
. These authors compared the linear three-dimensional instability results with the spatial structure of the experimental fields, showing qualitative agreement for the main flow characteristics. Furthermore, they also show that modifications of the spanwise boundary conditions can cause significant alterations of the flow field due to nonlinear effects. Finally, Meseguer-Garrido et al. (Reference Meseguer-Garrido, de Vicente, Valero and Theofilis2014) present a systematic study of the onset of the first instability, varying the Reynolds number, the incoming boundary-layer thickness and the cavity aspect ratio.
Experimental studies on the formation of three-dimensional structures in the start-up flow inside a lid-driven cavity (LDC) can be found in Migeon, Texier & Pineau (Reference Migeon, Texier and Pineau2000), Guermond et al. (Reference Guermond, Migeon, Pineau and Quartapelle2002) and Migeon, Pineau & Texier (Reference Migeon, Pineau and Texier2003). These three-dimensional modulations inside the cavities are also documented by de Vicente et al. (Reference de Vicente, Rodriguez, Theofilis and Valero2010) and Gonzalez et al. (Reference Gonzalez, Ahmed, Kühnen, Kuhlmann and Theofilis2011) who considered complex cross-sectional shapes.
1.2. Short-wave asymptotic analysis
Instabilities in open and closed cavities, as well as in separated flows, are interpreted as centrifugal instabilities. Centrifugal short-wave instabilities were first considered by Bayly (Reference Bayly1988) who used the geometrical optics approximation and Flouquet theory to extend the classical Rayleigh theory for centrifugal instabilities to general inviscid planar flows. Bayly proposed diagonalizing the convective operator of the linearized Euler equations (LEes) and constructing linear asymptotic eigenmodes in the limit of large spanwise wavenumber using a Wentzel–Kramers–Brillouin–Jeffreys (WKBJ) expansion, localized on the closed orbit characterized by the maximum Floquet exponent. Later, the same author showed qualitative agreement between the results obtained with the linearized Navier–Stokes equations and the asymptotic predictions (Bayly Reference Bayly and Caflisch1989). Lifschitz & Hameiri (Reference Lifschitz and Hameiri1991) investigated the asymptotic instability features, considering the initial value problem for the LEes and for the linearized equations of gas dynamics. Their more general approach was able to include both exponential and algebraic growth in time. Many efforts have since been made to quantitatively link the short-wave asymptotics and the normal-mode analysis: Sipp, Lauga & Jacquin (Reference Sipp, Lauga and Jacquin1999) showed agreement between the optimal streamline (i.e. the streamline where the inviscid growth rate is maximum) and the spatial distribution of the unstable eigenmodes, and between the inviscid and viscous amplification rate of an elliptic instability. Gallaire, Marquillie & Ehrenstein (Reference Gallaire, Marquillie and Ehrenstein2007) examined the centrifugal instability of the separated region behind a bump and were able to make a composite estimation of the growth rate taking into account the viscous effects (see also Landman & Saffman Reference Landman and Saffman1987) and the short-wave inviscid asymptotic limit. Recently, Giannetti (Reference Giannetti2015) applied a WKBJ approach to investigate the nature of the secondary instability arising in the periodic wake of a cylinder.
In this context, the main goal of the present work is to characterize the instabilities of the flow past an open cavity, develop an asymptotic approach to understand the instability mechanisms and finally relate the results of the global and local asymptotic analysis. The specific aims of the work, defining the outline of the article, are to:
-
(i) provide an accurate estimation of the critical Reynolds number and the spanwise wavenumber of the first three-dimensional bifurcation in incompressible open cavity flows;
-
(ii) determine the instability core by means of the adjoint-based structural sensitivity analysis;
-
(iii) investigate the sensitivity of the leading instability to base flow modifications induced by a perturbation of the inflow profile or wall blowing/suction;
-
(iv) provide a quantitative prediction of the onset of the instability by means of the short-wave asymptotic theory;
-
(v) show that the inviscid structural sensitivity (i.e. structural sensitivity based on direct and adjoint eigenmodes computed at a Reynolds number higher than that of the base flow) is able to accurately predict the particle orbit that provides the main contribution to the instability;
-
(vi) suggest a generalization of the expression used to calculate the instability growth rate from the Floquet exponent to predict the frequency of the time-dependent modes emerging at supercritical conditions.
2. Theoretical framework
2.1. Geometrical configuration and base flow
We investigate the stability and sensitivity of the flow over a spanwise-uniform square open cavity exposed to a uniform stream. The geometry, the frame of reference and the notation adopted in this work are all displayed in figure 1. The origin of the Cartesian reference system is located on the left edge of the cavity with
$x$
,
$y$
and
$z$
denoting the streamwise, wall-normal and spanwise directions. The fluid motion is described by the unsteady incompressible Navier–Stokes equations,
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{u}}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}P+\frac{1}{\mathit{Re}_{BF}}{\rm\nabla}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$
$\boldsymbol{u}$
is the velocity vector with components
$\boldsymbol{u}=(u,v,w)$
and
$P$
is the reduced pressure. Equations (2.1) are made dimensionless using the cavity depth
$D$
as the characteristic length scale and the velocity of the incoming uniform stream
$U_{\infty }$
as the reference velocity. The Reynolds number is thus defined as
$Re_{BF}=U_{\infty }D/{\it\nu}$
(here the subscript
$BF$
means base flow Reynolds number) with
${\it\nu}$
the fluid kinematic viscosity. To ease comparisons, we have chosen the same boundary conditions as Sipp & Lebedev (Reference Sipp and Lebedev2007) and Barbagallo, Sipp & Schmid (Reference Barbagallo, Sipp and Schmid2009). The system of differential equations (2.1) is closed by the following Dirichlet boundary conditions at the inflow
$\partial \mathscr{D}_{in}$
and stress-free conditions at the outflow
$\partial \mathscr{D}_{out}$
: 
$$\begin{eqnarray}\displaystyle \boldsymbol{u}=1\boldsymbol{\cdot }\boldsymbol{e}_{x},\quad x\in \partial \mathscr{D}_{in};\quad P\boldsymbol{n}-\mathit{Re}^{-1}(\boldsymbol{{\rm\nabla}}\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{n}=\mathbf{0},\quad x\in \partial \mathscr{D}_{out}, & & \displaystyle\end{eqnarray}$$
where
$e_{x}$
is the unit vector in the direction of the
$x$
axis,
$n$
is the normal vector (in this case, it is the vector perpendicular to the outlet of the computational domain). Symmetric conditions (i.e.
$\partial _{y}u=0$
and
$v=0$
) are imposed at the free-stream upper boundary of the computational domain
$\partial \mathscr{D}_{ext}$
and no-slip conditions
$\boldsymbol{u}=\mathbf{0}$
at the solid walls
$\partial \mathscr{D}_{w}$
. Note that a free-slip condition with zero tangential stress (i.e.
$\partial _{y}u=0$
and
$v=0$
) is used on the walls close to the inflow and outflow
$\partial \mathscr{D}_{fs}=\partial \mathscr{D}_{fs}^{1}\cup \partial \mathscr{D}_{fs}^{2}$
.
Figure 1. Flow configuration, frame of reference and computational domain
$\mathscr{D}$
. The main features of the flow are also sketched in the figure, i.e. the boundary layer developing over the walls, the shear layer above the cavity and the recirculation forming inside the cavity of length
$L$
and depth
$D=L$
.
2.2. Global stability analysis
The flow linear instability is studied with a classical normal-mode analysis. The analysis relies on the existence of a steady solution about which infinitesimal perturbations are superimposed. The velocity and pressure fields are decomposed into a two-dimensional base flow,
$\boldsymbol{Q}_{b}(x,y)=(\boldsymbol{u}_{b},P_{b})^{\text{T}}=(u_{b},v_{b},0,P_{b})^{\text{T}}$
, and a three-dimensional disturbance flow denoted by
$\boldsymbol{q}^{\prime }(x,y,z,t)=(\boldsymbol{u}^{\prime },P^{\prime })^{\text{T}}=(u^{\prime },v^{\prime },w^{\prime },P^{\prime })^{\text{T}}$
of small amplitude
${\it\epsilon}$
. Introducing this decomposition into (2.1) and linearizing the equations governing the disturbance evolution, we obtain the two systems describing the spatial structure of the base flow and the behaviour of generally unsteady perturbations. In particular, the base flow is governed by the steady version of (2.1), whereas the perturbation field is described by the linearized unsteady Navier–Stokes equations (LNSE)
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{u}^{\prime }}{\partial t}+\unicode[STIX]{x1D647}\{\boldsymbol{u}_{b}(\mathit{Re}_{BF}),\mathit{Re}_{STB}\}\boldsymbol{u}^{\prime }=-\boldsymbol{{\rm\nabla}}P^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{\prime }=0, & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D647}$
(
$\mathit{Re}_{STB}$
indicates the Reynolds number used for stability computations) 
$$\begin{eqnarray}\unicode[STIX]{x1D647}\{\boldsymbol{u}_{b},\mathit{Re}_{STB}\}\boldsymbol{u}^{\prime }=\boldsymbol{u}_{b}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\prime }+\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{b}-\frac{1}{\mathit{Re}_{STB}}{\rm\nabla}^{2}\boldsymbol{u}^{\prime }.\end{eqnarray}$$
As the base flow is homogeneous and stationary in the spanwise direction, a generic perturbation can be decomposed into Fourier modes of spanwise wavenumber
$k$
. The three-dimensional perturbations are expressed as
$$\begin{eqnarray}\boldsymbol{q}^{\prime }(x,y,z,t)={\textstyle \frac{1}{2}}\{(\hat{u} ,\hat{v},{\hat{w}},\hat{P})(x,y)\,\,\exp [\text{i}kz+{\it\gamma}t]+\text{c.c.}\},\end{eqnarray}$$
where
${\it\gamma}={\it\eta}+\text{i}{\it\omega}$
is the complex growth rate and
$\text{c.c}$
. stands for complex conjugate. The real part
${\it\eta}$
of
${\it\gamma}$
represents the temporal growth rate of the perturbation and the imaginary part
${\it\omega}$
its frequency. For
${\it\eta}>0$
, the flow is unstable whereas for
${\it\eta}<0$
it is stable. Introducing the ansatz (2.6) in the LNSE (2.3)–(2.4), we obtain the generalized eigenvalue problem
$$\begin{eqnarray}\unicode[STIX]{x1D63C}\hat{\boldsymbol{q}}+{\it\gamma}\unicode[STIX]{x1D63D}\hat{\boldsymbol{q}}=0,\end{eqnarray}$$
in which
$\hat{\boldsymbol{q}}=(\hat{u} ,\hat{v},{\hat{w}},\hat{P})^{\text{T}}$
and
$\mathscr{A}$
is the complex linearized evolution operator. The operators
$\mathscr{A}$
and
$\mathscr{B}$
, have the following expressions:
$$\begin{eqnarray}\unicode[STIX]{x1D63C}=\left(\begin{array}{@{}cccc@{}}\mathscr{C}-\mathscr{M}+\partial _{x}u_{b} & \partial _{y}u_{b} & 0 & \partial _{x}\\ \partial _{x}v_{b} & \mathscr{C}-\mathscr{M}+\partial _{y}v_{b} & 0 & \partial _{y}\\ 0 & 0 & \mathscr{C}-\mathscr{M} & \text{i}k\\ \partial _{x} & \partial _{y} & \text{i}k & 0\end{array}\right),\quad \unicode[STIX]{x1D63D}=\left(\begin{array}{@{}cccc@{}}1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\end{array}\right),\end{eqnarray}$$
where
$\mathscr{M}={\mathit{Re}_{STB}}^{-1}(\partial _{x^{2}}+\partial _{y^{2}}-k^{2})$
and
$\mathscr{C}=u_{b}\partial _{x}+v_{b}\partial _{y}$
describe the viscous diffusion of the perturbation and its advection by the mean flow. The boundary conditions associated with the eigenproblem (2.7) are derived from those used for the base flow, i.e.
$$\begin{eqnarray}\displaystyle & \hat{\boldsymbol{u}}=\mathbf{0},\quad \text{on}~\partial \mathscr{D}_{in}\cup \partial \mathscr{D}_{w}~\text{(inlet and wall)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \hat{P}\boldsymbol{n}-{\mathit{Re}_{STB}}^{-1}(\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{u}})\boldsymbol{\cdot }\boldsymbol{n}=\mathbf{0},\quad \text{on}~\partial \mathscr{D}_{out}~\text{(outlet)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \partial _{y}\hat{u} =\hat{v}={\hat{w}}=0,\quad \text{on}~\partial \mathscr{D}_{ext}\cup \partial \mathscr{D}_{fs}~\text{(free stream and free-slip boundary)}. & \displaystyle\end{eqnarray}$$
Finally, we note that the complex-conjugate pairs
$({\it\eta}+\text{i}{\it\omega};\hat{\boldsymbol{q}})$
and
$({\it\eta}-\text{i}{\it\omega};\hat{\boldsymbol{q}}^{\ast })$
are both solutions of the eigenproblem (2.7) with the boundary conditions (2.9) for a real base flow
$\boldsymbol{Q}_{b}$
. Thus, the eigenvalues are complex conjugates and the spectra are in a symmetric plane with respect to the real axis in the
$({\it\eta},{\it\omega})$
.
2.3. Determination of the instability core: structural sensitivity
In this section we present the structural sensitivity analysis following the framework in Pralits, Brandt & Giannetti (Reference Pralits, Brandt and Giannetti2010). The underlying idea is the concept of ‘wavemaker’, introduced by Giannetti & Luchini (Reference Giannetti and Luchini2007) to identify the location of the core of a global instability. The wavemaker is the region in the flow where variations in the structure of the problem provide the largest drift of a specific eigenvalue. We first consider the perturbed eigenvalue problem
$$\begin{eqnarray}\displaystyle & {\it\gamma}^{\prime }\hat{\boldsymbol{u}}^{\prime }+\unicode[STIX]{x1D647}\{\boldsymbol{u}_{b},\mathit{Re}_{STB}\}\hat{\boldsymbol{u}}^{\prime }=-\boldsymbol{{\rm\nabla}}\hat{P}^{\prime }+{\it\delta}\unicode[STIX]{x1D643}(\hat{\boldsymbol{u}}^{\prime },\hat{P}^{\prime }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{u}}^{\prime }=0, & \displaystyle\end{eqnarray}$$
${\it\delta}\unicode[STIX]{x1D643}$
is the generalized structural perturbation. It is assumed to be a momentum force localized in space and proportional to the local velocity perturbation through a
$(3\times 3)$
coupling matrix
${\it\delta}\unicode[STIX]{x1D648}_{0}$
and a Dirac delta function: 
$$\begin{eqnarray}{\it\delta}\unicode[STIX]{x1D643}(\hat{\boldsymbol{u}}^{\prime },\hat{P}^{\prime })={\it\delta}\unicode[STIX]{x1D648}(x,y)\boldsymbol{\cdot }\hat{\boldsymbol{u}}^{\prime }={\it\delta}(x-x_{0},y-y_{0}){\it\delta}\unicode[STIX]{x1D648}_{0}\boldsymbol{\cdot }\hat{\boldsymbol{u}}^{\prime }.\end{eqnarray}$$
Neglecting higher-order terms, variations of the eigenvalue
${\it\delta}{\it\gamma}$
and of the corresponding eigenfunction
$({\it\delta}\hat{\boldsymbol{u}},{\it\delta}\hat{P})$
satisfy the following expressions:
$$\begin{eqnarray}\displaystyle & {\it\gamma}{\it\delta}\hat{\boldsymbol{u}}+\unicode[STIX]{x1D647}\{\boldsymbol{u}_{b},\mathit{Re}_{STB}\}{\it\delta}\hat{\boldsymbol{u}}=-\boldsymbol{{\rm\nabla}}{\it\delta}\hat{P}+{\it\delta}\unicode[STIX]{x1D648}\boldsymbol{\cdot }\hat{\boldsymbol{u}}-{\it\delta}{\it\gamma}\hat{\boldsymbol{u}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\it\delta}\hat{\boldsymbol{u}}=0. & \displaystyle\end{eqnarray}$$
$\hat{\boldsymbol{g}}^{+}(x,y)=(\,\hat{\boldsymbol{f}}^{+},\hat{m}^{+})$
$$\begin{eqnarray}\displaystyle & -{\it\gamma}\hat{\boldsymbol{f}}^{+}+\boldsymbol{u}_{b}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{f}}^{+}-\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{b}\boldsymbol{\cdot }\hat{\boldsymbol{f}}^{+}+\displaystyle \frac{1}{\mathit{Re}_{STB}}{\rm\nabla}^{2}\hat{\boldsymbol{f}}^{+}+\boldsymbol{{\rm\nabla}}\hat{m}^{+}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{f}}^{+}=0. & \displaystyle\end{eqnarray}$$
$\mathscr{D}$
, accounting for the boundary conditions and introducing the sensitivity tensor 
$$\begin{eqnarray}\unicode[STIX]{x1D64E}(x_{0},y_{0};\mathit{Re}_{BF},\mathit{Re}_{STB})=\frac{\hat{\boldsymbol{f}}^{+}(x_{0},y_{0})\hat{\boldsymbol{u}}(x_{0},y_{0})}{\displaystyle \int _{\mathscr{D}}\hat{\boldsymbol{f}}^{+}\boldsymbol{\cdot }\hat{\boldsymbol{u}}\text{d}S},\end{eqnarray}$$
we can express the eigenvalue drift due to the local feedback as
$$\begin{eqnarray}{\it\delta}{\it\gamma}(x_{0},y_{0})=\frac{\displaystyle \int _{\mathscr{D}}\hat{\boldsymbol{f}}^{+}\boldsymbol{\cdot }{\it\delta}\unicode[STIX]{x1D648}\boldsymbol{\cdot }\hat{\boldsymbol{u}}\text{d}S}{\displaystyle \int _{\mathscr{D}}\hat{\boldsymbol{f}}^{+}\boldsymbol{\cdot }\hat{\boldsymbol{u}}\text{d}S}=\frac{\hat{\boldsymbol{f}}^{+}\boldsymbol{\cdot }{\it\delta}\unicode[STIX]{x1D648}_{0}\boldsymbol{\cdot }\hat{\boldsymbol{u}}}{\displaystyle \int _{\mathscr{D}}\hat{\boldsymbol{f}}^{+}\boldsymbol{\cdot }\hat{\boldsymbol{u}}\text{d}S}=\unicode[STIX]{x1D64E}:{\it\delta}\unicode[STIX]{x1D648}_{0}=\mathop{\sum }_{ij}\unicode[STIX]{x1D61A}_{ij}{\it\delta}\unicode[STIX]{x1D614}_{0ij}.\end{eqnarray}$$
Different norms of the tensor
$\unicode[STIX]{x1D64E}$
can be used to build a spatial map of the sensitivity. The spectral norm is chosen here to study the worst possible case.
3. Numerical approach
3.1. Base flow calculation
The numerical computation of the base flow has been performed using a finite element method. The variational formulation of the Navier–Stokes equations (2.1) is implemented in the software package FreeFem
$++$
(Hecht (Reference Hecht2012);http://www.freefem.org) using classical
$P2{-}P1$
Taylor–Hood elements for the spatial discretization. The resultant nonlinear system of algebraic equations, along with the boundary conditions, is solved by a Newton–Raphson procedure: given an initial guess
${\boldsymbol{w}_{b}}^{(0)}$
, the linear system
$$\begin{eqnarray}NS(\mathit{Re}_{BF},{\boldsymbol{W}_{b}}^{(n)})\boldsymbol{\cdot }{\boldsymbol{w}_{b}}^{(n)}=-rhs^{(n)}\end{eqnarray}$$
is solved at each iteration step using the MUMPS (multifrontal massively parallel sparse direct solver) (Amestoy et al. Reference Amestoy, Duff, Koster and L’Excellent2001, Reference Amestoy, Guermouche, L’Excellent and Pralet2006) for the matrix inversion. The base flow is then updated as
$$\begin{eqnarray}\boldsymbol{W}_{b}^{(n+1)}={\boldsymbol{W}_{b}}^{(n)}+{\boldsymbol{w}_{b}}^{(n)}.\end{eqnarray}$$
The initial guess is chosen to be the solution of the Stokes equations and the process is continued until the
$L^{2}$
-norm of the residual of the governing equations becomes smaller than
$10^{-12}$
. To test the implementation and convergence, we used three different meshes M1, M2 and M3 (see table 1). These are generated by the bidimensional anisotropic mesh generator (Bamg) that is part of the Freefem
$++$
package. The base flow computations are also validated using a variant of the second-order finite-difference code described in Giannetti & Luchini (Reference Giannetti and Luchini2007). A typical steady flow over the open cavity is depicted in figure 2.
Table 1. Comparison of the results obtained with the present implementation and those reported by Sipp & Lebedev (Reference Sipp and Lebedev2007) for the same configuration. The eigenfrequency
${\it\omega}$
and the growth rate
${\it\eta}$
have been calculated for the first two-dimensional unstable eigenmode at
$\mathit{Re}_{BF}=\mathit{Re}_{STB}=4140$
;
$n_{d.o.f.}$
and
$n_{t}$
indicate the total number of degrees of freedom of the linearized problem and the number of triangles for each of the unstructured meshes used.
Figure 2. Visualization of the steady two-dimensional base flow for the Reynolds number
$\mathit{Re}_{BF}=1370$
at which a three-dimensional instability is first observed. The background shading displays the streamwise velocity whereas solid lines indicate the streamlines inside the cavity.
3.2. Eigenvalue solver and adjoint field
Once the base flow is determined, the system (2.7) is used to perform the stability analysis. After spatial discretization, the governing equations and their boundary conditions (2.9) are recast in the following standard form:
$$\begin{eqnarray}[\unicode[STIX]{x1D63C}(\mathit{Re}_{STB},\boldsymbol{W}_{b}(\mathit{Re}_{BF}))+{\it\gamma}\unicode[STIX]{x1D63D}]\boldsymbol{\cdot }\boldsymbol{w}=\mathbf{0},\end{eqnarray}$$
where
$\boldsymbol{w}$
is the right (or direct) eigenvector. As methods based on the
$QR$
decomposition are not feasible for solving large-scale problems like those associated with the matrix
$\unicode[STIX]{x1D63C}$
obtained for our problem, we adopt an efficient matrix-free iterative method based on the Arnoldi algorithm (Arnoldi Reference Arnoldi1951). We use the state-of-the-art ARPACK package (Lehoucq et al.
Reference Lehoucq, Maschhoff, Sorensen and Yang2007), with implicit restarts to limit memory requirements. The solution of the linear system built by the Arnoldi iterations on the Krylov subspace is obtained with the same sparse solver (Amestoy et al.
Reference Amestoy, Duff, Koster and L’Excellent2001, Reference Amestoy, Guermouche, L’Excellent and Pralet2006) as used for the base flow calculations. The adjoint modes are computed as left eigenvectors of the discrete system derived from the discretization of the linearized equations and the sensitivity function is then computed by the product of the direct and the adjoint fields. The right (direct) and left (adjoint) eigenvectors are normalized by requiring
$$\begin{eqnarray}\max _{x,y\in \mathscr{D}}\{|\hat{u} (x,y)|\}=1,\quad \int _{\mathscr{D}}\hat{\boldsymbol{f}}^{+}\boldsymbol{\cdot }\hat{\boldsymbol{u}}\text{d}S=1.\end{eqnarray}$$
Figure 3. Stability analysis of the open cavity flow at
$\mathit{Re}=4140$
, where an unstable two-dimensional mode first emerges. (a) Real and (b) imaginary part of the eigenvalue
${\it\gamma}$
versus the spanwise wavenumber
$k$
. The different branches are numbered for future reference.
The code is validated against the results reported by Sipp & Lebedev (Reference Sipp and Lebedev2007). These authors investigate the stability of a Newtonian fluid in the same geometrical configuration and report the first instability of a two-dimensional eigenmode to occur at
$\mathit{Re}=4140$
. In table 1 we present the comparison between our results and the results in Sipp & Lebedev (Reference Sipp and Lebedev2007) for different meshes. In these tests, 50 eigenvalues were obtained, with an initial Krylov basis of dimension 150; the convergence criterion for the Arnoldi iterations is based on a tolerance of
$10^{-9}$
. To independently check the accuracy of the results we a posteriori computed the residual
$\max _{i}|(\unicode[STIX]{x1D608}_{ij}+{\it\gamma}\unicode[STIX]{x1D609}_{ij})w_{j}|$
: this turns out to be always below
$10^{-9}$
for the results reported in this paper, typically less than
$10^{-12}$
for the least stable modes. The majority of the computations presented in the following are obtained using mesh M2. Henceforth whenever
$\mathit{Re}_{BF}=\mathit{Re}_{STB}$
we will simply use
$\mathit{Re}$
.
4. Linear stability results
4.1. Three-dimensional versus two-dimensional instability
As the cavity is typically considered an example of centrifugal instability, we expect the first bifurcation to be characterized by the appearance of steady three-dimensional modes of relatively short wavelength in the spanwise direction (Albensoeder, Kuhlmann & Rath Reference Albensoeder, Kuhlmann and Rath2001). To verify this, we scan the
$k$
-axis seeking for unstable modes at the Reynolds number
$Re=4140$
where a two-dimensional mode first becomes unstable (Sipp & Lebedev Reference Sipp and Lebedev2007). The results in figure 3 clearly show that eight unstable branches can be found for this value of
$Re$
where the most unstable mode has wavenumber
$k=22$
and represents a steady disturbance (
${\it\omega}=0$
in figure 3
b). The flow over an open cavity is therefore characterized by a first bifurcation to a steady three-dimensional configuration.
The full eigenvalue spectrum at
$k=22$
is shown in figure 4 where seven unstable modes appear (four branches with respect to
$k$
), placed symmetrically with respect to the real axis. It is interesting to highlight that the steady mode is the most unstable one, while the others are characterized by frequencies which are integer multiples of a fundamental one and arise (as explained in § 7) as results of resonances on a particular streamline.
Figure 4. Eigenvalue spectrum for
$Re=4140$
and
$k=22$
. The numbers in the figure relates to the branches identified in figure 3.
In figure 5 we display the modulus of both the direct and adjoint eigenfunctions corresponding to the eigenvalues denoted by 1, 2, 3, and 4 in figure 4. The velocity perturbations are most evident in the circular region inside the cavity, with a tail in the shear region just above the downstream wall. The adjoint modes, indicating the region in the flow most receptive to forcing in the momentum equations, have a similar structure, except for the thin region close to the upstream tip of the cavity, where instability is most efficiently triggered. The unstable modes are spatially localized in the same region. The secondary flow generated by the leading instability can be described as a flat roll lying within the square cavity. The different resonances are associated with periodic oscillations again concentrated in the region inside the cavity.
Figure 5. Contour plots of the absolute value of the direct and adjoint eigenfunctions of modes 1, 2, 3 and 4 as denoted in figure 4 for
$Re=4140$
and
$k=22$
. Global mode 1 is stationary (
${\it\omega}=0$
), mode 2 has
${\it\omega}\approx 0.3$
, mode 3 has
${\it\omega}\approx 0.6$
and mode 4 has
${\it\omega}\approx 0.9$
.
To document the appearance of this two-to-three-dimensional bifurcation, we determine the critical Reynolds number at which the instability first occurs: as shown in figure 6 the critical value is about
$\mathit{Re}_{cr}\approx 1370$
and the first mode to become unstable is associated with a wavenumber
$k\approx 13.4$
. All modes whose growth rate is reported in figure 6 have zero frequency; this stationary instability will be analysed in detail in the rest of the paper.
Figure 6. Leading eigenmode for cavity flow with
$L/D=1$
. Growth rate of the most unstable mode versus the spanwise wavenumber
$k$
for the Reynolds numbers indicated. The circular frequency is zero for all modes displayed. Thus, the first instability of the flow over an open cavity is a three-dimensional steady mode.
4.2. Structural sensitivity of the first bifurcation
We study the characteristics of the bifurcation by first showing the spatial structure of the fluctuation of the least stable mode at
$Re=1370$
,
$k=13.4$
. Like the modes at higher Reynolds number, the mode is localized along the external streamlines of the recirculation region inside the cavity (see figure 7). The level of fluctuations is largest in the streamwise component, the cross-stream and spanwise ones being respectively about 61 % and 88 % of the streamwise fluctuations. The adjoint of the critical mode is displayed in figure 7(d–f): its spatial structure closely resembles that of the direct mode, with a strong localization along the circular streamlines inside the cavity. As noted above, the direct mode presents a second region of noticeable amplitude near the downstream tip of the cavity and in the shear region just downstream of it, whereas the amplitude of the adjoint mode is not negligible near the upstream tip.
Figure 7. Contour plots of the streamwise (a direct; d adjoint), wall-normal (b direct; e adjoint) and spanwise (c direct; f adjoint) component of the direct and adjoint mode close to the critical Reynolds number
$Re=1370$
,
$k=13.4$
.
The structural sensitivity of this mode is displayed in figure 8. This quantity indicates the regions in the flow where a feedback forcing proportional to the local perturbation velocity most alters the eigenvalue, in other words the wavemaker of the instability. The sensitivity, product of the direct and adjoint mode, is largest inside the cavity, with no significant contributions from the regions of strong shear above it. It is interesting to note that the wavemaker is similar to that computed for a lid-driven square cavity (Giannetti, Luchini & Marino Reference Giannetti, Luchini, Marino, Schlatter and Henningson2010; Haque et al. Reference Haque, Lashgari, Giannetti and Brandt2012).
Figure 8. Structural sensitivity, the core of the instability, at the neutral conditions,
$\mathit{Re}=1370$
,
$k=13.4$
, for the flow over an open square cavity. The streamline where the asymptotic analysis predicts the maximum inviscid growth rate is also depicted (see § 7).
5. Instabilities in cavity flows
5.1. Open cavity and lid-driven cavity
The sensitivity analysis performed in the previous section clearly shows that the core of the three-dimensional instability leading to the first bifurcation in a square open cavity is highly localized in space and completely contained inside the cavity. This is in contrast with the first 2D instability (Sipp & Lebedev Reference Sipp and Lebedev2007), arising at
$Re=4140$
, which is more similar to a wake-type instability (Sipp Reference Sipp2012; Yamouni et al.
Reference Yamouni, Sipp and Jacquin2013) and localized downstream near the second tip of the cavity.
Examining the results obtained from the stability analysis, it is clear that the external flow plays little role in the generation mechanism of the three-dimensional instability. It is thus reasonable to assume that the configuration studied here is subject to the same type of instabilities as those appearing in a LDC and discussed by Albensoeder et al. (Reference Albensoeder, Kuhlmann and Rath2001) and Albensoeder & Kuhlmann (Reference Albensoeder and Kuhlmann2006) among others. In the open cavity the shear layer detaching from the upstream corner has the same role as that of the lid in the formation of the vortical motion inside the LDC configuration. The velocity of the lid is uniform while in the present study the fluid velocity along the line connecting the two upper corners is not constant: the velocity starts from zero at the left corner, increases, reaches a maximum and then decreases again, vanishing at the downstream corner. In addition, the vertical velocity is different from zero, its magnitude being smaller than that of the horizontal component, however. As a result the flow field, although qualitatively similar, also has important differences with respect to that occurring in the lid-driven problem.
For the LDC the critical Reynolds number for the first bifurcation has been calculated independently by Theofilis (Reference Theofilis2000) and Albensoeder et al. (Reference Albensoeder, Kuhlmann and Rath2001) (see also Ding & Kawahara Reference Ding and Kawahara1998; Shatrov, Mutschke & Gerbeth Reference Shatrov, Mutschke and Gerbeth2003). The numerical three-dimensional linear stability analysis of Albensoeder et al. (Reference Albensoeder, Kuhlmann and Rath2001) covers a wide range of cavity aspect ratios and presents the corresponding unstable modes, which appear to be qualitatively different when varying the cavity aspect ratio. These authors explain the centrifugal instability mechanism in terms of the perturbation energy budget and the criterion proposed by Sipp & Jacquin (Reference Sipp and Jacquin2000).
As suggested by Brés & Colonius (Reference Brés and Colonius2008) (see § 4.2 of their paper), if we introduce a (base flow) Reynolds number
$Re_{av}$
based on the cavity depth
$D$
and on the average velocity
$\tilde{U}$
along the line connecting the two opposite corners, the critical Reynolds number for the first instability of the open cavity flow becomes
$$\begin{eqnarray}\mathit{Re}_{av}=\frac{\tilde{U} D}{{\it\nu}}\approx 490,\end{eqnarray}$$
which is around 38 % lower than the value found by Albensoeder et al. (Reference Albensoeder, Kuhlmann and Rath2001). Despite this difference in the value of the critical Reynolds number, the spanwise wavenumbers at which the instability first occurs are comparable, being
$k_{lid}\approx 15.4$
in the LDC case and
$k\approx 13.4$
in the present configuration. These qualitative similarities, both in terms of base flows and modes, suggest that the same kind of instability is acting in the two configurations.
5.2. Link between open cavity flows
Brés & Colonius (Reference Brés and Colonius2008) performed DNS of open cavity flows for several
$\mathit{Re}_{{\it\delta}^{\ast }}$
, where
${\it\delta}^{\ast }$
is displacement thickness, to investigate the effect of this parameter on the instability properties. In the present configuration the shear layer starts developing at
$x^{bl}=-0.4$
(we recall that the origin of our frame of reference is located on the left edge of the cavity) leading to a displacement thickness at the upstream edge of (
$\mathit{Re}=1370$
)
$$\begin{eqnarray}{\it\delta}^{\ast }(x^{bl})=\int _{0}^{0.5}\frac{u_{b}(x^{bl},0.5)-u_{b}(x^{bl},y)}{u_{b}(x^{bl},0.5)}\text{d}y\approx 0.029.\end{eqnarray}$$
Thus, the critical Reynolds number based on this boundary-layer thickness is equal to
$\mathit{Re}_{{\it\delta}^{\ast }}=U_{\infty }{\it\delta}^{\ast }(x^{bl})/{\it\nu}\approx 39.7$
.
As discussed in Brés & Colonius (Reference Brés and Colonius2008), the critical conditions for the instability are only weakly affected by the starting position of the upstream laminar boundary layer when the critical Reynolds number is appropriately defined. The main idea is that the dynamics inside the cavity is approximately driven by the average velocity between the two edges rather than by the shear-layer thickness, thus strengthening the connection to the LDC flow. Analysing the limiting case arising when the incoming free-stream velocity
$U_{\infty }$
is simply linear, i.e. when a Couette profile of velocity
$U_{\infty }(y)=\boldsymbol{y}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{y}}$
is imposed at the inlet, we found that the first bifurcation occurs at
$\mathit{Re}^{Couette}\approx 20\,200$
for modes of spanwise wavenumber
$k\approx 13.0$
. Interestingly, the critical Reynolds number scaled with the cavity depth and average velocity,
$Re_{av}\approx 470$
, is in a very good agreement with the value obtained previously, (5.1), although the boundary-layer thickness of Couette flow is infinity. In the light of this result we confirm that the averaged Reynolds number
$\mathit{Re}_{av}$
is a relevant parameter to predict the onset of instability for open cavity flows.
6. Structural sensitivity to a velocity-based linear feedback
The so-called sensitivity to base flow variations is a concept introduced by Bottaro, Corbett & Luchini (Reference Bottaro, Corbett and Luchini2003) and Marquet, Sipp & Jacquin (Reference Marquet, Sipp and Jacquin2008) within the global framework. In this analysis a small structural velocity-based perturbation acts at the base flow level: the effect of the base flow modifications on the leading eigenvalue of the stability problem allows us to study the different mechanisms that can suppress or enhance the instability. The spatial structure of the so-called adjoint base flow can be used to identify the features of the base flow that provide the main contribution to the instability dynamics and the regions where to locate effective passive control devices. In other words, this modification of the structure of the Navier–Stokes operator causes a variation of the base flow which in turn produces a drift of the leading eigenvalue
${\it\gamma}={\it\eta}+\text{i}{\it\omega}$
.
For the sake of brevity, only the main ingredients are outlined here; an extensive and detailed derivation can be found in Marquet et al. (Reference Marquet, Sipp and Jacquin2008) and Pralits et al. (Reference Pralits, Brandt and Giannetti2010). Using a formalism based on control theory, the eigenproblem (2.7) represents the state equation, the state vector is composed of the global mode
$\hat{\boldsymbol{q}}$
and the complex eigenvalue
${\it\gamma}$
, and the base flow
$\boldsymbol{Q}_{b}$
is the control variable. As in Pralits et al. (Reference Pralits, Brandt and Giannetti2010) we express the eigenvalue drift
${\it\delta}{\it\gamma}$
as
$$\begin{eqnarray}{\it\delta}{\it\gamma}={\it\delta}{\it\eta}+\text{i}{\it\delta}{\it\omega}=\frac{\displaystyle \int _{\mathscr{D}}(\hat{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{f}}^{+}-\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{u}}\boldsymbol{\cdot }\hat{\boldsymbol{f}}^{+})\boldsymbol{\cdot }{\it\delta}\,\boldsymbol{u}_{b}\text{d}S}{\displaystyle \int _{\mathscr{D}}\hat{\boldsymbol{f}}^{+}\boldsymbol{\cdot }\hat{\boldsymbol{u}}\text{d}S},\end{eqnarray}$$
where
${\it\delta}\boldsymbol{u}_{b}$
is a generic modification of the base flow. The relation (6.1) provides the effect of a specified velocity distribution implying a dedicated computation for each specific actuation, e.g. wall blowing/suction. The optimal boundary velocity distribution, instead, can be directly found as follows (see Meliga & Chomaz Reference Meliga and Chomaz2011):
$$\begin{eqnarray}{\it\delta}{\it\gamma}=\frac{\displaystyle \int _{\partial \mathscr{D}_{w,i,fs}}\left[\left(m_{b}^{+}\boldsymbol{n}+\displaystyle \frac{1}{\mathit{Re}}\boldsymbol{n}^{\text{T}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{f}_{b}^{+}\right)\boldsymbol{\cdot }{\it\delta}\,\boldsymbol{U}_{w,i,fs}\,\right]\text{d}l}{\displaystyle \int _{\mathscr{D}}\hat{\boldsymbol{f}}^{+}\boldsymbol{\cdot }\hat{\boldsymbol{u}}\text{d}S}\end{eqnarray}$$
where
$m_{b}^{+}$
and
${\hat{\boldsymbol{f}}_{b}}^{+}$
are the adjoint base flow pressure and the three-dimensional adjoint base flow velocity field and the subscripts
$w,i,fs$
indicate the boundaries (inlet, wall or free-slip) on which we calculate the integral
$\text{d}l$
is the length element along
$\partial D_{w,i,fs}$
. The adjoint base flow field
$\boldsymbol{Q}_{b}^{+}$
must satisfy the following set of linear equations (Pralits et al.
Reference Pralits, Brandt and Giannetti2010):
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{b}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{f}_{b}^{+}-\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{b}\boldsymbol{\cdot }\boldsymbol{f}_{b}^{+}+\frac{1}{\mathit{Re}}{\rm\nabla}^{2}\boldsymbol{f}_{b}^{+}+\boldsymbol{{\rm\nabla}}m_{b}^{+}=\hat{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{f}}^{+}-\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{u}}\boldsymbol{\cdot }\hat{\boldsymbol{f}}^{+}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{f}_{b}^{+}=0, & \displaystyle\end{eqnarray}$$
$m_{b}^{+}\boldsymbol{n}-Re^{-1}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{f}_{b}^{+}=-(\boldsymbol{u}_{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{f}_{b}^{+}+(\hat{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{n})\hat{\boldsymbol{u}}^{+}$
at
$\partial \mathscr{D}_{out}$
and zero-velocity conditions at the solid walls and at the inlet.
Figure 9. Sensitivity of the first bifurcation to streamwise and wall-normal mean velocity modifications at the inflow. The profile shown would provide the largest possible stabilization of the first bifurcation.
The sensitivity of the instability with respect to the incoming flow is examined first. Figure 9 shows the sensitivity to both the streamwise and the wall-normal components of the inflow velocity profile, where the profiles shown would provide the optimal decrease of the instability growth rate. The
$x$
-component is found to be always negative and attains significant values only near the wall. This fact is not surprising because the base flow modifications have effect only if related to the shear layer that drives the core vortex inside the cavity. Perturbations in the free stream do not affect the flow at the edge and inside the cavity, the regions where the instability is triggered. Negative modifications of the inlet velocity profile cause stabilization due to decrease of the momentum inside the shear layer. The effect of the wall-normal component is related to the same mechanism, decreasing of the total streamwise momentum at the cavity tip by normal advection.
In view of an active control of the first bifurcation, we depict the wall-normal component of the sensitivity along the cavity walls in figure 10; this corresponds to the optimal blowing/suction profiles giving the largest stabilization. We see that the sensitivity is vanishing along the first free-slip boundary just downstream of the inflow (
$y=0;-1.2<x<-0.4$
), while a combination of blowing and suction is found to be optimal on the wall upstream of the cavity (
$y=0;-0.4<x<0$
). Inside the cavity, on both lateral walls, a stabilizing normal component is directed in the streamwise direction. The optimal blowing and suction on the lower wall would create a flow opposite to the vortex inside the cavity, thus trying to quench it. The analysis, finally, shows that it is not possible to significantly modify the instability by applying control on the downstream wall (lowest sensitivity magnitude).
Figure 10. Sensitivity to wall-normal blowing and suction of the leading-mode growth rate at
$\mathit{Re}=1370$
and
$k=13.4$
. Positive (negative) values indicate blowing (suction) and lead to stabilization (destabilization). The profiles shown on each wall provide the optimal decrease of the instability growth rate.
7. Asymptotic inviscid stability theory
The spatial distribution of the structural sensitivity (that is spatially concentrated around a streamline inside the cavity) suggests the possibility of using the local theory to describe the evolution of the instability and provide more quantitative evidence for the mechanism from which it arises. An appealing approach in this context is offered by the short-wavelength approximation (WKBJ) developed by Bayly (Reference Bayly1988).
This approach is briefly outlined here; for a more detailed presentation the reader is referred to Lifschitz & Hameiri (Reference Lifschitz and Hameiri1991), Lifschitz (Reference Lifschitz1994) and references therein. The solution of the linearized Navier–Stokes equations is sought in the form of a rapidly oscillating and localized wave-packet evolving along the Lagrangian trajectory
$\boldsymbol{X}(t)$
and characterized by a wave-vector
$\boldsymbol{k}(t)=\boldsymbol{{\rm\nabla}}{\it\phi}(\boldsymbol{X},t)$
and an envelope
$\boldsymbol{a}(\boldsymbol{X},t)$
such that
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}(\boldsymbol{X},t)=\text{e}^{\text{i}{\it\phi}(\boldsymbol{X},t)/{\it\epsilon}}\boldsymbol{a}(\boldsymbol{X},t,{\it\epsilon})=\text{e}^{\text{i}{\it\phi}(\boldsymbol{X},t)/{\it\epsilon}}\mathop{\sum }_{n}\boldsymbol{a}_{n}(\boldsymbol{X},t){\it\epsilon}^{n}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p(\boldsymbol{X},t)=\text{e}^{\text{i}{\it\phi}(\boldsymbol{X},t)/{\it\epsilon}}b(\boldsymbol{X},t,{\it\epsilon})=\text{e}^{\text{i}{\it\phi}(\boldsymbol{X},t)/{\it\epsilon}}\mathop{\sum }_{n}b_{n}(\boldsymbol{X},t){\it\epsilon}^{n+1}, & \displaystyle\end{eqnarray}$$
${\it\epsilon}\ll 1$
and
$\boldsymbol{X}={\it\epsilon}\boldsymbol{x}$
is a slowly varying variable. In the limit of vanishing viscosity (
$\mathit{Re}\rightarrow \infty$
) and large wavenumbers (
$\Vert \boldsymbol{k}\Vert \rightarrow \infty$
), the theory provides the leading-order term for the growth rate associated with a localized perturbation. This is obtained by integrating the following set of ordinary differential equations (ODEs): 
$$\begin{eqnarray}\displaystyle & {\displaystyle \frac{\text{D}\boldsymbol{k}}{\text{D}t}}=-\unicode[STIX]{x1D647}^{t}(\boldsymbol{X})\boldsymbol{k}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\displaystyle \frac{\text{D}\boldsymbol{a}}{\text{D}t}}=\left({\displaystyle \frac{2\boldsymbol{k}\boldsymbol{k}^{\text{T}}}{|\boldsymbol{k}|^{2}}}-\unicode[STIX]{x1D644}\right)\unicode[STIX]{x1D647}(\boldsymbol{X})\boldsymbol{a}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}{\displaystyle \frac{\text{D}\boldsymbol{X}(t)}{\text{D}t}}=\boldsymbol{u}_{b}(X(t),t).\end{eqnarray}$$
In the equations above
$\unicode[STIX]{x1D647}=\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{b}$
is the base flow velocity gradient tensor and
$\unicode[STIX]{x1D644}$
the identity matrix. Since the flow under investigation is steady, the Lagrangian trajectory corresponds to the streamlines of the base flow. Three initial conditions have to be assigned to solve the problem above:
$\boldsymbol{k}(t=0)=\boldsymbol{k}_{0}$
,
$\boldsymbol{a}(t=0)=\boldsymbol{a}_{0}$
and
$\boldsymbol{x}(t=0)=\boldsymbol{x}_{0}$
. The last condition imposes the Lagrangian origin of the streamline and thereby entirely identifies it.
Lifschitz & Hameiri (Reference Lifschitz and Hameiri1991) proved that a sufficient condition for inviscid instability is that the system (7.3)–(7.5) has at least one solution for which
$\Vert \boldsymbol{a}(t)\Vert \rightarrow \infty$
as
$t\rightarrow \infty$
. This theory has been successfully applied in the past to study elliptic, hyperbolic and centrifugal instabilities of two-dimensional stationary base flows (Sipp et al.
Reference Sipp, Lauga and Jacquin1999; Godeferd, Cambon & Leblanc Reference Godeferd, Cambon and Leblanc2001). In order to characterize the instability mechanism arising inside the cavity with this local theory, the self-excited nature of the instability must be properly accounted for. In this context, a central role is played by closed Lagrangian trajectories (closed streamlines in our case), i.e. orbits described by material points which return to their initial positions after a given time
$T$
(the period of revolution of a material particle). These closed trajectories play a special role in the dynamics of the instability: on the closed orbits, local instability waves propagate and feed back on themselves leading to a self-excited unstable mode.
To apply the theory, both (7.3) and (7.4) must be integrated along the closed orbits existing inside the cavity. Since the base flow is steady and the streamlines are closed, (7.3) is a linear ODE with periodic coefficients whose general solution can be written in terms of Floquet modes. In particular, the solution can be found by building the fundamental Floquet matrix
$\unicode[STIX]{x1D648}(T)$
, solution of the system
$$\begin{eqnarray}{\displaystyle \frac{\text{D}\unicode[STIX]{x1D648}}{\text{D}t}}=-\unicode[STIX]{x1D647}^{t}(\boldsymbol{X})\unicode[STIX]{x1D648}\quad \text{with}~\unicode[STIX]{x1D648}(0)=\unicode[STIX]{x1D644},\end{eqnarray}$$
and extracting its eigenvalues and the corresponding eigenvectors. Using these eigenvectors as initial conditions, it is possible to retrieve the temporal evolution of
$\boldsymbol{k}$
during a lap around the closed streamline. Equation (7.3) admits three independent solutions related to the three eigenvectors of the fundamental Floquet matrix
$\unicode[STIX]{x1D648}(T)$
. However, since the base flow is two-dimensional, there exists for each orbit one eigenvalue equal to one, with the corresponding eigenvector remaining constant in time and orthogonal to the base flow. In other words, since the third column of
$\unicode[STIX]{x1D647}$
and the third line of
$\unicode[STIX]{x1D647}^{t}$
are zero, the transverse component of
$\boldsymbol{k}$
remains constant as time evolves. In contrast, the in-plane components evolve under the action of the deformation tensor. Once (7.3) is solved, the amplitude
$\boldsymbol{a}$
can be found by integrating (7.4). One can use any linear combination of the Floquet modes from (7.6) to set the specific
$\boldsymbol{k}$
in (7.4).
Since we are trying to determine a self-excited mode, we need only to consider solutions of (7.3) that are periodic in time, i.e. solutions such that
$\boldsymbol{k}(0)=\boldsymbol{k}(T)$
. Moreover Bayly (Reference Bayly1988), Lifschitz & Hameiri (Reference Lifschitz and Hameiri1991) and Sipp & Jacquin (Reference Sipp and Jacquin2000) have shown that centrifugal and hyperbolic instabilities attain their maximum growth rate for modes characterized by purely transverse wavenumbers. Therefore, only eigenvectors orthogonal to the base flow will be considered in the following analysis. Solutions of (7.4) associated with a
$\boldsymbol{k}$
orthogonal to the plane of motion are usually termed pressureless modes (see also Godeferd et al.
Reference Godeferd, Cambon and Leblanc2001). With this choice, (7.4) reduces to an ordinary linear differential equation with periodic coefficients. According to Floquet theory, its solution can be written in terms of Floquet modes
$$\begin{eqnarray}\boldsymbol{a}(t)=\bar{\boldsymbol{a}}(t)\exp ({\it\sigma}t),\end{eqnarray}$$
where
$\bar{\boldsymbol{a}}(t)$
is a periodic function (with the same period
$T$
as the material point moving along the selected closed streamline) and
$\text{Re}\{{\it\sigma}\}={\it\sigma}_{r}$
is the growth rate of the perturbation. In order to make a quantitative comparison with the eigenvalues predicted by the global analysis, we have to compute the values of
${\it\sigma}$
in (7.7) for each closed orbit inside the cavity. To this end, we parameterize each streamline, and the corresponding growth rate
${\it\sigma}$
, with the distance along the horizontal line connecting the centre of the vortex to the left-hand wall of the cavity (see figure 11).
As for (7.3), the fundamental Floquet matrix
$\unicode[STIX]{x1D63C}$
corresponding to (7.4) is built by integrating the system
$$\begin{eqnarray}\displaystyle & {\displaystyle \frac{\text{D}\unicode[STIX]{x1D63C}}{\text{D}t}}=\left({\displaystyle \frac{2\boldsymbol{k}\boldsymbol{k}^{\text{T}}}{|\boldsymbol{k}|^{2}}}-\unicode[STIX]{x1D644}\right)\unicode[STIX]{x1D647}(\boldsymbol{X})\unicode[STIX]{x1D63C}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D63C}(0)=\unicode[STIX]{x1D644}, & \displaystyle\end{eqnarray}$$
${\it\mu}_{i}(x_{0})$
and the corresponding eigenvectors of
$\unicode[STIX]{x1D63C}(T)$
are then easily extracted.
Figure 11. Streamlines for the flow inside the cavity at
$\mathit{Re}=1370$
. The asymptotic inviscid stability theory identifies three streamlines along which three different WKBJ modes present their maximum growth rate: - - - -, orbit of the global maximum growth rate pertaining to the steady mode
${\it\sigma}_{\infty }$
(see figure 12); — ⋅ —, orbit of the unsteady mode related to
${\it\sigma}_{2}$
;
$\cdots \cdots$
, orbit of the steady mode related to
${\it\sigma}_{3}$
. The evolution of a particle along the streamline - - - - is also depicted. The revolution period of this streamline is
$T=18.3$
. The horizontal line connecting the centre of the vortex to the left-hand wall of the cavity is used in the present work to parameterize the streamlines.
As mentioned above, since the base flow is two-dimensional and the wave-vector
$\boldsymbol{k}$
is orthogonal to the
$x{-}y$
plane, we expect one eigenvalue of
$\unicode[STIX]{x1D63C}$
to be 1. The other two, for the incompressibility constrain, must multiply to 1, i.e.
${\it\mu}_{1}(x_{0})~{\it\mu}_{2}(x_{0})=1$
. The Floquet exponent
${\it\sigma}(x_{0})$
of the perturbation on the selected orbit
${\it\psi}_{0}$
is obtained from the Floquet multiplier
${\it\mu}(x_{0})$
of
$\unicode[STIX]{x1D63C}$
by the simple relation
$$\begin{eqnarray}{\it\sigma}^{\{n\}}({\it\psi}_{0})={\it\sigma}_{r}({\it\psi}_{0})+\text{i}{\it\sigma}_{i}^{\{n\}}({\it\psi}_{0})={\displaystyle \frac{\log ({\it\mu})}{T({\it\psi}_{0})}}+\text{i}{\displaystyle \frac{2n{\rm\pi}}{T({\it\psi}_{0})}}\quad \text{with}~n\in \mathbb{N},\end{eqnarray}$$
where
$T({\it\psi}_{0})$
is the period of revolution.
The growth rate of each WKBJ mode is simply given by the real part of
${\it\sigma}^{\{n\}}$
. The frequency is related to the imaginary part and is not unique. According to the formula (7.10), modes with the same growth rate (at leading order) but different frequencies are admissible: in particular the admissible frequencies are integer multiple of the frequency of revolution along the same streamline.
7.1. Asymptotic estimate of the first bifurcation
The numerical computations of the asymptotic stability are performed on the same base flow fields as used for the global stability analysis. Several numerical methods are available to solve the system of ODEs (7.5)–(7.8) along with their initial conditions. We chose a fourth-order Runge–Kutta method: starting from the points located on the horizontal line, connecting the centre of the vortex to the left-hand wall of the cavity (see figure 11), the algorithm marches along the orbits ensuring the spatial periodicity of each streamline. In the figure we also report the position of a material point along its trajectory at equal time intervals to give a visual impression of the local velocity along the streamline.
The asymptotic eigenpairs have been computed with several discretizations and only the eigenvalues with an accuracy of four significant digits are presented. In figure 12, we show the real and imaginary part of the eigenvalues obtained with the WKBJ approximation as function of the
$x$
coordinate defining the different orbits.
Figure 12. WKBJ growth rate
${\it\sigma}_{r}$
(
$\ast$
) and eigenfrequency
${\it\sigma}_{i}$
(○) at
$\mathit{Re}_{BF}=1370$
. Here the parameter
$x$
is the physical coordinate showed in figure 11. The lines denote the maximum of the different branches: - - - -, maximum inviscid growth rate
${\it\sigma}_{r\,\infty }$
; — ⋅ —, maximum related to second branch
${\it\sigma}_{r\,2}$
;
$\cdots \cdots$
, third maximum
${\it\sigma}_{r\,3}$
.
The asymptotic analysis reveals three maxima of the growth rates
${\it\sigma}_{r\,\infty }(={\it\sigma}_{r\,1})$
, which is also the global maximum,
${\it\sigma}_{r\,2}$
and
${\it\sigma}_{r\,3}$
. The first and the third branches (
${\it\sigma}_{1}$
and
${\it\sigma}_{3}$
) are characterized by zero-frequency eigenvalues, while the second branch (
${\it\sigma}_{2}$
) is associated with unstable oscillations with frequency of
${\approx}0.2$
. As further discussed below, the closed streamline of maximum growth rate
${\it\sigma}_{r\,\infty }$
is located within the wavemaker of stationary unstable global mode.
However, despite this agreement, the viscous correction term and the correction term relative to finite wavenumber effects need to be taken into account for a correct prediction of the instability, see Landman & Saffman (Reference Landman and Saffman1987) and Gallaire et al. (Reference Gallaire, Marquillie and Ehrenstein2007).
In figure 13(a) we report the growth rate of the unstable mode computed on the base flow at
$\mathit{Re}_{BF}=1370$
when increasing the Reynolds number in the linearized stability equations,
$\mathit{Re}_{STB}$
, and the growth rate obtained by integrating along the closed orbits with the corrections discussed in appendix A,
$$\begin{eqnarray}s={\it\sigma}({\it\psi}_{0})-\frac{A}{k}-\frac{k^{2}}{\mathit{Re}_{BF}}.\end{eqnarray}$$
The value of
$A$
above is not estimated by a least square fitting as in previous studies, but computed analytically using the information provided by the local adjoint and direct field on the streamline. The values obtained with this procedure are reported in table 2. Figure 13(a) shows that the scaling provided by the global stability analysis estimates correctly the asymptotic growth rate
${\it\sigma}_{r\,\infty }$
. The corresponding optimal spanwise wavenumber
$k$
is depicted in 13(b) as a function of
$\mathit{Re}_{STB}^{1/3}$
. The spanwise wavenumber, like the maximal growth rate, follows the correct scaling laws,
${\it\sigma}_{r}\propto \mathit{Re}_{STB}^{-1/3}$
and
$k_{opt}\propto \mathit{Re}_{STB}^{1/3}$
(Bayly Reference Bayly1988; Sipp et al.
Reference Sipp, Lauga and Jacquin1999).
Figure 13. Global and asymptotic stability results. (a) Viscous growth rate
${\it\eta}$
(○, ▫
$=\mathit{Re}_{STB}=300\,000$
) at
$\mathit{Re}_{BF}=1370$
and asymptotic estimate of growth rate
${\it\sigma}_{r}$
(
$\ast$
- - - -) according to the correction in appendix A. We also depict the regression line (——) related to the global growth rates. (b) Optimal global spanwise wavenumber
$k_{opt}$
(○) and prediction from the asymptotic theory (
$\ast$
) as a function of
$\mathit{Re}_{STB}^{1/3}$
. The predicted optimal spanwise wavenumber is simply obtained by finding the maximum of the scaling law (A 6), i.e.
$k=(\mathit{Re}_{STB}A/2)^{1/3}$
.
Table 2. Convergence of parameters arising in the asymptotic estimation (see appendix A) of the viscous growth rate;
${\rm\Delta}t_{orbit}$
is the step used to discretize the critical orbit. (Here,
$\mathit{Re}_{BF}=1370$
).
Finally, we focus our attention on the spatial distribution of the structural sensitivity fields computed with the maximum
$\mathit{Re}_{STB}$
considered (equal to 300 000). Figure 14(a) shows the agreement between the critical streamline (i.e. the streamline
${\it\psi}$
where the inviscid growth rate is maximum) and the sensitivity map. At large (stability) Reynolds numbers
$\mathit{Re}_{STB}$
, therefore, the sensitivity analysis indicates that the instability core is located on the orbit with maximum growth rate.
Figure 14. Comparison between the optimal streamline of WKBJ branches and sensitivity maps related to (global) eigenvalues: (a)
$0.36+\text{i}0.00$
; (b)
$0.14+\text{i}0.20$
; (c)
$0.07+\text{i}0.00$
. Parameter settings:
$\mathit{Re}_{BF}=1370$
,
$\mathit{Re}_{STB}=300\,000$
and
$k=93.5$
.
The global analysis performed at
$\mathit{Re}_{BF}$
also provides information about the sub-critical branches arising in the asymptotic computations. We depict the structural sensitivity extracted from the global analysis of these two sub-critical WKBJ eigenmodes in figure 14(b,c). As for the leading eigemode we observe an excellent correspondence between the sensitivity spatial map and the two critical orbits. Interestingly, we note also the agreement between the frequency of mode
${\it\sigma}_{2}$
(see figure 14
b) and the frequency predicted by the WKBJ analysis. From a physical point of view, this matching can be associated with the fact that these eigenmodes are of centrifugal nature, i.e. inviscid, and therefore the inviscid structural sensitivity is able to isolate accurately the regions where each of the three instability branches presents the main contribution to the instability mechanism.
7.2. Asymptotic results for
$\mathit{Re}_{BF}=4140$
As previously discussed, when we consider the stability to three-dimensional perturbations at supercritical Reynolds numbers, we find several unstable branches (see figure 3). If we consider the spanwise wavenumber
$k=22$
and
$\mathit{Re}_{BF}=4140$
, we observe the occurrence of several harmonics of the fundamental leading eigenvalue. As shown in figure 3(b), these modes are characterized by a quantized eigenfrequency,
${\it\omega}\approx 0.32n$
with
$n$
integer. To show that the asymptotic analysis is also able to accurately predict the frequency of these harmonics, we carry out computations for the base flow at Reynolds number
$\mathit{Re}_{BF}=4140$
and report the results in figure 15, using the same conventions used for the onset of the bifurcation, for
$\mathit{Re}_{BF}=1370$
, in figure 12. We first need to identify the closed streamlines and then calculate the instability properties along the orbit. We observe again three local maxima of the asymptotic growth rate, corresponding to two steady and one time-dependent modes.
Figure 15. Asymptotic results for
$\mathit{Re}_{BF}=4140$
. See figure 12 for details.
Figure 16. Revolution period of cavity orbits as a function of the abscissa
$x$
(upper half). Contour plot of the base flow modulus (
$\Vert \boldsymbol{u}_{b}\Vert _{2}$
) at
$\mathit{Re}_{BF}=4140$
and the three streamlines corresponding to the local maxima of the inviscid growth rate in figure 15 (lower half).
The variation of the revolution period
$T$
as a function of the coordinate
$x$
, defining the different orbits, is depicted in the upper half of figure 16, while the corresponding orbits inside the cavity are displayed in the lower half. The main result we present here is that the period of the higher harmonics of the zero-frequency leading mode is selected by the period of revolution along the streamline of maximum growth rate. In table 3 we show that the frequencies obtained from the global stability analysis and displayed in figure 3 do indeed correspond to the frequency computed by the local analysis (7.10). It is interesting to note that the global mode frequencies are uniquely related to the revolution period of a Lagrangian particle transported along the orbit. Thus, we conclude that the different frequencies of the multiple unstable branches are obtained as multiples
$n$
of the period of revolution along the critical (most unstable) orbit; the data in the table show an error lower than 5 %.
Table 3. Comparison of the results obtained using global stability analysis with those provided by the asymptotic analysis. We selected the orbit that has the maximum inviscid growth rate using a base flow characterized by
$\mathit{Re}_{BF}=4140$
. The WKBJ frequencies are calculated according to (7.10).
8. Final remarks
In this work, we study the instability of the flow past an infinitely wide open square cavity. First, we identify the critical Reynolds number (
$\mathit{Re}_{BF}=1370$
) at which the first bifurcation occurs. This instability drives the flow from a steady two-dimensional to a steady three-dimensional configuration characterized by a relatively short modulation in the spanwise direction, the spanwise scale of the modulation being of about 0.47 cavity depths.
The spatial structure of the direct and adjoint eigenmodes is examined to describe the features of the flow past the bifurcation. The direct mode is concentrated inside the cavity in a circular region with a tail on the shear region just above the downstream wall; the adjoint mode has a similar structure except for a small region near the upstream edge of the cavity where the flow is most receptive to momentum forcing. The overlapping of these two fields provides information about the instability mechanism (the so-called wavemaker) and is concentrated within the square cavity, suggesting that the generation of the instability mechanism is spatially concentrated around a closed streamline inside the cavity, around the core vortex.
We examine different types of cavity flows (i.e. characterized by different boundary conditions) and identify the corresponding critical Reynolds number at which the first bifurcation occurs. As proposed by Brés & Colonius (Reference Brés and Colonius2008), we confirm that the mean velocity computed along the line connecting the two opposite edges allows us to roughly estimate the critical conditions at which the first bifurcation arises. The critical value of a Reynolds number based on this averaged velocity, the cavity depth and the fluid viscosity is found to be
$\mathit{Re}_{av}\approx 470$
.
The sensitivity to base flow modifications is then considered to study the mechanisms that can suppress or enhance the instability. We follow here the approach by Meliga & Chomaz (Reference Meliga and Chomaz2011) and compute the optimal linear velocity distribution at the walls and at the inlet of the computational domain able to stabilize the flow. The resulting blowing/suction profiles show that each modification (when possible) is aimed at decreasing the total momentum of the cavity core vortex (identified above as the core of the instability).
The WKBJ approximation is then introduced to predict the first instability and its characteristics as suggested by Bayly (Reference Bayly1988). Considering the asymptotic stability along the closed streamlines inside the cavity, we find three different branches of unstable orbits (two stationary branches and an unsteady branch) and select the three critical orbits
${\it\psi}$
whose corresponding growth rates are local maxima (
${\it\sigma}_{r\,\infty },{\it\sigma}_{r\,2},{\it\sigma}_{r\,3}$
). The asymptotic values of the growth rate and of the spanwise wavenumber of the unstable modes show very good agreement with the global stability analysis once the correction for finite Reynolds number and spanwise length scale are applied to the inviscid asymptotic result. The three critical orbits detected by the asymptotic analysis are also shown to overlap with the structural sensitivity map of unstable modes at low viscosities, large
$\mathit{Re}_{STB}$
(we refer to this field as inviscid structural sensitivity). This procedure allows us to identify the spatial region where the core of the inviscid mechanism of the instability is located.
To identify a frequency selection mechanism for the time-dependent sub-leading unstable global modes emerging at supercritical conditions, we consider the stability of the flow at
$\mathit{Re}_{BF}=4140$
where the global analysis shows the occurrence of four branches of unstable modes characterized by frequencies that are multiples of a fundamental value
${\it\omega}_{0}$
. We show that the value of
${\it\omega}_{0}$
corresponds to the period of revolution of Lagrangian fluid particles along the closed orbit of largest growth rate in the asymptotic limit. We thus conclude that the asymptotic theory is able to predict accurately the global stability results, enabling us to estimate the critical conditions leading to the instability. Furthermore, the inviscid structural sensitivity, discussed here, is a general concept that can be used whenever the instability is of inviscid type.
Appendix A. Construction of pressureless modes
In what follows we briefly recall the theory related to the dynamics of asymptotic modes (Bayly Reference Bayly1988). First of all, we express the evolution of the perturbation using the normal-mode ansatz
$$\begin{eqnarray}\displaystyle [\boldsymbol{u}^{\prime },P^{\prime }]=[\hat{\boldsymbol{u}},\hat{P}]\exp \{\text{i}kz+st\}. & & \displaystyle\end{eqnarray}$$
The main idea is to use the eigenpairs of the fundamental Floquet matrix
$\unicode[STIX]{x1D648}(T)$
to build a vector basis
$\boldsymbol{f}_{i}$
for the representation of the modes along the orbit:
$$\begin{eqnarray}\hat{\boldsymbol{u}}(\boldsymbol{x})=\mathop{\sum }_{i=1}^{3}\hat{u} _{i}(\boldsymbol{x})\boldsymbol{f}_{i}(\boldsymbol{x}).\end{eqnarray}$$
This basis diagonalizes the nonlinear operator
$\boldsymbol{u}_{b}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(\,)+(\,)\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{b}$
and can be computed as
$\boldsymbol{f}_{i}=\text{e}^{-{\it\sigma}t}\unicode[STIX]{x1D648}(t)\boldsymbol{e}_{\boldsymbol{i}}$
. Considering the limit of
$\Vert \boldsymbol{k}\Vert \rightarrow \infty$
, we re-scale the WKBJ eigenmode as
$$\begin{eqnarray}[\hat{\boldsymbol{u}},\hat{\boldsymbol{v}},\hat{\boldsymbol{w}},\hat{\boldsymbol{p}}](\tilde{{\it\Psi}})=[\hat{\boldsymbol{U}},\hat{\boldsymbol{V}}/k,\hat{\boldsymbol{w}}/\sqrt{k},\hat{\boldsymbol{p}}/k\sqrt{k}](\tilde{{\it\Psi}}),\end{eqnarray}$$
where the new streamfunction
$\tilde{{\it\Psi}}=\sqrt{k}({\it\psi}-{\it\psi}_{0})$
allows us to magnify the region near the critical orbit
${\it\psi}_{0}$
. Introducing the scaling (A 3) into the LEes, we get the equation of a quantum harmonic oscillator (see e.g. Bender & Orszag Reference Bender and Orszag1978)
$$\begin{eqnarray}\hat{\boldsymbol{U}}^{\prime \prime }(\tilde{{\it\Psi}})+\left[\frac{A}{\mathscr{J}({\it\psi}_{0})}-{\it\lambda}^{2}\tilde{{\it\Psi}}^{2}\right]\hat{\boldsymbol{U}}=0,\end{eqnarray}$$
with
$\hat{\boldsymbol{U}}(\pm \infty )=0$
,
${\it\lambda}^{2}=-{\it\sigma}^{\prime \prime }({\it\psi}_{0})/(2\mathscr{J})$
and
$$\begin{eqnarray}\displaystyle \mathscr{J}=\frac{1}{T({\it\psi}_{0})}\int _{0}^{T({\it\psi}_{0})}(\,\boldsymbol{f}_{1}^{\dagger }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\psi})\left[{\it\sigma}({\it\psi}_{0})+\frac{\text{d}}{\text{d}t}\right]\{\,\boldsymbol{f}_{1}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\psi}\}\text{d}t. & & \displaystyle\end{eqnarray}$$
In (A 4), the constant
$A$
is the parameter that governs the scaling of the eigenvalue
$s$
(i.e.
$s={\it\sigma}({\it\psi}_{0})-A/k$
) and the adjoint vector
$\boldsymbol{f}_{i}^{\dagger }$
is normalized as
$\boldsymbol{f}_{i}^{\dagger }\boldsymbol{\cdot }\boldsymbol{f}_{j}={\it\delta}_{ij}$
.
A better quantitative estimate of the viscous growth rate can be achieved using the viscous correction introduced by Landman & Saffman (Reference Landman and Saffman1987) (see also Gallaire et al. Reference Gallaire, Marquillie and Ehrenstein2007). The composite estimation thus reads
$$\begin{eqnarray}s={\it\sigma}({\it\psi}_{0})-\frac{A}{k}-\frac{k^{2}}{\mathit{Re}_{BF}}.\end{eqnarray}$$
