Keywords: choc, tuyères, FSS, DDES, stabilité globale, champ moyen shock-wave, nozzles, FSS, DDES, global stability analysis, meanflow

Les interactions onde de choc/couche limite à l'intérieur de tuyères sur-détendues produisent des grandes régions détachées entraînant des forces non-axisymétriques appelées charges latérales. Le mécanisme générateur (mettant en jeu le décollement, les couches de cisaillement et les disques de Mach) est auto-entretenu avec des échelles de temps et d'espace relativement bien séparées de celles de la turbulence. Dans de tels cas, une approche partiellement résolue pour la turbulence (e.g. DDES) est envisageable afin de réduire les coûts de calcul. Cette thèse vise à étudier la nature potentiellement globalement instable de cette instationarité à l'aide de simulations numériques hybrides dans le cas d'une tuyère sur-détendue et par une analyse de stabilité linéaire 3D. La géométrie considérée est une tuyère TIC, en régime free shock seperation (FSS) et fonctionnant à 3 differents nombres de Mach de jet M j = [1.83, 2.09, 2.27]. Les calculs non linéaires confirment les résultats expérimentaux: la densité spectral d'énergie des perturbations de pression proche paroi est répartie sur 2 pics à fréquence intermédiaire (St = 0.2 -0.3) et deux bosses à basse et haute fréquence, respectivement à St < 1 et St ≈ 1. Pour un Mach de jet M j = 1.83, l'énergie aux fréquences St = 0.2 et St = 0.3 sont comparables. A M j = 2.09, le pic d'énergie à St = 0.2 domine tandis que les deux finissent par disparaître pour un Mach de jet M j = 2.27. Une PSD calculée pour différentes composantes azimutales de la perturbation de la pression proche paroi montre une claire séparation azimutale pour toutes les contributions frequentielles. En particulier, alors que le pic à St = 0.3 a une double contribution m = 2-3 à M j = 1.83 et il est exclusivement de type m = 2 à M j = 2.09. Le pic à St = 0.2 a, quant à lui, une symétrie m = 1, avec une signature persistante à l'intérieur de la tuyére à M j = 2.09. Dans un second temps, une analyse de stabilité globale est effectuée autour du champ moyen DDES à M j = 2.09. Une telle analyse donne un mode instable à St = 0.2 avec une symetrie azimutale de type m = 1. Ce mode se développe à partir du point de décollement et se localise au niveau de la couche de cisaillement externe.

List of Tables

 [START_REF] Jacquin | Experimental study of shock oscillation over a transonic supercritical profile[END_REF]) is also reported. . . . . . . . . N P R = 12 (bottom). On the left, contours of static pressure in terms of ratio with static ambient pressure p = p * /p * a . The colorscale is centered at p = 1 (i.e. where pressure p * matches ambient pressure p a ) to highlight the shock system. x s is the separation point, while x d = (x d1 , x d2 , x d3 ) the locations of the 1st, 2nd and 3rd Mach Disk. All these characteristic positions are listed in Table 5.2. On the right, contours of streamwise velocity u. The solid black line is the sonic line. . . . . . . . . . . . . . . 5.5 Comparison between experiments and RANS calculations in terms of wall pressure and separation point location along the same generatrix for different NPRs. p 0 and L are the combustion chamber total pressure and the divergent length, respectively. . . . . .

List of Figures

RANS solution, N P R = 9. Contours of turbulent to molecular viscosity ratio (µ T /µ).

Turbulent viscosity is too much high to permit the shear layers destabilisation. . . . . The iso-surfaces of the (a) streamwise, (b) radial and (a) azimuthal velocity mode component are plotted for the levels ûx = ±1.0 × 10 -6 , ûr = ±1.0 × 10 -7 and ûθ = ±1.0 × 10 -7 , respectively. White and black iso-surfaces for positive and negative values, respectively. The azimuthal m = 1 feature is well visible in the shape of the mode. . . 5.24 LURANS: three-dimensional view of the stable global mode at (σ, St) = (-0.077, 0.298).

The iso-surfaces of the (a) streamwise, (b) radial and (a) azimuthal velocity mode component are plotted for the levels ûx = ±1.0 × 10 -6 , ûr = ±1.0 × 10 -7 and ûθ = ±1.0 × 10 -7 , respectively. White and black iso-surfaces for positive and negative values, respectively. The azimuthal m = 2 feature is well visible in the shape of the mode. . .

Introduction

Flow separation is an essential feature in the dynamics of most aerodynamic configurations in aeronautics or space. In rocket engines, thrust is produced by the expulsion of an exhaust fluid that has been accelerated to supersonic conditions through a convergent-divergent nozzle. This acceleration goes along with pressure reduction and the latter could be so strong that the flow separates and an oblique shock occurs inside the nozzle, yielding performance losses. In this case, nozzles are called over-expanded. In order to get the optimum performance over most of the flight trajectory, nozzles are designed to avoid any shock occurrence before exit. Here, where the low flow pressure is forced to rapidly equilibrate ambient one, shock is unavoidable and it is followed by a series of other shocks and expansion fans forming the so-called "shock diamond" structures shown in Figure 1.

Organisation of the manuscript

In Chapter 1, fundamental notions about supersonic nozzles are provided. Secondly, a detailed analysis of the state of art concerning shock wave boundary layer interactions (SWBLI) and their connection with the main unsteady phenomena detected in over-expanded nozzles are described. Finally, phenomenology as well main results about linear stability analysis for laminar and turbulent flows are reported. Chapter 2 is dedicated to the description of nonlinear and linear governing equations, as well fundamental notions about stability analysis. In Chapter 3, description of the numerical strategies adopted in order to solve nonlinear, linear and stability problems is provided. Governing

Supersonic nozzles

The maximum possible thrust provided by a nozzle depending on its contour and the specific operating conditions, efficient design of nozzles has become a critical aspect of modern launchers. The main operating parameter is the ratio between the combustion chamber total pressure and the ambient static pressure (p 0 /p a ), commonly known as Nozzle Pressure Ratio (NPR). While there is an upper limit to the velocity of a convergent nozzle by the optimum adjustment of NPR, and that limit is the speed of sound, it has been proven by the Swedish engineer De Laval that it is possible to further accelerate the flow up to supersonic conditions by adding to the convergent duct a suitably designed divergent section. For a given exit section, there is however only a single value of NPR at which the 1.1. SUPERSONIC NOZZLES exhaust flow can be isentropically expanded up to the static ambient pressure (p e = p a ) † still being supersonic. This value is called design NPR (NPR d ) and it is related to the corresponding exit design Mach number (M d ) as

NPR d = 1 + γ -1 2 M 2 d γ (γ-1)
.

(1.1)

At this value of NPR, the flow is said to be adapted and the maximum exit velocity is reached. However, the ambient static pressure (resp. chamber total pressure) being greater (resp. lower) than the design ambient static pressure (resp. design chamber total pressure), the flow is abruptly compressed by an inner and/or outer shock wave up to the ambient pressure. In this case, nozzles are over-expanded and the resulting flow is subsonic or slightly supersonic. On the other hand, if NPR is greater than NPR d , there are outer expansions and nozzles are called under-expanded (see Figure 1.1). At values of NPR different from design, one can use the design Mach number as reference. In this perspective, people generally prefer using the jet Mach number (M j ) rather than the NPR, which is defined as

M j = NPR (γ-1)/γ -1 2 (γ -1) 1/2 .
(1.2) † pe is the nozzle exit static pressure

GEOMETRY DESIGN FOR ROCKET NOZZLES

In this framework, over-expanded nozzles are defined for M j < M d , while for under-expanded ones

M j > M d .
Although nozzles are designed for adapted conditions for most of their trajectory, at ground level As previously mentioned, for fixed nozzle geometry and combustion chamber conditions, there is only one ambient pressure value for which a flow isentropically and supersonically expands until the exit section. This pressure is called adaptation pressure and corresponds to the design condition for the ambient pressure (p a,d ). Conversely, the lowest ambient pressure value for which flow stays isentropically subsonic until exit section, even if sonic at throat, is called limit pressure (p l ). NPR corresponding to this condition is the minimum one guaranteeing a sonic condition at the throat. For NPR values greater than p 0 /p l and lower than p 0 /p usu (i.e. p usu /p 0 < p/p 0 < p l /p 0 ) the occurrence of a normal shock inside the divergent duct leads to a downstream subsonic solution. At the ambient pressure value of p usu the normal shock occurs just before the nozzle exit. For NPR values greater than p 0 /p usu but lower than p 0 /p a,d (i.e. p a,d /p 0 < p/p 0 < p usu /p 0 ), the solution would stay supersonic after the passing of an oblique shock. Shock-wave boundary layer interactions occurring in actual over-expanded nozzles however lead to flow separation and much more complex configurations, strictly related to NPR and nozzle contour. This scenario will be highlighted in the next chapter.

(p 0 = p 0,

Geometry design for rocket nozzles

In supersonic nozzles design, the shape of divergent must be chosen very carefully. In fact, if the nozzle contour is not proper, although optimal operating conditions normally ensure the correct behaviour of the system, shock waves may occur inside the duct. For this reason, several nozzle contours have been developed. ‡ p 0,d and p a,d are combustion chamber total pressure and ambient static pressure at design conditions, i.e. adapted conditions, respectively.

Truncated Ideal Contour (TIC) nozzles

The most common technique providing a shock-free and isentropic flow, taking into account the three dimensional character of the problem, is the Method of Characteristics (MOC). As one can remark by observing Figure 1.3a, a standard Ideal Contour geometry made by means of MOC however demands the nozzle to be very long and therefore not suitable for space applications.

Considering that thrust contribution in the last part of the nozzle is very low due to the small wall slope, a practical solution is truncating the contour, i.e. using a truncated ideal contoured nozzle (TIC) [START_REF] Frey | The tictop nozzle: a new nozzle contouring concept[END_REF].

Thrust Optimised Contour and Thrust Optimised Parabolic nozzles

The Method of Characteristics is based on the Riemann invariants as compatibility equations in order to solve the Euler equations inside the nozzle. They do not however provide an optimal nozzle contour in terms of thrust. Using a variational optimisation approach, Rao (1958) [START_REF] Rao | Exhaust nozzle contour for optimum thrust[END_REF] however showed that it is possible to obtain another set of compatibility equations providing a nozzle contour capable to produce the maximum thrust at a fixed exit Mach number. Nozzles designed with a such method are called Thrust Optimised Contour (TOC) nozzles. However, the shock pattern inside the nozzle drastically changes. In contrast with a TIC nozzle, a TOC nozzle produces an internal shock at the (a) The solid black line with red filled circles represents the wall points coming from MOC method. The red and blue net is the caractheristic one (i.e. kernel).

x/y t throat where the contour curvature changes more rapidly as shown in Figure 1.4b.

Later, Rao (1960) [START_REF] Rao | Approximation of optimum thrust nozzle contours[END_REF] proved that the contour designed with his variational method could be approximated with a skewed parabola, without introducing significant performance loss. In this case, nozzles are called Thrust Optimised parabolic (TOP) nozzles. As a TOC nozzle, a TOP produces an From Nasuti & Onofri (2009) [START_REF] Nasuti | Shock structure in separated nozzle flows[END_REF] internal shock originating at the discontinuity where the contour is approximated by a parabolic curve [START_REF] Frey | The tictop nozzle: a new nozzle contouring concept[END_REF].

Shock wave-boundary layer interactions in Supersonic nozzles

When flow travels towards less pressured regions, it could achive the null velocity gradient condition near wall and separates. In the divergent duct of supersonic nozzles, the flow accelerates against an adverse pressure gradient due to the supersonic expansion and causes the boundary-layer thickening.

In over-expanded nozzles, compression waves can originate from the deflected mean flow and eventually coalesce into an oblique shock inducing the flow separation, as shown in Figure 1.5.

Depending on nozzle contour, different shock patterns may appear in over-expanded jet flows (Figure 1.6):

• Mach reflection, where the incident shock reflects on the axis of the nozzle at high NPRs. It can be observed in the plume of all nozzle types but only during flight;

• Mach disk, with the occurring of an incident shock reflecting on a Mach disk. The meeting point among the incident shock, the reflected shock and the Mach disk is called Triple Point. The Mach disk is characterised by a strong shock with a downstream subsonic region;

• Cap-shock pattern: when the nozzle is optimised for maximum thrust (TOC, TOP), an internal shock originates at throat. This shock hits the Mach disk forming complex shock patterns inside the nozzle. In this case a subsonic zone with large recirculation bubbles, called trapped vortex, is observed downstream of the Mach disk. Moreover, pressure fluctuations have been observed to propagate in the opposite direction of the local mean velocity, leading to the hypothesis that an acoustic resonance exists in the separated region.

This phenomenon could be the source of the low-frequency displacement of separation point [START_REF] Kistler | Fluctuating wall pressure under a separated supersonic flow[END_REF].

SHOCK WAVE-BOUNDARY LAYER INTERACTIONS IN SUPERSONIC NOZZLES

Free Shock and Restricted Shock Separation in Over-expanded Nozzles

The complex separation mechanism, inducing self-sustained oscillations, has soon been suspected to be related to the asymmetric and unsteady forces (or side-loads) that may damage nozzles during nozzles first experience a decrease followed by an increase of the pressure separation with increasing M j (see Figure 1.10b). On the contrary, in higher M d (M d > 1.8) nozzles, pressure separation decreases all the way to the exit with increasing M j (see Figure 1.10c). Particularly, the dimensionless minimum wall pressure position is shown to scale with the ratio M j /M d like an hyperbolic tangent [START_REF] Zaman | Shock-induced boundary-layer separation in round convergent-divergent nozzles[END_REF].

The discovery made by Nave & Coffey increased the interest about such self-sustained unsteadiness in nozzle flows so that many experiments were conducted on divergent diffusers flows for almost two decades.

Self-sustained shock oscillation in divergent diffusers

The dynamics complexity as well as the limited numerical resources of time first led to the launch of several experimental campaigns on two-dimensional diffusers. One of the first experimental investigations aiming at the characterisation of self-induced fluctuations present in nozzle flows with reattaching flow separation have been conducted by [START_REF] Sajben | A high-speed schlieren investigation of diffuser flows with dynamic distortion[END_REF] [START_REF] Sajben | A high-speed schlieren investigation of diffuser flows with dynamic distortion[END_REF]. In this campaign, several two-dimensional small-angled divergent diffusers, with area ratio ranging from 1.84 to 2.5, at transonic flow conditions have been studied. In this scenario, time-mean flow properties are shown to be correlated to the shock-wave Mach number (M sw ) which, based on the spatial minimum static wall pressure p w min , reads

SHOCK WAVE-BOUNDARY LAYER INTERACTIONS IN SUPERSONIC NOZZLES

M sw =      (p w min p 0 -γ-1 γ -1   2 γ -1    1/2 (1.3)
and represents a measure of time-mean conditions just upstream of the shock. Particularly, it was found that for M sw < 1.3 flow separation is due to the adverse pressure gradient in the subsonic flow and occurs downstream of the shock, while for M sw > 1.3 the pressure jump in the flow, caused by a stronger lambda-pattern shock, is sufficiently large to immediately separate as shown in Figure 1.11.

For this reason the first separation is called Pressured-Gradient-Induced Separation (PGIS), while the second one Shock-Induced Separation (SIS). In the latter case, highspeed shadowgraphs have shown a complex dynamic behaviour driven by turbulent structures, upstream propagating shocklets and shock oscillations. the presence of a low-frequency unsteadiness with or without outer forcings [START_REF] Bogar | Structure of self-excited oscillations in transonic diffuser flows[END_REF]. [START_REF] Hsieh | Numerical simulation and comparison with experiment for self-excited oscillations in a diffuser flow[END_REF] [26] numerically reproduced the same experimental test and observed the expected self-sustained shock oscillation. Similarities in terms of both unsteady and time-average flow properties with the supersonic flow over forward facing step arose and with them the idea that the concerned unsteadiness could be due to a similar phenomenon of acoustic resonance. In fact, it would not be the first time that acoustics play a crucial role in supersonic nozzles, given that acoustic tones had already been detected in under-expanded nozzles.

Screech tones in under-expanded nozzles

In high NPR over-expanded nozzles and under-expanded nozzles, when the jet Mach number (M j ) far exceeds unity, a feedback close-loop mechanism known as screech settles. The latter is governed by three main features of the jet exhausted from the nozzle: the stationary repetitive or "cellular" pattern, called shock cells, arising from the train of oblique shock waves and expansion fans; the sinuous oscillation of the jet and the acoustic waves emanating from the turbulent flow (see Figure 1.12). This phenomenon, first described by Powell (1958) [START_REF] Powell | On the noise emanating from a two-dimensional jet above the critical pressure[END_REF] and then by Raman (1998) [START_REF] Raman | Advances in understanding supersonic jet screech: review and perspective[END_REF] is composed of different phases. During the first phase, that is inside the jet, there is the birth of turbulent structures in the mixing layer. Subsequently, these structures, conveyed in the mixing layer, pass through the shock cells. Interactions then take place between these structures and the network of shock cells and acoustic waves are created. These waves spread in all directions, especially upstream, where they then excite mixing layer at the outlet of the nozzle. This excitation is at the origin of disturbances in the mixing layer which will in turn produce turbulent structures.

Tam (1986) [START_REF] Tam | Proposed relationship between boroadband shock associated noise and screech tones[END_REF] realised that screech frequencies (f ) collapse on the same curve inversely decreasing with the fully expanded jet Mach number (M j ), when the fully expanded exit diameter (D j ) and velocity (U j ) are used as the length and the velocity scale of the Strouhal number (St) respectively, as

U j = γRT 0 M j 1 + γ -1 2 M 2 j , D j = D e    1 + γ -1 2 M 2 j 1 + γ -1 2 M 2 d    γ+1 4(γ-1) , St = f D j U j ≈ 0.67 M 2 j -1 -1 2 1 + 0.7M j 1 + γ -1 2 M 2 j -1 2 T A T 0 -1 2 -1 , 1.

SHOCK WAVE-BOUNDARY LAYER INTERACTIONS IN SUPERSONIC NOZZLES

Figure 1.12: Resonant screech loop mechanism, by Raman (1998) [START_REF] Raman | Advances in understanding supersonic jet screech: review and perspective[END_REF].

with D e the nozzle exit diameter, T 0 the combustion chamber temperature and T A the ambient jet (1986) [START_REF] Tam | Proposed relationship between boroadband shock associated noise and screech tones[END_REF] for screech frequency.

SHOCK WAVE-BOUNDARY LAYER INTERACTIONS IN SUPERSONIC NOZZLES

Transonic resonance in divergent diffusers and conical nozzles

Even though self-sustained shock oscillations in small divergent diffusers were suspected to be due to a standing wave produced by the interaction of an upstream and a downstream propagating travelling wave [START_REF] Hsieh | Numerical simulation and comparison with experiment for self-excited oscillations in a diffuser flow[END_REF], there were no solid theory satisfactorily proving the experiments evidences about that. For this reason, Zaman et al. (2002) [START_REF] Zaman | Investigation of a'transonic resonance'with convergent-divergent nozzles[END_REF] launched an experimental campaign on several small convergent-divergent nozzles and diffusers with different geometries and operating conditions. The unsteadiness detected, called transonic resonance, have been proven to be due to a phenomenon different from screech. While the former takes place in two or three stages at low M j with a frequency (f tr ) increasing almost linearly with M j , the latter takes place in only one stage at high M j with a frequency (f sc ) decreasing linearly with M j , as shown in Figure 1.14. As for Zaman, the inner shock and the external fluctuating pressure behave like a dual oscillator whose fundamental resonance frequency (f tr ) results in a good agreement with experiments if compared with a one-quarter acoustic standing wave as

f tr = a 1 -M 2 e 4L (1.4)
where M e is the Mach number in the separated region at exit, a the local speed of sound and L is approximated with the distance between the shock foot and the nozzle exit. However, this model is less and less predictive for higher-NPR nozzle flows. The reason for that, according to Wong (2005) [START_REF] Wong | Theoretical prediction of resonance in nozzle flows[END_REF], is that a standing-wave resonance can occur only in low supersonic or transonic flows while it tends to be attenuated in a high supersonic ones. By a perturbative quasi-one-dimensional model coupled with a dual-oscillator concept, the damping effect is shown to be due to the interaction and correlation between compressibility and velocity in the subsonic domain downstream of the shock [START_REF] Wong | Theoretical prediction of resonance in nozzle flows[END_REF].

Unsteadiness dynamics in modern over-expanded nozzles

As technology advances, many experimental studies have been conducted on supersonic or hypersonic cold subscaled rocket nozzle flows. As for unsteady dynamics, Torgren (2002) [START_REF] Torngren | Correlation between outer flow and internal nozzle pressure fluctuations[END_REF] side loads are described: the former due to random pressure fluctuations and the latter to FSS/RSS transition. In the first case, flow inside the separated region is intermittent as for the case of supersonic flow over forward-facing step: pressure perturbation can be correctly modelled with a low and high-frequency part, correlated with shock displacement and shear layer-boundary layer interaction respectively. The first attempt to find a modal separation in over-expanded nozzle dynamics was made by Baars et al. (2012) [START_REF] Baars | Wall pressure unsteadiness and side loads in overexpanded rocket nozzles[END_REF]. In this work, a TOP nozzle experiencing a FSS regime reveals in the energy spectrum of wall pressure fluctuations two distinct peaks without any outer forcing. The first is related to a low-frequency unsteadiness (0.2 kHz) due to the shock displacement, while the second covers a large band of high frequencies (1.0-10.0 kHz) and its energy grows with increasing distance downstream. By computing a PSD for each azimuthal mode, the entire spectrum appears to be distributed over the first two Fourier modes. Particularly, the 0 th mode energy decreases along the streamwise direction, while the 1 st mode energy (the only mode responsible for side loads) remains relatively constant [START_REF] Baars | Wall pressure unsteadiness and side loads in overexpanded rocket nozzles[END_REF]. The experiments carried out by one with a recirculation bubble inside the separated region at nozzle lip and another with a large trapped vortex behind the Mach disk. URANS methods turned out to be very efficient in capturing unsteadiness in many compressible flows such as airfoil buffetting [START_REF] Brunet | Computational study of buffet phenomenon with unsteady rans equations[END_REF][START_REF] Deck | Numerical simulation of transonic buffet over a supercritical airfoil[END_REF][START_REF] Thiery | Numerical prediction of shock induced oscillations over a 2d airfoil: Influence of turbulence modelling and test section walls[END_REF] as well in evaluating low-frequency oscillations and the resulting side loads found in subscale rocket models [START_REF] Deck | Numerical simulation of side loads in an ideal truncated nozzle[END_REF][START_REF] Deck | Unsteady side loads in a thrust-optimized contour nozzle at hysteresis regime[END_REF][START_REF] Moríñigo | Three-dimensional simulation of the self-oscillating flow and sideloads in an over-expanded subscale rocket nozzle[END_REF].

However, they fail when wall turbulence and self-sustained oscillations get close in term of frequencies (less than one order of magnitude). On the other hand, using a Direct Numerical Simulation (DNS) or (1997) [START_REF] Spalart | Comments on the feasibility of les for wings, and on a hybrid rans/les approach[END_REF]. The interest for the latter methodology, permitting a computational resource saving for the RANS-treated boundary layer resolution, has highly increased in recent years. Deck (2009) [START_REF] Deck | Delayed detached eddy simulation of the end-effect regime and side-loads in an overexpanded nozzle flow[END_REF] and [START_REF] Martelli | Flow dynamics and wall-pressure signatures in a high-reynolds-number overexpanded nozzle with free shock separation[END_REF] performed a DDES calculation of an over-expanded TIC nozzle experiencing a FSS regime. Simulations, validated with experimental measurements, confirmed for a prescribed NPR the existence of a low-frequency breathing mode as well an higher frequency secondary contribution. The azimuthal PSD analysis revealed a m = 0 azimuthal mode related to shock oscillation at low-frequency as well a m = 1 azimuthal mode responsible for side-loads associated to the higher frequency. The latter is compatible with screech correlation proposed by Tam et al. [START_REF] Tam | Proposed relationship between boroadband shock associated noise and screech tones[END_REF] as well with the findings by Jaunet et al. [START_REF] Jaunet | Wall pressure and external velocity field relation in over-expanded supersonic jets[END_REF]. A feedback-loop model involving the turbulent shear layers, triple point and second

SHOCK WAVE-BOUNDARY LAYER INTERACTIONS IN SUPERSONIC NOZZLES

Mach disk is proposed to explain the side-loads mode.

The insight that nozzle low-frequency shock oscillations were due to a standing wave produced by an upstream and a downstream propagating travelling wave [START_REF] Hsieh | Numerical simulation and comparison with experiment for self-excited oscillations in a diffuser flow[END_REF] is nowadays the commonly accepted theory [START_REF] Zaman | Investigation of a'transonic resonance'with convergent-divergent nozzles[END_REF][START_REF] Wong | Theoretical prediction of resonance in nozzle flows[END_REF][START_REF] Martelli | Flow dynamics and wall-pressure signatures in a high-reynolds-number overexpanded nozzle with free shock separation[END_REF] and known as transonic resonance. On the contrary, the nature of high frequency unsteadiness is not yet clear and it is often associated to screech [START_REF] Powell | On the noise emanating from a two-dimensional jet above the critical pressure[END_REF][START_REF] Tam | Proposed relationship between boroadband shock associated noise and screech tones[END_REF][START_REF] Raman | Advances in understanding supersonic jet screech: review and perspective[END_REF][START_REF] Jaunet | Wall pressure and external velocity field relation in over-expanded supersonic jets[END_REF][START_REF] Martelli | Flow dynamics and wall-pressure signatures in a high-reynolds-number overexpanded nozzle with free shock separation[END_REF].

Another question raised is whether the flow dynamics in a such unsteadiness has a linear or nonlinear origin. In the former case, analysis tools like Dynamic Mode Decomposition (DMD) [START_REF] Schmidt | Guide to spectral proper orthogonal decomposition[END_REF] or Global Stability Analysis could be used to extract information about its nature and source. The mathematical formulation of the flow dynamics is based on the well known Navier-Stokes equations. Current computational resources do not however permit to solve such partial differential nonlinear equations for highly turbulent flow regimes. The objective of this chapter is to provide an overview of the different strategies adopted in such cases. Particularly, a full description of RANS/LES methods will be provided by detailing the strategy adopted in this work of thesis. Finally, a presentation of the linear stability framework will be given.

Nonlinear Compressible Navier-Stokes equations

The instantaneous continuity, momentum and energy equation for a three-dimensional, compressible, unsteady flow in Cartesian coordinates can be written as:

2.1. NONLINEAR COMPRESSIBLE NAVIER-STOKES EQUATIONS                            ∂ρ * ∂t * + ∂ ∂x * j ρ * u * j = 0 ∂ (ρ * u * i ) ∂t * + ∂ ∂x * j ρ * u * j u * i + p * δ ij -τ * ij = 0 ∂ (ρ * E * ) ∂t * + ∂ ∂x * j ρ * E * u * j + p * u j + q * j -u * i τ * ij = 0 (2.1)
where the superscript * denotes dimensional quantities and Einstein summation for repeated indices is used. The flow quantities of interest are the density ρ * , the streamwise, wall-normal and spanwise velocity components u * i in the cartesian coordinate directions x * i (i = 1, 2, 3), the static pressure p * and the total specific energy E * . The latter is obtained by adding the flow specific internal energy e * to the specific kinetic one as

E * = e * + u * i u * i 2 . ( 2.2)
The viscous stresses due to the friction between the fluid and the surface of an element are provided by the viscous stress tensor τ * ij . It describes the stress component affecting a plane perpendicular to the i-axis, along the direction of the j-axis. For most practical problems, the fluid is assumed to be Newtonian and the components of the viscous stress tensor are defined by

τ * ij = µ * L ∂u * j ∂x * i + ∂u * i ∂x * j - 2 3 ∂u * k ∂x * k δ ij , (2.3)
where µ * L denotes molecular dynamic viscosity and represents the resistance for unit area of a fluid to a unit deformation by shear stress.

The heat transfer due to temperature gradient is provided by q * j and generally reads a form of Fourier's law of heat conduction, i.e.,

q * j = -κ * ∂T * ∂x j (2.4)
with κ * standing for the thermal conductivity coefficient and T * for the absolute static temperature.

For a three-dimensional problem the above equations result in a system of five equations for the five conservative variables ρ * , ρ * u * , ρ * v * , ρ * w * , and ρ * E * with seven unknown primitive variables,

NONLINEAR COMPRESSIBLE NAVIER-STOKES EQUATIONS

namely: ρ * , u * , v * , w * , E * , p * and T * . Therefore, to close the system we need two more additional equations, which have to be thermo-dynamic relations between the state variables, e.g. a function linking pressure to density and temperature, and another one expressing internal energy in terms of pressure and temperature.

Assuming that the working fluid is a thermally perfect gas, these closing are respectively the thermal and caloric equations of state given respectively by

p * = ρ * R * T * (2.5a) e * = c * v T * , ( 2.5b)
where R * = R * /M * is the perfect gas constant divided by the molar mass of the fluid, commonly set at 287.1 J kg -1 K -1 for air, c * v the specific heat capacity at constant volume linked to specific heat capacity at constant pressure c * p by the specific heat ratio γ =

c * p c * v , set at 1.4 for diatomic gas † .
Molecular dynamic viscosity µ * L , for a perfect gas, strongly depends on temperature but only weakly on pressure. For this reason fluid is often imposed to obey Sutherland's law, i.e.

µ * L = µ * L,ref T * T * ref 3/2 T * ref + T * S T * + T S , ( 2.6)
where µ * L,ref = 1.716 × 10 -5 N m -2 s is the dynamic viscosity at the reference temperature T * ref = 273.15 K reference temperature and T * S = 110.4 K is the Sutherland temperature. In the case of gases, the thermal conductivity coefficient κ * depends on temperature in a similar way.

For this reason, the relationship

κ * = c * p µ * L P r (2.7)
is usually accepted for air, for which the Prandtl number P r is assumed to be constant and set at the value P r = 0.72.

Dimensionless Navier-Stokes equations

In the present work, the dimensionless compressible Navier-Stokes equations are solved. They are obtained by expressing each dimensional variable as the product of a corresponding dimensionless

NONLINEAR COMPRESSIBLE NAVIER-STOKES EQUATIONS

quantity and a reference one, and injecting it in Eqs.(2.1). However the reference quantities are not those associated to free stream conditions (U ∞ , ρ ∞ , T ∞ ) as for external aerodynamics problems, but the total ones (a * 0,∞ , ρ * 0,∞ , T * 0,∞ ). For the sake of simplicity, we will refer to them omitting the symbol * . By using the dimensional group (a ∞ 0 , ρ ∞ 0 , T ∞ 0 , L ref ), with L ref the characteristic length of the problem, the dimensionless variables are defined as

ρ = ρ * ρ ∞ 0 , u i = u * i a ∞ 0 , T = T * T ∞ 0 , x i = x * i L ref , ( 2.8
)

µ L = µ * L µ ∞ L,0 Re ∞ 0 = 1 Re ∞ 0 (T ) 3 2 1 + T * s /T ∞ 0 T + T s /T ∞ 0 , (2.9) p = p * ρ ∞ 0 (a ∞ 0 ) 2 = p * p ∞ 0 γ , ( 2.10
)

τ ij = τ * ij ρ ∞ 0 (a ∞ 0 ) 2 = τ * ij p ∞ 0 γ , ( 2.11
)

q j = -c p µ L P r ∂T ∂x j , ( 2.12
)

t = t * a ∞ 0 L ref , ( 2.13)
with

p = ρT γ , c p = 1 γ -1 , c v = c p γ , e = c v T, a ∞ 0 = γRT ∞ 0 . (2.14)
The Reynolds number is based on the dimensional total quantities and is therefore defined as

Re ∞ 0 = ρ ∞ 0 a ∞ 0 L ref µ ∞ 0 , (2.15)
that, in case of external aerodynamics, is related to the free stream Reynolds number (Re ∞ ) and Mach number (M ∞ ) by the following expression

Re ∞ 0 = p ∞ 0 √ γ RT ∞ 0 L ref µ ∞ 0 = Re ∞ M ∞ T ∞ 0 T * ∞ 1 γ-1 1 + T * s /T ∞ 0 1 + T * s /T * ∞ (2.16)

COMPUTATIONAL APPROACHES TO TURBULENT FLOWS

with

T ∞ 0 T * ∞ = 1 + γ -1 2 M 2 ∞ , Re ∞ = ρ ∞ U ∞ L ref µ L,∞ , M ∞ = U ∞ a ∞ .
(2.17)

Computational approaches to turbulent flows

Solving numerically the Navier-Stokes equations is possible by means of Direct Numerical Simulation (DNS), which requires that the whole range of spatial and temporal scales of the flow must be resolved. As stated first by Richardson, and then assessed by Kolgomorov, the more turbulent a flow is, the greater the number of stages of dissipating eddies which composes it. According to the theory of energy cascade, in turbulent flows the most energetic eddies, whose size is comparable to the characteristic flow scale, are unstable and break up, transferring their energy to smaller eddies. The smallest eddies, at which energy is dissipated, are of the order of the so-called Kolmogorov scales (η, u η , τ η ), whose universal features permits us to derive the following important relations relative to the largest scales of the problem (l 0 , u 0 , τ 0 ) in space, velocity and time respectively [START_REF] Pope | Turbulent flows[END_REF]

l 0 η ≈ Re 3 4 , u 0 u η ≈ Re 1 4 , τ 0 τ η ≈ Re 1 2 .
(2.18)

Eq.(2.18) tell us that simulating the evolution of the system over a single time unit would require a number of volumetric elements scaling with Re 9 4 for a number of iterations scaling with Re 1 2 , with a consequent number of operations approximately of the order of Re 3 . These conclusions lead us to the impossibility, with the current computational resources, to tackle the vast majority of practical engineering problems by means of DNS.

To this day, three different approaches have been proposed to overcome this limit: the Reynolds-Average Navier-Stokes (RANS) equations, the Large-Eddy Simulations (LES) and Detached Eddy Simulations (DES). While in a RANS approach all turbulence is modelled by solving an additional equation, the idea behind LES is the spectral separation of larger turbulent scales from the smaller scales, where the former are resolved directly and the latter approximated by a model. Over the past few years, DES has become popular thanks to a greater accuracy compared to RANS combined with a smaller computational cost than LES.

COMPUTATIONAL APPROACHES TO TURBULENT FLOWS

Nonlinear Compressible URANS equations

The first approach to tackle turbulent flows was presented by Reynolds in 1883 [START_REF] Reynolds | Xxix. an experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels[END_REF]. The methodology is based on the decomposition of the flow variables into a time-average and a fluctuating part and solving the governing equations (2.1) only for the former. The obtained equations are called RANS (Reynolds Averaging Navier-Stokes) equations. However, in the case of non-constant density, such as high-Mach compressible flows in supersonic jet nozzles, it is more convenient to employ a time averaging (Reynolds averaging) for density and pressure, while a mass-weighted averaging (Favre averaging) in addition to the time averaging for the other variables. Mass-averaging the flow variables can be done by multiplying the primitive ones by the density, averaging this product and dividing by the average density as

u = ρu ρ v = ρv ρ w = ρw ρ T = ρT ρ e = ρe ρ E = ρE ρ = e + u i u i 2 + u i u i 2 K (2.19)
After that, assuming the perfect gas hypothesis

p = ρR T e = c v T γ = c p c v R = c p -c v (2.20)
and defining the following fluctuating quantities

u = u + u v = v + v w = w + w ρ = ρ + ρ p = p + p T = T + T e = e + e (2.21)
for which, the following relations are met

b = 1 T T 0 b dt = 0 b = 1 T T 0 b dt = 0 ρb = 1 T T 0 ρb dt = 0 (2.
                         ∂ρ ∂t + ∂ ∂x j (ρ u j ) = 0 ∂ ∂t (ρ u i ) + ∂ ∂x j ρ u i u j + pδ ij -τ tot ij = 0 ∂ ∂t ρ E + ∂ ∂x j ρ u j E + u j p + q tot j -u i τ tot ij = 0 (2.23)
where

                                                     τ tot ij = τ lam ij + τ turb ij τ lam ij = µ L ∂ u i ∂x i + ∂ u j ∂x j - 2 3 δ ij ∂ u k ∂x k τ turb ij = -ρ u i u j + µ L ∂u i ∂x i + ∂u j ∂x j - 2 3 δ ij ∂u k ∂x k =τ ij τ lam ij ≈ -ρ u i u j q tot j = q lam j + q turb j q lam j = q j = -c p µ L P r ∂ T ∂x j q turb j = ρ u j h (2.24)
In order to close the system (2.23)-(2.24), we need to provide:

• six components of the Favre-averaged Reynolds-stress tensor τ turb ij ;

• three components of the turbulent heat-flux vector q turb j .

A common way to write the term τ turb ij is using the Boussinesq's Eddy-Viscosity hypothesis, i.e.

-ρ u i u j ≈ µ T ∂ u i ∂x i + ∂ u j ∂x j - 2 3 δ ij ∂ u k ∂x k - 2 3 ρKδ ij (2.25)
while the term q turb j is usually approximated as

ρ u j h ≈ -c p µ T P r T ∂ T ∂x j (2.26)
where µ T stands for the turbulent eddy viscosity and P r T for the turbulent Prandtl number, which is assumed to be constant over the whole flowfield and set at 0.9 for air.

COMPUTATIONAL APPROACHES TO TURBULENT FLOWS

The Boussinesq's hypothesis permits to close the system after determining µ T by means of one or twoequation turbulence model. In the former case, due to the lack of any information about K modelling, the last term in Eq.(2.25) is ignored.

Spalart Allmaras (SA) one-equation model

One of the most used turbulence models in RANS calculations is the Spalart-Allmaras model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF].

It has been derived as a simple convection-diffusion equation for the pseudo turbulent viscosity variable

ν where the source term is obtained as a balance between a production term P and a destruction term D associated to same quantity, as

∂ (ρ ν) ∂t + ∂ (ρ νu j ) ∂x j - 1 σ ∂ ∂x j (ρν + ρ ν) ∂ ν ∂x j = c b1 (1 -f t2 ) Sρ ν + c b2 σ ∂ ν ∂x j ∂ρ ν ∂x j P -ρ c w1 f w - c b1 κ 2 f t2 ν d w 2 D (2.27)
where the presence of density ρ is justified by the will to extend the above equation to compressible flows [START_REF] Allmaras | Modifications and clarifications for the implementation of the spalart-allmaras turbulence model[END_REF].

The various constants in the model are standard-defined as

c b1 = 0.1355, c b2 = 0.622, σ = 2/3, κ = 0.41, c v1 = 7.1, c w2 = 0.3, c w3 = 2, (2.28
)

c t1 = 1, c t2 = 2, c t3 = 1.1, c t4 = 0.5

COMPUTATIONAL APPROACHES TO TURBULENT FLOWS

which build the following quantities for the model

χ = ν/ν, f v1 = χ 3 χ 3 + C 3 v1 , f v2 = 1 - χ 1 + χf v1 , f t2 = c t3 e (-c t4 χ 2 ) , f w = g 1 + c 6 w3 g 6 + c 6 w3 1 6 , S = Ŝ + ν κ 2 d 2 w f v2 , Ŝ = 2 Ω ij Ω ij , Ω ij = 1 2 ∂ u i ∂x j - ∂ u j ∂x i , ( 2.29
)

c w1 = c b1 /κ 2 + (1 + c b2 )/σ, r = min ν Sκ 2 d 2 w , 1 and g = r + c w2 r 6 -r .
The closure problem depends on turbulent viscosity µ T , which is computed by µ T = f v1 ρ ν.

To improve the robustness and the near-wall numerical behaviour of the model, which is modelled by the functions f v1 , f v2 , f w e r, the Edwards Modification [START_REF] Edwards | Comparison of eddy viscosity-transport turbulence models for three-dimensional, shock-separated flowfields[END_REF] has been developed. This new version resembles the "standard" version (SA), except for the lack of the trip term f t2 and the following two modified variables

S = W 1 2 1 χ + f v1 , r = tanh ν Sκ 2 d 2 w tanh(1.0) , (2.30) W = ∂ u i ∂x j + ∂ u j ∂x i ∂ u i ∂x j - 2 3 ∂ u k ∂x k 2 or 4 3 ∂ u i ∂x i 2 + ∂ u i ∂x j + ∂ u j ∂x i 2 i =j (2.31)
The turbulence model chosen for all calculations in this work is the Spalart-Allmaras One-Equation Model with Edwards Modification.

COMPUTATIONAL APPROACHES TO TURBULENT FLOWS

Fundamentals of LES approach

As mentioned in §2.2, in LES a low-pass filtering is applied to the flowfield variables so that the filtered field could be computed on a coarser grid than that required by a DNS calculation [START_REF] Pope | Turbulent flows[END_REF]. The filtering operation, commonly applied in wavenumber space, leads to a system that needs to be closed.

One or more equations are thus required to model the so-called residual stress tensor τ r ij concerning all the scales above the filter cut-off frequency. These equations form a Sub-Grid-Scale (SGS) model.

The simplest SGS model, first proposed by Smagorinsky et al. (1965) [START_REF] Smagorinksy | Numerical results from froma nine-level general circulation model of the atmosphere[END_REF] and used by [START_REF] Deardoff | A numerical study of three-dimensional turbulent channel flow at large reynolds numbers[END_REF] [54], builds the residual stress tensor through a linear eddy viscosity model and a mixing length hypothesis as

τ r ij = -2ν r Ω ij , ν r = (C s ∆) 2 Ŝ (2.32)
where C s is a constant to be calibrated and ∆ the grid spacing. Therefore, any LES calculation made on a sufficiently refined grid would tend to the solution provided by a DNS, in contrast with the RANS methodology, whose best converged solution would be at most as good as the turbulence model [START_REF] Spalart | Comments on the feasibility of les for wings, and on a hybrid rans/les approach[END_REF]. However, although more accurate than a RANS calculation, the computational cost of a LES remains prohibitive in terms of degree of refinement and time-step requirements, especially for configurations at high Reynolds numbers and for wall-bounded flows. For this reason, the interest in hybrid RANS-LES approaches, permitting major saving, both in term of computational resource and mamory footprint, for the boundary layer resolution, has highly increased in recent years.

Hybrid RANS/LES methodologies

Detached Eddy Simulation (DES)

The first formulation of an hybrid methodology RANS/LES, called Detached Eddy Simulation

f d = 1 -tanh (8r d ) 3 , r d = ν ∂u i ∂x j ∂u i ∂x j κ 2 d 2 w , ( 2.34)
d = d w -f d • max (0, d w -C DES ∆) , ∆ = ∆ max = max (∆x i ) . ( 2
∆ ω = 3 i=1 N 2 i ∆x j ∆x k i =j =k (2.36)
where N i is the unit vector giving the projection of the cell face along the direction of vorticity ω. An easier to implement version has been proposed by Pont (2015) [START_REF] Pont | Self adaptive turbulence models for unsteady compressible flows Modèles de turbulence auto-adaptatifs pour la simulation des écoulements compressibles instationnaires[END_REF] as

∆ ω = N f i=1 ω • S i 2||ω|| (2.37)
However, an exclusive use of ∆ ω in Eq.(2.35) is not recommended since it leads to reduce the eddy viscosity and worsening the MSD effect. For this reason, Deck (2011) [START_REF] Deck | Recent improvements in the zonal detached eddy simulation (zdes) formulation[END_REF] suggests to adopt a zonal approach to shield the boundary layer by the using ∆ max as well as to ensure a rapid switch into the LES mode through ∆ ω as

∆ = ∆ max f d < f d0 ∆ ω f d > f d0 (2.38)
where f d0 is an user-defined variable in the range (0.75, 0.99) and usually set at f d0 = 0.8, after calibrating on flat-plate boundary layer calculations.

The combination of Deck 's zonal approach and the DDES by This modified version of EDDES has already been successfully adopted by [START_REF] Martelli | Flow dynamics and wall-pressure signatures in a high-reynolds-number overexpanded nozzle with free shock separation[END_REF] [START_REF] Martelli | Flow dynamics and wall-pressure signatures in a high-reynolds-number overexpanded nozzle with free shock separation[END_REF] in the case of unsteady FSS over-expanded nozzles.

COMPUTATIONAL APPROACHES TO TURBULENT FLOWS

Correction of RANS model damping term behaviour

As mentioned by [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF], standard DES could mis-interpret the LES region as a situation of "wall proximity" and excessively lower the eddy viscosity through f v1 , f v2 and f t2 functions of the SA model [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF]. To prevent the erroneous activation of near wall functions, [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF] introduced the factor Ψ, build as

Ψ 2 = min 100, 1 -c b1 c w1 κ 2 f * w (f t2 + (1 -f t2 ) f v2 ) f v1 max (10 -10 , 1 -f t2 ) (2.39)
where f * w = 0.424 is the limit value of f w in Eqs. (2.29). In this way, the corrected DDES formulation in Eq.(2.35) changes in

d = d w -f d • max (0, d w -ΨC DES ∆) (2.
f v1 = 1, f v2 = 0, f t2 = 0, f w = 1 (2.42)
However, though this solution remains well posed in the zonal frame (Chauvet, Deck) it may corrupt the boundary layer solution when used with an EDDES formulation [START_REF] Mockett | A comprehensive study of detached-eddy simulation[END_REF][START_REF] Deck | Recent improvements in the zonal detached eddy simulation (zdes) formulation[END_REF]. For this reason, it has been decided to adopt the solution provided by Eq.(2.41) in the rest of this work.

GLOBAL STABILITY FRAMEWORK

Global Stability framework

The notion of instability is related to the tendency of a physical system to depart from its equilibrium as it is perturbed. Different definitions of stability exist in the literature. In the rest of this work, we will restrict ourselves to that of linear (in)stability. An equilibrium is linearly unstable if injecting an infinitesimal perturbation causes the system to depart exponentially fast from this equilibrium.

On the other hand, if the system eventually relaxes towards this equilibrium, the latter is said to be linearly stable. In the field of fluid mechanics, such equilibria are known as base flows. Depending on the properties of this base flow, various types of stability analyses could be conducted. For instance, base flows which slowly evolves in the advection direction are traditionally studied through the prism of local stability. By relying on a parallel flow assumption, the linearized equations governing the spatio-temporal dynamics of infinitesimal perturbations can be studied in the Fourier domain thus greatly facilitating their analysis and massively reducing the computational cost of such analyses. Alternatively, if such simplifying assumptions do not hold, one needs to study the properties of this base flow via the framework of global stability analysis. The latter approach has been considered in this work. In the rest of this section, a brief overview of global stability analysis for laminar and turbulent flows is given in section 2.3.1 and 2.3.1 while the precise set of equations considered are presented in section 2.3.3.

Laminar flows

Depending on the properties of the investigated flow (e.g. whether it is laminar or turbulent), different modelling strategies could be employed. For low Reynolds numbers flows, one classically models their dynamics using the standard (in-)compressible Navier-Stokes equations. In this framework, the canonical example of a globally unstable flow is the two-dimensional flow past a circular cylinder wherein, for Re > 48, the flow exhibits the well-known Bénard-von Kàrmàn vortex street (see Despite these remarkable agreements for low Reynolds number flows, as one moves away from the stability threshold, the qualitative and quantitative agreements between the predictions of stability analyses and the experimental or numerical observations deteriorate rapidly. This increasing discrepancy is classically attributed to the fact that, as the Reynolds number increases, the base flow and the mean flow rapidly differ. In this context, [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF] [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF] has shown in his seminal work that, for the cylinder flow, a stability analysis conducted around the mean flow (rather than the base flow) yields better agreements even relatively far from the stability threshold. In particular, such analysis predicts the existence of a marginally stable mode whose eigenfrequency and flow features are in good agreements with experimental and numerical observations. Although a proper theoretical justification is still lacking, similar results have been obtained by Sansica et al. (2018) [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF] for the three-dimensional supersonic flow past a sphere. It needs to be emphasized however that, by definition, such analysis requires a priori knowledge of the mean flow, thus limiting its applicability to situations wherein this mean flow can either be computed from long time average obtained from direct numerical simulations

GLOBAL STABILITY FRAMEWORK

or high-resolution experimental data.

It needs to be noted that the mean flow is not a solution to the Navier-Stokes equations. As such, the validity of stability analysis around the mean flow lies on the well-posedness of neglecting the influence of the Reynolds stresses in the perturbation equation. As to circumvent this issue, Mantic-Lugo et al.

(2014) [START_REF] Mantič-Lugo | Self-consistent mean flow description of the nonlinear saturation of the vortex shedding in the cylinder wake[END_REF] proposed an extension of the linear stability framework known as a self-consistent model.

In this new framework, the coupling between the "instantaneous" mean flow and the perturbations is accounted for by explicitly taking into account the associated Reynolds stresses. Note however that, for the sake of tractability, it is assumed that this perturbation can be well approximated by consid-

ering

Turbulent flows

For high-Reynolds number flows which are typically turbulent, the concept of base flow is strongly dependent on the modelling strategy adopted. One could use for instance an approach wherein turbulence is entirely modelled such as the RANS/URANS formulation. In this framework, the fixed point of the URANS equations is nothing but the solution of the RANS ones. As such, the concept of base flow is well posed and one can study whether this solution is linearly stable or not. Such an approach has first been proposed by [START_REF] Crouch | Predicting the onset of flow unsteadiness based on global instability[END_REF] [START_REF] Crouch | Predicting the onset of flow unsteadiness based on global instability[END_REF]. In their work, the authors considered a 2.3. GLOBAL STABILITY FRAMEWORK two-dimensional turbulent transonic flow past a NACA0012. After having computed the base flow (solution to the RANS equations), both the URANS equations and turbulence model have been linearized

in the vicinity of this particular solution and its stability investigated using classical tools from global stability analysis. This analysis reported the existence of an unstable mode whose behaviour was found to be in good agreement with results from the literature and observations from nonlinear calculations.

Shortly after, Sartor et al. (2015) [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF] used the same framework to conduct stability, receptivity and sensitivity analyses for a two-dimensional OAT15A airfoil. In the same vein, a fully three-dimensional stability analysis on swept and unswept wings have been performed recently [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF][START_REF] Paladini | Transonic buffet instability: From two-dimensional airfoils to threedimensional swept wings[END_REF][START_REF] Timme | Global instability of wing shock-buffet onset[END_REF].

When turbulence is only partially modelled, as in LES or Hybrid methods, or even unmodelled (as in direct numerical simulations), the fixed point of the equations is simply the laminar solution. As discussed previously, conducting a stability analysis around this particular solution is however of very limited practical value for high-Reynolds number flows. If one aims at studying the properties of turbulent fluctuations, the mean flow thus represents the most coherent statistical object as they live in the same region in phase space (in contrast to the laminar solution). It must be noted however that, for highly turbulent flows, the nonlinear terms in the perturbation's equations may not be neglected so easily. Unfortunately, the computational cost associated to converging the expected value of these nonlinear terms is tremendously high and has been, up to now, essentially limited to the incompressible limit. Various approaches have been proposed in the literature to bypass these computations. Most of them relies on the use of the so-called resolvent operator [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF][START_REF] Lesshafft | Resolvent-based modeling of coherent wave packets in a turbulent jet[END_REF][START_REF] Symon | A tale of two airfoils: resolvent-based modelling of an oscillator versus an amplifier from an experimental mean[END_REF][START_REF] Symon | Mean and unsteady flow reconstruction using data-assimilation and resolvent analysis[END_REF]. Despite their attractiveness, such methods however raise numerous theoretical and practical questions yet unanswered and far beyond the scope of the present thesis. For the sake of clarity, the discussion about the connection between the resolvent operator and the Reynolds stresses in turbulent flows is postponed.

∂q ∂t + ∇ • F (q) = 0 (2.43)
where q = (ρ, ρu, ρE, ρ ν) T is the RANS state vector and F represents the flux vector, including the turbulent component coming from the turbulence model. The operator F can be split into the sum of a linear operator L and a nonlinear quadratic one N , so that Eq.(2.43) can be written as

∂q ∂t + ∇ • [L (q) + N (q, q)] = 0. (2.44)
If one considers the system perturbed around its steady RANS solution (i.e. base flow) q = q b , Eq.(2.44) reads ∂ (q + q ) ∂t + ∇ • L (q) + L q + N (q, q) + N q, q + N q , q + N q , q = 0, (2.45) where L (q) and N (q, q) depend only on base flow, L (q ) only on perturbations, while N (q, q ) = N (q , q) are cross-terms involving both base flow and perturbations. Focusing on the effects of a low-amplitude perturbation (q 1) evolving on top of a desirable base flow results in neglecting the nonlinear terms so that we obtain

∂q ∂t + ∇ • L q + N q, q + N q , q = 0 (2.46a) ∇ • [L (q) + N (q, q)] = 0. (2.46b)
The equations for the perturbation can be recast in matrix form as

∂q ∂t -J q = 0, (2.47)
where q = ρ , (ρu) , (ρE) , (ρ ν) Assuming a normal mode decomposition, the asymptotic behaviour of a perturbation around the base flow can be described by q = q (x) exp λt +c.c., which injected in Eq.(2.47) permits to reduce it to the classical eigenvalue problem

(J -Iλ) q (x) = 0, (2.48)
with λ = σ + iω. A base flow q will be considered as being linearly unstable if at least one of the eigenvalues of the associated linearized operator has a positive real part σ. The dynamical behaviour (in terms of frequency) is given by the imaginary part of the corresponding eigenvalue while the spatial structure of the corresponding global instability is described by the eigenmode q (x). The numerical strategy adopted to solve Eq.(2.48) will be detailed in section 3.5.

Chapter 3

Numerical method Laboratory. The tool can manage single-block as well multi-blocks cartesian grids, that can be easily generated in the CGNS standard format, enabling to handle complex geometries. Phoenix is completely written in Fortran90 and adopts the MPI libraries for parallel computing. Although the code has been equipped with a variety of spatial discretizations and temporal integration schemes from the beginning [START_REF] Goncalves | Numerical simulations of a transport-aircraft configuration[END_REF], its parallelisation is recent and only ultimated during the last 3 years. During this period, an important work has been conducted in order to generalise the data communications among different processors in case of multi-block geometries, that are crucial when dealing with structured meshes for axisymmetric geometries (i.e. nozzle flows). Last but not least, in order to keep the same spatial precision across overlapping MPI boundaries, a modular user-defined updating of ghost cells has also been implemented.

In the present chapter, a detailed description of the numerical strategies adopted to solver the linear and nonlinear governing equations will be given.

Nonlinear solver Conservative form of nonlinear compressible URANS equations

In the finite volume framework, a conservative form of dimensionless compressible URANS equations must be formulated. This is possible by defining a state vector q = [ρ, m, ρE, ρ ν] T , with m i = ρu i the momentum flow quantity in the i-th direction. Integration of Eqs.(2.23) and (2.24) over a cell with volume dΩ and face dS as

∂ ∂t Ω q dΩ + ∂Ω (F c -F v ) dS R(q) = 0, (3.1)
where R(q) is the residual vector composed by F c and F v , convective and viscous flux vector respectively. These vectors, after considering the projection dS i = n i dS of the face dS along the i-th direction, result in

F c =                   n i m i m 1 ρ n i m i + n 1 p m 2 ρ n i m i + n 2 p m 3 ρ n i m i + n 3 p E + p ρ n i m i                   F v =             0 n 1 τ tot 11 + n 2 τ tot 12 + n 3 τ tot 13 n 1 τ tot 21 + n 2 τ tot 22 + n 3 τ tot 23 n 1 τ tot 31 + n 2 τ tot 32 + n 3 τ tot 33 n 1 Θ 1 + n 2 Θ 2 + n 3 Θ 3             (3.2)

NONLINEAR SOLVER

where

Θ i = τ tot ij m j ρ + c p µ L P r L + µ T P r T ∂T ∂x i (3.3)
and the viscous stress tensor τ tot ij reading

                       τ tot 11 = 2 3 (µ L + µ T ) 2 ∂u 1 ∂x 1 - ∂u 2 ∂x 2 + ∂u 3 ∂x 3 τ tot 22 = 2 3 (µ L + µ T ) 2 ∂u 2 ∂x 2 - ∂u 1 ∂x 1 + ∂u 3 ∂x 3 τ tot 33 = 2 3 (µ L + µ T ) 2 ∂u 3 ∂x 3 - ∂u 1 ∂x 1 + ∂u 2 ∂x 2                        τ tot 12 = τ tot 21 = (µ L + µ T ) ∂u 1 ∂x 2 + ∂u ∂x τ tot 13 = τ tot 31 = (µ L + µ T ) ∂u 1 ∂x 3 + ∂u ∂x τ tot 23 = τ tot 32 = (µ L + µ T ) ∂u 2 ∂x 3 + ∂u ∂x . (3.4)
Primitive variables are computed from the conservative ones as follows

                 u i = q i+1 q 1 p = (γ -1) q 5 - 1 2 q 2 2 + q 2 3 + q 2 4 q 1 T = pγ q 1 .
(

Afterwards, one can collect all the flux terms in the residual vector R(q) as

R(q) =                 q 2 dS 1 + q 3 dS 2 + q 4 dS 3 f 11 dS 1 + f 12 dS 2 + f 13 dS 3 f 12 dS 1 + f 22 dS 2 + f 23 dS 3 f 13 dS 1 + f 23 dS 2 + f 33 dS 3 fe 1 dS 1 + fe 2 dS 2 + fe 3 dS 3                
where

                             f 11 = q 2 • q 2 q 1 + p -τ tot 11 f 12 = q 3 • q 2 q 1 -τ tot 12 f 13 = q 4 • q 2 q 1 -τ tot 13 f 22 = q 3 • q 3 q 1 + p -τ tot 22 f 23 = q 4 • q 3 q 1 -τ tot 23 f 33 = q 4 • q 4 q 1 + p -τ tot 33 , (3.6)
with

       fe 1 = q 5 + p -τ tot 11 q 2 -τ tot 12 q 3 -τ tot 13 q 4 /q 1 + qc tot 1 fe 2 = q 5 + p -τ tot 22 q 3 -τ tot 21 q 2 -τ tot 23 q 4 /q 1 + qc tot 2 fe 3 = q 5 + p -τ tot 33 q 4 -τ tot 31 q 2 -τ tot 32 q 3 /q 1 + qc tot 3 and qc tot i = -c p µ L P r L + µ T P r T ∂T ∂x i . (3.7)
The closure problem depends on turbulent viscosity µ T , which is initialised as usual in literature at 0.1µ L and then updated at each time step after solving the turbulence model equation.

LINEAR SOLVER

Conservative form of Spalart Allmaras model

As for URANS equations, when a finite volume approach is chosen, the integral form rather than the differential equation (2.27) of the Spalart-Allmaras turbulence model is preferred,

∂ ∂t Ω ρ ν dΩ + ∂Ω F T c -F T v dS R T - Ω Q T dΩ = 0, (3.8)
where Ω represents the control volume, ∂Ω its surface, and dS a surface element of Ω.

The turbulent convective flux F T c is defined as

F T c = ρ ν (n j u j ) , ( 3.9)
while the turbulent viscous flux F T v is given by

F T v = n j τ T jj (3.10)
where n j is the j-th component of the unit normal vector.

The normal viscous stresses reads

τ T jj = 1 σ [ρ (ν + ν)] ∂ ν ∂x j (3.11)
and the source term in Eq.(3.8) is

Q T = c b1 (1 -f t2 ) Sρ ν + c b2 σ ∂ ν ∂x j ∂ρ ν ∂x j -ρ ν c w1 f w - c b1 κ 2 f t2 ν d 2 .
(3.12)

In the present work, the Edwards modification is adopted in the integral form of the turbulence model.

Linear solver

The linear solver adopted in this work has been developed at Dynfluid Laboratory [START_REF] Guiho | Instabilities in oblique shock wave/laminar boundarylayer interactions[END_REF] and integrated in Phoenix in order to use the same numerical strategies [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF]. In the next paragraph, as done for nonlinear URANS equations, an accurate description of the linear N-S equations will be provided.

LINEAR SOLVER

Conservative form of the linearised compressible URANS Equations

Starting from the conservative form of the compressible URANS equations (3.1), and by linearising the latter about the base flow q = ρ, m, ρE, ρ ν T , it is possible to obtain the following Linearised

Compressible URANS (LCURANS) equations

∂ ∂t Ω q dΩ + ∂Ω F c -F v dS R(q ) = 0, (3.13)
where q = [ρ , m , (ρE) ] T is the state vector of the perturbation (in conservative form) and R(q ) the residual vector composed by F c and F v , the linearised convective and viscous flux vector respectively.

It is important to note that the linear solver has been developed in analogy with the nonlinear one so that the conservative quantities are the independent variables while primitives are functions of the former. In this perspective, primitive variables perturbations are obtained as linear expansion about q of the corresponding formula based on conservative variables.

The residual vector can be rearranged as its equivalent nonlinear ones in the following way

R (q) =                 q 2 dS 1 + q 3 dS 2 + q 4 dS 3 f 11 dS 1 + f 12 dS 2 + f 13 dS 3 f 12 dS 1 + f 22 dS 2 + f 23 dS 3 f 13 dS 1 + f 23 dS 2 + f 33 dS 3 fe 1 dS 1 + fe 2 dS 2 + fe 3 dS 3                
where

                                 f 11 = q 2 q 1 2q 2 -q 2 q 1 q 1 + p -τ 11 f 12 = f 21 = 1 q 1 q 2 q 3 + q 3 q 2 -q 2 q 3 q 1 q 1 -τ 12 f 13 = f 31 = 1 q 1 q 2 q 4 + q 4 q 2 -q 2 q 4 q 1 q 1 -τ 13 f 22 = q 3 q 1 2q 3 -q 3 q 1 q 1 + p -τ 22 f 23 = f 32 = 1 q 1 q 3 q 4 + q 4 q 3 -q 3 q 4 q 1 q 1 -τ 23 f 33 = q 4 q 1 2q 4 -q 4 q 1 q 1 + p -τ 33 , ( 3.14)
with the total stress tensor perturbation τ ij † , that can be split in its laminar and turbulent contribution as

τ ij = τ lam + τ turb , ( 3.15)
3.2. LINEAR SOLVER which are defined as

τ lam = µ L ∂u i ∂x i + ∂u j ∂x j - 2 3 δ ij ∂u k ∂x k + µ L ∂u i ∂x i + ∂u j ∂x j - 2 3 δ ij ∂u k ∂x k (3.16
)

τ turb = µ T ∂u i ∂x i + ∂u j ∂x j - 2 3 δ ij ∂u k ∂x k + µ T ∂u i ∂x i + ∂u j ∂x j - 2 3 δ ij ∂u k ∂x k , ( 3.17)
with the following primitive variables perturbations

                       ρ = q 1 u i = 1 q 1 m i - m i q 1 q 1 p = (γ -1) q 5 + 1 2q 1 q 1 q 1 q 2 2 + q 3 2 + q 4 2 -2 (q 2 q 2 + q 3 q 3 + q 4 q 4 ) T = γ q 1 p - p q 1 q 1 .
(3.18)

The most sensible part of the linearised equations is the energy equation, that is directly linked to the state equation through the heat flux perturbation q c = q c 1 , q c 2 , q c 3 T, † in the following way

fe i = A i -B i + qc i (3.19) with                        τ tot ij = τ lam + τ turb qc i = -c p µ L P r L ∂T ∂x i + µ L P r L ∂T ∂x i -c p µ T P r T ∂T ∂x i + µ T P r T ∂T ∂x i = qc lam + qc turb A i = q 5 u i + q 5 u i + pu i + p u i B i = τ tot ij u j + τ ij u j (3.20)

Conservative form of Linearised SA model with Edward's modification

The choice to use the SA model is related to its differentiability. In fact, as done for LURANS equations, also the turbulence for perturbations has been obtained by linearising Eq.(3.8) about the base flow q. In the following formula, for the sake of simplicity, the superscript has been neglected for all perturbative terms. Conversely, all terms with overbar symbol refer to base flow quantities.

The linearised version of SA model has the same form as Eq.(3.8), with the linearised turbulent convective flux F T c defined as 3.2. LINEAR SOLVER

F T c = q 6 m i q 1 + q 6 q 1 m i - m i q 1 q 1 n i , ( 3.21)
and the linearised turbulent viscous flux F T v given by

F T v = 1 σ (µ + q 6 ) ∂ ∂x i q 6 q 1 + (µ + q 6 ) ∂ ∂x i q 6 q 1 - q 6 q 1 • q 1 q 1 n i (3.22)
where n i is the i-th component of the unit normal vector and q 6 is the turbulent conservative variable ρ ν. The source term in Eq.(3.8), including the Edward's modification, is linearised as

Q T =c b1 Sq 6 + Sq 6 + c b2 σ ∂q 6 ∂x i ∂ ∂x i q 6 q 1 - q 6 q 1 • q 1 q 1 + ∂q 6 ∂x i ∂ ∂x i q 6 q 1 + c w1 d 2 f w q 6 ν + f w νq 6 + q 6 νf w . (3.23)
with the pseudo-turbulent viscosity perturbation defined as

ν = 1 q 1 q 6 - q 6 q 1 q 1 (3.24)
and the other auxiliary functions reading

χ = 1 µ q 6 - q 6 µ µ , f v 1 = 3C 3 v1 χ 2 χ 3 + C 3 v1 2 χ, ( 3.25
) (3.28)

W = 8 3
∂u i ∂x i ∂u i ∂x i + 2 ∂u i ∂x j + ∂u j ∂x i i =j ∂u i ∂x j + ∂u j ∂x i i =j , ( 3.26)
S = W 1 2 1 2 µ q 6 + f v 1 W W + µ q 6 - µ q 6 2 q 6 + f v 1 , g = r • 1 + C w2 6r 5 -1 , (3.27) r = δ tanh(1.0)    4e -2•δ e -2•δ + 1 2    ν ν - S S , δ = ν Sκ 2 d 2 , f w = g 1 + c 6 w3 g 6 + c 6
The closure problem depends on the following equation linking turbulent pseudo-viscosity perturbation and turbulent viscosity one

µ t = f v 1 q 6 + f v 1 q 6 .
(3.29)

SPATIAL DISCRETISATION

Spatial discretisation

Phoenix is a cell-centered solver. Therefore, after initialising the conservative variables at each cell, the residual vector is discretised in order to integrate Eq. Phoenix adopts the first and the third strategy in the case of Central Schemes with Artificial Dissipation and in the case of Upwind Schemes, respectively. In this work the central scheme used to discretise the URANS equations is a third-order Directional Non Compact (DNC) scheme, while a second order TVD scheme based on the upwind Roe's scheme has been chosen for the discretisation of the turbulence equation. The same numerical strategies are adopted in solving both nonlinear and linear URANS equations.

Spatial discretisation of URANS equations

In the case of hexahedral structured grids, the fluxes balance

∂Ω FdS = ∂Ω (F c -F v ) dS in
Eq.(3.1) or Eq.(3.13) at each cell is based on the following finite difference approach

∂Ω FdS = R I+1/2,J,K -R I-1/2,J,K R I + R I,J+1/2,K -R I,J-1/2,K R J + R I,J,K+1/2 -R I,J,K-1/2 R K (3. 30
)
where I ± 1/2, J ± 1/2, K ± 1/2 represent the face in the positive/negative i-th direction, respectively.

From now on, for the sake of simplicity, let us consider only the contribution to fluxes balance along 3.3. SPATIAL DISCRETISATION the I direction.

If we introduce the spatial operators

µ () I+1/2 = 1 2 () I+1 + () I (3.31a) δ () I+1/2 = () I+1 -() I , ( 3.31b)
a second order scheme reads

R I = µR I+1/2 -µR I-1/2 = δµR I + O ∆x 2 1 , (3.32)
with the contributions in the other directions R J , R K built in the same way.

After explicitly correcting the dispersive error, one can obtain the following fourth order DNC scheme

R I = δ I - 1 6 δ 2 µR I + O ∆x 4 1 . (3.33)
However, both schemes in Eq.(3.32) and Eq.(3.33) are not dissipative and need artificial dissipation to stabilise convective fluxes F c (or F c ). For this reason, people generally prefer separating the viscous contribution from the convective one. In this work, the same spatial scheme is adopted both for F c (or F c ) and F v (or F v ) so that

R I = R c I + R v I , ( 3.34a)
R c I = δ I - 1 6 δ 2 µF c I , ( 3.34b)
R v I = δ I - 1 6 δ 2 µF v I . (3.34c)
In order to stabilise the above scheme, the author chose to use a Scalar Dissipation Scheme, also known as JTS (Jameson-Schmidt-Turkel) scheme, which was firstly proposed by Jameson et al. (1981) [START_REF] Jameson | Numerical solution of the euler equations by finite volume methods using runge kutta time stepping schemes[END_REF] to solve the Euler equations. In this solution the artificial dissipation d I+1/2 consists in -

d I+1/2 = |µΛ I+1/2 c | • I+1/2 (2) δq I+1/2 +O(∆x 1 ) q I+1 -q I + -|µΛ I+1/2 c | • I+1/2 (4) q I+2 -3q I+1 + 3q I -q I-1 δ 3 q I+1/2 +O(∆x 3
ν I = |p I+1 -2p I + p I-1 | p I+1 + 2p I + p I-1 ∝ | ∂ 2 p ∂x 2 ∆x 2 + O(∆x 4 )| § ; - I+1/2 (2)
= k (2) max ν I+1 , ν I ; -

I+1/2 (4) = max 0, k (4) - I+1/2 (2) =O(1)
where a is the sound velocity.

The pressure-based sensor

I+1/2 (4)
is used to switch off the third order term δ 3 q I+1/2 at shocks and suppress any high-order schemes oscillations near discontinuities. Additionally, the sensor

I+1/2 (2)
switches off the first-order term δq I+1/2 in smooth parts of the flowfield, in order to minimise the dissipation. The contributions in the other directions d J+1/2 , d K+1/2 are built in the same way so that the overall artificial dissipation d at each cell is computed as

d = d I+1/2,J,K -d I-1/2,J,K d I + d I,J+1/2,K -d I,J-1/2,K d J + d I,J,K+1/2 -d I,J,K-1/2 d K , ( 3.36)
and added to the convective discretised flux R c I so that Eq.(3.34) becomes

R I = R c I + d I + R v I , ( 3.37a)
d I = δµd I . (3.37b)
Since artificial dissipation in Eq. (3.35) results in a blend of third-order term and a first-order term, which is enabled only in strong pressure gradient region, it is traditionally accepted to be a third order accurate term. Therefore a fourth order DNC scheme with a such artificial dissipation is commonly accepted as a globally third order scheme.

As for viscous fluxes R v , the τ ij cell-centred value is computed after evaluating the velocity gradient components ∂u i /∂x j , as well as the temperature gradient ones ∂T /∂x i , by means of Green's Theorem

as ∂A ∂x i = 1 Ω ∂Ω AdS i ≈ 1 Ω N F m=1 A m n i,m dS (3.38) ‡
In the linear solver Λc is computed using the baseflow variables as Λc = ui dSi + a dS. § In the linear solver the pressure sensor is based on the baseflow pressure p to avoid any division by zero.

SPATIAL DISCRETISATION

with N F number of cell faces and A m quantity face values. In the present work, the latter is discretised by a 2nd order accurate finite difference approach as Eq.(3.32), so that the above formula reduces to the following discretised version

∂A ∂x i ≈ 1 Ω M δ µA M dS M i M = I, J, K. (3.39)

Spatial discretisation of turbulence transport equation

Eq.(3.8), as its linear counterpart, includes quantities depending only on flowfield variables and others depending on gradients. In the present work, the former are evaluated at each cell-center while the latter are computed using Eq.(3.39). After computing each single term, the turbulent residual R T is discretised by a 2nd order TVD scheme.

A TVD scheme consists in a central scheme obtained as an average of the convective fluxes combined with an additional limited dissipation term preventing the generation of new extrema in the flow solution [START_REF] Blazek | Computational fluid dynamics: principles and applications[END_REF]. In this work, the central scheme adopted for turbulent vector flux R T (including the viscous and convective part) is a 2nd order scheme (see Eq.(3.32)). In order not to downgrade the precision, a second order dissipation term is adopted. For the sake of simplicity, if we call ρ ν with q and consider only the contribution on the face (I -1/2, J, K), the aforementioned dissipation

d I-1/2
is obtained by a preprocessing approach in terms of Roe's solver [START_REF] Chakravarthy | Development of upwind schemes for the euler equations[END_REF] as

d I-1/2 = 1 2 |µλ I-1/2 | q I-1/2 R -q I-1/2 L dS, ( 3.40)
where q R and q L are the approximations of q on the right and left sides of the I-th cell and computed as q In the present work, the upwinding term in Eq.(3.41) is approximated as

I-1/2 R -q I-1/2 L = δq I-1/2 - 1 2 ψ δq I-1/2 , δq I+1/2 + ψ δq I-1/2 , δq I-3/2
minmod δq I-3/2 , minmod δq I-1/2 , δq I+1/2 , ( 3.44)
and the resulting 2nd order TVD scheme is

R T | I = δµR T | I + δµd I . (3.45)
The dissipation term, as its equivalent in Roe's scheme and many other upwind schemes, would produce a non-physical expansion shocks failing to satisfy the entropy condition at sonic point. In order to solve this problem, the modulus of the convective eigenvalue can be modified using the following Harten's entropy correction

|Λ T c | =        |Λ T c |, |Λ T c | > c Λ T c 2 + c 2 2c , |Λ T c | ≤ c (3.46)
where c is a small value, commonly set equal to some fraction of the local speed of sound a. In the present work c is set at 5% of velocity of sound a.

Time integration

In this work, an implicit loosely coupled algorithm is used to integrate the (linear/nonlinear) URANS equations and (linear/nonlinear) Spalart-Allmaras turbulence equation separately [START_REF] Park | Implementation of κ -ω turbulence models in an implicit multigrid method[END_REF]. Historically, this approach has been preferred due to the fact that the turbulence model equations are solved only to obtain the eddy viscosity as well as to its flexibility in using different turbulence model

TIME INTEGRATION

Time Stepping (DTS) method has been used to solve steady and unsteady problems respectively. In the next paragraph, a detailed description of the above time-integration methods will be provided.

Steady Flows: Local time stepping (LTS) technique

URANS equations

The or Eq.(3.13), after a linearisation in time, we obtain the following implicit scheme at each cell center

1 ∆t L ΩI + ∂ R ∂q n ∆q = -R n , ( 3.47)
where

R n = R n I + R n J + R n K and ∆q = q n+1 -q n .
It is worthwhile noting that if ∆t L tends to infinity, the above scheme reduces to a Newton's method, which is known to have a quadratic convergence property. However this property is guaranteed only if the Jacobian matrix ∂ R/∂q is exact.

Unfortunately, explicit expression for the Jacobian matrix resulting from the exact linearisation of any second or higher order numerical flux is extremely expensive. A solution for compressible structured solver has been proposed by [START_REF] Luo | A fast, matrix-free implicit method for compressible flows on unstructured grids[END_REF] [START_REF] Luo | A fast, matrix-free implicit method for compressible flows on unstructured grids[END_REF] under the name of Matrix-Free Method. This algorithm is based on the simplification of the Jacobian matrix with ∂ R * /∂q, where R * is the first order Rusanov numerical inviscid flux defined as

R * = δµR * I + δµR * J + δµR * K (3.48) with R * I+1/2 = 1 2 F c I+1 + F c I -|Λ I+1/2 c | q I+1 -q I = µ F c I+1/2 - |µΛ I+1/2 c | 2 δq I+1/2 . ( 3
∆q i Ω ∆t L + |µΛ I+1/2 c | 2 + |µΛ I-1/2 c | 2 D I = -R I i + |µΛ I+1/2 c | 2 ∆q I+1 i + |µΛ I-1/2 c | 2 ∆q I-1 i -δ µdF I c,i dS I (3.50)
where dF c,i is the exact differential of the i-th component of convective flux F c . To take in account the contribution of viscous term F v that is not simultaneously diagonalizable, an approximation of viscous spectral radius Λ v is added in the diagonal operator D I [START_REF] Luo | A fast, matrix-free implicit method for compressible flows on unstructured grids[END_REF][START_REF] Park | Implementation of κ -ω turbulence models in an implicit multigrid method[END_REF] as

D I = Ω ∆t L + |µΛ I+1/2 c | 2 + |µΛ I-1/2 c | 2 + µΛ I+1/2 v dS I+1/2 µΩ I+1/2 + µΛ I-1/2 v dS I-1/2 µΩ I-1/2 (3.51)
where

Λ I v = µ I L P r L + µ I T P r T γ ρ I . ( 3.52)
Under the assumption that

dF c ≈ F c (q + ∆q) -F c (q), (3.53)
Eq.(3.50) can be iteratively resolved to find ∆q and advance the solution in time.

Turbulence transport equation

An Euler implicit scheme, similar to that described in the previous paragraph, is employed for the time integration of the discretised turbulence transport equation. However, when dealing with turbulence models, the source term becomes so dominant that flow variables changes rapidly in space and in time at much smaller time scales than those of RANS equations. This process significantly increases the stiffness of Eq.(3.8) and the usual remedy to that is treating the source term Q T † † in an implicit way. For this purpose, the linearisation of the source term about the current time n results in

Q T n+1 = Q T n + ∂Q T ∂q 6 n q n+1 6 -q n 6 (3.

54)

In the linear solver

Λ I v = µL I P rL + µT I P rT γ ρ I . † †
The lack of symbol means that the term is well defined at each cell without the need to be approximated.

TIME INTEGRATION

with q 6 = ρ ν. If one uses Eq.(3.54) to make Q T implicit in Eq.(3.8) and adopts the same linearisation in time for the turbulent flux balance R T , an equivalent Euler implicit scheme for the turbulence model can be obtained as

1 ∆t L + ∂Q T ∂q 6 n Ω + ∂R T ∂q 6 n ∆q 6 = -R T -Q T Ω n , ( 3.55)
with ∂R T /∂q 6 = ∂F T c /∂q 6 + ∂F T v /∂q 6 dS. In the case of one or two turbulence transport equations, the implicit terms ∂Q T /∂q 6 and ∂F T c /∂q 6 are diagonal and easy to obtain analytically. Whilst, the turbulent viscous counterpart ∂F T v /∂q 6 is usually estimated with its spectral radius Λ T v as

∂F T v ∂q 6 ≈ Λ T v = 2 (µ L + µ T ) ρ ‡ ‡ (3.56)
and added to the diagonal operator as well, giving no reason to use a Matrix-Free approach.

In this work, Eq. and advance the turbulent variable q 6 in time.

Unsteady Flows: Implicit Dual Time Stepping (DTS) technique

URANS equations

The efficiency of Matrix-Free method in solving steady problems can be exploited to solve unsteady

Ω ∆t 3 2 q n+1 -2q n + 1 2 q n-1 = -R n , (3.57)
requires the solution of an equivalent pseudo-steady problem as

Ω ∂q * ∂t * = -R(q * ) (3.58) ‡ ‡ In the linear solver Λ T v = 2 (µL + µT ) ρ .

TIME INTEGRATION

with q * approximation at the time t * of q n+1 and the unsteady residual R(q * ) reading Matrix-Free approach, one can come to the following DTS implicit scheme

R(q * ) = R(q * ) + 3Ω 2∆t q * + Ω ∆t -2q n + 1 2 q n-1 . ( 3
           R l = R(q l ) + 3Ω 2∆t q l + Ω ∆t -2q n + 1 2 q n-1 1 ∆t L + 3 2∆t ΩI + ∂ R * ∂q q l+1 -q l = -R l .
(3.60)

Before advancing the solution in each cell q of the physical time step ∆t, the system above is marched in the fictitious time ∆t L by a LTS technique for a number of iterations l so that the unsteady residual R l drops at least of one order of magnitude. After that, the solution q n-1 is updated with q n as well as q n with the converged pseudo-steady solution q l+1 . As one can state by comparing Eq.(3.60) with Eq.(3.47), the pseudo-steady implicit formula differs from the the classical Euler one only by the diagonal term 3/ (2∆t) ΩI. This means that, compared to an explicit scheme, the DTS implicit scheme has a higher computational cost for each physical iteration. However, the advantage of DTS method is the possibility of advancing the solution with a grid resolution-independent physical time step, which is much higher of that required by an explicit scheme on the same grid. Given that the greatest level of grid refinement is usually demanded near wall, this numerical property of DTS scheme perfectly blends with hybrid formulations, which have been developed to save computational cost inside boundary layer where a LES approach is not strictly required.

Turbulence transport equation

An analogous version of Eq.(3.60) can be obtained for Eq.(3.55) as

           R 6 l = R 6 (q l 6 ) + 3Ω 2∆t q l 6 + Ω ∆t -2q n 6 + 1 2 q n-1 6 1 ∆t L + ∂Q T ∂q 6 + 3 2∆t Ω + ∂R T ∂q 6 q l+1 6 -q l 6 = -R 6 l .
(3.61)

TIME INTEGRATION

with q 6 = ρ ν and

R 6 = R T -Q T .
Giving that Eq.(3.61) has the same form of Eq.(3.55), the former is numerically solved as the latter to advance the turbulent variable q 6 in time.

Steady Flows: Selective Frequency Damping (SFD)

In contrast with laminar flows, turbulent flows require an higher grid refinement level or, as it happens for URANS methods, an additional turbulence transport equation in order to accurately resolve or model wall turbulence. As mentioned in section 3.4.1.2, a turbulence transport equation includes stiff source terms which, together with the grid refinement-sensitiveness from which the Local Time Stepping (LTS) technique suffers, does not permit to converge the steady state with reasonable computational resources or even with a reasonable precision by this method. In order to bypass this drawback, people generally adopt a multigrid technique, which consists in converging at best the steady solution on a rough grid and then interpolating the obtained solution on a more refined grid and iterate this process. Another expedient, which does not require any differently refined grids, is the Selective Frequency Damping (SFD).

The SFD approach consists in considering the continuous URANS equations as a nonlinear system q = f (q), where q is the vector of conservative quantities (including the pseudo-turbulent viscosity q 6 = ρ ν)

and f (q) = -∇ • F (q) the RHS vector (including any turbulence model term) with appropriate initial and boundary conditions for q. Therefore, the base flow q b is a steady solution, given for q = f (q b ) = 0. In the case of linear instability, if the system is very small perturbed about its steady solution, any solution q will quickly depart from q b . In this perspective, computing the base flow means stabilising the above system.

According to the Linear Control Theory, one can design a feedback control loop to stabilise an unstable system. The SFD method, by using exactly this strategy, applies a low-pass frequency filter to the system, obtaining a filtered solution q, in the following way:

q = f (q) -χ (q -q) q = (q -q) /∆ (3.62)
where ∆ is inverse of the cut-off frequency (ω c ) and χ is the control coefficient. Usually, ω c is set at an half of the dominant disturbance frequency, while χ twice the growth rate.

TIME INTEGRATION

The version implemented in this work is that proposed by Jordi et al. (2014) [START_REF] Jordi | Encapsulated formulation of the selective frequency damping method[END_REF], which avoids any modification of the original solver. According to this, the proportional filter is activated on the output q j of the j-th DTS iteration. In this view, Equations (3.62) can be written in a compact way as

d dt q q = -χI χI 1 ∆ I -1 ∆ I q * q (3.63)
and exactly integrated over one time step ∆t, yielding

q j+1 q j+1 = 1 1 + χ∆ I + χ∆eI χ∆(I -eI) I -eI χ∆I + eI q * q j (3.64)
where e = e -(χ+1/∆)∆t . In this way, at the j-th iteration of the unsteady solver, the solution q j is updated and filtered, until q n ≈ q n + , with imposed by a stop criterion. This method has been successfully tested by Richez et al. (2016) [START_REF] Richez | Selective frequency damping method for steady rans solutions of turbulent separated flows around an airfoil at stall[END_REF] in finding the steady RANS solution in a turbulent separated flow past an airfoil.

Even though this formulation is the more convenient when dealing with complicated solver because it only needs an output correction, the control parameters setting (∆, χ) is more tricky and if they are not carefully chosen a steady solution can be even not reached. Particularly, Casacuberta et al.

(2018) [START_REF] Casacuberta | Effectivity and efficiency of selective frequency damping for the computation of unstable steady-state solutions[END_REF] have shown that when flow exhibits only one instability mode (σ, ω), the steady solution is always obtained for the following choice of parameters

χ = |σ| + ω 2 , (3.65a) ∆ = 2 |σ| -ω . (3.65b)
It is evident that in order to successfully converge towards a steady state with such a method, an a priori knowledge of the instability is required. However, while ω can be easily obtained by a nonlinear calculation, only a linear stability analysis about the base flow would provide σ. In order to solve this ambiguity, Casacuberta [START_REF] Casacuberta | Effectivity and efficiency of selective frequency damping for the computation of unstable steady-state solutions[END_REF] suggested the so-called flow unleash method, which consists in first driving the solution at low enough residual and restarting the calculation from that point. In this way, instabilities would grow and solution will depart from the converged steady state with a growth

TIME INTEGRATION

rate σ * close to σ, which is employed to update the value of (∆, χ) and activate the SFD back. The more accurate the initial steady solution is, the less the process must be iterated. An application of this method for turbulent flows will be provided in section 4.2.1.

EIGENVALUE SOLVER

Eigenvalue solver

The unsteady Reynolds-averaged Navier-Stokes equations linearized in the vicinity of a given base flow read ∂q ∂t + ∇ • L q + N q, q + N q , q = 0.

(3.66)

For the sake of simplicity, we will however use hereafter the following notation

∂q ∂t -J q = 0, (3.67)
where J denotes the Jacobian operator of the URANS equations and q = ρ , (ρu) , (ρE) , (ρ ν) T is the state vector (in conservative form) of the system. For a fully three-dimensional flow, this Jacobian operator can be represented as a kn × kn matrix (k denoting the number of conservative variables and n the number of grid points in the domain). For the flow configurations investigated in the present work, computing the leading eigenvalues and eigenvectors of such a large matrix is currently intractable using standard direct solvers. Consequently, the time-stepper formulation of the instability problem coupled with an iterative eigenvalue solver (e.g. the Arnoldi method), a technique originally proposed by Edwards et al. (1994) [START_REF] Edwards | Krylov methods for the incompressible navier-stokes equations[END_REF] and popularized by Bagheri et al. (2009) [START_REF] Bagheri | Matrix-free methods for the stability and control of boundary layers[END_REF], is used in this thesis.

The solution to Eq. (3.67) can be written as q (τ ) = e τ J q (0)

where M = e τ J is a matrix known as the exponential propagator of the system, mapping an initial condition at time t = 0 to the corresponding solution at time t = τ . From a practical point of view, this matrix-vector product can be easily computed by marching forward in time the linearized URANS equations using q (0) as the input vector. It is then relatively easy to compute the leading eigenvalues and eigenvectors λ M q = Mq of this exponential propagator using a standard iterative eigenvalue solver. Note that the eigenvalues of M are related to those of the Jacobian operator J through the relation

λ J = log λ M τ
while their eigenvectors are identical. These equivalences come from the very definition of the matrix exponential.

EIGENVALUE SOLVER

Although it has become relatively standard nowadays, let us explain the basic idea behind the Arnoldi algorithm used to compute these eigenpairs for the sake of completeness and clarity. It can be considered an advanced variant of the simpler power iteration method. Given the exponential propagator M and a normalised initial vector u 0 = q 0 / q 0 , the power iteration constructs the following m-

dimensional Krylov subspace K m K m (M, u 0 ) = {u 0 , Mu 0 , . . . , M m-1 u 0 }. (3.68)
It can be shown that as m increases, the last vector of this sequence slowly converges toward the leading eigenvector of M while the corresponding eigenvalue can be recovered from the associated Rayleigh quotient. Moreover, raising the exponential propagator to a given power k simply amounts in practice to integrate the linearized equations foward in time from time t = 0 to time t = kτ where τ is our sampling period. It should be emphasized however that the Krylov sequence K m becomes increasingly ill-conditioned as m increases.

In order to overcome this drawback, the Arnoldi method combines the power iteration with a modified

Gram-Schmidt process. The Arnoldi factorisation reads

MV m = V m H m + βv m+1 e T m+1 (3.69)
with V m an orthonormal set of vectors, H m a m×m upper Hessenberg matrix, and βv m+1 the residual which, by construction, is orthogonal to V m . Figure 3.1 visualises the shape of the matrices involved in the process. It can be shown that the pair (λ H , V m ŷ), where λ H and ŷ are the eigenvalues and the corresponding eigenvectors of H m , provides a good approximation for the eigenpairs of the matrix M.

The Arnoldi procedure to be carried out in order to construct iteratively these different matrices is described step by step in Table 3.1. 3.5. EIGENVALUE SOLVER Table 3.1: Arnoldi algorithm step-by-step Set an initial arbitrary vector u 0 and normalize it such that it has unitary norm. v 1 = q 0 q 0 The new Krylov vector v k is computed from v k-1 by integrating forward in time the linearized equations.

v k = Mv k-1
A Gram-Schmidt iteration is then used to orthonormalize v k with respect to the previous Krylov vectors. This is accomplished in two steps. h j,k-1 are the terms in the upper diagonal Hessenberg matrix.

f or j = 1 → k -1 h j,k-1 = v j v k v k = v k -h j,k-1 v j end
The subdiagonal term of the Hessenberg matrix is the norm of the new orthogonal vector.

h k,k-1 = v k Finally v k is normalised to get the new orthogonal direction in the Krylov subspace.

v k = v k h k,k-1
Repeat from step 2 until k is equal to the desired dimension m of the Krylov subspace.

while k ≤ m k = k + 1 goto step 2 end
Solve the small m × m eigenvalue problem using LA-PACK [START_REF] Anderson | LAPACK Users' guide[END_REF] [µ H , ŷ] = eig(H)

Approximate the eigenvectors of the exponential propagator M from the orthonormal Krylov basis V m and its eigenvalues from those of the Hessenberg matrix.

λ = µ H and q = V m ŷ
Recover the eigenvalues of the Jacobian matrix J from those of the propagation matrix M and the relative residual r k associated at each qk .

Boundary conditions

Navier-Stokes as well as URANS equations require initial conditions as well as boundary conditions.

The latter include both physical conditions that solutions must respect (i.e. symmetry, solid walls, etc.) and artificial ones derived from the fact that any numerical simulation can consider only a part of the real physical domain (i.e. farfield for inflow/outflow, periodic boundaries, continuity over overlapping MPI boundaries, etc.) [START_REF] Blazek | Computational fluid dynamics: principles and applications[END_REF]. In the next paragraphs, a detailed description of boundary conditions used in this work will be provided.

Farfield non-reflecting boundary condtions

In order to avoid any outgoing disturbances be reflected back into the computational domain at boundaries, non-reflective boundary conditions based on characteristics are adopted. Depending on the sign of the eigenvalues of the convective flux Jacobians ∂F c /∂q, disturbances are conveyed outside or inside the domain along the characteristics. For example, if one considers the entire set of

eigenvalues Λ c = [u • n, u • n, u • n, u • n + a, u • n -a]
T , in the case of subsonic inflow there is only one positive eigenvalues and consequently an outgoing characteristic. Conversely, for subsonic outflow the situation reverses. The approach chosen in this work is that of Whitfield [START_REF] Whitfield | Three-dimensional unsteady euler equations solution using flux vector splitting[END_REF], which is based on the one-dimensional Euler equations. If we denote with the supscript ∞ the free-stream condition and with S that about a point just inside the domain, primitive variables at boundary B must be specified as follows

SUPERSONIC INFLOW =              p B = p ∞ ρ B = ρ ∞ u B = u ∞ SUPERSONIC OUTFLOW =              p B = p S ρ B = ρ S u B = u S (3.70) 3.6. BOUNDARY CONDITIONS SUBSONIC INFLOW =                            p B = 1 2 [p ∞ + p S -ρ S a S (u ∞ • n -u S • n)] ρ B = ρ ∞ + p B -p ∞ a 2 S u B = u ∞ -n p ∞ -p B ρ S a S (3.71) SUBSONIC OUTFLOW =                      p B prescribed ρ B = ρ S + p B -p S a 2 S u B = u S + n p S -p B ρ S a S (3.72)
Free-stream variables are imposed in terms of Mach number M thanks to the practical reference quantities (see section 2.1.1 for further details). As for the linear solver, the above formula are kept in the same for baseflow variables or linearised for perturbations, as in Crouch (2007) [START_REF] Crouch | Predicting the onset of flow unsteadiness based on global instability[END_REF]. Particularly, for the case of subsonic outflow, a null pressure perturbations (p = 0) is prescribed at boundary [START_REF] Guiho | Instabilities in oblique shock wave/laminar boundarylayer interactions[END_REF].

Injection boundary conditions

A non-reflecting boundary condition for the Euler equations consisting in the specification of the total pressure, total temperature, and two flow angles is adopted.

Let us consider a boundary with prescribed upstream total quantities (p * 0 , T * 0 ). The static pressure at this boundary follows

p * = 1 + γ -1 2 M 2 -γ γ-1 p * 0 , ( 3.73)
with the superscript * referring to dimensional quantities. Rewriting the above equation for dimensionless variables defined in §2.1.1 leads to

p = 1 γ RP I 1 + γ -1 2 M 2 -γ γ-1
(3.74)

BOUNDARY CONDITIONS

where RP I represents the ratio between the prescribed total pressure p * 0 and the reference total pressure p ∞ 0 . By using the isentropic relations, one can obtain the following formula linking the Mach number M and velocity u at the same station

M 2 = ||u|| 2 RT I 1 - γ -1 2
||u|| 2 RT I , ( 3.75)
with RT I defined as the ratio between the prescribed total temperature T * 0 and the reference total temperature T ∞ 0 . If Eq.(3.74) and Eq.(3.75) are coupled with the dimensionless state equation, the following system of equations is obtained

                                 b =         1 + γ -1 2 1 RT I u • n n i θ i 2 1- γ -1 2 • RT I u • n n i θ i 2         -1 p = RP I γ b γ γ-1 ρ = p • γ RT I • b (3.76)
where n i and θ i are the i-th component of the unit normal vector and the i-th flow direction cosine, respectively.

In order to close the above system of equations, in case of subsonic inflow at this boundary, Whithman [START_REF] Whitham | Linear and nonlinear waves[END_REF] suggests to solve the left running Characteristic Equation linearised about an interior point S as

dp -ρ S a S d (u • n) = 0 on C -: dx dt = -a S , (3.77)
and is easily integrated as follows

p S -ρ S a S (u • n) S = p -ρ S a S (u • n) . (3.78)
Inserting the above equation in the system (3.76), we obtain the following one that can be solved by a Newton-Raphson:

3.6. BOUNDARY CONDITIONS                                              b =         1 + γ -1 2 1 RT I u • n n i θ i 2 1- γ -1 2 • RT I u • n n i θ i 2         -1 p = RP I γ b γ γ-1 ρ = p • γ RT I • b p = ρ S a S (u • n) + p S -ρ S a S (u • n) S . (3.79)
As for linear solver, a null perturbations (q = 0) is imposed at boundary, while the same formulation is kept for base flow variables [START_REF] Guiho | Instabilities in oblique shock wave/laminar boundarylayer interactions[END_REF].

Adiabatic Wall Function boundary condition

When a sufficiently refined grid in the vicinity of walls costs too much in terms of computational resources, a wall function is adopted. In the present work, the laminar and logarithmic regions of boundary layer are identified by the following conditions

       U + s = y + y + < 11. 13
U + s = 1 κ log y + + 5.25 y + > 11.13 (3.80)
with:

U τ = τ w ρ w , U + s = V (1) t U τ y + = y (1) U τ ρ w µ w (3.81)
where the superscript (1) denotes quantities related to the first cell near wall and κ the Von-karman constant. V t is the tangent velocity, while ρ w , τ w and ν w are density, shear stress and dinamic viscosity at wall, respectively.

The first step is obtaining the value of temperature at wall T w , is accomplished by imposing the adiabatic flux condition at wall as 3.6. BOUNDARY CONDITIONS

T w = T (1) +        µ (1) + µ (1) t c p µ (1) P r + µ (1) t P r t        V (1) t 2 2 . (3.82)
Consequently, dynamic viscosity at wall µ w is computed through T w by the Sutherland's law as

µ w = µ (1) T w T (1)
1 2

1 + S/T (1) 1 + S/T w .

(3.83)

Imposing the null pressure gradient condition inside the boundary layer permits to compute ρ w as follows

dP dy = 0 =⇒ ρ w = ρ (1) T (1) T (1) w . ( 3.84)
In order to tackle compressible boundary layer formulation, the following Van Driest transformation [START_REF] Huang | Van driest transformation and compressible wall-bounded flows[END_REF] is needed

V * t = A arcsin V (1) t A (3.85)
where

A = 2c p µ (1) P r + µ (1) t P r t T w µ (1) + µ (1) t . (3.86)
At this stage, depending on whether the point considered is inside the linear region or the logarithmic one, Eq.(3.87) rather than Eq.(3.88) is iteratively solved to compute U τ (1) µ w + 5.25.

V * t U τ = ρ w U τ y (1) µ w , (3.87) V * t U τ = 1 κ ln ρ w U τ y
(3.88)

BOUNDARY CONDITIONS

Finally, stress tensor τ ij and heat flux qc i are simply obtained via τ w as

τ (1) ij = τ w ij = τ w (t i • n j + t j • n i ) , (3.89) qc (1) i = τ w V
(1) t

• n i , (3.90)

qc w i = 0. (3.91)

Spalart Allmaras boundary condition

As for the turbulent transport equation, the only variable to be specified at boundaries is q 6 = ρν.

A null-value is specified at walls, while when adiabatic wall functions are adopted, the value for the first cell near wall (denoted with the superscript (1) ) is based on a the following mixing-length formulation for turbulent viscosity µ t µ

(1)

t = ρ (1) κy (1) 2
                                                                                         C 1 = √ 3C b 256C b + 27µ 4 t 1 2 C b = C 3 v µ 3 µ t C 0 = C v µµ t =⇒ C 3 0 = C b µ 2 t C 2 = C 1 18 + 0.5C 3 0 C 3 = C 1 18 -0.5C 3 0 C 4 = 4 C 1 3 2 -C 1 3 3 rc 4 = µ 2 t + C 4 1 2 ρν = 1 4      (µ t + rc 4 ) + √ 2 rc 4 µ 2 t -2C 1 3 3 rc 4 + 2C 1 3 2 rc 4 + µ 3 t 1 2 rc 1 2 4      (3.94) Chapter 4
Validation of numerical strategy

2D laminar flow past a cylinder

The chosen laminar case is the 2D quasi-incompressible flow past a cylinder. The cartesian structured 1-block grid used for all computations is shown in Figure 4.1 and consists in (n x , n y ) = (241, 145)

points in the streamwise and wall-normal direction, respectively. The computational domain has been first divided in 4 portions, each of them independently meshed and finally connected in a single C-grid.

The flow parameters used in all computations are listed in Table 4.1.

Free-stream Mach number M ∞ = 0.1 Free-stream stagnation temperature T 0,∞ = 288 K

Reynolds number

Re ∈ [START_REF] Kistler | Fluctuating wall pressure under a separated supersonic flow[END_REF][START_REF] Riou | Improvement of delayed-detached eddy simulation applied to separated flow over missile fin[END_REF] Table 4.1: Flow parameters for the 2D quasi-incompressible cylinder flow.

No slip boundary condition with adiabatic temperature and pressure extrapolation is prescribed at the cylinder walls. Farfield boundaries are treated with a non-reflective boundary condition [START_REF] Whitfield | Three-dimensional unsteady euler equations solution using flux vector splitting[END_REF] at

imposed Mach number M = M ∞ .
The numerical parameters used for all calculations are listed in Table 4.2.

2D LAMINAR FLOW PAST A CYLINDER

Linear Stability around base flow

A Linear Stability Analysis has been performed around the base flow. In accordance with nonlinear and linear calculations, for unstable base flows, the resulting eigenspectrum (see Figure 4.10a)

shows a single complex-conjugate pair of modes with positive-real part eigenvalue σ, which additionally imaginary part ω = 2πSt is compatible with the observed limit cycle frequency. The streamwise velocity component of this unstable mode has a packet shape and localised in the cylinder wake (see As expected, linear stability analysis is less and less predictive moving away from the stability threshold, as shown in table 4.4. On the other hand, as described by Barkley [66], far from stability threshold a linear stability analysis around the meanflow would return a marginally stable mode matching the limit cycle frequency. However, this is valid only because this particular flow has no higher harmonics [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF], being a monochromatic oscillation [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF] (see Figure 4.5b).

Linear Stability around mean flow

For the cases at Re = 60-70, a Linear Stability Analysis around the meanflow has been performed.

The meanflow has been previously computed by time-averaging the unsteady solution, once limit cycle settled, for 60 periods at a convergence residual about 10 -6 . In this case, Global Stability Analysis returns a marginally stable mode at a frequency very close to the nonlinear one, as shown in Figure 4.12. The comparison in terms of (σ, St) between linear stability around the base flow and around the meanflow in Table 4.5 confirms the better frequency capture by the latter.

Re (σ, St) BF (σ, St) M F (σ, St) N L ∆St BF % ∆St M F

2D TURBULENT BUFFETING ON OAT15A

2D turbulent buffeting on OAT15A

The turbulent case chosen to validate the nonlinear and linear URANS solver is the 2D transonic turbulent flow past a OAT15A airfoil. The cartesian structured 1-block grid used for all computations is shown in

Base flow calculation: Steady RANS solution

As mentioned in section 3.4.3, turbulence models contain stiff terms lowering the capability of Local Time Stepping (LTS) techniques in converging steady state solutions. For this reason, in this case the author adopted the Selective Frequency Damping (SFD) method in order to obtain the base flow (see section 3.4.3 for further details). Particularly, the control coefficients (∆, χ) have been selected by using the flow-unleash method and due to the need for an a priori knowledge of the unsteady solution, the latter has been computed from scratch. After having obtained a first URANS solution, as shown in Figure 4.14a, the former has been initially forced towards the steady state by a LTS method. Once we reached a first approximation for base flow, the latter is used to initialise another unsteady calculation, which is run until saturation. During this trip, an approximation for instability growth rate σ * is extracted by the slope of lift coefficient trend in logarithmic scale, while frequency is taken in term of angular velocity equal to that of the observed limit cycle ω * † . Using formula in same order of LTS, this "head start" of SFD leads it to get earlier the steady solution.

(a) Lift coefficient trend under flow-unleash method. Firstly, LTS method is activated on the saturated limit cycle in order to get an approximated steady solution. Secondly, instabilities are left to naturally develop and limit cycle is reached again. Finally, SFD is properly activated to stabilise the unsteady solution.

2D TURBULENT BUFFETING ON OAT15A

Nonlinear unsteady calculation: Unsteady RANS (URANS) solution

The converged steady RANS solution, obtained by filtering the unsteady URANS one, is used as initialisation for a new unsteady nonlinear calculation to get the instability growth rate. In order to accelerate the development of instabilities, flow has been perturbed with a pulse localised at the shock foot. As expected, base flow at α = 4.5 • is unstable and flow evolves with a growth rate σ N L towards a limit cycle with a frequency f N L , expressed as Strouhal number A flow visualisation by numerical schlieren of saturated limit cycle is reported in Figure 4.17. Unsteadiness seems mainly due to a lambda-pattern shock displacing back and forth along the airfoil surface. Frequency f N L at α ∞ = 4.5 • is about 67 Hz and buffet onset is at α ∞ ≈ 4.0 • . Previous calculations and experiments have found no significant changes in buffet frequency with the AoA (Angle of Attack) α ∞ . However, as shown in Table 4.8, even if the buffet onset detected in the current simulation is slightly delayed when compared with experiments (Jacquin et al. 2009 [START_REF] Jacquin | Experimental study of shock oscillation over a transonic supercritical profile[END_REF]) or with more accurate turbulence description simulations (Deck 2005 [35], Z-DES), it seems in quite good agreement with analogous previous calculations.

St N L = f N L • c/U ∞ (see
A meanflow has been computed by time-averaging the unsteady solution, once limit cycle settled, for observed in this turbulent flow is quite different from that occurring in the laminar flow past a cylinder.

Particularly, while the laminar cylinder flow DFT spectrum shows a monochromatic oscillation with a single harmonic even far from its stability threshold, in the turbulent one there are a second and third harmonic still visible. As shown in Figure 4.16b, at α ∞ = 4.5 • the energy related to the second harmonic of the signal is about ten per cent of that for the fundamental. This means that in this flow, nonlinear effects are so strong that a linear stability analysis performed around the meanflow may not be frequency-predictable according to Sipp & Lebedev (2007) [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF] and Turton et al. (2015) [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF]. This scenario will be highlighted in the last paragraph of this chapter. However, it is interesting to note that the weight of higher harmonics in the power spectrum increases the further away we move from the buffet onset, which current nonlinear calculations returns to be located at α ∞ ≈ 4.0 • . Actually, as shown in

Linear Stability around the base flow

A Linear Stability Analysis has been performed around the base flow. In accordance with nonlinear calculation, the same stability threshold at α ∞ ≈ 4.0 • is found. The resulting eigenspectrum (see

Linear Stability around the meanflow

For the case at α ∞ = 4.5 • , a Linear Stability Analysis around the meanflow has been performed.

The turbulent meanflow has been previously computed by time-averaging the turbulent unsteady solution, once limit cycle settled, for 120 periods. In this case, two different global stability analyses 4.2. 2D TURBULENT BUFFETING ON OAT15A have been carried out: one using LNS equations and another using LURANS equations. As for stability analysis around meanflow using LNS equations, one can state by observing Figure 4.23b that such an analysis returns a mode that is not marginally stable, in contrast with laminar flows, and does not match the frequency of nonlinear limit cycle. As previously mentioned, it is in the author opinion that this apparent inconsistency is due to the greater importance of nonlinearities present in the flow when operating far from stability threshold and that are visible in the frequency spectrum as higher harmonics of the fundamental oscillation. On the other hand, as shown in Figure 4.24b, a linear stability analysis performed using LURANS equations around the meanflow returns no unstable modes. In this chapter, we will finally tackle the nonlinear and linear analysis of under-expanded nozzle flows. First, a description of the computational setup adopted will be provided, including reference quantities for dimensionless N-S equations, meshing strategies, boundary conditions as well initialisation. Secondly, the above strategies will be validated on a steady RANS calculation that will be used to initialise the unsteady DDES calculations. In section 5.3, the TIC nozzle considered, experiencing a FSS unsteadiness and operating at 3 different jet Mach numbers M j = [1.83, 2.09, 2.27] (corresponding to N P R = [6, 9, 12]), will be analysed by means of nonlinear DDES calculations. The results will be compared with avalaible experimental data obtained at Pprime facility and published by Jaunet et al. (2017) [START_REF] Jaunet | Wall pressure and external velocity field relation in over-expanded supersonic jets[END_REF]. The section 5.5 will be dedicated to the presentation of Global stability Analysis for the case at M j = 2.09. First, a linear analysis performed around the steady RANS solution will be tackled and second a more consistent linear stability analysis around the DDES meanflow will be 5.1. NUMERICAL STRATEGY detailed. Conclusions will be drawn by the analysis of the stability spectrum and the spatial shape of an unstable global mode compatible in terms of frequency and symmetry with nonlinear calculations.

α ∞ (σ, St) BF (σ, St) M F (σ, St) N L ∆St BF % ∆St M F % 4.

Numerical Strategy

The nozzle geometry here considered is a Truncated Ideal Contour (TIC), obtained by means of

2D axisymmetric Method of characteristics technique with boundary layer correction and based on

Delery's strategies [START_REF] Délery | Aérodynamique interne, tuyères et arrières corps[END_REF]. The nozzle is finally truncated by imposing a wall exit angle θ = 3 • . The resulting shape is shown in Figure 5.1 and geometrical properties are listed in Table 5.1.

design Mach M d = 3.5 throat radius R t exit radius R e /R t = 2.55 divergent length L/R t = 9.63 convergent length l/R t = 3.68 curvature radius ρ = R c /R t = 10 wall exit angle θ = 3 • Table 5.1: TIC nozzle features
The reference quantities chosen are the free-stream stagnation ones (a * 0,∞ , ρ * 0,∞ , T * 0,∞ ), where the superscript * denotes dimensional quantities. Under the assumption of isentropically resting flow at free stream conditions, T * 0,∞ and p * 0,∞ are set equal to nozzle combustion chamber temperature (260 K) and standard sea-level pressure (0.1 MPa) respectively, while the nozzle throat radius (R t ) is chosen as reference length L ref .

The Reynolds number, based on the free-stream stagnation quantities, is therefore defined as

STEADY RANS SOLUTION

Re ∞ 0 = ρ * 0,∞ a * 0,∞ L ref µ * 0,∞ = 5.0 × 10 5 .
(

The three-dimensional computational domain includes the nozzle and the farfield region, discretised with 9 blocks. The nozzle and the external region obtained from its extrusion until the outlet is made Concerning the boundary conditions, wall function boundary condition with adiabatic temperature and pressure extrapolation are prescribed at nozzle walls (see section 3.6.3 for further details). Nozzle inflow is assumed to be subsonic and primitive variables (ρ, p, u) are found by coupling the linearised outgoing characteristic equation, 1D isentropic relations and state equation with imposed total temperature T 0 = T 0,∞ , total pressure p 0 = N P Rp a and flow angle α = 0 [START_REF] Whitham | Linear and nonlinear waves[END_REF][START_REF] Blazek | Computational fluid dynamics: principles and applications[END_REF] (see section 3.6.2 for further details). The resulting system of equations is iteratively solved by a Newton method. Outlet is assumed to be subsonic with prescribed pressure p a = p 0,∞ , while farfield lateral boundaries and inflow are treated with a non-reflective boundary condition [START_REF] Whitfield | Three-dimensional unsteady euler equations solution using flux vector splitting[END_REF] at imposed free stream Mach number M ∞ = 0.05 (see section 3.6.1 for further details).

Steady RANS solution

Grid resolution is based on a convergence study for the separation location predicted by RANS simulations. These steady solutions, which were computed using a Local Time Stepping (LTS) method, are shown in Figure 5.4. The resulting shock system appears to be composed of a first Mach disk (located at x d1 ), originating from the oblique shock coming from the separation point, and a series

NPR d = 1 + γ -1 2 M ideal d 2 γ (γ-1) = 18
(5.2)

x/y t As mentioned in section 3.4.3, LTS method suffers from poor convergence especially in RANS calculations with highly refined grids. Just to give an idea to the reader, the infinity norm residual Unfortunately, URANS methods did not permit to destabilise the shear layers and consequently no turbulent eddies were detected. As a consequence, unsteadiness could not be triggered. The reason for is the excessive turbulent viscosity produced as shown in Figure 5.6. In fact, previous calculations 5. [START_REF] Wired | We have ignition! nasa tests new rocket engine in mojave desert[END_REF]. UNSTEADY DDES SOLUTION performed on turbulent shear layers have shown that the turbulent to molecular viscosity ratio (µ T /µ) must be less than 1000. Otherwise, the development of the instabilities are damped in these regions.

It would interesting to investigate whether the same pathology is found with other turbulence models (e.g. κ -, κ -ω SST). In this work we choose to decrease the level of modelled turbulent viscosity, by adopting a DDES formulation for turbulence.

Statistical features

The flow topology of the mean pressure field obtained by time-averaging the unsteady DDES solution is shown in Figure 5.9. It seems that inter-Mach disks distances increases with the pressure ratio N P R, so that at high N P Rs the shock waves-expansion fans pattern is more stretched. The flow features of this mean flow are quite similar as those of the RANS solution, but if one observes in detail 5.3. UNSTEADY DDES SOLUTION Table 5.5 can easily realise that all characteristic positions are slightly shifted upstream. Actually, it is important to stress that mean flow computed is not a solution of neither N-S equations nor URANS equations, while the RANS solution is the base flow of URANS equations. Therefore, it is evident that, even though both these solutions appear globally quite similar, they are deeply different. The mean streamwise velocity field u is shown in Figure 5.10. Two different annular regions are distinguishable: the inner one limited by the internal shear layer, and the second one between the internal and the external shear layer. Particularly, the former appears to be strictly supersonic while the latter is subsonic just downstream a Mach disk and accelerates up to supersonic conditions just before the following Mach disk. This is possible thanks to external annular region that, "bending" towards the internal one because of expansion fans, reduces the inner region flow section. As one can state by observing this figure with special attention, the mean flow seems to be axisymmetric but not perfectly: this happens because, as we will see further, the unsteady solution is very rich in terms of frequency spectrum with a large amount of energy localised at low frequency. Consequently, since a mean flow is statistically converged only when all its frequency content is converged, the higher is the energy present at low frequencies, the higher is the time necessary to converge this solution. For this reason, in case of highly unsteady flows, the time required is so huge that it is even not affordable in terms of computational resources. For information, the signal length over which the unsteady solution was averaged is listed in Table 5.5 in terms of periods of a single-frequency wave oscillating at 2 kHz † . Just to give an idea to the reader, computing one period of such frequency wave, with a physical time-step ∆t = 2.35 × 10 -6 s, takes 16 hours on 308 processors of Occigen (CINES supercomputer).

Streamlines are added to the mean streamwise velocity field u in As for perturbations, Figure 5.13 gives interesting insights. Standard deviation (σ) shows a sharp peak localised at separation point due to shock oscillations as showed in previous works about nozzles experiencing a FSS regime [START_REF] Ostlund | Side-load phenomena in highly overexpanded rocket nozzles[END_REF][START_REF] Martelli | Flow dynamics and wall-pressure signatures in a high-reynolds-number overexpanded nozzle with free shock separation[END_REF]. The pressure signal at the streamwise location of this peak is As one can state observing Figure 5.14c, at N P R = 6 and N P R = 9 pressure jumps between two levels, one corresponding to a point located outside the separated region and the other inside. This is the result of separated region intermittency, already observed early by [START_REF] Kistler | Fluctuating wall pressure under a separated supersonic flow[END_REF] [START_REF] Kistler | Fluctuating wall pressure under a separated supersonic flow[END_REF] for supersonic flows over a forward-facing step (see Figure 5.14a) and later by Ostlund et al. (2004) [START_REF] Ostlund | Side-load phenomena in highly overexpanded rocket nozzles[END_REF] (see Figure 5.14b) in nozzle flows experiencing a FSS regime. On the other hand, at N P R = 12 this behaviour seems to be not present and replaced by a simpler harmonic oscillation.

UNSTEADY DDES SOLUTION

After the sharp peak, pressure RMS attenuates at nozzle exit, as already seen in the case of airfoil buffeting (see Figure 4.18b). By observing Figure 5.13, this attenuation seems to occur increasingly monotonic at higher NPRs. Nevertheless, while the DDES calculations for N P R = 9, 12 match quite well with experiments, the numerical simulation for N P R = 6 shows a greater unsteadiness at exit compared with the corresponding experimental data.

Spectral Analysis

Analysis of unsteady wall pressure perturbation p is performed by means of Power Spectral Density (PSD), referred here with pp * . Frequency is expressed in terms of Strouhal number

(St = f D j /U j ),
which is based on the fully expanded jet velocity (U j ) and the fully expanded jet diameter (D j ), given by

U j = γRT 0 M j 1 + γ -1 2 M 2 j , D j = D e    1 + γ -1 2 M 2 j 1 + γ -1 2 M 2 d    γ+1 4(γ-1) (5.4)
where T 0 is the stagnation temperature, D e is the nozzle exit diameter and M j is the fully expanded jet Mach number defined as

M j = NPR (γ-1)/γ -1 2 (γ -1) 1/2 .
(5.5)

In order to better recognise the dominant frequency ranges, PSDs are shown in a pre-multiplied form and normalised by the dynamic pressure of the fully expanded jet γP j M 2 j , with P j = p a . PSDs are computed using Welch's method [START_REF] Solomon | Psd computations using welch's method[END_REF], i.e. subdividing the pressure signal into N block blocks with 50 % overlap that are then individually Fourier-transformed. The PSD spectrum is finally obtained by averaging the DFT energy spectra of the various segments, in order to minimise the variance of the power spectral density estimator. The PSD parameters adopted are listed in Table 5 The question that now raised is the following: could this intense signature in the frequency-azimuthal wavenumber plane be the trace of a global instability? In the next chapter, we will try to give an answer to such a question. The result here shown seems to confirm this tendency, with a linear stability analysis returning global modes, which are consistent in terms of spatial shape and frequency, but stable.

Therefore, in order to resolve at least a part of turbulent spectrum, a partial resolved turbulence approach is required. In this way, the spectral filtering of turbulence energy, provided by a LES approach, permits to separate the large scale fluctuating part from the small one, and modelling only § In the case of airfoil buffeting, the low-frequency self-sustained oscillation observed (f ≈ 100 Hz) is two order of magnitudes lower than turbulent scales (f ≈ 10 kHz) 5.5. LINEAR STABILITY the latter. In this context, the cut-off wave number k c is given by the grid size ∆, as k c = π/∆ and the resulting filtered LES equations have the same mathematical form of URANS equations. In this view, one can consider to develop the turbulent field φ into several ranks of fluctuating parts φ

[k] as φ = φ + ∞ n=1 φ [k n-1 ,kn] (5.6)
and each one defined with the following partial integration of their generalised Fourier transform

φ [k n-1 ,kn] (ξ) = kn k n-1 φ e ikξ dk, (5.7)
where φ is the Fourier transform of the spatially filtered variable φ (ξ) and k n is a series of partitioning wave numbers. While URANS approach aims at modelling the entire turbulent spectrum φ [0,∞] , a LES approach only the smaller scales φ [kc,∞] . In this scenario, a DDES shifts from URANS to LES, when there is a sufficient grid refinement (∆ << d w ), permitting to resolve the large turbulent scales.

Nevertheless, it should be stressed that when turbulence is fully resolved (i.e. DNS), or only partially modelled (i.e. LES, DDES), base flow stays laminar, while the corresponding unsteady solution lies in the chaotic turbulent saddle. In such cases, even if obtaining a steady state solution could be possible, the latter would be inappropriate if one is interested in studying the turbulent fluctuation dynamics around it. In this context, mean flow represents the most coherent statistical object living in the same saddle of turbulent perturbations. In this scenario, could a global stability analysis, performed around the DDES mean flow, return any marginally-stable or even unstable global mode compatible with the unsteadiness at St ≈ 0.2 and St ≈ 0.3?

Linear Stability around the meanflow

In contrast with laminar flow, it is commonly accepted that, the nonlinear effects cannot be so easily neglected when studying the dynamics of perturbations on top of otherwise turbulent flows. For this reason, we will not talk about linearised equations but rather about nonlinear disturbance ones.

In this framework, let us assume that the spatial structure of the mean flow q is known. If we consider a turbulent solution perturbed around the mean flow which is based on a Reynolds decomposition, starting from a compact form of Navier-Stokes equations (2.1) analogous to that written for URANS equations (2.44), it is possible to obtain the following equations for turbulent perturbations (5.8a) and mean flow (5.8b) ∂q ∂t + ∇ • L (q) + N (q, q) + L q + N q, q + N q , q + N q , q = 0 (5.8a) ∇ • L (q) + N (q, q) + N (q , q ) = 0, (

where N (q , q ) is the so-called Reynolds Stress tensor.

The NonLinear Disturbance Navier-Stokes (NLDNS) equations around the mean flow are obtained by solving the system of Eqs.(5.8a)-(5.8b) as ∂q ∂t + ∇ • L q + N q, q + N q , q = ∇ • N (q , q ) -N q , q , (5.9)

where q = ρ , m , (ρE) T represents the N-S conservative variables perturbations and the different terms can be explicitly expressed in the following way NLDNS equations around the mean flow q can be written in the following compact form ∂q ∂t -J q = ∇ • N (q , q ) -N q , q f , (

                           ∂ρ ∂t + ∂m j ∂x j = 0 ∂m i ∂t + ∂ ∂x j m j u i + m j u i + p δ ij + τ lam = ∂ ∂x j f ij -f ij ∂ (
where J is the Jacobian operator of compressible N-S equations about the mean flow and turbulent fluctuations act as an unsteady forcing term f on the momentum and energy equations.

The comparison between LURANS equations and NLDNS ones around the N-S mean flow permits to draw some important considerations. First, both are obtained from two different approaches for turbulence: URANS equations in (2.44) are obtained for a statistical average q of the turbulent solution (∂ q /∂t = 0), while Eq.(5.8b) is valid for the long term time-averaged solution q. Moreover, in the URANS approach, as described in the previous paragraph, a full modelling of the turbulent energy spectrum, expressed by the term N (q , q ) is adopted through the introduction of one or two additionally turbulence equations. In this way, the Reynolds Stress tensor is no more unknown and it is reintroduced in the flux vector, as it is for turbulent viscosity in the new extended state vector q = [ρ, ρu, ρE, ρν] T . These assumptions permits to derive a new set of equations (i.e. URANS), admitting a steady turbulent solution (i.e. RANS) which is a fixed point and the base flow of the equations. Anyway, the statistical representation of turbulent solution and the turbulent energy spectrum modelling only explains the steady term absence in f , while the unsteady term N (q , q ) is due to the nonlinear turbulent fluctuations, that in case of infinitesimal amplitude perturbations are neglected.

However, when dealing with complex and high turbulent flows, requiring a partially or full resolved approach for turbulence, the forcing term f exists but it is difficult or even unfaisable to take it in account because of the large computational cost demanded. In these cases, there are two alternative 5.5. LINEAR STABILITY options: assuming f negligible and evaluating this hypothesis a posteriori by solving LNS equations and Linear Stability around the turbulent mean flow [START_REF] Mettot | Quasi-laminar stability and sensitivity analyses for turbulent flows: prediction of low-frequency unsteadiness and passive control[END_REF] or considering this forcing f as the sum of harmonic components, because of the presence of N (q , q ), without modelling it [START_REF] Beneddine | Unsteady flow dynamics reconstruction from mean flow and point sensors: an experimental study[END_REF]. The second approach leads to Resolvent analysis and assumes turbulent fluctuations to be the optimal forcing in exciting convective (linearly globally stable) instabilities.

In the laminar flow past a cylinder bifurcating into a single frequency limit cycle (Barkley [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF]), a global stability analysis performed around the mean flow results in a better capturing of the unsteadiness if compared with an analogous one performed around the base flow. In this case, Eq.(5.15 comes to Eq.(2.47) and a marginally stable mode (σ = 0) with a frequency matching that of the limit cycle observed in nonlinear calculations. Such condition has been called by [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF] [70] Real Zero Imaginary Frequency (RZIF) and it has been proven to be always verified for monochromatic oscillations. Moreover, the resolvent operator is not defined for any Real Zero Imaginary Frequency (RZIF) [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF] global mode. In fact, as already mentioned by Beneddine et al. (2016) [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF] this means that, if a resolvent analysis is performed on a marginally stable case, the resolvent operator should display a theoretically infinity response near this frequency and a consequent clear separation of the singular values near this frequency. However, what happens in non-single frequency oscillations? The answer to this question has been provided by [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF] [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF]. In such case, a global stability analysis around the mean flow returns a non-RZIF mode, i.e. globally unstable and with a frequency not matching the nonlinear one. This discrepancy is proven to be larger and larger with the increasing of higher harmonics strength.

The experimental results carried out on the over-expanded nozzle flow considered in this work have shown a self-sustained unsteadiness which emerge over the entire spectrum without higher harmonics.

This observation, confirmed by nonlinear calculations performed by means of Delayed Detached Eddy Simulations (DDES), is in the author's opinion a reasonable assumption to neglect the forcing term f as a first approximation.

A global stability analysis around a time-averaged DDES solution at N P R = 9 is tackled. As mentioned in section 5.3.2, the frequency spectrum is so rich that statistically converging all the frequency content is too much expensive in terms of computational resources. Consequently, the lack of a perfect convergence in the mean flow computation can induce the occurrence of spurious modes in the linear stability analysis. For this reason, the long-term time averaged DDES solution has been also averaged In contrast with the stability spectrum for LURANS equations, that provided by stability analysis around the mean flow shows a series of unstable modes even at low-frequency. This unstable modes are probably due to the activity related to the low-frequency hump observed in the PSD spectra and required further analyses in the future. The first main conclusion that we can draw is that LURANS Stability is consistent with the URANS nonlinear dynamics. In fact, no global unstable modes exist as well no unsteadiness are detected in the URANS nonlinear calculation. Anyway, this result opens three alternative scenarios. According to the the first, the unsteadiness at St = 0.2 and St = 0.3 are related to two effectively global unstable modes, but which result stable because of a modelling deficiency in the turbulence models adopted. In this scenario, the eigenfrequency of these mode is not enough separated from the turbulent unresolved scales and turbulent viscosity, tempting to close this resolution gap, is overproduced and stabilises the underlying self-sustained mechanism. In this case, it would be interesting to investigate whether 5.6. PRELIMINARY CONCLUSIONS ON LINEAR DYNAMICS another turbulence model (e.g. κ -ω SST), maybe limiting the turbulent viscosity production, would permit to reproduce the above mentioned unsteadiness in nonlinear URANS calculation as well provide global unstable modes in LURANS stability.

Preliminary conclusions on linear dynamics

Alternatively, in a second scenario, these modes could be related to convective instabilities, which are linearly globally stable and need a forcing to be triggered. In this case, URANS, and thus also LURANS, methods are known not to permit the development of convective mechanisms but only self-sustained ones occurring at frequencies well-separated from the turbulent scales. In this scenario, a strategy based on the resolvent or receptivity analysis of the turbulent mean flow (even URANS) could be attractive tools. Finally, in the case of a third possible scenario, if the mechanisms involved is nonlinear, not even such tools could be useful to extract any information [START_REF] Sansica | Instability and low-frequency unsteadiness in a shockinduced laminar separation bubble[END_REF].

Concerning the interest in global stability analysis around the DDES mean flow, one can easily assume that, the mean flow being not so different from the base flow, as shown in Figure (5.5b) and (5.12b), such an analysis permits to get rid of any turbulent viscosity effect. In this way, a linear stability analysis using LNS equations around the mean flow takes in account turbulence only in terms of distortion. In this approach, a global unstable mode at St ≈ 0.2 is effectively provided and its spatial shape is very similar to that of the corresponding stable mode in LURANS spectrum. Nevertheless, the weakness of such an analysis lies in the absence of the turbulent fluctuation forcing term f in Eq. (5.15). That being said, considering this term necessarily leads to three possible options: exactly computing each term, consuming large-time and memory resources, and solving Eq.(5.15) without the possibility to extract any global mode, because the system of equations is not homogeneous; modelling this forcing term and thus slipping back to a URANS-like method, with consequently above mentioned drawbacks; considering f as an harmonic forcing and turning towards a resolvent analysis.

CONCLUSION

In order to demonstrate the existence or not of any unstable global mode related to such unsteadiness, a linear stability analysis for the nozzle flow at N P R = 9 is tackled. Two globally stable modes, compatible in terms of frequency and azimuthal symmetry, are present in the stability spectrum obtained from the linearised unsteady RANS equations. This result is consistent with the corresponding URANS nonlinear equations, for which no unsteadiness was reproduced. The fact that these modes are exactly localised at shear layers gives different scenarios about the understanding of this phenomenon.

As already mentioned, triggering the shear layer instabilities/unsteadiness plays a crucial role in reproducing correctly the nonlinear unsteady dynamics. First, the shear layer stabilisation could be due to the adopted turbulence model over-producing turbulent viscosity. In this case, another turbulence model could eventually provides a self-sustained dynamics. Second, these modes could be effectively linearly globally stable, because related to convective instabilities and particularly sensitive to (possibly very low amplitude) external perturbations. Finally, the mechanism involved could be simply nonlinear.

In order to reduce the effect of the modelled turbulent viscosity, as a first approximation, a quasilaminar linear stability analysis was performed around the DDES solution averaged both in time and in the azimuthal direction. A globally unstable mode at St ≈ 0.2 is effectively provided and its spatial shape is very similar to that of the corresponding stable mode in the eigenspectrum of the linearised URANS equations. This result seems to corroborate the hypothesis according to which the stability of the modes at St = 0.2 and St = 0.3 is linked to the amplitude of the turbulent viscosity provided by the Spalart-Alamaras model, even though the presence of turbulent fluctuation forcing term is required, in order to be mathematically well-posed. It should be stressed that injecting the turbulent forcing term would change the mathematical nature of the problem and into a linear affine problem for which the eigenvalues and eigenvectors are not defined.

The physical phenomenon related to these unsteadiness is not clear, even though the frequencies related to the peak at m = 1, at different values of N P R, follow the Tam's correlation for screech frequency, as already reported in Jaunet et al. 2017 [START_REF] Jaunet | Wall pressure and external velocity field relation in over-expanded supersonic jets[END_REF]. Nevertheless, the authors concluded that this mode cannot be due to a screech-like phenomenon because no acoustic contribution is observed.

In support of this assessment, the m = 1 global mode at St = 0.2 obtained in this work appears not to emit acoustic radiations. This is in stark contrast with the m = 2 mode at St = 0.3. This scenario seems to be in contrast with that proposed by Martelli et al. (2020) [START_REF] Martelli | Flow dynamics and wall-pressure signatures in a high-reynolds-number overexpanded nozzle with free shock separation[END_REF].

PERSPECTIVES

Perspectives

The results produced in this work are not entirely conclusive. They nonetheless provide several insights and raise numerous questions that need answering.

As for the nonlinear dynamics, 1a) analysing the nozzle flow at N P R = 9 with other turbulence models (e.g. κ -ω SST ) is a priority, in order to corroborate the hypothesis according to which URANS methods are not capable to describe this type of self-sustained dynamics.

2a) The nozzle flow at N P R = 12 exhibits flow features as well as frequency content very different when compared with the other cases (i.e. N P R = 6, N P R = 9). At this regime, the pressure gradient is greater than the other ones, requiring a greater protection for the boundary layer in order to avoid any LES filtering occurrence in this region. This is possible by increasing a particular coefficient of RANS/LES transitioning, which necessarily delays the shear layer triggering. Although the comparison with the corresponding experimental data gives confidence in the simulation, the very low activity recorded in the shear layers, if compared with the other N P R conditions, requires more investigations maybe adopting another grid, forcing the shear layers, or testing other RANS/LES treatments (e.g. Z-DES).

3a) In order to complete the nonlinear analysis, one could investigate other N P R conditions and particularly verify whether the condition of isentropic flow is reached at a value close to that predicted by the MOC analysis (N P R = 18, see Figure 5.3). In this kind of parametric study, one would analyse the behaviour of the frequency content at increasing value of N P R up to under-expanded conditions, when classical screech phenomena is supposed to settle.

4a) It would be interesting to study other nozzle contours (e.g. Thrust Optimised Parabolic, Thrust Optimised Contour), characterised by different shock system topologies.

As for the linear dynamics, 1b) for the nozzle flow at N P R = 9, it is mandatory to verify the consistency of the linearised URANS equations when a different turbulence model is adopted (e.g. κ -ω SST ). In case both 6.2. PERSPECTIVES URANS and LURANS finds a self-sustained mechanism at St = 0.2 and/or St = 0.3 compatible in terms of azimuthal symmetry, we could conclude that this unsteadiness is definitely related to a global linear instability. Otherwise, the mechanism is very likely due to a convective instability or a nonlinear mechanism; 2b) If LURANS returns globally stable modes at St = 0.2 and St = 0.3 with other turbulence models (and in accordance with URANS), it would be interesting to tackle a receptivity analysis (e.g. Pseudo-spectrum [START_REF] Trefethen | Spectra and pseudospectra: the behavior of nonnormal matrices and operators[END_REF], Resolvent) to assess the hypothesis that these modes are effectively globally linearly stable and require a forcing to be excited.

3b) It would be interesting to extend the same analysis done for the TIC nozzle flow at N P R = 9

to other N P R regimes and nozzle contour (e.g. TOP, TOC).

3. 1

 1 Arnoldi algorithm step-by-step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Flow parameters for the 2D quasi-incompressible cylinder flow. . . . . . . . . . . . . . 4.2 Numerical parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Base flow calculation: comparison in terms of Drag Coefficient (C D ) with calculations from Giannetti & Luchini (2007) [1] (GL) . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Comparisons between global linear stability analysis around the base flow (BF) and nonlinear calculations (NL). The Strouhal number and growth rate of the nonlinear calculations (NL) are also reported as well the percent error in Strouhal prediction (∆St % ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Comparisons between global linear stability analysis performed on base flow (BF) and meanflow (MF). The Strouhal number and growth rate of the nonlinear calculations (NL) are also reported as well the percent error in Strouhal computation (∆St % ). . . . 4.6 Flow parameters for the 2D transonic turbulent buffeting . . . . . . . . . . . . . . . . 4.7 Numerical parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Comparison in terms of buffet onset α crit and corresponding frequency f N L with previous numerical calculations and experiments. The percent error in frequency computation (∆f % ) with experiments (Jacquin et al. 2009

4. 9

 9 Comparisons between global linear stability analysis performed on base flow (BF) and meanflow (MF). The Strouhal number and growth rate of the nonlinear calculations (NL) are also reported as well the percent error in Strouhal computation (∆St % ). . . . LIST OF TABLES 5.1 TIC nozzle features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 RANS calculation: comparisons in terms of separation point (x s ), location of 1st Mach Disk (x d1 ), location of 2nd Mach Disk (x d2 ) and location of 3rd Mach Disk (x d3 ) at different NPR conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Numerical parameters for nonlinear unsteady calculations . . . . . . . . . . . . . . . . 5.4 RANS/LES transitioning parameters adopted . . . . . . . . . . . . . . . . . . . . . . . 5.5 DDES mean flow calculations: comparisons in terms of mean separation point (x s ) and mean locations of the first 3 Mach disks (x d1 , x d2 , x d3 ), at different NPR conditions. Inter-Mach disk distances (∆d) are also reported, as well the signal length in terms of periods of a wave signal at 2 kHz (T f =2 kHz ), corresponding to St = 0.2. L is the nozzle divergent length. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 PSD computation parameters. The frequency at 2 kHz corresponds to a Strouhal number St = 0.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1Figure 1 . 2 )

 12 Figure 1.2), the dotted line when a normal shock is expected at the nozzle exit (p usu pressure condition in Figure 1.2), and the chain-dashed line when the flow is adapted (p a,d pressure condition in Figure 1.2). . . . . . . . . . . . . . . . . . . . . . . . . . .

5. 7

 7 Instantaneous DDES solution, N P R = 9. Contours of turbulent to molecular viscosity ratio (µ T /µ). The maximu value of turbulent viscosity is localised at the separation point, with a lower value at the shear layers. . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Unsteady DDES solution. Contours of density gradient (||∇ρ||) at different N P R conditions: N P R = 6 (top), N P R = 9 (middle), N P R = 12 (bottom). Logarithmic scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Mean DDES solution at different NPR conditions: N P R = 6 (top), N P R = 9 (middle), N P R = 12 (bottom), longitudinal view. Contours of mean pressure p = p * /p * a . The colorscale is centered at p = p * /p a = 1 (i.e. where mean pressure p matches ambient pressure p a ) to highlight the shock system. . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Contours of the mean streamwise velocity u at different NPR conditions: N P R = 6 (top), N P R = 9 (middle), N P R = 12 (bottom). The continuous black and red lines are the sonic line and 0-value streamwise velocity isolines, respectively. The transverse slice are taken just downstream the first Mach disk. . . . . . . . . . . . . . . . . . . . LIST OF FIGURES 5.12 Comparison between experiments and DDES in terms of mean wall pressure along the same generatrix for different NPRs. Solid lines, DDES; filled (empty) circles, measurements from experimental campaign of 2016 (2019). Black, red and blue color stands for N P R = 6, 9, 12 data, respectively. p 0 is the combustion chamber total pressure and L the divergent length. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.13 Comparison between experiments and DDES in terms of wall pressure RMS along the same generatrix for different NPRs. Solid lines, DDES; filled (empty) circles, measurements from experimental campaign of 2016 (2019). Black, red and blue color stands for N P R = 6, 9, 12 data, respectively. p 0 is the combustion chamber total pressure and L the divergent length. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.14 Intermittency inside the separated region. . . . . . . . . . . . . . . . . . . . . . . . . . 5.15 Nozzle flow at N P R = 12. Zoom visualisation for the Premultiplied PSD spectrum of the main m azimuthal mode for the wall pressure perturbation p in Figure 5.19e. . . . 5.16 Premultiplied PSD of wall pressure perturbation p for different azimuth positions at a fixed streamwise location. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.17 Premultiplied PSD of wall pressure perturbation p for different streamwise locations at azimuth position θ = 0 • . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.18 Premultiplied PSD of wall pressure perturbation for different streamwise locations at azimuth position θ = 180 • . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.19 Premultiplied PSD of main m azimuthal mode wall pressure perturbation p m at a fixed streamwise position. Green, blue and yellow denotes spactra for m = 0, 1, 2 respectively. 5.20 Comparison bewteen Tam's correlation (1986) [12] and 1 st azimuthal mode peak frequency at different values of NPR (M j ) . . . . . . . . . . . . . . . . . . . . . . . . . . 5.21 Contours of premultiplied PSD for main azimuthal mode of wall pressure perturbation (St p m p * m ) in space and frequency domain. The black lines denote the separation point (SP ) and the location of the first Mach Disk (M D), respectively. . . . . . . . . . . . . LIST OF FIGURES 5.22 Nozzle Flow at N P R = 9. Eigenspectrum for Global Stability Analysis around the RANS base flow. Two stable modes are visible at St ≈ 0.2 and St ≈ 0.3, denoted with a red and a blue empty circle, respectively. . . . . . . . . . . . . . . . . . . . . . . . . 5.23 LURANS : three-dimensional view of the stable global mode at (σ, St) = (-0.052, 0.200).

5. 25

 25 LURANS stability: contours of the streamwise velocity mode component ûx . The global mode at (σ, St) = (-0.077, 0.298) seems to emanate some kind of acoustic radiation. . 5.26 Nozzle Flow at N P R = 9. Eigenspectrum for Global Stability Analysis around the mean flow. An unstable mode (σ, St) = (0.06, 0.22) matching the unsteadiness at St N L ≈ 0.2 is denoted with an empty red circle. Conversely, another mode at St = 0.3, which can be probably associated to the nonlinear activity with symmetry m = 2, is stable. . . . 5.27 LNS around DDES meanflow: three-dimensional view of the unstable global mode at (σ, St) = (0.06, 0.22). The iso-surfaces of the (a) streamwise, (b) radial and (a) azimuthal velocity mode component are plotted for the levels ûx = ±1.0 × 10 -4 , ûr = ±7.0 × 10 -5 and ûθ = ±7.0 × 10 -5 , respectively. White and black iso-surfaces for positive and negative values, respectively. The azimuthal m = 1 feature is well visible in the shape of the mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1: Shock diamonds structure during 5M15 engine firing up test in Mojave Desert. From Wired [3]

Figure 2 :

 2 Figure 2: Ariane 5 take off . From www.arianespace.com [4].

1

 1 equations and numerical methods previously described are validated on a laminar and turbulent flow, respectively in Chapter 4. In Chapter 5, RANS, URANS and DDES nonlinear calculations are shown in the case of over-expanded nozzle flows and compared with corresponding experiments. Second, a Global Stability Analysis around RANS base flow and DDES mean flow is described. In Chapter 6, general conclusions and perspectives are provided . Supersonic nozzles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.2 Geometry design for rocket nozzles . . . . . . . . . . . . . . . . . . . . . . . . 27 1.2.1 Truncated Ideal Contour (TIC) nozzles . . . . . . . . . . . . . . . . . . . . . . 28 1.2.2 Thrust Optimised Contour and Thrust Optimised Parabolic nozzles . . . . . . 28 1.3 Shock wave-boundary layer interactions in Supersonic nozzles . . . . . . . 30 1.3.1 Free Shock and Restricted Shock Separation in Over-expanded Nozzles . . . . 33 1.3.2 Self-sustained shock oscillation in divergent diffusers . . . . . . . . . . . . . . . 34 1.3.3 Screech tones in under-expanded nozzles . . . . . . . . . . . . . . . . . . . . . . 37 1.3.4 Transonic resonance in divergent diffusers and conical nozzles . . . . . . . . . . 39 1.3.5 Unsteadiness dynamics in modern over-expanded nozzles . . . . . . . . . . . . 39

Figure 1 . 1 :

 11 Figure 1.1: Nozzle flow regime during flight at p 0 = const with decreasing value of p a from bottom to top.

  d and p a > p a,d ) or during the engine start-up (p 0 < p 0,d and p a = p a,d ), they are usually over-expanded ‡ . For sake of simplicity, let us consider a fully operational engine (p 0 = p 0,d ) at different value of ambient static pressure p. A typical inviscid one-dimensional supersonic nozzle should behave as shown in Figure 1.2.

1. 2 .Figure 1 . 2 :

 212 Figure 1.2: Nozzle flow at different NPR = p 0 /p for a fully operational engine (p 0 = p 0,d ): adapted flow for NPR = p 0 /p a,d ; sonic flow throat pressure p * ; oblique shock at exit for NPR = p 0 /p usu ; supersonicsubsonic flow with inner or outer shock occurrence for p 0 /p l < NPR < p 0 /p usu ; pure subsonic flow with NPR < p 0 /p l . p l and p a,d are called limit pressure and adaptation pressure, respectively.

  Zoom of the characteristic net (i.e. kernel). Right and left characteristic curve are denoted with red and blue solid line, rispectively. The red filled circle point detected at the intersection of each couple of caracteristic is uniquely determinated.

Figure 1 . 3 :

 13 Figure 1.3: Ideal Contour nozzle with sharp corner obtained by a 2D Method of Characteristics (MOC), M e = 3.5.

1. 3 .

 3 SHOCK WAVE-BOUNDARY LAYER INTERACTIONS IN SUPERSONIC NOZZLES (a) Rapresentation of kernel and nozzle contour in TIC nozzle. (b) Rapresentation of kernel and nozzle contour in TOC nozzle. An internal shock is well visible.

Figure 1 . 4 :

 14 Figure 1.4: Main flow differences for isentropically inviscid calculation between TIC and TOC nozzles. From Nasuti & Onofri (2009)[START_REF] Nasuti | Shock structure in separated nozzle flows[END_REF]

1. 3 .Figure 1 . 5 :Figure 1 . 6 :

 31516 Figure 1.5: Phenomenology of shock separation. At the incipient separation point, wall pressure increases sharply from p i to a plateau pressure p p . After that, pressure slightly increases until exit up to p e . From Aghababaie & Theunissen (2015)[START_REF] Aghababaie | Modeling free shock separation induced side loads in overexpanded rocket nozzles[END_REF]

Figure 1 . 7 :

 17 Figure 1.7: Supersonic flow over forward-facing steps from Matheis & Hickel (2015)[START_REF] Matheis | On the transition between regular and irregular shock patterns of shock-wave/boundary-layer interactions[END_REF]. A restricted ricirculation region is visible bewteen separation point and the forward facing step (left). Near separation, wall pressure rapidly increases from p s to a plateau value p p and then more slowly once more up to p r where flow reattaches (right).

  start-up engine. In order to investigate the nature of this unsteadiness and to characterise flow separation,[START_REF] Nave | Sea level side loads in high-area-ratio rocket engines[END_REF] [START_REF] Nave | Sea level side loads in high-area-ratio rocket engines[END_REF] performed several experiments on full scale and cold-flow subscale models of a J-2S engine. The results of such experiments revealed for the first time the conical and highly unsteady pattern of two different separation: the Free Shock Separation (FSS) without flow reattachment and the Restricted Shock Separation (RSS) with reattachment. A phenomenological sketch of the flowfield and the corresponding wall pressure trend for both cases is shown in Figure1.8-1.9.

Figure 1 . 8 :

 18 Figure 1.8: Free-shock separation in overexpanded rocket nozzles, wall pressure profile, and phenomenology:--, compression waves/shock,and ---, boundary/shear-layer edge. From Frey et al. (2000)[START_REF] Frey | Restricted shock separation in rocket nozzles[END_REF]

1. 3 .

 3 SHOCK WAVE-BOUNDARY LAYER INTERACTIONS IN SUPERSONIC NOZZLES (a) Separation pressure ratio p min w /p a = NPR • p min w /p 0 trend with NPR (p c /p a = p 0 /p a ) for FSS and RSS regimes. From Nave & Coffey (1973) [21] Good job, you got it! (b) Wall pressure ratio (p w /p a ) trend with M j for a M d = 1.6 convergent-divergent nozzle. From Zaman et al. (2015) [22]. (c) Wall pressure ratio (p w /p a ) trend with M j for a M d = 2.8 convergent-divergent nozzle. From Zaman et al. (2015) [22].

Figure 1 . 10 :

 110 Figure 1.10: Comparison between a FSS-RSS transitioning nozzle (top) and two different design Mach number M d supersonic nozzles (bottom). Two different minimum separation pressure ratio trends with M j are visible.

Figure 1 . 11 :

 111 Figure 1.11: Flow separation at different values of M sw , by Sajben et al. (1981) [10].

  temperature. The above mentioned correlation for screech frequency with M j in terms of Strouhal number St is shown in Figure 1.13.

Figure 1 . 13 :

 113 Figure 1.13: Strouhal number (St) trend with jet Mach number (M j ) according to Tam's correlation(1986) [START_REF] Tam | Proposed relationship between boroadband shock associated noise and screech tones[END_REF] for screech frequency.

Figure 1 . 14 :

 114 Figure 1.14: Difference between screech and transonic tones, from Zaman et al. (2002). The dashed line represents the condition when the flow is just choked (p l pressure condition in Figure 1.2), the dotted line when a normal shock is expected at the nozzle exit (p usu pressure condition in Figure 1.2), and the chain-dashed line when the flow is adapted (p a,d pressure condition in Figure 1.2).

a

  Wall-Resolved Large Eddy Simulation (WRLES) for such turbulent flows (Reynolds number ranging between 10 5 and 10 7 ) is difficult or even impossible with the current technology. The compromise adopted by Olson & Lele (2013) [40] was lowering the Reynolds number of their LES calculation. A low-frequency shock oscillation has been detected at a frequency slightly different from experiments unlike shock position and amplitude oscillations, which change significantly as strongly dependent on incoming boundary layer. Two alternative solutions are the use of a Wall-Modelled Large Eddy Simulation (WMLES) [41] or a Detached Eddy Simulation (DES) first proposed by Spalart et al.

22 )

 22 by replacing the averaged flow primitive variables in Eq.(2.21) in Eqs. (2.1) and subsequently time averaging the latter, one can achieve the following Favre-and Reynolds-Averaged Navier-Stokes equations 2.2. COMPUTATIONAL APPROACHES TO TURBULENT FLOWS

(

  DES97), was made by[START_REF] Spalart | Comments on the feasibility of les for wings, and on a hybrid rans/les approach[END_REF] [START_REF] Spalart | Comments on the feasibility of les for wings, and on a hybrid rans/les approach[END_REF] who proposed to make the destruction term of the canonical Spalart-Allmaras model (2.27) grid-sensitive, by replacing the near wall distance d w with d as d = min (d w , C DES ∆) (2.33) with ∆ = max (∆x i ) and C DES a constant calibrated for the case of isotropic homogeneous turbulence and usually set at C DES = 0.65. The above equation permits to keep the RANS behaviour inside the 2.2. COMPUTATIONAL APPROACHES TO TURBULENT FLOWS boundary layer, where d w ∆, as well as acting as a Smagorinsky SGS model where d w ∆ if the grid is adequate. Despite its simple heuristic, DES97 nonethless has some major drawbacks, namely: -the incursion of LES mode inside the boundary layer; -the delay in the development of shear layer instabilities; -the erroneous activation of near-wall damping terms in LES mode. Delayed Detached Eddy Simulation (DDES) As mention in Spalart et al. (1997), the RANS preservation of the boundary layer is only assured with a grid size exceeding the boundary layer thickness δ. Unfortunately this can happen despite the user's foresight in the case of boundary layer thickening and separation, leading to the activation of the DES limiter d = C DES ∆ inside the boundary layer and a consequent eddy viscosity reduction [55]. This phenomenon, commonly known as MSD (Model Stress Depletion), has been shown by Menter and Kuntz (2004) [56] to be the cause of the further upstream displacement of the separation point predicted by RANS (using the same model) in the case of high-loaded airfoils. For this reason the authors named this phenomenon grid-induced separation. To protect the boundary layer from the incursion of LES mode, Spalart et al. (2006) [55] proposed a new version of DES called DDES (Delayed Detached Eddy Simulation) including a shielding function f d , which is supposed to be 0 in the RANS region and 1 in the LES one, reading

  Spalart et al. (2006) comes down to an enhanced version of the Extended Delayed Detached Eddy Simulation (EDDES) proposed by Riou et al. (2009)[START_REF] Riou | Improvement of delayed-detached eddy simulation applied to separated flow over missile fin[END_REF] which originally required ∆ = 3 Π 3 i=1 ∆x i . However as highlighted by Deck (2011)[START_REF] Deck | Recent improvements in the zonal detached eddy simulation (zdes) formulation[END_REF], this modification is not negligible because it makes ∆ time-flow dependent and not only grid sensitive.

40 )

 40 An analogous derivation for the near-wall corrected version of SA model with Edward's modification (SAE) has been shown by[START_REF] Mockett | A comprehensive study of detached-eddy simulation[END_REF] [START_REF] Mockett | A comprehensive study of detached-eddy simulation[END_REF] asΨ 2 SAE = min  100, max f v1 , 10 -10 max (χ, 10 -10 ) An alternative near-wall treatment, proposed by Chauvet et al. (2007), Riou et al. (2009) and Deck (2011), suggests to fully disable the activation of f v1 , f v2 and f t2 functions by imposing

Figure 2 .

 2 Figure 2.1).Numerous studies have been conducted on this particular flow configuration[START_REF] Jackson | A finite-element study of the onset of vortex shedding in flow past variously shaped bodies[END_REF][START_REF] Hill | A theoretical approach for analyzing the restabilization of wakes[END_REF][START_REF] Zebib | Stability of viscous flow past a circular cylinder[END_REF][START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF]. Close to Re = 50, the predictions of linear stability analyses are consistent with the experimental observations and numerical simulations : the flow pattern associated to the unstable pair of modes is representative of the von Kàrmàn vortex street and the predicted frequency is in good agreement with the

Figure 2 . 1 :

 21 Figure 2.1: Incompressible flow past a cylinder at Re = 55. Streamwise velocity contours.

  only the leading eigenmode of the linearized Navier-Stokes operator. The resulting system of equations (i.e. an algebraic nonlinear equation for the mean flow forced by the perturbation's Reynolds stresses and the linearized equations for said perturbation) are solved iteratively as the amplitude of the perturbation is increased. For small amplitudes, the "instantaneous" mean flow and associated perturbation reduce to the classical base flow and eigenmode. As the amplitude of the perturbation increases, the "instantaneous" mean flow slowly gets distorted into the classical mean flow while the leading eigenmode increasingly resembles the marginal eigenmode obtained from a stability analysis around the mean flow. Further theoretical and empirical justifications of mean flow stability analysis have been proposed over the years. For instance,[START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF] [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF] have shown, by means of a weakly nonlinear expansion about the stability threshold, that such an analysis is valid only if the 0-th (fundamental) harmonic is much stronger than the second harmonic. Otherwise, as for the case of open cavity flow, linear stability analysis around the mean flow returns an unstable mode with a frequency different from nonlinear flow. Turton et al. (2015)[START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF] have provided strong evidence that stability analysis around the mean flow is valid only for monochromatic (single-frequency) oscillations.

T

  represents the RANS conservative variables perturbations and J the Jacobian operator of compressible URANS equations about the base flow. It is worthwhile to 2.3. GLOBAL STABILITY FRAMEWORK remark that LURANS equations in Eq.(2.47) differ from Linearised Navier-Stokes (LNS) equations only for the presence of turbulent stress tensor τ turb ij perturbations in the momentum equation and produced by the corresponding linearised version of turbulence model. Further details will be given in section 3.2.1.

( 3 . 1 )

 31 or Eqs.(3.13) with their corresponding turbulence transport equations, for nonlinear and linear equations respectively. Particularly, in the cell-centred framework, three different approaches for the approximation of the discretised residual vector at each face of a cell are available: -averaging the fluxes computed from the cell-centred values to the left and to the right of the cell face; -averaging the flow variables computed from values at the centroids of the grid cells to the left and to the right of the cell face; -computing the fluxes from flow quantities interpolated separately to the left and to the right side of the cell face.

1 )-

 1 Λ c = u i dS i + a dS, spectral radius of the convective flux Jacobian ∂F c /∂q ‡ ;

  limiter ψ defined byψ(x, y) = minmod (x, y) = sign(x) • max [0, min (|x|, y • sign(x))] . (3.42)Like in the case of RANS equations, λ I represents the derivative of the convective flux F T c with respect to conservative variable q as 3

  by simply adding separate routines to a pre-existing Navier-Stokes code[START_REF] Liu | A strongly coupled time-marching method for solving the navier-stokes and κ-ω turbulence model equations with multigrid[END_REF]. Particularly, in this work the turbulence transport equation is marched one time step for each Navier-Stokes iteration[START_REF] Park | Implementation of κ -ω turbulence models in an implicit multigrid method[END_REF]. A Local Time Stepping (LTS) Technique or a Selective Frequency Damping (SFD)[START_REF] Åkervik | Steady solutions of the navier-stokes equations by selective frequency damping[END_REF][START_REF] Jordi | Encapsulated formulation of the selective frequency damping method[END_REF] has been used to solve steady problems, with the latter method implemented by the author. Conversely, a Dual ¶ In the linear solver, Λ T c I

  idea of Local Time Stepping (LTS) technique is advancing the solution in each control volume with a local maximum time step ∆t L obeying to a CFL condition imposed all over the computational domain. If we use an Euler implicit time-integration for the spatially discretised version of Eq.(3.1)

. 49 )

 49 After replacing ∂ R/∂q with ∂ R * /∂q in Eq.(3.47) and several mathematical manipulations, one can obtain the following tridiagonal system with the I direction terms for the i-th conservative variable3.4. TIME INTEGRATION

( 3 . 6 -q n 6

 366 55) is first discretised in space and then marched in time with frozen flow variables immediately after Eqs.(3.50)-(3.53) are solved [85]. The tridiagonal system deriving from the discretisation in space of Eq.(3.55) is solved by a Thomas Algorithm to obtain ∆q 6 = q n+1

  ones in Dual Time Stepping (DTS) framework as proposed by Luo et al. (2001) [90]. This approach, based on the following second-order time accurate version of Eq.(3.1)

. 59 )

 59 If Eqs.(3.58)-(3.59) are spatially discretised and iteratively solved by an Euler implicit scheme with a

Figure 3 . 1 :

 31 Figure 3.1: Matrix representation of the Arnoldi decomposition. From Bucci, Phd thesis (2017)[START_REF] Bucci | Subcritical and supercritical dynamics of incompressible flow over miniaturized roughness elements[END_REF]

  and r k = |βy k (m)| 3.6. BOUNDARY CONDITIONS

Contents 4 . 1 94 4. 1 . 1 95 4. 1 . 2 97 4. 1 . 3 103 4. 2 . 1 104 4. 2 . 2 106 4. 2 . 3

 41941195129713103211042210623 2D laminar flow past a cylinder . . . . . . . . . . . . . . . . . . . . . . . . . Base flow calculation: steady Navier-Stokes solution . . . . . . . . . . . . . . . Nonlinear unsteady calculation: unsteady Navier-Stokes solution . . . . . . . . Linear unsteady calculation: LNS solution . . . . . . . . . . . . . . . . . . . . . 98 4.1.4 Linear Stability around base flow . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.1.5 Linear Stability around mean flow . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.2 2D turbulent buffeting on OAT15A . . . . . . . . . . . . . . . . . . . . . . . Base flow calculation: Steady RANS solution . . . . . . . . . . . . . . . . . . . Nonlinear unsteady calculation: Unsteady RANS (URANS) solution . . . . . . Linear unsteady calculation: LURANS solution . . . . . . . . . . . . . . . . . . 109 4.2.4 Linear Stability around the base flow . . . . . . . . . . . . . . . . . . . . . . . . 110 4.2.5 Linear Stability around the meanflow . . . . . . . . . . . . . . . . . . . . . . . 111As previously mentioned, the code used in this work is Phoenix[START_REF] Goncalves | Numerical simulations of a transport-aircraft configuration[END_REF], developed at Pprime Institute in collaboration with Dynfluid Laboratory. Although the code had been already validated for 2D laminar flows (Guiho et al. 2016[START_REF] Guiho | Instabilities in oblique shock wave/laminar boundarylayer interactions[END_REF]) and for simple three-dimensional laminar flows (Sansica et al. 2018[START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and hopf bifurcations[END_REF]), during the last three years a large amount of work has been spent in terms of data communications among different processors, in order to tackle 3D turbulent flow past more complex geometries. These modifications included the numerical strategies in nonlinear and linear solver as well as the implementation of new methods for the search of fixed points (i.e. SFD). For this reason, a full validation of the new version of the code was required before tackling the main objective of the present thesis, i.e. the linear and nonlinear analysis of turbulent nozzle flows. The 2D cases of laminar flow past a cylinder and transonic turbulent flow past an airfoil have been selected to validate the code 4.1. 2D LAMINAR FLOW PAST A CYLINDER and they will be presented in this chapter. In the considered laminar flow, Navier-Stokes equations are solved, while, in the turbulent one, the former are coupled with a Spalart-Allamaras turbulence model in a common RANS/URANS framework, as described in sections 3.1-3.2 for linear and nonlinear calculations. A 3rd order DNC scheme with a JST Scalar Dissipation Scheme is adopted for N-S equations, while the SA turbulence model is discretised with a 2nd order TVD scheme with Harten's entropy correction (see sections 3.3.1-3.3.2 for more details). Particularly, in the turbulent case a Ducros shock sensor[START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF] is added to the standard pressure sensor and a Quadratic Constitutive Relation (QCR2000)[START_REF] Spalart | Strategies for turbulence modelling and simulations[END_REF] is adopted for stress tensor formulation instead of classical Boussinesq's approximation.Finally, an implicit loosely coupled algorithm is used to integrate the N-S equations and turbulence equation separately by a Dual-Time Stepping (DTS) Technique. Further details are present in section 3.4.2.

Figure 4 . 1 :≤ 10 - 8 .

 41108 Figure 4.1: Computational grid used in the Nonlinear, Linear and Stability calculations for the 2D laminar cylinder flow.

Figure 4 . 2 :

 42 Figure 4.2: Quasi-incompressible flow past a cylinder at Re = 55. Convergence history in base flow computing by Local Time Stepping method.

Figure 4 . 4 and

 44 Figure 4.4 and Table 4.3, respectively.

Figure 4 . 3 : 54 Table 4 . 3 :Figure 4 . 4 :

 43544344 Figure 4.3: Base flow calculation: streamwise velocity contours for Re = 55. White lines are isolines at null velocity and represents the boundary of the wake bubble.

Figure 4 . 5 :

 45 Figure 4.5: Quasi-incompressible flow past a cylinder at Re = 55. Growth-rate (σ N L ) and limit cycle frequency (St N L ) for a nonlinear calculation.

Figure 4 . 6 :

 46 Figure 4.6: Quasi-incompressible flow past a cylinder at Re = 46 and Re = 47. Lift coefficient trend with time for a stable (4.6a) and unstable (4.6b) baseflow.

( a )

 a Lift coefficient trend with time, for perturbation. The red dotted line slope represents the perturbation growth-rate (σ L ). (b) Amplitude spectrum. The only peak is located at St L = 0.120.

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Quasi-incompressible flow past a cylinder at Re = 55. Growth-rate (σ L ) and frequency (St L ) instability, for a linear calculation.

Figure 4 .

 4 Figure 4.10b). When the same calculation is carried out at different Reynolds numbers, growth rate is seen shifting from negative values, for stable baseflows, to positive ones, for unstable baseflows, as shown in Figure 4.11a. The critical value of Re for which the global mode becomes unstable is about Re c = 46.6 as it occurs in nonlinear flow. Instabilities frequency grows with increasing value of Re as reported in Figure 4.11b. Both σ and St are in good agreement with previous equivalent calculations solved via incompressible solvers [1].

4. 1 .

 1 Figure 4.10: Quasi-incompressible flow past a cylinder at Re = 55. Linear Stability Analysis around base flow.

Figure 4 . 11 :

 411 Figure 4.11: Linear Stability Analysis around base flow. Red triangles correspond to the current simulation, black empty circles to data from Giannetti & Luchini (2007) [1].

% 60 ( 4 Table 4 . 5 :

 60445 0.046,0.120) (0.000,0.134) (0.045,0.135) 11.1 0.7 70 (0.072,0.120) (0.000,0.143) (0.070,0.145) 16.1 1.Comparisons between global linear stability analysis performed on base flow (BF) and meanflow (MF). The Strouhal number and growth rate of the nonlinear calculations (NL) are also reported as well the percent error in Strouhal computation (∆St % ).

  (a) Linear Stability Analysis around base flow. (b) Linear Stability Analysis around meanflow

Figure 4 . 12 :

 412 Figure 4.12: Quasi-incompressible flow past a cylinder at Re = 60. Comparison between Linear Stability Analysis around baseflow and around meanflow.

Figure 4 . 13 .Figure 4 . 13 :

 413413 Figure 4.13: Computational grid used in the Nonlinear, Linear and Stability calculations for the 2D turbulent OAT15A airfoil flow.

Eqs. ( 3 .

 3 [START_REF] Beneddine | Global stability analysis of underexpanded screeching jets[END_REF], ∆ and χ are computed and SFD is activated on the saturated unsteady solution. A base flow converged with a precision of 10 -8 is obtain in about 250, 000 Dual Time Step (DTS) iterations and it is shown in Figure4.15. As shown in Figure4.14b, SFD permits to converge more rapidly than LTS during the first iterations and although after this first residual fall down convergence speed is the 4.2. 2D TURBULENT BUFFETING ON OAT15A

  (b) Residuals trend in time for LTS and SFD methods in logarithmic scale. The black and red line correspond to the action of LTS and SFD, respectively. you found my trick! Good job! I know, I know, sorry That 's the last time I joke, I promise

Figure 4 . 14 :

 414 Figure 4.14: Transonic turbulent flow past an OAT15A airfoil at α ∞ = 4.5 • . SFD action under a flow-unleash method is used to get the steady solution.

Figure 4 . 15 :

 415 Figure 4.15: Base flow calculation for a 2D transonic turbulent flow past an OAT15A airfoil at α ∞ = 4.5 • : streamwise velocity contours. The solid black line is the sonic line.

Figure 4 .

 4 [START_REF] Rao | Approximation of optimum thrust nozzle contours[END_REF]).(a) Lift coefficient trend with time, starting from base flow. The solid red line slope represents the growth-rate (σ N L ). I need it, try to understand (b) DFT power spectrum computed on the saturated limit cycle signal. A peak at St N L = 0.0631 is well visible as its second and third harmonics.

Figure 4 . 16 :

 416 Figure 4.16: Transonic turbulent flow past an OAT15A airfoil at α ∞ = 4.5. Growth-rate (σ N L ) and limit cycle frequency (St N L ) for a nonlinear calculation.

4. 2 .Figure 4 . 17 :

 2417 Figure 4.17: Unsteady solution for a 2D transonic turbulent flow past an OAT15A airfoil at α ∞ = 4.5 • : instantaneous density gradient (||∇ρ||) contour in logarithmic scale. A strong lambda shock is well visible on the airfoil surface.

Figure 4 .

 4 18a with 4.15, the lambda shock is remarkably smoothed due to the averaging of unsteadiness in the shock foot. Inside the supersonic bubble, which appears more flattened along the vertical direction, it is possible to detect the upstream front of shock excursion, while the downstream one is located near the sonic line. Pressure Root Mean Square (RMS) σ has been computed run-time for different streamwise wall location x and compared with an analogous 2D URANS calculation done by[START_REF] Deck | Numerical simulation of transonic buffet over a supercritical airfoil[END_REF] inFigure 4.18b. A good agreement is found especially in the shock foot mean position, located where the continuous line reach its maximum. However, in terms of flow frequency response, comparing Figure 4.16b with Figure 4.5b, the limit cycle 4.2. 2D TURBULENT BUFFETING ON OAT15A (a) Meanflow streamwise velocity contours. The black line is the sonic line. Sorry, I need that Sorry, I need that Sorry, I need that It is more elegant, isn't it? (b) Comparison in terms of pressure RMS (σ) between current unsteady calculation (empty circles) and an analogous 2D URANS calculation (continuous line) performed by Deck (2005) [35] for different streamwise wall locations. c is the coord length.

Figure 4 . 18 :

 418 Figure 4.18: Transonic turbulent flow past an OAT15A aurfoil at α ∞ = 4.5 • . Statistical properties for 120 periods time-averaged unsteady solution.

Figure 4 .Figure 4 . 19 :Figure 4 . 20 :

 4419420 Figure 4.19: Unsteady solution for a 2D transonic turbulent flow past an OAT15A airfoil at α ∞ = 4.0 • : a peak at St N L = 0.0623 is well visible with no too much energetic harmonics.

Figure 4 . 21 :

 421 Figure 4.21: Transonic turbulent flow past an OAT15A aurfoil at α ∞ = 4.5 • . Consistency between linear and nonlinear calculation.

Figure 4 .

 4 Figure 4.22a) shows a single mode with positive-real part eigenvalue σ, which additionally imaginary part ω = 2πSt is compatible with the observed limit cycle frequency. The spatial structure of this unstable mode is localised in the shock and it is particularly important near the foot and in the reflected shock wave originating from the triple point (see Figure 4.22b).

4. 2 .Figure 4 . 22 :

 2422 Figure 4.22: Transonic turbulent flow past an OAT15A airfoil at α ∞ = 4.5 • . Linear Stability Analysis around base flow.

Figure 4 . 23 :

 423 Figure 4.23: Transonic turbulent flow past an OAT15A airfoil at α ∞ = 4.5 • . Comparison between Stability Analysis around base flow and meanflow. The latter using LNS equations.

Figure 4 . 24 : 5 Over-expanded nozzle dynamics Contents 5 . 1

 424551 Figure 4.24: Transonic turbulent flow past an OAT15A aurfoil at α ∞ = 4.5 • . Comparison between Stability Analysis around base flow and meanflow. The latter using LURANS equations.

Figure 5 . 1 :

 51 Figure 5.1: TIC nozzle axisymmetric geometry, M d = 3.5.

of 5

 5 blocks meshed together with an O-H topology to avoid any singularity about the longitudinal axis (see Figure 5.2c). The remaining farfield region is composed by four additional O-grid blocks. Domain boundaries are 4 and 2 times the overall nozzle length (L * = l + L) away from the nozzle exit in the streamwise and radial direction, respectively. For the O-H structured blocks of the domain, the grid resolution is n x = 813 in the longitudinal direction (697 of which are within the nozzle and 117 in the exterior part) and n r = 177 (n θ = 89) in the radial (azimuth) direction for each of the 4 blocks surrounding the core nozzle. The 4 O-grid blocks in the exterior part have 117, 309 and 352 points along the longitudinal, radial and azimuth directions, respectively. The resulting grid, consisting in 69 millions of cells, is shown in Figure 5.2.

5. 2 .Figure 5 . 2 :

 252 Figure 5.2: Computational grid used in the RANS, URANS and DDES calculations.

Figure 5 . 3 :

 53 Figure 5.3: Nozzle Contour. TIC nozzle geometry (continous line), Ideal Contour nozzle obtained by MOC (filled circles line), symmetry axis (dashed line). y t corresponds to the throat radius R t .

Figure 5 . 4 :

 54 Figure 5.4: RANS solution at different NPR conditions: N P R = 6 (top), N P R = 9 (middle), N P R = 12 (bottom). On the left, contours of static pressure in terms of ratio with static ambient pressure p = p * /p *a . The colorscale is centered at p = 1 (i.e. where pressure p * matches ambient pressure p a ) to highlight the shock system. x s is the separation point, while x d = (x d1 , x d2 , x d3 ) the locations of the 1st, 2nd and 3rd Mach Disk. All these characteristic positions are listed in Table5.2. On the right, contours of streamwise velocity u. The solid black line is the sonic line.

Figure 5 . 5 :

 55 Figure 5.5: Comparison between experiments and RANS calculations in terms of wall pressure and separation point location along the same generatrix for different NPRs. p 0 and L are the combustion chamber total pressure and the divergent length, respectively.

Figure 5 . 6 : 3 :Figure 5 . 7 :

 56357 Figure 5.6: RANS solution, N P R = 9. Contours of turbulent to molecular viscosity ratio (µ T /µ). Turbulent viscosity is too much high to permit the shear layers destabilisation.

Table 5 . 4 :

 54 Flowfield visualisationA numerical Schlieren is shown in Figure5.8. The flow at N P R = 6 and N P R = 9 shows quite similar features with instabilities apparently originating in the vicinity of the separation point, growing when advecting in the streamwise direction and eventually triggering where the reflected shock coming from the triple point hits the external shear layer. From this point to downstream, more and more turbulent eddies seem to detach and to propagate upstream, possibly closing the loop. Particularly, those structures detaching from the external shear layer near the second Mach disk seem to interact with the nozzle lip, from which some kind of radiation appears to emanate into the external environment. Moreover, compression waves similar to shocklets (already observed in Chen et al. (1994)[START_REF] Chen | Numerical investigation of separated nozzle flows[END_REF]) seem to radiate from the internal shear layer to the external one. Particularly, at N P R = 6 this phenomenon is well visible in the region between separation and the triple point.The case at N P R = 12 is quite different from the previous ones. The external shear layer seems to be less unsteady with no turbulent eddies developing inside the nozzle with most of turbulent motion being localised in the region near the second Mach disk. In this regard, as the previous regimes 123 5.3. UNSTEADY DDES SOLUTION examined, the first Mach disk seems to oscillate almost rigidly with the separation shock, in contrast with the more complex motion of the second Mach disk. Similarly, the internal shear layer seems to be significantly less unsteady than the external one, which is probably due to a slight delay in triggering the turbulent shear layers. In this regard, it should be stressed that, as mentioned in section 2.2.3, this is a very sensitive aspect of DDES formulation which depends on the local grid refinement (∆ max ), local vorticity (∆ ω ) and the distance of the closest cell to wall (d w ), with C DES the only easilytunable parameter. Particularly, the author observed that for low values of C DES , delay in triggering shear layer instabilities decreases but boundary layer becomes less and less protected from the LES incursion. In order to reduce this inconvenient, as proposed by Ashton (2017)[START_REF] Ashton | Recalibrating delayed detached-eddy simulation to eliminate modelled-stress depletion[END_REF], the classical shielding function f d is opportunely modified as followf d = 1 -tanh (c 1 r d ) 3(5.3)with c 1 and C DES values listed in Table5.4. RANS/LES transitioning parameters adopted Particularly, the RANS/LES parameters adopted for the case at N P R = 6 and N P R = 9 are the same used by Martelli et al. (2020)[START_REF] Martelli | Flow dynamics and wall-pressure signatures in a high-reynolds-number overexpanded nozzle with free shock separation[END_REF] in a similar DDES calculation. 5.3. UNSTEADY DDES SOLUTION (a) Unsteady visualisation for N P R = 6 (b) Unsteady visualisation for N P R = 9 (c) Unsteady visualisation for N P R = 12

Figure 5 . 8 :

 58 Figure 5.8: Unsteady DDES solution. Contours of density gradient (||∇ρ||) at different N P R conditions: N P R = 6 (top), N P R = 9 (middle), N P R = 12 (bottom). Logarithmic scale.

Figure 5 . 9 :

 59 Figure 5.9: Mean DDES solution at different NPR conditions: N P R = 6 (top), N P R = 9 (middle), N P R = 12 (bottom), longitudinal view. Contours of mean pressure p = p * /p * a . The colorscale is centered at p = p * /p a = 1 (i.e. where mean pressure p matches ambient pressure p a ) to highlight the shock system.

Figure 5 . 11 .

 511 Two different recirculation regions are well visible near the shock separation and the nozzle lip at N P R = 6 and N P R = 9, as already observed by Martelli et al. (2020) [45] in an analogous over-expanded TIC experiencing a FSS regime. Particularly, at N P R = 9 the recirculation region near the shock separation seems vanished, while at N P R = 12 that near the nozzle lip moved upstream. Moreover, in contrast with classical TIC flow topology, at this N P R condition a trapped vortex appears downstream the second † The frequency at 2 kHz corresponds to a Strouhal number St = 0.2.

5. 3 .Figure 5 . 10 :Figure 5 . 11 :Figure 5 . 12 :

 3510511512 Figure 5.10: Contours of the mean streamwise velocity u at different NPR conditions: N P R = 6 (top), N P R = 9 (middle), N P R = 12 (bottom). The continuous black and red lines are the sonic line and 0-value streamwise velocity isolines, respectively. The transverse slice are taken just downstream the first Mach disk.

( a )Figure 5 . 13 :

 a513 Figure 5.13: Comparison between experiments and DDES in terms of wall pressure RMS along the same generatrix for different NPRs. Solid lines, DDES; filled (empty) circles, measurements from experimental campaign of 2016 (2019). Black, red and blue color stands for N P R = 6, 9, 12 data, respectively. p 0 is the combustion chamber total pressure and L the divergent length.

5. 3 .Figure 5 . 14 :

 3514 Figure 5.14: Intermittency inside the separated region.

5. 3 .

 3 UNSTEADY DDES SOLUTION separation shock. The slight peak position shifting between simulations and experiments is very likely due to the fact that DDES predicts a separation point slightly upstream compared with experiments. By means of PSD, the unsteady pressure perturbation has revealed certain signatures in frequency. First Baars et al. (2012) and later Jaunet at al. (2017) for the same TIC nozzle here considered, have shown also a precise azimuth selection for the peaks present in the PSD spectrum. For this reason, a PSD of the first three pressure perturbation azimuthal modes is presented in Figure 5.19. For all cases, DDES reproduces the same distribution provided by experiments in terms of azimuth selection for each frequency component. Particularly, at N P R = 6 the component at St = 0.3 is associated almost exclusively to the mode m = 2, while the St = 0.2 component and the low-frequency hump (St ≈ 0.05) are related to the mode m = 1. A peak at St = 0.1 is also visible for the mode m = 0. The high-frequency hump (St ≈ 1) contributes in terms of energy to all azimuth modes, confirming to be related to turbulent structures. When N P R increases up to 9, the azimuth selection for all frequencies keeps constant and the St = 0.2 mode becomes dominant. The over-estimation of lowfrequency hump energy in DDES calculation is clustered in the m = 0 mode, confirming a probable deficiency in reproducing the axysimmetric time-behaviour in this case. The same problem is detected at N P R = 12: while experiments gives a low frequency peak energy almost homogeneously distributed over the first three azimuthal modes, DDES over-estimates the energy for the m = 0 mode. In any case, numerical simulation seems to globally reproduce the same mechanism provided by the experiments. Particularly, at N P R = 12 a very little peak between St = 0.15 and St = 0.2 with symmetry m = 1 is also captured by the numerical simulation (see Figure 5.15).

Figure 5 . 15 :

 515 Figure 5.15: Nozzle flow at N P R = 12. Zoom visualisation for the Premultiplied PSD spectrum of the main m azimuthal mode for the wall pressure perturbation p in Figure 5.19e.

5. 4 .Figure 5 . 16 :Figure 5 . 17 :Figure 5 . 18 :Figure 5 . 19 :Figure 5 . 20 : 12 Figure 5 . 21 :

 451651751851952012521 Figure 5.16: Premultiplied PSD of wall pressure perturbation p for different azimuth positions at a fixed streamwise location.

Figure 5 . 22 :

 522 Figure 5.22: Nozzle Flow at N P R = 9. Eigenspectrum for Global Stability Analysis around the RANS base flow. Two stable modes are visible at St ≈ 0.2 and St ≈ 0.3, denoted with a red and a blue empty circle, respectively.

  Figures 5.23-5.24.

  Figure 5.23: LURANS : three-dimensional view of the stable global mode at (σ, St) = (-0.052, 0.200). The iso-surfaces of the (a) streamwise, (b) radial and (a) azimuthal velocity mode component are plotted for the levels ûx = ±1.0 × 10 -6 , ûr = ±1.0 × 10 -7 and ûθ = ±1.0 × 10 -7 , respectively. White and black iso-surfaces for positive and negative values, respectively. The azimuthal m = 1 feature is well visible in the shape of the mode.

5. 5 .

 5 LINEAR STABILITY Both modes are localised at shear layers in a wave-packet shape and show a clear azimuthal symmetry. In fact, the mode at St ≈ 0.2 exhibits an azimuthal m = 1 symmetry, while that St ≈ 0.3 an m = 2 one, both of them perfectly in accordance with the peaks present in the PSD spectra related to nonlinear DDES (Figure 5.19c) and experiments (Figure 5.19d). This result is very interesting as it leads to the conclusion that these modes are already present in the LURANS stability spectrum, but they are stable. Moreover, Figure 5.25 seems to indicate that the physical nature of this two modes is different, with the global mode at St ≈ 0.3 emanating some kind of acoustic radiation in the exterior environment. (a) Longitudinal view of the stable global mode at (σ, St) = (-0.052, 0.200). (b) Longitudinal view of the stable global mode at (σ, St) = (-0.077, 0.298).

Figure 5 . 25 :

 525 Figure 5.25: LURANS stability: contours of the streamwise velocity mode component ûx . The global mode at (σ, St) = (-0.077, 0.298) seems to emanate some kind of acoustic radiation.

5. 5 .

 5 LINEAR STABILITY in the azimuth direction. The linearised N-S equations around the DDES mean flow are then solved forward in time with a sampling period ∆t s = 6.35 × 10 -5 s and a Krylov space dimension K = 85 was adopted. These choice of parameters ensures a 8 snapshots-discretisation par period and a signal length of approximately 11 periods of a single-frequency wave oscillating at 2 kHz (corresponding to St = 0.2). The stability spectrum is shown in Figure 5.22.

Figure 5 . 26 :Figure 5 . 27 :

 526527 Figure 5.26: Nozzle Flow at N P R = 9. Eigenspectrum for Global Stability Analysis around the mean flow. An unstable mode (σ, St) = (0.06, 0.22) matching the unsteadiness at St N L ≈ 0.2 is denoted with an empty red circle. Conversely, another mode at St = 0.3, which can be probably associated to the nonlinear activity with symmetry m = 2, is stable.

  The linear stability analyses performed first around the RANS base flow and secondly around the DDES mean flow have given interesting insights in the process of understanding of the unsteady activity observed in over-expanded TIC nozzle at N P R = 9.The azimuth-selected unsteadiness at St = 0.2 and St = 0.3 observed in DDES calculations and experiments are found both in the LURANS stability spectrum and LNS stability spectrum. Particularly, in the LURANS spectrum both global modes are stable, but are consistent in term of frequency and azimuthal symmetry with the unsteadiness expected. The global modes at St = 0.2 and St = 0.3 are localised at the inner and external shear layer in a packet shape. Conversely, for the linear stability analysis performed around the DDES mean flow, the global mode at St ≈ 0.2 is unstable but shows the same spatial structure of that provided by LURANS equations.

  According toDeck (2011) [START_REF] Deck | Recent improvements in the zonal detached eddy simulation (zdes) formulation[END_REF], though the shielding function f d aimed to protect the attached boundary layer from MSD, its combination with the use of ∆ max as subgrid scale has been proved to delay the formation of instabilities in free shear layers. A first solution has been proposed by[START_REF] Chauvet | Zonal detached eddy simulation of a controlled propulsive jet[END_REF][START_REF] Chauvet | Zonal detached eddy simulation of a controlled propulsive jet[END_REF] who introduced a flow-dependent length-scale ∆ ω

.

[START_REF] Deck | Numerical simulation of transonic buffet over a supercritical airfoil[END_REF]

In this strategy, r d acts as a sensor distinguishing the points inside from those outside the boundary layer, according on the values of local turbulent viscosity and velocity gradient norm. Particularly, the inner points are those with moderate velocity gradient values and must be treated in a RANS way, thus with a greater value of turbulent viscosity. Conversely, the points outside the boundary layer, treatable in a LES way, are those with high velocity gradients (e.g. shear layers) and thus with less levels of modelled turbulence (i.e. turbulent viscosity) required. Therefore, with function f d build to 2.2. COMPUTATIONAL APPROACHES TO TURBULENT FLOWS get a smooth transition between RANS regions and LES ones, d naturally switches between d w and C DES ∆ according to f d . Zonal Detached Eddy Simulation (ZDES) and Extended Delayed Eddy Simulations (EDDS)

Table 4 .

 4 4: Comparisons between global linear stability analysis around the base flow (BF) and nonlinear calculations (NL). The Strouhal number and growth rate of the nonlinear calculations (NL) are also reported as well the percent error in Strouhal prediction (∆St % )

	4.1. 2D LAMINAR FLOW PAST A CYLINDER
	Re	(σ, St) L BF	(σ, St) N L	∆St BF %
	48 (0.005,0.117) (0.005,0.120)	2.5
	50 (0.013,0.118) (0.013,0.123)	4.1
	55 (0.030,0.119) (0.030,0.130)	9.2
	60 (0.046,0.120) (0.045,0.135)	11.1

Table 4 .

 4

	6: Flow parameters for the 2D transonic turbulent buffeting
	physical time-step	0.05
	DTS iterations CFL 10
	DTS iterations	100
	JST k 2 coefficient	0.5
	JST k 4 coefficient	0.01
	Harten coefficient	0.05
	Table 4.7: Numerical parameters

Table 4 . 8

 48

: Comparison in terms of buffet onset α crit and corresponding frequency f N L with previous numerical calculations and experiments. The percent error in frequency computation (∆f % ) with experiments (Jacquin et al. 2009

[START_REF] Jacquin | Experimental study of shock oscillation over a transonic supercritical profile[END_REF]

) is also reported. 120 periods. As one can state by comparing

Table 4 .

 4 9: Comparisons between global linear stability analysis performed on base flow (BF) and meanflow (MF). The Strouhal number and growth rate of the nonlinear calculations (NL) are also reported as well the percent error in Strouhal computation (∆St % ).

	5 (0.054,0.055) (0.173,0.097) (0.065,0.062)	11.3	56.4

Table 5 .

 5 2: RANS calculation: comparisons in terms of separation point (x s ), location of 1st Mach Disk (x d1 ), location of 2nd Mach Disk (x d2 ) and location of 3rd Mach Disk (x d3 ) at different NPR conditions.||L|| ∞ = max ||L|| ρ ∞ , ||L|| ρu ∞ ,||L|| ρE ∞ dropped by three orders of magnitude down to 10 -1 in 30,000 iterations. The problem is that, in this case, computing 1,000 LTS iterations takes approximately 20 hours on 308 processors of Occigen (CINES supercomputer). Hence, obtaining a sufficiently converged base flow calculation ||L|| ∞ < 10 -5 is not feasible in terms of computational resources through this strategy. Nevertheless, even if only a residual ||L|| ∞ ≈ 10 -1 is obtained, a good agreement in terms of wall pressure trend (p w ) and separation point (x s ) with experiments is obtained, as shown in Figure

	5.2. STEADY RANS SOLUTION
	N P R x s /L x d1 /L x d2 /L x d3 /L
	6	0.332 0.524 0.846 1.080
	9	0.530 0.737 1.147 1.495
	12	0.643 0.901 1.417 2.024
	5.5.
	(a) Separation point location (x s ) at different	(b) Wall pressure trend at different NPRs.
	NPRs. Black stars, experiments; emply black	Solid lines, RANS; filled (empty) circles, mea-
	circles, RANS simulations.		surements from experimental campaign of 2016
	sorry, i need that		(2019). Black, red and blue color stands for
	sorry, i need that		N P R = 6, 9, 12 data, respectively.

  5.3. UNSTEADY DDES SOLUTIONN P R x s /L x d1 /L x d2 /L x d3 /L ∆d 1-2 /L ∆d 2-3 /L

					signal length
	6	0.287 0.493 0.877 1.161	0.384	0.284	68 T f =2 kHz
	9	0.426 0.669 1.116 1.526	0.447	0.410	78 T f =2 kHz
	12	0.576 0.805 1.324 1.869	0.519	0.545	66 T f =2 kHz

Table 5 .

 5 5: DDES mean flow calculations: comparisons in terms of mean separation point (x s ) and mean locations of the first 3 Mach disks (x d1 , x d2 , x d3 ), at different NPR conditions. Inter-Mach disk distances (∆d) are also reported, as well the signal length in terms of periods of a wave signal at 2 kHz (T f =2 kHz ), corresponding to St = 0.2. L is the nozzle divergent length.

Table 5 .

 5 .6.

	N P R N block N snapshots/block ∆f [Hz] N periods f =2 kHz
	6	6	1475	24	83
	9	6	1440	25	80
	12	4	1408	38	53

6: PSD computation parameters. The frequency at 2 kHz corresponds to a Strouhal number St = 0.2.

  + (ρE) u j + pu j + p u j -q lam -(u i τ lam + u i τ lam ) = g ij = (ρE) u j + p u j -u i τ lam + u i µ L + u i µ L

	5.5. LINEAR STABILITY
						q lam =	c p P r L	µ L	∂T ∂x j	+ µ L	∂T ∂x j	,	(5.12)
				f ij = m j u j -µ L	∂u j ∂x i	+	∂u i ∂x j		-	2 3	∂u k ∂x k	δ ij	and	(5.13)
															∂u j ∂x i	+	∂u i ∂x j	-	2 3	∂u k ∂x k	δ ij +	c p P r L	µ L	∂T ∂x j	. (5.14)
	ρE) ∂t	+	∂ ∂x j	ρEu j ∂ ∂x j	g ij -g ij
																	(5.10)
	with
			τ lam = µ L	∂u j ∂x i	+	∂u i ∂x j	-	2 3	∂u k ∂x k	δ ij + µ L	∂u j ∂x i	+	∂u i ∂x j	-	2 3	∂u k ∂x k	δ ij ,	(5.11)

† The gas constants are based on the following standard coefficients: R * ≈ 8.314 J K -1 mol -1 for ideal gas, M * = 0.028 96 kg mol -1 and c * p = 1005 J kg -1 K -1 for dry air at standard temperature and pressure (STP).

X X X X X X X X X X X X X X X X

2.3. GLOBAL STABILITY FRAMEWORK

† For the sake of simplicity, the superscript tot has been neglected for perturbations

† This choice is done under the assumption that limit cycle frequency is close to that of any perturbation around the base flow.

LIST OF FIGURES

UNSTEADY DDES SOLUTION

Spectra for different azimuthal positions at a fixed streamwise location are provided by Figure 5.16.

For the case at N P R = 6, DDES calculations match quite well experiments confirming the presence of a low-frequency hump (St < 0.1), two different peaks at St = 0.2 and St = 0.3, with the latter more energetic than the former, and a high-frequency hump around St ≈ 1. High-frequencies energy seems to be slightly over-estimated by numerics. Increasing NPR value up to 9, the peak at St = 0.2 becomes dominant compared to that at St = 0.3 and strongly emerges over the entire spectrum. Numerical calculations well reproduce this behaviour, except for the presence of an energetic low-frequency hump that is not captured experimentally. When NPR is increased up to 9, both experiments and numerical calculation present a drastic change in the frequency content: intermediate frequency peaks as well the high frequency hump disappear and all the pressure perturbation energy is localised at a low-frequency hump, which is non necessarily homogeneous along the azimuth direction. There is a discrepancy in terms of frequency and energy between DDES simulation and experiments, that in the author's opinion is mainly due to the shortness of numerical signal.

In order to analyse the unsteady behaviour of wall pressure perturbation along the nozzle contour, premultiplied spectra for different streamwise locations on two different generatrix are provided by Figure 5.17-5.18. At N P R = 6, both experiments and DDES show the peak at St = 0.3 to increasingly prevail that at St = 0.2 away from separation point, until the former merges with the high-frequency hump, which become more and more dominant downstream. Conversely, the low-frequency hump seems to gradually disappear. The over-estimation by numerical simulations of frequencies around St ≈ 1 seems to affect only turbulence structures near nozzle exit (x/L = 0.757 -0.852). As for N P R = 9, this numerical energy over-estimation seems to be lower. In this case, both intermediate frequency peaks (St = 0.2 -0.3) are dominant in the short-medium distance from separation point. Particularly, the energy of the peak at St = 0.2 increases in the streamwise direction except for near-separation and nozzle exit locations. Actually, at separation point a low-frequency hump dominates the spectrum, as provided by the experimental measurement at x/L = 0.481. However, while experiments observe an abrupt drop for the energy of this hump in the streamwise direction, DDES returns not so important energy loss for the same hump. On the other hand, energy in high-frequency hump increases in the streamwise direction both for experiments and numerical simulations and become dominant near nozzle exit. At N P R = 12, both DDES and experiments depict the low-frequency hump energy decreasing in the streamwise direction, giving birth to the suspect that flow motion be essentially due to the

PRELIMINARY CONCLUSIONS ON NONLINEAR DYNAMICS

The spatial and frequency evolution of PSD for each single azimuthal mode at different N P R conditions is shown in Figure 5.21. If one considers the evolution of the 0 th mode (m = 0) at increasing values of N P R, it seems that at N P R = 6 this mode appears mainly localised at low frequency (St ≈ 0.1) for locations ranging from the shock separation to the first Mach disk (x/L < 0.5) and at high frequency (St ≈ 1) for near nozzle exit positions, the latter related to developed turbulence. When N P R is increased up to 9, the same mode is energetic all over the nozzle at very low frequency (St ≈ 0.01), while at N P R = 12 it is localised near separation at moderate low frequency (St ≈ 0.1). The 1 st azimuthal mode (m = 1) appears as a signature peaking at St = 0.2 already at N P R = 6, exclusively for positions ranging from the shock separation to the first Mach disk, while moving in the streamwise direction it seems to greater and greater contribute to the energy of high frequency hump. At N P R = 9, the above mentioned signature spreads all over the nozzle domain, providing a lower contribution to the high frequency hump for near nozzle exit position. Finally, at N P R = 12 the persistent signature previously observed vanishes except for a weak energy peak at St ≈ 0.15, which is still visible all over the nozzle. At this condition, the mode is localised at moderate low-frequency in the region between separation and the first Mach disk. The 2 nd azimuthal mode (m = 2) exhibits a spatial distribution at different N P Rs similar to the 1 st one. Actually it emerges all over the domain already at N P R = 6 and its energy peak shifts towards the nozzle exit at N P R = 9. At N P R = 12, the frequency content changes as the 1 st mode and the signature previously observed vanishes. At this condition, the energy appears localised near the shock separation at 0.01 < St < 0.1.

In conclusion, the energy peak of PSD for the 1 st and 2 nd azimuthal modes at St = 0.2 and St = 0.3 appears to be located near 1st Mach disk position and somewhere between the 1st and 2nd Mach disks, respectively. Moreover, by observing Figure 5.20, if one compares the screech correlation formula proposed by [START_REF] Tam | Proposed relationship between boroadband shock associated noise and screech tones[END_REF] [START_REF] Tam | Proposed relationship between boroadband shock associated noise and screech tones[END_REF] with the frequency of the 1 st azimuthal mode peak at the different NPR (M j ) there is a good match. Nevertheless, as already told by Jaunet et al. (2017) [START_REF] Jaunet | Wall pressure and external velocity field relation in over-expanded supersonic jets[END_REF], even though this result could confirms that the 1 st azimuthal mode is possibly due to a screech-like phenomenon, further investigations regarding potential acoustic radiation should be carried out to assess this hypothesis.

Preliminary conclusions on nonlinear dynamics

In this section, the author has shown the capability of Delayed Detached Eddy Simulations (DDES) in numerically reproducing the behaviour of an over-expanded TIC nozzle experiencing a FSS regime at Cosimo TARSIA MORISCO Nonlinear dynamics and linear stability analysis of over-expanded nozzle flows Résumé : Les interactions onde de choc/couche limite à l'intérieur de tuyères sur-détendues produisent des grandes régions détachées entraînant des forces non-axisymétriques appelées charges latérales. Le mécanisme générateur (mettant en jeu le décollement, les couches de cisaillement et les disques de Mach) est auto-entretenu avec des échelles de temps et d'espace relativement bien séparées de celles de la turbulence. Dans de tels cas, une approche partiellement résolue pour la turbulence (e.g. DDES) est envisageable afin de réduire les coûts de calcul. Cette thèse vise à étudier la nature potentiellement globalement instable de cette instationarité à l'aide de simulations numériques hybrides dans le cas d'une tuyère sur-détendue et par une analyse de stabilité linéaire 3D. La géométrie considérée est une tuyère TIC, en régime free shock seperation (FSS) et fonctionnant à 3 differents nombres de Mach de jet Abstract : Shock wave/boundary layer interactions in over-expanded rocket nozzles are responsible for large detached regions resulting in non-axisymmetric forces called side-loads. The mechanism at stake is self-sustained and involves separation, shear layers and Mach disks. In such cases, an hybrid approach for turbulence is required to mitigate the computational cost. This thesis aims at investigating the possibly globally unstable nature of this unsteadiness by means of a Delayed Detached Eddy Simulations (DDES) on an over-expanded nozzle and comparing it with a fully-3D linear stability analysis. The geometry considered is a TIC nozzle, experiencing a FSS unsteadiness and operating at 3 different jet Mach number M j = [1.83, 2.09, 2.27]. Nonlinear calculations confirm the experimental outcomes: power spectral densities for wall perturbations is distributed over 2 peaks at intermediate frequency (St = 0.2 -0.3) and two humps at low-frequency (St < 1) and high frequency (St ≈ 1), respectively. Particularly, at M j = 1.83 the peak at St = 0.2 competes with that at St = 0.3, prevails on the latter at M j = 2.09 and finally vanishes as the other at M j = 2.27. A PSD computed for different azimuthal components of wall pressure perturbation show a clear azimuthal separation for all the contribution mentionned above. Particularly, while the peak at St = 0.3 has a double contribution m = 2 -3 at M j = 1.83 and exclusively a m = 2 symmetry at M j = 2.09, the peak at St = 0.2 has constantly a m = 1 symmetry, which behaves has a persistent signature inside the nozzle at M j = 2.09. Consequently, a global stability analysis is performed on the DDES mean flow at M j = 2.09. Such analysis returns an unstable mode at St = 0.2, characterized by a m = 1 azimuthal symmetry, which develops from the separation point and is localised at the external shear layer. Keywords : shock-wave, nozzles, FSS, DDES, global stability analysis, meanflow.