2.1 Open-loop control and synchronization studies

The PVC is a self-excited oscillator. The oscillatory motion is nearly harmonic with a characteristic natural frequency. When the system is harmonically perturbed or forced by open-loop control, a phenomenon called frequency synchronization may occur where the system stops oscillating at its natural frequency and starts to oscillate at the forcing frequency instead. For this, the forcing frequency needs to be sufficiently close to the natural frequency of the system, and the forcing amplitude needs to be of sufficient magnitude. In this work, only $1:1$ forced synchronizations are of interest. This means that the PVC and the forcing are unidirectionally coupled oscillators: the forcing influences the PVC, but not vice versa; and the synchronized frequency of the PVC is equal to the forcing frequency instead of a (reciprocal) integer multiple of the forcing frequency (Balanov et al. Reference Balanov, Janson, Postnov and Sosnovtseva2008).

In order to force the PVC with the goal of reaching synchronization, eight ZNMF actuation units are used, capable of exciting azimuthal modes of $m=0$ to $\pm 4$. Each unit consists of one loudspeaker with a rated input power of $15~\text{W}$ connected by small channels to the duct, as sketched in figure 3. Each stroke period of the speaker creates a ZNMF jet at the slot exit of the actuator channel with a rectangular cross-section of $1~\text{mm}\times 22~\text{mm}$. Addressing each speaker with a different phase lag at a specific time allows periodic excitation. The actuation is applied at the azimuthal wavenumber of $m=1$, equal to the azimuthal wavenumber of the PVC. The input actuation signal $A_{k}$ of speaker $k$ is, therefore, given by

(2.6)$$\begin{eqnarray}A_{k}=A_{f}\sin \left(2\unicode[STIX]{x03C0}\left[f_{f}t+\frac{k-1}{8}\right]\right),\end{eqnarray}$$

with $A_{f}$ being the forcing amplitude and $f_{f}$ being the forcing frequency.

All speakers are calibrated at no-flow conditions with a microphone placed at $r=0$ in the actuation plane for the relevant frequency range in order to produce repeatable sound pressure levels. Moreover, preliminary hot-wire measurements at no-flow conditions show that the transfer function of input voltage to output peak velocity is linear within the utilized frequency and amplitude ranges. Nonetheless, the amplitudes in the actuation experiments are still given in terms of voltage since the peak velocities at no-flow conditions are likely to be overestimated compared to flow conditions, thus not being directly proportional. Therefore, forcing and synchronization amplitudes are simply quantified by the input voltage of the speakers.

Five actuation positions are tested for the open-loop control studies. Measured from the exit plane of the duct, these are $x_{a}/D=-2$, $-1.5$, $-1$, $-0.75$ and $-0.5$. Actuation is applied within the range of $\pm 10\,\%$ of the natural frequency. The maximum frequency increments are $\pm 2.5\,\%\,f_{n}$ or lower between the selected forcing frequencies $f_{f}$. One of these five actuation positions is exemplarily illustrated in figure 1. The experimental parameters are summarized in table 1.

Table 1.

Overview of experimental parameters.

A viable criterion as to whether synchronization has been reached can be defined by considering the phase angle difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}$ between the instantaneous phase of the PVC, $\unicode[STIX]{x1D711}_{PVC}$, and the phase of the forcing, $\unicode[STIX]{x1D711}_{f}$:

(2.7)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D711}=\unicode[STIX]{x1D711}_{PVC}-\unicode[STIX]{x1D711}_{f}+\unicode[STIX]{x0394}\unicode[STIX]{x1D711}_{0},\end{eqnarray}$$

with $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}_{0}$ as the initial phase difference at $t=0$, arbitrarily set to $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}_{0}=0$ here. Synchronization is established when

(2.8)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D711}\approx 0,\quad \text{for }t\rightarrow \infty ,\end{eqnarray}$$

i.e. when the PVC and forcing frequency are approximately equal for the time of measurement of 10 s. Exact equality is not possible in this highly turbulent set-up. The forcing amplitude where synchronization is reached is called the ‘synchronization amplitude’.

3.1 Empirical modes by proper orthogonal decomposition

The PVC is a coherent structure and, therefore, typically described with the triple decomposition technique. To extract dominant coherent velocity fluctuations from a turbulent velocity field, POD is a well-tested technique (Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993). In the POD, an ensemble of $N$ velocity fields is projected onto an orthogonal $N$-dimensional basis that maximizes the turbulent kinetic energy for any subspace. Thus, POD provides an optimal set of spatial modes with which the energy-containing structures can be described in the most efficient way with as few modes as possible. The POD of a velocity field $\boldsymbol{u}$ is determined with

(3.2)$$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)=\overline{\boldsymbol{u}}(\boldsymbol{x})+\mathop{\sum }_{j=1}^{N}a_{j}(t)\unicode[STIX]{x1D731}_{j}(\boldsymbol{x})+\boldsymbol{u}_{res}(\boldsymbol{x},t)\end{eqnarray}$$

by minimizing the residual $\boldsymbol{u}_{res}$. The $a_{j}(t)$ are the POD time coefficients, while the $\unicode[STIX]{x1D731}_{j}(\boldsymbol{x})$ are the spatial POD modes. In the case of SPIV, $N$ is the number of snapshots. The POD modes $\unicode[STIX]{x1D731}_{j}$ provide the spatial shape of the mode and the time coefficients $a_{j}$ provide the time-dependent amplitude of the modes. Furthermore, the mean specific turbulent kinetic energy of the associated modes can be calculated via $\unicode[STIX]{x1D706}_{j}=\overline{a_{j}^{2}}$. The POD modes are energy ranked such that $\unicode[STIX]{x1D706}_{1}\geqslant \unicode[STIX]{x1D706}_{2}\geqslant \cdots \geqslant \unicode[STIX]{x1D706}_{N}$. In the case of the PVC, the global mode is typically captured by two modes with the highest energy content at similar energy level. Examining the phase portraits of the corresponding time coefficients reveals whether the modes are indeed periodic. When two such modes occur, the coherent velocity fluctuation of the PVC can be reconstructed by

(3.3)$$\begin{eqnarray}\widetilde{\boldsymbol{u}}(\boldsymbol{x},t)=\text{Re}\left\{\overline{\sqrt{a_{1}^{2}+a_{2}^{2}}}\,(\unicode[STIX]{x1D731}_{1}(\boldsymbol{x})+\text{i}\unicode[STIX]{x1D731}_{2}(\boldsymbol{x}))\text{e}^{-2\unicode[STIX]{x03C0}\text{i}\,f_{PVC}t}\right\},\end{eqnarray}$$

where $\text{Re}$ is the real part and $f_{PVC}$ is the frequency of the global mode. Since the SPIV measurements are not time-resolved with regard to the PVC oscillations, the frequency cannot be obtained via a spectral analysis of the time coefficients. Instead, the frequency is determined from the time-resolved pressure measurements.

3.2 Eigenmodes by global direct linear stability analysis

The LSA in turbulent flows is used to obtain modes of coherent velocity fluctuation. In incompressible flows, the governing equations for the hydrodynamic LSA are derived from the incompressible Navier–Stokes equations and the incompressible continuity equation. The triple decomposition ansatz (3.1) is substituted into both equations and both are time-averaged and phase-averaged. By subtracting the time-averaged set of equations from the phase-averaged set of equations, the governing equations for the coherent velocity fluctuations are obtained (Reynolds & Hussain Reference Reynolds and Hussain1972):

(3.4)$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\widetilde{\boldsymbol{u}}}{\unicode[STIX]{x2202}t}+(\widetilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\,\overline{\boldsymbol{u}}+(\overline{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\widetilde{\boldsymbol{u}} & = & \displaystyle -\,\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}\widetilde{p}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D708}(\unicode[STIX]{x1D735}+\unicode[STIX]{x1D735}^{\top })\widetilde{\boldsymbol{u}})\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70F}_{R}+\underbrace{\unicode[STIX]{x1D70F}_{N}}_{{\approx}0}),\end{eqnarray}$$

(3.5)$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\widetilde{\boldsymbol{u}}=0.\end{eqnarray}$$

The term $\unicode[STIX]{x1D70F}_{N}$ describes the nonlinear interactions of the perturbation with its higher harmonics. Under the assumption that these interactions are very weak, this term is neglected in the following. The term $\unicode[STIX]{x1D70F}_{R}=\langle \boldsymbol{u}^{\prime \prime }\boldsymbol{u}^{\prime \prime }\rangle -\overline{\boldsymbol{u}^{\prime \prime }\boldsymbol{u}^{\prime \prime }}=\widetilde{\boldsymbol{u}^{\prime \prime }\boldsymbol{u}^{\prime \prime }}$ describes the fluctuation of the stochastic Reynolds stresses due to the passage of a coherent perturbation (Reynolds & Hussain Reference Reynolds and Hussain1972). This term has to be modelled in order to close equation (3.4). In the context of the LSA in swirling flows with a PVC, it is well established to use Boussinesq’s eddy viscosity model as a closure (Reau & Tumin Reference Reau and Tumin2002a,Reference Reau and Tuminb; Crouch, Garbaruk & Magidov Reference Crouch, Garbaruk and Magidov2007; Paredes et al. Reference Paredes, Terhaar, Oberleithner, Theofilis and Paschereit2016; Rukes et al. Reference Rukes, Paschereit and Oberleithner2016; Tammisola & Juniper Reference Tammisola and Juniper2016; Kaiser et al. Reference Kaiser, Poinsot and Oberleithner2018). This approach is used here as well. It is assumed that the coherent fluctuations of the turbulent kinetic energy are small enough to be negligible (Reynolds & Hussain Reference Reynolds and Hussain1972). With that assumption, the Reynolds stresses are expressed as

(3.6)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{R}=\widetilde{\boldsymbol{u}^{\prime \prime }\boldsymbol{u}^{\prime \prime }}=-\unicode[STIX]{x1D708}_{t}(\unicode[STIX]{x1D735}+\unicode[STIX]{x1D735}^{\top })\widetilde{\boldsymbol{u}}.\end{eqnarray}$$

The unknown eddy viscosity is calculated from the known velocity field of the experiment. As the turbulence of the swirling jet is highly anisotropic (Rukes et al. Reference Rukes, Paschereit and Oberleithner2016), the approach in (3.6) yields six independent eddy viscosities. A reasonable compromise among the six eddy viscosities can be achieved by using a least-squares fit over all resolved stochastic Reynolds stresses (Ivanova, Noll & Aigner Reference Ivanova, Noll and Aigner2013):

(3.7)$$\begin{eqnarray}\unicode[STIX]{x1D708}_{t}=\frac{\langle -\overline{\boldsymbol{u}^{\prime \prime }\boldsymbol{u}^{\prime \prime }}+{\textstyle \frac{2}{3}}k^{\prime \prime }\unicode[STIX]{x1D644},\overline{\unicode[STIX]{x1D64E}}\rangle _{F}}{2\langle \overline{\unicode[STIX]{x1D64E}},\overline{\unicode[STIX]{x1D64E}}\rangle _{F}},\end{eqnarray}$$

where $\langle \,\cdot \,,\cdot \,\rangle _{F}$ is the Frobenius inner product, $k^{\prime \prime }$ is the turbulent-stochastic kinetic energy, $\unicode[STIX]{x1D644}$ is the identity tensor and $\overline{\unicode[STIX]{x1D64E}}={\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}+\unicode[STIX]{x1D735}^{\top })\overline{\boldsymbol{u}}$ is the mean strain-rate tensor. The eddy viscosity is then simply added to the kinematic viscosity in (3.4) to form an effective viscosity $\unicode[STIX]{x1D708}_{eff}=\unicode[STIX]{x1D708}+\unicode[STIX]{x1D708}_{t}$.

Figure 4.

Eddy viscosity $\unicode[STIX]{x1D708}_{t}$ normalized with molecular viscosity $\unicode[STIX]{x1D708}$, calculated via least-squares fit (3.7), $Re=20\,000$.

Figure 4 shows the eddy viscosity obtained with (3.7). The eddy viscosity is particularly high in the inner and outer shear layers of the swirling jet, with increasing eddy viscosity in the downstream direction. This is adequate to describe the downstream dissipation of the coherent structures that occur due to the interaction with small-scale turbulence. Owing to measurement noise, calculation of the discretized gradients in (3.7) is exposed to that same noise. In order to ameliorate the noisy gradient fields, a box blur low-pass filter is used with a discrete kernel size of $3$ corresponding to every neighbouring point, prior to the calculation of the eddy viscosity. Very small negative values only occur at $4\times 10^{-3}\,\%$ of the nodes of the entire domain and are confined to the boundary layer flow inside the duct. All of these negative values are rigorously set to zero.

The linearized Navier–Stokes and continuity equation for the coherent fluctuation become

(3.8)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\widetilde{\boldsymbol{u}}}{\unicode[STIX]{x2202}t}+(\widetilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\overline{\boldsymbol{u}}+(\overline{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\widetilde{\boldsymbol{u}}=-\,\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}\widetilde{p}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D708}_{eff}(\unicode[STIX]{x1D735}+\unicode[STIX]{x1D735}^{\top })\widetilde{\boldsymbol{u}}), & \displaystyle\end{eqnarray}$$

(3.9)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\widetilde{\boldsymbol{u}}=0. & \displaystyle\end{eqnarray}$$

These equations can be rewritten as

(3.10)$$\begin{eqnarray}\boldsymbol{{\mathcal{B}}}\frac{\unicode[STIX]{x2202}\widetilde{\boldsymbol{q}}}{\unicode[STIX]{x2202}t}-\boldsymbol{{\mathcal{A}}}\,\widetilde{\boldsymbol{q}}=0,\end{eqnarray}$$

where $\boldsymbol{{\mathcal{A}}}$ and $\boldsymbol{{\mathcal{B}}}$ represent the operators of (3.8) and (3.9) and $\widetilde{\boldsymbol{q}}=[\widetilde{\boldsymbol{u}},\widetilde{p}]^{\top }$ is the combined vector of the velocities and pressure.

The global LSA examines flows that are inhomogeneous in two or three spatial dimensions. These are typically named biglobal and triglobal LSA, respectively (Theofilis Reference Theofilis2003). Within the scope of the PVC, a biglobal analysis suffices due to the homogeneity in the azimuthal direction. Equation (3.10) is solved with a normal-mode ansatz in cylindrical coordinates:

(3.11)$$\begin{eqnarray}\widetilde{\boldsymbol{q}}(\boldsymbol{x},t)=\text{Re}\{\hat{\boldsymbol{q}}(x,r)\text{e}^{\text{i}(m\unicode[STIX]{x1D703}-\unicode[STIX]{x1D714}t)}\},\end{eqnarray}$$

where $m$ is the azimuthal wavenumber and $\unicode[STIX]{x1D714}$ is the complex angular frequency. Since the PVC is a single-helical mode, the azimuthal wavenumber is set to $m=1$. Discretization and rearrangement lead to a generalized eigenvalue problem with $\unicode[STIX]{x1D714}$ as the eigenvalue, which can be written in the form of

(3.12)$$\begin{eqnarray}\unicode[STIX]{x1D63C}\hat{\boldsymbol{q}}=\unicode[STIX]{x1D714}\unicode[STIX]{x1D63D}\hat{\boldsymbol{q}}\end{eqnarray}$$

in which $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}$ are the discretized operators of (3.10), including (3.11). Solving (3.12) provides the eigenmodes $\hat{\boldsymbol{q}}$, each accompanied with one complex eigenvalue $\unicode[STIX]{x1D714}$, respectively.

It consists of a real part $\text{Re}(\unicode[STIX]{x1D714})$ that corresponds to the angular frequency of the mode and an imaginary part $\text{Im}(\unicode[STIX]{x1D714})$ that corresponds to the growth rate of the mode. The mode is stable when $\text{Im}(\unicode[STIX]{x1D714})<0$, marginally stable when $\text{Im}(\unicode[STIX]{x1D714})=0$ and unstable when $\text{Im}(\unicode[STIX]{x1D714})>0$. In the eigenspectrum, an oscillator mode at limit cycle, such as the PVC, is expected to be represented by an eigenvalue isolated from any continuous eigenvalue branch and ideally marginally stable ($\text{Im}(\unicode[STIX]{x1D714})=0$) since the instability neither grows nor decays (Barkley Reference Barkley2006). With that criterion and the known frequency from the experiment, the PVC mode can be identified and written as

(3.13)$$\begin{eqnarray}\widetilde{\boldsymbol{u}}(\boldsymbol{x},t)=\text{Re}\{\hat{\boldsymbol{u}}(x,r)\text{e}^{\text{i}(\unicode[STIX]{x1D703}-\text{Re}(\unicode[STIX]{x1D714})t)}\}.\end{eqnarray}$$

For discretizing and solving the generalized eigenvalue problem (3.12) a Fortran code by Paredes (Reference Paredes2014) is employed that uses the Arnoldi algorithm. The domain is decomposed into three rectangular subdomains with conforming interfaces (Demaret & Deville Reference Demaret and Deville1991). A high-order finite difference scheme with non-uniformly distributed nodes is used for the discretization of the axial and radial directions, which is able to provide a converged solution at comparatively low mesh resolution (Paredes et al. Reference Paredes, Hermanns, Clainche and Theofilis2013). Figure 5 shows the grid for the final resolution with a total number of $70\,104$ grid nodes. The three subdomains are marked in terms of colour. The grid spacing in the axial and radial directions is particularly small inside the duct (domain $1$) and in the vicinity of the duct exit (lower part of domains $2$ and $3$) to accurately capture the PVC mode. The grid spacing is then expanded in the downstream and far-field directions.

Figure 5.

Mesh grid of the discretized domain for the global LSA including the boundaries $\unicode[STIX]{x1D6E4}$; the three rectangular subdomains with conforming interfaces are marked in terms of black/grey shade.

Homogeneous Neumann boundary conditions for velocity and pressure are imposed at the inlet $\unicode[STIX]{x1D6E4}_{in}$. This is necessary since the coherent perturbations are low but not zero at the upstream boundary, as observed in the experiments (figure 10). For the same reason, homogeneous Neumann boundary conditions for the velocity and pressure are set at the axial outlet $\unicode[STIX]{x1D6E4}_{out,x}$ in the direct LSA case. Here, the computational domain is truncated at a position with remaining non-zero perturbations. In the adjoint LSA case, homogeneous Dirichlet conditions are set at the axial outlet $\unicode[STIX]{x1D6E4}_{out,x}$ since no adjoint perturbations are allowed to enter the domain. At the radial outlet $\unicode[STIX]{x1D6E4}_{out,r}$, homogeneous Dirichlet boundary conditions are set. For all walls $\unicode[STIX]{x1D6E4}_{wall}$, homogeneous Dirichlet boundary conditions are imposed for the velocity due to the no-slip and no-penetration conditions. For the pressure, no physical boundary conditions exist. However, a compatibility condition from the governing momentum equations can be derived (Theofilis, Duck & Owen Reference Theofilis, Duck and Owen2004). The homogeneous Dirichlet conditions for the velocity are substituted into the linearized Navier–Stokes equation for the coherent perturbation (3.8). Knowing that the eddy viscosity is zero per definition on the wall and taking the inner product with the unit normal vector $\boldsymbol{n}$ at the respective boundaries, rearrangement leads to

(3.14)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}n}=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}u_{n}}{\unicode[STIX]{x2202}n^{2}},\quad \boldsymbol{x}\in \unicode[STIX]{x1D6E4}_{wall}.\end{eqnarray}$$

Assuming $\unicode[STIX]{x2202}^{2}u_{n}/\unicode[STIX]{x2202}n^{2}\approx 0$ provides homogeneous Neumann conditions for the pressure on the walls. Since the global LSA is conducted in a cylindrical coordinate system, the del operator $\unicode[STIX]{x1D735}$ in (3.8) and (3.9) exhibits a singularity at $r=0$. To ensure smoothness and boundedness on the centreline, the boundary conditions are set in agreement with Khorrami, Malik & Ash (Reference Khorrami, Malik and Ash1989). Homogeneous Neumann boundary conditions are set for $\hat{v}$ and ${\hat{w}}$ on the axis $\unicode[STIX]{x1D6E4}_{axis}$, whereas homogeneous Dirichlet conditions are set for $\hat{u}$ and $\hat{p}$. All boundary conditions are summarized in table 2.

Table 2.

Boundary conditions of the global LSA.

3.3 Receptivity and structural sensitivity by global adjoint linear stability analysis

The global adjoint LSA is employed to estimate the receptivity of the PVC. This provides a theoretical prediction where the PVC is most efficiently controlled by periodic forcing. Furthermore, the adjoint LSA enables one to calculate the structural sensitivity, which reveals the location of the wavemaker of the PVC that is responsible for the continuous self-excitation of the instability.

The discretized operator $\unicode[STIX]{x1D63C}$ in (3.12) is typically non-normal, especially in shear flows. This non-normality causes the eigenvectors of the system to be strongly non-orthogonal (Sipp et al. Reference Sipp, Marquet, Meliga and Barbagallo2010). When an arbitrary vector is expanded in the basis of eigenvectors, a dual basis is required due to this non-orthogonality (Salwen & Grosch Reference Salwen and Grosch1972; Tumin & Fedorov Reference Tumin and Fedorov1984; Tumin Reference Tumin1996; Oden & Demkowicz Reference Oden and Demkowicz2017). In the case of the global LSA, the dual of the direct eigenbasis is the basis of adjoint eigenvectors. These are obtained from the generalized adjoint eigenvalue problem (Luchini & Bottaro Reference Luchini and Bottaro2014)

(3.15)$$\begin{eqnarray}\unicode[STIX]{x1D63C}^{\text{H}}\hat{\boldsymbol{q}}^{+}=\unicode[STIX]{x1D714}^{+}\unicode[STIX]{x1D63D}^{\text{H}}\hat{\boldsymbol{q}}^{+},\end{eqnarray}$$

where $(\cdot )^{\text{H}}$ denotes the Hermitian transpose and $(\cdot )^{+}$ denotes an adjoint quantity (note that $\unicode[STIX]{x1D714}^{+}=\unicode[STIX]{x1D714}^{\text{H}}$). In the case of a converged adjoint solution, the adjoint eigenvalue $\unicode[STIX]{x1D714}_{j}^{+}$ of the adjoint eigenmode $\hat{\boldsymbol{q}}_{j}^{+}$ is the complex conjugate of the direct eigenvalue $\unicode[STIX]{x1D714}_{j}$. By that, the pair of mutual direct and adjoint eigenmodes, $\hat{\boldsymbol{q}}_{j}$ and $\hat{\boldsymbol{q}}_{j}^{+}$, are identified.

The bases of the direct and adjoint eigenvectors form a biorthogonal system (Oden & Demkowicz Reference Oden and Demkowicz2017). Let the Hermitian inner product on a complex vector space be defined by $\langle \unicode[STIX]{x1D648}\boldsymbol{x},\boldsymbol{y}\rangle =(\unicode[STIX]{x1D648}\boldsymbol{x})^{\text{H}}\boldsymbol{y}=\langle \boldsymbol{x},\unicode[STIX]{x1D648}^{\text{H}}\boldsymbol{y}\rangle$. Considering the generalized direct eigenvalue problem of (3.12) and the generalized adjoint eigenvalue problem of (3.15), the biorthogonality condition can be easily derived and reads (Luchini & Bottaro Reference Luchini and Bottaro2014)

(3.16)$$\begin{eqnarray}(\unicode[STIX]{x1D714}_{k}-(\unicode[STIX]{x1D714}_{j}^{+})^{\text{H}})\langle \hat{\boldsymbol{q}}_{j}^{+},\unicode[STIX]{x1D63D}\hat{\boldsymbol{q}}_{k}\rangle =0.\end{eqnarray}$$

If $k=j$, the expression inside the round brackets vanishes since $\unicode[STIX]{x1D714}_{k}=(\unicode[STIX]{x1D714}_{j}^{+})^{\text{H}}$. If $k\neq j$, it directly follows that $\langle \hat{\boldsymbol{q}}_{j}^{+},\unicode[STIX]{x1D63D}\hat{\boldsymbol{q}}_{k}\rangle =0$. In other words, all direct and adjoint eigenvectors are pairwise orthogonal with regard to $\unicode[STIX]{x1D63D}$, except if $k=j$. This property of biorthogonality can be used (a) to expand the solution of the inhomogeneous linearized Navier–Stokes equation for the coherent fluctuation (equation (3.10) including a source term on the right-hand side) in order to characterize the receptivity of an instability, and (b) to expand the solution of the perturbed generalized eigenvalue problem (equation (3.12) including small perturbations of the discretized operator $\unicode[STIX]{x1D63C}$) in order to localize the wavemaker of an instability.

For (a) characterizing the receptivity to periodic open-loop forcing, a momentum source on the right-hand side of (3.10) is introduced. The source term is assumed to be helical in the azimuthal direction, to be harmonic in time and not to interact with the homogeneous system. Hence, the eigenspectrum is not modified by the source and the source can be interpreted as an open-loop forcing. The initial value problem (3.10) is then modified and written as

(3.17)$$\begin{eqnarray}\boldsymbol{{\mathcal{B}}}\frac{\unicode[STIX]{x2202}\widetilde{\boldsymbol{q}}}{\unicode[STIX]{x2202}t}-\boldsymbol{{\mathcal{A}}}\,\widetilde{\boldsymbol{q}}=\hat{\boldsymbol{f}}\text{e}^{\text{i}(\unicode[STIX]{x1D703}-\unicode[STIX]{x1D714}_{f}t)},\end{eqnarray}$$

with $\hat{\boldsymbol{f}}$ being the spatial structure of the forcing and $\unicode[STIX]{x1D714}_{f}$ being the angular forcing frequency. The solution is derived making use of the biorthogonality condition in (3.16). The derivation is not shown here for brevity and the interested reader is referred to Chandler (Reference Chandler2011) for a complete derivation.

Let $\hat{\boldsymbol{q}}_{1}$ be the eigenmode of the PVC and $\hat{\boldsymbol{q}}_{0}$ be the initial condition. It is assumed that the system already oscillates with the PVC eigenmode in the natural state, prior to applying open-loop forcing, i.e. $\hat{\boldsymbol{q}}_{0}=\hat{\boldsymbol{q}}_{1}$. When oscillating close to limit cycle, the decay rate is much smaller than the decay rate of all the remaining modes, $|\text{Im}(\unicode[STIX]{x1D714}_{1})|\ll |\text{Im}(\unicode[STIX]{x1D714}_{j\neq 1})|$. Also, the forcing is assumed to be of constant amplitude, thus $\text{Im}(\unicode[STIX]{x1D714}_{f})=0$. For $t$ becoming very large, the discrete solution of (3.17) can then be approximated by

(3.18)

The resulting spatial mode shape of the system is not altered by the forcing $\hat{\boldsymbol{f}}$ and is the PVC mode $\hat{\boldsymbol{q}}_{1}$ itself. The temporal response is both at the PVC frequency $\unicode[STIX]{x1D714}_{1}$ and at the forcing frequency $\unicode[STIX]{x1D714}_{f}$ and proportional to $|R_{f}|$. The $|R_{f}|$ value increases when the spatial structure of the forcing $\hat{\boldsymbol{f}}$ approaches the spatial structure of the adjoint PVC mode $\hat{\boldsymbol{q}}_{1}^{+}$. Additionally, $|R_{f}|$ increases when the forcing frequency $\unicode[STIX]{x1D714}_{f}$ approaches the PVC frequency $\unicode[STIX]{x1D714}_{1}$. According to the Cauchy–Schwarz inequality, $|R_{f}|\propto |\langle \hat{\boldsymbol{q}}_{1}^{+},\hat{\boldsymbol{f}}\rangle |\leqslant |\hat{\boldsymbol{q}}_{1}^{+}||\,\hat{\boldsymbol{f}}|$ holds. Therefore, the magnitude of the adjoint mode can be physically interpreted as the receptivity of the system to a given open-loop periodic forcing. Note that, for $\hat{\boldsymbol{f}}=\mathbf{0}$, the original PVC mode as an eigensolution of (3.10), i.e. $\widetilde{\boldsymbol{q}}=\hat{\boldsymbol{q}}_{1}$, is retrieved.

It has to be noted that the receptivity to open-loop forcing defined by (3.18) is only valid in the sense it was mathematically introduced. The most severe restriction is the assumption that open-loop actuation only acts as a source term in the coherent momentum equation, without interacting with the PVC mode itself and without modifying the mean flow. As will be shown later, this assumption is not valid in general, and nonlinear interaction and mean-flow modifications do indeed occur in the process of synchronization. The adjoint theory has, therefore, its limits within this assumption, and will only give a gradient-like estimation but not a complete description of the synchronization process.

In order to locate the wavemaker of the PVC instability, the region of strongest intrinsic feedback is sought. The feedback can be modelled as a source term in (3.12), similar to (3.17) but this time proportional to the perturbation $\widetilde{\boldsymbol{q}}$ itself such that

(3.19)$$\begin{eqnarray}\boldsymbol{{\mathcal{B}}}\frac{\unicode[STIX]{x2202}\widetilde{\boldsymbol{q}}}{\unicode[STIX]{x2202}t}-\boldsymbol{{\mathcal{A}}}\,\widetilde{\boldsymbol{q}}=\boldsymbol{{\mathcal{C}}}\,\widetilde{\boldsymbol{q}}.\end{eqnarray}$$

This results in a perturbed generalized eigenvalue problem due to a structural perturbation of the linear operator $\boldsymbol{{\mathcal{A}}}$ . This is equivalent to a modification of the base or mean flow and results in a change of the eigenvalue. Considering only spatially localized structural perturbations, the region where the resulting eigenvalue change is the highest is the region where the intrinsic feedback is the strongest. This eigenvalue change can be calculated at first order, and a bound of this change at every position is obtained by utilizing the Cauchy–Schwarz inequality. The upper bound is typically called ‘structural sensitivity’ and, in $L^{2}$-norm, reads (Giannetti & Luchini Reference Giannetti and Luchini2007)

(3.20)$$\begin{eqnarray}\unicode[STIX]{x1D6EC}(x,r)=|\hat{\boldsymbol{u}}^{+}(x,r)||\hat{\boldsymbol{u}}(x,r)|.\end{eqnarray}$$

The position of the wavemaker, which drives the self-excitation of the instability, corresponds to the region where the structural sensitivity is significantly high.

4.1 Mean flow and frequency spectra of the swirling jet

Figure 6 shows the non-dimensionalized mean flow for all three velocity components as contours overlaid with in-plane streamlines. The axial component displays the strong jet discharging into the unconfined ambience with a slow decay downstream of the duct exit. The surrounding air is entrained by the jet into the outer shear layer, which becomes particularly visible on inspecting the transverse mean velocity component. The transverse component also reveals the strong initial cross-sectional spreading of the jet in the radial direction when the flow exits the duct. From $x/D=0.19$ to $1.71$ a prominent breakdown bubble can be made out in which recirculation occurs. It entails a wake-like velocity deficit downstream. The jet encloses the bubble, and an inner shear layer between bubble and jet is formed. Inside the confined duct, a long extending region of axial velocity deficit in the upstream direction exists as well. This significantly reduced velocity around the centreline is related to the fact that the global mode of the PVC also reaches far upstream into the duct, as will be discussed further below. To the knowledge of the authors, this important property has not been observed in any similar swirling jet set-up. However, this is definitely not a necessary condition for the onset of a PVC instability. For example, in Rukes et al. (Reference Rukes, Paschereit and Oberleithner2016), upstream of the breakdown bubble, the axial centreline velocity has indeed an excess instead of a deficit compared to the axial velocities at greater radial positions. The out-of-plane component of the mean flow shows the high rate of azimuthal rotation in the positive $\unicode[STIX]{x1D703}$-direction inside the duct. Downstream of the duct exit, the out-of-plane momentum decays and is spread in the radial direction due to the cross-sectional expansion, similar to the axial momentum.

Figure 6.

Normalized mean flow: (a) axial, (b) transverse and (c) out-of-plane; $Re=20\,000$.

In the following, the PVC mode is examined in detail. In figure 7, the power spectra are considered for all resolvable azimuthal modes of $m=-2$, $-1$, $0$, $1$ and $2$, as a function of the non-dimensional frequency, the Strouhal number $St=fD/u_{0}$. The spectra are based on the time-resolved pressure measurements, as explained in § 2. For $m=1$, a dominant peak at $St\approx 0.7$ is evident for all three Reynolds numbers, corresponding to dimensional frequencies of $f=60~\text{Hz}$, $78~\text{Hz}$ and $115~\text{Hz}$, respectively. These are the frequencies of the global, single-helical PVC mode. The first harmonic of the PVC manifests in the double-helical mode $m=2$ at double the natural PVC frequency, $St\approx 1.4$, with lesser amplitude. A slight residual peak also appears at the natural PVC frequency for $m=2$. This is due to an imperfect calibration of the pressure transducers. For the counter-rotating azimuthal mode of $m=-1$, no distinct peak exists. For $m=-2$, the same spectrum as for $m=2$ shows up. This is due to the indeterminate phase at the limiting Nyquist wavenumber. However, it can be assumed that the rotational direction is in the positive $\unicode[STIX]{x1D703}$-direction.

Another smaller peak is apparent for $m=1$. It occurs at $St\approx 0.4$ and is possibly related to another global mode, which has also been observed in previous works (Terhaar et al. Reference Terhaar, Reichel, Schrödinger, Rukes, Paschereit and Oberleithner2013; Sieber, Paschereit & Oberleithner Reference Sieber, Paschereit and Oberleithner2017). In the present configuration, this mode appears to be subcritical and it occurs in the data due to stochastic forcing. Since the peak is around two orders of magnitude lower than the PVC peak, the impact on the PVC dynamics is presumably low. At $m=0$, a small peak is visible at $St\approx 0.35$. It is presumably related to a weak pumping motion of the vortex breakdown bubble that has also been observed in previous studies (Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011).

Figure 7.

Power spectral density over Strouhal number $St$ for azimuthal modes $m=-2$ to $2$ and varied Reynolds numbers $Re$.

4.2 The precessing vortex core and mechanisms of its formation

For $Re=20\,000$, figure 8 shows that the first POD mode pair contains much more kinetic energy than the rest of any single mode or mode pair. This applies to the coefficients inside and outside of the duct. Furthermore, the phase portrait of the time coefficients in figure 9 demonstrates the harmonic nature of the first POD mode pair, which is linked to the dominant peak observed in the power spectrum of figure 7. It also explicitly shows that the harmonic nature of the PVC is clearly present inside the duct. Since the SPIV measurements are not time-resolved, the POD modes cannot be spectrally decomposed, in contrast to spectral POD (as e.g. proposed in Sieber, Paschereit & Oberleithner (Reference Sieber, Paschereit and Oberleithner2016) or Towne, Schmidt & Colonius (Reference Towne, Schmidt and Colonius2018)). However, the energy optimality of the POD combined with the energy separation of the PVC relative to all other modes ensures that the leading mode pair will mainly contain spectral parts that are related to the PVC and nothing else significant. This has been demonstrated in previous studies, such as in Terhaar et al. (Reference Terhaar, Reichel, Schrödinger, Rukes, Paschereit and Oberleithner2013) and Tammisola & Juniper (Reference Tammisola and Juniper2016).

Figure 8.

The POD mode energy $\unicode[STIX]{x1D706}_{j}$ normalized with total turbulent kinetic energy $\sum _{k}\unicode[STIX]{x1D706}_{k}$,$Re=20\,000$ (●, external domain with $x/D>0$; ▫, internal domain with $x/D<0$).

Figure 9.

Phase portrait of normalized POD time coefficients for the first two most energetic POD modes, $Re=20\,000$ (●, external domain with $x/D>0$; ▫, internal domain with $x/D<0$).

Figure 10.

Normalized POD mode: (a) axial, (b) transverse and (c) out-of-plane; $Re=20\,000$.

Figure 11.

Normalized LSA mode: (a) axial, (b) transverse and (c) out-of-plane; $Re=20\,000$.

The spatial mode shape of the mode pair, according to (3.3), is shown in figure 10 at arbitrary phase angle. The mode features non-zero velocity fluctuations in the transverse and out-of-plane component at the jet axis. For a kinematic reason, this is solely possible for azimuthal modes of $m=1$ (Khorrami et al. Reference Khorrami, Malik and Ash1989). This further establishes that the structure of this POD mode pair is associated with the dominant peak identified in the pressure spectrum at $m=1$ (figure 7).

The PVC can be conceived as a global precession instability with which the convective Kelvin–Helmholtz instabilities synchronize. Features of both instabilities are clearly visible in the contour plots. The spatial fluctuation of the axial and transverse coherent velocity downstream of the duct exit indicates the helical shear layer vortices of the Kelvin–Helmholtz instabilities. The amplitude of the vortices decays in the downstream direction after reaching a maximum around $x/D\approx 0.5$. The out-of-plane and transverse coherent velocities along the centreline describe the precession motion of the vortex core (Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011). Most strikingly, for all three components, it is evident that the PVC mode extends far upstream into the duct, with non-zero amplitudes up to the boundary of the resolved measurement domain. The spatial amplification of the mode is already initiated here. Similar to the centreline velocity deficit of the mean flow inside the duct, to the knowledge of the authors, this important property has not been observed before.

In the following, the term $(\widetilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\overline{\boldsymbol{u}}$ in the Navier–Stokes equations for the coherent part (3.4) is examined. This term is directly related to the production of coherent perturbations (Sipp et al. Reference Sipp, Marquet, Meliga and Barbagallo2010) and it gives insight into the mechanisms that initiate and generate the PVC. Figure 12 shows the magnitude of the production of coherent perturbations for each component separately. A clear spatial separation can be seen between the bulk of production of axial and the bulk of production of transverse as well as out-of-plane perturbations. The former are produced in the outer and inner shear layers of the swirling jet. While the outer shear layer is only present outside of the duct, the inner shear layer associated with the mean axial velocity deficit on the centreline is already forming inside the duct. The production of axial perturbations is attributed to the transverse velocity gradient of the mean axial velocity $\unicode[STIX]{x2202}\overline{u}/\unicode[STIX]{x2202}y$, which is strong in the shear layers in the vicinity of the duct exit. In contrast, the transverse velocity gradients of the mean transverse velocity $\unicode[STIX]{x2202}\overline{v}/\unicode[STIX]{x2202}y$ and mean out-of-plane velocity $\unicode[STIX]{x2202}\overline{w}/\unicode[STIX]{x2202}y$ are particularly strong on the centreline upstream of the duct exit. This gradient is responsible for the high production of both transverse and out-of-plane perturbations in this region. The significant bulk of production starts around $x/D\approx -0.8$. This region is responsible for the initiation of the precession motion. The impact of the absolute instability region of the wavemaker, which will be considered further below, apparently reaches this far upstream into the duct. Comparing the maximum values of production, it further becomes clear that the most significant contribution appears on the centreline inside the duct. This suggests that a forcing applied in these regions of high production should have the largest effect on the PVC dynamics. This will, indeed, be confirmed by the global adjoint LSA in § 4.4.

Figure 12.

Coherent momentum term associated with production of perturbations: (a) axial, (b) transverse and (c) out-of-plane; $Re=20\,000$.

4.3 Theoretical prediction of the precessing vortex core

The eigenvalues of the global direct LSA for all three Reynolds numbers are displayed in figure 13. The entire spectrum exhibits stable eigenvalues only. The identified PVC modes are clustered around the experimentally measured Strouhal number of $St\approx 0.7$. They are discrete and do not belong to any continuous branch. With increasing mesh resolution, monotonic convergence of the selected modes is not achieved. Instead, the selected eigenvalues converge to an average value and then start to ‘oscillate’ within a fixed limit when the mesh resolution is further increased. The lack of asymptotic convergence of the eigenvalues can be explained by the limited spatial resolution of the SPIV data. When the mesh of the global LSA is finer than the grid of the SPIV, the changes of the eigenvalues become small but seemingly ‘random’ since the addition of grid nodes no longer provides any additional information. These ‘random’ changes are attributed to discretization errors caused by the employed interpolation scheme coupled with the choice of the finite difference scheme in the LSA. Therefore, an uncertainty is associated with the eigenvalues. The error bars in figure 13 designate the standard deviation of the oscillations for both the non-dimensional frequency $\text{Re}(St)$ and growth rate $\text{Im}(St)$ for the selected eigenvalues representing the PVC.

In all three cases, the growth rate is close to the stability limit but below zero. However, it is more than an order of magnitude smaller than the frequency, which means that their time scales are well separated. Furthermore, the influence of inaccuracies in the mean-flow measurements and eddy viscosity estimation further affect the prediction of the true growth rate. Therefore, within uncertainty, the eigenmode of the PVC can be assumed to be marginally stable.

Other discrete modes emerge at $St\approx 0.4$, isolated from the continuous branches. This mode is probably related to the small additional peak that is observed in the pressure spectrum as outlined in § 4.1. This mode may be an actual subcritical mode that is excited by stochastic forcing (Farrell & Ioannou Reference Farrell and Ioannou1993). All other remaining modes are identified as numerical, non-physical modes since their corresponding eigenvalues are much more sensitive with regard to discretization errors than the physical modes (Sipp et al. Reference Sipp, Marquet, Meliga and Barbagallo2010). Furthermore, their decay rates increase with increasing mesh resolution, i.e. they become more stable.

Figure 13.

Eigenvalue spectrum for azimuthal wavenumber $m=1$ with real part of Strouhal number $\text{Re}(St)$ denoting the non-dimensional frequency, imaginary part of Strouhal number $\text{Im}(St)$ denoting the non-dimensional growth rate; selected PVC mode (filled markers ●, ▪, ▴); error bars designating standard deviation of frequency and growth rate; stability limit (horizontal dashed line).

Table 3 lists the explicit eigenvalues of the PVC obtained by the global direct LSA with their corresponding standard deviations in non-dimensional form. The LSA frequency converted to a dimensional quantity is compared to the experimentally measured frequency (as obtained by the pressure measurements, see § 4.1). For all three cases, the relative error between LSA and experimental frequency is always below $7\,\%$ within the standard deviation range.

Table 3.

Eigenvalues of the global direct LSA and comparison of LSA frequency $\text{Re}(f_{LSA})$ and experimental frequency $f_{exp}$ (see § 4.1); $\pm$ denotes standard deviation of the converged solution that starts ‘oscillating’ around an average value when mesh resolution is further increased.

The PVC eigenmode is shown in figure 11 for all three components at $Re=20\,000$. In comparison to the POD mode in figure 10, a very good match exists, although there are some minor discrepancies. The amplitude of the precession motion around the centreline is underestimated inside and outside of the duct, compared to the amplitude of the Kelvin–Helmholtz vortices in the outer shear layer. However, apart from that, the overall spatial shape of the mode is captured very well, including the quantitative growth and decay of the helical Kelvin–Helmholtz vortices. Furthermore, the initiation of the precession motion inside the duct is also very well predicted by the model.

Summarizing, the relative error of the predicted frequency is very small. Furthermore, the spatial eigenmode shape is in good agreement with the experimentally obtained POD mode. The results demonstrate the validity of the global LSA in the context of turbulent swirling jets.

Figure 14.

Normalized magnitude of the adjoint LSA mode vector indicating receptivity to open-loop forcing (a) and normalized structural sensitivity indicating sensitivity to mean-flow modifications (b); the closed red line denotes regions of high structural sensitivity ${>}0.92$ ($0.21<x/D<0.43$); the red triangle (▴) denotes stagnation point of breakdown bubble ($x/D=0.19$, $y/D=0$), $Re=20\,000$.

4.4 Receptivity and structural sensitivity

Figure 14(a) shows the magnitude of the adjoint PVC mode. From the start of the measured domain around $x/D\approx -1.75$ downstream to $x/D\approx -1.5$, the receptivity is very low and close to zero. Hence, the PVC should only respond very weakly or not at all when actuation is applied in these regions, at least when the introduced perturbations quickly decay in the downstream direction. In this case, advected residual perturbations do not survive long enough to couple with the PVC in a more receptive region. A quick decay will occur if there are no convective instabilities. Owing to the observed trends, even lower receptivities can be expected upstream of the resolved measurement domain. In the downstream direction, the highest receptivities are located from $x/D\approx -1$ to $-0.4$. In this region, manipulation of the PVC should be most efficient via open-loop actuation. It coincides exceptionally well with the region of maximum production of coherent perturbations (see figure 12). This also explains why the receptivity is primarily accumulated close to the centreline while being almost zero close to the walls of the duct. An open-loop forcing introduced in that region would be subjected to a strong convective amplification. The largest convective growth occurs when the actuation is introduced at the start of the production around $x/D\approx -0.8$. This point overlaps with the maximum of the adjoint magnitude being reached here as well. Further downstream, the receptivity decreases slowly up to the core of the breakdown bubble. When actuation is applied further downstream, the decreasing spatial length of convective growth can be interpreted to be responsible for that trend. Thereafter, the receptivity quickly drops to zero.

The adjoint mode compares well with the adjoint mode obtained by Qadri et al. (Reference Qadri, Mistry and Juniper2013), who also studied the PVC in a swirling flow for a Reynolds number two orders of magnitude lower. In their work, the amplitude of the adjoint mode reaches its maximum upstream of the breakdown bubble and concentrated around the centreline. This is also valid for the adjoint mode shown in figure 14. In their work, it is interpreted that the high receptivity is attributed to conservation of angular momentum. Any vorticity perturbation introduced in this region of high receptivity is subjected to a compression of the fluid element in the axial direction and a stretching in the radial and azimuthal directions (due to deceleration of the flow closely upstream of the breakdown bubble). Correspondingly, the radial and azimuthal vorticity increases while the axial vorticity decreases. This is basically the same mechanism as suggested above and related to the strong production of coherent fluctuations as shown in figure 12.

Based on structural perturbation theory, the structural sensitivity (figure 14b) shows where a perturbation of the linear stability operator induces the strongest change to the mode. In other words, the structural sensitivity quantifies where in the flow the mode is most sensitive to a change of the mean flow. This is in a way equivalent to the adjoint mode which quantifies where in the flow the mode is most sensitive to an introduction of an open-loop source term.

The structural sensitivity is typically interpreted to reveal how strong the internal coupling and feedback is between receptivity of the mode and self-excitation of the mode. In the regions where the highest values are reached, the so-called ‘wavemaker’ can be located. For the PVC, it is situated inside the breakdown bubble closely downstream of its stagnation point. The position roughly agrees with the determination of the wavemaker of similar swirling jet configurations from local LSA as well as global LSA. In Rukes et al. (Reference Rukes, Paschereit and Oberleithner2016) and Tammisola & Juniper (Reference Tammisola and Juniper2016), the wavemaker is located closely upstream of the stagnation point, whereas in Kaiser et al. (Reference Kaiser, Poinsot and Oberleithner2018) the wavemaker region is approximately situated at the stagnation point itself.

5.1 Modifications to the mean flow

In this section, the changes to the mean flow that are caused by the forcing and the associated mechanisms leading to synchronization are discussed.

Examining the mean-flow fields for different forcing frequencies at the actuator position $x_{a}/D=-0.5$ reveals only small mean-flow changes. Figure 16 indicates that the mean axial velocity on the centreline is slightly increased for below-natural forcing whereas it is slightly decreased for above-natural forcing. This is associated with small shifts of the breakdown bubble stagnation point. Below-natural forcing displaces the stagnation point downstream, above-natural forcing upstream. The mean azimuthal velocities do not change significantly. Thus, the changes to the mean flow are marginal in the responsive regime.

Figure 16.

Mean velocity profiles at different axial positions $x/D$ for baseline (black line), below-natural forcing $f_{f}/f_{n}=0.9$ (blue line) and above-natural forcing $f_{f}/f_{n}=1.08$ (red line) for actuator position $x_{a}/D=-0.5$ when synchronization is established, for axial component (a) and out-of-plane component (b), $Re=15\,000$.

In contrast, figure 17 clearly shows that, for the most upstream actuator position $x_{a}/D=-2$ in the non-responsive regime, the mean axial centreline velocity significantly increases at a given position for both below- and above-natural forcing. This is associated with the breakdown bubble shrinking and its stagnation point being displaced downstream. Furthermore, the maximum of the mean out-of-plane velocity is decreased in both cases.

Figure 17.

Mean velocity profiles at different axial positions $x/D$ for baseline (black line), below-natural forcing $f_{f}/f_{n}=0.9$ (blue line) and above-natural forcing $f_{f}/f_{n}=1.08$ (red line) for actuator position $x_{a}/D=-2$ when synchronization is established, for axial component (a) and azimuthal component (b), $Re=15\,000$.

In figure 18(a) the swirl number of the forced cases $S$ normalized on the natural swirl number of the baseline case $S_{n}$ is displayed as a function of normalized forcing frequency $f_{f}/f_{n}$ for actuator positions $x_{a}/D=-2$ (non-responsive regime) and $x_{a}/D=-0.5$ (responsive regime). Strikingly, the asymmetry from the tuning diagram (figure 15) is also clearly present here for $x_{a}/D=-2$. For below-natural forcing, the swirl number decreases almost continuously when the frequency shift $|\,f_{n}-f_{f}|$ is increased. On the other hand, for above-natural forcing, the swirl number is significantly reduced for every forcing frequency. The swirl number does not converge towards the natural swirl number when the frequency shift goes to zero. Furthermore, for the above-natural branch, the swirl numbers are all close to the critical swirl number. The critical swirl number defines the bifurcation point that quantifies the threshold at which a steady PVC can occur for the first time when increasing the swirl number. Although the critical swirl number of the natural unmodified mean flow may not be directly comparable to the critical swirl number of the forced modified mean flow, this observation suggests that the natural PVC is suppressed in these cases.

Figure 18.

Swirl number modifications $S/S_{n}$ for synchronization as a function of normalized forcing frequency $f_{f}/f_{n}$ (a) and as a function of synchronization amplitude $A_{s}$ (b); actuator positions $x_{a}/D=-2$ and $x_{a}/D=-0.5$.

In figure 18(b), in semilog scaling, the swirl number is shown as a function of synchronization amplitude. For $x_{a}/D=-2$, the swirl number decreases approximately linearly with increasing synchronization amplitude for each Reynolds number, respectively. The stronger the actuation, the lower the swirl number becomes. With this observation, the asymmetry between below- and above-natural forcing in the non-responsive regime can be easily explained now. The asymmetry is due to the fact that the swirl number is a function only of the synchronization amplitude in the non-responsive regime. For below-natural forcing, the swirl number reduction favours the goal of reducing the frequency of the PVC. The rotation rate of the jet, and thus the swirl number, is decreased until the natural frequency associated with the swirl number matches the forcing frequency. For above-natural forcing, the swirl number reduction acts adversely to the goal of increasing the frequency of the PVC. Ultimately, a threshold of forcing amplitude needs to be exceeded in order to reach the critical swirl number and suppress the natural PVC. At this point the global oscillations are completely dictated by the actuation itself. What is then seen as a ‘PVC’ are only Kelvin–Helmholtz vortices that still exist in the outer shear layer and that are excited at the dominating forcing frequency. Furthermore, the precession of the vortex core is artificially induced by the helical excitation of the actuator. This helical forcing substitutes the precession motion that is otherwise initiated through self-excitation in the region of absolute instability in the baseline case.

For the actuator positions of the responsive regime, it can be assumed that a strongly nonlinear interaction occurs between actuation mode and natural PVC mode due to the high receptivity at the location of actuation (see figure 14). Additionally, the trends of swirl number modification as a function of forcing frequency are clearly different from the non-responsive regime (see figure 18a). For above-natural forcing, the swirl number decreases linearly. For below-natural forcing, the swirl number even increases with forcing frequencies close to the natural frequency, before linearly decreasing with a larger frequency shift. Hence, in the responsive regime, there is an impact on the mean flow, even though small. Inspecting the swirl number as a function of synchronization amplitude (figure 18b), another qualitative difference from the non-responsive regime can be made out. Not only is the swirl number a function of synchronization amplitude but it also depends on whether the forcing is below or above the natural frequency, constituting two different branches. This suggests that there is a significant dynamical response, and not only an alteration of the mean flow. Synchronization is presumably reached by a classical phase lock-in mechanism via an inverse saddle-node bifurcation (Balanov et al. Reference Balanov, Janson, Postnov and Sosnovtseva2008). Mean-flow modifications occur but are weak. They are responsible for the persisting slight asymmetry between the lower and upper branches. Again, below-natural forcing has a favouring effect while above-natural forcing has an adverse effect.

Figure 19.

Changes to the stagnation-point position $x_{SP}/D$ (a) and changes to the axial breakdown bubble length $L_{bubble}$ normalized with baseline bubble length $L_{bubble,n}$ as functions of normalized forcing frequency $f_{f}/f_{n}$ (b); actuator positions $x_{a}/D=-2$ and $x_{a}/D=-0.5$.

Further evidence on the differences between responsive and non-responsive regimes can be determined when changes to the recirculation bubble are evaluated (see figure 19). For the non-responsive regime, the bubble is displaced in the downstream direction in every case. Since one main effect of the actuator is to dissipate rotational energy, which leads to a swirl number decrease, the downstream displacement can be attributed to a delay of the vortex breakdown when the frequency shift of the forcing is increased. Accordingly, the length of the breakdown bubble in the axial direction (measured from upstream stagnation point to downstream stagnation point on the centreline) is reduced. Oberleithner et al. (Reference Oberleithner, Paschereit, Seele and Wygnanski2012) shows that the size of the breakdown bubble correlates with the size of the region of absolute instability. Thus, one main effect of the actuation in the non-responsive regime can be hypothesized to be a reduction of this region of absolute instability. For above-natural forcing in the non-responsive regime, the actuation is strong enough that it leads to a complete suppression of the instability. In the responsive regime, the bubble is barely displaced and the length only changes insignificantly. Therefore, the region of absolute instability likely remains untouched. This suggests that the natural mechanism is not altered, apart from the slight swirl number modifications.

The mean-flow modifications occurring in the non-responsive regime clearly suggest that synchronization is not achieved by a direct interaction between actuation and PVC mode. Instead, the manipulation of the swirl number is responsible for synchronization. Either the decrease of the natural PVC frequency (below-natural forcing) or the suppression of the natural PVC (above-natural forcing) leads to synchronization. Strictly speaking, in the above-natural case, there is only a substitution of the helical PVC by the helical actuation and no synchronization takes place, since there is only one oscillator remaining in the system.

5.2 Global linear stability of the mean flow at synchronization

In the previous section (§ 5.1), it was concluded that synchronization in the non-responsive regime is dictated by mean-flow modifications. In order to further support this finding, a global LSA is conducted on the forced, synchronized mean-flow fields. In the non-responsive regime, the global LSA should predict the correct PVC frequencies for below-natural forcing (swirl number is reduced until the modified natural frequency matches the forcing frequency) and incorrect frequencies for above-natural forcing (PVC is suppressed). Note that the LSA can only be conducted for the cases with the actuator positioned at $x_{a}/D=-2$ since these are the only cases where a sufficient extent of the internal duct flow was measured. This upstream domain is required for valid results of the global LSA.

Figure 20.

Global LSA on synchronized mean flows with predicted PVC frequency (obtained from LSA) compared to measured PVC frequency (obtained by experiment) (a) and compared to normalized swirl number (b); $x_{a}/D=-2$.

Figure 20 shows the results of the global LSA. In figure 20(a), the frequency obtained by the LSA is compared to the frequency obtained by the experiment. For below-natural forcing, the predicted PVC frequency approximately coincides with the measured PVC frequency. The predicted frequencies are closely gathered around the dashed identity line. Ideally, they would collapse exactly onto the identity line. However, asymptotic convergence of the solutions is not achieved and the numerical values are coupled to an uncertainty, as discussed in § 3.2. Furthermore, this effect is even more pronounced for the forced cases. As will be shown in § 5.3, the coherent fluctuations are significantly large due to the forcing at the most upstream position of the SPIV resolved domain. Therefore, the setting of the inlet boundary conditions to homogeneous Neumann conditions (see § 3.2) is questionable. For above-natural forcing, the predicted frequencies are not in line with the measured frequencies. While the measured PVC frequencies are above the natural PVC frequency, the global LSA predicts frequencies below. This demonstrates that in this case there is no actual synchronization due to swirl number modification but due to suppression. This fact is further clarified in figure 20(b). The predicted LSA frequencies for above-natural forcing are closely gathered around the critical swirl number. This shows that the PVC is suppressed. For below-natural forcing, the predicted frequencies ($S/S_{n}>0.8$) increase with increasing swirl number. This is in line with the general relationship of the PVC frequency scaling proportionally to the swirl number. Both observations support that, in the non-responsive regime, synchronization is solely reached by mean flow and, thus, swirl number manipulation.

5.3 Modifications to the coherent flow

In order to gain further insight into the response of the PVC, the modifications to the coherent fluctuations are now investigated. This will reveal whether the PVC is significantly influenced by the forcing or not.

At first, the total coherent kinetic energy (KE) of the PVC is used to characterize its changes caused by the forcing. It is calculated by

(5.1)$$\begin{eqnarray}\widetilde{K}(x)={\textstyle \frac{1}{2}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\infty }|\widetilde{\boldsymbol{u}}|^{2}r\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

Figure 21 shows the total coherent KE for the natural baseline case and for a representative below- and above-natural forced case at both $x_{a}/D=-2$ and $-0.5$, respectively. In the natural case, the coherent KE is almost zero at the most upstream position $x/D=-1.75$ and the amplitude is slowly convectively amplified inside the entire duct. Outside of the duct, further substantial amplification of the PVC occurs due to the generation of coherent vortices by Kelvin–Helmholtz instabilities and saturation is reached at $x/D\approx 0.5$.

When the actuator forces at $x_{a}/D=-2$ (non-responsive regime), the introduced perturbations increase the total coherent KE in the vicinity of the actuator. However, these perturbations decay downstream. Hence, there are no convective instabilities that are amplified at this far-upstream location. For below-natural forcing, the coherent KE decays until it reaches the level of the natural case at $x/D\approx -0.8$, before the coherent KE starts to grow again. This turning point coincides with the start of major production of azimuthal perturbations observed in figure 12. Downstream of this position at $x/D\approx -0.8$, the convective amplification is very similar to the natural case, however, with slightly higher saturation amplitude. The similarity corroborates again that the intrinsic natural instability mechanism is barely manipulated in the case of below-natural forcing in the non-responsive regime.

For above-natural forcing, the coherent KE generated by the actuator also decays downstream until around $x/D\approx -1$. Thereafter, the energy stays at almost constant level. When the flow exits the duct, amplification of coherent fluctuations still occurs as in the natural case. However, the amplification is only related to Kelvin–Helmholtz vortices in the outer shear layer and the saturation amplitude is significantly reduced compared to the natural saturation amplitude. This observation supports the notion that the natural PVC instability mechanism is suppressed since no convective amplification occurs inside the duct, which is in contrast to the natural and below-natural case. The precession motion is not self-excited but upheld by the continuous introduction of helical perturbations by the actuator.

For actuator position $x_{a}/D=-0.5$ (responsive regime), below-natural forcing does not alter the coherent energy significantly. For above-natural forcing, the axial position where saturation sets in is only slightly shifted downstream. In both cases, the saturation energy is equal to the baseline case. This supports the hypothesis that in the responsive regime the synchronization is primarily achieved by a nonlinear modal interaction between actuation and PVC. A significant alteration of the natural instability mechanism does not occur.

Figure 21.

Total coherent kinetic energy $\widetilde{K}$ as a function of axial coordinate $x/D$ for natural PVC and synchronized PVC; actuator positions $x_{a}/D=-2$ and $-0.5$; $Re=15\,000$.

Now, the alterations to the production mechanism responsible for the spatial growth of the PVC mode are inspected in more detail. In figure 22 the production term is shown for below- and above-natural forcing when the actuator is placed at $x_{a}/D=-2$. Only the out-of-plane component is considered as representative for the precession motion (see § 4.2). For below-natural forcing, high production occurs in the vicinity of the actuator but decreases in the downstream direction, reaching almost zero around $x/D\approx -1$. Downstream, the production increases again. A very similar distribution to the natural baseline case is evident, with high levels of production in the inner and outer shear layers outside of the duct. Furthermore, regions of high production are found in accordance with the baseline case close to the centreline but slightly shifted downstream. Thus, the naturally acting production mechanism is not seriously altered.

In contrast, for above-natural forcing, the natural production mechanism is suppressed, as can be observed. The suppression is associated with the subcritical swirl number mentioned above. The production in the inner and outer shear layers is reduced. The maximum production occurs upstream where the actuator is situated and only decreases in magnitude downstream. In contrast to the below-natural case, there is no recovery and an increase of production in the downstream direction does not happen. Hence, the bulk of production of helical perturbations originates from the actuation itself and it serves as a substitute for the self-excited perturbations that are naturally generated in the baseline case. These substitute perturbations from the actuator are then simply convected downstream (not shown), artificially sustaining the precession motion.

Figure 22.

Impact of forcing at $x_{a}/D=-2$ (non-responsive regime) on coherent out-of-plane momentum term associated with production of perturbations, below-natural (a), $f_{f}/f_{n}=0.9$, and above-natural (b), $f_{f}/f_{n}=1.08$, forcing; natural case plotted side by side with forced case; $Re=15\,000$.

Figure 23.

Impact of forcing at $x_{a}/D=-0.5$ (responsive regime) on coherent out-of-plane momentum term associated with production of perturbations, below-natural (a), $f_{f}/f_{n}=0.9$, and above-natural (b), $f_{f}/f_{n}=1.08$ forcing; natural case plotted side by side with forced case, $Re=15\,000$.

Figure 23 displays the production of coherent fluctuations in the responsive regime for below- and above-natural forcing, respectively, compared side by side with the natural case. In both cases, the natural mechanism is barely touched. The almost symmetric response shows that mean-flow modifications play a minor role. Additionally, the fact that the production mechanism is not significantly altered suggests that synchronization cannot be attained by mean-flow modification but, instead, must be reached via nonlinear interactions between actuation and PVC mode. This poses a strong contrast to the alterations of the production mechanism in the non-responsive regime.

5.4 Comparison of experimental and theoretical receptivities

Since the impact of open-loop control on the flow field has been investigated in detail, the experimental and theoretical receptivities can now be compared and discussed. The theoretical receptivity was defined to be proportional to the magnitude of the adjoint mode according to § 3.3. Consistently, the experimental receptivity is now defined to be proportional to the gradient of synchronization amplitude $A_{s}$ with respect to the frequency shift $\unicode[STIX]{x0394}f=|\,f_{n}-f_{f}|$, i.e. $\unicode[STIX]{x0394}A_{s}/\unicode[STIX]{x0394}f$. A larger gradient means that higher amplitudes are necessary to achieve a defined frequency shift, indicating that the flow is less receptive at the current actuator position compared to another actuator position where the gradient is lower.

The results of the synchronization diagram (figure 15) generally agree with the receptivity obtained via the adjoint PVC mode in § 4.4. The non-responsive regime corresponds to the adjoint PVC mode being close to zero upstream of $x/D\approx -1.5$, suggesting extremely low receptivity. Therefore, the introduced perturbations decay before the PVC responds to the open-loop actuation and only the modifications to the mean flow have a lasting effect. The mean-flow modifications dictate the new frequency of the system, either by a favourable reduction of the swirl number (below-natural forcing) or by suppression of the natural PVC in conjunction with an artificial substitution of the precession motion through the actuation mode (above-natural forcing).

The responsive regime can be associated with the regions of high magnitude of the adjoint PVC mode, indicating maximum receptivity in the region of $x/D\approx -1$ to $-0.4$. Likewise, the synchronization amplitudes decrease from $x/D=-1$ to $-0.75$, suggesting an increasing receptivity of the PVC in the downstream direction. Within the uncertainty range, the flow seems to have a similar receptivity for actuator positions $x_{a}/D=-0.75$ and $-0.5$. In these regions of high receptivity, the introduced perturbations couple with the natural production mechanism of the PVC. The perturbations are convectively amplified and overwhelm the natural dynamics, overwriting the natural frequency of the system.

The results of the synchronization studies as well as the adjoint mode quantifying the receptivity coincide with the experimental lock-in results obtained by Kuhn et al. (Reference Kuhn, Moeck, Paschereit and Oberleithner2016). With a comparable set-up, but the Reynolds number being an order of magnitude lower, a similar slightly asymmetrical response of the PVC is observed, requiring lower synchronization amplitudes for below-natural forcing compared to above-natural forcing. This is observed for a large number of axial positions. In the case where the actuator is traversed downstream of the breakdown bubble stagnation point, the synchronization amplitudes steeply increase. This is in line with the magnitude of the adjoint mode (figure 14) decreasing downstream of the stagnation point, indicating a declining receptivity. Likewise, for the actuator located upstream of the stagnation point within the vicinity of the nozzle exit, the synchronization amplitudes stay almost constant, indicating that the receptivity of the PVC does not change much here. This can also be observed for the adjoint mode where the magnitude is approximately constant for a larger axial extent in the corresponding region.

5.5 Synchronization mechanisms

In the following, the physical mechanisms that lead to synchronization in the responsive and non-responsive regimes are summarized. The key mechanisms are visualized in the flowchart of figure 24. The forcing can have two effects on the flow: it either modifies the mean flow or it nonlinearly (directly) interacts with the PVC. The primary mean-flow modification is the reduction of the swirl number by dissipation of mean azimuthal momentum. For sufficiently weak forcing, the resulting swirl number stays above the critical swirl number and the PVC frequency is decreased. This leads to frequency pulling for below-natural forcing and to frequency pushing for above-natural forcing. When the pulling becomes strong enough, synchronization can be reached. As for the pushing, synchronization can never be reached. Instead, the forcing needs to be strong enough such that the swirl number falls below the critical swirl number and the PVC is suppressed. Then, the only remaining ‘global mode’ is artificially generated by the continuous actuation, substituting the PVC. The nonlinear interaction, on the other hand, leads to a frequency pulling in any case. When the forcing is strong enough, phase lock-in sets in and synchronization is attained.

Figure 24.

Schematic of physical mechanisms leading to synchronization or substitution of the PVC.

In the non-responsive regime, only mean-flow modifications are responsible for synchronization. This is associated with a region of virtually zero receptivity that prevents a direct nonlinear interaction. In the responsive regime, a superposition of mean-flow modification and nonlinear interaction occurs. With higher receptivity of the PVC, synchronization is increasingly achieved by nonlinear interaction and mean-flow modifications start to play a less important role.