Numerical Analysis of Side-loads Reduction in a Sub-scale Dual-bell Rocket Nozzle

A calibrated delayed detached eddy simulation of a sub-scale cold-gas dual-bell nozzle flow at high Reynolds number and in sea-level mode is carried out at nozzle pressure ratio NPR = 45.7. In this regime the over-expanded flow exhibits a symmetric and controlled flow separation at the inflection point, that is the junction between the two bells, leading to the generation of a low content of aerodynamic side loads with respect to conventional bell nozzles. The nozzle wall-pressure signature is analyzed in the frequency domain and compared with the experimental data available in the literature for the same geometry and flow conditions. The Fourier spectra in time and space (azimuthal wavenumber) show the presence of a persistent tone associated to the symmetric shock movement. Asymmetric modes are only slightly excited by the shock and the turbulent structures. The low mean value of the side-loads magnitude is in good agreement with the experiments and confirms that the inflection point dampens the aero-acoustic interaction between the separation-shock and the detached shear layer.


Introduction
The dual-bell nozzle is a kind of one-step altitude compensating nozzle that may represent a possible alternative to replace conventional bell nozzles in future launcher first-stage engines. The main feature of this advanced concept is the particular shape of the divergent section, composed by two bells, namely the base and the extension, separated by an 1 3 inflection point, as shown in Fig. 1 (a). It has two main operating conditions: i) sea-level mode and ii) high-altitude mode. In the first mode, it operates with a low area ratio, with a controlled and symmetrical flow separation at the inflection point, thus avoiding the onset of dangerous side loads that can induce structural damages to the engine. According to Schmucker (1984), side loads are caused by a non-symmetrical shock-induced flow separation and the degree of this non-symmetric behavior is proportional to the inverse of the wall-pressure gradient magnitude. The inflection point in the dual-bell contour induces a Prandtl-Meyer expansion and, consequently, an infinite wall-pressure gradient (in the inviscid sense) thus zeroing the side-load magnitude. In the high-altitude mode, the flow attaches to the extension and the nozzle works with a higher area ratio, thus increasing the thrust coefficient. The parameter governing the transition between the two operating modes is the nozzle pressure ratio (NPR), i.e. the ratio between the stagnation pressure in the chamber and the ambient static pressure p 0 ∕p a .
The dual-bell nozzle was first proposed by Cowles and Foster (1949) and all the studies carried out since then highlighted three main critical issues: i) the transition between the two operating modes, ii) the detached flow unsteadiness in sea-level mode, and iii) hot flow behavior and cooling effects. The first experimental test campaign was performed by Horn and Fisher (1993) using cold-gas sub-scale dual-bell nozzles. They found that both the constant-pressure and increasing-pressure extension exhibited good transition characteristics with a transition time of the order of 30 ms. Several numerical studies (Wong and Schwane 2002;Nasuti et al. 2005) confirmed a quick transition behavior. In particular, Nasuti et al. (2002Nasuti et al. ( , 2005 found that the dual-bell extension with a positive pressure gradient can ensure a transition faster than the one with a constant pressure profile (CP). On the other hand, this profile can ensure higher thrust performances and for this reason a trade-off is necessary. After the experiments by Horn and Fisher, the transition phenomenon was studied through an extensive test campaign conducted at the German Aerospace Research Center (DLR) Stark 2009, 2010), that investigated the effect of different geometrical parameters on the transition process and on the side-loads generation in cold-gas subscale nozzles. They found that both modes (sea-level and high-altitude) were associated to a level of side loads lower than the ones suffered by a comparable truncated ideal contour (TIC) nozzle. Nevertheless, the transition was characterized by a short-time high-peak side load, that could jeopardize the nozzle structure.
A peculiar characteristic of the switching between the two operating modes was indicated by Nasuti et al. (2005), who found that at the beginning of the transition the separation front moves into an inflection region, located immediately downstream of the inflection point and characterized by a negative wall-pressure gradient (due to viscosity) like conventional nozzles. As a consequence, the dual bell can start to experience non-axial forces before the full transition takes place. The length of the inflection region was found to be of the order of the throat radius and, as shown by Martelli et al. (2007), it depends on multiple factors, including the boundary-layer thickness at the end of the base, the Prandtl-Meyer expansion fan at the inflection point and the design wall-pressure gradient imposed to the extension during its construction with the method of characteristics (Nasuti et al. 2005). Génin and Stark (2011) confirmed experimentally the existence of the inflection region and that during the transition of the separation point in this zone the level of the side loads is similar to that suffered by a conventional nozzle. Verma et al. (2015) studied the flow unsteadiness when the separation front was located in the inflection region by analyzing the spectral content of the wall-pressure signature. They related the onset of side loads during the transition in the inflection region to the high level of unsteadiness suffered by the flow in this regime. Another important aspect considered in literature was the effect of the launcher-base flow on the nozzle internal flow behavior. In this context, the investigation of the complex interaction of a dual-bell exhaust flow with the unsteady external flow is of particular interest, especially when the nozzle works in sea-level mode. As observed and demonstrated (Perigo et al. 2003;Torngren 2002;Wong 2005), an external pressure fluctuation can cause an internal amplification of the pressure oscillations, which can be dangerous for the side-loads generation. Verma et al. (2014) suggested that the external unsteady perturbation could lead to a transition/re-transition cycle, generally known as flipflop effect. Loosen et al. (2019) numerically studied the interference of a turbulent wake coming from a generic space launcher with a dual-bell nozzle exhaust flow in sea-level mode. It was found that in supersonic freestream conditions, the presence of the outer flow leads to a premature transition, thus reducing the transition NPR.
Most of the studies available in literature refer to cold-gas flow investigations. In real flight conditions however, the dual-bell nozzle would work with hot exhaust flow and a film cooling technique could be used to protect the nozzle structure from thermal fatigue. Génin et al. (2013) performed one of the first experimental test campaign in order to study the nozzle behavior with inert hot air flow. The presence of the inflection point induces thermal loads during the transition process and the generated heat flux raises in the inflection region. Martelli et al. (2009) studied the effect of a film cooling injected in the base nozzle near the inflection point by means of numerical simulations. The emerging results highlight the efficiency of the film which is strongly influenced by the expansion fan at the inflection point. Unfortunately, from the study also emerges that the secondary gas increases the inflection region extension, leading to an increased risk of side-loads generation.
From the analysis of the available literature, it is clear that the development of side loads inside the nozzle is one of the main critical aspects that must be taken under consideration for the development and the realisation of a dual-bell nozzle. This manuscript focuses on the numerical investigation of the development of lateral loads inside an over-expanded dual-bell nozzle working in the sea-level mode. From a numerical point of view, the best approach to resolve the flow-separation dynamics and the shock-wave/boundary-layer interaction would be the large eddy simulation (LES) technique. But performing a wall-resolved LES of this high Reynolds numbers flow ( Re ≈ 10 7 ), requires an impractically computational effort. On the other hand, modeling the global effect of the turbulent scales as done in the URANS approach could suppress the important flow processes leading to the formation of the aerodynamic unsteady loads. In this scenario, a possible and valid choice is represented by the use of a hybrid RANS/LES methodology (Spalart et al. 2006;Spalart 2009;Fröhlich and von Terzi 2008;Weinmann et al. 2014), whose rationale is to simulate attached boundary layers in RANS mode, while separated shear layers and turbulent recirculating zones are solved by the LES mode. Among the different methods, in this work we chose to adopt the detached eddy simulation (DES) technique (Spalart et al. 1997) and, in this framework, few test cases of separated dual-bell nozzle flows are reported in the literature. In particular, Proschanka et al. (2012) performed a numerical study of a cold-gas dual-bell nozzle in sea-level mode, founding three types of pressure fluctuations: one symmetric and two asymmetric, the latter being associated with side-loads generation.
In this work, a calibrated version of the delayed detached eddy simulation (DDES) of a sub-scale cold-gas dual-bell nozzle in the first operating mode is carried out. The geometry is inspired by the experimental work of Verma et al. (2015) and it is characterized by a high nozzle-exit Mach number, closer to a real application than the dual bell used by Proschanka et al. (2012). The simulated NPR is experimentally characterized by a low level of side loads, and the separation front is located at the inflection point, at the very beginning of the inflection region. The present analysis focuses on the investigation of the wall-pressure signature and its spectral content, in particular in the azimuthal wave-number frequency plane, with the main aim of identifying a correlation between the forced flow separation and the energy level of the non-symmetrical azimuthal mode. By performing experimental and numerical analysis of TIC nozzles, Jaunet et al. (2017) and Martelli et al. (2020) have recently proposed that the side loads could be triggered by an aeroacoustic feed-back loop (screech-like phenomenon) involving the shock-cell structure and the separated shear layer present inside the nozzle. The scope of this study is to show how the dual-bell inflection point changes the receptivity of the separation shock to the acoustic disturbances coming from the interaction between the supersonic shear layer and the second Mach disk, thus altering the aeroacoustic feed-back loop observed in the TIC. This is demonstrated through an in-depth analysis of the flow organization, wall-pressure signature and side-loads level in the dual bell nozzle, that are compared with those relative to the TIC nozzle studied by Martelli et al. (2020). It is worthwhile to underline that, in real applications, conventional bell nozzles (like the TIC) operate with overexpanded flows without internal flow separation. Nevertheless, during the ignition transient they do suffer lateral forces, which must be foreseen during the design phase. Therefore, a deep physical understanding of the effectiveness of the inflection point in reducing the sideloads level could be important to develop more performant rocket nozzles. The paper is organised as follows. First, the methodological approach and the computational setup are presented in Sects. 2 and 3. After providing an overview of the flowfield organisation in Sect. 4.1, the main features of the wall-pressure signature are analysed in Sects. 4.2 and 4.3 by means of spectral analysis (through Fourier transform in time and wavenumber-frequency decomposition) and compared with experimental data. Then the aerodynamic side loads generated by the nozzle are inspected and compared with experiments and numerical data of a TIC nozzle in Sect. 4.4. Conclusions are finally given in Sect. 5.

Methodology
In the present work we consider the governing equations for a compressible, viscous and heatconducting gas, that in three dimensions can be written in conservation form as 1 3 where is the density of the flow, u i denotes the velocity component in the i-th coordinate direction ( i = 1, 2, 3 ), E is the total energy per unit mass and p is the thermodynamic pressure. ij and q j denote the total stress tensor and total heat flux, that are given by the sum of a viscous and a turbulent contribution, according to where the Boussinesq hypothesis is applied through the introduction of the eddy viscosity t (defined by the turbulence model), S * ij is the traceless strain-rate tensor and the kinematic viscosity, depending on temperature T through the Sutherland's law. The molecular Prandtl number Pr is constant and equal to 0.72. Similarly, a constant value for the turbulent Prandtl number Pr t = 0.9 has been considered, as commonly assumed in DES/LES of high-speed flows (Garnier et al. 2002).
The numerical methodology employed in this work is a calibrated version of the delayed detached eddy simulation method (Spalart et al. 2006), an efficient and powerful hybrid RANS/LES technique well suited to capture the unsteadiness of high-Reynolds number turbulent flows. The current implementation is based on the Spalart-Allmaras (SA) turbulence model, which solves a transport equation for a pseudo eddy viscosity w here d is the model length scale, f w is a near-wall damping function, S is a modified vorticity magnitude, and , c b1 , c b2 , c w1 are model constants. The pseudo eddy viscosity ̃ is directly linked to the eddy viscosity t through t =̃f v1 , where the correction function f v1 is used to guarantee the correct near-wall boundary-layer behavior. In the DDES approach the turbulence model automatically switches between a pure RANS mode, active in flow regions with attached boundary layers, to a pure LES mode, active in flow regions detached from the wall, where the computation can directly resolve the large scale, energy containing eddies. This objective is achieved by defining the length-scale d in (3) as where d w is the distance from the nearest wall, is the subgrid length-scale that controls the wavelengths resolved in LES mode and C DES is a calibration constant set equal to 0.65 in the original model. On the basis of previous calibration studies performed on DDES of a sub-scale rocket nozzle Martelli et al. (2019Martelli et al. ( , 2020) the constant C DES has been set equal to 0.20 and the function f d here employed is where U i,j is the velocity gradient and k the von Karman constant. The f d function, that is built in such a way that its value is 0 in boundary layers and 1 in LES regions, represents the main difference between the DDES strategy and the original DES approach (Spalart et al. 1997), denoted as DES97. It guarantees that attached boundary layers are always treated in RANS mode, even in the case of extremely fine grids, thus allowing to alleviate the well-known phenomenon of modeled stress depletion, which in turn can lead to gridinduced separation Spalart et al. (2006). The sub-grid length scale is specified according to the formulation proposed by Deck (2012), and it depends on the flow itself, through f d as 1∕3 . The main idea of this formulation is to take advantage of the f d function to switch between max , needed to shield the boundary layer, and vol , needed to ensure a rapid destruction of modelled viscosity to unlock the Kelvin-Helmholtz instability and accelerate the passage to resolved turbulence in the separated shear layer.

Flow Solver
The simulations have been performed by means of an in-house, compressible, finite-volume, structured solver, widely employed in the past to investigate the dynamics of turbulent, separated flows in transonic and supersonic rocket nozzles, involving complex shockwaves/boundary-layer interactions (Martelli et al. 2017(Martelli et al. , 2019(Martelli et al. , 2020Memmolo et al. 2018). In the flow regions away from the shock, the spatial discretization consists of a centered, second-order, energy consistent scheme, that makes the numerical method extremely robust without the addition of numerical dissipation (Pirozzoli 2011). This feature is particularly useful in the flow regions treated in LES mode, where in addition to the molecular, the only relevant viscosity should be that provided by the turbulence model. Strong compressions in the flow are identified by means of the Ducros shock sensor (Ducros et al. 1999), that is used to switch the discretization of the convective terms of the governing equations to third-order Weighted Essentially Non Oscillatory reconstructions for cell-faces flow variables. The viscous fluxes are evaluated through compact, second-order central-difference approximations. A low-storage, third-order Runge-Kutta algorithm (Bernardini and Pirozzoli 2009) is used for time advancement of the semi-discretized ODEs' system. The code is written in Fortran 90, it uses domain decomposition and it fully exploits the message passing interface (MPI) paradigm for the parallelism.

Test Case Description
The present investigation is carried out on the dual-bell geometry experimentally tested at DLR (Verma et al. 2015), characterized by a throat radius r t = 0.01 m, a wall inflection angle .3 , an extension area ratio E = r 2 e ∕r 2 t = 27.1 , with the extension displaying a constant wall-pressure profile (CP). The nozzle pressure ratio is fixed at the value 45.7, since at this NPR the nozzle flow experiences a symmetric shock separation located exactly at the inflection point. The other parameters of the simulation were chosen to reproduce the operating conditions of the experimental campaign. In particular, the total temperature T 0 and the static ambient temperature T a have been set equal to 300 K. The nozzle Reynolds number is where 0 is the density, a 0 the speed of sound and 0 = (T 0 ) is the molecular viscosity taken at the stagnation-chamber condition. The three-dimensional computational domain includes the external ambient (see Fig. 2) that extends up to 150 r t in the longitudinal direction and 80 r t in the radial one from the symmetry axis. The boundary conditions are imposed as follows: total temperature, total pressure and flow direction are enforced at the nozzle inflow, while on the outflow boundary at the end of the external domain a nonreflecting boundary condition is implemented. A back pressure equal to the ambient pressure p a is imposed at the other boundaries. The no-slip adiabatic condition is prescribed for the nozzle walls. The mesh resolution was selected following a preliminary sensitivity study carried out on steady-state axisymmetric RANS computations, for which convergence of the separation location was obtained. The RANS solution was also used to initialise the three-dimensional DDES simulation. The development of turbulent structures and the passage from modelled to resolved turbulence was triggered by adding random perturbations to the streamwise velocity field only at the initial time of the 3D computation. Those perturbations have a maximum magnitude of 3% of the inflow velocity and a flat spectrum. We point out that, in our approach, no resolved turbulence upstream of the separation is enforced. The computational domain is composed by 136 structured blocks in the x-y plane, each block discretized with 22 × 256 cells. The 3D mesh includes 256 cells in the azimuthal direction, for a total number of approximately 196 million cells. The computation was run with a time step t = 4.65 ⋅ 10 −8 s and a relatively long time span was simulated T = 0.0172 s, which guarantees coverage of frequencies down to at least f min ≈ 58 Hz. A total of 80 full three-dimensional fields have been collected at time intervals of 2.15 × 10 −4 s for post-processing purposes. Furthermore, samples of the pressure field at the wall and in an azimuthal plane have been recorded at shorter time intervals of 2.34 × 10 −6 s to guarantee sufficient resolution for the frequency analysis. In the following, the results for the dual bell are compared with those recently obtained for a truncated ideal contour (TIC) nozzle with free-shock separation (FSS) operating at NPR = 30.35, experimentally tested at the University of Texas at Austin and simulated with the same methodology here exposed (Martelli et al. 2020). The mesh resolution of the TIC computation was similar to that here employed for the dual-bell nozzle, being the two cases characterized by a similar Reynolds number. Basic geometric properties of the TIC nozzle are a throat radius of r t = 0.019 m, an exit radius of r e = 0.117 m and a throat-to-exit length of L = 0.351 m. A summary of the main geometrical parameters of both the dual-bell and TIC nozzle is reported in Table 1. An extensive analysis on the flow unsteadiness and wall-pressure signature for the TIC nozzle can be found in our previous paper (Martelli et al. 2020). Here, novel results for that geometry are presented concerning the analysis of the aerodynamic side loads, compared in Sec. 4.4 with those generated in the dual-bell nozzle.

Flowfield Organisation
The main features of the flow pattern inside the dual-bell nozzle are reported in Fig. 3, which shows a longitudinal x − r plane with the iso-contours of the mean Mach-number field, obtained by averaging in time and in the azimuthal direction. At the selected NPR, the flow is separated and anchored at the inflection point, generating a free shear layer. A shock system arises comprising the conical separation shock, a Mach-disk reflection and a second conical shock, which re-directs the shear layer in a direction almost parallel to the nozzle axis. The dual-bell extension is characterized by a considerable subsonic turbulent recirculating region. The flow downstream of the Mach disk is initially subsonic, then it experiences an expansion process across a fluid-dynamic throat, that is visible from the sonic line, and again accelerates to a supersonic velocity. Then a new shock appears to balance the pressure of the flow to the ambient level. It is also interesting to observe the presence of a recirculating region near the nozzle axis downstream the Mach disk. This region is fed by the vorticity produced by the shock curvature, which causes an entropy gradient in the radial direction and hence vorticity according to the Crocco's theorem. The unsteadiness of the flow is highlighted by Fig. 4 (a), where contours of the densitygradient magnitude are displayed in a longitudinal x − r plane. The developing turbulent structures of the two co-annular supersonic shear layers are well visible. The external shear layer originates from the jet detachment from the wall and the internal one originates from the triple point of the Mach reflection. The two layers merge downstream in the external ambient, where the Mach waves irradiated by the vortices are well visible. It is also possible to observe fine turbulent structures emitted by the Mach disk, coherently with the observation reported above on the entropy gradient generation due to the shock bending. The Q-criterion is generally used to identify the tube-like structures from a qualitative point of view (Hunt et al. 1988). In this work a modified definition of the Q-criterion has been adopted to include the effects of compressibility (Pirozzoli et al. 2008). An isosurface of the Q-criterion is shown in Fig. 4 (b), coloured according to the value of the streamwise Fig. 4 Visualisation of (a) the instantaneous density-gradient magnitude (numerical Schlieren) in an x-y plane and of (b) turbulent structures through an iso-surface of the Q criterion, coloured with the local value of the longitudinal velocity 1 3 velocity component. It shows that the Kelvin-Helmotz instability of the initial part of the shear layer is not characterised by the coherent toroidal vortices which can be found in incompressible flows. Instead, we observe that oblique modes dominate the initial part of the shear layer, then leading to the generation of small-scale three-dimensional structures, well resolved by the present DDES approach. The differences in the initial shear layer development observed with respect to the typical pattern of low-speed flows can be attributed to the large local convective Mach number ( M c ≈ 1 ) at the beginning of the detached shear layer, which changes the shear-layer instability process (Sandham and Reynolds 1991). A similar behavior was observed by Martelli et al. (2020) in a conventional TIC nozzle with flow separation and by Simon et al. (2007) in a supersonic cylindrical base flow.

Analysis of the Wall-Pressure Signature
The spatial evolution of the mean wall pressure ( p w ∕p a ), averaged in time and in the azimuthal direction, is presented in Fig. 5(a), compared with the reference experimental data (Verma et al. 2015). The mean wall pressure decreases until the separation point, which is located at the inflection point, then the oblique shock causes a sudden increase up to a plateau value close to the ambient pressure level. The numerical data are in good agreement with the experimental ones. In particular, the right values of wall pressure in the turbulent recirculating zone indicates that the LES-branch simulation is able to correctly capture the momentum exchange between the main jet and the separated flow, dominated by the high convective Mach number of the shear layer (Pantano and Sarkar 2002). The standard deviation of the wall-pressure signal ( w ∕p w ) is shown in Fig. 5(b) and compared with the corresponding experimental values. The distribution of w ∕p w along the longitudinal axis is typical of shock-wave/turbulent boundary layer interactions (Dolling and Or 1985), characterised by a dominant sharp peak at the separation-shock location and followed downstream by a lower and almost constant level in the zone where the turbulent shear layer develops. The numerical standard deviation agrees with the experimental behavior, especially in the turbulent recirculating zone. It seems that the sharp peak predicted by the DDES at the shock location is not captured by the experimental measurements, probably due to the probe spacing. Indeed, the flow separation is fixed at the inflection point, therefore the pressure fluctuations are confined in a very narrow spatial range. The dual-bell standard deviation of the wall-pressure fluctuations is also compared with that of the TIC nozzle in Fig. 6. It is possible to note that the TIC distribution is characterised by a higher value in correspondence of the peak and by a more marked drop downstream of the shock location. Then the energy level of the pressure fluctuation gradually increases during the development of the shear layer. The dual-bell nozzle instead has a more flat trend downstream of the peak.
To provide a global picture of the pressure energy distribution along the nozzle wall and in the frequency domain, contours of the premultiplied wall-pressure spectra G(f ) f ∕ 2 are reported in Fig. 7(a) with respect to the longitudinal coordinate x∕r t and frequency f. From the spectral map, it is possible to observe the presence of two regions characterised by high fluctuation energy. The first zone of interest is located near the separation-shock average location ( x∕r t ≈ 6.3 , see Fig. 3) and is characterised by a broad bump in the low-medium frequency range and very narrow in the longitudinal direction, that is clearly associated to the signature of the unsteady shock motion. The peak of this bump is at a frequency of ∼ 830 Hz and its footprint is well visible in the spatial direction until the end of the nozzle. According to Baars et al. (2012) and Martelli et al. (2020) the low-frequency peak can be attributed to an acoustic resonance, which can be described by the one-quarter acoustic standing wave model (Wong 2005) in open-ended pipes. The resonance frequency f ac can be expressed as follows, where a ∞ = 345m/s is the ambient speed of sound and L = l e = 0.0804 m is the distance between the inflection point and the nozzle lip. This distance should correspond to the pipe length in the 1D acoustic model and is an opened-end empirical correction parameter (Proschanka et al. 2012). According to Proschanka et al. (2012) this parameter can vary between 0.5L and 0 giving a range of frequency between 715 Hz and 1073 Hz, which include the numerical peak frequency of 833 Hz. It should also be considered that the presence of a large turbulent recirculating region, the unsteady detached shear layer and the area ratio variation in the second bell increase the difference of the present situation with the acoustic pipe 1D model. The second region characterized by high levels of pressure fluctuations is located downstream of the separation shock, in a high-frequency range between f = 10 4 and 10 5 Hz. The spectral map highlights that this second (broad) peak at high-frequency gradually appears moving in the downstream direction, and it is associated with the development of the separated shear layer, from which pressure disturbances are emitted through vortices that are convected downstream at high speed. For a better quantification of the evolution of the wall-pressure frequency content along the longitudinal axis, Fig. 7(b) shows the premultiplied spectra at six axial equispaced probes. The peaks associated to the shock movement have very high amplitudes at ≈ 833 Hz and they persists for almost all the stations considered. Starting from the first probe the energy content of the shock peak tends to reduce along the x-axis whereas the energy contribution of the separated shear layer rises in amplitude at higher frequencies. Table 2 reports the peak frequencies of the spectra (attributed to the acoustic resonance) at the three different axial stations corresponding to the location of the experimental probes (black dashed lines in Fig.7(a)). There is a rather good agreement between the numerical results (833 Hz) and the experiments (between 775 Hz and 787 Hz), since the deviation from the experimental values falls within the uncertainty of the numerical spectra ( f min = 58 Hz). In addition Verma et al. (2015) shows that the acoustic tone persists all along the second bell also in the experiment. Figure 8 shows the frequency spectra of the wall-pressure fluctuations corresponding to the TIC nozzle, as also reported in Martelli et al. (2020). In this case, to provide a further experimental-numerical validation, premultiplied spectra are reported at three axial stations corresponding to positions where experimental data from dynamic Kulite transducers are available. The main quantitative difference is that the footprint of the shock oscillation decreases much rapidly in the longitudinal direction with respect to the dual-bell case. Furthermore, the TIC wall-pressure signature presents some non-negligible energy in the intermediate frequency range (around 1 kHz) visible in both the numerical and experimental spectra, whereas the dual-bell nozzle show a very small energy content in the medium frequency range (around ≈ 2300 Hz). As previously speculated by Jaunet et al. (2017) and then confirmed by Martelli et al. (2020), this intermediate-frequency peak can be attributed to a screech-like mechanism with an antisymmetric-type mode occurring inside the nozzle. This aspect is deeply discussed in the following section.

Azimuthal Decomposition of the Pressure Field
The proper way to correlate the frequency and energy content of the wall-pressure signature with the aerodynamic side loads is to evaluate the Fourier azimuthal wavenumberfrequency spectrum, defined as whit R pp indicating the space-time correlation function of the wall-pressure fluctuations: where x and are the spatial separations in the streamwise and azimuthal directions, is the time delay, and ⟨ ⟩ denotes averaging with respect to the azimuthal direction (exploiting homogeneity) and time. The symbol m indicates the mode number in the azimuthal direction, also called Fourier-azimuthal wavenumber and it is noteworthy to remember that  the asymmetric mode m = 1 is the only one capable of providing a contribution to the side loads, as also shown in Fig. 9 where the real part ( m > 0 ) of the first five Fourier azimuthal modes is reported. Following Baars et al. (2012), we introduce the total resolved energy (x) for each longitudinal station of the nozzle as: where (m) (x) is the variance of the m th time-dependent Fourier-azimuthal mode coefficient. On the base of these definitions it is possible to show in Fig. 10 the eigenspectra of the Fourier modes as fractions of the resolved energy per mode (x, m): for two axial stations for both nozzles. The axial distances are evaluated from the separation location and normalized with throat radius for comparison purposes. Starting with the dual-bell nozzle case, Fig. 10(a), it is evident that the energy of the zeroth mode is by far the dominant one along the nozzle. The energy of the mode m = 1 is nearly negligible near the separation location and shows an increment at the downstream station, but always remaining a small fraction (1.8% at most) of the total resolved energy. The higher modes appear to be non-negligible (but always less than 1%) only at the downstream station. The picture highlighted by Fig. 10(b) for the TIC nozzle is rather different. Indeed, it is possible to appreciate a non-negligible decrease of the energy of the zeroth mode along the nozzle and an important growth of the first mode, that at (x − x sep )∕r t = 7.03 provides a 30% contribution to the total resolved energy. The energy of the higher modes is still negligible near the shock location, giving a minimal contribution only at the downstream station. To evaluate the frequency content of the time-dependent Fourier azimuthal coefficients, their power spectral densities are computed for both nozzles. Figure 11 shows the spectra for the symmetric ( m = 0 ) and asymmetric ( m = 1 ) modes at the same two axial stations of Fig. 10. To compare the two different nozzles, a Strouhal (St) number is introduced as in Tam et al. (1986): where U j is the fully expanded jet velocity  16r t (a, b) and 7.03r t (c, d). The black solid line refers to the dual-bell nozzle while the red solid one to the TIC. For visualization purposes, different scales are applied to the various figures T 0 is the stagnation temperature, the specific heat ratio, R the air constant, M j is the fully adapted Mach number, which is a function of the nozzle pressure ratio through the isentropic relation. The length scale D j is computed as a function of M j , the design Mach number M d and the nozzle exit diameter D through the mass flux conservation Bearing in mind the picture of the wall-pressure power spectral densities (Figs. 7(a) and 8(a)), we see in Fig. 11(a) that the high energy peaks characterising the shock region, at 833 Hz ( St = 0.057) for the dual-bell nozzle and at 315 Hz ( St = 0.038) for the TIC nozzle, belong to the symmetric mode. In the dual-bell case, the zeroth mode maintains its importance along the nozzle, while an energy decrease can be observed in the TIC case, Fig. 11(c), coherently with the contours in Fig. 8(a). As far as the asymmetric mode is concerned, we see in Fig. 11(b) and (d) that the dual bell shows some energy around ≈ 2300 Hz ( St = 0.16) and that this peak persists downstream. The same picture appears for the TIC case, but with a much higher energy contribution: near the shock it appears a clear tone at approximately 1 kHz ( St = 0.12). Moving downstream, the peak at St = 0.12 is still the most important but now some energy appears, as in the dual-bell case, at higher frequencies (around St ≈ 0.6), indicating the presence of the asymmetric mode in the turbulent shear layer. As already shown by Baars et al. (2012), Jaunet et al. (2017) and Martelli et al. (2020) it seems clear that the flow separation in a TIC nozzle is characterised by a shock movement which is mainly symmetric (piston-like oscillation) but with some energy present in the antisymmetric mode. This last mode, together with the higher ones, characterises also the detached shear layer, as shown by its footprint on the wall-pressure signature. This scenario is quite different from that observed in the dual-bell case, where the energy of the anti-symmetric mode is almost negligible with respect to that of the symmetrical one. All the peaks frequencies and St numbers of Fig. 11 are reported in Table 3. It is also adopted a modified Strouhal number St * = f L U j , where L is the distance between the mean separation point and the nozzle exit plane as defined before. From the data reported in the table it can be seen that L is, as expected, a better length scale for the acoustic resonance and in fact, the peak St * 's for the m = 0 mode are very similar for the two nozzles. On the contrary, it seems that D j is a better length scale for the antisymmetric mode and now the St's are closer.
The very low energy content of the first azimuthal mode in the dual-bell nozzle, with respect to the TIC case, is an indication that the dual bell working in sea-level mode at the present NPR should not develop a significant level of side loads. Jaunet et al. (2017) and   Martelli et al. (2020) speculated the existence of a screech-like mechanism inside conventional nozzles, with free shock separation, associated to this antisymmetric mode. Martelli et al. (2020) proposed a path for the feedback loop, involving the downstream propagation of hydrodynamic instabilities along the shear layers and in the nozzle core downstream the Mach disk and the upstream traveling of the acoustic waves, generated by the interaction of the vortices with the shock cells, in the subsonic turbulent recirculating region inside the nozzle. Indeed, according to Powell (1953a), Powell (1953), every aeroacoustic resonance in high-speed jets can be decomposed into four main processes: i) the downstream propagation of energy through hydrodynamic instabilities, ii) an acoustic generation process in which the downstream perturbations are converted into upstream perturbations, iii) the upstream propagation of these disturbances, iv) the generation of new hydrodynamic instabilities through the forcing of a sensitive point excited by the upstream acoustic waves (the nozzle lip in external screech). The last process is generally called receptivity process (Mitchell 2019). In the TIC nozzle, this process involves the forcing of the separation line and of the separation shock, which induces the generation of a new hydrodynamic instability. In the dual-bell nozzle the antisymmetric mode is only slightly excited due to the presence of the inflection point that anchors the separation shock and constrains its movement. In particular, the inflection point induces the shock to move symmetrically, thus modifying the efficiency of the receptivity process of the feedback loop, which hampers the generation of the first azimuthal mode. A quantitative analysis and comparison of the side-loads generation in the dual-bell and TIC nozzles is reported in the following section.

Aerodynamic Loads Distribution
The wall-pressure distribution is integrated along the nozzle wall (longitudinal and azimuthal direction) to obtain the time distribution of the side-loads components ( F y and F z ) and magnitude ( |F| = √ F 2 y + F 2 z ), which are reported in Fig. 12. The thrust and side-loads components are calculated by integration of the wall-pressure field as

3
The figures highlight that the side-loads components oscillate around the zero mean value and have a randomic distribution in time. Consistently with the previous discussion on the energy of the asymmetric m = 1 mode, the magnitude of the aerodynamic side loads in the TIC nozzle is found to be remarkably higher than that observed in the dual-bell nozzle. The experimental side-loads time history is not available from Verma et al. (2015), therefore, the numerical results are compared with the experimental data of Génin and Stark (2011), where a test campaign on a similar subscale dual-bell nozzle (named DB2) was performed. This nozzle has a slightly different area ratio of the second bell ( e = 24 ) and a different inflection angle ( i = 5 • ). The simulation results show an averaged value of the side loads of 3.26 N, which well agrees with the value of 3.42 N obtained by the experiments. The time evolution of the side-loads vector direction ( cos( ) , sin( ) ) with respect to the z axis, is reported in Fig. 13 as a function of time. This figure shows that the side loads have not a preferential direction in the space, similarly to the results obtained by Deck and Guillen (2002), who performed a numerical campaign to predict side loads in a sub-scale TIC nozzle. The probability density functions of F y (t) , F z (t) and |F|(t) are shown in Fig. 14(a-c) respectively for both nozzles. It clearly appears that the side-loads components are two independent normal random variables with zero mean value and the same standard deviation, thus indicating that their distribution is Gaussian. Indeed, the numerical distribution well fit the theoretical Gaussian one. As a consequence, the probability distribution of |F|, Fig. 14(c), is a Rayleigh distribution, as shown by Dumnov (1996), who studied the experimental statistical distribution of side loads in a sub-scale nozzle.
To better characterize the origin of the aerodynamic side loads, and identify the nozzle regions that most contribute to their generation, the axial distribution of the side-loads magnitude averaged in time is shown in Fig. 15 for the dual-bell and the TIC geometry. These profiles have been obtained by integrating the wall-pressure, at each streamwise location, along the azimuthal direction and averaging in time. As far as the dual bell is concerned, Fig. 15(a) shows a peak in correspondence of the separation-shock location, followed by a gradual increase until the end of the second bell, associated with the development of the turbulent shear layer. This behaviour is in close agreement with the observations made on the spatial distribution of the first Fourier azimuthal mode. The axial distribution of the side loads in the TIC nozzle, Fig. 15(b), further confirms the larger intensity of the lateral force for this type of nozzle. In this case, a strong peak is found in the proximity of the shock location, whereas a flat contribution is observed in the divergent section of the nozzle. It is worth to notice that these trends are qualitatively similar to that obtained by Deck and Guillen (2002). The Fourier spectra of the lateral forces F y and F z for both nozzles are reported in Fig. 16. It is clearly visible that both components in the TIC nozzle are characterized by a predominant peak at a Strouhal number approximately equal to 0.12, corresponding as expected to the peak frequency of the first Fourier azimuthal mode (see Table 3). The dual-bell nozzle shows a picture a little bit different: the spectrum of F y is characterized by a rather large bump with the highest peak centered at S t = 0.16 , that is the characteristic peak frequency of the anti-symmetric mode, while the spectrum of F z shows a bump split in two peaks at S t = 0.1 and S t = 0.16 . The difference in the behavior is most probably due to the fact that the anti-symmetric mode is less energetic and coherent than the one characterizing the TIC nozzle. Finally, the dimensionless longitudinal component of the force F x ∕(p c A t ) , corresponding to the thrust coefficient, is reported in Fig. 17(a) for both geometries. The two nozzles are characterised by similar levels of thrust oscillation. In particular, for the DB nozzle the fluctuation of the thrust coefficient from its mean value reaches a maximum departure of 1.37% and a standard deviation value of 0.74%. For the TIC nozzle these fluctuations have a standard deviation of 0.42% with a maximum departure of 1.02%. The value of the thrust coefficient fluctuation varies depending on several factors, including the operating regime. This behaviour can be observed, for example, in the study of Proschanka et al. (2012) on a sub-scale dual-bell nozzle, where the standard deviation of the thrust coefficient, calculated with respect to the NPR, is around 2.8%. This is valid both in the up-ramping and down-ramping process. A similar value is obtained by Moríñigo and Salvá (2006) for a sub-scale J-2S TOC nozzle, where the deviation of the St * . 17 a Time history of the thrust coefficient F x ∕(p 0 A t ) for the dual-bell and TIC nozzles and (b) corresponding premultiplied power spectral densities. The black solid line denotes the dual-bell nozzle, the red solid line refers to the TIC nozzle thrust oscillation amplitude is around 2%. The spectra of the thrust coefficients reported in Fig. 17(b) show that all the energy is concentrated at the frequency corresponding to the acoustic resonance given by the piston-like shock movement, St * = 0.10 and 0.13 for the dual-bell and TIC nozzle, respectively. This implies that, in the design phase, the level of vibrations at lift-off and during the first part of the trajectory (first operating mode), that are detrimental for the engine structure and the payload, must be carefully addressed.

Conclusions
A calibrated delayed detached eddy simulation of a dual-bell nozzle working in an overexpanded condition with flow separation anchored at the wall-inflection point has been carried out. The main goal of the investigation was to figure out the effectiveness of the wall discontinuity as a flow-separation control device. The study has been conducted by analyzing the spectral content of the wall-pressure signature, with the purpose of evaluating the aerodynamic loads and comparing them with those experienced by a truncated ideal contour (TIC) nozzle (Martelli et al. 2020). The latter is a well-known conventional nozzle, that suffers the onset of side loads when operates with an internal flow separation. The dual-bell nozzle geometry has been specified following the experimental work of Verma et al. (2015), whose data have been used to assess the accuracy of the simulation. The analysis of the unsteady wall-pressure signals showed a good agreement of the mean wall pressure and the standard deviation of the pressure fluctuations along the nozzle with the experimental data, although the standard deviation peak at the shock location seemed to be over-predicted. The spectral analysis highlighted the presence of an acoustic tone at ≈0.8 kHz, associated with the shock motion. This tone appeared to persist all along the nozzle wall, as also confirmed by the experimental trend, and it was found to be an acoustic resonance, its frequency being associated with a one-quarter standing wave that develops between the mean separation shock location and the nozzle lip. Also the TIC nozzle was characterised by a strong tone in the low frequency range, even if its energy decreases along the nozzle wall. Both nozzles also showed the presence of oscillation energy in the high frequency range (between 10 and 100 kHz), originated by the turbulent detached shear layer. The analysis carried out in the wavenumber-frequency space revealed that the lowfrequency peaks of both geometries are associated to the zeroth azimuthal Fourier mode, that is symmetric, and its energy represents the dominant contribution to the total fluctuation energy. The second important contribution came from the first (antisymmetric) mode, which is the only one to be non-symmetric and capable to trigger lateral forces. In particular, the TIC nozzle was characterised by a tone in the intermediate frequency range, ( ≈ 1 kHz, St = 0.12), as already shown in literature (Baars et al. 2012;Jaunet et al. 2017;Martelli et al. 2020). The first Fourier azimuthal mode is instead only slightly excited in the dual-bell nozzle with a peak at approximately 2300 Hz (St = 0.16). Its energy is almost negligible with respect to the zeroth mode and much lower with respect to the TIC case. Therefore, the triggering of the antisymmetric mode seems to be suppressed by the presence of the inflection point, which mainly forces a symmetric shock movement, altering the receptivity process of the separation line to the upstream travelling acoustic disturbances. The evaluation of the aerodynamic lateral forces revealed a qualitatively similar distribution for both geometries, in particular they follow the Rayleigh distribution as found by Dumnov (1996). The mean magnitude of the lateral force in the dual bell well agrees with the experimental value of Génin and Stark (2011) and it is an order of magnitude lower than the side loads found in the TIC nozzle. On the other hand, the fluctuation of the thrust coefficient is similar in both cases and of the order of 1%. This aspect should be considered in the design phase and carefully evaluated to avoid annoying vibrations on the payload and engine structure. To this purpose a more extensive comparison spanning a wider range of parameters in terms of nozzle pressure ratios, wall-inflection angles and second-bell design wall-pressure gradients, could be certainly interesting for a future study.