Nonlinear dynamics of vortex pairing in transitional jets

Abstract

This study investigates the onset of linear instabilities and their later nonlinear interactions in the shear layer of an initially-laminar jet using a combination of stability analysis and data from high-fidelity flow simulations. We provide a complete picture of the vortex-pairing process. Hydrodynamic instabilities initiate the transition to turbulence, causing the shear layer to spread rapidly. In this process, the shear layer rolls up to form vortices, accompanied by the exponential growth of the fundamental frequency. As the fundamental frequency grows, it gains energy from the mean flow. Subsequently, as it saturates and begins to decay, the fundamental vortices start to pair. During this vortex pairing process, the subharmonic vortex acquires energy both linearly from the mean flow and nonlinearly through a reverse cascade from the fundamental. The process concludes when the subharmonic vortex eventually saturates. Similarly, two subharmonic vortices merge to form a second subharmonic vortex. Our results confirm Kelly (1967)’s hypothesis of a resonance mechanism between the fundamental and subharmonic, which supplies energy to the subharmonic. In this multi-tonal, convective-dominated flow, we clarify the ambiguity surrounding the fundamental frequency by demonstrating that the spatially most amplified frequency should be considered fundamental, rather than the structure associated with the spectral energy peak. For the initially-laminar jet considered here, the fundamental frequency corresponds to the fourth largest spectral peak, highlighting the important distinction between the energetically and dynamical significance of a tone. Despite its low energy, the fundamental frequency is dynamically dominant as it determines all other spectral peaks and supplies energy to the subharmonics through a reverse energy cascade.

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’s exit. In particular, the flow field is sensitive to parameters such as momentum thickness, fluctuation level, and whether the boundary layer is laminar or turbulent. Several experimental (Hill Jr et al., 1976; Hussain & Zedan, 1978a, b; Husain & Hussain, 1979; Bridges & Hussain, 1987; Zaman, 1985, 2012; Fontaine et al., 2015) and numerical studies (Bogey & Bailly, 2005, 2010; Bogey et al., 2012; Kim & Choi, 2009; Brès et al., 2018) 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, 1985; Bridges & Hussain, 1987; Bogey & Bailly, 2010). Zaman (1985) showed that the initially-laminar jet exhibits a 4 dB increase in the radiated noise compared to its turbulent counterpart. Bogey & Bailly (2010) 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, 1974; Winant & Browand, 1974; Ho & Huang, 1982; Metcalfe et al., 1987; Moser & Rogers, 1993) and jet flows (Becker & Massaro, 1968; Zaman & Hussain, 1980; Hussain & Zaman, 1980; Meynart, 1983). It significantly contributes to turbulent mixing (Brown & Roshko, 1974), the production of Reynolds stresses (Zaman & Hussain, 1980), entrainment (Winant & Browand, 1974), and triggers the transition to turbulence (Ho & Huang, 1982; Moser & Rogers, 1993). Given its importance, vortex pairing has been the subject of numerous studies. In jets, it was first visualized by Becker & Massaro (1968). In the seminal study by (Crow & Champagne, 1971), 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 (1980) and Hussain & Zaman (1980) reveal that pairing occurs at two distinct frequencies: one around $\mbox{\it St}_{\theta}=f\theta/U\approx 0.012$, termed the shear layer mode, and the other around $\mbox{\it 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 (1980) forced the 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 (1982) observed multiple pairing when the flow was forced at the subharmonic frequency.

The nonlinear mechanisms behind the growth of the subharmonic during the vortex pairing process have been the focus of many studies (Kelly, 1967; Monkewitz, 1988; Paschereit et al., 1995; Husain & Hussain, 1995). Kelly (1967) proposed a resonance mechanism that provides energy to the subharmonics. Monkewitz (1988) used weakly nonlinear spatial theory to further support Kelly (1967)’s resonance mechanism. They demonstrated that the fundamental mode needs to reach a critical amplitude before both the fundamental and subharmonic can phase-lock, resulting in a energy transfer to the subharmonic. Experiments by Hajj et al. (1992) show that the subharmonic gains energy through this resonance mechanism. Works by Husain & Hussain (1995) and Cho et al. (1998) demonstrate that vortex pairing can be either enhanced or attenuated by controlling the phase difference between the fundamental and the subharmonic. Mankbadi (1985) employed an energy-integral method to show that the subharmonic gains energy both from the fundamental and the mean flow. Paschereit et al. (1995) 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.

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 and 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 vs 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 (1964, 1965) on a hyperbolic tangent profile. Here, we apply LST to the temporally averaged mean flow of the jet. This approach corresponds to applying the parallel flow assumption locally to the zero-frequency component, which has been utilized in cylinder wakes (Pier, 2002; Barkley, 2006) and jets (Suzuki & Colonius, 2006; Gudmundsson & Colonius, 2011; Schmidt et al., 2017). The mathematical framework of spectral property orthogonal decomposition (SPOD) dates back to the work of Lumley (1967, 1970). SPOD identifies the most energetic coherent structures at each time scale. Early applications of SPOD include the work of Glauser et al. (1987); Glauser & George (1992); Delville (1994). More recently, this method has attracted significant interest following its application to large flow databases by Schmidt et al. (2018) and the establishment of its relationship to other methods by Towne et al. (2018). Since then, SPOD has become a mainstay of physical exploration, in particular for identifying the different modal and non-modal instabilities (Schmidt et al., 2018; Nogueira et al., 2019; Pickering et al., 2020).

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 analyze the energy transfer between different scales, including the Karman-Howarth equation (Von Karman & Howarth, 1938; Danaila et al., 1999; Hill, 2001) and bispectrum analysis (Lii et al., 1976; Kim & Powers, 1979; Herring, 1980). The former is based on structure functions, and the latter on third-order statistics. is the frequency domain representation of third-order moments. Recently, a modal decomposition based on bispectral analysis was developed by Schmidt (2020). 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, 2016; Cho et al., 2018; Lee & Moser, 2019; Gomé et al., 2023) 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 et al., 2021; Jin et al., 2021), Fourier modes (Nekkanti et al., 2023), dynamic mode decomposition (DMD) modes (Kinjangi & Foti, 2023), and optimal mode decomposition (OMD) modes (Biswas et al., 2022) to estimate the spectral energy budget.

The paper is organized as follows: In §2, the methodologies of LST, SPOD, and SPOD-based spectral analysis budget 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 the linear stability theory. This LST determines the spatial growth rates of the corresponding frequency. Here, we use it to identify the most unstable frequency. LST is also referred to as local linear theory, and throughout this paper, we will use the terms linear stability theory 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 decompostion, $\mathbf{q}\pqty{x,r,\theta,t}=\bar{\mathbf{q}}\pqty{x,r}+\mathbf{q}^{\prime}\pqty{x,r,\theta,t}$, where $\bar{\mathbf{q}}$ is the mean flow and $\mathbf{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

$$\mathbf{q}^{\prime}(x,r,\theta,t)=\tilde{\mathbf{q}}(r)e^{i\pqty{\alpha x+m\theta-\omega t}},$$ (1)

where the streamwise wavenumber $\alpha$ is complex, and $\tilde{\mathbf{q}}$ is the radial profile. Linearizing the Navier-Stokes equation about the base flow yields the equation,

$$\pqty{i\omega\mathbf{I}+\mathbf{L}}\tilde{\mathbf{q}}=0,$$ (2)

where $\mathbf{I}$ is the identity matrix and $\mathbf{L}$ is the linearized 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 (2007). 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 using the methodology of Maia et al. (2021, 2022). Solving this eigenvalue problems yields monochromatic amplification rates, $\alpha_{i}$; For $\alpha_{i}<0$, the disturbances will grow exponentially downstream, whereas for $\alpha_{i}>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 optimized 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, 1967). For the mathematical derivation and computational details, refer to Towne et al. (2018) and Schmidt & Colonius (2020). The SPOD modes and eigenvalues are obtained by solving the eigenvalue problem

$$\int_{\Omega}\bm{S}\pqty{\mathbf{x},\mathbf{x}^{\prime},f}\mathbf{W}\pqty{\mathbf{x}^{\prime}}\bm{\Phi}\pqty{\mathbf{x}^{\prime},f}\differential\mathbf{x}^{\prime}=\bm{\Phi}\pqty{\mathbf{x},f}\bm{\lambda}(f),$$ (3)

where $\bm{S}$ is the cross-spectral density matrix, $\mathbf{W}$ the positive definite matrix that accounts for component-wise and numerical quadrature weights, $\bm{\Phi}$ the SPOD modes and $\bm{\lambda}$ the eigenvalues. The modes $\bm{\phi}^{(i)}(\mathbf{x},f)$ and eigenvalues $\lambda^{(i)}(f)$ are sorted by energy, where $i$ is the mode number index. At each frequency, the SPOD modes are orthogonal in space,

$$\int_{\Omega}\bm{\phi}^{(i)}(\mathbf{x},f)\mathbf{W}\pqty{\mathbf{x}}\bm{\phi}^{(j)}(\mathbf{x},f)\differential\mathbf{x}=\delta_{ij},$$ (4)

where $\delta_{ij}$ is the Kronecker delta function. The Fourier modes of each flow realization are expanded as

$$\mathbf{\hat{q}}(x,f)=\sum\limits_{i}a^{(i)}(f)\bm{\phi}^{(i)}(\mathbf{x},f).$$ (5)

Recently, Nekkanti & Schmidt (2021) proposed a convolution approach that computes time-continuous expansion coefficients

$$a^{(i)}(f,t)=\int_{\Delta T}\int_{\Omega}\quantity({\bf\it\phi}^{(i)}(\mathbf{x},f))^{*}\mathbf{W}(\mathbf{x})\textbf{q}(\mathbf{x},t+\tau)e^{-i2\pi f\tau}\differential x\differential\tau,$$ (6)

which facilitates obtaining the time-continuous Fourier modes

$$\mathbf{\hat{q}}(\mathbf{x},f,t)\approx\sum\limits_{i}a^{(i)}(f,t){\bf\it\phi}^{(i)}(\mathbf{x},f).$$ (7)

2.3 SPOD-based spectral energy budget

The optimality of the SPOD expansion, equation (7), 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,

$$\frac{\partial\hat{k}}{\partial t}=\mathcal{R}\bqty{-\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}_{ij}^{*}\hat{s}_{ij}}}_{\mathcal{D}}-\frac{\partial}{\partial x_{j}}\pqty{\overline{\hat{u}_{j}^{*}\hat{p}}}+\frac{2}{Re}\frac{\partial}{\partial x_{i}}\pqty{\overline{\hat{u}_{j}^{*}\hat{s}_{ij}}}},$$ (8)

where, $\hat{s}_{ij}=1/2\pqty{\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$. 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 scale-specific production term

$$\mathcal{P}\pqty{f,m}=-\mathcal{R}\overline{\bqty{\hat{u}_{j}^{*}\pqty{f,m}\hat{u}_{i}\pqty{f,m}\frac{\partial\bar{u}_{j}}{\partial x_{i}}}},$$ (9)

the scale-specific dissipation term

$$\mathcal{D}\pqty{f,m}=-\frac{2}{Re}\mathcal{R}\overline{\bqty{\pqty{\frac{\partial\hat{u}_{i}^{*}}{\partial x_{j}}\pqty{f,m}+\frac{\partial\hat{u}_{j}^{*}}{\partial x_{i}}\pqty{f,m}}\pqty{\frac{\partial\hat{u}_{i}}{\partial x_{j}}\pqty{f,m}+\frac{\partial\hat{u}_{j}}{\partial x_{i}}\pqty{f,m}}}},$$ (10)

and the scale-specific nonlinear transfer term

$$\mathcal{T}_{\text{nl}}\pqty{f,m}=-\mathcal{R}\overline{\bqty{\hat{u}_{j}^{*}\pqty{f,m}\widehat{u_{i}\frac{\partial u_{j}}{\partial x_{i}}}\pqty{f,m}}}.$$ (11)

The term $\widehat{u_{i}\frac{\partial u_{j}}{\partial x_{i}}}$ in (11) entails the contributions from all other frequencies and azimuthal wavenumbers to frequency $f$ and azimuthal wavenumber $m$. Following Cho et al. (2018), 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 $\widehat{u_{i}\frac{\partial u_{j}}{\partial x_{i}}}$ using the discrete convolution to obtain

$$\mathcal{T}_{\text{nl}}\pqty{f_{3},m_{3}}=-\mathcal{R}\bqty{\overline{\hat{u}_{j}^{*}\pqty{f_{3},m_{3}}\sum\limits_{\begin{subarray}{c}f_{1}+f_{2}=f_{3}\\ m_{1}+m_{2}=m_{3}\end{subarray}}\hat{u}_{i}\pqty{f_{1},m_{1}}\frac{\partial\hat{u}_{j}}{\partial x_{i}}\pqty{f_{2},m_{2}}}}.$$ (12)

To obtain a bispectral representation of the energy transfer, we further split the sum in (12) into individual frequency and wavenumber components that are triadically compatible

$$\mathscr{t}_{\text{nl}}\pqty{f_{1},f_{2},m_{1},m_{2}}=-\mathcal{R}\bqty{\overline{\hat{u}_{j}^{*}\pqty{f_{1}+f_{2},m_{1}+m_{2}}\hat{u}_{i}\pqty{f_{1},m_{1}}\frac{\partial\hat{u}_{j}}{\partial x_{i}}\pqty{f_{2},m_{2}}}},$$ (13)

such that $\mathcal{T}_{\text{nl}}(f_{3},m_{3})=\sum_{\begin{subarray}{c}f_{1}+f_{2}=f_{3}\\ m_{1}+m_{2}=m_{3}\end{subarray}}\mathscr{t}_{\text{nl}}\pqty{f_{1},f_{2},m_{1},m_{2}}$.

At each frequency, the dominant flow structure is represented by the leading SPOD mode. To characterize the production and nonlinear kinetic energy transfer of these structures, we use a rank-1 approximation of the velocity field, following equation (7). 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>0$ do not play a significant role in the dynamics of interest. Substituting equation (7) into equations (9), (10), (12), and (13) we get

$$\mathcal{P}^{\text{rank-1}}\pqty{f_{1}}=-\overline{a^{*}\pqty{f_{1}}{a}\pqty{f_{1}}}\phi_{j}^{*}\pqty{f_{1}}\phi_{i}\pqty{f_{1}}\frac{\partial\bar{u_{j}}}{\partial x_{i}}=-\lambda(f_{1})\phi_{j}^{*}\pqty{f_{1}}\phi_{i}\pqty{f_{1}}\frac{\partial\bar{u}_{j}}{\partial x_{i}},$$ (14)

$$\mathcal{D}^{\text{rank-1}}\pqty{f_{1}}=-\lambda(f_{1})\pqty{\frac{\partial\phi_{i}^{*}}{\partial x_{j}}\pqty{f_{1}}+\frac{\partial\phi_{j}^{*}}{\partial x_{i}}\pqty{f_{1}}}\pqty{\frac{\partial\phi_{i}}{\partial x_{j}}\pqty{f_{1}}+\frac{\partial\phi_{j}}{\partial x_{i}}\pqty{f_{1}}},$$ (15)

$$\mathscr{t}_{\text{nl}}^{\text{rank-1}}\pqty{f_{1},f_{2}} =-\mathcal{R}\pqty{\overline{a^{*}\pqty{f_{3}}a\pqty{f_{1}}a\pqty{f_{2}}}\phi_{j}^{*}\pqty{f_{3}}\phi_{i}\pqty{f_{1}}\frac{\partial\phi_{j}}{\partial x_{i}}\pqty{f_{2}}},$$ (16)

$$\mathcal{T}_{\text{nl}}^{\text{rank-1}}(f_{3}) =-\mathcal{R}\pqty{\sum\limits_{f_{1}+f_{2}=f_{3}}\overline{a^{*}\pqty{f_{3}}a\pqty{f_{1}}a\pqty{f_{2}}}\phi_{j}^{*}\pqty{f_{3}}\phi_{i}\pqty{f_{1}}\frac{\partial\phi_{j}}{\partial x_{i}}\pqty{f_{2}}}.$$ (17)

The term $\mathcal{P}^{\text{rank-1}}(f_{1})$ signifies energy transfer from frequency $f_{1}$ to the mean flow. $\mathcal{D}^{\text{rank-1}}(f_{1})$ represents energy dissipation at $f_{1}$, $\mathscr{t}_{\text{nl}}^{\text{rank-1}}(f_{1},f_{2})$ denotes nonlinear interactions between $f_{1}$ and $f_{2}$, and $\mathcal{T}_{\text{nl}}^{\text{rank-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 of subsonic jets are conducted at a Reynolds number of $\Rey=U_{j}D/\nu=50000$. Here, $U_{j}$ is 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 large-eddy simulations are performed using the compressible flow solver “Charles” developed at Cadence, formerly Cascade Technologies (Brès et al., 2017, 2018) and the reader is referred to Brès et al. (2017, 2018) 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 et al. (2022). This experimental nozzle geometry is meshed using the same strategy as by Brès et al. (2018), 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

Figure 1: Comparison of the initially-laminar and turbulent jets: instantaneous fluctuating streamwise velocity field of ($a$) initially-laminar jet, ($b$) turbulent jet; 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.

We compare the initially laminar and turbulent jets using instantaneous visualizations and RMS of the fluctuating streamwise velocity in figure 1. The 95% (white solid line) and 5% (white dashed line) 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 figure 1($a$) and ($b$). The hallmark of the untripped jet is the initially-laminar shear layer that is easily distinguished from the turbulent shear layer by its growth rate; the laminar shear layer exiting the nozzle has a significantly lower spreading rate untill $x\approx 0.8$, where its growth 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 panels ($c$)-($g$). At the nozzle exit in figure 1 ($c$), the RMS profile of the turbulent jet peaks at $\gtrsim 10\%$ of the freestream 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 will later show that the rapid growth of the shear layer is associated with exponential growth of the hydrodynamics instability modes supported by the initially laminar shear layer. We explore the nonlinear dynamics of the jet that lead this rapid growth and explain previous observations by (Zaman & Hussain, 1980; Kim & Choi, 2009; Bogey & Bailly, 2010).

Figure 2: Experimental validation of turbulent jet LES and comparison of initially-laminar jet LES with literature (Zaman & Hussain, 1980; Bogey & Bailly, 2010): ($a$) mean and ($b$) RMS of the streamwise velocity on the centerline. 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, turbulent jets, LES, and experiments. The mean and RMS streamwise velocities on the centerline are shown in figure 2 ($a$) and ($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 centerline, 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 (2010) 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 (1980) and Bogey & Bailly (2010). This hump and elevated RMS were previously associated with vortex-pairing by Kim & Choi (2009); Bogey & Bailly (2010); Bogey et al. (2012). The authors also demonstrated that both these phenomena strongly depend on the initial shear-layer thickness, with the 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 centerline (Bogey & Bailly, 2010). In the remainder of this paper, we go beyond this phenomenological description and analyze the underlying nonlinear mechanism in detail. To this end, we use spectral modal decomposition techniques and local linear theory.

3.2 Shear layer instability

Figure 3: Premultiplied radially integrated PSD along $x$ for initially laminar ($a$,$b$) and turbulent jets ($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 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, $\mbox{\it 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 will later confirm this using local stability theory 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 $\mbox{\it St}_{\theta}=f\theta/U_{j}=0.0132$, which closely matches that of the shear layer mode observed in forced jet experiments Zaman & Hussain (1980). For all three frequencies, the PSD of the $u_{r}$ component is greater than the $u_{x}$ component. The PSDs for the turbulent jet in figure 3 ($c$,$d$) are significantly lower than the initially-laminar jet. For the turbulent jet, the first tonal component is observed at $x\approx 2.6$ and $\mbox{\it 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 are absent in the fully turbulent jet.

Figure 4: Total RMS velocities and selected contributing frequency components along $x$: ($a$) $u_{x}^{rms}$ on the centerline, $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 $\mbox{\it 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 centerline ($r=0$) in ($a$), the streamwise velocity on the lipline ($r=0.5$) in ($b$), and the radial velocity on the lipline in ($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 centerline, 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 to 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 the centerline. 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 are in agreement with the research of Hajj et al. (1992), which demonstrates that the simultaneous decay of the fundamental frequency and growth of the subharmonic is a result of the energy transfer from the fundamental to the subharmonic through a resonance mechanism (Monkewitz, 1988). Later in §3.4, we will investigate this energy transfer using the spectral kinetic energy equation.

Figure 5: Comparison of TKE and amplification rate predicted from linear stability theory 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 equation (19), for the initially-laminar jet ($c,e$) and turbulent jet ($d,f$). The solid and dashed lines in ($a,b$) represent the TKE computed from data and LST using equation (20), 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$).

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 deducted from data and theory. We define the local amplitude of $\hat{\mathbf{q}}$ as

$$A(x,m,\omega)=\sqrt{\int_{r}\hat{\mathbf{q}}^{*}\hat{\mathbf{q}}r\differential r}=\sqrt{k^{\text{LES}}},$$ (18)

where $k^{\text{LES}}$ is the turbulent kinetic energy computed from the data. Using equation (18) and the normal mode ansatz, equation (1), allows us to compute the empirical growth rate as

$$\alpha_{i}^{\rm{emp}}=-\frac{1}{A}\frac{\differential A}{\differential x}.$$ (19)

The TKE can also be predicted from the local linear theory by integrating and squaring the theoretical amplification rate

$$k^{\rm{LST}}=\pqty{\int\alpha(x)dx}^{2}.$$ (20)

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 normalize $k^{\rm{LST}}$ by the maximum $k^{\text{LES}}$ of each frequency.

Figure 5 shows the TKE ($a,b$) and spatial growth rate ($c$-$f$) for the axisymmetric component $m=0$ in both jets. The spatial growth rate is computed from local linear theory in panels ($c,d$) and empirically using equation (19) in panels ($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 $\mbox{\it St}=1.76$, 0.88, and 0.44 in the initially laminar jet and for $\mbox{\it St}=0.56$ in the turbulent jet, is shown in figure 5($a$) and ($b$), respectively. The TKE estimated from the linear stability theory ($\text{k}^{\rm{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 figure 5($c$) and ($d$), which show the spatial growth rates of the initially laminar and turbulent jets, respectively. Figure 5($c$) reveals that $\mbox{\it 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 $\mbox{\it St}=1.76$ is the most unstable frequency and occurs most upstream at the nozzle’s exit confirms our previous interpretation of $\mbox{\it St}\approx 1.76$ as the fundamental frequency. The absence of a peak at $\mbox{\it St}=0.88$ and 0.44 indicates that these frequencies are not a result of 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 $\mbox{\it St}\approx 0.90$. The empirical spatial growth rates in $(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, and second, it supports the physical interpretation that the observed amplitude disturbances of individual frequency components represent hydrodynamic waves.

Figure 6: SPOD eigenspectra with a focus on the shear layer untill 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 ultra harmonics $\mbox{\it St}=0.66$ and 1.32, and the dotted line in ($b$) to the most energetic frequency of the turbulent jet.

As both jets are statistically stationary flows, we use SPOD to extract the spatiotemporal coherent structures. To emphasize 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 $\mbox{\it 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 $\mbox{\it St}=1.76$, which has lower energy compared to the two subharmonics. Additionally, the spectra reveal the presence of ultra and superharmonics at $\mbox{\it St}\approx 2.7$ and 3.5, respectively. We have confirmed from the spatial linear theory that $\mbox{\it 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 $\mbox{\it St}=0.66$ and $1.32$ (denoted by gray dotted lines). Later, in figure 12, we will 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\leq\mbox{\it St}\leq 1$ (Schmidt et al., 2018). Overall, figure 6 indicates that the presence of discrete tones leads to a much higher fluctuation level than the turbulent case.

Figure 7: Leading SPOD modes of the fundamental and four of its subharmonic frequencies: ($a$,$b$) $\mbox{\it St}=1.76$; ($c$,$d$) $\mbox{\it St}=0.88$; ($e$,$f$) $\mbox{\it St}=0.44$; ($g$,$h$) $\mbox{\it St}=0.22$; ($i$,$j$) $\mbox{\it St}=0.11$. The left column represents the streamwise velocity component, $u_{x}$ and the right 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 $\mbox{\it 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 materializes as a Kelvin-Helmholtz wavepacket with compact support in the region $0.8\leq x\leq 1.2$ around the lip line. At a lower frequency of $\mbox{\it St}=0.44$, the leading SPOD mode exhibits larger spatial support within the range $1\leq x\leq 6$ and $0\leq r\leq 1$. Notably, it exhibits high amplitude in $2\leq x\leq 3$, where the streamwise velocity is concentrated on the centerline 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 centerline 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 $\mbox{\it St}=0.88$ is at $(x,r)\approx(1.5,0.5)$, the corresponding to the peak of the blue curve in figure 4 ($c$).

3.3 Vortex pairing

Figure 8: Time traces exemplifying two successive vortex pairing events are visualized 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 visualized 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 $\mbox{\it 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 panels 8($b$-$d$). This vortex pairing results in the formation of the $\mbox{\it St}=0.88$ vortex. Another instance of vortex pairing is evident in panels 8($c$-$e$). Next, the two $\mbox{\it St}=0.88$ vortices, denoted by blue rectangles in ($e$), undergo pairing to form the $\mbox{\it St}=0.44$ vortex. Panels 8($e$-$f$) demonstrate 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.

Figure 9: $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\leq t\leq 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 (2021). 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 realization 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 realizations. 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$.

3.4 SPOD-based spectral energy transfer

We now investigate the spectral energy balance of turbulent kinetic energy 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. The nonlinear energy transfer term could alternatively be estimated using bispectral mode decomposition (BMD, Schmidt, 2020). 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 $\mbox{\it 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 equation (14). Positive production indicates the energy gain from the mean flow, and negative production represents energy loss to the mean flow. Figure 10($a$) shows that $\mbox{\it St}=1.76$ is the earliest to gain energy from the mean flow, and as it saturates, $\mbox{\it St}=0.88$ begins its growth. Subsequently, $\mbox{\it St}=0.88$ attains its global maximum at $x=1.25$, which also corresponds to the location of the global minimum of $\mbox{\it St}=1.76$. Similarly, the global maximum of $\mbox{\it St}=0.44$ and the global minimum of $\mbox{\it 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 (1982) and was confirmed by Hajj et al. (1992, 1993). Our findings are in agreement with these 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 $\mbox{\it 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 linear stability theory (denoted by dashed lines, also see figure 5). The correspondence is expected as the positive production is associated with the amplification rate, and the negative production is associated with the decay rate of each frequency.

Figure 10: Production and nonlinear energy transfer terms for $\mbox{\it 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 the linear stability theory, of the corresponding frequency. The black dotted line indicates the onset of nonlinear interactions.

The net nonlinear energy transfer, computed using equation (17), is shown in figure 10($b$). The net nonlinear energy transfer is always negative for $\mbox{\it St}=1.76$. On the other hand, for $\mbox{\it 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\leq x\leq 2.3$ exhibit a similar shape but are opposite in sign, indicating that nonlinear interactions result in energy loss for $\mbox{\it St}=0.88$ and energy gain for $\mbox{\it St}=0.44$. This phenomenon is expected for a vortex pairing process, as the energy is transferred from $\mbox{\it St}=0.88$ to its subharmonic, $\mbox{\it 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 $\mbox{\it St}=0.88$ starts to grow (see figure 10($a$)), implying that the nonlinear interactions are triggered due to 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.

Figure 11: Spatial fields of production ($a$,$d$,$g$), dissipation ($b$,$e$,$h$), and net nonlinear energy transfer ($c$,$f$,$i$).

Figure 11 shows the spatial fields of production, dissipation, and net nonlinear transfer at $\mbox{\it St}=1.76$, 0.88, and 0.44. These fields are computed using equations (14), (15), and (17), respectively. As in figure 10, the production fields in panels ($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 localized in the region about the lip line. 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 panels ($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.

Figure 12: Nonlinear energy transfer using SPOD: transfer term bispectrum for ($a$) entire domain and ($b$) shear layer subdomain, $\Omega_{r}$; spatial fields for the triads ($c$) (1.76,-0.88,0.88) and ($d$) (0.88,-0.44,0.44) are compared to the radially integrated TKE of $\mbox{\it St}=1.76,0.88$, and 0.44.

Figure 12($a$) and ($b$) show the triadic energy transfer obtained by integrating equation (16) over the entire spatial domain and a domain $\mathrm{\Omega}_{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 panels ($a$) and ($b$), positive values (red color) denote energy transfer to frequency $\mbox{\it St}_{3}$ from frequencies $\mbox{\it St}_{1}$ and $\mbox{\it St}_{2}$, while negative values (blue color) indicate the extraction of energy from $\mbox{\it St}_{3}$ by $\mbox{\it St}_{1}$ and $\mbox{\it 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 from $\mbox{\it St}=0.88$ and -0.44 to $\mbox{\it St}=0.44$. On the other hand, the triad $(0.44,0.44,0.88)$ has a negative value, indicating that $\mbox{\it St}=0.44$ extracts energy from $\mbox{\it 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 $\mbox{\it St}=1.32$. The occurrence of $\mbox{\it 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 $\mbox{\it St}=1.32$. The nonlinear energy transfer analysis reveals that the ultraharmonic $\mbox{\it 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 $\mbox{\it St}=1.76$ is dynamically significant but not energetically dominant. Therefore, to shed light on the nonlinear interactions of the fundamental frequency, we narrow our focus to the first two jet diameters. Figure 12($b$) now illuminates the triads (1.76,-0.88,0.88) and (0.88,0.88,1.76). These triads show a similar energy transfer behavior characteristic of the vortex-pairing process. 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 ($\mbox{\it St}=1.76$) to its first subharmonic ($\mbox{\it St}=0.88$), and then from the first subharmonic to the second subharmonic ($\mbox{\it St}=0.44$). Our results hence support the hypothesis of a parametric resonance mechanism proposed by Monkewitz (1988).

The spatial fields of $\mathscr{t}_{\text{nl}}^{\text{rank-1}}$ for the dynamically dominant and energetically significant triad are shown in figure 12($c$) and ($d$), respectively. The TKE of the fundamental, subharmonic, and second subharmonic frequencies are overlaid on the contours of the $\mathscr{t}_{\text{nl}}^{\text{rank-1}}$ field. For the (1.76, -0.88, 0.88) triad, the $\mathscr{t}_{\text{nl}}^{\text{rank-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 panel ($d$), demonstrating that the nonlinear energy transfer of the triad (0.88, -0.44, 0.44) is localized to the region of subharmonic decay and second subharmonic growth. In summary, SPOD-based transfer analysis allowed us to systematically catalog 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 focussed on the hydrodynamics instabilities (Michalke, 1965; Kelly, 1967), phase locking of vortices (Monkewitz, 1988; Husain & Hussain, 1995), and specific aspects of energy transfer (Mankbadi, 1985; Paschereit et al., 1995) during the vortex pairing process. This study reveals the complete physical picture of the vortex-pairing process, establishes a clear concept for identifying the fundamental frequency, and provides a framework to characterize the energy transfer between the most energetic coherent structures.

The analysis was conducted for two LES of initially laminar and turbulent jets at $\Rey=50000$. 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, its shear layer’s growth rate is more gradual compared to that of the initially-laminar jet. Local linear theory identifies the fundamental frequency, i.e., the frequency with the largest spatial growth rate, as $\mbox{\it St}=1.76$ for the initially-laminar jet. At this frequency, the shear layer rolls up into vortices. The tones at $\mbox{\it 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 of the fundamental frequency. 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 utilized to estimate spatial growth rates in turbulent flows. Conversely, linear stability theory can provide an accurate prediction of growth rates of coherent structures despite the presence of nonlinear dynamics. The growth rate spectrum estimated empirically is also broadband, which re-emphasizes that the peaks in PSD at first and second subharmonics are due to nonlinear effects.

Different visualization 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 clearly inferred from LST. 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 them mainly to its subharmonic frequencies, which result from nonlinear interactions. This is a manifestation of the importance of distinguishing between dynamical relevance and energetic significance, as highlighted, among others, by Schmid (2010).

SPOD-based spectral TKE analysis was performed to characterize 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 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. In a spatial sense, as two fundamental vortices merge, there is a backscatter of energy from the fundamental to its subharmonic. This process repeats itself for higher subharmonics. Few studies (Hajj et al., 1992, 1993) argue that the dominant interaction during vortex pairing is between the subharmonic and the fundamental frequency, while other studies (Mankbadi, 1985; Paschereit et al., 1995) argue that the subharmonic gains most of its energy from the mean flow and the fundamental-subharmonic interaction only acts as a catalyst. Our findings show that both views are correct in their own right. 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.

The entire process can be summarized in the following steps:

1. The hydrodynamic instabilities initiate the transition into turbulence, causing the shear layer to grow rapidly.

2. 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.

3. The saturation and subsequent decay of the fundamental frequency mark the onset of vortex pairing.

4. 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.

5. Processes analogous to steps (ii) and (iii) then repeat to create higher subharmonics.

Appendix A Nonlinear energy transfer using SPOD and BMD

Here, we qualitatively compare the nonlinear energy transfer estimated from SPOD and BMD. BMD is 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 the triadic interactions by maximizing the integrated point-wise bispectrum

$$b(f_{k},f_{l})=E\bqty{\int_{\Omega}\hat{\mathbf{q}}_{1}^{*}(\mathbf{x},f_{k})\circ\hat{\mathbf{q}}_{2}^{*}(\mathbf{x},f_{l})\circ\hat{\mathbf{q}}_{3}(\mathbf{x},f_{k}+f_{l})d\mathbf{x}}.$$ (21)

For further details on computing the mode bispectrum, the reader is referred to Schmidt (2020). We tailor BMD to estimate the nonlinear energy transfer by maximizing the following point-wise integral

$$b_{\rm{nl}}^{\rm{BMD}}(f_{k},f_{l})=E\bqty{\int_{\Omega}\hat{\mathbf{u}}_{i}^{*}(\mathbf{x},f_{k})\circ\frac{\partial\hat{\mathbf{u}}_{j}}{\partial\mathbf{x}_{i}}^{*}(\mathbf{x},f_{l})\circ\hat{\mathbf{u}}_{j}(\mathbf{x},f_{k}+f_{l})d\mathbf{x}}.$$ (22)

Figure 13: Comparison of nonlinear energy transfer: ($a$) SPOD; ($b$) BMD.

Figure 13 shows the nonlinear energy transfer estimated using SPOD (equation 16) and BMD (equation 22). Qualitatively similar trends are observed, particularly, the direction of energy transfer is 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 to that of the SPOD-based $\mathscr{t}_{\rm{nl}}$, this is because the BMD in comparison to SPOD requires more data for convergence.