2.1 Weakly nonlinear analysis

The parameter $S$ in (2.8) is chosen as the control parameter that is varied to change the characteristics of the flow. We presume that there exists a critical swirl number, $S_{c}$, beyond which the flow transitions from a stable steady state for $S<S_{c}$ to an unsteady state for $S>S_{c}$. We determine the characteristics of this unsteady state for small deviations of $S$ from $S_{c}$ when compared to the flow oscillation amplitude, as follows:

(2.9)$$\begin{eqnarray}S=S_{c}+\unicode[STIX]{x1D6E5}_{s}\unicode[STIX]{x1D716}^{2},\end{eqnarray}$$

where $\unicode[STIX]{x1D716}$ is the normalized amplitude of the unsteady state for $S>S_{c}$ and $\unicode[STIX]{x1D6E5}_{s}\sim O(1)$. Thus, the matrix ${\mathcal{B}}^{s}$ in (2.8) can be written as

(2.10)$$\begin{eqnarray}{\mathcal{B}}^{s}={\mathcal{B}}+\unicode[STIX]{x1D6E5}_{s}\unicode[STIX]{x1D716}^{2}{\mathcal{B}}_{1},\end{eqnarray}$$

where ${\mathcal{B}}=\text{diag}(1,S_{c},1,0)$ and ${\mathcal{B}}_{1}=\text{diag}(0,1,0,0)$.

We also presume that the unsteady solution can be characterized by a combination of variations over a ‘fast’ time scale $t_{1}=t$ and a ‘slow’ time scale $t_{2}=\unicode[STIX]{x1D716}^{2}t$. Thus, $\tilde{\boldsymbol{q}}$ is now expressed as an asymptotic series in terms of $\unicode[STIX]{x1D716}$ as follows:

(2.11)$$\begin{eqnarray}\displaystyle \tilde{\boldsymbol{q}}(r,\unicode[STIX]{x1D703},z,t_{1},t_{2}) & = & \displaystyle \boldsymbol{q}_{o}(r,z)+\unicode[STIX]{x1D716}\boldsymbol{q}_{1}(r,\unicode[STIX]{x1D703},z,t_{1},t_{2})+\unicode[STIX]{x1D716}^{2}\boldsymbol{q}_{2}(r,\unicode[STIX]{x1D703},z,t_{1},t_{2})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D716}^{3}\boldsymbol{q}_{3}(r,\unicode[STIX]{x1D703},z,t_{1},t_{2})+\cdots \,,\end{eqnarray}$$

where we have assumed that the contributions to $\tilde{\boldsymbol{q}}$ at each order depend on both time scales. Each of the nonlinear operators in (2.8) can be expanded in powers of $\unicode[STIX]{x1D716}$ as follows:

(2.12)$$\begin{eqnarray}\displaystyle & {\mathcal{N}}\{\boldsymbol{q}\}={\mathcal{N}}\{\boldsymbol{q}_{o}\}+\unicode[STIX]{x1D716}{\mathcal{N}}\{\boldsymbol{q}_{1}\}+\unicode[STIX]{x1D716}^{2}{\mathcal{N}}\{\boldsymbol{q}_{2}\}+\cdots \,, & \displaystyle\end{eqnarray}$$

(2.13)$$\begin{eqnarray}\displaystyle & {\mathcal{N}}^{s}\{\boldsymbol{q}\}={\mathcal{N}}^{s}\{\boldsymbol{q}_{o}\}+\unicode[STIX]{x1D716}{\mathcal{N}}^{s}\{\boldsymbol{q}_{1}\}+\unicode[STIX]{x1D716}^{2}{\mathcal{N}}^{s}\{\boldsymbol{q}_{2}\}+\cdots \,, & \displaystyle\end{eqnarray}$$

(2.14)$$\begin{eqnarray}\displaystyle & {\mathcal{N}}^{ss}\{\boldsymbol{q}\}={\mathcal{N}}^{ss}\{\boldsymbol{q}_{o}\}+\unicode[STIX]{x1D716}{\mathcal{N}}^{ss}\{\boldsymbol{q}_{1}\}+\unicode[STIX]{x1D716}^{2}{\mathcal{N}}^{ss}\{\boldsymbol{q}_{2}\}+\cdots \,. & \displaystyle\end{eqnarray}$$

Substituting (2.9) and (2.11)–(2.14) into (2.8), yields the following equation at $O(1)$:

(2.15)$$\begin{eqnarray}\displaystyle {\mathcal{N}}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{o}+S_{c}{\mathcal{N}}^{s}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{o}+S_{c}^{2}{\mathcal{N}}^{ss}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{o}={\mathcal{L}}_{T}\boldsymbol{q}_{o}+S_{c}{\mathcal{L}}_{T}^{s}\boldsymbol{q}_{o}. & & \displaystyle\end{eqnarray}$$

Comparing the above with (2.4) shows that $\boldsymbol{q}_{o}$ is just the time-averaged state, $\bar{\boldsymbol{Q}}$, at $S=S_{c}$.

The $O(\unicode[STIX]{x1D716})$ terms yield the linearized evolution equation for $\boldsymbol{q}_{1}$,

(2.16)$$\begin{eqnarray}\displaystyle & & \displaystyle {\mathcal{B}}\frac{\unicode[STIX]{x2202}\boldsymbol{q}_{1}}{\unicode[STIX]{x2202}t_{1}}+{\mathcal{N}}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{1}+{\mathcal{N}}\{\boldsymbol{q}_{1}\}\boldsymbol{q}_{o}+S_{c}\left({\mathcal{N}}^{s}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{1}+{\mathcal{N}}^{s}\{\boldsymbol{q}_{1}\}\boldsymbol{q}_{o}\right)\nonumber\\ \displaystyle & & \displaystyle \quad +\,S_{c}^{2}\left({\mathcal{N}}^{ss}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{1}+{\mathcal{N}}^{ss}\{\boldsymbol{q}_{1}\}\boldsymbol{q}_{o}\right)={\mathcal{L}}_{T}\boldsymbol{q}_{1}+S_{c}{\mathcal{L}}_{T}^{s}\boldsymbol{q}_{1}.\end{eqnarray}$$

Equation (2.16) is the linearized evolution equation for flow unsteadiness from the time-averaged state at $S=S_{c}$. For convenience of discussion in this paper, we rewrite (2.16) using a single operator to represent all terms with spatial derivatives as follows:

(2.17)$$\begin{eqnarray}\displaystyle \left({\mathcal{B}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t_{1}}+{\mathcal{L}}\right)\boldsymbol{q}_{1}=0, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{L}}$ is the linear spatial operator acting on $\boldsymbol{q}_{1}$ (see (A 8) in appendix A). Note that (2.17) does not have any term that depends on $t_{2}$. Therefore, the general solution to (2.17) can be written as a superposition of harmonically oscillating azimuthal Fourier modes. However, in this paper, we assume that the flow unsteadiness is characterized by a single neutrally stable Fourier mode as follows:

(2.18)$$\begin{eqnarray}\boldsymbol{q}_{1}(r,z,\unicode[STIX]{x1D703},t_{1},t_{2})=A_{1}(t_{2})\hat{\boldsymbol{q}}_{1m}(r,z)\text{e}^{\text{i}m\unicode[STIX]{x1D703}-\text{i}\unicode[STIX]{x1D714}_{g,m}t_{1}}+(A_{1}(t_{2})\hat{\boldsymbol{q}}_{1m}(r,z))^{\ast }\text{e}^{-\text{i}m\unicode[STIX]{x1D703}+\text{i}\unicode[STIX]{x1D714}_{g,m}t_{1}},\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{g,m}$ is a real number and the superscript, ‘*’, here and in the rest of this paper, denotes the complex conjugate. The mode shape, $\hat{\boldsymbol{q}}_{m}(r,z)=[\hat{u} _{r},\hat{u} _{\unicode[STIX]{x1D703}},\hat{u} _{z},\hat{p}]^{\text{T}}$, is the spatial mode shape of the $m$th Fourier mode that is determined by the following:

(2.19)$$\begin{eqnarray}\displaystyle (-\text{i}\unicode[STIX]{x1D714}_{g,m}{\mathcal{B}}+{\mathcal{L}}_{m})\hat{\boldsymbol{q}}_{1m}=0, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{L}}_{m}$ is the linearized spatial operator corresponding to the $m$th azimuthal mode. The solutions to (2.19) satisfy the following boundary conditions. At the centreline ($r=0$), $\hat{\boldsymbol{q}}_{1m}(r,z)$ must satisfy kinematic compatibility conditions as follows (Batchelor & Gill (Reference Batchelor and Gill1962)):

(2.20)$$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}c@{}}\hat{u} _{r}=\hat{u} _{\unicode[STIX]{x1D703}}=0\\ {\displaystyle \frac{\text{d}\hat{u} _{z}}{\text{d}r}}={\displaystyle \frac{\text{d}\hat{p}}{\text{d}r}}=0\end{array}\right\}\quad \text{if }m=0, & \displaystyle\end{eqnarray}$$

(2.21)$$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}c@{}}\hat{u} _{z}=\hat{p}=0\\ \hat{u} _{r}+\text{i}m\hat{u} _{\unicode[STIX]{x1D703}}=0\\ {\displaystyle \frac{\text{d}\hat{u} _{r}}{\text{d}r}}=0\end{array}\right\}\quad \text{if }|m|=1, & \displaystyle\end{eqnarray}$$

(2.22)$$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}c@{}}\hat{u} _{r}=\hat{u} _{\unicode[STIX]{x1D703}}=0\\ \hat{u} _{z}=\hat{p}=0\end{array}\right\}\quad \text{if }|m|>1. & \displaystyle\end{eqnarray}$$

Since the flow is unconfined, the following boundary conditions must be satisfied at the far-field boundary:

(2.23a,b)$$\begin{eqnarray}\hat{\boldsymbol{q}}_{m}\rightarrow 0,\quad r\rightarrow \infty .\end{eqnarray}$$

Thus, equations (2.19)–(2.23) represent a global temporal eigenvalue problem where $\unicode[STIX]{x1D714}_{g,m}$ is the temporal eigenvalue and $\hat{\boldsymbol{q}}_{1m}(r,z)$ is the associated global hydrodynamic mode.

The solution for $A_{1}(t_{2})$ in (2.18) must be derived from higher-order asymptotic relations as follows. At $O(\unicode[STIX]{x1D716}^{2})$, we have the following equation for $\boldsymbol{q}_{2}(r,z,\unicode[STIX]{x1D703},t_{1},t_{2})$:

(2.24)$$\begin{eqnarray}\displaystyle \left({\mathcal{B}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t_{1}}+{\mathcal{L}}\right)\boldsymbol{q}_{2} & = & \displaystyle -{\mathcal{N}}\{\boldsymbol{q}_{1}\}\boldsymbol{q}_{1}-S_{c}{\mathcal{N}}^{s}\{\boldsymbol{q}_{1}\}\boldsymbol{q}_{1}-S_{c}^{2}{\mathcal{N}}^{ss}\{\boldsymbol{q}_{1}\}\boldsymbol{q}_{1}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6E5}_{s}({\mathcal{N}}^{s}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{o}+2S_{c}{\mathcal{N}}^{ss}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{o}-{\mathcal{L}}_{T}^{s}\boldsymbol{q}_{o}).\end{eqnarray}$$

The general solution to (2.24) is given by,

(2.25)$$\begin{eqnarray}\displaystyle \boldsymbol{q}_{2}(r,z,\unicode[STIX]{x1D703},t_{1},t_{2}) & = & \displaystyle \unicode[STIX]{x1D6E5}_{s}\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}(r,z)+|A_{1}|^{2}(t_{2})\hat{\boldsymbol{q}}_{A^{\ast }A}(r,z)\nonumber\\ \displaystyle & & \displaystyle +\,\left(A_{1}^{2}(t_{2})\hat{\boldsymbol{q}}_{AA}(r,z)\text{e}^{\text{i}2m\unicode[STIX]{x1D703}-\text{i}2\unicode[STIX]{x1D714}_{g,m}t_{1}}+{A_{1}^{\ast }}^{2}(t_{2})\hat{\boldsymbol{q}}_{AA}^{\ast }(r,z)\text{e}^{-\text{i}2m\unicode[STIX]{x1D703}+\text{i}2\unicode[STIX]{x1D714}_{g,m}t_{1}}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\left(A_{2}(t_{2})\hat{\boldsymbol{q}}_{1m}(r,z)\text{e}^{\text{i}m\unicode[STIX]{x1D703}-\text{i}\unicode[STIX]{x1D714}_{g,m}t_{1}}+A_{2}^{\ast }(t_{2})\hat{\boldsymbol{q}}_{1m}^{\ast }(r,z)\text{e}^{-\text{i}m\unicode[STIX]{x1D703}+\text{i}\unicode[STIX]{x1D714}_{g,m}t_{1}}\right)\!.\quad\end{eqnarray}$$

The first two terms and the third term, written as a pair of complex conjugate terms, are a particular solution to (2.24). The fourth term, also written as a pair of complex conjugate terms, is the solution to the homogeneous equivalent of (2.24). The functions $\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}(r,z)$, $\hat{\boldsymbol{q}}_{AA^{\ast }}(r,z)$ and $\hat{\boldsymbol{q}}_{AA}(r,z)$ all satisfy the same boundary conditions given in (2.20)–(2.23) and are given by solutions to the following three equations.

The function $\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}(r,z)$ is given by the following:

(2.26)$$\begin{eqnarray}{\mathcal{L}}_{0}\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}=-{\mathcal{N}}_{0}^{s}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{o}-2S_{c}{\mathcal{N}}_{0}^{ss}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{o}+{\mathcal{L}}_{T,0}^{s}\boldsymbol{q}_{o}.\end{eqnarray}$$

Note that the source terms in (2.26) are a function of $\boldsymbol{q}_{o}$, i.e. the base flow at $S_{c}$, alone. Therefore, $\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}$ quantifies to leading order in $S-S_{c}$, the modification of $\boldsymbol{q}_{o}$ with increasing $S$. Therefore, we will refer to $\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}$ as the intrinsic base-flow modification function. It is important to note that, in general, equation (2.26) shows that $\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}$ does not depend on any characteristics of the linear instability mode.

The function $\hat{\boldsymbol{q}}_{A^{\ast }A}(r,z)$ is given by the following:

(2.27)$$\begin{eqnarray}\displaystyle {\mathcal{L}}_{0}\hat{\boldsymbol{q}}_{A^{\ast }A} & = & \displaystyle -{\mathcal{N}}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{1m}^{\ast }-{\mathcal{N}}\{\hat{\boldsymbol{q}}_{1m}^{\ast }\}\hat{\boldsymbol{q}}_{1m}-S_{c}({\mathcal{N}}_{-m}^{s}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{1m}^{\ast }+{\mathcal{N}}_{m}^{s}\{\hat{\boldsymbol{q}}_{1m}^{\ast }\}\hat{\boldsymbol{q}}_{1m})\nonumber\\ \displaystyle & & \displaystyle -\,S_{c}^{2}({\mathcal{N}}_{-m}^{ss}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{1m}^{\ast }+{\mathcal{N}}_{m}^{ss}\{\hat{\boldsymbol{q}}_{1m}^{\ast }\}\hat{\boldsymbol{q}}_{1m}).\end{eqnarray}$$

Note that the source terms on the right of (2.27) are independent of $\boldsymbol{q}_{o}$, showing that the function $\hat{\boldsymbol{q}}_{A^{\ast }A}$ quantifies to leading order in $S-S_{c}$, the non-oscillatory, time-averaged distortion that the linear instability imposes on the base flow at $S_{c}$. Therefore, we will refer to $\hat{\boldsymbol{q}}_{A^{\ast }A}$ as the base-flow distortion function.

The function $\hat{\boldsymbol{q}}_{AA}(r,z)$ quantifies the spatial distribution of the amplitude associated with the first harmonic of the linear instability and is given by the following:

(2.28)$$\begin{eqnarray}(-\text{i}2\unicode[STIX]{x1D714}_{g,m}{\mathcal{B}}+{\mathcal{L}}_{2m})\hat{\boldsymbol{q}}_{AA}=-{\mathcal{N}}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{1m}-S_{c}{\mathcal{N}}_{m}^{s}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{1m}-S_{c}^{2}{\mathcal{N}}_{m}^{ss}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{1m}.\end{eqnarray}$$

At $O(\unicode[STIX]{x1D716}^{3})$ we obtain,

(2.29)$$\begin{eqnarray}\displaystyle \left({\mathcal{B}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t_{1}}+{\mathcal{L}}\right)\boldsymbol{q}_{3} & = & \displaystyle -{\mathcal{B}}\frac{\unicode[STIX]{x2202}\boldsymbol{q}_{1}}{\unicode[STIX]{x2202}t_{2}}-\unicode[STIX]{x1D6E5}_{s}{\mathcal{B}}_{1}\frac{\unicode[STIX]{x2202}\boldsymbol{q}_{1}}{\unicode[STIX]{x2202}t_{1}}-{\mathcal{N}}\{\boldsymbol{q}_{1}\}\boldsymbol{q}_{2}-{\mathcal{N}}\{\boldsymbol{q}_{2}\}\boldsymbol{q}_{1}\nonumber\\ \displaystyle & & \displaystyle -\,S_{c}\left({\mathcal{N}}^{s}\{\boldsymbol{q}_{1}\}\boldsymbol{q}_{2}+{\mathcal{N}}^{s}\{\boldsymbol{q}_{2}\}\boldsymbol{q}_{1}\right)\nonumber\\ \displaystyle & & \displaystyle -\,S_{c}^{2}\left({\mathcal{N}}^{ss}\{\boldsymbol{q}_{1}\}\boldsymbol{q}_{2}+{\mathcal{N}}^{ss}\{\boldsymbol{q}_{2}\}\boldsymbol{q}_{1}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6E5}_{s}[{\mathcal{N}}^{s}\{\boldsymbol{q}_{1}\}\boldsymbol{q}_{o}+{\mathcal{N}}^{s}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{1}+2S_{c}({\mathcal{N}}^{ss}\{\boldsymbol{q}_{1}\}\boldsymbol{q}_{o}+{\mathcal{N}}^{ss}\{\boldsymbol{q}_{o}\}\boldsymbol{q}_{1})-{\mathcal{L}}_{T}^{s}\boldsymbol{q}_{1}].\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Physically relevant, non-trivial solutions to (2.29) must be bounded in time and space. As is common, the norm that is used to quantify the magnitude of $\boldsymbol{q}_{3}$ is generated by the following definition of the inner product between two functions:

(2.30)$$\begin{eqnarray}\langle \boldsymbol{q}_{1},\boldsymbol{q}_{2}\rangle =\int _{0}^{Z}\int _{0}^{R}\int _{0}^{2\unicode[STIX]{x03C0}}(\boldsymbol{q}_{1}^{H})^{\text{T}}\boldsymbol{q}_{2}r\,\text{d}\unicode[STIX]{x1D703}\,\text{d}r\,\text{d}z,\end{eqnarray}$$

where the superscript ‘$H$’ denotes the transpose complex conjugate operation.

Thus, the boundedness condition on $\boldsymbol{q}_{3}$ requires that the oscillatory contributions with frequency, $\unicode[STIX]{x1D714}_{g,m}$, to the source term in (2.29) must be orthogonal to the corresponding adjoint eigenfunction of the operator on the left, associated with the inner product defined in (2.30). These adjoint eigenfunctions can be written as $\boldsymbol{q}_{1m}^{\dagger }(r,z,\unicode[STIX]{x1D703},t_{1})=\hat{\boldsymbol{q}}_{1m}^{\dagger }(r,z)\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}_{g,m}^{\dagger }t_{1})}$, where $\unicode[STIX]{x1D714}_{g,m}^{\dagger }$ and $\hat{\boldsymbol{q}}_{1m}^{\dagger }(r,z)$ are given by solving the following equation:

(2.31)$$\begin{eqnarray}(\text{i}\unicode[STIX]{x1D714}_{g,m}^{\dagger }{\mathcal{B}}+{\mathcal{L}}_{m}^{\dagger })\hat{\boldsymbol{q}}_{1m}^{\dagger }=0,\end{eqnarray}$$

where ${\mathcal{L}}_{m}^{\dagger }$ is the adjoint linear operator that satisfies the following relation, $\langle \hat{\boldsymbol{q}}_{1m}^{\dagger },{\mathcal{L}}_{m}\hat{\boldsymbol{q}}_{1m}\rangle =\langle {\mathcal{L}}_{m}^{\dagger }\hat{\boldsymbol{q}}_{1m}^{\dagger },\hat{\boldsymbol{q}}_{1m}\rangle$. This relation, along with (2.19) and (2.31), gives $\unicode[STIX]{x1D714}_{g,m}^{\dagger }=\unicode[STIX]{x1D714}_{g,m}^{\ast }$ (see Schmid & Henningson Reference Schmid and Henningson2001).

Thus, substituting (2.18) and (2.25) into (2.29) and invoking the boundedness condition, yields the evolution equation for $A_{1}(t_{2})$ as follows (Landau & Lifshitz Reference Landau and Lifshitz1959):

(2.32)$$\begin{eqnarray}\displaystyle \frac{\text{d}A_{1}}{\text{d}t_{2}}=\unicode[STIX]{x1D6E5}_{s}B_{A}A_{1}-N_{A}A_{1}|A_{1}|^{2}. & & \displaystyle\end{eqnarray}$$

Equation (2.32) is referred to as the Stuart–Landau amplitude evolution equation. The coefficient in the linear term on the right in (2.32), $B_{A}$, is given by the following inner product:

(2.33)$$\begin{eqnarray}\displaystyle B_{A} & = & \displaystyle \frac{1}{\langle \hat{\boldsymbol{q}}_{1m}^{\dagger },{\mathcal{B}}\hat{\boldsymbol{q}}_{1m}\rangle }\langle \!\hat{\boldsymbol{q}}_{1m}^{\dagger },[\!\text{i}\unicode[STIX]{x1D714}_{g,m}{\mathcal{B}}_{1}\hat{\boldsymbol{q}}_{1m}-{\mathcal{N}}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}-{\mathcal{N}}\{\boldsymbol{q}_{\unicode[STIX]{x1D6E5}}\}\hat{\boldsymbol{q}}_{1m}\nonumber\\ \displaystyle & & \displaystyle -\,{\mathcal{N}}_{0}^{s}\{\hat{\boldsymbol{q}}_{1m}\}\boldsymbol{q}_{o}-{\mathcal{N}}_{m}^{s}\{\boldsymbol{q}_{o}\}\hat{\boldsymbol{q}}_{1m}\nonumber\\ \displaystyle & & \displaystyle -\,S_{c} ({\mathcal{N}}_{0}^{s}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}+{\mathcal{N}}_{m}^{s}\{\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}\}\hat{\boldsymbol{q}}_{1m}\nonumber\\ \displaystyle & & \displaystyle -\,2{\mathcal{N}}_{0}^{ss}\{\hat{\boldsymbol{q}}_{1m}\}\boldsymbol{q}_{o}-2{\mathcal{N}}_{m}^{ss}\{\boldsymbol{q}_{o}\}\hat{\boldsymbol{q}}_{1m} )\nonumber\\ \displaystyle & & \displaystyle -\,S_{c}^{2}({\mathcal{N}}_{0}^{ss}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}+{\mathcal{N}}_{m}^{ss}\{\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}\}\hat{\boldsymbol{q}}_{1m})\nonumber\\ \displaystyle & & \displaystyle +\,{\mathcal{L}}_{T,m}^{s}\hat{\boldsymbol{q}}_{1m}\! ]\!\rangle \!.\end{eqnarray}$$

Thus, it is clear that the growth of flow oscillations, determined by $B_{A}$ in linear term in (2.32), is in turn determined by $\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}$, i.e. the intrinsic base-flow modification for $S>S_{c}$, which therefore, is the cause of the onset of flow oscillations.

The coefficient of the nonlinear term in (2.32), $N_{A}$, is given by the following inner product:

(2.34)$$\begin{eqnarray}\displaystyle N_{A} & = & \displaystyle \frac{1}{\langle \hat{\boldsymbol{q}}_{1m}^{\dagger },{\mathcal{B}}\hat{\boldsymbol{q}}_{1m}\rangle }\langle \!\hat{\boldsymbol{q}}_{1m}^{\dagger },[\!{\mathcal{N}}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{A^{\ast }A}+{\mathcal{N}}\{\hat{\boldsymbol{q}}_{1m}^{\ast }\}\hat{\boldsymbol{q}}_{AA}\nonumber\\ \displaystyle & & \displaystyle +\,{\mathcal{N}}\{\hat{\boldsymbol{q}}_{A^{\ast }A}\}\hat{\boldsymbol{q}}_{1m}+{\mathcal{N}}\{\hat{\boldsymbol{q}}_{AA}\}\hat{\boldsymbol{q}}_{1m}^{\ast }\nonumber\\ \displaystyle & & \displaystyle +\,S_{c} ({\mathcal{N}}_{0}^{s}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{A^{\ast }A}+{\mathcal{N}}_{2m}^{s}\{\hat{\boldsymbol{q}}_{1m}^{\ast }\}\hat{\boldsymbol{q}}_{AA}\nonumber\\ \displaystyle & & \displaystyle +\,{\mathcal{N}}_{m}^{s}\{\hat{\boldsymbol{q}}_{A^{\ast }A}\}\hat{\boldsymbol{q}}_{1m}+{\mathcal{N}}_{-m}^{s}\{\hat{\boldsymbol{q}}_{AA}\}\hat{\boldsymbol{q}}_{1m}^{\ast })\nonumber\\ \displaystyle & & \displaystyle +\,S_{c}^{2} ({\mathcal{N}}_{0}^{ss}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{A^{\ast }A}+{\mathcal{N}}_{2m}^{ss}\{\hat{\boldsymbol{q}}_{1m}^{\ast }\}\hat{\boldsymbol{q}}_{AA}\nonumber\\ \displaystyle & & \displaystyle +\,{\mathcal{N}}_{m}^{ss}\{\hat{\boldsymbol{q}}_{A^{\ast }A}\}\hat{\boldsymbol{q}}_{1m}+{\mathcal{N}}_{-m}^{ss}\{\hat{\boldsymbol{q}}_{AA}\}\hat{\boldsymbol{q}}_{1m}^{\ast })\!]\!\rangle \!.\end{eqnarray}$$

The expression for $N_{A}$ in (2.37) can be decomposed into two individual components such that $N_{A}=N_{A1}+N_{A2}$. The component, $N_{A1}$, is determined by base-flow distortion as follows:

(2.35)$$\begin{eqnarray}\displaystyle N_{A1} & = & \displaystyle \frac{1}{\langle \hat{\boldsymbol{q}}_{1m}^{\dagger },{\mathcal{B}}\hat{\boldsymbol{q}}_{1m}\rangle }\langle \!\hat{\boldsymbol{q}}_{1m}^{\dagger },[\!{\mathcal{N}}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{A^{\ast }A}+{\mathcal{N}}\{\hat{\boldsymbol{q}}_{A^{\ast }A}\}\hat{\boldsymbol{q}}_{1m}\nonumber\\ \displaystyle & & \displaystyle +\,S_{c}({\mathcal{N}}_{0}^{s}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{A^{\ast }A}+{\mathcal{N}}_{m}^{s}\{\hat{\boldsymbol{q}}_{A^{\ast }A}\}\hat{\boldsymbol{q}}_{1m})\nonumber\\ \displaystyle & & \displaystyle +\,S_{c}^{2}({\mathcal{N}}_{0}^{ss}\{\hat{\boldsymbol{q}}_{1m}\}\hat{\boldsymbol{q}}_{A^{\ast }A}+{\mathcal{N}}_{m}^{ss}\{\hat{\boldsymbol{q}}_{A^{\ast }A}\}\hat{\boldsymbol{q}}_{1m})\!]\!\rangle \!.\end{eqnarray}$$

The component, $N_{A2}$, is the component due to harmonic generation as follows:

(2.36)$$\begin{eqnarray}\displaystyle N_{A2} & = & \displaystyle \frac{1}{\langle \hat{\boldsymbol{q}}_{1m}^{\dagger },{\mathcal{B}}\hat{\boldsymbol{q}}_{1m}\rangle }\langle \!\hat{\boldsymbol{q}}_{1m}^{\dagger },[\!{\mathcal{N}}\{\hat{\boldsymbol{q}}_{1m}^{\ast }\}\hat{\boldsymbol{q}}_{AA}+{\mathcal{N}}\{\hat{\boldsymbol{q}}_{AA}\}\hat{\boldsymbol{q}}_{1m}^{\ast }\nonumber\\ \displaystyle & & \displaystyle +\,S_{c}({\mathcal{N}}_{2m}^{s}\{\hat{\boldsymbol{q}}_{1m}^{\ast }\}\hat{\boldsymbol{q}}_{AA}+{\mathcal{N}}_{-m}^{s}\{\hat{\boldsymbol{q}}_{AA}\}\hat{\boldsymbol{q}}_{1m}^{\ast })\nonumber\\ \displaystyle & & \displaystyle +\,S_{c}^{2}({\mathcal{N}}_{2m}^{ss}\{\hat{\boldsymbol{q}}_{1m}^{\ast }\}\hat{\boldsymbol{q}}_{AA}+{\mathcal{N}}_{-m}^{ss}\{\hat{\boldsymbol{q}}_{AA}\}\hat{\boldsymbol{q}}_{1m}^{\ast })\!]\!\rangle \!.\end{eqnarray}$$

Next, we rewrite (2.32) in terms of a single time variable, $t$, by introducing $A(t)=\unicode[STIX]{x1D716}A_{1}(\unicode[STIX]{x1D716}^{2}t)+\unicode[STIX]{x1D716}^{2}A_{2}(\unicode[STIX]{x1D716}^{2}t)$ into (2.32), which then yields the following equation at leading order:

(2.37)$$\begin{eqnarray}\displaystyle \frac{\text{d}A}{\text{d}t}=(S-S_{c})B_{A}A-N_{A}A|A|^{2}. & & \displaystyle\end{eqnarray}$$

Thus, the asymptotic solution to (2.8), accurate to $O(\unicode[STIX]{x1D716}^{3})$, is given by

(2.38)$$\begin{eqnarray}\displaystyle \tilde{\boldsymbol{q}}(r,z,\unicode[STIX]{x1D703},t) & = & \displaystyle \boldsymbol{q}_{o}(r,z)+(S-S_{c})\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6E5}}(r,z)+|A(t)|^{2}\hat{\boldsymbol{q}}_{A^{\ast }A}(r,z)\nonumber\\ \displaystyle & & \displaystyle +\,A(t)\hat{\boldsymbol{q}}_{1m}(r,z)\text{e}^{\text{i}m\unicode[STIX]{x1D703}-\text{i}\unicode[STIX]{x1D714}_{g,m}t}+A^{\ast }(t)\hat{\boldsymbol{q}}_{1m}(r,z)^{\ast }\text{e}^{-\text{i}m\unicode[STIX]{x1D703}+\text{i}\unicode[STIX]{x1D714}_{g,m}t}\nonumber\\ \displaystyle & & \displaystyle +\,(A(t))^{2}\hat{\boldsymbol{q}}_{AA}(r,z)\text{e}^{\text{i}2m\unicode[STIX]{x1D703}-\text{i}2\unicode[STIX]{x1D714}_{g,m}t}+(A^{\ast }(t))^{2}\hat{\boldsymbol{q}}_{A^{\ast }A^{\ast }}(r,z)\text{e}^{-\text{i}2m\unicode[STIX]{x1D703}+\text{i}2\unicode[STIX]{x1D714}_{g,m}t}\nonumber\\ \displaystyle & & \displaystyle +\,O((S-S_{c})^{3/2}).\end{eqnarray}$$

Except for the first two terms, the remaining terms in (2.38) are contributions that become relevant when coherent unsteadiness is established in the flow for $S>S_{c}$. Note that this first happens when the linear hydrodynamic mode, $\hat{\boldsymbol{q}}_{1m}$, becomes neutrally stable, i.e. $\unicode[STIX]{x1D714}_{g,m}$ is purely real. In this case, the nature of flow unsteadiness for $S>S_{c}$ depends on the characteristics of $A(t)$. These can be determined by setting $A(t)=D(t)\text{e}^{\text{i}\unicode[STIX]{x1D719}(t)}$ in (2.37), where $D(t)$ and $\unicode[STIX]{x1D719}(t)$ are real-valued functions that represent the time evolution of the amplitude and argument of $A(t)$. Thus, equating real parts on both sides of (2.37) yields an evolution equation for $D(t)$ as follows:

(2.39)$$\begin{eqnarray}\frac{\text{d}D}{\text{d}t}=(S-S_{c})B_{Ar}D-N_{Ar}D^{3},\end{eqnarray}$$

where $B_{Ar}$ and $N_{Ar}$ are the real parts of the coefficients $N_{A}$ and $B_{A}$, respectively. Thus, the unsteady characteristics of the flow are determined by establishing the behaviour of $D(t)$ as $t\rightarrow \infty$ for a given non-zero initial perturbation, $D_{i}=D(0)$.

Assuming that the flow is nominally stable for $S<S_{c}$, the following types of unsteady flow behaviour can be identified for $S>S_{c}$ based on the signs of $B_{Ar}$ and $N_{Ar}$ as follows.

(i) Stable limit cycle ($B_{Ar}>0$, $N_{Ar}>0$). For any value of $D_{i}>0$, equation (2.39) shows that for $S>S_{c}$ the first term causes $D(t)$ to increase exponentially, while the second term limits this growth. Eventually, $D(t)\rightarrow D_{LC}$ as $t\rightarrow \infty$ where $D_{LC}$, determined by letting $\text{d}D/\text{d}t\rightarrow 0$ in (2.39), is given by

(2.40)$$\begin{eqnarray}D_{LC}=\sqrt{\frac{(S-S_{c})B_{Ar}}{N_{Ar}}}.\end{eqnarray}$$

The change in oscillation frequency of the flow from $\unicode[STIX]{x1D714}_{g,m}$ baseline at $S=S_{c}$ with increasing $S$ can be determined using (2.39) by choosing $A(t)=D_{LC}\text{e}^{\text{i}\unicode[STIX]{x1D714}_{c}t}$ and comparing imaginary parts on both sides to yield:

(2.41)$$\begin{eqnarray}\unicode[STIX]{x1D714}_{c}=(S-S_{c})\left[B_{Ai}-B_{Ar}\left(\frac{N_{Ai}}{N_{Ar}}\right)\right].\end{eqnarray}$$

Thus, equation (2.38) shows that the angular frequency, $\unicode[STIX]{x1D714}_{LC}$, of the limit-cycle oscillation is given by

(2.42)$$\begin{eqnarray}\unicode[STIX]{x1D714}_{LC}=\unicode[STIX]{x1D714}_{g,m}-(S-S_{c})B_{Ai}+(S-S_{c})B_{Ar}\left(\frac{N_{Ai}}{N_{Ar}}\right).\end{eqnarray}$$

Note also that for $S<S_{c}$, both terms on the right of (2.39) cause $D(t)\rightarrow 0$ as $t\rightarrow \infty$ for any value of $D_{i}$. This type of flow behaviour corresponds to a supercritical Hopf bifurcation in the flow state at $S=S_{c}$, i.e. the emergence of a stable limit cycle for $S>S_{c}$ due to a linear hydrodynamic mode of the system becoming neutrally stable. Thus, the spatial region where the integrand in (2.33) is non-zero can now be identified as the nonlinear wavemaker region, i.e. the region of the flow that contributes to the growth of coherent oscillations. Note that we use the term ‘nonlinear’ to emphasise that this definition of the wavemaker is different from that defined from linear, structural sensitivity considerations (Hill Reference Hill1992; Giannetti & Luchini Reference Giannetti and Luchini2007; Juniper & Pier Reference Juniper and Pier2015; Tammisola & Juniper Reference Tammisola and Juniper2016). Further, we refer to the spatial regions where the integrand in (2.34) is non-zero as the nonlinear wavelimiter region because $N_{A}$ determines the impact that the nonlinear terms have in determining the final saturated oscillation amplitude of the limit cycle. Note that (2.34) shows that amplitude saturation is caused by the modification induced on the time-averaged flow for $S>S_{c}$ as well as interaction between the linear eigenmode and its first harmonic.

(ii) Subcritical flow ($B_{Ar}>0$, $N_{Ar}<0$). In this case, equation (2.39) suggests that for $S>S_{c}$, $D(t)\rightarrow \infty$ as $t\rightarrow \infty$ for any value of $D_{i}$. This suggests that $S=S_{c}$ represents a catastrophic loss of stability. As such, the flow dynamics cannot be quantitatively described by the present leading-order nonlinear analysis. On the other hand, for $S<S_{c}$, an unstable limit cycle with an amplitude $D_{LC}=\sqrt{(S-S_{c})B_{Ar}/N_{Ar}}$ exists. Note that for $D_{i}<D_{LC}$, $\text{d}D/\text{d}t<0$ and therefore, $D(t)\rightarrow 0$ as $t\rightarrow \infty$. For $D_{i}>D_{LC}\,\text{d}D/\text{d}t>0$ and therefore, $D(t)\rightarrow \infty$ as $t\rightarrow \infty$. Thus, coherent flow unsteadiness is triggered when $D>D_{LC}$ for $S<S_{c}$.

(iii) Stable flow ($B_{Ar}\leqslant 0$, $N_{Ar}>0$). For this case, equation (2.39) shows that $D(t)\rightarrow 0$ as $t\rightarrow \infty$ for any value of $D_{i}$.

Note that other possibilities for $N_{Ar}$ and $B_{Ar}$ are not considered as they do not result in steady flow behaviour for $S<S_{c}$. We next use these results to gain insight into the emergence of globally self-excited helical oscillations in a variable swirl jet experiment that will be described in the next section.

4.1 Time-averaged characteristics

Figure 2. Variation of the time-averaged axial velocity field with $S$ ($Re=59\,000$) with streamlines. The solid black curve shows the time-averaged recirculation zone. The broken curves show the nominal location of the inner (black) and outer (magenta) shear layers in the time-averaged flow.

The evolution of the time-averaged flow field for $S=0.67{-}1.83$ is shown in figure 2. In each of the figures, the colour field shows $\bar{U}_{z}$. Streamlines in white are overlaid to visualize the direction of the flow. The broken black and magenta lines denote the time-averaged location of the shear layers in the flow. These curves are determined from the radial position of the inflection points in the $\bar{U}_{z}$ profile at each $z$ location. The thick black contours correspond to $\bar{U}_{z}=0$ and denote the boundaries of the recirculation zones. The same scheme is used to label these time-averaged flow field features in all field plots in the rest of this paper.

Figure 2 shows that increasing $S$ leads to increased jet spreading due to the formation of the breakdown bubble, whose boundaries are captured by the thick black $\bar{U}_{z}=0$ contour. Note that the absence of a bubble for the $S=0.67$ case is due to the fact that for this case, the bubble is still intermittent and hence does not appear in the time-averaged result presented in figure 2(a). This change in the time-averaged flow field structure is accompanied by the emergence of a coherent flow oscillation in the jet. We next present results that characterize the nature of the flow oscillation in the next subsection.

4.2 Unsteady dynamics

We characterize the global spatial structure of flow oscillations using the spectral proper orthogonal decomposition (SPOD) on the time-resolved flow field measurements from sPIV. The SPOD calculates the optimal energy-ordered but spectrally resolved basis modes that can be used to reconstruct the unsteady flow field (Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018). Thus, SPOD is a superior method for identifying narrowband oscillations when compared to the traditional proper orthogonal decomposition (POD), which does not discern modes based on frequency. As such, highly tonal but low-energy oscillations that are not necessarily separated from other unsteady motions by the POD are unambiguously identified by the SPOD. All results from SPOD in this paper have been determined by dividing the time record of sPIV measurements into 18 ensembles with 256 points per ensemble and 50 % overlap between successive ensembles for the spectral estimation step.

Figure 3(a–e) shows the mode energy spectra from SPOD for $S=0.67{-}1.83$. The mode energy spectra show the energy spectra for each mode, where mode 1, the most energetic mode, is highlighted in red and the remainder of the modes fade from black to grey in descending order of total mode energy. As $S$ is increased, the peak frequency, which appears in mode 1 between 800 and 1100 Hz for each case, increases as well. Additionally, the energy of the mode increases in amplitude and decreases in spectral width as the coherence of the PVC increases. At the highest two swirl numbers shown, $S=1.43$ and 1.83, some of the PVC energy appears in the second mode (mode 2) as well, where a smaller, more broadband peak is present beneath the mode 1 spectrum. Additionally, the first harmonic of the PVC appears in mode 1 at the higher swirl number cases, a result of the very high amplitude coherent oscillation.

Another feature that appears in mode 1 of the SPOD for swirl numbers in the range $S=0.95{-}1.15$ (see figure 3c,d) is two additional peaks near 500 Hz and 300 Hz. A continuous wavelet transform analysis of the velocity fluctuations at several points in the annular jet near the nozzle exit for these two values of $S$ showed that these peaks are due to an intermittent turbulent burst that occurred at frequencies of 300 Hz and 500 Hz. Further, the pressure transducer signals show no evidence of these oscillations (Manoharan Reference Manoharan2019). Note, however, that the energy associated with these fluctuations is an order of magnitude less than that associated with the coherent oscillation. As such, they are not the focus of the discussion in this paper.

Figure 3. Variation of the energy spectra with $S$ from spectral proper orthogonal decomposition of the time-resolved flow fields determined from sPIV measurements ($Re=59\,000$). The red curve shows the spectrum for the most energetic mode. The spectra for the remaining modes are shown in grey scale from high (dark) to low (light) total energy. The values of $S$ associated with each result are shown within each plot.

Figure 4. Variation of the horizontal in plane velocity component modes determined from spectral proper orthogonal decomposition at the PVC frequency, $f_{PVC}$, with $S$ ($Re=59\,000$). The value of $S$ associated with each result is shown on the title of each of the plots. The solid black curve shows the time-averaged recirculation zone boundary. The broken curves show the nominal location of the inner (black) and outer (magenta) shear layers in the time-averaged flow.

Figure 4(a–f) shows the normalized mode shapes corresponding to the mode 1 peak frequency at each value of $S$. Each figure shows the real part of the horizontal velocity component in the plane of light sheet, $Re(\hat{u} _{x})$, associated with each SPOD mode. These mode shapes show the evolution of the spatial dynamics of the coherent oscillation as a function of $S$. Note first that the mode at $S=0.67$ shows a mode shape with strong oscillations on the centreline, which suggests that the coherent oscillation induces a precession of the axial vortex. Therefore, we identify the coherent oscillation associated with the frequency corresponding to the peak in figure 3(a–e) as the precessing vortex core (Syred Reference Syred2006). We will show in a forthcoming section of this paper that the PVC is a consequence of the formation of the breakdown bubble in the flow, using the results from the theoretical analysis presented in § 2. It is useful at this point to emphasize the distinction between the PVC oscillations shown in figure 4 and helical spiral breakdown studied by Meliga et al. (Reference Meliga, Gallaire and Chomaz2012) and Qadri et al. (Reference Qadri, Mistry and Juniper2013). In these two studies, the helical rollup of the shear layers is due to the instability of helical modes that have significant oscillation amplitudes downstream of the breakdown bubble. The PVC oscillation on the other hand is associated with the precession of the vortex breakdown bubble and results in shear layer rollup that originates at the upstream end of the breakdown bubble as figure 4 and other prior experimental studies show (Escudier & Keller Reference Escudier and Keller1985; Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011).

With increasing $S$, figure 4(a–e) shows that the centreline oscillation moves upstream into the nozzle. This is consistent with the prior observation of Escudier & Keller (Reference Escudier and Keller1985) in a very similar experimental set-up. Also, the region where the centreline oscillations are concentrated lies at the upstream end of the bubble, consistent with observations in prior studies with experimental set-ups similar to ours (Escudier & Keller Reference Escudier and Keller1985; Escudier Reference Escudier1988; Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011). Note also that the oscillations along the OSL become progressively stronger relative to those on the centreline with increasing $S$ due to the development of helical oscillations of the shear layers in the flow downstream of the bubble.

We next estimate the azimuthal mode number associated with the PVC oscillation and the critical swirl number, $S_{c}$, at which it emerges from a spatial cross-spectral analysis of the sPIV results as follows. Figure 5(a) shows the spatial variation of $\overline{u_{z}^{\prime }u_{z}^{\prime }}$ determined from the sPIV measurements at $S=1.83$. The broken black and magenta curves show the nominal position of the ISL and the OSL respectively. Points 1 and 2, positioned at ($r/D=\pm 0.4$, $z/D=0.25$) (see figure 5a), are chosen as reference points for cross-spectral analysis (Bendat & Piersol Reference Bendat and Piersol2011). All cross-spectrum values are computed by dividing the time record of signals at points 1 and 2 into 19 ensembles with 500 samples each and with 50 % overlap between successive ensembles. Figure 5(b) shows a typical result for the variation of the magnitude of the cross-spectrum, $|C_{12}|$, between $u_{z}^{\prime }$ signals at these two points ($S=1.83$). Note that $|C_{12}|$ shows a strong peak (${\sim}0.9$) at 1090 Hz, which is close to the value of the peak frequency in the corresponding SPOD spectra for mode 1 – see figure 3(e) – as may be expected. In the rest of this paper, we define the characteristic PVC oscillation frequency ($f_{PVC}$) as the frequency at which the maximum value of $|C_{12}|$ occurs.

Figure 5. Variation of (a) $\overline{u_{z}^{\prime }u_{z}^{\prime }}$ in $r{-}z$ plane. (b) Coherence ($|C_{12}|$) as a function of frequency for $S=1.83$. The broken curves show the nominal location of the inner (black) and outer (magenta) shear layers in the time-averaged flow for this case.

Figure 5(b) also shows a second peak at the first harmonic, $f_{1h}=2f_{PVC}$, consistently with the corresponding spectrum from the SPOD (figure 4e). The lower value of $|C_{12}|\sim 0.6$ at $f_{1h}$ (see figure 5b) is because of smaller amplitudes of velocity oscillation at $f_{1h}$ and therefore, lower signal-to-noise ratio when compared with $f_{PVC}$. These results are typical and similar results are observed at other spatial locations for this case and for other cases with $S\geqslant 0.67$. The Fourier mode number, $m$, associated with the oscillations at $f_{PVC}$ and $f_{1h}$, can be determined from the phase, $\unicode[STIX]{x1D719}_{12}$, of $C_{12}$ at these frequencies for each $S$. Figure 6(a) shows the variation $\unicode[STIX]{x1D719}_{12}$ with $S$, determined at $f_{PVC}$ and $f_{1h}$. The data in figure 6 correspond to values of $S$ at which $|C_{12}|\geqslant 0.7$ at $f_{PVC}$ and $|C_{12}|\geqslant 0.5$ at $f_{1h}$. Note that $\unicode[STIX]{x1D719}_{12}\sim \unicode[STIX]{x03C0}$ for all $S$, suggesting that the PVC oscillation shown in figure 4 corresponds to $m=1$ helical oscillations, as may be expected. The corresponding result for $f_{1h}$ shows that $\unicode[STIX]{x1D719}_{12}\sim 0$, which suggests that the first harmonic corresponds to an $m=2$ oscillation. These results are typical and, although not shown here, similar variations for $\unicode[STIX]{x1D719}_{12}$ with $S$ are observed at other downstream locations as well. Also, note that this difference in $m$ values associated with the PVC and its first harmonic is consistent with the general asymptotic solution for stable limit-cycle oscillations determined from the theoretical analysis, equation (2.38).

Figure 6. Variation of the phase of the cross-spectrum between axial velocity signals at points 1 and 2, ($\unicode[STIX]{x1D719}_{12}$) at $f_{PVC}$ and $f_{1h}=2f_{PVC}$, with $S$ ($Re=59\,000$).

Figure 7. Variation of (a) $f_{PVC}$ and (b) $|\hat{u} _{z}|^{2}$ as a function of $S$. The solid lines show linear fits through the experimental values. Panel (c) shows the SPOD spectrum for the $S=0.61$ case. The broken vertical line corresponds to $f=723~\text{Hz}$, determined from the fit in (a).

We next determine the critical swirl number, $S_{c}$, at which the coherent PVC oscillation emerges as follows. Figures 7(a) and 7(b) show the variation of $f_{PVC}$ and $|u_{z}^{\prime }|^{2}$ at $f_{PVC}$, respectively, with $S$. The data in figure 7(b) have been determined from the component at $f_{PVC}$, of the continuous wavelet transform, $\tilde{u} _{z}$, of the velocity time series at point 1 (see figure 5a). All non-zero data shown in figure 7(a,b) correspond to values of $S$ for which $|C_{12}|\geqslant 0.7$ at their respective $f_{PVC}$ values. Linear least-squares fits through the data, determined using non-zero values for $f_{PVC}$ and $|u_{z}^{\prime }|^{2}$, have been superposed in each case. The upstream movement of the breakdown bubble can induce changes to $|u_{z}^{\prime }|^{2}$ at point 1 because its distance from the region where the PVC oscillation amplitude is maximum, varies as the bubble moves upstream with increasing $S$. However, figure 7(a,b) qualitatively show a nearly linear variation of $f_{PVC}$ and $|u_{z}^{\prime }|^{2}$ with $S$. This is consistent with the results in (2.40) and (2.42) and hence suggest that the jet undergoes a supercritical Hopf bifurcation. Requiring $|\tilde{u} _{z}|^{2}=0$ at $S=S_{c}$ yields, from the linear fits shown in figure 7(a,b), $S_{c}\sim 0.61$ and $f_{PVC}=723~\text{Hz}$ at $S_{c}$. Figure 7 shows the spectrum determined from SPOD of the sPIV measured velocity fields at $S=0.61$. The broken vertical line in figure 7 corresponds to $f=723~\text{Hz}$. Note that this result shows no evidence of a coherent oscillation when compared with the results in figure 3. This result, together with those in figure 3, additionally confirms that the bifurcation in the flow state does indeed occur at $S=0.61$.

It is interesting to note here that Oberleithner et al. (Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011) report additional evidence for the occurrence of supercritical Hopf bifurcation from experiments wherein the PVC is axially forced near the jet exit using an arrangement of loudspeakers. Their results show that the forcing amplitude needed to synchronize the oscillation frequency of the PVC with that of the forcing scales linearly with $|\,f-f_{PVC}|$, where $f$ is the forcing frequency. These results suggest that the underlying PVC oscillation is a stable limit cycle and is therefore produced by a supercritical Hopf bifurcation process, at least in the present class of variable swirl jet experiments.

We confirm the above facts in this paper by determining a series of results from the weakly nonlinear theory presented in § 2. The time-averaged base-flow field, $\boldsymbol{q}_{o}$, at $S=S_{c}$, needed to compute these results is determined from sPIV measurements at $S=0.61$. Further, we adopt a numerical approach to solve the various equations presented in § 2. All of these details are discussed in the next section.

5.1 Linear and leading-order nonlinear solutions

The global linear stability problem given by (2.19), along with (2.20)–(2.22) and (5.1)–(5.2), are solved using the pseudospectral collocation method (Boyd Reference Boyd2000). Accordingly, all dependent variables are written in terms of Chebyshev polynomials of the second kind using function values defined on a two-dimensional mesh of Gauss–Lobatto points in computational space. The parameters in (C 1) and (C 2) are tuned to ensure sufficient collocation point density in physical space in the regions where shear layers of the time-averaged flow field are located and at the domain centreline. The radial and the axial derivatives in (2.19) are replaced by equivalent discrete differentiation matrices (Boyd Reference Boyd2000), which yields a generalized matrix eigenvalue problem as follows:

(5.3)$$\begin{eqnarray}A\hat{\boldsymbol{q}}_{1m,d}=\text{i}\unicode[STIX]{x1D714}_{g,m}B\hat{\boldsymbol{q}}_{1m,d},\end{eqnarray}$$

where $B$ and $A$ are matrices representing the discrete equivalent of continuous operators ${\mathcal{B}}$ and ${\mathcal{L}}_{m}$, respectively. Also, $\unicode[STIX]{x1D714}_{g,m}$ is the complex eigenvalue and $\hat{\boldsymbol{q}}_{1m,d}$ the discrete equivalent of the eigenfunction $\hat{\boldsymbol{q}}_{1m}$ in (2.19). The boundary conditions are written in discrete form and are imposed by replacing the rows corresponding to points on the boundaries in $B$ and $A$. The eigenvalue problem given in (5.3) is solved using the Arnoldi iterative solver ‘eigs’ available in MATLAB.

The solutions for the adjoint mode are obtained using the discrete adjoint approach outlined in Juniper & Pier (Reference Juniper and Pier2015) as follows. First, the inner product, equation (2.30), is discretized using Gauss–Chebyshev quadrature in computational space as $\langle \hat{\boldsymbol{q}}_{1,d},\hat{\boldsymbol{q}}_{2,d}\rangle =\hat{\boldsymbol{q}}_{1,d}^{H}M\hat{\boldsymbol{q}}_{2,d}$, where $M$ is the matrix generated by the weights of the quadrature rule and the superscript ‘$H$’ denotes the transpose conjugate. The matrices of the discrete adjoint problem corresponding to (5.3) are now given by $A^{\dagger }=M^{-1}A^{H}M$ and $B^{\dagger }=M^{-1}B^{H}M$. Thus, the adjoint mode, $\hat{\boldsymbol{q}}_{1m,d}^{\dagger }$, is computed using the ‘eigs’ iterative solver as well.

Figure 8. Variation of $E=|\unicode[STIX]{x1D714}_{g,-1}^{k+1}-\unicode[STIX]{x1D714}_{g,-1}^{k}|$ as a function of number of collocation points using the swirling base flow from computations reported in Pradeep (Reference Pradeep2018) of the $Re=200$ and $S=1$ case discussed in Meliga et al. (Reference Meliga, Gallaire and Chomaz2012). The superscript $k$ indexes computations in increasing order of $N_{I}$.

We validate our linear stability analysis solver against the results reported by Meliga et al. (Reference Meliga, Gallaire and Chomaz2012) for the most unstable $m=-1$ mode of an unconfined swirling jet undergoing a spiral-type vortex breakdown at $S=1.0$ and $Re=200$. The base flow for the validation case is obtained from flow simulations performed by Pradeep (Reference Pradeep2018). These simulations are performed on a on a cuboidal computational domain that extends 15 jet inflow diameters in the streamwise direction and 20 jet inflow diameters in the transverse direction. The domain was discretized using 481 and 241 points in the streamwise and transverse directions respectively. These solutions were obtained using the INCOMPACT3D solver (Laizet & Li Reference Laizet and Li2011). The solution captures unsteady flow features described by Meliga et al. (Reference Meliga, Gallaire and Chomaz2012) for this case (see Pradeep (Reference Pradeep2018) for details).

Linear stability computations are performed on a computational mesh comprising $N$ points each in the radial and axial directions in the computational domain. The quasi-steady flow solutions are interpolated onto the linear stability computational mesh at each axial location using cubic spline interpolation along the radial direction. Figure 8 shows the variation of $E=|\unicode[STIX]{x1D714}_{g,-1}^{k+1}-\unicode[STIX]{x1D714}_{g,-1}^{k}|$ with the number of collocation points in the streamwise direction ($N_{I}$). Note that the vertical axis uses a logarithmic scale. The superscript $k$ indexes the computation in the order of increasing $N_{I}$. The value $\unicode[STIX]{x1D714}_{g,-1}$ is the prediction from our numerical solutions. Note that $E\sim 1\times 10^{-4}$ for $N_{I}=120$. The curve fit through the data, shown in figure 8, shows that $\unicode[STIX]{x1D714}_{g,-1}^{k}$ converges exponentially with grid refinement, as can be expected for the present pseudospectral method. Table 1 shows the values of $\unicode[STIX]{x1D714}_{g,-1}$ for $N_{I}=120$ and the final converged value reported by Meliga et al. (Reference Meliga, Gallaire and Chomaz2012). The value of $N_{DOF}$ in table 1 refers to the total number of degrees of freedom for three-dimensional state vector in each study. The small difference between these results (${\sim}2.5\,\%$) validates the implementation of the linear global stability solver used in the present study.

Table 1. Most unstable $m=-1$ mode eigenvalues for the $Re=200$, $S=1.0$ case reported in Meliga et al. (Reference Meliga, Gallaire and Chomaz2012). Present results are determined using the swirling base flow from computations reported in Pradeep (Reference Pradeep2018). The value of $N_{DOF}$ refers to the total number of degrees of freedom for three-dimensional state vector.

Once the linear stability solutions have been determined, the leading-order nonlinear solutions are determined from (2.28)–(2.26). These are solved by replacing spatial derivatives occurring in the operators on the left with the discrete differentiation matrices and solving the resulting linear system using LU decomposition. The Stuart–Landau coefficients are then determined using (2.34) and (2.33), using the discretized inner product.

We apply the numerical methods described in this section to the time-averaged base-flow state at $S=0.61$ as determined from sPIV measurements in order to gain insight into the fundamental physical mechanism that results in the emergence of PVC oscillations. Note that the fact that the optimal locations for collocation points do not correspond to the centres of the interrogation windows used while processing the sPIV measurements. Accordingly, the base-flow state on the computational mesh is determined using a combination of interpolation and extrapolation procedures as will be described next.

5.2 Base flow

Figure 9(a–c) shows the spatial variation of the time-averaged velocity components, $\bar{U}_{z}$, $\bar{U}_{\unicode[STIX]{x1D703}}$ and $\bar{U}_{r}$ determined using sPIV at $S=0.61$. Three-point moving-average filtering of these time-averaged fields is used to remove small length scale oscillations due to uncertainties in the experimental data. The spatial variations of $\bar{U}_{z}$ and $\bar{U}_{\unicode[STIX]{x1D703}}$ (figure 9a,b) show expected symmetry characteristics across the flow centreline. However, the spatial variation of $\bar{U}_{r}$ (figure 9c) is not perfectly anti-symmetric about $r=0$, as one may expect. This is because the presence of turbulence in the flow along with intermittent vortex breakdown observed for this case, both of which cause the rotational axis of the flow to ‘wobble’ incoherently about the geometric jet centreline; details of this process are described in Clees et al. (Reference Clees, Lewalle, Frederick and O’Connor2018). Ideally, given a long enough time record, one can expect that the impact of these incoherent flow centreline oscillations would average out, yielding an axisymmetric $\bar{U}_{r}$ field (Syred, O’Doherty & Froud Reference Syred, O’Doherty and Froud1994). However, the fact that figure 9(c) shows that this is not the case suggests that the intensity of turbulence fluctuations in the present experiment is comparable to $\bar{U}_{r}$ and that additional flow field samples are necessary to achieve good statistical confidence in the data. We next describe how base-flow velocity fields for stability analysis computations are determined from these sPIV measurements.

Figure 9. Time-averaged velocity fields (a) $\bar{U}_{z}$ (b) $\bar{U}_{\unicode[STIX]{x1D703}}$ (c) $\bar{U}_{r}$ ($S=0.61$, $Re=59\,000$). The broken curves show the nominal location of the inner (black) and outer (magenta) shear layers in the time-averaged flow.

Note that the collocation points used for computations do not necessarily lie at the same locations in physical space as the interrogation window centres used in sPIV processing. Also, the field of view of the cameras used in the experiments limits the availability of sPIV data to a region corresponding to $0\leqslant r<1.5$ around the geometric centreline of the flow. Therefore, we use a combination of interpolation and extrapolation using sPIV data to determine base flow $\bar{U}_{\unicode[STIX]{x1D703}}$ and $\bar{U}_{z}$ fields for $1.5<r\leqslant 50$ as follows. First, least-squares fits to the time-averaged $\bar{U}_{\unicode[STIX]{x1D703}}$ and $\bar{U}_{z}$ sPIV data of model flow profiles for a swirling jet are determined using the model flow profiles for a swirling jet suggested by Oberleithner et al. (Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011). Next, individual cubic spline fits to the radial profiles of $\bar{U}_{\unicode[STIX]{x1D703}}$ and $\bar{U}_{z}$ at each $z$ are determined. A uniform interpolation is then constructed by smoothly merging the cubic spline fits with the fits to model flow profiles at each $z$ in a small region centred at $r=1.5$ using a hyperbolic tangent function. Figure 10(a,b) shows typical radial profiles of $\bar{U}_{\unicode[STIX]{x1D703}}$ and $\bar{U}_{z}$ at $z=1$ determined from the fitting procedure. The time-averaged sPIV data are overlaid on these profiles for comparison. These results are typical and a similar quality of agreement with sPIV data is seen at other $z$ locations as well. The values of $\bar{U}_{\unicode[STIX]{x1D703}}$ and $\bar{U}_{z}$ at the collocation points on the computational mesh are determined using these fits.

Figure 10. Extrapolated time-averaged velocity fields at $z/D=1$ (a) $\bar{U}_{z}$ (b) $\bar{U}_{\unicode[STIX]{x1D703}}$ ($S=0.61$, $Re=59\,000$).

We determine the $\bar{U}_{r}$ field from the measured $\bar{U}_{z}$ field using the fact that the time-averaged flow must be divergence free rather than use the experimentally measured result shown in figure 9(c). This ensures that our stability analysis is anchored around a base flow that satisfies the time-averaged mass balance condition. This constraint on the time-averaged flow field yields,

(5.4)$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}+\frac{1}{r}\right)\bar{U}_{r}=-\frac{\unicode[STIX]{x2202}\bar{U}_{z}}{\unicode[STIX]{x2202}z},\end{eqnarray}$$

along with $\bar{U}_{r}=0$ at $r=0$ (centreline). The spatial derivatives in (5.4) are discretized using the same Chebyshev pseudospectral scheme described in the previous subsection. The resulting linear system of equations for $\bar{U}_{r}=0$ at the collocation points in the domain is solved using an LU decomposition. Figure 11 shows the $\bar{U}_{r}$ field computed thus. For all computations involving $\bar{U}_{r}$, the field shown in figure 11 is used instead of the corresponding time-averaged sPIV field (see figure 9c).

Next, the turbulent viscosity, $\unicode[STIX]{x1D708}_{T}$, field is determined as follows. Reynolds stress components in the region $0<r\leqslant 1.5$ are determined on the sPIV mesh from the sPIV measurements. These values are then interpolated onto the computational mesh using the same interpolation procedure as was used for the time-averaged velocity. The value of each of the Reynolds stresses is allowed to smoothly vanish for $r\geqslant 1.5$ using a hyperbolic tangent function. The $\unicode[STIX]{x1D708}_{T}$ field is then determined using (B 3) at each point. Figure 12 shows the spatial variation of $\unicode[STIX]{x1D708}_{T}$ normalized by the value of the molecular viscosity, $\unicode[STIX]{x1D708}$. As may be expected, $\unicode[STIX]{x1D708}_{T}$ is concentrated between the two shear layers and increases with axial distance downstream. Note that (B 3) represents the solution to an unconstrained least-squares minimization problem. Also, the experimental uncertainty in the measured velocity values, translate into uncertainty in the Reynolds stresses and mean strain-rate fields. Both of these reasons cause $\unicode[STIX]{x1D708}_{T}$ to become negative in some regions of the flow. We find that these values appear predominantly in regions where turbulence intensities are small. Accordingly, we set to $\unicode[STIX]{x1D708}_{T}=0$ at all points where $\unicode[STIX]{x1D708}_{T}<0$ in order to avoid the generation of spurious eigenvalues in the linear stability analysis (Rukes et al. Reference Rukes, Paschereit and Oberleithner2016).

Figure 11. Base-flow field for $\bar{U}_{r}$ determined from time-averaged sPIV profiles for $\bar{U}_{\unicode[STIX]{x1D703}}$ and $\bar{U}_{z}$ using the divergence free constraint ($S=0.61$, $Re=59\,000$). The broken curves show the nominal location of the inner (black) and outer (magenta) shear layers in the time-averaged flow.

Figure 12. Turbulent viscosity, $\unicode[STIX]{x1D708}_{T}$, field determined from sPIV data at $S_{c}=0.61$ ($Re=59\,000$). The broken curves show the nominal location of the inner (black) and outer (magenta) shear layers in the time-averaged flow.

All curve fitting and linear equation solutions discussed in this subsection is performed using functions available as part of MATLAB. The results presented in this paper use meshes that have between 60 and 65 points in the axial direction and 80 and 100 points in the radial direction. These limits represent the finest meshes possible such that there are no more than two collocation points per cell of the sPIV data. For all cases, the number and location of collocation points in physical space were refined until the change in the magnitude of all physically relevant eigenvalues was less than $1\times 10^{-5}$ between successive grid refinements.