Beverley Mckeon

Etienne Mémin

Lutz Lesshafft

Paola Cinnella

Carlo Cossu

Pedro Jean-Lou Carlos

Damien Holly

Ye-Bonne Nicolas

Arnold

Anthony

Nicolo, Colin, Tristan, Samir, Robin, Mathieu Simon Benjamin

for the mixed periodic/broadband fluctuations on turbulent flows

Chapter 5. Mean-flow analysis

v

Résumé long

La modélisation des écoulements turbulents aérodynamiques à l'aide des simulation numérique est une tâche difficile, principalement en raison des questions liées à la modélisation de la turbulence. C'est ce qu'on appelle le problème de fermeture, visant à modéliser l'influence de la partie non résolue du spectre dans la partie résolue. Ces modèles ont commencé à être proposés il y a plusieurs décennies et c'est encore un sujet de recherche actif à ce jour. Ils diffèrent fondamentalement en fonction de la portion du spectre modélisée. Par exemple, les simulations numériques directes (DNS) tentent de résoudre le spectre entier, là où rien n'est modélisé. C'est bien sûr la méthode la plus coûteuse et, dans de nombreuses applications industrielles, elle n'est pas réalisable avec les ordinateurs actuels. Les simulations par grands tourbillons (LES) visent à utiliser le maillage sur lequel le problème est discrétisé comme un filtre, de sorte que les fluctuations non résolues sont celles qui ne peuvent pas être résolues sur ce maillage. Ce type de modélisation peut encore être très coûteux. En effet, si des couches limites sont présentes, le maillage doit y être bien affiné pour que la simulation soit précise. Un autre type de modélisation est le modèle RANS (Reynolds-Averaged Navier-Stokes), où seul l'écoulement moyen dans le temps est résolu et la totalité des fluctuations est modélisée. C'est le modèle le plus simple à résoudre et aussi le plus utilisé à l'industrie. Cependant, il produit souvent des solutions imprécises, surtout en présence de décollemet. Il pourrait donc être avantageux d'incorporer des données de plus haute fidélité afin de corriger ces défauts de la solution et, en conséquence, d'augmenter sa valeur prédictive. C'est l'objectif de l'assimilation de données. Dans ce travail, nous traiterons de l'assimilation de données variationnelles où on un paramètre de réglage lié au modèle est ajusté de telle sorte que la solution du modèle soit aussi proche que possible des données de haute fidélité. Dans la première partie de ce manuscrit, nous aborderons ce problème en utilisant un modèle RANS particulier, le modèle Spalart-Allmaras. Ceci veut dire que nous evisageons de reconstruire le champ moyen seulement. Deux types de paramètres d'accord seront utilisés, à savoir une force volumique et un paramètre liè directement au modèle de Spalart-Allmaras. Nous étudierons leurs performances en fonction du type de donnés à être assimilés, allant des champs de vitesse denses (provenant éventuellement d'une Particule Image Velocimetry bien résolue), aux mesures ponctuelles telles que la pression de paroi et la vitesse, obtenues par des tubes de Pitot. La configuration étudiée ici est une marche descendente à Re = 28275.

Dans la deuxième partie du manuscrit, nous étudierons le rôle de la reconstruction des fluctuations à l'aide d'une analyse du Resolvent des équations linéarisés sur le champ moyen. Il s'agit d'une technique linéaire où la forme des modes, à une fréquence donnée, peut être prédite, sous certaines hypothèses assez générales. Cependant, comme cette technique est linéaire, l'énergie de ces modes (ou leur amplitude, en d'autres termes) ne peut pas être déterminée. Pour cette raison, nous utiliserons des données externes pour les régler. Ces données, cette fois-ci, doivent être résolus en temps, ce qui diffère de la première partie où seules des données moyennes (moyennes) ont été utilisées, mais elles peuvent être extrêmement peu nombreuses, généralement quelques sondes ponctuelles. Le principal inconvénient de cette approche est la nécessité de connaître a priori le champ moyen. Pour cette raison, une possibilité serait d'appliquer les techniques présentées vi dans la première partie du manuscrit pour obtenir ce champ moyen et d'appliquer ensuite l'analyse de Resolvent pour récupérer la fluctuation. Cependant, comme nous disposons d'un modèle pour la fluctuation, nous pourrions incorporer ses effets dans la modélisation du champ moyen, ce qui éviterait d'avoir besoin des données du flux moyen et permettrait de reconstruire le flux moyen et la fluctuation à partir de quelques données ponctuelles résolues dans le temps. Cette technique sera appliquée à deux flux de transition (sur un cylindre de section carrée, Re = 100, et une marche descendente, Re = 500) et un flux turbulent (sur un cylindre carré, Re = 22000), où un modèle de turbulence (de nouveau Spalart-Allmaras) est résolu avec le flux moyen et sa fluctuation par l'analyse de Resolvent.

Dans une troisième partie, nous nous concentrerons sur l'analyse du résolvant lorsqu'elle est appliquée au flux autour du cylindre carré à Re = 22000. Pour cette configuration, nous observons deux phénomènes distincts. Le premier est le détachement de vortex périodique, typique des écoulements sur les corps peu profilés. L'autre phénomène est le mode de Kelvin-Helmholtz, produit au niveau des couches de cisaillement qui se déplacent selon la fréquence de détachement de vortex. Pour pouvoir saisir cette caractéristique, nous proposons un cadre pour l'opérateur résolvant dans lequel nous examinons la fluctuation à partir des équations linéarisées autour du détachement ve cortex. Cette linéarisation est naturelle si nous considérons la triple décomposition où le champ moyen, ainsi que les modes de détachement de vortex sont séparés des phénomènes de fréquence à large bande, représentés par le Kelvin-Helmholtz. Cela conduit à un opérateur de résolution de type Floquet, qui est calculé ici de manière quasi-stationnaire.

Chapter 1 Introduction 1 Context and Motivations

Experimental and numerical fluid dynamics have been used for several years as tools to investigate and understand the physics of fluid flow phenomena. Separately, they provide both useful and different insights on those phenomena. From the point of view of the numerical scientist or engineer, we usually have complete control on the numerical simulation we run in the sense that we can design the geometry of the physical settings, apply different boundary or initial conditions and investigate several regimes of non-dimensional numbers, such as the Reynolds and Mach numbers, etc. However, all this freedom has a drawback since, in order to produce a solution that represents a real physical experiment, we need to precisely know those actual quantities, such as boundary and initial conditions. Furthermore, apart from those aspects, the numerical model itself is often seen as a major source of discrepancy between the resulting numerical solution and what happens in real flow phenomena, especially in turbulent flow simulations, as it will be discussed later. On the other hand, experimental scientists and engineers have developed several tools to probe fluid flows, generating data from real experiments that are often regarded as more valuable than numerical ones. However, this experimental approach suffers from several drawbacks as well. For example, the data extracted this way are often insufficient for describing, in a complete manner, the fluid flow. For example, having the wall-pressure only does not provide direct and precise information of what goes on the bulk of the flow. Another related issue is the fact that the extracted data is often contaminated by external noise and errors intrinsic to the measurement, leading to an uncertainty on those quantities. For those reasons, those two techniques are regarded as complementary and have been extensively used in parallel as cross-validation of one another. For example, having a numerical solution that matches well with the experimental data gives us confidence to extrapolate this solution to regions where the probes could not reach, enlarging this way our understanding of the physical phenomenon. Although this approach provides several insights, it becomes less useful when the numerical solution does not match well the experimental data. In those cases, the numerical simulation (including the model, boundary and initial conditions, etc), becomes questionable and may lose its predictive power. It is interesting then to modify the model such that its solution matches the experimental data, improving this way its quality and, hopefully, its predictability in regions far from the probe locations. For this purpose, the data-assimilation procedure has been designed as a mathematically well-defined manner to include the experimental data in the model, correcting it. In the literature, several approaches in several different contexts have been developed in order to solve this problem. In what follows, we provide a brief description of the most common approaches, identified in two (maybe overlapping) groups, Kalman Filters and Variation data-assimilation, followed by some of their applications in the context of Fluid Mechanics. We will then precise what is the goal of the present thesis, for which we present a brief state-of-the-art, followed by an overview of the thesis.

Data-Assimilation Techniques

Kalman Filters (KF)

The data-assimilation concept goes back to the theory of modern control of linear dynamical systems (see [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]), where an estimation of the whole state of the system from limited partial measurements is needed in order to act on the system with a control law, stabilizing it. One can show that, for deterministic linear systems (under the condition of observability), one can always design a (not unique) dynamical system, called observer, mimicking the dynamics of the original one, such that its output tends towards the true state as time progresses. However, if noise is present, the performance of such observer can be drastically deteriorated and the error committed in the prediction of the state can be large. The Kalman filter is the design of such observer so that its prediction is as close as possible to the actual state. The final form of this filter can be seen as a prediction step, where the model of the dynamical system is applied on the a priori state, and a correction step, where the external data (or measure) is taken into account through the filter, which is constructed based on the covariance matrix of the state vector. We remark that, in this procedure, not only the state vector is propagated in time, but also its corresponding covariance matrix. For this reason, this approach is mainly used when the system has only a few degrees of freedom, as it is often the case in many engineering application, not involving Partial Differential Equations (PDE). However, in cases where the dimension of the state vector is prohibitively large, the propagation of the covariance matrix becomes unfeasible. For this purpose, the Kalman filter has been modified with an approximation of this matrix, consisting in the propagation of a cloud of points (or ensemble) in the phase space, instead of only one point. The evolution of this cloud provides, if the number of points is sufficiently large, an estimation of the statistics. In those circumstances, this filter is called Ensemble Kalman Filter (EnKF, see [START_REF] Evensen | Data-Assimilation: yhe ensemble Kalman filter[END_REF]). This concept was first introduced by [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte-Carlo methods to forecast error statistics[END_REF] in the context of oceanography.

In fluid mechanics, this technique has been applied for several purposes: the estimation of velocity field from Particle-Tracking Velocimetry snapshots in a planar-jet configuration at Re = 2000 [START_REF] Suzuki | Reduced-order Kalman-filtered hybrid simulation combining particle tracking velocimetry and direct numerical simulation[END_REF]), the recovery of low-Reynolds wall-bounded turbulent flows from wall measure [START_REF] Colburn | State estimation in wall-bounded flow systems. part 3. the ensemble Kalman filter[END_REF]), estimation of Mach number, angle of attack and eddy-viscosity field at high Reynolds number flows in complex flow configuration [START_REF] Kato | A data assimilation methodology for reconstructing turbulent flows around aircraft[END_REF].

Variational data-assimilation (VDA)

The variational data-assimilation technique is based on the theory of optimal control (see, for example [START_REF] Lions | Optimal control of systems Governed by partial differential equations[END_REF]) where we act on a system such that its output satisfies 3. Goal of the Present Work 3 as closely as possible some target task. In the data-assimilation context we dispose in general of some tuning parameter, that plays the role of the control variables, which is related to the model that we effectively solve (which can be initial/boundary conditions, the spatial distribution of an unknown forcing term, etc). This tuning parameter will be found such that the model's solution is as close as possible to the external reference data. One of the earliest applications of this technique was on the weather forecast problem [START_REF] Talagrand | Variational assimilation of meteorological observations with the adjoint vorticity equation[END_REF], [START_REF] Liu | An ensemble-based four-dimensional variational data-assimilation scheme. Part i: technical formulation and preliminary test[END_REF]) under the name of 3D/4D-Var where the tuning parameters were the initial condition of the dynamical model such that some cost, function of the error between the model's prediction and the external data, is minimal. In this context, we can see that the main difference between this approach and Kalman Filters is that, VDA can be posed as an optimal problem on the initial condition where EnKF can be view as a way to correct an a priori condition with the information of the external data. Indeed, in general, VDA problems can be viewed as optimization under constraint ones and, when the number of tuning parameters becomes large, gradient-descent methods need to be employed so that the optimization procedure reaches a local minimum efficiently. This gradient is often obtained with the aid of the Lagrangian formalism where the adjoint problem needs to be solved. This adjoint problem represents in general a supplementary numerical tool (other than the model's resolution) and may pose some difficulties.

Several examples can be presented in the context of fluid mechanics and aerodynamics, such as estimation of boundary and/or initial conditions from PIV data using DNS in low Reynolds numbers configurations (Gronskis A. [2013], [START_REF] Yegavian | Performance assessment of PIV super-resolution with adjoint-based data assimilation[END_REF]), estimation of fluid flow from a sequence of sparse snapshots [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF]), as a tool for determining pollutant sources in urban areas, where Large Eddy Simulations (LES), coupled with Proper Orthogonal Decomposition (POD), are used in the data-assimilation scheme [START_REF] Mons | Data assimilationbased reconstruction of urban pollutant release characteristics[END_REF]), and to determine the fluid flow on a POD-based lower order system [START_REF] D'adamo | Variational assimilation of POD low-order dynamical systems[END_REF]).

We remark that, although those two techniques differ in nature, some hybrid ones may exist as well. For example, in order to bypass the need for an adjoint solver, necessary for VDA, we may consider instead the idea of ensembles, where, in order to compute the sensitivity of the cost functional with respect to the control parameters, a cloud of points is propagated, each one of those points presenting a different outcome. This technique is called Ensemble-based Variational assimilation (EnVar). It was first designed by [START_REF] Liu | An ensemble-based four-dimensional variational data-assimilation scheme. Part i: technical formulation and preliminary test[END_REF] in the context of weather prediction. More recently, [START_REF] Mons | Reconstruction of unsteady viscous flows using data assimilation schemes[END_REF]), performed some studies on this technique, in comparison with EnKF and VDA in the context of Fluid Mechanics.

Goal of the Present Work

At this point, data-assimilation has been presented in a very general manner, together with a few applications it may have in Fluid Mechanics. In what follows, we intent to present a more specific overview of the precise problem we will consider in this work. This particularization will come mainly from the type of problem we will consider as numerical model and on the measures that will be the most compatible with the problem.

In most of the engineering applications, for example in aerodynamics, the timeaveraged (mean) flow is sufficient for the global understanding of the fluid flow, allowing us to design mechanical components, to evaluate loads on structures, etc. Indeed, in most of industrial wind tunnel facilities, the nature of the data extracted is usually timeaveraged. Examples of such data are the mean wall-pressure or friction distributions or mean-velocity fields, obtained with either point-wise probes (such as Pitot) or Particle Image Velocimetry (PIV) techniques. From the numerical point of view, a very common tool used in industry is the Reynolds-Averaged Navier-Stokes (RANS) equations, whose solution represents the mean-flow. Since those equations are steady in nature, their solution does not require a fine time-stepping procedure, minimizing the computational effort to solve them. This makes them very popular in industry, where the procedure of engineering design often demands easily obtainable solutions. The RANS equations are formally derived from the Navier-Stokes ones, by applying the Reynolds decomposition of the fluid flow variables, separating them into a mean-flow component and a fluctuation around it. Due to nonlinearities, the mean-flow variables in RANS equations will depend on the statistics of the fluctuation, given by the Reynolds-stress tensor. Those statistics cannot be explicitly computed since the fluctuation is not known, once we solve only the mean-flow component. Several models exist that attempt to provide an expression of this tensor in terms of mean-flow quantities and possibly of supplementary variables. This is the so-called closure problem. Generally those models are build using physical considerations such as eddy-viscosity, mixing-lengths, scaling laws, boundary layer theory, etc. This physical reasoning, in general, is capable of fixing somewhat the structure of the model, but may leave a set of parameters or coefficients to be specified. The specification of those parameters is done by considering some benchmark physical configurations, providing a calibration for them. This means that, for those configurations, the models are designed to work well. However, if we consider other configurations far away from those design points, typically with strong separation regions, the model may lose its accuracy when compared with experimental or higher-fidelity numerical data.

The main goal here is then to modify those models (or propose others) in a way that their solution matches the best some higher-fidelity data. We can see that this kind of problem can easily be posed in the logic of Variational Data-Assimilation (VDA) by identifying a tuning parameter, related with the Reynolds-stress tensor itself or with the closure model. This is the main goal of this work. In what follows, we present the brief state-of-the-art of mean-flow data-assimilation, together with a common application of it, the mean-flow stability analysis (in a broad sense), which will also be explored in the present work. This analysis is capable to provide an approximation for the nonlinear fluctuation around this mean-flow.

Mean-flow data-assimilation

The mean-flow data-assimilation, as presented in the previous section, has been considered in the past few years mainly from the point of view of variational data-assimilation, where some parameters related to the Reynolds-stress tensor are tuned such that the solution of RANS equations matches the best with external data. [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF] proposed a model consisting in the laminar, steady Navier-Stokes, where the correction parameter is a volume-force modeling directly the divergence of the Reynolds-stress tensor. They showed that, in a laminar / transitional context (Re = O(10 2 )) that this model performs well to assimilate synthetic Particle Image Velocimetry (PIV) data. The same model was applied in a higher-Reynolds number context (Re = O(10 4 ), [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF]).

They reported some difficulties related to the use of the laminar equations in a higher Reynolds number configuration. To circumvent those difficulties, a series of increasing Reynolds number data-assimilation needed to be performed, which can be costly. Other approaches, more appropriate for high-Reynolds number flows, rely on a RANS turbulence model. For example [START_REF] Li | A data-driven adaptive Reynolds-averaged Navier-Stokes κ -ω model for turbulent flow[END_REF] tuned the κ -ω model's parameters such that its solution is the closest to reference ones. [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF], on the other hand, proposed a field inversion method where a spatially-dependent parameter related to the eddy-viscosity production term is tuned. In their work, the output solution was used for machine learning techniques such that the considered model becomes improved by learning the modifications caused by the data-assimilated procedure. We remark that, although those techniques are based on Variational Data-Assimilation, ensemble-based assimilation can also be used in a steady (mean) flow simulations contexts. [START_REF] Kato | Approach for uncertainty of turbulence modeling based on data assimilation technique[END_REF] used ensemble techniques by varying the parameters (constants) in the Spalart-Allmmaras model, and a further Kalman filter, where the external measure was applied, leading to the converged optimal set of parameters. We remark that this kind of approach is somewhat equivalent to variational assimilation, where if the statistics are generated in an iterative procedure, with a large enough ensemble, VDA and KF solutions will be similar [START_REF] Iglesias | Ensemble Kalman methods for inverse problems[END_REF]).

Beyond mean-flow: stability analysis

In this paragraph we make a few remarks on the mean-flow and how it can be used to investigate aspects of the nonlinear fluid flow other than mean velocity-field, mean wall-pressure, etc. It was shown that, stability analysis (in a broad sense) using the mean-flow, instead of the steady-state solution, could provide an approximation of the unsteady nonlinear fluctuation around this mean-flow. Those studies have been formalized by [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF], who applied this technique to the flow around a circular cylinder at low Reynolds number. He showed that the nonlinear solution corresponds to a periodic limitcycle, that can be decomposed in a mean-flow and on a fluctuation around it. The stability analysis of this mean-flow leads to an eigenvalue whose real part was zero and whose imaginary part matched very well with the nonlinear frequency. Later, Sipp and Lebedev [2007a] and [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF] provided conditions for the validity of this kind of analysis in the context of oscillator flows (flows whose dynamics in not dependent on external forces) and showed that, under this conditions, the eigenmode provided an approximation for the first harmonics of the fluctuation, represented by the periodic limit-cycle. An extension of this analysis for flows whose eigenspectrum does not indicate unstable or neutral mode (for example, "noise amplifiers") was investigated by [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF], [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]. They proposed the use of Singular-Value decomposition of the Resolvent operator (transfer function of the mean-flow linearized system) rather than the eigenvalue/eigenvector analysis. They showed that, at a given frequency, if the leading Singular Value is much higher that the others, the nonlinear fluctuation could be approximated by the leading (left) Singular Vector at that frequency. They argued that this condition is verified in most flows presenting a strong instability mechanism. We remark however that, since this analysis is linear, it is not capable to determine alone the amplitudes (or energies) of those modes at each frequency. For this reason, most of the studies aiming to obtain the actual fluctuation needed time-resolved data, which is usually point-wise information, from which we can extract, for each frequency, the energy of the mode, together with its phase. Furthermore, if this low-rank representation of the Resolvent operator holds, in principle only one (or a few, for robustness) time-resolved probe is sufficient for determining those amplitudes. We remark that even if this kind of measure is not time-averaged ones, the very low quantity of those probes still leaves this kind of technique feasible. For this reason, in the past few years, the determination of the mean-flow by variational data-assimilation and its further Resolvent-based analysis is increasingly popular technique, since it allows one to obtain an approximation of the full flow field from mean-flow measurements and very few time-resolved probes (see Symon et al. [2019a], [START_REF] Symon | Mean and unsteady flow reconstruction using data-assimilation and resolvent analysis[END_REF], He et al. [2019a]).

Overview of the Thesis

In this work, we will investigate a few techniques for mean-flow data-assimilation, together with mean-flow analysis. More precisely, in the chapters 2 and 3, we will consider the modification of the RANS equations for mean-flow data-assimilation purposes, by tuning some parameters of the model. In the chapters 4, 5 and 6, we will deal with mean-flow analysis and Resolvent-based data-assimilation. A detailed description of each chapter individually is provided:

• Chapter 2: In this chapter, we will discuss the mean-flow data-assimilation with the aid of the RANS equations, completed by a turbulence model, the Spalart-Allmaras. More precisely, we will introduce two different tuning parameters related to a correction of those equations. The first will be a volume force, following [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF], modeling the corrective force supplementary to the one coming from the eddy-viscosity provided by the SA model. We note that the high-Reynolds number issues reported by [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF] will not occur in our case since we dispose of a backup model (SA), providing enough eddy-viscosity such that the unforced solution (the original RANS-SA one) is well posed. This way, we avoid the need for iterative assimilation for increasing Reynolds numbers. The second one is a corrective source term acting in the SA equation (similar to [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF]), modifying how the production or destruction of the turbulent eddy-viscosity was designed. The data to be assimilated will be the mean-velocity field, either in the whole domain, to investigate controlability aspects of the proposed model, or at partial/sparse velocity measurements, as we would have using point-wise Pitot probes. By doing so we assess how those models perform in function on the quantity of information. The flow configuration will be a Backward-Facing Step (BFS) at Re = 28275, investigated in [START_REF] Dandois | Numerical simulation of active separation control by a synthetic jet[END_REF] with Direct Numerical Simulations (DNS), which will be used as ground-truth.

• Chapter 3: In this chapter we essentially investigate the same models as in the Chapter 2 but changing now the nature of the data to be assimilated. Here, we consider wall-information (pressure and friction). In a first time, we derive the adjoint equations for the laminar regime (without turbulence model), as their boundary conditions will be non standard, due to the wall-measurements. This formalism will also be, at first, tested in a laminar configuration. Then, we perform those data-assimilation on the turbulent case presented on Chapter 2, studying how the 4. Overview of the Thesis 7 different models manage to assimilate this type of data, together with a few studies on the number of wall-probes necessary for convergence.

• Chapter 4: The mean-flow Resolvent analysis, as presented in §3.2, has been extensively used for the recovering of the nonlinear fluctuation. Most of the studies consider the mean-flow as a given quantity, both directly from experiments (or higherprecision numerical simulations) or from previous mean-flow data-assimilation. However, we remark that, if we know an approximation for the fluctuation, we can also have an approximation of the mean-flow. Indeed, from this fluctuations, we can compute the Reynolds-stress tensor and thus the forcing term on the RANS equations, that we could use to compute the mean-flow itself. For this reason, we propose here to skip the a priori knowledge of the mean-flow, proposing instead a nonlinear model composed by the mean-flow equation coupled with the Resolvent modes through the Reynolds-stress tensor. This model is not closed since, as mentioned before, the amplitudes of each mode (for each frequency) is unknown and will be tuned such that some quantity of the reconstructed flow matches the best with some external ones. Note that this measure can, in principle, be the same few time-resolved point-wise probes used by [START_REF] Beneddine | Unsteady flow dynamics reconstruction from mean flow and point sensors: an experimental study[END_REF] or Symon et al. [2019a], since their frequency spectra provide, individually for each frequency, the point-wise energy of the fluctuation. Moreover, those quantities may provide the information about which frequencies our model should have. For this reason, from the point of view of the input information needed, this as an extremely efficient data-assimilation technique. This technique will be applied on transitional flows: a squared-section cylinder, an example of flows whose instabilities saturate into a periodic limit cycle (oscillator flows), and a backward-facing step, an example of flows driven by external noise (noise amplifiers).

• Chapter 5: In this chapter we investigate the mean-flow stability analysis on flows whose fluctuation is composed by a quasi-periodic component and a broadband one. The application case is a squared-section cylinder at Re = 22000, whose quasiperiodic component is in fact represented by the periodic vortex-shedding and its broadband fluctuation is partially composed by Kelvin-Helmholtz (KH) structures located at the shear layer generated from the detachment point. To obtain both structures from mean-flow analysis (obtained from DNS), we decompose the flow variables using the triple decomposition as proposed by [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF], where the Koopman formalism is used, leading to two different equations, one for each one of them. The periodic vortex-shedding mode is obtained with a classical mean-flow stability. We investigate here the need of taking in consideration an eddy-viscosity, modeling the broadband fluctuation's nonlinear interactions and computed from the DNS statistics. The Kelvin-Helmholtz structures, on the other hand were obtained with a Resolvent analysis based on the linearized equations around the periodic component, instead of the mean-flow, similarly to what is done in Floquet stability theory. By using this formalism, we are able to capture the oscillation of those KH structures with the periodic vortex-shedding, which is also observed in DNS.

• Chapter 6: In this chapter, we intend to extend the Resolvent-based data-assimilation technique (presented in chapter 4) to the turbulent flow around the square cylinder, same configuration presented in chapter 5. This is a much more challenging task since now we have both periodic and broadband components to be taken into account. Numerically, this is a costly task and will be circumvented here with the use of a URANS model. The chosen one is the Spalart-Allmaras. This model provides naturally an eddy-viscosity, necessary for the resolution of the vortex-shedding mode and for the mean-flow equation, modeling on part of the Reynolds-stress tensor. The final model is composed by the vortex-shedding mode, whose amplitude is tuned similarly to what was done in the third chapter, coupled with the mean-flow equations.

Introduction

Experimental measurements are commonly used to validate numerical simulations. Yet, experimental measurements and numerical simulations can also be more intimately combined to generate an improved picture that makes a compromise between numerical modeling and experimental measurements. In this view, data-assimilation aims at establishing a mathematically well-defined procedure to take advantage of the knowledge of few experimental measurements to complete or correct an incomplete or imprecise numerical model. Data assimilation historically emerged in the field of meteorological forecasts [START_REF] Lorenc | Analysis methods for numerical weather prediction[END_REF], [START_REF] Liu | An ensemble-based four-dimensional variational data-assimilation scheme. Part i: technical formulation and preliminary test[END_REF]), where the lack of accurate models, uncertain initial and boundary conditions yield poor predictions. Several data-assimilation algorithms were designed that can be classified into two categories: direct-adjoint optimization techniques, also called 3D/4D Var, where correction terms or boundary conditions in the numerical model are tuned to minimize the measurement mismatch between the experiment and the numerical solution Le [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspect[END_REF], and ensemble based-techniques, where a state estimate is corrected using the measurement mismatch [START_REF] Evensen | Data-Assimilation: yhe ensemble Kalman filter[END_REF]. In other areas of fluid mechanics, the use of these kind of techniques has been increasing in the past years. The motivations for this were multifold: estimation of initial or inlet boundary conditions in open flows Gronskis et al. [2013], [START_REF] Mons | Reconstruction of unsteady viscous flows using data assimilation schemes[END_REF], interpolation of velocity fields between a sequence of images [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF], [START_REF] Yang | Enhanced ensemble-based 4DVar scheme for data assimilation[END_REF], identification of pollutant release in urban area [START_REF] Mons | Data assimilationbased reconstruction of urban pollutant release characteristics[END_REF] and even as a theoretical tool to investigate aspects of the decay of large scales in homogeneous isotropic turbulence [START_REF] Mons | Is isotropic turbulence decay governed by asymptotic behavior of large scales? An eddy-damped quasi-normal Markovian-based data assimilation study[END_REF]. In aerodynamics, the recovery of time-averaged quantities is in many applications sufficient. For example, in industrial wind-tunnels, most of the measurements consist in time-averaged pressure distributions and forces and from a numerical point of view, engineers commonly solve Reynolds-Averaged Navier-Stokes equations with turbulence models to predict the mean features of the flow around a body. In so far, a time-averaged approach is more suited and the data assimilation problem should rather be considered in the framework of inverse problems (see [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF] for the direct-adjoint optimization approach and [START_REF] Iglesias | Ensemble Kalman methods for inverse problems[END_REF], [START_REF] Kato | Approach for uncertainty of turbulence modeling based on data assimilation technique[END_REF] for the ensemble-based framework). Concerning the direct-adjoint optimization approach, [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF] started with a low-Reynolds cylinder flow (Re = O(10 2 )) 10

Chapter 2. Data-Assimilation for minimal Corrections of RANS equations exhibiting vortex shedding. They tuned a volume-force (modeling the force associated to the Reynolds-stress) acting in the steady Navier-Stokes equations such that its corresponding solution best matches velocity measurements, mimicking a real experimental situation where such measurements are provided by a Particle Image Velocimetry setup. The same procedure was applied at a higher Reynolds number (Re = O(10 4 )) by [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF] on a flow around an idealized airfoil. This time, additional difficulties related to the well-posedness of the steady Navier-Stokes equations at such high Reynolds number had to be faced. They overcame this problem by performing several optimization steps for increasing Reynolds numbers, which can be costly. Another possible approach at high Reynolds-numbers is to consider turbulence models. This was, for example, employed by [START_REF] Li | A data-driven adaptive Reynolds-averaged Navier-Stokes κ -ω model for turbulent flow[END_REF], where a set of model-related coefficients (in their work, RANS-k -ω) have been tuned such that the solution of the model is as close as possible to some given partial higher-fidelity-data. Yet, such an approach is strongly constrained by the structure of the model and does not allow a high-enough flexibility to adjust the model to regions exhibiting different behaviors. Concerning the ensemble approach, [START_REF] Kato | Approach for uncertainty of turbulence modeling based on data assimilation technique[END_REF], [START_REF] Kato | A data assimilation methodology for reconstructing turbulent flows around aircraft[END_REF] generated an ensemble by varying abritrarily the coefficients of the Spalart-Allmaras model and a Kalman filtering method than led to the tuning of optimal values of these coefficients with respect to given measurements. The difficulty here is again that the Reynolds stresses need to comply with the structure of the Spalart-Allmaras model (in particular the Boussineq approximation) which is known to be very constraining. [START_REF] Duraisamy | Turbulence modeling in the age of data[END_REF] also employed data-assimilation techniques to tune production terms in RANS models to mean-flow data obtained by DNS or experiments. They then used the assimilated eddy-viscosity quantities to design new turbulence models with deep-learning techniques.

In this opening chapter, we extend the time-averaged direct-adjoint optimization approach introduced by [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF], [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF] to higher Reynolds numbers by using as a baseline model RANS equations combined with turbulence models instead of the steady (laminar) Navier-Stokes equations. Contrarily to [START_REF] Li | A data-driven adaptive Reynolds-averaged Navier-Stokes κ -ω model for turbulent flow[END_REF], the tuning parameters are space-dependent functions related to a correction of the modelled Reynolds-stress tensor. In the present work, we choose the Spalart-Allmaras model (SA) [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] for its numerical simplicity, even though other models could also be used. Within this model, we explore two correction possibilities.

The first consists in a volume-force acting in the momentum equations. This force is supposed to stand for a correction of the viscous force induced by the eddy-viscosity produced by the SA model. This therefore allows the turbulent stresses not to be constrained by the Boussinesq assumption, which is known to be only well adapted in shear dominated regions. Since this correction is not attached to the turbulence model itself, this approach may be considered as an extension of [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF] and [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF], the difference being the consideration of a background eddy-viscosity that evolves with the data-assimilation process according to a turbulence model (in this work, SA). This modeling is better suited for turbulent flows. For example, the uncorrected solution (for which the correction term is null), which is needed to initialize the optimization algorithm, now corresponds to the RANS-SA solution. Such an initial solution may easily be obtained with standard algorithms and is already much closer to the targeted mean-flow than an initial condition obtained with a laminar model. The assimilation procedure therefore only needs to compensate the small imperfections of the initial turbulence model, while using a laminar model for turbulent flows implies compensating for the whole discrepancy between the laminar steady solution and the targeted mean-flow.

The second correction possibility is a source term acting on the SA equation that governs the eddy-viscosity. This time the control parameter acts within the turbulence model, correcting the production and dissipation terms of the used turbulence model. It is worth-mentioning that a similar technique was employed by [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF] and [START_REF] Parish | A paradigm for data-driven predictive modeling using field inversion and machine learning[END_REF] to identify correction terms in turbulence models: they tuned a space-dependent function modulating the strength of the turbulent-viscosity production term. More recently, this very same procedure was used by [START_REF] He | Data assimilation and resolvent analysis of turbulent flow behind a wall-proximity rib[END_REF] to assimilate some PIV measurements. We believe however that an additive source term in the equation governing the eddy viscosity is more general and robust.

This chapter is organized as follows. First ( §2), we will describe the physical configuration of interest, by showing the reference solution (obtained by DNS) and the corresponding (uncorrected) RANS-SA solution, pointing out the differences between them, motivating the need for data-assimilation. This configuration consists of a rounded backward-facing step [START_REF] Dandois | Numerical simulation of active separation control by a synthetic jet[END_REF]), for which the SA model provides a bad prediction of the re-circulation length. Then ( §2), we will introduce the two models for dataassimilation mentioned above, together with the direct-adjoint optimization procedure to tune the correction fields. Finally ( §4), we will present the data-assimilation results on the backward-facing step configuration, by comparing the assimilated flowfields obtained with the two approaches to the true mean-flow. We will also assess the performance of the reconstruction as a function of the sparsity of the used velocity measurements.

Flow Configuration and numerical solutions

Before introducing the mathematical models and the formalism for data-assimilation, we briefly introduce the physical configuration on which we will apply this technique, together with the reference (DNS) and a priori (RANS-SA) solutions. We will in particular highlight and discuss the errors of the (RANS-SA) solution, motivating the data-assimilation procedure.

The chosen configuration, is a rounded Backward-Facing Step (BFS) at Reynolds number Re = 28275, based on the height of the step and on the inflow velocity. Those quantities will also be used to make all variables dimensionless in the following. This configuration was considered in [START_REF] Dandois | Numerical simulation of active separation control by a synthetic jet[END_REF] and consists of a smooth expansion, generating a re-circulation region that detaches at x d = 0.53 and reattaches at x r = 3.93, those positions being measured after the start of the expansion. The geometry of the rounded step is given by the equation:

y(x) = [sin(aπx) -aπx]/(2π) + 1, 0 ≤ x ≤ 2/a (2.1)
where a = 0.703. As detailed in [START_REF] Dandois | Numerical simulation of active separation control by a synthetic jet[END_REF], the Direct Numerical Simulation (DNS) was fed at the inlet with fluctuations generated from a time-dependent zero-pressure-gradient turbulence simulation (see [START_REF] Lund | Generation of turbulent inflow data for spatially-developing boundary layer simulations[END_REF]). This allows to describe a situation where the incoming boundary layer is fully turbulent. The resulting mean-flow can be seen in Figure 2.1. We now turn our attention to the solution given by the RANS equations coupled with the Spalart-Allmaras model. Here we consider the incompressible Navier-Stokes equations, describing the velocity ũ and pressure p (the • notation referring to quantities obtained by a model), together with an extra equation governing an

ũ • ∇ũ + ∇p -∇ • ((ν + ν t (ν))∇ s ũ) = 0 (2.2) ∇ • ũ = 0 (2.3) ũ • ∇ν -∇ • (η(ν)∇ν) = P (ν, ∇ũ) + D(ν, ∇ũ) + C(∇ν) (2.4) = s(ν, ∇ν, ∇ũ),
where η(ν) is the diffusivity of ν, P (ν, ∇ũ) the production term, D(ν, ∇ũ) the destruction term and C(∇ν) the cross diffusion term, all of them being lumped into a single source term s(ν, ∇ν, ∇ũ). More information about the original baseline model (especially the functions ν t (ν), P (ν, ∇ũ), D(ν, ∇ũ) and C(ν, ∇ũ)) and the chosen implementation can be found in [START_REF] Crivellini | A Spalart-Allmaras turbulence model implementation in a discontinuous galerkin solver for incompressible flow[END_REF] and [START_REF] Oliver | A High-Order and Adaptive and Discontinuous Galerkin Finite Element Method for the Reynolds-Averaged Navier-Stokes Equations[END_REF], as well as in Appendix A. To properly solve this set of equations, we need to specify boundary conditions, especially the inflow velocity profile (a sketch of the flow domain, with the remaining boundary conditions is provided in figure 2.2). We will consider here that the inflow boundary layer is known and mimics the inflow conditions of the DNS (this constraint can in principle be removed by tuning or assimilating this inflow boundary condition with respect to some measurements more downstream; yet, this is not the focus of this work and has not been studied). One of the issues here concerns the value of ν, which should reflect the value of the shearstress in the DNS. For this variable, we will keep it simple and impose a constant value (ν/ν) ∞ with a fast-decay to zero very close to the wall. We have made a parametric study for values of this quantity ranging from (ν/ν) ∞ = 3 (recommended for a fully-turbulent boundary layer, see Allmaras et al. [2012]) to (ν/ν) ∞ = 0. The stream-wise velocity, the normalized eddy-viscosity ν t /ν, the wall-pressure and friction coefficients of the RANS-SA solution for (ν/ν) ∞ = 3 are respectively shown in Figures 2.3 (a,b,c,d). The dependency of the incoming boundary layer (at the step location) as a function of (ν/ν) ∞ is shown in Figure 2.3 (e), indicating a small sensitivity with respect to the choice of (ν/ν) ∞ . In table 2.1, we can see that the re-circulation length is strongly overestimated by the RANS-SA model and that the choice of (ν/ν) ∞ strongly affects these overestimated results. This overestimation may be due to the fact that the Spalart-Allmaras model has been initially tuned for higher Reynolds numbers and the gradients are not sufficiently strong in the present shear-layer to strengthen the eddy-viscosity (and reduce the re-circulation length). For this reason, throughout this work, we fix that value to be (ν/ν) ∞ = 3, which is consistent with the fact that the DNS itself was performed with a fully-turbulent incoming boundary-layer. As mentioned above, the RANS-SA solution over-predicts the reattachment length by 70%. This also reflects in a strong mismatch in the pressure and friction coefficients (Figure 2.3 (c,d)). These discrepancies motivate the use of data- assimilation techniques to improve the recovery of the right flow features. In this study, we will consider the time-averaged DNS as the reference solution from which measurements will be extracted to feed the assimilation process.

Models and Data assimilation procedure

In this section, we will discuss the models we consider for the data-assimilation process.

We recall that the time-averaged u and pressure p quantities are governed by the RANS equations:

u • ∇u + ∇p -∇ • (ν∇ s u) = -∇ • u ⊗ u = -∇ • τ = f , (2.5)
where τ is the Reynolds-stress tensor, and the notation u ⊗ v stands for the tensor

[u i v j ] i,j=1,••• .
This tensor contains all the nonlinear interactions of the non-resolved fluctuations u . For this reason, those equations are not closed and several hypotheses / approximations / models exist in the literature to close the problem, that is provide a closure τ = τ (u, p, + supplementary variables ). It is thus clear that quantities related to this closure can legitimately be used as tuning functions for our data-assimilation procedure. In the following, we will describe two possible tuning functions, both of them relying on the SA-turbulence model. The first one corresponds to adding a volume-force in the momentum equations and the second one to adding a volume source in the ν equation.

Model 1: source term in momentum equation

The first model considered here is given by a volume force on the momentum equations. However, using the Navier-Stokes equations with solely a molecular viscosity as a basemodel and compensating the measurement mismatch through only a volume-force (as done in [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF]) may be prone of difficulties at high Reynolds numbers for several reasons: [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF] had to perform several expensive optimization steps at increasing Reynolds numbers to deal with a case at Re = O(10 4 ), the correction term may be very large since the uncorrected laminar steady state is usually very far from the turbulent mean, and in some cases, this laminar state may not even exist or there might be lots of such steady states [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF]. All these reasons render the assimilation process difficult and expensive. In any case, it seems a good idea to have an assimilation process with an initial guess that is easy to obtain and not too far from the targeted solution so that only a small discrepancy needs to be compensated by the optimization process.

In order to keep the physical idea of a volume-force in the momentum equations and to improve the robustness of the assimilation process, we choose to use existing turbulence 3. Models and Data assimilation procedure models in connection with a momentum force in the RANS equations to serve as a model for the optimization. This means that the Reynolds stress tensor will be divided into a modeled and a tuned part:

f = -∇ • τ ≈ -∇ • (2/3κI -ν t ∇ s ũ) + fx .
(2.6)

The first term on the right-hand-side of equation (2.6) corresponds to the Boussinesq term, stating that part of the Reynolds-stresses can be modeled by an eddy-viscosity ν t (the term involving the turbulent kinetic energy κ = u • u /2 being incorporated in the pressure). The second term is the correction term fx , which can in principle account for non-Boussinesq effects or other modelling errors. With this approximation, instead of (2.2), we consider the following forced momentum equations:

ũ • ∇ũ + ∇ (q + 2/3κ) p -∇ • ((ν + ν t )∇ s ũ) = fx , (2.7)
where the turbulent kinetic energy p = q + 2κ/3 has been lumped into the pressure term.

The eddy viscosity ν t present in these equations is given by a turbulence model, here the Spalart-Allmaras model ν t (ν) where ν is governed by (2.4).

Model 2: source term in ν equation

The second tuning function considered corresponds to a source term in the equation governing ν. Equation (2.4) is replaced by ũ • ∇ν -∇ • (η(ν)∇ν) = s(ν, ∇ν, ∇ũ) + fν .

(2.8)

The correction fν only has an indirect effect on the velocity ũ and pressure p fields, since fν only modifies the source-production-diffusion term s acting on ν, which then modifies the eddy-viscosity ν t that appears in the momentum equations (2.3). The Reynolds stress tensor therefore remains constraint by the Boussinesq assumption, which is nevertheless a valid assumption in a great variety of flows, especially in shear-dominated flows. This model is sufficiently flexible to allow for example correction of a wrong free shear-layer expansion ratio, as is the case in figure 2.3 (a). Model 2 is nevertheless far more constrained than model 1, since with model 1 we are able to encompass all corrections obtained by model 2, while the reverse is not true. Yet, this should not just be seen as drawback for model 2, since the associated solution will be forced to remain in a subspace that might in some cases be more physical (or representative of the exact DNS solution) than the solution obtained with a much more free model, such as model 1. One has to keep in mind that depending on the number of measurements that are available, the inverse problem that aims at determining the tuning function that yields the best measurement match is more or less severely under-determined: for example, in presence of sparse (resp. many) measurements, it might be better to pick model 2 (resp. model 1) that is more constrained (resp. freer).

Data Assimilation procedure

In this section, we describe the data-assimilation procedure employed. Let m be a set of higher-fidelity or experimental measurements that correspond to information extracted

Chapter 2. Data-Assimilation for minimal Corrections of RANS equations from the flow and M(•) the measurement operator that allows to extract the corresponding measure from our simulation result ũ. In this chapter, since we are dealing only with velocity measurements, this operator will act on the velocity field, ũ, yielding m = M(ũ) ∈ M (other measure operators will be considered in the next chapter). M is the measurement space, whose norm is given, generically, by ||•|| M . The data-assimilation problem can now be cast into an optimization one, for which we want to tune the forcing terms (either fx or fν ) such that the cost functional:

J (q) = 1 2 ||M(q) -m|| 2 M (2.9)
is minimal, with the following constraints:

ũ • ∇ũ + ∇p -∇ • ((ν + ν t (ν))∇ s ũ) = fx (2.10) ∇ • ũ = 0 (2.11) ũ • ∇ν -∇ • (η(ν)∇ν) -s(ν, ∇ν, ∇ũ) = fν ,
(2.12)

where here, for conciseness, we consider a general framework in which both fx and fν could simultaneously be considered in the optimization process. To solve this optimization problem, following [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF], [START_REF] Parish | A paradigm for data-driven predictive modeling using field inversion and machine learning[END_REF], [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF] or [START_REF] Mons | Is isotropic turbulence decay governed by asymptotic behavior of large scales? An eddy-damped quasi-normal Markovian-based data assimilation study[END_REF], we choose a gradient-based algorithm. For the descent algorithm, we pick the low memory version of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, which is far more efficient than a simple steepest gradient-descent one. It requires in particular the computation of the cost functional gradient with respect to the forcing parameter, ∇ fx J or ∇ fν J . For this, we resort to a Lagrangian formalism, in which the state is augmented with a set of Lagrange multiplier (or adjoint) variables (ũ † , p † , ν † ) and we look for critical points of the Lagrangian functional:

L([ũ, p, ν], [ũ † , p † , ν † ], [ fx , fν ]) =J (ũ) + ũ † , ũ • ∇ũ + ∇p -∇ • ((ν + ν t (ν))∇ s ũ) -fx Ω + p † , ∇ • ũ Ω + ν † , ũ • ∇ν -∇ • (η(ν)∇ν) -s(ν, ∇ν, ∇ũ) -fν Ω , (2.13) where (v 1 , v 2 ) Ω = Ω v 1 • v 2 dΩ
represents the inner product related to the classical (real valued) energy-norm. Setting the variation of this Lagrangian with respect to the adjoint variables to zero, yields the governing equations (2.10, 2.11, 2.12). Setting its variation with respect to the direct variables, we obtain the adjoint equations:

ũ † • (∇ũ) T -ũ • ∇ũ † -∇ • ((ν + ν t )∇ s ũ † ) -∇p † (2.14) +ν † ∇ν + ∇ • (ν † ∂ ∇ũ s) = - ∂M ∂ ũ † (M(ũ) -m) ∇ • ũ † = 0 (2.15) -ũ • ∇ν † -∇ • (η∇ν † ) + (∂ ν η)∇ν † • ∇ν + (∂ ν ν t )∇ũ † : ∇ s ũ -(∂ ν s)ν † + ∇ • (ν † ∂ ∇ν s) = 0, (2.16) 4. Results

17

where the boundary conditions for this adjoint variables are given by:

N o -slip/Inf low : ũ † = 0, ν † = 0 (2.17)

Outf low : n • ũ † ⊗ ũ + (ν + ν t )∇ s ũ † + I p † + (∂ ∇ũ s)ν † = 0 (2.18) n • ũν † + η∇ν † -ν † (∂ ν η)∇ν -(∂ ∇ν s)ν † = 0 (2.19) Slip : ∂ y ũ † x = 0, ũ † y = 0, η∂ y ν † = (∂ ∂y ν s)ν † .
(2.20)

Although here we derive the strong equations for the adjoint problem, together with its boundary conditions, in practice we use the discrete adjoint, meaning the transpose of the Jacobian matrix. This procedure is fully justified in a Finite-Element Framework, since it provides an actual FEM discretization of the adjoint equations. For further details on this issue, see Appendix B 4. Taking now the variation of the Lagrangian with respect to the forcing terms (either fx or fν ), we have:

∂L ∂ fx , δ fx Ω = -ũ † , δ fx Ω → ∇ fx J = -ũ † (2.21) ∂L ∂ fν , δ fν Ω = -ν † , δ fν Ω → ∇ fν J = -ν † . (2.22)
At this point we have all the necessary ingredients for the BFGS optimization algorithm.

A brief description of the iterative procedure is depicted in figure 2.4. Roughly speaking, this algorithm is more elaborate than usual gradient-descent methods, since the descent direction is obtained by inverting an approximation of the Hessian matrix, constructed from previous gradient evaluations. This procedure allows one to capture the second-derivative behaviour of the cost functional J , which, formally provides a second-order convergence of the optimization (further details on its implementation with a finite-element solver in Appendix B). It is interesting to notice that, in the case where we use the fx tuning parameter, even if the original gradient is modified by the Hessian matrix, the finally chosen descent direction is still divergence-free, just as in [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF] for the adjoint velocity field ũ † . However, since we now have a background eddy-viscosity, the total force in 2.6 that acts on the momentum equations is not anymore divergence-free. We remark as well that, at any iteration of the optimization procedure, the forcing term fx may have a non zero projection onto the space of all forces compatible with a Boussinesq hypothesis. This means that the eddy-viscosity ν t = ν t (ν) does not hold anymore the 'total amount' of eddy-viscosity, a part of which can as well be contained in fx . This slightly hides the physical meaning of the sole term ν t .

Results

We will now analyze the reconstruction results of the backward-facing step flow for two different measurement operators. Firstly, we will consider the complete mean velocity field as known in order to evaluate the ability of the two models to accurately reconstruct the velocity, pressure and Reynolds stress forcing. This is interesting because the nature of the forcing terms in the two models is different, which induces different constraints on the reconstructed fields, so that the ensemble of reachable velocity-pressure fields is different between the two models. For example, the mean-flow generated by the DNS / Figure 2.4: Sketch of the BFGS algorithm for a general problem, having f as a control parameter and q as state. The adjoint variable is denoted by q † . LES solutions might in fact not be fully compatible with the Boussinesq approximation of model 2. Secondly, we will use only few point-wise velocity measurements such as those given by a Pitot probe or a hot-wire measurement that is moved along a line to extract cross-stream velocity profiles. Such a situation is actually much more realistic when considering industrial wind-tunnels, since optical measurements such as Particle-Image-Velocimetry (PIV) measurements remain difficult to implement in such facilities. We actually expect different behaviors between model 1 and model 2 when only few measurements are available since the reconstruction constraints are different between the two models. [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF] and [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF] have already explored point-wise measurements in the case of a laminar model at lower Reynolds numbers: in particular, [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF] showed that a measurement operator involving a spatial averaging of the measure over a finite region improved the quality of the reconstruction by smoothing the gradients of the reconstructed velocity field in the vicinity of the measurements. We will apply here a similar approach.

Complete velocity field measurements

We consider the measure as being the whole exact (mean) velocity field m = u. This field was extracted from the DNS results of [START_REF] Dandois | Numerical simulation of active separation control by a synthetic jet[END_REF] and then interpolated on our mesh, which is coarser than the one from the DNS. The measure is thus defined on the same mesh as the solution itself. For this reason, the measurement space M ends up being the velocity space itself. Its norm also simplifies into || • || M = || • || Ω , the norm induced by the kinetic energy inner product (without the pressure component). We then have:

- ∂M ∂ ũ † (M(ũ) -m) = -(ũ -u) (2.23)
Under those conditions, the optimization of the tuning fields fx and fν proceeds as shown in figure 2.5. Model 1 with the optimization of the volume-force fx manages to reconstruct the time-averaged DNS solution extremely accurately due to the large flexibility of this model: the decrease in the cost functional J n /J 0 is around 10 -4 , which is very good. Note also that this decrease is achieved in less than 200 iterations, which is much quicker than the 10 3 iterations required by [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF] to achieve similar results. This gain is due to the choice of far-superior l-BFGS optimization algorithm, which provides an approximation of the Hessian matrix based on previous gradient evaluations, improving the descent direction estimation. When using model 2, the optimization procedure reaches a plateau after n ≈ 10 iterations. In this case, the reduction of the cost functional is of the order of J n /J 0 ≈ 0.13. This is less accurate than with model 1. Yet, it is coherent with the fact that the control parameter fν only indirectly affects the velocity field ũ (through was reported by [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF], who explained this discrepancy by the fact that the adjoint fields (and thus the forcing term fx ) are divergence-free by construction if no pressure measurements are considered in the cost functional. This constraint ends up preventing the pressure to be fully reconstructed, since we have no information on the potential part of the forcing term. It is interesting to notice however that, despite this fact, the wall-pressure (figure 2.8 (a)) is well represented. To explain this, we turn our attention to the exact RANS equations (2.5) and the approximated ones (2.7). We can see that subtracting the two equations yields:

∇(p -p) + (ũ • ∇ũ -u • ∇u) -∇ • (ν∇ s (ũ -u)) = fx + ∇ • (ν t ∇ s ũ) -f . (2.24)
To evaluate the difference between the approximated wall-pressure p and the exact one p, we evaluate this equation on the tangential direction at the solid wall. From the noslip condition (satisfied both by the approximation and by the DNS), we have ũ = u = f = 0. We also have that the viscous term ∇ • (ν t ∇ s ũ) has to vanish at the wall, since ∇ • (ν t ∇ s ũ) = ν t ∇ • ∇ s ũ + ∇ν t • ∇ s ũ, the first term being zero since ν t = 0 at the walls and the second being also zero since, with the Spalart-Allmaras model ∇ν t = ∂ ν ν t ∇ν, which is zero because ∂ ν ν t | ν=0 = 0. Moreover, by imposing the divergence-free property of both the reference velocity field and the approximated one and projecting the result onto the wall-tangential direction t we have:

∂ t (p -p) -ν∆(ũ t -u t ) = ∂ t (p -p) -ν∂ nn (ũ t -u t ) = fx • t.
(2.25)

Since, at any iteration of the optimization procedure, the forcing term fx is a combination of adjoint variables that satisfy the no-slip boundary condition, we have fx = 0 at the wall.

If we apply the approximation ũ ≈ u (true at the end of successful optimization process), we have that the variation of the pressure difference along the wall is approximately null ∂ t (p -p) ≈ 0, justifying the good estimation of the wall-pressure. We remark that this approximation is valid only if, at the wall, the reconstructed tangential velocity field exhibits normal-derivatives that are close to the reference ones ν∂ nn ũt ≈ ν∂ nn u t . This condition is therefore only met if the quality of the reconstruction is very good near the wall, as is the case here with complete velocity measurements. Also, as mentioned before, this approximation for the pressure obtained at the wall does not hold anymore in the bulk of the flow (see figure 2.9). This is due to the fact that the evolution of fx is restricted to a solenoidal vector space. In order to study this error in the pressure field, following Foures et al. [2014], (2.24) shows that the pressure difference is related to the optimal forcing fx , the actual DNS forcing f and the eddy-viscosity term through:

∇(p -p) ≈ fx + ∇ • (ν t ∇ s ũ) -f . (2.26) (a) (b) (c) (d) (e)
Figure 2.9: Assimilated pressure (a) p (thick solid lines), reconstructed one pφ + φ (thin solid lines) and DNS one p (dashed lines). The divergence-free forcing terms fx + fs coming from the model (b) is also compared with the one coming from the DNS f s (c). The turbulent kinetic energy 2/3κ (computed from DNS statistics) (d) and the combination of the potentialsφ + φ (e).

The forcing term fx is solenoidal and the pressure error p -p appears as a scalar potential.

In order to identify in this equation solenoidal and scalar-potential parts, we perform a Helmholtz decomposition of the DNS forcing term f and of the eddy-viscosity term related to turbulence model:

f = f s + ∇φ (2.27) ∇ • (ν t ∇ s ũ) = fs + ∇ φ (2.28)
where f s and fs are the solenoidal functions and φ and φ are the scalar potentials. We therefore have: (2.30) where on the left-hand-side we have only terms stemming from the reference itself (thus the notation •) while those on the right-hand-side are linked to the optimal reconstruction (thus the notation •). If we plot the reconstructed pressure p + φ -φ, we can see that it matches very well with the actual pressure p (figure 2.9 (a)). The "divergence-free" component of the DNS forcing term f s compares now much better with the optimal one fx + fs (figure 2.9(b,c)). Also, in (2.7), we have seen that p = q+2/3κ, so that φ-φ ≈ 2/3κ if we assume that q = p, which is also verified (see figure 2.9 (d,e)).

p -φ ≈ p -φ (2.29) f s ≈ fx + fs ,

Sparse velocity measurements

Sparse velocity measurement operator

Now that we have performed and analysed the complete-velocity field measure, we turn our attention to a more realistic case where fewer (but still) velocity measurements are 4. Results

23

considered in the assimilation process. Instead of a 'spatial average' of the velocity field around the measurement location, we will consider a measure operator under the form:

M(ũ) = {ũ(x i m )} i=1,••• ,N .
(2.31)

Here x i m are N points where the velocity measurements are performed. Typically, these N measurements may correspond to cross-stream profiles of N y points at N x stream-wise stations, such that N = N x × N y . Experimentally, this would correspond a Pitot probe that is moved in the cross-stream direction at N x different stations, a rather common measure in experimental fluid dynamics. This way, the cost functional reads:

J = 1 2 N i=1 (ũ(x i m ) -m i ) 2 , (2.32)
where

m i = u(x i m )
. With this definition of the measurement operator, the right-hand-side of the adjoint equations becomes:

- ∂M ∂ ũ † (M(ũ) -m) = - N i=1 (ũ(x i m ) -m i )δ x i m (2.33)
where δ xm is the Dirac-mass centered at x m . With this right-hand-side defined, we can perform the data-assimilation for both models.

Reconstruction results using model 1 and model 2

Results of the reconstruction using model 1 and model 2 are presented in figure 2.10 for N x = 3 and N y = 10. From this figure, we can see that the volume-force model produces a solution which, although it respects the measurements (J n /J 0 ≈ 10 -6 ), does not correspond well to the reference u. We can see from figure 2.10 (a) that the isocontours of the assimilated velocity are not smooth and non physical. For example, the wall-friction shown in figure 2.10 (e) shows a (rather large) secondary re-circulation bubbles, a feature that does not exist neither in the RANS-SA solution nor in the DNS. We can explain these observations by noticing that the optimal forcing term fx is given by a linear combination of adjoint fields with right-hand-sides given by (2.33). Since those equations are excited on the momentum equation by a sum of Dirac masses, the velocityadjoint variables are 'peaked' around the measurements, which contaminates the solution, leading to those undesired oscillations. On the other hand, this does not occur in the fxdata-assimilation, since the ν-adjoint is not forced directly by the Dirac masses, leaving the gradients smooth. Furthermore, one can see that the result for model 2 is very similar to the complete velocity field assimilation. If we look at the wall-pressure or wall-friction plot (figure 2.10 (d,e)), we can see that the results in the sparse case present the same general features as in the complete case, including the (rather small) re-circulation bubble presented before. Model 2 therefore appears to be more robust than model 1 in the sense that, seemingly, for any type of velocity measurements (dense, sparse), the results of the reconstruction will be almost the same. This is confirmed when looking at the global errors: Table 2.2: Global error e Ω and convergence results for partial-velocity data-assimilations 2.2 we can also see that, in the cases where few measurements are provided (N y = 5, for example), the error made with model 2 is smaller than with model 1. This means that, whenever we have only few measurements, it is preferable to use the more robust model 2 than the more flexible model 1 as a means to compensate the lack of knowledge of the actual reference flow.

e Ω = 1 V Ω Ω (ũ -u) 2 dΩ, (2.34) (a) (b) (c) (d) (e)
e Ω J n /J 0 f x f x , γ 2 = 1 f ν f x f x , γ 2 = 1 f ν Total Measure ∼ 10 -3 - 0.035 ∼ 10 -4 - 12.3 % N y = 5 0.051 0.042 0.036 ∼ 10 -6 0.1 % 8.7 % N x = 3 N y =

Optimization with gradient penalization of model 1

However, in the case of denser velocity measurements, we may still prefer to use the more flexible model 1, since we have shown a more drastic reduction in e Ω for model 1 than for model 2. Yet, we need to improve the reconstruction associated to model 1 by trying to damp the oscillations seen in the solution (see figure 2.10 (a)). For this, we suggest to penalize the forcing term by modifying the cost functional J following:

J = J J 0 + γ 2 2 Ω |∇ f | 2 dΩ (2.35)
We remark that, by dividing the original cost functional by its value at the first iteration, we are normalizing the penalization with a unity value, making the penalization term independent of the measurement and only dependent on the parameter γ 2 . This parameter should not be taken too small (the penalization needs to be effective) and not too big (the optimization should not get stuck in the vicinity of the RANS-SA solution). It tells the algorithm to favor a solution that is smoother, since the tuning term f in the momentumequations will exhibit smaller gradient values. This is a physically sound constraint since we know that the Reynolds-stress forcing is very smooth. The results for this modified data-assimilation procedure with model 1 are shown in figure 2.10 (c) for γ 2 = 1. We can see that the resulting solution is smoother and still matches well with the measurements (J n /J 0 1%), suggesting that the kind of penalisation employed (here derivatives of fx ) still allows enough freedom for the tuning field to match the measurements, while constraining the solution in a smooth subspace. Furthermore, except for the very sparse setup (N x = 3, N y = 5), the penalized algorithm even manages to give a better reconstruction of the flow away from the measurement points (see table 2.2). Indeed, for all such test-cases, the penalized model 1 solution always manages to outperform both the non-penalized model 1 reconstruction and the model 2 reconstruction.

Understanding the 'rigidity' of model 2 through observability Gramian analysis

This section is devoted to the understanding of the rigidity of model 2, pointed out in the last section. One possible way to understand this observation is to examine the linearized optimization problem, say, around the RANS-SA solution. In this linear framework, a variation of the control parameter δ f (here, either fx or fν ) induces a variation of the state variable q = (ũ, p, ν):

δ q = ∂R ∂ q -1 P δ f = (∂ qR) -1 P δ f (2.36)
where R(q) is the nonlinear residual composing the RANS-SA equations and P being an operator that transforms the variation of the control parameter into the actual forcing term of the Jacobian. This variation in the state induces a variation on the observable data according to:

δ m = ∂M ∂ q δ q = (∂ qM)δ q = (∂ qM)(∂ qR) -1 P δ f = Aδ f (2.37)
To be able to identify which forcing terms induce the most energetic variation on the measure we may optimize the gain:

G(δ f , δ m) = ||δ m|| 2 M ||δ f || 2 = δ m, δ m M (δ f , δ f ) Ω , (2.38)
which is equivalent to solving the Singular-Value decomposition of the operator A = (∂ qM)(∂ qR) -1 P :

A † Aδ fi = λ 2 i δ fi , (2.39)
where λ 2 i are the ranked squared singular-values (λ 2 i ≥ λ 2 i+1 ). The values λ i represent the measurement variations along the unit optimal measurement directions δ mi (=λ -1 i Aδ fi ) associated to the unit optimal forcing directions δ fi . If all those singular-values are of the same order, we can say that any measurement δ m is equally reachable (or, more suitably, 'observable'). If we have a strong separation of singular-values, that is λ 2 0 λ 2 i for some index i, we have measurement states that cannot be reached ('observed') with the chosen model. This can be seen by considering that, for a given δ f , the resulting perturbation on the measure is given by:

δ m = Aδ f = A i α i δ fi = i α i λ i δ mi = λ 0 i α i λ i λ 0 δ mi . (2.40)
Since δ mi is unit norm, this shows that the amplitude α i required for the forcing term to reach measurement δ mi scales as λ 0 /λ i , which can be very large, and therefore not achievable. In figure 2.11, we have plot the quantities (λ i /λ 0 ) 2 for model 1 and model 2, taking as measurement operator the full velocity field M(q) = ũ. We can see that the fx model 1 is much more flexible than the fν model 2 so that more measurement states δ m can be reached. As for the fν model 2, we can see that the separation is very strong and increases very rapidly for i > 1.

Conclusion

In this chapter, we have introduced two models for the reconstruction of mean-flow features from mean-velocity measurements. The first one consists in the perturbation of the momentum equations by a volume-force in a similar fashion as in [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF], the difference being that we have a background eddy-viscosity turbulence model, making the numerical procedure more robust and well-posed for high Reynolds number flows. The second one consists in the correction of the equation governing the eddy-viscosity (here the Spalart-Allmaras model). The physical configuration on which those models are tested is a smooth backward-facing step. We can show that, if the whole velocity field is considered as measured (like in a PIV setup), the fx model 1 produces a solution that matches exactly the reference, showing the high flexibility of the model. Moreover, we could also show that the optimal forcing parameter matches almost perfectly the projection of the real forcing term onto the space of solenoidal functions, in a similar way as [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF]. As for the fν model 2, we showed that, if a complete velocity field measurement is available, it is not capable to reproduce exactly the reference state. The model is too 'rigid' and lots of measurements are not accessible with this model. If now only sparse point-wise velocity measurements are available, we show that the first model, despite its high flexibility to exactly recover the measurements, may lead to noisy / unphysical state reconstructions. For this reason, an additional penalisation on the gradients of the tuning field needed to be considered to obtain smooth solutions. As for the second model, it turned out to be more rigid in the sense that, independently of the quantity of available measurements, the optimal state was approximately the same. This behaviour could be analysed by performing a singular value decomposition of the linear operator between the forcing space and the measurement space. We showed that the controlability/observability of the second model was actually restricted to a rather small subspace. However, such a rigidity may be helpful whenever the number of measurements is very low since its optimal solution is almost independent of the number of measurements. Hence, model 2 can provide in such cases a better reconstruction of the flow field than the much more flexible model 1.

Chapter 3

Wall-measure assimilation for RANS equations 1 Introduction

In the second chapter, we dealt with the case where the information or data to be assimilated was considered to be sparse velocity-field information. The measure of this kind of quantity from experiments can be easy to do if not many probes are used, employing, for example, hot-wire or Pitot probes. However if one needs a great number of velocity measurements, more refined experimental techniques may be needed (for example, Particle Image Velocimetry, [START_REF] Raffel | Particle image velocimetry: a practical guide[END_REF]), the drawback in those cases being the requirement of more sophisticated experimental apparatus. Instead, in this chapter, we evaluate other quantities that are more easily obtained experimentally, for example, the wall-pressure and friction. We remark that those measures are more indirect, regarding the velocity field itself, and their assimilation's performance needs to be investigated. It is interesting to note that, from a practical point of view, the assimilation of, say, the wall-friction, can be useful even in cases where sufficiently bulk velocity field information is available since it provides an easy way to express the localized behavior of the flow near the wall, making the optimization procedure to impose its correct behavior there.

From a mathematical point of view, in an incompressible framework, the derivation of the gradient for this cost functional is less straightforward than the velocity measure one. The reasons for that are related with a seemingly over-constrained set of boundary conditions for the adjoint problem. Since this issue is not so much discussed in the literature (even if similar physical problems have already been addressed, with the correct adjoint equations, see [START_REF] Bewley | Skin friction and pressure: the "footprints" of turbulence[END_REF]), we will dedicate the next section to the exposure of the mathematical problem and to the proper derivation of the adjoint equations on a laminar case, together with a few numerical results on the transitional flow around a circular cylinder (as in [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF]), where the control parameter will be a volume-force acting on momentum equations. Then, we will apply this procedure in the case of the turbulent Backward-Facing Step (BFS), already introduced in the previous paragraph, using for this purpose the two models discussed there. As we will see, this case is much more challenging since, depending on the model used and on the quantity we measure, the model will provide a solution that corresponds more or less to the reference. For instance, if the wall-friction information for assimilation using model 2 ( fν ), the optimal solution will produce a re-circulation bubble that is much thinner than the actual 29 Chapter 3. Wall-measure assimilation for RANS equations one, over-estimating the pressure. This chapter can thus be viewed as an extension of the first, since the overall goals and methods will essentially be the same.

2 Wall-measure for laminar RANS equations

Theoretical developments on laminar equations

In this section, we introduce the theoretical aspects of the adjoint derivation for the wallinformation measure operator. For the sake of clarity, we choose to place ourselves first in the case where the flow is governed by the laminar RANS equations. We remark that all technicalities related to the adjoint boundary conditions, that will allow us to assimilate wall-information are present in the same manner if we consider a turbulence model. The exact RANS equations, starting point of the analysis, are given by:

u • ∇u + ∇p -ν∆u = -∇ • u ⊗ u = f , ∇ • u = 0, (3.1)
which differ from the steady solution equation by the presence of the right-hand side term f , coming from the nonlinear interactions of the fluctuation u . Following Foures et al.

[2014], [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF], we will use this right-hand side term as tuning parameter for our assimilation procedure. Its goal will be then to find the forcing term f , satisfying:

ũ • ∇ũ + ∇p -ν∆ũ = f , ∇ • ũ = 0, (3.2)
such that the following cost functional:

J Γ (ũ, p) = 1 2 ||M Γ (∇ũ, p) -m Γ || 2 M Γ (3.3)
is minimal. M Γ is now a (rather general) measure operator that takes as input the information of the gradient of velocity field and pressure at (a portion of) the wall Γ and gives the measure output on the measure space M Γ :

M Γ : (∇ũ, p)| Γ → M Γ (∇ũ, p) ∈ M Γ (3.4)
and m Γ = M Γ (∇u, p) represents the measure applied on the actual reference solution.

We note that, since we suppose the measure operator acts on a solid wall Γ, on which the no-slip boundary conditions are assumed, it does not make sense to measure the velocity field. However, in order to measure the wall-friction, its derivatives are needed. The measure space is equipped with a norm || • || M Γ that will quantify the distance from the measure of our model and the actual target data to be assimilated m Γ . Just as before, in order to derive the gradient, we consider the Lagrangian:

L = J Γ + (ũ † , ũ • ∇ũ + ∇p -ν∆ũ -f ) Ω + (p † , ∇ • ũ) Ω (3.5)
The variations of this Lagrangian with respect to the adjoint variables will, again, lead to the set of direct equations (3.2). Its variations with respect to the state (ũ, p) will lead to the adjoint equations, which will be slightly different from those found in the previous chapter. The reason for this, will be the nature of the measure operator, that takes as input boundary terms. As we will see, this fact will change the forcing term of the adjoint equations from a volume force to a non-homogeneous boundary condition. To properly derive them, we will proceed with the formal derivation of the Lagrangian:

∂L ∂(ũ, p) , δ(ũ, p) Ω = M Γ (∇ũ, p) -m Γ , ∂M Γ ∂(∇ũ, p) • (∇δ ũ, δ p) M Γ + ũ † , ũ • ∇δ ũ + δ ũ • ∇ũ + ∇δ p -ν∆δ ũ Ω + p † , ∇ • δ ũ Ω = 0 (3.6)
We remark that, for any linear operator that acts on the boundary such as:

A Γ : u Γ → A Γ u ∈ M Γ (3.7)
we can define its adjoint as:

(m, A Γ u) M Γ = A † Γ m, u Γ (3.8) where (v 1 , v 2 ) Γ = Γ v 1 • v 2 dΓ.
With this in mind, and performing the appropriated integration by parts, we can derive the adjoint equations to be:

ũ † • (∇ũ) T -ũ • ∇ũ † -∇p † -ν∆ũ † = 0, ∇ • ũ † = 0 (3.9)
We remark that, since the measure operator acts now on the boundary of the domain, the adjoint equations are no longer forced by a volume term. Those equations must however be completed with a set of boundary conditions. For the inflow Γ I , outflow Γ O and slip Γ S boundaries the boundary conditions are exactly the ones obtained by [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF]. However, since now we are performing some measurements on the solid wall Γ, the classical no-slip boundary condition does not hold anymore. To re-establish those conditions there, we have to impose that all boundary terms, coming from integration by parts, vanish (having in mind that δ ũ = 0 at Γ). Those terms can be written as:

  -νn ⊗ ũ † + ∂M Γ ∂∇ũ † (M Γ (∇ũ, p) -m Γ ), ∇δ ũ  Γ +   n • ũ † + ∂M Γ ∂ p † (M Γ (∇ũ, p) -m Γ ), δp   Γ = 0, ∀(δ ũ, δ p) (3.10)
At a first glance, this seems impossible since, if we take ∇δ ũ = 0, we see that we need to impose that the normal adjoint velocity has to be driven by the measure operator on the pressure:

n • ũ † = - ∂M Γ ∂ p † (M Γ (∇ũ, p) -m Γ ) (3.11)
Inversely, if we take δ p = 0, we have that:

νn ⊗ ũ † = ∂M Γ ∂∇ũ † (M Γ (∇ũ, p) -m Γ ) (3.12)
which not only demands the imposition of all the adjoint velocity on the wall but also imposes some restrictions on the measure operator itself, since the above equation is a tensor-equation. We can see then that, instead of having a tensor-like equation, we should

Chapter 3. Wall-measure assimilation for RANS equations have a scalar one, such that, together with 3.11, a complete set of boundary conditions for ũ † is given. To try to achieve that, we will simplify the first term of equation 3.10 using the fact that, along the wall, the test function must not change (in fact it is identically zero). This condition can be written as t • ∇δ ũ = 0, where t is any tangent vector (t

• n = 0).
In what follows, the choice for t will be fixed to be t = (-n y , n x ) (in a 2D configuration). With this arbitrary choice for the tangent direction, we can rewrite this condition as:

t • ∇δũ x = t x ∂ x δũ x + t y ∂ y δũ x = -n y ∂ x δũ x + n x ∂ y δũ x = 0 t • ∇δũ y = t x ∂ x δũ y + t y ∂ y δũ y = -n y ∂ x δũ y + n x ∂ y δũ y = 0
For a general point on the solid wall (for which n x n y = 0), we can solve this set of equations for ∂ x δũ y and ∂ y δũ x . Those extra scalar equations establish two links between the components of the tensor ∇δ ũ, reducing the possible variations of it:

∇δ ũ = ∂ x δũ x ∂ x δũ y ∂ y δũ x ∂ y δũ y = ∂ x δũ x n x /n y ∂ y δũ y n y /n x ∂ x δũ x ∂ y δũ y (3.13)
Using those conditions on equation (3.10), we now can re-derive a new equation for the boundary conditions for ũ † , taking again δ p = 0. Since we use now a supplementary vector-like equation (t • ∇δ ũ = 0), we manage to reduce the previous equation from a tensor-like one to a vector-like one. However, we still observe contradiction with respect to the boundary condition (3.11). For this reason, at this point, we also suppose that ∇ • δ ũ = 0 at Γ. This hypothesis may seem arbitrary, since the test function δ ũ does not need, in general, to satisfy one of the volume-differential equations (continuity) of our model, specially at the boundary Γ. However, one may make this argument rigorous by changing the space of our test functions such that it verifies ∇ • δ ũ = 0 at Γ. With this extra condition, the quantity ∇δ ũ becomes:

∇δ ũ = 1 -n x /n y n y /n x -1 ∂ x δũ x = D∂ x δũ x (3.14)
Replacing this assumption on the first term of equation 3.10, we have:

  -νn ⊗ ũ † : D + ∂M Γ ∂∇ũ : D † (M Γ (∇ũ, p) -m Γ ), ∂ x δũ x   Γ (3.15)
Where the notation A : B stands for i,j A i,j B i,j . We can see that we now have a scalar equation that will lead to only one extra boundary condition for the adjoint, to be imposed with (3.11). If we simplify the term:

-νn ⊗ ũ † : D = -ν∂ x δu x ũ † x n x - ũ † y n y = ν(n x n y ) -1 ∂ x δu x t • u † (3.16)
we notice that, for the boundary-integral term to vanish, we need to impose the tangential velocity from the measurements on the viscous friction together with the normal velocity from the pressure at the wall: With those boundary conditions, equations (2.1) are well posed. Lastly, by taking the variation of the Lagrangian with respect to the forcing term (in a similar way as before), we obtain the gradient ∇ f J Γ = -ũ † . With those tools, we are ready to perform wallmeasure data-assimilation. We will first apply this technique on a laminar/transitional flow, described in the next paragraph.

t • ũ † = - 2(n x n y ) ν ∂M Γ ∂∇ũ : D † (M Γ (∇u, p) -m Γ ) n • ũ † = - ∂M Γ ∂ p † (M Γ (∇ũ, p) -m Γ ) (3.17) (a) (b) (c) (d)

Physical Description of the Considered Laminar Flow

The physical configuration that will be treated here is the same as the one in [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF], the flow around a cylinder at Re = 150. The corresponding steady solution (figure 3.1 (a)) is unstable and the time-dependent solution deviates from it, saturating into a nonlinear limit cycle, whose time-average u x can be depicted in figure 3.1 (b), together with the corresponding divergence of the Reynolds-stress tensor f (figure 3.1 (c,d)). The re-circulation bubble's length is much shorter for the mean-flow velocity field (L DN S = 1.55, when L Steady = 9.75). The goal of the optimizations described on this section will be to get as close as possible to those fields. The reference mean-flow state was obtained through a two-dimensional Direct Numerical Simulation of equations (5.1). This simulation was performed in the software FreeFem++, using P 1 b , P 1 b , P 1 discretization (see [START_REF] Hecht | New development in freefem++[END_REF]) on a mesh of around 70000 triangles with local refinement on the cylinder's surface and on the wake region. This simulation will also provide the necessary data to perform the assimilation in the following sections.

Results

In this section we will evaluate how our optimization procedure works for wall-pressure and wall-friction information assimilation.

Wall-Pressure Assimilation

To account for the wall-pressure, the generic measure operator introduced before needs to be particularized. Since here we do not perform any measure on the velocity-gradientrelated quantities, this operator takes the form:

M Γ (∇ũ, p) = p (3.18)
We remark that, by this choice of measure-operator, we implicitly suppose that the information of wall-pressure of the reference state (DNS presented in the last section) is given on the same computational mesh as the solution. This hypothesis may hold in practice when several pressure probes are available on the surface of the cylinder. However, here, for the laminar case, since our goal is to validate the assimilation procedure, we will not be interested in coarsening the spatial distribution of the measure. For this reason, the measure space M Γ will be the set of all the functions with support on the boundary. A natural norm for this space is given by || • || Γ = ( Γ (•) 2 dΓ) 1/2 . We can thus infer that:

∂M Γ ∂∇ũ = 0, ∂M Γ ∂ p = I ⇒ ∂M Γ ∂ p † = I (3.19)
And, consequently, we have that, in a simplified manner:

t • ũ † = 0, n • ũ † = -(M Γ (∇ũ, p) -m Γ ) = -(p -p) (3.20)
With those considerations, the optimization procedure leads to the results presented in figure 3.2 and 3.3. We can see that the wall-pressure reconstruction manages to exactly capture the actual reference pressure coming from the DNS. Moreover, it does so in a very small number of iterations (Figure 3.2), such that the cost functional ((J Γ ) n /(J Γ ) 0 ) converges in a nearexponential manner. We can safely say then that the data-assimilation procedure, based on the theory presented in the last section, was successful for the pressure measurement. Furthermore, we can see that, by measuring the pressure, one also reduces the error in the wall-friction (3.2). Looking now at the reconstructed velocity field (3.3 (a)), we can see that, although the re-circulation bubble is now shorter (L ≈ 5.15), it is still far away from the actual one from the DNS (L DN S = 1.55), showing that wall-pressure information may be insufficient for reconstruction of the reference state. Furthermore, we can see that the spatial structure of the forcing term (3.3 (b)) does not compare with the DNS one (Figure 3.1 (c,d)) and is localized around the cylinder due to the inhomogeneity of the adjoint's boundary condition.

Wall-Friction Assimilation

We turn now our attention to the wall-friction measure. In this case, the measure operator will be given by:

M Γ (∇ũ, p) = νt • ∇ s ũ • n = ν(n ⊗ t + t ⊗ n) : ∇ũ (3.21)
As before, the norm of the measure space is given by || • || Γ = ( Γ (•) 2 dΓ) 1/2 , leading to:

∂M Γ ∂ p = 0, ∂M Γ ∂∇ũ = ν(n ⊗ t + t ⊗ n) ⇒ ∂M Γ ∂∇ũ : D † = ν(n ⊗ t + t ⊗ n) : D (3.22)
The boundary conditions for the adjoint in this case are:

t • ũ † = -(M Γ (∇ũ, p) -m Γ ), n • ũ † = 0 (3.23)
The results for this optimization, presented in figures 3.4 and 3.5 are close to the pressure-measurements one. First of all, as before, we can see that the convergence of the cost functional is very fast and the reconstructed wall-friction matches perfectly with the one coming from the DNS. The wall-pressure distribution, although not considered in the cost functional, is also better represented. Furthermore, the reconstructed velocity field presents a re-circulation bubble of L ≈ 5.20, which is very similar with the one previously shown. Moreover, the forcing term, although slightly different from the one obtained for the pressure measurement, is still concentrated near the wall, also a symptom of the inhomogeneous boundary conditions for the adjoint variable, now on the tangential velocity field. 3 Wall-measure for turbulent RANS

Once the data-assimilation is now validated, we proceed to the turbulent case. Although now the equations are more complex than in the laminar case, due to the consideration of the Spalart-Allmaras model, the volume-adjoint equations remain the same as the ones presented in the previous chapter, with a null volume-forcing term:

ũ † • (∇ũ) T -ũ • ∇ũ † -∇ • ((ν + ν t )∇ s ũ † ) -∇p † (3.24) +ν † ∇ν + ∇ • (ν † ∂ ∇ũ s) = 0 (3.25) ∇ • ũ † = 0 (3.26) -ũ • ∇ν † -∇ • (η∇ν † ) + (∂ ν η)∇ν † • ∇ν + (∂ ν ν t )∇ũ † : ∇ s ũ -(∂ ν s)ν † + ∇ • (ν † ∂ ∇ν s) = 0, (3.27)
The boundary conditions for those equations on the solid wall are derived in a similar manner as before. We note that, since the variation of the turbulent variable at the wall is strictly null in the Lagrangian formalism, the adjoint boundary equation does not differ from those presented for the laminar case (equation (3.17)), the difference being the imposition of a homogeneous Dirichlet boundary condition for the variable ν † = 0 on the solid wall. We remark that, for both models (model 1 and 2), the adjoint equations are the same. Furthermore, the gradient expressions, in function of the model, are the same ones given in the previous chapter.

Results for turbulent backward-facing step

We now apply those data-assimilation techniques on the turbulent Backward-Facing Step, already presented in the last chapter. As we will see, contrarily to the laminar case where the final optimal state seems to be almost independent on the measurement of wallpressure or friction, the final solution now can greatly depend on the measure performed and on the model (1 or 2). We will first present the results of the wall-friction measurement in detail, that seems to be the most difficult quantity to be assimilated. Then we will present the wall-pressure assimilation that, in a way, resembles the partial bulk-velocity assimilation, presented in the previous chapter.

Wall-Friction Information

We turn our attention now to the wall-friction measure. As before, we will consider a complete knowledge of the wall-friction at every point of our computational mesh, limited 3.8 (a,b,c)) are distant from the DNS, specially the one at x = 4.5, where, in order to match the wall-friction, the velocity profile creates an inflection point and a negative-velocity zone near the wall, leaving the impression that the re-circulation bubble is actually bigger than it is. This problem is due to the particular form of the forcing term fx , which is non null at the wall and very localized there (see figure 3.8 (d)). We note that this characteristics of the forcing term is not representative of what happens with the reference solution, since the actual term f = -∇ • τ is null at the wall. This may be at the origin of this anomalous behavior of the wall-friction data-assimilation using the model 1. We note however that, even with this theoretical inconsistency on fx , one can still produce a solution that is less distorted. For example, in the laminar case presented before, where the Reynolds number was much lower, the velocity field did not seem to produce such behavior. This is consistent with the fact that the viscosity was much higher, allowing the forcing term and its corresponding solution to be smoother. To simulate those conditions in this high-Reynolds number flow, we propose the penalization of the gradients of the forcing term on the cost functional, in the same way presented in the second chapter. We make clear, however, that this proposition does not at all solve the actual theoretical inconsistency problem discussed, Chapter 3. Wall-measure assimilation for RANS equations as the inhomogeneous boundary conditions will still exist. It only attenuates one symptom of it at high Reynolds number flows. We can see that this strategy manages to correct the tendency of the velocity profiles (3.8), while keeping a good representation of the wall-friction. However, in terms of global velocity-error, the error increases slightly (see table 3.1). From this table, we can also see that neither of the proposed models managed to recover well the actual velocity field from the wall-friction. For this reason, we will set aside this measure and focus on the wall-pressure in the next section.

Ω (J Γ ) n /(J Γ ) 0 f x f x , γ 2 = 1 f ν f x f x , γ 2 =

Wall-Pressure Information

We consider now the complete information of the pressure distribution along the wall as our measure, restricted to the interval x ∈ (0, 10). The main results are presented in figures 3.9 and 3.10. By performing such data-assimilation for the two models, we can see that, in terms of the reduction of the cost-functional, the model 1 manages to converge precisely to the exact pressure, whereas the second one, although it produces a very good agreement in most of the wall-region, fails to capture the exact pressure at x ∈ (0.5, 1). It is as well at this region where the corresponding wall-friction (figure 3.9 (d)) oscillates the most, provoking a (rather small) secondary bubble. This is coherent with the relation between the pressure variation along the wall and the concavity of the wall-normal velocity-error (deduced in the second chapter) ∂ t (p -p) = ν∂ nn (ũ tu t ), allowing an (adverse) error in wall-pressure to gradually modify the inflection of the normal velocity profile, causing the separation. For the model 1, a (small) secondary bubble also occurs around x ≈ 1.9. We can also see that both models make similar reconstructions for the re-circulation bubble It is interesting to note as well that, for the model 1, the fact that fx • n = 0 at the wall did not seem to affect the solution in the same way it did for the wall-friction measure. Indeed, the u y velocity profiles obtained from data-assimilation of both models (3.11) are close to each other and do not present the "anomalous" behavior as in the case of wallfriction measure. The explanation for this fact is given by the Helmholtz decomposition of the forcing term fx , similarly to those used in the second chapter. Indeed, we have:

fx = f s x + ∇φ (3.28)
where f s x is solenoidal, satisfying f s x • n = 0 at the solid wall Γ and φ is a scalar potential. We note that, although the forcing term fx is solenoidal (in the bulk), since it holds nonnull value at Γ, it can produce a scalar potential different from zero φ = 0. Indeed, by taking the divergence of this equation, we have that φ is given by a Laplace equation -∆φ = 0 in Ω and n • ∇φ = fx • n = 0 at the wall. This fully determines φ. Since φ is a scalar potential, its contribution, which holds the non-null component of fx • n, can be merged into the pressure, leaving unchanged the actual forcing acting on the velocity field alone. For this reason, the fact that fx • n = 0 does not seem to be an important issue. It is important to note that the present argument does not apply for the previous case theoretical aspects of the wall-pressure measure, we now proceed to the sparse-measure of the pressure. Here, we suppose that the we dispose of N w measurements m i . Those measurements are supposed to be localized at some discrete locations on the wall, given by the points x i . For numerical simplicity of the imposition of the adjoint boundary conditions, we will assume that those measures will correspond to the following quantity:

mi = Γ p(x)g x i dΓ ≈ Γ p(x)δ x i dΓ = p(x i ) (3.29)
where g x i is a normalized Gaussian function centered at x i . The typical standard-deviation of those Gaussian functions are of a few mesh-elements. The boundary conditions of the adjoint problem are given by:

ũ † • n = - Nw i=1 g x i ( mi -m i ), ũ † • t = 0 (3.30)
We also suppose that those measure points x i will be evenly distributed in the interval x ∈ (0, 10). The results of this data-assimilation are summarized in figure 3.12, where the wall-pressure and wall-friction are plotted, together with a vertical gray line, indicating the measure locations. We can see that the pressure is indeed satisfied at the measure locations for both models 1 and 2. We can also see that both models present a reluctance to reconstruct well the pressure around x ∈ (0.5, 1). This was already observed by the dense wall-pressure measure for the model 2. Thus, this seems to be a difficult region for the pressure assimilation, probably because of the presence of the detachment point. We can also see from figure 3.12 that the resulting wall-friction is globally similar to the one obtained by the full information of the pressure. However, it is interesting to note that, due to the sharp/localized nature of the gradient, given by the measure, the forcing term will also be sharp and may lead to some oscillation on the wall-friction, especially for the model 1, where the forcing term acts more directly on the momentum equations. This issue can be attenuated, again, with the aid of the penalization term of the gradients of fx . This leads to smoother wall-frictions, but not necessarily to better global velocityfield reconstructions (see table 3.2). From the point of view of the global velocity fields reconstruction, we can see that the model 2 produces systematically a better solution than the first model. We can see that, although the model 1 is more adaptable to fit data (as presented throughout the second chapter), the information of the pressure is not sufficient for it to reconstruct a correct solution. On the other hand, since the second model is more rigid than the first one, we can, again, take advantage of this to produce better solutions, for example, with N w ≥ 10. We can see as well that, for the model 2, the value of the error is approximately the value obtained with velocity data-assimilation, presented in the second chapter.

Ω (J Γ ) n /(J Γ ) 0 f x f x , γ 2 = 1 f ν f x f x , γ 2 = 1 f ν
Total Measure 0.057 -0.037 ∼ 10 -5 -6%

N w = 5 0.061 0.061 0.045 ∼ 10 -5 ∼ 10 -5 ∼ 10 -5 N w = 10 0.057 0.057 0.038 ∼ 10 -5 ∼ 10 -4 0.7% N w = 20 0.057 0.057 0.036 ∼ 10 -5 0.2% 0.2% RANS-SA 0.094 -

Combined Measure: Wall-Pressure and Velocity

In the last paragraph, we investigated how the data-assimilation procedure performs under sparse wall-pressure assimilation. We could see that, for the models 2, the final assimilated solution was almost independent on the number of measure points. This solution was also very similar to the one obtained with velocity data-assimilation. For this reason, the combined measure will not be performed with this model. However, since the model 1 is fully controlable, it may be interesting to understand the compromise between the pressure or velocity assimilation. For this purpose, we perform two assimilations, one with minimal information (N x = 3, N y = 5 and N w = 5) and other with maximum information (N x = 6, N y = 20 and N w = 20). We can see from table 3.3 that the data assimilation results do not change much with the addition of the wall-pressure measure. This suggests that, from the point of view of the global error e Ω , the wall-pressure does not include much extra information with respect to the velocity measurement. If, however, wallpressure probes are preferred instead of bulk velocity ones, we may look for minimizing the amount of the last ones. For this purpose, we could look for optimally place those velocity probes on top of the wall-pressure ones such that the information provided by them is not redundant (see, for example, Mons et al. [2017a]), allowing us to explore directions not observed by wall-pressure probes. This will be considered however out of the scope of the present work. Pressure/Velocity Velocity alone Pressure alone

N x = 3, N y = 5 N w = 5 0.041 0.042 0.061 N x = 6, N y = 20 N w = 20
0.027 0.027 0.057

Conclusions

We could see from the presented wall-measure data assimilation that indeed this information is not in general sufficient to reconstruct correctly the velocity field. In particular, the wall-friction assimilation produced some unexpected results. For example, the second model produced a much thinner re-circulation bubble, over-estimating the wall-pressure.

Since the control parameter is completely associated with the turbulence model, we may expect that, by considering other models, we could obtain more appropriated results. As for the model 1, we could notice strong gradients of the forcing term near the wall, leading to a solution that oscillates there, creating a seemingly larger re-circulation bubble. This is shown to be related to some inconsistencies related to a inhomogeneous wall-boundary condition for the adjoint, leading to non-null values for forcing term, which is not the case of the actual forcing term coming from the Reynolds-stress tensor. A possible way to attenuate this effect is to penalize the gradients of the forcing term, smoothing it and consequently the final solution as well.

As for the wall-pressure assimilation, the results obtained are much more coherent with those obtained with the velocity-field, presented on chapter 2. The model 1 does not manage to reconstruct well the velocity field, and we conclude that the wall-pressure is not sufficient information for the capture of the velocity field. The second model, on the other hand, even with few pressure probes, manages to reconstruct the velocity field in the same manner as with the velocity-field data-assimilation. This issue was already explained in the second chapter, where the rigidity of this model was exposed and its use with very few measure was proposed.

Chapter 4

Resolvent-based Data-Assimilation 1 Introduction

The reconstruction of a time-averaged (mean) flow together with the fluctuation around it is a topical issue within data-assimilation in fluid mechanics. In the previous chapters, we explored the mean-flow data-assimilation only through a variational minimization algorithm, where some tuning parameter of the considered model for the mean flow is adjusted so that the resulting approximation of the mean-flow is as close as possible to some input data, coming from experiments of higher-fidelity numerical simulations. However, one common drawback of this approach is the fact that, in order for the final assimilated solution to represent well the reference one, a sufficiently large number of probes is needed. This issue is intimately related to the large number of tuning parameters considered then since, if few measurements are provided, the gradient of those spatiallydistributed tuning parameters may be sharp and may demand regularization (such as the gradient penalization employed with sparse velocity measurements, present in the second chapter).

Although those techniques can lead to correct mean-flow reconstructions, they remain purely steady approaches, where mean-flow information is used to reconstruct the meanflow itself. In the case where one desires to obtain a time-resolved assimilated field, one has to resort to unsteady data-assimilation techniques. Classically, those methods belong to two (possibly overlapping) categories: 3/4DVar (Le [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspect[END_REF]) and Ensemble Kalman filters [START_REF] Evensen | Data-Assimilation: yhe ensemble Kalman filter[END_REF]). Those methods have been extensively applied in the context of meteorology [START_REF] Lorenc | Analysis methods for numerical weather prediction[END_REF], [START_REF] Liu | An ensemble-based four-dimensional variational data-assimilation scheme. Part i: technical formulation and preliminary test[END_REF]) and more recently in fluid mechanics [START_REF] Mons | Reconstruction of unsteady viscous flows using data assimilation schemes[END_REF], Gronskis A. [2013], D'Adamo et al. [2007]), [START_REF] Meldi | Augmented prediction of turbulent flows via sequential estimators[END_REF]. This kind of data-assimilation requires some obvious extra numerical elements that were not needed before, such as a time-stepper code, and possibly a time-adjoint code (for 4DVar) and/or a statistical treatment of an ensemble of points in the phase space (for Ensemble Kalman filter).

Another approach proposed in this chapter, more suited for statistically-steady flows, and which does not require many extra numerical tools, relies on the analysis of the linear Navier-Stokes operator around the mean-flow. The origin of such analysis lies in the work of [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF], where he observed that, in the case of the supercritical flow around a circular cylinder, the stability analysis of the linearized Navier-Stokes operator on the time-averaged (mean) flow (instead of the steady solution, or the base-flow) gives rise to a mode whose eigenvalue is purely imaginary and its frequency matches well with the one coming from the nonlinear limit cycle. The theoretical explanation of this fact was further explored by Sipp and Lebedev [2007a] with the aid of weakly-nonlinear analysis near the bifurcation. They also noticed that the mode coming from this analysis should approximate the first harmonic of the nonlinear limit-cycle. Later, this lead [START_REF] Mantic-Lugo | A self-consistent model for the saturation dynamics of the vortex shedding around the mean flow in the unstable cylinder wake[END_REF] to develop a self-consistent model describing the nonlinear saturation of such flow, where the steady solution evolves to the mean-flow through the action of the Reynolds-stress tensor, that can be reconstructed with the mode coming from the stability analysis. All the previously mentioned references concerned oscillator flows (see [START_REF] Chomaz | Global instabilities in spatially developing flows: Non-normality and nonlinearity[END_REF]). The extension of the mean-flow analysis to flows whose spectral content is arbitrarily spread in frequency has been conducted by [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF], who showed that, under some fairly general assumptions, the Fourier mode of the nonlinear fluctuation around the mean-flow at a given frequency is proportional to the leading mode of the Singular Value Decomposition of the Resolvent operator around the meanflow, evaluated at that frequency. This means that with the knowledge of the mean-flow one can estimate the spatial structure of the fluctuation, at a given frequency, up to a (complex) multiplicative constant, modeling the energy and the phase of that mode. This constant could be estimated with some extra (sparse) information of the fluctuation field, typically some pointwise time-resolved signal. This procedure was validated for turbulent flow over a backward-facing step, in which case the mean-flow was obtained through a time-resolved simulation (ZDES, see [START_REF] Deck | Zonal-detached-eddy simulation of the flow around a high-lift configuration[END_REF]). The main drawback of this approach is of course the need for the a priori knowledge of the mean-flow. To overcome this issue, we propose here to combine the estimation of the fluctuation field from the Resolvent analysis in order to reconstruct the Reynolds-stress tensor and the mean-flow equation involving this Reynolds-stress tensor. By doing so, we obtain a nonlinear model whose unknown parameters are the amplitudes of the Resolvent modes at each considered frequency. For this reason, one could consider this model as an a priori reduced order model for the unsteady flow, which is then used for data-assimilation. In addition, only very sparse information is needed for the reconstruction of not only the mean-flow but also the fluctuation around it.

In this chapter, this procedure will be applied to simple model problems of transitional laminar flows. The first configuration will be the flow around a square cylinder (Re = 60, 100), falling in the category of oscillator flows, for which the frequency spectrum exhibits an important peak on the vortex-shedding frequency and smaller ones on its honlinear interactions. This means that only one Resolvent mode is needed at one frequency (and thus only one parameter) for the reconstruction procedure. We will show that the velocity at some point in the flow then suffices for the estimation of the mean-flow and the fluctuation field. The second configuration is the two-dimensional backward-facing step (Re = 500), a typical example of noise amplifier [START_REF] Hervé | A physics-based approach to flow control using system identication[END_REF]. For this flow, the spectral content of the fluctuation will strongly depend on the external noise that triggers the dynamics. We will consider two different cases, one for which the noise is characterized by two incommensurable frequencies and the other for which it is a broadband white noise. The frequency content of the flow, which can be inferred from time-resolved signals at some points of the flow, can then be used as an input signal for the data assimilation process or as a reference signal to assess the quality of the reconstruction.

This chapter proceeds as follows. First ( §2), we present the model combining the mean-flow and the perturbation equations. In particular, we show how to solve for this problem for a given set of Resolvent amplitudes. Then, we successively present the results of the data-assimilation procedure to the cylinder flow case ( §3) and then to the backward facing step flow ( §4).

Theory

This section is devoted to the derivation of the model that will be used in the dataassimilation procedure. We consider an incompressible flow, whose velocity and pressure fields are given by the Navier-Stokes equations :

∂ t u + u • ∇u + ∇p -Re -1 ∆u = 0, ∇ • u = 0 (4.1)
where Re = U ∞ L/ν is the Reynolds number based on an arbitrary velocity U ∞ and a characteristic length L, which are used to make all the variables non-dimensional. Those equations can be recast into a more compact form:

B∂ t q + L 0 q + N (q, q) = 0 (4.2)
where the state vector q is composed by the concatenation of the velocity and pressure fields q = (u, p) and the operators B, L 0 (Stokes operator) and N (q 1 , q 2 ) (nonlinear convection, sometimes also denoted by N (u 1 , u 2 ) since it is independent on the pressure fields) are given by:

B = I 0 0 0 , L 0 = -Re -1 ∆() ∇() ∇ • () 0 , N (q 1 , q 2 ) = u 1 • ∇u 2 + u 2 • ∇u 1 0 (4.3)
Since we focus on statistically steady regimes, we may decompose the instantaneous state q(x, t) into a mean and a fluctuating component, as q(x, t) = q(x) + q (x, t) (4.4)

where the time-averaged operator is defined as

(•) = lim T →+∞ 1 T T 0 (•)(t)dt. (4.5)
Introducing the decomposition (4.4) into the Navier-Stokes equations (4.2) and taking the time-average yields the mean flow equation:

N (q, q) + L 0 q = -N (q , q ) = P f (4.6)
which differs from the fixed point solution equation by the presence of the right-hand side term f = -u • ∇u (acting on the equations through a prolongation P = (I, 0) T ) that expresses the time-averaged effect of the fluctuating velocity induced by the quadratic non-linearity (forcing term coming from the Reynolds stress tensor). The goal here is to model this term by approximating the fluctuation u by singular vectors of the Resolvent operator [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]). To do so, we subtract the mean flow equation (4.6) from the instantaneous equations (4.2), leading to:

B∂ t q + L q q = P f (4.7)
where the operator L q = L 0 + N q = L 0 + N (q, •) + N (•, q) is the linearized Navier-Stokes operator and f = -u • ∇u + u • ∇u represents the nonlinear fluctuations, acting as a forcing term. We then try to find solutions of this equation in the frequency space by using the ansatz q = qe iωt + c.c. and f = f e iωt + c.c. (where "c.c." stands for complex conjugate), leading to : (iωB + L q )q = P f (4.8)

This equations represents an (infinite-dimensional) input/output relation between the forcing term f and its perturbation q. In our approach we would like to simplify this equation, retaining only the dynamics that maximizes the gain in kinetic energy:

G(ω) = ||û|| 2 || f || 2 = (P q, P q) ( f , f ) = (R † (q; ω)R(q; ω) f , f ) ( f , f ) (4.9)
where R(q, ω) is the so-called Resolvent operator, defined as R(q, ω) = P T (iωB + L q ) -1 P and its adjoint R † (q, ω) given by the general definition (R † a, b) = (a, Rb), ∀a, b. Here we consider the energy inner product (a, b) = Ω a * • bdΩ.

We can show that the extraction of such dynamics is equivalent to the Singular Value Decomposition of the Resolvent operator R(q, ω), i.e:

R † (u; ω)R(u; ω) fk = µ 2 k (ω) fk , || fk || = 1 (4.10)
where µ 2 k are the Singular Values (ordered such that µ 2 0 ≥ µ 2 1 ≥ • • • ≥ 0), fk are the Right Singular Functions or the optimal forcing ( f0 being the optimal forcing term). The Left Singular Functions are defined as ûk = µ -1 k R(u; ω) fk . Alternatively, one can write the SVD problem for the variables ûk directly. To do so, we apply the operator µ -1 k R(u; ω) on equation (4.10), leading to:

R(u; ω)R † (u; ω)û k = µ 2 k (ω)û k , ||û k || = 1 (4.11)
If the nonlinear term f (ω) is expanded in the orthonormal basis constituted of the Right Singular Vectors fk and if we apply the Resolvent operator on it, we obtain an expansion of û(ω) in terms of the Left Singular Vectors:

û(x, ω) = +∞ k=0 µ k (ω)û k (x, ω) fk (ω), f (ω) (4.12)
This expression indicates that, when the forcing Fourier modes f (ω) is not precisely known, the perturbation û cannot be exactly determined. However, when the flow presents a strong convective unstable behavior, one dominant mode is expected such that µ 0 µ k , k > 0. In this case, a good approximation of the Fourier mode of the fluctuation is given by the first term in (4.12):

û(x, ω) ≈ µ 0 (ω)û 0 (x, ω) f0 (ω), f (ω) = A(ω)û 0 (x, ω) (4.13)
where the (complex) function A(ω) holds the information of the energy and phase of the mode at frequency ω. At this point, this approximation of the Fourier representation of the fluctuation is still not suited for numerical purposes as it is not yet 'discrete' in frequency space. In the following, we consider two distinct cases: first a simple case where the frequency content is mostly 'peaked' (oscillator flows) around one frequency and then, a more complex one, where it is more broadband (amplifier flow).

Time-periodic flows

Many flows present a periodic behaviour such as the flows around bluff bodies. In those cases the steady-state solution is unstable and the instability develops into a saturated periodic limit-cycle. Furthermore, it is not uncommon that most of the energy of this periodic limit-cycle is concentrated at the fundamental frequency ω 0 (see Turton et al. [2015a]). For those flows, a single-frequency description of the fluctuation may be acceptable:

u ≈ Aû 0 e iω 0 t + c.c. (4.14)
where the parameter A is now a single complex scalar. Based on this, we can approximate the Reynolds-stress tensor and the force acting on the mean flow equation following f = -2|A| 2 Re{û 0 • ∇û * 0 }. Considering the mean flow equation 4.6, an approximation ũ of the mean-flow u may be obtained through:

N (q, q) + L 0 q = P f = -2|A| 2 Re{N (û 0 , û * 0 )} (4.15a) R(ω 0 )R † (ω 0 )û 0 = µ 2 0 (ω 0 )û 0 (4.15b)
where R(ω 0 ) = R(q; ω 0 ). At this point, we remark that the system formed by the mean flow equation and the Resolvent analysis (4.15) is a closed set of nonlinear equations, if one knows the frequency of the nonlinear signal ω 0 and its energy |A| 2 . Arguably, already from a data-assimilation point of view, one can know a priori the frequency of the flow from a time-resolved signal from the flow, for instance a hot-wire or a point-wise wall-pressure measurement. An additional measure (amplitude of the oscillations, length of re-circulation bubble, strength of backflow, ...) can then be used to determine the energy of the fluctuation |A| 2 , since, for every value of |A| 2 one can solve for the nonlinear solution of 4.15 and therefore tune this parameter so that the solution matches to the actual measure. The output of this procedure is the full flow-field u composed of the reconstructed mean-flow ũ and the fluctuation u given in equation (4.14). The strength of the following approach lies in the simplicity of the optimization problem, since it only involves a single parameter contrarily to the procedure described in [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF], where complete spatial fields had to be determined.

Broadband Flows

For broadband flows , the previous approximation is no longer suited. In those cases we have thus to consider a multiple (N ) frequency approach {ω j } j=1,••• ,N , for which (A j = A(ω j )):

u = N j=1 A j ûj 0 e iω j t + c.c. (4.16)
This approximation is only valid if the dominant frequencies of the unsteadiness are located within a frequency range where the assumption of separation of the first singular value holds. This is usually the case in flows dominated by shear-driven unsteadiness, such as jets, wakes, shear-layers, etc.

If we consider this approximation in the mean flow equation 4.6, we have the following equation governing the mean flow ũ:

N (q, q) + L 0 q = P f = -2 N j=1 |A j | 2 Re{N (û j 0 , ûj, * 0 )} (4.17a) R(ω j )R † (ω j )û j 0 = µ 2 0 (ω j )û j 0 (4.17b)
The choice of the discrete frequencies {ω j } will depend on the knowledge of the physics of the problem, and can be determined from few point-wise measurements located at highly energetic regions of the flow. Moreover, similarly as before, the energy of each frequency |A j | 2 can be tuned using some measurements of the flow field through the model 4.17. The major difference here lies in the fact that the optimization procedure now involves N parameters

(|A 0 |, |A 1 |, • • • ) ∈ R N to be tuned.
The output of this model is as before the reconstructed mean flow ũ and the corresponding fluctuations at the sampling frequencies ω j , given by 4.16.

Numerical resolution of models

In this paragraph we provide some elements on how we chose to solve the system of nonlinear equations given by 4.15 (or 4.17). One possible way is to iteratively solve the mean flow for a given SVD mode and then solve a new SVD problem with the new mean flow and keep alternating in this way up to the point where the solution remains fixed. This method was employed by Mantic-Lugo et al. [2015] and Mantic-Lugo and Gallaire [2016] on a similar model than ours. However, it was observed that for relatively high values of |A|, the method needs a quite low under-relaxation factor, making the convergence quite slow. Here we propose an alternative, which is to regard the set of equations 4.15 as a whole nonlinear system and apply the Newton method to it. To do so, first of all, we note that the SVD eigensystem (equation 4.11) is not yet suited for this purpose, since, to extract the complete Jacobian matrix, we need to calculate, for example, the variations of those equations with respect to the state q, which does not show itself explicitly. To pose the same problem in a more suitable way, we write an eigensystem for the variables (q 0 , â0 ):

λ(iω 0 B + L q)q 0 = Bâ 0 , q * ,T 0 B q0 = 1 (iω 0 B + L q) * â0 = B q0 (4.18)
such that û0 = P T q0 , f0 = µ 2 0 P T â0 and λ = µ 2 0 . The variable â0 is defined as â0 = (iω 0 B +L q) -1, * B q0 . We remark here that, for the sake of clarity, we keep the 'continuous' notation introduced in previous paragraphs, even if now those objects are discrete, since we are dealing with the numerical method. For example, q, û0 , etc, are now vectors and B, L q, etc, matrices.

Together with those equations, we need to solve the mean-flow equation as well (for simplicity, we only show the single-frequency case given by equation 4.15). Its lineariza-3. An "Oscillator" Flow -Squared-Section Cylinder 51 tion, together with the linearization of the rewritten eigensystem leads to:

      L q |A| 2 N q * 0 (•) + |A| 2 N q0 (•) * 0 0 λN (•) q0 λ(iω 0 B + L q) (iω 0 B + L q)q 0 -B 0 q * ,T 0 B(•) 0 0 N * ,T (•) â0 -B 0 (iω 0 B + L q) * ,T            δ q δ q0 δλ δâ 0      =       -L 0 q -N (q, q) -|A| 2 N (q 0 , q * 0 ) -|A| 2 N (q * 0 , q0 ) Bâ 0 -λ(iω 0 B + L q)q 0 (1 -q * ,T 0 B q0 )/2 B q0 -(iω 0 B + L q) * ,T â0       (4.19)
The iterative solution of this system of equations and the update of the solution leads to the aforementioned Newton method. It was observed that the convergence of such method was often difficult and the reason for this is not fully understood. However, we could observe that, by changing slightly the equations, the convergence of the method was improved. First of all, we no longer suppose that λ is a real parameter, and we will let it free to evolve in the complex plane as the solution converges. Moreover, since now we have added a new variable to the system, namely the imaginary part of λ, we need a supplementary equation for the Jacobian of the system to be a square matrix. Note that if, instead of having an SVD problem, we had an eigenvalue one to be solved, involving a non-normal operator (such as, for example Mantic-Lugo et al. [2015]), the eigenvalue should not be considered to be real and this extra equation had to be considered. For this purpose, we remark that, from equations 4.10, 4.11 or 4.18, the Singular Modes are defined up to an arbitrary phase e iφ . We construct then the extra equation that fixes this phase, which can be, for example the imposition that the imaginary part of some point-wise quantity of the mode is null. This strategy leads to a nonlinear system that converges much better than the previous one. We remark that, in order to have a valid solution of the SVD system, the imaginary part of λ needs to be verified to be null once the residual has been decreased. Indeed, for all computations here presented, the ratio between imaginary and real parts of λ is of order 10 -12 . This linear system is then solved with an GMRES iterative solver, preconditioned with the lower triangular matrix (lower Gauss-Seidel), meaning that the equation for δ q are solved first, followed by the equation for δ q0 and δâ 0 .

An "Oscillator" Flow -Squared-Section Cylinder

The two-dimensional flow around a squared-section cylinder is a typical example of oscillator flows. For Reynolds number larger than Re ≈ 50 (see [START_REF] Sohankar | Low-reynolds-number flow around a square cylinder at incidence: study of blockage, onset of vortex shedding and outlet boundary condition[END_REF] for discussions around the critical value), it naturally develops an unsteady behaviour, corresponding to a periodic limit cycle, characterized by a single frequency and its harmonics. Below, we will first ( §3.1) present the details of the configuration and the numerical implementations. Then ( §3.2), we will consider a single time resolved measurement in the wake of the cylinder to reconstruct both the mean-flow and the perturbation dynamics. Finally ( §3.3), the case of other measurements solely based on mean-flow characteristics will be explored.

Configuration and numerical implementation

A sketch of the computational domain and the inflow boundary condition is provided in Figure 4.1. The inflow boundary is located 15 diameters away from the center of the cylinder and the velocity profile imposed is uniform. Two Reynolds numbers Re = DU ∞ /ν have been considered here, Re = 60, 100. The side walls are located 20 diameters away from the cylinder, on which we impose a slip condition. The outflow boundary condition is imposed 30 diameters downstream of the cylinder and reads (pI -Re -1 ∇u)•n = 0. Both the Direct Numerical Solution and the models are solved on the same two-dimensional mesh. All codes are second order in space and are based on finite elements. The DNS solver employs a second-order semi-implicit temporal scheme for time-advancement. All implementations have been performed in the FreeFEM++ open source code [START_REF] Hecht | New development in freefem++[END_REF]).

In Table 4.1, we show various quantities of interest to characterize the mean-and unsteady features of the flowfield: the re-circulation bubble's length L (for the steady solution and the time-averaged unsteady solutions), the fundamental frequency ω 0 and the amplitude |A DN S | of the Fourier mode at the fundamental frequency in the unsteady simulation. Results are given for 3 meshes of different density. All of them are strongly refined around the square and in its wake and are coarsened in the freeflow region. We can see that mesh 2 provides a good agreement with the finer mesh 3 for all observed quantities. Mesh 2 will therefore be chosen as the default mesh in all following computations. In particular, the model will be computed on this mesh.

Single point velocity measurement

Let us first assume that we only know the time evolution of (for example) the cross-stream velocity component at a single point in the wake of the cylinder, here x m = (x, y) = (2, 0). This quantity is represented in figures 4.2(a,d) for the two Reynolds numbers. The frequency content of this signal may be extracted by considering harmonic averages [START_REF] Arbabi | Study of dynamics in post-transient flows using Koopman mode decomposition[END_REF]: where F T (ψ)(ω) represents the Fourier transform on the finite time interval (0, T ):

ψ(ω) = lim T →+∞ 1 T F T (ψ)(ω), ( 4
F T (ψ)(ω) = T 0 e -iωt ψ(t)dt. (4.21)
The quantity ψ(ω) is finite and independent of the value of T . The amplitude |ǔ y (x m , ω)|, represented as a function of ω in figures (b,e), indicates that the fundamental mode at ω 0 strongly dominates those of the harmonics pω 0 with p ≥ 2. This suggests that a reconstruction with a single frequency may be sufficient. The fundamental frequencies of the periodic limit cycles are equal respectively to ω 0 = 0.77 and 0.91 for Re = 60 and Re = 100, while the amplitudes correspond to |ǔ y (x m , ω 0 )| = 0.09 and 0.21.

Since there is only one free parameter |A| in the model 4.15, we may explore its reconstruction ability by evaluating how the solution depends on |A| and comparing it with the actual measurement data shown in figures 4.2(a,b,d,e). From a numerical point of view, this does not necessarily represent a high computational effort since one may proceed by progressively increasing the value of |A| and restart each new computation from the solution corresponding to the highest available |A| (if no solution is available, we may start from |A| = 0, the steady solution). By doing so and by changing the value |A| slowly within the nonlinear solver, each nonlinear solver converges in a few iterations. This convergence deteriorates when the gap between two values of |A| becomes large. For example, if we desire to solve the model for a relatively high value of |A| from the steady solution (|A| = 0), we may need, for example, 10 iterations to solve it for |A| = 1, 5 or 20 iterations for |A| = 2.2.

Once those solutions obtained, in order to compare the model's solution in function of |A| with the available information, we apply the operator (•) to the fluctuation predicted by model (4.14): ). This is due to the fact that the model only considers the fundamental harmonic and lacks a representation of the higher order harmonics. The effect of the latter modes on the mean-flow distortion is taken into account by a slight overestimation of the amplitude of the fundamental Fourier mode.

|ǔ y (ω 0 )| = lim T →+∞ 1 T T 0 e -iω 0 u (x m , t)dt = |A||û 0 (x m )| (4.
The reconstructed solution of the model 4.15 at the crossing point is respectively shown in figure 4.3 (resp. 4.4) for Re = 60 (resp. 100). We can see that the reconstructed streamwise component of the mean-flow (a), the reconstructed streamwise component of the fundamental mode at frequency ω 0 (b) and the reconstructed streamwise meancomponent of the Reynolds stress tensor induced force (c) compare extremely well with their DNS counterparts (b,d,f). This is confirmed if we focus on the reconstruction of these quantities on the symmetry line y = 0: figure 4.5 shows that for both considered Reynolds numbers the streamwise profiles of mean-velocity, first-harmonic and force induced by Reynolds stress tensor closely follow the results obtained by the DNS.

Other measurements

Other quantities can be considered as inputs to tune the perturbation amplitude |A|. For example, this parameter can also be identified with mean quantities such as the mean-recirculation length, the mean-drag or the mean-velocity field. In figure 4.6, we have assessed the quality of the reconstruction considering the above mentioned mean-measurements.

We can see that the mean measurements based on the re-circulation length and meandrag vary monotonically (either increasing or decreasing) with |A|, showing that the stronger the fluctuations, the smaller the re-circulation bubble and the higher the average drag. As before, if we know the actual value of these mean-measurements, we can find a particular perturbation amplitude |A| that matches the measurement value. It turns out that all of these predictions for |A| are consistent and also correspond to the particular value of |A| for which the error on the full mean-velocity field is minimal. However, since these are mean-quantities, it is not possible to predict the fundamental frequency. To asses We can see that the re-circulation bubble's length (2nd row) and the mean drag (3rd row), at a fixed frequency, exhibit a behavior that is similar to the results shown before with ω = ω 0 : all curves are monotonic and they all eventually cross the actual value from the DNS. This means that we cannot uniquely determine a frequency from this data. If we observe now the global error in the mean-velocity field (1st row), we can see an elongated 'valley' in the frequency direction. This suggests that a precise knowledge of the chosen frequency is not crucial for the mean-flow estimation. The minima in these 'valleys' (black solid symbols) lie on points that do not necessarily correspond to the DNS frequency ω 0 (red solid symbols), since these minima are reached for (|A|, ω) = (1.68, 1.01) at Re = 60 and (|A|, ω) = (2.2, 1.08) at Re = 100. Note finally, that as discussed above, the amplitude parameter |A| that enables a fit to the measured DNS quantities (black solid symbols and thick solid lines) always slightly overestimate the amplitude of the Fourier mode in the DNS. This is due to the fact that the model lacks a representation of the higher order harmonics, whose effect on the mean-flow is taken into account by a slight overestimation of the amplitudes of the fundamental Fourier mode. At this point, we remind that a similar model was proposed by Mantic-Lugo et al.

[2015] for the flow around a circular cylinder, where instead of using the Resolvent formalism the classic stability analysis was used for the approximation of the harmonics. However, instead of taking external data to provide a criterion for the choice of |A|, their model stopped for a value of |A| such that the eigenmode was critically stable. This condition is justified whenever the effect of the second harmonics on the first is negligible. However, even if this method is similar to the here presented one, it is almost exclusive to "oscillator" flows, where its owns dynamics triggers the nonlinear saturation. For this reason, in the next section we show the generality of the proposed formalism on a backward-facing step flow, where the fluctuation is triggered by external noise, amplified through linear mechanisms [START_REF] Marquet | Global sustained perturbations in a backward-facing step flow[END_REF]).

Noise Amplifier Flows -Backward-facing step

In this section, we consider the (more challenging) backward facing step flow configuration. After a brief presentation of the configuration ( §4.1), we assess the reconstruction procedure in two distinct noise environements. The first case ( §4.2) corresponds to a quasi-periodic flow with two incommensurable frequencies, while the second one ( §4.3) refers to a broadband white noise excitation.

Configuration

The configuration chosen for representing the noise amplifier type flow is the Backwardfacing step described in [START_REF] Barkley | Three-dimensional instability in flow over a backward-facing step[END_REF], [START_REF] Hervé | A physics-based approach to flow control using system identication[END_REF]. The Reynolds number, based on the height of the step H and the maximum velocity U ∞ of the Poiseuille flow imposed at the inlet, is fixed at Re = U ∞ H/ν = 500. The inflow boundary condition is located at 5H upstream of the step and the outflow is located at 50H downstream, where, again, the boundary condition (pI -Re -1 ∇u)•n = 0 is imposed. The other boundaries of the domain correspond all to solid walls on which a no-slip boundary condition is imposed. A sketch of the computational domain is given in Figure 4.8.

For the present flow conditions, there exists a steady-state solution that is globally stable (see [START_REF] Barkley | Three-dimensional instability in flow over a backward-facing step[END_REF], [START_REF] Blackburn | Convective instability and transient growth in flow over a backward-facing step[END_REF]) but which exhibits strong transient mechanisms. The latter can be triggered by considering upstream stochastic noise, whose amplitude is tuned so that the dynamics of the perturbations is nonlinear and develops a

f (x, t) = 0 φ(x) w(t), (4.23)
where w(t) contains the temporal dependency (and thus frequency content) of the noise and φ(x) = e -|x-x G | 2 /σ 2 is a Gaussian function centered at x G = (-0.5, 0.25) and of width σ = 0.1 (for further details, see [START_REF] Hervé | A physics-based approach to flow control using system identication[END_REF]). The precise choice of w(t) will be given below and will enable us to consider different flow environments. In table 4.2, we provide some brief mesh-convergence results for the steady solution and the case of two incommensurable frequencies. Two meshes with different refinement levels have been generated. Both meshes are characterized by a higher refinement near the step location x = (0, 0) and a smooth coarsening towards the outflow (downstream of the main re-circulation bubble, which ends around x ≈ 10). We can see that for both meshes all observed quantities are almost identical (up to ∼ 3%). For this reason, we found reasonable to keep the coarser mesh (Mesh 2) as the default mesh for all subsequent computations.

4.2 Upstream noise with two incommensurable frequencies: ω = 0.48 and ω = 0.70

In this section, we investigate the case where the dynamics of the backward-facing step flow is triggered with two distinct frequencies. The first frequency is the one that maximizes the Resolvent gain of the base-flow, ω 1 = 0.48. The second one is chosen such that it lies in the zone where the Resolvent mode is still important and such that it is incommensurable with ω 1 . The second frequency is chosen to be ω 2 = 0.70. The time-signal presented in the external noise can thus be rewritten as:

w(t) = E 1 e iω 1 t + E 2 e iω 2 t + c.c. (4.24)
where the amplitudes of the forcing are equal to E 1 = E 2 = 0.5. Some characteristics of the unsteady flow are presented in figure 4.9 and in table 4.2. Compared to the base-flow, the mean-flow (figure 4.9 (a)) exhibits a shorter main recirculation bubble (L = 7.5 versus 10.9 for the base-flow) and also a smaller secondary bubble at the top wall (x d = 5.8 and x r = 10.4 versus 8.7 and 17.5 for the base-flow). The spectra of the cross-velocity ǔy , probed at y = 0.2 and x = 1, 4 and 15 show that as the perturbation is convected downstream, higher-order harmonics which are different ω 1 and ω 2 . The reconstructed flowfield (mean-flow and Fourier modes) for the optimal parameters (|A 1 |, |A 2 |) = (0.77, 0.44) corresponding to the minimum of ē and their exact DNS counterparts are presented in figure 4.11. We can see that, even if the energy of the mode for ω 1 is over-predicted, its spatial structure is well represented by the resolvent mode. This is also the case for ω 2 , especially in the high-amplitude region close to the bubbles. Comparing the model's Reynolds stress-tensor (u x u y component) with the one of the reference flow, we see that they have similar spatial structure and amplitude, confirming that the contribution of harmonics that are not considered in the model may be accounted for by overestimating the amplitudes of the modes retained in the model.

If we now consider measures involving the fluctuation amplitudes at given points, figures 4.10(e-j) indicate that the reconstructed fluctuation amplitude at frequency ω 1 (resp. ω 2 ) exhibits nearly a linear behaviour as a function of |A 1 | (resp. |A 2 |) and does only weakly depend on |A 2 | (resp. |A 1 | ), especially when x ≥ 10 (see figs (g-j)). Deviations from this behavior observed with the measurements at x = 4 may be attributed to the fact that, at those considered frequencies, the resolvent modes are only energetic if the re-circulation bubble's length L is sufficiently short. This can be seen through the similarity of the plots in figures 4.10 (b) and (e). Considering the measurements at x = 15 (see figs. i,j), the fluctuation amplitude at ω 1 obtained in the DNS (thick solid line in fig. i) is best-recovered for |A 1 | ≈ 0.4, while the best-fit for the amplitude at frequency ω 2 is obtained for |A 2 | ≈ 0.2 (see fig. j). Those amplitudes are not far from the DNS ones (red solid symbols). However, the mean-flow is not well represented with these parameter values (see figs a,b,c,d). This means that one cannot, with the two frequencies model, well represent at the same time the quantities related to the mean-flow and its fluctuations. This nevertheless points out an interesting feature: if one is interested only in the mean-flow and Reynolds-stress tensor reconstruction, one may consider a model with few frequencies and tune the few parameters by mean-flow measurements (such as L, x r or x d ). This may be considered as a very efficient way to perform mean-flow dataassimilation with extremely sparse data. However, if, in addition to the mean-flow, one desires to have a good estimate of the fluctuation, more frequencies need to be included in the model. This will be shown in the next section dealing with the BFS forced with white-noise external perturbation.

Upstream noise with broadband white noise forcing

We consider now the case where the dynamics is triggered with a (time) white-noise w(t). Since the noise-amplification mechanism associated to the linearized Navier-Stokes equations in the case of the present configuration is a broadband low-frequency band, we expect the output fluctuations to exhibit the same feature. This is in stark contrast with respect to the cylinder flow and the two-frequencies backward-facing step flow cases discussed before. The harmonic averages used above, that extract Koopman modes, now converge to zero for all frequencies. Following [START_REF] Arbabi | Study of dynamics in post-transient flows using Koopman mode decomposition[END_REF], a more suited quantity for the frequency analysis consists in considering the two-time auto-correlation function. For a given point x in space, this function can be given by:

T x,x (t) =E[u (x, τ )u (x, τ + t)] = lim T →+∞ 1 T T 0 u (τ )u (τ + t)dτ.
(4.25) Since no Koopman modes (i.e. quasi-periodic component) exist in the present configuration, this function T x,x (t) tends to 0 for large t, so that it is square integrable. Hence, (4.26) where the Fourier modes may be obtained from:

T x,x (t) = 1 2π +∞ -∞ T ω x,x e iωt dω,
T ω x,x = lim T →+∞ 1 T E[|F T (u )(ω)| 2 ]. (4.27)
For t = 0, we obtain that the time-averaged energy of the signal is equal to: (4.28) which establishes that T ω x,x is also the Power Spectral Density (PSD) of the signal. Coming back to the description of the flow unsteadiness, if we probe the resulting unsteady flow (presented in figure 4.12) at several flow locations, we can see that, as expected, the further we go downstream, the stronger and the more dominant are the low-frequencies in the PSD, shown in figure 4.12 (a-d).

T x,x (0) = lim T →+∞ 1 T T 0 u (τ ) 2 dτ = 1 2π +∞ -∞ T ω x,x dω = 1 π +∞ 0 T ω x,x dω,

Objective functional

We now turn to the model used for the reconstruction of the flow using partial measurements. The auto-correlation function in the case where the fluctuation field is approximated by eq. ( 4.16) reads:

T x,x (t) = N j=1 |A j | 2 |û j 0 (x)| 2 e -iω j t + e iω j t . (4.29)
For t = 0, the time-averaged energy of the signal therefore reads:

T x,x (0) = N j=1 2|A j | 2 |û j 0 (x)| 2 (4.30)
We therefore note that the PSD in the reference signal is dense (see ((4.28)) and figure 4.13 (a)), while the PSD in the model (see ((4.30)) and figure 4.13 (b)) is discrete. For the model, that contains a finite number of discrete Fourier modes, to be representative of the whole energy contained in the DNS, the two PSDs need to match. To establish this x,x in frequency domain (in equation (4.28)). To do so, we discretize the frequency space by considering {Ω j } j=1,N +1 , where ω j = (Ω j + Ω j+1 )/2 falls into the center of the interval (Ω j , Ω j+1 ). In our case, we fix the first frequency to be Ω 0 = 0 and Ω N +1 = 2, since Fourier modes corresponding to frequencies higher than Ω N +1 = 2 are weak and located upstream only (figure 4.12), near the location of the external forcing. Moreover, we pick a uniform frequency discretization such that ∆ω = Ω j+1 -Ω j for all j. This choice is in accordance with the actual PSDs of the DNS signal, in which only one 'bump' at low frequencies is visible. However, in the case where several 'bumps' in frequency would exist, a non-uniform distribution of Ω j s could be considered. The energy of the signal at point x contained within the frequency band (Ω j , Ω j+1 ) is:

I x,ω j = Ω j+1 Ω j T ω x,x dω = 2π|A j | 2 |û j 0 (x)| 2 (4.31)
With this, one can, as previously, compare the energy of the fluctuation at a measure point x m predicted by the model to the one measured. The measure of the error associated to this quantity is given by:

J ω j ,xm = 2π|A j | 2 |û j 0 (x)| 2 I x,ω j -1 2 . (4.32)
We note that this new objective functional to be minimized, contrarily to the meanflow related ones (such as L, x r , x d ), should be capable to discriminate efficiently between the different frequency contributions associated to each {|A j |} j=1,N . Indeed, following the recommendations of [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF], we will choose the measurement points at regions where the Resolvent modes are energetic to maximize the robustness of the reconstruction. The spatial distribution of the energy of the Resolvent modes is shown in figure 4.14. We can see that for high frequencies (ω ≥ 1), the resolvent modes are concentrated near the step location above the main re-circulation bubble, where for low frequencies the Resolvent modes tend to be more spread towards the outflow. For this reason, we choose to consider for data-assimilation the frequency content of different points of the domain. For high frequencies, we will consider the signal at x = 4, while for low frequencies, we consider two different possible points x = 10, 15. This leads to two cost functions to quantify the error for the data-assimilation procedure: Figure 4.15: Convergence of the amplitudes |A j | and the cost functional J n /J 0 for the cost functional J 1 (solid lines) and J 2 (dashed lines). Case N = 5.

J 1 = ω j <1 J ω j ,x=10 + ω j ≥1 J ω j ,x=4 , J 2 = ω j <1 J ω j ,x=15 + ω j ≥1 J ω j ,x=4 (4.33) (a) (b)

Optimization procedure

Now that the cost functionals have been precisely defined, we can proceed with the optimization procedure. We note that this white-noise case is richer in frequency domain.

For this reason, it is interesting to evaluate how the assimilated solution will change as a function of the number of considered frequencies N . For that purpose, we will consider the cases N = 2, 5, 10. We can see now that a parametric study (similar to the one performed for the two-frequencies case) is no longer feasible. Hence, we need to resort to an automatized optimization algorithm, where the cost functional J(|A j |) is minimized. For simplicity, here we choose a gradient-free algorithm, which provides good performance if the number of frequencies is not too large. We picked here the COBYLA algorithm, a simplex-based algorithm, where the optimal condition is evaluated within a N -simplex evolving in the parameter space. For further details on this algorithm, see [START_REF] Powell | A direct search optimization method that models the objective and constraint functions by linear interpolation[END_REF].

To speed up the convergence for the cases N = 5 and 10, we first perform the optimization for the case N = 2, which is quick. Then, we initialize the

Reconstruction results

In figure 4.15 we represent the convergence of the optimization process for the case N = 5 for two cost functionals described. We can see that starting from the case N = 2 provides indeed a good initial guess, especially for J 1 , since the optimization did not drift apart too much from the initialization and the optimisation did not take more than 20 iterations to converge up to a precision of J n /J 0 ≈ 10 -2 . For the cost functional J 2 , the convergence needed slightly more iterations, mainly because of a poor estimation of the initial amplitude for ω 2 = 1. Furthermore, we observe that the estimation of the amplitude for ω 0 = 0.2 and ω 1 = 0.6 is considerably lower than those predicted for J 1 . The rate of convergence of the cost functional J n /J 0 is almost exponential for both cost functionals considered. In table 4.3, we have compared, for the 2 cost functionals and 3 numbers of frequencies N , some overall features of the reconstructed flow with those of the DNS. Those features are the detachment and reattachment points, the mean-flow error and the the overall amplitude of the fluctuations, defined as

J 1 J 2 J 1 J 2 J 1 J 2 J 1 J 2 J 1 J 2 N =
A = 1 T T 0 Ω |u | 2 , which is equal to A = j 2|A j | 2
for our model. We could observe that, by refining the frequency discretization (increasing N ), the value of A increased for both cost functionals. This increase of this value lead to a mean-flow that was more saturated, leading to a decrease of almost all detachment / attachment points, with exception of x f for the cost functional J 1 . Indeed, this variable seemed less sensitive to the increase in A. Globally, the cost functional that met the best value for A was J 2 , whereas for best mean-flow approximation was J 1 . This is explained, in part, by the mismatch between the resolvent modes and the actual Fourier (Spectral Proper Orthogonal Decomposition, SPOD, as it will be explained shortly) modes. Indeed, in Figure 4.14, it is seen that the Resolvent modes only approximately reconstruct the Fourier modes and that significant discrepancies may exist. For example, for ω = 0.2, at the point used to tune the amplitudes (x = 10 or 15), the amplitude of the normalized Resolvent mode is either lower (x = 10) or higher (x = 15) than the normalized Fourier mode. In particular, tuning the Resolvent amplitude with the value of the Fourier mode at x = 10 will lead to an over-prediction of the fluctuation level and thus higher values of A.

We turn now our attention to the comparison of the assimilated fields with the ones coming from the DNS, namely the mean-flow and the Fourier modes. The definition of the Fourier modes is here based on the more general two-point, two-time correlation tensor (see [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]), whose components are the (time) correlation between the state vector at two distinct points in space:

T i,j (t) = E[u (x i , τ )u (x j , τ + t)].
Since there are no Koopman modes, all components of the correlation tensor tend to 0 Chapter 4. Resolvent-based Data-Assimilation when t becomes large, so that one may define the Fourier transform:

T ω i,j = F ∞ (T i,j )(ω). (4.34)
The diagonal entries (for which x i = x j ) of this tensor correspond to the PSDs defined earlier. This tensor can also be expressed in terms of the Fourier transforms of u as:

T ω i,j = lim T →+∞ 1 T E[F T (u )(x i , ω)F * T (u )(x j , ω)]. (4.35)
The explicit computation of this tensor from the DNS data is performed in the following way: we divide our total time-signal into N b 'bins' (with 50% of overlap), each one with a (sufficiently large) time-duration T , and we perform the average of each Fourier transform of each 'bin':

T ω ≈ 1 T 1 N b N b k=1 F k T (u )(ω)F k, * T (u )(ω) (4.36)
To extract the most coherent structure over all bins, ie the Fourier mode, we look for the largest eigenvalue of T ω :

T ω y ω = (λ ω 0 ) 2 y ω 0 , ||y ω 0 || = 1. (4.37)
The Fourier mode is then:

ǔSPOD (ω) = λ ω 0 y ω 0 . (4.38)
In figure 4.16 we compare the reconstructed fields for J 2 and N = 5 with the DNS. We can see that the mean-flow is well represented, as well as the Fourier modes. This shows the validity of the presented data-assimilation procedure, where, from point-wise time-resolved sensors, we manage to estimate the whole steady and unsteady features of a given flowfield.

Conclusions and Perspectives

In this chapter, we presented a data-assimilation technique based on Resolvent modes around the mean-flow. The technique was based on a model composed of a mean-flow equation coupled with Resolvent modes corresponding to left singular vectors of the Resolvent operator. The model is not closed since the amplitudes of those modes are unknown and need to be tuned with the help of measurements. This technique is essentially different from what was presented in the second and third chapters since the model has typically few degrees of freedom to be tuned, corresponding to the energies of each mode. This feature allows one to assimilate also very sparse data, without, for example, the need for supplementary penalization or regularization. This data can come from various types of measurements, such as mean-flow characteristics (length of separation) or point-wise time-resolved information (easily obtainable with, for example, a Pitot probe). Furthermore, the optimization / tuning procedure is much simpler, not involving necessarily an adjoint solver. This technique was applied on an oscillator flow, whose frequency content was essentially monochromatic. For this reason, the data-assimilation results showed that we are indeed capable of recovering not only the mean-flow but also the fluctuation around it. This technique was also applied on a backward-facing step where the dynamics is driven by an external forcing. Two cases were considered: the first where the forcing spectrum was discrete in frequency and composed of two incommensurable frequencies and, the second, continuous, with white-noise broadband forcing. In the first case, even if the data from the DNS showed more than those two peaks in frequency, we chose to use a two-frequencies model. With this model, we could only obtain partial information about the complete flow, meaning that we could either accurately reconstruct the unsteady features (the amplitudes of the actual Fourier modes) or the mean-flow-related quantities. A natural way of completing the model and solving those issues would be to add higher-order harmonics as done in Symon et al. [2019a]. However, this would lead to a much more complex model than the one proposed by them, since in their work the mean-flow equation was not solved, and the mean-flow was estimated from variational data-assimilation and frozen for all subsequent Resolvent analysis. The second case dealt with white-noise external perturbation, that lead to a dense (or broadband) frequency content. This is an essentially different test-case since before, the nature of the fluctuation was discrete in frequency (given by the Koopman modes, see [START_REF] Arbabi | Study of dynamics in post-transient flows using Koopman mode decomposition[END_REF]), and so was our model. For this reason, a connection between this broadband fluctuation and our model needed to be established. This connection is based on the integral of the Power-Spectral Density, which is a finite quantity either for the actual reference flow and for our model. By enforcing our model to satisfy those quantities at certain point locations, where the velocity field was measured, we were able to obtain a solution that reproduced this time both mean-flow quantities and the fluctuation as well. The reconstruction results show an increase of the fluctuation energy content when the frequency-refinement is increased. We could also see that the assimilated solution depends somewhat on the spatial location of the measurements. This comes from a mismatch between the SPOD modes and the Resolvent ones. However, as pointed out by [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF], those two modes are the same if the projection coefficient of the actual nonlinear forcing onto the Left-Singular Vectors (Resolvent-forcing) are uncorrelated between the optimal mode and its suboptimal ones. This suggests that, even in this broadband case, nonlinear interactions may be produced and may be important.

(c) (d) (e) (f) (g) (h) (i) (j) (k) (l)

Chapter 5

Mean-flow analysis for the mixed periodic/broadband fluctuations on turbulent flows 1 Introduction

The mean-flow stability analysis was formally introduced for transitional flows, particularly the flow around a circular cylinder, where the unsteady behavior was produced by a global instability that saturates into a periodic limit cycle. In this case, [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF] showed that the eigensystem formed by the linearized Navier-Stokes equations around the time-averaged (mean) flow lead to an eigenmode whose eigenvalue was (almost) purely imaginary, predicting very well the nonlinear frequency of the periodic limit-cycle for a range of Reynolds numbers near the bifurcation. Further applications [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF]) and theoretical elucidations (Sipp and Lebedev [2007b], [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF]) lead to the conclusion that, if the nonlinear forcing acting on the, say, first harmonic equation was negligible, this mean-flow analysis could reproduce the main features of the nonlinear limit-cycle, namely the frequency (from the eigenvalue) and the structure of this harmonic (from the eigenmode). In the laminar/transitional regime this concept has widely been successfully employed, for example for control purposes [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF]) or for designing models for the self-sustained dynamics of periodic limit-cycles [START_REF] Mantič-Lugo | Self-consistent mean flow description of the nonlinear saturation of the vortex shedding in the cylinder wake[END_REF]). In a turbulent context, for example in bluff-body configurations, this concept has also been applied with success in determining the frequency and the spatial structure of the vortex-shedding mode (see [START_REF] Mettot | Quasi-laminar stability and sensitivity analyses for turbulent flows: prediction of low-frequency unsteadiness and passive control[END_REF], [START_REF] Meliga | Sensitivity of 2-D turbulent flow past a D-shaped cylinder using global stability[END_REF].

Although the mean-flow stability analysis (based on the solution of the eigensystem of the linearized equations) was successful for determining dynamics in the context of oscillator flows, whose corresponding fluctuation is characterized by a discrete frequency distribution, noise-amplifier-like flows that hold a broadband-frequency representation, do not, in general, exhibit a neutral/unstable mode [START_REF] Ehrenstein | On two-dimensional temporal modes in spatially evolving open flows: the flat-plate boundary layer[END_REF], Ehrenstein and Gallaire [2008]). Rather, they may exhibit several (quasi-neutral) stable modes. In those cases, a more suitable formalism is given by the Resolvent operator that establishes the linear input/output transfer function from the linearized equations around the mean-flow. The Singular-Value Decomposition of this operator provides then a framework [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF]) where the most energetically-amplified input/output dynamics are evidenced. Furthermore, it was shown [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]) that, if those Singular-Values are well separated from one-another, the flow dynamics can be described by a low-rank approximation of the (infinite dimensional) Resolvent operator. This technique has been extensively applied in several contexts: to highlight coherent turbulent mechanisms [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], [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF]), as a tool for data-assimilation (Symon et al. [2019a]) and control [START_REF] Leclercq | Linear iterative method for closed-loop control of quasiperiodic flows[END_REF]), to mention a few.

In this chapter, we apply those concepts to the flow around a squared-section cylinder at Re = 22000 on the reference solution coming from a DNS run. This flow is characterized by a periodic component at a well defined frequency, corresponding to the vortex-emission, a common feature of flows around bluff-bodies, and a broadband-frequency one. This last contribution holds both a coherent part, corresponding to the well-defined structure of the Kelvin-Helmholtz instabilities happening at the shear-layer of the separation region, and an uncoherent/turbulent part. A sketch of such signal is depicted in figure 5. 1 (a).

In what follows, we wish to recover both coherent components, namely the vortexshedding and Kelvin-Helmholtz, through mean-flow stability or Resolvent analysis. To do so, we apply one version of the triple decomposition to the complete fully-developed turbulent signal. This decomposition was first introduced by Reynolds and [START_REF] Reynolds | The mechanics of an organized wave in turbulent shear flow. part 3. Theoretical models and comparisons with experiments[END_REF] and used to extract a periodic-wave component from the background turbulence. Further work, employed other versions of this decomposition in order to extract (not necessarily periodic) coherent structures [START_REF] Hussain | Coherent structures-reality and myth[END_REF], [START_REF] Hussain | Coherent structures and turbulence[END_REF]). Since then, several studies have been carried out using this framework in numerous contexts, for example, in turbulence modeling [START_REF] Speziale | Analytical methods for the development of reynolds-stress closures in turbulence[END_REF], [START_REF] Palkin | Scrutinizing URANS in shedding flows: the case of cylinder in cross-flow in the subcritical regime[END_REF]) or in experimental data-treatment [START_REF] Lyn | A laser-doppler velocimetry study of ensemble-averaged characteristics of the turbulent near wake of a square cylinder[END_REF]). More recently, [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF] proposed another version of this decomposition based on the Koopman operator formalism where the flow is separated in a quasi-periodic component (holding possibly multiple well-defined frequencies) from a broadband one. In this work, the latter will be preferred since its rigorous application leads automatically to two separate equations, one for the (quasi-)periodic component and the second one for the broadband component.

The first one corresponds to a nonlinear set of equations, where the influence of the broadband fluctuation (coherent and uncoherent) is given classically through the Reynolds-stress tensor, which is modeled by an eddy-viscosity (Boussinesq hypothesis), leading to a similar formalism as in URANS. This eddy-viscosity is then computed from treating the statistics from a DNS run. Upon further decomposition of the periodic component into mean-flow and periodic fluctuation, we can establish the mean-flow equations, where the eignvalue analysis is performed to recover the vortex-shedding mode.

The second one can be viewed as forced linear equation, where the forcing term comes from the quadratic nonlinearities of the broadband component itself and the linear operator corresponds to the linearized Navier-Stokes equations around the periodic component, in a similar way as in Floquet stability theory. In order to establish the input/output relation (Resolvent) in this Floquet-like framework, we consider that both the forcing term and the perturbation depend not only on the imposed frequency but also on the vortex-shedding time scale. By doing so, we are able to recover the dependency of the high-frequency Kelvin-Helmholtz modes on the phase of the slow and periodic movement of the shear layer, given by the vortex-shedding frequency. It is interesting to note that this phenomenon is observed in experiments (see [START_REF] Brun | Coherent structures and their frequency signature in the separated shear layer on the sides of a square cylinder[END_REF]) and also in the DNS performed in this work. For example, if we probe the unsteady flow at some point within the oscillatory motion of the shear-layer (figure 5.1 (b)), we can see that, according to the phase of the periodic movement, we can observe the Kelvin-Helmholtz structures or not at all. This chapter will be structured as follows: in the section §2, we will introduce and apply the concepts of the triple decomposition, giving rise to the equations governing periodic and broadband components. Then, in section §3 we will present the DNS in more detail and the mean-flow analysis deduced in the previous section. Finally, we provide some conclusions and discussions.

Triple Decomposition and Mean-Flow analysis

The goal of this section is to provide the theoretical framework for Resolvent analysis in general turbulent flows. We suppose here that those flows are governed by the incompressible Navier-Stokes equations:

∂ t u + u • ∇u + ∇p -∇ • (ν∇ s u) = 0, ∇ • u = 0, (5.1)
where ∇ s is a shortcut for ∇ + ∇ T . In general, a turbulent flow may hold a contribution whose frequency content is discrete (quasi-periodic part) and a contribution that is continuously spread over a frequency band. In so far, following [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF], we may decompose any fully-developed ergodic velocity / pressure field into a quasi-periodic component, u and p , and a broadband one u and p such that:

u(x, t) = u (x, t) + u (x, t), p(x, t) = p (x, t) + p (x, t).
(5.

2)

The extraction of the quasi-periodic component u (t) of any given signal u(t) can be formally expressed in terms of the Koopman operator theory. According to [START_REF] Arbabi | Study of dynamics in post-transient flows using Koopman mode decomposition[END_REF], [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF], there exists a countable set of frequencies {ω k } such that the harmonic averages, defined as (5.3) converge to a finite non-zero quantity. Typically this set of discrete frequencies is composed of few fundamental (incommensurable) frequencies, that arise, for example, from global instability mechanisms (absolute instabilities, [START_REF] Huerre | Absolute and convective instabilities in free shear layers[END_REF], [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: A linearized approach[END_REF]), and their nonlinear interactions (see, for example [START_REF] Leclercq | Linear iterative method for closed-loop control of quasiperiodic flows[END_REF]).

û(ω k ) = lim T →+∞ 1 T T 0 e -iω k t u(t) dt,
The resulting modes û(ω k ) are called Koopman modes. The quasi-periodic component is then reconstructed as the Koopman Mode Decomposition (KMD):

u (t) ≡ ω k e iω k t û(ω k ) + c.c., (5.4)
where, here we can see that the notation • is used to group together the extraction of the quasi-periodic component of u described by equations ( 5.3) and (5.4). The remainder of the original signal u = uu holds a broadband-like distribution in frequency, as the corresponding harmonic averages of û (ω) converge to zero for any frequency ω.

As stated by [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF] the Koopman approach provides an extension of the triple decomposition concept, that was initially introduced by Reynolds and [START_REF] Reynolds | The mechanics of an organized wave in turbulent shear flow. part 3. Theoretical models and comparisons with experiments[END_REF] and based on phase-average of the time-signals. We remark here that this decomposition does not necessarily separate the flow field into a coherent and an uncoherent part, since there are many well-known physical instabilities that manifest in the form of coherent structures over a continuous broadband frequency range. A typical example of this is the convective Kelvin-Helmholtz instability occurring in shear-layers (see [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]). In the following, we aim at modeling not only the coherent part linked to the quasi-periodic component but also the coherent contribution associated to the broadband component.

To proceed, we write the equations respectively governing the quasi-periodic component u and the broadband fluctuation u :

∂ t u + u • ∇ u +∇ p -∇ • (ν∇ s u ) = -∇ • u ⊗ u , (5.5) ∂ t u + u • ∇u + u • ∇ u +∇p -∇ • (ν∇ s u ) = f , ( 5.6)
where f = ∇ • ( u ⊗ uu ⊗ u ) represents the nonlinear interactions of the fluctuation u . The derivation of those equations relies on some properties of the operatior • that extracts the quasi-periodic component of a given flowfield such as: f g = f g + f g and f = 0. A more complete description of those properties is given in [START_REF] Reynolds | The mechanics of an organized wave in turbulent shear flow. part 3. Theoretical models and comparisons with experiments[END_REF] and are verified to be true even in the Koopman framework [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF].

In the following, we will separetely consider the modelling of the quasi-periodic component u and the broadband component u .

Modelling of quasi-periodic component

In equation (5.5) governing the quasi-periodic component, the term u ⊗ u corresponds to the Reynolds stress tensor stemming from the broadband contribution u . Even though u may involve a part related to coherent structures, we model this tensor using the Boussinesq hypothesis:

-u ⊗ u + 2 3 κ I ≈ ν t ∇ s u , ( 5.7)
where ν t = ν t (t) is an eddy viscosity that is, in general, quasi-periodic as well, and thus the slightly unusual notation for it around the brackets • . This approximation 2. Triple Decomposition and Mean-Flow analysis 75 results in the following equation:

∂ t u + u • ∇ u + ∇ p -∇ • ((ν + ν t )∇ s u ) = 0, (5.8)
where, in an abuse of notation, the turbulent kinetic energy κ = u • u /2 has been incorporated in the pressure term p . This equation governs the quasi-periodic part of the flowfield. If a turbulence model is provided f ( ν t , u ) = 0, then this equation is closed and may be solved for by timestepping. In general, such an equation may exhibit various types of solutions: fixed-points, periodic limit-cycle solutions linked for example to Hopf bifurcations of these fixed-points, or even multi-frequency quasi-periodic solutions linked to secondary instabilities of the periodic limit-cycles. Hence, classical stability analysis may be performed around fixed points [START_REF] Crouch | Origin of transonic buffet on aerofoils[END_REF][START_REF] Meliga | Sensitivity of 2-D turbulent flow past a D-shaped cylinder using global stability[END_REF], Mettot et al., 2014a] or Floquet type stability analysis around limit cycles. Yet, one has to keep in mind that only quasiperiodic solutions to this equation are "strictly" relevent.

Modelling with mean-flow marginal global modes

As explained in [START_REF] Mettot | Quasi-laminar stability and sensitivity analyses for turbulent flows: prediction of low-frequency unsteadiness and passive control[END_REF], a mean-flow approach may be favoured to characterize the unsteady quasi-periodic component of u from the sole knowledge of the time-averaged mean-flow ū. Indeed, from [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF], it is known that the Koopman mode corresponding to the fundamental harmonic of a time-periodic limit-cycle associated to equation (5.8) corresponds to a marginal eigenvalue of the linearized operator around the mean ū. There are actually two conditions for this result to hold: either the limitcycle should be purely harmonic (only the fundamental) or the amplitudes of all Koopman modes (Fourier modes) be weak. We thus further decompose the quasi-periodic component into its time-average (denoted with (•)) and its corresponding fluctuation (denoted with (•)), completing the triple decomposition. We remark that the eddy-viscosity ν t , having the same time-dependency as u , may have, in general, a fluctuation associated to it. Introducing the decompostion

u (x, t) = u(x) + ũ(x, t), p (x, t) = p(x) + p(x, t), ν t (x, t) = ν t (x) + νt (x, t)
(5.9) in the phase-average equation (5.8) yields the equations governing the mean-flow (with in particular additional Reynolds stresses stemming from the quasi-periodic component ũ):

u • ∇u + ∇p -∇ • ((ν + ν t )∇ s u) = -ũ • ∇ũ + ∇ • (ν t ∇ s ũ),
(5.10) and the linearized equations around the mean-flow:

∂ t ũ + u • ∇ũ + ũ • ∇u + ∇p -∇ • ((ν + ν t )∇ S ũ) -∇ • ( νt ∇ s u) = 0, (5.11)
If a turbulence model is at hand, we may explicitly relate the fluctuations νt to variations of the closure model. If not, we neglect this term, which we will assume in the following to make it simple. We finally look for eigenvalues of these linearized equations in the form: ũ = ûe λt , (5.12) so that (with a slight abuse of notations since the pressure term has been omitted)

λû + L(u, ν t )û = 0, (5.13)
where the operator L(u, ν t ) represents the linearized Navier-Stokes operator for a given mean-flow u and a given mean eddy-viscosity ν t and λ = σ+iω is the eigenvalue, dictating how the mode will decay (with σ) and oscillate (with ω). We can see that, in a natural manner, the mean-eddy-viscosity is incorporated in this operator. This is in accordance with [START_REF] Meliga | Sensitivity of 2-D turbulent flow past a D-shaped cylinder using global stability[END_REF], where, to obtain the vortex-shedding mode, the eddy-viscosity, coming from the URANS (Spalart-Allmaras) model, was taken into account. In [START_REF] Mettot | Quasi-laminar stability and sensitivity analyses for turbulent flows: prediction of low-frequency unsteadiness and passive control[END_REF], it was argued that in most of the cases, the frequency of the eigenmode was not really sensitive to νt , so that a crude approximation νt = 0 could also be considered in some cases. According to [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF], if one of the two conditions mentioned above are met, then the eigenvalue λ of (5.13) should be marginal with a frequency equal to the fundamental frequency of the limit cycle. Also, the eigenvector should correspond to the Koopman mode at the fundamental frequency. Note also that considering eigenvalues of the linearized operator (rather than resolvent modes) complies with the discrete nature of the Koopman modes.

How to obtain ν t

The mean eddy-viscosity νt may be obtained through different ways. If a turbulence model is used, then a straightforward time-average of ν t may be performed. If the fluctuation field u is obtained by a DNS, then by taking the time-average of equation 5.8, we obtain:

-u ⊗ u + 2 3 κI = ν t ∇ s u + νt ∇ s ũ. (5.14)
Applying the hypothesis νt = 0 mentioned above leads to the following mean eddyviscosity:

ν t = ∇ s u : (-u ⊗ u + 2/3κI) ∇ s u : ∇ s u .
(5.15)

Note that a DNS computation generally straightforwardly provides the mean-flow u, the periodic fluctuation ũ (by performing harmonic-averages) and the statistics u ⊗ u. These statistics may be related to the u ⊗ u through:

u ⊗ u = u ⊗ u -u ⊗ u -ũ ⊗ ũ.
The hypothesis (ν t = 0) may suffice as a first-order approximation, that leads to qualitatively correct global modes. However, a finer definition of νt could be found either by a more detailed study of the snapshots from the DNS (if available) or through a data-assimilation technique (see [START_REF] Camarri | Investigation of passive control of the wake past a thick plate by stability and sensitivity analysis of experimental data[END_REF]).

Resolvent analysis for broadband component

For the recovery of the broadband component, a frequency domain input-output approach, such as the resolvent formalism, is more suited. However, we remark that the linearized operator involved in the equations modeling the broadband component u , equation (5.6), is quasi-periodic, since the field around which the equations are linearized, u (t), is quasiperiodic:

∂ t u + L( u (t))u = f . (5.16)
Here, for the sake of clarity, we will suppose it is only periodic. This formalism is similar to the one appearing in Floquet stability theory, where we suppose that the solution may be written in the form u = ǔ(t)e λt , where ǔ(t) is a T -periodic mode that grows/decays with the Floquet exponent λ. Here, instead of analyzing the stability of the operator (as done in Floquet stability theory), we study its input/output relation. To do so, and similarly to Floquet stability analysis, we look for solutions of the forcing term and its corresponding response under the separated form: f = f (t)e iωt + c.c., u = ǔ(t)e iωt + c.c.

(5.17)

where ǔ(t), f (t) are T-periodic fields and ω is a given frequency. Here we suppose that the frequency ω does not coïncide with a Floquet exponent eigenvalue λ = iω of the periodic operator ∂ t + L( u (t)), which is the case if stable periodic limit-cycles are considered. We may then rewrite equation (5.16) as:

∂ t ǔ + iω ǔ + L( u (t))ǔ = f . (5.18)
This equation establishes the input/output relation mentioned before, since, for to a given T-periodic forcing term f , corresponds a T-periodic response ǔ, upon integration of 5.18. We remark however that only a T -periodic solution to this equation makes sense. One way of recovering this periodic solution can be to resort to feedback control techniques, where an extra forcing term is considered, penalizing the non-periodic component (similarly to [START_REF] Shaabani-Ardali | Time-delayed feedback technique for suppressing instabilities in time-periodic flow[END_REF]). Another way is to enforce periodicity by construction. This is the basic idea of, for example, the Floquet-Fourier-Hill Method (FFHM, see [START_REF] Deconinck | Computing spectra of linear operators using the Floquet-Fourier-Hill method[END_REF]), where the solution ǔ, the forcing term f and the operator L( u (t)) are expanded as Fourier series, i.e.: (5.19) where ω 0 = 2π/T and (ǔ n , fn , L n ) are the different harmonics. With this, equation (5.18) can be recast in the form:

ǔ(t) = n ǔn e inω 0 t , f (t) = n fn e inω 0 t , L( u (t)) = n L n e inω 0 t ,
(iωI + L)U = F, L =           . . . . . . . . . . . . . . . • • • L 0 + iω 0 I L -1 L -2 • • • • • • L 1 L 0 L -1 • • • • • • L 2 L 1 L 0 -iω 0 I • • • . . . . . . . . . . . . . . .           , ( 5.20)
where the vectors U, F represent collections of the harmonics of the perturbation and forcing terms. The matrix L is the so-called Hill matrix [START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF]) and plays the same role as the operator L(u) in the case of mean-flow analysis. We can see thus that, in this case, there exists a generalization for the Resolvent operator, where:

R(ω) = (iωI + L) -1 .
(5.21) Equation (5.20) formally establishes the existence of the T-periodic response ǔ(t) if ω does not coincide with any of the eigenvalues of L (since in such a case the matrix iωI + L is invertible). The gain to be optimized under the constraint U = RF should be viewed

Chapter 5. Mean-flow analysis for the mixed periodic/broadband fluctuations on turbulent flows as an energy-ratio between the fluctuation and the forcing term, both integrated over the period T . A possible form of such a gain is given by: (5.22) where the last equality is the application of Parseval's theorem and || • || denotes the usual space energy norm

G(ω) = 1 T T 0 ||ǔ(t)|| 2 dt 1 T T 0 || f (t)|| 2 dt = n ||ǔ n || 2 n || fn || 2 ≡ |||U||| 2 |||F||| 2 ,
|| • || 2 = Ω | • | 2 dΩ.
We can see that this energy gain is a natural extension of the usual steady one. This system is costly to solve, as it requires inversions of the operator iωI + L, if the FFHM is preferred or many forward/backward time-integrations, if time-methods are preferred. In what follows, we propose two simplifications, namely the steady-flow u ≈ u and the quasi-static approximations.

Steady flow approximation: u ≈ u

In the case u ≈ u, we may avoid the complexity of equations (5.18) or (5.20)and recover the more usual Resolvent formalism. In this approximation, the modes f , ǔ would evolve on top of the mean-flow and their time-dependency would be eliminated: ∂ t ǔ = 0. We recover then the classical Resolvent operator:

ǔ(ω) = R(u, ω) f (ω), R(u, ω) = (iωI + L(u)) -1 .
(5.23)

Following [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF], we are interested in the pair ( f , ǔ) that maximizes the energy gain:

G(ω) = ||ǔ|| 2 || f || 2 . (5.24)
This is equivalent to studying the Singular-Value Decomposition of R(

•) ≡ R(ω, u)(•) = i µ i ûi ( fi , •)
, where the modes ( fi , ǔi ) are obtaining by solving:

RR † fi = µ 2 i fi , R † Rǔ i = µ 2 i ǔi , ||ǔ i || = || fi || = 1 (5.25)
Note that the approximation u ≈ u may be valid as soon as the forcing fi and response ǔi (or more precisely the wavemaker region associated to ( fi , ǔi ), see [START_REF] Brandt | Effect of base-flow variation in noise amplifiers: the flat-plate boundary layer[END_REF]) are located in a region where u ≈ u.

Quasi-static Approximation

Another possible simplification of equations (5.18) is possible when the evolution timescale of the periodic component u is slow in comparison with the excitation frequency ω.

In those cases, we may postulate that the modes ǔ, f also evolve on the slow time-scale. This hypothesis can be formally translated by the introduction of the slow time-scale τ = εt that is involved in the fields ( u (τ ), ǔ(τ ), f (τ )). Therefore ∂ t ǔ = ε∂ τ ǔ and focusing on the leading-order equation, yields:

iω ǔ + L( u (τ ))ǔ = f ⇒ ǔ(τ ) = H τ f (τ ).
(5.26)

We can see that the resulting equation is nothing but the Resolvent operator linearized around the instantaneous phase-averaged flow H τ ≡ R(ω, u (τ )). This means that, in 3. Mean-Flow analysis on DNS data 79 this model, the modes f and ǔ coincide with the classical Resolvent analysis based on the high frequency ω and snapshots of u (τ ) taken along the slowly evolving quasiperiodic component. For this reason, we will refer to this approach as the quasi-static approximation. Furthermore, the gains from the SVD of the operator

H τ G(ω, τ ) = ||ǔ(τ )|| 2 || f (τ )|| 2 , ( 5.27)
that classically are solely dependent on the frequency ω are now dependent on τ as well and denoted by µ τ i (ω).

Remarks on the choice of averaging procedure and eddyviscosity

In the previous sections, we proposed a theoretical decomposition of the flow variables based on the Koopman formalism, leading to two separate equations, one for the quasiperiodic component u and one for the broadband one u . We remark, however, that the broadband component u , may itself be composed of a coherent and an uncoherent part u = u co + u unco , the first one related to spatio-temporal coherent structures driven by instability mechanisms, the second one related rather to chaotic/turbulent motion. In equation (5.6) governing u co , the effect of this latter uncoherent broadband term may be modelled by an additional eddy-viscosity term, say ν t which is a priori different from ν t . Such a term then appears in equations (5.16) and (5.20) to take into account the effects of the broadband uncoherent part u unco on u co . Such a modelling strategy for a broadband uncoherent part was already discussed by [START_REF] Morra | On the relevance of Reynolds stresses in resolvent analyses of turbulent wall-bounded flows[END_REF], who found, in a case where no quasi-periodic component existed ( u ≈ u), that considering an eddyviscosity term in the linearized equations for the fluctuations improved the prediction of the Resolvent modes. Furthermore, such an eddy viscosity term has been extensively used, see for example, [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], to model broadband large scale coherent structures in shear-flows. However, in few cases [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]], it appears sufficient to consider solely the molecular viscosity (as done in the previous section, 3 of this chapter) to reconstruct the coherent broadband fluctuation field. This may stem for example from a lower Reynolds number so that u ≈ u co , while at higher Reynolds numbers the coherent part is always accompanied by an uncoherent part. The modelling strategy described may be sketched in figure 5.2.

Mean-Flow analysis on DNS data

In this section, we illustrate the theory in the case of a squared-section cylinder at moderately high Reynolds number. First, we will present the Direct Numerical Simulation (DNS) results, exhibiting a well defined periodic component, corresponding to the vortexshedding phenomenon, and a broadband part, comprising in particular small Kelvin-Helmholtz type vortices arising in the shear-layer that issues from the detachment point. We will then analyse the fluctuation dynamics of the periodic and broadband components with the Resolvent based techniques presented in the previous section.

Direct Numerical Simulation solution

The configuration corresponds to the flow around a squared-section cylinder at Reynolds number Re = U ∞ Dρ/µ = 22000, based on the cylinder's diameter. It corresponds to a benchmark study case for numerical methods in turbulent flows [START_REF] Trias | Turbulent flow around a square cylinder at reynolds number 22,000: A dns study[END_REF], [START_REF] Nakayama | On the influence of numerical schemes and subgrid-stress models on large eddy simulation of turbulent flow past a square cylinder[END_REF], [START_REF] Verstappen | Fourth-order dns of flow past a square cylinder: First results[END_REF]) and the assessment of turbulence models [START_REF] Bosch | Simulation of vortex shedding past a square cylinder with different turbulence models[END_REF], [START_REF] Bao | Numerical prediction of aerodynamic characteristics of prismatic cylinder by finite element method with Spalart-Allmaras turbulence model[END_REF], [START_REF] Minguez | Experimental and high-order LES analysis of the flow in near-wall region of a square cylinder[END_REF]).

The DNS code solves the three-dimensional compressible Navier-Stokes equations, with a standard Sutherland law for the dynamic viscosity. The code is run at a low inflow Mach number, M = 0.1, representing a near-incompressible flow regime. The spatial discretisation used in the solver corresponds to a finite-volume method based on a modification of the AUSM+(P) scheme, which is second-order accurate. The timeintegration corresponds to the second-order accurate backward scheme of Gear, with a non-dimensional time step (based on the inflow velocity and in the diameter of the cylinder) of 3 × 10 -3 , which is around 100 times smaller than the highest observed frequency in the flow corresponding to the Kelvin-Helmholtz instability. For further details, see [START_REF] Dandois | Large-eddy simulation of laminar transonic buffet[END_REF] and references herein. We use the FastS code, developed by ONERA, which is a highly optimised flow solver for high performance computing clusters.

The time-averaged flow is shown in figure 5.4 (b). A re-circulation length of L = 0.97 (measured from the center of the square) is obtained, which is in close agreement with the literature (for example, L = 1.04 by [START_REF] Trias | Turbulent flow around a square cylinder at reynolds number 22,000: A dns study[END_REF]). The main physical unsteady feature of the flow corresponds to the vortex shedding phenomenon, as also observed in other bluff-body configurations. The frequency of this phenomenon can be seen in, for example, figure 5.3 (a), which represents the Fourier Transform of the cross-stream velocity signal at a point on the symmetry axis in the wake of the cylinder. We can see that the peak in frequency is sharp and corresponds to the frequency ω 0 ≈ 0.83 (St ≈ 0.133, which is in a good agreement with [START_REF] Rodi | Comparison of LES and RANS calculations of the flow around bluff bodies[END_REF], [START_REF] Bosch | Simulation of vortex shedding past a square cylinder with different turbulence models[END_REF], [START_REF] Trias | Turbulent flow around a square cylinder at reynolds number 22,000: A dns study[END_REF]). The spectrum also exhibits smaller sharp peaks at the harmonics, due to the nonlinear interactions associated to the fundamental mode at ω 0 . Coming back to the triple decomposition framework introduced in the last section, we will consider this phenomenon as being discrete in frequency, and therefore accounted for in u . The structure of the Koopman mode at ω = ω 0 , given by the harmonic averaging procedure, is shown in 5.4 (b). The flowfield also holds a contribution at higher frequencies, shown as a 'bump' in figure 5.3 (b), which is relative to the frequency spectrum of the stream-wise velocity near the upstream corner of the square. This broadband bump corresponds to small-scale convective Kelvin-Helmholtz instabilities (see [START_REF] Brun | Coherent structures and their frequency signature in the separated shear layer on the sides of a square cylinder[END_REF], [START_REF] Trias | Turbulent flow around a square cylinder at reynolds number 22,000: A dns study[END_REF]) evolving in the shear-layer produced by the separation. This instability manifests itself over a broad range of frequencies ω KH ≈ 10 -20 approximately. For this reason, we will consider this phenomenon belonging to the broadband component u .

Mean-Flow Stability Analysis: periodic component

In this section, we perform the stability analysis on the mean-flow operator to recover the periodic component of the flow. As mentioned before, this analysis takes into account a mean-eddy-viscosity ν t which models the broadband fluctuations and will be computed from the DNS statistics. This quantity is given by equation 5.15. We can see from this definition that its spatial distribution is dependent on the numerator prod = ∇ s u : (-u ⊗ u + 2/3κI), also referred to as the production term, and the denominator dest = ∇ s u : ∇ s u, referred to as the destruction term. The destruction (plotted in 5.5 (a)) is a positive quantity and is only dependent on the gradients of the mean-flow. We can see that this quantity is strong at the shear-layers and decays exponentially away from the cylinder. The production term (plotted in 5.5 (b)) also follows a similar spatial dependency. We remark that at some regions of the flow near the corners of the cylinder and at the far wake, the production is negative. At those regions, in the computation of the eddy-viscosity, the production value is set to zero. This is done to avoid ill-posedness problem for the operator L(u, ν t ). The eddy viscosity is then given by the ratio between the production and the destruction. However, at regions where the destruction term is small, the value of ν t can be large. To avoid this issue, we use instead a denominator under the form √ dest 2 + ε, where ε is a regularizing parameter. We plot the resulting ν t for ε = 0.01 and ε = 10 in figure 5.5 (c,d). We can see that for ε = 0.01, the eddy-viscosity is indeed large, leading to local Reynolds numbers of order 10. For this reason, we will assume, if not mentioned otherwise, that the quantity ν t was computed with ε = 10, although some tests will be performed with ε = 0.01 as well.

With this definition of the eddy-viscosity we are now in means of performing the stability analysis for the recovery of the periodic component. In figure 5.6 we provide those results for three different strategies: the laminar one (assuming ν t = 0, as in [START_REF] Mettot | Quasi-laminar stability and sensitivity analyses for turbulent flows: prediction of low-frequency unsteadiness and passive control[END_REF]), the frozen eddy-viscosity (as presented in section 2) and a third one which takes into account the fluctuation of νt . To do so, we assume that the eddy-viscosity is governed by a URANS turbulence model (here the Spalart-Allmaras, see [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF]) and suppose that the mean-flow and the eddy-viscosity computed from the DNS are near to the ones issuing from a URANS run for this model. This hypothesis, at least for the mean-flow u, could be supported since mean-flow URANS results have been reported as been comparable with experiments and with the present DNS, for example in [START_REF] Meliga | Experimental and theoretical sensitivity analysis of turbulent flow past a square cylinder[END_REF], where their re-circulation bubble was around L = 1.11.

We can see that, for the laminar analysis (figure 5.6 (a,b)), although the frequency is very well predicted (ω = 0.80), the unstable mode is far away from the neutral line (σ = 0.38), and it corresponds to a mode whose spatial distribution is localized at a thin shear-region around the square and almost nonexistent towards the outflow. However, we can see that, taking into account the eddy-viscosity ν t in the analysis (figure 5.6 (c,d)), we are able to recover the overall features of the DNS Koopman mode (figure 5.4 (b)). Its eigenvalue is much closer to the neutral line (σ = 0.03), presenting now a lower frequency than before (ω = 0.69). In order to investigate if this effect comes from a lack of modeling of the eddy-viscosity fluctuation, we perform as well the stability analysis considering the Spalart-Allmaras model (figure 5.6 (e,f)), as described before. We can see that, indeed, the predicted frequency now is bigger than before (ω = 0.75), suggesting that, indeed, the fluctuation of νt can have a positive effect on the eigenmode. We remark that in a URANS formalism, this fluctuation can be straightforwardly taken into account by the consideration of the complete Jacobian of the nonlinear system. Lastly, another question we may ask ourselves is the effect of the chosen eddy-viscosity between the ones computed with the two regularization parameters introduced before ε = 10/0.01. To answer this, we provide the modes coming from the eddy-viscosity computed with ε = 0.01 on the eigen-spectra in figure 5.6 (c,e) by a red dot (the case ε = 10 was already discussed). We can see that those modes tend to have a lower frequency and they tend to be closer to the neutral axis, which is coherent with the fact that this eddy-viscosity is stronger than the previous one. However, even with this (quite important) change in the eddy-viscosity choice, the modification of the mode was not too important. This suggests that the region most sensible to the quality eigenvalue is not at the wake of the cylinder, but at the shear layers, where the destruction term dest is very important and the regularization procedure is not active.

Resolvent analysis: broadband component

We now turn our attention to the reconstruction of the broadband component. First, we will provide the results for the quasi-static approximation, then we will compare them with the steady-flow one. As mentioned in section 2, this approach is valid if the frequency ω is much higher than the one coming from the periodic flow, ω 0 = 0.83. Thus, we place ourselves in the interval ω ∈ (5, 35). The Resolvent gains for ω > 35 quickly decay to zero and no more energetic input/output dynamics is recovered. In this interval, the Resolvent gains are given in figure 5.7 (a-d) for four different times (or phases) of the baseline period, corresponding to τ /T = 0, 1/4, 2/4 and 3/4, where the initial time was set such that the y-velocity of the periodic flow at (0.5, 1) is maximal. From those plots, we can see a clear "bump" on the frequency range ω ≈ 10 -20, which corresponds well to the region in frequency where the probes on the DNS accused the Kelvin-Helmholtz instabilities. We can see from the gains as well that the two first modes have similar energies while the third and fourth modes are several orders of magnitude less energetic. This suggests that they do not participate actively in the nonlinear fluctuation of the actual signal. It is interesting to note that this "bump" corresponding to the leading energy gain is also present in the steady-flow analysis, u ≈ u, (blue/red curves on figure 5.7 (e), although a little bit less energetic. We can say then that the present analysis, considering the time-periodic Resolvent operator corresponds to a refinement of the classical mean-flow procedure.

It is interesting to note that the baseline periodic flow has a special space/time symmetry, which reflects itself on the modes. Indeed, the fact that the Koopman mode (figure 5.4 (b)) is anti-symmetric with respect to the mean-flow implies that at a given phase τ , the field u (τ ) is the mirror-symmetry of the field u (τ + T /2). This means that the same physical mechanisms observed at a given phase τ will also be present in the phase τ + T /2, meaning also that the same Resolvent modes found at τ should be found as well at τ + T /2 as their mirror-symmetric counterpart. Indeed this phenomenon is observed. We can see that the upper modes, given for the frequencies ω = 10/20, (marked in red in figure 5.7 (f-i,k-n)) are mirror-symmetry version of the lower (blue) ones after a semi- period T /2. This is also true for their corresponding Resolvent gains, corresponding to the same colour code as the modes. As for the static approach, we can see that, since the mean-flow is symmetric, the two upper/lower modes do not exist anymore. Instead, we have two identical modes, having the same Resolvent gain, corresponding this time to a symmetric or an anti-symmetric structure. Their anti-symmetric counterpart only is plotted in figure 5.7 (j,o), also for ω = 10/20. We compare now those modes, obtained with the quasi-static approximation, with (span-wise-averaged) snapshots of the actual DNS flow at two different phases, namely τ /T = 0 and 1/2. This is achieved in figure 5.8. We can see that, indeed, the broadband fluctuation, represented there by the pressure p = p -p , depends on the phase of the periodic velocity field. Furthermore, this fluctuation responds almost instantaneously with changes of u , validating the quasi-static approximation in this case. An interesting observed feature is the fact that at the phase τ /T = 0, we can identify some structures on the pressure snapshot that have the same wavelength as the mode obtained for ω = 20 and maybe some larger structures corresponding to the lower frequencies towards the rightmost corner. This is coherent with the energy gains observed by the upper mode at that phase (red curve in figure 5.7 (a)), that presents a peak near ω = 20. As for the phase τ /T = 1/2, we can see that, most of the fluctuation on the shear layer is comparable with the lower frequency ω = 10 instead. This is also coherent with the energy gains (red curve in figure 5.7 (c)), which presents higher values for lower frequencies and almost no amplification for frequencies such as ω = 20. Also, by comparing the fluctuations with the modes coming from the steady-flow approximation, we can see that they do not match so well, as those modes do not oscillate with the shear-layers.

Conclusion

In this chapter, we performed flowfield reconstructions, based on linearized Navier-Stokes equations, for a turbulent flow holding both a periodic large-scale low-frequency vortex shedding component and a broadband small-scale high-frequency component (Kelvin-Helmholtz type instabilities) evolving on top of the vortex-shedding mode. The triple decomposition [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF] separates in an accurate and rigorous manner these two components. The approaches employed for flow reconstruction are different for both components.

For the periodic component, a mean-flow stability analysis with an eddy-viscosity modelling had to be taken into account. The reconstruction accuracy even improved when augmenting the Jacobian operator with a linearized Spalart-Allmaras model governing the eddy-viscosity fluctuations.

For the broadband component, input-output linearized Navier-Stokes equations around the vortex-shedding component naturally step in. This leads to optimal forcing and response modes whose structure depends on the phase of the vortex-shedding cycle. We note that with the triple-decomposition of the flow-field used here, no eddy-viscosity is taken into account for those fluctuations; however, we believe that, if finer statistics of the fluctuations were provided, an eddy-viscosity associated with the uncoherent/chaotic component of the broadband spectrum could be taken into account and could improve reconstruction accuracy. The general reconstruction problem involves an expensive inputoutput problem to solve. We proposed two approximations, for which the complex inputoutput problem simplifies into a classical resolvent analysis: the quasi-static and the steady-state approximations. It was shown that the quasi-static approach gave finer results than the steady-state approach, and that the structure of the quasi-static resolvent modes very accurately compared with DNS snapshots at different phases of the vortex shedding cycle.

Chapter 6

Resolvent-based Data-Assimilation for the turbulent flow around a squared-section cylinder 1 Introduction

In the fourth chapter, we explored the data-assimilation based on a model composed by the mean-flow equation, coupled with the mean-flow resolvent modes, modeling the fluctuation. In that context, the Reynolds number was small and the full complexity of a general turbulent flow could not be assessed or discussed. Furthermore, in those cases, we studied flows having either a quasi-periodic fluctuation, in which case only one or two well-defined frequencies were taken into account by the model, or a (coherent) broadband one, for which a frequency-discretization was employed. Generally, for turbulent flows, we potentially have both behaviors together, as discussed in the fifth chapter. For this reason, in the present chapter, we attempt to extend the idea of resolvent-based dataassimilation to turbulent flows. More precisely, we wish to apply those ideas on the same physical configuration considered in the fifth chapter, namely the turbulent flow around a squared-section cylinder, where the same DNS is used as the reference solution. To do so, we need a model composed by the mean-flow equation coupled with unsteady modes coming from mean-flow analysis, whose amplitudes will be used as tuning parameters. However, we can see now that designing a model for this configuration is a much more challenging task than before. Indeed, in order to obtain the correct vortex-shedding structure, necessary for the proper reconstruction of the Reynolds-stress tensor, we need an eddy-viscosity, as shown in fifth chapter, which is a spatially-dependent function that models the broadband nonlinear interactions. One could argue that this quantity could be computed from the resolvent modes modeling the broadband component, in a similar way we did for the white-noise forced BFS, presented in the fourth chapter. This would however represent a difficult task since now the broadband signal is richer, presenting not only coherent structures, as the Kelvin-Helmholtz modes but also uncoherent ones. Furthermore, if we look to capture as well the time-dependency of those modes with the vortex-shedding frequency (present in the quasi-static approximation), this time scale needs also to be discretized, increasing even more the complexity of the model. For those reasons, we choose instead to rely on an URANS turbulence model, here the Spalart-90 Chapter 6. Resolvent-based Data-Assimilation for the turbulent flow around a squared-section cylinder

Allmaras. This is useful since it naturally provides a physically plausible eddy-viscosity. Moreover, as we will see, this model predicts correctly the periodic component from the DNS. This chapter will be structured as follows: in the section §2, we will design the model, based on URANS equations and its resolvent analysis, used in the data-assimilation context. Then, in section §3 we will present the data-assimilation results, followed by conclusions.

Development of the model for Data-Assimilation

In this section we will present an extension of the model used for data-assimilation process discussed in the fourth chapter in the laminar/transitional case. As discussed before, the inclusion of broadband fluctuations in the model, coupling them with the mean-flow equation through the Reynolds-stress tensor and with the vortex-shedding mode through the eddy-viscosity leaves the model too complex. For this reason, we will drastically simplify the problem by considering a URANS turbulence model. The chosen model will the the Spalart-Allmaras model, as introduced in the second chapter. The solution of this model, together with its steady version (RANS) will be presented in the next paragraph, together with a model for the data-assimilation.

RANS and URANS Spalart-Allmaras Solutions

The starting point of our data-assimilation for this turbulent case will be the evaluation of the RANS-SA solution, presented in figure 6.1. This helps to quantify how far we are from the actual DNS solution by solving a spatial-only problem. This solution (U, P, Ñ ) is supposed to represent, in a way, the mean-flow U ≈ u. The RANS-SA equations can be summarized in:

U • ∇U + ∇P -∇ • (ν + ν t ( Ñ ))∇ s U = 0 (6.1) ∇ • U = 0 (6.2) U • ∇ Ñ -∇ • (η( Ñ )∇ Ñ ) = s( Ñ , ∇ Ñ , ∇U), (6.3)
where s = s( Ñ , ∇ Ñ , ∇U) is the nonlinear source term for the turbulent variable. For a better description of it, see Appendix A. We can see however that this model fails to capture its overall characteristics, presenting a much larger re-circulation bubble L ≈ 3.8. This is actually a known fact that RANS modeling fails to predict the mean-flow in bluffbody-like configuration. This comes from the fact that this kind of model was designed for capturing small-scale turbulent motion. However, in this configuration, an important large-scale motion exists, the vortex-shedding. Indeed, this flow is unstable and, if we allow it to evolve in time (in an URANS manner), it will produce a much shorter meanre-circulation bubble (see [START_REF] Bosch | Simulation of vortex shedding past a square cylinder with different turbulence models[END_REF]). Those URANS equations are given by: The solution of those equations corresponds to a saturated periodic limit-cycle, whose mean-flow (both velocity and eddy-viscosity) and nonlinear fluctuation, given by the first and second harmonics are presented in figure 6.2. The mean-flow is similar to the one coming from the DNS, with L ≈ 1.4. The nonlinear frequency of this solution happens to be very close to the one coming from the DNS (ω ≈ 0.83). Also, the structure of the first harmonic is very similar as well to the DNS one. Those results are in accordance to the work of [START_REF] Bosch | Simulation of vortex shedding past a square cylinder with different turbulence models[END_REF]. We remark that most of the energy of the fluctuation is concentrated on the first harmonic (at ω 0 ), the second one having a norm 6.7 times weaker than the first one. We can see thus that this URANS framework provides a good basis for capturing the vortex-shedding mode and for modeling the broadband structures. For this reason, it will be used as guideline for the construction of the model in the following. This model will focus on coupling the mean-flow equations with the vortex-shedding mode only, in a very similar way we did for the laminar square-cylinder case, present in the fourth chapter. However, before doing so, we investigate the predictability power of the resolvent formalism in this case.

∂ t U + U • ∇U + ∇P -∇ • (ν + ν t ( Ñ ))∇ s U = 0 (6.4) ∇ • U = 0 (6.5) ∂ t Ñ + U • ∇ Ñ -∇ • (η( Ñ )∇ Ñ ) = s( Ñ , ∇ Ñ , ∇U) (6.6)

URANS Mean-Flow resolvent Analysis

In order to assess the effectiveness of the resolvent analysis in describing the vortexshedding mode, we perform here two distinct ones: the frozen-eddy-viscosity, where only the momentum and continuity equations are taken in the Jacobian (where the fluctuation is taken under the form (U , P ) = ( Û, P )e iωt ), and the variable eddy-viscosity one, where the SA model is linearized and included in the Jacobian as well (and the fluctuation is taken under the form (U , P , Ñ ) = ( Û, P , N )e iωt ). Those analyses arise when we decompose the time-dependent variables U(t), P (t)(, Ñ (t)) into their mean-flow U, P , N and fluctuation U , P (, Ñ ), in a similar way we did in the last chapter. Note that for the variable eddy-viscosity case, the linearization of the model needs to be done in order to close the Jacobian matrix, which will now contain a variation Ñ in the state vector. We remark that, for this data-assimilation model, we will stay consistent with the fourth chapter modeling, meaning that we will perform resolvent analysis for the recovery of the vortex-shedding mode rather than a classical stability analysis. The reason for this is that, in a way, the resolvent approach can be more general, allowing one to impose the exact frequency coming from a time-resolved probe, avoiding the possible mismatch of frequencies presented in the last chapter. Furthermore, in configurations that are dominated by only one globally unstable mechanism, those approaches lead to similar results (see [START_REF] Symon | Non-normality and classification of amplification mechanisms in stability and resolvent analysis[END_REF]). We can see from figure 6.3 that the gains produced from both approaches are very close together and have a maximal value for nearly the same frequency. The maximal frequency for the frozen ν t approach is slightly lower than the one coming from the DNS, as it was the case for the DNS-based eddy-viscosity analysis discussed in the last chapter. However, despite this small difference, we can see that the structure of the modes (both computed at ω 0 = 0.83) is virtually the same. This means that, for modeling the vortexshedding purposes, the frozen eddy-viscosity may suffice. Turning now our attention to the eddy-viscosity fluctuation, in the case where it is modeled with the linearized SA equation, we can see that its structure is not corresponding to the one coming from the URANS simulations (see figure 6.4). Moreover, it is interesting to notice as well that, in the URANS computations, the second harmonics (2ω 0 ) of this quantity is of the same order of magnitude as the first one. This shows that the dynamics of the eddy-viscosity is more complex than the one for the velocity field only. With this in mind, we are ready to design our model for data-assimilation.

Data-Assimilation Model

In order to mimic the behavior of the URANS equations in a minimal way, the proposed model will be composed by the mean-flow moment equation, coupled with the SA model for the mean turbulent variable, together with the resolvent mode, modeling the vortexshedding. The complete nonlinear model is given by:

U • ∇U + ∇P -∇ • (ν + ν t (N ))∇ s U = -2|A| 2 Re{ Û * 0 • ∇ Û0 } (6.7) ∇ • U = 0 (6.8) U • ∇N -∇ • (η(N )∇N ) = s(N , ∇N , ∇U) (6.9) R(U, N , ω 0 )R † (U, N , ω 0 ) Û0 = µ 2 0 Û0 ,
where the fluctuation of the velocity field is approximated by: U = A Û0 e iω 0 t + c.c. (6.10) and the operator R = R(U, N , ω) corresponds to the linearized Navier-Stokes operator around the mean-flow U with the (frozen) eddy-viscosity N . We notice that the nonlinearities in the mean-momentum equation (and in the SA equation for N ) coming from the fluctuation of the eddy-viscosity are neglected since we do not have a good estimation for it, as discussed in the last paragraph. Moreover, the second-harmonics of the velocity field is also neglected, since its energy is 6.7 times lower for the URANS simulation and 8.6 times lower for the DNS. The only non-linearity taken into account in this model, other than those already existent for the RANS-SA model itself, is the advection in the momentum equation, similarly to what was proposed in the fourth chapter. Note that, although we base our model on a URANS model, the amplitude of the vortex-shedding mode is not necessarily the same as the one coming from the URANS run, and will be tuned using external data. This gives the model and external parameter for it to adapt to the DNS solution in a way that the URANS model itself cannot. In the next section, we will present some results on the dependency of this model's solution with respect to the parameter |A|, also tuning the energy of the mode. We remark that, since the model consists of a system of equations very similar to the one presented in the fourth chapter, we use the same strategy to solve it numerically. This means we solve it with a Newton method, once we rewrite the system as a set of nonlinear equations. It is worthmentioning that, in order accelerate the resolution of the model, only a half-domain mesh was employed, where, at the center-line, symmetry conditions were used for the mean-flow and anti-symmetry conditions were used for the vortex-shedding mode.

Results

In this section, we study how the solution of the previously introduced model varies with |A|, in the same manner we did for the laminar case, where a parametric study is employed, shown in figure 6.5 (continuous black lines). We can see that, by increasing the value of |A|, the mean-flow related quantities represent better the actual DNS ones. The re-circulation bubble becomes shorter and the mean-flow error (with respect to the DNS mean-flow) also decreases and creates a minimum for |A| ≈ 2.18. We can see as well that the quantities related to the fluctuation do not match the DNS ones necessarily at the same value of |A|. For example, if we look at the magnitude of the resolvent mode at a point close to the cylinder (|A|| Û0 |), we see that the model does not manage to react fast when |A| is small. This issue was already discussed in the fourth chapter and it is related to the fact that, at this range of parameter, the bubble is still too long and the resolvent mode does not have important values close to the cylinder. However, for points far from the cylinder, this quantity is almost linear. We can see that even quantities related to the resolvent mode can hide a mean-flow dependency.

Another interesting aspect of the model would be to study how the eddy-viscosity, produced naturally by the SA model, evolve with |A|. We can see from figure 6.6 that this quantity becomes weaker with an increase in |A|, meaning that the model adjusts itself in order to produce less eddy-viscosity whenever the saturation mechanism coming from the vortex-shedding mode becomes important. Note that the effect of the eddyviscosity on the mean-flow is similar to the nonlinear forcing term, coming from the vortex-shedding mode: an increase on both of them tends to reduce the re-circulation bubble. In our model, this interplay between them is somewhat fixed and given by the RANS-SA equation. However, this relation does not necessarily capture the actual physical behavior. Indeed, for all the mean-flow related quantities, the predicted value of |A| was always much higher |A| ≈ 2.3 than the DNS one |A| DN S = 1.58 or the URANS one |A| U RAN S = 1.35. It seems that the eddy-viscosity produced by our model is not sufficient to produce a short re-circulation bubble for lower values of |A|. Indeed, if we compare the eddy-viscosity in figure 6.6 (d) (for which the predicted re-circulation bubble is not too far from the URANS one) with the mean-eddy-viscosity from URANS in figure 6.2 (b), we can see that we under-predict this quantity. To investigate how a different value of eddy-viscosity would produce a different solution, we modify the model into a simpler one, where we no longer solve the SA equations. Instead, we take the eddy-viscosity to be fixed with |A| and equal to the RANS-SA one, which is much higher (see figure 6.6 (a)). The results for this model are also presented in figure 6.5 with a dashed line. We can see that the mean flow quantities are predicted in the proper range of |A| ≈ 1.6, where the mean-flow error is even smaller with this simplified model. Furthermore, the resolvent mode-related quantities are also better represented at the same range of values |A| ≈ 1.5 -1.8. This suggests that the SA equation, taking into account only mean-flow quantities, may not be a good approximation and should also include nonlinear terms on the fluctuation. Examples of such terms are the nonlinear advection U • ∇ Ñ , a nonlinear diffusion ∇ • (∂ Ñ η Ñ )∇ Ñ and the (second-order Taylor expanded) nonlinear source term, exemplified by ∂ Ñ ∂ Ñ s( Ñ ) 2 . We note that taking into account those terms is a challenging task since we do not have a proper approximation for Ñ from the mean-flow resolvent analysis. Moreover, the last two terms rely on the truncation of a Taylor expansion, that may itself represent a limited description of the whole dynamics, especially if the fluctuation is large.

As for the reconstructed fields (mean-flow and vortex-shedding mode), the results are presented in figure 6.7 for values of |A| for which the mean-flow error is minimal. We can see that the vortex-shedding modes, although they represent the vortex-emission, as in the DNS or URANS simulations, their structure can be a little "discontinuous" along some lines in the flow. This may be due to the fact that only the first harmonics is modeled and its nonlinear interactions with the second harmonics could correct locally the spatial distribution of the mode. We remark that this possible lack of the second harmonics in Iso-contours are saturated with mean-URANS eddy-viscosity color scale.

our model does not necessarily explain the over-prediction of the parameter |A| discussed before since the ratio between the first and second harmonics energies is or order ≈ 7, leading to an over-prediction of 14% for the quantity |A|, against ≈ 45% for the model.

Conclusion

This chapter dealt with the resolvent-based data-assimilation for the flow around a square cylinder at Reynolds number Re = 22000, similarly to what was done in the fourth chapter, where the vortex-shedding mode is coupled to the mean-flow equations, using a URANS turbulence model to provide the turbulent eddy-viscosity. Although this model managed to reproduce the overall features of the flow, it was less successful than the one presented in the fourth chapter in the prediction of mean-flow quantities together with unsteady ones. This comes from the fact that the URANS model did not produce the necessary amount of eddy-viscosity such that, for example, the re-circulation bubble predicted the correct one at the right range of vortex-shedding mode amplitude. Since the fully nonlinear URANS equations actually predicted better the DNS one, we believe that the model proposed is over-simplified and should take into account more nonlinear terms. However, this is a challenging task because all the remaining terms depend on fluctuations of Ñ , which is poorly reconstructed from the resolvent analysis. As a possible remedy for this could be the use of a tuning source term in the SA equation, modeling all those remaining terms, similarly to what was done in the second and third chapters of this manuscript, where some other external data is considered in the assimilation. Chapter 7

Conclusions and Perspectives

In this manuscript, we investigated some modeling techniques for data-assimilation of aerodynamic flows. Those techniques can be divided into two large groups. The first one was the recovery of the time-averaged (mean) flow through minimal modifications in the RANS (Reynolds-Averaged Navier-Stokes) equations, performed by an optimization, minimizing the error committed by the solution and some mean-flow external data. This was the subject of chapters 2 and 3. The second one relied on the mean-flow analysis, mainly based on the Resolvent formalism, subject of the subsequent chapters. For the data-assimilation purposes, we were able to build a reduced order model given by the mean-flow equation coupled with the Resolvent modes that provide a representation of the fluctuation, from which we can compute the Reynolds-stress tensor, both in a transitional and turbulent scenarios. Also, in a slightly different context from data-assimilation, we provided a framework from which both periodic and broadband components could be recovered from linear analysis. In the following we will provide a quick overview of the methods and results discussed in the manuscript, together with ideas for future work for each of them separately.

RANS-based Data-Assimilation

From the point of view of the first one (RANS-based data-assimilation), we investigated two different models. The first (model 1) consisted by the tuning parameter on the momentum equations, which is an extension of what was done by [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF] to a turbulent regime. The second one (model 2) consisted by a source term in the turbulence model (Spalart-Allmaras, in the present case), correcting it. We investigated how the choice of tuning parameter and the choice of the measure operator would change the optimized procedure and the final solution. For example, we observed that, if the whole velocity field information is provided, the model 1 manages to adapt itself to represent correctly the whole velocity-field everywhere. For the model 2 however, we did not manage to capture the whole velocity-field and the overall error decreased of one order of magnitude only. This shows that this model is not fully controlable. However, when sparse velocity field is provided, the model 1 may produce some nonphysical behavior related to the sharpness of the gradients, and further penalization on the derivatives of the forcing term had to be considered to smooth it, leading to a more physical solution. This was shown not to be the case for the model 2, where the gradients were smoother and the reconstructed final solution, although did not necessarily match the provided data, could provide a better solution when very few measurements were provided. Moreover, we observed that the model 2 produced almost the same solution independently on the number of velocity measure points. We further explained this fact with an Observability Gramian analysis, showing that this model presents a much more rigid behavior, where only a few directions in the space of forcing terms are easily reachable/observable. This rigidity was also observed when, for example, we considered wall-pressure measurements, for example. In this case, we could show that, by considering enough wall measure points, the final solution tended towards the same one produced by the velocity-field information. However, for the model 1, even with complete wall-pressure information, we did not manage to outperform the model 2, showing that, in this case, wall-pressure is not enough.

With this in mind, we can propose a few directions for future studies. Firstly, the Spalart-Allmaras model showed to be rigid when we considered a corrective forcing term and this behavior could modify the final assimilated solution. For this reason, it would be interesting to explore in a future work how other turbulence models would perform. If we consider, say, the k -ω model instead of the SA, we could observe a different behavior. For example, the model could present itself being a little less rigid than the SA, opening more room for the solution to adapt to the data, if sufficient information is available. Alternatively, even if this model ends up being as rigid as the SA, it would be interesting to evaluate if it manages to produce an assimilated solution that is closer to the reference one than the one with SA. This would mean that, even if the Observability Gramian analysis still produces a large separation between the leading Singular Values and the remaining ones (meaning that the model is still poorly observable), the leading Singular Vectors may "point" towards a better direction, on which the DNS field may lie. This would favor the choice for this model instead of SA, since it appears to be more physical.

Another interesting direction would be the consideration of turbulence models that do not rely on Boussinesq hypothesis for the data-assimilation. This was already the goal of the model 1, where a free volume force was used as tuning parameter. However, the restriction of this tuning parameter within such turbulence model would perhaps limit the action of the actual forcing term to a more physical direction, improving the assimilated solution when few (or localized) measurements are provided, always retaining the (to be verified) fully controlability. This kind of model would be useful in bluff-body-like configurations that do not satisfy the Boussinesq hypothesis and for which the volumeforcing term was not capable to correctly assimilate, for example, the wall-pressure (as in the cylinder flow presented in the third chapter). Examples of such models are the Reynolds-stress tensor models (RSM, see for example [START_REF] Hanjalić | A reynolds stress model of turbulence and its application to thin shear flows[END_REF]), where the Reynolds-stress tensor is not replaced by an eddy-viscosity and will be modeled by a supplementary set of PDEs, typically one for each component of the tensor and another one for the dissipation rate. In this case, the forcing/tuning parameter would be, for example, a volume term in each of those extra equations.

Mean-flow analysis and Data-Assimilation

As for the Resolvent-based data-assimilation (presented in the fourth chapter), inspired on the idea that the Resolvent modes can approximate the fluctuation around the mean-flow 2. Mean-flow analysis and Data-Assimilation 101 [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF], [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]), we proposed a model composed by the mean-flow equation coupled with leading Resolvent modes through the Reynoldsstress tensor. This model represents a nonlinear set of equations, whose only unknown parameters are the energies (or amplitudes) of each of those modes, one for each modeled frequency. Since the unknowns are intimately related to the spectral (frequency) content of the fluctuation, those amplitudes can be tuned by assimilating, for example, few pointwise time-resolved probes. From an experimental point of view, this represents a gain, with respect to the RANS-based data-assimilation, in the amount of information needed to recover, not only the mean-flow but also the fluctuation around it. This procedure was first applied in three different laminar/transitional configurations: squared-section cylinder and the backward-facing step forced by two incommensurable frequencies and with (time) white-noise. Generally, what was observed in those cases is that the success of the reconstructed solution of this model, together with its robustness with respect to the choice of point-wise data, can be assured if the mean-flow Resolvent modes do manage to reconstruct the main features of the fluctuation and if all of the relevant frequencies are taken into account. This tends to be the case when the reference flow does not present strong nonlinearities. An example of flows having an important amount of nonlinearities is the two-frequencies backward-facing step where, even if the external forcing contains only two frequencies, the fluctuation contains more frequencies arising from nonlinear interactions between the two fundamental ones. The presented model did not consider those nonlinear interactions, limiting the quality of the reconstruction of the final solution.

In those cases, we may employ a finer model that captures those nonlinear terms. An example of such approach is given by Symon et al. [2019a], where the fluctuation is solved in a nonlinear manner for a given (frozen) mean-flow, assimilated with similar techniques presented in the first two chapters of this manuscript. For this reason, the application of this procedure is much more challenging in our case, where the mean-flow is also an unknown of the problem. The inclusion of nonlinearities in the model of the fluctuation can thus be considered as subject of future work.

As for the turbulent case, corresponding to the square cylinder at Re = 22000, we considered first, in the fifth chapter, the mean-flow analysis (in a broad sense) for the recovery of both periodic mode (vortex-shedding) and the coherent broadband component (Kelvin-Helmholtz) that characterize this configuration and other flows around bluff bodies. We showed that those components could be isolated using one version of the triple decomposition, proposed by [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF]. We could recover the periodic component with the usual mean-flow stability analysis, where, in order to model the broadband component's effect on the quasi-periodic one, we considered an eddy-viscosity computed from the statistics of a DNS. We showed as well that the dependency of the broadband component (here exemplified by the Kelvin-Helmholtz instability) on the quasi-periodic component could be captured by the Resolvent operator, where the flow around which the equations were linearized was the quasi-periodic one, in a similar manner as in Floquet stability theory. Those modes showed themselves predictive when compared with DNS snapshots, where the dependency of the Kelvin-Helmholtz structures on the oscillation of the shear-layer is observed. For future work, it would be interesting to test this periodic Resolvent technique on a flow configuration where the frequency of the periodic motion is closer to the broadband structures one and the complete resolution of the modes with, for example, the Floquet-Fourier-Hill technique would be required.

In the sixth chapter, we applied the Resolvent-based data-assimilation, tested for lam-inar flows, in a turbulent one, consisting also in the square cylinder, same configuration as in the fifth chapter. Since this configuration is much more complex than those discussed in the fourth chapter, we proposed the use of a URANS model, from which we could use the eddy-viscosity, necessary for the prediction of the vortex-shedding mode and for the mean-flow equation saturation. The model consisted then of the mean-flow equations (momentum equations, together with the Spalart-Allmaras one) coupled through the nonlinear advection term with the resolvent mode modeling the vortex-shedding. This means that this model did not resolve the coherent broadband component, since it would leave the model much more complex. We observed that the solution of this model, although it captured the overall feature of the reference flow, did not predict well the range of amplitudes of the fluctuation of the reference flow (DNS or URANS). We could also show that the balance between this amplitude and the eddy-viscosity was not quite the one present in the URANS solution (and arguably in the DNS one). We believe the reason for that was an over-simplification of the model where, the only nonlinear term considered for the fluctuation was the nonlinear advection. The reason why the consideration of other nonlinear terms in the model is difficult is that the remaining ones are all dependent on the fluctuation of the eddy-viscosity, which is poorly represented by the mean-flow analysis when compared to the URANS one. Furthermore, the first harmonics of this quantity is of the same order of the second one, suggesting that there may be strong nonlinear interactions for this variable and a more complex dynamics for it needs to be used. One possible direction for future work could be the design of a such more refined model for this variable and its further use in the mean-flow equations through nonlinear terms. Another possible solution would be the use of of a tuning parameter such as a source term in the mean SA equation (similarly to what was done in chapters 2 and 3) that corrects this erroneous balance between the eddy-viscosity production and the vortex-shedding amplitude and its further determination with an optimization procedure. An obvious drawback of this approach would be the need for more data to be assimilated since now we tune a spatially-dependent function as well.

χ < 18.4. Moreover, it has a singular behavior for χf v2 = -S(κ 2 d 2 )/ν. A . As we will see, the phenomenon of negative r (including its sigularity) may be a symptom of other anomolous behavior. Indeed, if we write the nonlinear production term, say, s P D = s P + s D = c b1 S ν -c w1 f w (ν/d) 2 as: .16) we can see that the production term inverts its signals for r < 0. Moreover the destruction term also inverts its signal if f w < 0, which corresponds to the region -1.185 < r < 0. This phenomenon may seem a small deviation of the "normal behavior" of the model, which tends to produce solutions for which r > 0. However, it may not be the case in the convergence process, which in this work, is based on the Newton method. Indeed, the term ∂s P D /∂ ν presents an oscillatory behaviour at the regions where r < 0 (see figure A.1). This behaviour deteriorates drastically the convergence of the Newton algorithm [START_REF] Crivellini | A Spalart-Allmaras turbulence model implementation in a discontinuous galerkin solver for incompressible flow[END_REF]). To correct this behavior, we limit its value to r max such that:

s P D = νχ d 2 c b1 κ 2 r -c w1 f w (A
r * = ν Sκ 2 d 2 , r(r * ) =     
r max r * < 0, r * 0 ≤ r * < r max , r max r * ≥ r max .

(A.17)

This function is not differentiable (nor continuous) if we see it as r = r(r * ). However, what wee need to have is a differentiable r = r(ν, u) or r = r(S(κ 2 d 2 )/ν, χ), as evoked before.

On can from figure A.1 see that, below the singular line, the value of r * is negative. Thus, in those regions, r = r max . The same is true if we approach the singular line from the positive r * region. We will have an increasing value of r * , which also translates into r = r max . This means that at the originally singular line, the function r is differentiable (and constant, r = r max ). Moreover, the only source of non differentiability is when r * = r max . At this region, the model is continuous, but not rigorously differentiable. So that the jump itself on the quantity ∂s/∂ ν is: .20) This jump function depends directly on the factor df w /dr| r * =r - max . However, for r > 2 the function f w stays at a plateau, and its derivative is almost zero. This means that, for a sufficiently large value of r max , the jump is almost nonexistent. In our work, this value is set to be r max = 10, as in [START_REF] Crivellini | A Spalart-Allmaras turbulence model implementation in a discontinuous galerkin solver for incompressible flow[END_REF].

d 2 ν ∂s ∂ ν = c w1 dfw dr r * =r - max r max 1 -r max f v2 + χ df v2 dχ χ (A
With this modification, the s P D function, together with its Jacobian are much more smoother and lead to a much better Jacobian matrix (see figure A.2). We remark that the production term s P was not in our code written in function of r, as presented before. For this reason, the negative-production problem was not solved by this approach. However, it was not observed to lead any convergence problems, as documented by A.1.

Stabilization procedure -SUPG

Contrarily to finite volume, finite differences and, more recently, discontinuous Galerkin methods, the Finite Element Method (FEM) does not provide a natural framework in which one can add some 'upwind' effect, which stabilizes the numerical scheme for high Reynolds number flows. However, this effect can be achieved by modifying the variational formulation through a method called Streamline-Upwind Petrov-Galerkin (SUPG). More precisely, it was observed that the modification of the test function of a classical variational formulation in a precise way leads to stabilization of the scheme [START_REF] Brooks | Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF]). Several forms of such modifications are proposed in the literature (see Franca et al. [1992], Franca and Frey [1992]). The main idea is to modify the test function of the FEM through the advection operator, such that the scheme becomes decentralized and stable. For the RANS equations, for simplicity, we only apply this modification of the test function on the advection operators. This leads to a simpler variational formulation than the ones from, for example, Franca et al. [1992]. However, this is extensively done in the context of unsteady simulations, without loss in precision (see, for example [START_REF] Bao | Numerical prediction of aerodynamic characteristics of prismatic cylinder by finite element method with Spalart-Allmaras turbulence model[END_REF], [START_REF] Haferssas | Efficient and scalable discretization of the navier-stokes equations with LPS modeling[END_REF]). The variational formulation employed is: R ([u, p, ν], [v, q, ν]) = R F EM ([u, p, ν], [v, q, ν]) (A.21)

+ Ω k Ω k τ N S SU P G u • ∇v (u • ∇u) (A.22) + Ω k Ω k τ SA SU P G u • ∇ν (u • ∇ν) (A.23)
where R F EM is the standard, unstable, FEM discretization. The parameters τ N S SU P G and τ SA SU P G are parameters tuning the amount of artificial diffusivity and should not be taken too large, sacrificing precision, nor too low, sacrificing stability. A good compromise is given by: where the parameters ξ N S (Re N S h ) and ξ SA (Re SA h ) are limiters, each function on the mesh Reynolds numbers such that they saturate for Re h → +∞, preventing it to over-diffuse.

Here we introduce two different parameters for the NS and SA equations individually since we will study their dependency on the final solution separately. Furthermore, this parameters must tend to a constant whenever |u| → 0, such that τ SU P G tends to a finite value. Several expressions for those functions exist in the literature, all of them providing similar results. The one we will use here is given by: where the definition of the mesh Reynolds numbers is given by: .26) where here we introduce two diffusivities, one for the NS equations µ and another one for the SA equations η. The reason for this is that, in what follows, we will study how their definitions can affect the solution, either by considering their laminar version (ν) or their turbulent versions (for example, µ = ν + ν t ). Finally, the parameter h T is the local element size. When the computational mesh contains high-aspect ratio elements, it is preferable, in order to avoid over-diffusion (see [START_REF] Mittal | On the performance of high aspect ratio elements for incompressible flows[END_REF]), to take it as the smallest length of the element. It was the case for the backward-facing step, treated in the second and third chapters of this manuscript. For the square cylinder case, since the mesh is essentially build with low aspect-ratio elements, this parameter is taken to be the its largest length. However, in order to evaluate how the SUPG works and how it converges in mesh refinement, it is interesting to modify this value to be equal to γ × h t , so that the amount of artificial eddy-viscosity will be decreased with a small value of γ. In what follows, we perform several tests using the models and numerical implementations described above.

Re N S h = |u|h T 2µ , Re SA h = |u|h T 2η (A

Numerical tests and Mesh Convergence

In this section, we are going to provide some mesh-convergence on three different problems related to the RANS-SA equations, applied on the square-section cylinder, considered in the fifth chapter. The first one consists on the Newton algorithm for the convergence of the RANS-SA equations. The second one consists on the Resolvent modes, in this case the vortex shedding one. The third is the solution of the model considered in the fifth chapter for data-assimilation. All those problems will be solved on three different meshes, displayed on figure A.3. The first one ( 30000 elements) is the coarsest one, having a relatively good refinement around the cylinder. The second ( 75000 elements) one is more refined around the cylinder and on the wake up until x < 5. The third one ( 100000 elements) is essentially the same as the second one, however, with a finer refinement up to the outflow boundary condition. We remark that all those meshes correspond to the half-domain y > 0. Their use when the whole domain was needed (for example in the fifth chapter) was done by constructing a full-domain mesh constructed by uniting the original mesh with its symmetric.

RANS-SA Solution

This first problem we will consider here is the steady-flow equations, the RANS-SA ones.

Those equations are solved with the Newton method, together with an increase of the Reynolds number. The strategy here is to give ourselves a few Reynolds numbers, starting with 100, for which a solution can be easilly found, and, starting from this solution, we converge a solution for Re = 500, then for Re = 2000, Re = 5000, Re = 10000 and finally Re = 22000.

From the point of view of the final solution, in a first time, we investigate how the final solution is modified with the mesh and with the diffusivity taken into account in the SUPG stabilization. This diffusivity, here denoted µ for the Navier-Stokes equations and ζ for SA equation, can either be considered as the laminar one ν or the turbulent ones µ = ν + ν t and ζ = ν + η. Those results are summarized in table A.1. We can see .3, where the spatial distribution of the energy of the mode is evaluated. We can see that whenever the mesh is coarse and the SUPG used is the full version γ N S h = 1, the energy tends to be more localized at the vicinity of the cylinder where if the mesh is refined and/or the parameter γ N S h is decreased, the energy tends to be more distributed at the outflow region of the domain. In particular, even for the coarsest mesh, we manage already to converge, if the parameter γ N S h is sufficiently small. Table A.3: Mesh convergence for SVD modes for some SUPG implementations and mesh refinements. The energy distribution of the Resolvent mode (given by κ Ω i = Ω i |û| 2 ) is studied in function of the SUPG implementation, for Ω 1 = (-0.75, 0.75) × (-0.75, 0.75), Ω 2 = (-0.75, 5)×(-0.75, 0.75)\Ω 1 and Ω 3 = Ω\(Ω 1 ∪Ω 2

Data-Assimilation model resolution

The final test for our model would be the resolution of the model proposed in the last chapter for data-assimilation, where not only the mean-flow equation is solved with the Resolvent vortex-shedding mode. This is thus some sort of combination between the two last mesh-convergence studies. To perform this test, we solve the model for several values of |A| and evaluate the model's solution in function of it. This is performed for the three meshes and for three parameters γ N S h = 0.05/0.2/1. Those results are plotted in figure drastically. Indeed, from figure B.3 we can see that, in order for the algorithm to correct those effects of the forcing term (that may extend to the solution), it takes generally 200 iterations for a level of convergence of 10 -7 , while for the mesh 1, it takes 10-20 iterations. This poses a real problem for the applications presented in this manuscript since we cannot afford to have structured and uniform meshes. Proposed Solution: Cholesky Decomposition It is interesting to note that, at the same time where the Discrete Formalism produced the best results so far, it produces some inconsistencies in the FEM formulation. On the other hand, the Continuous formalism did not produce any result, but, if we look closely, we find that the inconsistency described does not take place, since the gradient definition does not contain the mass-matrix. To make this Continuous Formalism applicable on a BFGS algorithm, we have to perform a change the variables, on which the Euclidean norm (performed inside BFGS) will make sense both physically and numerically. Generally, this change of variables reads F = L T F . The Euclidean inner-product of two arbitrary vector is given the by:

( F1 ) T F2 = F T 1 LL T F 2 = F T 1 MF 2 = (f 1 , f 2 ) Ω (B.18)
where the last two terms correspond to the imposition that this quantity represents the

Figure 2

 2 Figure 2.1: DNS solution: streamwise mean-velocity field. Dashed lines represent negative values

Figure 2 . 2 :

 22 Figure 2.2: Sketch of computational domain, together with its boundary donditions.

Figure 2 . 5 :

 25 Figure 2.5: Convergence history J n /J 0 as a function of the iteration count n when optimizing fx (solid line) or fν (dashed line).

  Figure 2.6: Complete velocity field assimilation: optimal velocity field when optimizing (a) fx (model 1) and (b) fν (model 2). The reference (DNS) is provided in dashed lines in both plots.

  Figure 2.8: Complete velocity field assimilation: (a) wall-pressure and (b) wall-friction. Red dots correspond to the RANS-SA baseline solution, blue dots to the time-averaged DNS, continuous black line to model 1 and black dashed line to model 2.

Figure 2 .

 2 Figure 2.10: Sparse-velocity measurement assimilation: assimilated velocity field ũ (solid lines) together with reference solution u (dashed lines) for (a) fx , (b) fν and (c) fx with penalized gradients. Green dots correspond to measurement locations x i m . Wall-pressure (d) and wall-friction (e) are plotted as well. Blue dots correspond to the DNS solution, continuous black line to model 1, black dashed line to model 2 and orange line to model 1 with γ-penalisation.

Figure 2

 2 Figure 2.11: Separation of Singular-Values (λ i /λ 0 ) 2 for the fx model 1 (black dots) and the fν model 2 (white dots), computed around the RANS-SA solution.

Figure 3

 3 Figure 3.1: (a) Unstable steady-state (longitudinal velocity) solution and also the starting point of our optimization procedure and (b) mean-flow and exact divergence Reynoldsstress tensor f x , f y (c,d)

  Figure 3.2: Results for pressure on the wall optimization: (a) convergence of the cost functional, (b) wall-pressure and (c) wall-friction in function of the angular coordinate of the cylinder's surface, starting at the stagnation point (θ = 0).

  Figure 3.3: Reconstructed velocity field (a) and Reynolds-stress tensor force (b) for wallpressure assimilation.

  Figure 3.5: Reconstructed velocity field (a) and Reynolds-stress tensor force (b) for wallfriction assimilation.

  Figure 3.6: Wall-friction assimilation: wall-pressure (b) and wall-friction (c) for model 1 (black curve), model 2 (black dots) and γ = 1 penalized model 1 (orange curve). DNS results present in blue dots. The color-code will be retained in all the plots in the following

  Figure 3.7: Wall-friction assimilation: u x iso-contours for model 1 (a), model 2 (b) and γ = 1 penalized model 1 (c).

  Figure 3.8: Wall-friction assimilation: profiles of u x (y) velocity at x = 3.5 (a), x = 4 (b) and x = 4.5 (c). Stream-wise component of the optimal forcing for (d) non-penalized case and (e) penalized case. Black lines represent two streamlines.

  Figure 3.11: Wall-pressure assimilation: profiles of u x (y) velocity at x = 3.5 (a), x = 4 (b) and x = 4.5 (c).

  Figure 3.12: Wall-pressure assimilation for three numbers of probes N w = 5, 10 and 20: wall pressure (a,c,e) and wall-friction (b,d,f). Blue points correspond to the DNS, solid lines correspond to model 1, dashed black lines to the model 2 and orange lines to the model 1 with fx -gradient penalization. The location of pressure measurements are indicated by gray lines.

Figure 4 . 1 :

 41 Figure 4.1: Sketch of the physical domain and inflow boundary condition for the square cylinder case. The re-circulation bubble's length is indicated in the figure and will be defined as L.

  Figure 4.2: Square cylinder flow at Re = 60 (a,b,c) and Re = 100 (d,e,f): (a,d) measured signal corresponding to the cross-velocity u y at (x, y) = (2, 0), (b,e) Discrete Fourier Transform of signal in (a,d), giving the modulus of the Fourier coefficient (|ǔ y |(ω 0 ) = 0.09/0.21) at the fundamental frequency ω 0 = 0.77/0.91 and (c,f) comparison of this quantity with the one given by model (4.15) for several values of |A|, the optimal one being |A| ≈ 1.71/2.23.

  Figure 4.3: Comparison between reconstructed unsteady flows (left column) and DNS (right column) for the square-cylinder configuration at Re = 60. Stream-wise component of (a,b) mean flow velocities (c,d) first-harmonic velocities Aû 0 and (e,f) mean force induced by the Reynolds stress tensor f = 2|A| 2 Re{û 0 • ∇û * 0 } with |A| = 1.71. The black curves delimit the re-circulation regions of the reconstructed mean flow. The frequency of the periodic flows is respectively ω 0 = 0.77 and the amplitude of the Fourier mode at the fundamental frequency is |A DN S | = 1.68.

  Figure 4.7: (ω, |A|) parametric study : (a,b) global error (minimal value marked by a black dot), (c,d) bubble's length and (e,f) mean drag. Results presented for Re = 60, 100, respectively. Reference (DNS) measured quantities marked in thick line and reference frequency and first-harmonic's amplitude (ω, |A|) represented by the red dot.

Figure 4 . 8 :

 48 Figure 4.8: Sketch of the physical domain and inflow boundary condition for the backwardfacing step case. The length of the main re-circulation bubble, the detachment and reattachment points are indicated in the figure and defined to be L, x d and x r , respectively.

  Figure 4.9: (a) Exact time-averaged flow (with L = 7.5, x d = 5.8 and x r = 10.4), (b) Fourier mode at ω = 0.48, (c) at ω = 0.70, (d) at ω = 0.70 -0.48 = 0.22 and spectrum of the cross-velocity-spectra field at (e) (x, y) = (1, 0.2), at (f) (x, y) = (4, 0.2) and (g) (x, y) = (15, 0.2).

  Figure 4.11: Reconstructed fields with optimal parameters (|A 1 |, |A 2 |) = (0.77, 0.44) minimizing ē : reconstructed (a) and exact (b) mean flows, reconstructed (normalized) Fourier modes at ω 1 = 0.48 (c) and (e) ω 2 = 0.7 and reference ones (d,f). Reynolds stress tensor u x u y for the model (g) and DNS (h). Identical colour-bars have been used for reconstructed and DNS fields to help comparison.

  Figure 4.13: Power Spectral Density (PSD) for a broadband flow (a). Its integral on five (N = 5) frequency intervals are compared with the discrete-in-frequency energies coming from the model (b).

Figure 4 .

 4 Figure 4.14: Spatial distribution of the energy, evaluated on the line y = 0.2, of the resolvent modes on the mean flow (solid lines) and Fourier modes (dashed lines) at several frequencies (both modes are normalized).

  |A j | coefficients for N = 5 by computing the Resolvent modes based on the mean-flow determined with N = 2 and by setting |A j | = I xm,ω j /2π/|û j 0 (x m )|, where x m is the measurement point. We repeat the same procedure to initialize the optimization for N = 10 based on the mean-flow obtained for N = 5.

Figure 4 .

 4 Figure 4.16: Comparison between reconstructed (a,c,e,f,i,k) and reference (b,d,f,h,j,l) fields: mean-flow (a,b) and (normalized) Fourier modes at ω = 0.2 (c,d), ω = 0.6 (e,f), ω = 1.0 (g,h), ω = 1.4 (i,j), ω = 1.8 (h,l).

Chapter 5 .

 5 Figure 5.1: Sketch of the frequency representation of a flow's signal containing the vortexshedding and the Kelvin-Helmholtz (coherent contributions) and the remaining (uncoherent) fluctuations (a). The time-signal from which this sketch was produced is given in (b) and corresponds to a point on the oscillating shear-layer.

Chapter 5 .

 5 Figure 5.2: Proposed scenario for a decomposition and modeling of the flow field.

  Figure 5.3: DNS: FFT of signal for the (a) cross-wise velocity component at (0, 2.5) (second green dot in figure 5.4) and (b) Power-Spectral Density for the stream-wise velocity on the point (-0.45, 0.63) (first green dot)

  Figure 5.5: Mean-eddy-viscosity computation: (a) production term (numerator in equation 5.15, prod = ∇ s u : (-u ⊗ u +2/3κI), (b) destruction term (denominator in equation 5.15), dest = ∇ s u : ∇ s u and (c,d) computed eddy viscosity ν t with penalisation terms equal to ε = 0.01/10, respectively. The eddy-viscosity computed with ε = 0.01 (c) is saturated with the maximal value of the eddy-viscosity computed with ε = 10 (d) and completed with iso-contours (dashed white lines).

Chapter 5 .

 5 Figure 5.6: Stability analysis: (a,c,e) eigen-spectra and (b,d,f) unstable mode for laminar viscosity approach (a,b), frozen eddy-viscosity (from figure 5.5 (d)) approach(c,d) and variable eddy-viscosity approach (e,f), where νt is modeled with the Spalart-Allmaras. Red dots correspond to the unstable mode computed from the ε = 0.01 mean-eddy-viscosity.

  Figure 5.8: Pressure p = p -p snapshots (a,d) at two different phases of the baseline period, corresponding to τ /T = 0 and 1/2 (same phase as in figure 5.7 (a,f,k) and (c,h,m)). The pressure of the modes p at those phases are also plotted for ω = 10 (b,e) and ω = 20 (c,f). Streamlines correspond to u . Steady-flow modes for ω = 10 and ω = 20 are also provided in (g,h) for comparison.

  Figure 6.1: RANS solution: U x velocity field (a) and eddy-visocisty Ñ /ν (b).

  Figure 6.3: Resolvent Analysis on (URANS) mean-flow: optimal (thick lines) and suboptimal gains (dashed lines) and leading mode at ω = 0.83 (grey vertical line) for (a,b) frozen ν t and (c,d) variable ν t .

  Figure 6.4: (a) resolvent mode for N (variable Ñ analysis). We compare it with the URANS first (b) and second (c) harmonics of Ñ .

  Figure 6.5: Solution of the model in function of |A|. Mean re-circulation bubble length (a), global error of velocity field with respect to the mean-DNS (b) and mean-URANS (c), amplitude of Fourier mode of the cross-velocity signals at (1.5, 0) (c) (2.5, 0) (d) and (3.5, 0) (e). Black continuous curves correspond to the model given by (6.7) and dashed ones to the comparative model where the ν t is frozen to be the RANS-SA one. DNS value of |A| ≈ 1.58 on blue vertical line and URANS value of |A| ≈ 1.35 on red vertical lines. Exact values of each appropriate quantities are given in dashed horizontal lines for DNS (blue) and URANS (red).

  Figure 6.6: Results for reconstruction procedure: evolution of eddy-viscosity ν t /ν for increasing values of |A|, |A| = 0 (a), |A| = 1.4 (b), |A| = 2.0 (c), |A| = 2.4 (d) . Iso-contours are saturated with mean-URANS eddy-viscosity color scale.

  Figure 6.7: Results for reconstruction procedure: mean-flow (a) and normalized vortexshedding mode (b) for the model given by (6.7). The same is presented for the model where the eddy-viscosity is fixed to be the RANS-SA one (c,d).

Figure A. 1 :

 1 Figure A.1: original SA parameters: r, s P D and its Jacobian ∂s P D /∂ ν as function of the two non-domensional groups χ and S(κd) 2 /ν. The dashed lines on the plots represent the set of parameters for which r is singular. Figure extracted from Crivellini et al. [2013].

  can calculate how much is the jump in ∂s/∂ ν at r * = r max :∂s ∂ ν r * =r + max

Figure

  Figure A.2: r-Modified SA model: s P D and its Jacobian ∂s P D /∂ ν. Figure extracted from Crivellini et al. [2013].

  Figure A.3: Three meshes used for numerical tests. Mesh 1 (a), 2 (b), 3 (c).

  Figure B.2: Discrete Formalism Optimization: mesh (a,d), optimal solution u * (b,e) and optimal forcing term f * (c,f).

Figure B. 3 :

 3 Figure B.3: Discrete Formalism Optimization: Convergence of the cost functional J n /J 0 .

Table 2

 2

	.1: Comparison of separation and reattachment points for RANS, LES, and DNS
	simulations for the smooth backward-facing step, Re = 28275
		(ν/ν) ∞ x s	x r
		3	0.50 6.67
	RANS -SA	3 • 10 -6 0.46 6.87
	(present studies)	0	0.40 7.35
	DNS Dandois et al. [2007]	-	0.53 3.93
	LES Dandois et al. [2007]	-	0.51 3.96

  which are provided in table 2.2. Indeed, for several values of N x , N y , the optimization based on model 2 leads to very similar values for e Ω . whereas for model 1, those values become smaller if more measurements are considered. From table

		10 0.037	0.031	0.037 ∼ 10 -6	0.3 %	8.7 %
		N y = 20 0.030	0.028	0.037 ∼ 10 -6	0.4 %	9.0 %
		N y = 5	0.042	0.035	0.037 ∼ 10 -6	0.1 %	11.9 %
	N x = 6	N y = 10 0.032	0.030	0.037 ∼ 10 -6	0.5 %	9.4 %
		N y = 20 0.027	0.027	0.038 ∼ 10 -6	1.1 %	13.4 %
	RANS-SA			0.094		-

Table 3

 3

.1: Global results of wall-friction data-assimilation e

  1 f ν

	Optimal Solution 0.060	0.065	0.080 0.6%	1.2%	13%
	RANS-SA	0.094		-

reconstructs the wall-friction (figure

3

.6 (b)), producing however a small overestimation of the wall-pressure. Looking close to what goes on near the wall, we see that the boundary layer profiles for this model (black solid lines in figure

Table 3 .

 3

2: Global results of wall-pressure data-assimilation e

Table 3

 3

.3: Global results of sparse wall-pressure/velocity data-assimilation for model 1 ( fx ) with penalization γ 2 = 1 e Ω

  22)We can now investigate which value of |A| provides a solution for which the quantity |A||û 0 (x m )| is closest to the measurement |ǔ y (ω 0 )|. To do so, we have represented in figure4.2 (c,f), the measurement of the model |A||û 0 (x m )| as a function of the parameter |A|. The dashed horizontal line indicates the value linked to the measurement. We can see that the prediction of the amplitude by the model increases monotonically with |A| and that it crosses the measurement value for a particular value of parameter |A|: |A| = 1.71

(resp. 2.23)

for Re = 60 (resp. 100). These values slightly overestimate the actual DNS values, since |A DN S | = 1.68 (resp. 2.13) for Re = 60 (resp. 100

Table 4 .

 4 2: Mesh dependency of steady and unsteady simulations for backward-facing step (|ADN S

	1	| and |A DN S 2	| represent the amplitudes of the Fourier modes ǔ(ω 1 = 0.48) and
	ǔ(ω 2 = 0.7) in the case of 2-frequencies external noise) -chosen mesh for the computations
	in boldface
			Steady Solution		DNS -2 frequencies
		# triangles	L	x d	x r	L	x d	x r	|A 1 |	|A 2 |
	Mesh 1	35000	10.875 8.688 17.508 7.455 5.765 10.390 0.448 0.282
	Mesh 2	17000	10.870 8.685 17.508 7.461 5.761 10.405 0.461 0.292
	mean-flow. In the present study, we pick:

Table 4 .

 4 3: Reconstructed parameters in function of number of frequencies considered N .

		L	x d	x r	e = ||ũ -u|| 2 A = 1 T	T 0 Ω |u | 2
	DNS	7.07	5.62	9.99	-	1.07

Table A .

 A 2: RANS-SA solutions for different implementations of SUPG stabilisation. The mesh element size is masked with a parameter γ N S

					h	and γ SA h	in the SUPG, tuning the
	diffusivity of the method.
	Mesh (# triangles) γ N S h	γ SA h	L	||ν||		Ω xν/ Ω ν # iterations
		1	1	3.800 0.213	8.497	5
	1 (29384)	0.2	1	3.876 0.225	8.496	8
		0.05	1	3.891 0.227	8.495	8
		0.05 0.05 3.898 0.231	8.516	12
		1	1	3.832 0.221	8.502	4
	2 (74563)	0.2	1	3.886 0.228	8.505	8
		0.05	1	3.896 0.230	8.505	19
		0.05 0.05 3.900 0.231	8.514	10
		1	1	3.832 0.221	8.501	7
	3 (99922)	0.2	1	3.886 0.228	8.503	7
		0.05	1	3.896 0.230	8.504	7
		0.05 0.05 3.900 0.231	8.513	10

it converges. This results are summarized in table

A

  ). Results presented for ζ = ν +ν t

		γ N S h	= 1/0.2/0.05
	Mesh (# triangles)	κ Ω 1	κ Ω 2	κ Ω 3
	1 (29384)	2.8/2.0/1.9 % 33.0/27.2/26.8 % 64.2/70.8/71.2 %
	2 (74563)	1.6/1.6/1.5 % 26.5/24.4/24.2 % 71.8/74.0/74.3 %
	3 (99922)	1.8/1.7/1.7 % 26.2/24.8/24.6 % 72.0/73.5/73.7 %

Chapter 2. Data-Assimilation for minimal Corrections of RANS equations

Results

Wall-measure for turbulent RANS

Chapter 3. Wall-measure assimilation for RANS equations

Theory

Noise Amplifier Flows -Backward-facing step

Triple Decomposition and Mean-Flow analysis

Chapter 5. Mean-flow analysis for the mixed periodic/broadband fluctuations on turbulent flows

Chapter 6. Resolvent-based Data-Assimilation for the turbulent flow around a squared-section cylinder

Acknowledgements

Appendix A

RANS -Spalart-Allmaras Numerical Resolution 1 Incompressible RANS

In this section we are interested in giving some details on the (incompressible) RANS equations, coupled with the Spalart-Allmaras model, together with its resolution in a Finite-Element code, FreeFem++ (see [START_REF] Hecht | New development in freefem++[END_REF]). They are given by:

The last equation in A.1 is a nonlinear Advection-Diffusion-Reaction equation for the auxiliary quantity ν, such that ν t = ν t (ν). In the next subsection, we introduce the original Spalart-Allmaras model such that it was first proposed [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF]). Then, we will comment on its modifications, in order to make it more stable and differentiable.

Original model

The dependency of ν t with respect to ν is given by:

Moreover, the nonlinear term s(ν, ∇ν, ∇u) is composed by a production, cross-diffusion and destruction terms, respectively given by:

each of those terms correspond to nonlinear functions on the state, which are further closed with the extra auxiliary parameters:

) The original model may have some issues when one tries to derive the Jacobian operator corresponding to the nonlinear model, such as non differentiable terms or numerical problems related to the ill-posed structure of the resulting matrix. In the next subsections, we will document all the modifications made and a brief explanation of them.

Modification in S

The simplest modification made was to replace the term S = |∇×u| by S = |∇ × u| 2 + M 2 -M , where M is a user defined parameter (typically M = 10 -4 ). This new definition of S makes it differentiable, retaining the property S > 0 and S = 0 iif ∇ × u = 0.

Negative ν modifications

In the original model, it is possible to observe some undesirable phenomena that may not occur in a continuous framework, but do occur in a discrete framework. One of them is the apparition of negative values of ν [START_REF] Crivellini | A Spalart-Allmaras turbulence model implementation in a discontinuous galerkin solver for incompressible flow[END_REF], [START_REF] Oliver | A High-Order and Adaptive and Discontinuous Galerkin Finite Element Method for the Reynolds-Averaged Navier-Stokes Equations[END_REF]) at the thin layer that transitions from the turbulent to the free-stream region. This phenomenon, coupled with a high-order representation may lead to highly oscillatory and unstable behavior [START_REF] Persson | Numerical investigation of the flow over an axisymmetric hill using LES, DES, and RANS[END_REF], [START_REF] Oliver | An unsteady adaptation algorithm for discontinuous galerkin discretizations of the RANS equations[END_REF]).

A differentiable way to overcome this problem was given by [START_REF] Oliver | A High-Order and Adaptive and Discontinuous Galerkin Finite Element Method for the Reynolds-Averaged Navier-Stokes Equations[END_REF], in the context of Discontinuous Galerkin Method discretization. His idea was to modify the original ν equation in such way that the energy of the turbulent variable ν2 /2 to reduce when ν in negative. For this purpose, let us derive the equation for this quantity:

Let us define our domain Ω and let us divide it into two parts, one of them containing positive Ω + and negative Ω -values of ν. If we assume that ∂Ω ∩ ∂Ω -= ∅, meaning that the negative values of ν only happen at the interior of Ω, then:

The second and third terms are null since, by continuity, ν| ∂Ω -= 0. By opening the source term s into its three terms, we have:

Since the convective fluxes are null, the above equation represents the "energy" of ν in the region where it is negative in a control volume Ω -moving with the flow. This means 1. Incompressible RANS 105 that, for the energy to be positive if ν < 0, we have to ensure the right-hand-side to be negative. For that matter, we start by choosing a modification on η:

ensuring its positivity and ensuring the first integral to be negative if ν < 0. The sign of term νs P = c b1 S ν2 is only dictated by the sign of S, that is positive (maybe unless 1 < χ < 18.4). For this reason, the following modification is used: .10) where: .11) We remark that this modification does not always ensure νs P < 0 when ν < 0, since g n has the "wrong" sign if -1/999 < χ < 0. However, this small interval is unavoidable if the differentiability of the model has to be retained, since s P | χ=0 = 0 and ∂ χ s P | χ=0 = c b1 Sν ≥ 0. Finally, the last term νs D = c w1 f w ν3 /d 2 is negative ("wrong" sign) if f w < 0, which happens for -1.185 < r < 0. Therefore, another modification is made:

Rewriting the whole nonlinear forcing term:

In addition, inspired by the η modification and from the fact that negative diffusivity coefficients may lead to an ill posed behavior, we also modify the turbulent viscosity ν t as: .14) This modification seems to be contrary to the general purpose of differentiability of the model, however it is not the case, since the derivative of f v1 ν approaches 0 when χ → 0.

Negative r modifications

We now discuss another issue with the original SA model, which is the possibility of having negative values of the originally positive definite quantity r, being related to a mixing length l such as r = l 2 /(κ 2 d 2 ) > 0. This phenomenon can be better understood by rewriting this function as:

We see that the function r is function of the two non dimensional groups χ, S(κ 2 d 2 )/ν, the first being responsible for negative r, since the function f v2 can be negative for 1.003 < that the global quantities we chose to look at (namely the recirculation bubble length L, the amount of eddy-viscosity produced ||ν|| and its baricenter Ω xν/ Ω ν) do not depend much on the choice of the diffusivity. Moreover, the number of iterations needed for converging from Re = 10000 to Re = 22000 is almost independent on this choice. We can see that the choice on the mesh seems to be more important than the choice on the diffusivity. For this reason, if explicitly said otherwise, we will retain the definitions µ = ν + ν t and ζ = ν. We now study how the effect of the extra diffusivity given by the version of SUPG implemented as function of the mesh. To do so, we mask the mesh-size h that appears in the computation of the SUPG functions by the factors γ N S h (in NS equations) and γ SA h (in SA equation). By doing so, the SUPG functions will "see" the element sizes smaller than they are, minimizing their effect and thus the numerical diffusion. We can see from table A.2 that, by taking smaller values of those γ parameters, the solution converges towards a longer bubble. It is interesting to note that even the coarser mesh manages to recover the same solution for the finer one, requiring only a sufficient small value for the γ parameters. Moreover, the number of iterations that the Newton method takes, say for the mesh 1 (coarser one), is similar to the one for the mesh 3 (finer one). This shows that, even if we decrease the numerical diffusivity, the convergence is not altered. We note that this is not true for the mesh 2 (intermediate one), where the discretization towards the outflow border is very coarse. This is related to some oscillations of the variable ν at the outflow region where it transitions from the farfield value to the turbulent one inside the bubble. Apparently, this is caused by a negative ν that persists even with the corrections mentioned in previous paragraphs. This effect is not fully understood and will be avoided by considering meshes sufficiently refined at that region and more refined towards the outflow.

Resolvent mode analysis

The Resolvent analysis (or the Stability analysis) is a much more different test case from the RANS-SA solution, since it has a connection with the Jacobian matrix and not with the residual itself and the SUPG stabilization may reflect differently for this problem. To assess this difference, we consider only the vortex-shedding mode, corresponding to the Resolvent mode on the RANS-SA solution, at frozen ν, at ω = 0.83 and investigate how A.4. We can see from those results that indeed the model's solution converges either with a mesh refinement or with a decrease of γ N S h , as observed before. Furthermore, the mesh chosen for this kind of model, presented in the fourth chapter, is the coarser one, since no observed quantity seemed to depend on the discretization.

Appendix B BFGS Algorithm and its use with a Finite-Element Solver 1 BFGS and Line search Methods

With the information of the gradient ∇ f J of our cost functional J with respect to the vector of control variables (here denoted f ), several optimization algorithms can be used to find the optimal conditions. They all rely on the building of a sequence of forcing terms f n such that f n → f * , where f * is the optimal forcing term, for which, ideally, ∇ f J | f * = 0. A simple gradient-based method (employed, for example, in [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean ow reconstruction[END_REF], [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF], etc) creates this sequence such that f n+1 = f n -ε n ∇ f J , where the parameter ε n is determined, for example, with a line-search algorithm. However, this is a first-order method since it relies solely on the information of the gradient. A well known example of second-order method is the Newton method for optimization, that constructs this sequence of forcing terms such that f n+1 -

, where H n = ∇ f ∇ f J is the Hessian matrix at iteration step n. This method is generally not so easy to apply since one does not, in general, have the Hessian matrix. However, some methods try to approximate it, retaining a good memory requirement (see [START_REF] Fletcher | A rapidly convergent descent methods for minimization[END_REF], [START_REF] Broyden | Quasi-Newton methods and their application to function minimization[END_REF], [START_REF] Fletcher | A new approach to variable metric algorithms[END_REF], [START_REF] Nocedal | Numerical optimization[END_REF], [START_REF] Bonnans | Numerical optimization: theoretical and practical aspects[END_REF]). In this context, we propose to use the limited-memory Broyden-Fletcher-Goldfarb-Shanno (l-BFGS) method. This method approximates the Hessian matrix from previous iterations computations in the following way (with H 0 = I):

where

is the difference of the gradient between two successive iterations and s n = f n+1 -f n = ε n p n (p n being the descent direction, p n = H -1 n ) the difference of control parameters. Since at every successive iteration this matrix is incremented by two sets of line/column, we can use the Woodbury identity to invert it:

Where, in order to respect some memory requirement, the storage of previous vectors s n , y n is limited to a certain number k of previous iterations. Together with the estimation of Appendix B. BFGS Algorithm and its use with a Finite-Element Solver the Hessian matrix, the method is coupled with a line search algorithm. This algorithm finds the parameter ε n such that the following (Wolfe) conditions are satisfied:

This roughly summarizes how the BFGS algorithm works. In the next section, we will discuss how to use it with a Partial Differential Equations (PDE) solver, in particular with a Finite-Element solver.

BFGS with Finite-Element Solver 2.1 A toy model for EDP solver -Laplace equation

In what follows, we will be interested in applying the optimization procedure on a very simple PDE, namely the Laplace equation (with homogeneous Cauchy boundary conditions):

We remark however that the presentation of the discretization of those equation, together with the optimization procedure are general and valid as well for nonlinear model, such as the ones considered in this manuscript. The discrete version of the continuous problem (B.4) is given here by the Finite-Element Method (FEM). This method arises from the variational formulation associated with the strong equations (B.4). To derive it, we multiply it by an arbitrary test function, integrate by parts and apply the boundary conditions, leading to:

In FEM, we give ourselves a set of base functions {φ i } that will be used to span our solution and the forcing term. We can say that:

where U i , F i are the coefficients that represent numerically the function u(x) and f (x).

Since the variational formulation has to be true for any v, we may choose v(x) = φ j (x), ∀j.

Each such choice of test function will lead to a line of the linear system:

where

is the so-called rigidity matrix and M i,j = Ω φ i φ j dΩ is the so-called mass matrix. We remark that the inner product associated with the continuous problem is given in its discrete version by the mass matrix:

Once the discretization of the continuous equations are set, we proceed to the associated optimization problem.

Optimization Problem

The optimization problem will consist of finding the optimal forcing term f such that the cost functional:

is minimal, where the relation between u and f is given by the equation (B.4). Here the reference state u r will consist on the solution of this Laplace equation imposing f r = tanh(10(x -0.5)) * sin(πx) * sin(πy), an arbitrary function. This problem can be viewed in two different manners: a continuous one and a discrete one.

Continuous Formalism:

The continuous version is essentially the one presented in the main body of this manuscript in chapters 2 and 3 and consists of the Lagrangian multiplier formalism posed for the continuous equations, where the adjoint system and the gradient are set. The Lagrangian of this system is given by:

By deriving this functional with respect to the state, adjoint and forcing term, we have:

The discrete version of those equations, necessary for numerical implementation, is given by:

where G c is the continuous gradient in this framework.

Discrete Formalism:

The discrete formalism arises when we treat the optimization problem as the minimization of the cost functional:

whose inputs are the coordinates of the vector U and its relation with F is given by the discrete system KU = MF . If we apply again the Lagrangian multiplier formalism, given by the Lagrangian:

If we regard this functional as being dependent on each component of the state U i , the adjoint U † i and the forcing term F i , we can derive it with respect to each of those components, leading to:

We notice that the discrete gradient G d differs from the continuous one by the factor of the mass-matrix. This can be understood by noticing that the discrete gradient should quantify the derivative of the cost functional with respect to the coefficients of the forcing vector F , whereas the continuous one informs how the cost functional varies with respect to an arbitrary variation of the forcing term f . Indeed, if we give ourselves a perturbation in the forcing term δf = i δF i φ i , a variation of the cost functional δJ will be given by:

In what follows, we will discuss the implications of the choise between those approaches on the numerical results.

Optimization Results

The goal of this section is to, through the optimization procedure, find the pair (u r , f r ), displayed in figure B.1. At first, we will naively apply the continuous and discrete formalism by providing BFGS either (J , G c ) or (J , G d ).

Continuous Formalism: By applying this formalism we were not able to converge to any solution. This was the case not only in this simple Laplace equation model, but also in any other PDE-constrained optimization problem considered in this manuscript. In fact, the algorithm was not able to exit the very first global iteration, even if the cost functional clearly decreased in the line-search procedure. Typically, the line-search method kept decreasing the parameter ε until the point it became small enough for the optimization procedure to stop. We are lead to believe that the main reason for that is a norm-inconsistency between our PDE-based model and the internal BFGS routines. Indeed, as mentioned in section 1, we can see that all the norms and inner products performed by the BFGS algorithm are based on the euclidean norm U T V = i U i V i , which does not have any physical or numerical meaning, if we regard U, V as finiteelement representations of their original functions. For this reason, we find necessary to test the discrete formalism.

Discrete Formalism: Contrarily to the continuous formalism, the discrete one always produces an optimal condition. The optimal solutions u * (figure B.2 (b,e)) correspond very well to the reference one. However, we can see that, in function of the mesh, typically if the mesh is structured or not, the forcing term may become very noisy. To understand this, we need to take a closer look to how the gradient information G d , provided by our FEM solver, will be re-fed to it in the form of a new forcing term F . Indeed, for the sake of clarity, if we place ourselves at the first global iteration, for which F 1 -F 0 = ε 0 G d,0 , we can see that the solution of the next iteration will be given by the system:

we can see that on the right-hand-side of this equation we have a term with a double mass-matrix product, acting on a FEM representation of the adjoint function U † 0 . This procedure produces a forcing term that is inconsistent with the FEM. If the mesh we consider is somewhat structured and uniform (as it is the case of the mesh 1), this inconsistency may produce some negligible effects, where the mass-matrix itself will be very well structured and will not establish locally any preferential direction. However, if the mesh is unstructured (mesh 2), this effect is quite visible. Furthermore, in those cases, the number of iterations for the convergence of the optimization algorithm may increase energy inner product of the two FEM functions f 1 , f 2 . We can see that the matrix L has to be given by the Cholesky decomposition of the mass-matrix: M = LL T . For this reason, the continuous gradient, written now in terms of this new variable becomes:

We note that, with this changes of variables, not only all the norms performed within BFGS will correspond to the physical (energy norm) ones, but also the FEM formulation will be consistent, since the forcing term on the system B.17 will be:

Indeed, with this kind of approach, we obtain better results (figure B.4). The are meshindependent and manage to reconstruct both the forcing term and the solution itself. Furthermore, the convergence of the cost functional B.5 is also (almost) mesh-independent. We notice as well that, for the model 1, where the tuning parameter is the volume force fx , the incompressibility of this forcing term, discussed in the second chapter, is not obvious since the Hessian matrix may modify the spatial structure of the gradients, leading to an increment in fx that is not solenoidal. However, this is not the case as we will see here. To be able to prove this, we need to show that the application of H -1 n+1 onto any divergencefree vector z (here, the gradient) is also divergence-free. To do so, we use induction, where we suppose it is true at some iteration n, i.e., DH -1 n z = 0 (true for n = 0 since H 0 = I), where D is the discrete-version of the divergence operator. If we suppose that Ds n = 0, we have that:

which is zero since both terms on the last equation are applications of DH -1 n onto divergence-free vectors.

We notice that the fact we used a linear transformation coming from the Cholesky decomposition does not invalidate this argument since the vector z can be considered as being already transformed by it and the matrix D modified to take into account for the Cholesky matrix L incorporated in z.

Adjoint Discretization in FEM

In this section, we briefly discuss the adjoint problem for a (sufficiently) general nonlinear Partial Differential Equation (PDE) and its further discretization in FEM. The main goal here is to prove that the transpose of the Jacobian matrix provides indeed possible discretization of the adjoint problem. For our purpose here, a nonlinear scalar advectiondiffusion equation under the form:

Here f is a forcing term similar to those used as tuning parameters in the body of this work. The boundary conditions associated with this problem are typically of Dirichlet or Neumann types:

The Finite-Element discretization of such equation considers its weak-form, given by:

This (semi-linear) form provides a nonlinear residual that needs to be zero for all test functions v so that u is a weak solution of our problem. In practice, it is achieved with a Newton method, which converges towards a solution by successive inversions of the Jacobian of the system, which corresponds to the linearization of this residual:

where a(δu, v) is a bilinear form that leads to the Jacobian matrix upon discretization.

Here, for simplicity, we suppose that the solution guess u already satisfies the Dirichlet boundary condition, meaning that the correction δu has to satisfy δu = 0 at Γ D . This imposition of the Dirichlet boundary conditions is done through, for example, elimination of columns and lines of the complete matrix corresponding to the degrees of freedom on which we impose u = 0. Another way, which is default in FreeFEM++ and used throughout this work, is the penalization of those degrees of freedom in the diagonal of the matrix, freezing them.

We proceed now to the adjoint problem. This problem appears whenever we desire to optimize some functional J under the PDE constraint. We consider the Lagrangian:

The variation of this Lagrangian with respect to u, and the further integration by parts of the viscous terms lead to:

where, at this point we can identify the bilinear form a † (u † , δu) = a(δu, u † ) in the formulation of the adjoint equations. Since this form is the same as the one leading to the Jacobian, only with changed arguments, we can see that the bilinear form a † (u † , δu) will be represented numerically by the transpose of the Jacobian matrix. However, to be able to fully justify the use of the transposed matrix as the adjoint one, we need to prove two things: (a) that the last boundary term appearing in (B.27) is zero and (b) that the adjoint variable should satisfy a Dirichlet boundary condition at Γ D . This last point presents itself numerically in the imposition of the boundary condition by line/columns elimination or penalization, and we need to ensure that the positions of the adjoint matrix where this procedure was done are the same as the ones in the Jacobian matrix. We remark that, if an homogeneous Dirichlet condition is verified at Γ D , the aforementioned boundary term is only non null at Γ N , where, by application of the linearized Neumann boundary conditions for δu, i.e., ν(u)n • ∇(δu) + ∂ u ν(u)(δu)n • ∇u = 0, we prove it is also null. This means that both points (a) and (b) can be assured by showing that u † = 0 at 4. Adjoint Discretization in FEM 121 Γ D . To do so, we continue the integration by parts in the adjoint derivation:

This determines the volume (strong) equation for the adjoint problem as:

and the boundary conditions. Indeed, at Γ N , by applying the linearized Neumann boundary condition, we have that, to the boundary terms to be zero, n • (au † + ν(u)∇u † ) = 0. At Γ D , since δu = 0, to the remaining boundary terms to be zero, we have that u † + 0. We have proven this way that the transpose of the Jacobian matrix provides indeed a FEM discretization of the adjoint equations. Abstract : In a first moment we will be interested in the recovery of the mean-flow quantities from partial or sparse information, ranging from point-wise velocity probes to wall-pressure and friction. This will be achieved by considering the Reynolds-Averaged Navier-Stokes (RANS) equations, completed with a model, here the Spalart-Allmaras. This kind of modeling has been conceived for some benchmark flow configurations and may lack generality, leading to erroneous predictions, especially when re-circulation is present. We modify this model with a tuning parameter such that its solution matches the best the aforementioned mean-flow data. The configuration considered was a Backward-Facing Step at Re = 28275, with actual data stemming from a DNS. Then, we turn our attention to linear mean-flow analysis and its use to predict the nonlinear unsteady fluctuation. In particular, we design a reduced-order model, composed by the mean-flow equation coupled with the resolvent modes, predicting the fluctuation for each existing frequency. The energies of those modes are used as tuning parameters for the dataassimilation procedure, that will take as input typically (very) few point-wise time-resolved information. This technique will be applied in transitional flows such as the one around a squared-section cylinder, a benchmark case for oscillator flows, and a backward-facingstep, a typical noise-amplifier flow. We then consider a turbulent case corresponding to the flow around a squared-section cylinder at Re = 22000, having both oscillator (periodic vortexshedding) and noise-amplifier-like characteristics (represented by the Kelvin-Helmholtz structures). Classical mean-flow stability analysis is used to recover the the vortex-shedding mode and a resolvent technique, based on the linearized equations around the periodic component, is used to recover the dependency of the Kelvin-Helmholtz modes with the vortex-shedding.

Universit é Paris-Saclay

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