Swirling Flow Effects on the Aeroacoustic Signature of an Aerospike Nozzle

High-speed aircraft or rocket propulsion systems use nozzles to generate thrust. Those nozzles cannot always be operated at design conditions; the supersonic jet is not expanded ideally, leading to very intense acoustic noise. In particular, a shock cell structure may arise due to a deviation from the ideally expanded jet conditions at the inlet of the supersonic nozzle. Generally, the aeroacoustic signature of such device is governed by the boundary conditions at their inlet. For instance, an increase in jet temperature can enhance the Mach wave and crackle noise radiation. Inhomogeneities in the velocity distribution might also lead to increased sound pressure levels (SPLs) in the far field. In aerospace applications, mechanical work is usually gained by implementing a turbine between the combustion chamber and the nozzle. While such turbines have been largely used in the aerospace industry for subsonic applications for many decades, concepts of supersonic turbines have emerged in the past few years [1–3]. In those configurations, the high-speed flow at the outlet of the turbine can feature a swirl component. This in turn might affect the aeroacoustic signature of the nozzle and of the whole engine.

Several noise components are present in supersonic jets. They are linked to the development of flow instabilities and the emerging shock cell structure in nonideally expanded cases. Turbulent mixing noise is present in both subsonic and supersonic jets. It is generated at the end of the potential core by large turbulent structures emitting noise at low frequencies [4]. The characteristic Strouhal number might vary slightly, but is usually around $St=fD/uj∼0.30$⁠, where $D$ is the jet diameter and $uj$ is the velocity of the ideally expanded jet [5].

Screech noise is a sound component generated by the interaction between vortical structures convected downstream and the quasi-periodic shock cell structure. It propagates in the upstream direction in the far field. It was first detected by Powell [6]. The generation mechanism is described as follows: in the first step, turbulent structures are convected downstream in the shear layer and interact with the quasi-periodic shock cell structure of the jet. The passage of vortices across the shock generates acoustic waves propagating in the upstream direction. Then, the upstream propagating waves reach the nozzle trailing edge and impinge on the latter, exciting the jet shear layer. This closes the resonant loop. However, two distinct mechanisms for the generation of such upstream propagating waves exist [7,8]. In the first mechanism, the interaction between the quasi-periodic shock cell structure and the turbulent structures convected downstream generates waves propagating in all directions, some of them propagating in the upstream direction. In the second mechanism, the resonant loop is closed by neutral modes of the jet. In this scenario, the supersonic jet exhausting the nozzle is acting as a waveguide and the feedback acoustic mode closing the resonant loop is guided upstream along the jet [9]. These two mechanisms can coexist in a supersonic jet [7].

Broadband shock-associated noise (BBSAN) is linked to the passage of vortical structures across the quasi-periodic shock cells of the nonideally expanded jet. The passage through a shock leads to vortex stretching and the emission of an acoustic wave. It has first been identified by Martlew [10]. Harper-Bourne and Fisher derived a model allowing to compute the central frequency of this noise component as a function of the observation angle [11]. It is in particular linked to the shock cell length and the convection velocity of the vortical structures in the shear layer.

Mach waves are generated when the convection velocity of the flow vortical structures in the jet shear layer is above the speed of sound. The Mach wave directivity angle is given by the model of Oertel [12] and is a function of the convective Mach number. The generation mechanism can be explained by the wavy wall analogy; the wavy wall consists of the turbulent structures in the jet [13]. It travels downstream at supersonic speed and generates compression waves propagating in the far field. For highly heated jets, crackle noise along the Mach wave propagation direction can also arise. It is generated in the jet and can be amplified by nonlinear wave steepening effects [14]. Crackle is observed as sudden bursts in the pressure signals [15].

Swirling conditions have been studied for jets, both numerically and experimentally [16–19]. Stability analyses of compressible, swirling mixing layers were also performed [20–22]. It was found that an additional swirling component significantly enhances the maximum amplification rate for compressible mixing layers and alters the stability of helical and axisymmetric jet modes [20]. Furthermore, it was found that vortex breakdown plays an major role in the dynamics of a swirling jet; vortex bursting produces high turbulence fluctuations spreading radially from the vortex core to the jet shear layer. The interaction between the shear layer and the high turbulence fluctuations leads to a faster decay of the coherent structures [23]. Subsequently, the spreading rate, the mass entrainment rate, and the growth rates of the shear layer momentum thickness are higher for a swirling jet than a nonswirling one [24]. The effect of the swirling flow on the aeroacoustics of supersonic jets has been looked at, however, from a theoretical or experimental point of view [17,25]. It was shown that swirling could reduce or eliminate BBSAN, leading to a reduction of the total acoustic power radiated [25]. Experimental studies suggested that an additional swirling motion would be responsible for the disappearance of screech tones [17]. In particular, the additional swirl component leads to a shortened shock cell structure, and the mixing rate is increased compared to nonswirling cases. Large Eddy simulations for a subsonic swirling jet at $M=0.9$ were performed [16]. It was found that the instability amplification rates in swirling jets are higher than that in the nonswirling case. Moreover, the jet spreading as well as the vorticity thickness are enhanced. However, the effect of swirling on the flow characteristics using numerical simulations and the consequent aeroacoustic signature for annular, supersonic jets has not been investigated. This study aims to fill this gap by examining how an additional swirling component influences the attributes of supersonic jets and the resulting aeroacoustic signature.

In Sec. 2, the geometry, the method, and the numerical setup are presented. The time-averaged flow characteristics as well as the swirl number and the jet spreading effects are presented in Sec. 3. Subsequently, the aeroacoustic results are discussed in Sec. 4.

Compressible implicit large eddy simulation (LES) of an annular supersonic jet exhausting a three-dimensional, axisymmetric, truncated aerospike nozzle configuration are performed. The simulations are carried out at a constant nozzle pressure ratio of 3 and for a temperature ratio of 1, without swirling effects, referred to as the baseline case, as well as for three different Swirl numbers $S$⁠, i.e., 0.10, 0.20, and 0.30, respectively. A grid sensitivity study for the cold aerospike nozzle jet without swirling was performed [26]. It was used to identify the necessary grid resolution able to represent the jet characteristics relevant for aeroacoustic computations. In the present investigation, the finest grid of 170 million cells was used.

The equivalent diameter for the annular section is $Deq=84.7$ mm. It aligns with the diameter that the annular chamber would have for an equivalent circular cross-sectional exit area. In the present study, all geometric parameters are normalized using this equivalent diameter. The aerospike nozzle has a total length of $2.08Deq$⁠. The annular channel has a width $0.16Deq$⁠. The nozzle lip has a width $0.14Deq$⁠. Throughout this article, velocities and Mach numbers will be scaled using the values of the jet velocity $uj=399$ m/s and Mach number $Mj=1.36$⁠, respectively, corresponding to the ideally expanded cold jet parameters without swirling. The dimensionless time-steps used for each simulation are given in Table 1. The introduction of an extra swirling component at the nozzle inlet plane results in smaller time-steps, thereby increasing the computational costs.

Large eddy simulations are performed using an in-house finite volume-based compressible flow solver. An explicit four-stage Runge–Kutta algorithm is applied for time integration and a second-order central difference scheme is used for spatial discretization. A modified artificial dissipation scheme after Jameson has been used for this study. It is applied at the end of every first stage of the Runge–Kutta method. It was validated for three acoustic benchmark cases: an acoustic pulse in a uniform flow, a shock propagation, and a vortex–shock interaction. The modified artificial dissipation scheme aims to capture shock waves and avoid Gibbs-like oscillations near shocks, as well as to relax subgrid-scale turbulent energy [28]. Detailed explanation about the scheme and its application in the LES context can be found in Ref. [29]. Further details about the artificial dissipation modification for the solver used in this study can be found in Ref. [30]. An implicit LES approach (also referred to as monotonically integrated large eddy simulation) is adopted. In other words, no explicit subgrid-scale model is used for the simulations. The physical dissipation produced by small-scale eddies is implicitly considered by taking advantage of truncation errors [31]. Viscosity is estimated using Sutherland’s law, and a fourth-order central difference is used to calculate viscous flux term. At $T=293K$⁠, the viscosity is $ν=1.789×10−5Pas$⁠. The validation of the aforementioned method using this solver has been successfully carried out across diverse nozzle configurations and operating conditions, extending up to jet temperature ratios of 4 [30,32–38]. The computational fluid dynamics results for both the flow and acoustic near and far fields are compared with particle image velocimetry data and aeroacoustic measurements for the various nozzle geometries considered. Very good levels of agreement between the experimental and simulation data are reached.

Figure 2 features the axisymmetric computational grid on the $yz$-section planes with relevant domain dimensions. The flow direction is the $z$-direction. The plane $z/Deq=0$ corresponds to the outlet of the annular section. The three-dimensional, axisymmetric domain has a diameter of $28Deq$ at this plane. The computational domain is of total length $65Deq$ and extends from $z/Deq=−8$ to $z/Deq=57$⁠. At the far-field outlet (⁠$z/Deq=57$⁠), the domain has a diameter of $44Deq$⁠. The domain size has been chosen in accordance with previous studies of supersonic jets with various nozzle geometries using the same solver [30,37,38]. Further details about the mesh are given in the next paragraph.

A total pressure of $pt=304,000Pa$ and a total temperature of $Tt=293K$ are applied as inlet boundary conditions for all the cases. The turbulent intensity is set to zero at the inlet plane, which corresponds to a laminar state. The velocity vector is prescribed at the inlet plane (⁠$z=−1.65Deq$⁠) to achieve the desired Swirl number $S$⁠. The static pressure, the static temperature, and the velocity are set to $p=101,325Pa$⁠, $T=293K$⁠, and $u=0m/s$ at the outlet and radial far-field boundaries, respectively. A secondary flow with a Mach number of $0.05$ is applied at the coflow inlet located at $z/Deq=−8$ and at a radial distance $r>1.06Deq$ (see Fig. 2). Adiabatic no-slip boundary conditions are applied in the annular chamber and at the aerospike body in a weak sense. Characteristic boundary conditions are used to avoid reflections [39] at the far-field boundaries. In addition, a mesh stretching rate of 5% is applied outside of the flow region (shown as a thin white line on Fig. 2) to further dissipate the acoustic waves generated in the near field and propagating in the far field. This makes the far-field boundaries reflection free. Lower stretching rates would not sufficiently dissipate the acoustic waves propagating toward the outer domain boundaries, and larger stretching rates would lead to reflections from the grid [40,41]. The same grid stretching strategy was applied in previous studies [26,30,32–38,42]. This meshing strategy leads to a structured grid of 170 million cells composed of 119 blocks. The mesh is designed to fulfill $y+∼1$ in the annular duct. In the region where the shear layer is formed (around $r∼Deq=0.9Deq$⁠), the smallest cell has a size of about $Δz2+Δr2∼Deq/400$⁠. Further downstream, the cells have a size of $Δz=Deq/150$⁠, $Δz=Deq/13$⁠, $Δz=Deq/8$ at $z=2Deq$⁠, $z=20Deq$⁠, and $z=40Deq$⁠. The average cell size in the radial direction is kept around $Δr∼Deq/105$ in the potential core. This grid resolution in the potential core guarantees the radial propagation of acoustic waves at high frequency with a point per wavelength number (⁠$PPW)=16$⁠.

The Ffowcs Williams–Hawkings (FWH) equation is solved to compute the aeroacoustic signature in the far field. It requires data acquisition on a permeable surface enclosing the relevant acoustic sources, that is, the shear layers and the shock cell structure for a supersonic jet. The surface in the downstream at $z/Deq=40$ is open to avoid pseudo-sound due to hydrodynamic fluctuations overlapping the cross section [43]. The corresponding domain enclosing the flow region has a total length of $42Deq$ and a diameter of $7.5Deq$ at the annular duct outlet (⁠$z/Deq=0$⁠) and $15Deq$ at $z/Deq=40$⁠, location at which the surface is open. The choice of these dimensions is based on previous simulations using the same solver. The mesh stretching is kept below 1% to guarantee an undisturbed wave propagation in that region. The permeable surface enclosing the relevant acoustic sources is depicted as a thin white line in Fig. 2 with the dimensions previously mentioned. Henceforth, this surface will be denoted as the FWH surface. The numerical simulations are performed for a total physical time of 65 ms with a transient of 15 ms. Data for the far-field acoustic analysis are exported on 448,400 nodes for a physical time of 50 ms every $2×10−6s$⁠. The four simulations required around 570,000 CPU hours each on average to get signals of the sufficient length.

Figure 3 features the time-averaged Mach numbers normalized with $Mj$⁠. It is averaged on the $xz$- and $yz$-planes.

In the baseline case without swirling, the aerospike nozzle jet is composed of two distinct parts: an annular and a circular part [26], located between $z/Deq=0$ and $z/Deq∼2.1$ and downstream of $z/Deq∼2.1$⁠, respectively. The jet in the baseline case is depicted in Fig. 3(a). For the cases with swirling, an annular shock cell structure is formed as well. The supersonic flow is nonattached in the direct vicinity of annular nozzle. The flow reattaches at a distance $z∼1.20Deq$ further downstream compared to the case without swirling with a flow reattachment around $z∼Deq$ (see Fig. 3(b)). A separation bubble featuring negative velocity is formed around $(r/Deq,z/Deq)∼0.65$⁠. Due to the higher shear in the inner annular mixing layer, the separation bubble is lengthened with the increasing Swirl number. The three first annular shock cells belong to the nonattached part of the jet in the swirling case, whereas only two belong to the nonattached part of the jet in the baseline case. The shock cells downstream (from $z=1.20Deq$⁠) that are reattached to the aerospike bluff body exhibit a much lower Mach number jump across them, making them less significant in terms of noise generation. The Mach number jumps across the shock cells are reported in Table 2 with an uncertainty of $±0.003$ due to the interpolation and the choice of monitoring lines. The observed shocks in the annular shock cell structure for the swirling cases are stronger with respect to the baseline case without swirling. The Mach number jumps across the shocks decrease marginally with the increasing swirl number. Additionally, the shock cell length reported in Table 3 with an uncertainty of $±0.01$ is also affected by the swirling motion. The first shock cell is longer compared to the case without swirling motion. The lengths of the second and third shock cells are comparable. Additionally, the three shock cells maintain consistent lengths across all swirl numbers. However, in the baseline case, the third shock cell is elongated, attributed to the flow reattachment on the aerospike bluff body occurring further upstream in comparison to the swirling cases.

A shock cell structure downstream of the aerospike bluff body is observed in the baseline case without swirling and in the case $S=0.10$⁠. The shock cell structure disappears for higher Swirl numbers. The introduction of swirling motion at the nozzle inlet plane intensifies shear in the azimuthal direction. Hence, the pressure adaption takes place more quickly compared to the cases without swirling, resulting in a decreased number of shock cells. In the baseline case, annular shock cells downstream of the aerospike bluff body coalesce into a circular shock cell structure. The jet regains a fully circular character in the baseline case at approximately $z∼4.5Deq$ at the location of the third shock cell. For the case $S=0.10$⁠, an annular shock cell structure with fewer shock cells is formed. In that case, the annular shock cells do not merge. The radial velocity gradient is observed not only at the outer shear layer around $r∼0.5Deq$ but also in the inner shear layer around $r∼0.1Deq$⁠. This dual presence contributes to increased dissipation and a subsequent reduction of the shock cell count downstream of the aerospike bluff body. Ten circular shock cells are observed in the baseline case without swirling on Fig. 3(a), while only four annular shock cells are observed for the cold case at a Swirl number $S=0.10$ up to a distance $z∼5Deq$ in Fig. 3(b) in the jet part located downstream of the aerospike bluff body (from $z∼2.1Deq$⁠). The Mach number jumps and the shock cell length in that part of the jet are reported in Tables 4 and 5. These were computed with an uncertainty of $±0.003$ for the Mach number jumps and $±0.01$ for the shock cell length due to the interpolation and the choice of monitoring lines. The correct Mach number jump across the fourth annular shock cell in the swirling jet with $S=0.10$ was too low to be correctly computed.

The Mach number jump across the first shock cell downstream of the aerospike bluff body is slightly lower for the case $S=0.10$ compared to the baseline case without swirling. The Mach number difference across the following shock cells decreases rapidly for the case with $S=0.10$⁠. Due to the increased shear, the swirling motion leads to a quicker pressure adaption, and as a consequence, the shocks are less strong. The Mach number difference across the first three shock cells in the circular part of the jet is $ΔM¯/Mj=0.907$ for the baseline case and $ΔM¯/Mj=0.058$ for the case with $S=0.10$⁠. Moreover, the shock cell length also decays faster for the swirling case compared to the baseline case without swirling motion as can be seen in Table 5. This leads to a shorter shock cell structure. On average, the first three shock cells display a length of $L¯/Deq=0.68$ for the baseline case without swirling and $L¯/Deq=0.57$ in the case $S=0.10$⁠. This decrease in shock spacing is consistent with the previous experimental observations [17]. The baseline case without swirling exhibits a shock cell length decay of 3%, which matches with André’s results obtained with a secondary flow of $M=0.05$ [44]. Shock strength affects sound generation; the passage of vortical structures through weaker shocks leads to reduced SPLs. Additionally, a decreased count of shock cells within the potential core leads to less interaction between the shock cell structure and the aforementioned vortical structures. Shorter shock cells lead to a shift of the central frequency of the BBSAN toward higher frequencies. The effect of the flow characteristics on the aeroacoustic signature will be characterized further in Sec. 4.

At Swirl Number $S=0.20$⁠, one annular shock cell of length $L/Deq=0.67$ is still noticeable. However, the Mach number jump across it is much lower compared to the baseline case and the case with $S=0.10$ (around $ΔM/Mj=0.020$⁠). The radial velocity gradient observed in the inner shear layer in the case with $S=0.10$ is more pronounced for in the case with $S=0.20$⁠, with a slight negative axial velocity. Hence, the annular character of the jet downstream of the aerospike bluff body is more pronounced with the increasing Swirl number. At Swirl number $S=0.30$⁠, the jet undergoes an expansion at the vicinity of the aerospike body tip and the shock cell structure disappears completely. In that case, the supersonic flow undergoes an acceleration toward the tip of the aerospike bluff body at $z=2.08Deq$⁠. Due to the cross-sectional change, the supersonic flow is further accelerated downstream of $z∼2.1Deq$⁠. The observed trend persists as swirl numbers further increase in subsequent simulations. For brevity, detailed results at these higher swirl numbers are not presented in this article. The flow acceleration along the aerospike nozzle is located further upstream for $S=0.30$ compared to $S=0.20$⁠.

The pressure fluctuations for the baseline case and the swirling case with $S=0.10$ are shown in Fig. 4. High levels of pressure fluctuations in jets are linked to increased sound pressure levels [38]. The baseline case without swirling features high levels of pressure fluctuations in the annular part of the jet (see Fig. 4(a)). In particular, high levels are found at the location of the junction between the first and the second shock cells. Lower fluctuation levels are found in the outer annular shear layer and at the location of the separation bubble, as well as at the location of the expansion fan downstream of the aerospike bluff body (at $z∼2.3Deq$⁠). High fluctuation levels are found at the end of the circular shock cell structure, around $z∼7Deq$⁠. This is linked to the flapping motion of the circular shock cell structure associated with screech. In the case with $S=0.10$⁠, high fluctuation levels are found further downstream. This indicates a change of dynamics in the motion of the annular shock cell structure. However, higher pressure fluctuation levels are found downstream of the aerospike bluff body. This is due to the formation of a swirling mixing layer between the attached annular jet and the air trapped in this central region. Similar pressure fluctuation fields were observed at higher swirl numbers but are not shown here for brevity. The central region around $r/Deq<0.25$ where the pressure fluctuation levels increases with the increasing swirl number.

At the outlet of the annular chamber, the high shear due to the additional azimuthal component leads to the formation of an inner and outer swirling mixing layer. The flow reattaches on the aerospike further downstream at around $z/Deq∼1.20$⁠. The inner annular mixing layer disappears. Downstream of the aerospike bluff body, a second expansion takes place and an annular mixing layer with azimuthal shear is formed. For the lowest Swirl numbers $S=0.10$ and $S=0.20$⁠, the recirculation region remains confined close to the tip of the aerospike bluff body. The axial velocity decreases, but no backflow is observed. In the case with $S=0.30$⁠, a strong adverse pressure gradient is generated in the central region located at a radius $r/Deq<0.15$⁠, leading to a backflow towards the aerospike bluff body. The recirculation zone in that region features stronger velocity gradients with increasing Swirl number. This region is characterized by wake-like conditions and might feature vortex breakdown with further increasing swirl number.

Figure 5(a) shows that the local swirl number is maximum at the aerospike body tip, before decreasing further downstream. For high swirl numbers, the annular character of the jet is preserved, with an outer and an inner swirling mixing layer with strong velocity and pressure fluctuations. The absence of shock cell structure in this region eliminates the generation of BBSAN and screech resulting from the interaction between vortical structures and the quasi-periodic shock cell structure. Hence, although the pressure fluctuations intensity increases with the increasing swirl numbers, those pressure fluctuations are linked to a free shear flow. In supersonic jets, this represents a secondary source of noise compared to phenomena like screech, Mach waves, or BBSAN. Moreover, the corresponding vortical structures have less coherence and decay faster in time and space, making them less-efficient quadrupole sources. As a consequence, lower SPL are expected for the increasing swirl numbers. Further details about acoustic generation mechanism is discussed in Sec. 4.

Since the simulation case is axisymmetric, the time-averaged flow data remains invariant under rotation. The integration along the azimuthal direction yields a factor $2π$ canceling out, which leaves us solely with the integration over the radius $r$⁠. $ρ$⁠, $uθ$⁠, and $uz$ correspond to the local density, the local azimuthal velocity, and the local velocity in the $z$-direction, respectively. These are a function of the radius $r$ and the axial distance $z$⁠. The integration is carried out over the radial domain $Ωr=[R(z);r0.05(z)]$⁠, where $R(z)$ corresponds to the radius of the aerospike body depicted in Fig. 1 and $r0.05(z)$ is defined as the radius in a cross-sectional $Σ$ located at the axial distance $z$⁠, where $M(r0.05(z),z)=0.05maxΣM(r,z)$⁠, with M being the local Mach number. In other words, $r0.05(z)$ depicts the outer shear layer of the annular jet and is a measure of the jet spreading. It is shown in Fig. 5. $r(z)¯$ is the average radius of the radial integration domain $r(z)¯=[R(z)+r0.05(z)]/2$⁠. For the circular part of the jet, the integration approach is the same as shown in Eq. (2); however, the integration domain is now $Ωr=[0;r0.05(z)]$ and the average radius is $r(z)¯=r0.05(z)/2$ due to the aerospike bluff body truncation. This change of average radius leads to a discontinuity for $s(z)$ at the interface between the annular and circular regions located at $z/Deq=2.08,$ which is noticeable in Fig. 5(a).

The local swirl number for the jet downstream of the aerospike bluff body is computed along the flow direction between $z/Deq=2.08$ and $z/Deq=8$⁠. This corresponds to the length of the potential core for the baseline case without swirling (see Fig. 3(a)). The results for all jets are presented in Fig. 5. The swirl number as defined in Eq. (2) as function of $z$ and normalized with the imposed Swirl number $S$ at the nozzle inlet plane is shown in Fig. 5(a). The separation between the annular and circular regions is marked by a black thin, dashed line. The local swirl number $s(z)$ is displaying a discontinuity at this interface, as previously mentioned, for the cases with swirling. The baseline case without swirling features a normalized swirl number equal to 1 for all $z$⁠, which is consistent with the expectations. The local swirl number on a cross section increases in the annular region of the jet, following the same pattern for all nonzero swirl numbers. The flow acceleration observed at the tip of the aerospike bluff body shown in Figs. 3(c) and 3(d) is hence mostly due to an enhancement of the swirling component. Downstream of the aerospike bluff body, the swirl spreads out due to viscosity effects. The normalized local swirl numbers fall approximately on the same line up to an axial distance $z∼5Deq$ for the three Swirl numbers, suggesting that the swirling behavior follows a self-similarity law in that region. Further downstream, only the lines for $S=0.10$ and $S=0.20$ fall on the same line. This might be due to the shortening of the potential core in the case $S=0.30$ as observed in Fig. 3(d). Figure 5(b) features the parameter $r0.05$⁠, which serves as a measure to characterize jet spreading. The jet almost does not spread in the annular part for the baseline case and across the three Swirl numbers. In this region, $r0.05$ is highest for $S=0.30$⁠. Downstream of the aerospike bluff body, the jet radially spreads more with the increasing Swirl number. At $z=8Deq$⁠, the jet is almost wider by $Deq$ for $S=0.30$ compared to the baseline case without swirling.

The increase in Swirl number $S$ affects the shock cell structure, the jet spreading, and the flow entrainment. As a consequence, the thrust performance is affected by the Swirl number change. The strong adverse pressure gradient and the resulting backflow in the case of $S=0.30$ might lead to a significant decrease in thrust compared to the other cases.

The two lowest swirl numbers feature a slight thrust decrease compared to the baseline. The case with $S=0.30$ features a significant drop in performance compared to the baseline case. This is due to the strong adverse pressure gradient in the region $r<0.20Deq$ with a more pronounced backflow and featuring wake-like conditions. The flow entrainment is enhanced with the increasing Swirl numbers. However, at lower Swirl numbers, the absence of the strong backflow allows for the sustained generation of high thrust.

Screech is a tonal noise component linked to a feedback mechanism characterized by the presence of upstream propagating waves in the jet shear layer. The upstream propagating waves can correspond to free propagating acoustic waves or neutral modes of the jet, the latter acting as a waveguide. The underlying physical mechanisms have been described in Refs. [13,45]. The screech tones propagate in the far field at upstream angles (⁠$θ=20$ deg). Such upstream propagating waves are detected by computing two-point space-time cross correlations in the jet shear layer. In this study, pressure data were recorded on four lines composed of 84 and 96 monitoring points equally spaced in the circular jet shear layer (⁠$r/Deq=0.40,0.55$⁠). The reference point for the cross correlation is fixed at $z/Deq=2.05$ at the junction between the circular and annular shear layers. The corresponding lines in the circular shear layer are lines $L2$ and $L3$ shown in Fig. 6, consisting of 84 and 96 monitoring points, respectively.

The flapping motion in radial direction of both the shock cell structure and the separation bubble trapped between the shock cell structure and the aerospike bluff body leads to sound radiation in the far field (as reported in Refs. [26,42]). Monitoring points are placed in the annular shock cell structure to identify its oscillations modes. They are considered at the location of the separation bubble trapped between the shock cell structure and the aerospike bluff body (point $P1$⁠), at the junction between the expansion cell and the first shock cell (point $P2$⁠), and at the junction between the first and the second shock cell (point $P3$⁠). Figure 6 features the three monitoring points. The exact location of the points varies slightly between the baseline case and the swirling case, since the length of the annular shock cells and the location of the flow reattachment on the aerospike bluff body differ. Velocity data at the these monitoring points are acquired for $50$ ms. The power spectral density (PSD) of the radial velocity is shown for the four cases in Fig. 8. The first, second, and third columns show the PSD results for points $P1$⁠, $P2$⁠, and $P3$⁠, respectively. The most important peaks are marked with a black, dashed line.

A more thorough analysis of the radial oscillations modes for the baseline case without swirling can be found in Ref. [26]. The main results are summarized in Table 7. The spectra for the baseline case without swirling are shown on the first row in Fig. 8. The results for the case with $S=0.10$⁠, $S=0.20$⁠, and $S=0.30$ are shown in the second, third, and fourth rows of Fig. 8, respectively. In the baseline case, the separation bubble (point $P1$⁠) features two main oscillations modes with a broadband character around $St∼0.52$ and $St∼1.25$⁠. The latter is divided into two subpeaks at $St∼1.21$ and $St∼1.36$⁠. Furthermore, secondary peaks are observed, among others a peak at $St=0.68$⁠. The first shock cell oscillates at a Strouhal number $St∼2.60$ as shown on the spectrum for point $P2$⁠. The lower oscillation levels and the higher radial frequency suggest a stiff motion of the first shock cell. The second shock cell (point $P3$⁠) consistently oscillates at the same frequency as the separation bubble $St∼1.21,$ which exhibits a broadband character. A second peak at $St=0.68$ is observed, which couples with one of the secondary peak of the separation bubble.

The same Strouhal number of $St=0.85$ is found for the radial oscillation of the separation bubble (point $P1$⁠) for all swirling cases. This suggests that the dynamic behavior of the separation bubble trapped between the shock cell structure and the aerospike bluff body is governed by the reattachment characteristics rather than the swirl number itself. The flow reattaches at the same axial distance for all swirling cases at $z∼1.20Deq$⁠, while it reattaches further upstream at a distance $z∼Deq$ for the baseline case without swirl. This Strouhal number is consistently observed for point $P3$ but at a lower amplitude compared to point $P2$⁠. A subharmonic of the main oscillation mode at $St=1.71$ is observed in the spectra for all the swirling cases. Further subharmonics of this main mode are observed in the spectra. These are especially visible for $P2$ in the case with $S=0.10$⁠. They are found around $St∼2.55$ and $St∼3.35$ at lower amplitudes compared to the main mode and the first subharmonic. In the case with $S=0.30$⁠, the azimuthal mode identified in the two-point cross correlation in Fig. 9 exhibits a Strouhal number $St=0.90$⁠. A coupling phenomenon between the azimuthal mode and the radial oscillation mode of the annular shock cell structure takes place since a peak at $St=0.91$ is observed only in the case with $S=0.30$⁠. The amplitude for the radial oscillation mode at $St=0.85$ in the case $S=0.30$ is lower compared to the cases with lower Swirl numbers $S=0.10$ and $S=0.20$⁠. The observed peaks at $St=0.85$ in the spectra of the two latter cases correspond to an overlap of the identified azimuthal mode (see Table 8) and the radial oscillation mode of the annular shock cell structure.

Moreover, a secondary peak is observed at a lower frequency around $St∼0.5$ for all the cases. The Strouhal number of this component slightly decreases with the increasing Swirl Number: $St=0.57$ for $S=0.10$⁠, $St=0.53$ for $S=0.20,$ and $St=0.48$ for $S=0.30$⁠. Further analyses are required to determine the physical origin of this frequency shift. In the baseline case without swirling, the peak observed at $St=0.52$ is linked to a coupling with the screech tones while it cannot be the case for the swirling cases due to the absence of screech. Further details can be found in Sec. 4.3. The peak at around $St∼0.5$ is observed at a lower amplitude in the PSD spectra for point $P2$⁠. Generally, the energy levels are lower for point $P2$ for all cases, indicating a lower level of fluctuations compared to the separation bubble and further downstream. Furthermore, the oscillation modes of the shock cell structure are affected by the swirling. Secondary peaks at Strouhal numbers $St=1.35$ and $St=1.40$ are observed for cases $S=0.10$⁠, $S=0.20$⁠, and $S=0.30$⁠, respectively, in the PSD spectra for point $P2$⁠.

The swirl imposed as a boundary condition at the nozzle inlet plane leads to an enhancement of azimuthal modes in the jet. Two-point space-time cross correlations of the pressure and the velocity between several lines located in the annular shear layer are performed to identify azimuthal modes and compute the convection velocity of the vortical structures. The relevant monitoring lines in the $yz$-plane used are shown in Fig. 6. Fifty lines each consisting of 25 points are placed in the annular shear layer. They cover an angular domain of $π$⁠. The first set of 25 lines corresponds to the rotation of line $L0$ with an angular increment of $7.5deg=π/24$ in the counterclockwise direction. The second set of 25 lines corresponds to the rotation of line $L1$ with an angular increment of $7.5deg=π/24$ in the counterclockwise direction. They cover the annular flow domain between $z/Deq=[0.18,2.08]$⁠. The method used to identify the azimuthal modes of the annular jet will be described in the following.

Let us call the lines covering the angular domain of $π$$L0nπ/24(z)$ and $L1nπ/24(z)$⁠, with $n$ varying from 0 to 24, and $L0$ and $L1$ corresponding to the lines located in the $yz$-plane (see Fig. 6). Let us fix a position $z0$ along the flow direction. The two-point space-time cross correlations are computed between the time series at point $Lijπ/24(z0)$ and time series along the line $Li(j+1)π/24(z)$ as well as between point $Lijπ/24(z0)$ and line $Li(j+2)π/24(z)$ with $i={0,1}$⁠. This allows to follow the convection of a vortical structures in the azimuthal direction. A space-time cross correlation matrix of the size of the time-lag and the vector $z$ of considered locations in the flow direction is obtained. The maximum cross correlation level and its location $z1$ are identified, indicating the most probable movement pattern of vortical structures in the annular shear layer in the azimuthal direction. This location $z1$ is used as an input for the computation of the cross correlation between the time series at point $Li(j+1)π/24(z1)$ and line $Li(j+2)π/24(z)$ or point $Li(j+2)π/24(z1)$ and line $Li(j+4)π/24(z)$⁠. A new space-time cross correlation matrix of the size of the time-lag and the vector $z$ of considered locations along the flow direction is obtained. The maximum cross correlation level and its location $z2$ are identified. The latter is used for the next cross correlation computation. This method is applied again for the following lines. This procedure allows to follow the movement pattern of vortical structures over the whole azimuthal range and not only between two individual lines. The resulting space-time two-point cross correlation over all the lines is obtained by summing up all the individual cross correlations between two single lines. Moreover, to ensure the absence of any other tonal components, the cross correlations underwent temporal filtering through multiple bandpass filters. This aimed to reveal potential periodic patterns not readily apparent behind more dominant periodic patterns. Only a single periodic pattern for the various swirling cases was detected. The results are shown in Fig. 9 for cross correlations between lines $L0jπ/24$ and $L0(j+2)π/24$⁠, with $j∈{0,…,22}$ for $S=0.30$⁠. The $x$-axis corresponds to the axial distance along the flow direction and the $y$-axis corresponds to the time-lag in seconds. Additional cross correlations following the same procedure described previously were performed between lines $L0jπ/24$ and $L0(j+1)π/24$ as well as $L0jπ/24$ and $L0(j+3)π/24$⁠, as well as for the set of azimuthal lines $L1$ at a larger radius (see Fig. 6). The results for these additional cross correlations are very similar to the results shown in Fig. 9 and are not shown here for brevity. The same procedure has been applied to the cases $S=0.10$ and $S=0.20$⁠. Almost identical cross correlations were obtained; only the time-lag between the maxima was larger.

Figure 9 displays periodic patterns with negative slopes. This indicates that vortical structures are convected in the counterclockwise azimuthal direction. The obtained time-lag for a periodic pattern is $Δτ=2.35×10−4s$ with an uncertainty of $±5×10−6s$⁠, which yields a Strouhal number $St=0.90±0.01$⁠. As mentioned, similar cross correlations are computed for the remaining cases. The results resemble the cross correlation shown in Fig. 9 with similar periodic patterns, but a slightly larger time-lag. The corresponding cross correlations for the other cases with swirling are not shown for brevity. The Strouhal numbers based on the time-lag observed between two periodic patterns on the two-point cross correlations are summarized in Table 8 for all four cases.

In the baseline case without swirling, the azimuthal mode at $St∼0.53$ is due to the coupling with the screech component generated downstream in the circular shock cell structure. Moreover, the slopes were almost horizontal suggesting a weak convection of vortical structures in the azimuthal direction [26]. This coupling between an azimuthal mode and the screech mode was also observed at higher temperature ratios for nonswirling jets [42]. The Strouhal number of the azimuthal mode slightly increases with the increasing swirl numbers but are in the same Strouhal number range (⁠$St∼0.85$⁠). The enhancement of azimuthal modes leads to sound propagation with a helical shape at upstream angles in the far field. Moreover, the PSD of the radial velocity presented previously (see Fig. 8) suggests a coupling phenomenon between the azimuthal mode and the radial motion of the annular shock cell structure.

The computation of the central frequency of the BBSAN requires the knowledge of the convection velocity (see Sec. 4.5). The axial convection velocity of the vortical structures was determined through cross correlation of the velocity along the flow direction. This holds particular significance in establishing the central frequency of the BBSAN component, especially in the presence of a shock cell structure, as will be shown in the following subsection. The specific lines along which these cross correlations were conducted are depicted in Fig. 6. For the annular shear layer, computations were carried out on lines $L1$ and $L0$⁠, while for the circular shear layer, the analysis was performed on line $L3$⁠. A summary of the results for both annular and circular shear layers is presented in Table 9.

The convection velocity is higher in annular shear layer than in the circular shear layer, aligning with earlier observations for the aerospike nozzle jet [42]. The computed value for the circular jet $uc=0.72uj$ is in good agreement with findings from previous simulations of cold jets [30]. Moreover, the convection velocity decreases with the increasing swirl number as suggested by Carpenter [25]. Finally, the convection Mach number remains below one, which is consistent with observations indicating the absence of Mach waves in the near field. The results for the convection velocity in the azimuthal direction are not presented for brevity.

The far-field aeroacoustic signature is computed using the Ffowcs Williams–Hawkings equation. Flow data are exported on 448,400 nodes located on the FWH surface (shown as a thin white line in Fig. 2) every $2×10−6$ s. The original formulation is found in Ref. [49]. This approach has been widely adopted for jets and validated against experimental data with the current solver [30,38]. A detailed formulation for jets is given in Ref. [50]. The observation points are considered on an arc of radius $60Deq$ centered on the aerospike body tip. The spectra are averaged on the $xz$- and $yz$-planes.

At downstream side angles $θ∼150$ deg, mixing noise produced by large turbulent flow structures propagates at low frequencies ranging in the neighborhood of $St∼0.25$ in agreement with previous observations for both subsonic and supersonic jets [5,52]. Mixing noise is still present in the swirling cases, but at lower amplitudes compared to the baseline case. This phenomenon could be attributed to a loss of coherence, which is characteristic for swirling jets [23,24]. Moreover, the high shear in the vicinity of the nozzle trailing edge outlet affects the stability of the vortical structures generated in that region, leading to a faster decay. In particular, it was found that swirling enhanced the maximum amplification rate for compressible mixing layers [20]. These factors impede the convection of coherent large structures further downstream and lead to the attenuation of mixing noise.

BBSAN is due to the passage of vortical structures through the quasi-periodic shock cell structure. The central frequency of this sound component is given by Eq. (5). Carpenter generalized the results developed by Howe and Ffowcs Williams for broadband shock-associated noise [25,46]. It was shown that an additional swirling component should not fundamentally modify the spectrum compared to the nonswirling cases. The formula to predict the central frequency in the baseline case without swirling given in Eq. (5) applies for the swirling cases as well. The cross correlation for the axial velocity in the annular and circular shear layers provides the convection velocity of the vortical structures (see Table 9). Since the axial convection velocity in the annular shear layer decreases with the increasing swirl number, the central frequency for BBSAN emerging in the annular part of the jet is shifted toward lower frequencies with the increasing swirl numbers. This can be observed in the spectra in Fig. 10. A shock cell structure is observed downstream of the aerospike bluff body only in the baseline case and for $S=0.10$⁠. The axial convection velocity also decreases in the jet part located downstream of the aerospike bluff body for the case $S=0.10$ compared to the baseline case. This leads to a slight shift towards lower Strouhal numbers.

Tonal components at the Strouhal numbers of the azimuthal modes identified in Sec. 4.3 are observed in the far-field spectra, at $St=0.85$⁠, $St=0.85$⁠, and $St=0.90$⁠, for the cases $S=0.10$⁠, $S=0.20,$ and $S=0.30$⁠, respectively. These components propagate at upstream angles (⁠$θ∼30$ deg) as well as downstream (⁠$θ∼140$ deg). Additionally, the SPL for the azimuthal mode increases with the increasing Swirl number. A tonal component at slightly lower SPL close to $St∼1.70$ is observed for the case with $S=0.10$ propagating at right angles (⁠$θ∼80$ deg). This corresponds to a subharmonic of the main radial oscillation mode or azimuthal mode. Similar subharmonics are found for the case $S=0.20$ around $St∼1.70$ and for the case $S=0.30$ around $St∼1.80$⁠. Further subharmonics of the main radial oscillation modes are found at lower amplitudes $St∼2.55$ and $St∼3.40$ for the case $S=0.10$⁠, propagating at angles $θ∼145$ deg. High SPL are also found along and between the central frequency lines of the BBSAN component using the shock cell lengths in both parts of the jet for $S=0.10$⁠. A further peak is identified in Fig. 8 at the location of the separation bubble and at the junction between the two first shock cells. These peaks were identified at $St=0.56$⁠, $St=0.50$⁠, and $St=0.48$ for $S=0.10$⁠, $S=0.20$⁠, and $S=0.30$⁠, respectively. The far-field spectra for these cases display peaks at these Strouhal numbers at upstream angles $θ∼30$ deg. The intensity decreases with the increasing swirl number such that the peak is not as distinct for $S=0.30$⁠. Additional secondary tonal components identified in the PSD of the radial velocity (see Fig. 8) are also observed in the far field. These propagate at right downstream angles between $θ∼90$ deg and $θ∼140$ deg. Further tonal sound components at low Strouhal numbers are observed in the spectrum for $S=0.30,$ which were not present in the $S=0.20$⁠. These might correspond to the tonal components observed in the power spectral density for point $P1$ depicted in Fig. 8 at around $St∼0.17$ in the case $S=0.30$⁠.

The aeroacoustic signature in the baseline case displays more complex features compared to the swirling cases due to the observed shock cell structure. Screech noise and BBSAN are observed, as well as mixing noise. Imposing a swirling component to the flow has a direct impact on the far-field aeroacosutic signature of the jets. It has been shown that the screech noise component disappears with increasing Swirl number. Instead, the azimuthal modes are enhanced, leading to sound radiation in the far field at higher Strouhal numbers as displayed in Table 8, data computed with the two-point cross correlations. The SPL in the far field linked to this azimuthal mode increases with the increase in the Swirl number. Mixing noise is weakened due to the increased mixing. Additionally, the supersonic flow reattaches on the aerospike further downstream, leading to a frequency shift for the radial oscillation modes of the separation bubble and the annular shock cell structure. The Strouhal numbers shift from $St∼1.25$ in the baseline case to $St∼0.50$ in the swirling cases. Globally, the SPL decreases with the increasing Swirl numbers.

The study investigates the swirling flow effects on the cold supersonic jet exhausting a three-dimensional axisymmetric aerospike nozzle. Three different Swirl numbers are considered imposed as inlet flow boundary conditions. The results were compared with a baseline case without swirl. It was shown that the introduced swirl reduces the length of the shock cell structure found downstream of the aerospike plug at baseline (without swirl) or even suppresses it, thus eliminating the screech component. The annular shock cell structure is still preserved with slightly longer shock cells for the swirling cases compared to the baseline case. The flow reattachment toward the aerospike bluff body occurs further downstream. In particular, this affects the radial oscillation mode of the separation bubble and the annular shock cell structure, leading to radial oscillation modes at lower Strouhal numbers compared to the baseline case. The additional swirling motion enhances the azimuthal modes of the jet. Two-point space-time cross correlations of pressure data in the annular shear layer allowed to identify the Strouhal numbers of the azimuthal modes of the swirling jets. Additionally, the jet spreading is increased with the increasing swirl numbers, which suggests an enhanced mixing and lower sound pressure levels. The far-field aeroacoustic signature is computed using the Ffowcs Williams–Hawkings equation. The obtained spectra display both broadband noise, associated with large vortical structures and the passage of vortices through the quasi-periodic shock cell structure, as well as tonal components. The frequency of those tonal components is in agreement with the Strouhal numbers observed in the PSD of the azimuthal and radial oscillation modes. Furthermore, mixing noise at low Strouhal numbers decreased in intensity due to the loss of coherence for vortical structures. Finally, swirling was found to decrease the global sound pressure levels compared to the baseline case while maintaining reasonable thrust performance.

This project was funded by the “INSPIRE” EU Project H2020-MSCA-ITN-2020, Marie Skłodowska-Curie Innovative Training Networks, Project No. 956803. The computations were performed on resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS, Project NAISS 2023/1-19) at the PDC Centre for High-Performance Computing (PDC-HPC) and at the National Supercomputer Center (NSC) in Sweden. The authors would like to thank Dr. Stefan Wallin and Dr. Peter Eliasson for the assistance of the code implementations.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.