1. Introduction
The hydrodynamic near-field and the far-field acoustics of jets are significantly influenced by the inflow conditions and the state of the boundary layer at the nozzle exit. In particular, the flow field is sensitive to parameters such as momentum thickness and fluctuation level, and whether the boundary layer is laminar or turbulent. Several experimental (Hill Jr et al. Reference Hill, Jenkins and Gilbert1976; Hussain & Zedan Reference Hussain and Zedan1978a , Reference Hussain and Zedanb ; Husain & Hussain Reference Husain and Hussain1979; Zaman Reference Zaman1985; Bridges & Hussain Reference Bridges and Hussain1987; Zaman Reference Zaman2012; Fontaine et al. Reference Fontaine, Elliott, Austin and Freund2015) and numerical (Bogey & Bailly Reference Bogey and Bailly2005; Kim & Choi Reference Kim and Choi2009; Bogey & Bailly Reference Bogey and Bailly2010; Bogey, Marsden & Bailly Reference Bogey, Marsden and Bailly2012; Brès et al. Reference Brès, Jordan, Jaunet, Le Rallic, Cavalieri, Towne, Lele, Colonius and Schmidt2018) studies have explored the effect of these parameters. For a jet with a laminar boundary layer at laboratory-scale Reynolds number, the flow emerging from the nozzle mixes with the surrounding fluid and transitions to turbulence within the first few jet diameters. This rapid mixing results in roll-up and pairing of vortices, consequently leading to an increase in the radiated noise (Zaman Reference Zaman1985; Bridges & Hussain Reference Bridges and Hussain1987; Bogey & Bailly Reference Bogey and Bailly2010). Zaman (Reference Zaman1985) showed that an initially laminar jet exhibits a 4 dB increase in the radiated noise compared with its turbulent counterpart. Bogey & Bailly (Reference Bogey and Bailly2010) demonstrated that initially laminar jets with larger momentum thickness exhibit stronger vortex pairing that increases the sound pressure levels in the sideline direction.
Vortex pairing is a main characteristic of mixing layers (Brown & Roshko Reference Brown and Roshko1974; Winant & Browand Reference Winant and Browand1974; Ho & Huang Reference Ho and Huang1982; Metcalfe et al. Reference Metcalfe, Orszag, Brachet, Menon and Riley1987; Moser & Rogers Reference Moser and Rogers1993) and jet flows (Becker & Massaro Reference Becker and Massaro1968; Hussain & Zaman Reference Hussain and Zaman1980; Zaman & Hussain Reference Zaman and Hussain1980; Meynart Reference Meynart1983). It significantly contributes to turbulent mixing (Brown & Roshko Reference Brown and Roshko1974), the production of Reynolds stresses (Zaman & Hussain Reference Zaman and Hussain1980), entrainment (Winant & Browand Reference Winant and Browand1974) and triggers the transition to turbulence (Ho & Huang Reference Ho and Huang1982; Moser & Rogers Reference Moser and Rogers1993). Given its importance, vortex pairing has been the subject of numerous studies. In jets, it was first visualised by Becker & Massaro (Reference Becker and Massaro1968). In the seminal study by Crow & Champagne (Reference Crow and Champagne1971), pairing was found to be more regular at low Reynolds numbers and became increasingly chaotic at higher Reynolds numbers. Detailed experimental studies on forced jets conducted by Zaman & Hussain (Reference Zaman and Hussain1980) and Hussain & Zaman (Reference Hussain and Zaman1980) reveal that pairing occurs at two distinct frequencies: one around
${St}_\theta = f \theta /U \approx 0.012$
, termed the shear-layer mode, and the other around
${St}_D = fD/U \approx 0.85$
, referred to as the jet-column mode. Here,
$f$
is the frequency,
$U$
is the jet velocity,
$D$
is the diameter and
$\theta$
is the momentum thickness. Kibens (Reference Kibens1980) forced a jet at the shear-layer instability frequency
$St_{\theta }$
and found that this forcing results in three successive vortex pairings, producing the subharmonic frequencies
$St_{\theta }/2$
,
$St_{\theta }/4$
and
$St_{\theta }/8$
. Similarly, in mixing layers, Ho & Huang (Reference Ho and Huang1982) observed multiple pairing when the flow was forced at the subharmonic frequency.
Vortex pairing is a nonlinear process that involves the merging of two smaller vortices into a larger vortex with half the frequency. The nonlinear mechanisms behind the growth of the subharmonic during this process have been the focus of many studies (Kelly Reference Kelly1967; Monkewitz Reference Monkewitz1988; Paschereit, Wygnanski & Fiedler Reference Paschereit, Wygnanski and Fiedler1995; Husain & Hussain Reference Husain and Hussain1995). Kelly (Reference Kelly1967) proposed a resonance mechanism that provides energy to the subharmonics. Monkewitz (Reference Monkewitz1988) used weakly nonlinear spatial theory to further support Kelly’s (Reference Kelly1967) resonance mechanism. It was demonstrated that the fundamental mode needs to reach a critical amplitude before both the fundamental and subharmonic can phase lock, resulting in an energy transfer to the subharmonic. Cohen & Wygnanski (Reference Cohen and Wygnanski1987) demonstrated that fundamental–subharmonic interactions occur only when both waves propagate at the same phase velocity, allowing sufficient time for the transfer of energy between the two waves. Experiments by Hajj, Miksad & Powers (Reference Hajj, Miksad and Powers1992) show that the subharmonic gains energy through this resonance mechanism. Works by Husain & Hussain (Reference Husain and Hussain1995) and Cho, Yoo & Choi (Reference Cho, Yoo and Choi1998) demonstrate that vortex pairing can be either enhanced or attenuated by controlling the phase difference between the fundamental and the subharmonic. Mankbadi (Reference Mankbadi1985) employed an energy-integral method to show that the subharmonic gains energy from both the fundamental and the mean flow. Paschereit et al. (Reference Paschereit, Wygnanski and Fiedler1995) estimated the energy transfer to the subharmonic wave based on the production terms. Their findings reveal that the subharmonic primarily gains its energy from the mean flow, with the fundamental wave acting as a catalyst. This catalytic role refers to the scenario where the fundamental promotes subharmonic growth while remaining largely unaffected itself. This concept is consistent with Monkewitz’s (Reference Monkewitz1988) theoretical framework and is, in fact, a common characteristic of resonant wave interactions in shear flows (see e.g. Wu, Stewart & Cowley Reference Wu, Stewart and Cowley2007).
Using classical linear theory and spectral modal analysis, this work answers the following questions: How does the growth of the shear layer differ between an initially laminar jet and a turbulent jet? What is the influence of vortex pairing on the development of the shear layer? What should be considered as the fundamental frequency in a multi-tonal flow? In particular, our study sheds light on the ambiguity of energetic versus dynamical significance in this context. We further demonstrate how the energy transfer during vortex pairing can be quantified.
Linear stability theory (LST) has been used in the past with great success, first by Michalke (Reference Michalke1964, Reference Michalke1965) on a hyperbolic tangent profile. Here, we apply LST to the temporally averaged mean flow of a jet. This approach corresponds to applying the parallel flow assumption locally to the zero-frequency component, which has been utilised in cylinder wakes (Pier Reference Pier2002; Barkley Reference Barkley2006) and jets (Suzuki & Colonius Reference Suzuki and Colonius2006; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Schmidt et al. Reference Schmidt, Towne, Colonius, Cavalieri, Jordan and Brès2017). Recent studies (Wu & Zhuang Reference Wu and Zhuang2016; Zhang & Wu Reference Zhang and Wu2020, Reference Zhang and Wu2022) have also examined the effects of non-parallelism and nonlinearity, finding both to be important. The mathematical framework of spectral proper orthogonal decomposition (SPOD) dates back to the work of Lumley (Reference Lumley1967, Reference Lumley1970). SPOD identifies the most energetic coherent structures at each time scale. Early applications of SPOD include the work of Glauser, Leib & George (Reference Glauser, Leib and George1987), Glauser & George (Reference Glauser and George1992) and Delville (Reference Delville1994). More recently, this method has attracted significant interest following its application to large flow databases by Schmidt et al. (Reference Schmidt, Towne, Rigas, Colonius and Brès2018) and the establishment of its relationship to other methods by Towne, Schmidt & Colonius (Reference Towne, Schmidt and Colonius2018). Since then, SPOD has become a mainstay of physical exploration, in particular for identifying different modal and non-modal instabilities (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Nogueira et al. Reference Nogueira, Cavalieri, Jordan and Jaunet2019; Pickering et al. Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020). While SPOD does not provide a priori insights, it helps reveal the underlying physics of turbulent flows in a statistical sense, particularly when combined with nonlinear analysis. In contrast, LST offers predictive capabilities, providing a priori insights. Therefore, in this study, we employ both LST and SPOD to characterise the vortex-pairing process.
Vortex pairing is a nonlinear process involving energy transfer between different frequencies or scales. This interscale energy transfer arises from the quadratic nonlinearity of the Navier–Stokes equations, leading to triadic interactions. Various approaches have been employed to analyse the energy transfer between different scales, including the Karman–Howarth equation (Von Karman & Howarth Reference Von Karman and Howarth1938; Danaila et al. Reference Danaila, Anselmet, Zhou and Antonia1999; Hill Reference Hill2001) and bispectrum analysis (Lii, Rosenblatt & Van Atta Reference Lii, Rosenblatt and Van Atta1976; Kim & Powers Reference Kim and Powers1979; Herring Reference Herring1980). The former is based on structure functions, and the latter is the frequency-domain representation of third-order moments. Recently, a modal decomposition based on bispectral analysis was developed by Schmidt (Reference Schmidt2020). Consistent with our data analysis approach, we quantify the production, dissipation and nonlinear transfer between the leading SPOD modes associated with different harmonics and the mean flow, based on the spectral turbulent kinetic energy (TKE) equation. Previous studies (Mizuno Reference Mizuno2016; Cho, Hwang & Choi Reference Cho, Hwang and Choi2018; Lee & Moser Reference Lee and Moser2019; Gomé et al. Reference Gomé, Tuckerman and Barkley2023) have investigated interscale energy transfer in turbulent channel flows using the spectral TKE equation. Additionally, researchers have employed various bases, such as resolvent modes (Symon, Illingworth & Marusic Reference Symon, Illingworth and Marusic2021; Jin, Symon & Illingworth Reference Jin, Symon and Illingworth2021), Fourier modes (Nekkanti et al. Reference Nekkanti, Nidhan, Schmidt and Sarkar2023a ), dynamic mode decomposition modes (Kinjangi & Foti Reference Kinjangi and Foti2023) and optimal mode decomposition modes (Biswas, Cicolin & Buxton Reference Biswas, Cicolin and Buxton2022), to estimate the spectral energy budget.
The paper is organised as follows. In § 2, the methodologies of LST, SPOD and SPOD-based spectral energy budget analysis are discussed. Results focusing on shear-layer instability, vortex pairing and spectral energy transfer are presented in § 3. The paper concludes with discussions and conclusions in § 4.
2. Methodology
2.1. Linear stability theory
We employ the spatial form of LST. This LST determines the spatial growth rates of the corresponding frequency. Here, we use it to identify the most unstable frequency. The LST is also referred to as local linear theory, and throughout this paper, we will use the terms LST and local linear theory interchangeably. In LST, the frequency
$\omega$
is assumed to be real, and the eigenvalue problem is solved for a complex
$\alpha$
. The real part of
$\alpha$
is the streamwise wavenumber and the imaginary part is its amplification rate. We start off with the Reynolds decomposition,
${q} (x,r,\theta ,t )=\bar {{q}} (x,r )+{q}^{\prime } (x,r,\theta ,t )$
, where
$\bar {{q}}$
is the mean flow and
${q}^{\prime }$
is the fluctuation. The basic assumption of the local linear theory that the flow is locally streamwise parallel, i.e. the flow is homogeneous in the streamwise direction. Using the normal mode ansatz, the fluctuation is expressed as
\begin{equation} {q}^{\prime }(x,r,\theta ,t) = \tilde {{q}}(r) {\textrm e}^{{\textrm i} \left (\alpha x + m \theta - \omega t\right )}, \end{equation}
where the streamwise wavenumber
$\alpha$
is complex and
$\tilde {{q}}$
is the radial profile. Linearising the Navier–Stokes equation about the base flow yields the equation
\begin{equation} \left ({\textrm i}\omega {I} + {L}\right )\tilde {{q}}= 0, \end{equation}
where
${I}$
is the identity matrix and
${L}$
is the linearised compressible Navier–Stokes operator. The domain is extended to the far field by mapping the original domain
$r \in [0,6]$
to
$r \in [0, \infty )$
by the mapping function suggested by Lesshafft & Huerre (Reference Lesshafft and Huerre2007). Dirichlet boundary conditions are used for a far-field boundary, and the radial direction is discretised using Chebyshev collocation points. Finally, the eigenvalue problem is solved using the methodology of Maia et al. (Reference Maia, Jordan, Cavalieri, Martini, Sasaki and Silvestre2021, Reference Maia, Jordan and Cavalieri2022). Solving this eigenvalue problem yields monochromatic amplification rates,
$\alpha _i$
. For
$\alpha _i \lt 0$
, the disturbances will grow exponentially downstream, whereas for
$\alpha _i \gt 0$
, they will decay downstream.
2.2. Spectral proper orthogonal decomposition
SPOD decomposes stationary flow data into monochromatic modes that represent the spatial flow structures optimised in terms of the flow’s energy. The eigenvalues corresponding to these modes represent their energy. We employ a SPOD algorithm based on Welch’s method (Welch Reference Welch1967). For the mathematical derivation and computational details, we refer the reader to Towne et al. (Reference Towne, Schmidt and Colonius2018) and Schmidt & Colonius (Reference Schmidt and Colonius2020). The SPOD modes and eigenvalues are obtained by solving the eigenvalue problem
\begin{equation} \int _{\varOmega } \boldsymbol{S}\left ({x},{x}',f\right ){W}\left ({x}'\right )\boldsymbol{\phi }\left ({x}',f\right ){\textrm{d}}x' = \boldsymbol{\phi }\left ({x},f\right ) \boldsymbol{\lambda }(f), \end{equation}
where
$\boldsymbol{S}$
is the cross-spectral density matrix,
${W}$
the positive-definite matrix that accounts for component-wise and numerical quadrature weights,
$\boldsymbol{\phi }$
the SPOD modes and
$\boldsymbol{\lambda }$
the eigenvalues. At each frequency, the eigenvalue decomposition of the cross-spectral density matrix
$\boldsymbol{S}$
yields a countably infinite number of SPOD modes,
$\boldsymbol{\phi }^{(i)}({x},f)$
, and corresponding eigenvalues,
$\lambda ^{(i)}(f)$
, where
$i$
is the mode number index. The modes are sorted by energy:
$\lambda ^{(1)}(f) \geqslant \lambda ^{(2)}(f) \geqslant \dots \geqslant 0$
. At each frequency, the SPOD modes are orthogonal in space:
\begin{equation} \int _{\varOmega } \boldsymbol{\phi }^{(i)}({x},f) {W}\left ({x}\right )\boldsymbol{\phi }^{(j)}({x},f) {\textrm{d}}x = \delta _{ij}, \end{equation}
where
$\delta _{ij}$
is the Kronecker delta function. Recently, Nekkanti & Schmidt (Reference Nekkanti and Schmidt2021) proposed a convolution approach that computes time-continuous expansion coefficients:
\begin{equation} a^{(i)}(f,t) = \int _{\Delta T} \int _{\varOmega } {\left ( {\phi }^{(i)}({x},f)\right )}^{*}{W}({x}){\textbf{q}}({x},t+\tau ) {\textrm e}^{-{\textrm i} 2\pi f \tau }{\textrm{d}} x\, {\textrm{d}} \tau , \end{equation}
which facilitates obtaining the time-continuous Fourier modes:
\begin{equation} {\hat {q}}({x},f,t) \approx \sum \limits _{i} a^{(i)} (f, t) {\phi }^{(i)} ({x},f). \end{equation}
2.3. SPOD-based spectral energy budget
The optimality of the SPOD expansion, (2.6), can be leveraged to quantify the nonlinear interactions between the most salient flow features and their global energy budget. We specifically focus on production and nonlinear energy transfer. The starting point is the spectral TKE equation:
\begin{equation} \frac {\partial \hat {k}}{\partial t} \!=\! {\mathcal{R}} \! \left [ -\bar {u}_i \frac {\partial \hat {k}}{\partial x_i} \underbrace {\vphantom {\displaystyle \sum _{i=1}^N}- \overline {\hat {u}_j^*\widehat { u_i\frac {\partial u_j}{\partial x_i}}}}_{\mathcal{T}_{nl}} \underbrace {\vphantom {\displaystyle \sum _{i=1}^N}-\overline {\hat {u}_j^*\hat {u}_i\frac {\partial \bar {u}_j}{\partial x_i}}}_{\mathcal{P}} \underbrace {\vphantom {\displaystyle \sum _{i=1}^N} \!-\! \frac {2}{{Re}} \overline {\hat {s}_{i j}^* \hat {s}_{i j}}}_{ \mathcal{D}} -\frac {\partial }{\partial x_j}\big (\overline {\hat {u}_j^*\hat {p}}\big ) \!+\! \frac {2}{{Re}} \frac {\partial }{\partial x_i} \big (\overline {\hat {u}_j^*\hat {s}_{ij}}\big ) \! \right ], \end{equation}
where
$\hat {s}_{ij} = 1/2 (\partial \hat {u}_i/\partial x_j + \partial \hat {u}_j/\partial x_i )$
is the spectral strain rate,
${\mathcal{R}}$
denotes the real part and
$\hat {k}(f,m)= \overline {\hat {u}_i^{*}(f,m)\hat {u}_i(f,m)}/2$
. Readers are referred to Bolotnov et al. (Reference Bolotnov, Lahey, Drew, Jansen and Oberai2010) and Cho et al. (Reference Cho, Hwang and Choi2018) for the derivation of the spectral TKE equation, with a detailed derivation provided in Appendix B for completeness. For brevity, we suppress the spatial dependence of all quantities. For statistically stationary flows, the left-hand side goes to zero. In this work, we focus on the spectral production term
\begin{equation} {\mathcal{P}}\left (f , m\right ) = -{\mathcal{R}}\overline {\left [\hat {u}_j^*\left (f,m\right )\hat {u}_i\left (f,m\right )\frac {\partial \bar {u}_j}{\partial x_i}\right ]}, \end{equation}
the spectral dissipation term
\begin{equation} \mathcal{D}\left (f , m\right ) = - \frac {2}{{Re}}{\mathcal{R}}\overline {\left [ \left (\frac {\partial \hat {u}_i^*}{\partial x_j}\left (f,m\right ) + \frac {\partial \hat {u}_j^*}{\partial x_i}\left (f,m\right )\right ) \left (\frac {\partial \hat {u}_i}{\partial x_j}\left (f,m\right ) + \frac {\partial \hat {u}_j}{\partial x_i}\left (f,m\right )\right )\right ]} \end{equation}
and the spectral nonlinear transfer term
\begin{align} \mathcal{T}_{{nl}}\left (f, m\right ) = - {\mathcal{R}}\overline {\left [\hat {u}_j^*\left (f , m\right )\widehat {u_i\frac {\partial u_j}{\partial x_i}}\left (f, m\right )\right ]}. \\[-18pt] \nonumber \end{align}
The term 
in (2.10) entails the contributions from all other frequencies and azimuthal wavenumbers to frequency
$f$
and azimuthal wavenumber
$m$
. Following Cho et al. (Reference Cho, Hwang and Choi2018), we isolate the energy transfer of individual triads, i.e. frequency and wavenumber triplets related by the resonance conditions,
$m_1 \pm m_2 \pm m_3=0$
and
$f_1 \pm f_2 \pm f_3=0$
, by splitting 
using the discrete convolution to obtain
\begin{equation} \mathcal{T}_{{nl}}\left (f_3 , m_3\right ) = - {\mathcal{R}}\left [\overline {\hat {u}_j^*\left (f_3,m_3\right )\sum \limits _{\substack {f_1 +f_2=f_3 \\ m_1 +m_2=m_3}}\hat {u}_i\left (f_1,m_1\right )\frac {\partial \hat {u}_j}{\partial x_i}\left (f_2,m_2\right )}\right ]. \end{equation}
To obtain a bispectral representation of the energy transfer, we further split the sum in (2.11) into individual frequency and wavenumber components that are triadically compatible:
\begin{equation} {\mathscr{t}}_{{nl}}\left (f_1,f_2 , m_1,m_2\right ) = - {\mathcal{R}}\left [\overline {\hat {u}_j^*\left (f_1+f_2,m_1+m_2\right )\hat {u}_i\left (f_1,m_1\right )\frac {\partial \hat {u}_j}{\partial x_i}\left (f_2,m_2\right )}\right ], \end{equation}
such that
$\mathcal{T}_{{nl}}(f_3,m_3)= \sum _{\substack {f_1 +f_2=f_3 \\ m_1 +m_2=m_3}} {\mathscr{t}}_{{nl}} (f_1,f_2 , m_1,m_2 )$
.
At each frequency, the dominant flow structure is represented by the leading SPOD mode. To characterise the production and nonlinear kinetic energy transfer of these structures, we use a rank-1 approximation of the velocity field, following (2.6). As the focus is on self-interactions of the axisymmetric component, we further suppress the azimuthal wavenumber and SPOD mode number dependence with the understanding that
$m_1=m_2=m_3=0$
and
$i=1$
. We have confirmed that higher azimuthal wavenumber components with
$m\gt 0$
do not play a significant role in the dynamics of interest. Substituting (2.6) into (2.8), (2.9), (2.11) and (2.12) we get
\begin{align} {\mathcal{P}}^{{\textit{rank}\text{-}1}}\left (f_1\right ) &= -\overline {a^*\left (f_1\right ){a}\left (f_1\right )} \phi _j^*\left (f_1\right ) \phi _i\left (f_1\right ) \frac {\partial \bar {u_j}}{\partial x_i} = -\lambda (f_1) \phi _j^*\left (f_1\right ) \phi _i\left (f_1\right ) \frac {\partial \bar {u}_j}{\partial x_i}, \end{align}
\begin{align} {\mathcal{D}}^{{\textit{rank}\text{-}1}}\left (f_1\right ) &= -\lambda (f_1) \left (\frac {\partial \phi _i^*}{\partial x_j}\left (f_1\right ) + \frac {\partial \phi _j^*}{\partial x_i}\left (f_1\right )\right ) \left (\frac {\partial \phi _i}{\partial x_j}\left (f_1\right ) + \frac {\partial \phi _j}{\partial x_i}\left (f_1\right )\right ), \end{align}
\begin{align} {\mathscr{t}}_{{nl}}^{{\textit{rank}\text{-}1}}\left (f_1,f_2\right ) &= - {\mathcal{R}}\left (\overline {a^*\left (f_1+f_2\right )a\left (f_1\right )a\left (f_2\right )} \phi _j^*\left (f_1+f_2\right )\phi _i\left (f_1\right )\frac {\partial \phi _j}{\partial x_i}\left (f_2\right )\right ), \end{align}
\begin{align} \mathcal{T}_{{nl}}^{{\textit{rank}\text{-}1}}(f_3) &= - {\mathcal{R}}\left (\sum \limits _{f_1 +f_2=f_3}\overline {a^*\left (f_3\right )a\left (f_1\right )a\left (f_2\right )} \phi _j^*\left (f_3\right )\phi _i\left (f_1\right )\frac {\partial \phi _j}{\partial x_i}\left (f_2\right )\right ). \end{align}
The term
${\mathcal{P}}^{{\textit{rank}\text{-}1}}(f_1)$
signifies energy transfer from frequency
$f_1$
to the mean flow,
${\mathcal{D}}^{{\textit{rank}\text{-}1}}(f_1)$
represents energy dissipation at
$f_1$
,
${\mathscr{t}}_{{nl}}^{{\textit{rank}\text{-}1}}(f_1, f_2)$
denotes nonlinear interactions between
$f_1$
and
$f_2$
and
$\mathcal{T}_{{nl}}^{{\textit{rank}\text{-}1}}(f_3)$
indicates net nonlinear energy transfer into
$f_3$
. In § 3.4, we use these equations to shed light on the energy transfer during the vortex-pairing process.
3. Results
This study investigates the nonlinear dynamics of the initially laminar shear layer of jets at laboratory-scale Reynolds numbers. To validate against experiments and contrast to jets with different shear layers, two large-eddy simulations (LES) of subsonic jets are conducted at a Reynolds number of
${Re} = U_j D/\nu = 50\,000$
. Here,
$U_j$
is the jet exit velocity,
$D$
is the diameter and
$\nu$
is the kinematic viscosity. In both simulations, the jet plume is turbulent. The main difference between the two cases is the state of the boundary layer inside the nozzle and, consequently, the initial free shear layer over the first few jet diameters downstream of the nozzle. In the first case, the boundary layer inside the nozzle is laminar but quickly transitions within the first jet diameter. In the second case, the boundary layer is tripped inside the nozzle and, hence, turbulent from the start. We refer to the first case as the initially laminar jet and the second case as the turbulent jet. Independent of the boundary-layer tripping, both jets exhibit a potential core length that extends over several jet diameters.
The LES are performed using the compressible flow solver ‘Charles’ developed at Cadence, formerly Cascade Technologies (Brès et al. Reference Brès, Ham, Nichols and Lele2017, Reference Brès, Jordan, Jaunet, Le Rallic, Cavalieri, Towne, Lele, Colonius and Schmidt2018), and the reader is referred to Brès et al. (Reference Brès, Ham, Nichols and Lele2017, Reference Brès, Jordan, Jaunet, Le Rallic, Cavalieri, Towne, Lele, Colonius and Schmidt2018) for further details on the numerical method and validation on jet flows. To ensure the accuracy of our simulations, we first validate them by comparing them against companion experiments conducted by Maia, Jordan & Cavalieri (Reference Maia, Jordan and Cavalieri2022). This experimental nozzle geometry is meshed using the same strategy as that of Brès et al. (Reference Brès, Jordan, Jaunet, Le Rallic, Cavalieri, Towne, Lele, Colonius and Schmidt2018), resulting in a total grid size of 16.6 million control volumes. The LES are performed at the experimental Reynolds number. The Mach number is artificially increased to
$M_j=0.4$
to avoid the very small explicit time steps associated with the incompressible limit. It has been confirmed by comparison with the experimental data that the effects of compressibility at this relatively small Mach number are negligible for the purpose of this study (see figure 2).
3.1. Turbulent jets
We compare the initially laminar and turbulent jets using instantaneous visualisations and root mean square (RMS) of the fluctuating streamwise velocity in figure 1. The 95 % (white solid) and 5 % (white dashed) contour lines of the mean streamwise velocity outline the potential core and the jet plume, respectively. We refer to the initial shear layer as the free shear layer between the potential core and the ambient free stream, extending over about five jet diameters downstream of the nozzle exit. The instantaneous streamwise velocity fluctuating fields shown in figure 1(a,b) reveal that the initially laminar jet exhibits a shorter potential core and larger jet width. The initial shear layer over the first two jet diameters is highlighted in the insets of figures 1(a) and 1(b). The hallmark of the untripped jet is the initially laminar shear layer that is easily distinguished from the turbulent shear layer by its spreading rate; the laminar shear layer exiting the nozzle has a significantly lower spreading rate until
$x\approx 0.8$
, where its spreading rate rapidly increases after it transitions to turbulence. This distinction between the laminar and turbulent portions of the shear layer is quantitatively confirmed by the RMS profiles in figure 1(c–g). At the nozzle exit in figure 1(c), the RMS profile of the turbulent jet peaks at
$\gtrsim 10\,\%$
of the free-stream velocity close to the nozzle wall. On the other hand, expectedly, the laminar shear layer of the untripped jet has zero RMS. By
$x=1$
, the RMS profile of the initially laminar jet already resembles that of the turbulent jet. Further downstream, at
$x=5$
, the RMS of the initially laminar jet surpasses that of the turbulent jet, and eventually, in the region of self-similarity, the RMS profiles become nearly identical again. We later show that the rapid growth of the shear layer is associated with exponential growth of the hydrodynamic instability modes supported by the initially laminar shear layer. We explore the nonlinear dynamics of the jet that leads this rapid growth and explain previous observations (Zaman & Hussain Reference Zaman and Hussain1980; Kim & Choi Reference Kim and Choi2009; Bogey & Bailly Reference Bogey and Bailly2010).
Comparison of the initially laminar and turbulent jets: instantaneous fluctuating streamwise velocity field of (a) initially laminar jet and (b) turbulent jet. The RMS of streamwise velocity at (c)
$x=0$
; (d)
$x=1$
; (e)
$x=2$
; (f)
$x=5$
; (g)
$x=15$
. The potential core and the jet width are indicated as lines of constant
$u_x$
at 95 % and 10 % of the jet velocity
$U_j$
, respectively.
Experimental validation of turbulent jet LES and comparison of initially laminar jet LES with the literature (Zaman & Hussain Reference Zaman and Hussain1980; Bogey & Bailly Reference Bogey and Bailly2010): (a) mean and (b) RMS of the streamwise velocity on the centreline. The intersection of the black dashed line at
$\overline {u_x}/U_j=0.95$
with the mean streamwise velocity defines the length of the potential core.
Figure 2 shows a four-way comparison between initially laminar jet, turbulent jet, LES and experiments. The mean and RMS streamwise velocities on the centreline are shown in figures 2(a) and 2(b), respectively. In the absence of the initially laminar jet for the considered nozzle geometry, the turbulent jet is used to validate the numerical setup. Good agreement is observed for the turbulent jet in terms of the mean and RMS streamwise velocities, but the LES of the initially laminar jet shows significant differences. Along the centreline, the mean flow of the initially laminar jet exhibits a dip at
$x \approx 2$
(see the inset of figure 2
a) and decays rapidly beyond the end of the potential core. The mean flow velocity profile of Bogey & Bailly (Reference Bogey and Bailly2010) exhibits the same phenomenon. The RMS of the initially laminar jet in figure 2(b) is notably higher than the fully turbulent case and exhibits a distinct hump at
$x \approx 2.8$
. This hump, to a certain degree, was also observed by Zaman & Hussain (Reference Zaman and Hussain1980) and Bogey & Bailly (Reference Bogey and Bailly2010). This hump and elevated RMS were previously associated with vortex pairing by Kim & Choi (Reference Kim and Choi2009), Bogey & Bailly (Reference Bogey and Bailly2010) and Bogey et al. (Reference Bogey, Marsden and Bailly2012). The authors also demonstrated that both these phenomena strongly depend on the initial shear-layer thickness, with a thicker shear layer growing faster and exhibiting enhanced vortex pairing. In contrast, a thinner shear layer will grow earlier but slower, resulting in lower
$u_x^{rms}$
on the centreline (Bogey & Bailly Reference Bogey and Bailly2010). In the remainder of this paper, we go beyond this phenomenological description and analyse the underlying nonlinear mechanism in detail. To this end, we use spectral modal decomposition techniques and local LST.
3.2. Shear-layer instability
Premultiplied radially integrated PSD,
${\int }_r {St}\boldsymbol{\cdot }{PSD}r\,{\textrm d}r$
, along
$x$
for initially laminar jet (a,b) and turbulent jet (c,d). (a,c) Streamwise velocity,
$u_x$
; (b,d) radial velocity,
$u_r$
. Dashed lines indicate the three tones of the initially laminar jet, and the dotted line corresponds to the most energetic frequency of the turbulent jet.
Figure 3 shows the premultiplied power spectral density (PSD) of the streamwise and radial velocities for the initially laminar and turbulent jets, respectively. The PSD is computed at each spatial location and integrated radially for each streamwise location. For the initially laminar jet, the PSDs of
$u_x$
and
$u_r$
show the prominence of three frequency components:
${St} = 1.76, 0.88$
and 0.44. The three peaks are located at (
$St$
,
$x) =(1.76, 1.0)$
, (0.88, 1.45) and (0.44, 2.4). The most upstream peak occurs at
$St = 1.76$
as a result of the convective Kelvin–Helmholtz instability of the initial shear layer. Notably, it is an integer multiple of the second and third peaks. This observation suggests that
$St = 1.76$
is the fundamental frequency. We later confirm this using LST and demonstrate in § 3.3 that the second and third peaks arise from vortex pairing, i.e. due to nonlinearity as opposed to hydrodynamic instability. The fundamental frequency translates to
${ St}_{\theta } = f\theta /U_j = 0.0132$
, which closely matches that of the shear-layer mode observed in forced-jet experiments (Zaman & Hussain Reference Zaman and Hussain1980). For all three frequencies, the PSD of the
$u_r$
component is greater than that of the
$u_x$
component. The PSDs for the turbulent jet in figure 3(c,d) are significantly lower than those for the initially laminar jet. For the turbulent jet, the first tonal component is observed at
$ x \approx 2.6$
and
${St}=0.56$
(exact frequency determined from figure 6), denoted by the dotted line. No peaks are observed at
$St=1.76$
and
$0.88$
, indicating that the dynamics associated with these frequencies is absent in the fully turbulent jet. Bridges & Hussain (Reference Bridges and Hussain1987) demonstrated that tripping the boundary layer suppresses the development of shear-layer instabilities and eliminates the associated vortex pairing.
Total RMS velocities and selected contributing frequency components along
$x$
: (a)
$u_x^{rms}$
on the centreline,
$r=0$
; (b)
$u_x^{rms}$
on the lipline,
$r=0.5$
; (c)
$u_r^{rms}$
on the lipline. The total RMS (black curve) is plotted on the left ordinate, and the remaining curves are plotted on the right ordinate. The red curve is the sum of five frequencies
${St}=1.76$
, 0.88, 0.44, 0.22 and 0.11.
Figure 4 shows the contribution of tonal frequencies to the total RMS of the initially laminar jet. This contribution is shown for the streamwise velocity on the centreline (
$r=0$
) in figure 4(a), the streamwise velocity on the lipline (
$r=0.5$
) in figure 4(b) and the radial velocity on the lipline in figure 4(
$c$
). The total RMS is plotted on the left ordinate, whereas the RMS of five frequencies and their sum (denoted by a red line) are plotted on the right ordinate. In all cases, the frequency
$St=0.44$
exhibits the maximum RMS. The lower frequencies are more significant on the centreline, whereas the higher frequencies are dominant on the lipline. In figure 4(a), the curve representing the sum of the five frequencies, while of a lower magnitude, closely resembles the shape of the total RMS. Upon comparing the sum curve with the individual frequencies, it becomes evident that the frequencies
$St=0.44$
(magenta) and
$St=0.22$
(cyan) correspond to the accumulation of RMS at
$x \approx 2.8$
and
$x \approx 6$
, respectively.
Along the lipline, the total RMS of streamwise and radial velocity peaks at
$x \approx 1.8$
. The curve representing the sum of five frequencies matches the shape of the total curve only up to
$x \lessapprox 1.8$
. This is because the flow downstream exhibits a broadband-like nature, necessitating more frequencies to capture the fluctuating dynamics. The higher frequencies,
$St=1.76,0.88$
and 0.44, are more prominent on the lipline than on the centreline. In figure 4(c), these three frequency components peak successively as the preceding frequency starts to decline. This observation is consistent with the expectation that lower-frequency components peak at more downstream locations where the shear layer is thicker. Our findings indicate energy transfer between the fundamental and subharmonic modes. Later, in § 3.4, we investigate this transfer using the spectral kinetic energy equation.
Next, we perform local linear stability analysis to understand the origin of the fundamental frequency. In particular, we seek a quantitative comparison between the empirical growth rates deduced from data and theory. We define the local amplitude of
$\hat {{q}}$
as
\begin{equation} A(x,m,\omega ) = \sqrt {\int _r \hat {{q}}^*\hat {{q}} r {\textrm{d}} r } = \sqrt {k^{\textit{LES}}}, \end{equation}
where
$k^{\textit{LES}}$
is the turbulent kinetic energy computed from the data. Using (3.1) and the normal mode ansatz (2.1) allows us to compute the empirical growth rate as
\begin{equation} \alpha _i^{{emp}} = - \frac {1}{A}\frac {{\textrm{d}} A}{{\textrm{d}} x}. \end{equation}
The TKE can also be predicted from the local linear theory by integrating and squaring the theoretical amplification rate:
\begin{equation} k^{{LST}} = \textrm {exp}\left \{2\int ^x\alpha (x) {\textrm{d}} x\right \} \end{equation}
In the absence of an amplitude in local linear theory, only a qualitative comparison of the TKE distribution can be made. Here, we choose to normalise
$k^{{LST}}$
by the maximum
$k^{\textit{LES}}$
of each frequency.
Comparison of TKE and amplification rate predicted from LST with empirical data: TKE of the dominant frequencies for (a) initially laminar jet and (b) turbulent jet. Amplification rate predicted from LST (c,d) and empirically from data (e, f), using (3.2), for the initially laminar jet (c,e) and turbulent jet (d,f). The solid and dash-dotted lines in (a,b) represent the TKE computed from data and LST using (3.3), respectively. The green line in (c) denotes the most unstable frequency at each streamwise location. The neutral stability curve is represented by the black line in (c–f). The rightward-pointing green triangle in (d) denotes the most unstable frequency,
$St= 0.9$
, for the turbulent jet.
Figure 5 shows the TKE (figure 5
a,b) and spatial growth rate (figure 5
c–f) for the axisymmetric component
$m=0$
in both jets. The spatial growth rate is computed from local linear theory in figure 5(c,d) and empirically using (3.2) in figure 5(e,f). The black line represents the neutral stability curve, while the green line denotes the frequency associated with the largest growth rate at each streamwise location. The TKE for
${St}=1.76$
, 0.88 and 0.44 in the initially laminar jet and for
${St}=0.56$
in the turbulent jet is shown in figures 5(a) and 5(b), respectively. The TKE estimated from LST (
$k^{{LST}}$
), denoted by the dotted-dashed lines, is also shown, demonstrating a good fit with the data. It is observed that the TKE peak of each frequency is approximately located on the neutral stability curve. This observation is in agreement with the fact that the local growth rate at the maximum is zero. As a visual aid, dashed lines indicating the locations of the maximum TKE extend into figures 5(c) and 5(d), which show the spatial growth rates of the initially laminar and turbulent jets, respectively. Figure 5(c) reveals that
${St} \approx 1.76$
is the most unstable frequency at the nozzle’s exit for the initially laminar jet. The coalescence of the observations that
${St}=1.76$
is the most unstable frequency and occurs most upstream at the nozzle’s exit confirms our previous interpretation of
${St} \approx 1.76$
as the fundamental frequency. The absence of a peak at
${St}=0.88$
and 0.44 indicates that these frequencies are not a result of purely linear hydrodynamic instability. Despite being spatially unstable in their own right, these tonal frequencies arise from the nonlinearity involving the fundamental frequency. On the other hand, for the turbulent jet (figure 5
d), the growth rates are much lower in magnitude, and the most unstable frequency is
${St} \approx 0.90$
. The differences in the most amplified frequencies for the laminar and turbulent jet are in accordance with the findings of Morris & Foss (Reference Morris and Foss2003). The empirical spatial growth rates in figure 5(e,f) are similar to those estimated by local linear theory. In particular, the most unstable frequencies
$St\approx 1.76$
and 0.9, for the initially laminar and turbulent jet, respectively, are well estimated. The notable similarity between the empirical and theoretical results reveals two key points: first, it confirms the validity of the assumption of local linear theory; second, it supports the physical interpretation that the observed amplitude disturbances of individual frequency components represent hydrodynamic instability waves. Note that the TKE predicted by LST deviates from the LES data beyond the neutral point, where nonlinear effects become significant. Beyond this location, nonlinear stability theory can potentially offer a more accurate estimate of the LES data (Zhang & Wu Reference Zhang and Wu2020, Reference Zhang and Wu2022). Nevertheless, the good agreement between the growth rates estimated from LST of the mean flow and those empirically derived from data up to the neutral point aligns with the findings of Pier (Reference Pier2002) and Barkley (Reference Barkley2006). They demonstrated that for flow past a cylinder, LST of the mean flow can accurately predict nonlinear vortex shedding. The LST of the time-averaged mean flow indirectly incorporates nonlinear interactions through mean flow deformations, which arises from the self-interaction of each mode with its conjugate. In weakly nonlinear regimes, the dominant contribution is at quadratic order, such as the self-difference interaction of the fundamental frequency that contributes to mean flow distortion.
Recent works, such as the Floquet analysis of Shaabani-Ardali, Sipp & Lesshafft (Reference Shaabani-Ardali, Sipp and Lesshafft2019) and the harmonic resolvent analysis of Padovan & Rowley (Reference Padovan and Rowley2022), use a time-varying base flow to capture first-order triadic interactions at specific frequencies, showing remarkable agreement between linear models and fully nonlinear direct numerical simulation. These methods account for the periodic structure induced by the fundamental wave, whereas LST captures only the indirect nonlinear effects associated with mean flow deformation. Similarly, the close correspondence between growth rates estimated from local linear theory and those derived from LES data in figure 5 supports this observation. Overall, these results indicate that while the linear instabilities initiate and drive the vortex-pairing process, the process itself is inherently nonlinear. As shown later in figure 10, the production term, associated with the linear effects, is energetically more significant, which explains why the linear model can effectively approximate the overall behaviour. Nevertheless, direct nonlinear interactions are dynamically important and essential for the vortex-merging process.
As both jets are statistically stationary flows, we use SPOD to extract the spatiotemporal coherent structures. To emphasise the dynamics in the initial shear layer, SPOD is computed using a weighting function that assigns zero weights to the region we wish to exclude. Specifically, zero weights are assigned to the areas outside the shear layer and beyond the end of the potential core. These weighting functions are shown in the top row of figure 6, where the black regions represent zero weights. Figure 6 shows the SPOD eigenspectra of the initially laminar and turbulent jets. Corresponding to the previously observed tones in figure 3, SPOD identifies the tonal peaks at
${St}=1.76$
, 0.88 and 0.44 for the initially laminar jet in figure 6(a). Interestingly, a broader peak is observed at the fundamental frequency
${St}=1.76$
, which has lower energy compared with the two subharmonics. Additionally, the spectra reveal the presence of ultra- and superharmonics at
${St} \approx 2.7$
and 3.5, respectively. We have confirmed from the spatial linear theory that
${St} = 3.5$
is indeed a higher harmonic, as there was no associated peak in the amplification rate. No distinct peaks are observed at the ultraharmonic frequencies
${St}=0.66$
and
$1.32$
(denoted by grey dotted lines). Later, in figure 12, we demonstrate that despite not being energetically important, they play an active role in triadic interactions. The SPOD spectrum of the turbulent jet is relatively broadband and exhibits a low-rank behaviour for
$ 0.3 \leqslant {St} \leqslant 1$
(Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018). Overall, figure 6 indicates that the presence of discrete tones leads to a much higher fluctuation level than the turbulent case.
SPOD eigenspectra with a focus on the shear layer until the end of the potential core: (a) initially laminar jet; (b) turbulent jet. The white-shaded area in the top row denotes the focus region of SPOD. Dashed lines indicate the three tones of the initially laminar jet. Dotted lines in (a) correspond to the ultraharmonics
${St}=0.66$
and 1.32, and the dotted line in (b) to the most energetic frequency of the turbulent jet.
Leading SPOD modes of the fundamental and four of its subharmonic frequencies: (a,b)
${St}=1.76$
; (c,d)
${St}=0.88$
; (e,f)
${St}=0.44$
; (g,h)
${St}=0.22$
; (i,j)
${St}=0.11$
. The left-hand column represents the streamwise velocity component
$u_x$
and the right-hand column represents the radial velocity
$u_r$
. The potential core and the jet width are indicated as lines of constant
$u_x$
at 95 % and 10 % of the jet velocity
$U_j$
, respectively.
The leading SPOD modes associated with
${St} = 1.76, 0.88, 0.44, 0.22$
and 0.11 are shown in figure 7. Due to the spreading of the shear layer, structures with lower frequencies are supported farther downstream than those associated with higher frequencies. At the fundamental frequency, the leading SPOD mode materialises as a Kelvin–Helmholtz wavepacket with compact support in the region
$0.8 \leqslant x \leqslant 1.2$
around the lipline. At a lower frequency of
${St} = 0.44$
, the leading SPOD mode exhibits larger spatial support within the range
$1 \leqslant x \leqslant 6$
and
$0 \leqslant r \leqslant 1$
. Notably, it exhibits high amplitude in
$2 \leqslant x \leqslant 3$
, where the streamwise velocity is concentrated on the centreline and the radial velocity on the lipline. Furthermore, the presence of this Kelvin–Helmholtz-type wavepacket structure is responsible for the accumulation of the RMS streamwise velocity on the centreline at
$x \approx 3$
. Similarly, the structures at other frequencies are also associated with the peaks of the RMS velocities in figure 7. For instance, the global maximum of the spatial structure of the radial velocity for
${St} = 0.88$
is at
$(x,r)\approx (1.5,0.5)$
, corresponding to the peak of the blue curve in figure 4(c).
3.3. Vortex pairing
(a–h) Time traces exemplifying two successive vortex-pairing events are visualised in terms of the azimuthal vorticity for the
$m=0$
component. The green, blue and magenta rectangles enclose the fundamental, subharmonic and second subharmonic vortices, respectively.
The vortex-pairing process is visualised in figure 8. Eight time instances of the vorticity
$\omega _{\theta }$
for the axisymmetric component
$m=0$
are shown. These snapshots follow two successive vortex-pairing events, where four fundamental vortices merge into two subharmonic vortices, eventually coalescing into a single vortex corresponding to the second subharmonic frequency. This process involves the formation of a larger vortex associated with half the frequency and half the wavenumber of the previous vortices. The highlighted vortical structures are unambiguously associated with their respective frequencies. This association is determined by considering their spatial support consistently with figures 4(c), 5(a) and 7(a–f). The first snapshot highlights two developing and two developed vortices enclosed in green rectangles, corresponding to a frequency of
${St} = 1.76$
. These vortices are formed due to the roll-up of the shear layer. The following three snapshots show the pairing of these vortices in the region
$1.2 \lesssim x \lesssim 2.1$
, denoted by the blue rectangle in figure 8(b–d). This vortex pairing results in the formation of the
${St} = 0.88$
vortex. Another instance of vortex pairing is evident in figure 8(c–e). Next, the two
${St}= 0.88$
vortices, denoted by blue rectangles in figure 8(e), undergo pairing to form the
${St} = 0.44$
vortex. Figure 8(e–f) demonstrates this process. In this process, the vortex at a more upstream location accelerates and catches up with the decelerating downstream vortex. Eventually, they wrap around each other and form a single vortex at
$x \approx 3.5$
. These consecutive vortex-pairing events are qualitative evidence of an inverse cascade, transferring energy from smaller to larger structures.
The
$x$
–
$t$
plots along the lipline showing the vorticity fluctuations,
$\omega ^{\prime }_{\theta }$
: (a) representative time interval; (b) conditional average of the SPOD-band filtered data about the spatial location
$x\approx 3.5$
. The successive vortex-pairing events shown in figure 8 are enclosed by the yellow box in (a). The green, blue and magenta lines correspond to the fundamental, subharmonic and second subharmonic vortices, respectively.
A different perspective on vortex pairing is presented in terms of the
$x$
–
$t$
plots in figure 9. Figure 9(a) shows the vorticity fluctuations,
$\omega _{\theta }^{\prime }$
, of the
$m=0$
component along the lipline in the time interval
$ 35 \leqslant t \leqslant 65$
. Merging lines with slopes greater than zero indicate the pairing of vortices. For instance, the vortex-pairing process in figure 8 is highlighted in yellow dashed lines. However, note that this represents a single event of vortex pairing chosen for its clarity in figure 8. In order to confirm that this vortex-pairing sequence is indeed a prevailing flow feature, we deploy a statistical perspective. To this end, we use conditional averaging and SPOD-band-pass filtering developed by Nekkanti & Schmidt (Reference Nekkanti and Schmidt2021). To isolate the two successive vortex-pairing events, we band-pass-filter the data using SPOD, retaining only the fundamental frequency and its first two subharmonics. Next, we select the location
$x=3.5, r=0.5$
and extract all local peaks that exceed 25 % of their global maximum. We obtain the
$i$
th realisation of the vortex-pairing sequence by collecting all the snapshots in the interval
$[t_0^{(i)}-20, t_0^{(i)}+10]$
, where
$t_0^{(i)}$
is the time instant of the
$i$
th local peak. Finally, the statistical representation of the vortex-pairing sequence is obtained by averaging over all realisations. The
$x$
–
$t$
plot along the lipline obtained from this conditional averaging is shown in figure 9(b). Evident here are the first and second vortex-pairing events at
$x\approx 1.5$
and
$x\approx 3$
, respectively. For illustrative purposes only, we use a polynomial curve fitting to demonstrate this process. The merging of two lines indicates the acceleration and deceleration of the vortices in the vortex-pairing event. We observe that the phase speed outside of the vortex-pairing events is constant. This is expected, as these result from the Kelvin–Helmholtz-type instability waves governed by the dispersion relation
$\omega /k = c_{ph}$
. Here, for the second subharmonic frequency,
$c_{ph}=0.55$
. The constant phase speeds of the fundamental wave and its subharmonics indicate phase locking, confirming the resonance condition of Cohen & Wygnanski (Reference Cohen and Wygnanski1987) that both waves must propagate at the same speed to interact and exchange energy.
3.4. SPOD-based spectral energy transfer
We now investigate the spectral energy balance of TKE as outlined in § 2.3. In particular, we use it to identify the net production and nonlinear transfer of TKE associated with the frequencies involved in the vortex-pairing process. This nonlinear transfer term quantifies the direct nonlinear effects. The nonlinear energy transfer term could alternatively be estimated using bispectral mode decomposition (BMD) (Schmidt Reference Schmidt2020). However, to maintain consistency with the modes employed in this analysis, we estimate the nonlinear energy transfer using the leading SPOD modes. Moreover, employing the SPOD modal basis ensures direct comparability among different terms such as production, nonlinear transfer and dissipation. In Appendix A, we tailor BMD to estimate nonlinear energy transfer and compare it with those estimated from SPOD. We find that both methods yield qualitatively similar results.
Figure 10 shows the production and nonlinear energy transfer at frequencies
${St} = 1.76,\ 0.88$
and 0.44. These terms are radially integrated and plotted as a function of the streamwise location. The production term is computed using (2.13). Positive production indicates energy gain from the mean flow, and negative production represents energy loss to the mean flow. Figure 10(a) shows that
${St} =1.76$
is the earliest to gain energy from the mean flow, and as it saturates,
${St} = 0.88$
begins its growth. Subsequently,
${St} = 0.88$
attains its global maximum at
$x = 1.25$
, which also corresponds to the location of the global minimum of
${St} = 1.76$
. Similarly, the global maximum of
${St} = 0.44$
and the global minimum of
${St}= 0.88$
are in close proximity. The saturation of the subharmonic frequency has been linked to the onset of vortex pairing by Ho & Huang (Reference Ho and Huang1982) and was confirmed by Hajj et al. (Reference Hajj, Miksad and Powers1992, Reference Hajj, Miksad and Powers1993). Our findings are in agreement with those studies. The production curves, in combination with figure 8, demonstrate this, where the location corresponding to the peak production of the subharmonic also marks the beginning of the merging process. Furthermore, the production of
${St} = 1.76,0.88$
and 0.44 becomes negative at
$x = 1.07$
, 1.64 and 2.64, respectively. These locations closely correspond to the location of the neutral stability points as predicted by LST (denoted by dashed lines; also see figure 5). The correspondence is expected as positive production is associated with the amplification rate and negative production is associated with the decay rate of each frequency.
(a) Production and (b) nonlinear energy transfer terms for
${St} = 1.76,\ 0.88$
and 0.44 integrated in
$r$
and as a function of streamwise location. Dashed lines indicate neutral stability points, predicted by LST, of the corresponding frequency. The black dotted line indicates the onset of nonlinear interactions.
The net nonlinear energy transfer, computed using (2.16), is shown in figure 10(b). The net nonlinear energy transfer is always negative for
${St}=1.76$
. On the other hand, for
${St}=0.88$
and 0.44, the energy is initially transferred into these frequencies and as the flow evolves downstream, energy is extracted from these frequencies. The blue and magenta curves in the region
$1.3 \leqslant x\leqslant 2.3$
exhibit a similar shape but are opposite in sign, indicating that nonlinear interactions result in energy loss for
${St} = 0.88$
and energy gain for
${St} = 0.44$
. This phenomenon is expected for a vortex-pairing process, as the energy is transferred from
${St}=0.88$
to its subharmonic,
${St} = 0.44$
. The black dotted line denotes the location for the onset of nonlinearity at
$x=0.72$
. This also corresponds to the location where
${ St}=0.88$
starts to grow (see figure 10
a), implying that the direct nonlinear effects trigger the growth of the subharmonic. Overall, these observations suggest that the spatially unstable fundamental grows linearly, gaining energy from the mean flow, followed by successive nonlinear growth and decay of its subharmonics. As pointed out earlier, this successive growth and decay is characteristic of the vortex-pairing process.
The production curves in figure 10 are significantly larger than the nonlinear transfer curves, indicating that subharmonic frequencies gain more energy from the mean flow than from their harmonics. Paschereit et al. (Reference Paschereit, Wygnanski and Fiedler1995) suggested that the subharmonic derives most of its energy from the mean flow, with fundamental–subharmonic interactions serving as a catalyst. However, they did not quantify these contributions. Our analysis bridges this gap by providing clear, quantitative evidence that the energy transfer from the mean flow to the subharmonic is energetically dominant. Additionally, both production and nonlinear transfer for the fundamental frequency are smaller in magnitude than those for the subharmonics, reaffirming that the fundamental is not energetically dominant. However, the nonlinear transfer curves show that quadratic interactions extract energy from the fundamental and supply energy to the subharmonic. This highlights the fundamental’s dynamical importance as it supplies energy to the subharmonic through these direct nonlinear interactions.
Figure 11 shows the spatial fields of production, dissipation and net nonlinear transfer at
${St}=1.76$
, 0.88 and 0.44. These fields are computed using (2.13), (2.14) and (2.16), respectively. As in figure 10, the production fields in figure 11(a,d,g) indicate that the fundamental frequency and its subharmonics initially gain energy from the mean flow and subsequently transfer energy back to the mean flow. These fields also reveal that this energy transfer from the mean flow is localised in the region about the lipline. Figure 11(b,e,h) shows that the spatial dissipation fields are significantly lower in magnitude and, hence, have minimal impact on the energy budget of the vortex-pairing process. The net nonlinear energy transfer fields, shown in figure 11(c,f,i), exhibit multilobe structures. These fields elucidate the spatial dependence of the nonlinear kinetic energy transfer. Figure 11 highlights that individual terms in the spectral TKE budget can exhibit high local values, yet their overall contributions may remain small when integrated across the entire spatial domain.
Nonlinear energy transfer using SPOD: transfer term bispectrum for (a) entire domain and (b) shear-layer subdomain
$\varOmega _r$
. Spatial fields for the triads (c) (1.76, −0.88, 0.88) and (d) (0.88, −0.44, 0.44) are compared with the radially integrated TKE of
${St}=1.76,\ 0.88$
and 0.44.
Figures 12(a) and 12(b) show the triadic energy transfer obtained by integrating (2.15) over the entire spatial domain and a domain
${\varOmega }_r$
that focuses on the initial shear layer within the first two jet diameters in
$x, r \in [0,2] \times [0,6]$
, respectively. In this bispectral representation, the abscissa represents the first frequency of the triad,
$St_1$
, corresponding to contributions from
$\hat {u}(St_1)$
. On the ordinate, the second frequency of the triad,
$St_2$
, is depicted, whose contribution arises from
${\partial \hat {u}}/{\partial x_j}(St_2)$
. The third (or resulting) frequency is determined by the sum of these two frequencies,
$St_3 = St_1 + St_2$
. Different combinations of frequencies
$(St_1,St_2)$
interact to generate the same frequency
$(St_1 + St_2$
= constant) along the diagonals of slope
$-1$
in this representation. In figure 12(a,b), positive values (red colour) denote energy transfer to frequency
${St}_3$
from frequencies
${St}_1$
and
${St}_2$
, while negative values (blue colour) indicate the extraction of energy from
${St}_3$
by
${St}_1$
and
${St}_2$
. The two triads with the highest positive and negative intensity are
$(0.88,-0.44,0.44)$
and
$(0.44,0.44,0.88)$
. These triads exhibit positive and negative energy transfer, respectively. However, both triads convey the same information. The triad
$(0.88,-0.44,0.44)$
has a positive value, signifying the transfer of energy to
${St}= 0.44$
. On the other hand, the triad
$(0.44,0.44,0.88)$
has a negative value, indicating that
${St} = 0.44$
extracts energy from
${St} = 0.88$
. These triads clearly suggest that the energy is transferred from the subharmonic to the second subharmonic frequency. Other significant triads include (0.44, 0.88, 1.32), (1.32, −0.44, 0.88) and (1.32, −0.88, 0.44), highlighting the presence of the ultraharmonic frequency
${St} = 1.32$
. The occurrence of
${St}=1.32$
is interesting in itself as its relevance was not apparent from the SPOD analysis, i.e. no distinct peak was found at
${St}=1.32$
. The nonlinear energy transfer analysis reveals that the ultraharmonic
${St}=1.32$
is created from the sum interaction of 0.44 and 0.88.
In figure 12(a), the triads (1.76, −0.88, 0.88) and (0.88, 0.88, 1.76) exhibit lower amplitudes. This is because the fundamental frequency
${St} = 1.76$
is dynamically significant but not energetically dominant. Therefore, to shed light on the direct nonlinear interactions of the fundamental frequency, we narrow our focus to the first two jet diameters. Figure 12(b) now illuminates the triad (0.88, 0.88, 1.76). This triad exhibits energy transfer behaviour characteristic of the vortex-pairing process, with a loss of energy at
${ St}=1.76$
and a gain at
${St}=0.88$
. These findings, in combination with observations from figures 5, 8 and 9, suggest that during the vortex-pairing process, the nonlinear interactions cause the energy to be initially transferred from the fundamental (
${St}=1.76$
) to its first subharmonic (
${St}=0.88$
), and then from the first subharmonic to the second subharmonic (
${St}=0.44$
). Our study, for the first time, quantifies the triadic energy transfer between the fundamental and its subharmonic waves.
The spatial fields of
${\mathscr{t}}_{{nl}}^{{\textit{rank}\text{-}1}}$
for the dynamically dominant and energetically significant triad are shown in figures 12(c) and 12(d), respectively. The TKE of the fundamental, subharmonic and second subharmonic frequencies are overlaid on the contours of the
${\mathscr{t}}_{{nl}}^{{\textit{rank}\text{-}1}}$
field. For the (1.76, −0.88, 0.88) triad, the
${\mathscr{t}}_{{nl}}^{{\textit{rank}\text{-}1}}$
field is concentrated in the region
$0.95 \lesssim x \lesssim 1.4$
, corresponding to where the fundamental frequency decays and the subharmonic frequency grows. A similar trend is evident in figure 12(d), demonstrating that the nonlinear energy transfer of the triad (0.88, −0.44, 0.44) is localised to the region of subharmonic decay and second subharmonic growth. In summary, SPOD-based transfer analysis allowed us to systematically catalogue the triadic energy transfer among the SPOD modes, thus establishing a direct link between energy flow analysis and the phenomenon of vortex pairing.
4. Discussion and conclusions
Many studies have investigated the vortex-pairing process in mixing layers, free shear layers and natural and forced jets. It is a nonlinear process that involves the merging of two smaller vortices into a larger vortex with half the frequency. Previous studies have focused on hydrodynamic instabilities (Michalke Reference Michalke1965; Kelly Reference Kelly1967), phase locking of vortices (Monkewitz Reference Monkewitz1988; Husain & Hussain Reference Husain and Hussain1995) and specific aspects of energy transfer (Mankbadi Reference Mankbadi1985; Paschereit et al. Reference Paschereit, Wygnanski and Fiedler1995) during the vortex-pairing process. This study presents a detailed quantitative analysis of the vortex-pairing process using high-fidelity simulations of initially laminar jets and provides direct data-driven verification of previous experimental and theoretical results, particularly Monkewitz’s (Reference Monkewitz1988) parametric resonance mechanism. While recent methods aim to characterise or optimise triadic interactions, including BMD (Schmidt Reference Schmidt2020), triadic orthogonal decomposition (Yeung, Chu & Schmidt Reference Yeung, Chu and Schmidt2024) and related approaches (Jin et al. Reference Jin, Symon and Illingworth2021; Biswas et al. Reference Biswas, Cicolin and Buxton2022; Kinjangi & Foti Reference Kinjangi and Foti2023; Freeman, Martinuzzi & Hemmati Reference Freeman, Martinuzzi and Hemmati2024), the present analysis focuses on nonlinear interactions of the statistically most energetic coherent structures.
The analysis was conducted for two LES of initially laminar and turbulent jets at
${Re}= 50\,000$
. The PSDs and local linear theory confirm the expectation that the hydrodynamic instabilities are more pronounced in the initially laminar jet. These instabilities cause the initially laminar jet to transition quickly into turbulence within the first two jet diameters. A comparison of the two jets shows that the initially laminar jet starts to develop from a more downstream location but grows more rapidly. This delayed but faster growth is triggered by the hydrodynamic instabilities. The boundary layer of the turbulent jet is tripped inside the nozzle, and the shear layer is already fully turbulent at the nozzle’s exit. Therefore, the spreading rate of its shear layer is more gradual compared with that of the initially laminar jet. Additionally, vortex pairing is present and strongly influences the dynamics for the initially laminar jet but it is not observed in the turbulent jet. Local linear theory identifies the fundamental frequency, i.e. the frequency with the largest spatial growth rate, as
${St}=1.76$
for the initially laminar jet. At this frequency, the shear layer rolls up into vortices. The tones at
${St}=0.88$
and 0.44 are not distinguished as distinct peaks in the local linear theory and are therefore a consequence of the nonlinear interactions between the fundamental frequency and the subharmonic frequency, and between the subharmonic and the second subharmonic frequency, respectively. A remarkable agreement was found between the spatial growth rates predicted from theory and empirically from data. This validates the parallel flow assumption and strongly suggests that the purely empirical approach can be utilised to estimate spatial growth rates in turbulent flows. Similarly, LST can provide an accurate prediction of growth rates of coherent structures, even when weak nonlinear effects are present, particularly in the initial shear layer.
Different visualisation techniques give a clear phenomenological understanding of the vortex-pairing process, wherein the accelerating upstream vortex and decelerating downstream vortex merge to form a larger vortex with twice the wavelength of the preceding ones. The process starts upstream with two vortices associated with the fundamental frequency merging to form a vortex associated with the subharmonic frequency. Subsequently, the two subharmonic vortices merge to result in a second subharmonic vortex. The second subharmonic frequency is energetically the most significant, while the fundamental frequency, despite its low energy, is dynamically the most important as it dictates the entire nonlinear dynamics. This dynamical significance can be inferred from (i) LST and (ii) the three-way resonance between the fundamental and its subharmonics, which determines the other frequencies observed further downstream. While LST predicts the fundamental frequency to have the highest amplification rate, it does not indicate the presence of any other peaks. As the flow is convectively dominated, the fundamental frequency primarily influences the dynamics, confining it mainly to the subharmonic frequencies, which result from direct nonlinear interactions. The fundamental and subharmonic waves propagate at the same phase speed, allowing sufficient time for interaction, consistent with the resonance condition proposed by Cohen & Wygnanski (Reference Cohen and Wygnanski1987). This is a manifestation of the importance of distinguishing between dynamical relevance and energetic significance, as highlighted, among others, by Schmid (Reference Schmid2010).
SPOD-based spectral TKE analysis was performed to characterise the energy transfer during the vortex-pairing process. The focus was on the interactions between the mean flow, fundamental frequency and subharmonic frequencies. The production and nonlinear transfer terms are major contributors to energy transfer, while the effect of dissipation is negligible. The energy transfer between non-zero frequencies and the mean flow is quantified using the production term, while the triadic energy transfer among different frequencies is measured using the nonlinear transfer term. The fundamental frequency gains energy from the mean flow, while its subharmonics gain energy from both the mean flow and their harmonic. Quantitatively, the energy gained from the mean flow is greater than from its harmonic. As two fundamental vortices merge, there is a backscatter of energy from the fundamental to its subharmonic. This process repeats itself for higher subharmonics. The backscatter of energy from the fundamental to the subharmonic occurs when the fundamental wave saturates, initiating direct nonlinear interactions. During this process, the subharmonic wave simultaneously gains energy both from the mean flow and through direct nonlinear interactions with the fundamental wave. A few studies (Hajj et al. Reference Hajj, Miksad and Powers1992, Reference Hajj, Miksad and Powers1993) show that the dominant interaction during vortex pairing is between the subharmonic and the fundamental frequency, and other studies (Mankbadi Reference Mankbadi1985; Paschereit et al. Reference Paschereit, Wygnanski and Fiedler1995) show that the subharmonic gains most of its energy from the mean flow and the fundamental–subharmonic interaction only acts as a catalyst. Our findings agree with all these studies and present a comprehensive analysis of the vortex-pairing process. In an energetic sense, the energy extracted by the subharmonic from the mean flow is more significant. However, from a dynamical perspective, the fundamental–subharmonic interaction is more significant.
This study provides a quantitative analysis that supports the previous findings outlined below from a data-driven perspective:
-
(i) The hydrodynamic instabilities initiate the transition into turbulence, causing the shear layer to grow rapidly.
-
(ii) Through exponential growth, the fundamental frequency attains significant amplitude and triggers the roll-up of the shear layer. As the fundamental frequency grows, it extracts energy from the mean flow.
-
(iii) The saturation and subsequent decay of the fundamental frequency mark the onset of vortex pairing.
-
(iv) As the vortex pairing continues, the subharmonic frequency acquires energy linearly from the mean flow and nonlinearly through backscatter from the fundamental frequency. The eventual saturation of the subharmonic frequency signals the completion of the vortex-pairing process.
-
(v) Processes analogous to steps (ii) and (iii) then repeat to create higher subharmonics.
Acknowledgements
The main LES calculations were performed at SDSC through UC@HPC. The authors would like to thank the anonymous reviewers for their insightful comments.
Funding
We acknowledge support from ONR under grant N000142312650. O.T.S. and A.N. are grateful for support from the Office of Naval Research under grant N00014-20-1-2311.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Nonlinear energy transfer using SPOD and BMD
Here, we qualitatively compare the nonlinear energy transfer estimated from SPOD and BMD. BMD is a modal decomposition technique that can be understood as an extension of classical bispectral analysis to multidimensional and multivariate data. It identifies the spatially coherent structures associated with triadic interactions by maximising the integrated point-wise bispectrum:
\begin{equation} b(f_k,f_l) = E\left [\int _\varOmega \hat {{q}}_1^{*}({x},f_k) \circ \hat {{q}}_2^{*}({x},f_l) \circ \hat {{q}}_3({x},f_k+f_l) {\textrm d}{x} \right ]. \end{equation}
For further details on computing the mode bispectrum, the reader is referred to Schmidt (Reference Schmidt2020). Recently, BMD has been applied to identify the triadic energy cascade in forced jets (Nekkanti et al. Reference Nekkanti, Maia, Jordan, Heidt, Colonius and Schmidt2022, Reference Nekkanti, Schmidt, Maia, Jordan, Heidt and Colonius2023b ) and to investigate the interplay between streaks and Kelvin–Helmholtz wavepackets (Nekkanti et al. Reference Nekkanti, Pickering, Schmidt and Colonius2025). Here, we tailor BMD to estimate the nonlinear energy transfer by maximising the following point-wise integral:
\begin{equation} b_{ {nl}}^{\textit{BMD}}(f_k,f_l) = E\left [\int _\varOmega \hat {{u}}_i^{*}({x},f_k) \circ \frac {\partial \hat {{u}}_j}{\partial {x}_i}^{*}({x},f_l) \circ \hat {{u}}_j({x},f_k+f_l) {\textrm d}{x} \right ]. \end{equation}
Comparison of nonlinear energy transfer: (a) SPOD; (b) BMD.
Figure 13 shows the nonlinear energy transfer estimated using SPOD (2.15) and BMD (A2). Qualitatively similar trends are observed. In particular, the direction of energy transfer is the same for both methods, i.e. all the significant triads exhibit the same sign. Additionally, the triad with highest intensity for both methods is (0.88, −0.44, 0.44), which is representative of energy transfer from subharmonic to second subharmonic. Note that the BMD-based nonlinear energy transfer exhibits more noise in comparison with that of the SPOD-based
${\mathscr{t}}_{{nl}}$
. This is because the BMD in comparison with SPOD requires more data for convergence.
Appendix B. Spectral TKE derivation
We start with the momentum equation:
\begin{align} \frac {\partial {u}_j}{\partial t}+{u}_i \frac {\partial {u}_j}{\partial x_i} &=-\frac {\partial {p}}{\partial x_j}+\frac {2}{R e} \frac {\partial {s}_{i j}}{\partial x_i}. \end{align}
By substituting the Reynolds decomposition,
$u =\bar {u} + u^{\prime }$
, into the momentum equation (B1) and taking its mean, we obtain the mean-momentum equation:
\begin{align} \frac {\partial \bar {u}_j}{\partial t}+\bar {u}_i \frac {\partial \bar {u}_j}{\partial x_i}+\frac {\partial }{\partial x_i} \overline {u_i^{\prime } u_j^{\prime }} &=-\frac {\partial \bar {p}}{\partial x_j}+\frac {2}{R e} \frac {\partial \bar {s}_{i j}}{\partial x_i}. \end{align}
Next, subtracting (B2) from (B1), we obtain the fluctuating-momentum equation:
\begin{align} \frac {\partial u_j^{\prime }}{\partial t}+\bar {u}_i \frac {\partial u_j^{\prime }}{\partial x_i}+u_i^{\prime } \frac {\partial u_j^{\prime }}{\partial x_i} &=-u_i^{\prime } \frac {\partial \bar {u}_j}{\partial x_i}+\frac {\partial }{\partial x_i} \overline {u_i^{\prime } u_j^{\prime }}-\frac {\partial p^{\prime }}{\partial x_j}+\frac {2}{R e} \frac {\partial s_{i j}^{\prime }}{\partial x_i}. \end{align}
For statistically stationary, azimuthally homogeneous flows, (B3) can be Fourier-transformed in both time and the azimuthal direction, such that
${u}^{\prime }(x,r,\theta ,t) = \sum \nolimits _{f} \sum \nolimits _m \hat {u}(x,r) {\textrm e}^{{\textrm i} (m \theta -2\pi f t )}$
, which results in the spectral momentum equation for each azimuthal wavenumber
$m$
and frequency
$f$
:
\begin{align} \left [\frac {\partial \hat {u}_j}{\partial t}+\bar {u}_i \frac {\partial \hat {u}_j}{\partial x_i}+ \widehat {u_i\frac {\partial {u}_j}{\partial x_i}} =-\hat {u}_i\frac {\partial \bar {u}_j}{\partial x_i}-\frac {\partial \hat {p}}{\partial x_j}+\frac {2}{R e} \frac {\partial \hat {s}_{i j}}{\partial x_i}\right ]_{m,\omega }. \\[-18pt] \nonumber \end{align}
Here, the term
is a nothing but a Fourier transform of the convective term. Finally, the spectral TKE equation at each azimuthal wavenumber
$m$
and frequency
$f$
can be obtained by multiplying (B4) with
$\hat {u}^{*}_j (m,f)$
and invoking continuity,
$\partial \hat {u}_j/ \partial x_j =0$
:
\begin{equation} \frac {\partial \hat {k}}{\partial t}= {\mathcal{R}}\! \left [\!-\bar {u}_i \frac {\partial \hat {k}}{\partial x_i} - \hat {u}_j^*\widehat { u_i\frac {\partial u_j}{\partial x_i}} -\hat {u}_j^*\hat {u}_i\frac {\partial \bar {u}_j}{\partial x_i} -\frac {2}{{Re}} \hat {s}_{i j}^* \hat {s}_{i j} -\frac {\partial }{\partial x_j}\left (\hat {u}_j^*\hat {p}\right ) +\frac {2}{{Re}} \frac {\partial }{\partial x_i}\! \left (\hat {u}_j^*\hat {s}_{ij}\! \right )\! \right ], \end{equation}
where
$\hat {k}=\hat {u}_j(m,f)\hat {u}_j(m,f)$
and
${\mathcal{R}}[\boldsymbol{\cdot }]$
denotes the real part.
Open access
