2.1. Geometrical configuration

The analysis presented here concerns the flow of a Newtonian fluid through a symmetric sudden expansion. With reference to figure 1, the flow is coming from the left at the inlet of a channel of height $d$ and length $L_{\textit{in}}=5d$ , followed by a sudden enlargement in a planar duct of height $D=3d$ , whose length $L_{\textit{out}}=50d$ has been chosen to recover a plane Poiseuille flow at the downstream end of the channel. Let $\textit{ER}=D/d$ be the ratio of the height of the channel downstream ( $D$ ) and upstream ( $d$ ) the geometrical discontinuity. In this study, a value of the expansion ratio $\textit{ER}=3$ has been considered (as, for instance, in Fearn et al. Reference Fearn, Mullin and Cliffe1990; Hawa & Rusak Reference Hawa and Rusak2000; Fani et al. Reference Fani, Camarri and Salvetti2012, among others).

Figure 1. Symmetric sudden expansion geometry. Dashed lines qualitatively represent the recirculation regions.

The origin of the Cartesian reference frame is placed at the midpoint of the inlet edge ( $\partial \mathcal{D}_{\textit{in}}$ ), denoting the streamwise, wall-normal and spanwise directions with $x,y$ and $z$ , respectively. The base flow is assumed to be two-dimensional and homogeneous in the spanwise direction $(z)$ .

The separation and reattachment points of the three recirculation regions (that characterise the secondary asymmetric flow field) are also qualitatively represented in figure 1.

2.2. Problem statement

The fluid motion is governed by the unsteady incompressible Navier–Stokes equations that, in non-dimensional form, can be written as follows:

(2.1a) \begin{align} \frac {\partial \boldsymbol{u}}{\partial t} + \left ( \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol{u} & = -\boldsymbol{\nabla }P + \frac {1}{\textit{Re}} {\nabla} ^2 \boldsymbol{u}, \end{align}

(2.1b) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} & = 0, \\[9pt] \nonumber \end{align}

where $\boldsymbol{u}$ is the velocity vector, $\boldsymbol{u}(\boldsymbol{x},t)=(u,v,w)^T$ , $P(\boldsymbol{x},t)$ is the reduced pressure field and $\boldsymbol{x}=(x,y,z)^T$ indicates the spatial coordinates, while $t$ denotes time.

Equations (2.1) are made dimensionless using the height of the channel before the expansion ( $d$ ) and the maximum (centreline) velocity ( $U_c$ ) of the incoming flow as reference scales for the length and for the velocity, respectively. The Reynolds number ( $\textit{Re}$ ) is thus defined as

(2.2) \begin{equation} \textit{Re}=\frac {U_c \ d}{\nu }, \end{equation}

with $\nu$ the fluid kinematic viscosity.

The boundary conditions complete the system of differential equations, (2.1): at the inlet ( $\partial \mathcal{D}_{\textit{in}}$ ), a fully developed Poiseuille velocity profile is prescribed, no-slip boundary conditions are imposed along the solid walls of the channel ( $\partial \mathcal{D}_w$ ), while a traction-free condition is enforced at the outlet ( $\partial \mathcal{D}_{\textit{out}}$ ) of the domain

(2.3) \begin{align} \left \{ \begin{aligned} \boldsymbol{u} &= U_c \left (1- \left ( \frac {2y}{d}\right )^2 \right ) \boldsymbol{e}_x, & \quad \forall \boldsymbol{x} &\in \partial \mathcal{D}_{\textit{in}} \quad \text{(inlet)}, \\[10pt] \boldsymbol{u} &= \boldsymbol{0}, & \quad \forall \boldsymbol{x} &\in \partial \mathcal{D}_{w} \quad \text{(wall)}, \\[10pt] P \boldsymbol{n} &- \textit{Re}^{-1} (\boldsymbol{\nabla }\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{n}) = \boldsymbol{0}, & \quad \forall \boldsymbol{x} &\in \partial \mathcal{D}_{\textit{out}} \quad \text{(outlet)}, \end{aligned} \right . \end{align}

where $\boldsymbol{e}_x$ is the unit vector in the direction of the $x$ -axis, while $\boldsymbol{n}$ denotes the surface-normal vector (in this case, it is the vector perpendicular to the outlet of the computational domain).

3.1. Methodology

3.1.1. Base flow and direct eigenvalue problem

The onset of the secondary instability is investigated by a classical linear stability analysis (LSA). Thus, the total flow field $ \boldsymbol{Q} = (\boldsymbol{u}, P)^T$ is assumed as the sum of a steady, two-dimensional base flow (or basic state) $\boldsymbol{Q}_b$ and a time-dependent, three-dimensional perturbation $\boldsymbol{q'}$ , i.e. a small deviation from the base flow of amplitude $\epsilon$ , so that

(3.1) \begin{equation} \boldsymbol{Q} (\boldsymbol{x},t) = \begin{pmatrix} \boldsymbol{u} \\ P \end{pmatrix} (\boldsymbol{x},t) = \boldsymbol{Q_b}(\boldsymbol{x}) + \epsilon \boldsymbol{q'}(\boldsymbol{x},t), \end{equation}

where

(3.2a) \begin{align} \boldsymbol{Q}_b(x,y)=(\boldsymbol{u}_b,P_b)^T=(u_b,v_b,0,P_b)^T, \\[-29pt] \nonumber \end{align}

(3.2b) \begin{align} \boldsymbol{q}'(x,y,z,t)=(\boldsymbol{u}',P')^T=(u',v',w',P')^T. \\[6pt] \nonumber \end{align}

Substituting the decomposition (3.1) into the system of differential equations, (2.1), and neglecting quadratic terms, two systems can be obtained that describe the spatial structure of the base flow and the evolution of perturbations.

The base flow is solution of the steady version of (2.1)

(3.3) \begin{equation} \left \{ \begin{aligned} &\left ( \boldsymbol{u}_b \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol{u}_b = -\boldsymbol{\nabla }P + \frac {1}{\textit{Re}_{\textit{BF}}} {\nabla} ^2 \boldsymbol{u}_b, \\[10pt] &\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_b = 0, \end{aligned} \right . \end{equation}

supplemented with the boundary conditions listed in (2.3), while the perturbed field is governed by the linearised, unsteady Navier–Stokes equations (LNSEs)

(3.4) \begin{equation} \left \{ \begin{aligned} &\frac {\partial \boldsymbol{u}'}{\partial t} + \boldsymbol{L} \{\boldsymbol{u}_b(\textit{Re}_{\textit{BF}}), \textit{Re}_{\textit{STB}}\} \ \boldsymbol{u}'+ \boldsymbol{\nabla }P' = \boldsymbol{0}, \\[10pt] &\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}' = 0, \end{aligned} \right . \end{equation}

where $\boldsymbol{L}$ stands for the linearised Navier–Stokes operator that in vector notation can be written as follows:

(3.5) \begin{equation} \boldsymbol{L} \{\boldsymbol{u}_b (\textit{Re}_{\textit{BF}}),\textit{Re}_{\textit{STB}} \} \ \boldsymbol{u}'= \boldsymbol{u}_b \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}'+\boldsymbol{u}'\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}_b-\frac {1}{\textit{Re}_{\textit{STB}}} {\nabla} ^2 \boldsymbol{u}'. \end{equation}

In the above systems (3.3–3.4), the Reynolds number (defined as in 2.2) used for base-flow computations (3.3) is denoted as $\textit{Re}_{\textit{BF}}$ , while the Reynolds number used in the LNSEs (3.4) is labelled as $\textit{Re}_{\textit{STB}}$ . This notation is intended to decouple the Reynolds number that serves as the control parameter for determining the steady base flow, from the Reynolds number that governs the evolution of perturbations (as required by the inviscid structural sensitivity framework discussed in § 5). Henceforth, when $\textit{Re}_{\textit{BF}}=\textit{Re}_{\textit{STB}}$ we simply use $\textit{Re}$ .

Exploiting the homogeneity of the base flow in the spanwise direction ( $z$ ), the perturbation ( $\boldsymbol{q'}$ ) can be expressed using the normal mode decomposition as

(3.6) \begin{equation} \boldsymbol{q'}(x,y,z,t)=\boldsymbol{\hat {q}}(x,y) e^{\gamma t + i k z} + \textrm{c.c.}, \end{equation}

where the term $\textrm{c.c.}$ stands for complex conjugate, while $\boldsymbol{\hat {q}}=(\boldsymbol{\hat {u}},\hat {P})^T=(\hat {u},\hat {v},\hat {w},\hat {P})^T$ denotes the global (direct) eigenfunction and $\gamma \ =\ \sigma +i \omega$ is the complex eigenvalue, that comprises the real growth/decay rate of the perturbation ( $\sigma$ ) and its oscillation circular frequency ( $\omega$ ). Together, eigenvectors and eigenvalues are called the eigenmodes or global modes.

Injecting the modal ansatz (3.6) into the LNSEs (3.4), the following two-dimensional (so-called biglobal) generalised (direct) eigenvalue problem (EVP) is obtained:

(3.7) \begin{equation} \textit{Direct EVP:} \ \left \{ \begin{aligned} &\gamma \ \boldsymbol{\hat {u}} + \boldsymbol{L}\{\boldsymbol{u}_b (\textit{Re}_{\textit{BF}}),\textit{Re}_{\textit{STB}} \} \ \hat {\boldsymbol{u}} +\boldsymbol{\nabla }\hat {P}= \boldsymbol{0}, \\[10pt] &\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat {u}} = 0. \end{aligned} \right . \end{equation}

For a given $k$ , the complex $\gamma$ can be computed by solving the stability problem (3.7), complemented with the following boundary conditions:

(3.8) \begin{equation} \left \{ \begin{aligned} &\boldsymbol{\hat {u}} = \boldsymbol{0}, \quad \text{on} \ \partial \mathcal{D}_{\textit{in}} \cup \partial \mathcal{D}_{w} \ \text{(inlet and wall)}, \\[10pt] &\hat {P} \boldsymbol{n} - \textit{Re}_{\textit{STB}}^{-1} (\boldsymbol{\nabla }\boldsymbol{\hat {u}} \boldsymbol{\cdot }\boldsymbol{n}) = \boldsymbol{0}, \quad \text{on} \ \partial \mathcal{D}_{\textit{out}} \ \text{(outlet)}. \end{aligned} \right . \end{equation}

If all the eigenvalues ( $\gamma$ ) are such that $\sigma$ is negative, then the base flow is stable with respect to small three-dimensional perturbations (i.e. the disturbance is damped). If, on the other hand, there exists at least one eigenvalue with positive real part, then the base flow is unstable (i.e. the disturbance is amplified), namely the perturbations ( $\boldsymbol{q'}$ ) grow with time and the flow will evolve away from its initial state ( $\boldsymbol{Q}_b$ ).

3.1.2. Adjoint eigenvalue problem

To determine the instability core by means of the adjoint-based structural sensitivity analysis (in § 4), the study of both the direct and adjoint modes is required. In particular, the adjoint global mode represents the receptivity of the direct mode to external forcing.

Following the framework introduced by Giannetti & Luchini (Reference Giannetti and Luchini2007), the adjoint perturbation $\boldsymbol{\hat {q}}^\dagger =(\boldsymbol{\hat {u}}^\dagger ,\hat {P}^\dagger )^T$ satisfies the following equations:

(3.9) \begin{equation} \textit{Adjoint EVP:} \ \left \{ \begin{aligned} &-\gamma \boldsymbol{\hat {u}^\dagger } + \boldsymbol{L}^{\dagger } \{\boldsymbol{u}_b (\textit{Re}_{\textit{BF}}),\textit{Re}_{\textit{STB}} \} \hat {\boldsymbol{u}}^{\dagger } +\boldsymbol{\nabla }\hat {P}^\dagger = \boldsymbol{0 }, \\[10pt] &\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat {u}}^\dagger = 0. \end{aligned} \right . \end{equation}

with corresponding boundary conditions and $\boldsymbol{L}^{\dagger }$ the adjoint linearised Navier–Stokes operator expressed as

(3.10) \begin{equation} \boldsymbol{L}^{\dagger } \{\boldsymbol{u}_b (\textit{Re}_{\textit{BF}}), \textit{Re}_{\textit{STB}} \} \boldsymbol{u}'^{\dagger }=\boldsymbol{u}_b \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}'^{\dagger }-\boldsymbol{\nabla }\boldsymbol{u}_b \boldsymbol{\cdot }\boldsymbol{u}'^{\dagger }+\frac {1}{\textit{Re}_{\textit{STB}}} {\nabla} ^2 \boldsymbol{u}'^{\dagger }. \end{equation}

More details on the derivation of the adjoint problem can be found in Giannetti & Luchini (Reference Giannetti and Luchini2007), while for a complete review of the role of the adjoint equations in stability analysis the reader is referred to Luchini & Bottaro (Reference Luchini and Bottaro2014).

3.2. Results

As is well known from the literature (Fani et al. Reference Fani, Camarri and Salvetti2012), the flow in a $1 \! : \! 3$ plane sudden expansion remains symmetric for Reynolds numbers below $\textit{Re}_{\textit{cr}}^{(I)}\approx 81$ . However, under supercritical conditions, while the flow remains two-dimensional and steady, it becomes asymmetric.

We perform a global LSA around this secondary flow to investigate the onset of the secondary bifurcation. The numerical results, shown in figure 2, reveal that the flow undergoes a steady ( $\omega =0$ ; see figure 2 b), three-dimensional ( $k \ne 0$ ) bifurcation at a critical Reynolds number of approximately $\textit{Re}_{\textit{cr}}^{(\textit{II})}\approx 306$ , associated with spanwise wavenumber $k\approx 0.65$ (figure 2 a), in agreement with previous numerical investigations (see e.g. Lanzerstorfer & Kuhlmann Reference Lanzerstorfer and Kuhlmann2012a ). The eigenvalue spectrum (at $\textit{Re}_{\textit{cr}}^{(\textit{II})}=306$ and $k=0.65$ ), illustrated in figure 3, shows that the eigenvalues are real or come in complex-conjugate pairs, as expected, and that the least stable eigenvalues are discrete.

Figure 2. (a) Growth rate $\sigma$ and (b) frequency $\omega$ of the least-damped eigenvalue $\gamma$ as a function of the spanwise wavenumber $k$ , for the cases: $\textit{Re}=304$ (dashed line) and the critical value $\textit{Re}=\textit{Re}_{\textit{cr}}^{(\textit{II})}=306$ (solid line). The red marker in (a) highlights the value of the wavenumber ( $k=0.65$ ) corresponding to the maximum amplification rate for the case $\textit{Re}=306$ .

Figure 3. Eigenvalue spectrum at $\textit{Re}=\textit{Re}_c^{(\textit{II})}=306$ and $k=0.65$ . The red marker, labelled with the letter A, indicates the least stable eigenvalue, with coordinates $(7\times 10^{-5},0)$ in the plane $(\sigma ,\omega )$ .

The structure of the base flow, computed at the critical Reynolds number, is depicted in figure 4. The flow separates from the lower wall, impinges on the upper edge of the domain and then reattaches to the bottom wall. This phenomenon is usually explained by a Coanda effect, that is a consequence of the pressure difference established in the cross-stream direction (see the early work by Bourque & Newman Reference Bourque and Newman1960). The resulting asymmetric flow field can have one of two specular orientations, corresponding to the main stream being diverted towards one wall of the channel or the other. In particular, there is an equal probability of the jet bending towards the upper or the lower side of the domain. The streamline plot, of one of the two possible asymmetric flows, shows that the flow field is characterised by three recirculation bubbles. Two of them separate immediately past the sharp corners, with the larger vortex referred to as the primary vortex, while the smaller one, on the opposite side of the channel, is denoted as the secondary vortex. Further downstream, on the same side of the secondary vortex, a third separated region arises (tertiary vortex), while at a streamwise distance of approximately $x=35$ , the flow recovers its parabolic profile.

Figure 4. Stationary, two-dimensional and asymmetric base flow computed at the critical Reynolds number ( $\textit{Re}=\textit{Re}_{\textit{cr}}^{(\textit{II})}=306$ ). The contour plot displays the modulus of the velocity field, while solid lines indicate the streamlines inside the domain.

The features of the bifurcation can be described by investigating the spatial structure of the direct and adjoint global eigenfunctions, computed at the critical conditions (figure 5). The absolute value of the fluctuating velocity field and its adjoint are depicted in figure 5(a,e). It can be observed that the direct mode (figure 5 a) is strictly confined within the primary recirculation region and around the tertiary vortex, while for the adjoint mode (figure 5 e) the main contribution of the perturbation is localised close to the primary vortex core, on the upper wall of the domain, and it decays moving toward the outlet of the channel.

Figure 5. Spatial structure of the direct and adjoint eigenfunctions at the critical Reynolds number ( $\textit{Re}=\textit{Re}_{\textit{cr}}^{(\textit{II})}=306$ ) associated with a spanwise wavenumber $k=0.65$ : (a) modulus, (b) streamwise $\hat {u}$ , (c) cross-stream $\hat {v}$ and (d) spanwise $\hat {w}$ velocity (direct) perturbed components. Similarly, (e–h) depict the corresponding adjoint eigenfunctions.

With reference to the direct eigenfunctions (figure 5 b–d), the $u$ -perturbation is the strongest, with the cross-stream ( $\hat {v}$ ) and spanwise ones ( $\hat {w}$ ) being respectively approximately $29.7\,\%$ and $17.3\,\%$ of the streamwise fluctuations. On the other hand, the adjoint field (figure 5 f–h) reveals that the maximum receptivity is reached near the upper edge of the domain and the downstream side of the primary vortex for the streamwise velocity component ( $\hat {u}^\dagger$ ), whereas in the cross-stream ( $\hat {v}^\dagger$ ) and spanwise ( $\hat {w}^\dagger$ ) directions, the fluctuations are most evident inside the primary recirculation region and the bottom wall, respectively (in particular, $\hat {v}^\dagger \approx 0.524 \ \hat {u}^\dagger$ and $\hat {w}^\dagger \approx 0.241 \ \hat {u}^\dagger$ ).

4.1. Methodology

To investigate the physical mechanism responsible for the instability, we compute the structural sensitivity, which allows us to localise the core of the global instability, i.e. the ‘active’ flow regions where the amplification of the perturbations and receptivity combine to trigger the instability. The fundamental idea of such approach was introduced by Giannetti & Luchini (Reference Giannetti and Luchini2007) in an attempt to extend to strongly non-parallel flows the concept of ‘wavemaker’, first introduced in the asymptotic analysis of slowly evolving flows by Huerre & Monkewitz (Reference Huerre and Monkewitz1990). Thus, following Giannetti & Luchini (Reference Giannetti and Luchini2007), the structural sensitivity tensor can be evaluated by combining the direct ( $\boldsymbol{\hat {u}}$ ) and adjoint ( $\boldsymbol{\hat {u}^\dagger }$ ) eigenvectors

(4.1) \begin{equation} \boldsymbol{S}(x_0,y_0;\textit{Re})=\frac {\boldsymbol{\hat {u}}^\dagger \ \boldsymbol{\hat {u}}}{\int _{\mathcal{D}} \boldsymbol{\hat {u}^\dagger } \boldsymbol{\cdot }\boldsymbol{\hat {u}} \ dV }. \end{equation}

The tensor $\boldsymbol{S}$ gives the shift in the eigenvalue produced by a velocity feedback localised at $(x_0,y_0)$ . Here, $\boldsymbol{\hat {u}}^\dagger \ \boldsymbol{\hat {u}}$ stands for the dyadic product of the two vectors.

A spatial sensitivity map is then traced by computing the norm of the tensor at each point of the domain (see Giannetti & Luchini Reference Giannetti and Luchini2007 for details). Large values of $\boldsymbol{S}$ identify the regions where the localised feedback produces the largest drift in the growth rate or in the frequency of the mode. As in Giannetti & Luchini (Reference Giannetti and Luchini2007), here, we use the Frobenius norm of this sensitivity tensor to identify where a modification in the linearised equations produces the greatest drift of the eigenvalue and thereby reveal the region of the flow that acts as the wavemaker.

4.2. Results

The spatial distribution of the sensitivity tensor, computed at the critical Reynolds number, is depicted in figure 6. It can be observed that the primary recirculation bubble is the most sensitive region of the flow, as the structural sensitivity field is highly localised therein and its peak spreads near the primary vortex centre. Additional information could be gained by inspecting the absolute value of the individual components of the sensitivity tensor (as proposed by Qadri et al. Reference Qadri, Mistry and Juniper2013, Reference Qadri, Chandler and Juniper2015). This analysis allows us to distinguish the direct and adjoint velocity components that take a relevant role in the feedback process. All the components depicted in figure 7 show spatial support within the primary recirculation bubble. It is evident that the spatial structure of the sensitivity tensor is essentially due to the three components related to the streamwise velocity of the direct global mode (figure 7 a,b,c). Thus, we can state that the feedback process consists of a loop in which the streamwise perturbations couple with the local receptivity of the flow in all the spatial directions. Although this analysis allows us to visualise the flow regions where feedback between the components of the velocity is strong, one can hardly assess the character of the instability simply based on the inspection of these maps. In particular, it is not clear whether the wavemaker is linked to the vortex centre (see e.g. figure 7 a,b) or whether the core of the instability is distributed over a closed streamline (see e.g. figure 7 c), that plays a crucial role in the feedback process. The question of the nature of the amplification process, operating within the primary recirculation zone, is addressed in the next section (§ 5).

Figure 6. Structural sensitivity map, showing the core of the instability at the critical conditions ( $\textit{Re}\ =\textit{Re}_{\textit{cr}}^{(\textit{II})}=\ 306, \ k=0.65$ ).

Figure 7. The absolute values of the components of the sensitivity tensor $\boldsymbol{S}$ at the critical Reynolds number ( $\textit{Re}=\textit{Re}_{\textit{cr}}^{(\textit{II})}=306$ ) associated with a spanwise wavenumber $k=0.65$ . The shading on all the plots is scaled from $0$ (white) to $2$ (red). Each frame spans $0\lt x\lt 20$ on the horizontal axis and $-1.5\lt y\lt 1.5$ on the vertical axis.

5.1. Methodology

The sensitivity analysis performed in the previous section (§ 4) clearly shows that the core of the three-dimensional instability, leading to the secondary bifurcation, is confined within the primary recirculation bubble. However, this observation is insufficient to unveil the underlying physical process, even if the relatively compact spatial support of the structural sensitivity field, whose peak spreads near the primary vortex centre, and the presence of locally closed elliptical streamlines in the core vortex, suggest that the instability could be associated with an inviscid mechanism.

Sipp et al. (Reference Sipp, Lauga and Jacquin1999) showed that different inviscid instabilities are characterised by spatial supports in different regions of the flow domain. Specifically, centrifugal and elliptic eigenmodes are concentrated along closed streamlines and in the vortex centres, respectively. However, a simple inspection of the sensitivity map (see figure 6) does not clarify whether the wavemaker is linked to the vortex centre or whether the core of the instability is distributed over a closed streamline, that plays a crucial role in the feedback process.

To shed light on the physical nature of the secondary bifurcation, and to assess whether it is of centrifugal or elliptic type, we aim at isolating the effect of the inviscid mechanism through a sensitivity-based protocol. The fundamental idea is to decouple the system of equations that govern the base flow from those that describe the disturbance evolution (as anticipated in § 2 with the notation $\textit{Re}_{\textit{STB}}$ and $\textit{Re}_{\textit{BF}}$ ) to change the viscosity effects solely in the global stability equations, while maintaining the base-flow field constant.

From a practical standpoint, we freeze the base flow at the critical Reynolds number ( $\textit{Re}_{\textit{BF}}=\textit{Re}_{\textit{cr}}^{(\textit{II})}=306$ ) and then repeat the structural sensitivity analysis while progressively increasing $\textit{Re}_{\textit{STB}}$ in the linearised equations. The resulting inviscid structural sensitivity tensor is thereby a function of both the base flow ( $\textit{Re}_{\textit{BF}}$ ) and the stability ( $\textit{Re}_{\textit{STB}}$ ) Reynolds number

(5.1) \begin{equation} \boldsymbol{S}^*(\textit{Re}_{\textit{BF}}, \textit{Re}_{\textit{STB}})=\frac {\boldsymbol{\hat {u}}^\dagger (\boldsymbol{u}_b(\textit{Re}_{\textit{BF}});\textit{Re}_{\textit{STB}}) \ \boldsymbol{\hat {u}} (\boldsymbol{u}_b(\textit{Re}_{\textit{BF}});\textit{Re}_{\textit{STB}})}{\int _{\mathcal{D}} \boldsymbol{\hat {u}^\dagger } \boldsymbol{\cdot }\boldsymbol{\hat {u}} \ dV }. \end{equation}

The norm of $\boldsymbol{S}^*$ can then be used to trace the inviscid structural sensitivity map. If the latter concentrates in a specific spatial region as $\textit{Re}_{\textit{STB}}$ goes to infinity, then one can conjecture that an inviscid instability mechanism is at play and determine the flow regions where it acts.

5.2. Results

The results given by the inviscid structural sensitivity analysis are plotted in figure 8(a). The graph shows the evolution of the growth rate ( $\sigma$ ) associated with the unstable global mode under consideration as a function of the spanwise wavenumber ( $k$ ) for several stability Reynolds numbers. In particular, numerical experiments were conducted at $\textit{Re}_{\textit{BF}}=306$ and $\textit{Re}_{\textit{STB}}= 612, 918, 1530, 3060,6120$ . It can be observed that, as $\textit{Re}_{\textit{STB}}$ increases, the curves reach their peaks (highlighted by the red markers in the figure) at progressively higher wavenumbers (this aspect is discussed in detail in the following section § 6). By observing the spatial structure of the inviscid sensitivity tensor (figure 8 b–d), it is evident that, as $\textit{Re}_{\textit{STB}}$ becomes larger, the wavemaker tends to focus in the primary vortex centre. Furthermore, the spatial support of these fields is concentrated around their maximum values. These results preclude a centrifugal instability from being at play and suggest that the amplification process could be associated with an elliptic mechanism. In particular, by inspecting the spatial distribution of the structural sensitivity field, computed at the maximum $\textit{Re}_{\textit{STB}}$ considered (equal to $6120$ ) (figure 8 d), it is evident that the critical mode is confined within the primary eddy and is strongest in its centre, which is a hallmark of an elliptic instability (see Bayly Reference Bayly1986; Pierrehumbert Reference Pierrehumbert1986; Waleffe Reference Waleffe1990; Kerswell Reference Kerswell2002, for a comprehensive description of the elliptic instability mechanism).

Figure 8. Inviscid structural sensitivity results at $\textit{Re}_{\textit{BF}}=\textit{Re}_{\textit{cr}}^{(\textit{II})}=306$ , while progressively increasing $\textit{Re}_{\textit{STB}}$ in the LNSEs only. (a) Growth rate $\sigma$ vs spanwise wavenumber $k$ . The red markers denote the maximum of the curves. The last two points are labelled with the letters $A$ and $B$ as they will be referenced in the subsequent section (§ 6). (b–d) Inviscid structural sensitivity maps at (b) $\textit{Re}_{\textit{STB}}=612$ , (c) $\textit{Re}_{\textit{STB}}=1530$ and (d) $\textit{Re}_{\textit{STB}}=6120$ .

In summary, our findings indicate that the classical structural sensitivity ( $\boldsymbol{S}$ ) effectively identifies the wavemaker within the primary recirculation bubble (see figure 6 of § 4). However, it is only the inviscid structural sensitivity field ( $\boldsymbol{S}^*$ ) that reveals that the instability core is concentrated around the centre of an elliptical vortex (see figure 8 d).

Therefore, our results complete the picture proposed by Lanzerstorfer & Kuhlmann (Reference Lanzerstorfer and Kuhlmann2012a ). Using an a posteriori energy transfer method, they observed that the flow dynamics changes with the expansion ratio. For large expansion ratios (i.e. $3.484 \leqslant \textit{ER} \leqslant 10$ ; note that this range has been redefined according to our definition of $\textit{ER}$ ) the critical mode passes from stationary to oscillatory and pure centrifugal and elliptical amplification processes are identified. However, for moderate expansion ratios ( $1.538 \leqslant \textit{ER} \leqslant 3.334$ ), which is the case of the present study where a stationary three-dimensional bifurcation is detected, the physical interpretation becomes challenging and the energy approach fails to predict the inviscid character of the instability. The authors argue that, in the moderate expansion ratios range, the flow becomes unstable due to streamline convergence within the downstream side of the primary vortex, in combination with flow deceleration near the reattachment point, and an amplification process due to shear stresses near and between the primary and tertiary vortex. In contrast, by using inviscid structural sensitivity analysis, we demonstrate that an elliptic instability is already present in the steady regime.

In this context, the sensitivity protocol acts as a magnifying glass, allowing us to clearly isolate the inviscid core of the instability, thereby completing the characterisation of the secondary instability also for moderate ERs. Furthermore, the sensitivity protocol presented in this work is a general concept that can be used to discern any inviscid-type instability. Thus, this method can effectively be used also in other flow configurations (see for instance the works of Citro et al. Reference Citro, Giannetti, Brandt and Luchini2015; Chiarini & Auteri Reference Chiarini and Auteri2023), to unveil the inviscid character of an instability, thereby enhancing the insight into the mechanism whereby complex three-dimensional motion can arise directly from two-dimensional coherent structures.

5.3. Results for a large expansion ratio

The geometric configuration examined in this study is known to exhibit different instability mechanisms depending on the value of the expansion ratio (Lanzerstorfer & Kuhlmann Reference Lanzerstorfer and Kuhlmann2012a ). Here, we demonstrate the effectiveness of the sensitivity-based protocol in detecting inviscid instabilities in a different configuration, specifically one characterised by $\textit{ER}=20$ . As in Lanzerstorfer & Kuhlmann (Reference Lanzerstorfer and Kuhlmann2012a ), we found that the flow undergoes an unsteady three-dimensional bifurcation at a critical Reynolds number of approximately $\textit{Re}_{\textit{cr}}^{(\textit{II})}\approx 42.9$ , associated with a spanwise wavenumber of $k\approx 0.13$ and an oscillation circular frequency of $\omega \approx 0.048$ .

To assess the physical mechanism responsible for the instability, we compute both the classical and the inviscid structural sensitivity fields. The results are depicted in figure 9. It can be observed that, at the critical Reynolds number, the structural sensitivity field appears to be spatially diffused (figure 9 a). However, by inspecting the inviscid sensitivity field computed at $\textit{Re}_{\textit{STB}}=100$ and $\textit{Re}_{\textit{STB}}=3000$ (figure 9 b–c), it is evident that as $\textit{Re}_{\textit{STB}}$ increases, the wavemaker tends to focus around a closed orbit, as expected for a centrifugal instability. Thus, in accordance with the energetic approach of Lanzerstorfer & Kuhlmann (Reference Lanzerstorfer and Kuhlmann2012a ), it is possible to attribute the nature of the instability to a centrifugal mechanism due to the localisation of the inviscid sensitivity field around a closed streamline that falls in a region in which the Rayleigh (Reference Rayleigh1917) criterion is valid (see figure 12 of Lanzerstorfer & Kuhlmann Reference Lanzerstorfer and Kuhlmann2012a ). In conclusion, we highlight that, not only is the inviscid structural sensitivity able to identify the centrifugal character of the instability associated with a closed streamline in the recirculating region, but it also accurately predicts the particle orbit that provides the main contribution to the instability.

Figure 9. Inviscid structural sensitivity results for the case $\textit{ER}=20$ , at $\textit{Re}_{\textit{BF}}=\textit{Re}_{\textit{cr}}^{(\textit{II})}=\ 42.9$ , while progressively increasing $\textit{Re}_{\textit{STB}}$ in the LNSEs only. Panels show (a) $\textit{Re}_{\textit{BF}}=\textit{Re}_{\textit{STB}}= \ 42.9$ , (b) $\textit{Re}_{\textit{BF}}= \ 42.9$ and $\textit{Re}_{\textit{STB}}=\ 100$ , (c) $\textit{Re}_{\textit{BF}}= \ 42.9$ and $\textit{Re}_{\textit{STB}}= \ 3000$ , plotted at the optimal spanwise wavenumber.

6.1. Methodology

To confirm the evidence, given by the inviscid structural sensitivity analysis (in § 5) of an elliptical instability mechanism, we pursue the investigation using a local theory based on the geometric optics (or short-wave asymptotic) approach. This theory has been successfully applied in the past to study elliptic, hyperbolic and centrifugal instabilities developing on two-dimensional stationary base flows. It is based on the fact that short-wave instabilities can be described as localised wave envelops that move along trajectories of fluid elements. The motion is governed by a system of characteristic equations along the corresponding trajectories, consisting of the eikonal equation for the wavenumber and the transport equation for the velocity amplitude. Thus, the burden of proof that a given flow is unstable is placed on the analysis of ordinary differential equations along trajectories of fluid elements rather than of the three-dimensional Euler equations for perturbations (Lifschitz Reference Lifschitz1994). In the following, this approach is briefly presented; for a detailed description, the reader is referred to the work of Lifschitz & Hameiri (Reference Lifschitz and Hameiri1991) in which the whole theory is thoroughly explained and applied.

Within the classical WKB approximation, the solution of the LNSEs (3.4) is sought in the form of a localised rapidly oscillating wave packet evolving along a closed Lagrangian trajectory, $\boldsymbol{X}(t)$ such that

(6.1a) \begin{align} \boldsymbol{u'}(\boldsymbol{X},t)= e^{i \frac {\phi (\boldsymbol{X},t)}{\epsilon }} \ \boldsymbol{a}(\boldsymbol{X},t,\epsilon )=e^{i \frac {\phi (\boldsymbol{X},t)}{\epsilon }} \ \sum _n \boldsymbol{a}_n(\boldsymbol{X},t) \epsilon ^n, \\[-28pt] \nonumber \end{align}

(6.1b) \begin{align} p'(\boldsymbol{X},t)= e^{i \frac {\phi (\boldsymbol{X},t)}{\epsilon }} \ b(\boldsymbol{X},t,\epsilon )=e^{i \frac {\phi (\boldsymbol{X},t)}{\epsilon }} \ \sum _n b_n(\boldsymbol{X},t) \epsilon ^{n+1}, \\[6pt] \nonumber \end{align}

where $\phi$ denotes the eikonal function, whose spatial derivatives represent the local wavenumber components $\boldsymbol{k}=\boldsymbol{\nabla }\phi (\boldsymbol{X},t)$ , while $\epsilon \ll 1$ is a small parameter expressing the smallness of the scale of the waves represented by (6.1).

Substituting the expansions (6.1a )–(6.1b ) into the perturbation equations (3.4) and grouping terms multiplied by the same power of $\epsilon$ , a hierarchy of equations is obtained. The leading-order approximation yields a set of ordinary differential equations (ODEs) that, in the inviscid $(\textit{Re} \rightarrow \infty )$ and short-wave $(||\boldsymbol{k}|| \rightarrow \infty )$ limits, provide the growth rate associated with a localised perturbation. Specifically, the system of equations describes the evolution of the wavenumber vector $\boldsymbol{k}$ and the amplitude $\boldsymbol{a}$ as follows:

(6.2) \begin{equation} \left \{ \begin{aligned} \frac {D \boldsymbol{k}}{\textrm{D}t} &= - \mathscr{L}^T (\boldsymbol{X}) \boldsymbol{k}, \\ \frac {D \boldsymbol{a}}{\textrm{D}t} &= \left ( \frac {2 \boldsymbol{k} \boldsymbol{k}^T}{|\boldsymbol{k}|^2} - \mathscr{I} \right ) \mathscr{L}(\boldsymbol{X}) \boldsymbol{a}, \end{aligned} \right . \end{equation}

where $\boldsymbol{X}$ is the position of the fluid particle moving along the streamlines of the steady base flow

(6.3) \begin{equation} \frac {\textrm{D}\boldsymbol{X}(t)}{\textrm{D}t}=\boldsymbol{u}_b(\boldsymbol{X}(t),t). \end{equation}

In the above equations, $D/\textrm{D}t$ represents the so-called material derivative (i.e. $D/\textrm{D}t=\partial /\partial t + \boldsymbol{u}_b \boldsymbol{\cdot }\boldsymbol{\nabla}$ ), $\mathscr{L}=\boldsymbol{\nabla }\boldsymbol{u}_b$ denotes the base-flow velocity-gradient tensor, $\mathscr{I}$ is the identity matrix and the superscript $T$ denotes the transpose. Lifschitz & Hameiri (Reference Lifschitz and Hameiri1991) proved that the flow is unstable if there exists a streamline along which the amplitude $\boldsymbol{a}(t)$ grows without bound as $t\rightarrow \infty$ . Since the inviscid structural sensitivity analysis (§ 5) suggested the occurrence of an elliptic amplification process, we evaluated the inviscid local stability equations (6.2) at the vortex centre of the primary recirculation bubble (thereby, following the work of Bayly Reference Bayly1986, the tensor $\mathscr{L}$ has been computed at this point).

Provided $\boldsymbol{k}(t)$ is periodic in time, the second equation of (6.2) is a Floquet problem for $\boldsymbol{a}(t)$ . As is well known (Bender & Orszag Reference Bender and Orszag2013), the general solution is a linear superposition of Floquet modes

(6.4) \begin{equation} \boldsymbol{a}(t)=\boldsymbol{\bar {a}}(t)e^{\gamma t}, \end{equation}

where $\boldsymbol{\bar {a}}(t)$ is a periodic function with the same period ( $T_p$ ) as $\boldsymbol{k}$ , while $\gamma$ is the corresponding complex Floquet exponent, to be determined by computing the monodromy matrix $\mathscr{A}(T_p)$ that satisfies

(6.5) \begin{equation} \frac {\textrm{D}\mathscr{A}}{Dt}=\left ( \frac {2 \boldsymbol{k} \boldsymbol{k}^T}{|\boldsymbol{k}^2|}-\mathscr{I}\right ) \mathscr{L}(\boldsymbol{X}) \mathscr{A} \ \ \ \text{with} \ \ \mathscr{A}(0)=\mathscr{I}. \end{equation}

Thus, the instability problem reduces to the calculation of the matrix $\mathscr{A}$ and its non-trivial eigenvalues ( $\gamma _i$ with $i=1,2,3$ ). Assuming that the real parts ( $\mu _i$ ) of such eigenvalues satisfy the condition $\mu _1\gt \mu _2\gt \mu _3$ , the asymptotic growth rate is given by the simple relation

(6.6) \begin{equation} \sigma _{\infty } = \frac {\log (\mu _1)}{T_p}, \end{equation}

where $T_p$ denotes the period of revolution of $\boldsymbol{k}$

To validate our numerical implementation of the asymptotic approach, we replicated Bayly’s analytical results (Bayly Reference Bayly1986) (as detailed in Appendix B) and subsequently used this numerical tool to compute the amplification rate ( $\sigma _{\infty }$ ) for the problem under investigation.

6.2. Results

The asymptotic amplification rate ( $\sigma _{\infty }$ ), computed in the centre of the primary vortex, is displayed in the fifth column of table 1. To make a quantitative comparison between the results obtained from the local and global stability analysis, we followed the work of Sipp et al. (Reference Sipp, Lauga and Jacquin1999). Therefore, considering a relation of the type

(6.7) \begin{equation} \sigma (k)=\sigma _0-\frac {\xi }{k}, \end{equation}

we computed an extrapolated value of the (global) growth rate in the limit of large wavenumbers, namely $\sigma (k\rightarrow \infty )=\sigma _0$ . In particular, the parameters of (6.7), i.e. $\sigma _0$ and $\xi$ (the latter called the eigenvalue convergence parameter), have been evaluated using the amplification rates ( $\sigma _A$ , $\sigma _B$ ) and the corresponding wavenumbers ( $k_A$ , $k_B$ ) given by the global LSA, for the last two cases investigated using the inviscid sensitivity protocol (indicated in figure 8 a of § 5 with A and B, while in table 1 they are denoted $\sigma _A$ and $\sigma _B$ ).

Table 1. Comparison between the amplification rates given by the LSA and that given by the short-wave asymptotics (WKB) in the centre of the primary recirculation bubble ( $\sigma _{\infty }$ ). The last column estimates the discrepancy between the two methods. $\xi$ and $\sigma _0$ are extrapolated from the results shown in the plane $(k,\sigma )$ (see figure 8a in § 5) using the points $A=(k_A,\sigma _A)=(1.77,0.05166)$ and $B=(k_B,\sigma _B)=(2.08,0.05922)$ and assuming a relation of the form $\sigma =\sigma _0-\xi /k$ (from Sipp et al. Reference Sipp, Lauga and Jacquin1999).

As reported in the sixth column of table 1, the discrepancy in percentage between the predictions provided by the two theories (global, $\sigma _0$ and local, $\sigma _{\infty }$ ) is approximately $10\,\%$ .

Figure 10 a shows the comparison between the global growth rates and those predicted through (6.7). A good agreement is observed between the two sets of data. In particular, the global stability analysis is more computationally expensive with respect to the geometric optics approach, since the computation of a global eigenpair involves the solution of a two-dimensional generalised EVP, while the WKB approach involves the simple integration of a tiny set of ODEs in one dimension.

Figure 10. (a) Comparison between the maximum amplification rates $\sigma$ (i.e. the y coordinate of the red markers of figure 8 a) given by the LSA (lines) and those predicted by the WKB method using the relation proposed by Sipp et al. (Reference Sipp, Lauga and Jacquin1999). (b) Spanwise wavenumber $k$ at maximum growth rate (i.e. the abscissa of the red markers of figure 8 a) as a function of $\textit{Re}_{\textit{STB}}$ . A logarithmic scale is used on the x-axis only for panel (b).

Figure 10 b shows how the global least stable spanwise wavenumber $k$ grows with the Reynolds number ( $\textit{Re}_{\textit{STB}}$ ). In particular, it increases almost linearly with the logarithm of $\textit{Re}_{\textit{STB}}$ . This is in agreement with the hypothesis of the asymptotic theory, suggesting that $k$ will tend to infinity in the limit of vanishing viscosity.