1. Introduction
Jet noise is an environmental challenge that the aerospace industry must face today and for years to come. High-speed aircraft and space launchers commonly feature multi-jet engines to supply the required levels of thrust, with twin-jet configurations being the most common. Twin-jet flow fields present a number of challenges in comparison to single jets, which have been widely studied in the past. When close to each other, the jets may interact at the hydrodynamic and acoustic levels, yielding complex physical mechanisms that remain to be understood. Early experiments studying twin-jet configurations (Bhat Reference Bhat1977; Kantola Reference Kantola1981) identified mechanisms that result in a reduction of the far-field noise with respect to an equivalent isolated jet, due in part to an acoustic shielding of one jet by the other. In contrast, recent experimental observations by Bozak & Henderson (Reference Bozak and Henderson2011) have revealed enhanced noise levels in the plane contained between the two jets, which are above those corresponding to two linearly superposed incoherent jets.
The development of jet-noise reduction and control technologies requires understanding and appropriate simplified models for the mechanisms of sound generation. This has motivated study of the dynamics of turbulent shear flow. In this regard, the importance of ordered motions (today known as coherent structures) in turbulent planar mixing layers and jets has been recognised since the pioneering experimental observations by Mollo-Christensen (Reference Mollo-Christensen1967), Crow & Champagne (Reference Crow and Champagne1971) and Brown & Roshko (Reference Brown and Roshko1974). Since then, the connection between the sound radiated by high-speed jets and the coherent structures has been the focus of significant research, as reviewed by Jordan & Colonius (Reference Jordan and Colonius2013). Coherent structures were found to be reminiscent of linear instability waves growing in harmonically forced supersonic jets, a finding that prompted the use of linear stability theory for their modelling (Crow & Champagne Reference Crow and Champagne1971; Crighton & Gaster Reference Crighton and Gaster1976; Michalke Reference Michalke1984; Tam & Morris Reference Tam and Morris1985) and established the term ‘wavepacket’ to refer to them. Since those studies, the presence of wavepackets in unforced high-speed jets and their relationship with mixing-noise emission has been demonstrated both from experimental analyses (Juvé et al. Reference Juvé, Sunyach and Compte-Bellot1980; Guj et al. Reference Guj, Carley, Camussi and Ragni2003) and large-eddy simulations (Cavalieri et al. Reference Cavalieri, Daviller, Comte, Jordan, Tadmor and Gervais2011), together with the finding that the radiated sound is highly directional both for subsonic and perfectly expanded supersonic jets (Crighton & Huerre Reference Crighton and Huerre1990; Tam Reference Tam1995; Cavalieri et al. Reference Cavalieri, Jordan, Colonius and Gervais2012). In parallel, several works consolidated the ability of linear stability calculations to model wavepackets correctly, both in subsonic (Suzuki & Colonius Reference Suzuki and Colonius2006; Piot et al. Reference Piot, Casalis, Muller and Bailly2006; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Cavalieri et al. Reference Cavalieri, Rodríguez, Jordan, Colonius and Gervais2013) as well as supersonic jets (Sinha et al. Reference Sinha, Rodríguez, Brès and Colonius2014). Although wavepackets constitute a relatively small fraction of the total fluctuation energy of the turbulent flow (Cavalieri et al. Reference Cavalieri, Rodríguez, Jordan, Colonius and Gervais2013), they nevertheless are known to be acoustically significant owing to their high spatio-temporal coherence with respect to the small turbulent scales (Jordan & Colonius Reference Jordan and Colonius2013).
Among the different linear stability theories employed for wavepacket modelling in single jets, the parabolised stability equations (PSE) have been broadly used. PSE were found to deliver a satisfactory agreement with both high-fidelity simulations and experiments for subsonic and supersonic jets at a small computational cost (Yen & Messersmith Reference Yen and Messersmith1999; Piot et al. Reference Piot, Casalis, Muller and Bailly2006; Ray, Cheung & Lele Reference Ray, Cheung and Lele2009; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Rodríguez et al. Reference Rodríguez, Sinha, Brès and Colonius2012; Cavalieri et al. Reference Cavalieri, Rodríguez, Jordan, Colonius and Gervais2013; Rodríguez et al. Reference Rodríguez, Sinha, Brès and Colonius2013; Sinha et al. Reference Sinha, Rodríguez, Brès and Colonius2014; Breakey et al. Reference Breakey, Jordan, Cavalieri, Nogueira, Léon, Colonius and Rodríguez2017; Sasaki et al. Reference Sasaki, Cavalieri, Jordan, Schmidt, Colonius and Brès2017). These works have also served to identify and address some limitations of the PSE model for jet noise, which may be overcome by more computationally expensive approaches such as the one-way Navier–Stokes equations (Towne & Colonius Reference Towne and Colonius2015; Towne et al. Reference Towne, Rigas, Kamal, Pickering and Colonius2022) or resolvent analysis (Garnaud et al. Reference Garnaud, Lesshafft, Schmid and Huerre2013; Jeun, Nichols & Jovanović Reference Jeun, Nichols and Jovanović2016; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Cavalieri, Jordan & Lesshafft Reference Cavalieri, Jordan and Lesshafft2019). On the one hand, difficulties have been found in predicting the wavepacket properties at low frequencies and low azimuthal wavenumbers, conditions at which recent studies based on resolvent analysis have found non-modal growth to be important through the Orr and lift-up mechanisms (Tissot et al. Reference Tissot, Lajús, Cavalieri and Jordan2017; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Lesshafft et al. Reference Lesshafft, Semeraro, Jaunet, Cavalieri and Jordan2019; Nogueira et al. Reference Nogueira, Cavalieri, Jordan and Jaunet2019; Pickering et al. Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020). Similarly, discrepancies have been encountered in the PSE description of wavepackets downstream of the region at which the growth of Kelvin–Helmholtz instabilities saturates, where locally parallel transient growth calculations and resolvent analyses reveal non-modal phenomena to be important and point to the role of nonlinear interactions in their activation (Jordan et al. Reference Jordan, Zhang, Lehnasch and Cavalieri2017; Tissot et al. Reference Tissot, Lajús, Cavalieri and Jordan2017; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Lesshafft et al. Reference Lesshafft, Semeraro, Jaunet, Cavalieri and Jordan2019). On the other hand, an underprediction of the amplitude of Mach-wave radiation has been observed, which can be attributed to the inherent limitations of the regularised marching problem, as described by Towne, Rigas & Colonius (Reference Towne, Rigas and Colonius2019). Despite these shortcomings, coherent structures in supersonic jets at the frequency range of interest in the study of mixing noise are dominated by Kelvin–Helmholtz instabilities (Pickering et al. Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020), which are well modelled by PSE owing to their modal convective nature. For supersonic mixing layers featuring Kelvin–Helmholtz waves with high (supersonic) phase speed, the error in the PSE representation of Mach-wave radiation is found to be small (Towne et al. Reference Towne, Rigas and Colonius2019). These conditions are satisfied for most wavepacket calculations considered in this work.
Compared with single round jets, studies of wavepackets and their modelling in twin-jet systems are scarce in the literature, mainly due to the increased complexity of the flow field. The twin-jet mean flow is no longer axisymmetric, which prevents the introduction of azimuthal Fourier modes in both the modelling approaches and the experimental postprocessing techniques, and in turn requires inhomogeneity in at least two spatial directions. From the modelling point of view, Sedel’nikov (Reference Sedel’nikov1967) was the first to derive dispersion relationships for multi-jet configurations employing vortex sheets, although no solutions were determined. After this, only three works addressed the parallel-flow linear instability problem for twin jets before the last decade (Morris Reference Morris1990; Du Reference Du1993; Green & Crighton Reference Green and Crighton1997), making use of vortex sheet and finite-thickness models based on bipolar coordinate systems to study the inviscid instability of two axially homogeneous parallel jets. Morris (Reference Morris1990) established the classification of twin-jet fluctuations into four possible families according to the natural symmetries of the system. More recently, cross-plane local linear stability theory (Rodríguez et al. Reference Rodríguez, Jotkar and Gennaro2018; Nogueira & Edgington-Mitchell Reference Nogueira and Edgington-Mitchell2021; Rodríguez et al. Reference Rodríguez, Stavropoulos, Nogueira, Edgington-Mitchell and Jordan2023), along with plane-marching parabolised stability equations (PM-PSE) (Rodríguez et al. Reference Rodríguez, Jotkar and Gennaro2018; Rodríguez Reference Rodríguez2021) have enabled characterisation of the three-dimensional structure of the Kelvin–Helmholtz instabilities associated with mixing noise in the round twin-jet system. Rodríguez et al. (Reference Rodríguez, Jotkar and Gennaro2018) computed wavepacket models for subsonic twin jets employing local stability analyses and plane-marching PSE on a tailored twin-jet mean flow, consisting of the linear superposition of two single-jet analytical velocity fields fitted from hot-wire measurements. Similarly, Rodríguez (Reference Rodríguez2021) obtained wavepacket models for perfectly expanded supersonic twin jets via PM-PSE, once again using a tailored twin-jet mean flow fitted from large eddy simulation (LES) calculations. Nogueira & Edgington-Mitchell (Reference Nogueira and Edgington-Mitchell2021) characterised Kelvin–Helmholtz instabilities in a supersonic twin-jet system operating at underexpanded conditions via a locally parallel Floquet stability analysis, featuring a twin-jet mean flow obtained from planar particle image velocimetry (PIV) measurements revolved around the nozzle axes. Lastly, Rodríguez et al. (Reference Rodríguez, Stavropoulos, Nogueira, Edgington-Mitchell and Jordan2023) revisited the locally parallel stability of perfectly expanded twin jets using a vortex sheet and a finite-thickness model based on a tailored analytical mean flow, revealing that the coupling between the fluctuation fields of the two jets favours flapping motions over helical ones. These works employed simplified mean-flow models that account for linear interaction between the two jets. However, experimental mean-flow measurements and high-fidelity simulations indicate that the interaction is nonlinear, particularly in the case of closely spaced jets.
The aforementioned modelling efforts have predicted wavepackets analogous to those in single round jets, although with observable deviations in their spatial structure from the exact azimuthal Fourier modes characteristic of single round jets. The few experimental works aimed at their identification and characterisation in round twin-jet systems are very recent and focused on conditions at which screech is dominant (Kuo, Cluts & Samimy Reference Kuo, Cluts and Samimy2017; Knast et al. Reference Knast, Bell, Wong, Leb, Soria, Honnery and Edgington-Mitchell2018; Bell et al. Reference Bell, Cluts, Samimy, Soria and Edgington-Mitchell2021; Nogueira & Edgington-Mitchell Reference Nogueira and Edgington-Mitchell2021; Nogueira, Stavropoulos & Edgington-Mitchell Reference Nogueira, Stavropoulos and Edgington-Mitchell2021; Wong et al. Reference Wong, Stavropoulos, Beekman, Towne, Nogueira, Weightman and Edgington-Mitchell2023). Therefore, while the previously described investigations have provided physically sound wavepacket models for mixing noise, an assessment of their validity is still missing. Significant work has also focused on the study of the coupling and control of twin rectangular jets under screech conditions (Raman, Panickar & Chelliah Reference Raman, Panickar and Chelliah2012; Esfahani, Webb & Samimy Reference Esfahani, Webb and Samimy2021; Yeung, Schmidt & Brès Reference Yeung, Schmidt and Brès2022; Samimy et al. Reference Samimy, Webb, Esfahani and Leahy2023; Jeun, Wu & Lele Reference Jeun, Wu and Lele2024; Karnam et al. Reference Karnam, Ahn, Mihaescu, Saleem and Gutmark2025). However, no quantitative comparison between wavepacket models and experimental measurements currently exist for such rectangular configurations either.
A reduced number of recent works employed LES to study the flow dynamics in round twin jet configurations. Brès et al. (Reference Brès, Bose, Ham and Lele2014) and Gao, Xu & Li (Reference Gao, Xu and Li2016) simulated twin jets at under-expanded conditions with exit Mach number 1.36 and centre-to-centre spacing of
$s=2.625D$
and
$2.25D$
, respectively; however, Brès et al. (Reference Brès, Bose, Ham and Lele2014) considered heated jets and fully turbulent boundary layers, while Gao et al. (Reference Gao, Xu and Li2016) simulated an isothermal jet with initially laminar shear layers. These differences lead to clear impact on the turbulent flow dynamics and radiated noise. Goparaju & Gaitonde (Reference Goparaju and Gaitonde2018) simulated a twin-jet system with
$s = 2D$
at perfectly expanded conditions, Mach number 1.23 and initially laminar shear layers. Finally, Muthichur et al. (Reference Muthichur, Vempati, Hemchandra and Samanta2023) simulated two configurations, with
$s = 1.1D$
and
$2D$
, Mach number 0.9, comparatively low Reynolds number and initially laminar shear layers. Similar to the case of experiments, the relative scarcity and sparsity of LES data for supersonic round twin jets has precluded a complete validation of the wavepacket models based on linear calculations.
This work contributes to the modelling of wavepackets in supersonic round twin jets in two new ways: (i) by employing a more accurate mean-flow representation based on the compressible RANS equations, which is used as an input for the plane-marching parabolised stability equations, and (ii) by performing the first quantitative comparisons of twin-jet wavepackets against coherent structures educed from experimental measurements. The three-dimensional RANS flows account for the nonlinear interactions between the mean flows of the two jets. By using them to inform the plane-marching PSE, the impact of the jet-to-jet mean flow interactions is incorporated in the wavepacket predictions. Here, RANS solutions are computed for perfectly expanded twin jets and validated by means of PIV measurements of the mean flow.
The validation of twin-jet wavepackets computed by PSE against experimental data is realised by application of spectral proper orthogonal decomposition (SPOD) to high-speed schlieren visualisations. Recent investigations focused on the study of screech in single (Edgington-Mitchell et al. Reference Edgington-Mitchell, Li, Liu, He, Wong, MacKenzie and Nogueira2022; Karnam, Saleem & Gutmark Reference Karnam, Saleem and Gutmark2023) and twin jets (Esfahani et al. Reference Esfahani, Webb and Samimy2021; Nogueira et al. Reference Nogueira, Stavropoulos and Edgington-Mitchell2021; Prasad et al. Reference Prasad, Gaitonde, Esfahani, Webb and Samimy2022; Wong et al. Reference Wong, Stavropoulos, Beekman, Towne, Nogueira, Weightman and Edgington-Mitchell2023; Karnam et al. Reference Karnam, Ahn, Mihaescu, Saleem and Gutmark2025) have educed organised structures at resonant frequencies by means of SPOD applied to schlieren measurements. In this case, the dynamics of the resonant mechanism feature a strong low-rank behaviour and a highly organised structure that is contained in the most energetic SPOD mode. The eduction of coherent structures from schlieren visualisations in perfectly expanded twin-jet systems is more challenging, as the schlieren images under non-screeching conditions are dominated by highly energetic, high-azimuthal wavenumber, small-scale vortical structures localised in the turbulent shear layer (Cavalieri et al. Reference Cavalieri, Rodríguez, Jordan, Colonius and Gervais2013). These structures dominate the fluctuation energy of the turbulence and complicate the use of techniques like proper orthogonal decomposition in obtaining a low-rank representation of wavepackets.
In this work, recent improvements (Prasad & Gaitonde Reference Prasad and Gaitonde2022; Padilla-Montero et al. Reference Padilla-Montero, Rodríguez, Jaunet and Jordan2024) that facilitate the extraction of mixing-noise-related coherent structures from schlieren images are adopted. These consist in performing spectral proper orthogonal decomposition on a filtered quantity derived from the schlieren images, instead of applying it to the schlieren fields directly. The quantity used is intimately related to the irrotational momentum potential field introduced by Doak (Reference Doak1989), which does not feature small-scale vortical fluctuations present in the schlieren visualisations. Through this methodology, coherent structures are extracted from the experimental datasets in the form of SPOD modes which can be compared with the PM-PSE wavepackets. To quantify the agreement between both, an alignment metric based on the projection of one SPOD field into one PM-PSE fluctuation field is proposed and evaluated for different frequencies and jet separations. Due to the line-of-sight integration implicit in the schlieren visualisations, only those fluctuations which are symmetric with respect to the plane containing both jets can be realised from the employed experimental set-up. The analysis and validation of wavepackets which are antisymmetric with respect to this plane are thus not considered in this work.
The remainder of this paper is organised as follows. Section 2 describes the twin-jet geometry and conditions considered for study. Section 3 presents the adopted modelling strategy for twin-jet wavepackets, formulating the governing equations as well as providing details on the numerical methodology employed for the computations. In § 4, the twin-jet experimental set-up is reported, including a description of the method employed for extracting coherent structures from the experimental data. Section 5 presents the obtained mean-flow solutions and wavepacket models. Particular characteristics of the twin-jet fluctuations for closely spaced jets are discussed, and their comparison against the experimentally educed structures is reviewed. Finally, concluding remarks are summarised in § 6.
2. Twin-jet configuration
Figure 1. Sketch representing the twin-jet configuration and the associated geometrical parameters:
$(a)$
cross-stream plane (constant
$x$
);
$(b)$
streamwise plane containing both jets axes at
$z = 0$
.
The twin-jet configuration considered is sketched in figure 1. The jets are generated by two round convergent–divergent nozzles. The centre of each of the nozzles is located along the
$y$
-axis, with
$y = 0$
denoting the symmetry plane between the jets. Each nozzle exit is located at
$x = 0$
and the spacing between the nozzles axes is denoted by
$s$
. The interior nozzle geometry follows a truncated ideal contour (TIC) profile with an exit diameter of
$D = 0.025$
m and an exit-to-throat area ratio of
$A_e / A_t = 1.225$
. The nozzle geometry is shown in figure 3.
The jets are generated under a fixed nozzle pressure ratio (NPR) corresponding to the perfectly expanded condition and a fixed total temperature of
$T_{0} = 300$
K. The nozzle pressure ratio, defined as the ratio of the total pressure in the reservoir
$p_{0}$
to the ambient pressure
$p_\infty$
, is expressed in terms of the isentropic jet-exit Mach number
$M_{\kern-1pt j}$
through the isentropic relation
$p_{0} / p_\infty = [1 + 0.5(\gamma - 1) M_{\kern-1pt j}^2 ]^{\gamma /(\gamma - 1)}$
, where
$\gamma = 1.4$
denotes the ratio of specific heats. The jets are unheated and the flow acceleration within the nozzles results in jet static temperatures lower than the ambient temperature. The operating condition closest to the perfectly expanded regime has been calibrated experimentally by means of flow visualisation and corresponds to
$M_{\kern-1pt j} = 1.54$
(NPR = 3.89). This operating condition is considered throughout the study. Two different nozzle separations are investigated:
$s/D = 1.76$
and
$s/D = 3$
. The first corresponds to the minimum distance allowed between the nozzles’ axes due to the size of the outer nozzle diameter. A strong interaction between jets is expected for
$s/D = 1.76$
, while a weak interaction is expected for
$s/D = 3$
. The system exhibits two symmetry planes: the plane that contains both jets (the
$xy$
plane located at
$z = 0$
) and the plane between the jets (the
$xz$
plane located at
$y = 0$
).
3. Modelling strategy
3.1. Mean flow
In the past, twin-jet mean-flow descriptions have been obtained using large eddy simulations (Brès et al. Reference Brès, Bose, Ham and Lele2014; Gao et al. Reference Gao, Xu and Li2016; Goparaju & Gaitonde Reference Goparaju and Gaitonde2018; Brès et al. Reference Brès, Yeung, Schmidt, Isfahani, Webb, Samimy and Colonius2021; Troyes & Vuillot Reference Troyes and Vuillot2022; Muthichur et al. Reference Muthichur, Vempati, Hemchandra and Samanta2023) or simplified approaches such as the axisymmetric revolution of twin-jet planar PIV measurements (Nogueira & Edgington-Mitchell Reference Nogueira and Edgington-Mitchell2021) or via a tailored twin-jet model, in which the twin-jet flow field is constructed by linearly superposing the velocity fields of two single jets. In the latter, the single-jet mean fields may come either from analytical models such as hyperbolic tangent profiles (Stavropoulos et al. Reference Stavropoulos, Mancinelli, Jordan, Jaunet, Weightman, Edgington-Mitchell and Nogueira2023; Rodríguez et al. Reference Rodríguez, Stavropoulos, Nogueira, Edgington-Mitchell and Jordan2023) or analytical Gaussian fittings to single-jet PIV measurements (Gudmundsson & Colonius Reference Gudmundsson and Colonius2011), hot-wire measurements (Rodríguez et al. Reference Rodríguez, Jotkar and Gennaro2018) or LES calculations (Rodríguez Reference Rodríguez2021).
One of the goals of this work is to model wavepackets in supersonic twin-jets by means of linear parabolised stability equations based on an accurate mean-flow representation. A mean-flow description that properly accounts for the nonlinear interactions between the jets at a mean-flow level is thus preferred. A LES calculation of the twin-jet configuration would provide one of the most accurate mean-flow predictions, as the actual jet dynamics and their mutual interactions are captured explicitly by LES and their imprint in the mean flow is retained by the time averaging operation. However, the computational cost of a properly resolved LES of the flows considered in this work, comprising complex geometries and very high Reynolds numbers, are considerable: of the order of thousands of cores and millions of CPU hours in a supercomputer (Brès et al. Reference Brès, Bose, Ham and Lele2014). Instead, three-dimensional compressible RANS calculations are used here to compute the mean flows, including the mean-flow jet interaction, with a good accuracy at a significantly reduced computational cost. The accuracy of the RANS solutions is assessed in § 5.1 by comparison to PIV measurements.
3.1.1. Governing equations
The twin-jet mean flow is modelled using the compressible RANS equations (also known as Favre-averaged Navier–Stokes equations), see e.g. Wilcox (Reference Wilcox2006):
\begin{align} \frac {\partial \bar {\rho }}{\partial t} + \frac {\partial }{\partial x_i} \left ( \bar {\rho } \tilde {u}_i \right ) = 0, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \left ( \bar {\rho } \tilde {u}_i \right )}{\partial t} + \frac {\partial }{\partial x_j} \left ( \bar {\rho } \tilde {u}_i \tilde {u}_j \right ) = -\frac {\partial \bar {p}}{\partial x_i} + \frac {\partial }{\partial x_j} \left ( \bar {\tau }_{\textit{ij}} - \overline {\rho u_i'' u_{\kern-1pt j}''} \right )\!, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial }{\partial t} \left ( \bar {\rho } \tilde {E} + \frac {\overline {\rho u_i'' u_i''}}{2} \right ) & + \frac {\partial }{\partial x_j} \left ( \bar {\rho } \tilde {u}_j \tilde {H} + \tilde {u}_j \frac {\overline {\rho u_i'' u_i''}}{2} \right ) = \frac {\partial }{\partial x_j} \left [ \tilde {u}_i \left ( \bar {\tau }_{\textit{ij}} - \overline {\rho u_i'' u_{\kern-1pt j}''} \right ) \right ]\\\nonumber & + \frac {\partial }{\partial x_j} \left ( \overline {\tau _{\textit{ij}} u_i''} - \frac {\overline {\rho u_i'' u_i'' u_{\kern-1pt j}''}}{2} \right ) - \frac {\partial }{\partial x_j} \left (-\bar {\kappa } \frac {\partial \tilde {T}}{\partial x_j} + c_p \overline {\rho T'' u_{\kern-1pt j}''} \right )\!, \\[-4pt] \nonumber \end{align}
where
$u_i$
is the velocity component along the
$i$
th direction (with
$i = 1,2,3$
),
$\rho$
is the density,
$p$
denotes pressure and
$T$
the temperature. For a given instantaneous flow variable
$q$
, the following decompositions are applied:
$q = \bar {q} + q' = \tilde {q} + q''$
. The symbol
$\bar {\boldsymbol{\cdot }}$
indicates time averaging, whereas
$\tilde {\boldsymbol{\cdot }}$
denotes Favre averaging, defined as
\begin{align} \tilde {q} = \frac {\overline {\rho q}}{\bar {\rho }}. \end{align}
The quantity
$\bar {\tau }_{\textit{ij}}$
is the viscous stress tensor, given by
\begin{align} \bar {\tau }_{\textit{ij}} = \bar {\mu } \left ( \frac {\partial \tilde {u}_i}{\partial x_j} + \frac {\partial \tilde {u}_j}{\partial x_i} \right ) - \frac {2}{3} \bar {\mu } \frac {\partial \tilde {u}_k}{\partial x_k} \delta _{\textit{ij}}, \end{align}
where
$\delta _{\textit{ij}}$
is the Kronecker delta. The quantity
$-\overline {\rho u_i'' u_{\kern-1pt j}''}$
denotes the Favre-averaged Reynolds stress tensor, which is modelled by means of the Boussinesq eddy-viscosity approximation. Here,
$\tilde {E}$
refers to the Favre-averaged total energy, defined as
\begin{align} \tilde {E} = c_v \tilde {T} + \frac {\tilde {u}_i \tilde {u}_i}{2}, \end{align}
where
$c_v$
represents the specific heat at constant volume, while
$\tilde {H}$
is the Favre-averaged total enthalpy, given by
$\tilde {H} = \tilde {E} + \bar {p}/\bar {\rho }$
. The quantity
$(1/2) \overline {\rho u_i'' u_i''} = \bar {\rho } k$
denotes the kinetic energy of the turbulent fluctuations per unit volume. The transport coefficients representing the dynamic viscosity and thermal conductivity are respectively denoted by
$\bar {\mu }$
and
$\bar {\kappa }$
. The evolution of viscosity with temperature follows Sutherland’s law and
$\bar {\kappa }$
is obtained from the assumption of a constant Prandtl number. Finally, the averaged perfect gas equation of state reads
\begin{align} \bar {p} = \bar {\rho } R_g \tilde {T}, \end{align}
where
$R_g = 287$
J kg−1 K−1 is the specific gas constant.
The two-equation shear stress transport (SST) turbulence model introduced by Menter (Reference Menter1994) is considered, together with the compressible mixing-layer correction proposed by Wilcox (Reference Wilcox2006, (5.83)). Modified values of the empirical constants of the turbulence model have been employed and are listed in table 1. These values were obtained by manual calibration following a linear interpolation between the original SST constants provided by Menter (Reference Menter1994) and the optimised constants recently obtained by Ozawa & Nonomura (Reference Ozawa and Nonomura2024) via data assimilation for a perfectly expanded jet. The selected values were found to yield a small mean squared error with respect to PIV mean-flow measurements with the employed solver. Appendix A provides a comparison between the mean-flow results obtained with the three different sets of empirical constants, which justifies the choice of parameters for this work.
Table 1. Employed values for the parameters of the Menter SST turbulence model. The same nomenclature as in Menter (Reference Menter1994) is followed.
The reference quantities employed for non-dimensionalisation are the nozzle-exit diameter
$D$
and the free stream (ambient) flow conditions. The flow velocity components are non-dimensionalised with the free stream speed of sound
$c_\infty$
; the density with the respective free stream value
$\rho _\infty$
; the pressure is made dimensionless with
$\rho _\infty c_\infty ^2$
and the temperature with
$(\gamma - 1) T_\infty$
. The transport properties
$\mu$
and
$\kappa$
are non-dimensionalised with their respective free stream values.
3.1.2. Numerical methodology for the RANS calculations
Figure 2. Computational domain employed for the RANS calculation of the twin-jet flow field, including boundary conditions:
$(a)$
cross-stream plane at a fixed streamwise location;
$(b)$
$xy$
symmetry plane plane at
$z = 0$
. The dashed lines mark the nozzle exit geometry.
The RANS calculations are performed using the CFD solver TAU, developed by the German Aerospace Centre (DLR) (Schwamborn, Gerhold & Heinrich Reference Schwamborn, Gerhold and Heinrich2006). The numerical solution employs a second-order finite-volume method with an upwind discretisation scheme based on Roe’s approximate Riemann solver. A steady RANS solution is targeted, which is approached using a backward Euler implicit scheme for time integration.
Figure 2 shows the computational domain employed for the three-dimensional RANS calculations. Only one quarter of the twin-jet system is simulated, taking advantage of the two symmetry planes inherent in the geometry. The size of the domain is
$L_x = 35D$
in the streamwise direction. In the
$z$
direction, a size of
$L_z = 8D$
is used at the domain inflow, which extends linearly until
$L_z = 14D$
at the domain outflow. Similarly, in the
$y$
direction,
$L_y = 9D$
is employed at the inflow, while
$L_y = 15D$
is used at the outflow. The inflow is located at the nozzle exit plane (
$x/D = 0$
), where the nozzle exit flow field is imposed as a Dirichlet boundary condition. The flow at the nozzle exit is computed by means of an axisymmetric compressible RANS calculation of the entire nozzle flow field, including the boundary layer developing on the nozzle inner walls. This solution is shown in figure 3, which shows contours of the streamwise velocity field inside the nozzle and in the vicinity of the nozzle exit, and in figure 4, which illustrates the mean velocity and temperature profiles at the nozzle exit (
$x/D = 0$
) and slightly downstream of it (
$x/D = 0.05$
and
$0.1$
). The profile corresponding to
$x/D = 0$
is the one imposed at the inflow of the twin-jet calculation.
Figure 3. Contours of mean streamwise velocity for the axisymmetric (single) jet calculation.
Figure 4. Mean flow profiles extracted at the nozzle exit (
$x/D = 0$
) and slightly downstream of the nozzle exit (
$x/D = 0.05$
and
$x/D = 0.1$
) from the axisymmetric (single jet) calculation:
$(a)$
streamwise velocity;
$(b)$
static temperature. The black dotted horizontal lines indicate the radial position of the nozzle lip walls.
The single-jet axisymmetric RANS calculation uses a two-dimensional computational domain with a hybrid mesh mainly consisting of triangular cells, together with a structured region of quadrilateral cells located at the nozzle boundary layer. The cell size is calibrated such that
$y^+ \approx 1$
to resolve the boundary layer developing on the nozzle walls, which feature an adiabatic thermal boundary condition. The resulting number of grid points across the boundary layer thickness is 50. A small co-flow of
$M_{\textrm {co}} = 0.01$
is included in the far-field to replicate the experimental conditions.
Table 2 reports the total and the ambient flow conditions employed for both the single-jet and the twin-jet RANS calculations. In addition to the jet-exit Mach number
$M_{\kern-1pt j} = u_{\kern-1pt j}/c_{\kern-1pt j}$
, the acoustic jet Mach number is also reported, defined as
$M = u_{\kern-1pt j}/c_\infty$
. Here,
$u_{\kern-1pt j}$
is the jet-exit velocity and
$c_{\kern-1pt j}$
is the speed of sound corresponding to the jet-exit temperature (
$T_{\kern-1pt j}$
), which is determined assuming an isentropic expansion. Two different Reynolds numbers are considered, the reference, acoustic Reynolds number used for non-dimensionalisation,
$ \textit{Re} = \rho _\infty c_\infty D / \mu _\infty$
, and the jet-exit Reynolds number,
$ \textit{Re}_j = \rho_{\kern-1pt j} u_{\kern-1pt j} D / \mu_{\kern-1pt j}$
, where
$\rho_{\kern-1pt j}$
and
$\mu_{\kern-1pt j}$
respectively denote the density and the dynamic viscosity of the jet at the nozzle exit, which are also determined based on the isentropic flow relations. The Prandtl number is defined as
$ \textit{Pr} = c_p \mu _\infty / \kappa _\infty$
, where
$c_p$
is the specific heat at constant pressure.
For the three-dimensional calculations, a fully unstructured tetrahedral mesh is used. It comprises successive refinement regions downstream of the nozzle exit following the growth of the mixing layers. A total of 154 million tetrahedral cells are employed for the case
$s/D = 1.76$
and 174 million cells for
$s/D = 3$
. The impact of the grid resolution of the RANS calculations on the predicted wavepackets is discussed in Appendix B. At the symmetry planes
$z = 0$
and
$y = 0$
of the domain, symmetry boundary conditions are imposed. At the domain outflow, a subsonic outflow boundary condition is used, which imposes the ambient pressure
$p_\infty$
and extrapolates the other flow variables from the interior of the domain. At the far-field boundaries, a far-field inflow/outflow boundary condition is imposed that evaluates boundary fluxes via a characteristics method.
3.2. Formulation of the plane-marching parabolised stability equations
The PSE constitute one approach for studying the growth of linear and nonlinear disturbances in shear flows, and are well adapted when there exists a slow divergence of the mean-flow properties in the streamwise direction. They have been shown to deliver results comparable to direct numerical simulations for convectively unstable laminar and transitional flows (Bertolotti, Herbert & Spalart Reference Bertolotti, Herbert and Spalart1992; Chang et al. Reference Chang, Malik, Erlenbacher and Hussaini1993). PSE has been extensively used as a model for the spatial evolution of coherent structures in single round turbulent jets (wavepackets), yielding good agreement with coherent structures educed either from experiments or large eddy simulations (Yen & Messersmith Reference Yen and Messersmith1999; Piot et al. Reference Piot, Casalis, Muller and Bailly2006; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Cavalieri et al. Reference Cavalieri, Rodríguez, Jordan, Colonius and Gervais2013; Rodríguez et al. Reference Rodríguez, Sinha, Brès and Colonius2013; Sinha et al. Reference Sinha, Rodríguez, Brès and Colonius2014; Sasaki et al. Reference Sasaki, Cavalieri, Jordan, Schmidt, Colonius and Brès2017). For perfectly expanded supersonic jets, screech does not occur and the sound generated at this condition is largely dominated by downstream-propagating Kelvin–Helmholtz waves. Owing to its lower computational cost and complexity compared with alternative techniques, PSE remains an attractive simplified modelling framework for this problem.
Following recent investigations (Rodríguez et al. Reference Rodríguez, Jotkar and Gennaro2018; Rodríguez Reference Rodríguez2021), wavepackets associated with mixing noise in perfectly expanded supersonic twin jets are modelled here using linear plane-marching parabolised stability equations (PM-PSE). This model is a direct extension of the classical PSE approach that accounts for mean flows featuring a strong inhomogeneity in the cross-stream planes, i.e.
$y$
and
$z$
(Broadhurst & Sherwin Reference Broadhurst and Sherwin2008).
Table 2. Flow conditions employed for the RANS calculations and resulting dimensionless parameters.
The plane-marching parabolised stability equations employed in this study are derived from the compressible Navier–Stokes equations. To formulate the PM-PSE linear stability problem, first, the instantaneous turbulent flow field is separated into a stationary mean-flow component
$\bar {\boldsymbol{q}}$
and an unsteady fluctuation component
$\boldsymbol{q}'$
:
\begin{align} \boldsymbol{q} (x,y,z,t) = \bar {\boldsymbol{q}} (x,y,z) + \boldsymbol{q}'(x,y,z,t), \end{align}
where
$\boldsymbol{q} = [u,v,w,p,T]^{\mathrm{T}}$
is the vector of primitive flow variables, with
$u$
,
$v$
and
$w$
representing the velocity component along
$x$
,
$y$
and
$z$
, respectively. Next, the fluctuating component is expressed as a sum of discrete modes in frequency:
\begin{align} \boldsymbol{q}'(x,y,z,t) = \sum _\omega \hat {\boldsymbol{q}}_\omega (x,y,z) \mathrm{e}^{-\mathrm{i}\omega t} + \textrm{c.c.}, \end{align}
where
$\omega = 2\pi f$
denotes the angular frequency and
$\textrm{c.c.}$
stands for the complex conjugate. The Fourier-transformed fluctuations in time
$\hat {\boldsymbol{q}}_\omega$
are then decomposed into a shape function (
$\tilde {\boldsymbol{q}}_\omega$
), which is assumed to undergo a slow variation along the streamwise direction
$x$
, and a wave function (
$A_\omega$
), which is allowed to vary rapidly with
$x$
:
\begin{align} \hat {\boldsymbol{q}}_\omega (x,y,z) = A_\omega (x) \ \tilde {\boldsymbol{q}}_\omega (x,y,z) = A_\omega (x_0) \exp \left (\mathrm{i} \int ^x _{x_0} \alpha _\omega (\xi ) \,\mathrm{d} \xi \right )\ \tilde {\boldsymbol{q}}_\omega (x,y,z), \end{align}
where the complex quantity
$\alpha _\omega = \alpha _r + \mathrm{i} \alpha _i$
is the streamwise wavenumber for frequency
$\omega$
, for which a slow variation with
$x$
is also required. The streamwise coordinate
$x_0$
is the location where the PSE integration is initialised. In the context of this work, this is typically a cross-section close to the nozzle exit. The quantity
$A_\omega (x_0)$
refers to the initial disturbance amplitude, which is arbitrary in a linear problem. Introducing (3.9) and (3.10) into the linearised Navier–Stokes equations and performing an order of magnitude analysis to neglect terms of order
$1/Re^2$
, the following linear system is obtained:
\begin{align} \boldsymbol{L} \frac {\partial \tilde {\boldsymbol{q}}_\omega }{\partial x} = \boldsymbol{R} \tilde {\boldsymbol{q}}_{\omega }, \end{align}
where the linear operators
$\boldsymbol{L}$
and
$\boldsymbol{R}$
depend on the mean flow quantities and their first and second spatial derivatives, the complex spatial wavenumber
$\alpha _\omega$
and real frequency
$\omega$
, and the non-dimensional parameters
$ \textit{Re}$
,
$ \textit{Pr}$
and
$\gamma$
. Additional details on the derivation of the operators
$\boldsymbol{L}$
and
$\boldsymbol{R}$
can be found from Rodríguez et al. (Reference Rodríguez, Jotkar and Gennaro2018).
The ansatz (3.10) allows the spatial disturbance growth to be accounted for by either the shape function or the wave function, making the solution non-unique. To solve this ambiguity, a normalisation is imposed to eliminate the exponential dependence from
$\tilde {\boldsymbol{q}}_\omega$
, which takes the form (Herbert Reference Herbert1997):
\begin{align} \iint _\varOmega \tilde {\boldsymbol{q}}^* _\omega \frac {\partial \tilde {\boldsymbol{q}}_\omega }{\partial x} \,\mathrm{d}y\mathrm{d}z = 0, \end{align}
where the superscript
$^*$
denotes the complex conjugate and
$\varOmega$
is the cross-stream spatial domain in which the shape functions are defined for a given streamwise location.
3.2.1. Initial condition: locally parallel linear stability analysis
The PM-PSE system takes the form of a downstream-marching problem in the streamwise coordinate
$x$
, which requires initial conditions at a given location
$x_0$
. To obtain a consistent set of initial conditions for the PM-PSE problem, a locally parallel stability problem is derived from the PSE formulation following that of Rodríguez et al. (Reference Rodríguez, Sinha, Brès and Colonius2013). This is achieved by assuming that
$\mathrm{d} \alpha _\omega / \mathrm{d} x = 0$
,
$\partial \hat {\boldsymbol{q}}_\omega / \partial x = \mathrm{i} \alpha _\omega \tilde {\boldsymbol{q}}_\omega$
, and that, for a given streamwise station
$x$
, the mean flow is locally parallel:
$\partial \bar {\boldsymbol{q}} / \partial x = 0$
. These assumptions yield a generalised eigenvalue problem of the form:
\begin{align} \mathrm{i} \alpha _\omega \boldsymbol{L} \tilde {\boldsymbol{q}}_\omega = \boldsymbol{R} \tilde {\boldsymbol{q}}_\omega . \end{align}
For a given cross-stream plane at an initial location
$x_0$
and a given frequency
$\omega$
, the solution of the eigenvalue problem (3.13) delivers a set of complex eigenvalues and their corresponding two-dimensional eigenfunctions (defined on the cross-stream plane). For unbounded mean flows like the present ones, classical linear stability theory (Mack Reference Mack1984; Balakumar & Malik Reference Balakumar and Malik1992) shows the existence of a fixed number of continuous branches, related to the uniform flow surrounding the jets, in addition to an indefinite number of discrete eigenmodes associated with localised mean-flow shear. For a single round jet, the different eigenmode families are described by Rodríguez et al. (Reference Rodríguez, Sinha, Brès and Colonius2013, Reference Rodríguez, Cavalieri, Colonius and Jordan2015). Identification of the discrete modes of interest (Kelvin–Helmholtz instabilities in this case) from the local stability spectrum is achieved by visual inspection of the corresponding two-dimensional amplitude functions.
3.2.2. Numerical methodology for the plane-marching PSE calculations
The numerical solution of both the marching problem and the local stability problem require the spatial discretisation of the two-dimensional linear operators
$\boldsymbol{L}$
and
$\boldsymbol{R}$
. Following the studies by Gennaro et al. (Reference Gennaro, Rodríguez, Medeiros and Theofilis2013) and Rodríguez & Gennaro (Reference Rodríguez and Gennaro2017), finite-difference stencils of variable number of points are used in both
$y$
- and
$z$
-directions. A seven-point stencil is used in the inner points. The stencil size is gradually reduced towards the boundary points to preserve the banded structure of the differentiation matrices. At the boundary points on each coordinate direction, four-point forward and backward stencils are used. This approach, employing high-order finite differences with banded differentiation matrices, leverages on the sparse storage and algebra employed, offering a good trade-off between convergence of results and computational cost. According to the mean-flow symmetries, only a quarter of the domain needs to be discretised. Consequently, a rectangular domain
$\varOmega = [0,y_{\infty }]\times [0,z_{\infty }]$
is used for the cross-stream planes. Two independent coordinate transformations (mappings) are used to concentrate grid points in the jet mixing regions. The following transformation is used in both
$y$
and
$z$
directions:
\begin{align} \eta = \eta _c + (\eta _\infty - \eta _c) \frac {\sinh \left (a(\xi - b) \right )}{\sinh \left (a(1-b)\right )}, \end{align}
where
$\eta$
is the mapped coordinate in the physical domain (
$y$
or
$z$
),
$\xi$
is the associated coordinate in the computational domain (
$\xi \in [0,1]$
);
$\eta _c$
is the location of the jet axis, where the discretisation is refined;
$\eta _\infty$
is the maximum coordinate and
$a$
is a real number that controls the intensity of the clustering of points. The remaining parameter
$b$
must be determined for each combination of
$\eta _c$
,
$\eta _\infty$
and
$a$
so that
$\eta = 0$
for
$\xi = 0$
. In the present calculations, a square domain
$\varOmega = [0,10] \times [0,10]$
is used, discretised with a spatial resolution of
$N_y \times N_z = 201 \times 201$
points. For each case, the values
$a = 5$
,
$y_c/D = s/(2D)$
and
$z_c/D = 0.5$
are used. Mesh convergence tests, not reproduced here, ensure the numerical consistency of the results using these parameters.
The Cartesian coordinate system allows for the use of standard finite differences for the independent differentiation on
$y$
and
$z$
, resulting in the differentiation matrices
$\mathcal{D}_y$
and
$\mathcal{D}_z$
for first-order derivatives, and
$\mathcal{D}_{\textit{yy}}$
and
$\mathcal{D}_{zz}$
for second-order derivatives. The same stencil is used for first- and second-order differentiation matrices, allowing for the control of the matrix structure required for efficiency of the sparse implementation. Explicit finite-difference formulae are employed for the second-order derivatives in
$\mathcal{D}_{\textit{yy}}$
and
$\mathcal{D}_{zz}$
. The cross-differentiation matrix is then obtained as
$\mathcal{D}_{yz} = \mathcal{D}_y \times \mathcal{D}_z$
.
For a given cross-stream plane, four possible sets of boundary conditions can be imposed, leading to four solution families (Morris Reference Morris1990; Rodríguez et al. Reference Rodríguez, Jotkar and Gennaro2018, Reference Rodríguez, Stavropoulos, Nogueira, Edgington-Mitchell and Jordan2023): (i) symmetry conditions with respect to both
$z = 0$
and
$y = 0$
(denoted by SS); (ii) symmetry with respect to
$z = 0$
and antisymmetry with respect to
$y = 0$
(denoted by SA); (iii) antisymmetry with respect to
$z = 0$
and symmetry with respect to
$y = 0$
(denoted by AS); (iv) antisymmetry conditions with respect to both
$z = 0$
and
$y = 0$
(denoted by AA). The symmetry imposition is accomplished through appropriate sets of symmetry and antisymmetry boundary conditions on the perturbation quantities, which are implemented through the corresponding Neumann and Dirichlet conditions. These are summarised as follows (refer to figure 2):
\begin{eqnarray} \begin{array}{ll} \mbox{Symmetry on } y=0, & \dfrac {\partial \hat {u}}{\partial y} = \dfrac {\partial \hat {w}}{\partial y} = \dfrac {\partial \hat {p}}{\partial y} = \dfrac {\partial \hat {T}}{\partial y} = \hat {v} = 0; \\ \mbox{Antisymmetry on } y=0, & \hat {u} = \hat {w} = \hat {p} = \hat {T} = \dfrac {\partial \hat {v}}{\partial y} = 0; \\ \mbox{Symmetry on } z=0, & \dfrac {\partial \hat {u}}{\partial z} = \dfrac {\partial \hat {v}}{\partial z} = \dfrac {\partial \hat {p}}{\partial z} = \dfrac {\partial \hat {T}}{\partial z} = \hat {w} = 0; \\ \mbox{Antisymmetry on } z=0, & \hat {u} = \hat {v} = \hat {p} = \hat {T} = \dfrac {\partial \hat {w}}{\partial z} = 0. \\ \end{array} \end{eqnarray}
The local stability problem is solved by means of a parallelised sparse implementation of the shift-and-invert Arnoldi algorithm (Arnoldi Reference Arnoldi1951), which employs the MUMPS library (Amestoy et al. Reference Amestoy, Duff, L’Excellent and Koster2001) for solution of the linear systems of equations. The implicit Euler scheme is used to march the PM-PSE solution downstream. The resulting linear system of equations at each streamwise station is also solved via a parallel implementation using MUMPS. To adjust the value of
$\alpha$
so that the normalisation condition (3.12) is satisfied, the solution of the linear system is iterated together with the following relation:
\begin{align} \alpha ^{(k+1)}_{\omega , j+1} = \alpha ^{(k)}_{\omega , j+1} - \frac {\mathrm{i}}{\Delta x} \frac {\iint _{\varOmega } \tilde {\boldsymbol {q}}_{\omega , j+1}^* (\tilde {\boldsymbol {q}}_{\omega , j+1} - \tilde {\boldsymbol {q}}_{\omega , j}) \,\mathrm{d}y \,\mathrm{d}z}{\iint _\varOmega \tilde {\boldsymbol {q}}_{\omega , j+1}^* \tilde {\boldsymbol {q}}_{\omega , j+1} \,\mathrm{d}y\, \mathrm{d}z}, \end{align}
where
$k$
is the iteration index,
$j+1$
refers to the solution at the next downstream streamwise position,
$j$
refers to the current streamwise station and
$\Delta x = x_{j+1} - x_j$
is the streamwise step size. The iteration continues until
$\alpha _\omega$
is converged up to a relative error
$|\alpha ^{(k+1)}_\omega -\alpha ^{(k)}_\omega |/|\alpha ^{(k+1)}_\omega | \lt 10^{-4}$
.
For all the calculations presented in this work, the locally parallel stability problem is solved at
$x_0 = 0.75 D$
. This is an arbitrary location relatively close to the nozzle lip where robust results are obtained. Computing the initial marching conditions in the close proximity of the nozzle lip would be influenced by the sharp gradients of the mean flow, as shown in figure 4, requiring a more refined cross-plane mesh. The streamwise step size
$\Delta x \approx 0.2 D$
is used for all cases. This value was determined empirically as the smallest one for which the PSE marching does not present instability issues for the most critical
$St=0.1$
cases. Stabilisation procedures, such as the one introduced by Andersson, Henningson & Hanifi (Reference Andersson, Henningson and Hanifi1998), are not employed.
4. Experimental set-up
Figure 5. Experimental set-up and elements of the high-speed schlieren measurement system.
High-speed schlieren visualisations and low-speed PIV measurements of the twin-jet flow field have been performed at the T200 wind tunnel facility of the PROMETEE platform of Institut Pprime (CNRS – Université de Poitiers – ISAE-ENSMA), France. The set-up is shown in figure 5. The flow is generated using a
$200$
bar compressed air network that can reach operational conditions up to an isentropic Mach number
$M_{\kern-1pt j} = 2$
for the twin-jet configuration. A heating system based on a series of tanks with heated nickel balls is used to increase and maintain the total temperature of the air reaching the nozzles. The left image of figure 5 shows the twin-jet system situated in an anechoic chamber.
Schlieren visualisations are obtained using a classical Z-type set-up, for which some components are also illustrated in figure 5. A continuous light source is generated by a 60 W LED, which passes through an aperture that prevents direct light from the source to enter the test section. Two parabolic mirrors, each 30 cm in diameter and with a 3 m focal length, produce a collimated light beam that travels through the test section in the
$z$
direction, according to the reference frame in figure 1. To form a Z-shaped optical path that allows for a more compact experimental set-up, two additional flat mirrors, each 12 cm in diameter, are incorporated. A vertical knife edge is positioned at the focal length of the second parabolic mirror, making the resulting light intensity field recovered in the schlieren images proportional to the streamwise density gradients in the flow.
A Phantom v2640 high-speed camera was used to record the images. For the current datasets,
$N_s = 30\,000$
instantaneous snapshots were recorded at a sampling frequency of
$f_s = 68$
kHz (
$\Delta t = 1.47 \times 10^{-5}$
s) and a spatial resolution of
$352 \times 512$
pixels. The Strouhal number corresponding to the recording Nyquist frequency is
$St_{\textit{Ny}} = f_s D / (2 u_{\kern-1pt j}) = 1.93$
(
$2 \Delta t u_{\kern-1pt j} / D = 0.52$
), which is high enough to resolve the coherent structures associated with mixing noise, which develop in the range
$St \approx [0.1, 1]$
for supersonic jets (Sinha et al. Reference Sinha, Rodríguez, Brès and Colonius2014). The total number of convective time units spanned by the recording is
$\tau c_{\textit{KH}} / D \approx 5440$
, where
$\tau = \Delta t N_s = 0.44$
s denotes the total recording time and
$c_{\textit{KH}} = 0.7u_{\kern-1pt j}$
is the estimated convective velocity of Kelvin–Helmholtz structures. Figure 6
$(a{,}b$
) illustrates two instantaneous schlieren snapshots respectively recorded for
$s/D = 1.76$
and
$s/D = 3$
. The mean schlieren fields resulting from averaging all the instantaneous snapshots are also shown in figure 6
$(c{,}d$
). Due to the high sensitivity of the perfectly expanded condition to imperfections in the nozzle geometry, weak residual shock and expansion waves remain in the flow. A small discrepancy is present between the two jets, which originates from small machining inaccuracies between the two nozzles. The weak asymmetry introduced by this difference between both jets has not been found to have a meaningful impact on the coherent fluctuations supported by the system and their eduction based on their symmetric and antisymmetric behaviour, and hence does not affect the validation process of the wavepacket models.
PIV measurements are performed to obtain a mean velocity field to validate the RANS calculations. The PIV measurement system consists of two side-by-side LaVision Imager LX 16 M cameras which record 500 instantaneous, high-resolution images (dual frames) at a sampling frequency of 4.4 Hz and at a resolution of 16 megapixels each, employing a time delay between frames of 1.5
$\unicode {x03BC}$
s. The flow is seeded with pressurised glycerine vapour in the settling chamber, prior to nozzle expansion, while the ambient region is seeded with a fog generator employing a glycerine–water mixture. A double-pulse Nd-YAG laser operating at 532 nm is used, which illuminates the seeded particles in the
$z = 0$
plane, enabling measurement of the mean-flow velocity in the symmetry plane containing the two jets. The PIV measurement window covers the region
$(x/D,y/D) \in [0,10]\times [-5,5]$
, extending downstream approximately up to the end of the potential core. After processing, the resulting instantaneous velocity fields have a resolution of
$695 \times 622$
pixels.
Figure 6. Schlieren snapshots of the twin-jet flow field at
$M_{\kern-1pt j} = 1.54$
:
$(a{,}c)$
$s/D = 1.76$
;
$(b{,}d)$
$s/D = 3$
;
$(a{,}b)$
instantaneous snapshots;
$(c{,}d)$
mean of all snapshots. The dark regions near the downstream boundary of the images correspond to the edges of the right parabolic mirror.
4.1. Eduction of coherent structures from schlieren visualisations
To validate the wavepackets modelled by plane-marching PSE, coherent structures are extracted by means of spectral proper orthogonal decomposition (SPOD) (Lumley Reference Lumley1970; Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018) applied to the high-speed schlieren visualisations. The coherent structures educed by SPOD in subsonic and perfectly expanded supersonic turbulent jets are known to represent wavepackets related to the Kelvin–Helmholtz instability of the shear layers and non-modal dynamics including Orr and lift-up mechanisms (Suzuki & Colonius Reference Suzuki and Colonius2006; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Tissot et al. Reference Tissot, Lajús, Cavalieri and Jordan2017; Jordan et al. Reference Jordan, Zhang, Lehnasch and Cavalieri2017; Sasaki et al. Reference Sasaki, Cavalieri, Jordan, Schmidt, Colonius and Brès2017; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Lesshafft et al. Reference Lesshafft, Semeraro, Jaunet, Cavalieri and Jordan2019; Nogueira et al. Reference Nogueira, Cavalieri, Jordan and Jaunet2019; Pickering et al. Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020).
The application of SPOD to schlieren visualisations of twin jets has proved successful in obtaining empirical descriptions of the screech phenomenon (Edgington-Mitchell et al. Reference Edgington-Mitchell, Li, Liu, He, Wong, MacKenzie and Nogueira2022), which involves a highly organised structure at localised tonal frequencies. The perfectly expanded twin-jet flow is convectively unstable for a broad band of frequencies, which leads to the formation of flow structures with a wide range of spatial scales. Schlieren visualisation emphasises the small-scale vortical structures, making the extraction of coherent structures difficult: a large number of snapshots is necessary to reach converged SPOD modes. For this reason, instead of directly using the schlieren images to educe the coherent structures, a different quantity derived from the schlieren field is used to build the cross-spectral density (CSD) matrix.
Doak’s momentum potential theory (Doak Reference Doak1989; Jordan, Daviller & Comte Reference Jordan, Daviller and Comte2013) can be employed to extract the irrotational momentum potential component associated with the line-of-sight integrated momentum density field (Prasad & Gaitonde Reference Prasad and Gaitonde2022), denoted by
$\varTheta$
, from the schlieren field, denoted by
$\sigma$
. As illustrated by Padilla-Montero et al. (Reference Padilla-Montero, Rodríguez, Jaunet and Jordan2024), for a given dataset, using
$\varTheta$
instead of
$\sigma$
to build the CSD matrix in SPOD allows for a more effective extraction of coherent structures, facilitating comparison with wavepacket models based on linear stability theory. The irrotational field acts as a filter that removes much of the small-scale turbulence fluctuations from the schlieren images, reducing the dimensional complexity of the dataset.
The same methodology described by Padilla-Montero et al. (Reference Padilla-Montero, Rodríguez, Jaunet and Jordan2024) is here employed for educing coherent structures. Doak’s momentum potential theory yields a relationship between the density fluctuation field (
$\rho '$
) and the momentum potential fluctuation field (
$\psi '$
) in the form of a linear Poisson equation (Doak Reference Doak1989):
\begin{align} \frac {\partial ^2 \psi '}{\partial x_i^2} = \frac {\partial \rho '}{\partial t}. \end{align}
Following Prasad & Gaitonde (Reference Prasad and Gaitonde2022), (4.1) can be differentiated with respect to
$x$
, yielding
\begin{equation} \frac {\partial ^2}{\partial x_i^2} \left ( \frac {\partial \psi '}{\partial x} \right ) = \frac {\partial }{\partial t} \left ( \frac {\partial \rho '}{\partial x} \right )\!, \end{equation}
and can be integrated along
$z$
(the line of sight of the schlieren measurements) to obtain
\begin{equation} \frac {\partial ^2 \varTheta '}{\partial x_i^2} = \frac {\partial \sigma '}{\partial t}, \end{equation}
where
$\sigma ' = \int (\partial \rho ' / \partial x)\,\mathrm{d}z$
denotes the schlieren fluctuation field and
$\varTheta ' = \int (\partial \psi ' / \partial x)\,\mathrm{d}z$
is the line-of-sight integrated streamwise derivative of the momentum potential fluctuation. Equation (4.3) can be solved using ad hoc boundary conditions on the truncated domain given by the schlieren images, allowing calculation of the
$\varTheta '$
field associated with the experimental measurements. Here, solution of the Poisson equation is incorporated into the SPOD algorithm and
$\varTheta '$
is computed in the spectral domain. Special care is taken to filter out spurious harmonic components that arise from solution of the Poisson equation in a truncated domain with unknown boundary conditions. The reader is referred to Padilla-Montero et al. (Reference Padilla-Montero, Rodríguez, Jaunet and Jordan2024) for details of the solution approach.
When (4.3) is applied to the experimental schlieren snapshots,
$\sigma$
represents the pixel intensity obtained directly from the schlieren experimental system without any calibration, i.e. it is a field proportional to
$\partial \rho / \partial x$
. The particular proportionality constant depends on the elements of the experimental set-up, such as the camera sensor, the knife edge and the mirrors. This does not have an impact on the validation methodology employed in this work as it is based on the structure of the schlieren fields and not on their magnitude.
To facilitate comparison of PM-PSE modes against experimentally educed structures, the symmetry of the system along the
$y = 0$
plane is also exploited in the treatment of the experimental datasets. Before solving the SPOD eigenvalue problem, two new sets of
$\varTheta '$
realisations are computed by splitting each
$\hat {\varTheta }$
snapshot by the line at
$y = 0$
, namely, a symmetric set
$\hat {\varTheta }_s = (\hat {\varTheta }_u + \hat {\varTheta }_l)/2$
and an antisymmetric one
$\hat {\varTheta }_a = (\hat {\varTheta }_u - \hat {\varTheta }_l)/2$
, with
$\hat {\varTheta }_u$
and
$\hat {\varTheta }_l$
respectively denoting the upper (
$y \gt 0$
) and lower (
$y \lt 0$
) halves of the original realisations. The SPOD eigenvalue problem is then solved for the symmetric and antisymmetric datasets separately. This decomposition is analogous to the
$D_2$
decomposition (Sirovich & Park Reference Sirovich and Park1990) employed for the study of rectangular twin jets (Yeung et al. Reference Yeung, Schmidt and Brès2022), which allows symmetric and antisymmetric fluctuations to be extracted without loss of generality. This permits a systematic comparison with the wavepackets computed using PM-PSE featuring symmetric/antisymmetric boundary conditions.
Each schlieren dataset is divided into 57 blocks of 1024 snapshots, with an overlap of 50 %, yielding frequency bins of
$\Delta f = 66$
Hz (
$\Delta St = 0.0038$
). A temporal Hamming window is used to reduce spectral leakage. The SPOD decomposition of the momentum potential fluctuation field yields SPOD modes in terms of
$\varTheta '$
. To compare the educed coherent structures with the PM-PSE wavepackets, SPOD modes expressed in terms of schlieren fluctuations are also of interest. These SPOD fields can be computed as an extended SPOD problem (see for example Freund & Colonius (Reference Freund and Colonius2002), Borée (Reference Borée2003), Sinha et al. (Reference Sinha, Rodríguez, Brès and Colonius2014), Souza et al. (Reference Souza, Rodríguez, Himeno and Medeiros2019), Himeno et al. (Reference Himeno, Souza, Amaral, Rodríguez and Medeiros2021), Karban et al. (Reference Karban, Bugeat, Towne, Lesshafft, Agarwal and Jordan2023)). In the results shown later in this work, two different flow-field variables are shown for each SPOD mode, namely, the
$\hat {\varTheta }$
field and the
$\hat {\sigma }_\varTheta$
field, which denotes the schlieren fluctuation field corresponding to the coherent structure educed based on the cross-spectral density of
$\varTheta '$
.
Figure 7. Comparison of mean streamwise velocity profiles along
$y$
(and at
$z = 0$
) for four different streamwise locations, between the RANS solution and the PIV mean flow for
$s/D = 1.76$
:
$(a)$
$x/D = 2$
;
$(b)$
$x/D = 4$
;
$(c)$
$x/D = 6$
;
$(d)$
$x/D = 8$
.
Figure 8. Comparison of mean streamwise velocity profiles along
$y$
(and at
$z = 0$
) for four different streamwise locations, between the RANS solution and the PIV mean flow for
$s/D = 3$
:
$(a)$
$x/D = 2$
;
$(b)$
$x/D = 4$
;
$(c)$
$x/D = 6$
;
$(d)$
$x/D = 8$
.
5. Results
5.1. Mean flow
First, the RANS mean flow employed for the PSE calculations is discussed and compared with the PIV measurements. Figures 7 and 8 show the mean streamwise velocity profiles obtained from the RANS computations at the symmetry plane
$z = 0$
. For each jet separation, the RANS and PIV fields are in good agreement. The RANS model accurately captures the mixing-layer thickness, the centreline velocity and the velocity increase in the region between the two jets. The profiles highlight the interaction between the jets that takes place at the merging region between the inner shear layers. In this region, the flow velocity increases above that of the equivalent isolated jet, making the velocity gradient (shear magnitude) in the inner mixing layers smaller with respect to that of the external one. This effect reflects the loss of axisymmetry, which is significantly stronger for
$s/D = 1.76$
. Small discrepancies are progressively found for
$x/D \gt 6$
, which manifest as a small overprediction of the external mixing layer thickness and of the velocity magnitude in the region between the jets.
Figure 9 shows the evolution of the inner and outer shear-layer boundaries of the twin jet for each spacing in both the RANS and PIV fields. This is defined by the
$\bar {u}/u_{\kern-1pt j} = 0.05$
velocity contour. As discussed already, the comparison against the PIV measurements reveals a small overprediction of the RANS mixing-layer thickness which increases further downstream. The axisymmetric (single jet) solution is also shown for comparison. In the twin-jet system, the spread rate of the outer mixing layer is reduced with respect to that of the isolated jet, while the spreading of the inner mixing layer is increased owing to entrainment of quiescent fluid in the region between the jets. These trends are consistent with the findings of previous work (Moustafa Reference Moustafa1994; Alkislar et al. Reference Alkislar, Krothapalli, Choutapalli and Lourenco2005; Goparaju & Gaitonde Reference Goparaju and Gaitonde2018).
Figure 9. Comparison of mean streamwise velocity contours corresponding to the external shear-layer boundary (
$\bar {u}/u_{\kern-1pt j} = 0.05$
) at
$z = 0$
between RANS and PIV:
$(a)$
$s/D = 1.76$
;
$(b)$
$s/D = 3$
. The single jet (axisymmetric) RANS solution is also added for comparison.
Figure 10. Comparison of mean streamwise velocity profiles along the top nozzle axis between the RANS solution and the PIV mean flow for
$(a)$
$s/D = 1.76$
and
$(b)$
$s/D = 3$
.
Figure 10 shows the streamwise velocity evolution along the top nozzle axis for each
$s/D$
, illustrating the decay of the centreline velocity near and after the end of the potential core. The effect is intimately related to the turbulent mixing and constitutes an additional validation of the numerical solution. The reciprocal quantity (
$u_{\kern-1pt j}/\bar {u}$
) and the ratio
$\sqrt {k}/\bar {u}$
along the nozzle axis are also reported in Appendix C.
The mean-flow results demonstrate the ability of the three-dimensional (3-D) RANS solution to accurately model the twin-jet mean flow, effectively accounting for the nonlinear mean-flow interaction between the jets. This overcomes the shortcomings of more simplified models employed in previous work, such as the tailored twin-jet mean flow constructed by the linear superposition of two isolated jets (Rodríguez et al. Reference Rodríguez, Jotkar and Gennaro2018; Rodríguez Reference Rodríguez2021; Rodríguez et al. Reference Rodríguez, Stavropoulos, Nogueira, Edgington-Mitchell and Jordan2023).
5.2. Twin-jet wavepackets modelled with PM-PSE
PM-PSE calculations have been performed in the frequency range
$St = [0.1,1]$
. For each frequency, four different disturbances have been marched downstream, denoted by SS0, SA0, SS1 and SA1 following the same nomenclature as Rodríguez et al. (Reference Rodríguez, Jotkar and Gennaro2018), Nogueira & Edgington-Mitchell (Reference Nogueira and Edgington-Mitchell2021), Stavropoulos et al. (Reference Stavropoulos, Mancinelli, Jordan, Jaunet, Weightman, Edgington-Mitchell and Nogueira2023) and Rodríguez et al. (Reference Rodríguez, Stavropoulos, Nogueira, Edgington-Mitchell and Jordan2023). The first letter refers to the symmetry with respect to
$z = 0$
, the second letter to the symmetry with respect to
$y = 0$
(as described in § 3.2.2) and the number denotes whether the disturbance is of toroidal nature (equivalent to
$m = 0$
in a single round jet) or of flapping nature (equivalent to
$m = 1$
in a single round jet). Rodríguez et al. (Reference Rodríguez, Stavropoulos, Nogueira, Edgington-Mitchell and Jordan2023) discusses the nature of these modes and demonstrates that the linear coupling of the fluctuation fields of the two jets favours flapping modes over the spiral modes typically observed in single round jets. Disturbances antisymmetric with respect to the symmetry plane at
$z = 0$
(AS1 and AA1) are not treated in this work as they cannot be characterised by the current schlieren visualisations. The relative importance of these disturbances with respect to those that are symmetric across
$z = 0$
is not assessed this study. However, it is important to note that there is no limitation in the proposed wavepacket model that prevents the calculation of AS and AA mode families by means of PM-PSE. For illustrative purposes, flapping modes AS1 and AA1 have been included in figures 11 and 12 for
$St = 0.4$
and
$s/D = 1.76$
. These figures are described further later. Higher modes (twin-jet analogous for
$m \gt 1$
), although being recovered by the two-dimensional linear stability theory, are not considered as they are more challenging to educe from the experimental visualisations.
Figure 11. Isosurfaces of the real part of the pressure fluctuation for PM-PSE modes
$(a)$
SS0,
$(b)$
SS1,
$(c)$
SA0,
$(d)$
SA1,
$(e)$
AS1 and
$(f)$
AA1 at
$St = 0.4$
,
$s/D = 1.76$
. Values are normalised with respect to the maximum absolute value of the real part of
$\hat {p}$
. Displayed isosurfaces correspond to
$\textrm {Re} \{ \hat {p} \} = 0.1$
(orange) and
$\textrm {Re} \{ \hat {p} \} = -0.1$
(blue). The projected filled contours correspond to the real part of the pressure fluctuation in the
$xy$
symmetry plane located at
$z = 0$
and the
$xz$
symmetry plane located at
$y = 0$
. The colourbars refer to the projected contours.
Figure 12. Contours of the real part of the pressure fluctuation for different PM-PSE modes at
$St = 0.4$
,
$s/D = 1.76$
:
$(a{,}d{,}g{,}j{,}m{,}p)$
symmetry plane at
$z/D = 0$
;
$(b{,}e{,}h{,}k{,}n{,}q)$
symmetry plane at
$y/D = 0$
;
$(c{,}f{,}i{,}l{,}o{,}r)$
nozzle mid-plane
$y/D = s/(2D)$
;
$(a{,}b{,}c)$
mode SS0;
$(d{,}e{,}f)$
mode SS1;
$(g{,}h{,}i)$
mode SA0;
$(j{,}k{,}l)$
mode SA1;
$(m{,}n{,}o)$
mode AS1;
$(p{,}q{,}r)$
mode AA1. Values are normalised with respect to the maximum absolute value of the real part of
$\hat {p}$
in the entire field. Only one quarter of the twin-jet system is displayed according to the two inherent symmetry planes. The black dashed lines denote the nozzle lip lines.
Figure 13. Contours of the real part of the pressure fluctuation for the different PM-PSE modes at
$St = 0.4$
,
$s/D = 3$
:
$(a{,}d{,}g{,}j)$
symmetry plane at
$z/D = 0$
;
$(b{,}e{,}h{,}k)$
symmetry plane at
$y/D = 0$
;
$(c{,}f{,}i{,}l)$
nozzle mid-plane
$y/D = s/(2D)$
;
$(a{,}b{,}c)$
mode SS0;
$(d{,}e{,}f)$
mode SS1;
$(g{,}h{,}i)$
mode SA0;
$(j{,}k{,}l)$
mode SA1. Values are normalised with respect to the maximum absolute value of the real part of
$\hat {p}$
in the entire field.
Figure 11 shows the three-dimensional structure of the PM-PSE wavepackets computed for
$St = 0.4$
and
$s/D = 1.76$
(closely spaced jets). For each mode, isosurfaces of the real part of the pressure fluctuation are shown, together with projected filled contours that correspond to the two symmetry planes at (
$y = 0$
) and (
$z = 0$
). The isosurfaces clearly illustrate the toroidal nature of modes SS0 and SA0, characterised by ring-like structures that grow downstream within the potential core and the mixing layers of each jet, while the projected contours serve to highlight the near pressure field associated with the wavepacket structure. In the range
$St = [0.4, 1]$
, the PSE wavepackets clearly exhibit the distinct acoustic beam of Mach-wave radiation; towards the lower
$St$
, the acoustic beam becomes gradually less evident. This can be attributed to different causes: (i) for lower
$St$
, the wavepackets may gradually produce a weaker Mach-wave radiation owing to a slower convection speed or elongated envelope; or (ii) the inherent limitations of PSE in capturing the Mach-wave radiation for decreasing
$St$
, as demonstrated by Towne et al. (Reference Towne, Rigas and Colonius2019). The three-dimensional structure of modes SS1 and SA1 reveals flapping fluctuations featuring a change of sign across the axis of each jet.
To better highlight some key features of the three-dimensional structure of the twin-jet wavepackets, figures 12 and 13 present filled contours of the real part of the pressure fluctuation at
$St = 0.4$
along three different cutplanes for
$s/D = 1.76$
and
$s/D = 3$
, respectively. In addition to both twin-jet symmetry planes, the
$xz$
plane located at
$y = s/(2D)$
(nozzle axis) is also shown for comparison. For both nozzle spacings, mode SS0 features similar amplitude levels and Mach-wave radiation signatures in both the
$z = 0$
and
$y = s/(2D)$
planes. In contrast, for mode SA0, while the wavepacket radiates in the
$z = 0$
plane up to the end of the displayed streamwise range (
$x/D = 10$
), no distinct acoustic beam is visible in the
$y = s/(2D)$
plane. Similarly, the flapping modes do not exhibit clear near-field signatures associated with acoustic radiation at this particular frequency.
For closely spaced jets (
$s/D = 1.76$
), significant differences are found in the streamwise evolution of the wavepackets between the inner and the outer mixing layers of each jet. This is clearly visible in the
$z/D = 0$
contour plots (figure 12
$a{,}d{,}g{,}j$
). For the toroidal fluctuations, this makes the ring-like structures appear tilted with respect to the axis of the jet, rather than perpendicular to it as would be expected in the axisymmetric case. For the flapping modes, it causes the fluctuations to not be in perfect counter-phase across the nozzle axis, preventing the amplitude at the jet axis from reaching a zero value. As shown later, these phase differences are associated with a mismatch in the phase-speed evolution
$c_{\textit{ph}}$
of the Kelvin–Helmholtz waves along the inner and outer mixing layers, which is attributed to the mean-flow interaction taking place in the region between the two jets.
To quantify the disparity between the part of the disturbances respectively evolving through the inner and outer mixing layers, figure 14 shows the phase difference in the pressure fluctuation of the four PM-PSE modes between the inner and outer lip lines of each nozzle (represented by the black dashed lines in the
$z/D=0$
plane in figures 12 and 13). Along each lip line
$y_{o,i}/D = s/(2D) \pm 1/2$
, where subscripts
$o$
and
$i$
refer to the outer and inner lip lines, the phase is defined as
$\phi _{o,i}(x) = \arg (\hat {p}_\omega (x,y_{o,i}))$
(
$\hat {p}_\omega$
is the fluctuation field in (3.10)). The phase distributions are referenced to their respective phase at the initial position of the PM-PSE marching (
$x_0/D = 0.75$
), so that both
$\phi _o$
and
$\phi _i$
have a value of zero at that point. For
$s/D = 1.76$
, all four modes accumulate a notable phase difference near the end of the potential core (
$x/D \approx 10$
), with mode SA1 exhibiting the largest phase shift between the inner and outer lip lines, which exceeds
$80^\circ$
.
Figure 15 presents the streamwise evolution of the phase speed of the pressure fluctuation for each mode along the lip lines. The phase velocity is computed as
$c_{\textit{ph}}(x) = \omega / k_x$
, where
$k_x$
is the streamwise derivative of
$\phi _{o,i}$
. The comparison between the phase-speed evolution along each line reveals a higher phase speed in the inner lip line with respect to the outer one. For modes SS0 and SS1, the phase-speed evolution correlates with the streamwise mean-flow velocity evolution at these positions. In particular, SS0 and SS1 exhibit an important phase-speed mismatch for
$x/D \gt 5$
, which is the streamwise range in which the mean-flow velocity along the inner lip line becomes significantly higher than along the outer lip line. Therefore, the higher streamwise velocity found in the inner mixing layers leads to a larger convective speed for the wavepackets in this region. However, the antisymmetric modes SA0 and SA1 present a significant discrepancy in
$c_{\textit{ph}}$
between both lines for the entire streamwise range. This is attributed to the combination of two effects: the higher mean-flow velocity in the inner mixing layer, as for the symmetric modes, but also the fact that the antisymmetric wavepacket structures contained in the inner mixing layer are confined between the nozzle axis and the symmetry plane. For low
$St$
, the spatial wavelength of the fluctuations becomes comparable to the spacing between the nozzle axis and the symmetry axis. This has an impact in the development of the antisymmetric structures, modulating the wavenumber in the inner mixing layer such that the phase speed is increased in this case.
The phase-speed mismatch between the inner and outer mixing layers is found to be much weaker for higher frequencies (
$St \gt 0.4$
), as the radial reach of the wavepacket structures across the plane between jets decreases with increasing frequency. For lower frequencies (
$St \leqslant 0.4$
), the antisymmetric flapping modes (SA1) remain the fluctuations which manifest the largest distortion owing to the phase-speed difference, in line with the results presented in figure 14
$(a)$
. More details on the evolution of SA1 modes with
$St$
are provided in § 5.4.3.
In the case of larger jet spacing (
$s/D = 3$
), the phase-speed mismatch and associated phase difference between the inner and outer mixing layers is much weaker, as shown in figures 14
$(b)$
and 15
$(e$
–
$h)$
. In this configuration, the structure of the modes is closer to the axisymmetric, single-jet behaviour, and the structure of the wavepackets features a perpendicular alignment with respect to the nozzle axis (see figure 13).
Figure 14. Streamwise evolution of the phase difference between the outer and inner lip lines of the pressure fluctuation for the different PM-PSE modes (
$St = 0.4$
):
$(a)$
$s/D = 1.76$
;
$(b)$
$s/D = 3$
.
$\phi _o$
and
$\phi _i$
respectively denote the phase along the outer and inner lip lines. The phase along each line is referenced to the initial streamwise station where the PSE marching is initialised.
Figure 15. Streamwise evolution of the phase speed of the pressure fluctuation along the outer and inner lip lines for the different PM-PSE modes at
$St = 0.4$
:
$(a$
–
$d)$
$s/D = 1.76$
;
$(e$
–
$h)$
$s/D = 3$
;
$(a{,}e)$
SS0;
$(b{,}f)$
SS1;
$(c{,}g)$
SA0;
$(d{,}h)$
SA1.
5.3. Coherent structures educed from the schlieren measurements
Figure 16. SPOD spectra obtained using the cross-spectral density of
$\varTheta$
fluctuations:
$(a{,}c)$
$s/D = 1.76$
;
$(b{,}d)$
$s/D = 3$
;
$(a{,}b)$
symmetric datasets;
$(c{,}d)$
antisymmetric datasets. Only the first ten SPOD modes are shown, ranked following a grey scale between mode 1 (black) and mode 10 (white).
$\lambda$
represents the spectral density (eigenvalue of the cross-spectral density matrix) associated to each SPOD mode for each frequency. The orange line denotes the sum of the spectral density of all 57 SPOD modes.
Figure 16 shows the SPOD spectra obtained for the symmetric and antisymmetric datasets for both jet separations, based on the cross-spectral density of the
$\varTheta$
fluctuations (see § 4.1). For all cases, the spectrum presents a broadband structure with highest energy contained in the range
$St \approx [0.1,1]$
, corresponding to the signature of the mixing noise. The decomposition is low-rank for most of the frequencies of interest, as indicated by the proximity of the line corresponding to mode 1 to the orange line, which represents the sum of the energy for all the modes. For frequencies below
$St = 0.6$
, more than 50 % of the total energy is contained in the first SPOD mode.
The structure of the first and second symmetric SPOD modes educed for
$St = 0.4$
and
$s/D = 1.76$
is shown in figures 17
$(a{,}c)$
and 18
$(a{,}c)$
. Both the
$\varTheta$
fluctuation field as well as the schlieren fluctuation field obtained from the SPOD problem based on
$\varTheta$
, denoted by
$\hat {\sigma }_{\varTheta }$
, are shown for each mode. SPOD mode 1 exhibits coherent structures that largely resemble toroidal fluctuations, while SPOD mode 2 features a clear flapping structure with low amplitude at the nozzle axes. The decomposition separates the two types of fluctuation into different SPOD modes, facilitating comparison with the PM-PSE models. The efficacy of the decomposition in separating toroidal and flapping modes is favoured by the fact that the two types of structure are naturally largely orthogonal to one another. In addition, both SPOD modes contain the associated Mach-wave radiation signature. These properties are observed over the entire frequency range of interest.
The peaks seen in the SPOD spectra at
$St \approx 0.55$
and
$St \approx 1.15$
are anomalies in the experimental data whose nature has not yet been clarified. Different tests have been conducted to shed light on their origin, in particular: schlieren measurements at different
$M_{\kern-1pt j}$
, different jet spacing and different total temperature, as well as microphone measurements at different
$M_{\kern-1pt j}$
. The results show that these peaks are not sensitive to the nozzle pressure ratio (or equivalently, to
$M_{\kern-1pt j}$
) or the nozzle spacing, and they are consistently found in both schlieren and microphone measurements taken at different times. Their frequency, however, is found to be sensitive to the jet total temperature. Their frequency ratio is found to be approximately 2.1. This might indicate that the second peak is an harmonic of the first one, as they are not tonal but centred around a small frequency band. At present, it is hypothesised that these anomalies are generated by structural vibrations in the twin-jet system, which weakly force the jets at those specific frequencies.
5.4. Comparison between PM-PSE wavepackets and SPOD modes
The comparison of PM-PSE wavepackets with the coherent structures educed from the schlieren measurements requires the calculation of a numerical schlieren field for the PM-PSE fluctuations. The corresponding PM-PSE schlieren fluctuation field is computed as
\begin{align} \hat {\sigma } (x,y) = \int _{-\infty }^{\infty } \frac {\partial \hat {\rho } (x,y,z)}{\partial x} \,\mathrm{d} z, \end{align}
where the density fluctuation is calculated from
$\hat {p}$
and
$\hat {T}$
through the linearised perfect gas equation of state, and the streamwise derivative of
$\hat {\rho }$
is evaluated by differentiating the PM-PSE ansatz (3.9) with respect to
$x$
. The PM-PSE
$\hat {\varTheta }$
fluctuation fields are then computed following the same procedure employed for the experimental datasets, i.e. solving the Poisson equation (4.3) in the spectral domain including a phase-speed filter for spurious waves that propagate with unphysically large supersonic speed (Padilla-Montero et al. Reference Padilla-Montero, Rodríguez, Jaunet and Jordan2024).
To quantify the structural similarity between the SPOD and PM-PSE modes, the following alignment factor is calculated as the normalised projection of the schlieren SPOD eigenvectors on the corresponding PM-PSE schlieren perturbations for a given frequency, as done in previous studies of single jets (Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Cavalieri et al. Reference Cavalieri, Rodríguez, Jordan, Colonius and Gervais2013; Sinha et al. Reference Sinha, Rodríguez, Brès and Colonius2014; Sasaki et al. Reference Sasaki, Cavalieri, Jordan, Schmidt, Colonius and Brès2017):
\begin{align} \chi _{\sigma _\varTheta } = \frac {\big|\big\langle \hat {\sigma }_{\varTheta }, \hat {\sigma }_{\textit{PSE}} \big\rangle \big|}{\big\langle \hat {\sigma }_{\varTheta }, \hat {\sigma }_{\varTheta } \big\rangle ^{1/2} \big\langle \hat {\sigma }_{\textit{PSE}}, \hat {\sigma }_{\textit{PSE}} \big\rangle ^{1/2}}, \end{align}
where the inner product between two fields
$\langle \hat {q}_1, \hat {q}_2 \rangle$
is defined as
\begin{align} \langle \hat {q}_1, \hat {q}_2 \rangle = \iint _\varOmega \hat {q}_1^* \hat {q}_2 \,\mathrm{d}x\, \mathrm{d}y, \end{align}
with
$\varOmega$
denoting the spatial domain over which the projection is evaluated, in this case, limited by the spatial window of the schlieren images. To avoid including the dark region corresponding to the mirror of the optical schlieren set-up, a horizontal window between
$x/D = 0.75$
and
$x/D = 8.3$
is employed here for the calculation of the inner product. Alternatively to
$\chi _{\sigma _\varTheta }$
, an alignment coefficient based on
$\hat {\varTheta }$
is also evaluated, given by
\begin{align} \chi _{\varTheta } = \frac {\big|\big\langle \hat {\varTheta }_{\textit{SPOD}}, \hat {\varTheta }_{\textit{PSE}} \big\rangle \big|}{\big\langle \hat {\varTheta }_{\textit{SPOD}}, \hat {\varTheta }_{\textit{SPOD}} \big\rangle ^{1/2} \big\langle \hat {\varTheta }_{\textit{PSE}}, \hat {\varTheta }_{\textit{PSE}} \big\rangle ^{1/2}}. \end{align}
By definition, the alignment factor takes values between 0 and 1. A value of 0 means that the fields compared are spatially orthogonal, while a value of unity means that they are parallel, and differ at most by a complex scale factor.
5.4.1. Small jet separation,
$s/D = 1.76$
Figure 17. Contours of the real part of the symmetric toroidal fluctuation for
$St = 0.4$
and
$s/D = 1.76$
:
$(a)$
schlieren field of SPOD mode 1 (
$\hat {\varTheta }$
-based CSD);
$(b)$
schlieren field of PM-PSE mode SS0;
$(c)$
$\hat {\varTheta }$
field of SPOD mode 1;
$(d)$
$\hat {\varTheta }$
field of PM-PSE mode SS0. Dashed lines indicate the nozzle lip lines.
Figure 18. Contours of the real part of the symmetric flapping fluctuation for
$St = 0.4$
and
$s/D = 1.76$
:
$(a)$
schlieren field of SPOD mode 2;
$(b)$
schlieren field of PM-PSE mode SS1;
$(c)$
$\hat {\varTheta }$
field of SPOD mode 2;
$(d)$
$\hat {\varTheta }$
field of PM-PSE mode SS1.
Figure 19. Alignment factors for
$s/D = 1.76$
:
$(a{,}b)$
alignment between schlieren fluctuation fields
$\hat {\sigma }$
;
$(c{,}d)$
between
$\hat {\varTheta }$
fields;
$(a{,}c)$
symmetric fluctuations with respect to
$xz$
plane;
$(b{,}d)$
antisymmetric fluctuations with respect to
$xz$
plane.
First, the case of closely spaced twin jets (
$s/D = 1.76$
) is discussed. A qualitative comparison between SPOD and PM-PSE fluctuations for symmetric modes is presented in figures 17 and 18. These figures show filled contour plots of the schlieren fluctuation field
$\hat {\sigma }$
as well as the
$\hat {\varTheta }$
fluctuation field for both SPOD and PM-PSE modes at
$St = 0.4$
. Figure 17 compares SPOD mode 1 against the PM-PSE mode SS0, i.e. the symmetric toroidal fluctuation. The experimentally educed structure and the modelled wavepacket are in excellent agreement, showing ring-like structures tilted with respect to the nozzle axis owing to the phase-speed mismatch between the inner and outer mixing layers, as discussed in § 5.2. The two representations of the fluctuation exhibit spatial wavenumber distributions that are in qualitative accordance both inside and outside of the jets, indicating that the Mach-wave radiation is also successfully captured by the PM-PSE model.
Figure 18 compares SPOD mode 2 with PM-PSE mode SS1. The results are also in good agreement in this case, with both representations showing symmetric flapping structures with similar wavenumbers and Mach-wave radiation patterns. The PM-PSE mode, however, does not yield a zero amplitude region at the nozzle axis. This is attributed to the phase-speed mismatch observed between the inner and outer mixing layers, which prevents the flapping mode structures from being perfectly in counter-phase with each other across the axis. In view of the results presented in § 5.2, this is considered a physical characteristic of closely spaced twin jets. The SPOD algorithm, however, extracts a flapping structure which resembles more closely the
$m = 1$
fluctuation encountered in single jets, featuring a chequerboard pattern with nearly zero amplitude on the nozzle axis. This discrepancy between the SPOD structure and the model, rather than implying a deficiency of PM-PSE in modelling flapping wavepackets, reflects a limitation of educing three-dimensional coherent structures from two-dimensional experimental images. The lack of information in the third dimension prevents SPOD from extracting a fluctuation that incorporates the effect of the line-of-sight integrated phase-speed mismatch.
Figure 20. Contours of the real part of the antisymmetric fluctuations for
$St = 0.4$
and
$s/D = 1.76$
:
$(a)$
schlieren field of SPOD mode 1;
$(b)$
schlieren field of PM-PSE mode SA0;
$(c)$
schlieren field of SPOD mode 2;
$(d)$
schlieren field of PM-PSE mode SA1;
$(e)$
$\hat {\varTheta }$
field of SPOD mode 1;
$(f)$
$\hat {\varTheta }$
field of PM-PSE mode SA0;
$(g)$
$\hat {\varTheta }$
field of SPOD mode 2;
$(h)$
$\hat {\varTheta }$
field of PM-PSE mode SA1.
Two alignment factors quantifying the agreement between PM-PSE and the experimentally educed structures for
$s/D = 1.76$
are shown in figure 19. Alignment factors are computed both for schlieren fluctuations (denoted by
$\chi _{\sigma _\varTheta }$
) and for
$\hat {\varTheta }$
fluctuations (denoted by
$\chi _\varTheta$
), as well as for symmetric and antisymmetric modes. For each case, four different alignments are evaluated:
-
(i) first SPOD mode against the toroidal PM-PSE mode, (SPOD1,SS0) and (SPOD1,SA0);
-
(ii) second SPOD mode against the toroidal PM-PSE mode, (SPOD2,SS0) and (SPOD2,SA0);
-
(iii) first SPOD mode against the flapping PM-PSE mode, (SPOD1,SS1) and (SPOD1,SA1);
-
(iv) second SPOD mode against the flapping PM-PSE mode, (SPOD2,SS1) and (SPOD2,SA1).
According to the splitting of toroidal and flapping fluctuations into SPOD modes 1 and 2, respectively, alignments (SPOD1,SS0) and (SPOD2,SS1) are expected to yield much higher values than (SPOD2,SS0) and (SPOD1,SS1). This is what is observed in figure 19
$(a{,}c)$
, except for
$St = 0.1$
, for which the streamwise wavelength of the wavepacket structure is comparable to the size of the schlieren window and the resulting SPOD modes are not accurate representations. The alignment factors based on
$\sigma _\varTheta$
are significantly lower than those based on
$\hat {\varTheta }$
due to the presence of small-scale vortical structures in the fluctuation fields. The obtained alignment factors for both symmetric fluctuations are high over most of the studied frequency range. Particularly in the interval
$St = [0.3, 0.7]$
, values of
$\chi _{\varTheta } \approx 0.97$
(
$\chi _{\sigma _\varTheta } \approx 0.85$
) are reached for toroidal fluctuations and
$\chi _{\varTheta } \approx 0.95$
(
$\chi _{\sigma _\varTheta } \approx 0.83$
) for flapping fluctuations, indicating an excellent agreement between the experimentally educed coherent structures and the PM-PSE wavepacket models. According to the SPOD spectra shown in figure 16, this is the range of frequencies at which the coherent structures have highest energy.
The alignment factors for the antisymmetric fluctuations are shown in figure 19
$(b{,}d)$
. A clear dominance of alignment factors (SPOD1,SA0) and (SPOD2,SA1) is also observed (as for the symmetric structures) for frequencies above
$St = 0.4$
. An anomalous behaviour is found at lower frequencies for which both PM-PSE modes SA0 and SA1 have a high alignment against SPOD mode 1. This is discussed in more detail in § 5.4.3. Figures 20 and 21 compare the first two antisymmetric SPOD modes against the PM-PSE modes SA0 and SA1 for
$St = 0.4$
and
$St = 0.6$
, respectively. For toroidal fluctuations, the comparison at these frequencies reveals the same findings as for the symmetric case, namely, a good agreement for mid-range frequencies, which progressively deteriorates for higher and lower frequencies. For antisymmetric flapping fluctuations, their distortion, owing to the phase-speed mismatch between the inner and outer mixing layers, is found to be stronger than for the symmetric counterpart. This is supported by figure 14
$(a)$
, where mode SA1 has the highest accumulated phase difference. This distortion reduces the agreement with SPOD modes, and shifts the region of optimal alignment towards higher frequencies, namely,
$St = 0.6$
as opposed to
$St=0.4$
as for the symmetric case. The agreement between antisymmetric fluctuations at
$St = 0.6$
is excellent, as shown in figure 21.
5.4.2. On the sources of misalignment between PM-PSE and SPOD modes
Figure 21. Contours of the real part of the antisymmetric fluctuations for
$St = 0.6$
and
$s/D = 1.76$
:
$(a)$
schlieren field of SPOD mode 1;
$(b)$
schlieren field of PM-PSE mode SA0;
$(c)$
schlieren field of SPOD mode 2;
$(d)$
schlieren field of PM-PSE mode SA1;
$(e)$
$\hat {\varTheta }$
field of SPOD mode 1;
$(f)$
$\hat {\varTheta }$
field of PM-PSE mode SA0;
$(g)$
$\hat {\varTheta }$
field of SPOD mode 2;
$(h)$
$\hat {\varTheta }$
field of PM-PSE mode SA1.
Figure 22. Contours of the real part of the symmetric fluctuations for
$St = 0.2$
and
$s/D = 1.76$
:
$(a)$
schlieren field of SPOD mode 1;
$(b)$
schlieren field of PM-PSE mode SS0;
$(c)$
schlieren field of SPOD mode 2;
$(d)$
schlieren field of PM-PSE mode SS1;
$(e)$
$\hat {\varTheta }$
field of SPOD mode 1;
$(f)$
$\hat {\varTheta }$
field of PM-PSE mode SS0;
$(g)$
$\hat {\varTheta }$
field of SPOD mode 2;
$(h)$
$\hat {\varTheta }$
field of PM-PSE mode SS1.
For frequencies above
$St = 0.7$
and below
$St = 0.3$
, the alignment factors for the symmetric fluctuations (SPOD1,SS0) and (SPOD2,SS1) progressively decrease, although values
$\chi _{\sigma _\varTheta } \geqslant 0.7$
are maintained throughout most of the studied frequency range. To illustrate the differences between the structures at other frequencies for which the alignment is lower, figures 22 and 23 compare the symmetric SPOD modes and PM-PSE modes for
$St = 0.2$
and
$St = 0.8$
, respectively.
For low frequencies (e.g.
$St = 0.2$
), the deterioration of the agreement is attributed to the following observations. First, the streamwise wavelength of the PM-PSE wavepackets appears to be larger than that of the corresponding SPOD modes, especially for the toroidal fluctuations. This might reflect a limitation of the PSE model at low frequencies. Past studies (Suzuki & Colonius Reference Suzuki and Colonius2006; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Cavalieri et al. Reference Cavalieri, Rodríguez, Jordan, Colonius and Gervais2013; Sinha et al. Reference Sinha, Rodríguez, Brès and Colonius2014) have reported difficulties in modelling low-frequency wavepackets in single jets by means of PSE. Recent investigations based on resolvent analysis (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Lesshafft et al. Reference Lesshafft, Semeraro, Jaunet, Cavalieri and Jordan2019; Pickering et al. Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020) have shown that non-modal effects are important for the linear dynamics at low
$St$
and
$m = 0$
, and these cannot be captured by the PSE.
Figure 23. Contours of the real part of the symmetric fluctuations for
$St = 0.8$
and
$s/D = 1.76$
:
$(a)$
schlieren field of SPOD mode 1;
$(b)$
schlieren field of PM-PSE mode SS0;
$(c)$
schlieren field of SPOD mode 2;
$(d)$
schlieren field of PM-PSE mode SS1;
$(e)$
$\hat {\varTheta }$
field of SPOD mode 1;
$(f)$
$\hat {\varTheta }$
field of PM-PSE mode SS0;
$(g)$
$\hat {\varTheta }$
field of SPOD mode 2;
$(h)$
$\hat {\varTheta }$
field of PM-PSE mode SS1.
Second, the streamwise range of linear growth and development of coherent structures for low
$St$
is larger than for higher frequencies, which reach amplitude saturation earlier upstream. Given the limited streamwise length of the schlieren measurement window, the coherent fluctuations educed by SPOD at low
$St$
exhibit an organised structure in a smaller portion of the domain. This is more evident for the symmetric flapping structure in SPOD2 at
$St = 0.2$
(figure 22
$c{,}g$
), which does not emerge as a clearly organised structure until
$x/D \gt 5$
.
Third, the wavelength of the fluctuations for low
$St$
becomes comparable to the size of the schlieren window. In this circumstance, the coherent structures extracted by means of SPOD are influenced by the domain truncation enforced by the size of the schlieren images. In particular, the size of the domain along
$y$
limits the eduction of Mach-wave radiation after a certain streamwise position, which for e.g.
$St = 0.2$
, is found to be very early upstream (
$x/D \approx 2$
). The last two shortcomings are inherent to the methodology employed here for extracting coherent structures from the experimental data. In addition to the aforementioned limitations, the PM-PSE flapping modes for low
$St$
are also highly distorted owing to the effect of the phase-speed mismatch described in § 5.2. For the reasons already discussed there, the SPOD modes do not capture this mismatch accurately.
For high frequencies (
$St \geqslant 0.7$
), the small decrease in the alignment factors is mainly attributed to limitations of the linear PM-PSE model. As an illustrative example, for
$St = 0.8$
(see figure 23
$f{,}h$
), the PM-PSE wavepacket grows rapidly after the nozzle exit up to
$x/D \approx 5$
and then progressively vanishes further downstream. The corresponding educed SPOD structures, in contrast, feature energetic coherent structures within the potential core extending up to the downstream boundary of the visualisation window, without a decrease in the amplitude. Previous works on single jets have also observed similar discrepancies between the predictions of the linear stability approaches limited to modal instabilities and experimental measurements in the downstream region (Suzuki & Colonius Reference Suzuki and Colonius2006; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Cavalieri et al. Reference Cavalieri, Rodríguez, Jordan, Colonius and Gervais2013; Sinha et al. Reference Sinha, Rodríguez, Brès and Colonius2014; Breakey et al. Reference Breakey, Jordan, Cavalieri, Nogueira, Léon, Colonius and Rodríguez2017). Studies based on a spatial model of the Orr mechanism (Tissot et al. Reference Tissot, Lajús, Cavalieri and Jordan2017), spatial transient-growth calculations (Jordan et al. Reference Jordan, Zhang, Lehnasch and Cavalieri2017) and resolvent analysis (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Lesshafft et al. Reference Lesshafft, Semeraro, Jaunet, Cavalieri and Jordan2019) show that non-modal mechanisms gradually become important in this region, and suggest that nonlinear interactions are key in their activation (Jordan et al. Reference Jordan, Zhang, Lehnasch and Cavalieri2017). The PM-PSE wavepackets exhibit an additional difference: the amplitude of the Mach-wave radiation relative to the amplitude of the wavepacket evolving inside the jet is found to be weaker for the PM-PSE modes than for SPOD modes. This discrepancy has also been observed in previous studies dealing with isolated jets (Sinha et al. Reference Sinha, Rodríguez, Brès and Colonius2014; Breakey et al. Reference Breakey, Jordan, Cavalieri, Nogueira, Léon, Colonius and Rodríguez2017) and has been attributed to inherent limitations of the parabolised stability equations (Towne et al. Reference Towne, Rigas and Colonius2019).
The aforementioned effects have a small impact for the intermediate frequency range,
$St \approx [0.3, 0.7]$
, which corresponds to the most energetic part of the spectrum according to the experimental data (see figure 16). Large alignment factors between the PM-PSE model and the experimentally educed structures are obtained for these frequencies. This result is consistent with the dominance of modal, Kelvin–Helmholtz instabilities in the intermediate frequency range highlighted by Pickering et al. (Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020) for single jets.
5.4.3. On the inconsistent agreement for antisymmetric flapping modes at low
$St$
Figure 24. Contours of the real part of the antisymmetric fluctuations for
$St = 0.2$
and
$s/D = 1.76$
:
$(a)$
schlieren field of SPOD mode 1;
$(b)$
schlieren field of PM-PSE mode SA0;
$(c)$
schlieren field of SPOD mode 2;
$(d)$
schlieren field of PM-PSE mode SA1;
$(e)$
$\hat {\varTheta }$
field of SPOD mode 1;
$(f)$
$\hat {\varTheta }$
field of PM-PSE mode SA0;
$(g)$
$\hat {\varTheta }$
field of SPOD mode 2;
$(h)$
$\hat {\varTheta }$
field of PM-PSE mode SA1.
Figure 25. Contours of the real part of the PM-PSE antisymmetric flapping fluctuations (SA1) for various
$St$
(
$s/D = 1.76$
):
$(a{,}d{,}g{,}j{,}m)$
pressure fluctuation at
$z = 0$
;
$(b{,}e{,}h{,}k{,}n)$
density fluctuation at
$z = 0$
;
$(c{,}f{,}i{,}l{,}o)$
schlieren fluctuation;
$(a{,}b{,}c)$
$St = 0.1$
;
$(d{,}e{,}f)$
$St = 0.2$
;
$(g{,}h{,}i)$
$St = 0.3$
;
$(j{,}k{,}l)$
$St = 0.4$
;
$(m{,}n{,}o)$
$St = 0.5$
.
Figure 26. Alignment factor between the schlieren (
$\sigma$
) and
$\varTheta$
perturbation fields of PSE modes SA0 and SA1 as a function of
$St$
:
$(a{,}b)$
$s/D = 1.76$
;
$(c{,}d)$
$s/D = 3$
;
$(a{,}c)$
alignment between
$\hat {\sigma }$
fields;
$(b,d)$
alignment between
$\hat {\varTheta }$
fields. For comparison, the alignment between modes SS0 and SS1 is also included.
The alignment factor for antisymmetric flapping modes at low frequencies (
$St \lt 0.4$
, see figure 19
$b{,}d$
) shows an inconsistent behaviour in which both PM-PSE modes SA0 and SA1 have a high alignment with SPOD mode 1. To illustrate this, figure 24 compares empirical (SPOD) and modelled (PM-PSE) antisymmetric fluctuations for
$St = 0.2$
. The PM-PSE schlieren fluctuation and the
$\varTheta$
fluctuation fields in this case better resemble a toroidal fluctuation than a flapping one. The structure of the line-of-sight integrated fields for SA1 at this frequency is similar to that of mode SA0 and to that of mode SPOD1. As a result, both the (SPOD1,SA0) and (SPOD1,SA1) factors yield similar alignments. In addition, the educed coherent structure in SPOD mode 2 is poorly organised for most of the measurement window at this frequency, leading to small values of
$\chi _{\sigma _\varTheta }$
for the projection (SPOD2,SA1).
The explanation of this phenomenon is linked to the deformation suffered by the antisymmetric flapping fluctuations when the jets are closely spaced and the value of
$St$
is low. This deformation is, on one side, caused by the phase-speed mismatch induced by the mean-flow differences between the inner and outer mixing layers of the twin-jet system, as described in previous sections. On the other side, flapping antisymmetric modes modelled by PM-PSE are confined between the nozzle axis and the symmetry plane at
$y = 0$
(where antisymmetry boundary conditions are imposed). Under the combination of low
$St$
and small jet spacing, the resulting Kelvin–Helmholtz wavelength is large compared with the distance between the nozzle axis and the symmetry axis at
$y = 0$
. As a consequence, the flapping wavepacket structures for
$s/(2D) \gt y/D \gt 0$
become squeezed in the vertical direction and elongated in the streamwise direction, leading to a phase-speed increase in this region and to fluctuations which, after line-of-sight integration, adopt a shape such as the one depicted in figure 24
$(d)$
.
To support this argument, figure 25 represents the evolution of PM-PSE mode SA1 as a function of
$St$
. This figure shows contours of the PM-PSE pressure fluctuation, density fluctuation and schlieren fluctuation for
$St = [0.1,0.5]$
in the
$xy$
plane at
$z = 0$
. As
$St$
is progressively decreased, the wavelength of the wavepackets becomes larger, the phase shift between the flapping structures above and below the nozzle axis deviates from the value of
$\pi$
radians for an exact
$m=1$
, and the shape of the structures below the nozzle axis progressively becomes more distorted. When looking at the corresponding schlieren PSE fluctuations, the chequerboard pattern characteristic of pure flapping motions progressively fades as
$St$
decreases and the line-of-sight integrated wavepacket signature shifts to structures that appear closer to toroidal motions tilted at an angle with respect to the
$y$
axis. To further illustrate this, the normalised projection of the schlieren field
$\hat {\sigma }$
and the
$\hat {\varTheta }$
field of PM-PSE mode SA0 onto the PM-PSE mode SA1 is shown in figure 26
$(a{,}b)$
. For
$St \lt 0.6$
, the alignment between the schlieren fields of modes SA0 and SA1 becomes progressively high, indicating that both structures resemble each other. This effect is much weaker for the projection between modes SS0 and SS1, which remains significantly lower for the entire frequency range.
Although the SPOD is able to decompose the experimental data into structures that strongly resemble toroidal and flapping fluctuations in the leading modes, there is no guarantee that the decomposition would strictly separate the structures into the actual,
$m$
-like twin-jet modes (SA0, SA1, SA2 etc.). Instead, it tends to build modes that are spatially orthogonal within the defined inner product. As a consequence, the individual SPOD modes may fail to capture the subtleties of the phase-speed mismatch and shape distortion observed in the toroidal and flapping antisymmetric PM-PSE modes, features which deteriorate the orthogonality and that, in turn, may be spread over different SPOD modes. This limitation associated with the use of schlieren (two-dimensional) visualisations to educe the empirical structures leads to the anomalous behaviour in the alignment factors reported in figure 19
$(b{,}d)$
.
5.4.4. Moderate jet separation,
$s/D = 3$
Figure 27. Alignment factors for
$s/D = 3$
:
$(a{,}b)$
alignment between schlieren fluctuation fields
$\hat {\sigma }$
;
$(c{,}d)$
between
$\hat {\varTheta }$
fields;
$(a{,}c)$
symmetric fluctuations with respect to
$xz$
plane;
$(b{,}d)$
antisymmetric fluctuations with respect to
$xz$
plane.
Figure 28. Contours of the real part of the symmetric toroidal fluctuation for
$St = 0.4$
and
$s/D = 3$
:
$(a)$
schlieren field of SPOD mode 1;
$(b)$
schlieren field of PM-PSE mode SS0;
$(c)$
$\hat {\varTheta }$
field of SPOD mode 1;
$(d)$
$\hat {\varTheta }$
field of PM-PSE mode SS0.
Figure 29. Contours of the real part of the antisymmetric toroidal fluctuation for
$St = 0.4$
and
$s/D = 3$
:
$(a)$
schlieren field of SPOD mode 1;
$(b)$
schlieren field of PM-PSE mode SA0;
$(c)$
$\hat {\varTheta }$
field of SPOD mode 1;
$(d)$
$\hat {\varTheta }$
field of PM-PSE mode SA0.
Figure 30. Contours of the real part of the symmetric flapping fluctuation for
$St = 0.4$
and
$s/D = 3$
:
$(a)$
schlieren field of SPOD mode 2;
$(b)$
schlieren field of PM-PSE mode SS1;
$(c)$
$\hat {\varTheta }$
field of SPOD mode 2;
$(d)$
$\hat {\varTheta }$
field of PM-PSE mode SS1.
Figure 31. Contours of the real part of the antisymmetric flapping fluctuation for
$St = 0.4$
and
$s/D = 3$
:
$(a)$
schlieren field of SPOD mode 2;
$(b)$
schlieren field of PM-PSE mode SA1;
$(c)$
$\hat {\varTheta }$
field of SPOD mode 2;
$(d)$
$\hat {\varTheta }$
field of PM-PSE mode SA1.
The comparison between SPOD and PM-PSE for
$s/D = 3$
is reported in figures 27–31. Figure 27 presents the alignment factors, which show good agreement between SPOD1 and toroidal PM-PSE modes, and between SPOD2 and flapping PM-PSE modes. For this jet spacing, the anomaly in the alignment induced by the deformation of antisymmetric flapping modes manifests only at very low frequencies (
$St \lt 0.2$
; see also figure 26
$c{,}d$
), which is consistent with the fact that the larger jet spacing does not constrain the development of antisymmetric flapping structures over a longer range of wavelengths. Similarly, since the phase-speed mismatch between the inner and outer mixing layers is much smaller for this nozzle spacing (see figures 14
$b$
and 15
$e$
–
$h$
), the distortion of symmetric flapping structures is also considerably weaker, and the evolution of the alignment factor as a function of
$St$
is close to the evolution of the alignment for toroidal modes. These observations are supported by figures 28 and 29, which show the symmetric and antisymmetric comparisons for toroidal fluctuations at
$St = 0.4$
, as well as figures 30 and 31, which display the respective comparison for flapping fluctuations. For
$s/D=3$
, the toroidal structures are almost perpendicular to the nozzle axes, and the modelled flapping wavepackets present a well-defined chequerboard pattern, respectively resembling the
$m = 0$
and
$m = 1$
modes found in isolated round jets. This reflects how, for
$s/D = 3$
, the interaction between jets is much weaker than for
$s/D = 1.76$
, with each jet behaving closer to an isolated axisymmetric case.
The magnitude of the alignment factors for
$s/D = 3$
is not significantly higher than for
$s/D = 1.76$
, especially near the best alignment frequency
$St \approx 0.4$
. This is a quantitative indicator that the PM-PSE wavepacket models for
$s/D = 1.76$
are as accurate as for
$s/D = 3$
and, therefore, are able to incorporate the impact of the mean-flow interactions occurring when the jets are closely spaced.
The alignment factors obtained for both values of
$s/D$
are comparable to those reported in previous studies dealing with single round jets (see for example Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Cavalieri et al. Reference Cavalieri, Rodríguez, Jordan, Colonius and Gervais2013; Sinha et al. Reference Sinha, Rodríguez, Brès and Colonius2014), demonstrating the ability of plane-marching PSE to successfully model wavepackets in perfectly expanded supersonic twin jets for a significant portion of the frequency spectrum, and especially for the most energetic part of the spectrum according to SPOD.
6. Conclusions
The generation of mixing noise by isolated supersonic turbulent jets is known to be related to coherent structures (wavepackets) which can be successfully described as Kelvin–Helmholtz instabilities supported by the mean flow (Crighton & Gaster Reference Crighton and Gaster1976; Tam & Burton Reference Tam and Burton1984; Yen & Messersmith Reference Yen and Messersmith1999; Piot et al. Reference Piot, Casalis, Muller and Bailly2006; Rodríguez et al. Reference Rodríguez, Sinha, Brès and Colonius2013; Sinha et al. Reference Sinha, Rodríguez, Brès and Colonius2014). The performance of mean-flow-based linear stability calculations in describing coherent structures in supersonic twin jets, however, has not yet been fully characterised. Models are available based on simplified descriptions of the three-dimensional mean flow (see for instance Rodríguez Reference Rodríguez2021), which have been successful in computing wavepacket structures, but which lack experimental validation.
This work contributes to the modelling of twin-jet coherent structures, important for mixing noise, in two ways: by providing a higher fidelity, yet affordable description of the twin-jet dynamics using three-dimensional RANS calculations combined with plane-marching parabolised stability equations, and by providing the first quantitative validation of wavepacket models based on linear stability theory against experimental measurements of supersonic twin jets.
The analysis focused on supersonic round twin-jets generated by convergent–divergent nozzles operating at perfectly expanded conditions. Two different nozzle spacings
$s/D$
were considered, consisting of a closely spaced case (
$s/D = 1.76$
) featuring a strong interaction between jets, and a moderately spaced configuration (
$s/D = 3$
), characterised by a weaker interaction between the jets. Experimental measurements were performed to provide validation data. In particular, PIV measurements of the mean velocity field were made to compare against the RANS calculations, allowing quantitative comparisons in the symmetry plane containing the two jets (
$y/D = 0$
), and high-speed schlieren visualisations were performed to extract coherent structures present in the twin-jet system by means of SPOD, enabling qualitative and quantitative comparison against the PM-PSE fluctuations.
The three-dimensional RANS calculations were performed using a second-order finite volume solver featuring the Menter SST turbulence model. Modified turbulence model constants were employed, calibrated by linear interpolation between the original values provided by Menter (Reference Menter1994) and the optimised values provided by Ozawa & Nonomura (Reference Ozawa and Nonomura2024) for a single supersonic jet. The resulting mean-flow solutions are found to be in good agreement with the PIV measurements, demonstrating the ability of the RANS model to properly account for the nonlinear mean-flow interaction between the jets. The RANS results illustrate important mean-flow differences that emerge between the inner and outer mixing layers of each jet when the nozzle separation is small, such as for
$s/D = 1.76$
.
Plane-marching PSE wavepackets were computed for a range of frequencies relevant for mixing noise (
$St = [0.1,1]$
). The computed structures exhibit toroidal and flapping fluctuations following the natural symmetries of the twin-jet system, together with a non-axisymmetric Mach-wave radiation signature. Owing to the schlieren set-up, only those fluctuations symmetric with respect to
$z = 0$
can be visualised. The analysis was thus focused on toroidal and flapping fluctuations symmetric with respect to
$z = 0$
, namely oscillation modes SS0, SS1, SA0 and SA1. The PM-PSE calculations for closely spaced twin jets reveal significant differences in the phase of structures between the inner and the outer shear layers, especially for lower frequencies (
$St \leqslant 0.4$
), where interaction between the two jets is strong in terms of both nonlinear (mean-flow) and linear (wavepacket) dynamics. For symmetric modes, these differences are found to be induced by a mismatch in the phase speed of the instabilities which correlates with differences in the streamwise mean flow velocity. For antisymmetric modes, in addition to the mismatch induced by differences in mean-flow velocity, changes in the phase speed are also caused by a spatial constraint of the wavepacket structures in the shear layers that develop between the two jets. With some limitations inherent to the PSE formulation (Towne et al. Reference Towne, Rigas and Colonius2019), the modelled wavepackets for supersonic jets directly capture their Mach-wave radiation, which dominates the near-field structure.
A methodology for obtaining comparisons of the PM-PSE wavepacket models against experimentally educed structures was also presented, making use of a recently developed approach that facilitates the extraction of coherent structures from schlieren images (Prasad & Gaitonde Reference Prasad and Gaitonde2022; Padilla-Montero et al. Reference Padilla-Montero, Rodríguez, Jaunet and Jordan2024). The technique derives the line-of-sight integrated streamwise derivative of the momentum potential field from the schlieren images and extracts coherent structures of the momentum potential field instead of the original schlieren field, by means of SPOD. Comparison of SPOD modes and PM-PSE wavepackets reveals good qualitative agreement for symmetric and antisymmetric toroidal structures, as well as for symmetric flapping structures for the frequency range considered. Quantitative comparisons are performed through the definition of an alignment coefficient which measures the normalised projection of one field into the other. The alignment factors show agreement to be especially high in the range of frequencies at which the energy of the SPOD spectra is maximum (
$St = [0.3,0.7]$
), providing a first validation of the modelled wavepacket properties.
The results presented herein show that the modelling strategy based on RANS and PM-PSE can be useful for the physical understanding of how the interaction between the two jets affects the wavepacket properties, and, in particular, for parametric studies exploring different jet spacings and nozzle shapes. Agreement with the experimental measurements confirms mean-flow-based linear stability calculations can correctly capture coherent structures present in the twin-jet flow at conditions where non-modal effects are not important. More sophisticated approaches based on linear stability equations which enable the recovery of non-modal instabilities, like resolvent analysis (Garnaud et al. Reference Garnaud, Lesshafft, Schmid and Huerre2013; Jeun et al. Reference Jeun, Nichols and Jovanović2016; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018) or one-way stability equations (Towne & Colonius Reference Towne and Colonius2015), can offer improved modelling capabilities over the PM-PSE employed here, but their application to fully three-dimensional flows like the twin-jet configuration remains more computationally expensive despite recent developments (Martini et al. Reference Martini, Rodríguez and Towne2021; Towne et al. Reference Towne, Rigas, Kamal, Pickering and Colonius2022). In scenarios involving shape optimisation via iterative processes, where reducing the computational cost is of primary importance, the PSE approach continues to be of high interest.
New sets of experimental data remain necessary for the validation of structures antisymmetric with respect to
$z = 0$
, as well as for assessing the accuracy of the growth rates predicted by the linear theory. Finally, the methodology presented is only strictly valid for perfectly expanded conditions. Coherent structures developing in supersonic twin jets operating at over- or underexpanded conditions may be dominated by screech resonances, which pose additional modelling challenges for linear stability theory owing to the strong mean-flow gradients induced by shock and expansion waves, and the presence of upstream travelling waves (Edgington-Mitchell et al. Reference Edgington-Mitchell, Li, Liu, He, Wong, MacKenzie and Nogueira2022).
Acknowledgements
The authors thank S. Girard and D. Eysseric for the technical design and support with the experimental campaign, as well as for their insightful discussions.
Funding
The work of I.P.-M. has received funding from the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 101063992. The work of D.R. has been carried out under the project MIADDTRANS (PID2024-157642 MB-I00) funded by MCIN/AEI/10.13039/501100011033 and European Union’s FEDER. This work has also received funding from the Government of the Community of Madrid within the multi-annual agreement with Universidad Politécnica de Madrid through the Program of Excellence in Faculty (V-PRICIT line 3).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available upon request.
Appendix A. Comparison of mean flow solutions for different sets of constants of the Menter SST model
Figure 32. Comparison of mean streamwise velocity profiles along
$y$
(and at
$z = 0$
) for four different streamwise locations, between the RANS solutions with different sets of SST model constants and the PIV mean flow (
$s/D = 1.76$
):
$(a)$
$x/D = 2$
;
$(b)$
$x/D = 4$
;
$(c)$
$x/D = 6$
;
$(d)$
$x/D = 8$
.
Table 3. Different sets of values considered for the parameters of the Menter SST turbulence model. The same nomenclature as in Menter (Reference Menter1994) is followed.
Table 4. Details of the mean-flow grids considered for RANS calculations (
$s/D = 1.76$
).
Figure 33. Comparison of mean streamwise velocity profiles along the top nozzle axis between the RANS solutions with different sets of SST model constants and the PIV mean flow for
$s/D = 1.76$
.
Figures 32 and 33 illustrate the mean flow comparison between the RANS solutions computed using three different sets of constants for the Menter SST model, namely, the default SST constants implemented in TAU, the optimised SST constants from Ozawa & Nonomura (Reference Ozawa and Nonomura2024) and the average between both, all against the PIV mean flow. The corresponding constants for each case are listed in table 3. For the problem under study, the default set of constants is found to yield an overly diffusive solution (potential core is too short), while the use of Ozawa & Nonomura (Reference Ozawa and Nonomura2024)’s constants is found to yield a not enough diffusive solution (potential core is too long). For this reason, a set of constants corresponding to the average between both has been employed in this work, which is able to accurately reproduce the measured mean flow.
Appendix B. Effect of RANS grid resolution on PM-PSE wavepackets
Figure 34. Comparison of alignment factors for the five different mean-flow grids reported in table 4,
$s/D = 1.76$
:
$(a{,}b)$
alignment between schlieren fluctuations;
$(c{,}d)$
alignment between
$\hat {\varTheta }$
fields;
$(a{,}c)$
symmetric toroidal fluctuation;
$(b{,}d)$
antisymmetric toroidal fluctuation.
The RANS calculations have been performed on five different grids to study the impact of the mean-flow grid on the PM-PSE calculations and, in particular, on the alignment between the PM-PSE and SPOD structures. The five grid resolutions considered are summarised in table 4. The grids have been obtained by successively reducing the element size by a factor of
$\sqrt {2}$
between two consecutive grids.
Figure 34 shows a comparison of the alignment factors obtained for wavepackets SS0 and SA0 computed using the five mean-flow grids for selected frequencies (
$St = 0.4, 0.6$
and
$0.8$
). Both the symmetric and antisymmetric disturbances exhibit convergence with respect to the mean-flow grid. The sensitivity of PM-PSE wavepackets to the mean-flow grid resolution is higher for higher frequencies, which is a consequence of the smaller spatial wavelength associated with the developing structures. The PM-PSE results presented throughout this work have been computed using the RANS solution obtained on the finest mean-flow grid (G5).
Appendix C. Evolution of mean-flow quantities along the nozzle centreline
Figure 35. Evolution of the reciprocal of the normalised velocity (
$u_{\kern-1pt j}/\bar {u}$
) along the nozzle axis as a function of
$x$
:
$(a)$
$s/D = 1.76$
;
$(b)$
$s/D = 3$
. The single jet RANS solution is also added for comparison.
Figure 36. Evolution of the ratio of the square root of turbulent kinetic to the mean streamwise velocity along the nozzle axis as a function of
$x$
:
$(a)$
$s/D = 1.76$
;
$(b)$
$s/D = 3$
.
Figure 35 shows the evolution of
$u_{\kern-1pt j}/\bar {u}$
along the nozzle axis for the mean-flow solutions of the two studied twin-jet configurations and of the single jet (axisymmetric) case. The twin-jet curves are the direct reciprocal of the curves shown in figure 10. Far downstream of the end of the potential core (
$x/D \gt 20$
), the evolution of the axial velocity for the axisymmetric jet solution is linear, which indicates that the RANS calculations are consistent and able to asymptotically recover the self-similar region of jet development. The twin-jet solutions also exhibit a near-linear behaviour, with the case
$s/D = 3$
behaving much closer to the axisymmetric case.
For reference, the axial velocity decay constant
$B_u$
associated with each of the RANS solutions is calculated following the definition of Hussein, Capp & George (Reference Hussein, Capp and George1994) and Pope (Reference Pope2000), resulting in
$B_u = 6.8$
for the axisymmetric case,
$B_u = 10.6$
for the twin-jet with
$s/D = 1.76$
and
$B_u = 7.5$
for
$s/D = 3$
. The self-similarity of the velocity profile for a single round jet in the fully developed region implies that
$u_{\kern-1pt j}/\bar {u}$
must grow linearly with
$x$
, but not that the proportionality constant
$B_u$
should be universal. The value
$B_u = 5.8$
reported by Hussein et al. (Reference Hussein, Capp and George1994) and Pope (Reference Pope2000) corresponds to incompressible, unheated jets. Compressibility effects are expected to lead to different values of
$B_u$
. For comparison, the values of the decay constant for other numerical datasets of compressible single jets that were obtained from LES calculations are:
$B_u \approx 6.8$
for the openly available LES dataset from Brès et al. (Reference Brès, Jordan, Jaunet, Le Rallic, Cavalieri, Towne, Lele, Colonius and Schmidt2018) and Towne et al. (Reference Towne2023) (
$M_{\kern-1pt j} = 0.9$
, isothermal);
$B_u \approx 10.5$
for the B118 LES dataset from Sinha et al. (Reference Sinha, Rodríguez, Brès and Colonius2014) (
$M_{\kern-1pt j} = 1.5$
, isothermal) and
$B_u \approx 7.4$
for the B122 LES dataset from Sinha et al. (Reference Sinha, Rodríguez, Brès and Colonius2014) (
$M_{\kern-1pt j} = 1.5$
, temperature ratio 1.74). This range of values is aligned with the values obtained in the RANS calculations presented in this work.
Figure 36 also shows the evolution of the ratio of the square root of the turbulent kinetic energy to the mean streamwise velocity along the nozzle centreline for the twin-jet RANS solutions and the PIV measurements. For the experimental data,
$k$
is computed from the root mean square of the velocity fluctuations in the PIV plane. A good agreement is obtained between the RANS mean-flow and the PIV data within the available experimental window.
Open access
