).

. These good agreements are explained by the fact that the base flow stability assumptions, namely the existence of a steady solution of the equations and the linear assumption, are perfectly relevant at criticality. However, in some situations, the transition from a steady to an unsteady state occurs much earlier than what is predicted by a linear base flow stability analysis. This is the case for a Poiseuille flow for instance, where the analysis predicts a critical Reynolds number Re c = 5772, while in experiments, the transition to a turbulent state have been observed as

Chapter 1 -Introduction

Stability analyses for the study of unsteady flows

Unsteady flows play a central role in a very large number of problems, ranging from natural phenomena to engineering applications. For instance, both military and civil aircrafts are using increasingly large and powerful engines. As a consequence, the resulting noise has become an important nuisance that scientists are continuously trying to reduce. A major part of this noise does not originate from the engine itself, but from the unsteady interaction between the exhaust jet and the ambient medium. Therefore, understanding the nature and the features of the unsteady dynamics of jets is of great interest, and a large body of work has been devoted to this question. the underlying physical mechanisms. Classically, the perturbations are assumed to be small, which yields the simplified linear problem detailed in section 1.2. This linear assumption may be justified in some situations, such as in flow control for instance, where one may try to reduce the perturbation amplitude as much as possible. But its physical relevancy may be questioned in other contexts. Neglecting the nonlinear nature of the Navier-Stokes equations yields results that might not correspond to any observed physical behavior. Yet, surprisingly, in several situations, such a linear stability analysis accurately predicts nonlinear dynamical features. The underlying reasons of success remain a rather eluding question, for which we only have a partial understanding. Through the study of several different cases, this thesis aims at giving a clearer picture of different linear stability approaches, their physical relevancy, the underlying assumptions and the information they are able to give about the dynamics of flows.

Stability about the base flow

Classically, a linear stability analysis consists of considering a base flow -a steady solution of the governing equations -and then studying the behavior of small disturbances about this field. Such a base flow stability analysis can be either global, by studying perturbations over a whole domain, or local, by focusing on a cross-stream profile and see if locally, some disturbance waves may get amplified. The following sections successively describe global and local approaches, and present several classical stability theories.

Global approaches

Global base flow stability analysis for self-sustained unsteadiness The global base flow theory has been first introduced by [START_REF] Pierrehumbert | The two-and three-dimensional instabilities of a spatially periodic shear layer[END_REF] for the study of a spatially periodic shear layer. Since then, this has become a widely used approach to address the stability of a large range of flows (see for instance the review of [START_REF] Theofilis | Global linear instability[END_REF]). The main idea is as follows: let us consider a general flow governed by the Navier-Stokes equations, which can be recast, for instance in a compressible framework, in the following compact form: ∂q ∂t = N (q), (1.1) with q a state vector which contains the set of conservative variables describing the flow, and N the Navier-Stokes operator. A classical linear stability analysis consists of assuming the existence of a steady state q b (the base flow, verifying N (q b ) = 0), and then studying the temporal evolution of perturbations q ′ about this base flow. Assuming that these perturbations are small, q ′ is governed by a linear equation that reads:

∂q ′ ∂t = L(q b )q ′ , (1.2)
with L(q b ) = ∂N ∂q q b the linearized Navier-Stokes operator about the base flow, also called the Jacobian of the system. A global linear stability analysis consists of searching modes of the form q ′ = qe λt , which reduces equation (1.2) to the eigenproblem L(q b )q = λq. The corresponding eigenvalues characterize the stability of q b : if an eigenvalue λ i has a positive real part σ i = ℜ(λ i ) (positive growth rate), then the base flow is unstable and will be subject, in the linear regime, to self-sustained oscillations growing exponentially over the whole domain, and oscillating at a frequency ω i , given by the imaginary part of λ i (ω i = ℑ(λ i )). The spatial structure of these oscillations is then described by the eigenvector qi .

Resolvent analysis for the study of pseudo-resonance Another global approach, known as a resolvent analysis, consists of studying the effect of an external forcing on a base flow instead of focusing on self-sustained oscillations. This was first introduced by [START_REF] Trefethen | Hydrodynamic stability without eigenvalues[END_REF] to explain the unsteady behavior that was sometimes observed for stable flows (for instance, a Poiseuille flow is known to be unsteady even for a sub-critical Reynolds number, see section 1.2.3). This observed unsteady state relates to the non-normality of the linearized operator L, which allows external perturbations to get amplified by very large factors through non-modal linear mechanisms, even if all the modes are stable (see [START_REF] Trefethen | Hydrodynamic stability without eigenvalues[END_REF], [START_REF] Schmid | Stability and transition in shear flows[END_REF]).

Resolvent analyses were introduced to address this behavior, and focus on the evolution of small disturbances under an external forcing. This requires to add a forcing term f to equation (1.2), which yields:

∂q ′ ∂t = L(q b )q ′ + f .

(1.3)

In the frequency domain, this equation may be recast as:

q = R f , (1.4)
with R = (iωI -L(q b )) -1 the so-called resolvent operator, I the identity operator, ω the frequency, q and f the Fourier transform of q ′ and f , respectively. Then, a resolvent analysis aims at finding the optimal forcing which maximizes the energy gain G( f ) = ∥R f ∥2 /∥ f ∥ 2 , with ∥ • ∥ a relevant Hermitian norm 1 , defined from an inner product denoted as ⟨•, •⟩. This is related to the Schmidt decomposition of R, that reads:

∀ f , R f = i∈N µ i ⟨ϕ i , f ⟩ψ i , ⟨ϕ i , ϕ j ⟩ = ⟨ψ i , ψ j ⟩ = δ ij , µ i ≥ µ i+1 , (1.5)
where the sets (ϕ i ) i∈N , (ψ i ) i∈N and (µ 2 i ) i∈N are the so-called optimal forcings (of unit norm), optimal responses (of unit norm) and optimal gains, respectively. These quantities are linked by the relation:

µ i ψ i = Rϕ i , (1.6)
with ϕ i the i -th optimal forcing that maximizes the energy gain G, equal to µ 2 i . The optimal forcings (ψ i ) i∈N and optimal gains (µ 2 i ) i∈N are given by the eigenvectors and eigenvalues of R † R, respectively, with R † the Hermitian adjoint of R, defined as the operator verifying ∀a, b: ⟨a, Rb⟩ = ⟨R † a, b⟩. The optimal responses (ψ i ) i∈N correspond to the eigenvectors of RR † , and can be deduced from (ϕ i , µ i ) i∈N with equation (1.6).

In a discrete framework, R would be a discrete matrix of size n, and the discrete counterpart of the Schmidt decomposition is called the singular value decomposition (s.v.d.) of R. The sets (ϕ i ) 0<i≤n , (ψ i ) 0<i≤n and (µ i ) 0<i≤n are referred to as left singular vectors, right singular vectors and singular values of R, respectively. More details about the implementation of a resolvent analysis are given in chapter 2. Fig. 1.2 Schematic evolution of a perturbation close to x = 0 in the space-time plane (x, t), for a one-dimensional base flow with advection from left to right, in the case where the base flow is (a): absolutely unstable, (b): convectively unstable.

Local approaches

Local stability analysis A local stability analysis is conceptually close to a global analysis, but focuses on infinitesimal disturbances superimposed to a given crosswise profile instead of the complete base flow. Similarly to a global stability analysis, it yields an eigenproblem, whose size is however drastically smaller than its global counterpart.

Historically, this local approach has been introduced prior to the global methods, and early local stability analyses can be found, for instance, in the works of [START_REF] Lin | On the stability of two-dimensional parallel flows[END_REF], [START_REF] Batchelor | Analysis of the stability of axisymmetric jets[END_REF], [START_REF] Watson | On spatially-growing finite disturbances in plane poiseuille flow[END_REF] or [START_REF] Crighton | Stability of slowly diverging jet flow[END_REF]. The full procedure and main equations are detailed in chapter 2, and the central assumption is that the flow is approximately parallel, such that all the streamwise derivatives of base flow quantities could be neglected. Thus, only slowly diverging flows may be treated by such an analysis.

Two types of locally unstable behaviors may then be observed, denoted as convective and absolute local instability, respectively [START_REF] Huerre | Absolute and convective instabilities in free shear layers[END_REF]. An absolute instability describes disturbances that would eventually grow over a whole region, and never get convected out (see figure 1.2(a)). Hence, an initial impulse at a given location may yield self-sustained oscillations over this region, and this situation relates to a global instability. On the other hand, a convective instability describes a situation where local disturbances grow exponentially while traveling in the fluid, but they never contaminate the full domain due to advection (see figure 1.2(b)). At first sight, convectively unstable flows could be considered as globally stable configurations, whose unsteadiness stems from non-modal mechanisms, that could be studied with a resolvent analysis. But that is not always true: some global instabilities are the result of an acoustic feedback loop between an upstream and a downstream location (more details can be found in chapter 4). Due to its local nature, a local analysis only focuses on one single streamwise location, and cannot account for such a downstream/upstream interaction. Therefore, this may lead to globally unstable flow that would not contain any absolutely unstable region.

Parabolised stability equations

In the past decades, a non-modal approach, introduced by [START_REF] Herbert | Parabolized stability equations[END_REF] and known as a Parabolized Stability Equations (PSE) analysis, has become widely used for the study of weakly non-parallel flows. The detailed PSE procedure can be found in chapter 2. For a given input frequency, this yields a global disturbance field, obtained by marching the parabolized stability equations downstream. This spatial march is classically initialized by the most unstable local mode at a given upstream location, given by a local stability analysis. Chapter 5 shows that this may approximate the optimal response modes given by a resolvent analysis. Note that contrary to a resolvent analysis, which is a fully global approach, the PSE marching procedure consists of solving local equations at successive streamwise locations. Therefore, a PSE analysis may be classified as a local approach, and its computational cost is close to that of a local stability analysis.

However, the PSE approach presents stronger limitations than a resolvent analysis. Indeed, the spatial marching procedure requires equations that avoid information to propagate upstream (as explained, for instance, in [START_REF] Towne | One-way spatial integration of hyperbolic equations[END_REF]). In the case of the PSE, in which some ellipticity remains [START_REF] Haj-Hariri | Characteristics analysis of the parabolized stability equations[END_REF], this is achieved by numerically damping all upstream-propagating waves, either by taking sufficiently large marching steps or by explicitly adding a damping term in the equations [START_REF] Andersson | On a stabilization procedure for the parabolic stability equations[END_REF]. The PSE are therefore designed to track one single downstream-going wave, and may significantly distort and damp all other waves. Therefore, any unsteadiness related to both hydrodynamics and acoustic waves, such as the acoustic feedback-loop behaviors mentioned previously, cannot be treated by PSE, and a global stability or resolvent analysis is then required. The case of screeching jets, studied in chapter 3, is an example of such a situation.

Characterizing a flow with a linear base flow analysis

The previous sections introduced several linear stability techniques, and the following paragraphs present a few results from the literature that highlights some of their limitations. We show that due to their linear nature, they are relevant only in some specific situations. Fig. 1.3 Eigenspectrum of the cylinder flow, computed from a linear base flow stability analysis, for Re = Re c (data extracted from [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF]). σ is the real part of the eigenvalues (the growth rate), and ω the imaginary part (the frequency). The red circled point is the marginally stable mode that becomes unstable for Re > Re c . early as for Re ≈ 1000 (see for instance [START_REF] Patel | Some observations on skin friction and velocity profiles in fully developed pipe and channel flows[END_REF]). This apparent contradiction has been explained by [START_REF] Trefethen | Hydrodynamic stability without eigenvalues[END_REF], and this behavior is related to non-modal linear mechanisms, that have been discussed in section 1.2.1. For such flows, studying the criticality is not physically relevant, and one should rather carry a resolvent or a PSE analysis. [START_REF] Mantič-Lugo | Too big to grow: self-consistent model for nonlinear saturation in open shear flows[END_REF] explained that as we go further from criticality, a linear base flow stability analysis is condemned to fail because of the increasingly strong influence of the nonlinearities. Indeed, at criticality, the linear assumption is justified, but beyond this critical point, its relevancy may quickly be invalidated. A typical example of this may be observed in the case of the cylinder, introduced in the previous section. For instance for Re = 100 > Re c ≈ 46, [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF] showed that starting a nonlinear simulation from the unstable base flow leads to oscillations whose frequency initially matches the stability prediction, but which quickly shift to a higher value. For the whole range of Reynolds numbers that he studied (up to 180), the stability analysis always under-predicted the frequency, and the discrepancy grew as he went further from criticality: for Re = 60, the prediction error (compared to a nonlinear simulation) was approximately 10%, and for Re = 180, this reached almost 50%. This is explained by the increasing importance of the nonlinear mechanisms in the frequency selection process, that is not accounted for by a linear stability about the base flow.

Prediction of the dynamics

However, claiming that a linear base flow stability analysis would systematically fail is erroneous, as several examples in the literature show otherwise. One example is the case of an open cavity flow: a base flow analysis not only predicts the critical Reynolds number for the onset of oscillations, but it also yields a good prediction of the frequency involved. For instance, [START_REF] Yamouni | Interaction between feedback aeroacoustic and acoustic resonance mechanisms in a cavity flow: a global stability analysis[END_REF] studied a cavity flow at a Reynolds number Re = 7500 (therefore far from criticality), for different Mach numbers. They carried out a base flow stability analysis, and found unstable eigenvalues (positive real part) whose imaginary part matches a characteristic frequency of the flow. This property, that we shall call the RPIF (real positive imaginary frequency) property, is therefore fulfilled by a few base flows. This shows that, in some situations, a linear base flow stability analysis yields relevant information about the nonlinear dynamics of a flow, even far from criticality.

This RPIF property may be understood from the work of [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF]. They performed a linear and a third-order weakly nonlinear analysis of the cylinder wake and the cavity flow, based on the small parameter ε = Re -1 c -Re -1 . From that, they showed the existence of a saturated limit cycle. Then, they computed a linear and a weakly nonlinear estimate of the corresponding non-dimensional frequency ω as a function of ε, that we shall Fig. 1.4 Comparison between a weakly nonlinear and a linear analysis (using the results of [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF]), showing the relative discrepancy between the weakly nonlinear estimate of the frequency ω n and the linear prediction ω l , with respect to the Reynolds number, for (a): the cylinder flow and (b): the cavity flow. The dotted line shows the value of the critical Reynolds number for each configuration, and the dashed line represents the asymptotic upper bound of the discrepancy. refer to as ω l and ω n , respectively. For the cavity flow, they found ω l = 7.5 + 3000ε,

(1.7)

ω n = 7.5 + 5000ε, (1.8)

and for the cylinder wake, they obtained ω l = 0.74 + 3.3ε, (1.9)

ω n = 0.74 + 34.3ε.

(1.10)

From these relations, they concluded that in both cases, nonlinearities strongly affect the frequency of the flow in the super-critical regime Re > Re c . But it should be noted that ε is a monotonous bounded function of Re, that verifies for Re > Re c : 0 ≤ ε(Re) < Re -1 c . Similarly, the relative discrepancy e between the weakly non-linear and the linear model, defined as e = |ω n -ω l |/ω n , is also a function of Re, which displays an asymptotic upper bound. For the cylinder, this bound is approximately equal to 83%, but for the cavity, it is less than 6% (see figure 1.4). This reveals that actually, the nonlinearities have only a moderate influence on the dynamics of the cavity flow, which explains that the RPIF property is verified in this case.

Mean flow stability for the study of nonlinear behavior

The previous section has revealed several limits of a linear base flow stability analysis, which stem from its linear nature. In some cases (e.g. the cylinder flow), the frequency selection process is strongly driven by nonlinear effects. Consequently, a linear base flow stability analysis is unable to correctly predict the dynamical behavior of the flow, and this yields situations where the RPIF property is not satisfied. However, there are many papers [START_REF] Mattingly | The stability of an incompressible two-dimensional wake[END_REF][START_REF]Global linear stability analysis of time-averaged flows[END_REF]Pier, 2002, for instance) showing that those nonlinear effects may be factored in by using the time-averaged field (the mean flow) instead of the base flow. The procedures are the exact same as those described in section 1.2, except that the equations are linearized about the mean flow.

A typical example for such an analysis is the cylinder case: [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF] has showed that a linear stability analysis about the mean flow yields a marginally stable eigenvalue (real part equal to zero) with a frequency (the imaginary part) close to that of the von Kármán street (see figure 1.5). A justification of this so-called "real zero imaginary frequency" (RZIF) property has been first studied 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]. They showed that in the vicinity of the bifurcation threshold, the RZIF property is approximately satisfied when the fundamental frequency of the flow dominates the other harmonics. Recently, [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF] more generally demonstrated that this is exactly satisfied for a flow exhibiting monochromatic oscillations. This RZIF criterion has even been considered by [START_REF] Mantič-Lugo | Self-consistent mean flow description of the nonlinear saturation of the vortex shedding in the cylinder wake[END_REF][START_REF] Mantič-Lugo | Too big to grow: self-consistent model for nonlinear saturation in open shear flows[END_REF] to build a self-consistent model of a cylinder flow that predicts the frequency of the vortex shedding for Reynolds numbers up to 110. Note that however, since the growth rate is zero, this cannot be used to study the onset of a given unsteadiness.

The generality of the RZIF property has been studied by [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF], who showed that in the case of thermosolutal convection, the mean field corresponding to standing waves does not satisfy it. They attributed this failure to the nature of the frequency spectrum of the flow, that is rather broad and not strongly peaked. Yet, several papers show that mean flow approaches may yield physically relevant results even for flows exhibiting a broadband spectrum originating from convective instabilities, especially for the prediction of spatially correlated structures linked to a nonlinear unsteadiness. For instance, [START_REF] Gudmundsson | Instability wave models for the near-field fluctuations of turbulent jets[END_REF] focused on turbulent round jets: using an array of microphones, they experimentally measured the pressure fluctuations outside the jet shear layer and showed that the pressure amplitude and phase of the data accurately matched predictions from a PSE analysis about the mean flow. The same conclusions were obtained by [START_REF] Oberleithner | Mean flow stability analysis of oscillating jet experiments[END_REF] in the case of a transitional jet, where the mean flow and unsteady structures were fully characterized by Particle Image Velocimetry measurements: at given frequencies, a local stability analysis around the mean velocity profiles led to the full reconstruction of the 1.4 Scope of the research problem Fig. 1.5 From [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF]: the thick line shows the evolution of the Strouhal number of the cylinder flow (the nondimensional frequency) with respect to the Reynolds number, computed from a direct numerical simulation. The vertical dotted line shows the critical Reynolds number for the onset of the oscillations. The dashed line represent the frequency prediction from a linear base flow stability analysis, and the triangles represent that of a linear mean flow stability analysis. spatial structure (in magnitude and phase) of the perturbation field. These are only a part of the examples showing that a mean flow stability analysis is able to predict some dynamical features of broadband-spectrum flows.

Scope of the research problem

The previous sections have introduced several examples of base flow and mean flow analyses, and while the procedure is very similar in one case or another, there seems to have no clear link between these two approaches. The present chapter has raised several unanswered questions, that need to be addressed in two steps. In one hand, for a few situations, the base flow is found to satisfy the RPIF property, even far from criticality. To our knowledge, the existing literature does not provide any clear explanation to identify a priori the situations where this would be the case.

In the other hand, regarding mean flow approaches, there are a few theoretical works that have addressed the RZIF property. But these studies does not apply to broadband-spectrum flows, and the situation where the frequency spectrum is not peaked is very little understood. As a results, clear conditions for the RZIF property are still wanting, and the general meaning of a mean flow analysis as well as its link with the actual nonlinear dynamics of a flow remain eluding.

Aim and structure of the thesis

This thesis aims at developing robust tools and methods, based on the stability theory, for the study of flows. This requires to improve the state-of-the-art understanding of the existing approaches, and thus, to address the points mentioned above. On account of this, the present work starts with a chapter detailing the main stability tools and numerical strategies that are used in the rest of the study ( §2). Then, the manuscript is divided in two main parts, dedicated to base flow and mean flow stability analyses, respectively. The first part starts with a chapter focusing on the screech phenomenon, an original case that has never been studied through a base flow analysis in the literature ( §3). We evaluate the relevancy of such an analysis in terms of frequency prediction by confronting stability results with the existing literature. These stability results are then analyzed to identify an underlying noise generation mechanism. This part dedicated to base flow analyses is then concluded by a discussion in chapter 4 about the RPIF property, that relies on the findings about the screech phenomenon. This aims at delimiting situations where a base flow analysis should or should not be used to characterize a given dynamics.

The second part is dedicated to the mean flow approach, and contains two chapters ( §5 and §6). Chapter 5 focuses on mean flow analyses for general flow configurations. We give mathematical conditions as well as physical interpretations regarding the mean flow stability approach. This serves as a basis to introduce a new predictive model for the characterization of a flow dynamics. A turbulent backward facing step flow simulation is used to validate the results of this chapter. Then, chapter 6 studies this new predictive model in an experimental context, by focusing on a transitional round jet experiment. In particular, the robustness of the model with respect to physical uncertainties is assessed in this chapter.

Note that three chapters of the present manuscript are based on articles ( §3, §5 and §6), whose content has been modified and adjusted for the present work. These chapters are introduced by some preliminary remarks, that include the complete reference for the corresponding article as well as the original abstract. In the following lines, we quickly detail the division of work between each co-author of the said articles. The experiment has been done by B. Leclaire (BL) and R. Yegavian (RY), with the help of G. Losfeld and C. Illoul (ackowledged in the paper). The PIV processing has been performed by R. Yegavian (RY). The PSE code is a Python porting (done by SB) of a Matlab code provided by LL (the same as that of paper 2). The full reconstruction procedure, including the processing of the velocity snapshots, has been done by SB. The study of the robustness of the method has been done by SB, in collaboration with RY. The paper has been written by SB with feedback from DS, BL, RY.

Remark

The end of chapter 2 contains more details about the contribution of each collaborator with respect to the development of codes.

Chapter 2 -Procedures and numerical strategies

In chapter 1, several stability analysis techniques have been introduced. In the following, we give more details about these techniques, in particular regarding their implementation. The last section presents the personal contribution for the implementation of each technique.

Global stability analysis

Linearization strategy

As explain in chapter 1, a linear global base flow/mean flow analysis requires a linearization of the governing equations of the system. This is also required for resolvent-based approaches. Several options co-exist for this linearization, that can be performed in a continuous or a discretized framework, by an analytic derivation or difference methods. Each method has advantages and drawbacks, that have been discussed by [START_REF] Peter | Numerical sensitivity analysis for aerodynamic optimization: A survey of approaches[END_REF] or [START_REF] Mettot | Stabilité linéaire, sensibilité et contrôle passif d'écoulements turbulents par différences finies[END_REF], for instance. For the present study, we have chosen a fully discrete framework (the equations are first discretized and then linearized), with a linearization performed by a finite differences method. As explained by [START_REF] Mettot | Stabilité linéaire, sensibilité et contrôle passif d'écoulements turbulents par différences finies[END_REF], this is the easiest approach, as all the complexity (boundary conditions, numerical scheme, etc.) is accounted for by the method itself, without any particular treatment. The result is however more easily prone to approximation errors, due to the finite differences method.

We now give details about the linearization procedure, based on the notations introduced in chapter 1: the equations to linearize are taken of the form

∂q ∂t = N (q).
(2.1)

In the following, the tilde-notation refers to discretized operators/vectors. We aim at computing

L(q 0 ) = ∂ Ñ ∂q q0 ,
where q0 is either a base flow or a mean flow. The approach consists of using a computational fluid dynamics (CFD) code (in our case, ONERA's finite-volume compressible code ElsA, see [START_REF] Cambier | The onera elsa cfd software: input from research and feedback from industry[END_REF]), able to compute Ñ (q) for any arbitrary q. We then consider the following first order approximation, valid for any q up to an accuracy of O(ε):

L(q 0 )q = 1 ε Ñ (q 0 + εq) -Ñ (q 0 ) , (2.2)
whose right-hand side may be computed by the CFD code. If q is chosen equal to the i-th unit vector (0, 0 . . . 0, 1, 0, . . . 0, 0) t , equation (2.2) yields the i-th column of L(q 0 ). Therefore, using successive unit perturbation vectors gives the columns of the Jacobian one by one. The number of required evaluations of equation (2.2) is then equal to the number of degrees of freedom of the discrete system.

Note that if q is a unit vector, q0 + εq differs from q0 only at one mesh cell. As a result, this unit-vectors strategy may be improved by taking advantage of the features of the CFD code ElsA, which has a compact stencil represented in figure 2.1. Indeed, Ñ (q 0 + εq) -Ñ (q 0 ) (which contains all elements of one column of the Jacobian) is different from zero only inside the stencil of the cell perturbed by the unit vector. Consequently, if two cells are sufficiently far away, i.e. their stencil do not overlap, then then could not "interfere" with one another. It is therefore possible to simultaneously perturb them, and then separate, in Ñ (q 0 + εq) -Ñ (q 0 ), their respective contribution (each one gives a full column of the Jacobian matrix). This is the cornerstone of the approach introduced by Mettot (2013), who explained that for a one-dimensional case, for a stencil of width n s , the perturbation of all flow variables at all grid points may be optimally achieved by defining perturbation vectors for each variable that act every n s points. The number of evaluations of equation (2.2) needed to get all elements of the Jacobian matrix is then n s n v , with n v the number of flow variables of the system of equations. This is the optimal procedure, in the sense that it is not possible to simultaneously perturb a larger number of cells without having an overlap of stencils. This yields, in one evaluation of Ñ (q 0 + εq) -Ñ (q 0 ), the maximal number of non-zero Jacobian terms.

Fig. 2.1 Stencil dependency of the ElsA CFD code: if a flow variable at the gray cell i, j is perturbed, only the flow variables at the cells within this stencil are affected. Note that several numerical scheme are implemented in ElsA, and this stencil corresponds to that of all the schemes used in this thesis (but this may differ for other schemes).

He then extended this approach to two-dimensional (2-D) problems, by saying that one should perturb every n s points in each direction, with n s the width of the 2-D stencil. In our case, the stencil width is n s = 5 (see figure 2.1), but perturbing the flow every five points in each direction results in the sub-optimal situation represented in 2.2(a), where Ñ (q 0 + εq) -Ñ (q 0 ) would still contain many zero elements. An optimal situation corresponds to that of figure 2.2(b), where the perturbations are arranged in a compact way. This improved strategy has been implemented for the present work. The number of required evaluations of equation (2.2) is then 13n v (13 is the number of cells in the stencil, see figure 2.1). This procedure is easily parallelisable using up to 13n v cores, since each evaluation is independent from the others. Fig. 2.2 Perturbation pattern for the efficient computation of the Jacobian matrix by finite differences, (a): sub-optimal strategy introduced by Mettot ( 2013), (b): optimal strategy used in this thesis. The perturbed cells are in black, and their stencil is represented in gray.

Note that the discretization parameter ε needs to be small enough to avoid approximation errors, but not too small to avoid round-off errors. This has been discussed by [START_REF] Knoll | Jacobian-free newton-krylov methods: a survey of approaches and applications[END_REF], and detailed again by [START_REF] Mettot | Stabilité linéaire, sensibilité et contrôle passif d'écoulements turbulents par différences finies[END_REF] for the present approach. We have found the computations to be robust with respect to this parameter, and for all treated cases, after adimensionalization of the base flow/mean flow, we have set ε = 10 -7 , and have ensured that ε = 10 -8 yields nearly identical results (in particular with respect to the eigenspectrum of the Jacobian). An example of convergence study with respect to ε is presented in section 2.2.

Computation of eigenspectra

The eigenspectrum computations, required for global stability analyses (see chapter 1) have been carried out by Krylov methods with a shift-and-invert strategy, using the open source library ARPACK [START_REF] Lehoucq | ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF]. This strategy requires a linear solver. In our case, we have used the direct LU solver MUMPS (see http://graal.ens-lyon.fr/MUMPS/). The solver is parallel, which makes it potentially very fast. However, the LU factorization step has high memory requirements, that may be non-affordable for large configurations (typically tri-dimensional configurations). Given the size of the problems treated in this thesis (only bi-dimensional configurations), this last aspect was not a limitation. To study larger configurations, one may consider using an iterative solver such as BICGSTAB (see for instance [START_REF] Mack | A preconditioned krylov technique for global hydrodynamic stability analysis of large-scale compressible flows[END_REF]).

Base flow computation

The CFD solver ElsA, used in this thesis, is able to compute a base flow qb (steady solution of (2.1)) using a backward-Euler scheme with a local time-stepping strategy, but this strategy is very time-consuming. An alternative solution is the use of a Newton method, which yields qb by iteratively computing qb,i+1 = qb,i -L(q b,i ) -1 qb,i . When initialized properly, this method is known to converge quadratically.

In our case, given the chosen approaches for the computation of Jacobian operators (see section 2.1.1) as well as for the inversions (LU solver MUMPS, see section 2.1.2), a more efficient strategy consists of skipping the computation of L(q b,i ) at each Newton iteration except the first one, and re-using L(q b,0 ) at each step. This quasi-Newton method requires a unique Jacobian computation and a unique LU factorization (performed during the first iteration), which are the most time-consuming operations in our procedure. This yields a much faster method, despite the larger number of iterations needed to reach convergence.

Note that for the quasi-Newton method to work, the initial guess should be rather close to the solution. Therefore, we used the solver ElsA to get a quick unconverged solution, that was then used to initialize the method.

Example case: a laminar impinging jet

The configuration presented in this section has been studied during a collaboration with O. Browne from the Universidad Politécnica de Madrid. Originally, the aim was a crossvalidation of numerical tools and a comparison of performances, with no physical study of the results. This collaboration has never been achieved, but the work that has been done yields a good illustration case for the global stability procedure, and provides a validation of some of the tools introduced above.

In the following, we present a laminar impinging jet configuration, inspired from that of [START_REF] Chiriac | A numerical study of the unsteady flow and heat transfer in a transitional confined slot jet impinging on an isothermal surface[END_REF]. The global stability procedure is illustrated on this case, and the results are compared with those of [START_REF] Chiriac | A numerical study of the unsteady flow and heat transfer in a transitional confined slot jet impinging on an isothermal surface[END_REF]. This illustrates each step of a global linear stability analysis. The information about the performances and level of convergence that we obtained is representative of the other stability results presented in this thesis.

Physical configuration and numerical parameters

We consider a laminar impinging jet flow. In the following, all the quantities are made non-dimensional by using the density, velocity, temperature and pressure at the nozzle exit and the nozzle height h. The exit Mach number is equal to M = 0.2, and the Reynolds number, based on the hydraulic diameter d = 2, is 660. The jet-to-plate spacing is 5. The computational domain as well as the boundary conditions are represented in figure 2.3. The lower confining wall is adiabatic while the target plate is modeled as an isotherm surface, at T w = 1.1. We consider no-slip conditions for every wall. For the outflow condition, we impose an output pressure p out = 1. To prevent any fluid to enter from these boundaries, two converging sections have been added in both ends of the domain, which ensure that all the fluid goes out (see figure 2.3). The mesh is composed of two blocks, one for the nozzle (denoted as block 1 in the following) and one for the rest of the domain (denoted as block 2 in the following), containing 125 × 120 and 299 × 1063 points, respectively. For the numerical scheme, we used a Jameson flux with scalar dissipation, and a second-order centered scheme for the diffusive fluxes.

Base flow computation

The computation of the base flow has been done by following the procedure described in section 2.1.3. The quasi-Newton method has been initialized with an unconverged base flow qb,0 computed by local time-stepping with the compressible solver ElsA (the level of convergence of this initial guess can be seen in figure 2.4). We have then computed the L2 and infinity norm of the residual Ñ (q b,i ) for each iteration i. The evolution is represented in figure 2.4, and roughly one order of magnitude is gained by iteration. The resulting base flow can be seen in figure 2.5. Note that the good convergence of the quasi-Newton method gives an indirect validation of the Jacobian computation tool (see section 2.1.1), since this method could not work without a rather accurate evaluation of this matrix.

Unstable global modes

A global linear base flow stability analysis has been carried out, following the procedures detailed in the previous sections. The resulting spectrum contains two unstable modes, at St = 0.48 (mode A) and St = 0.55 (mode B) (see figure 2.6). The spatial structure of the modes is shown in figures 2.7 and 2.8, respectively.

As said previously, this impinging jet is close to the configuration of [START_REF] Chiriac | A numerical study of the unsteady flow and heat transfer in a transitional confined slot jet impinging on an isothermal surface[END_REF], who performed incompressible simulations with a similar geometry, for different Reynolds numbers ranging from 250 to 750. The main difference is the compressible nature of our simulation, but the compressible effects should however be rather limited since M = 0.2. For Re = 660, they found that the jet was unsteady, with a characteristic Strouhal number St = 0.55. This compares favorably with the Strouhal numbers of the unstable modes we found (see figure 2.6).

Critical Reynolds number of the unstable modes

Another point, studied by [START_REF] Chiriac | A numerical study of the unsteady flow and heat transfer in a transitional confined slot jet impinging on an isothermal surface[END_REF], may be compared with stability results. From their simulations, they showed that the flow becomes unsteady for a critical Reynolds number between Re = 585 and Re = 610. To compare with this, we have performed a stability analysis for several configurations that cover a range of Re from 560 to 660. As it can be seen in figure 2.9, mode A becomes unstable for Re between 580 and 600, and mode B, between 560 and 580. These values are in good agreement with the results of [START_REF] Chiriac | A numerical study of the unsteady flow and heat transfer in a transitional confined slot jet impinging on an isothermal surface[END_REF], especially for mode A. Note that some differences may appear due to the compressible nature of our computations. Also, the existence of two modes instead of a single one is something that has not been studied here. Further analysis would be required to determine their origin and their role in the dynamics of the flow.

Convergence of the unstable eigenvalues

Mesh convergence To ensure that the two unstable modes have reached convergence with respect to the mesh, the same base flow analysis has been performed with a denser grid, Fig. 2.9 Evolution of the growth rate of mode A (circles) and mode B (triangles) with respect to the Reynolds number. containing 183 × 153 points for block 1 and 439 × 1363 points for block 2. The variation of the modulus of the eigenvalue has been found equal to 0.12% and 0.07% for mode A and mode B, respectively. The variation of the frequency (imaginary part of the eigenvalue) was 0.09% for mode A and 0.02% for mode B, and for the growth rate (real part of the eigenvalue), we found 6.6% for mode A, and 2.7% for mode B. Note that this stronger relative variation of the growth rate represents a small absolute variation, since the growth rate value is very small (σ ∼ 10 -2 ). Based on these results, the modes have been considered to be converged.

Epsilon convergence

The Jacobian matrix, needed to perform the global stability analysis, is computed by finite differences. This implies the use of a small parameter of discretization ε (see section 2.1.1), that was set to ε = 10 -7 for the results presented above. The same analysis has been performed with ε = 10 -8 . The variation of the eigenvalue modulus was then 0.05% and 0.09% for mode A and mode B, respectively. The variation of the frequency (imaginary part of the eigenvalue) was 0.04% for mode A and 0.05% for mode B. For the growth rate (real part of the eigenvalue), we found 3% for mode A and 3.1% for mode B. Based on these results, the modes have been considered to be converged.

Performances

In the following, we give some information about the numerical procedure with respect to memory and computational time. This configuration (initial mesh) contains 1.3 million degrees of freedom (approximately 333000 mesh cells and four variables: density, streamwise and cross-stream momentum and total energy). The computation and storage of the Jacobian required approximately 2 minutes using 52 cores (processors Intel Xeon X5675). The quasi-Newton method (without the initial Jacobian computation) took approximately 30 seconds for the initial Jacobian LU factorization, then 5 seconds per iteration (using 52 cores as well). The matrix contained approximately 45 millions of non-zero elements, and its LU factorization required 28 gigabytes of memory (approximately 540 megabytes per core).

Resolvent analysis

As explained in chapter 1, resolvent analyses require the computation of the Jacobian operator as well as a linear solver and a code able to compute eigenvalues. This can be done using the procedures previously described in sections 2.1.1 and 2.1.2. The following details some other numerical aspects for the implementation of a resolvent analysis tool.

As introduced in chapter 1, a resolvent analysis aims at maximizing the gain G( f ) defined as:

G = ⟨R f , R f ⟩ ⟨ f , f ⟩ (2.
3)

The definition of the inner product is a modeling choice, and while we often consider an inner product corresponding to the energy norm (see chapter 5), any other choice is possible.

In the following, we consider a general situation with two different inner products for the numerator and the denominator, denoted with the subscripts a and b, respectively. This yields

G = ⟨R f , R f ⟩ a /⟨ f , f ⟩ b .
From a discrete point of view, these inner products can be performed by introducing the corresponding symmetric semi-definite positive mass matrices Q a and Q b (for the definition of these matrices for the energy norm, see [START_REF] Mettot | Stabilité linéaire, sensibilité et contrôle passif d'écoulements turbulents par différences finies[END_REF]). The discrete counterpart of equation (2.3) is then:

G = f * R * Q a R f f * Q b f , (2.4)
where the tilde notation refers to discretized quantities, and R * and f * are the transconjugate of R and f , respectively. This corresponds to a generalized Rayleigh quotient, and the optimal gains µ i are given by the largest eigenvalues of the generalized eigenvalue problem R * Q a R f = µ 2 Q b f , which can be solved using the ARPACK/MUMPS strategy described in section 2.1.2. The optimal forcings of unit norm ϕ i are then the corresponding normalized eigenvectors, and the corresponding optimal responses are given by ψ i = Rϕ i . These responses may then be normalized using the inner product defined by Q a .

Remark As explained in [START_REF] Mettot | Stabilité linéaire, sensibilité et contrôle passif d'écoulements turbulents par différences finies[END_REF], this can be further modified to include a restriction matrix, such that the forcing is restrained to a given spatial region for a given set of conservative variables. This has not been used for the present work, and is therefore not detailed here.

Local approaches

In the following, we give details related to the local approaches introduced in chapter 1, namely the PSE and the local stability analysis. We present the main equations and give information about the numerical strategy adopted. Note that the theory for each approach is introduced using a base flow, but the exact same could be done using a mean flow instead.

Parabolized stability equation analysis

A parabolised stability equation (PSE) analysis is a non-modal approach for the study of weakly-nonparallel flows [START_REF] Herbert | Parabolized stability equations[END_REF]. It has been briefly introduced in section 1.2, and is further discussed in chapters 5, where we show that a PSE analysis may be used to approximate the optimal response mode ψ 1 given by a resolvent analysis. In the following, we introduce the PSE equations in an incompressible framework and present the numerical strategy adopted.

PSE in Cartesian coordinates

Let us consider a weakly non-parallel base flow of the form q = (u x (x, y), u y (x, y), p(x, y)), where u x , u y and p are the streamwise velocity, crossstream velocity and the pressure, respectively. A PSE analysis about this base flow consists of considering small harmonic perturbations (u ′ x , u ′ y , p ′ ) at a frequency ω, which are governed by the linearized Navier-Stokes equations, that reads for an incompressible flow:

∂ x u ′ x + ∂ y u ′ y = 0, (2.5) ∂ t u ′ x + u x ∂ x u ′ x + u y ∂ y u ′ x + ∂ x u x u ′ x + ∂ y u x u ′ y = -∂ x p ′ + 1 Re (∂ xx u ′ x + ∂ yy u ′ x ), (2.6) ∂ t u ′ y + u x ∂ x u ′ y + u y ∂ y u ′ y + ∂ x u y u ′ x + ∂ y u y u ′ y = -∂ y p ′ + 1 Re (∂ xx u ′ y + ∂ yy u ′ y ), (2.7)
where Re is the Reynolds number. Note that p is not needed since it does not appear in the equations. The perturbation (u ′ x , u ′ y , p ′ ) is then separated into amplitude functions (ũ x , ũy , p) and an exponential function α. For instance, u ′

x is taken of the form:

u ′ x = ũx (x, y) exp( x x 0 α(ξ)dξ -iωt).
(2.8)

Substituting this form into {2.5;2.6;2.7} yields the PSE equations:

∂ x ũx = -αũ x -∂ y ũy ,
(2.9)

∂ x ũy = 1 u x iω -αu x -u y ∂ y -∂ y u y + 1 Re (∂ y + α 2 ũy - 1 u x ∂ y p,
(2.10)

∂ x p = (iω -u y ∂ y -u x + 1 Re (∂ yy + α 2 )ũ x + (u x ∂ y -∂ y u y )ũ y -αp.
(2.11)

Note that to obtain these equations, one has to consider the slowly varying base flow assumption, which consists of setting ∂ x , u y ∼ O(1/Re), and then neglecting all the terms of order ∼ 1/Re 2 and higher. The previous equations may be recast in the compact matrix form:

∂ x q = Lq,
(2.12)

with q = (ũ, ṽ, p) T , and L the operator corresponding to equations {2.9;2.10;2.11}.

Starting at a certain upstream location x 0 , from an initial value for the shape functions and for α, equation (2.12) can be marched downstream to compute the fluctuation field for x > x 0 : the state vector qj+1 = q(x j+1 ) at a downstream location x j+1 may be computed from the value qj directly upstream from it, by solving the first-order approximation:

(I -∆xL)q j+1 = qj , (2.13) with ∆x = x j+1 -x j the marching step, and I the identity matrix. The initial conditions (q 0 , α 0 ) are usually computed from a local spatial linear stability analysis (detailed in section 2.4.2) using the base flow velocity profile at x 0 . Note that an additional condition must be imposed to remove the ambiguity of the decomposition in (2.8), which allows the streamwise development of the wave to be either absorbed by the shape function or by α.

The auxiliary condition usually considered is: (2.14) where the superscript * denotes the complex conjugate. This relation ensures that most of the streamwise variation is absorbed by α. This provides an update algorithm for α j+1 = α(x j+1 ) (see for instance [START_REF] Gudmundsson | Instability wave models for the near-field fluctuations of turbulent jets[END_REF]):

q * ∂ x qdy = 0,
α n+1 j+1 = α n j+1 - i ∆x (q n j+1 ) * (q n j+1 -qj )dy / (q n j+1 ) * qn j+1 dy (2.15)
By iterating between (2.13) and (2.15) until convergence, the solution can be advanced from x j to x j+1 .

In practice, this procedure may not converge for spatial steps ∆x that are too small. For the present study, we therefore used the stabilizing procedure described in [START_REF] Andersson | On a stabilization procedure for the parabolic stability equations[END_REF], which consists of rewriting equation (2.13) as:

(I -∆xL -sL) qj+1 = (I -sL)q j , (2.16)
where s is a positive scalar which value depends on ∆x and the local value of α (see [START_REF] Andersson | On a stabilization procedure for the parabolic stability equations[END_REF]). With this procedure, arbitrarily small marching step ∆x may be used.

Axisymmetric PSE in cylindrical coordinates

Using similar notations than those of section 2.4.1, the cylindrical PSE for an axisymmetric base flow reads:

∂ x ũx = -αũ x -(∂ r + 1 r ) ũr ,
(2.17)

∂ x ũr = 1 u x iω -αu x -u r ∂ r -∂ r u r + 1 Re (∂ rr + α 2 + ∂ r r - 1 r 2 ) ũr - 1 u x ∂ r p, (2.18) ∂ x p = (iω -u r ∂ r -u x + 1 Re (∂ rr + α 2 + ∂ r r ))ũ x + (u x ∂ r + u x r -∂ r u x )ũ r -αp, (2.19)
with u x and u r the axial and radial base flow velocity, receptively, and ũx , ũr , p the shape function of the axial velocity fluctuation, the radial velocity fluctuation and the pressure fluctuation, respectively. Then, the only difference from the procedure in Cartesian coordinates is the auxiliary condition that becomes rmax 0 q * ∂ x qrdr = 0.

(2.20) Equation (2.15) then becomes

α n+1 j+1 = α n j+1 - i ∆x (q n j+1 ) * (q n j+1 -qj )rdr / (q n j+1 ) * qn j+1 rdr . (2.21)
Numerical strategy For the present study, the discretization of the PSE equations is based on second-order centered finite differences. Equation (2.16), as well as its axisymmetric counterpart, have been solved using the LAPACK library (implemented in the python package Numpy). The integrations involved in equations (2.15) and (2.21) have been carried out numerically with a trapezoidal rule.

Linear local spatial stability analysis

Linear local spatial stability analysis has been introduced in chapter 1, and such an analysis is usually required to initialize a PSE procedure. In the following, we present its main equations and assumptions, and give details about the numerical schemes that have been used for this thesis.

Cartesian coordinates

The 2D incompressible equations of a local stability analysis may be derived by assuming a parallel base flow of the form q = (u x (y), u y = 0, p(y)). Considering small fluctuations (u ′ x , u ′ y , p ′ ) around the base flow gives the following linearized equations:

∂ x u ′ x + ∂ y u ′ y = 0, (2.22) ∂ t u ′ x + u x ∂ x u ′ x + u y ∂ y u ′ x = -∂ x p + 1 Re (∂ xx u ′ x + ∂ y u ′ x ),
(2.23)

∂ t u ′ y + u x ∂ x u ′ y = -∂ y p ′ + 1 Re (∂ xx u ′ y + ∂ yy u ′ y ).
(2.24)

We then assume every fluctuating quantities a ′ of the classical form a ′ (x, y, t) = â(y) exp(αx -iωt), yielding the final set of equations:

αû x + ∂ y ûy = 0, (2.25) - 1 Re ∂ yy + u x α - 1 Re α 2 ûx + ∂ y u x ûy + αp = iωû x , (2.26) - 1 Re ∂ yy + u x α - 1 Re α 2 ûy + ∂ y p = iωû y .
(2.27)

By setting q = (û x , ûy , p) t and A 1 , A 2 , A 3 , A 4 as

A 1 =     0 0 0 -1/Re 0 0 0 -1/Re 0     , A 2 =     1 0 0 u x 0 α 0 u x 0     , A 3 =     0 ∂ y 0 -(1/Re)∂ yy ∂ y u x 0 0 -(1/Re)∂ yy 0     , A 4 =     0 0 0 1 0 0 0 1 0     ,
equations {2.25;2.26;2.27} may be recast as

α 2 A 1 q + αA 2 q + A 3 q = iωA 4 q.
(2.28)

A linear spatial stability analysis consists of computing α and q for a given fixed real value of ω. This can be achieved by considering the augmented eigenvalue problem for α:

0 -A 3 + iωA 4 I -A 2 X Y = α I 0 0 A 1 X Y , (2.29)
with I the identity operator. The pair (α, q = Y ) formed from an eigenvalue and eigenvector of (2.29) is also a solution of equation (2.28).

Axisymmetric flow in cylindrical coordinates

The local stability analysis for an axisymmetric base flow q = (u x (r), u r = 0, p(r)) may be carried by following the exact same procedure than that of section 2.4.2, applied on the following equations

αû x + (∂ r + 1 r )û r = 0, (2.30) ω - 1 Re (∂ rr + ∂ r r ) + u x α - 1 Re α 2 ûx + ∂ r u x ûr + αp = 0, (2.31) ω - 1 Re (∂ rr + ∂ r r - 1 r 2 ) + u x α - 1 Re α 2 ûr + ∂ r p = 0. (2.32)
Numerical strategy For the present study, the equations have been discretised using second-order centered finite differences for the derivatives, then the eigenvalue problem (2.29) and its axisymmetric counterpart have been solved using the ARPACK library [START_REF] Lehoucq | ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF], with a shift-and-invert strategy. For each streamwise location, the computation of the eigenvalue may be initialized by the value of α computed just upstream from the current location. For the initial guess for the most upstream position, the values usually come from physical argument (by knowing the frequency and approximate wavelength of the searched mode).

Personal contributions

In this section, for each of the techniques presented above, we detail the work that has been personally done as a part of this thesis, and we present the collaborators as well as their own contribution.

Linearization procedure (section 2.1.1) The numerical code was based on a previous code created by C. Mettot, a former ONERA PhD student, which has been validated during his PhD [START_REF] Mettot | Stabilité linéaire, sensibilité et contrôle passif d'écoulements turbulents par différences finies[END_REF]. It presented several limitations with respect to the mesh configurations that could be treated (not detailed here), and has been modified to handle any mesh configuration. Moreover, the new optimal strategy for the Jacobian computation, described in section 2.1.1, has been implemented. In addition, as detailed in section 2.1.1, the code requires the use of an external CFD software (ElsA). The integration of ElsA has been improved: the inputs and outputs for this external software are now treated automatically (input parameters, input/output files, etc.). The code has also been parallelized: the matrix computation is based on multiple independent residual evaluations, which are now handled simultaneously on all the available cores. Several other optimizations, not detailed here, have been added to reduce the computational time. A new user interface has been developed with the help of F. Sartor.

Eigenspectrum computation (section 2.1.2) This computation is based on existing libraries (see section 2.1.2), that are called by a simple FORTRAN code, developed by D. Sipp, and slightly modified for a better integration with the Jacobian computation code.

Quasi-Newton code (section 2.1.3) This has been entirely developed and implemented as a part of the present work.

Resolvent analysis code (section 2.3) This tool is based on a previous code developed by C. Mettot. The main contributions correspond to some bug-fixing related to forcing restrictions, and improvements of the user interface.

PSE code (section 2.4.1)

The PSE code developed for this thesis is a porting to Python of a Matlab code provided by L. Lesshafft. The main modification concerns the spatial differentiation schemes, that were initially based on a Chebyshev collocation, and that are now performed by second-order finite differences. This has been done in order to get a more robust code with respect to bias and noise in the input data, motivated by the willing of using it on experimental data (see chapter 6). The counterpart of this robustness is the need of denser spatial grids in order to achieve a similar accuracy. On account of this, the code has been optimized in order to reduce the computational time.

Local stability code (section 2.4.2) The local stability code has been entirely developed and implemented as a part of the present work, with the help of N. Bonne (PhD student at ONERA), who suggested the "augmented eigenproblem" strategy (see section 2.4.2). A compressible version has also been developed (not presented here).

Part I

Base flow stability analysis

Chapter 3 -The screech phenomenon

Preliminary comments

The results of this section have been published in the European Journal of Mechanics B/fluids as "Global stability analysis of underexpanded screeching jets", S. Beneddine, C. Mettot and D. Sipp, 49, pp. 392-399 (2014). The content of the article has been modified to conform with the present manuscript. The original abstract of the article is the following:

This article deals with a global stability analysis of the screech phenomenon. We have shown that a laminar underexpanded supersonic cold jet can exhibit globally unstable modes. A closer look at the structure of these modes shows that they present upstream propagating waves, which is known to be a major component of the screech phenomenon. Furthermore, we find a good agreement between the frequency of the eigenmodes and existing empirical formulas for the prediction of screech frequency. We have then studied the influence of two key parameters on the linear stability of the flow, the jet pressure ratio (JPR) and the nozzle lip thickness, which are known to play an important role in the screech phenomenon. Finally, a careful study of the structure of the unstable modes shows that the upstream propagating acoustic waves of those modes are generated by supersonic phase velocity disturbances, a well-known sound generation mechanism.

Motivation

Chapter 1 introduced the RPIF property, that refers to a base flow which is linearly unstable, with unstable mode(s) displaying a frequency close to that of the actual nonlinear dynamics of the corresponding flow. This property is verified far from criticality for a few configurations, such as the cavity flow. But this is not always the case, and a classical counter-example is the flow past a cylinder.

In the following, we focus on screeching jets and show that this case satisfies the RPIF property. Beside yielding new results about the screech phenomenon, absent from the existing literature, this introduces a new case that is then used for a general discussion in chapter 4 about the RPIF property.

Introduction

The study of imperfectly expanded supersonic jets is an active subject of research since such flows can be found in a broad variety of industrial applications. The most common one is military aircrafts, whose engines often operate at off-design conditions. In this chapter, we focus on underexpanded jets, where the flow pressure at the nozzle exit is higher than the ambient pressure. This mismatch in pressure induces the apparition of a complex quasiperiodic "shock-cell" structure. The jet periodically overexpands and re-converges, attempting to match the ambient pressure, and consequently, forms a standing wave pattern. As a result, shocks and expansion fans appear periodically, creating the so-called shock-cells. Despite the fact that those flows are highly nonlinear, it is possible to predict the gross features of such jets, such as the shock-cell length, with a good agreement with experimental data [START_REF] Pack | A note on prandtl's formula for the wave-length of a supersonic gas jet[END_REF][START_REF] Prandtl | Stationary waves in a gaseous jet[END_REF][START_REF] Tam | A multiple-scales model of the shock-cell structure of imperfectly expanded supersonic jets[END_REF].

One of the important features of supersonic jets is that they can generate strong noise, a point which has been intensively studied in the past decades. A large number of articles on the subject have been written since the first work of [START_REF] Lighthill | On sound generated aerodynamically. i. general theory[END_REF]. One may refer for instance to the review of [START_REF] Tam | Jet noise: since 1952[END_REF] for further details. It is now known that the noise of shock-containing supersonic jets has three components : the broadband shock-associated noise resulting from the interaction of instability waves and the shocks, the turbulent mixing noise generated by the turbulent fluctuations, and the screech tones, which are the subject of this chapter. More information about the two first noise components can be found in [START_REF] Tam | Jet noise: since 1952[END_REF].

The screech phenomenon was first studied by [START_REF] Powell | On the mechanism of choked jet noise[END_REF]. He observed that, under certain conditions, supersonic imperfectly expanded jets can produce very loud discrete frequency tones, the so-called screech tones. This phenomenon can be so intense that in real flight conditions, it can damage the structure of an aircraft. The first observation of such damages was made by the British Aircraft Corporation in the 1960s, where in-flight measurements showed that screech was responsible of minor cracking on VC 10 aircrafts [START_REF] Hay | In-flight shock cell noise[END_REF][START_REF] Raman | Advances in understanding supersonic jet screech: review and perspective[END_REF]. Such concerns do not affect most of modern commercial engines though, and usually, screech tones are observed only with military aircrafts.

In one of his papers, [START_REF] Tam | Jet noise: since 1952[END_REF] referred to screech as "the least understood, least predictable component of supersonic jet noise". Indeed, many questions, such as the prediction of the amplitude of the noise, or its sensitivity to the surrounding environment, have remained unanswered. However, the dominant physical mechanism is known and has been described by [START_REF] Powell | On the mechanism of choked jet noise[END_REF] as a feedback loop between the shocks and the nozzle lip: instability waves developing in the shear layer interact with the shocks, giving birth to acoustic waves propagating upstream. When those waves reach the nozzle lip, they are reflected and excite the shear layer, giving birth to new embryo perturbations that undergo the same process, closing the resonant loop.

There is an abundant literature available on the topic of screeching jets: as mentioned in [START_REF] Raman | Supersonic jet screech: Half-century from powell to the present[END_REF], from Powell's first observation to now, more than 200 papers have been published. An extensive bibliography and a detailed review on screech can be found in the article of [START_REF] Raman | Advances in understanding supersonic jet screech: review and perspective[END_REF]. But despite this large amount of studies, our knowledge of the phenomenon remains mainly qualitative. The only real quantitative prediction available is the frequency of the tones [START_REF] Powell | On the mechanism of choked jet noise[END_REF]. This lack of understanding is the reason why screech is still an active field of research. But, to our knowledge, screech has never been studied in the light of a linear global stability analysis.

Recent works have shown that stability theory appears as a very successful framework for sound prediction in jets. We can cite for example [START_REF] Lesshafft | Aerodynamic sound generation by global modes in hot jets[END_REF] on global modes in adapted subsonic hot jets, the work of [START_REF] Ray | Sound generated by instability wave/shock-cell interaction in supersonic jets[END_REF] on broadband shock-associated noise in supersonic underexpanded jets using parabolized stability equations (PSE, see section 2.4.1 for more details), or [START_REF] Nichols | Global mode decomposition of supersonic jet noise[END_REF], who performed a global mode decomposition on supersonic adapted jets. The case of screech presents one strong particularity: as briefly detailed above, one of the key features of a screeching jet is that it presents upstream propagating acoustic waves that play a major role in the instability process. Consequently, as explained in chapter 1, a stability analysis based on PSE is unable to capture such upstream propagating structures and cannot therefore be used here. On the contrary, a global stability analysis, in which both the cross-stream and stream-wise directions of the perturbation are solved for, is able to capture upstream-propagating waves and may also handle more precisely the non-parallelism induced by the shock-cell structures. In the present chapter, we aim at finding an underexpanded supersonic jet that is globally unstable, and analyze the link between the unstable structures and the screech phenomenon.

The outline of this chapter is as follows. After a presentation of the studied jets and corresponding equations ( §3.3) and the numerical strategy ( §3.4), we focus on an underexpanded jet configuration that is marginally unstable, and relate the features of the unstable global mode to screech ( §3.5). Then, we assess the effects of two key parameters, the jet pressure ratio ( §3.6) and the lip thickness ( §3.7). In a last section ( §3.8), we focus on the noise generation mechanism associated with the unstable global modes.

Physical configuration and linear global stability analysis

We focus on two-dimensional cold jets of air surrounded by a co-flow. In this study, the jet pressure ratio (JPR), defined as the ratio between the static jet pressure and the ambient pressure, and the lip thickness are the two parameters that will be varied. All other parameters of the configuration are fixed: the Mach number of the co-flow is 0.5, while the Mach number of the jet is 1.02. The height of the nozzle exit is equal to 3 mm. The stagnation temperature of both the jet and the co-flow is T 0 = 288 K, the Prandtl number is 0.72 and the viscosity follows a Sutherland law, with standard coefficients for air. The static pressure of the ambient air is set to 3000 Pa. The static pressure of the jet, and thus the Reynolds number (based on the jet velocity, the height of the nozzle, and the static density/temperature of the jet) depend on the JPR. The value of the height of the nozzle and the static pressure of the ambient air have been chosen such that, for all the studied JPR in this paper, the order of magnitude of the Reynolds number is 10 3 , ensuring that the flow is in a laminar transitional situation. To simplify the study, we have imposed adiabatic slip conditions on the walls of the nozzle, so that the boundary layer thicknesses (inside and outside the nozzle) are zero at the nozzle exit. The effect of the boundary layer thickness is not studied here.

The flow dynamics is modeled using the compressible 2D Navier-Stokes equations, that can be recast in the following compact form:

∂q ∂t = N (q), (3.1)
where q = (ρ, ρu, ρv, ρE) T designates the variables describing the flow (density, streamwise momentum, cross-stream momentum, total energy) and N (q) designates the conservation of mass, momentum, and energy equations. From now on, we consider that all quantities are made nondimensional using the jet velocity, the height of the nozzle, the static pressure and density of the ambient air.

We then follow the formalism and notations of chapter 1, such that the linear global stability analysis about the base flow q b reduces to the following eigenproblem:

L(q b )q = λq, (3.2)
where L(q b ) correspond to the linearization of N about the base flow q b .

Numerical strategy

All the computations were performed using the compressible ElsA solver developed at ONERA [START_REF] Cambier | The onera elsa cfd software: input from research and feedback from industry[END_REF]. The spatial discretization of the governing equations is based on a second-order finite-volume conservative formulation. We use a Jameson numerical flux with scalar dissipation with Martinelli correction for the convective fluxes, and a second-order centered scheme for the diffusive fluxes.

The computational domain has 810 nodes in the streamwise direction, and 435 in the cross-stream direction, and covers an area of (165 × 44). Sponge regions, necessary to damp spurious reflections on the boundaries of the domain have been added. We impose respectively supersonic and subsonic injection conditions for the jet and the co-flow. Similarly to the work of [START_REF] Prudhomme | Investigation of supersonic underexpanded jets using adaptive unstructured finite elements[END_REF], if the flow reaches the downstream boundary with a subsonic speed, we impose its pressure to be equal to the static pressure of the co-flow, if it is supersonic, no conditions are imposed. Figure 3.2 shows the sponges, the physical domain of interest and the boundary conditions. Based on these numerical parameters, the base flow, the Jacobian operator L(q b ) and the solutions of equation (3.2) have been computed following the procedures described in chapter 2.

Marginally unstable configuration and associated mode

Physical configuration and base flow computation

We consider a configuration as described in section 3.3, with a JPR equal to 1.12 and a lip thickness h l = 0.063. With those parameters, the Reynolds number is 3020.

Figure 3.3 shows the computed density field. We can see the expected shock-cell structures. We have checked that this solution is spatially converged by considering meshes twice as dense in the streamwise and cross-stream direction, respectively. We have in particular compared the pressure distribution along the central line of the jet for the different computed fields. The maximum relative error, defined as

max x∈[0,xmax] |p(x) -p ref (x)| p ref (x) ,
where p is the pressure distribution along the central line corresponding to the denser meshes, and p ref that of the original mesh, is inferior to 5%.

Finally, to check the validity of our computed base flow, we have compared the length of the computed shock-cells with Tam's formula [START_REF] Tam | A multiple-scales model of the shock-cell structure of imperfectly expanded supersonic jets[END_REF]:

l s = 2 M 2 j -1 M j 1 + γ-1 2 M 2 j 1 + γ-1 2 (γ+1)/2(γ-1)
, (3.3) Fig. 3.3 Left: density field of the computed base flow (JPR=1.12, h l = 0.063). One can notice the expected shock-cell structure. Right: zoom on the exit of the nozzle, on which we have represented the main features of underexpanded jets. The figures display 15 contours equally spaced between 0.9 and 1.4 (nondimensional density unit).

where γ is the heat capacity ratio, l s the nondimensional length of the shock-cells, and M j , the fully expanded Mach number, defined as (see [START_REF] Berland | Numerical study of screech generation in a planar supersonic jet[END_REF]):

M j = 2 γ -1 1 + γ -1 2 M 2 e JPR (γ-1)/γ -1 . (3.4)
For the studied configuration, the formula gives l s = 0.99. Using the computed base flow, we can measure the length of the cells, for instance by evaluating the distance between two local pressure maxima along the center-line of the jet. Table 3.1 shows that the results are in very good agreement with equation (3.3).

Shock-cell number 1 2 3 4 5 Shock-cell length (computed base flow) 0.960 0.978 0.983 0.978 0.995

Relative error (w.r.t. formula (3.3)) 3.0% 1.3% 0.7% 1.3% 0.5% Table 3.1 Comparison between the measured shock-cell lengths from the computed base flow (for the five first shock-cells) and the theoretical values l s = 0.99 (equation (3.3)), for the case JPR=1.12, h l = 0.063. Fig. 3.4 Spectrum of the Jacobian matrix, for the marginally unstable configuration (JPR=1.12, h l = 0.063). Circled red dot: the marginally unstable eigenvalue, vertical red line: frequency prediction (Powell's formula [START_REF] Powell | On the mechanism of choked jet noise[END_REF]).

Global mode decomposition

Figure 3.4 shows the least stable eigenvalues corresponding to the base flow presented in the above section. We can see that there exists one marginally unstable global mode and a series of damped global modes.

Figure 3.5 shows the spatial structure of the marginal global mode. One can clearly notice that the density field presents upstream propagating waves in accordance with the screech phenomenon. We can observe that the radiated sound field is out-of-phase on either side of the jet, in accordance with most of the works carried out on rectangular and planar jets (see the review of [START_REF] Raman | Advances in understanding supersonic jet screech: review and perspective[END_REF]). However, it has also been shown that rectangular jets can sustain a weak symmetric mode [START_REF] Lin | Symmetrical oscillation modes in choked-jet edge tones and screech from rectangular nozzles[END_REF][START_REF] Raman | Advances in understanding supersonic jet screech: review and perspective[END_REF], but several papers suggest that, in the strictly two-dimensional case, this mode should not appear [START_REF] Lin | Symmetrical oscillation modes in choked-jet edge tones and screech from rectangular nozzles[END_REF][START_REF] Umeda | Instability of astrophysical jets. ii numerical simulation of two-dimensional choked underexpanded slab jets[END_REF] To confirm that this mode is indeed related to screech tones, we have compared its frequency with empirical formulas. [START_REF] Powell | On the mechanism of choked jet noise[END_REF] provides a simple formula to predict the frequency of the tones, given by

f = u c l s (1 + u c /c a ) , (3.5)
where f is the frequency, l s the length of the shock-cells, u c the convective speed of the disturbance, and c a the ambient sound velocity. This formula indicates that one period of screech corresponds to the sum of the time needed for a disturbance to reach the first shock (t 1 = l s /u c ) and the time needed for the generated acoustic waves to travel back from this shock to the nozzle lip (t 2 = l s /c a ). This leads to a period for each screech cycle equal to

T = t 1 + t 2 = l s (u c + c a )/(u c c a )
, from which one can easily deduce Powell's formula. The convective velocity u c is generally taken equal to a • u j , with u j the jet velocity, and a a scalar around 0.5 to 0.7 [START_REF] Raman | Advances in understanding supersonic jet screech: review and perspective[END_REF]. Equation (3.5) holds when the ambient air is at rest. In our case, since we have a co-flow, the convection speed of the perturbation needs to be modified according to:

u c = a(u j -u a ) + u a ,
where u a is the ambient velocity. The acoustic feedback is also impacted by the co-flow, and sound waves propagate back to the nozzle with a speed equal to c a -u a instead of c a . This yields the following final relation for the frequency prediction:

f = (a(u j -u a ) + u a ) l s 1 + a(u j -u a ) + u a c a -u a -1 . (3.6)
In the present case, this formula yields ω = 2πf ≈ 1.96. We can see that the frequency of the unstable mode is in reasonable agreement with this value (18 % of error, see figure 3.4).

A part of the discrepancy may stem from the fact that global stability analysis performed around base flows may provide frequencies that are different from those observed on saturated limit-cycles, since nonlinear mechanisms are not accounted for (see section 1.2.3). But we can also point out that predictions based on (3.5) overestimate frequencies for low values of M j (and thus, low values of JPR) with respect to experimental data (see for example Raman (1998)), in agreement with the present results. Therefore, we expect better consistency between the global mode frequency and equation (3.6) for higher values of JPR (see next section). But those first results show that both the structure of the mode and its frequency support the idea that screech is linked to a global instability of the jet.

Influence of the Jet Pressure Ratio

The influence of the JPR on the screech phenomenon is still an open question, in particular there is no simple relation between the JPR and the intensity of the tones [START_REF] Raman | Advances in understanding supersonic jet screech: review and perspective[END_REF]).

Yet, we know that this parameter plays a major role in the phenomenon since it has a strong influence on the length of the shock-cells l s and on the convective velocity u c . An increase of the JPR leads to a larger value of M j , and therefore larger values of u c and l s (see equation (3.3)), the resulting effect being a decrease of the frequency.

Physical configurations

To assess the effect of the JPR parameter on the eigenvalues, we consider seven configurations with a fixed value of the lip thickness h l = 0.063 but with JPR values ranging from 1.12 to 1.215, and compute the associated spectrum. The studied values of JPR and the corresponding Reynolds numbers are reported in table 3.2.

JPR 1.12 1.13 1.14 1.15 1.16 1.18 1.215 Re 3020 3050 3070 3100 3120 3170 3270 Table 3.2 Studied values of JPR and corresponding Reynolds numbers (h l = 0.063).

Results

We have reported in figure 3.6 the full spectrum for JPR=1.12 (circular symbols) and the unstable eigenvalues for higher values of the JPR (other symbols). We can see that there are two branches of unstable modes (red and blue symbols), that will be called modes A and B in the following, respectively. The marginal global mode described in the previous section is an A-mode. The spatial structure of the new unstable mode (mode B) can be seen in figure 3.7. Similarly to mode A, it exhibits upstream propagating waves in accordance with the screech phenomenon. The two modes exhibit a frequency decreasing when the JPR increases, confirming their link with the screech phenomenon. Figure 3.8 displays the two mode branches in terms of growth rate σ and frequency ω, with respect to the JPR. we observe that the two modes display a frequency in much closer agreement with Powell's formula as the JPR is increased. One last point about those results can be commented: the growth rate, for which we have no theoretical results, presents an interesting behavior. Indeed, we can observe that mode A, the first to become unstable, stabilizes again when we increase the JPR, while the other mode has, within the studied range of JPRs, a simpler behavior, with a growth rate increasing when the JPR becomes higher. This shows that mode A is unstable only for a very narrow range of JPRs, meaning (3.3)), and the frequencies of mode A (red) and B (blue). Right: evolution of the growth rate σ of mode A (red) and mode B (blue) with respect to the JPR that mode B is certainly more relevant for the study of screech, which is known to be a phenomenon occurring for a wider range of JPRs.

Influence of the Nozzle lip thickness

The lip thickness is another parameter known to have a major influence on screech tones. As mentioned before, the mechanism of generation of the tones is based on a resonant feedback loop with upstream propagating acoustic waves that are reflected on the lip and excite the mixing layer. The thicker the lip is, the larger the reflecting surface is. The experiment of [START_REF] Ponton | The effects of nozzle exit lip thickness on plume resonance[END_REF] reported that increasing the nozzle lip thickness tends to significantly increase screech tone amplitudes. Also, we know that screech ceases to exist when the fully expanded Mach number M j is large enough [START_REF] Raman | Advances in understanding supersonic jet screech: review and perspective[END_REF], but it appears that even in such non-screeching jets, it can be reactivated by adding a thickener on the lip [START_REF] Raman | Cessation of screech in underexpanded jets[END_REF]. The frequency is also slightly affected by this parameter, and tends to increase when h l is larger, as it has been shown for instance by [START_REF] Kim | Acoustic properties associated with nozzle lip thickness in screeching jets[END_REF]. But, to our knowledge, there is no prediction formula accounting for this parameter. The influence of the lip thickness is still one of the points that is not entirely understood in the screech phenomenon.

To study the impact of the lip thickness on the global stability, we consider five different lip thicknesses (see table 3.3), with a fixed value of the JPR=1.215 (Re=3270).

h l 0.0302 0.0392 0.0630 0.0728 0.0965 Table 3.3 Values of h l that have been considered for the parametric study (JPR=1.215).

Fig. 3.9 Eigenvalues of mode A (red) and B (blue) for different nozzle lip thicknesses h l (JPR=1.215). Darker colors correspond to higher thicknesses. Remark: the eigenvalues corresponding to the two higher lip thicknesses are very close, making hard to distinguish one from the other.

Results

Figure 3.9 shows the evolution of the eigenvalues corresponding to mode A and mode B. Mode A is stable for the three thinner lips, but is unstable with the two thicker lips. A similar behavior is observed in mode B, which has a growth rate increasing when h l becomes larger. As expected, the thicker the lip is, the more unstable the modes are. The frequency is also slightly altered and becomes higher for thicker lips, in agreement with the paper of [START_REF] Kim | Acoustic properties associated with nozzle lip thickness in screeching jets[END_REF] who pointed out that an increase of lip thickness causes a decrease in screech wavelength and therefore an increase of the frequency. One can also notice that the spectrum displays a saturation-like behavior: we can see that both the frequency and the growth rate are practically not altered when we compare the two smaller or the two larger thicknesses. This behavior is reminiscent of the experiments led by [START_REF] Raman | Advances in understanding supersonic jet screech: review and perspective[END_REF], where the strength of the oscillations of the jet increases when the lip becomes thicker, until a certain point beyond which no appreciable change is seen anymore. But this last remark should be considered with caution, as the link between growth rate and noise intensity is not obvious, and the general physical interpretation of the growth rate in the nonlinear dynamics is very unclear (see chapter 4).

Noise generation mechanism

We have identified two unstable modes related to screech, and we have observed that they exhibit both upstream and downstream propagating acoustic waves. Here, we focus on the sound generation mechanism responsible for the birth of those waves.

One of the well-known mechanism of sound generation is the presence of supersonic phase velocity disturbances [START_REF] Tam | Jet noise: since 1952[END_REF]. Global stability analysis provides the temporal and spatial linear behavior of the conservative variables of modes A and B. We have then computed the associated vorticity field for each mode, using a second-order centered scheme for spatial derivations of the velocity (see figure 3.10). Then, the local streamwise phase velocity of the resulting field is given by:

v x ϕ = - ω ∂ϕ/∂x , with ϕ = arg( ω), (3.7)
where ω is the complex vorticity field of the mode. We consider the case JPR=1.18 and h l = 0.063. In the presence of a co-flow, one should consider the relative phase velocity with respect to the ambient velocity u a . Therefore, in our case, we search for disturbances which verify v x ϕ -u a > c a , with c a the ambient sound velocity. Since the ambient Mach number is M a = 0.5, this is equivalent to v x ϕ /c a > 1.5.

As seen in figure 3.11, for the two modes, several streamwise locations verify this relation, confirming that some sound generation is caused by the supersonic phase velocity mechanism. Both modes display evenly distributed supersonic phase velocity locations. Such observation is reminiscent of the work of [START_REF] Panda | An experimental investigation of screech noise generation[END_REF], where the convective velocity of coherent fluctuations inside the shear layer, measured experimentally, presented such a feature. Note that here, the supersonic disturbances that we identified have a positive phase velocity. Therefore, they are likely to be responsible for the downstream-propagating part of the acoustic waves observed in the modes.

As a remark, we want to point out that equation (3.7) is meaningless where the vorticity is too small or equal to zero. Indeed, the computed vorticity values need to be significantly greater than machine precision or numerical errors due to the finite difference scheme used to compute vorticity. The computed phase velocity has no physical meaning in low-vorticity areas. Therefore, for both modes, figure 3.11 shows the phase velocity only where we have ensured that vorticity is high enough.

Concluding remarks

The screech phenomenon has been studied through a linear global stability analysis. We have seen that, under certain circumstances, underexpanded cold jets can be globally unstable. The associated unstable modes have a structure reminiscent of the screech phenomenon. In particular, they all present upstream propagating acoustic waves, whom frequency is in good agreement with the empirical prediction formula of [START_REF] Powell | On the mechanism of choked jet noise[END_REF]. We have performed a parametric study on two parameters, the JPR and the lip thickness, and the behavior of the frequency of the unstable modes is consistent with Powell's formula and existing experimental results. Therefore, the base flows considered in this study satisfy the RPIF property. In such a situation, the results of a linear base flow stability analysis may be used to further study the dynamics of the flow. For instance, we have been able to study a sound generation mechanism: the unstable modes contain supersonic phase velocity disturbances, which is a mechanism known to generate acoustic waves. These successful findings about underexpanded jets are further investigated in the next chapter, dedicated to a more general discussion about base flow stability approaches.

Remark As mentioned in the preliminary remarks, this chapter is based on a published article, which includes supplementary materials available online. Two videos show an animation of each unstable mode, and two other videos show the vorticity of the modes: one can see that the supersonic regions mentioned in 3.8 correspond to locations where the structures accelerate abruptly.

Chapter 4 -Discussion about linear base flow stability

RPIF property and acoustic feedback loop

The previous chapters have introduced two flow configurations satisfying the RPIF property even far from criticality, namely cavity flows ( §1) and screeching jets ( §3). Similarly, the impinging jet case, briefly studied in section 2.2, also verifies this property. These three examples are in stark contrast with the flow past a cylinder, where a similar analysis inaccurately estimates the vortex shedding frequency of the wake, except very close to criticality (see chapter 1). An explanation of those opposite results, already presented in chapter 1, is that for some flows, the frequency selection process is weakly impacted by nonlinearities, which explains why a linear analysis accurately predicts the nonlinear frequency.

In general, determining a priori if an unsteadiness comes from a linear mechanism is a complex question, but the present study reveals one recurrent mechanism: in the three cases where the RPIF property is satisfied, the unsteady behavior stems from an acoustic feedback loop, schematically illustrated in 4.1 (see [START_REF] Rossiter | Wind tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds[END_REF] for the cavity, [START_REF] Powell | On the mechanism of choked jet noise[END_REF] for the screeching jets and [START_REF] Ho | Dynamics of an impinging jet. part 1. the feedback phenomenon[END_REF] for the impinging jet, for instance).

In these flow configurations, the situation is as follows: a shear layer amplifies disturbances through the Kelvin-Helmholtz mechanism, that are convected downstream by the fluid. They eventually generate acoustic waves by interacting with a wall or a shock. Finally, these waves propagate downstream where they excite the shear layer, thus closing the loop. This results Fig. 4.1 Schematic illustration of the acoustic feedback mechanism: perturbations are convected downstream and amplified inside the shear layer. They eventually interact with a wall or a shock, for instance, which creates acoustic waves that eventually got reflected on a surface and excite the shear layer. This kind of self-sustained mechanism, occurring in cavities or screeching jets for instance, is nearly not affected by nonlinearities, yielding situations where the RPIF property is verified.

in an oscillating behavior at discrete frequencies, that only depends on the geometry, the convective and sound velocity, which are, in the present cases, well described by a linear model. While these observations cannot be extended to any flow involving a feedback loop, it still reveals a category of flows where the RPIF property seems to be often observed. Moreover, these examples show that the high amplitude of fluctuations is not a sufficient criterion to invalidates the results from a linear stability analysis.

Base flow analysis using RANS equations

The previous findings also apply to a different base flow stability approach, not used for the present work, but that relies on similar concepts than a classical base flow analysis. A linear So far, we have considered equation (4.1) to be the Navier-Stokes equations, but it could be otherwise. For instance, one may consider the unsteady Reynolds-Averaged Navier-Stokes (RANS) equations. A base flow would then be a solution of the steady RANS equations, and comparing a RANS base flow analysis with respect to an unsteady RANS simulation is analog to confronting a classical base flow analysis to a direct numerical simulation.

The RANS equations describe the large scale dynamics, with the influence of the small scales being accounted for by an additional viscosity term (eddy viscosity). By analogy with the previous section, we expect that nonlinearities in the RANS equations would not affect the frequency selection process when the large scale unsteadiness originates from an acoustic feedback loop. This is indeed the case for the cavity configuration, studied by Mettot et al. (2014a): figure 4.2 shows the good agreement between the eigenspectrum computed from a RANS base flow analysis of a high Reynolds number cavity flow (Re = 860000) and a pressure spectrum computed from an unsteady RANS simulation of the same configuration. Such a RANS approach may be helpful to carry out a base flow stability analysis of high Reynolds number flows, where a steady solution of the Navier-Stokes equations may be very challenging to compute (it may not even exist). Note that the physical relevancy of the unsteady RANS equations is not discussed here, and the results of such an analysis are conditioned by the intrinsic ability of the unsteady RANS equations to describe the dynamics.

Conclusion

While a linear base flow analysis is often relevant to study the onset of some unsteady behavior, the situation is more complex for the prediction of unsteady features, such as the dominant frequency of an oscillating behavior. Chapters 1 and 3 have shown that a few flows do verify the RPIF property because they become unsteady through linear mechanisms.

In particular, flows involving an acoustic feedback loop are likely to satisfy this condition.

Several strategies may then be considered, depending on the Reynolds number Re. If Re is not too high, then a classical base flow analysis may be carried out, similarly to what has been done in the case of the screech. This analysis can then be used, for instance, to study the influence of complex parameters, such as the nozzle lip thickness for the screeching jet. Finally, as shown for the screech phenomenon, the unstable global modes may help to identify, for instance, a noise generation mechanism.

If the Reynolds number is too high, finding a steady solution of the Navier-Stokes equations may be rather challenging, if not impossible. Then, one may proceed by computing a RANS base flow (a steady solution of the RANS equations), and then performing the stability analysis on the RANS equations. This approach could be used to study the onset of large scale unsteadiness. This has been done, for instance, by [START_REF] Crouch | Predicting the onset of flow unsteadiness based on global instability[END_REF]. He predicted the attack angle that yields the onset of the buffeting phenomenon on a NACA0012 airfoil at Re = 10 7 and M = 0.76. This could also be used to study the large scale dynamics, under the condition that the unsteady behavior is nearly not affected by nonlinear terms of the RANS equations, such that the RPIF property is verified. Finally, the cavity example also reveals that comparing the growth rate of co-existing unstable modes gives no reliable indication regarding their relative importance in the nonlinear dynamics. For instance, the dominant eigenvalue of the cavity corresponds to the third harmonics, that is not particularly energetic in the unsteady RANS simulation (see figure 4.2). Therefore, regarding the growth rate, the only relevant information about the nonlinear dynamics relates to its sign and the onset of an unsteady behavior. The value itself gives no quantitative information about the nonlinear regime.

The next part focuses on mean flow stability analyses, that may be used in several situation where a base flow analysis would not be adapted. This includes situations where the RPIF property is not satisfied, or when one cannot compute a steady solution of the Navier-Stokes or RANS equations. It may also be used when the unsteady RANS equations are not able to describe a given dynamics of interest.

Part II

Mean flow stability analysis

Chapter 5 -Conditions for validity of mean flow stability analysis

Preliminary comments

This chapter is based on an article, published in the Journal of Fluid Mechanics as "Conditions for validity of mean flow stability analysis", S. Beneddine, D. Sipp, A. Arnault, J. Dandois and L. Lesshafft, 798, pp. 485-504 (2015). The content has been modified to conform with the present manuscript. The original abstract of the article is the following:

This article provides theoretical conditions for the use and meaning of a stability analysis around a mean flow. As such, it may be considered as an extension of the works by McKeon & Sharma (J. Fluid Mech., vol. 658, 2010, pp. 336-382) to non-parallel flows and by Turton et al. (Phys. Rev.E,vol. 91 (4), 2015, 043009) to broadband flows. Considering a Reynolds decomposition of the flow field, the spectral (or temporal Fourier) mode of the fluctuation field is found to be equal to the action on a turbulent forcing term by the resolvent operator arising from linearization about the mean flow. The main result of the article states that if, at a particular frequency, the dominant singular value of the resolvent is much larger than all others and if the turbulent forcing at this frequency does not display any preferential direction toward one of the sub-optimal forcings then the spectral mode is directly proportional to the dominant optimal response mode of the resolvent at this frequency. Such conditions are generally met in the case of weakly non-parallel open flows exhibiting a convectively unstable mean flow. The spatial structure of the singular mode may in these cases be approximated by a local spatial stability analysis based on parabolised stability equations (PSE). We have also shown that the frequency spectrum of the flow field at any arbitrary location of the domain may be predicted from the frequency-evolution of the dominant optimal response mode and the knowledge of the frequency spectrum at one or more points. Results are illustrated in the case of a high-Reynolds number turbulent backward facing step flow.

Introduction

The previous chapter has focused on base flow stability analyses, which were found particularly relevant in some specific situations where the nonlinearities do not affect significantly the frequency, typically when the instability mechanism is an acoustic feedback loop. But many flow configurations do not fall into this category, and mean flow stability analyses are an alternative solution for the study of such flows, that has turned out to be extremely successful for the prediction of the frequency in many situations (see chapter 1).

The few theoretical results from the literature, that investigate the relevancy of such an approach, mainly focus on self-excited systems which present a strong dominating frequency [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] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF]). These articles address the so-called RZIF (real zero imaginary frequency) property, that has been introduced in section 1.3. The work of [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF] has showed that this RZIF property may not occur for flows with a rather broad spectrum. Nonetheless, many examples from the literature show that a mean flow analysis is still able to characterize spatio-temporal features of flows presenting a broadband spectrum [START_REF] Gudmundsson | Instability wave models for the near-field fluctuations of turbulent jets[END_REF], [START_REF] Oberleithner | Mean flow stability analysis of oscillating jet experiments[END_REF], see section 1.3). On account of this, mathematical foundations for the mean flow approach are still wanting.

In this chapter, we aim at providing, in the general case of a flow field that presents either a broadband or a peaked spectrum, mathematical insights to justify the efficiency of a mean flow stability approach to predict the spatio-temporal features of a flow field. Note that our objective is not to tackle the turbulence closure problem, as in the stochastic structural stability theory (SSST) of [START_REF] Farrell | Dynamics of streamwise rolls and streaks in turbulent wall-bounded shear flow[END_REF]: the mean flow is here considered as known a priori (for example, from experimental measurements) and we aim at establishing the link between the spatial structures of the underlying fluctuation field (the coherent structures) and the linearized equations about this mean flow. The theoretical developments are expressed in a general framework, without consideration of the turbulent/laminar or amplifier/oscillator nature of the flow, and are based on the optimal response modes computed from a singular value decomposition (s.v.d.) of the resolvent operator around the mean flow. The use of the s.v.d. in hydrodynamic stability has first been considered as a way to study energy amplification due to non-modal mechanisms in transitional flows [START_REF] Schmid | Stability and transition in shear flows[END_REF]; [START_REF] Trefethen | Hydrodynamic stability without eigenvalues[END_REF], see section 1.2.1). Non-modal analysis with turbulent mean flows has been introduced by [START_REF] Butler | Optimal perturbations and streak spacing in wall-bounded turbulent shear flow[END_REF], [START_REF] Chernyshenko | The mechanism of streak formation in near-wall turbulence[END_REF], Del Álamo and Jimenez (2006), [START_REF] Cossu | Optimal transient growth and very large-scale structures in turbulent boundary layers[END_REF], [START_REF] Pujals | A note on optimal transient growth in turbulent channel flows[END_REF] and [START_REF] Moarref | Model-based design of transverse wall oscillations for turbulent drag reduction[END_REF]. [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF] showed that the first optimal response computed from an s.v.d. around the mean flow dominates the velocity field of the full flow field when the first singular value is significantly larger than all others. However, they restricted their analysis to a turbulent pipe configuration, which is invariant (homogeneous) in the streamwise direction. Hence, the mean flow was constant in the streamwise direction and all fluctuating quantities were Fourier-transformed in this direction. Therefore, only a local stability analysis was required. The present study is an extension of their work, where we consider more general configurations, in particular open flows which are not invariant in the streamwise direction, such as a backward-facing step or jet configurations. In such cases, a global stability analysis is required, in which all operators and variables depend on the streamwise coordinate in a non-specified (general) manner. As a consequence, the convective non-normality [START_REF] Marquet | Direct and adjoint global modes of a recirculation bubble: lift-up and convective non-normalities[END_REF] becomes a dominant physical mechanism, which is responsible for the downstream location of the optimal response modes and the upstream location of the optimal forcing modes. Finally, we also aim at elucidating the link between the global stability results and those provided by local stability approaches, such as spatial stability or PSE (parabolised stability equations, see [START_REF] Herbert | Parabolized stability equations[END_REF]) analyses.

The present chapter is divided in four sections. First, we highlight the role of the resolvent operator around the mean flow for the determination of the temporal spectral (Fourier) mode of the fluctuation field ( §5.2). Then, we introduce a rank 1 approximation of the resolvent operator to link this spectral mode to the dominant optimal response, and relate this result to local and global stability analyses ( §5.3). This relation yields models for the prediction of the frequency spectrum at arbitrary points of a flow ( §5.4). The last part is dedicated to an application of these results to a turbulent backward facing step ( §5.5).

Resolvent-based equation for the Fourier mode of the fluctuation field

We consider the incompressible, homogeneous Navier Stokes equations:

∂ t u = -u•∇u + ν∇ 2 u -∇p, ∇•u = 0 (5.1)
with u, ν, p respectively standing for the velocity, kinematic viscosity and pressure of the fluid. If we consider a Reynolds decomposition of the flow variables, the mean quantities (denoted by an overline) verify:

-u•∇u -u ′ •∇u ′ + ν∇ 2 u -∇p = 0, ∇•u = 0, (5.2)
and the fluctuating quantities, denoted by a prime:

∂ t u ′ = -u•∇u ′ -u ′ •∇u + ν∇ 2 u ′ -∇p ′ -u ′ •∇u ′ + u ′ •∇u ′ (5.3) ∇•u ′ = 0.
(5.4)

Defining the turbulent forcing term f ′ = u ′ •∇u ′ -u ′ •∇u ′ , and the linear operators L, B and P by:

L = -u•∇() -()•∇u + ν∇ 2 () -∇() ∇•() 0 , B = 1 0 0 0 , P = 1 0 ,
equations {(5.3),(5.4)} can be rewritten in the compact form:

B∂ t u ′ p ′ = L u ′ p ′ + P f ′ .
(5.5)

The temporal Fourier transform of equation (5.5) gives, for each frequency ω:

iωB û p = L û p + P f , (5.6)
with û, p, f as respectively the Fourier transform of u ′ , p ′ and f ′ . Finally, by using the restriction operator R = (1 0) to eliminate the pressure term, equation (5.6) can be rearranged as:

û = R(ω) f , (5.7)
where R(ω) = R(iωB -L) -1 P is the resolvent around the mean flow of the system. Note that (iωB -L) has to be non-singular for the resolvent to be defined. This equation is the cornerstone of the analysis that will follow. It highlights the central role that the linear operator R plays in the natural dynamics of a flow. Note that f is not an external forcing, and cannot be seen as such, since it depends on the perturbation field. In this study, we avoid any closure problem by assuming that the mean flow is known a priori, for example by experimental data (similarly to [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF]). In this situation, equation (5.7) can be seen as a linear system governed by R relating the spectral mode û to a specific forcing f .

Relation between the spectral mode and the dominant optimal response

This section focuses on the study of the linear operator R introduced in the previous section.

For a given harmonic forcing φ and corresponding harmonic response ψ = R φ, we define the gain function µ 2 ( φ) as:

µ 2 ( φ) = ∥ ψ∥ 2 ∥ φ∥ 2 = ⟨R φ, R φ⟩ ⟨ φ, φ⟩ = ⟨R † R φ, φ⟩ ⟨ φ, φ⟩ , (5.8)
with ∥ • ∥ a relevant Hermitian norm, ⟨•, •⟩ the corresponding scalar product, and R † the adjoint of R, defined as the operator R † that verifies ∀a, b; ⟨a, Rb⟩ = ⟨R † a, b⟩. For example, we may consider the energy norm defined as ∥ ψ∥ 2 = ⟨ ψ, ψ⟩ = ( ψ• ψ)dx, where ()•() is the local Hermitian inner product.

As explained in section 1.2.1, the set of optimal forcings of unit norm (ϕ j , j ≥ 1), that maximize µ 2 , can be obtained by solving R † Rϕ j = µ 2 j ϕ j . The solution gain values (µ 2 j , j ≥ 1) and the corresponding forcings (ϕ j , j ≥ 1) are sorted such that (µ 2 1 ≥ µ 2 2 ≥ . . . ), and so the optimal forcing ϕ 1 corresponds to the harmonic forcing which verifies µ 2 1 = µ 2 (ϕ 1 ) = sup φ µ 2 ( φ).

Note that the operator R † R is hermitian, consequently ⟨ϕ i , ϕ j ⟩ = δ ij , meaning that (ϕ j , j ≥ 1) is an orthonormal basis of the forcing space. The projection of f in (5.7) onto this basis then yields: û = R(ω) j≥1 ⟨ϕ j , f ⟩ϕ j .

(5.9)

For each optimal forcing ϕ j , the corresponding optimal response of unit norm ψ j is defined by ψ j = µ -1 j R(ω)ϕ j . Similarly to the set of optimal forcings, the set of optimal responses of unit norm (ψ j , j ≥ 1) defines an orthonormal basis of the response space (⟨ψ i , ψ j ⟩ = δ ij ). Using this basis, equation (5.9) can be rewritten as:

û = j≥1 ψ j µ j ⟨ϕ j , f ⟩.
(5.10)

This final relation is related to the singular value decomposition (s.v.d.) of the operator R, with (µ j , j ≥ 1) the so-called singular values of R, and (ψ j , j ≥ 1) and (ϕ j , j ≥ 1) respectively the left and right singular vectors of R. In the case where, for a given frequency ω, the first singular value µ 1 (ω) is significantly larger than the others, a rank 1 approximation of equation (5.10) yields:

û(x, ω) ≈ ψ 1 (x, ω)µ 1 (ω)⟨ϕ 1 (ω), f (ω)⟩ (5.11)
The approximation error reads:

||û -ψ 1 µ 1 ⟨ϕ 1 , f ⟩|| 2 ||û|| 2 = 1 - µ 2 1 |⟨ϕ 1 , f ⟩| 2 i≥1 µ 2 i |⟨ϕ i , f ⟩| 2 .
(5.12)

In the following, the condition µ 1 ≫ µ j≥2 will be called the dominant singular value (DSV) condition. Note that from equation (5.12), this DSV condition is not strictly sufficient for the approximation to hold. A more rigorous one is:

µ 2 1 |⟨ϕ 1 , f ⟩| 2 ≫ i≥2 µ 2 i |⟨ϕ i , f ⟩| 2 .
(5.13) However, under the reasonable assumption that f does not display any preferential direction toward one of the sub-optimal forcings ϕ j≥1 , this stricter condition reduces to the DSV condition. Note that the accuracy of the approximation is ensured in a global sense (the error is evaluated with a 2-norm), which yields a good approximation in the high-energy regions. In contrast, this might be locally inaccurate, especially in low-energy regions, which are hopefully, in most applications, of lesser interest. Building approximations that hold at any point would require infinity-norm based approaches, which have been considered in hydrodynamic stability only very recently [START_REF] Foures | Localization of flow structures using-norm optimization[END_REF].

We have explicitly indicated the variable-dependency of each term of (5.11), highlighting that the only term on the right-hand side depending on the spatial position x is ψ 1 , since µ 1 is a function of the frequency only, and the result of the scalar product ⟨ϕ 1 (ω), f (ω)⟩ is space independent. This implies that for the frequencies where the DSV hypothesis holds, û is directly proportional to the dominant optimal response ψ 1 . In other words, this result states that the spatial structure of the unsteadiness at these frequencies may be closely related to the dominant singular mode around the mean flow. From the literature (for instance [START_REF] Dergham | Stochastic dynamics and model reduction of amplifier flows: the backward facing step flow[END_REF]), the DSV hypothesis seems fulfilled when a flow presents one dominant convective instability mechanism (for example, the Kelvin-Helmholtz mechanism in a shear layer). This point is discussed in further detail in section 5.5.5, and the rest of the present study is dedicated to such flows. This excludes flows displaying a range of frequencies for which the DSV condition is not verified, due, for example, to the existence of several strong instability mechanisms. These cases could be treated by a higher rank approximation over this particular frequency range. The present results justify why a mean flow stability analysis may accurately predict unsteady spatial features of a flow.

The findings of the present study can also relate to the global stability modes (eigenvectors of the linearized Navier-Stokes operator) around a mean flow. [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF] showed that if the flow field is monochromatic, then the mean flow satisfies the RZIF property, i.e. it exhibits a marginally stable global mode at the frequency of the flow. In the context of resolvent analysis, monochromatic oscillations at a frequency ω 0 yield a monochromatic turbulent forcing f ′ = u ′ •∇u ′ -u ′ •∇u ′ at a frequency 2ω 0 . This is not consistent with the energy content of the response, that is concentrated at ω 0 . This problem can only be solved if the resolvent operator degenerates such that the response is infinite for ω 0 , which happens when the RZIF property is fulfilled. The situation is likely to be similar for periodic behaviors with a very strong fundamental frequency and weak harmonics. The forcing term is then expected to display most of its energy for ω = 2ω 0 , and for the flow response to be consistent with oscillations at ω = ω 0 , the resolvent operator should display a very strong response in the very vicinity of ω 0 . Note that in this situation, the first optimal gain function µ 2 1 will naturally be much higher than all others for this frequency, which would justify the validity of a mean flow stability analysis.

Predictive models for the local frequency spectra

Equation (5.11) is a function of the frequency and of space. Except in the very particular situation where the RZIF property is fulfilled (see the previous section), predicting the frequency spectrum of the flow at different points requires the determination of all the terms of the equation, including f . Two models are proposed in this section, which will provide a prediction of the frequency spectrum. By defining an amplitude function Λ(ω) ≡ µ 1 ⟨ϕ 1 , f ⟩, equation (5.11) may be re-written as û(x, ω) ≈ Λ(ω)ψ 1 (x, ω).

(5.14)

The function Λ is directly related to the energy spectral density e(ω) of the fluctuation field integrated over the whole domain Ω. Indeed, Parseval's theorem yields:

1 2 Ω ∞ 0 u ′ • u ′ dt dx = 1 4π Ω ∞ 0 û • ûdω dx = ∞ 0 ∥û∥ 2 4π dω.
(5.15)

Consequently, e(ω) = ∥û(ω)∥ 2 /4π, and since ψ 1 is of unit norm, equation (5.14) gives:

e(ω) = 1 4π ∥û(ω)∥ 2 ≈ 1 4π |Λ(ω)| 2 .
(5.16)

A first simple assumption consists in considering that e is proportional to the dominant optimal energetic gain function µ 2 1 , which can be achieved by posing Λ = µ 1 (model 1). Note that from the definition of Λ, this implies that ⟨ϕ 1 , f ⟩ = 1 for all ω, which is a very strong assumption since the projection of the turbulent forcing term on ϕ 1 must depend on the frequency and is not equal to 1 in general. Hence, we expect that model 1 may provide (at best) only a qualitative shape of the frequency spectrum. Moreover, this model presents one strong limitation that stems from considering Λ as a positive real-valued function (the actual function Λ is complex-valued). By doing so, we lose all phase-related information, and can only predict the amplitude of the spectra. The main advantage of this model is that it does not require any knowledge of û. Results in section 5.5.6 confirm that this model provides good qualitative insights.

A more precise and rigorous model can be built by assuming a partial local knowledge of û (model 2). For example, from equation (5.14), the function Λ may be defined using one component ûi of û, as Λ(ω) = ûi (x, ω)/ψ 1,i (x, ω). Considering that ûi (ω) is known at a given point x 0 yields Λ(ω) = ûi (x 0 , ω)/ψ 1,i (x 0 , ω). This approach consists in rescaling the response mode for each frequency, such that its amplitude and phase match those of û at a given point. Since the spatial shape of this response mode and û are close, This provides a reconstruction of the spectra over the whole domain. Contrary to model 1 that only predicts a qualitative shape of the spectra, model 2 provides a quantitative prediction of the amplitude, but requires the knowledge of the spectrum at at least one point x 0 . The issue of the choice of x 0 is addressed in section 5.5.6, where the two models are compared in the turbulent backward facing step case.

Application to a turbulent backward facing case

Computation of the mean flow

We consider the case of a two-dimensional backward facing step. From now on, we consider all the quantities to be non-dimensionalized by the step height h and the free-stream velocity u ∞ .

We performed an unsteady 3D simulation by using the compressible CFD solver ElsA [START_REF] Cambier | The onera elsa cfd software: input from research and feedback from industry[END_REF], at low Mach number M = 0.09 (the flow field may therefore be considered as homogeneous and incompressible) and Re = 57500. These parameters correspond to those of an experiment that has been carried in ONERA wind tunnels [START_REF] Gallas | Rapport annuel d'activité du pr hybexcit[END_REF].

The solid walls of the step are modeled by adiabatic no-slip conditions. We impose a fluid injection upstream of the step on the left boundary of the computational domain, such that the boundary layer thickness is δ = 0.35 at the step location (the momentum thickness is θ = 0.041 and the displacement thickness is δ * = 0.056). Periodic boundary conditions are used in the spanwise direction. All other boundary conditions are treated by imposing an output pressure. The computational domain dimensions are 38 in length (23 downstream from the step edge, 15 upstream), 4 in width and 10 in height (downstream from the step). The mesh contains approximately 19 million cells (block-structured mesh, L×H×W: 597 × 185 × 161 for the downstream block, x > 0, and 91 × 79 × 161 for the upstream block, x < 0), and satisfies y + ≤ 1 at the wall boundaries. We then performed a Zonal Detached Eddy Simulation [START_REF] Deck | Zonal-detached-eddy simulation of the flow around a high-lift configuration[END_REF] for which the region upstream of the step edge is described by a steady RANS model, and the region downstream by an unsteady LES model. A 2D mean flow was obtained by averaging the results over 310 convection time units, and over the spanwise direction.

To assess the validity of the mean flow, three methods have been used to compute the streamwise reattachment location x r : (a) the location where the mean streamwise velocity U is equal to 0 at the first cell away from the wall, (b) the location of the dividing streamline of the mean flow, and finally (c) the location where the wall-shear stress ∂U/∂y = 0. These three methods yield x r = 6.45 ± 0.2%, which corresponds to the results from the literature for comparable Reynolds number and boundary layer thickness [START_REF] Adams | Effects of the separating shear layer on the reattachment flow structure part 2: Reattachment length and wall shear stress[END_REF][START_REF] Driver | Unsteady behavior of a reattaching shear layer[END_REF], for instance). The results of (a) can be seen in figure 5.1. A secondary reattachment point is located at x = 0.88, which is consistent with the work of [START_REF] Spazzini | Unsteady behavior of back-facing step flow[END_REF], who found a secondary separation at x = 1.1 for a Reynolds number of 16000. They also showed that this length tends to decrease as the Reynolds number increases. Figure 5.1 displays a third corner recirculation region near (x, y) = (0, 0), with a size of 0.053. As a comparison, this value is close to that found by [START_REF] Le | Direct numerical simulation of turbulent flow over a backwardfacing step[END_REF] at a lower Reynolds number (0.042 for Re = 5100).

The unsteady features computed from the simulation also present a good agreement with the literature. The spectrum at (x = 2, y = 1.5) is shown in figure 5. free-stream velocity and the momentum thickness at the step location. This is the frequency expected for the Kelvin Helmholtz instability (see for example [START_REF] Chun | Control of turbulent separated flow over a backward-facing step by local forcing[END_REF]). Further downstream, near to the reattachment, the flow presents another energetic frequency.

The associated Strouhal number F + based on the free-stream velocity and the reattachment length is around 0.6 (ω = 0.7) (see figure 5.2(b)), which is consistent with the values that can be found in other studies (values between 0.5 and 1, see for example [START_REF] Dandois | Numerical simulation of active separation control by a synthetic jet[END_REF]).

Singular values and DSV hypothesis

Since the configuration is invariant in the spanwise direction, the fluctuation field can be expressed via its spatio-temporal Fourier coefficients, depending on the frequency ω and on the spanwise number β. The resolvent operator therefore also depends on β and we expect that the Fourier mode û(ω, β) is proportional to

ψ 1 (ω, β) if µ 1 (ω, β) ≫ µ 2 (ω, β).
In the following, we will focus on two-dimensional perturbations (β = 0), for which the spanwise component u ′ 3 is null. An s.v.d. of the resolvent matrix about the 2-D mean flow has been performed following the procedures described in chapter 2. Figure 5.3 presents the three highest singular values over a wide range of frequencies. It shows that the first singular value µ 1 is several orders of magnitude larger than all others, as expected from a flow presenting a strong convective instability mechanism (see section 5.5.5). Consequently, following the results of section 5.3, at all these frequencies, the spectral mode û is expected to be proportional to the dominant optimal response ψ 1 . Fig. 5.3 Comparison of the first three energetic gain functions in the case of two-dimensional perturbations (β = 0): one can see that the first singular value µ 1 is significantly higher than the others, which ensures the validity of the leading singular value assumption (section 5.3).

Computation of the Fourier modes

Checking the proportionality between Fourier and dominant resolvent modes requires the computation of û, from the finite-length discrete temporal data given by the simulation. Spectral properties of such a signal, that is not square integrable, is usually computed from a statistical point of view, by dividing the time series of spanwise-averaged velocity snapshots u ′ (x, y, t, β = 0), computed from the simulation, into N b bins. These bins are then processed to extract a statistically converged spectral information out of the flow. The final spectrum may be obtained by two approaches: (a) the amplitude of the spectrum is computed from the r.m.s. of all the bins, which corresponds to the classical Welch algorithm:

|û welch | =   1 N b N b k=1 |û k | 2   1/2
, (5.17) (b) the amplitude and phase of the spectrum are computed by performing a POD-filtering of these bins. Method (b) aims at isolating the most energetic, spatially correlated spectral field from each bin. This processing presents one strong advantage over method (a): it yields not only the modulus of the modes, but also the phase. It has been used for instance by [START_REF] Gudmundsson | Instability wave models for the near-field fluctuations of turbulent jets[END_REF], and consists of computing, for each bin k and for each frequency ω, the cross-spectral tensor T (k,ω) defined as: (5.18) where x i,j are the discretisation points. Then, one may compute the ensemble-averaged tensor T ω 0 :

T (k,ω) ij = ûk (ω, x i ) • ûk (ω, x j ),
T ω 0 = 1 N b N b k=1 T (k,ω) , (5.19)
and finally solve the eigenvalue problem:

T ω 0 y ω = λy ω .

(5.20)

By considering the unit-norm eigenvector y ω m associated to the highest eigenvalue, one may obtain the POD-filtered spectral mode ûpf (ω) = A ω y ω m , where the amplitude A ω is defined as:

A ω = 1 N b N b k=1 |⟨û k (ω), y ω m ⟩|.
(5.21)

First optimal response vs. Fourier mode

For the present study, we have N b = 30, with a 50% overlap between the bins. Each bin has a frequency resolution of ∆ω = π/10. To ease the comparison, each field has been normalized such that its maximal amplitude is 1. Profiles extracted from these fields at three streamwise locations (x = 2, x = 3.5, x = 5) are respectively compared in figures 5.4(d), 5.4(e) and 5.4(f). These figures show that for ω = 2.1, the agreement between the optimal response and the spectral mode of the flow is significantly improved when POD-filtering is used. The mode obtained from method (a) displays moderately energetic fluctuations near the wall, that are not found in the resolvent mode. Since this discrepancy is absent from the POD-filtered mode, it means that this is associated to non-spatially correlated motion near the wall. Similar results were obtained for all the tested frequencies of the present study. This is reminiscent of the work of [START_REF] Gudmundsson | Instability wave models for the near-field fluctuations of turbulent jets[END_REF], who made the same observation in the downstream region of a turbulent jet, and attributed this to some uncorrelated turbulent motion. This shows that the optimal response mode only describes the coherent motion of the flow, which fortunately accounts for most of the energy of the fluctuations, despite their turbulent nature.

Method (b) also gives the phase structure of the spectral modes, and for all the frequencies investigated, it displays a strong similarity with the phase of the optimal response. As an example, figure 5.5 illustrates this statement again for ω = 2.1. Note that from equation (5.14), û and ψ 1 are approximately equal up to a complex multiplicative constant, therefore their phase fields are equal up to an additive constant. To ease the comparison, the phase of the optimal response has been shifted such that it matches the one of the spectral mode at an arbitrary position, here x = (3, 1.5). Note that despite a very good overall agreement, a careful comparison of figures 5.4(a)/(b) and figures 5.5(a)/(b) shows that the strongest discrepancy between the first optimal response and the spectral mode appears in low-energy regions, where the spectral mode is not smooth (for instance the left bottom corner in figure 5.5(a)). Increasing the simulation time should yield better agreement, but there is no guarantee of the capability of the rank 1 approximation to finely describe these low-energy regions. This issue is addressed in section 5.5.6.

This section introduced two practical definitions of the spectral mode, and method (b) has yielded a better agreement with the optimal response. Therefore, in the following, the remaining discussions of the chapter concerning the spectral mode all refer to the mode computed with method (b), which corresponds to the spectral description of the spatially correlated motion of the flow.

Link with local stability analysis

The question of the link between a resolvent analysis and a local stability analysis has been addressed by [START_REF] Sipp | Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer[END_REF]. They showed that a flow presenting a strong convective instability, such as a backward facing step flow, displays respectively an upstream/downstream location of the optimal forcings/responses. In such a situation, the spatial growth of the first optimal response compares favorably with the growth rate computed from a local spatial stability analysis: by defining the streamwise energy density d ψ1 (x, ω) of ψ 1 (ω) as

d ψ1 (x, ω) = ymax y min ψ 1 (x, y, ω) • ψ 1 (x, y, ω) dy, (5.22)
then the energy growth σ ψ1 = (∂d ψ1 /∂x)/2d ψ1 at a given frequency is close to the spatial amplification rate obtained from a local spatial stability analysis. Figure 5.6(a) compares these two quantities in the backward facing step case for ω = 2.1, and shows that they display a good agreement in the region where the optimal forcing is weak (the local spatial stability results have been obtained by following the procedure detailed in section 2.4.2). This link between a local stability analysis and the dominant optimal response also justifies the argument that when a convective instability mechanism of a flow is strong, then the dominant optimal forcing ϕ 1 will produce a particularly energetic response (see [START_REF] Garnaud | The preferred mode of incompressible jets: linear frequency response analysis[END_REF]; [START_REF] Sipp | Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer[END_REF]). Therefore, the associated gain µ 2 1 will be high. Thus, the configurations presenting a strong convective instability mechanism are good candidates to fulfill the DSV assumption required for the results of section 5.3. This conjecture is supported by the literature, which shows that such flows indeed present a clear separation of singular values [START_REF] Boujo | Sensitivity and open-loop control of stochastic response in a noise amplifier flow: the backward-facing step[END_REF][START_REF] Dergham | Stochastic dynamics and model reduction of amplifier flows: the backward facing step flow[END_REF][START_REF] Sartor | Unsteadiness in transonic shock-wave/boundarylayer interactions: experimental investigation and global stability analysis[END_REF].

We will now show that not only the energy growth rates of the spatial most amplified k + mode and of the dominant optimal response mode compare well, but also the shape of the modes, at least in the locally unstable region (where the local mode is amplified). In this region (x < 3 for ω = 2.1, see figure 5.6(a)), the shape of the most amplified local k + mode at a given streamwise location is close to the leading optimal response profile at the same location (figure 5.6(b)). However, in the region where the mode is spatially damped, this agreement may strongly deteriorate due to the fact that the spatial structure of the leading optimal response mode may no longer be dominated by one single spatial mode. In the backward facing step flow, this leads to a strong discrepancy between the first optimal response and the local eigenmode in the spatially damped regions of the mode, as illustrated in figure 5.6c for x = 3.1 and ω = 2.1. Note that this failure of the local spatial stability analysis was observed for all the tested frequencies of the present study. Fig. 5.6 Comparison between dominant optimal response mode (black continuous line) and most amplified k + mode obtained with local spatial stability analysis (red dashed line) for ω = 2.1: (a) streamwise energy growth of the optimal response vs. spatial growth rate computed from a local spatial stability analysis, (b) amplitude profile of the optimal response mode at x = 0.75 vs. local spatial mode at the same location (streamwise velocity component), (c) amplitude profile of the optimal response mode at x = 3.1 vs. local spatial mode at the same location (streamwise velocity component). The local stability profiles are normalized such that their maxima match those of the optimal response profile. This problem can be overcome by solving the parabolised stability equations (see section 2.4.1), initialized just downstream of the step (amplified region) by the most amplified local k + mode. This procedure yields an accurate reconstruction of the perturbation field in both the amplified and damped region of the mode. Figures 5.7(a) and 5.7(b) illustrate the similarity between ψ 1 and the PSE field, for ω = 2.1 (to ease the comparison, both fields have been normalised such that their maximum amplitude is 1, and such that their phase match at an arbitrary point, here (x = 3, y = 1.5)). This agreement was observed for all the investigated frequencies. This shows that under the condition that the resolvent fulfill the DSV property (or that the flow displays a strong local spatial growth rate), solving the PSE equations around the mean flow may provide a good prediction of the spatial structure of the unsteadiness.

Note that, even if here it gives excellent results, PSE is not guaranteed to work for any configuration since it may damp any behavior not consistent with a single downstream-going spatial wavelength (see section 1.2.1).

As a conclusion, the link between the resolvent analysis and the local spatial stability analysis shows that the findings of the present study also hold when based on quantities obtained from a local stability analysis: if a mean-flow exhibits a strong convective instability at some frequency (the spatial growth rate is strictly positive in some region of the flow), then the structure of the spectral mode at that frequency should correspond to the structure obtained with a PSE analysis based on the mean flow.

Comparison of the predictive models for the local frequency spectra

The ability of model 1 and model 2 (defined in section 5.4) to predict the frequency spectrum of the flow field at given points is assessed by comparing the resulting estimated energy spectrum with results obtained from the unsteady simulation data.

The question of the choice of the matching points for model 2 needs to be addressed beforehand (see section 5.4). Since model 2 consists of predicting the frequency spectra at every point of the domain by matching, at all frequencies, the optimal response mode and the spectral mode at given spatial locations, it may lead to poor predictions if these locations are ill-chosen. In section 5.5.2, we showed that, for a given frequency, the energetic regions were well predicted by the dominant optimal response. But the prediction was not so accurate in low-energy regions. This probably stems from the inability of the rank 1 approximation to capture the low-energy regions of the spectral mode. Hence, at a given frequency, low-energy regions should be avoided for the choice of the matching point. This issue is illustrated in figure 5.8(a): for an arbitrary frequency ω = 5, the amplitude profiles at x = 1 of û and ψ 1 are rescaled such that their amplitudes match at x 0 = (1, 1.4), which is a point displaying a low level of energy at this particular frequency. This results in two quite different profiles. However, when the rescaling is based on a point whose energy is high, for instance x 1 = (1, 0.9), the overall agreement is strongly improved (figure 5.8(b)). A robust way to compute the function Λ of model 2 is to consider several points, such that at least one point is located in an energetic area of the flow at all the frequencies of interest, and then solving the over-determined set of n equations {Λ(ω) = ûi (x k , ω)/ũ 1,i (x k , ω), 1 ≤ k ≤ n} by Fig. 5.8 Illustration of the issue of the choice of the matching point for the determination of Λ. The three figures display profiles at x = 1 extracted from û (continuous line) and the optimal response (red dashed line), for ω = 5. Figures (a) and (b) respectively compare the resulting profiles when the optimal response is scaled to match the amplitude of the low-energy point (x = 1, y = 1.4) and the higher-energy point (x = 1, y = 0.9). Figure (c) shows the result after a scaling that minimizes the square of the error at those two points. least squares, where n is the number of points considered. Figure 5.8(c) shows the result of least squares fitting based on the two points x 0 = (1, 1.4) and x 1 = (1, 0.9). The resulting profile is close to what is obtained when x 0 is not considered. In brief, the least-squares technique chooses automatically the best matching point at each frequency (several matching points are generally mandatory due to the fact that the energetic regions of the spectral modes are not the same at all frequencies).

For the present study, we used two points for the determination of Λ, located at x 0 = (4, 1.5) and x 1 = (7, 0.1). These points have been chosen in order to correctly capture both the high frequency behavior near the edge and the low frequency unsteadiness in the reattachment region. Figure 5.9 presents the results for the prediction of the spectra of the streamwise component û1 of the velocity, at five points x a = (3, 1.5), x b = (5, 3), x c = (6, 1.5), x d = (8, 0.1) and x e = (9, 0.1). Note that since model 1 is unable to predict the amplitude of the spectra, it has been rescaled such that its maximum matches the one predicted by model 2. Model 1 gives a reasonable order of magnitude for the energetic frequencies of the flow, and in particular it is able to qualitatively discriminate if a point displays either a high or a low frequency content. On the other hand, model 2 led to a much more accurate prediction of those frequencies. The overall shape of the spectra is well-predicted, as well as the amplitude that remarkably agrees with the simulation data for all five tested points.

Concluding remarks

The resolvent operator around the mean flow naturally appears in the full nonlinear equations that govern the perturbation field around the mean flow. If the mean flow exhibits a strong convective instability, the resolvent is expected to fulfill the DSV condition for the frequencies dominated by the instability. Consequently, by assuming that the turbulent forcing does not display any preferential direction toward one of the sub-optimal forcings, the resolvent can be approximated by using the first singular value and the associated optimal response and forcing modes. This approximation yields a mathematical link between the actual dynamics of the flow (the spectral mode) and the dominant optimal response mode. In the case of weakly non-parallel flows, we have shown that the DSV assumption generally holds at some frequency if a strong convective instability mechanism affects the mean flow at that frequency. In such a case, the spatial structure of the spectral mode should be very close to the spatial structure obtained by a PSE analysis around the mean flow. This article also focuses on the development of methods to predict the frequency spectrum at every point of a flow. It presents a model (model 2) based on partial information of the flow, in order to reconstruct the complete fluctuation field at any point. In experimental applications, it is often easy to accurately obtain the mean flow as well as an accurate spectrum at a few points in the flow, for example wall-measurements such as pressure or skin-friction. In the context of compressible flows, the experimental jet setup from [START_REF] Gudmundsson | Instability wave models for the near-field fluctuations of turbulent jets[END_REF] uses an array of microphones measuring the pressure fluctuations just outside the jet shear layer. The approach that is developed in the present study could give a reconstruction of the local spectra at any point from this data and from a measurement of the mean flow. The extent of the region in which mean flow measurements are required also depends on the method chosen to compute the dominant optimal response modes. This region is smaller when using the PSE, since the computation is restricted to the spatial support of those modes. In contrast, a much larger domain is required in the streamwise direction for a global stability computation in order to adapt to the upstream and downstream boundary conditions. Finally, time-resolved data for frequency spectra are only needed at a few points; the mean flow may be obtained by less sophisticated means, such as 5 hole probes or traditional (non time-resolved) Particle Image Velocimetry measurements. The methodology presented here is well-suited for the exploitation of experimental results, and the next chapter is dedicated to the study of this approach in an experimental context.

Chapter 6 -Mean flow stability for flow reconstruction: an experimental study

Preliminary comments

The present chapter focuses on the previous findings about a mean flow stability analysis. These are studied in an experimental context in order to yield an alternative approach to existing measurement methods.

This chapter is based on an article submitted to the Journal of Fluid Mechanics as "Timeresolved reconstruction of a round jet from point-wise measurements and mean flow stability analysis", S. Beneddine, R. Yegavian, D. Sipp and B. Leclaire. The content has been modified to conform with the present thesis. The original abstract of the article is the following:

This article presents a reconstruction of the unsteady behavior of a round jet at a Reynolds number equal to 3300, from the sole knowledge of the time-averaged flow field and one pointwise unsteady measurement. The reconstruction approach is an application of the work of [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF] and relies on the computation of the dominant resolvent modes of the flow, using a PSE analysis. To validate the procedure, the unsteady velocity field of the jet has been characterized by time-resolved particle image velocimetry, yielding an experimental reference. We first show that the dominant resolvent modes are proportional to the experimental Fourier modes, as predicted by [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]. From these results, it is then possible to fully reconstruct the unsteady velocity and pressure fluctuation fields, yielding a flow field that displays a good agreement with the experimental reference. Finally, it is found that the robustness of the reconstruction mainly depends on the location of the point-wise unsteady measurement, that should be within energetic regions of the flow, and this robustness as well as the quality of the reconstruction can be greatly improved by considering a few point-wise measurements instead of just one.

Introduction

A large quantity of work has been devoted to the improvement of measurement methods, which are essential tools for the study of physical mechanisms. Nowadays, most of the conventional methods in fluid mechanics, such as hot wire probes or microphones, are able to give an accurate point-wise time-resolved characterization of a given physical quantity of a flow. When the need for flow characterization goes beyond a single point, it is possible to use arrays of sensors, but this has material limitations toward the spatial extent and resolution of the characterization, and may be too intrusive to get a global flow field measurement.

Alternatively, a single sensor may be displaced to a large number of locations, yielding an arbitrarily dense set of measurements, which are however uncorrelated due to their nonsimultaneity. Another quite widely used option consists of using Particle Image Velocimetry (PIV), that classically yields the two or three instantaneous displacement components, by acquiring and processing two images of the seeded flow separated by a very short time interval. Due to technical constraints related to illumination (usually with pulsed lasers) and camera imaging, the technique can still be considered as characterized by a trade-off between accuracy and temporal resolution. Indeed, low frame rate lasers have a high energy per pulse that guarantees a high signal-to-noise ratio in the images and a good measurement accuracy, but such traditional PIV systems cannot resolve the unsteady flow behavior. On the other hand, high frame rate systems (Time-Resolved PIV, TR-PIV) can characterize frequencies up to 10 kHz, however at the cost of a much lower signal-to-noise ratio, possibly hindering the accuracy unless some specific advanced processing are used (see for instance [START_REF] Jeon | Fluid trajectory evaluation based on an ensembleaveraged cross-correlation in time-resolved piv[END_REF]; [START_REF] Lynch | A high-order time-accurate interrogation method for time-resolved piv[END_REF]; [START_REF] Yegavian | Lucas-kanade fluid trajectories for time-resolved piv[END_REF]). Besides, it is worthwhile noticing that this highest measurable frequency remains one or two orders of magnitude lower than that of a hot-wire probe for instance, which can make a difference in the context of high-speed flows.

In view of these limitations, reconstructing the time-resolved flow field based on quantities measurable by point-wise sensors and/or traditional, low frame rate PIV can be of great interest. Several reconstruction techniques exist to rebuild global information from point-wise measurements, and among them, the stochastic estimation (SE) is one of the most widely used in fluid mechanics. Initially introduced by [START_REF] Adrian | Conditional eddies in isotropic turbulence[END_REF] as a way to extricate the coherent structures in a turbulent flow, this technique has been extensively used to get the instantaneous least-mean-square error estimate of the velocity at various locations, with the sole information of the velocity at a few other points (see for instance [START_REF] Adrian | Conditional eddies in isotropic turbulence[END_REF]; [START_REF] Cole | Applications of stochastic estimation in the axisymmetric sudden expansion[END_REF]; [START_REF] Guezennec | Stochastic estimation of coherent structures in turbulent boundary layers[END_REF]; [START_REF] Stokes | Multi-point measurement techniques used in the study of separated flows[END_REF]; [START_REF] Tung | Higher-order estimates of conditional eddies in isotropic turbulence[END_REF]). This requires to have access to simultaneous unsteady measurements at points of interest. The method may also yield the estimation of the pressure, but requires the use of a higher-order SE model than for the velocity estimation (see [START_REF] Hudy | Stochastic estimation of a separated-flow field using wall-pressure-array measurements[END_REF]; [START_REF] Murray | Estimation of the flowfield from surface pressure measurements in an open cavity[END_REF]; [START_REF] Naguib | Stochastic estimation and flow sources associated with surface pressure events in a turbulent boundary layer[END_REF]), for which the linear SE (LSE) gives satisfactory results. Other more advanced methods, relying on similar techniques, have been introduced for instance by [START_REF] Tu | Integration of non-time-resolved piv and time-resolved velocity point sensors for dynamic estimation of velocity fields[END_REF]. They elaborated a three-steps estimation of the unsteady field, that uses down-sampled TR-PIV snapshots and point sensors. Their method is based on a variant of LSE, coupled with proper orthogonal decomposition, Kalman smoothing and Kalman filtering. They obtained satisfactory results for the estimation of the wake behind a flat plate at a Reynolds number of 3600, but to the cost of a rather heavy processing of the data. Moreover, similarly to classical SE, this requires a simultaneity of the PIV acquisition and the sensor data.

The SE as well as any global estimation technique based on local measurements naturally relies on a certain degree of spatial correlation in the flow field and the existence of coherent structures. It is now generally admitted that even fully-turbulent flows present such structures, which has been addressed in chapter 5 from a stability point of view. We have showed that these structures relate to resolvent modes about the time-averaged flow field (mean flow).

Based on that, chapter 5 has demonstrated that it is possible to estimate the frequency spectrum at any point of a flow from the knowledge of the mean flow and a few point-wise measurements. The analysis was based on a rank-one approximation of the resolvent. We have also showed that cheaper techniques such as a PSE analysis may be used to approximate the dominant resolvent mode. A similar work has been conducted by Gòmez et al. (2016a), who built a reduced-order model of a three-dimensional lid-driven cavity at Re = 1200, based on resolvent modes. Their model yielded a flow reconstruction that accurately compared with direct numerical simulation (DNS) results. Alike the work presented in chapter 5, their input data were the mean flow and a few local unsteady measurements. More recently, Gòmez et al. (2016b) used the same model to estimate aerodynamic forces from point-wise data, and they once again successfully compared their results to DNS.

These theoretical articles claim that the use of such a model could be, in some situations, an interesting alternative to TR-PIV. Yet, to our knowledge, the only experimental demonstration that stability techniques could be used to rebuild a flow field from point-wise measurements is a proceeding paper by [START_REF] Sasaki | Real-time modelling of wavepackets in turbulent jets[END_REF]. They considered an experimental turbulent round jet, and measured pressure fluctuations near the shear layer with microphone rings, located at several streamwise locations. They assumed the existence of a linear transfer function between upstream and downstream measurements, and qualitatively showed that a PSE analysis of the mean field may well-predict this function. However, they do not provide any theoretical background explaining this, and their study relied on a very different formalism than that of chapter 5, which introduced a more appropriate framework to understand PSE/resolvent-based models. Moreover, they presented the accuracy of the reconstruction in a qualitative way, and did not address the robustness of the reconstruction. Therefore, if used in a new experimental context, its reliability would be unknown, and any new practical application would lead to highly uncertain results.

Scope of the chapter

The following sections present a temporal reconstruction method, inspired from the model introduced in chapter 5, and investigate its accuracy and robustness in an experimental context. The case that is considered for this study is a transitional round jet at a Reynolds number Re = 3300. We focus on the reconstruction of the flow field from the sole knowledge of the mean flow (which can be measured for instance by classical non-time-resolved PIV or by a large number of point-wise probings) and a few local unsteady measurements. This reconstruction is then used to assess the impact of experimental uncertainties on the results. Note that this article does not aim at predicting the future dynamics of the jet, since the reconstruction is performed only for the time interval of the input data. This puts the present study in stark contrast with estimation techniques such as that of Guzmán Iñigo et al. ( 2014), who were able to predict the linear dynamics of a flow from local measurements. Contrary to the present work, their model relies on a prior learning of the dynamics, based on a knowledge of the time-resolved flow field.

This chapter is divided in four main sections. The first part briefly reintroduce the reconstruction procedure described in chapter 5 ( §6.2). The second part is dedicated to a physical description of the jet, through the analysis of TR-PIV measurements, which will be considered as the reference to assess the performance of the approach ( §6.3). We then apply the procedure to rebuild the time-resolved flow field of the jet from only the mean flow and a single point-wise measurement ( §6.4). Finally, we focus on the robustness of the method and present the main guidelines to get a robust reconstruction in an experimental context ( §6.5).

Reconstruction procedure

The following briefly reintroduces the model of chapter 5, in the context of a temporal reconstruction. The unsteady velocity field U (x, t) of a given flow can be equivalently represented in the frequency domain by its Fourier modes Û (x, ω). For any given frequency ω, the Fourier mode of a flow is expected to be proportional to the dominant resolvent mode about the mean flow, given that the resolvent operator fulfill the DSV condition (see chapter 5).

From a physical point of view, this DSV property relates to the existence of a strong instability mechanism. In the case of the present jet flow, the proportionality between Fourier and dominant resolvent mode occurs for all energetic frequencies (see section 6.4.1), which may be explained by the global dynamics of the jet that is driven by the Kelvin-Helmholtz instability mechanism. Using the notations of the previous chapter, this proportionality may be formalized, for instance in an axisymmetric framework and for the axial velocity component u x , as ûx (ω, x, r) ≈ Λ(ω)ũ ω x (x, r), (6.1)

with ûx the Fourier transform of the measured axial velocity, ũω x the axial velocity component of the dominant resolvent mode at the frequency ω and Λ the amplitude function (see chapter 5). By assuming the knowledge of ûx at a given point (x 0 , r 0 ), we get an evaluation of Λ as Λ(ω) = ûx (ω, x 0 , r 0 )/ũ ω x (x 0 , r 0 ), (6.2) yielding the prediction of ûx at any point of the domain (with equation (6.1)). Then, the temporal unsteady field may be reconstructed by inverse Fourier transform. Note that, as mentioned in the previous chapter, the resolvent modes may be computed by an s.v.d. of the resolvent, or alternatively, by a PSE analysis. The latter presents several advantages (discussed in chapter 5) that makes it more interesting in the present experimental context. The final procedure that we followed is graphically illustrated in figure 6.1.

The formalism introduced here presents one strong difference with that of chapter 5, which has introduced different methods to compute a statistical representation of ûx from the simulated data. Here, since we aim at performing a temporal reconstruction, which involves an inverse Fourier transform, the Fourier mode has to come from a straightforward FFT of the time signal u x (t), without any bin segmentation nor POD processing. As shown in section 6.4.1, such FFT Fourier modes are found to be indeed proportional to the resolvent modes, which makes the reconstruction possible. The generality of this is discussed in chapter 7.

The previous explanations arbitrarily focus on the axial velocity for the description of the procedure. But the same relations hold for other flow quantities, and as shown in the previous chapter, the resulting Λ would not depend on which quantity is used. We could therefore use another velocity component, the pressure, or any other quantity that is linearly dependent on these variables (such as the wall shear stress) to compute Λ, and this Λ may then be used to reconstruct all the flow quantities. This aspect is demonstrated in the case of the round jet studied here, in section 6.4.2 where both the radial velocity and the pressure fluctuation field are reconstructed from the sole knowledge of the axial velocity at one point.

Note that Λ can be set equal to zero over a given frequency range, as this would simply result in a frequency-filtering of the reconstruction. This is generally useful, since the unsteady input measurement, used for the computation of Λ, is a discrete signal sampled in time. A Fast Fourier transform (FFT) gives its Fourier transform for a set of discrete values ω i that depends on the sampling and duration of the signal. Resolvent modes can be computed for every ω i , but even if the model could accurately reconstruct all frequencies, it is likely that a large number of them would have a very weak contribution to the dynamics (for instance the very low or very high frequencies). Setting Λ = 0 for low-energy frequencies may spare non-useful computations. This may also be necessary for flows where the DSV condition holds for a limited frequency range: one may then focus on this range only, and get an accurate filtered reconstruction of the unsteady field.

Characterisation of the application case

Experimental set-up and data processing

The experimental configuration studied in the present work corresponds to that of [START_REF] Yegavian | Lucas-kanade fluid trajectories for time-resolved piv[END_REF]. We have focused on a cold round jet at a Reynolds number Re ≈ 3300 (based on the exit diameter of the nozzle D = 12 mm, the jet exit velocity U e = 4.0 m/s and the air viscosity at T = 15 • C). From now on, all quantities are made non-dimensional by using U e and D. The flow dynamics has been characterized in a diametral plane of the jet by TR-PIV measurements, using the experimental set-up and parameters presented in [START_REF] Yegavian | Lucas-kanade fluid trajectories for time-resolved piv[END_REF]. Note that the coincidence of the laser sheet with a diametral plane has been ensured using precision devices, such that these planes can be considered nearly perfectly parallel (up to 3 • 10 -3 radians) and distant by less than 9 • 10 -3 non-dimensional length units.

In the following, as depicted in figure 6.2, the streamwise direction of the jet is denoted as x, while the cross-stream axis of the laser plane is y. The velocity has been measured in the streamwise plane (x, y), using 10000 snapshots taken at a frequency of 10 kHz (corresponding to 30 snapshots per non-dimensional time unit). Particles have been illuminated by a 2 mm thick laser sheet, using a Litron LDY303HE laser that provides an energy of 5 mJ per pulse. The snapshots have been processed using the FOLKI-PIV software, based on a classical two-frame estimation technique [START_REF] Champagnat | Fast and accurate piv computation using highly parallel iterative correlation maximization[END_REF]. We used Gaussian interrogation windows of 19×19 pixels (corresponding to a size of 0.09×0.09 in non-dimensional units) with a standard deviation σ = 4 pixels. Given the seeding density of our set-up (approximately 0.05 particle per pixel), the choice of this interrogation window size yields a good trade-off between noise and spatial resolution, in particular for the computation of the mean flow required for our reconstruction technique. This is further discussed in section 6.5.3.

Note that such a seeding density is however low for planar PIV, and thus yields noisy instantaneous velocity fields with the classical two-frame approach. However, [START_REF] Yegavian | Lucas-kanade fluid trajectories for time-resolved piv[END_REF] showed that an accurate estimation may be obtained in these conditions when using the Lucas Kanade Fluid Trajectory (LKFT) algorithm described in their paper, similar to other advanced time-resolved algorithms such as FTC [START_REF] Lynch | A high-order time-accurate interrogation method for time-resolved piv[END_REF] or FTEE Jeon et al. (2014). Figure 6.3 compares a snapshot obtained with the two approaches, and the noise appears indeed strongly reduced when using this alternative PIV processing. We therefore have two sets of velocity estimations for the present study. The two-frame noisy fields were used for the reconstruction procedure, both for the computation of the mean flow (figure 6.4) and for extracting the local input signal necessary for the reconstruction (figure 6.5).

On the other hand, the LKFT snapshots have been used as a comparison to evaluate the Fig. 6.2 Picture of the experimental set-up, displaying the round jet, the high-speed laser and camera. In the laser sheet plane, the streamwise direction is denoted x, and the cross-stream direction, y. . 6.3 Comparison between the instantaneous streamwise velocity field at t = 50 obtained from (a): classical two-frame processing and (b): the LKFT algorithm. Velocities are made non-dimensional using the jet exit velocity. The point-wise measurement and mean flow used in the reconstruction will be extracted or computed from the classical processing, therefore hereafter referred to as 'reconstruction set'. The LKFT processing will be considered as the objective to attain, i.e. the 'reference set'. quality of the reconstruction. In the following, we will refer to these sets of snapshots as the reconstruction set and the reference set, respectively. Using such a reconstruction set shows that the procedure is rather robust even in sub-optimal experimental conditions. Moreover, to our knowledge, the experimental studies based on stability theory seldom use raw PIV results. For instance, the PSE analysis performed by [START_REF] Gudmundsson | Instability wave models for the near-field fluctuations of turbulent jets[END_REF] uses a mean flow computed from a Gaussian fitting of PIV measurements. As demonstrated in the next sections, the present approach may be successfully used without any fitting of the data acquired from classical measurement techniques, standard algorithms and a rather poor experimental set-up.

Fig

Fig. 6.6 Instantaneous fluctuation field about the mean flow for an arbitrary time t = 50 of (a) the streamwise velocity and (b) the cross-stream velocity (reference set). The streamwise velocity is symmetric w.r.t. y = 0 while the cross-stream velocity is anti-symmetric, as expected from a round jet, known to be dominated by axisymmetric fluctuations.

Characterization of the unsteady behavior of the jet

The unsteady behavior of the jet near field can be qualitatively observed both in figures 6.3 and 6.6, where we see typical oscillations in the shear layer due to the Kelvin-Helmholtz instability mechanism. A quantitative characterization of the dynamics has been made by computing frequency spectra at several locations in the jet from the reference velocity field.

To obtain statistically converged spectra, the temporal velocity signals have been processed following the classical Welch's algorithm: the time series have been divided into 49 bins of 600 snapshots with a 66 % overlap, the final spectra being obtained by averaging the spectra of each bin. Close to the nozzle, they display a clear peak around a Strouhal number St = 0.76 (see figure 6.7(a)). The boundary layer momentum thickness near the nozzle has been measured as θ ≈ 0.023, and therefore, the Strouhal based on θ is around 0.009, which is consistent with the existing work in the literature for low Reynolds number jets (see for instance [START_REF] Gutmark | Preferred modes and the spreading rates of jets[END_REF]). At later stages of the development of the shear layer, the spectra become more broadband, and the energy shifts to lower frequencies. The dominant frequency becomes close to St = 0.38, corresponding to a sub-harmonic of the Kelvin-Helmholtz frequency (see figure 6.7(b)). This can be attributed to downstream vortex pairing, as explained in [START_REF] Yegavian | Lucas-kanade fluid trajectories for time-resolved piv[END_REF].

The global dynamics of round jets is known to be dominated by fluctuation modes of azimuthal wavenumbers m = 0 and m = 1, the m = 0 mode being dominant where the shear layer thickness is small with respect to the diameter (see for instance [START_REF] Davoust | Dynamics of m = 0 and m = 1 modes and of streamwise vortices in a turbulent axisymmetric mixing layer[END_REF]). In the present article, the investigated jet is measured rather near to the nozzle exit, and θ is around 0.023 near the nozzle, explaining the apparent axisymmetry of the fluctuations that can be clearly observed in the PIV snapshots (symmetric streamwise velocity field and antisymmetric cross-stream velocity field, see figure 6.6). Consequently, for the reconstruction procedure used in the rest of the study, we consider an axisymmetric framework, which has been validated a posteriori by the agreement between the PIV results and the axisymmetric PSE analysis (section 6.4.1).

The axisymmetric assumption requires to post-process the PIV velocity fields in order to accurately determine the streamwise direction and the location of the symmetry axis of the jet. To this end, the camera was carefully oriented to be approximately aligned with this axis. The small remaining misalignment has been corrected by computing, for each streamwise location x 0 , the center y c (x 0 ) of the corresponding mean streamwise velocity profile u(x 0 , y) defined as

y c (x 0 ) = min y 1 0 (u(x 0 , y + ỹ) -u(x 0 , y -ỹ)) 2 dỹ, ( 6.3)
The streamwise direction has been computed by linear regression over the computed points y c , and finally, the velocity components have been corrected with respect to this new orientation. The origin of the new system of coordinates is chosen at the center of the nozzle exit. Figure 6.4 presents the mean streamwise velocity field in this new frame of reference. Note that since the new fields are slightly rotated with respect to the original PIV snapshots, the discrete velocity values have been evaluated on a new grid centered around the axis of the jet, using a third-order spline interpolation. The new grid covers the domain (0.5 ≤ x ≤ 2.75,-1 ≤ y ≤ 1) and contains 225 × 200 points in the streamwise and cross-stream direction, respectively.

Unsteady flow field reconstruction

Prediction of the Fourier modes from the mean flow

In this section, we aim at predicting for any arbitrary frequency ω 0 the spatial structure of the velocity Fourier mode Û (x, ω 0 ) from the sole knowledge of the mean flow, computed by time-averaging the reconstruction set of PIV snapshots. As explained in section 6.2, this may be achieved by computing the dominant resolvent modes with a PSE analysis, by following the procedure detailed in chapter 2. Note that the PSE analysis has been performed in an axisymmetric framework, using the axisymmetric mean flow defined by the upper half of the complete two-dimensional PIV mean flow (see figure 6.4). In the following, we therefore switch from Cartesian coordinates (x, y) to cylindrical coordinates (x, r), and the axial and radial velocity, denoted as u x and u r , are taken equal to the streamwise velocity component u and the cross-stream velocity component v of the upper half of the PIV domain, respectively. One may alternatively use the lower half of the domain (but u r has to be taken equal to -v), but in our case, this second choice led to similar results, and therefore is not presented here.

Note that in such an axisymmetric configuration, it is also possible to use both the upper and lower part of the domain, which virtually gives twice as many snapshots to produce a mean flow that would be better converged. This may be useful when the number of available snapshots is rather low, but this was not the case for the present study, explaining why we did not use such processing.

The Fourier modes have been computed by an FFT of the reference set of PIV snapshots. Figures 6.8 and 6.9 compare the axial velocity of the dominant resolvent mode and the Fourier mode (modulus and phase comparison of the fields, respectively), for St = 0.76 (the Kelvin-Helmholtz frequency, see 6.3.1). Figures 6.10 and 6.11 show the same comparison for the radial velocity, for St = 0.38 (the Kelvin-Helmholtz sub-harmonic). For the sake of comparison, the modulus of every mode has been normalized such that its maximum is 1, and the phases such that they are equal to zero at an arbitrary location x = 2.25, r = 0.3. These figures illustrate for the two velocity components, and for two different characteristic frequencies of the flow, that, as claimed in section 6.2, the dominant resolvent mode is approximately proportional to the Fourier mode (their modulus is approximately equal up to a multiplicative constant and their phase, up to an additive constant). However this agreement strongly deteriorates in low-energy parts of the flow. This is particularly striking in figure 6.11 for instance. We see that close to r = 0, as well as for x < 1.7, the two fields present a strong discrepancy. While it is known that the dominant modes may not reproduce well the actual dynamics in the low-energy parts of the flow (see chapter 5), here the discrepancy could also be attributed to the difficulty of measuring a signal with such a low fluctuating energy. Fluctuations in these regions are instead dominated by measurement noise, which makes a relevant comparison impossible. Similar overall agreement has been found for all the other investigated frequencies, both for the axial and radial velocity components. This agreement justifies a posteriori the validity of the axisymmetric hypothesis. Note that here, contrary to the turbulent backward facing step case of chapter 5, we do not need to use POD-filtering to remove any uncorrelated turbulent fluctuations from the FFT, and still get an accurate agreement with the PSE fields. This is due to the moderate Reynolds number of the jet, which leads to a highly organized behavior. For turbulent flows, as seen previously, the situation is more complex, and this is further discussed in the general conclusion of the thesis (chapter 7).

The PSE analysis also yields the prediction of the pressure modes, for which we do not have any experimental comparison. Figure 6.12 presents the real part of the pressure mode for the Kelvin-Helmholtz frequency St = 0.76 and the sub-harmonic St = 0.38, where one can see alternated positive and negative pressure regions along the shear layer, as classically observed in shear layers. We also observe that the downstream structure that appears in figure 6.12(b) is approximately twice as large as the structures of figure 6.12(a), consistently with the assumption that frequency St = 0.38 is related to downstream vortex pairing.

Time-resolved reconstruction of the snapshots

Following the procedure of section 6.2, we have computed the amplitude function Λ using the axial velocity u x at x 0 = 2.5, r 0 = 0.4 (the input signal can be seen figure in 6.5). This choice of input point is discussed in more detail in section 6.5.1. The choice of the axial velocity as the input quantity was arbitrary; using the radial velocity led to similar results (not presented here). The snapshots are finally reconstructed by performing an inverse FFT at every spatial point. The Λ function has been computed for 0.2 < St < 1.4, which contains most of the energy of the flow, see figure 6.7, following relation (6.2). Outside of this range, we set Λ = 0, which filters the low-energy part of the spectra. In this frequency range, we computed 400 resolvent modes. This number stems from the frequency resolution of the FFT of the input signal, a mode being computed for every frequency within the range considered. Note that the final reconstructed set of snapshots has therefore the same time-sampling and duration as the input signal used for the computation of Λ.

The reconstructed snapshots display a good agreement with the reference set, as can be seen in figure 6.13, which compares an axial velocity snapshot with the reference PIV field. The size and the location of energetic structures are well reconstructed. The agreement is also good in lesser-energy locations. This can be seen in figure 6.14, which compares the temporal evolution of the axial velocity from the reconstruction and the reference set at the point (x = 2.0, r = 0). The oscillations of both signals are well in phase, and their amplitude is very close. Note that the comparison cannot be made in very low-energy parts of the flow, such as the near-axis and most upstream zones in figure 6.11. This stems from the fact that the reference set is not accurate in such regions (the signal-to-noise ratio is low) and the model is not designed to reconstruct low energy behavior anyway (see chapter 5).

In order to provide a more quantitative comparison, we introduce an instantaneous global measure of the error, denoted as e(t), defined as:

e(t) = 1 N p Np i=1 u i x (t) -u i x,0 (t) 2 , (6.4)
with N p the number of discrete points where the velocity is known, u i x and u i x,0 the i-th discrete streamwise velocity value of the reconstructed field and of the reference field, respectively. This quantity corresponds to the root mean square over the whole domain of the velocity Fig. 6.14 Comparison of the axial velocity computed from the TR-PIV snapshots (reference set, dashed line) and the reconstructed signal (red continuous line) at x = 2.0, r = 0 (a point presenting an intermediate level of energy), for an arbitrary time range 50 < t < 100. Fig. 6.15 Time evolution of the global axial velocity error (continuous line) for an arbitrary temporal range 50 < t < 100. The dashed line represents the mean error over the full time range (0 < t < 333), which is approximately equal to 0.038. error, expressed in non-dimensional unit. One should note that both the reconstruction errors and the PIV measurement errors contribute to the value of e(t), therefore it would not be zero even in the case of a perfect reconstruction. Figure 6.15 shows the temporal evolution of this error, and we see that the discrepancy between the reference and the reconstruction does not vary much with time (around 0.04), ensuring that the global quality of the reconstruction is approximately constant over the full time range considered.

As explained in section 6.2, the amplitude function Λ computed from only one flow variable (here the axial velocity) yields the reconstruction of all other fluctuating quantities (pressure and radial component of velocity). Figure 6.16 compares a reconstructed radial velocity snapshot with the corresponding reference field, and the agreement is once again favorable. Figure 6.17 shows the temporal evolution of e, and the level of error is again rather steady over time, with values close to that of the axial velocity (approximately 0.03). For the pressure reconstruction, we do not have experimental results to serve as reference, but the accurate reconstruction of the radial velocity is a strong argument in favor of the quality of the pressure reconstruction. Indeed, nothing distinguishes these two variables in our approach. Moreover, the resulting pressure fluctuation field p ′ (x, r) is reminiscent of what is expected in such a jet (see figure 6.18(a)), with alternated positive and negative pressure regions that grow in size and amplitude when moving downstream. Fig. 6.16 Comparison of the radial velocity of (a) the reconstructed field and (b) the reference PIV fields, at t = 75. Fig. 6.17 Time evolution of the global radial velocity error (continuous line) for an arbitrary temporal range 50 < t < 100. The dashed line represents the mean error over the full time range (0 < t < 333), which is approximately equal to 0.031. Another remarkable feature of the reconstructed fields is their smoothness, that makes them easily differentiable. This is of high importance for the computation of derived quantities such as the vorticity, which is sometimes difficult to accurately compute from PIV measurements, especially in a time-resolved context where the fields display a stronger noise. As an illustration, figure 6.19 compares a vorticity snapshot computed from the PSE reconstruction, the reference PIV set and the reconstruction PIV set. The fields have been obtained from differentiation of the velocity snapshots, based on second-order centred finite differences, with a stencil length equal to 0.08 (twice the length of the interrogation window used for the PIV estimation). The reconstructed vorticity compares very favorably with the reference. The vorticity directly derived from the PIV reconstruction set illustrates the kind of results that are obtained from PIV when the level of noise is too high. This strongly deteriorates the estimation of the derivatives, while the same level of noise in the input data for the PSE reconstruction has no similar impact on the reconstructed field. Note that the vorticity is one among many other derived quantities that may be determined from the present reconstruction. For instance, figure 6.18(b) shows the mean square pressure fluctuation field p ′2 , and any mean square fluctuation or Reynolds stress term could be similarly computed, and exploited for further physical analysis.

Robustness of the reconstruction method

Influence of the choice of input point

The reconstruction model of this work is based on that presented in chapter 5, which has explained that the point used for the computation of the function Λ should be in a high-energy region of the flow to yield an accurate reconstruction. This requirement stems from the fact that the dominant resolvent modes reproduce accurately the Fourier modes in these regions, but we explained that there is no guarantee that this agreement would be as good in low-energy regions. This is confirmed in our case, and it can be seen in figures 6.8 to 6.11, where we see that the strongest discrepancy between dominant resolvent modes and Fourier modes appears in low-energy regions. But as mentioned in section 6.4.1, in our case, this could be attributed to a poor signal-to-noise ratio in low-energy regions of the flow, a common problem in experimental studies. This issue can be clearly seen in figures 6.9(a) and 6.11(a) for instance, where the non-energetic upstream part of the domain is strongly noised, yielding to a significant discrepancy with the dominant resolvent fields.

To assess the sensitivity of the results with respect to the input point, we have considered seven additional points along the shear layer for the determination of Λ (in total, four points above the shear layer, four points below, see figure 6.20 for the position of the points). The reconstruction is based on the streamwise component of the velocity only, as was done in section 6.4.2. For each of these points, we have computed the global error E defined as

E = (1/T max ) Tmax 0 e(t)dt (6.5)
with T max = 333 the duration of the TR-PIV acquisition, and e the instantaneous error (see equation (6.4)). The resulting values can be seen in table 6.1, and we observe two clear tendencies: the points located downstream yield a smaller error than the ones located upstream, and the points above the shear layer display a larger error than the ones below. This is fully consistent with the findings of chapter 5, which explained that input points located in high-energy regions yields smaller error. Indeed, here the downstream region contains more energy than the upstream one, where the Kelvin-Helmholtz instability has not fully developed yet. Also, concerning the axial velocity, the region inside the jet has been found more energetic than the region outside. This can be observed for instance in figure 6.8 for St = 0.76.

For a given flow configuration, the location of high energy region may usually be determined by a prior physical knowledge of the qualitative dynamics of the flow. It may therefore be easy to predict good locations for the input sensor in most situations. However, it is likely that the flow dynamics would involve several characteristic frequencies, or even a range of frequencies, and that these frequencies would be related to different locations of the flow. In such a case, it is not possible to find a single location that would be energetic for all these frequencies. As explained in the previous chapter, this may be treated by considering n points for the determination of Λ, leading to an over-determined set of n equations

{Λ(ω) = ûx (ω, x k , r k )/ũ ω x (x k , r k ), 1 ≤ k ≤ n}
that may be solved by least squares. The resulting Λ would be weakly impacted by low-energy points, such that for every frequency, only high-energy points would contribute to its value. The reconstruction would also be likely to be more accurate because it would be based on more input data. Finally, multiplying the number of input points may be useful when there is not any prior knowledge of the energetic regions, or when this knowledge is not accurate enough, since only one of the sensor needs to be well-positioned.

We have tested this approach by considering the best and worst points together (A and H, see table 6.1). The result, close to that obtained with A only, shows that the reconstruction is not degraded by the poorly chosen point H. Therefore, in an experimental context, an accurate prediction would be achieved by using several point-wise measurements as long as at least one of them is located in a high-energy part of the flow for each frequency of interest. We also tested the four points above the shear layer (E to H) together, as well as all points together, and table 6.1 shows that when Λ is based on a set of points, the results appear always more accurate than the single-point reconstruction based on the best point of the set. This is of high importance in an experimental context, when using a sensor at some specific high-energy points may be too intrusive (for example at points A to D). This gives more flexibility to the method, since lower-energy points might also be used, as long as they are sufficiently numerous. In particular, one may position a rather large number of sensors downstream from the region of interest and get a satisfactory reconstruction, even if this region is not the most energetic one.

Impacts of an inaccurate knowledge of the input sensor position

The present procedure aims at being used with unsteady data obtained from point-wise sensors. However, in practice, the spatial positions of these sensors can be subject to some uncertainties. In this section, the impact of such an inaccurate knowledge of the input sensor position on the quality of the unsteady reconstruction is studied by considering an input sensor positioned at (x 0 + δx, r 0 + δr), but which would be erroneously assumed to be located at (x 0 , r 0 ). In the procedure, this amounts to select a different input signal than that of figure 6.5, such that the computed amplitude function is no longer defined by equation (6.2), but becomes equal to the following biased amplitude function

Λ b = ûx (ω, x 0 + δx, r 0 + δr) ũω x (x 0 , r 0 ) . (6.6)
For the reconstruction, we found that such a misplacement mainly results in a streamwise translation of the structures of the fields, and to a lesser extent, to a moderate change in their overall amplitude and shape. This overall effect on the reconstruction can be seen in figure 6.21 for δx = δr = 0.05, δx = δr = 0.1 and δx = δr = -0.2 (Λ is computed from the axial velocity and (x 0 , r 0 ) = (2.5, 0.4)). In physical units, this corresponds respectively to a misplacement of 0.85, 1.7 and 3.4 mm, i.e. from realistic up to over-estimated experimental positioning errors. When compared with an unbiased reference snapshot (figure 6.21(a)), these biased snapshots seems to all present a phase shift. Beside this dephasing, the resulting field displays the expected physical features (alternated positive and negative structures along the shear layer that grow in size and amplitude in the streamwise direction). This demonstrates that despite a possibly significant misplacement, the reconstruction does not degenerate but rather keeps a certain physical relevancy. This may be explained by the fact that here we consider flows displaying coherent structures, which present by definition a strong spatial correlation. Therefore, the energy content of the frequency spectrum will not suddenly change by considering a point slightly misplaced (the spectrum will display a similar shape with the same dominant frequencies). The phase of the spectrum may however change significantly, which explains why the main observable effect of a misplacement is a dephasing of the fields.

Despite this weak overall impact on the physical features of the reconstruction, it may still be interesting to minimize the sensitivity with respect to a misplacement by a relevant choice of input data. To address this issue, let us decompose the biased amplitude function Λ b as

Λ b = Λ + δΛ, (6.7)
with Λ = ûx (ω, x 0 , r 0 )/ũ ω x (x 0 , r 0 ) the unbiased amplitude function defined by equation (6.2), and δΛ a spurious term that reads δΛ = δ ûx ũω

x (x 0 , r 0 ) , (6.8)

with δ ûx = ûx (ω, x 0 + δx, r 0 + δr) -ûx (ω, x 0 , r 0 ). For the reconstruction to be weakly impacted by the misplacement of the sensor, δΛ has to be small with respect to Λ, which would ensure that each Fourier mode is rebuilt with a good accuracy. The spurious effect of the misplacement may therefore be evaluated by the ratio b = |δΛ/Λ| = |δ ûx /û x |. Finally, using the proportionality between the Fourier modes and the dominant resolvent modes yields the following expression of b b = ũω x (x 0 + δx, r 0 + δr) -ũω x (x 0 , r 0 ) ũω

x (x 0 , r 0 ) . (6.9)

This expression only involves the dominant resolvent modes, and may therefore be computed by only knowing the mean flow. Assuming that the misplacement (δx, δr) is small, b may be linearized as: b = |b x (ω, x 0 , r 0 )δx + b r (ω, x 0 , r 0 )δr|, (6.10) with b x and b r respectively the axial and radial sensitivity coefficient defined as (6.11) Small values of |b x | and |b r | correspond to a low sensitivity with respect to misplacements of the sensors, but the reciprocal is not true, as b may be small even for large values of |b x | and |b r | (errors along x and r may compensate each other). In addition, the value of this coefficient does not give any indication regarding the type of effect of a misplacement, which could be a simple phase shift (as evidenced in figure 6.21), or a more significant distortion of the fields. Therefore, this coefficient only gives qualitative guidelines to minimize the impact of misplacements.

b x = ∂ x ũω x /ũ ω x , b r = ∂ r ũω x /ũ ω x .
First, it should be noted that these coefficients depend on the input physical quantity considered: here b x and b r are defined from u x , but using for instance the radial velocity u r as input data for the reconstruction yields b x = ∂ x ũω r /ũ ω r and b r = ∂ r ũω r /ũ ω r . Consequently, the sensitivity with respect to misplacements may be reduced by considering an input physical quantity whose Fourier modes do not present strong spatial gradients. In the case of the jet, while the axial velocity modes display some abrupt variations across the shear layer, as it may be seen for instance in figure 6.9, the radial velocity or pressure modes do not display such high-gradient regions (see figures 6.10 and 6.12). The reconstruction may be more robust if based on a local record of one of these two quantities, especially if the sensors are designed to be near the shear layer. Equations (6.10;6.11) give also insight about the best locations for the input sensors. Due to the division by the local amplitude of the resolvent modes in equation (6.11), b x and b r are expected to be small in high-energy regions, as long as the local spatial gradients do not become too strong. Moreover, the spatial gradients are expected to be rather small in the direct neighborhood of a local energy maxima. Figure 6.22 compares the quantity |b x | + |b r | (computed from the radial velocity modes) and |ũ ω r | for St = 0.76, and we see that high-energy regions are indeed rather weakly sensitive. While the considerations related to the sensors misplacements completely differ from that of section 6.5.1, we are here led to a similar conclusion: the input sensors should be preferably positioned in energetic regions.

Finally, this rather good correspondence between high-energy and low-sensitivity points is interesting if multiple input points are considered for the computation of Λ. The procedure that is used for the determination of Λ is then based on a least squares minimization that rules out the low-energy points, for each single frequency of the reconstruction (see section 6.5.1). It is therefore expected to also rule out the high-sensitivity points, which increases the overall robustness of the method.

Sensitivity with respect to the mean flow measurements

The PIV measurements need to be accurate enough to yield a proper mean flow. In addition to the experimental conditions, the key parameter that influences the quality of the PIV velocity estimation is the size of the interrogation window: a large interrogation window contains more particles, which reduces the measurement noise, but it tends to smooth down the spatial gradients of the flow.

In the previous sections, we have used a mean flow that yields an accurate reconstruction, confirming a posteriori that the choice of the size of the interrogation windows was appropriate (see section 6.3.1). Choosing a smaller window increases the level of noise, but this does not Fig. 6.23 Comparison between mean flow A (continuous line) and mean flow B (dashed line), the latter corresponding to a PIV processing with a larger interrogation window. Figures (a) and (b) respectively display u x and ∂u x /∂r at x = 2, where u x is the mean axial velocity.

have any significant impact on the final mean flow, as long as the number of snapshots used for the time-averaging is high enough for convergence to be reached. However, taking an overly large interrogation window yields a mean flow that displays inaccurate, biased spatial gradients. This may be observed in figure 6.23, which compares the original mean flow (mean flow A) with a new mean flow (mean flow B), obtained with interrogation windows twice as large in each direction (Gaussian window, 37×37 pixels, σ = 8). The overall effect of this enlargement is an under-estimation of the spatial derivatives of the mean flow. For St = 0.38 (the dominant frequency in the downstream zone of the flow), the dominant resolvent modes computed from mean flow A or B have been found hardly distinguishable. However, as we go to higher frequencies, an increasingly strong discrepancy appears when considering one mean flow or another. This discrepancy is therefore particularly important for high frequencies, as it can be seen for instance in figure 6.24: for the Kelvin-Helmholtz frequency St = 0.76, the mode computed from mean flow B is abnormally energetic in the upstream region of the jet. The impact on the final reconstruction may be observed in figure 6.25 (reconstruction based on the axial velocity at (x = 2.5,r = 0.4)): while the large low-frequency structures are correctly reconstructed downstream, some spurious high-frequency structures appear upstream of the flow.

These results show that the quality of the reconstruction is conditioned by the accuracy of the mean flow measurement, and one crucial aspect is the correct evaluation of the spatial gradients of the mean flow. In our case, it is found that when these gradients are erroneously evaluated, the final reconstruction presents abnormal levels of energy for high-frequency structures only. Note that this kind of preoccupation mainly concerns configurations where the seeding density would be particularly low, or where some parts of the flow contain really few particles. These two situations would be the only ones justifying the use of windows so large that it would have a significant impact on the PSE analysis.

Concluding remarks

This study shows that, in the case of a round jet at Re = 3300, the sole knowledge of the mean flow and the unsteady behavior of one velocity component at a given point is enough to yield a reconstruction of all the variables of the flow field, including the pressure fluctuations. The reconstruction procedure relies on the ability of a PSE analysis on the mean flow to yield the dominant resolvent modes, which gives an approximation of the spatial structure of the Fourier modes of the flow. Such an analysis is computationally low-demanding and easy to implement. It is also particularly well adapted to the study of experimental configurations since it can be used even when the mean flow is known on a rather small region.

The reconstruction quality is conditioned by the choice of the input measurement location: for the reconstruction to be accurate, it has to be located in a rather energetic area of the flow. Moreover, high-energy locations are likely to be robust with respect to a small misplacement of the sensors (but in our case the overall effect of a misplacement was rather limited). In some situations, determining a priori these optimal locations may be difficult. Hopefully, as demonstrated in the last section of the article, this issue may be solved by increasing the number of inputs: the reconstruction technique is then more robust and accurate even if the sensors are not optimally positioned. This gives more flexibility for the location of the inputs. For instance, with enough sensors, one may avoid intrusive regions, and place all the inputs downstream of the flow. The accuracy of the reconstruction is also naturally related to the quality of the mean flow measurement. But the present study proves that a sufficiently accurate mean flow may be obtained from classical two-frame PIV even in the case of difficult experimental conditions (e.g., here, with a low seeding density). The noise present in the PIV snapshots should not impact the reconstruction, since it cancels out in the mean flow as soon as enough snapshots have been acquired and convergence is achieved. Standard PIV parameters yield, in our case, a good reconstruction, and the main precaution that emerged concerns the interrogation window used for the PIV processing, as it may lead to bias present in the mean flow. If it is excessively large, then the final mean flow may present erroneous spatial derivatives, which had in our reconstruction an impact on high-frequency structures.

The present work is a successful example of the use of our reconstruction method, but it does not ensure the generality of this approach. For higher Reynolds number flows, while chapter 5 has showed that the statistical spectra associated to the coherent part of the unsteadiness may be predicted, a temporal reconstruction would be a more challenging task, that is discussed in the general conclusion of the thesis (chapter 7).

Note that here, we have only reconstructed the m = 0 fluctuations, that are dominant for the present jet. In the case where an m = 1 mode would be of interest, then the PSE analysis could be performed for m = 1 disturbances, to get the associated resolvent mode. Then, the only requirement would be to separate, in the unsteady point-wise data, the m = 0 and m = 1 contributions. This may be done by considering a circular array of sensors, and using classical azimuthal decomposition techniques (see for instance [START_REF] Kopiev | Azimuthal decomposition of turbulent jet noise and its role for diagnostics of noise sources[END_REF]). Then, the m = 0 contribution of the input data would serve to compute an amplitude function Λ 0 for the axisymmetric PSE mode, and the m = 1 one, for the computation of its counterpart Λ 1 for the m = 1 mode. This approach may be generalized to other situations where several modes would come into play.

Part III

General conclusion and future work

Chapter 7 -Conclusion

Summary of the work

The purpose of this thesis was to develop reliable tools for the characterization of a general flow dynamics, based on linear stability theory. A large body of work has already focused on that goal, which gave birth to a manifold of linear stability methods. Unfortunately, as showed in chapter 1, none of these linear approaches has the ability to fully characterize any general flow. Instead, each one presents its own limitations, and a clear understanding of these limitations is still missing.

On account of this, the present work addressed the conditions of validity of these methods. To this end, chapter 2 provided efficient implementation strategies for the main linear stability techniques. Then, the work was carried out in two steps, by studying successively the base flow and the mean flow approach.

The base flow approach was first addressed in chapter 3, by focusing on a global stability analysis of screeching jets. This revealed unstable modes, with a frequency close to that of screech tones. Moreover, the stability analysis correctly accounted for the influence of two key parameters, namely the jet-to-ambient pressure ratio and the nozzle lip thickness. Finally, a careful study of the unstable modes revealed supersonic disturbances, which is a well-known noise generation mechanism that is responsible for acoustic waves emission in the jet.

These successful results served as a discussion basis in chapter 4, which addressed the general ability of a linear base flow analysis to predict the nonlinear frequency of a flow. The main finding was that for oscillating flows, when the frequency selection process originates from an acoustic feedback loop (encountered for instance in impinging jets, cavity flows, screeching jets), then the nonlinearities weakly impact the dynamics. In this situation, a linear base flow analysis would be physically relevant even far from criticality. More specifically, the eigenspectrum of the Jacobian is expected to fulfill the so-called RPIF (real positive imaginary frequency) property, which refers to the existence of an eigenvalue whose real part is positive (unstable eigenvalue) and imaginary part matches the nonlinear frequency of the flow.

In the second part of the manuscript, we investigated the mean flow approach, with a focus on resolvent analyses. From the formalism introduced in chapter 5, the Fourier mode of a flow was found equal to the product of the resolvent operator by a turbulent forcing term. We then considered a singular value decomposition of this operator. This yielded a set of optimal responses and optimal forcings, which were used as projection basis for the Fourier mode and the turbulent forcing, respectively. The resulting relation may then be simplified by assuming that the first singular value is much larger than all other. From a physical point of view, this dominant singular value (DSV) assumption relates to the existence of a strong convective instability in a flow. Under this condition, the first optimal response was found to be approximately proportional to the Fourier mode. Note that for this to be true, the turbulent forcing should not be preferentially oriented toward sub-optimal forcings. This was the only modeling assumption of our study. These findings, that addressed the prediction of spatial features of an unsteady flow, are consistent with the work of [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 [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF], who focused on the prediction of the frequency of self-oscillating flows. They showed that for a (quasi)monochromatic dynamics, the mean flow verifies the so-called RZIF property (real zero imaginary frequency), i.e. a global stability analysis of this field yields a marginally stable mode whose frequency matches the frequency of the flow. In a resolvent framework, such monochromatic oscillations could only appear if the first singular value tends to infinity in the vicinity of this frequency, which is equivalent to the RZIF property. In addition, the DSV condition would then be automatically fulfilled for the frequency of the flow oscillations. Therefore, a mean flow analysis would predict not only the frequency, but spatial features of the self-oscillating behavior too.

In the more general case where the RZIF property is not fulfilled, there are two options to predict the frequency spectrum at any point: one needs to either model the turbulent forcing term, or use additional information about the dynamics. On account of that, we introduced two models. The first one considered a modeling of the forcing as a response-independent white noise. The second one was based on a few point-wise unsteady measurements from the flow. The models were tested in the case of a turbulent backward facing step simulation, which was found to fulfill the DSV condition for all considered frequencies. The results from the first model were qualitatively acceptable, but quantitatively inaccurate. This was due to the overly-simple modeling of the turbulent forcing. On the other hand, the second model gave excellent results, and proved to well predict the dominant frequency and the level of energy at all the tested locations. Finally, we showed that in the backward facing step case, the optimal response may be well-approximated by a PSE analysis. This gives an inexpensive alternative that may replace the s.v.d. for both models.

The second model, based on point-wise measurements, is very interesting for experimental studies. Indeed, the mean flow and point-wise data are generally much easier and much cheaper to measure than the unsteady field over a whole region. Therefore, in chapter 6, we focused on the use of this model in an experimental context, and we reconstructed the unsteady flow field of a transitional round jet from the sole knowledge of the mean flow and a single point-wise measurement. This reconstruction favorably compared with an experimental reference, given by TR-PIV measurements. Finally, we addressed the robustness of the model with respect to experimental uncertainties. We gave a few guidelines to ensure robust results, and the main finding was that the measurements should be in high-energy regions. Moreover, using a few point-wise measurements instead of a single one increases both the accuracy and the robustness of the reconstruction.

Meaning of a mean flow stability analysis

The present work showed that a mean flow stability analysis has a clear physical meaning from the resolvent point of view. From this standpoint, and despite its name, a mean flow stability analysis is not a stability analysis: the stability of the mean flow has no physical meaning, and such an analysis cannot yield any kind of information with respect to the onset of an unsteadiness. However, this work pinpointed the central role of the resolvent about the mean flow in the nonlinear dynamics of a flow. For flows displaying a strong convective instability mechanism, this operator is expected to fulfill the DSV property. This leads to a particular situation where the flow would display high-energy spatially correlated structures that can be easily computed from the resolvent. Regarding the frequency prediction, this cannot be achieved by a classical mean flow analysis in a general case, and it requires to properly model the turbulent forcing term, or to use additional information about the dynamics (such as point-wise measurements).

However, in the particular case of monochromatic and quasi-monochromatic flow oscillations, the situation is different: for the system to display a consistent behavior, the resolvent has to degenerate such that the response is infinite for the frequency of the flow. We showed that this yields a very particular situation where the frequency may be predicted without any modeling of this turbulent forcing term, by performing a simple global stability analysis about the mean flow.

Temporal reconstruction and convergence of spectra

Chapter 6 has demonstrated that the techniques developed in this thesis may yield an accurate flow reconstruction. The studied case was a jet at Re = 3300. Such a flow is often classified as a noise amplifier, in the sense that its dynamics is rather broadband, and does not stem from self-sustained oscillations. This may be observed with the spectra presented in the previous chapter for the characterization of the jet (figure 6.7). These spectra were obtained by the classical Welch algorithm to compute the power spectral density. This is a common way to get a converged spectral information from non-periodic flows, that cannot be obtained from a straightforward FFT, since the signal is not square integrable. However, the reconstruction itself is based on a straightforward FFT. But the non-convergence of the local spectra is not an issue to rebuild a signal. Indeed, an FFT of a sampled signal u(x, y, t i ) gives spatial Fourier modes ûn (x, y) for a discrete set of frequencies (ω n ), such that we have:

u(x, y, t i ) = N n=0 ûn (x, y)e iωnt i + c.c, ( 7.1)
valid for all discrete times t i of the sampled signal. The spectrum of an FFT at a given point gives the relative amplitude and phase of all the FFT modes. Regardless of the convergence of this spectrum, it still gives the exact linear combination of these modes that reproduces the original signal. If the length of the signal increases, the number of modes increases too, and therefore the phase and amplitude of the modes have to change to reproduce the extended signal. As a consequence, the spectra obtained from an FFT does not converge. Nonetheless, the finite-length signal is always perfectly reproduced by this superposition of Fourier modes, by construction.

Therefore, in the case where the resolvent modes would be perfectly proportional to the FFT Fourier mode, then the flow field may also be perfectly reconstructed from these modes. This would be true without any consideration regarding the local convergence of the spectra. The real issue is therefore the convergence of the Fourier modes, obtained from an FFT, rather than the local spectra. In the case of the jet, this convergence occurs, and the reconstruction is therefore possible (see chapter 6).

But the jet is a moderate Reynolds number flow, and for turbulent flows, the question is more complex. In the case of the turbulent backward facing step, the most energetic behavior was associated to ω = 0.7. The corresponding Fourier mode, computed from a simple FFT, yields highly organized structures, that is easily converged with a simple FFT, as long as the input signal is long enough to account for this low-frequency behavior. For this frequency, the resolvent mode is very close to this FFT mode (see figure 7.1), and therefore the situation is similar to that of the jet, where a reconstruction should be possible.

At a higher, less-energetic frequency, the situation is different. The Fourier mode, computed from an FFT, displays a spatially organized region, that is well converged when the length is increased, but the downstream part presents an apparent noise, that does not converged as the length of the input signal is changed (figure 7.2). Similarly to what is shown in chapter 5, this noisy part is absent from the resolvent mode. The nature of these fluctuations is still unknown, and may be further studied in future works.

Unsteady reconstruction of turbulent flows

The previous results show that only the reconstruction of highly-organized behavior is possible from the rank-one approach developed in part II. For non-turbulent flows, this yields a straightforward reconstruction, similar to that of chapter 6. For turbulent flows, the previous findings show that some high-energy behavior may correspond to such spatially correlated structures, which can therefore be accurately reconstructed. However, some lower-energy behaviors, typically associated to higher frequencies, can only be partially reconstructed. Indeed, at least in the case of the backward facing step, above a certain frequency, the Fourier modes display a coherent behavior superimposed to some uncorrelated fluctuations. Only the former can be reconstructed, but such a partial reconstruction would require to isolate, in the point-wise input data, this coherent part. This could be done by performing a proper orthogonal decomposition of the input point-wise signals, similarly to what has been done by [START_REF] Gudmundsson | Instability wave models for the near-field fluctuations of turbulent jets[END_REF], but the number of input sensors should be high enough for such a technique to be used.

Alternatively, one could perform a reconstruction based on a rank-n model, using the suboptimal resolvent modes. In theory, with enough modes, the reconstruction can be arbitrarily accurate, but the cost is, again, a higher number of input sensors, in order to compute an amplitude function Λ i for each of these sub-optimal modes.

Future work

The following section presents several ideas of further developments, for the base flow and the mean flow approach, respectively.

Base flow approach

Chapters 3 and 4 defined some limits of validity for a linear base flow analysis. We pinpointed the cases where such an analysis may yield relevant results, even far from criticality. These chapters were based on a stability analysis of underexpanded jets, which revealed that the screech phenomenon originates from a global instability of the jets. This result is very interesting from a physical point of view, especially in the context of control, but the aim of the chapter was such that we did not explore this aspect. Nonetheless, based on the existing tools described in chapter 2, it would be rather straightforward to perform a sensitivity analysis of the unstable modes with respect to an external forcing (see [START_REF] Mettot | Stabilité linéaire, sensibilité et contrôle passif d'écoulements turbulents par différences finies[END_REF]). This would bring new results for the problem of screech suppression, that has been intensively studied during the past decades (see for instance the work of [START_REF] Norum | Screech suppression in supersonic jets[END_REF] or more recently the article of [START_REF] Ramakrishnan | Screech suppression of supersonic jet noise[END_REF]).

Mean flow approach

The present study led to conditions for the validity of a mean flow stability analysis, in the case where the flow dynamics is governed by the Navier-Stokes equations. Since the approach is rather general, an interesting development would be the extension to an augmented set of equations, that would include equations from other physical or chemical mechanisms. From this, one may then study reactive fluids, for instance. It would be also interesting to see if the formalism introduced in chapter 5 could be extended to multi-phased flows, or even coupled solid/fluid systems. In both cases, the notion of mean field needs to be addressed, in particular in region where the nature of the medium changes during time. A shared mathematical framework is also needed to describe both media, in particular for fluid/solid systems, which are classically described from an Eulerian and a Lagrangian point of view, respectively. Moreover, the state of each medium may be described by different variables, which adds to the complexity of the problem.

Other developments may be done regarding the reconstruction model described in chapter 5.

In the present work, we only considered 2-D data. But in the case where one would be able to produce a three-dimensional mean flow measurement, then the exact same procedure may be used to rebuild a full 3-D time-resolved flow field. This could be of high interest since the existing methods to measure a flow field in a volume are more complex and costly than their 2-D counterparts. Note that several flows, such as the backward facing step of chapter 5, have a 3-D dynamics but a 2-D mean flow. Therefore, this mean flow could be measured from 2-D PIV, and a few point-wise measurements may yield the full 3-D unsteady flow field.

In chapters 5 and 6, the PSE analysis is presented as an inexpensive alternative to an s.v.d. of the resolvent, that may be useful to treat 3-D configurations (too expensive for an s.v.d. with the numerical tools presented in this thesis). But it can also be argued that the PSE analysis presents some limitations that a resolvent analysis do not have. For instance, acoustic waves are often impossible to properly capture by PSE (see chapter 1). Other spatial marching methods, that may be superior to PSE in some situations, could be used, such as the one-way equations introduced by [START_REF] Towne | One-way spatial integration of hyperbolic equations[END_REF]. This may replace the PSE analysis in the present reconstruction model if needed. Its computational cost is between that of a PSE analysis and an s.v.d. of the resolvent.

The reconstruction method is based on a rank-one approximation of the resolvent operator, that is only valid if the DSV property is fulfilled. But if that is not the case, then a higher-rank approximation may be used. This requires to compute sub-optimal responses, which cannot be obtained from a PSE analysis. Therefore, it would be interesting to modify the PSE procedure, and to use the adjoint PSE equations such that, from a random initialization, we compute the optimal response by adjoint optimization. We may then compute sub-optimal responses by imposing this optimization procedure in a sub-space that would exclude the first optimal response. With such a technique, one may be able to perform a higher-rank reconstruction without the need of a costly resolvent analysis.

Note that in chapters 5 and 6, the reconstruction is performed in the frequency domain. Gòmez et al. (2016a) performed a similar reconstruction in the time domain. They assumed that the flow field is equal to a Fourier series involving the modes computed with the resolvent. The coefficients of the series are then determined by minimizing the error at at few points, from a few point-wise measurements. An alternative approach would be to minimize the error from a PIV snapshot instead of a point-wise measurement: instead of using an information at one point for a large set of different times, we use the information at one time for a large number of points. We may then use classical (non time-resolved) PIV, which yields the mean flow and a set of snapshots. The minimization can be performed over the whole set of snapshot to make the reconstruction more robust, which gives a time-resolved field from a non-time-resolved set of snapshots. Note that from a physical point of view, a snapshot contains structures of various sizes, associated to various frequencies. Therefore, a single snapshot contains phase and amplitude information about all the energetic frequencies. In the approach developed in this thesis, the snapshots are only used for the computation of the mean flow, and therefore a lot of information is lost in the process. It may be also possible to use both point-wise measurements and the snapshots to increase the robustness of the reconstruction.

This snapshot-based approach could also be used with Particle Image Tracking (PTV, see for instance [START_REF] Maas | Particle tracking velocimetry in three-dimensional flows[END_REF]) instead of PIV. PTV is based on a tracking of individual particles, and yields an instantaneous velocity measurement at the location of these particles, that may be more accurate than PIV measurements. However, this velocity measurement is usually sparser than with PIV, as the density of particles needs to be low enough to get an accurate tracking. Following an approach similar to that mentioned above, PTV results may be used to reconstruct a dense time-resolved velocity field.

Finally, one may attempt to reconstruct a flow field from the sole knowledge of a mean flow, without any other extra data about the dynamics. To this end, the modeling of the turbulent forcing term should be further analyzed. In this thesis, the only modeling that we proposed was extremely simple, and could only recover some qualitative trends of the unsteady dynamics (we could only identify high frequency or low frequency regions). Undoubtedly, a better modeling would yield an improved reconstruction, but this requires to further study this turbulent term.

A.2.1 Cas d'étude: le phénomène de screech

Les jets supersoniques sont qualifiés de sous-détendus lorsque la pression à la sortie de la tuyère est supérieure à celle du milieu ambiant. Le ratio de pression jet/milieu ambiant (JPR, Jet Pressure Ratio) caractérise la structure spatiale du jet, qui s'organise en une succession de cellules de chocs (voir figure A.2).

Sous certaines conditions, un bouclage acoustique peut se former: des perturbations naissent dans la couche de mélange, et sont convectées vers les chocs où elles interagissent et génèrent des ondes acoustiques. Ces dernières remontent l'écoulement, jusqu'à être reflétées par la lèvre de la tuyère, ce qui excite la couche de mélange et génère de nouvelles perturbations qui subissent le même processus. Ce mécanisme donne naissance à une dynamique oscillatoire, associée à l'apparition d'un bruit tonal de très forte intensité se propageant vers l'amont du jet, connu sous le nom de screech. Ce phénomène a été décris pour la première fois par [START_REF] Lighthill | On sound generated aerodynamically. i. general theory[END_REF][START_REF] Powell | On the mechanism of choked jet noise[END_REF].

Une formule analytique de prédiction de fréquence, basée sur le calcul du temps nécessaire pour une perturbation de parcourir une cellule de choc et pour une onde acoustique de revenir au niveau de la tuyère, a été introduite par [START_REF] Powell | On the mechanism of choked jet noise[END_REF]. Cette formule prédit assez bien la fréquence du bruit de screech. Dans cette thèse, il est montré que contrairement au cas du cylindre décris précédemment, une analyse globale de stabilité linéaire la prédit avec une bonne précision également. En effet, l'analyse révèle deux modes instables (mode A et B), dont la fréquence est très proche de ce qui est obtenu à partir de la formule de [START_REF] Powell | On the mechanism of choked jet noise[END_REF] (voir figure A.3). De plus, les écarts existants correspondent qualitativement à des observations expérimentales présentes dans des travaux précédents (surestimation de la fréquence par la formule de [START_REF] Powell | On the mechanism of choked jet noise[END_REF], voir par exemple Raman (1998)).

A.2.2 Interprétation et discussion des résultats

Les résultats obtenus sur le cas du screech, très différents de ceux d'un écoulement de cylindre (voir figure A.1) révèlent que malgré la grande amplitude des fluctuations, le mécanisme de sélection de fréquence (le bouclage acoustique) est peu influencé par des dynamiques non-linéaires. Il en résulte qu'une analyse linéaire prédit correctement la fréquence en jeu.

Ce mécanisme de bouclage acoustique est présent dans d'autres configurations, telles que les écoulements de cavités et les jets impactants. Des analyses de stabilité linéaire similaires semblent aboutir également, dans ces cas, à une bonne prédiction de fréquence, même loin de la criticalité (voir par exemple les travaux de Mettot et al. (2014b)). Il semblerait donc que le bouclage acoustique soit un mécanisme récurrent où les non-linéarités jouent un rôle marginal dans la sélection de fréquence d'un écoulement.

Bien que la généralité de ce résultat n'est ici pas démontrée, il est facile de l'interpréter d'un point de vue physique. La fréquence est déterminée principalement par la géométrie de l'écoulement (par exemple dans le cas de la cavité, par la longueur de la cavité), expliquant qu'un modèle linéaire aboutisse à une bonne prédiction de cette dernière. Outre cette interprétation physique, l'un des points centraux de la présente discussion est qu'elle prouve l'existence de situations où une analyse de stabilité linéaire est pertinente, même loin de la criticalité. Pour la prédiction de fréquence, l'amplitude des fluctuations n'est pas un critère suffisant pour trancher sur la validité d'une analyse de stabilité linéaire.

A.3 Analyse linéaire du champ moyen

Une analyse linéaire du champ de base présente néanmoins un certains nombre de limites, détaillées précédemment. Dans certains cas, ces limitations peuvent être contournées en étudiant la stabilité linéaire du champ moyen (voir Pier (2002), [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF]). Les conditions de validité de cette approche sont encore mal comprises, et ont fait l'objet de la seconde partie de la thèse.

A.3.1 Cadre théorique

Le présent manuscrit démontre qu'un opérateur linéaire nommé résolvante, qui fait directement intervenir la linéarisation des équations de Navier-Stokes autour du champ moyen, est au coeur de la dynamique non-linéaire d'un écoulement. En conséquence, il est possible d'extraire des informations sur la dynamique d'une configuration à partir de cet opérateur linéaire, qui peut être calculé par la seule connaissance du champ moyen. Cela fournit un cadre permettant de comprendre en quoi une analyse de stabilité linéaire du champ moyen peut caractériser un écoulement.

En particulier, une décomposition en valeur singulière (s.v.d.) de la résolvante fournit de nombreuses informations sur la dynamique de l'écoulement. Dans le cas particulier où la première valeur singulière de la résolvante présente une nette séparation avec les autres (hypothèse DSV, Dominant Singular Value), il est alors possible de montrer que les modes de Fourier sont directement proportionnels à la première réponse optimale (premier vecteur 1997)), dont le lien avec les analyses de résolvante est bien établi dans la littérature [START_REF] Sipp | Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer[END_REF].

Cela pose cependant la question de l'interprétation physique de l'hypothèse DSV. Une des grandes conclusions des travaux présents dans cette thèse lie cette dernière à l'existence d'un mécanisme d'instabilité convectif dominant (tel que le mécanisme de Kelvin-Helmholtz dans une couche de mélange). Cette situation est typiquement rencontrée dans un jet simple flux, ou dans un écoulement de marche descendante, étudié dans le paragraphe suivant.

A.3.2 Prédiction des modes de Fourier: validation sur un écoulement de marche

Ces résultats ont été vérifiés sur des données de simulation d'une marche descendante turbulente (Re = 57000), à partir de laquelle un champ moyen a été calculé (figure A.4). Il a été ensuite possible d'effectuer une analyse de résolvante pour une large gamme de fréquence. Comme attendu, la propriété DSV est vérifiée sur l'ensemble de cette gamme (figure A.5).

Cependant, malgré cette large séparation, les modes de Fourier, évalués par des techniques d'estimation spectrale de type méthode de Welch, ne sont proportionnels aux réponses optimales que dans certaines zones de l'écoulement. La figure A.6 montrent que même si l'accord est globalement bon, proche de la paroi, la différence entre réponse optimale et mode de Fourier normalisé est significative. Une approche différente pour le calcul de ces modes de Fourier, basée sur un filtrage POD (Proper Orthogonal Decomposition) des

A.3.3 Prédiction des spectres locaux

La partie précédente montre qu'il est possible de prédire, à une constante de proportionnalité près, la structure des modes de Fourier (la partie cohérente uniquement, obtenue par filtrage POD). En utilisant le spectre en quelques points de l'écoulement, il est possible de calculer pour chaque fréquence cette constante de proportionnalité, ce qui donne accès aux spectres locaux associés aux fluctuations.

Pour vérifier cette approche, ces constantes de proportionnalité (qui définissent une fonction d'amplitude Λ) ont été calculées à partir de deux points x 0 = (4, 1.5) et x 1 = (7, 0.1), ce qui a permis ensuite une prédiction en tout point. La figure A.7 présente les résultats de la prédiction des spectres associés à la vitesse longitudinale en cinq points: x a = (3, 1.5), x b = (5, 3), x c = (6, 1.5), x d = (8, 0.1) et x e = (9, 0.1). On y voit que la forme générale des spectres, la fréquence dominante et le niveau associé sont remarquablement bien prédits.

A.4 Conclusion

Une analyse linéaire de stabilité du champ de base permet de prédire et comprendre le caractère stationnaire ou instationnaire d'un écoulement. Cette thèse a cependant montré que contrairement à l'idée reçu, une telle analyse peut aussi caractériser fréquentiellement le comportement instationnaire de certains écoulements loin de la criticalité, malgré son caractère linéaire. En effet, un certain nombre de mécanismes de sélection de fréquence, souvent associés à un bouclage acoustique, sont très peu influencés par des dynamiques non-linéaire, indépendamment de l'amplitude des fluctuations associées. Pour ces cas, il est possible de prédire, avec une bonne précision, le fréquence des oscillations de l'écoulement.

Pour traiter les autres situations, un certains nombre de travaux antérieurs ont montré qu'une étude de la stabilité linéaire du champ moyen aboutissait à une caractérisation précise de certains écoulements. La présente thèse fournit un cadre mathématique qui permet de comprendre la pertinence de cette approche ainsi que sa généralité. Ce cadre met en évidence le rôle que joue la résolvante autour du champ moyen dans la dynamique d'un écoulement. Il est montré que lorsque cet opérateur présente une forte séparation de valeurs singulières, ce qui correspond à l'existence d'un mécanisme d'instabilité fort, alors les modes de Fourier de l'écoulement sont proportionnels aux modes de résolvante dominants. Ce résultat fournit des conditions mathématiques et physiques pour l'utilisation et le sens de diverses méthodes d'analyse du champ moyen, telles qu'une analyse PSE. De plus, cela permet de mettre en place un modèle de prédiction du spectre fréquentiel en tout point d'un écoulement, à partir d'une ou de quelques mesures ponctuelles et du champ moyen. L'ensemble de ces résultats a été illustré et validé sur un cas de marche descendante turbulente. Enfin, cela a été exploité dans un cadre expérimental, afin de reconstruire le comportement instationnaire d'un jet rond transitionnel, à partir de la seule connaissance du champ moyen et d'une mesure ponctuelle.

L'étude a révélée que, sous certaines précautions expérimentales, la reconstruction est très précise et robuste. Cette problématique est traitée dans cette thèse en s'intéressant dans un premier temps aux approches classiques de stabilité, qui étudient l'évolution de petites perturbations autour d'une solution stationnaire -un champ de base -des équations de Navier-Stokes. Pour cela, le phénomène du screech -un bruit tonal que peuvent causer les jets sous-détendus -est étudié d'un point de vue de la stabilité linéaire. Les résultats obtenus montrent que la dynamique non-linéaire du phénomène est correctement prédite par une analyse linéaire de stabilité du champ de base. Une confrontation avec d'autres analyses similaires montre qu'un tel résultat n'est pas toujours observé. Cependant, lorsque les oscillations auto-entretenues d'un écoulement sont provoquées par un bouclage acoustique, comme c'est le cas entre autres pour le screech, l'écoulement de cavité ou encore les jets impactants, alors les non-linéarités ont une faible influence sur le phénomène de sélection de fréquence. Cela explique la capacité d'une analyse linéaire à caractériser ces écoulements, même dans le régime non-linéaire. Abstract: Linear stability theory has been intensively used over the past decades for the characterization of unsteady flow behaviors. While the existing approaches are numerous, none has the ability to address any general flow. Moreover, clear validity conditions for these techniques are often missing, and this raises the question of their general reliability.

In this thesis, this question is addressed by first considering the classical stability approach, which focuses on the evolution of small disturbances about a steady solutiona base flow -of the Navier-Stokes equations. To this end, the screech phenomenon -a tonal noise that is sometimes generated by underexpanded jets -is studied from a linear stability point of view. The results reveal that the nonlinear dynamics of this phenomenon is well-predicted by a linear base flow stability analysis. A confrontation with other similar analyses from the literature shows that such a satisfactory result is not always observed. However, when a self-sustained oscillating flow is driven by an acoustic feedback loop, as it is the case for the screech phenomenon, cavity flows and impinging jets for instance, then the nonlinearities have a weak impact on the frequency selection process, explaining the ability of a linear analysis to characterize the flow, even in the nonlinear regime.

Another alternative approach, based on a linearization about the mean flow, is known to be successful in some cases where a base flow analysis fails. This observation from the literature is explained in this thesis by outlining the role of the resolvent operator, arising from a linearization about the mean flow, in the dynamics of a flow. The main finding is that if this operator displays a clear separation of singular values, which relates to the existence of one strong convective instability mechanism, then the Fourier modes are proportional to the first resolvent modes. This result provides mathematical and physical conditions for the use and meaning of several mean flow stability techniques, such as a Parabolised Stability Equations analysis of a mean flow. Moreover, it leads to a predictive model for the frequency spectrum of a flow field at any arbitrary location, from the sole knowledge of the mean flow and the frequency spectrum at one or more points. All these findings are illustrated and validated in the case of a turbulent backward facing step flow. Finally, these results are exploited in an experimental context, for the reconstruction of the unsteady behavior of a transitional round jet, from the sole knowledge of the mean flow and one point-wise measurement. The study shows that, after following a few experimental precautions, detailed in the manuscript, the reconstruction is very accurate and robust.

Université Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

Fig. 1

 1 Fig. 1.1 Satellite observation of a Von-Karman vortex street behind the Juan Fernandez Island, by B. Cahalan, NASA GSFC. This type of fluid behavior may be explained by a stability analysis.

  The codes used to produce the results have been developed by S. Beneddine (SB), based on 1.5 Aim and structure of the thesis an existing code from C. Mettot. All the computations have been performed by SB. The paper has been written by SB with feedback from D. Sipp (DS). Paper 2 ( §5): "Conditions for validity of mean flow stability analysis", S. Beneddine, D. Sipp, A. Arnault, J. Dandois and L. Lesshafft, Journal of Fluid Mechanics 798, pp. 485-504 (2015). The simulation has been performed by A. Arnault under the supervision of J. Dandois. The rest of the codes and computations have been done by SB, except for the PSE code that is a Python porting (done by SB) of a Matlab code provided by L. Lesshafft (LL). The mathematical developments and models have been done by SB with feedback from DS. The paper has been written by SB with feedback from DS and LL. Paper 3 ( §6): "Time-resolved reconstruction of a round jet from point-wise measurements and mean flow stability analysis", S. Beneddine, R. Yegavian, D. Sipp and B. Leclaire, submitted to Journal of Fluid Mechanics.

  Fig. 2.3 Impinging jet configuration: geometry of the domain and boundary conditions.

Fig. 2 . 4

 24 Fig. 2.4 Evolution of the residual norm with respect to the number of quasi-Newton iterations, for the impinging test case.

Fig. 2 . 5

 25 Fig. 2.5 Density field of the impinging jet base flow, computed with the quasi-Newton method.

Fig. 2

 2 Fig. 2.6 Spectrum of the impinging jet (Re = 660). The red dots represent the two unstable eigenvalues that has been found (St = 0.48 and St = 0.55).

Fig. 2 . 7

 27 Fig. 2.7 Density field of mode A (St = 0.48) for Re = 660.

  Fig. 3.1 Scheme of the physical configuration of the flow (nondimensional quantities). The two variable parameters are the JPR and the lip thickness h l .

  Fig. 3.2 Computational domain and boundary conditions. Sponge regions are in grey and the physical domain of interest is in white.

Fig. 3 . 5

 35 Fig. 3.5 Real part of the density field of the marginally unstable mode (JPR=1.12, h l = 0.063). One can notice the upstream propagating waves, typical from screech. The structures in the red rectangle propagate downstream. The black arrows show the direction of propagation of the acoustic waves.

Fig. 3

 3 Fig.3.6 Spectrum of the JPR=1.12 configuration (circles), with the evolution of mode A (red) and B (blue) when we increase the JPR, h l = 0.063. Darker colors correspond to higher JPRs.

Fig. 3

 3 Fig. 3.8 Left: Comparison between the analytically derived evolution of the frequency with respect to the JPR (based on equation (3.6), using equations (3.4) (continuous line) and (3.3)), and the frequencies of mode A (red) and B (blue). Right: evolution of the growth rate σ of mode A (red) and mode B (blue) with respect to the JPR

Fig. 3 .

 3 Fig. 3.10 Real part of the vorticity field of the two unstable modes (JPR=1.18, h l = 0.063). Top: mode A, bottom: mode B.

Fig. 3 .

 3 Fig. 3.11 Spatial distribution of the quantity v x ϕ /c a , with v x ϕ the local phase velocity of the vorticity field, and c a = 1.04 the ambient sound velocity, for the case JPR=1.18, h l = 0.063. Only regions where vorticity is high enough with respect to numerical errors are shown.

Fig. 4

 4 Fig. 4.2 From Mettot et al. (2014a): cavity flow at a Reynolds number Re = 860000. Comparison between the unstable (σ > 0) eigenvalues obtained from a RANS base flow analysis (left) and the pressure spectrum at a given point in the shear layer, computed by time marching the RANS equations (right). Both results are based on the k -ω model of Wilcox (Roe scheme). The dashed lines indicates the Rossiter peak and harmonics observed in the experiment by[START_REF] Forestier | The mixing layer over a deep cavity at high-subsonic speed[END_REF]. The eigenvalues matching these frequencies are indicated by the red dots.

Fig. 4

 4 Fig. 4.3 Schematic procedure to identify the situations where a given kind of linear base flow stability analysis is relevant. Note that the results of a RANS base flow stability analysis are however conditioned by the ability of the RANS equations to model a given behavior.

Fig. 5 . 1

 51 Fig. 5.1 Streamwise velocity of the mean flow, computed by averaging in time and in the spanwise direction the results of the 3D unsteady simulation. Contours are in equal increments from -0.2 to 1. Negative contours are dashed, the thickened contours correspond to a zero velocity.

  Fig. 5.2 Spectra of the streamwise (spanwise-averaged) velocity computed from the unsteady ZDES simulation at two locations: (a) x = 2, y = 1.5, (b) x = 7.0, y = 0.1. St θ is the Strouhal number based on the free-stream velocity and the momentum thickness at the step location, and F + the Strouhal number based on the free-stream velocity and the reattachment length.

  Figure 5.4 compares the amplitude of the streamwise component of û computed with method (a) (figure 5.4(c)) and method (b) (figure 5.4(a)), with the optimal response ψ 1 (figure 5.4(b)), for ω = 2.1 (the frequency associated with the Kelvin-Helmholtz mechanism, see the peak in figure 5.2(a)).

Fig. 5 . 4

 54 Fig. 5.4 Comparison between the amplitude of the normalized streamwise velocity field of (a) the spectral mode computed by POD-filtering, (b) the optimal response, (c) the spectral mode computed by r.m.s., for ω = 2.1 and two-dimensional perturbations β = 0. The three vertical dashed lines represent the locations where profiles were extracted. Figures(d), (e), (f) compare the pod-filtered profile (thick continuous line), the optimal response profile (red dashed line), and the spectral r.m.s. profile (thin continuous line) respectively for x = 2, x = 3.5 and x = 5.

Fig. 5

 5 Fig. 5.7 Comparison between the real part of the streamwise velocity of (a) the dominant optimal response and (b) the PSE field, for ω = 2.1 (normalized in phase and amplitude). Each figure displays ten equally spaced contours ranging from -1 to 1.

Fig. 5 . 9

 59 Fig. 5.9 Comparison of the estimate of the streamwise velocity spectra, from model 1 (dotted line) and model 2 (red dashed line), with the simulation (continuous line), at five points: (a): x 1 = (3, 1.5), (b): x 2 = (5, 3), (c): x 3 = (6, 1.5), (d): x 4 = (8, 0.1), (e): x 5 = (9, 0.1).

  Fig. 6.1 Graphical illustration of the reconstruction procedure.

Fig. 6

 6 Fig.6.4 Mean streamwise velocity field, computed from the reconstruction set of snapshots. Point A (x 0 = 2.5, y 0 = 0.4) corresponds to the location where the unsteady signal used for the reconstruction has been extracted (also extracted from the reconstruction set). The nozzle position (schematically represented on the left of the figure) has been used to set the origin of the reference frame.

Fig. 6 . 5

 65 Fig. 6.5 Streamwise velocity at point A (x 0 = 2.5, y 0 = 0.4) vs. time, extracted from the reconstruction set. The reconstruction presented in section 6.4.2 has been done solely based on this signal and the mean flow shown in Figure 6.4.

Fig. 6

 6 Fig.6.7 Frequency spectrum of the streamwise velocity at (a) x = 1.0, y = 0.3 and (b) x = 1.5, y = 0.25, computed from the PIV results (reference set). The upstream spectrum (a) displays a clear peak at St = 0.76, linked to the Kelvin-Helmholtz mechanism, while further downstream, (b) displays a broader spectrum center around St = 0.38. This latter frequency is related to the typical downstream vortex pairing that occurs in such flows.

Fig. 6 . 8

 68 Fig. 6.8 Comparison of the normalized modulus of the axial velocity of (a): the Fourier mode computed by an FFT of the TR-PIV measurements and (b): the dominant resolvent mode computed from the experimental mean flow only (St = 0.76). Figures(c), (d) compare profiles from the Fourier mode (dashed line) and the dominant resolvent mode (continuous red line), extracted at r = 0.4 and x = 2.0, respectively.

Fig. 6

 6 Fig. 6.9 Comparison of the normalized phase of the axial velocity of (a): the Fourier mode computed by an FFT of the TR-PIV measurements and (b): the dominant resolvent mode computed from the experimental mean flow only (St = 0.76). Figures (c), (d) compare profiles from the Fourier mode (dashed line) and the dominant resolvent mode (continuous red line), extracted at r = 0.4 and x = 2, respectively. The phase profiles have been unwrapped.

Fig. 6 .

 6 Fig. 6.10 Comparison of the normalized modulus of the radial velocity of (a): the Fourier mode computed by an FFT of the TR-PIV measurements and (b): the dominant resolvent mode computed from the experimental mean flow (St = 0.38). Figures(c), (d) compare profiles from the Fourier mode (dashed line) and the dominant resolvent mode (continuous red line), extracted at r = 0.5 and x = 2.5, respectively.

Fig. 6 .

 6 Fig. 6.11 Comparison of the normalized phase of the radial velocity of (a): the Fourier mode computed by an FFT of the TR-PIV measurements and (b): the dominant resolvent mode computed from the experimental mean flow (St = 0.38). Figures (c), (d) compare profiles from the Fourier mode (dashed line) and the dominant resolvent mode (continuous red line), extracted at r = 0.5 and x = 2.5, respectively. The phase profiles have been unwrapped.

Fig. 6 .

 6 Fig. 6.13 Comparison of the axial velocity of (a): the reconstructed field and (b): the TR-PIV field (reference set) at t = 75.

Fig. 6 .

 6 Fig. 6.18 (a): Reconstructed pressure fluctuation field at t = 75, (b): mean square pressure fluctuation field p ′2 .

Fig. 6 .

 6 Fig.6.20 Definition of the eight points considered to study the influence of the input point location on the accuracy of the reconstruction. A is the original point considered in section 6.4.2.

Fig. 6 .

 6 Fig. 6.21 Comparison of the streamwise velocity of (a): the unbiased reconstructed field (no misplacement of the input sensor), (b): biased reconstructed field with δx = δr = 0.05, (c): biased reconstructed field with δx = δr = 0.1, (d): biased reconstructed field with δx = δr = -0.2, for t = 10 and (x 0 , r 0 ) = (2.5, 0.4). The reconstruction is performed using the axial velocity component.

Fig. 6 .

 6 Fig. 6.22 (a): sensitivity map with respect to misplacements of the input sensors, for the radial velocity u r used an an input (the quantity displayed is |b x | + |b r |), (b): radial velocity modulus of the normalized dominant resolvent mode. For both figures, St = 0.76. The figures show a relative correspondence between low-sensitivity and high-energy regions.

Fig. 6 .

 6 Fig. 6.24 Comparison between the normalized modulus of the dominant resolvent mode (radial velocity) at St = 0.76 computed from (a): mean flow A (reference mean flow) and (b): mean flow B (larger interrogation window in the PIV processing).

Fig. 6 .

 6 Fig. 6.25 Comparison between a reconstructed radial velocity snapshot computed from (a): mean flow A (reference mean flow) and (b): mean flow B (larger interrogation window in the PIV processing).

Fig. 7

 7 Fig. 7.1 Backward facing step treated in chapter 5: comparison of the real part of the streamwise velocity of (a) the Fourier mode obtained by a straightforward FFT based on the full ZDES data of chapter 5, (b) the resolvent mode computed from the mean flow, for ω = 0.7 (most energetic frequency of the flow). The figures display ten equally spaced contours from the minimal to the maximal value of each field.

Fig. 7 . 2

 72 Fig. 7.2 Backward facing step treated in chapter 5: comparison of the real part of the streamwise velocity of (a) the Fourier mode obtained by a straightforward FFT based on the full ZDES data of chapter 5 (simulated time: 310 non-dimensional time units), (b) a similar FFT mode obtained from a truncated signal corresponding to a simulated time of 280, (c) the resolvent mode computed from the mean flow, for ω = 11. The figures display ten equally spaced contours from the minimal to the maximal value of each field.

Fig

  Fig.A.2 Gauche: Densité du champ de base d'un jet sous-détendu calculée pour JPR=1.12, faisant apparaître les cellules de chocs. Droite: zoom sur la sortie de la tuyère, présentant les principaux éléments d'une cellule de choc. Les figures présentent 15 contours également répartis entre 0.9 et 1.4 (unités adimensionnées par la densité du milieu ambiant). Le nombre de Reynolds du jet est de 3030.

Fig. A. 3

 3 Fig. A.3 Comparaison de la fréquence obtenue à partir d'une analyse globale de stabilité linéaire du champ de base (rouge: mode A, bleu: mode B) et de la prédiction donnée par la formule de Powell (1953), pour différents JPRs.

Fig. A. 4

 4 Fig. A.4 Champ moyen (vitesse longitudinale) calculé en moyennant temporellement et suivant la direction transverse une simulation 3D. Les contours sont également espacés entre -0.2 et 1 (unités adimensionnées par la vitesse à l'infini). les contours négatifs sont en pointillés, le contour épais correspond à une vitesse nulle.

Fig. A. 5

 5 Fig. A.5 Comparaison des trois premières valeurs singulières de la résolvante pour des perturbations 2D (β = 0): la première valeur singulière µ 1 est beaucoup plus grande que les autres, i.e. l'hypothèse DSV est vérifiée pour la gamme de fréquence présentée ici.

Fig. A. 6

 6 Fig. A.6 Comparaison du module des champs de vitesse longitudinal normalisés de (a) le mode de Fourier calculé par filtrage POD, (b) la réponse optimale, (c) le mode de Fourier calculé par méthode de Welch, pour ω = 2.1 (pulsation temporelle adimensionée) et perturbations 2D (β = 0). les trois lignes verticales pointillées représentent les positions des profils qui ont été extraits. Les figures (d), (e), (f) comparent le profil du mode avec filtrage POD (ligne continue épaisse), le profil de la réponse optimale (ligne rouge pointillée), le profil du mode obtenu par la méthode de Welch (ligne continue fine), pour x = 2, x = 3.5 et x = 5, respectivement.

Fig. A. 9

 9 Fig. A.9 Comparaison de la vitesse axiale instantanée de (a): le champ reconstruit (b): le champ TR-PIV (champ de référence).

Titre:

  Caractérisation de comportement d'écoulement instationnaire par analyse de stabilité linéaire. Mots clés: Champ de base, Champ moyen, Stabilité, Résolvante, Structures Cohérentes, Reconstruction d'écoulement.Résumé: Au cours des dernières décennies, la théorie de la stabilité a été intensivement utilisée pour caractériser le comportement instationnaire d'écoulements. Cela a donné naissance à un grand nombre d'approches, mais malheureusement chacune d'entre elles semble présenter ses propres limitations. De plus, leurs conditions de validité sont encore très mal connues, ce qui soulève la question de la fiabilité de ce genre de méthodes dans un cas général.

  Une autre approche, consistant à étudier la stabilité linéaire du champ moyen, a montré de bons résultats dans certaines configurations qui ne peuvent être correctement étudiées par une analyse linéaire du champ de base. Cela est justifié dans cette thèse en mettant en évidence le rôle que joue la résolvante autour du champ moyen dans la dynamique d'un écoulement. Il est montré que lorsque cet opérateur présente une forte séparation de valeurs singulières, ce qui correspond à l'existence d'un mécanisme d'instabilité fort, alors les modes de Fourier de l'écoulement sont proportionnels aux modes de résolvante dominants. Ce résultat fournit des conditions mathématiques et physiques pour l'utilisation et le sens de diverses méthodes d'analyse du champ moyen, telles qu'une analyse d'équations de stabilité parabolisées (Parabolised Stability Equations). De plus, cela permet de mettre en place un modèle de prédiction du spectre fréquentiel en tout point d'un écoulement, à partir d'une ou de quelques mesures ponctuelles et du champ moyen. L'ensemble de ces résultats est illustré et validé sur un cas de marche descendante turbulente. Enfin, cela est exploité dans un cadre expérimental, afin de reconstruire le comportement instationnaire d'un jet rond transitionnel, à partir de la seule connaissance du champ moyen et d'une mesure ponctuelle. L'étude montre que, sous certaines précautions expérimentales, la reconstruction est très précise et robuste. Title: Characterization of unsteady flow behavior by linear stability analysis. Keywords: Base flow, Mean flow, Flow Stability, Resolvent, Coherent Structures, Flow Reconstruction.

  These different situations are summarized in figure 4.3.Another point may motivate the use of a linear base flow analysis: when a base flow verifies the RPIF property, then the resulting unstable modes indicate what phenomenon may be generated through linear mechanisms. For instance in the case of the cavity, this gives insights about the harmonics that can be seen in figure4.2. Such harmonics are generally considered to stem from nonlinear interactions. But as one can see in figure4.2 in the cavity case, some unstable eigenvalues match harmonic frequencies. This reveals that the harmonics may also be generated through purely linear mechanisms.

	Unsteady flow behavior
	Feature of interest ?	Dynamical feature (e.g. freq.)
	Onset of unsteadiness (e.g. Rec)
		Linear (eg. acoustic feedback)	Freq. selection process ?
	High	Re ?	Low	Nonlinear
	Perform a RANS linear	Perform a classical linear	Base flow approach
	base flow stability analysis	base flow stability analysis	is not adapted

  Comparison of the global error for different input point/set of points.

		A is the original point considered in section
	6.4.2.
	Point(s) error E	x	r
	A	0.038 2.5 0.4
	B	0.047	2 0.4
	C	0.062 1.5 0.4
	D	0.092	1 0.4
	E	0.059 2.5 0.6
	F	0.061	2 0.6
	G	0.069 1.5 0.6
	H	0.117	1 0.6
	A and H	0.037	-	-
	E to H	0.046	-	-
	All points	0.035	-	-
	Table 6.1

2 Analyse linéaire du champ de base

  Fig. A.1[START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF]: la ligne épaisse représente l'évolution du nombre de Strouhal (fréquence adimensionnée) du sillage du cylindre, en fonction de nombre de Reynolds. La courbe pointillée représente la prédiction provenant d'une analyse de stabilité globale linéaire du champ de base. Les triangles correspondent à la même analyse effectuée autour du champ moyen.La littérature scientifique présente des cas pour lesquels une analyse linéaire de stabilité du champ de base ne permet pas une prédiction correcte de la fréquence d'oscillation d'un écoulement. Un exemple classique de cela est l'écoulement derrière un cylindre, qui présente une bifurcation de Hopf super-critique pour un nombre de Reynolds critique Re c = 46. Dès lors que l'on s'éloigne de cette valeur pour aller vers des nombres plus grand, l'analyse linéaire sous-estime significativement la fréquence (voir figure A.1).Cette situation n'est cependant pas systématique, et la littérature scientifique contient un certain nombre de contre-exemples. Dans cette thèse, il est proposé l'étude d'un phénomène connu sous le nom de screech, pour lequel une analyse linéaire du champ de base prédit avec une bonne précision la fréquence de l'écoulement, alors même que les perturbations ont une grande amplitude. Cette étude, présentée dans le paragraphe suivant, a servi de base pour une discussion plus générale sur la pertinence physique d'une telle analyse.

solution stationnaire. Dans une telle situation,

[START_REF] Pier | On the frequency selection of finite-amplitude vortex shedding in the cylinder wake[END_REF]

puis

[START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF]

ont proposé de mener des analyses de stabilités linéaires, reposant sur des approches classiques, mais utilisant le champ moyen (moyenne temporelle du champ instationnaire) au lieu du champ de base. Leurs résultats concernant l'étude d'un écoulement de cylindre révèlent que cette approche fournit une excellente prédiction de la fréquence d'oscillation du sillage, même loin de la criticalité (voire figure A.1). Ce résultat très intéressant soulève cependant de nombreuses questions: outre la fréquence, quelles quantités peuvent-être ainsi caractérisées ? Quelle est la généralité de l'approche ? Quelle sont ses conditions de validité ? La réponse à ces interrogations est encore mal connue, et constitue la problématique traitée dans la seconde partie de la thèse.

A.

for a given flow behavior, the choice of this norm depends on the physical quantity of interest (kinetic energy, acoustic energy, etc.). This is briefly discussed in chapter

Procedures and numerical strategies

Remerciements

Appendix A -Résumé en français

A.1 Problématique scientifique

La très grande majorité des applications en mécanique des fluides fait intervenir des comportements instationnaires qu'il est nécessaire de caractériser (aéronautique, ingénierie, météorologie, etc.). Au cours des dernières décennies, la théorie de la stabilité linéaire a été intensivement utilisée pour cela. Le concept central de cette dernière consiste en l'étude de l'évolution de petites perturbations autour d'un état stationnaire donné (un champ de base). Il est ainsi possible de montrer que certains comportements instationnaires sont le résultat du caractère instable de ces perturbations, qui peuvent être amplifiées à des niveaux importants par des mécanismes d'instabilité.

Cela correspond historiquement à l'approche classique pour analyser la stabilité d'un écoulement, qui présente plusieurs limites. Tout d'abord, même si elle permet d'étudier les bifurcations d'un état stationnaire vers un état instationnaire, et donc de prédire l'apparition d'instationnarités, une hypothèse centrale de ce type d'étude consiste à considérer des perturbations de petite amplitude, ce qui permet une linéarisation des équations. Par conséquent, dans le cas d'un écoulement présentant une bifurcation, il n'y a aucune garantie sur la pertinence physique de l'approche loin de la criticalité, où les fluctuations ont une grande amplitude. Cette question est traitée dans la première partie de la thèse, dédiée aux analyses linéaire de stabilité du champ de base.

Une autre limite vis-à-vis de l'étude du champ de base est qu'il n'y a aucune garantie théorique sur l'existence d'un tel champ. En effet, un écoulement est régi par les équations de Navier-Stokes, et dans un cas général, il est possible que ces équations n'aient pas de