robustness . . . . .

iii Nomenclature ρ density q velocity vector u, v, w velocity components in the aerodynamic reference frame ω vorticity vector ω = ∇ × q t time p static pressure

p i stagnation pressure p i = p 1 + γ-1 2 M 2 γ γ-1 T temperature T i stagnation temperature T i = T 1 + γ-
C p pressure coefficient: C p = p -p ∞ 1 2 ρ ∞ u 2 ∞ S ref ϕ
velocity potential: q = ∇ϕ for a potential flow S rate-of-stress tensor: S = 1 2 ∇q + ∇q T Ω vorticity tensor:

Ω = 1 2 ∇q -∇q T Q Q-criterion: Q = 1
CD drag coefficient CD = D 1 2 ρ ∞ u 2 ∞ S ref n
normal vector pointing outside the flow domain i unitary freestream direction vector V control volume S a surface of the body S e outer surface of the fluid volume S d downstream wake plane S w surface for the integration of wave drag S ′ w downstream plane of S w S wd downstream wake plane of the streamtube enclosing the shock S v surface for the integration of viscous drag S vd downstream wake plane of the streamtube enclosing the body and its boundary layer S cd complementary downstream wake plane S cd = S d \ (S wd ∪ S vd ) V w volume enclosed within S w V wd volume downstream of Although all these phenomena are interlinked, superimposed in the flow field and interacting with each other, a drag breakdown method attempts to artificially isolate each contribution. It can of course only be a modeling of the intricate flow physics, and leads to an utopian decomposition.

V w V v volume enclosed within S v V c complementary volume V c = V \ (V w ∪ V wd ∪ V v ) Subscript ∞ freestream

Problem Statement

The numerical methods for drag prediction and breakdown have been developed and extensively used for the past years. They have proved their reliability and accuracy on many steady cases, that is to say in cruise conditions, in high-lift conditions assuming a steady flow, around rotors in quasistatic conditions, etc. The performance in those conditions has therefore almost reached its highest level. To this day, no reliable method for the breakdown of drag into phenomenological components for unsteady flows only has been developed. Losses are however observed at the boundary of the flight envelop, at take-off or landing for example, when unsteadiness of the flow develops and cannot be neglected anymore. Those losses can imply important reduction of the overall performance, which can no more be left aside. Another point to keep in mind is the performance assessment of complex innovative systems, such as CRORs, which by nature cannot be evaluated in a quasi-static frame. The lack of information concerning the drag and its sources of such breakthrough designs brakes their development.

Many authors have looked into experimental and numerical drag breakdown methods for steady flows over the years. A detailed state-of-the-art review can be found in the first part of this dissertation. The contributions from the main authors are summarized here.

Van der Vooren [START_REF] Van Der Vooren | Drag/Thrust Analysis of a Jet-Propelled Transonic Transport Aircraft: Definition of Physical Drag Components[END_REF] has for example developed a phenomenological drag breakdown formulation based on thermodynamic considerations. He succeeded in formulating a method which allows to evaluate the drag due to shock waves, the drag due to viscous interaction within the boundary layers and wakes, and the drag due to wingtip vortices in the case of finite-span wings. Destarac [START_REF] Destarac | Far-Field Drag in Transonic Potential Flow: Analysis and Optimisation[END_REF] then enhanced the method by identifying and defining a spurious drag component, which is mainly due to the artificial diffusion introduced by numerical schemes. Their formulation has proved to be reliable and robust, and is now widely used in the designing phase by the industrial sector. It is however restricted to steady flows. Moreover, the induced drag component is defined by default as the remaining part of drag once the wave and viscous components defined. It does therefore not perfectly correspond to the phenomenological source which is the wingtip vortex. Finally, the formulation requires the use of an irreversible axial velocity which can become ill-defined in certain zones of the flow.

A first attempt in generalizing this formulation to unsteady flows has been proposed by Gariépy [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF]. He included additional terms coming from the unsteadiness in the equations for conservation of momentum and mass in a non-inertial reference frame. He also adopted an alternative expression for the irreversible axial velocity which had been developed by Méheut [START_REF] Méheut | Drag-Breakdown Methods from Wake Measurements[END_REF], so as to get rid of the ill-definition problem. His resulting unsteady drag breakdown method, which takes all unsteady terms into account, has been successfully tested on three unsteady test cases. He however chose to assign the unsteady additional contributions to a new drag component called unsteady drag, which does not correspond to a phenomenological source. The other drag components are also deprived from terms which would ensure a synchronization in time by taking the propagation in time and space into account.

Another family of formulations, which rely on replacing the pressure term by velocity and its derived variables such as vorticity, has also been extended to unsteady flows. The original steady version has been developed by Wu, Ma, and Zhou [START_REF] Wu | Vorticity and Vortex Dynamics[END_REF] with contributions from Noca, Shiels, and Jeon [START_REF] Noca | A Comparison of Methods for Evaluating Time-Dependent Fluid Dynamic Forces on Bodies, Using only Velocity Fields and their Derivatives[END_REF]. A breakdown of the drag for compressible flows has been suggested by Mele and Tognaccini [START_REF] Mele | Aerodynamic Force by Lamb Vector Integrals in Compressible Flow[END_REF]. Finally, Marongiu and Tognaccini [START_REF] Marongiu | Far-Field Analysis of the Aerodynamic Force by Lamb Vector Integrals[END_REF] proposed a breakdown in the unsteady case. These formulations may be promising for a drag breakdown in the unsteady case, since the physics involved is richer, and should allow defining the induced component directly rather than by default. However, the few breakdown attempts have proved rather unsatisfactory, when compared to other drag breakdown methods. The physics involved is also quite complex, and difficult to link to phenomenological observations. The maturity of such methods can therefore be questioned.

As a conclusion, we can summarize this short literature review as follows: on the one hand, several authors have contributed to the development of accurate and robust drag breakdown methods in the steady case. A first attempt in a generalization to unsteady flows has been suggested, but leads to the definition of drag components which are not phenomenological. On the other hand, formulations based on the velocity vector have been developed, with extension to compressible and unsteady flows. Despite their rich and promising physical content, they do not allow a convincing drag breakdown for the time being. Hence this study. It consists in developing a far-field drag prediction method aiming at a phenomenological breakdown of drag for unsteady flows.

Thesis Outline

The approach consists in developing an unsteady formulation starting from Van der Vooren's formulation. The first step consists in identifying the difficulties which we are likely to encounter. The second step is to propose a new proof of Van der Vooren's formulation which is directly generalizable to unsteady flows, the main contribution of this work concerning the volume splitting. The assumptions restricted to steady cases are also pinpointed at this step. The formulation is then rigorously extended to unsteady flows, including unsteady contributions into each drag component, wave, viscous and unsteady. Finally, the resulting method is discussed, and flaws identified. (Chapter I)

The next step consists in trying to improve the flaws detected at the end of the previous chapter. Concerning the robustness of the formulation, an alternative expression for the irreversible axial velocity found in the literature is studied and compared on steady cases. An unsteady version of the criterion used for the definition of the wave volume is also evaluated on two unsteady test cases. The physical background for the definition of the unsteady induced drag component is then looked into. Acoustic effects are highlighted and quantified, and a further breakdown proposed. Finally, the method is summarized, along with good practice recommendations. (Chapter II)

The new unsteady drag breakdown method is then applied to natural unsteady cases. A vortex shedding and a buffet case are chosen in order to assess increasingly complex cases. The wave component adds indeed to the breakdown for the buffet case. Orders of magnitude and synchronization in time are looked into in order to validate the results. Comparisons with another unsteady formulation are also performed, as well as comparisons with steady breakdown on the averaged flow field. (Chapter III)

The next step consists in assessing the formulation on mobile test cases: a pitching airfoil in both inviscid and viscous flows. The same method is used for the validation of the breakdown, with the addition of the investigation of the influence of the reduced frequency. Comparisons with steady computations at each angle of attack are finally carried out. (Chapter IV)

Finally, the unsteady drag breakdown is applied to complex cases. The aim is here to demonstrate the capabilities of the formulation. The first complex case is a 3-D pitching case. A wing tip vortex develops and evolves along with the pitching motion, allowing to better understand the breakdown of the induced drag component. A buffet case simulated with ZDES (Zonal Detached Eddy Simulations) is the last test case of this study. The rich content in terms of both spatial and time resolution gives promising results as far as the drag breakdown is concerned. URANS and ZDES results are also compared. (Chapter V)

Communications

The work contained in this dissertation has been presented in three international conferences. The formulation obtained at the end of Chapter I has also been published in the AIAA Journal, and the one obtained at the end of Chapter II is in preparation.

Conference Proceedings

[85] Toubin, H. and Bailly, D., "Far-Field Drag Decomposition for Unsteady Flows," 49th 3AF International Symposium of Applied Aerodynamics, Lille, France, 2014.

Presentation of the Main Existing Methods of Numerical Drag Prediction

The aim of this state-of-the-art review is to present the historical background which has led to the drag breakdown formulations used nowadays, and to introduce innovative axes recently explored by peers.

The first section concerns the review of the historical formulations for thermodynamic drag breakdown, starting from 1925. The formulations proposed by Betz, Jones, Oswatitsch, Maskell and Van der Vooren are presented, focusing on the assumptions used and the domain of application of each method. A set of illustrative examples of the capabilities of drag breakdown, using Van der Vooren's formulation, are also presented.

The formulations using asymptotic developments are deliberately left aside to focus in the second section on the formulations based on the velocity vector. These formulations have the benefit of having a very rich physical content. Breakdown into induced and profile components have been suggested and tested in the incompressible case. However, achieving a decomposition in the compressible case appears much more difficult.

There exist a few attempts of generalization to unsteady flows. Unsteady versions of the formulations based on the velocity vector are first looked into, with both incompressible and compressible versions. The breakdown in this case is however an even harder issue. The only unsteady formulation which allows a phenomenological breakdown, proposed by Gariépy as a generalization of Van der Vooren's theory, is finally presented, along with one of its application case: a pitching profile.

Historical presentation of the main thermodynamic methods

The far-field philosophy was inspired by Von Kármán, who suggested to apply the conservation laws on a control volume surrounding a body. The effort exerted on the body skin can thus be equivalently evaluated by integrations in the fluid domain. The drag can then be broken down into phenomenological components, such as lift-induced drag, which has been identified at the very beginning of aerodynamics for aircraft.

This section gives a review of the historical developments concerning thermodynamic formulations, starting from the steady theory which has led to Van der Vooren's formulation, now widely used, towards a first attempt in generalizing to unsteady flows, by Gariépy. Applications demonstrating the capabilities of far-field drag breakdown are also presented.

Steady methods

The thermodynamic breakdown of drag has first been developed for steady flows. The first formulation for the drag breakdown has been proposed in 1925. Other formulations have then been developed, but they only gave an expression for the profile drag. The first precise formulation for the induced drag dates back to 1972.

Betz

Far-field drag breakdown was first introduced by Betz [START_REF] Betz | A Method for the Direct Determination of Wing-Section Drag[END_REF] in 1925. His aim was to evaluate the profile drag of an aircraft using experimental measurements. Induced drag was indeed at the time approximated by formulas based on elliptic loading and aspect ratio. The profile drag was then evaluated as the difference between the total drag measured or computed and this approximated induced drag. The estimation of profile drag therefore lacked accuracy. The idea of Betz's formulation is to compute the profile drag from wake measurements. His formulation is valid for incompressible flows only.

In order to do so, he assumes that the velocity on a wake plane S d is aligned with the free stream velocity, and that the variations of stagnation temperature are negligible:

• v = w = 0 • T i = T i∞
The resulting expression is the following:

D Betz vw = S d (p t∞ -p t ) dS - ρ 2 S d (u ′ -u)(2u ∞ -u ′ -u) dS (1.1)
where p t denotes the total pressure p t = p + 1 2 ρu 2 , p t∞ the total pressure upstream of the body, and u ′ the velocity of a potential flow which would be identical to the real flow outside the vortical region. The second integral therefore limits to the vortical region. Betz admits that this second term is non trivial to compute, but shows, using orders of magnitude arguments, that it is negligible compared to the first one (less than 1 20 ).

Jones

Jones [START_REF] Jones | Measurement of Profile Drag by the Pitot-Traverse Method[END_REF] proposed another formulation in 1936. It is still valid for incompressible flows only, although extensions to compressible flows have been proposed later [START_REF] Bicknell | Determination of the Profile Drag of an Airplane Wing in Flight at High Reynolds Number[END_REF]. He uses more restrictive assumptions in order to get rid of the second integral in Betz's formulation which could not be easily evaluated by wake measurements. He assumes that on the wake plane S d the static pressure is homogeneous and equal to the reference static pressure, and that the velocity is aligned with the free stream velocity:

• p = p ∞ • v = w = 0
Jones obtains a formulation which depends only on the dynamic and static pressure for the estimation of profile drag:

D Jones vw = p d∞ S d 2 p * d -p * (1 -p * d ) dS (1.2)
where p * d and p * are the dimensionless dynamic and static pressure: p * d = p d p d∞ and p * = p p d∞ , the dynamic pressure being defined as p d = 1 2 ρu 2 = p tp. The surface S d can now be limited to the wake boundary of the streamtube enclosing the body, the boundary layer, and the wake. In practice, Jones advises to move the wake plane close to the trailing edge of the body. He indeed observes that the formulation does not depend on the location of the wake plane, except in what he calls dead zones. These dead zones are regions where p * d becomes negative, hence the formulation becomes undefined. These zones may appear very close behind a bluff body.

Jones applied his formulation to a flight test case. He used the airplane Hart K.1442 and performed the drag evaluation using static-pressure tubes. He compared the drag prediction using two variants of his formulation, and Betz's formulation, as shown in Figure 2. We can see that Betz's and Jones's formulations are equivalent, except when the integration is performed too close to the trailing edge. Jones also investigated the effect of smooth or fabric wings, proving that the profile drag was indeed reduced using smooth surfaces.

Oswatitsch

Oswatitsch [START_REF] Oswatitsch | Gas Dynamics[END_REF] developed in 1956 a formulation for the profile drag which differs from the previous ones by its nature. He indeed chose to use thermodynamic considerations instead of assuming no transverse velocity on a wake plane. Expressing the velocity and pressure as a function of the thermodynamic variables H and s, and neglecting second order terms and assuming isenthalpic flow, he achieves the expression for the profile drag:

D Oswatitsch vw = T ∞ u ∞ S d ρu∆s dS (1.3)
This definition allows avoiding using restrictive assumptions on a wake plane. Méheut has shown during his Ph.D. thesis [START_REF] Méheut | Évaluation des Composantes Phénoménologiques de la Traînée d'un Avion à partir de Résultats Expérimentaux[END_REF] that all three formulations are equivalent at first order.

Maskell

In 1972, Maskell [START_REF] Maskell | Progress towards a Method for the Measurement of the Components of Drag for a Wing of Finite Span[END_REF] achieved a further step into drag breakdown by proposing an expression for the induced drag as well. He starts from an expression for the total drag, obtained using the definition of the total pressure, and keeps the cross components of the velocity vector. His expression for the profile drag remains very similar to Betz's formulation, except for the introduction of a blocking velocity u b which allows taking the blocking of the wind tunnel into account in the experimental case.

D Maskell vw = S d (p t∞ -p t ) dS + ρ 2 S d (u ′ -u) u ′ + u -2(u ∞ + u b ) dS (1.4)
The expression of the induced drag is obtained from the difference between an expression of the total drag and this profile drag. An alternative expression as a function of the stream function ψ, the velocity potential φ, the axial vorticity ζ = ∂w ∂y -∂v ∂z , and the source term σ = ∂v ∂y + ∂w ∂z :

D Maskell i = ρ S d (v 2 + w 2 ) dS = ρ S d (ψζ -φσ) dS (1.5)
This second expression allows reducing the integration surface S d to the downstream boundary of the streamtube enclosing the body, the boundary layer, and the wake.

Another important contribution was brought by Kusunose [START_REF] Kusunose | A Wake Integration Method for Airplane Drag Prediction[END_REF][START_REF] Kusunose | Extension of Wake Survey Analysis Method to cover Compressible Flows[END_REF] who extended their application domain to compressible flows, suggesting to breakdown the profile drag into wave and viscous components using distinct integration surfaces. All these formulations were designed in order to address experimental results probed in a wake plane. The objective was therefore to obtain the best accuracy level with the smallest measuring plane. The development of CFD methods has led to a completely different approach, with data available anywhere in the flow. Methods better suited for numerical results have therefore been developed, Van der Vooren's formulation being the most advanced and practical one.

Van der Vooren and Destarac

Van der Vooren first addressed the thermodynamic breakdown of far-field drag in 1990 with Slooff [91], in the framework of a state-of-the-art review. His method is like a crossover between the assumptions on a wake plane and the thermodynamic considerations of Oswatitsch. The main contribution is the breakdown of the profile drag into wave and viscous components. It was however not directly applicable to industrial cases due to discrepancies between the near-field and far-field values as observed earlier by Yu, Chen, Samant and Rubbert [START_REF] Yu | Inviscid Drag Calculations for Transonic Flows[END_REF], and Steger and Baldwin [START_REF] Steger | Shock Waves and Drag in the Numerical Calculation of Compressible, Irrotational Transonic Flows[END_REF]. This discrepancy issue was solved by the notion of spurious drag first introduced by Lock [START_REF] Lock | Prediction of the Drag of Wings at Subsonic Speeds by Viscous / Inviscid Interaction Techniques[END_REF] and formalized by Destarac in 1993 [START_REF] Destarac | Far-Field Drag in Transonic Potential Flow: Analysis and Optimisation[END_REF].

Destarac and Van der Vooren then associated to implement a drag breakdown method in a code at ONERA. The method was also extended to powered configurations and published in [START_REF] Van Der Vooren | Drag/Thrust Analysis of a Jet-Propelled Transonic Transport Aircraft: Definition of Physical Drag Components[END_REF]. The formulation for non propelled aircraft presented hereafter is drawn from this article.

Theoretical formulation for non-propelled aircraft

The formulation starts from the far-field equation, obtained applying the conservation of mass and momentum in a control volume V whose boundary is S A ∪ S F ∪ S D (see Figure 3):

S A (p -p ∞ )(i • n) -(τ x • n) dS = - S F ∪S D ρ(u -u ∞ )(q • n) + (p -p ∞ )(i • n) -(τ x • n) dS (1.6)
The near-field drag can be recognized on the left-hand side, with its pressure and friction components. The right-hand side defines the far-field drag. The upstream and lateral surfaces S F are chosen at infinity, so that the integration limits to S D . The formulation aims at breaking down this total drag definition. The principle consists in assuming a flow free of vortices, i.e. a case where only profile drag is exerted on the body. In this case, the flow on a wake plane S D far from the sources is such that:

• τ x = 0 • v = w = 0 • p = p ∞
Under these assumption, the axial velocity defect uu ∞ can be expressed as:

∆u = u ∞ 1 + 2∆H u 2 ∞ - 2 (γ -1) M 2 ∞ e γ-1 γ ∆s r -1 -u ∞ (1.7)
The notation ∆u is linked to the notation u irr used in this manuscript by the simple relation:

∆u = u irr -u ∞ .
The total far-field drag, which is equal to the profile drag D vw in the case of a flow free of vortices, can now be derived using the former assumptions:

D = D vw = - S D ρ∆u(q • n) dS (1.8)
This wake surface integral is then turned into a volume integral over the volumes V w and V v surrounding the sources of drag, as defined in Figure 3, using the divergence theorem and assuming truly inviscid flow outside V w ∪ V v . No other justification is provided by the authors. The profile drag can then naturally be broken down into a wave and a viscous contribution from a splitting of the integration volumes:

D w = - Vw ∇ • (ρ∆uq) dV
(1.9)

D v = - Vv ∇ • (ρ∆uq) dV (1.10)
The complementary part defines the induced drag, which is also written as a volume term, so that the integration on the body skin coming from the divergence theorem can be written as a function of the near-field components:

D i = - Vw∪Vv ∇ • ρ(u -u ∞ -∆u)q + (p -p ∞ )i -τ x dV -D p -D f (1.11)
This expression was obtained using the property of vector

f = -ρ(u -u ∞ )q -(p -p ∞ )i + τ x : ∇ • f = 0 (1.12)
The balance between near-field and far-field total drag is theoretically ensured.

Numerical deviations from the theory and practical refinements

Besides this formulation which allows a breakdown into the three phenomenological drag components with a correct balance of the total drag, the main contribution of this method is the notion of spurious drag.

Two causes of numerical deviations are identified in the article:

• the spurious contribution of ∇ • (ρ∆uq) in the near-field outside V w ∪ V v , which is due to the artificial creation of entropy.

• the numerical and physical dissipation of the trailing vortices further outside of V w ∪ V v which causes conversion of apparent induced drag into apparent viscous drag.

These spurious phenomena can be accounted for by defining a spurious drag component D sp and by enlarging the integration volume for the induced drag. Denoting V sp the volume outside V w ∪ V v :

D sp = - Vsp ∇ • (ρ∆uq) dV
(1.13)

D i = - V ∇ • ρ(u -u ∞ -∆u)q + (p -p ∞ )i -τ x dV -D p -D f (1.14)
The balance between near-field and far-field can then be written as:

D p + D f = D w + D v + D i + D sp (1.15)
This formulation has known many extensions and refinements over the years. Among the extensions:

• propelled aircraft [START_REF] Van Der Vooren | Drag/Thrust Analysis of a Jet-Propelled Transonic Transport Aircraft: Definition of Physical Drag Components[END_REF] • rotating frames with ALE simulations (hovering helicopter rotors [START_REF] Verley | Torque Analysis based on the CFD "Far-Field" Philosophy for Rotors in Hover[END_REF] and counter rotating open rotors [START_REF] Méheut | Thrust and Torque Far-field Analysis of Propeller and Counter Rotating Open Rotor Configurations[END_REF])

The practical refinements, such as the attribution of the friction contribution to the profile drag instead of the induced drag, the use of a one vector formulation, or the use of integration surfaces instead of volumes [START_REF] Destarac | Far-field / Near-field Drag Balance and Applications of Drag Extraction in CFD[END_REF][START_REF] Destarac | Drag Extraction from Numerical Solutions to the Equations of Fluid Dynamics: The Far-Field Philosophy[END_REF][START_REF] Hue | Computational Drag Prediction of the DPW4 Configuration Using the Far-Field Approach[END_REF], are detailed in Chapter I. Other Ph.D. dissertations dedicated to far-field drag extraction inspired from Van der Vooren's theory are also listed in the bibliography [START_REF] Amant | Calcul et Décomposition de la Traînée Aérodynamique des Avions de Transport à partir de Calculs Numériques et d'Essais en Soufflerie[END_REF][START_REF] Arntz | Civil Aircraft Aero-thermo-propulsive Performance Assessment by an Exergy Analysis of High-fidelity CFD-RANS Flow Solutions[END_REF][START_REF] Bergé | Études et Développement de Méthodes d'Évaluation de la Traînée Aérodynamique Appliquées à l'Installation Motrice d'Avions Civils de Transport[END_REF][START_REF] Esquieu | Évaluation de la Traînée d'un Avion de Transport à partir de Calculs Numériques de Mécanique des Fluides[END_REF][START_REF] Méheut | Évaluation des Composantes Phénoménologiques de la Traînée d'un Avion à partir de Résultats Expérimentaux[END_REF][START_REF] Verley | Évaluation du Couple "Champ Lointain" d'un Rotor d'Hélicoptère en[END_REF].

Examples of application on steady non-propelled cases

There have been many applications of Van der Vooren's formulation, implemented and refined in the ONERA code family ffd. Far-field drag prediction is used to assess the performance of new designs, and to optimize one or several drag components through shape deformations or flow control. Here is a representative set of applications which allow to illustrate the capabilities of the formulation.

Study of the wave drag reduction using bumps on a profile

The first example is a study of the wave drag reduction using local curvature modifications slightly downstream of the shock wave location on a supercritical profile, in the framework of a cooperation between ONERA and Airbus Germany [START_REF] Destarac | Far-field / Near-field Drag Balance and Applications of Drag Extraction in CFD[END_REF]. The aim was to assess the improvement in terms of performance suggested by Ashill and Fulker [START_REF] Ashill | A Novel Technique for Controlling Shock Strength of Laminar Flow Airfoil Sections[END_REF], and to evaluate which drag component was affected by such a shape deformation.

Two bumps have been designed (see Figure 4(a)). The performance of the three profile shapes have then been assessed both experimentally in the ONERA T2 wind tunnel and numerically at Mach number 0.762. The results of the drag extraction are presented in Figure 4(b). We can see that the wave drag component is effectively decreased of 20% for the bump 1 and 50% for the bump 2, whereas the viscous drag remains approximately constant. There is also a good agreement with the experimental results, with an error of approximately 10 drag counts, and very similar relative gains in terms of total drag. Figure 4 -Study of the wave drag reduction thanks to bumps using Van der Vooren's formulation (extracted from [START_REF] Destarac | Far-field / Near-field Drag Balance and Applications of Drag Extraction in CFD[END_REF])

Study of the induced drag reduction using unconventional wingtips designs

The second example is the reduction of the induced drag by modifying the wingtip of a wing, extracted from a paper by Grenon [START_REF] Grenon | Numerical Study of Unconventional Wing Tip Devices for Lift-Induced Drag Reduction[END_REF]. Induced drag can indeed be reduced by the use of devices at the tip of wings, such as winglets. The levels of drag reduction are however generally small and require accurate drag prediction. The effect of the wingtip devices on the other drag components must also be carefully looked into. Phenomenological breakdown provides insight into the balance between the reduction of one component and the possible penalty for another.

Three wingtip devices have been tested: two winglets of different heights and a spiroid, the overall span being kept constant (see Figure 5(a)). The four configurations have been numerically simulated in an inviscid flow using the Euler equations. The results of the application of Van der Vooren's formulation are presented in Figures 5(b), and 5(c). We can see that all three wingtip devices allow to slightly reduce the induced drag (Figure 5(b)), the long winglet 2 being the most efficient. The wave drag component (Figure 5(c)) remains approximately constant with both winglets designs. It is however slightly increased with the spiroid, as a small shock wave appears on the spiroid itself. The induced drag reduction achieved by the use of the spiroid is unfortunately canceled by the wave drag penalty.

The successful phenomenological breakdown teaches us that the performance of the spiroid could be improved by redesigning its shape to avoid the formation of the small shock wave.

Optimization of competing drag components

The third example is drawn from an article by Yamazaki, Matsushima and Nakahashi [START_REF] Yamazaki | Aerodynamic Design Optimization Using the Drag-Decomposition Method[END_REF]. It consists in the aerodynamic shape optimization of the planform of the DLR F6 wing configuration, using a Taylor expansion of Van der Vooren's formulation developed by Tognaccini [START_REF] Tognaccini | Methods for Drag Decomposition, Thrust-Drag Bookkeeping from CFD Calculations[END_REF].

Yamazaki et al. first show an interesting result on a preliminary study: they show that taking the far-field expression of drag instead of the near-field, which includes the spurious drag, allows to better reduce the drag on the profile. They then tackle the issue of planform design optimization with a multi-objective approach: they search a trade-off between induced and wave drag components using a genetic algorithm. They obtain a Pareto set, as depicted in Figure 6.

The designs corresponding to several points of the Pareto front are also presented in Figure 7. This optimization method allows to isolate the effect of each drag coefficient, and to obtain very rich design information. The designer can then choose which physical drag component, or which corresponding design, he wants to promote, based on the targeted mission of feasibility considerations.

Performance evaluation of a wing-body configuration

The last example is the performance assessment of the wing-body configuration from the 5th edition of the AIAA Drag Prediction Workshop [START_REF] Hue | Fifth Drag Prediction Workshop: Computational Fluid Dynamics Studies Carried Out at ONERA[END_REF]. Figure 8 shows the integration surfaces for the wave and viscous drag components computed by the ONERA code ffd72. The formulation is integrated on these surfaces determined using physical criteria (see Chapter I, Section 2.4.3.3).

The results of the far-field drag extraction are presented in Figures 9(a), and 9(b). We can see in Figure 9(a) that the three phenomenological drag components increase with the angle of attack, while the friction part of the viscous component remains approximately constant. The evolution of the share of each component, represented in Figure 9(b) is very interesting: the share of the viscous drag decreases while the share of the wave drag increases, and the share of the induced drag remains almost constant.

There are many other examples of application [START_REF] Van Der Vooren | Drag/Thrust Analysis of a Jet-Propelled Transonic Transport Aircraft: Definition of Physical Drag Components[END_REF][START_REF] Méheut | Thrust and Torque Far-field Analysis of Propeller and Counter Rotating Open Rotor Configurations[END_REF][START_REF] Destarac | Far-field / Near-field Drag Balance and Applications of Drag Extraction in CFD[END_REF], in particular contributions to the Drag Prediction Workshops organized by AIAA, where ONERA with the ffd72 code compares often very well with the other participants [START_REF] Hue | Computational Drag Prediction of the DPW4 Configuration Using the Far-Field Approach[END_REF][START_REF] Hue | Fifth Drag Prediction Workshop: Computational Fluid Dynamics Studies Carried Out at ONERA[END_REF]. An interesting application is also a mesh refinement method based on the spurious drag coefficient developed by Yamazaki [START_REF] Yamazaki | Drag Decomposition-Based Adaptative Mesh Refinement[END_REF]. The formulation has therefore proven its accuracy and robustness on steady cases, but is however not applicable as such on unsteady cases. Another flaw of the formulation is the definition by default of the induced drag component.

There exists a second family of formulations, which rely on asymptotic expansions, along the same line of thoughts as Oswatitsch, and which are valid only for wake integrations [START_REF] Giles | Wake Integration for Three-Dimensional Flowfield Computations: Theoretical Development[END_REF][START_REF] Hunt | Wake Integration for Three-Dimensional Flowfield Computations: Applications[END_REF]. The idea was also used by Chao and Van Dam [START_REF] Chao | Wing Drag Prediction and Decomposition[END_REF] and Tognaccini [START_REF] Tognaccini | Methods for Drag Decomposition, Thrust-Drag Bookkeeping from CFD Calculations[END_REF]. All these formulations however lack accuracy and are not easily generalizable to unsteady flows due to the many restrictive assumptions made.

First unsteady generalization of Van der Vooren's formulation

Gariépy tackled during his Ph.D. thesis [START_REF] Gariépy | Méthodes d'Extraction et de Décomposition de la Traînée pour des Écoulements Permanents et Transitoires[END_REF] the issue of a generalization to unsteady flows of Van der Vooren's formulation. He proposed in 2013 [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF] a first generalization of the far-field drag breakdown to unsteady flows. His method consists in computing the additional terms due to the unsteadiness of the flow and to allocate them to a new drag component: the unsteady drag D uns . The other drag components are slightly modified compared to Van der Vooren's definition. An alternative expression for the axial velocity defect, first introduced by Méheut in a study of experimental drag breakdown [START_REF] Méheut | Drag-Breakdown Methods from Wake Measurements[END_REF], is in particular used.

Theoretical developments

Gariépy starts with writing the far-field equation in the relative reference frame (see Figure 10):

∂ ∂t V ρ(u -u ∞ ) dV + ∂V ρ(u -u ∞ )(q r • n) + (p -p ∞ )(i • n) -(τ x • n) dS + V a x dm = 0 (1.16)
where q r denotes the relative velocity, and a x the relative acceleration of the control volume.

He then applies the same reasoning as Van der Vooren and Destarac for the profile drag, assuming a flow free of vortices to define the irreversible part of drag. A divergence with the steady theory is the use of an alternative expression for the irreversible axial velocity defect ∆u. As Jones [START_REF] Jones | Measurement of Profile Drag by the Pitot-Traverse Method[END_REF] had first foreseen, this expression can indeed be undefined in certain regions of the flow for given values of local stagnation pressure (see Chapter I, Section 4.1.1 for the detailed proof). Since this situation can occur more frequently for unsteady flows, Gariépy chooses to use an alternative expression of the axial velocity defect using the reversible axial velocity u * suggested by Méheut [START_REF] Méheut | Drag-Breakdown Methods from Wake Measurements[END_REF], and studied in [START_REF] Gariépy | A New Axial Velocity Defect Formulation for a Far-Field Drag Decomposition Method[END_REF]:

∆u * = u -u * = u -u ∞ 1 + 2∆H u 2 ∞ - 2 (γ -1) M 2 ∞   p p ∞ γ-1 γ -1   - v 2 + w 2 u 2 ∞ (1.17)
He then uses the divergence theorem in order to turn the surface integral into a volume integral. The next step consists in splitting the control volume into regions enclosing the sources of drag, the exact same way Van der Vooren proceeded. He also chooses to allocate the contribution outside these sources to a spurious drag component, arguing that ∆u * should be theoretically zero outside V w and V v :

D w = - Vw ∇ • (ρ∆u * q r ) dV (1.18) D v = - Vv ∇ • (ρ∆u * q r ) dV (1.19) D sp = - V \(Vw∪Vv) ∇ • (ρ∆u * q r ) dV (1.20)
Gariépy then argues that the stagnation enthalpy variation is solely due to reversible processes and linked to a time derivative, so that its contribution should be assigned to the unsteady drag component, along with the time derivative and the acceleration of the control volume terms. To do so, he defines a reversible axial velocity u * ,s assuming a steady flow:

u * ,s = u ∞ 1 - 2 (γ -1) M 2 ∞   p p ∞ γ-1 γ -1   - v 2 + w 2 u 2 ∞ (1.21)
and computes the reversible axial velocity caused by unsteady fluctuations as u * ,u = u *u * ,s . He can therefore define the unsteady drag component as:

D uns = - V ∇ • ρu * ,u q r dV - ∂ ∂t V ρ(u -u ∞ ) dV - V a x dm (1.22)
The induced drag component is defined similarly to the steady theory by default as the complementary part in the total drag. Using again the divergence theorem, Gariépy defines the induced component as a volume integral:

D i = - V ∇ • ρu * ,s q r + (p -p ∞ )i -τ x dV -D p -D f (1.23)
The formulation can be summarized as:

D p + D f = D w + D v + D i + D uns + D sp (1.24)

Results on a pitching case

Gariépy presents results on three unsteady cases: a motionless profile with a variation in time of the reference Mach number, a profile pitching between -0.5 and +0.5 deg at the reduced frequency of 0.04, and a profile pitching between -5 and +5 deg at the reduced frequency of 0.1. The latter is summarized hereafter.

The test case is a NACA0012 airfoil pitching in a flow at Mach 0.3 and Reynolds number 6.6 million. Concerning the definition of the integration volume, the authors use the same physical criteria as at ONERA, but the extension of the volumes is said to vary during the convergence process so as to minimize the gap between near-field and far-field estimations. The resulting volumes are not shown in the article.

The drag components are extracted over one period. Figure 11(a) shows the evolution of the total drag coefficients, computed using the far-field and near-field expressions. The drag coefficients reach a maximum when the angle of attack is at +5 and -5 deg. There is a very good agreement between the near-field and the far-field estimations.

The results of the unsteady drag breakdown are presented in Figure 11(b). Only the viscous and unsteady components are visible. There is no mention of the induced component, although both curves do not add up to give the far-field drag curve. The wave drag component is zero since the flow is subsonic and therefore the volume for the integration of wave drag is zero. We can see that the viscous drag evolves very little as a function of the angle of attack. This is due to the fact that α remains small so that there is no flow separation. The variations in drag come therefore mainly from the variations of the unsteady component.

Gariépy's formulation will be discussed in Chapter I, Section 4.3. The results of his formulation, which has been implemented in the ONERA code ffd72 will also be compared in Chapters III, and IV.

A radically different family of far-field methods, the family of formulations based on the velocity vector, or vortex methods as they are sometimes called, seem promising in terms of physical content. The induced drag components in particular, which was defined by default in the previous methods, may have a true definition, at least in the steady case. They are summarized in the next section. In 1989, Wu and Wu [START_REF] Wu | Vortical Sources of Aerodynamic Force and Moment[END_REF] proposed to replace the integral of pressure in the far-field equation by an integral of velocity related variables. This manipulation is especially interesting for experimental studies where the velocity field is easily acquirable contrarily to the pressure. As far as drag breakdown is concerned, it also allows to introduce the vorticity vector, as well as the Lamb vector, which seem promising for defining the induced drag component. The way to a whole new family of drag formulations was open.

Formulation for steady incompressible flows

The formulation as first expressed by Wu and Wu is described in [START_REF] Wu | Vortical Sources of Aerodynamic Force and Moment[END_REF]. They use integration by parts to replace the surface pressure integration by volume integral using velocity and vorticity components. In this first paper, the authors also use enthalpy and entropy to suggest a breakdown into friction, wake vortices, and pressure variations components.

In 1997, Noca derived in his Ph.D. thesis [START_REF] Noca | On the Evaluation of Time-Dependent Fluid-Dynamic Forces on Bluff Bodies[END_REF] another expression of the force for unsteady, incompressible flows, which was later rederived and included in Wu's book [START_REF] Wu | Vorticity and Vortex Dynamics[END_REF]. The unsteady extension of Noca will be described in Section 2.5.1. The proof is similar to Wu's first method: it uses mathematical identities for the pressure and the impulse-momentum which allow replacing the surface pressure integral by a volume integral of velocity and vorticity combinations, with the appearance of the spatial dimension as a parameter.

The expression for the force F (drag and lift components), for steady incompressible cases, is:

F = - ρ N -1 V x × (∇ × l) dV - µ N -1 Se x × n × (∇ × ω) dS + µ Se ω × n dS (2.1)
where l = ω × q denotes the Lamb vector, µ the dynamic viscosity, x the position vector, and N the spatial dimension.

The main issue consists now in breaking down this expression into phenomenological drag components.

Breakdown into induced and profile components for steady incompressible cases

No breakdown is directly achievable with Equation 2.1. Wu, Ma, and Zhou therefore proposed an alternative expression:

F = -ρ V l dV - ρ N -1 ∂V x × (n × l) dS + F Se (2.2)
where F Se denotes the sum of the last two surface integrals:

F Se = - µ N -1 Se x × n × (∇ × ω) dS + µ Se ω × n dS (2.3)
At high Reynolds numbers, F Se has been shown to be negligible, and the authors argue that the volume integral yields the lift and induced drag, while the surface integral yields the profile drag:

D i = ρ V (vω z -wω x ) dV (2.4) D vw = - ρ N -1 ∂V (yl y + zl z ) dS (2.5)
Yang, Zhang, An, and Wu [START_REF] Yang | Steady Vortex Force Theory and Slender-Wing Flow Diagnosis[END_REF] applied in 2007 the steady incompressible formulation on a slender delta-wing. They first notice that the Lamb-vector field is highly localized in the flow field, within the near wall zone of the boundary layer, and has a very complicated behavior, with high peaks, as depicted in Figure 12. The breakdown into induced and profile component is very difficult to validate, owing to the fact that no alternative value is provided in the paper. The results are presented in Figures 13(a Marongiu, Tognaccini, and Ueno also addressed this issue in 2010 [START_REF] Ueno | Lift Induced Drag Computation by Lamb Vector Integration[END_REF]. They first proved that Equation 2.4 is equivalent to Prandtl's lifting line theory [START_REF] Prandtl | Application of Modern Hydrodynamics to Aeronautics[END_REF] for a steady inviscid incompressible flow around a large aspect ratio wing at small angles of attack. Applications on 2-D and 3-D test cases were also proposed, but the comparisons with Maskell's formulation [START_REF] Maskell | Progress towards a Method for the Measurement of the Components of Drag for a Wing of Finite Span[END_REF] for the induced drag and Tognaccini's formulation [START_REF] Tognaccini | Methods for Drag Decomposition, Thrust-Drag Bookkeeping from CFD Calculations[END_REF] for the profile drag were not very accurate. The profile drag estimation by Equation 2.5 was described as less robust.

In 2013, they applied the breakdown on more test cases, with more accurate comparisons [START_REF] Marongiu | Lift and Lift-Induced Drag Computation by Lamb Vector Integration[END_REF]. The authors show in particular that the integration domain for the Lamb vector can be reduced to the boundary layer region without implying a loss of accuracy. Figure 14 shows that the induced drag computed by Equation 2.4 is very close to the induced drag using Maskell's formulation (Equation 1.5). The volume expression however seems less robust when the downstream extension of the integration volume is reduced towards the trailing edge.

The breakdown into induced and profile drag in the steady incompressible case can therefore be considered satisfactory. Civil aircraft applications however require to study compressible cases as well. The breakdown may not be as easily achieved in this case.

Extension to steady compressible flows

Compressible extension of the formulation has been addressed by Wu, Ma, and Zhou [START_REF] Wu | Vorticity and Vortex Dynamics[END_REF]. They simply add the density gradient terms which had been neglected. The expression for the total force in the compressible case is the following:

F = - V ρl - 1 2 q 2 ∇ρ dV - 1 N -1 ∂V x × n × ρl - 1 2 q 2 ∇ρ dS + F Se (2.6)
In 2007, Zhu, Beanman, and Graham [START_REF] Zhu | Prediction of Drag and Lift of Wings from Velocity Fields[END_REF] used a similar steady compressible formulation to compute the total drag of a 3-D configuration, without any breakdown. The total drag is accurately predicted by the formulation.

Breakdown in the steady compressible case

A first attempt of breakdown has been proposed by Mele and Tognaccini in 2014 [START_REF] Mele | Aerodynamic Force by Lamb Vector Integrals in Compressible Flow[END_REF]. They manipulate Equation 2.6 in order to gather the density gradients terms into a same force component:

F l = - V ρl dV (2.7)
F ρ = 1 N -1 V x × ∇ × 1 2 q 2 ∇ρ dV (2.8)
F S = - 1 N -1 ∂V x × (n × (ρl)) dS (2.9)
The component on the outer surface S e of the surface term F S is shown to be negligible for the test cases considered. An alternative breakdown is also proposed, where density gradient terms are kept in the surface component F S , but the results seem less promising.

This compressible formulation is applied on a 2-D viscous test case, a NACA0012 airfoil, for several Mach numbers. The results of the drag breakdown are summarized in Figure 15. We can see that the coefficients CD l and CD ρ compensate each other until M ∞ = 0.7, above which the values are said to lack accuracy. The drag due to compressibility effects is actually a thrust. The total drag is accurately predicted by CD S . Therefore, this breakdown method does not allow to extract a wave component yet, since compressibility effects produce thrust according to this decomposition.

For completeness, we can also cite Chang, Su and Lei [START_REF] Chang | On Aerodynamic Forces for Viscous Compressible Flows[END_REF] who proposed an extension to compressible flows of the projection theory introduced by Quartapelle and Napolitano [START_REF] Quartapelle | Force and Moment in Incompressible Flows[END_REF]. Their formulation however provides only the pressure contribution, and requires the knowledge of proper harmonic functions.

Generalization to unsteady flows

Extension to unsteady flows of the velocity-based formulations presented in the previous section has been proposed as early as 1997 by Noca [START_REF] Noca | On the Evaluation of Time-Dependent Fluid-Dynamic Forces on Bluff Bodies[END_REF]. Unsteady compressible flows have also been addressed very recently, with a few hints provided regarding a possible breakdown of drag into phenomenological components.

Noca

Noca studied the extension of the vortex force theory to unsteady flows during his Ph.D. His thesis [START_REF] Noca | On the Evaluation of Time-Dependent Fluid-Dynamic Forces on Bluff Bodies[END_REF] is remarkably well detailed and pedagogical. His formulation has been summarized in Wu's book [START_REF] Wu | Vorticity and Vortex Dynamics[END_REF], and can be written in the incompressible case as:

F = - ρ N -1 V x × ∂ω ∂t dV -ρ V l dV - ρ N -1 Se x × (n × l) dS + F Se + F Sa (2.10)
where F Se is still defined by Equation 2.3 and F Sa is the force exerted on the body surface S a :

F Sa = ρ N -1 Sa x × σ a dS (2.11)
where σ a is the boundary vorticity flux due to the acceleration of the body surface. Noca applied this formulation to compute the total drag of the unsteady flow around a circular cylinder performing a transverse oscillation. The results of the total drag prediction computed from experimental data are provided in [START_REF] Noca | A Comparison of Methods for Evaluating Time-Dependent Fluid Dynamic Forces on Bodies, Using only Velocity Fields and their Derivatives[END_REF], as depicted in Figure 16. No breakdown is however suggested.

Wu

As mentioned earlier, Wu also addressed the issue of extension to unsteady flows in his book. In 2007, Wu, Lu, and Zhuang studied the behavior of the different terms in Equation 2.10 in the case of the Von Kármán instability downstream of a circular cylinder in [START_REF] Wu | Integral Force Acting on a Body due to Local Flow Structures[END_REF]. They plotted the field distribution of each of the source terms in the formulation, and showed that some of them allowed capturing different flow structures. The interpretation is however not easy, due to the complexity of the flow considered and of the patterns observed for the source terms.

Marongiu

Marongiu also addressed this issue during his Ph.D. thesis [START_REF] Marongiu | On the Aerodynamic Force of the Oscillating Airfoils[END_REF]. He focused on the impact of turbulence on the formulation. He also wrote an equivalent of Equation 2.10 in 2-D:

CD l = 2 V v ω dV
(2.12)

CD uns = -2 V y ∂ω ∂t dV (2.13)
CD S = -2 Se y ω q • n dS (2.14) CD Sa = 2 Sa y n × (a -l) dS (2.15)
where a is the fluid acceleration. Marongiu and Tognaccini applied this unsteady formulation to a pitching profile in [START_REF] Marongiu | Far-Field Analysis of the Aerodynamic Force by Lamb Vector Integrals[END_REF]. The drag component CD Sa computed at the skin is very low. The comparison between CD S and CD l allows comparing the impact of the bound vorticity and the free vorticity respectively. The bound vorticity (vorticity created in the boundary layer) leads to a lift contribution which is in phase with the angle of attack, and a drag contribution shifted towards higher values (see Figure 17(a)), whereas the free vorticity (vorticity shed in the wake) gives out of phase results (see Figure 17(b)). The unsteady contribution CD uns is however not analyzed.

Xu

An extension to compressible flows of Equation 2.10 is easily achievable, similarly to the steady case. It was proposed by Wu, Ma, and Zhou [START_REF] Wu | Vorticity and Vortex Dynamics[END_REF]: When analyzing the drag experienced by a wavy cylinder, Xu, Chen, and Lu [START_REF] Xu | Large-Eddy Simulation of the Compressible Flow past a Wavy Cylinder[END_REF] used Equation 2.16 to breakdown the drag. They propose the following breakdown:

F = - 1 N -1 V x × ∇ × ∂ρq ∂t dV - V ρl - 1 2 q 2 ∇ρ dV - 1 N -1 Se x × n × ρl - 1 2 q 2 ∇ρ dS + F Se + F Sa (2.16)
F l = - V ρl dV - 1 N -1 ∂V x × (n × (ρl)) dS (2.17)
F ρ = V 1 2 q 2 ∇ρ dV + 1 N -1 ∂V x × n × 1 2 q 2 ∇ρ dS (2.18)
F uns = - 1 N -1 V x × ∇ × ∂ρq ∂t dV (2.19)
The results on their test case is the following: F Se due to viscous effects gives a negligible contribution. F Sa contributes to about half the total drag. The compressibility effects represented by F ρ contribute to 3,500 d.c., and the vortex force F l gives a drag coefficient of 2,400 d.c. The compressible effects are therefore more important in this test case, which is confirmed by another article. No estimation for the unsteady drag is provided. This breakdown is however not equivalent to a viscous/wave/induced breakdown, since the vortex force thus defined contributes to several of these drag components, not only to the induced one. The compressibility effects are also not limited to the wave contribution.

Other contributions to the unsteady generalization

There have also been many contributions from the experimental field, with authors using this family of unsteady formulations with time-resolved PIV data [START_REF] Devoria | On Calculating Forces from the Flow Field with Application to Experimental Volume Data[END_REF][START_REF] Sällström | Force Estimation from Incompressible Flow Field Data Using a Momentum Balance Approach[END_REF][START_REF] Sterenborg | Assessment of PIV-based Unsteady Load Determination of an Airfoil with Actual Flap[END_REF][START_REF] Villegas | Evaluation of Unsteady Pressure Fields and Forces in Rotating Airfoils from Time-Resolved PIV[END_REF]. They have however shown no interest in the breakdown of drag. In 2014, Liu et al. [START_REF] Liu | Longitudinal-Transverse Aerodynamic Force in Viscous Compressible Complex Flow[END_REF] proposed a breakdown of the drag starting from the velocity vector based formulation, in terms of longitudinal and transverse components, which are respectively related to compressibility and irreversible thermodynamics, and shearing. This breakdown is however difficult to link with the usual phenomenological components.

Niu and Chang also proposed an extension of their theory to unsteady flows in [START_REF] Niu | How Do Aerodynamic Forces of the Pitching Rigid and Flexible Airfoils Evolve?[END_REF], which again requires the computation of proper harmonic functions. The breakdown suggested in terms of acceleration of the body, velocity of the body, pressure force due to vorticity in the flow field, and surface vorticity and friction, gives little information about the phenomena at stake.

Chapter Summary

A historical review of the thermodynamic formulations for the breakdown of drag has been presented. The advantages and limits of such methods have been highlighted, and Van der Vooren's formulation has stood out as the most advanced and practical one. A first generalization attempt, proposed by Gariépy, has also been extensively analyzed, showing deficiency in the phenomenological breakdown. Another family of formulations, based on the velocity vector, has then been investigated. The literature review has shown that unsteady compressible formulations exist, although the few breakdown suggestions, even in the steady case, are not equivalent to the classical wave/viscous/induced decomposition aimed at in this study. These methods therefore do not seem mature enough yet. This study will therefore focus on a new generalization of Van der Vooren's theory to unsteady flows.

Chapter I

Development of an Unsteady Formulation starting from Van der Vooren's Formulation

The aim of this first chapter is to extend Van der Vooren's formulation [START_REF] Van Der Vooren | Drag/Thrust Analysis of a Jet-Propelled Transonic Transport Aircraft: Definition of Physical Drag Components[END_REF] to unsteady flows. This formulation was chosen for its robustness and precision for steady cases (see Chapter Presentation of the main existing methods of numerical drag prediction). It was also the choice of Gariépy [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF]. His formulation however contained a drag coefficient which is not considered a phenomenological component, and the synchronization in time does not seem to be correctly accounted for inside each drag coefficient, since steady formulas have been used. A new careful generalization is therefore required.

In order to generalize the steady formulation to unsteady flows, the obstacles and difficulties to avoid are first identified. The propagation in time and space is particularly looked into.

A new proof of Van der Vooren's theory is then derived by carefully noting all the assumptions required. The first step consists in obtaining the far-field equation. The thermodynamic breakdown is then carried out following the line of reasoning of Van der Vooren, in particular assuming a flow free of vortices in a wake plane. The main contribution of this work comes from the volume splitting, which allows to breakdown the profile drag. We chose to work in a streamtube in order to rigorously identify the assumptions made when reducing the integration surfaces.

Once all the assumptions and mechanisms of the demonstration identified, the formulation can be rigorously extended to unsteady flows. The far-field equation is first generalized to unsteady flows, with the emergence of additional terms. The attribution of these terms into the various phenomenological components has led to several attempts before resulting in the final formulation which has been published in [86]. The key point is to obtain unsteady contributions into each drag component, and to take the wakes into account.

The formulation is then discussed. It is first compared to Gariépy's formulation. The robustness and the physical background are also carefully looked into.

1 What is difficult about a generalization to unsteady flows?

The first step in the generalization of a method is to identify the obstacles. Before trying to obtain the unsteady formulation, let us take a little time to locate the difficulties we shall encounter.

1.1 To account for all the additional terms due to unsteadiness and relate them to phenomenological components

It may seem obvious but in order to be certain not to forget any additional term due to the unsteadiness of the flow, one should start back from the balance of mass and momentum equations in order to obtain the correct far-field equation. Trying to add time derivative terms in the expression of the total far-field drag would lead to a wrong starting point.

Furthermore, these additional terms, once correctly accounted for, must be attributed to existing drag components, or to new ones. There are a few studies in the literature concerning unsteady drag in incompressible multiphase cases, which can be applied to moving airships for example. A so-called virtual mass force, first identified by Prandtl [START_REF] Prandtl | Essentials of Fluid Dynamics[END_REF], and a Basset force [START_REF] Basset | Treatise on hydrodynamics[END_REF] have been defined as resulting from the acceleration of the body.

However, nothing was found in the literature for transonic cases which we are interested in. The first and only team who has tackled the question of unsteady drag breakdown [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF] has chosen to create a so-called unsteady drag coefficient, which is not satisfactory on the physical and phenomenological point of view.

The difficulty will therefore be to identify the physical source of this additional drag. Is it an unsteady contribution into each drag component, such as the contribution of a moving shock in a buffet case? Are there new phenomena which cause drag when the flow is unsteady, such as acoustics? Can it be considered a phenomenological component of drag?

To avoid applying the steady theory as such

We have now settled that the demonstration for the steady theory should be adapted to unsteady flows from the very beginning, that is to say the derivation of the far-field equation. Now another obstacle may arise in the following steps of the demonstration. The assumptions used in the steady theory, which can have been used without being clearly mentioned in the original paper from Van der Vooren and Destarac [START_REF] Van Der Vooren | Drag/Thrust Analysis of a Jet-Propelled Transonic Transport Aircraft: Definition of Physical Drag Components[END_REF], may indeed be inapplicable to unsteady flows. All the assumptions made will therefore have to be pinpointed and adapted to unsteady flows.

Another point is that we need to ponder in the unsteady framework, when we are rather used to thinking in a steady framework. It may for example seem unnatural that a snapshot at a given time of the flow around an object can be sufficient to give the force applied by the fluid on the body at the exact same instant. This is however what the momentum and mass balance equations tell us. The information which allows to account for propagation delays and history of the flow is provided by the time derivative terms as we will see later.

To take the propagation delays into account

The last obstacle consists in taking correctly into account the propagation in time and space, which imply delays in the transmission of perturbations. The history of the flow must also somehow be taken into account to correctly assess the force experienced at a given instant by the body. Remember that the principle of far-field methods consists of integrating phenomena occurring in the flow over a control volume.

A concrete example will make it more clear. Let us consider an airfoil in a flow field, with a shock wave on the upper surface. The viscous effects are assumed to be one order higher than those due to the shock wave and its motion, and are therefore neglected. A streamtube encloses the wave, as depicted in The perturbation information will go on propagating downstream until dissipating in the flow field or leaving the boundaries of our control volume. Now the idea of far-field methods consists in integrating flow field quantities on surfaces or volumes surrounding the sources of drag. In steady cases, it is usual to choose the integration surface for the shock wave close to the wave (the order of a few grid cells). The perturbation will accordingly cross and leave this surface quite soon, and the later propagation will not be accounted for. It however plays a role in the force experienced by the profile at every instant. A means to take this history and propagation process is therefore required when considering unsteady flows.

Derivation of a directly generalizable proof of Van der Vooren's steady formulation

Van der Vooren's formulation allows to split the drag into three phenomenological components: wave, viscous and induced drag, using a thermodynamic breakdown. A spurious contribution is also assessed. It accounts for a part of the numerical drag due to numerical dissipation and can be of great interest to characterize the quality of a mesh or to optimize a configuration with a poor mesh quality. This steady formulation has proven its robustness and precision in many industrial cases and has been developed conjointly by Van der Vooren and Destarac at ONERA for about 20 years. A historical review of the existing drag prediction methods is given in the Chapter Presentation of the main existing methods of numerical drag prediction.

The idea in this section is to give a new derivation of Van der Vooren's theory [START_REF] Van Der Vooren | Drag/Thrust Analysis of a Jet-Propelled Transonic Transport Aircraft: Definition of Physical Drag Components[END_REF], since some steps are missing in the original demonstration and restrictive assumptions regarding the extension to unsteady flows are used. There are three main steps in the method: the derivation of the farfield equation from conservation equations, the thermodynamic breakdown into induced and profile components, and the volume splitting of the profile component into a wave and a viscous component.

The approach adopted here keeps the same process, but differs as far as the volume splitting is concerned. All steps will however be demonstrated in details in order to ensure that no hypothesis is forgotten. Practical refinements will then be presented.

Derivation of the far-field equation in the steady case

The far-field theory consists in computing the aerodynamic force from the flow field analysis instead of the integration of the local stress on the body. The balance between the two approaches relies on the conservation of fluid momentum. It requires no further hypothesis and is thus valid for unsteady compressible flows.

In an inertial reference frame, the conservation of mass in a control volume V with a closed surface ∂V (see Figure I.2) is:

∂V ρ (q • n) dS = 0 (I.2.1)
We can multiply this equation by the free stream velocity q ∞ which is constant:

∂V ρq ∞ (q • n) dS = 0 (I.2.2)
We also have that the integration of the reference pressure p ∞ over the closed surface ∂V is zero:

∂V p ∞ n dS = 0 (I.2.3)
The conservation of momentum in the same control volume also gives:

∂V ρq (q • n) dS + ∂V pn dS - ∂V (τ • n) dS = 0 (I.2.4)
Finally, subtracting Equations I.2.2, and I.2.3 to Equation I.2.4, we get:

∂V ρ (q -q ∞ ) (q • n) dS + ∂V (p -p ∞ ) n dS - ∂V (τ • n) dS = 0 (I.2.5)
In the steady case, this equation gives the conservation of flux momentum at the boundaries ∂V , a closed surface.

The drag is obtained by projecting this vector equation along the free stream velocity direction:

∂V ρ (u -u ∞ ) (q • n) dS + ∂V (p -p ∞ ) (i • n) dS - ∂V (τ x • n) dS = 0 (I.2.6) Figure I.

-Control volume, surfaces and normal vectors

The frontier of the domain ∂V is then split into the body surface S a and the outer surface S e (∂V = S a ∪ S e , see Figure I.2). The body has to be motionless in the inertial reference frame with a steady flow, so that the normal velocity at the skin is zero. Rearranging the terms, we get:

Sa (p -p ∞ ) (i • n) -(τ x • n) dS = Se -ρ (u -u ∞ ) (q • n) -(p -p ∞ ) (i • n) + (τ x • n) dS (I.2.7)
The left-hand side of the equation is the near-field drag D nf , which can be broken down into pressure drag D p and friction drag D f .

D nf = D p + D f (I.2.8) D p = Sa (p -p ∞ ) (i • n) dS (I.2.9) D f = Sa -(τ x • n) dS (I.2.10)
The right-hand side defines the far-field drag, which can be written:

D ff = Se (f • n) dS (I.2.11) with f = -ρ (u -u ∞ ) q -(p -p ∞ ) i + τ x .
This equation represents the flux of momentum through the outer surface and the forces applied on S e . It also gives a property of the vector f which will be used in the following paragraphs, but is valid only for steady flows:

∇ • f = 0 (I.2.12)
No assumption other than classical fluid mechanics flow description has been made so far. Additional terms must however be added in these equations in the unsteady case. Let us first carry on the demonstration in the steady case.

The physical decomposition of drag is then obtained in two steps. The first step consists in breaking down the vector f into two distinct parts, reversible and irreversible. This notion of reversibility is to be understood in the thermodynamic sense, involving assumptions on the entropy variations, rather than in the temporal sense (no flow phenomenon is strictly reversible, in the temporal sense). The second step consists in splitting the control volume V into several volumes attached to the distinct drag sources, in order to breakdown the so-called irreversible drag into wave and viscous contributions. The principles of the decomposition for steady flows are presented in the next sections.

Thermodynamic breakdown

The aim here is to breakdown the far-field expression of drag obtained in Equation I.2.11 into phenomenological components. It was noted that the drag due to the wing-tip vortices, i.e. the induced drag was linked to reversible phenomena. On the contrary, shock waves and viscous interactions are irreversible phenomena, in the sense that entropy is created. The notion of irreversibility commonly used is not strictly rigorous, since enthalpy is also created and taken into account for the profile drag, which is the sum of the wave and viscous drag components. It should rather be called entropic and enthalpic rather than irreversible. We will however stick to the usual terminology for the sake of simplicity.

Breakdown of vector f

Let us first reduce the integration domain for the total far-field drag. Upstream of the body, we have u = u ∞ , p = p ∞ and τ x = 0. On the lateral surfaces sufficiently far or aligned with the streamlines, q • n = 0 and τ x • n = 0, so that

D ff = S d (f • n) dS (I.2.13)
Let us assume that we are in the irreversible case. On a wake plane S d sufficiently far downstream from the drag sources, with a flow submitted to only irreversible processes, it is assumed that:

• p = p ∞ the pressure goes back to the reference pressure

• v = w = 0 the velocity is parallel to the reference velocity q ∞
The assumption concerning the pressure can be justified by taking the y and z components of Equation I.2.5, assuming that the y and z components of the viscous stress tensor are negligible on the wake plane and applying the second assumption (v = w = 0).

We denote u irr the axial velocity on S d under these assumptions. This so-called irreversible axial velocity will be computed in the next section. We can now introduce the assumptions into f in order to obtain the irreversible part of f :

f irr = -ρ (u irr -u ∞ ) q + τ x (I.2.14)
The complementary is the reversible part:

f rev = f -f irr = -ρ (u -u irr ) q -(p -p ∞ ) i (I.2.15)
The profile and induced drags can then be defined as:

D vw = S d (f irr • n) dS = S d -ρ (u irr -u ∞ ) q + τ x • n dS D i = S d (f rev • n) dS = S d -ρ (u -u irr ) q -(p -p ∞ ) i • n dS (I.2.16) (I.2.17)
We now have expressions for the irreversible drag, or profile drag, and the reversible drag, or induced drag. An element of the equations is however missing: the expression of the irreversible axial velocity u irr , which is the axial velocity measured on the wake plane when assuming an irreversible flow.

Derivation of the irreversible axial velocity

The computation of u irr is derived from a general expression of the axial velocity as a function of the thermodynamic variables.

Derivation of the general expression of the full axial velocity

The aim is here to express the full axial velocity as a function of the thermodynamic variables: the stagnation enthalpy H and the entropy s. Owing to the definition of the stagnation enthalpy:

∆H = H -H ∞ = h + q 2 2 -h ∞ + u 2 ∞ 2 (I.2.18)
So that:

q 2 = u 2 ∞ + 2∆H -2 (h -h ∞ ) (I.2.19)
Reminding that q 2 = u 2 + v 2 + w 2 and h = C p T :

u 2 = u 2 ∞ + 2∆H -2C p (T -T ∞ ) -v 2 + w 2 (I.2.20) Dividing by u 2 ∞ : u 2 u 2 ∞ = 1 + 2∆H u 2 ∞ - 2C p T ∞ u 2 ∞ T T ∞ -1 - v 2 + w 2 u 2 ∞ (I.2.21)
We now want to include the entropy in the expression of u. Let us recall that the entropy can be expressed as:

∆s = C p ln   T T ∞ p ∞ p γ-1 γ   (I.2.22)
We can therefore write:

e ∆s Cp = T T ∞ p ∞ p γ-1 γ (I.2.23)
Injecting this expression for T T ∞ into Equation I.2.21 yields:

u 2 u 2 ∞ = 1 + 2∆H u 2 ∞ - 2C p T ∞ u 2 ∞   e ∆s Cp p p ∞ γ-1 γ -1   - v 2 + w 2 u 2 ∞ (I.2.24)
Using the expressions of the heat capacity

C p = rγ γ -1 , the sound velocity a 2 ∞ = γrT ∞ and the Mach number M ∞ = u ∞ a ∞ : u 2 u 2 ∞ = 1 + 2∆H u 2 ∞ - 2 (γ -1) M 2 ∞   e γ-1 γ ∆s r p p ∞ γ-1 γ -1   - v 2 + w 2 u 2 ∞ (I.2.25)
The general expression of the axial velocity is therefore:

u = u ∞ 1 + 2∆H u 2 ∞ - 2 (γ -1) M 2 ∞   e γ-1 γ ∆s r p p ∞ γ-1 γ -1   - v 2 + w 2 u 2 ∞ (I.2.26)
Note that no hypothesis was made in this step. The equations used are simply the definitions of stagnation enthalpy and entropy and are valid in all cases. The derivation of the irreversible axial velocity is now straightforward.

Application of the assumptions to obtain the irreversible expression

The irreversible velocity is defined as the axial velocity measured on the wake plane S d where we apply the assumptions of an irreversible flow, i. The irreversible axial velocity u irr can now be computed plugging these assumptions into Equation I.2.26:

u irr = u ∞ 1 + 2∆H u 2 ∞ - 2 (γ -1) M 2 ∞ e γ-1 γ ∆s r -1 (I.2.27)
This axial velocity depends only on the entropy and the stagnation enthalpy variations from the reference state, hence the name irreversible. Note that this expression is valid only on the wake plane S d , where the entropy and stagnation enthalpy variations, created across the shock waves and within the boundary layers and viscous wakes, are probed. Outside the wake, there are no entropy or enthalpy variations, so that this irreversible velocity becomes the reference velocity u ∞ .

Volume splitting using streamtubes

The irreversible drag (or profile drag) has now to be further broken down into a wave and a viscous drag. This splitting is possible thanks to the location of the sources of drag in distinct zones of the flow, as depicted in Figure I.4. The shock-boundary layer interaction zone is however an exception: both phenomena are there intricate and strongly linked with one another. A choice must therefore be made. It was arbitrarily chosen by Destarac [START_REF] Destarac | Far-field / Near-field Drag Balance and Applications of Drag Extraction in CFD[END_REF] to attribute this contribution to the viscous drag component. The local splitting is then achieved successively considering each source of irreversible drag isolated. The drag computation carried out by Van der Vooren and Destarac [START_REF] Van Der Vooren | Drag/Thrust Analysis of a Jet-Propelled Transonic Transport Aircraft: Definition of Physical Drag Components[END_REF] involves the fact that ∇ • f = 0 for steady flows. Since it is no more true for unsteady cases, the expressions of the drag coefficients were derived using a new methodology relying on a balance over streamtubes. This new argument is presented hereafter.

Wave drag

Let us first assume an inviscid flow with an isolated shock, i.e. without any boundary layers, streamwise vorticity, or wakes. The total drag coefficient is then equal to the wave drag coefficient and the wake integration is reduced to the downstream surface of the streamtube enclosing the shock S wd (see Figure I.5):

D = D w = S wd (f irr • n) dS (I.2.28)
The viscous stress tensor is negligible on the wake plane S wd . The key point is to now legitimize the restriction of the integration to S wd . This restriction is mandatory to split into wave and viscous components. Since we are working in the streamtube enclosing the shock wave, we know that the flow is isentropic and isenthalpic outside of it, and that neither entropy nor stagnation enthalpy was created there. The irreversible velocity is thus equal to the reference velocity on S d \ S wd , as noted in Section 2.2.2.2, giving a zero contribution outside of S wd .

We then try to move the wake surface of integration, S wd to S ′ w , i.e. closer to the source of wave drag. The aim is to be able to process complex configurations with several shock waves or bodies intersecting the wake of the shock, and to improve the numerical accuracy since the cells are usually much smaller closer to the aircraft than in the far-field. Applying the divergence theorem in V wd (volume located between S wd and S ′ w , i.e. downstream of the shock wave as depicted in Figure I.5), given that the lateral contribution is zero since q • n = 0 on a streamline, we get:

D w = S ′ w -ρ (u irr -u ∞ ) (q • n) dS - V wd ∇ • ρ (u irr -u ∞ ) q dV (I.2.29)
Using the continuity equation for steady flows ∇ • (ρq) = 0:

∇ • ρ (u irr -u ∞ ) q = ρ q • ∇ (u irr -u ∞ ) (I.2.30)
Since the reference velocity is constant, it can be taken out of this equation. The irreversible axial velocity depends only on entropy and enthalpy, so that the chain rule gives:

ρ q • ∇u irr = ρ ∂u irr ∂s q • ∇s + ρ ∂u irr ∂H q • ∇H (I.2.31)
It is assumed that the flow is isentropic and isenthalpic downstream of the shock so that both terms on the right hand side are zero. The wave drag is then:

D w = S ′ w -ρ (u irr -u ∞ ) (q • n) dS (I.2.32)
It is as if the entropy and stagnation enthalpy were probed on this surface close to the shock, but one must keep in mind that it really is an integration on a wake plane far from the sources of drag. The local axial velocity is on no account equal to u irr this close to the shock wave. Note that even if this property is here used only for the wave drag, it is also true for the viscous and the induced drag, so that the choice of the integration volumes has theoretically no impact on the results for the drag breakdown, so long as the flow remains isentropic and isenthalpic on the wake plane.

When several shock waves are present on a body, it becomes necessary to turn the wake integral into closed surfaces enclosing the shock waves. Another benefit of using closed surface integrals is the numerical robustness. A closed surface is for example easier to define by the outer boundary of a volume determined using physical criteria.

The wake surface integral is therefore always transformed into a closed surface integral denoted S w as defined in Figure I.6. This manipulation is directly valid since q • n = 0 on the lateral surface of the streamtube, and u irr = u ∞ on the upstream plane (there is no variation of entropy and stagnation enthalpy upstream of the source). The final expression for the wave drag is then: Note that the choice of the small surface S w enclosing the shock wave must be such that the flow is indeed isentropic and isenthalpic downstream of it. The wake surface integral can indeed be moved as far upstream as these assumptions are valid.

D w = Sw -ρ (u irr -u ∞ ) (q • n) dS (I.2.33)
Finally, it must be noted that Van der Vooren original formulation [START_REF] Van Der Vooren | Drag/Thrust Analysis of a Jet-Propelled Transonic Transport Aircraft: Definition of Physical Drag Components[END_REF] used the volume form of this equation. The surface form is however much more robust since the integration limits to a surface enclosing a part of the fluid. The ffd72 code does indeed use the surface form described here.

Viscous drag

We then apply the same reasoning to the viscous drag, assuming a subsonic flow without vortices. Denoting S vd the downstream surface of the streamtube enclosing the body, its boundary layer, and its wake, we get for the viscous drag:

D v = S vd (f irr • n) dS (I.2.34)
Here, the viscous stress tensor could also be neglected on the wake plane S vd similarly to Van der Vooren's article [START_REF] Van Der Vooren | Drag/Thrust Analysis of a Jet-Propelled Transonic Transport Aircraft: Definition of Physical Drag Components[END_REF]. Destarac has however shown that its inclusion into the viscous drag component gave more reliable results [START_REF] Destarac | Far-field / Near-field Drag Balance and Applications of Drag Extraction in CFD[END_REF]. We therefore choose to include it in the unsteady case as well. We could similarly to the wave drag move the surface S vd as close to the trailing edge as possible, that is to say as long as the flow is isentropic and isenthalpic. The viscous stress tensor should nevertheless be neglected in order to do so. We judge that the surface manipulation is not mandatory in the viscous component case, since the downstream boundary could not be moved very close to the trailing edge, given that the assumption of isentropic flow is not valid very close to the body. It actually represents a refinement level that we consider premature here. It could however constitute a further development of the unsteady formulation.

The wake surface is turned into a closed surface S v (see Figure I.6) with the same arguments as in the shock wave case:

D v = Sv -ρ (u irr -u ∞ ) q + τ x • n dS (I.2.35)
Again, in the original formulation by Van der Vooren, the volume form of this equation is used, but the surface one is used in practice at ONERA.

Another justification of the use of streamtubes

The use of streamtubes comes rather naturally when we want to isolate each drag source. Another physical justification arises when we try to compute the force applied by a shock wave on a body in an inviscid flow. The total drag, equal in that case to the wave drag, is indeed expressed as the integral of the vector f on the skin of the body in the near-field vision:

D = D w = - Sa (f • n) dS (I.2.36)
A balance of the vector f over the boundary of the streamtube enclosing the body Ω a gives:

∂Ωa (f • n) dS = 0 (I.2.37)
In 2-D, the boundary of Ω a can be split into S au , Σ 2 , S ad , Σ 3 , and S a (as defined in Figure I.7). A similar reasoning can be applied in 3-D.

Sau (f • n) dS - Σ 2 (f • n) dS + S ad (f • n) dS + Σ 3 (f • n) dS + Sa (f • n) dS = 0 (I.2.38)
Since the flow is inviscid without vortices, the pressure and the velocity on S au and S ad are equal to their reference values p ∞ and u ∞ , and the friction is negligible. The contributions on the upstream and downstream boundaries of Ω a are therefore zero. The same balance can be made on the streamtube Ω 2 , with zero contributions everywhere, so that we can deduce that the contribution on Σ 3 is also zero. Finally: Let us now apply the same reasoning on the streamtube Ω w enclosing the shock wave:

D w = - Σ 2 (f • n) dS (I.2.39)
∂Ωw (f • n) dS = 0 (I.2.40)
Splitting the frontier:

Swu (f • n) dS + Σ 1 (f • n) dS + S wd (f • n) dS + Σ 2 (f • n) dS = 0 (I.2.41)
For the same reason as before, the contribution on the upstream boundary S wu is zero. We can again apply the same balance on the streamtube above, Ω 1 , deducing that the contribution over Σ 1 is also zero. Finally:

D w = S wd (f • n) dS (I.2.42)
If S wd is chosen such that the assumptions p = p ∞ and v = w = 0 apply, we obtain the starting point from Section 2.3.1, which is:

D w = S wd (f irr • n) dS (I.2.43)
This proof also leads to a physical interpretation of the wave drag: due to the creation of entropy across the shock wave and conservation of mass flux, the streamtube downstream of the shock is enlarged, resulting in an obstruction of the streamtube enclosing the body itself. A pressure force thus arises which is through the fluid applied on the body skin. Its longitudinal component is what we call the wave drag.

Derivation of the final steady formulation

Now that all elements are defined, the final formulation can be obtained. A "raw" version is first given. Refinements were carried out over the years and resulted into a robust and reliable postprocessing tool used by industry, the ONERA code ffd72.

A first "raw" formulation

As far as induced drag is concerned, it is defined as the complementary drag, i.e. the integral of the reversible part of f . This definition by default rather than by using physical assumptions is the main flaw of Van der Vooren's theory. Other physical phenomena can indeed be taken into account in this drag component.

The last step consists in computing the far-field drag as the sum of the three drag components. The resulting formulation is theoretically equivalent to the one exposed in Van der Vooren and Destarac's paper [START_REF] Van Der Vooren | Drag/Thrust Analysis of a Jet-Propelled Transonic Transport Aircraft: Definition of Physical Drag Components[END_REF], the only difference being the use of surface integrals rather than volume ones.

D w = Sw -ρ (u irr -u ∞ ) (q • n) dS D v = Sv -ρ (u irr -u ∞ ) q + τ x • n dS D i = Se -ρ (u -u irr ) q -(p -p ∞ ) i • n dS D ff = D w + D v + D i (I.2.44) (I.2.45) (I.2.46) (I.2.47)
The formulation was implemented in the early version of the drag extraction codes family ffd. Deviations from the theory were then observed, and refinements were proposed by Destarac [START_REF] Destarac | Far-field / Near-field Drag Balance and Applications of Drag Extraction in CFD[END_REF] so as to result in a reliable and accurate tool.

Numerical deviations from the theory

The far-field drag as defined in this formulation is theoretically equal to the near-field drag. There are in practice numerical deviations from the theory.

Spurious drag

The first deviation is the spurious drag due to the assumptions made when deriving the equations. This physical spurious drag component, denoted D phy sp , can be evaluated from the neglected term, that is to say the irreversible part of f in the volume outside S w and S v , denoted V sp :

D phy sp = Vsp (∇ • f irr ) dV + Se \S vd (τ x • n) dS (I.2.48)
The spurious drag can also originate from numerical effects. The quality of the mesh has indeed an influence on the far-field drag computed from integrations in the fluid domain, whereas the nearfield drag depends only on the mesh quality at the skin. Moreover, the numerical scheme can induce additional spurious diffusion which can also imply further discrepancy between far-field and near-field values. Moreover, errors may arise from the boundary conditions chosen. In 2-D, vortex type boundary conditions can be applied to reduce these errors. Note finally that the near-field drag is also affected by numerical errors. These errors are however often small compared to the errors in the far-field coefficient, and will not be accounted for in this method.

An estimation of the numerical spurious drag D num sp can be defined. Its exact expression cannot be directly obtained. We can however compute the total spurious drag as the balance between the near-field and the far-field, neglecting the errors on the near-field coefficient:

D sp = D nf -D ff (I.2.49)
The numerical spurious drag is then defined as:

D num sp = D nf -D ff -D phy sp (I.2.50)
This is a strong benefit of the far-field formulation: it allows extracting the spurious drag, which stands for a part of the numerical errors (with the restriction that the near-field drag must be sufficiently accurately computed). This spurious drag can be useful to evaluate the quality of a mesh, the accuracy of a numerical scheme, or in an optimization process to optimize a physical component rather than the numerical errors, as illustrated by Destarac [START_REF] Destarac | Far-Field Drag in Transonic Potential Flow: Analysis and Optimisation[END_REF].

Vortex decay

Another discrepancy between the theoretical and the numerical behavior of the flow concerns the decay of the trailing vortices downstream from finite lifting surfaces. In computations, mainly because of the coarseness of the grid in the downstream far-field, and to the spurious diffusion due to the numerical scheme, numerical smoothing dominates over physical dissipation and causes the trailing vorticity to decay. Induced drag expressed through a surface integral as in Equation I.2.46, will decrease as the downstream extension of S e increases. Such losses of apparent induced drag (which are also a loss of crossflow kinetic energy) due to a spurious dissipative phenomenon will be matched by the production of spurious viscous drag. There is neither loss nor production of total drag, but a transfer of one form of drag to another.

Practical refinements of the theoretical formulation

Many refinements of the raw formulation were carried out over the years at ONERA in order to develop a robust tool.

The single vector formulation

Destarac [START_REF] Destarac | Far-field / Near-field Drag Balance and Applications of Drag Extraction in CFD[END_REF] relied on the fact that for steady flows, we have ∇ • f = 0 to express all drag coefficients as a function of only one part of vector f = f irr + f rev . He observed that the results were losing accuracy and robustness when using only f irr , without being able to explain the phenomenon. It would seem that the physical equation ∇ • f = 0 is numerically not accurately verified. Moreover, the numerical behavior of f irr and f rev is different. The numerical errors in ∇ • f seem to mainly originate from f irr .

As a result, Destarac chose to use a formulation using only the reversible part f rev . The expressions of each drag coefficient have first to be turned into volume integrals to do so, using the divergence theorem:

D w = Vw (∇ • f irr ) dV (I.2.51) D v = Vv (∇ • f irr ) dV - Sa (f irr • n) dS (I.2.52
)

D i = V (∇ • f rev ) dV - Sa (f rev • n) dS (I.2.53)
The body surface integrals simply give the near-field drag coefficients (since q • n = 0 at the skin):

Sa (f irr • n) dS = Sa -ρ (u irr -u ∞ ) q + τ x • n dS = Sa (τ x • n) dS = -D f (I.2.54) Sa (f rev • n) dS = Sa -ρ (u -u irr ) q -(p -p ∞ ) i • n dS = - Sa (p -p ∞ ) (i • n) dS = -D p (I.2.55)
Now, using ∇ • f = 0 and f = f irr + f rev , we can write:

D w = - Vw (∇ • f rev ) dV (I.2.56) D v = - Vv (∇ • f rev ) dV + D f (I.2.57) D i = V (∇ • f rev ) dV + D p (I.2.58)
We can then apply again the divergence theorem to go back to surface integrals which allow integrating on a smaller domain:

D w = - Sw (f rev • n) dS (I.2.59) D v = - Sv (f rev • n) dS - Sa (f rev • n) dS + D f (I.2.60
)

D i = Se (f rev • n) dS + Sa (f rev • n) dS + D p (I.2.61)
So that finally:

D w = - Sw (f rev • n) dS D v = - Sv (f rev • n) dS + D nf D i = Se (f rev • n) dS (I.2.62) (I.2.63) (I.2.64)
This one-vector formulation has the advantages of being both faster to compute and more accurate. In particular, the viscous stress coefficient, a second order derivative subject to accuracy difficulties, is required only at the skin and not in the far-field.

The decomposition of the viscous drag coefficient

When the one-vector formulation is used, a further breakdown of the viscous drag coefficient is achievable. The friction drag has indeed appeared in Equation I.2.60 and can be removed to define the viscous pressure drag coefficient.

D vp = - Sv (f rev • n) dS + D p (I.2.65)
So that:

D v = D vp + D f (I.2.66)
The viscous pressure drag component is the part of the viscous drag which is reducible, for instance by minimizing flow separation through shape modifications. The friction drag component is indeed only due to the viscous stress occurring at the skin, and is almost impossible to reduce, except by using riblets on the skin of the aircraft [START_REF] Bechert | Experiments on Drag-Reducing Surfaces and their Optimization with an Adjustable Geometry[END_REF].

Criteria used for the definition of the integration surfaces

The theory as derived here requires the definition of streamtubes as integration surfaces for the wave and viscous drag components. The induced surface is in the raw version defined as the outer surface of the control volume. In practice, it has been observed that smaller surfaces imply better robustness. It also has consequences on the spurious drag coefficient, which has been observed to be smaller with smaller surfaces. Moreover, streamtubes are not easily constructed in a postprocessing tool. Physical criteria are consequently used in practice, leading to small discrepancies with the theoretical definitions. The corresponding surfaces will be denoted with a superscript p for practice.

Wave drag criterion

We have seen with our definition that the downstream face of the streamtube enclosing the shock could be moved close to the shock wave. Since the surface is quite small around the shock wave, and since the contributions over lateral surfaces will remain small whatever their shape, a physical criterion of shock detection is in practice used to define the integration surface.

This criterion was first used by Tognaccini [START_REF] Tognaccini | Methods for Drag Decomposition, Thrust-Drag Bookkeeping from CFD Calculations[END_REF]. It allows to detect the presence of a shock wave from the comparison between the local velocity in the direction of the pressure gradient and the speed of sound: if this velocity is greater than the speed of sound, then the cell must be inside the shock integration surface. In practice, a coefficient c shock close to 1 is added. It is chosen equal to 0.95 in most cases but can be tuned by the user.

q • ∇p ≥ c shock a ∇p ⇒ cell ∈ S p w (I.2.67)
Defaults values and the number of cell layers added in ffd72 can be found in the Chapter Methods used in this work .

Viscous drag criteria

In the theoretical demonstration, the surface chosen for the viscous drag was the streamtube enclosing the body and its boundary layer. The downstream extension can be reduced similarly to the shock wave case, provided that the viscous dissipation is negligible in the wake. Again, the streamtube is replaced by a physical contour based on viscous stress and turbulent considerations. The greatest contribution to drag comes from the downstream plane again, since the lateral surfaces are almost aligned with the streamlines where the integration of the normal components cancel out. The contribution of the viscous stress and momentum flux on the lateral boundaries is actually very small and considered negligible compared to the contribution of the downstream boundary.

The first physical criterion consists in comparing the local viscous stress tensor to the mean value at the skin τ skin . This allows detecting the boundary layers. If the local value is stronger than a small percentage of the value at the skin, then the cell must be inside the viscous integration surface. This percentage c skin close to 10% is again used for tuning.

τ ≥ c skin τ skin ⇒ cell ∈ S p v (I.2.68)
The second criterion regards the turbulent level in the flow in the vicinity of the skin, in order to detect the boundary layers as well as the wakes. The eddy viscosity is compared to the dynamic viscosity: if it is stronger, then the cell must be located within the viscous surface. A tunable coefficient c turb close to 1 is also used.

µ t ≥ c turb µ ⇒ cell ∈ S p v (I.2.69)
Cell layers are also added for robustness, and the downstream extension can be tuned. See Chapter Methods used in this work for details and default values.

Induced drag criteria

In the theory, the integration surface for the induced drag is the outer surface of the control volume. However, as noted in Section 2.4.2, it was observed that the induced drag prediction was more accurate when integrating closer to the trailing edges of the body due to the numerical vortex decay. Geometrical criteria are therefore used to define a box around the airplane, with a downstream extension a given distance from the wing tips. The frontiers visible inside the surfaces are simply due to the frontiers between the different grid types, since the case is a chimera mesh. The box for the induced drag integration stops at a given distance from the most downstream position of the wings trailing edge. The tail of the aircraft is also included in the surface, with a margin chosen as a given distance from the skin. All these refinements put together lead to the formulation as implemented in the ONERA code ffd72. This postprocessing software is presented in details in the Chapter Methods used in this work .

Examples of integration surfaces

The equations can now be written as: with the integration surfaces as described earlier. The drag coefficients are defined as the nondimensional version of the drag components described here:

D w = S p w ρ (u -u irr ) q + (p -p ∞ ) i • n dS D vp = S p v ρ (u -u irr ) q + (p -p ∞ ) i • n dS + D p D v = D vp + D f D i = S p i -ρ (u -u irr ) q -(p -p ∞ ) i • n dS D ff = D w + D v + D i D sp = D nf -D ff (I.
CD = D 1 2 ρ ∞ u 2 ∞ S ref (I.2.76)
where S ref is the reference surface, usually chosen as the planform area of the wing of the aircraft. This formulation has proved its robustness and accuracy over the years. Many extensions have also been added, such as propulsive configurations, relative rotating reference frames, or the computation of sensitivities for optimization. It is however restricted to steady flows.

Generalization of Van der Vooren's formulation to unsteady flows

The extension to unsteady flows of Van der Vooren's theory is not as straightforward as it seems. Where some authors [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF] chose to define an unsteady drag coefficient concentrating all unsteady contributions, we have finally decided to allocate unsteady terms to each drag component, in order to take the early drag creation and propagation in the fluid into account. The main difficulty tackled here is indeed the synchronization of the fluid phenomena and the loading experienced by the body at the same instant.

In order to do so, we have derived a new proof of Van der Vooren's formulation and carefully located all the assumptions which have been made:

• The unsteady terms in the balance of mass and momentum equations were removed.

• During the thermodynamic breakdown step, under the hypothesis of an only irreversible flow, it was assumed that the pressure went back to the reference pressure and the velocity was parallel to the reference velocity in a downstream wake plane.

• At the volume splitting step, it was assumed that the flow was not affected by the presence of the shock wave or the boundary layer outside the streamtube enclosing the source.

• When moving the integration surfaces, it was assumed that the flow was isentropic and isenthalpic downstream of a shock wave.

The generalization is now made achievable. We only have to adapt all the assumptions to unsteady flows, starting with the first step: the derivation of the far-field equation.

Implementation of the additional unsteady terms in the far-field equation

The derivation of the far-field equation follows the same principle as in the steady case, with additional terms into each conservation equation due to the unsteadiness. Note that all physical variables are now time dependent. We will therefore write the conservation equations at a given instant t.

In an inertial reference frame, the conservation of mass in a control volume V with a closed surface ∂V is:

V ∂ρ ∂t dV = - ∂V ρ (q • n) dS (I.3.1)
We can multiply this equation by the free stream velocity q ∞ which is constant:

V ∂ρq ∞ ∂t dV = - ∂V ρq ∞ (q • n) dS (I.3.2)
As in the steady case, the integration of the reference pressure p ∞ over the closed surface ∂V is zero:

∂V p ∞ n dS = 0 (I.3.3)
The conservation of momentum in the same control volume also gives:

V ∂ρq ∂t dV = - ∂V ρq (q • n) dS - ∂V pn dS + ∂V (τ • n) dS (I.3.4)
Finally, subtracting Equations I.3.2, and I.3.3 to Equation I.3.4, we get:

V ∂ρ (q -q ∞ ) ∂t dV = - ∂V ρ (q -q ∞ ) (q • n) dS - ∂V (p -p ∞ ) n dS + ∂V (τ • n) dS (I.3.5)
This equation gives the balance between the variation of fluid momentum in a control volume V with the aerodynamic force and the momentum flux at the boundaries ∂V , a closed surface, at instant t.

The drag is obtained by taking the component parallel to the free stream velocity:

V ∂ρ (u -u ∞ ) ∂t dV = - ∂V ρ (u -u ∞ ) (q • n) dS - ∂V (p -p ∞ ) (i • n) dS + ∂V (τ x • n) dS (I.3.6)
Splitting the frontier of the domain ∂V into the body surface S a and the outer surface S e (∂V = S a ∪ S e , see Figure I.2) and rearranging the terms, we get:

Sa (p -p ∞ ) (i • n) -(τ x • n) dS = Se -ρ (u -u ∞ ) (q • n) -(p -p ∞ ) (i • n) + (τ x • n) dS - Sa ρ (u -u ∞ ) (q • n) dS - V ∂ρ (u -u ∞ ) ∂t dV (I.3.7)
The left-hand side of the equation still defines the near-field drag at instant t in the exact same way as in the steady case:

D nf (t) = Sa (p -p ∞ ) (i • n) -(τ x • n) dS (I.3.8)
The right-hand side defines the far-field drag at instant t, which can be written:

D ff (t) = Se (f • n) dS - Sa ρ (u -u ∞ ) (q • n) dS - V ∂ρ (u -u ∞ ) ∂t dV (I.3.9)
The first surface term represents the flux of momentum through the outer surface and the forces applied on S e similarly to the steady case. The second surface term is the variation of longitudinal momentum due to the motion of the body. It is zero in steady cases in the inertial reference frame. The volume term accounts for the time dependence as well as the propagation in time of momentum.

The notion of far-field can be considered abusive, given that it now includes a contribution at the skin and in the fluid volume which extends from the skin to the far-field. The designation of far-field drag is however deliberately kept as in opposition to near-field drag which consists only of an integration at the skin.

These equations also give a property of the vector f which will be used in the following sections:

∇ • f = ∂ρ (u -u ∞ ) ∂t (I.3.10)
Note that this result is quite different from the steady case, where we had ∇ • f = 0. This justifies the choice of deriving the splitting step using streamtubes rather than the property of the vector f . Before trying to obtain the breakdown in the unsteady case, let us remind the main objective of our work: we will focus on ensuring the synchronization between the integration in the far-field and the effort experienced by the body at the same instant. The volume time derivative term will be playing a key role.

Derivation of the four components unsteady formulation

The unsteady formulation is achieved using the new demonstration of the steady theory given in Section 2. The presentation of this breakdown method is an opportunity of giving the rigorous demonstration.

Like for the steady case, the wave drag D w is first defined for a flow where we assume that only one shock occurs. This approach is generalized to the viscous drag D v . The unsteady induced drag, D ui , is defined as the complementary part. It can obviously enclose other unsteady phenomena and is not only due to vorticity variations, hence the denomination unsteady induced.

Unsteady wave drag expression

Let us first consider the wave drag of an isolated normal shock moving in a perfect fluid. It is assumed that there are no other processes within the flow like boundary layers or vortices. We will again work in the streamtube enclosing the shock as defined for the steady case (see Figure I.9). We choose for the downstream boundary S wd such that p = p ∞ and v = w = 0. We obtain the same irreversible velocity as in Equation I.2.27.

If the flow is supersonic, there is no information propagating upstream of the shock wave, so that the integration can be reduced to the downstream part of the streamtube. In the subsonic case, however, there can be information ascent, e.g. acoustic waves. We assume that the upstream contribution is negligible. The neglected part will be assigned to the unsteady induced drag component as we will see in Section 3.2.3. As a conclusion, we limit the integration to the downstream part of the streamtube, that is to say V w ∪ V wd . The time-derivative term now adds to the steady terms. The wave drag is indeed by definition caused by entropy creation and volume effects downstream of the shock. The total drag, which is equal to the wave drag, is therefore:

D = D w = S wd -ρ (u irr -u ∞ ) (q • n) dS - Vw ∪V wd ∂ρ (u -u ∞ )
∂t dV (I.3.11) The wave drag is the sum of the flux on S wd of loss of axial momentum due to the irreversible process (shock) and the variation in time of axial momentum in the streamtube.

Like for steady flows, the flux integrated on S wd is transformed into the same flux but integrated on S ′ w (see Figure I.9). The aim is here to probe the entropy due solely to the shock wave. The entropy probed on S wd contains indeed a viscous component which we are trying to get rid of. The application of the divergence theorem in V wd gives:

D w = S ′ w -ρ (u irr -u ∞ ) (q • n) dS - Vw ∂ρ (u -u ∞ ) ∂t dV - V wd ∂ρ (u -u ∞ ) ∂t + ∇ • ρ (u irr -u ∞ ) q dV (I.3.12)
The last term, I, can be written:

I = - V wd ∂ρ (u -u irr ) ∂t dV - V wd ∂ρ (u irr -u ∞ ) ∂t + ∇ • ρ (u irr -u ∞ ) q dV (I.3.13)
We can expand the last two terms J under the integral as:

J = ρ ∂ (u irr -u ∞ ) ∂t + ρq • ∇ (u irr -u ∞ ) + (u irr -u ∞ ) ∂ρ ∂t + (u irr -u ∞ ) ∇ • (ρq) (I.3.14)
The last two terms cancel out due to the continuity equation ∂ρ ∂t + ∇ • (ρq) = 0: A local balance of internal energy in the absence of heat flux gives:

J = ρ d (u irr -u ∞ ) dt =
ρ de dt = ∇ • (σ • q) -∇ • σ • q (I.3.17)
On the other hand, a local balance of kinetic energy gives:

ρ d dt u 2 2 = ∇ • σ • q (I.3.18)
so that ∇ • σ • q cancels out. Given the definition of the stress tensor σ = -p1 + τ , it can be written that:

∇ • (σ • q) = -∇p • q + ∇ • (τ • q).
The definition of the time derivative yields for the pressure:

dp dt = ∂p ∂t + ∇p • q (I.3.19)
so that the term ∇p • q cancels out as well. As a result, the time derivative of stagnation enthalpy can be expressed as simply:

dH dt = 1 ρ ∂p ∂t + ∇ • (τ • q) (I.3.20)
Assuming that the viscous effects are negligible in the wake of a shock wave, we get:

dH dt = 1 ρ ∂p ∂t (I.3.21)
The only term left to compute is the partial derivative of the irreversible axial velocity with respect to stagnation enthalpy:

∂u irr ∂H = 2 u 2 ∞ u ∞ 2 1 + 2∆H u 2 ∞ - 2 (γ -1) M 2 ∞ e γ-1 γ ∆s r -1 = 1 u irr (I.3.22)
Finally:

I = - V wd ∂ρ (u -u irr ) ∂t + 1 u irr ∂p ∂t dV (I.3.23)
Concerning the surface contribution, it can be turned into a closed surface integral S w similarly to the steady case, assuming that u irr = u ∞ on the upstream surface, and q • n = 0 on the lateral surfaces aligned with the streamtube.

As a result, the unsteady wave drag definition is:

D w = Sw -ρ (u irr -u ∞ ) (q • n) dS - Vw ∂ρ (u -u ∞ ) ∂t dV - V wd ∂ρ (u -u irr ) ∂t + 1 u irr ∂p ∂t dV (I.3.24)
It must be noted that for a steady flow this equation gives the classical steady wave drag. The unsteady wave drag is the sum of three elements. The two first stand for the instantaneous creation of drag due to the entropy creation across the shock wave in the thin volume V w . The first integral is the flux of momentum and the second one the volume creation. The third element takes into account the propagation and the past effects of entropy production downstream of the shock. The fluid is assumed inviscid for the latter, so that it only stands for the propagation of the past creation of entropy in the volume V w . It is obvious that irreversibilities take place downstream of the shock. The drag due to these physical phenomena is small compared to D w . This contribution may appear in the viscous drag coefficient if the volume V v is sufficiently large, or else in the induced drag coefficient.

Unsteady viscous drag expression

We now consider a profile isolated in a flow, without shock waves and reversible processes. We work again in the streamtube enclosing the body, its boundary layer and its wake. We again neglect the upstream part of the streamtube and work in the volume V v (see Figure I.10). The time-derivative term also adds to the steady term:

D v = Sv -ρ (u irr -u ∞ ) q + τ x • n dS - Vv ∂ρ (u -u ∞ ) ∂t dV (I.3.25)
One could also move the wake surface a little upstream similarly to the shock case, the same volume term would appear, but we do not judge it necessary here since irreversibilities occur within the wake itself, so that the wake integral would remain large. In order to reduce numerical errors due to the dissipation of the entropy downstream of the body, the control domain can however be bound at a given distance of the trailing edge, as long as a sufficient part of the wake is included.

Unsteady induced drag expression

As for the steady case, the induced drag is here defined as the complementary part of the drag. Since it will likely enclose other unsteady phenomena, we decided to call it unsteady induced drag D ui rather than simply induced drag. The notion of lift induced is kept in order to highlight the strong connection between this component and the variations of lift.

We define the complementary volume

V c = V \ (V w ∪ V wd ∪ V v )
D ui = S d -ρ (u -u irr ) q -(p -p ∞ ) i • n dS + S cd -ρ (u irr -u ∞ ) (q • n) dS - Vc ∂ρ (u -u ∞ ) ∂t dV (I.3.26)
The second term has the same structure as the original expression of wave drag in Equation I.3.11, given that S cd is the downstream boundary of V c . The idea comes naturally to simplify this expression using the same process as in the wave drag case, that is to say to move the downstream integration S cd towards the upstream infinity S u , the same volume terms arising since the flow has the same properties in the complementary volume ds dt = 0 and dH dt = 1 ρ ∂p ∂t . We can thus write: Similarly to the wave drag case, this manipulation enables to reduce the numerical errors by ensuring the isentropic feature of the flow.

S cd -ρ (u irr -u ∞ ) (q • n) dS - Vc ∂ρ (u -u ∞ ) ∂t dV = Su -ρ (u irr -u ∞ ) (q • n) dS - Vc ∂ρ (u -u irr ) ∂t + 1 u irr ∂p ∂t dV (I.3.27)
The integration on the upstream surface S u gives a zero contribution since u irr = u ∞ there. Using the usual assumption that the contributions on the lateral surfaces are zero, the wake integral can be turned into a closed surface integral over the outer surface. The final expression for the induced drag component is therefore:

D ui = Se -ρ (u -u irr ) q -(p -p ∞ ) i • n dS - Vc ∂ρ (u -u irr ) ∂t + 1 u irr ∂p ∂t dV (I.3.28)
If the fluid is inviscid without any shock, then V d becomes the whole volume and thus the whole drag is taken into account. It must be noted that this expression of induced drag takes into account, within the volume term, phenomena which can be part of the irreversible drag, such as propagation of viscous interactions outside the viscous volume. We will focus in the applications on analyzing its behavior.

Final decomposition

The spurious drag is defined as the balance between the near-field and the far-field approaches, similarly to the steady case. The final decomposition is then: Note that if we remove all the time-derivative terms, we obtain the steady formulation as expressed in Section 2. This unsteady formulation has been published in [START_REF] Toubin | Far-Field Drag Decomposition for Unsteady Flows[END_REF], [84] and [86]. The definition of the surfaces and volumes is shown in Figure I.10. Again, the drag coefficients are defined as the non-dimensional version:

D w = Sw -ρ (u irr -u ∞ ) (q • n) dS - Vw ∂ρ (u -u ∞ ) ∂t dV - V wd ∂ρ (u -u irr ) ∂t + 1 u irr ∂p ∂t dV D v = Sv -ρ (u irr -u ∞ ) q + τ x • n dS - Vv ∂ρ (u -u ∞ ) ∂t dV D ui = Se -ρ (u -u irr ) q -(p -p ∞ ) i • n dS - Vc ∂ρ (u -u irr ) ∂t + 1 u irr ∂p ∂t dV D m = Sa -ρ (u -u ∞ ) (q • n) dS D ff = D w + D v + D i + D m D sp = D nf -D ff (I.
CD = D 1 2 ρ ∞ u 2 ∞ S ref (I.3.35)
Since the whole volume is here taken into account, the spurious drag should come from both the numerical errors and the assumptions made during the demonstration. There are three:

• The flow was assumed isentropic in the complementary volume and in the wake of the shock.

• The viscous contribution in the equation for the stagnation enthalpy variation was neglected in the complementary volume and in the wake of the shock.

• The longitudinal viscous stress vector was neglected outside the viscous volume.

The neglected terms can consequently be considered as a physical definition of the spurious drag, to which the numerical errors have to be added. Another definition of the spurious drag is then:

D sp = D phy sp + D num sp (I.3.36) D phy sp = - Vc ∪V wd ρ ∂u irr ∂∆s d∆s dt + ∂u irr ∂∆H ∇ • (τ • q) dV + Se \S vd (τ x • n) dS (I.3.37) D num sp = D nf -D ff -D phy sp (I.3.38)

Criteria used in practice for the integration volumes definition

We have seen in the theory that the equations are valid on streamtubes enclosing the shock wave or the body and its boundary layer. In practice, such a surface is not trivial to compute, especially since the postprocessing tool takes only the conservative variables in the fluid domain as inputs. For unsteady turbulent flow simulated by DES methods in particular, such a streamtube can even be illdefined if it is intersected by turbulent structures. An almost equivalent, more practical, definition has thus been found.

The idea is to find a physical criterion which allows to make the distinction between the shock wave and its wake and the boundary layer and the viscous wake. The variables which allow this distinction are the entropy and the vorticity:

• Entropy is created across the shock wave and is advected downstream in the wake of the shock, following the streamlines. Assuming a regular structure for the shock wave, we can say that there is no or very little vorticity in the wake of a shock.

• Entropy and vorticity are created in the boundary layers by the friction between fluid particles. Both are advected along the wake as well.

Critical values are user defined for the entropy and the vorticity. The steady criteria are also used since refinements were made over the years to ensure the continuity of the surfaces and a sufficient distance from the body (with cell layers added for example). The deviatoric viscous stress tensor is also used to avoid holes in the volumes. The routine for the definition of the new surfaces is the following:

• V w : same criterion as in the steady case (see Section 2.4.3.3).

• V wd : five criteria are used. In the following, -sign stands for a logical and. A cell is tagged as inside V wd if:

-outside V w -outside V p v (volume inside S p v , defined in the steady case in Section 2.4.3.3) -∆s ≥ c entropy -ω ≤ c vorticity -τ ≤ c f riction τ m • V v : three criteria are used. A cell is tagged as inside V v if: -outside V w -inside V p v or -∆s ≥ c entropy or -∆s ≥ c entropy -ω ≥ c vorticity -τ ≥ c f riction τ m
• V c : no refinement on the induced volume is used. The surface S e is the outer surface of the control domain, and V c is the complementary volume of the whole control domain. The downstream extension of the control volume can however be tuned.

Examples of the resulting integration volumes are shown in Figure I.12. We can see that the wake of the shock wave is effectively included in V wd . The viscous volume extends also further downstream than in the steady case. Finally, the complementary volume used for the integration of the unsteady induced drag is as large as the control volume.

When there is no shock wave, V p v is used for the unsteady viscous drag. Other refinements similar to the steady formulation, such as a one vector formulation or the definition of a viscous pressure component, could be made. The method is however not considered mature enough to address such advanced improvements. This new unsteady formulation, although derived with as much care as possible, may lack robustness. It is therefore crucial to discuss its validity and the assumptions made during its derivation. A main flaw, which is the physics behind the unsteady induced drag component, also requires attention. Finally, it is natural to compare the formulation with the only other one available in the literature.

Robustness of the formulation

Let us now look into the robustness of the formulation. Two main possible issues are identified: the domain of definition of the irreversible axial velocity and the physical criteria used for the definition of the integration volumes.

Domain of definition of the irreversible axial velocity

The first source of loss of robustness is the use of the irreversible axial velocity u irr . As first pointed out by Méheut [START_REF] Méheut | Drag-Breakdown Methods from Wake Measurements[END_REF], this variable can indeed be undefined in certain zones of the flow. Its expression comes from the application of a flow free of vortices on a wake plane far from the sources. The irreversible velocity is there always defined. However, we have seen in Section 2.3 that this velocity was in practice probed on surfaces closer to the sources of drag. It is therefore required to check whether its expression is still valid there. We need to check whether the expression under the square root, denoted Γ irr , can become negative:

Γ irr = 1 + 2∆H u 2 ∞ - 2 (γ -1) M 2 ∞ e γ-1 γ ∆s r -1 (I.4.1)
The entropy and stagnation enthalpy variations can also be written as a function of the stagnation temperature and pressure:

∆s = C p ln   T i T i∞ p i∞ p i γ-1 γ   (I.4.2) ∆H = C p (T i -T i∞ ) (I.4.3)
So that:

Γ irr = 1 + 2C p u 2 ∞ (T i -T i∞ ) - 2 (γ -1) M 2 ∞   T i T i∞ p i∞ p i γ-1 γ -1   (I.4.4)
We have seen earlier that

C p u 2 ∞ = 1 (γ -1) M 2 ∞ T ∞
, which gives:

Γ irr = 2 (γ -1) M 2 ∞   (γ -1) M 2 ∞ 2 + 1 T ∞ (T i -T i∞ ) - T i T i∞ p i∞ p i γ-1 γ + 1   (I.4.5)
The definition of the stagnation temperature gives

T i∞ = T ∞ 1 + (γ -1) M 2 ∞ 2
, so that:

Γ irr = 2 (γ -1) M 2 ∞   T i∞ -T ∞ T ∞ + T i -T i∞ T ∞ - T i T i∞ p i∞ p i γ-1 γ + 1   (I.4.6) Γ irr = 2 (γ -1) M 2 ∞   T i T ∞ - T i T i∞ p i∞ p i γ-1 γ   (I.4.7)
The expression for the irreversible velocity u irr is hence defined if:

Γ irr ≥ 0 ⇔ p i∞ p i γ-1 γ ≤ T i∞ T ∞ = p i∞ p ∞ γ-1 γ ⇔ p i ≥ p ∞ (I.4.8)
This formulation is therefore undefined if p i < p ∞ . Recall that this expression has been devised on a wake plane far from the sources. It is thus not surprising that it may become undefined in other zones of the flow. This situation can happen in the regions where the fluid is detached, in the core of strong vortices, in the boundary layers if M > M ∞ locally and downstream of strong shocks, which means inside almost any kind of drag source.

Another physical explanation can be found from Section 2.3.3. The ratio between the upstream and downstream boundaries of the streamtube S wu and S wd can indeed be expressed as:

S wu S wd = ρ d ρ u u d u u = e γ-1 γ ∆s r 1 + 2∆H u 2 ∞ - 2 (γ-1)M 2 ∞ e γ-1 γ ∆s r -1 (I.4.9)
p i becoming smaller than p ∞ corresponds to an entropy variation becoming large, so that the downstream boundary S wd becomes large as well. We can therefore interpret the non definition of the irreversible axial velocity defect as a complete obstruction of the streamtube enclosing the body, as depicted in Figure I. [START_REF] Brunton | Unsteady Aerodynamic Forces on Small-Scale Wings: Experiments, Simulations and Models[END_REF].

For steady cases, the zones where u irr is undefined can be easily enclosed within S p v or S p w , resulting in small discrepancies on the drag coefficients.

The physics for unsteady flows may however lead to issues. The unsteady flow could indeed be such that the non-definition zones become much larger and evolve in time. Including them into S p v would then lead to absurd contours, or even impossible separation between the viscous and the wave surfaces. Not including them involves to choose a default value to these cells, leading to a loss of accuracy. An alternative expression would therefore be the best solution, as long as the accuracy and robustness is kept.

Physical criteria used for the definition of the integration volumes

The criteria used for the definition of the integration domains S w and S p v recalled in Section 2.4.3.3 and still used in the unsteady case are valid for steady flows only. Let us check whether it is still the case for unsteady flows.

Wave drag criterion

The wave drag criterion compares the local velocity in the direction of the pressure gradient to the speed of sound. In unsteady flows, we can encounter a moving shock wave. In that case, the criterion should be applied in the shock frame of reference. Indeed, in the inertial reference frame, a shock wave traveling in the opposite direction of the fluid motion can exist even if the velocity of the fluid is less than the speed of sound. Figure I.14 shows an example of a shock wave moving upstream in a tube. The relative velocity, in the shock frame of reference, is: q r = q -(-q shock ) = q + q shock can be larger than a even though q < a. This shock wave might therefore be undetected: the criterion used to define S w lacks robustness for unsteady flows.

Viscous drag criteria

The viscous drag criteria rely on viscous stress and turbulence level considerations. The velocity is not involved in the definitions, so that the criteria do not depend on the reference frame chosen. They are therefore valid for unsteady flows as such.

Physical background for the definition of the unsteady induced drag

It must be understood that the breakdown method, both thermodynamically and in terms of location, remains an artificial decomposition. All phenomena are actually linked and intricate in the fluid. If the wave and viscous components can be considered correctly defined, the definition of the induced component by default is a main flaw. In the unsteady case in particular, it must enclose several types of sources of drag, due to the unsteadiness.

This force remains however strongly associated with variations of lift, so that we chose to keep the notion of lift induced drag in the name of this drag component. It includes the force caused by reversible phenomena at first order, such as vortices and acoustics, as well as second order irreversible phenomena, such as the propagation of viscous interactions outside the viscous volume or variations of momentum upstream of the shock wave. These second order irreversible contributions are global in the fluid domain. This is why we chose to keep them in this induced drag component.

Comparison with Gariépy's formulation

The only other unsteady breakdown method found in the literature is Gariépy's formulation [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF]. The authors introduce a new definition of the axial velocity defect [START_REF] Gariépy | A New Axial Velocity Defect Formulation for a Far-Field Drag Decomposition Method[END_REF] within the steady decomposition and suggest to include the time derivatives into an unsteady drag coefficient. It is detailed in the Chapter Presentation of the main existing methods of numerical drag prediction. The comparison with our formulation requires to write both of them in the same reference frame. Gariépy had indeed chosen to derive the equations in the mobile reference frame. His formulation, translated in the inertial reference frame, gives, with similar notations:

D G w = Sw -ρ (u -u * rev ) (q • n) dS (I.4.10) D G v = Sv -ρ (u -u * rev ) q + τ x • n dS (I.4.11) D G i = Se -ρ (u rev -u ∞ ) q -(p -p ∞ ) i • n dS (I.4.12) D G uns = Se -ρ (u * rev -u rev ) (q • n) dS - V ∂ρ (u -u ∞ ) ∂t dV (I.4.13) D G sp = Vc ∪V wd ∇ • -ρ (u -u * rev ) q + τ x dV (I.4.14) D m = Sa -ρ (u -u ∞ ) (q • n) dS (I.4.15) D G ff = D G w + D G v + D G i + D G uns + D G sp + D m (I.4.16)
Here is a list of differences between the two formulations:

• Use of an unsteady drag component: unsteady terms are gathered in D G uns rather than distributed in the wave and viscous drag components. As a result, the synchronization with the force exerted on the skin is not ensured, and the drag coefficients are likely to be dependent on the extension of the integration surfaces.

• Theoretical expression for the spurious drag: the authors choose an analytic definition of the spurious drag, based on ∇ • f irr being allegedly zero. However, we have seen that it was not the case. This definition of the spurious drag may therefore lead to abnormally large values of artificial drag.

• No accounting for the wake of the shock in the wave drag component: we have seen in the derivation that the assumption of irreversible flow was valid only on a wake plane far from the sources. Applying the assumption on a surface closer to the shock wave requires to take the wake of the shock into account. Gariépy does not take it into account, so that his definition of wave drag lacks time history properties.

• Use of another definition of the axial velocity defect: an alternative expression first introduced by Méheut [START_REF] Méheut | Drag-Breakdown Methods from Wake Measurements[END_REF] is used here. Its robustness and accuracy needs validation on steady cases. A comparative study of both expressions is carried out in Chapter II.

Both unsteady formulations will be compared on simple unsteady test cases in Chapters III, and IV.

Chapter Summary

The obstacles of the generalization to unsteady flows have been identified. The steady formulation chosen has been carefully redemonstrated in order to locate all the assumptions made. The formulation has then been generalized to unsteady flows. The result is a breakdown method which includes unsteady contribution into each phenomenological drag component. Finally, flaws in the robustness and physical background have been discussed. Axes to improve these two points will now be explored.

Chapter II

Study of Improvement Axes for the Robustness and the Physical Background

This chapter aims at improving the robustness and the physical background of the formulation developed in the previous chapter. It has indeed been noted that the irreversible axial velocity could be undefined in certain regions of the flow, inducing loss of robustness. The criterion used for the definition of the wave volume can also lack accuracy with a moving shock wave. Robustness is however of crucial importance in a generalization to unsteady flows. Complex unsteady flows around innovative designs or at the boundary of the flight envelop are indeed the targeted applications. Another flaw of the unsteady formulation, which comes from the steady method itself, is the lack of physical background for the induced drag, which is derived by default, as the complementary part of drag once the other components defined.

The robustness is first looked into, with both issues that are the definition of the irreversible axial velocity and the criterion used for the wave drag integration. A possible improvement found in the literature for the first is the use of an alternative expression for the irreversible axial velocity proposed by Méheut. The derivation of the alternative expression is recalled in order to analyze its validity in the frame of the drag breakdown method. The theoretical accuracy is then studied, especially when the integration surfaces are moved. Finally, the formulation using this alternative expression is compared to the classical one on a whole set of steady test cases. The aim here is to assess the accuracy of the drag prediction on simple and well known steady test cases, and to infer the results for unsteady flows.

The second cause of robustness loss is then tackled. An unsteady version of the wave criterion found in the literature to improve the robustness of the computation of the wave surface is looked into, implemented and tested on unsteady test cases, with and without a moving shock. The effect of a filter, which is mandatory with the unsteady wave criterion, is also investigated.

The last axis of improvement concerns the physical background. The source terms of the unsteady induced drag are analyzed, the aim being to identify contributions of other sources of drag. These contributions may come from either the irreversible sources, wave or viscous, whose second order effects may have propagated outside the streamtubes, or new sources existing only in unsteady flows, such as acoustics. The share of each contribution also needs to be quantified.

Finally, a definitive formulation is proposed, along with good practice recommendations drawn from this study.

Study of an alternative expression for the irreversible axial velocity to improve the robustness

Méheut [START_REF] Méheut | Drag-Breakdown Methods from Wake Measurements[END_REF] introduced another expression for the irreversible axial velocity, since the non definition of u irr in certain regions of the flow may cause robustness trouble. The applications targeted in his paper were only drag prediction from experimental data, with measurements in a wake plane about one chord downstream of the model. The use of this expression closer to the body in numerical simulations could however imply an accuracy loss, so that a careful comparative study is required.

Derivation of the expression developed by Méheut

The definition of the new u irr is obtained by tackling the problem the other way around: Méheut assumes that the flow is reversible, which means isentropic and isenthalpic according to Van der Vooren's theory, in order to obtain a reversible axial velocity u rev as a function of p, v and w. He can therefore assume that, on a plane S d (as depicted in Figure II.1), there has been neither entropy nor stagnation enthalpy production:

• ∆s = 0 the entropy is equal to its reference state • ∆H = 0 the stagnation enthalpy is equal to its reference state These assumptions are actually valid on any plane downstream of a body. There is no need to use a wake plane far from the sources as in Van der Vooren's theory.

The general expression for the axial velocity as a function of the thermodynamic variables obtained in Equation I.2.26 is valid for any flow:

u = u ∞ 1 + 2∆H u 2 ∞ - 2 (γ -1) M 2 ∞   e γ-1 γ ∆s r p p ∞ γ-1 γ -1   - v 2 + w 2 u 2 ∞ (II.1.1)
Simply applying the assumptions of reversibility gives the expression of a reversible axial velocity u rev :

u rev = u ∞ 1 - 2 (γ -1) M 2 ∞   p p ∞ γ-1 γ -1   - v 2 + w 2 u 2 ∞ (II.1.2)
This velocity depends only on the pressure variations probed on the wake plane, and on the cross velocity components. Now we can define the alternative irreversible axial velocity as:

u † irr = u -u rev + u ∞ (II.1.3)
The reversible and irreversible parts of f are then defined using this alternative expression u † irr , or equivalently the reversible axial velocity defined:

f † rev = -ρ(u -u † irr )q -(p -p ∞ ) i = -ρ (u rev -u ∞ ) q -(p -p ∞ ) i (II.1.4) f † irr = -ρ(u † irr -u ∞ )q + τ x = -ρ (u -u rev ) q + τ x (II.1.5)
The rest of the breakdown method is identical. The resulting induced and profile drag coefficients are denoted D † i , D † w and D † v :

D † i = S d f † rev • n dS (II.1.6) D † w = S wd f † irr • n dS (II.1.7) D † v = S vd f † irr • n dS (II.1.8)

Domain of definition of the reversible axial velocity

The motivation behind the use of an alternative expression of the axial velocity was the question of the domain of definition. Let us now check whether this reversible axial velocity is indeed defined everywhere in the fluid, i.e. whether the term under the square root is always positive:

Γ rev = 1 - 2 (γ -1) M 2 ∞   p p ∞ γ-1 γ -1   - v 2 + w 2 u 2 ∞ ≥ 0 (II.1.9)
Since the flow is now assumed isentropic and isenthalpic, we know that the stagnation pressure is constant:

p i = p i∞ (II.1.10)
Using the isentropic laws, we can therefore write:

p 1 + γ -1 2 M 2 γ γ-1 = p ∞ 1 + γ -1 2 M 2 ∞ γ γ-1 (II.1.11) So that: p p ∞ γ-1 γ = 1 + γ -1 2 M 2 ∞ 1 + γ -1 2 M 2
(II.1.12)

Adding -1:

p p ∞ γ-1 γ -1 = γ -1 2 1 + γ -1 2 M 2 M 2 ∞ -M 2 (II.1.13)
Multiplying by -

2 (γ-1)M 2 ∞ therefore yields: - 2 (γ -1) M 2 ∞   p p ∞ γ-1 γ -1   = - 1 1 + γ -1 2 M 2 1 - M 2 M 2 ∞ (II.1.14)
Owing to the fact that

M 2 M 2 ∞ = q 2 u 2 ∞ a 2 ∞ a 2
, and using the isentropic law

a 2 ∞ = a 2 1 + γ-1 2 M 2 : - 2 (γ -1) M 2 ∞   p p ∞ γ-1 γ -1   = - 1 1 + γ -1 2 M 2 + q 2 u 2 ∞ (II.1.15) Adding 1 -v 2 +w 2 u 2 ∞ : 1 - 2 (γ -1) M 2 ∞   p p ∞ γ-1 γ -1   - v 2 + w 2 u 2 ∞ = 1 - 1 1 + γ -1 2 M 2 + q 2 u 2 ∞ - v 2 + w 2 u 2 ∞ (II.1.16) Factoring 1 1+ γ-1
2 M 2 finally yields:

Γ rev = 1 1 + γ -1 2 M 2 γ -1 2 M 2 + γ -1 2 M 2 u 2 u 2 ∞ + u 2 u 2 ∞ (II.1.17)
which is always positive.

We have therefore proved that the reversible axial velocity is defined everywhere.

Study of its theoretical validity

Now remember that in order to obtain the final formulation, the integration surfaces are moved upstream closer to the source of drag. In this step, it is assumed that the flow is isentropic and isenthalpic downstream of the sources. This step is mandatory to justify why we can integrate on a small surface around the shock wave, but also explains why the choice of the integration surfaces, as long as they enclose the sources, has a small impact on the drag coefficients. This point is of course very important since the method must give the same results with any position of the wake integration surface.

Here though, the alternative irreversible axial velocity does not depend only on ∆s and ∆H anymore. We can indeed write:

u † irr = u ∆s r , ∆H u 2 ∞ , p p ∞ , v 2 + w 2 u 2 ∞ -u rev 0, 0, p p ∞ , v 2 + w 2 u 2 ∞ + u ∞ (II.1.18)
Now let us expand at the first order this expression of u † irr for small variations of normalized entropy ∆s * = ∆s r , stagnation enthalpy

∆H * = ∆H u 2 ∞
, pressure p * = p p∞ , and cross velocity q * vw = v 2 +w 2 u 2 ∞ . The perturbations are assumed to be of the same order of magnitude. Using Taylor's theorem, we get:

u † irr ∼ u(0, 0, 1, 0) + ∂u ∂∆s * ∆s * + ∂u ∂∆H * ∆H * + ∂u ∂p * p * + ∂u ∂q * vw q * vw -u rev (0, 0, 1, 0) - ∂u rev ∂p * p * - ∂u rev ∂q * vw q * vw + u ∞ + O(∆s * ) + O(∆H * ) + O(p * ) + O(q * vw ) (II.1.19)
We have u(0, 0, 1, 0) = u rev (0, 0, 1, 0) = u ∞ . Furthermore, we can assume that at first order ∂u ∂p

≃ ∂u rev ∂p and ∂u ∂(v 2 + w 2 ) ≃ ∂u rev ∂(v 2 + w 2 )
. We thus obtain a first order approximation of this alternative expression of the irreversible axial velocity:

u † irr ∼ u ∞ + ∂u ∂∆s ∆s + ∂u ∂∆H ∆H ∼ u irr (II.1.20)
As a conclusion, both expressions for the irreversible axial velocity are equivalent at first order . However, it is no more the case at the second order, since the second order crossed terms will no more cancel out, and the assumption that the derivatives of u and u rev are close will no more hold. Both expressions should accordingly have similar behaviors if the variations in the fluid remain small, but it can be expected divergences when they become larger. We will verify this allegation shortly hereafter.

Analysis of the variant suggested by Gariépy

Gariépy and Trépanier assert in their article [START_REF] Gariépy | A New Axial Velocity Defect Formulation for a Far-Field Drag Decomposition Method[END_REF] that ∆H should be taken into account in u rev , which is equivalent at first order to withdrawing it from u † irr . It may be of limited consequences for steady test cases. It might however become very large for unsteady cases, since the stagnation enthalpy varies strongly versus time. We will denote it u * rev and add a * superscript to the corresponding variables.

u * rev = u ∞ 1 + 2∆H u 2 ∞ - 2 (γ -1) M 2 ∞   p p ∞ γ-1 γ -1   - v 2 + w 2 u 2 ∞ (II.1.21) u * irr = u -u * rev + u ∞ (II.1.22)
They try to show that stagnation enthalpy variations are solely caused by reversible effects, expressing the material derivative of H as a function of the time derivative of pressure: dH dt = 1 ρ ∂p ∂t . Since the variations of pressure on a wake plane are related to reversible processes, they infer that the variations of stagnation enthalpy probed on the wake plane are solely due to reversible effects.

If this statement might be true, there is no theoretical motivation for the variations of stagnation enthalpy in the whole fluid domain, as it is computed by the formulation, to be due to reversible effects only. The time derivative of pressure should in particular not be linked to reversible processes only. This is why we chose to keep this term in our evaluation of u † irr . The variation of stagnation enthalpy due to the unsteadiness of a shock should for example appear inside the wave drag.

All three formulations will now be compared for several steady applications.

Comparison of the three expressions on several steady test cases

The three expressions derived for the irreversible axial velocity defect, u irr , u † irr and u * irr , are now compared on a whole set of steady test cases, with a growing complexity. The aim is to assess the robustness and accuracy of the alternative expressions, especially when moving the downstream extension of the integration surfaces. For all test cases, second order central discretization with Jameson artificial viscosity and the turbulence model of Spalart-Allmaras, when the flow is viscous and turbulent, are used.

Airfoil in a transonic inviscid flow: assessment of CD w

The first test case consists in isolating the wave drag component. The inviscid transonic flow around a 2-D airfoil is therefore studied. The wave drag component can therefore be compared to the near-field drag evaluation.

Description of the test case

The test case is a NACA0012 profile. The mesh is a very regular 1,048,276 elements grid adapted for Euler computations, courtesy from Vassberg and Jameson [START_REF] Vassberg | In Pursuit of Grid Convergence for Two-Dimensional Euler Solutions[END_REF], which can be seen in Figure II.2(a). The mesh extends to 150 chord lengths around the airfoil. The aerodynamic conditions are: M ∞ = 0.8, and α = 0 deg.

Convergence study

The convergence of the CFD computation can be analyzed in

Analysis of the flow field resulting from the simulation

A visualization of the flow field can be found in Figure II.3(a). We can see that a strong shock wave is present on both the upper and the lower sides, at about mid-chord.

Application of the drag extraction method

The quantity u irr is defined everywhere in this case, so that there is no issue with the integration surface. The wave surface S w computed by ffd72, using the physical criteria with three additional cell layers, encloses the shock wave, as we can see in

Analysis of the drag breakdown results

The previous wave surface is then progressively extended downstream by varying the number of cell layers. This number is also reduced to zero in order to analyze the behavior of all three expressions for the drag coefficient. The resulting curves are visible in Figure II. [START_REF] Basset | Treatise on hydrodynamics[END_REF]. We can see that the classical expression CD w , which uses u irr , rapidly converges to a value very close to the near-field value (83.59 d.c.).

On the contrary, CD † w and CD * w , which use u † irr and u * irr , converge towards a different value (about 77 d.c.). Moreover, the behavior of CD * w very close to the source of drag is very different from the other two expressions, with strong discrepancies. The discrepancy with CD † w can only be explained by numerical enthalpy variations very close to the shock wave which are not taken into account in this expression. As noted in the theoretical study, the expressions are equivalent at first order only, so that discrepancies close to the source of drag could be expected.

Airfoil in a subsonic viscous flow: assessment of CD v

The second drag component looked into is the viscous drag. To do so, a 2-D airfoil in a viscous subsonic flow is studied.

Convergence study

The convergence curves for the residuals converge slowly but decrease of about 6 orders, as we can see in Figure II.5(b). The near-field drag coefficient is also converged to the thousandth of drag count after 10,000 iterations. A grid convergence study has also been carried out in this case and can be found in Appendix B, Section B.2.

Analysis of the flow field resulting from the simulation

A visualization of the flow field can be found in Figure II.8. The boundary layer developing around the airfoil and the wake are clearly visible.

Application of the drag extraction method

The viscous surface is computed using the physical criteria described earlier. We can see in They are easily enclosed within the viscous surface, without needing to deform it, so that the drag results should not be impacted.

Analysis of the drag breakdown results

The downstream extension of the viscous volume is varied from 0 to 5 chords downstream of the trailing edge in Figure II.8. We can see that all three formulations converge towards the same value, which is consistent with the near-field value CD nf = 90.64 d.c. The two variants of the reversible expression CD † v and CD * v give almost superimposed curves. There is actually a very small discrepancy between both curves when the convergence is reached. We can however see that the classical expression CD v converges more rapidly. We can conclude that the use of u irr allows to bring the integration surface much closer to the trailing edge than the use of u † irr or u * irr . This confirms the equivalence at first order only highlighted in the theoretical study. The third test case is a 3-D wing in a subsonic inviscid flow. The aim is here to isolate the third drag component: the induced drag.

Description of the test case

The wing is an elliptic wing based on the NACA0012 profile, with an aspect ratio λ = 8. The grid has been adapted from the 2-D grids of Vassberg and Jameson by Destarac and are available online [92]. The mesh is about 8 million elements, with an extension of 150 chord lengths around the wing, and is represented in Figure II.9(a). The aerodynamic conditions are a Mach number M ∞ = 0.2, and α = 6 deg.

Convergence study

The residuals of density and stagnation energy decrease of more than 8 orders after 4,000 iterations, as shown in Figure II.9(b). The near field drag coefficient is also converged to the thousandth of drag count, so that the convergence is considered satisfactory. A grid convergence study is also provided in Appendix B, Section B.3.

Application of the drag extraction method

In this case, there is no zone where u irr is undefined. The induced volume consists in considering the whole fluid domain cut at a given downstream extension. The discrepancy is even stronger with the alternative expressions CD † i and CD * i , which give very similar results. Again, the equivalence at first order only can explain these discrepancies between the various far-field formulations in this zone close to the trailing edge.

Wing in a transonic viscous flow: assessment of all three drag components

This fourth test case is a synthesis case, which allows to evaluate the behavior of all three drag components when all phenomena are present.

Description of the test case

The wing is a rectangular NACA0012-based wing with an aspect ratio of 8. The mesh is around 1 million elements and is shown in

Analysis of the drag breakdown results

The downstream extension of the surfaces was set to vary in Figures II. [START_REF] Campobasso | Turbulent Navier-Stokes Analysis of an Oscillating Wing in a Power-Extraction Regime using the Shear Stress Transport Turbulence Model[END_REF] The results of the synthesis test case confirm the observations made in the previous test cases: the classical expression CD using u irr as defined by Van der Vooren is quite reliable, depending very little on the distance to the body, even very close to the sources of drag. The expressions CD † and CD * , using u † irr and u * irr , give less satisfactory results: they are strongly dependent on the integration volumes, especially close to the trailing edge. These expressions are therefore valid for experimental data measured far from the body but less reliable for numerical results. These observations are consistent with the theoretical remarks made in Section 1.3. Another comment is that there is very little difference between CD † and CD * in the steady case, since the variation of stagnation enthalpy is then very small. It will no more be the case for unsteady flows.

The conclusion of this first study is that the integration volumes and surfaces must be chosen very carefully when using u † irr or u * irr . Although these conclusions are valid only for steady flows, we draw some recommendations for unsteady flows: we chose to keep u irr in our unsteady formulation, and the downstream extension will also be chosen sufficiently far from the body when using Gariépy's formulation.

As noted in Chapter I, Section 4.1.2, the physical criterion used for the determination of the surface enclosing the shock wave may lose robustness in the case of a mobile shock wave or a mobile body. A literature review [START_REF] Ziniu | Review of Shock Wave Detection Method in CFD Post-Processing[END_REF] indicates that there exist filters as well as an unsteady version of the wave criterion in order to account for transient phenomena. They have been implemented in the ffd72 code, and assessed on two unsteady cases, with and without shock waves.

Expression of the unsteady criterion

The criterion for a mobile shock wave has first been introduced by Lovely and Haimes [START_REF] Lovely | Shock Detection from Computational Fluid Dynamics[END_REF], and has been used by Gariépy [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF]. The idea consists in adding an unsteady term to the steady physical criteria as expressed in Equation I.2.67. Lovely and Haimes note that, in the unsteady case, a transient shock wave corresponds to a material derivative of pressure such that, in the inertial reference frame:

dp dt ≥ a ∇p (II.2.1)
The definition of the material derivative yields:

dp dt = ∂p ∂t + q • ∇p (II.2.2)
The unsteady version of the criterion can therefore be expressed as:

q • ∇p + ∂p ∂t ≥ c shock a ∇p ⇒ cell ∈ S uns w (II.2.3)
Removing the unsteady term ∂p ∂t allows retrieving the steady criterion. Since the computation of a time derivative requires computational effort and storage, Lovely and Haimes then transform the time derivative of pressure into a divergence term. Gariépy also chose to use this expression with a divergence. However, since we are already computing time derivatives for the formulation itself, we choose to keep the original expression of Equation II.2.3 in order to avoid errors in the computation of the gradient.

Evaluation on an unsteady subsonic test case

The unsteady wave criterion is first tested on an unsteady case which is free of shock. The aim is here to make sure that the criterion does not wrongly detect cells as belonging to a shock surface. The test case is presented in details in Chapter III, Section 1. It consists of a NACA0012 profile placed in a flow field at low Mach number M ∞ = 0.2 and at very high angle of attack α = 20 deg. A natural instability develops as the flow separates, forming two vortices which are periodically emitted and advected along the wake. There is no shock wave, so that the criterion should stay inactive everywhere in the fluid. Figure II.16 shows however that many cells activate the unsteady wave criterion, whereas the steady criterion is never activated. This is quite problematic for a low Mach number case as this one.

Evaluation on an unsteady transonic test case

The unsteady wave criterion is then tested on a buffet case. The test case is presented in details in Chapter III, Section 2. The supercritical profile OAT15A is placed in a flow field at Mach number M ∞ = 0.73 with an angle of attack α = 4.5 deg. Due to the interaction between the shock wave and the boundary layer, a natural instability develops: the shock wave starts oscillating on the upper side of the airfoil.

During the motion of the shock wave in the upstream direction, the situation described in Chapter I, Section 4.1.2 may occur: the shock wave may go undetected if the velocity of the fluid becomes less than the speed of sound due to the upstream motion, the relative velocity being still as strong. It is therefore a good test case.

The computation of the integration volumes and surfaces by the ffd72 code using both steady and unsteady wave criteria however gives very different results. Without filtering, the unsteady criterion tags many cells which do not belong to the wave surface as such, as depicted in Figure II.17 In order to solve this problem, Lovely and Haimes propose several filtering techniques. The conclusion of their study is that a simple filter on the pressure gradient is the best solution, since shock waves correspond to zones with high pressure gradients. The main issue is however to find the correct threshold value. They devise a sophisticated method using the derivative of the curve representing the number of cells for each value of the pressure gradient.

Filtering

We decided to initially test the filter on the pressure gradient using an arbitrary threshold value:

∇p ≤ ∇p threshold ⇒ cell / ∈ S uns,filtered w (II.2.4)
For the subsonic case, since we want all tagged cells to be cut out of the wave surface, the threshold value simply has to be chosen large enough. The transonic case is a more challenging issue. Several threshold values have been manually tested. The best value was ∇p threshold = 10 6 Pa.m -1 . The results evolving with respect to time are presented in Appendix C, Figure C. [START_REF] Huang | Vortex Shedding and Shear-Layer Instability of Wing at Low Reynolds Numbers[END_REF]. We can see that cells around the leading edge are still tagged despite the filtering. A larger threshold value actually removes cells from the surface around the shock wave rather than at the leading edge. Despite this observation, the surface using the filtered unsteady criterion S uns,filtered w is very close, even slightly larger, to S p w as the shock wave is moving upstream, i.e. when the criterion is expected to be better suited, as in Figure II.18(b). On the contrary, as the shock wave moves downstream, there are more discrepancies between both surfaces, and S uns,filtered w becomes even much smaller, as at time step 400 in Figure II.18(a). The foot of the shock wave is also excluded from the wave surface at most time steps.

Conclusions on the validity of the unsteady wave criterion

As a conclusion, although the filtered unsteady wave criterion gives acceptable results for the buffet case when the shock wave is moving in the upstream direction, the other results are not satisfactory. In particular, it is quite problematic that wave surfaces should be identified in a subsonic flow. In practice, we would advise not to use this unsteady criterion, except if the steady version fails at detecting a moving shock wave, which has not yet occurred in this work.

On the theoretical point of view, the fact that the criterion without filter selects almost all cells in the fluid domain is actually quite astonishing. It is in fact as if only the filtering on the pressure gradient allowed defining the wave surface. Questions about the theoretical validity of the unsteady criterion are therefore raised.

Study of the physical interpretation of the volume term in the unsteady induced drag component

We have seen in Chapter I, Section 4.2 that the unsteady induced drag component, defined by default as the remaining component once the wave and viscous components defined, could enclose phenomena other than lift induced ones. A careful physical study is therefore required in order to be able to correctly interpret the behavior of this drag component.

Link between surface and volume terms

The aim here is to better understand the role played by the volume term in the unsteady induced drag component. Let us imagine that the induced drag is defined only by the surface term in Equation I. We can similarly as before reduce the integration surface to the wake plane S d . The induced drag can therefore be expressed as:

D i = S d -ρ (u -u irr ) q -(p -p ∞ ) i • n dS (II.3.1)
We want the induced drag computed on the wake plane S ′ d to be equal to this expression. A correction term denoted D corr needs to be added:

D i = S ′ d -ρ (u -u irr ) q -(p -p ∞ ) i • n dS + D corr (II.3.2)
We can now compute this corrective term as the difference between both surface terms:

D corr = S d -ρ (u -u irr ) q -(p -p ∞ ) i • n dS - S ′ d -ρ (u -u irr ) q -(p -p ∞ ) i • n dS (II.3.3)
Applying the divergence theorem, and using the fact that the contributions on the lateral surfaces are once again zero, we can express this corrective term as a volume integral:

D corr = - V d ∇ • -ρ (u -u irr ) q -(p -p ∞ ) i dV (II.3.4)
Expanding the integrand I:

I = -∇ • ρuq - ∂p ∂x + ρq • ∇u irr + u irr ∇ • ρq (II.3.5)
The x-component of the Navier-Stokes equations, neglecting the viscous term, implies that:

∂ρu ∂t + ∇ • ρuq = - ∂p ∂x (II.3.6)
On the other hand, the definition of the material derivative of u irr gives:

du irr dt = ∂u irr ∂t + q • ∇u irr (II.3.7)
Plugging Equations II.3.6, and II.3.7 into Equation II.3.5, we get:

I = ∂ρu ∂t + ρ du irr dt -ρ ∂u irr ∂t + u irr ∇ • ρq (II.3.8)
The second term has already been computed in Chapter I, Section 3.2.1:

ρ du irr dt = 1 u irr ∂p ∂t (II.3.9)
The third term can be written:

ρ ∂u irr ∂t = ∂ρu irr ∂t -u irr ∂ρ ∂t (II.3.10)
The continuity equation ∂ρ ∂t + ∇ • ρq = 0 allows to cancel out the last term with the last term from Equation II.3.8. We finally get:

D corr = - V d ∂ρ(u -u irr ) ∂t + 1 u irr ∂p ∂t dV (II.3.11)
This expression corresponds exactly to the volume term in the definition of the unsteady induced drag. We can therefore deduce that this term allows accounting for the non dependence on the integration surface, as well as for the synchronization in time.

Now we can analyze separately the surface and the volume terms. The surface term indeed accounts for the reversible effect measured on the wake plane, whereas the volume term allows to remove the dependency on the position of the wake surface. This volume term is in practice computed over the whole complementary volume V c . It may therefore enclose other sources of drag. We will now focus on analyzing its physical content.

Acoustic effects

One of the first ideas to better understand what the volume term may enclose has been to plot the field distribution of its source term. It has been done with all unsteady test cases.

The first unsteady test case is the vortex shedding case presented in the previous section and in details in Chapter III, Section 1. The source term is plotted in the flow field in Figure II.20. The color scale was chosen so as to highlight the phenomena occurring in the far field. Acoustic waves can be recognized in red and blue, mainly centered around the trailing edge of the airfoil. They also propagate around the airfoil as time goes by. The shape of the circles is slightly deformed by the reference velocity, which is relatively small (M ∞ = 0.2). The distance between two waves in the direction of the upcoming flow has been checked to be equal to the speed of sound corrected by the reference velocity, multiplied by the time separating two successive waves. The wake of vortices also seems to behave like a secondary acoustic source. A positive and a negative wave are emitted respectively over and under the airfoil, so that the acoustic field is very similar to the field emitted by an acoustic dipole. The same kind of flow visualization has been applied to the second unsteady test case, the buffet case (details can be found in Chapter III, Section 2). The waves propagating around the airfoil are also quite specific of acoustic phenomena, as depicted in Figure II.21. The circles are here much more deformed than in the vortex shedding case, due to the fact that the reference velocity is larger (M ∞ = 0.73). Again, the distance between two waves corresponds to the distance traveled by a perturbation propagating at the speed of sound in the field moving at the reference velocity during the time separating two waves. The alternating positive and negative waves are also very similar to the acoustic field around a dipole source.

The source term of the volume integral, ∂ρ(uu irr ) ∂t + 1 u irr ∂p ∂t , therefore seems to be a very good indicator of acoustic effects. Further investigation has revealed that both terms were equally strong and both showed the same kind of pattern in the flow field. Only the term -∂ρu irr ∂t is negligible

∂ρ(u -u irr ) ∂t + 1 u irr ∂p ∂t = ∇ • ρ (u irr -u) q + ∂p ∂x = ∂ρ(u -u ∞ ) ∂t + ∇ • ρ (u irr -u ∞ ) q (II.3.12)
The behavior of each term in every formulation is however always the same: strong variations in time and correlation between the two of them.

Although the presence of acoustic effects has been proven, it may be interesting to know whether acoustics can be a source of drag and, if so, what is the order of magnitude of such effects, in particular whether it is negligible with respect to the other phenomena or not.

Breakdown of the unsteady induced drag component

We have seen in the previous section that acoustics could be a source of unsteady drag worth accounting for, and that its contribution was enclosed within the volume term of the unsteady induced drag component. Other sources of drag can be enclosed within this volume term: second order irreversible phenomena can appear, such as the propagation of viscous interactions outside the viscous volume, or variations of momentum upstream of a shock wave. These second order irreversible contributions are global in the fluid domain.

The idea in this section is to achieve a breakdown of the volume term, in order to part the unsteady induced drag component into a contribution related to acoustics, and a component related to lift variations, and global contributions in the fluid domain. The idea consists in retrieving the Riemann invariant for a potential flow R + = u + 2a γ-1 . We can write the volume term I as:

I = ∂ρ(u -u irr ) ∂t + 1 u irr ∂p ∂t = ρ ∂(u -u irr ) ∂t + (u -u irr ) ∂ρ ∂t + 1 u irr ∂p ∂t (II.3.13)
Assuming an isentropic flow, which is a valid assumption in the considered volume, we can write: ∂p ∂t = a 2 ∂ρ ∂t . Therefore:

I = ρ ∂u ∂t -ρ ∂u irr ∂t + u -u irr + a 2 u irr ∂ρ ∂t (II.3.14)
Unsteady Bernoulli theorem says:

γ γ -1 p ρ + q 2 2 + ∂ϕ ∂t = a 2 ∞ γ -1 + u 2 ∞ 2 (II.3.15)
Deriving with respect to time, assuming a 1-D flow and applying the isentropic law for the pressure, yields:

a 2 ρ ∂ρ ∂t + u ∂u ∂t + ∂ 2 ϕ ∂t 2 = 0 (II.3.16)
So that the derivative of density can be written as:

∂ρ ∂t = - ρ a 2 u ∂u ∂t + ∂ 2 ϕ ∂t 2 (II.3.17)
The Bernoulli theorem as written in Equation II.3.15 also yields:

∆H = - ∂ϕ ∂t (II.3.18)
So that the time derivative of the irreversible axial velocity can be written:

∂u irr ∂t = - 1 u irr ∂ 2 ϕ ∂t 2 (II.3.19)
As a result, Equation II.3.14 yields:

I = ρ(u irr -u) a 2 ∂ 2 ϕ ∂t 2 + u + a 2 u irr ∂u ∂t (II.3.20)
The time derivative of the Riemann invariant is expressed as:

∂R + ∂t = ∂u ∂t + 2 γ -1 ∂a ∂t (II.3.21)
The time derivative of the Bernoulli theorem, written with the speed of sound, yields:

2 γ -1 ∂a ∂t = - u a ∂u ∂t - 1 a ∂ 2 ϕ ∂t 2 (II.3.22)
The time derivative of the Riemann invariant can therefore be expressed as:

∂R + ∂t = - 1 a ∂ 2 ϕ ∂t 2 + (u + a) ∂u ∂t (II.3.23)
Finally, we get for the volume term I:

I = ρ(u irr -u) a - ∂R + ∂t + 1 a - 1 u irr ∂u ∂t (II.3.24)
We can now isolate an acoustic contribution, expressed as a function of the Riemann invariant, in the volume term of the unsteady induced drag component. It may also enclose other reversible propagation effects. The resulting drag coefficient is therefore called propagation and acoustics drag component D pa :

D pa = Vc ρ(u irr -u) a ∂R + ∂t dV (II.3.25)
The residual term is then added to the surface term of the unsteady induced drag, so as to define an induced drag component D i :

D i = Se -ρ (u -u irr ) q -(p -p ∞ ) i • n dS - Vc ρ (u -u irr ) 1 a - 1 u irr ∂u ∂t dV (II.3.26)
Although the theoretical content of this term is difficult to analyze, it exhibits interesting features, such as positiveness, good orders of magnitude, and independence on the integration volume chosen, if sufficiently far from the body. Remember that the whole volume term was the one responsible for the independence on the integration surface. We have split this term in two, and kept only the residual term. It appears however that this residual can ensure the independence once sufficiently far from the body. It should be equivalent to the whole volume term in this zone.

Taking all these reflections and observations into account, we have devised a robust drag breakdown method which is summarized here.

Final formulation with five components

Introducing the breakdown of the unsteady induced drag derived in Section 3.3 into the unsteady formulation developed in Chapter I yields: This formulation has been implemented in the ffd72 code. It has also been published in [88].

D w = Sw -ρ (u irr -u ∞ ) (q • n) dS - Vw ∂ρ (u -u ∞ ) ∂t dV - V wd ∂ρ (u -u irr ) ∂t + 1 u irr ∂p ∂t dV D v = Sv -ρ (u irr -u ∞ ) q + τ x • n dS - Vv ∂ρ (u -u ∞ ) ∂t dV D i = Se -ρ (u -u irr ) q -(p -p ∞ ) i • n dS - Vc ρ (u -u irr ) 1 a - 1 u irr ∂u ∂t dV D pa = Vc ρ(u irr -u) a ∂R + ∂t dV D m = Sa -ρ (u -u ∞ ) (q • n) dS D ff = D w + D v + D i + D pa + D m D sp = D nf -D ff (II.

Good practice recommendations

The good practice recommendations drawn from the studies carried out in this chapter are summarized in Table II Ideas to improve the robustness of the unsteady drag breakdown method have been explored. The contemplated solutions have proved disappointing, while the solutions used in the steady method actually raise no issue with all the test cases considered so far. The use of the alternative irreversible axial velocity as well as the unsteady wave criterion should therefore be used only when the classical tools fail. The knowledge of the physical background of unsteady induced drag has also been improved, highlighting in particular acoustic effects. A new breakdown of this drag component has been accordingly proposed. The resulting final formulation, as well as good the practice recommendations, can now be applied to unsteady test cases.

Chapter III

Assessment of the Wave, Viscous, and Acoustic Drag Components on Naturally Unsteady Cases

This chapter aims at assessing the wave, viscous, and propagation and acoustics drag components which are expected to appear in 2-D naturally unsteady cases. Since no other data is available for the breakdown of drag, only the total far-field drag can be validated against the near-field drag. The behavior of the components will therefore have to be analyzed in terms of order of magnitude and synchronization in time with the flow features. The decomposition will also be compared to Gariépy's.

The first test case is the shedding of vortices downstream of an airfoil at high angle of attack. The instability is natural, and the flow is inviscid, so that the only expected drag components are the viscous, and the propagation and acoustics drag coefficients. The conditions of the simulation are first introduced, the convergence is then evaluated, as well as the resulting unsteady flow field. The drag breakdown method can then be applied. The behavior of the resulting drag coefficients is then analyzed and compared to the ones obtained applying Gariépy's formulation, and steady methods.

The second test case is a buffeting instability on the upper side of a supercritical airfoil at transonic speed. A wave drag component will therefore add to the viscous and propagation and acoustics components. The flow simulation is carried out and converged, before being post-processed. The evolution in time of the drag coefficients resulting from the breakdown is then analyzed, in terms of orders and magnitude as well as synchronization in time with the flow features. Gariépy's formulation and steady methods are also compared. Conclusions regarding the validity of the far-field breakdown can finally be drawn.

Application to a vortex shedding case

The first unsteady case is the natural instability which develops downstream of an airfoil at very high angles of attack and low speed: similarly to the wake behind a bluff body, positive and negative free vorticity is alternatively shed in the wake. This test case is representative of the flow which can be experienced by the blades of a rotor, for which the local angle of attack can reach such high values. The drag components playing a role in this test case should be the viscous and the propagation and acoustics components.

Quick literature review

The idea for this simple test case has been found in a Ph.D. dissertation about unsteady optimization [START_REF] Rumpfkeil | Airfoil Optimization for Unsteady Flows with Application to High-Lift Noise Reduction[END_REF]. This test case has been chosen since it is one of the most simple: a motion less 2-D airfoil, and no shock waves. Only viscous and acoustic effects are present in this case. The features of this type of flow are: a strong separation with emission of positive and negative vortices, which are shed in the wake. In this Ph.D. dissertation, the flow is simulated by URANS simulations at Reynolds number 2 million.

Such a simple test case has not been often studied. An experimental study by Mesquita [START_REF] Mesquita | Experimental Analysis of Airfoil for High Angle of Attack[END_REF] show that an average drag of around 2,000 d.c. is experienced by the NACA0012 airfoil when the angle of attack reaches 20 deg at Reynolds number 1 million. Studies of the frequency of the vortex shedding experimentally measured for several Reynolds number and angles of attack have been carried out by Huang and Lin [START_REF] Huang | Vortex Shedding and Shear-Layer Instability of Wing at Low Reynolds Numbers[END_REF], and Lee and Huang [START_REF] Lee | Frequency Selection of Wake Flow behind a NACA 0012 Wing[END_REF], showing that the Strouhal number at 20 deg is around 0.2 for small Reynolds numbers (up to 30,000). No data for larger Reynolds numbers is however provided.

A numerical study of the vortex shedding downstream of a NACA0012 profile at a low Reynolds number (50,000), and for two smaller angles of attack (9.25 and 12 deg), with DNS simulations, has also been performed by Rodriguez et al. [START_REF] Rodríguez | Direct Numerical Simulation of a NACA0012 in Full Stall[END_REF], highlighting the flow structures appearing on the upper side of the airfoil.

Description of the test case

The considered airfoil is a NACA0012 profile. The aerodynamic conditions are a very high angle of attack α = 20 deg, a low Mach number M ∞ = 0.2, and a Reynolds number Re = 2 × 10 6 based on the free stream velocity and the chord. The turbulence model Menter kω SST was chosen for its better capabilities on detached flows (see Appendix A for details), and the numerical scheme is AUSM-P. The mesh is a 270,000 elements 2-D mesh (see Figure III.1(a)) whose extension is 20 chord lengths in every direction, and adapted for a wake with a high skew angle.

Convergence study

The time step is 5 × 10 -6 s, which corresponds to 3,200 steps by period. A steady computation has first been carried out, from which the unsteady simulation has been run. Several periods have been simulated before extracting the flow field. The convergence to a periodic state is ensured looking at the lift vs drag curve shown in The induced drag is rather small in this case, but not negligible. It was expected from the experience with steady flows to be almost zero for a 2-D case. We can try to interpret its evolution in time by comparing it with the evolution of the square of the lift coefficient, as has been done in Figure III.6(a). There is a strong correlation between both coefficients, similarly to a law of the form CD th i = CL 2 πλ which is used for steady elliptical loading in 3-D cases. We can however observe a small time delay between both coefficients, which may indicate a more complex law in this case. It could also be due to second order irreversible phenomena which could not be separated from the rest. The conclusion is therefore that the variations in time of the lift coefficient induce small variations of the induced drag component as defined with this formulation, whereas it was zero for a steady case where the lift was constant.

The propagation and acoustics drag component is also small, but not negligible. This can be explained by the fact that the reference Mach number is small, so that the compressible effects are low. We have already seen that the acoustic propagation did exist in this case, as visualized in

Comparison with Gariépy's formulation

Gariépy's formulation, as adapted in Chapter I, Section 4.3, has been implemented and applied on this test case. The corresponding results are presented in Figure III.7. The far-field drag is here closer to the near-field drag since the spurious drag coefficient is included in its definition. The spurious drag component, on the other hand, is much stronger: at most 100 d.c. or 5% of the total drag. The viscous drag is smaller, around 1500 d.c., whereas the unsteady drag component is large and varies strongly. The induced drag is almost constant, with a small negative value which can be considered Concerning the synchronization in time, it was found during the study that the downstream extension of the integration volumes had a very strong impact on the phase of the drag coefficients. In Figures III.8(a), and III.8(c), we can see for example the evolution of the viscous and the unsteady drag coefficients respectively as the downstream extension of the volume is set to vary. The results given in Figure III.7 correspond to the largest downstream extension, for which the spurious drag was smallest. We can see that the phase as well as the amplitude vary strongly. The synchronization in time is therefore not ensured with this formulation. In Figures III.8(b), and III.8(d), the same study has been carried out with our definitions of the viscous and induced drag components. We can see that there is this time very little dependence on the volume, and that the phase, i.e. the synchronization in time, is in particular very stable.

Comparison between steady and time-averaged unsteady results

The idea is now to compute the averaged drag coefficients, as well as the standard deviations, and to compare the results with the steady formulation applied to the averaged flow field. The results are expected to be different. The aim is to prove that a steady evaluation cannot be considered reliable.

This comparison is presented in Table III.1. The total drag is quite similar for both cases. The viscous drag component is however quite different between both cases. The steady formulation predicts no induced drag while it is non negligible in the unsteady case. Propagation and acoustics drag is also present in the unsteady case. The unsteady analysis is also richer since it provides the time evolution of the coefficients as well as their maximal values. 2 Application to a buffet case simulated by a URANS method

The second unsteady test case consists in studying the drag coefficients of an airfoil subject to buffeting. Buffet is a natural instability which can appear on a supercritical airfoil at transonic regimes: the interaction between the shock wave and the boundary layer causes the development of an instability which triggers the motion of the shock wave along the upper side of the airfoil. The expected drag components are the wave, viscous, and propagation and acoustics coefficients.

Quick literature review

In transonic flow conditions, the shock wave/turbulent boundary layer interaction and flow separations induce flow instabilities called buffet. This unsteady behavior of the flow may provoke structure vibrations named buffeting, which can have an important influence on the aerodynamic behavior of the aircraft. Design standards for aeroplanes limits the intensity of the buffeting phenomenon, and so, the prediction of the onset of the buffet is a problem of outstanding importance, which cannot be solved easily in wind tunnel experiments because of the dynamic behavior of the structure of models. A good way to predict the onset of the buffet is the numerical approach. A numerical drag breakdown analysis could also provide an interesting analysis tool to better predict and control this instability.

Buffet has been extensively studied for the past 30 years. First experimentally, as in the work of Benoit and Legrain [START_REF] Benoit | Buffeting Prediction for Transport Aircraft Applications based on Unsteady Pressure Measurements[END_REF], or Caruana et al. [START_REF] Caruana | Buffet and Buffeting Active Control[END_REF], then numerically. The first numerical studies focused on the prediction of the apparition of the instability itself, using steady computations [START_REF] Chung | Prediction of Transonic Buffet Onset for an Airfoil with Shock Induced Separation Bubble using Steady Navier-Stokes Solver[END_REF]. The unsteady behavior was however not predicted with such methods. Unsteady simulations for 2-D airfoils have then been initiated [START_REF] Wang | A Computational Study of Unsteady Turbulent Buffet Aerodynamics[END_REF], encountering difficulties when using RANS models. The development of the instability was indeed ill-predicted compared to the experimental data, as well as the pressure coefficient curves. The simulations were particularly sensitive to the turbulence model chosen. Brunet [START_REF] Brunet | Computational Study of Buffet Phenomenon with Unsteady RANS Equations[END_REF] finally achieved a satisfactory unsteady simulation of a 2-D case using the k -ω turbulence model from Wilcox [START_REF] Wilcox | Reassessment of the Scale-Determining Equation for Advanced Turbulence Models[END_REF], although the instability appeared at α = 4.5 deg rather than 3.5 deg in the experiments. He also studied the influence of the numerical scheme, and obtained good comparisons with the experimental data by Jacquin et al. [START_REF] Jacquin | Experimental Study of Shock Oscillation over a Transonic Supercritical Profile[END_REF]. The influence of the test section geometry has also been assessed by Thiéry and Coustols [START_REF] Thiéry | URANS Computations of Shock Induced Oscillations over 2-D Rigid Airfoils: Influence of Test Section Geometry[END_REF], highlighting the major effects of upper and lower walls as well as side walls with corner flow separations on the dynamics of the shock motion.

Description of the test case

This test case consists of an OAT15A profile under buffet conditions from Brunet [START_REF] Brunet | Computational Study of Buffet Phenomenon with Unsteady RANS Equations[END_REF]. The chord is 0.23 m, and the mesh is 2-D with around 150,000 elements (see Figure III.9(a)), extending about 10 m upstream, downstream and laterally. The turbulence model is kω Wilcox and the numerical scheme is Roe, chosen for their favorable comparison with experimental results. The Mach number is M ∞ = 0.73, α = 4.5 deg and Re = 3 × 10 6 based on the free stream velocity and the chord.

Convergence study

The study was carried out over one period, with a time step of 1.7× 10 -5 s, which corresponds to 1,000 steps by period. The unsteady computation was converged over several periods in order to reach the full periodicity and avoid transient phenomena.

Analysis of the drag breakdown results

The drag extraction was carried out over one period. The resulting time evolutions are shown in Figure III.13. The near-field drag varies between 500 and 800 drag counts. The total far-field drag (pink) increases while the shock moves downstream and decreases while the shock moves upstream. It is in good agreement with the near-field drag. The spurious drag (orange), computed as the balance between near-field and far-field, is indeed at most 20 drag counts or 2% of the total drag. As far as the order of magnitude of the drag components are concerned, we can see that the drag is mainly due to viscous, wave, and acoustic effects, as was expected. The viscous drag, in green, varies strongly, which is in agreement with the strong variations in the flow from attached to separation states. The wave drag, in red, varies much less. Although the shock waves is in motion above the upper side of the airfoil, the Mach number does not vary, so that these small variations are logical. We have also computed the theoretical expression of the irreversible axial velocity for a normal shock in motion, and found that the expression varies very little when the velocity of the shock varies. On the other hand, the propagation and acoustics drag, in light blue, varies strongly around zero. The flow is this time transonic, however the variations were not expected to be this strong. The induced drag, in dark blue, is very small compared to the other components, and remains almost positive during the period, as

Comparison with Gariépy's formulation

Gariépy's formulation results are presented in Figure III.15. The far-field drag in pink is here again very close to the near-field drag (in black), owing to the fact that the spurious drag is included in its definition. The latter (in orange) is however very strong, at most 15% of the total drag, which is quite close to physical components such as the wave or the unsteady ones. It also varies strongly in time. As suggested in Chapter I, Section 4.3, the reason for that is probably that the term integrated is not theoretically zero but encloses some physical drag. The order of magnitude of the wave and viscous drag components are rather similar to the ones obtained with our formulation. The induced drag is quite small and constant in time, which could be expected for a 2-D airfoil. The unsteady drag coefficient (light blue) on the contrary is strong with positive and negative values. We can also note that the sum of Gariépy's induced and unsteady components is roughly consistent with the sum of our induced and propagation and acoustics component. However, his unsteady drag coefficient encloses more than just unsteady contributions due to the buffet phenomenon. The synchronization in time raises however questions: it was indeed found that the choice of the integration surfaces, in particular the downstream extension, had again a very strong impact on the results. You can see in . For the induced drag, the extension has to be increased further in the wake in order to reach a convergence, as depicted in Figure III.16(f). This can be explained by the use of the residual volume term only rather than the whole volume term, as was already commented in Chapter II, Section 3.3.

Comparison between steady and time-averaged unsteady results

We can finally compare the unsteady results with steady results obtained using the steady formulation. The idea is to both validate the unsteady formulation and conclude whether the steady formulation is enough to predict the average drag, a question raised by designers. Table III.2 compares the averaged results. The first line stands for the average of the unsteady results presented above. The second line stands for the steady formulation applied to the average unsteady flow. The mean unsteady induced drag is close to zero, which is the value expected with a steady formulation. The wave drag is stronger for the steady calculation, resulting in a total drag higher for the steady case. The viscous drag is however lower. The spurious drag is stronger for the steady case. It can be explained by the fact that the criteria for the shock detection are not so effective on the mean flow. This comparison confirms that a study of the unsteady flow field with the steady formulation does not allow to predict the drag coefficients correctly. 3 Conclusions regarding the validity of the method

Conclusions can be drawn from these first two unsteady cases, although no rigorous validation can be achieved concerning the phenomenological breakdown.

• The total drag is accurately predicted compared to the near-field values.

• The orders of magnitude of the phenomenological components seem consistent and correct.

• The drag components do not depend on the downstream extension of the integration volumes, except for the induced and propagation and acoustics components which require large downstream extensions in order to converge towards a stable curve.

• The behavior in time is correlated to physical phenomena occurring in the flow, with time delays which can be attributed to relaxation times.

• The propagation and acoustics component is not negligible. Acoustic phenomena propagating in the flow have been highlighted, but the behavior in time cannot be straightforwardly evaluated.

• Gariépy's formulation has three main flaws: the unsteady drag component, which is not related to any physical phenomenon, is large with strong contributions, the drag components depend on the integration volumes, so that the synchronization in time cannot be assessed, and the spurious component is large, reaching the same order of magnitude as the phenomenological drag components.

• A steady drag evaluation performed on the mean flow does not allow to retrieve the mean values of the drag components observed with the unsteady method.

This first mobile test case is an airfoil oscillating in pitch in a subsonic and inviscid flow. The aim is here to focus on the induced and propagation drag components.

Quick literature review

Pitching cases have been extensively studied, whether in order to study dynamic stall on helicopters rotors, the moderate oscillations at high frequencies of a structure in aeroelasticity cases, or the aerodynamic efficiency of flapping flight. Most study focus on very large variations of the angle of attack, and solve the viscous RANS equations to model the very complex behavior of realistic flow fields. We are however focusing here on inviscid flows and small variations of the angle of attack, in order to validate our drag coefficients. Little material can be found for such simple cases.

The theoretical framework has been established by Theodorsen [START_REF] Theodorsen | General Theory of Aerodynamic Instability and the Mechanism of Flutter[END_REF] for potential flows. Pitching in an inviscid incompressible flow is expected to behave consistently with this theory. Hunsaker and Phillips [START_REF] Hunsaker | Propulsion Theory of Flapping Airfoils, Comparison with Computational Fluid Dynamics[END_REF] have compared numerical inviscid results with Theordorsen's theoretical results for the pitching of an airfoil between -5 and +5 deg in an inviscid and incompressible flow. Good agreement has been obtained for the averaged and time-dependent forces at high reduced frequency (above 0.1). Another study by Yang, Luo, and Liu [START_REF] Yang | Subsonic Flow over Unstalled Pitching Airfoil Computed by Euler Method[END_REF] compared Euler compressible simulations with experimental data for the pitching of a NACA0012 profile between -5 and +5 deg at the reduced frequency 0.2, focusing on the influence of the mean angle of attack of the pitching motion. They obtained a good agreement with the experimental data, as well as with the incompressible linear theory for the general trends.

Description of the test case

The airfoil is a NACA0012 profile, in an inviscid flow at Mach number M ∞ = 0.3. The motion is an oscillation in pitch around the first quarter of chord point. The mesh is rotated around this point as a block, the angle of attack oscillating between -5 and +5 deg at the reduced frequency k = 0.1, following the equation:

α(t) = 5 sin(2M ∞ kt) (IV.1.1)
The mesh is a very fine mesh, courtesy of Vassberg and Jameson [START_REF] Vassberg | In Pursuit of Grid Convergence for Two-Dimensional Euler Solutions[END_REF]. Several levels of refinements were tested but trouble was encountered in the vicinity of the skin, which could lead to loss of accuracy for the computation of the motion drag coefficient with coarse grids. These troubles at the skin had already been observed with our solver for Euler flows. A correction of the pressure at the skin has also been applied. The very fine mesh chosen has about 4 million elements (2048 × 2048). It extends up to 150 chord lengths around the airfoil. A visualization of the 1 node over 4 grid (512 × 512) is proposed in Figure IV. 1(a). With the benefit of insight, the effects on the motion drag component were probably very small compared to the computational effort required with such a grid. We would recommend for other studies to use coarser grids for such inviscid cases.

Convergence study

In order to further enhance the quality of the simulation near the body, the computation was first converged in a steady state with the profile fixed, before triggering the motion. A pressure correction was also applied on the boundary condition. Due to the size of the mesh, the number of time steps The curve has the shape of a butterfly, starting from the black dot (α = 0 deg), following the arrows, reaching a drag maximum before the maximum lift, then the angle of attack decreases again until zero (grey dot on the figure), and the behavior is symmetric for the negative angles. We can also note that there are two regions where the drag is negative, i.e. where the pitching motion creates thrust. Finally, there exists a point where both lift and drag are zero, but shifted compared to the angle of attack. We can explain this phenomenon by the fact that it corresponds to a quasi-steady flow at α = 0 deg delayed by the relaxation time. The shape of the curve is also very similar with the results from Yang, Luo, and Liu [START_REF] Yang | Subsonic Flow over Unstalled Pitching Airfoil Computed by Euler Method[END_REF] previously cited, although the reduced frequency is slightly higher in their study (0.2).

Analysis of the flow field resulting from the simulation

The flow field can be seen in We can see the periodicity between the positive and negative angles, but the flow is not symmetrical between the ascent and the descent at the same angle of attack. The flow field at α = 0 deg is in particular not symmetrical. This is due to the relaxation time. Note also that the flow in the vicinity of the skin is not very well predicted, despite our efforts to improve the convergence.

Application of the drag extraction method

The only volume used here is the fluid domain for the integration of induced and propagation and acoustics drag. The downstream extension is chosen as far as possible. There is no cell where the irreversible axial velocity u irr is undefined in this case. The postprocessing method is simply applied, taking the angle of attack at each time step into account. A study of the influence of the downstream extension of the integration volumes has been carried out, similarly to the applications in Chapter III, showing that the formulation gave drag curves which are independent on the extension chosen. The far-field drag, in pink, matches the near-field drag rather well. The spurious drag, in orange, which oscillates between -9 and 3 drag counts, is a little stronger than in the other unsteady cases. This could be due to the lack of accuracy in the vicinity of the skin or to the large time step chosen. The motion drag, in grey, takes values between -1 and 5 drag counts. It is largest when the angle of attack is around 0 deg, and smallest when the angle of attack is around +5 or -5 deg. This evolution is consistent since the velocity of the pitching airfoil follows the same variations (largest velocity at α = 0 deg and zero velocity at α = ±5 deg).

The induced drag varies between -10 and 30 drag counts. Its variations are quite in phase with the total drag, i.e. with a slight advance of phase compared to the angle of attack. In Figure IV.4, we can see that the variations are again strongly correlated with the variations of the square of lift, with a small time shift. The propagation and acoustics drag component is here very strong and oscillates between -20 and 65 drag counts. The synchronization in phase is identical to the induced drag's.

Comparison with Gariépy's formulation

Gariépy's formulation has also been applied to this test case. Figure IV.5 shows the corresponding results. Gariépy's formulation was written in the relative reference frame as described in the Chapter Presentation of the main existing methods of numerical drag prediction. It has been rewritten in the inertial reference frame for comparison in Chapter I, Section 4.3, so that the same motion drag is also present here. The far-field drag in pink is in good agreement with the near-field drag. The spurious drag (orange) is very small, which is logical given its definition and the fact that the flow is inviscid.

The induced drag, in dark blue, has the same order of magnitude as our formulation's, but its phase is shifted. It has been shown in Chapter III that the phase with this formulation was dependent on the extension of the integration volumes chosen. It is also true for this test case, which explains the

Influence of the reduced frequency

The reduced frequency is then decreased, from k = 0.1 to k = 0.05 and k = 0.01. The aim here is to evaluate the behavior of the drag components as the flow is closer from a steady solution. The results are also compared to the steady computations at each angle of attack. The influence of the reduced frequency on the polars of lift vs drag is shown in Figure IV.6(a). The steady drag is zero for all angles of attack since the flow is inviscid and subsonic. As the reduced frequency decreases, the relaxation of the flow is expected to be fast enough for the field to be considered quasi-static. The total drag should therefore converge to zero. On the other hand, the lift coefficient for the steady computations shows variations. We can see that the polars converge towards the steady line, the amplitude of the variations of the near-field drag being drastically reduced. This phenomenon can also be visualized in Figure IV.6(b). We can see that the curve of lift with respect to time is very little affected by the change of reduced frequency (the phase shift is reduced), whereas the drag sees its amplitude decrease, and its phase shifted to the left.

The evolution of the induced drag coefficient depending on the reduced frequency can also be seen in Figure IV.6(c). The amplitude of oscillation varies very little between k = 0.1 and k = 0.05, and is then drastically reduced for k = 0.01, whereas the phase is weakly affected. The same study has also been carried out on the propagation and acoustics component in Figure IV.6(d). Both amplitude and phase are affected for this coefficient, which is mainly responsible for the evolution of the total drag. The same behaviors are indeed observed for CD nf and CD pa .

The drag breakdown results for both reduced frequencies are also presented in Figures IV.7(a), and IV.7(b). We can see there that the order of magnitude of the induced and propagation and acoustics drag components are very close for the smaller reduced frequencies. The spurious drag coefficient is The averaged results for the drag breakdown are now compared to steady results, obtained from applying the steady formulation to the airfoil fixed at zero angle of attack. This comparison is given in Table IV.1. We can see that the steady formulation predicts no drag, except a tenth of drag count due to spurious effects, whereas the averaged total drag is more than 30 d.c. This can be explained by the energy supplied to the airfoil to perform the pitching motion, although the motion drag is very low. The main contribution to drag in average is the propagation and acoustics drag. This can be explained by the fact that we have an inviscid, subsonic flow, so that inertia, relaxation, and vorticity effects are preponderant, and are likely enclosed in this drag component. Now that the motion, induced, and propagation components have been assessed in an inviscid case, we can move on to a viscous case. The same airfoil follows the same pitching motion, in a viscous flow. The viscous drag component will therefore be added.

Quick literature review

As it has already been mentioned in Section 1, pitching cases have been extensively studied for high angles of attack and dynamic stall conditions, in order to reproduce animal flapping flights, which researchers attempt to mimic for the propulsion of micro air vehicles. Panda and Zaman [START_REF] Panda | Experimental Investigation of the Flow Field of an Oscillating Airfoil and Estimation of Lift from Wake Surveys[END_REF] have for example investigated the flow field around such airfoils using flow visualization techniques, and estimated the lift from wake measurements for several motions. Lee and Gerontakos [START_REF] Lee | Investigation of Flow over an Oscillating Airfoil[END_REF] have carried out similar studies, with the qualitative measurement of lift, drag, and moment coefficients. Finally, Gibertini et al. [START_REF] Gibertini | Experimental Method for Drag Measurement of an Oscillating Airfoil in Dynamic Stall Condition[END_REF] have also used far-field techniques to compute the drag from experimental wake measurements of an oscillating profile in dynamic stall conditions. No breakdown is however suggested.

Numerical studies have also been performed. We can cite Young and Lai [START_REF] Young | Mechanism Influencing the Efficiency of Oscillating Airfoil Propulsion[END_REF], and Xie et al. [START_REF] Xie | Computational Analysis of Propulsion Performance of Modified Pitching Motion Airfoils in Laminar Flow[END_REF], who attempted to assess the propulsion performance of such motions depending on the motion characteristics or shape of the airfoil. An unsteady optimization of the shape of the airfoil using adjoint methods has even been performed by Mani and Mavriplis [START_REF] Mani | An Unsteady Discrete Adjoint Formulation for Two-Dimensional Flow Problems with Deforming Meshes[END_REF], for small pitching amplitudes but transonic Mach numbers. As a matter of fact, little can be found in the literature concerning non-stalled motions (small amplitudes for the oscillating motion). Only Gariépy [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF] has chosen this kind of configuration for the validation of his drag breakdown method.

Description of the test case

It is the same configuration as the last case of the paper [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF], which has been described in the Chapter Presentation of the main existing methods of numerical drag prediction, Section 1.2.2. The pitching angle is oscillating between -5 and +5 deg at the reduced frequency k = 0.1, following the equation:

α(t) = 5 sin(2M ∞ kt) (IV.2.1)
The Mach number is M ∞ = 0.3 and the Reynolds number based on the free stream velocity and the chord is Re = 6.6 × 10 6 . Jameson numerical scheme and the turbulence model of Spalart-Allmaras are used. The grid used here is a 530,000 elements 2-D mesh as shown in Figure IV.8(a). The size of the computational domain is 20 chord lengths in the downstream direction and 25 chord lengths in the other directions.

This grid may seem very coarse compared to the one used for the inviscid case. As mentioned in the previous section, the very fine grid for the Euler case had been chosen to solve troubles at the skin with our solver. RANS simulations do not exhibit such issues, so that a reasonably fine mesh could be chosen. The extension chosen is also smaller, but experience with steady cases tells us that it is usually widely enough.

Convergence study

The time step chosen is 3.42 × 10 -5 s, which corresponds to 9,000 steps by period. The convergence is quickly reached after three or four periods and the drag extraction was carried out over the Again, we can see that the flow is periodic but not symmetric between the positive and negative angles of attack during ascent or descent, due to the relaxation time.

Application of the drag extraction method

The integration volume at the same instant for the viscous drag is presented in ) and are of course enclosed within the viscous surface without implying deformations. They never expand further around the airfoil, so that it was not judged necessary to show their evolution in time for this test case.

Analysis of the drag breakdown results

The curves of evolution of the drag as a function of time are presented in one positive peak and one negative peak, with a small positive phase shift. This is consistent with the paper [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF]. The far-field drag (pink) is again in very good agreement with the near-field drag (black): both curves are almost superimposed. The spurious drag in orange is indeed very small, less than a drag count. We can first focus on the orders of magnitude.

D m = Sa -ρ (u -u ∞ ) (q • n) dS (IV.2.2)
This very low value can be explained by the fact that the velocity on the upper side projected on the normal vector cancels out with the lower side, as depicted in Figure IV.11(a). The density is also present in the formula, but the flow is almost incompressible (M ∞ = 0.3), and the frequency is high so that the relaxation time prevents it from evolving too much between both sides. Density variations are also one order higher than velocity variations. The viscous drag in green varies around 100 drag counts. This value is very close to the shift of total drag between the inviscid and the viscous case, so that the order of magnitude seems correct. The induced drag in dark blue is positive with relatively small values, i.e. below 30 drag counts. It is not negligible, but remains rather small. The propagation and acoustics component, in light blue, varies more strongly, with positive and negative values. Acoustic waves are not so easily visible in this case, probably due to the size of the fluid domain.

As far as the synchronization in time is concerned, the viscous drag shows two peaks which are quite in phase with the angle of attack. The variations are not so strong since there is no separation of the flow and the angle of attack remains small. The induced drag variations are in advance of phase compared to the angle of attack. The comparison with the square of the lift coefficient, shown in A study of the influence of the grid refinement level can also be found in Appendix B, Section B.4.

Comparison with Gariépy's formulation

Difficulties were encountered with the choice of the integration volumes with Gariépy's formulation. The domain had to be reduced to 10 chords in the upstream and vertical directions, and three in the downstream direction, in order to match the far-field with the near-field. The discrepancy is probably due to the quality of the mesh. The computation on a coarser mesh confirmed that the balance between far-field and near-field required an even smaller integration domain with a coarser grid. Our formulation proved to be quite robust to any choice of domain extension with both mesh qualities. Another difficulty occurs when trying to retrieve the results shown in the paper [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF]. The choice of the integration volumes has indeed a great influence on the results and they were not specified in the original paper. The results presented in Figure IV.12 are the best achieved in this study, i.e. the closest to the results from [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF]. The viscous drag in green is stronger, varying little around 100 drag counts instead of 85. The induced drag in dark blue is not zero but remains positive. No result for the induced drag was mentioned in the article [START_REF] Gariépy | Generalization of the Far-Field Drag Decomposition Method to Unsteady Flows[END_REF], although the curves of viscous and unsteady components do not add up to give the total far-field drag curve. The phase shift of the viscous and induced drag coefficients depends on the downstream extension chosen and can therefore not be analyzed. The unsteady drag (light blue) varies between -60 and 70 drag counts instead of -15 and 110. The phase is also shifted. The spurious drag (orange) is small, at most 2 drag counts.

Influence of the reduced frequency

The reduced frequency has then be decreased to k = 0.05 and k = 0.01. The unsteady simulations have been run starting from the oscillatory state obtained with the reduced frequency k = 0.1. The same time step has been kept in order to keep reasonable restitution times, so that the convergence level may be less satisfactory as the reduced frequency is decreased. Steady computations at each angle of attack have also been performed, the drag coefficients being predicted by the steady formulation. The results for each drag coefficient, as well as the lift coefficient, are presented in The viscous drag coefficient is also relatively little affected by the change of reduced frequency, as we can see in Figure IV.13(d). The phase is almost constant, and the order of magnitude is very close to the steady results, with a gap of about 5 d.c. between unsteady and steady results. We can therefore consider that the viscous drag coefficient is correctly accounted for, since a change of frequency does not impact the viscous phenomena. The amplitude of the induced drag coefficient is strongly reduced (see Figure IV.13(e), giving even small negative values around the horizontal position. This may be due to the convergence level which is not as good as for the highest reduced frequency. The phase is however not affected by the change of frequency, contrarily to the inviscid case. The phase change observed for the total drag actually comes from the propagation and acoustics coefficient, as depicted in Figure IV.13(f). The evolution of this drag component is actually very similar to the inviscid case.

The drag breakdown results for the two reduced frequencies are also given in Figures IV.14(b), and IV.14(c). We can see how the induced and propagation and acoustics components are attenuated whereas the viscous component becomes the main contribution to drag.

The averaged unsteady results were also compared to the steady formulation, applied to the steady flow field around the airfoil at 0 deg angle of attack. The comparison could also be made for each angle of attack. One can see in Table IV.2 that the total average drag is again largely under estimated by the steady evaluation. The difference comes from both the induced and the propagation and acoustics drag components. Note that the average of CD pa is almost zero. The spurious drag is in both cases quite low. 3 Conclusions regarding the validity of the method Again, no rigorous validation of the phenomenological drag components could be made, but strong arguments supporting the validity of the method have been highlighted. The conclusions drawn in the previous applications (Chapter III, Section 3) still apply in the mobile cases. Additional conclusions can be added:

• The order of magnitude of the viscous component is consistent with the drag jump between inviscid and viscous cases.

• The motion drag is almost negligible in these cases, which can be explained by the nature of the pitching motion.

• The propagation and acoustics drag component is not negligible in these low speed mobile cases either, with negative values and an evolution in phase with the total drag.

• The behavior of the drag components with respect to the reduced frequency is consistent: they converge to quasi-static results.

• The viscous effects simply add to the inviscid behavior without changing the general behavior of the system.

• Steady evaluations at each angle of attack do not allow to correctly predict the instantaneous drag.

The unsteady formulation has been applied to pitching cases, both inviscid and viscous. The motion drag component has been validated, and observed to be negligible for such motions. The correct behavior of the viscous drag has been confirmed on the viscous case. The breakdown between induced and propagation and acoustics components is again not easily explained. Arguments towards a satisfactory estimation have however been given, such as the influence of the reduced frequency, the correlation with the lift coefficient, or the comparison with steady results. Although not strictly validated, the unsteady method is now considered ready to use on complex unsteady cases.

Chapter V

Application of the Unsteady Formulation to Complex Cases

This chapter aims at demonstrating the capabilities of the unsteady formulation on complex cases. The aimed applications are indeed flows with complex turbulent features, such as a realistic buffet case, the separated flow downstream of a spoiler, or the flow around a counter-rotating open rotor. The first step towards a realistic complex application consists in moving on to 3-D cases. The first complex test case chosen is therefore a 3-D pitching case. The second step consists in taking the time and spatial complexity of the flow into account, so that the second test case is a ZDES simulation. Since the complexity is addressed for the first time in this chapter, and due to the cost of the numerical computations, no advanced study, such as the ones carried out in the previous chapters, will be done here. The resulting evolution of the drag coefficients will only be qualitatively commented.

The first complex case is the pitching of a 3-D wing in an inviscid flow. Beyond the demonstration that the code is effective on a 3-D case, it is hoped that the induced drag component will be better understood for this test case where a wing tip vortex is present and evolves in time along with the pitching motion.

The second case is a three dimensional buffet case simulated by the ZDES model. This method allows a fine representation of a large scope of turbulent scales, and is therefore much more realistic than URANS simulations. The aim is here to evaluate the behavior of the formulation on such rich results. The drag results will also be qualitatively compared to the ones obtained previously with URANS simulations.

Application to a pitching wing in an inviscid flow

The third test case is a 3-D version of the 2-D inviscid case presented in Chapter IV, Section 1. A vortex will now appear and evolve in time with the flapping of the wing, so that the induced drag component should be affected.

Quick literature review

If little could be found in the literature concerning non-stalled 2-D cases, as pointed out in Chapter IV, even less is available concerning non-stalled 3-D inviscid cases. Most studies focus indeed on stall behavior in viscous flows. Good examples of such studies can be found in the work of Brunton et al. [START_REF] Brunton | Unsteady Aerodynamic Forces on Small-Scale Wings: Experiments, Simulations and Models[END_REF], with experimental, numerical, and modeling results and comparisons, or in the work of Campobasso et al. [START_REF] Campobasso | Turbulent Navier-Stokes Analysis of an Oscillating Wing in a Power-Extraction Regime using the Shear Stress Transport Turbulence Model[END_REF].

The paper which is the closest from our study was proposed by Neef and Hummel [START_REF] Neef | Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications[END_REF]. They numerically study a rectangular wing based on a NACA0012 profile in an inviscid flow, experiencing a plunging/pitching motion with a twisting motion at the tip. They get only thrust over the period, if the mean angle of attack is zero, but which is reduced by the presence of induced drag compared to the 2-D case.

Description of the test case

The wing considered is an elliptic wing based on the NACA0012 profile, with an aspect ratio λ = 8. The Mach number is still M ∞ = 0.3, and the pitching motion is again the same:

α(t) = 5 sin(2M ∞ kt) (V.1.1)
The mesh is an 8 million nodes grid as shown in Figure V.1(a). It is much coarser than the mesh for the 2-D case but was already very costly for an unsteady application. Lack of accuracy is therefore expected close to the skin.

Convergence study

The time step chosen is here 7.7 × 10 -5 s, corresponding to 4,000 time steps by period. The convergence is considered reached after 3 periods since the curve of lift vs drag in Figure V.1(b) is a closed loop. The same general behavior as for the 2-D cases is observed: starting at α = 0 deg on the black dot, the angle of attack increases when following the arrows. A maximum of drag is reached, before a maximum of lift. The angle of attack is then decreasing, reaching a minimum of drag, before passing in the horizontal position again at the grey dot. The symmetric behavior is observed for the negative angles. The shape of the butterfly is however slightly different from the viscous 2-D case: the thrust zones are smaller and the wings are more stretched.

Analysis of the flow field resulting from the simulation

Application of the drag extraction method

There is neither viscous nor wave volumes in this case, and the volume for the induced drag integration is the whole fluid domain. There is also no cell where the irreversible axial velocity is undefined.

Application of the drag extraction method

The drag breakdown method is applied the exact same way as in URANS cases: every 50 time steps, the integration volumes are computed using the physical criteria, and the integration of the equations defining the drag components is carried out. The integration surfaces and volumes at time step 90,000 are presented in Figure V.8(a). The shock is well enclosed in the wave volume, its wake is also well detected, as well as the boundary layer and the turbulent wake. In Appendix C, Figure C.43, you can find the evolution in time of these volumes over one period.

The criteria are however put to their limits, and a few cells are in the wrong zones (orange cells under the viscous volume or in front of the shock wave for example in Figure V.8(a)). These few cells must however not account for large contributions to the drag coefficients and their influence is neglected in this study. The behavior of the criteria should however be studied in details before moving on to industrial cases.

Comparison with URANS results

Remember that the URANS simulation in Chapter III, Section 2, has been performed at an angle of attack of 4.5 deg, since URANS methods fail to predict buffet at smaller angles of attack. Furthermore, the physics simulated in the boundary layer by both models is completely different. The results can therefore only be compared as trends.

Figure V.11 shows the comparison over one pseudo-period of the ZDES results on the left and the URANS results on the right. We can see that the trends are quite similar in terms of phase between the components and with the motion of the shock wave. ZDES gives more detailed evolutions in time when the flow is separated, i.e. between time steps 80,000 and 95,000. The viscous drag coefficient in particular is especially similar. The average wave drag is smaller in the ZDES case, which is consistent with the angle of attack being smaller. The amplitude of oscillation is also slightly smaller. The induced drag has a similar behavior in both cases. The propagation drag also has a similar behavior with a slightly stronger amplitude in the URANS case. The conclusion of this comparison is that the trends are remarkably similar between both simulation methods, and richer with the ZDES computations when the flow is detached. The induced drag takes here values between -10 and 25 drag counts. It is compared to the square of lift coefficient in Figure V.12. The correlation seems less strong in this case, with a peak of lift not observed on the induced drag coefficient. The same time shift as in the others cases is also observed.

Chapter Summary

The new unsteady formulation has been applied to two complex cases, and has demonstrated robustness and efficiency. The first test case, a 3-D complex case, has allowed to confirm the seemingly correct behavior of the induced drag component defined in Chapter II. The second test case, the ZDES simulation of a buffeting instability, has demonstrated the capabilities of the formulation on such a rich flow field. The information level which can be drawn from such drag prediction results is also quite interesting. Finally, the comparison with URANS results has confirmed the good behavior of the unsteady drag prediction method.

Comparisons with the only other unsteady drag breakdown formulation have also been carried out. The conclusions confirm the correct behavior of the wave, viscous, and motion components, and the absence of synchronization in time for the other formulation. The induced drag component has shown unexpected results, but which have been partially explained and linked with physical phenomena. We have however been unable to validate the behavior of the propagation and acoustics drag. Finally, all validations are qualitative only, since no quantitative result can be found for the unsteady drag components. This point remains an important limit to this study.

Concerning the demonstration on complex cases

Finally, the last objective was to demonstrate the capabilities of the method on complex cases. Two cases have been studied: a 3-D inviscid pitching case, and a buffet case simulated with the ZDES method. The method has proved to be robust to such difficulties. These two complex cases have also allowed confirming the correct behavior of the induced drag component, as well as the satisfactory comparison between URANS and ZDES results on similar test cases. The number of application cases and their complexity remains however limited. This general discussion leads to many perspectives, which will be described in the next section.

Conclusion and Perspectives

Synthesis

Several steps have led to the fulfillment of the previously described objectives. Each step has brought out important results hereafter summarized:

• A new proof of Van der Vooren's steady formulation has been devised. This rigorous demonstration should remain a starting point for other attempts in the development of unsteady far-field drag formulations.

• An unsteady far-field drag breakdown formulation has been developed in a rigorous manner. Two versions are available, one with four components, and one with five components.

• New drag components have been identified and isolated, namely a motion drag linked to the motion of a body, and a propagation and acoustics drag linked to acoustic effects. Other possible sources of unsteady drag have also been discussed.

• A ready-to-use postprocessing tool has been implemented. A Python routine has also been written in order to compute the evolution in time of the drag coefficients over a given period. Good practice recommendations have finally been formulated.

• The estimation of the wave, viscous, induced, and motion drag components has been validated on simple cases, fixed and mobile. The orders of magnitude as well as the synchronization in time with the flow phenomena has indeed be judged correct. The propagation and acoustics component probably encloses several phenomena and is therefore difficult to validate.

• The results have been compared to the results using Gariépy's formulation, the only other unsteady drag breakdown method available, showing that the same orders of magnitude were obtained. The synchronization in time could however not be compared since it depended on the extension of the integration volumes with this other formulation, unlike with ours.

• The method has been validated with several turbulence models, and several fidelity levels of flow simulation, proving its robustness for all methods considered.

• The unsteady drag breakdown formulation has allowed to better understand the phenomena at stake for all test cases considered, as well as the drag creation processes. It has also confirmed that steady estimation could not be used, even for the estimation of the average drag levels. Finally, new postprocessing perspectives have been open, in particular for the spectral analysis of a flow using a global variable.

Perspectives

Limits to these positive results have been identified in the General Discussion. They lead towards a potential for improvement, as well as widening perspectives:

Concerning the theoretical formulation and the physical understanding

• An interesting perspective would be to try to combine this formulation with the ones based on the velocity vector summarized in the first part of this thesis, in order to improve the theoretical definition of the induced drag component, and get rid of the definition by default of this component.

• The criterion used for the definition of the integration surface for the computation of the wave component could be generalized to unsteady flows in a better way. Nothing was found in the literature for such developments.

• The propagation and acoustics drag component could be further explored, for example by numerically computing the drag in the presence of an acoustic source only. A further decomposition of the term, by isolating contributions rather related to non-acoustic propagation, could also be attempted. An analytic evaluation of the drag created by an acoustic dipole is being studied and will be published in a future paper.

Concerning the numerical implementation

• The loop currently used for the computation of drag at each time step could be directly included within the postprocessing code, and even within the CFD code, so that no storage of the flow field would be required.

• The time derivatives could be directly extracted from the data used in the CFD code, in order to be able to apply the method to more complex motions, with deformations of the mesh.

Concerning the aerospace applications

• Cases with unsteady experimental data for the drag components, if achievable, would be a means to finally validate the unsteady breakdown. Experimental data of pitching cases isolating each drag components (similarly to the steady cases used in this study) may be a first step towards this validation.

• The method could be applied to the detached flow downstream of a spoiler. Simulations with the ZDES method are already available at Onera for this test case. It may be interesting for the aircraft designers to better understand the role of the viscous and induced drag components in this configuration.

• The unsteady formulation could also be useful for aeroelasticity problems, such as flutter. The understanding of this complex instability driven by the interaction between the stiffness of the structure and the aerodynamic loading may indeed benefit from such drag breakdown analysis.

• The unsteady loading experienced by flapping wings is another promising application case. Flapping flight is indeed still being extensively studied, and only clues have been raised so far.

• Unsteady optimization of one of the drag component may also be achievable with this formulation. Optimization of the total unsteady drag is already being carried out by other research teams. We can easily imagine to use the wave drag for example in an optimization problem. The objective could be a target time evolution of the coefficient, or to minimize the average value or even the peak values. Adjoint methods would however not be a realistic path, since the computation of the sensitivities with this unsteady formulation would require a huge effort.

• Finally, the most promising application case for the aircraft designer in 2015 is the performance assessment of a Counter-Rotating-Open-Rotor (CROR). It is indeed crucial to finally have access to reliable data concerning the performance of such a breakthrough configuration, before the technology leap can be seriously contemplated for a commercial use. This tool could indeed provide information concerning the location of the sources of drag, the interactions between the front and the rear propellers, and could help enhance the global performance.

The value commonly used for RANS models is P r t = 0.9. Thus, classical RANS models based on the Boussinesq assumption try to express µ t with the mean field to close the system. Many ways have been developed in the literature: algebraic models, models with transport equations of one or more variables, etc.

In the following paragraph, the Spalart-Allmaras and kω turbulence models are presented in more details.

A.1.1.2 The Spalart-Allmaras turbulence model

The Spalart-Allmaras model [START_REF] Spalart | A One-Equation Turbulence Model for Aerodynamic Flows[END_REF] is a one-equation turbulence model based on the Boussinesq assumption, which focuses on the determination of µ t as a function of the mean field. The model proposes to solve a transport equation for the quantity ν:

∂ ν ∂t +u j ∂ ν ∂x j = c b1 (1-f t2 ) S ν -c w1 f w - c b1 K 2 f t2 ν d 2 + 1 σ ∂ ∂x j (ν + ν) ∂ ν ∂x j + c b2 ∂ ν ∂x i ∂ ν ∂x i (A.6)
The different terms of the right hand side have been determined heuristically and calibrated from increasingly complex flows. Details can be found in [START_REF] Spalart | A One-Equation Turbulence Model for Aerodynamic Flows[END_REF].

This model is successfully used in the industry, but often tends to predict too massive separations, due to an overestimation of the turbulent viscosity in the vortex center. The turbulent eddy viscosity is then expressed as µ t = ρk ω . The transport equation for the rate of viscous dissipation ω can however lack robustness far from the wall, so that Menter [START_REF] Menter | Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications[END_REF] proposed a correction: he suggests to replace this transport equation by the transport equation for the viscous dissipation ǫ far from the walls. The Shear Stress Transport (SST) correction also allows to bound the value of µ t so as to avoid overestimation of shearing stress in presence of positive pressure gradients.

A.1.1.4 Unsteady RANS

The RANS equations are actually time-dependent. If a resolved time step is used to simulate an unsteady flow, one can therefore access to an unsteady flow field. This is often referred to as Unsteady RANS (URANS), in order to insist on the unsteady feature of such simulations. The resulting timedependent flow field must however be carefully looked into, since the turbulence models used to close the equations are valid only as long as the time over which these changes in the mean flow occur is large compared to the time scales of the turbulent motion containing most of the energy.

31 I 1

 311 1.1.2 Jones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.1.3 Oswatitsch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.1.4 Maskell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.1.5 Van der Vooren and Destarac . . . . . . . . . . . . . . . . . . . . . . . . 12 1.2 First unsteady generalization of Van der Vooren's formulation . . . . . . . . . . . 18 1.2.1 Theoretical developments . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.2.2 Results on a pitching case . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2 Presentation of the formulations based on the velocity vector . . . . . . . . . . . . . . . 25 2.1 Formulation for steady incompressible flows . . . . . . . . . . . . . . . . . . . . . . 25 2.2 Breakdown into induced and profile components for steady incompressible cases . 25 2.3 Extension to steady compressible flows . . . . . . . . . . . . . . . . . . . . . . . . 28 2.4 Breakdown in the steady compressible case . . . . . . . . . . . . . . . . . . . . . . 28 2.5 Generalization to unsteady flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.5.1 Noca . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.5.2 Wu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.5.3 Marongiu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.5.4 Xu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.5.5 Other contributions to the unsteady generalization . . . . . . . . . . . . Development of an Unsteady Formulation starting from Van der Vooren's Formulation 33 What is difficult about a generalization to unsteady flows? . . . . . . . . . . . . . . . . 35 1.1 To account for all the additional terms due to unsteadiness and relate them to phenomenological components . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 1.2 To avoid applying the steady theory as such . . . . . . . . . . . . . . . . . . . . . 35 1.3 To take the propagation delays into account . . . . . . . . . . . . . . . . . . . . . 35 2 Derivation of a directly generalizable proof of Van der Vooren's steady formulation . . . 37 2.1 Derivation of the far-field equation in the steady case . . . . . . . . . . . . . . . . 37

  u irr axial velocity under irreversible flow assumptions u rev axial velocity under reversible flow assumptions u † irr irreversible axial velocity defined from u rev u * rev enthalpy-corrected axial velocity under reversible flow assumptions u * irr irreversible axial velocity defined from u * rev ∆u irreversible axial velocity defect defined by Van der Vooren and Destarac [90]: ∆u = u irru ∞ D p pressure drag D f friction drag D nf near-field drag D vw profile drag D vw = D w + D v D w

Figure 1 -

 1 Figure 1 -Visualization thanks to favorable atmospheric conditions of the phenomenological sources of drag

Figure 2 -

 2 Figure 2 -Evolution of the profile drag coefficients with respect to the position of the wake plane by Jones using two variants of his formulation, and Betz's formulation (extracted from [42])

Figure 3 -

 3 Figure 3 -Control volume and boundaries as defined by Destarac and Van der Vooren (extracted from [90])

  (a) Visualization of the bump shapes on the profile (b) Results of the drag breakdown

Figure 5 -Figure 6 -

 56 Figure 5 -Study of the induced drag reduction thanks to unconventional wingtips using Van der Vooren's formulation (extracted from [34])

Figure 7 -Figure 8 -Figure 9 -Figure 10 -

 78910 Figure 7 -Examples of designs obtained from the multi-objective optimization of a wing using a Taylor expansion of Van der Vooren's far-field breakdown method (extracted from [105])

  (a) Far-field and near-field drag components (b) Viscous and unsteady drag components

Figure 11 -

 11 Figure 11 -Results of the drag prediction of the DPW5 case using Van der Vooren's formulation (extracted from [31])

Figure 12 -

 12 Figure 12 -Contours of the z component of the Lamb vector on a sectional plane of the delta-wing (extracted from [107])

  ), and 13(b).

Figure 13 -Figure 14 -

 1314 Figure 13 -Results of the drag breakdown using the steady incompressible velocity-based formulation (extracted from [107])

Figure 15 -

 15 Figure 15 -Drag coefficients for several Mach numbers computed by the steady compressible velocity-based formulation (extracted from [59])

Figure 16 -

 16 Figure 16 -Total drag computed by the unsteady incompressible velocity-based formulation applied on experimental data with respect to time (extracted from [65])

( a )Figure 17 -

 a17 Figure 17 -Evolution of the drag coefficients defined by Marongiu with respect to the angle of attack (extracted from [52])

  Figure I.1(a) at time t = t i . Now let us consider a perturbation of this shock wave, say that it suddenly moves downstream at time t = t i + ∆t as sketched in Figure I.1(b). This perturbation of the shock wave induces a perturbation of the flow which is shed downstream in the wake of the wave, inside the streamtube, such as in Figure I.1(c).

  Figure I.1 -Schematic representation of the propagation in the wake of a shock wave of a perturbation due to the motion of the wave

  e. a flow free of vortices (see Figure I.3): p = p ∞ and v = w = 0.

Figure I. 3 -

 3 Figure I.3 -Schematic representation of the computation of the irreversible axial velocity on a wake plane S d

Figure I. 4 -

 4 Figure I.4 -Schematic representation of the streamtubes enclosing the shock wave and the body with its boundary layer

Figure I. 5 -

 5 Figure I.5 -Streamtube enclosing the shock wave with the surfaces used in the demonstration

Figure I. 6 -

 6 Figure I.6 -Integration surfaces in the steady case

Figure I. 7 -

 7 Figure I.7 -Surfaces and normal vectors used to demonstrate the force applied by a shock wave on a body

Figure

  Figure I.8 shows an example of integration surfaces computed on a wing-body-tail configuration.The frontiers visible inside the surfaces are simply due to the frontiers between the different grid types, since the case is a chimera mesh. The box for the induced drag integration stops at a given distance from the most downstream position of the wings trailing edge. The tail of the aircraft is also included in the surface, with a margin chosen as a given distance from the skin.

Figure I. 8 -

 8 Figure I.8 -Integration surfaces for the CRM with tail [36]: S p w (red), S p v (green) and S p i (blue)

Figure I. 9 -

 9 Figure I.9 -Streamtube enclosing the unsteady shock wave with the surfaces used in the demonstration

Figure I. 10 -

 10 Figure I.10 -Integration volumes and surfaces in the unsteady case

  and the complementary downstream surface S cd = S d \ (S wd ∪ S vd ) (see Figure I.11).

Figure I. 11 -

 11 Figure I.11 -Complementary volumes and surfaces used for the demonstration of the unsteady induced drag component

Figure I. 12 -

 12 Figure I.12 -Integration volumes for the CRM with tail[START_REF] Hue | Fifth Drag Prediction Workshop: Computational Fluid Dynamics Studies Carried Out at ONERA[END_REF] with the unsteady criteria: V p w (red), V wd (orange), V v (green) and V c (blue)

Figure I. 13 -

 13 Figure I.13 -Schematic representation of the obstruction process due to the increase of the entropy variation across a shock wave

Figure I. 14 -

 14 Figure I.14 -Schematic representation of a shock wave in a tube traveling in the upstream direction

Figure II. 1 -

 1 Figure II.1 -Schematic representation of the computation of the reversible axial velocity on a wake plane S d

  Figure II.2(b). The simulation reaches convergence after 2,000 iterations. The near-field drag coefficient is then CD nf = 83.52 d.c. and varies less than a thousandth of drag counts. This value for the total drag is quite close from the one achieved by Vassberg and Jameson with their finest grid (83.42 d.c.). A grid convergence study has been carried out and is presented in Appendix B, Section B.1.

Figure II. 2 -

 2 Figure II.2 -Mesh and convergence curves for the steady Euler transonic case

  Figure II.3(b).

  Figure II.3 -Flow field and integration surface for the steady viscous subsonic case

Figure II. 4 -

 4 Figure II.4 -Evolution of the wave drag coefficient with respect to the downstream extension of the integration surface for the steady Euler transonic case

  Figure II.7(a) that the boundary layer and the wake are completely enclosed within the surface. There exist cells where u irr is undefined in this case. They are located in the first cell layer around the skin, on the upper and the lower side of the airfoil, leading edge and trailing edge excluded (see Figure II.7(b)).

  Figure II.5 -Mesh and convergence curves for the steady viscous subsonic case

  Figure II.7 -Integration surface and cells where u irr is undefined for the steady viscous subsonic case

Figure II. 8 - 1 .

 81 Figure II.8 -Evolution of the wave drag coefficient with respect to the downstream extension of the integration surface for the steady viscous subsonic case

Figure II. 9 -

 9 Figure II.9 -Mesh and convergence curves for the steady Euler subsonic 3-D case

Figure II. 10 -

 10 Figure II.10 -Mach contours, isosurface of Q criterion colored by the vorticity magnitude, and stream ribbons for the steady Euler subsonic 3-D case

  Figure II.12(a). It extends up to 30 chord lengths around the wing in every direction. The aerodynamic conditions are: M ∞ = 0.8, α = 2.5 deg and Re = 2.7 × 10 6 based on the free stream velocity and the reference chord.

Figure II. 11 -

 11 Figure II.11 -Evolution of the wave drag coefficient with respect to the downstream extension of the integration surface for the steady Euler subsonic 3-D case

  Figure II.12 -Mesh and convergence curves for the steady viscous transonic 3-D case

  (a), and II.15(b), allowing a comparison between the drag coefficients obtained with the three expressions of the axial velocity defect. For the viscous and induced drag coefficients, the downstream extension of the integration surface expands quickly in the wake, whereas for the wave drag coefficient the surface is slowly increased with a certain amount of cell layers. The converged values of drag are around 104 d.c. for the wave drag, 154 d.c. for the viscous drag and 27 d.c. for the induced drag. The lift is 0.256, so that the theoretical induced drag for an elliptic load is CD th i = CL 2 πλ = 26 d.c.

( a )

 a Figure II.15 -Evolution of the wave, viscous and induced drag coefficients with respect to the downstream extension of the integration surfaces for the steady viscous transonic 3-D case

Figure II. 16 -

 16 Figure II.16 -Cells activating the unsteady criterion for the subsonic vortex shedding case, without filtering

  . The steady criterion on the other hand seems to give satisfactory results. The time evolution of the S p w surfaces is presented in Appendix C, Figure C.33.

  Lovely and Haimes also observed these falsely tagged cells in the whole fluid domain. They explained it by the errors in the time derivative which become of the same order of magnitude as the time derivative itself. The zones in red in Figure II.16 and Figure II.17 also correspond to zones where the acoustics may play a role, as we will see in Section 3.2.

Figure II. 17 -

 17 Figure II.17 -Cells activating the unsteady criterion for the buffet case, without filtering

  Figure II.18 -Cells activating the steady criterion (red), and cells activating the unsteady criterion for the buffet case, with filtering (black)

  3.28. The definition should not depend on the choice of the wake integration plane S d . Let us consider a wake plane S ′ d located further downstream. The volume between these two planes is denoted V d , as in Figure II.19.

Figure II. 19 -

 19 Figure II.19 -Schematic representation of the effect of a change in the wake surface on the induced drag component

Figure II. 20 -

 20 Figure II.20 -Visualization in the flow field of the volume term of the unsteady induced drag component for the vortex shedding case

Figure II. 21 -

 21 Figure II.21 -Visualization in the flow field of the volume term of the unsteady induced drag component for the buffet case

  Figure III.1(b) which is an almost closed loop.

  Figure III.1 -Mesh and convergence curve for the vortex shedding case

Figure III. 2 -

 2 Figure III.2 -Flow field visualization taken at time step 2,000 for the vortex shedding case

Figures

  Figures retracing the evolution in time of this vorticity field are presented in Appendix C, Figure C.27. In these animated figures, we can see how the pockets of vorticity detach from the airfoil and interact with each other.

Figure III. 5 -

 5 Figure III.5 -Evolution of the drag coefficients with respect to time for the vortex shedding case

  Figure III.6(b). The evolution in time of this drag component is however very difficult to validate.

  Figure III.6 -Assessment of the induced and acoustic drag components for the vortex shedding case

Figure III. 7 -

 7 Figure III.7 -Evolution of the drag coefficients defined by Gariépy with respect to time for the vortex shedding case

  Figure III.8 -Effect of the downstream extension on the drag coefficients for the vortex shedding case (in distance from the leading edge, recalling that the chord is 1)

  Figure III.9(b) shows that the periodic regime is almost reached, since the loop of lift as a function of drag is almost closed.

  Figure III.9 -Mesh and convergence curve for the buffet case

Figures III. 10

 10 Figures III.10(a), and III.10(b) show the instantaneous flow field taken at two different time steps. The time evolution of the vorticity field can be seen in Appendix C, Figure C.31. A strong shock appears on the upper side which oscillates in time as a consequence of the interaction between the shock wave and the boundary layer. A separation occurs right downstream of the shock when it is in an upstream position as seen in Figure III.10(b).

  Figures III.10(a), and III.10(b) show the instantaneous flow field taken at two different time steps. The time evolution of the vorticity field can be seen in Appendix C, Figure C.31. A strong shock appears on the upper side which oscillates in time as a consequence of the interaction between the shock wave and the boundary layer. A separation occurs right downstream of the shock when it is in an upstream position as seen in Figure III.10(b).

Figure III. 10 -

 10 Figure III.10 -Flow field visualization taken at time steps 500 and 725 for the buffet case

  (a), and III.11(b). Recall that the integration surface for the induced drag is the outer surface, and V c is the complementary of the volumes shown in these figures in the whole domain. The evolution along the period is also presented in Appendix C, Figure C.33. It confirms that the criteria are capable to trace the shock wave and the evolution of the boundary layer.

Figure III. 11 -

 11 Figure III.11 -Integration volumes taken at time steps 500 and 725 for the buffet case

Figure III. 12 -

 12 Figure III.12 -Cells where u irr is undefined taken at time steps 500 and 725 for the buffet case

  expected. Now let us focus on the synchronization in time. The wave drag coefficient is strongest around time step 350, and lowest around time step 750. The evolution of the C p curves at the skin is presented in Appendix C, Figure C.32. We can see that the most downstream position occurs around time step

Figure III. 13 -

 13 Figure III.13 -Evolution of the drag coefficients with respect to time for the buffet case

  Figure III.14 -Assessment of the induced and acoustic drag components for the buffet case

Figure III. 15 -

 15 Figure III.15 -Evolution of the drag coefficients defined by Gariépy with respect to time for the buffet case

  Figure III.16(a) the effect of the downstream extension on the wave drag coefficient, in Figure III.16(c) the effect on the viscous one, and in Figure III.16(e) on the unsteady one. The order of magnitude is almost not affected but the phase shifts as the downstream extension changes. In Figure III.15, the downstream extension for the viscous volume was chosen of 3 chords (or 0.69 m in Figure III.16(c)), 11 for the unsteady volume, and S w for the wave surface. This choice has been made in order to get the least spurious drag. Our definitions of the viscous and wave drag coefficients on the contrary give quite stable results (see Figures III.16(b), and III.16(d))

  Figure III.16 -Effect of the downstream extension on the drag coefficients for the buffet case (distance in m from the leading edge, recalling that the chord is 0.23 m)

(a) 1

 1 Figure IV.1 -Mesh and convergence curve for the inviscid pitching case

Figure IV. 2 .

 2 The evolution in time of the Mach contours can also be seen in Appendix C, Figure C.36. The C p curves at the skin for the unsteady mobile computation and the steady computation at the same angle of attack is also provided in Appendix C, Figure C.37.

Figure IV. 2 -

 2 Figure IV.2 -Flow field visualization taken at time step 500 (5 deg) for the inviscid pitching case

Figure

  Figure IV.3 shows the results of the drag breakdown over a period. The total drag, in black, exhibits two positive peaks, which correspond to the largest angles of attack α = 5 deg and α = -5 deg. We can first note that the drag coefficients exhibit a negative phase shift compared to the angle of attack, whereas the lift coefficient, in red in the bottom figure, has a small positive phase shift.The far-field drag, in pink, matches the near-field drag rather well. The spurious drag, in orange, which oscillates between -9 and 3 drag counts, is a little stronger than in the other unsteady cases. This could be due to the lack of accuracy in the vicinity of the skin or to the large time step chosen. The motion drag, in grey, takes values between -1 and 5 drag counts. It is largest when the angle of attack is around 0 deg, and smallest when the angle of attack is around +5 or -5 deg. This evolution is consistent since the velocity of the pitching airfoil follows the same variations (largest velocity at α = 0 deg and zero velocity at α = ±5 deg).The induced drag varies between -10 and 30 drag counts. Its variations are quite in phase with the total drag, i.e. with a slight advance of phase compared to the angle of attack. In Figure IV.4, we can see that the variations are again strongly correlated with the variations of the square of lift, with a small time shift. The propagation and acoustics drag component is here very strong and oscillates between -20 and 65 drag counts. The synchronization in phase is identical to the induced drag's.

Figure IV. 3 -

 3 Figure IV.3 -Evolution of the drag coefficients with respect to time for the inviscid pitching case

  Figure IV.6 -Influence of the reduced frequency on the drag and lift coefficients for the viscous pitching case

  Figure IV.7 -Drag breakdown results for two smaller reduced frequencies for the inviscid pitching case

  Figure IV.8 -Mesh and convergence curve for the viscous pitching case

2. 4

 4 Analysis of the flow field resulting from the simulation A snapshot of the flow field is shown in Figure IV.9(a). It is very close to the inviscid case, with the addition of the boundary layer and the viscous wake. There is no flow separation even at the highest angles of attack. The evolution in time of the Mach contours is given in Appendix C, Figure C.38. The C p curves at the skin for the unsteady mobile computation and the steady computation at the same angle of attack is also provided in Appendix C, Figure C.39.

  Figure IV.9(b). The evolution in time of this volume is also presented in Appendix C, Figure C.40. It matches the boundary layer and the wake nicely along the pitching motion. The cells where u irr is undefined are concentrated very close to the skin (see Figure IV.9(c)

  Figure IV.10. The total drag over one period has two peaks, one for each maximal angle (-5 deg and +5 deg), with again a negative phase shift. It varies between 65 and 190 drag counts. The lift, on the other hand, has

  Figure IV.9 -Flow field and integration surfaces taken at time step 2250 (5 deg) of the viscous pitching case

Figure IV. 10 -

 10 Figure IV.10 -Evolution of the drag coefficients with respect to time for the viscous pitching case

  Figure IV.11 -Assessment of the motion and induced drag components for the viscous pitching case

  Figure IV.11(b), confirms the strong correlation with a phase shift, which can again be explained by the effect of second order irreversible phenomena. The negative values of the propagation and acoustics component correspond to a thrust occurring as the airfoil is moving upwards or downwards.

Figure IV. 12 -

 12 Figure IV.12 -Evolution of the drag coefficients defined by Gariépy with respect to time for the viscous pitching case

Figure IV. 13 .

 13 We can first see in FigureIV.13(a) that the polar of lift as a function of drag converges towards the steady curve. The magnitude of drag variations is in particular strongly attenuated. The lift coefficient is very little affected by the change of frequency, as we can see in FigureIV.13(b). Its phase shift with the angle of attack is reduced, and the order of magnitude converges towards the steady

  Figure IV.13 -Influence of the reduced frequency on the drag and lift coefficients for the viscous pitching case

  Figure V.1 -Mesh and convergence curve for the 3-D pitching case

  The flow field is shown in Figure V.2 with a visualization in a cross plane on the left and a 3-D one on the right. The Mach contours in a cross plane (Figure V.2(a)) are very similar to the 2-D ones, except in the vicinity of the skin where unexpected jumps are visible. It is the same issue as in the 2-D case, but a very fine grid could not be used in this case for computational cost reasons. A wing tip vortex is visible in Figure V.2(b) with the help of the Q-criterion. As the wing pitches, the vortex sequentially vanishes, then reappears with the vorticity oriented in the opposite direction, before vanishing again. This evolution in time is described in Appendix C, Figure C.41.

Figure V. 2 -

 2 Figure V.2 -Flow field visualizations for the 3-D pitching case

  Figure V.5 -Mesh visualizations for the ZDES buffet case

  Figure V.7 -Visualization of the instantaneous flow field at time step 90,000 for the ZDES buffet case

Figure V. 9 -

 9 Figure V.9 -Evolution of the drag coefficients with respect to time over ten periods for the ZDES buffet case

Figure V. 11 -

 11 Figure V.11 -Comparison of ZDES and URANS results over one pseudo-period for the buffet case

Figure V. 12 -

 12 Figure V.12 -Induced drag coefficient and square of the lift coefficient for the ZDES buffet case

A. 1 . 1 . 3 k

 113 ω model Wilcox[START_REF] Wilcox | Reassessment of the Scale-Determining Equation for Advanced Turbulence Models[END_REF] suggested a two-equations turbulence model. Coupled transport equations for the turbulent kinetic energy k and the rate of viscous dissipation ω are added to the model. They can be written as:

1 4 2 36 I. 2

 4362 Visualization thanks to favorable atmospheric conditions of the phenomenological sources of drag . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Evolution of the profile drag coefficients with respect to the position of the wake plane by Jones using two variants of his formulation, and Betz's formulation . . . . . . . . . 11 3 Control volume and integration surfaces as defined by Destarac and Van der Vooren . 13 4 Study of the wave drag reduction thanks to bumps using Van der Vooren's formulation 15 5 Study of the induced drag reduction thanks to unconventional wingtips using Van der Vooren's formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 6 Results of the multi-objective optimization of a wing using a Taylor expansion of Van der Vooren's far-field breakdown method . . . . . . . . . . . . . . . . . . . . . . . . . 18 7 Examples of designs obtained from the multi-objective optimization of a wing using a Taylor expansion of Van der Vooren's far-field breakdown method . . . . . . . . . . . 19 8 Integration surfaces for the DPW5 case computed by the ONERA code ffd72 . . . . . 19 9 Results of the drag prediction of the DPW5 case using Van der Vooren's formulation 20 10 Schematic representation of the moving and rotating control volume as defined by Gariépy 21 11 Results of the drag prediction of the DPW5 case using Van der Vooren's formulation 23 12 Contours of the z component of the Lamb vector on a sectional plane of the delta-wing 26 13 Results of the drag breakdown using the steady incompressible velocity-based formulation 27 14 Evolution of the induced drag coefficient with respect to the downstream extension by Lamb vector integral compared to Maskell's formulation . . . . . . . . . . . . . . . . . 27 15 Drag coefficients for several Mach numbers computed by the steady compressible velocitybased formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 16 Total drag computed by the unsteady incompressible velocity-based formulation applied on experimental data with respect to time . . . . . . . . . . . . . . . . . . . . . . . . 30 17 Evolution of the drag coefficients defined by Marongiu with respect to the angle of attack 31 I.1 Schematic representation of the propagation in the wake of a shock wave of a perturbation due to the motion of the wave . . . . . . . . . . . . . . . . . . . . . . . . . . . Control volume, surfaces and normal vectors . . . . . . . . . . . . . . . . . . . . . . . 38 I.3 Schematic representation of the computation of the irreversible axial velocity on a wake plane S d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 I.4 Schematic representation of the streamtubes enclosing the shock wave and the body with its boundary layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 I.5 Streamtube enclosing the shock wave with the surfaces used in the demonstration . . 43

  state

	Superscript
	p	practice (integration surfaces and volumes used in practice in the postprocessing tool)
	Acronyms
	CFD	Computational Fluid Dynamics
	DNS	Direct Numerical Simulation
	RANS	Reynolds Averaged Navier Stokes
	URANS Unsteady Reynolds Averaged Navier Stokes
	DES	Detached Eddy Simulation
	ZDES	Zonal Detached Eddy Simulation
	d.c.	drag count (10 -4 )

  By assumption, the flow is isentropic in V wd so ds dt = 0 in V wd . The second term is however non-zero, contrarily to the steady case. Starting from the definition of the stagnation enthalpy H =

									ρ	ds dt	∂u irr ∂s	+	dH dt	∂u irr ∂H	(I.3.15)
	e +	p ρ	+	u 2 2	, one can write a local balance of stagnation enthalpy:
					ρ	dH dt	= ρ	de dt	+	dp dt	+ ρ	d dt	u 2 2	(I.3.16)

  .1.

	Problem	General case	Special case	Special recommendation
	u irr or u † irr ?	u irr	u † irr if the zone of non	use a larger downstream ex-
			definition deforms the	tension
			viscous surface
	which wave crite-	steady criterion	unsteady wave criterion	remove spurious cells by hand
	rion?		if the steady one fails
	downstream	more than 3 chord
	extension?	lengths

Table II.1 -Good practice recommendations for the robustness of the method

Table III

 III

.1 -Comparison between steady and averaged unsteady results for the vortex shedding case

  Table III.2 -Comparison between steady and averaged unsteady results for the buffet case

	formulation/flow	CD nf	CD w	CD v	CD i	CD pa	CD sp	CD ff
	average/unsteady	620.96	218.73	391.75	11.18	-11.46	10.50	610.91
	(±	standard	(± 101.99) (± 35.65) (± 122.04) (± 13.69) (± 101.37) (± 6.75) (± 102.66)
	deviation)
	steady/average	640.47	302.03	341.36	-	-	-56.22	696.69

  Table IV.1 -Comparison between steady and averaged unsteady results for the inviscid pitching case

	formulation/flow	CD nf	CD m	CD i	CD pa	CD sp	CD ff
	average/unsteady	31.97	2.50	11.44	21.54	-3.50	35.47
	(± standard deviation) (± 52.27) (± 2.04) (± 16.18) (± 32.83) (± 3.89) (± 48.45)
	steady/average (0 deg)	0.08	-	-	-	0.08	0.00

  Table IV.2 -Comparison between steady and averaged unsteady results for the viscous pitching case

	formulation/flow	CD nf	CD m	CD v	CD i	CD pa	CD sp	CD ff
	average/unsteady 127.56	-0.04	99.65	18.47	9.69	-0.22	127.78
	(±	standard	(± 44.94) (± 0.07) (± 3.76) (± 10.95) (± 30.58) (± 0.31) (± 45.01)
	deviation)
	steady/average	96.39	-	96.34	-	-	0.05	96.35
	(0 deg)

Historical presentation of the main thermodynamic methods

Application to a vortex shedding case

Application to a buffet case simulated by a URANS method

Application to a pitching airfoil in an inviscid flow

Application to a pitching wing in an inviscid flow Figure V.4 -Induced drag coefficient and square of the lift coefficient for the 3-D pitching case

Analysis of the flow field resulting from the simulation

The flow field is presented in Figure II. [START_REF] Brunton | Unsteady Aerodynamic Forces on Small-Scale Wings: Experiments, Simulations and Models[END_REF]. The shock wave on the upper side is quite strong. The Q-criterion defined by Hunt [START_REF] Hunt | Eddies, Stream, and Convergence Zones in Turbulent Flows[END_REF] allows locating vortices in a flow field. An isosurface of this criterion is plotted in Figure II.13, highlighting the wingtip vortex in this case.

Application of the drag extraction method

The flow field is post-processed after 10,000 iterations. Figure II.14(a) shows an example of the integration surfaces for a given downstream extension. The cells where u irr is undefined, shown in black in Figure II.14(b), are located very close to the skin and a little downstream of the shock. They are entirely enclosed within the viscous surface, so that there is no trouble for using u irr .

Application of the drag extraction method

Once the flow field obtained from the numerical simulation and stored every 50 time steps, the postprocessing method can be applied. The first step consists in computing the integration volumes using the criteria described in Chapter I and Chapter II, with a downstream extension of 20 chord lengths. They are computed at each extraction of drag. An example of the viscous volume is shown in Figure III.3. The induced surface is there again the outer surface and the integration volume is the complementary of the viscous volume. [START_REF] Gariépy | A New Axial Velocity Defect Formulation for a Far-Field Drag Decomposition Method[END_REF]. In spite of the relatively high number of cells concerned, we can see that the automatic inclusion within the viscous surface does not deform it. There is therefore no need to use u † irr .

Figure III.4 -Cells where u irr is undefined taken at time step 2,000 for the vortex shedding case

Analysis of the drag breakdown results

The final formulation derived in Chapter II, Section 4.1 can now be applied. The drag extraction is carried out over one period. The resulting time evolution curves are shown in Figure III.5. The total drag varies between 1825 and 1975 drag counts. The total far-field drag (pink) is in good agreement with the near-field drag (black). Accordingly, the spurious drag in orange is quite small, around 1% of the total drag.

As far as the breakdown into phenomenological components is concerned, a validation is not achievable. We can only try to evaluate the order of magnitudes and the behavior in time compared to the evolution of the flow features. The main contribution to drag is here the viscous drag (in green) with around 1800 d.c., as was expected. Both induced and propagation and acoustics drag components, in dark and light blue respectively, are very small in comparison, around 50 d.c. each.

The synchronization between the far-field drag components and the phenomena occurring in the flow can also be discussed. The viscous drag is the strongest around time step 2400. If we look at the

Chapter Summary

The unsteady formulation has been successfully applied to two naturally unsteady test cases, a vortex shedding and a buffet case. The wave and viscous components have been analyzed in terms of both orders of magnitude and synchronization in time. Their behavior is considered validated. The strong unsteady induced drag, as defined in Chapter I, is surprising but can be explained by the breakdown suggested in Chapter II into acoustic effects and lift-induced effects. The comparison with steady results has also proved the benefit of such methods. The formulation will now be assessed on mobile cases, with the apparition of new drag components.

Chapter IV

Assessment of the Motion, Induced, and Propagation Drag Components on Mobile Cases

The formulation can now be applied to mobile cases, in order to validate the motion, induced, and propagation and acoustics drag components. Again, there is no other data available concerning the evolution in time of the drag coefficients. Their behavior will therefore be assessed in terms of order of magnitude and synchronization in time. Gariépy's results, as well as steady computations at each position will also be used for comparison.

The first mobile test case is an airfoil pitching around its axis in an inviscid and subsonic flow. The numerical configuration is first described, then the evolution in time of the flow field is presented. The drag breakdown method is then applied. The behavior of the resulting drag components is extensively analyzed, with a study of the influence of the reduced frequency. Comparisons with Gariépy's formulation are also performed, with a study of the dependence on the integration volume chosen. The results are also compared to steady computations at each angle of attack.

The second mobile test case consists of the exact same pitching motion, in a viscous flow. The viscous drag component will therefore add up, allowing us to confirm the breakdown into viscous and induced components. The same studies are performed, in terms of influence of the reduced frequency, extension of the integration volumes, and comparison with another formulation and steady evaluations. From these two studies are drawn conclusions concerning the validity of the breakdown method.

Analysis of the drag breakdown results

The results of the drag extraction are visible in Figure V.3. The results for the total drag look similar to the Euler 2-D case although the negative or thrust zones are smaller and the curves have a smoother variation. The same negative phase shift with the angle of attack is observed. The near-field drag takes values between -12 and 101 drag counts. The far-field drag is again in good agreement, with a spurious drag between -3 and 5 drag counts. The lift is smaller than in the 2-D case, but exhibits the same small positive time shift with respect to the angle of attack. The phase shifts observed for both near-field lift and drag coefficients are consistent in trend with the results from [START_REF] Neef | Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications[END_REF]. The motion drag is almost zero. As far as the induced drag coefficient is concerned, the downstream extension of the integration volume has to be reduced since the wing tip vortex is quickly dissipated by numerical effects. In this case, the corresponding dissipated drag is actually attributed to the spurious component. Closer to the body, the definition chosen for the induced component is however dependent on the volume chosen, especially for the phase. A compromise was found at 15 chords downstream of the body, where the phase is almost converged but the induced drag is not too dissipated. The corresponding induced drag takes values between -10 and 80 drag counts. The strongest values correspond to instants when the wing tip vortex is the strongest. The evolution of the induced drag coefficient is also compared to the theoretical steady expression of the induced drag as a function of the square of lift in Figure V.4. The same strong correlation, with a negative time shift, is again observed. There is also a shift in the order of magnitude which can be attributed to the numerical dissipation of the wingtip vortex.

The propagation drag is not negligible, varying between -10 and 20 drag counts, but the dominant drag component is now the induced one, which was expected for a 3-D case.

A study of the influence of the grid refinement level can also be found in Appendix B, Section B.5.

Application to a buffet case simulated by the ZDES method

The aim of this last test case is to demonstrate the capabilities of the formulation on a complex case and compare the results to those obtained in Chapter III, Section 2 with URANS simulations. The ZDES, from Zonal Detached Eddy Simulation, is a powerful simulation method developed by Deck [START_REF] Deck | Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF]. The characteristics and benefits of the method are exposed in Appendix A, Section A.1.3.2.

Quick literature review

We have seen in Chapter III, Section 2 that the buffet phenomenon had been studied by several authors with URANS simulations. Although the comparisons of global variables with experimental data was satisfactory (rms, pressure distribution), we have seen that the angle of attack needed to be increased to 4.5 deg instead of 3.5 to trigger the instability, and that the fine temporal features are not predicted with such methods (typically, only the first frequency of a spectrum is predicted). Authors such as Fan et al. [START_REF] Fan | Hybrid LES/RANS Simulation of a Shock Wave/Boundary Layer Interaction[END_REF] have therefore attempted to use DES methods to predict the detailed unsteady behavior of the flow. The agreement with the experimental results was however disappointing, and the angle of attack also had to be increased to 4 deg. Deck then used the ZDES method to study the buffet over an airfoil [START_REF] Deck | Numerical Simulation of Transonic Buffet over a Supercritical Airfoil[END_REF], comparing the results to experimental data by Jacquin et al. [START_REF] Jacquin | Experimental Study of Shock Oscillation over a Transonic Supercritical Profile[END_REF]. He obtained the instability at 3.5 deg, and very good agreement with the experiments, despite a relatively coarse grid compared to today's standards. These fine computations allow retrieving the rich temporal content of the flow which can help us better understand this complex instability.

Note finally that the buffet over a 3-D configuration is much more complex: in 2-D, the motion of the shock wave over the airfoil is periodic, whereas no single characteristic frequency can be observed in 3-D, and most of time there is a broad band frequency distribution which presents high instability levels. This is due to the the strong interaction between the turbulence of the separation downstream of the shock and the shock itself on the upper side of wings. One can easily conclude that 2-D and 3-D buffet are completely different phenomena since the 2-D one is a kind of resonant flow and the 3-D one can be more considered as an important excitation of the separation which is very energetic and pilots the shock location and motion. 3-D buffet has been studied with ZDES simulations by Brunet and Deck [START_REF] Brunet | Zonal-Detached Eddy Simulation of Transonic Buffet on a Civil Aircrfat Type Configuration[END_REF].

We will focus on 2-D buffet only in this preliminary test case.

Description of the test case

The case is the transonic buffet over the supercritical OAT15A profile with chord 0.23 m at M ∞ = 0.73, α = 3.5 and Re = 12.14 × 10 6 . The time step is 5 × 10 -7 s, so that there are about 29,000 time steps in one period. The mesh is a 2 million nodes grid shown in Figure V.5.

The DES zones are located around and downstream of the shock above the upper side of the airfoil and in the wake, as depicted in Figure V. 6(a). This corresponds to the ZDES mode 2. More details can be found in [START_REF] Deck | Numerical Simulation of Transonic Buffet over a Supercritical Airfoil[END_REF].

Convergence study

In Figure V.6(b), we can see the polars of lift with respect to drag over ten periods. Variations in the lift and drag values do occur between two periods of movement of the shock wave. The position of the shock wave is actually not quite periodic. The loops keep however the same general shape, so that the convergence is considered reached for the ten periods considered. The time derivatives, computed by finite differences between two time steps as described in Appendix A, Section A.2, were not accurate with only a one time step difference. This can be explained by the fact that the time step is very small for ZDES computations, so that the flow evolves very little in between and numerical errors become the same order of magnitude as the derivative itself. A distance of 10 time steps has therefore been chosen for the computation of the time derivatives, which is a good compromise for the accuracy.

Analysis of the drag breakdown results

Figure V.9 shows the results of the drag extraction over ten pseudo-periods. We can first notice that a secondary frequency appears when the flow is detached, and vanishes as the flow reattaches. The drag breakdown will be analyzed in further details over one period.

Spectral analysis

A spectral analysis has also been performed in order to investigate the possibilities offered by such rich information on the evolution of the drag coefficients, and to compare the frequencies obtained for the drag coefficients with the ones measured in the flow field. In Figure V.10(a), we can see a power spectral density function of the wave drag coefficient. This can be compared to Figure V.10(b), where the power spectral density function of the pressure fluctuations obtained by Deck [START_REF] Deck | Numerical Simulation of Transonic Buffet over a Supercritical Airfoil[END_REF] is plotted: the first four peaks have the same frequency. We can conclude that the wave drag coefficient is strongly correlated with the pressure fluctuations measured on the upper skin of the airfoil.

General Discussion

The objective of this thesis was to develop a formulation for the phenomenological breakdown of drag for unsteady flows, to implement it in a postprocessing code, to validate it on simple test cases, and finally to apply it to complex cases. All steps of the method have now been carried out. The fulfillment of the objectives can now be discussed, before concluding this work.

Concerning the theoretical formulation

The first objective was to obtain a formulation as rigorous and robust as possible, and to ensure the synchronization in time and the independence on the integration volumes for each phenomenological drag component. To do so, a new proof of Van der Vooren's formulation has been derived, and then been carefully generalized to unsteady flows. The main flaw of the demonstration is however the definition by default of the unsteady induced drag. This starting point will always remain shaky. The physical background of the induced drag component has then been improved, with the highlighting of acoustic effects. A propagation and acoustics drag component has been isolated. Doubts however remain concerning the validity of the breakdown, as well as the physical background of the remaining term which has been attributed to the induced drag component. In order to improve the robustness of the method, an alternative expression of the ill-defined irreversible axial velocity has been studied and compared on several test cases, proving to be less robust and not as accurate. An unsteady extension of the criterion used for the definition of the wave surface has also been studied, with disappointing results. These attempts to improve the robustness have therefore been dropped, and the robustness remains an issue.

Concerning the numerical implementation

The second objective was to implement this new formulation in a postprocessing code. It has been done as an extension of the Onera code ffd72. There are however limits for the resulting postprocessing tool: the code has to be run entirely for each time step, so that there is a great waste in computational cost. The computation of the time derivatives, defined for the time being as a finite difference between two instant of the flow at a given grid node, has been corrected in order to account for rigid body motions, but cannot be used for deforming meshes. There is no theoretical restriction for that, it is only a computational issue.

Concerning the validation on test cases

The third objective consists in validating the formulation on simple test cases. Two motionless and two simple pitching cases have been chosen to achieve this goal. The orders of magnitude, synchronization in time, influence of the viscous effects, and influence of the frequency have been investigated.

Appendices

Appendix A

Numerical tools used for the applications

This appendix quickly describes the tools used for the applications, from the modeling and simulation of aerodynamics, to the codes used.

A.1 Modeling of aerodynamics

If the viscosity is neglected, the governing equations for the flow field can be simplified to the Euler equations, with no turbulence involved. However, if the viscosity cannot be neglected, turbulent features appear for aerospace applications operating at high Reynolds numbers. A turbulent flow is dominated by vorticity dynamics, which explains the chaotic character of the flow when the non-linear mechanisms are preponderant. As a result, a completely turbulent flow presents many different spatial and temporal scales which are correlated.

In the Kolmogorov cascade theory, the major part of the kinetic energy is carried by the biggest structures, which interact with the mean motion and whose kinetic energy is too high to be dissipated by the viscosity. That energy is transferred to the smaller structures by a mechanism of eddy stretching, up to a certain scale of eddy, the scale of Kolmogorov, permitting the viscous dissipation. This smallest length-scale for eddy structures is noted η. We can distinguish three ranges of scales:

• the integral scales, enclosing the biggest structures, carrying the major part of the energy and interacting with the mean flow

• the inertial scales which for the mechanism of cascade happens

• the dissipative scales

The resolution of the Navier-Stokes equations up to the smallest scales implies many constraints of enforcement. The computational cost of such a Direct Numerical Simulation (DNS) becomes rapidly prohibitive and must be restricted to simple cases at low Reynolds numbers, up to the thousand. For example, Spalart et al [START_REF] Spalart | Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[END_REF] indicate that a DNS on a complete aircraft would require about 10 16 mesh points. The current computing capabilities impose to model turbulent interactions. Two main approaches are currently used, the Reynolds Averaged Navier-Stokes (RANS) approach and the Large Eddy Simulation (LES) approach. These two approaches are both based on a separation of scales and are described in the next sections.

The main idea to solve the Navier-Stokes equations without calculating all the spectra of turbulent structures is to separate variables in the following manner :

with f the solved part of the variable and f ′ the unsolved part of the variable. It represents a filtering of the equations by the operator . For compressible flows, it is useful to introduce the Favre decomposition:

Additional terms appear when this decomposition is introduced into the Navier-Stokes equations, among which the turbulent stress tensor τ tur ij which requires modeling.

A.1.1 RANS approach

In the RANS approach, the scale separation operator is the statistical average. In practice, it is assimilated to a temporal average operator, by doing the ergodicity hypothesis. Thus, the solved part of the field is the mean field and the unsolved part of the field is the time-varying field. The turbulent stress tensor becomes:

In the various RANS models, it only remains to model the Reynolds stress tensor to close the system of averaged Navier-stokes equations.

A.1.1.1 The Boussinesq hypothesis

Most of common RANS models are based on the concept of turbulent viscosity developed by Boussinesq in 1897 [START_REF] Boussinesq | Théorie de l'Écoulement Tourbillonnant et Tumultueux des Liquides[END_REF]. The turbulent viscosity enables to link the Reynolds stress tensor to the mean field by the equation:

It can be observed the analogy with the kinematic viscosity, but it is important to notice that the kinematic viscosity is a property of the fluid whereas the turbulent viscosity is a property of the flow. The turbulent conduction is evaluated thanks to the turbulent Prandtl number:

A.1.2 LES Approach

In the RANS approach, all the turbulent scales are modeled, which provides very low computing costs. On the contrary, in the LES approach, a big part of these scales is solved. The largest structures interacting with the mean flow, as well as a part of structures in the inertial range are solved, whereas the smallest dissipative ones are modeled. It is assumed that, in the inertial range and beyond, the turbulence is locally isotropic, justifying the modeling of the smallest turbulent length scales. The scale-separator operator is submitted to a low-pass filter of width ∆ (in the wave number space).

This width is dependent on the grid refinement. To capture a certain length scale of turbulent structures, it is required to get a sufficiently refined mesh, at least as small as the length scale. In some ways, the grid acts as a filter. Many different filters exist in the literature, the most common of them being:

• The "box" or top-hat filter, providing a spatial average on a given length. It is a cardinal sinus in the spectral space.

• The "spectral or sharp cutoff" filter, representing a cardinal sinus in the physical space but is an ideal low-pass filter in the spectral space (and owns the idempotence property).

• The Gaussian filter whose representations in the physical and spectral space is a Gaussian distribution.

Once the filtering applied with the Favre decomposition, the term τ tur ij has to be modeled. In the LES approach, τ tur ij can be noted τ SGS ij , where SGS means Sub-Grid Scales, as the modeling deals with scales lower than ∆. It exists various ways to model this tensor that have been classified by Sagaut [START_REF] Sagaut | Large-Eddy Simulation for Incompressible Flows[END_REF].

In the functional models all the tensor τ SGS ij is modeled in the aim to reproduce the behavior of the sub-grid turbulence and its interactions with solved scales. Such an approach requires a very fine analysis of the physics of turbulent interactions and is largely based on the Kolmogorov's cascade. A strong hypothesis often used is that the effect of the subgrid scales on the resolved scales is essentially dissipative, the action of subgrid scales can be summarized by the energy equilibrium with the resolved scales. It is then assumed that subgrid scales have locally an isotropic behavior which is not necessarily the case in fully turbulent flows other than isotropic homogeneous turbulence. In this frame, it yields to introduce the notion of subgrid viscosity µ sgs =< ρ > ν sgs to model the energy transfer mechanisms from the resolved to the subgrid scales, by analogy with the molecular mechanisms used in the RANS approach:

A closure equation to evaluate ν sgs is now required. Many subgrid viscosity models have been developed in the literature, among which the Smagorinsky model will be used later. It can be summarized as such:

• Only the forward energy cascade process is taken into account, i.e. the action of subgrid scales on the resolved scales (backward process) is neglected.

• A characteristic length l 0 and a characteristic time t 0 are sufficient for describing the subgrid scales.

• There exists a total separation between the subgrid and resolved scales.

• The flow is in constant spectral equilibrium, there is no accumulation of energy at any frequency.

That leads to the following expression of the subgrid viscosity ν sgs :

The constant C s can be determined by different manners but is often estimated assuming an equilibrium between the production and dissipation terms in isotropic homogeneous turbulence. The common value is C s = 0.18, but it can be tuned according to the application, between 0.1 and 0.2.

A.1.3 Hybrid RANS/LES approaches

The LES approach requires a computing cost which is still prohibitive at high Reynolds numbers with the current calculation capabilities. In order to reduce this CPU cost, the idea has been advanced to couple the RANS and the LES approaches to solve complex flows. Thus, the hybrid approach combines the efficiency of the RANS and the high resolution level of the LES. The hybrid RANS/LES methods take profit of the similarity of filtered equations form of the two approaches which only differ by the scale separation operator. Despite the mechanisms modeled in the tensor τ SGS ij and τ RAN S ij are very different, the common formalism enables to switch from one modeling to the other. Switching from one model to another one is equivalent to change the effective filtering, or in other words, the characteristic length scale. All the challenge is to manage the transfer from one model to the other, as they do not represent the same physics.

A very large number of hybrid RANS/LES strategies have been developed over the last two decades, as SDM, LNS, VLES, XLES, DES, SAS, DDES, ZDES,... The reader is invited to report to the review of Sagaut, Deck and Terracol [START_REF] Sagaut | Multiscale and Multiresolution Approaches in Turbulence[END_REF] to get a quite exhaustive insight on these methods.

In the following, only the Detached Eddy Simulations, and its extension Zonal Detached Eddy Simulation, methods are described.

A.1.3.1 Detached Eddy Simulation

One of the most popular hybrid approach for simulating industrial flows is the Detached Eddy Simulation (DES), originally developed by Spalart [START_REF] Spalart | Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[END_REF]. Considering that RANS modeling is very efficient close to the wall to simulate an attached boundary layer, but exhibits many difficulties to provide realistic solutions of massively separated flows, and that LES requirements are prohibitive in the boundary layer, the main idea of this approach lies in the modification of the transport equation of the turbulent variables depending on the wall-distance, imposing thus a RANS model close to the wall and a LES model further. In the originated DES of Spalart et al. [START_REF] Spalart | Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[END_REF], the initial turbulent viscosity of the RANS model, close to the wall, is modified toward a subgrid viscosity of Smagorinsky type in the free stream and flow separated areas, thanks to the introduction of a new characteristic length d, which replaces the wall distance d w in the pseudo-viscosity transport equation. This length scale is defined as:

with ∆ = ∆ max = max(∆ x , ∆ y , ∆ z ) the characteristic mesh length and C DES a constant settled to 0.65.

The definition of d based on ∆ max induces a switching from the RANS modeling to the LES modeling according to the local mesh refinement. It allows the use of RANS grid standards close to the wall (which is a very big advantage of this approach) and to satisfy LES grid standards far from the solid walls. Nevertheless the transition area between the two models remains uncertain when d w ≈ ∆, the "grey zone" as called by the authors of DES [START_REF] Spalart | Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[END_REF]. In particular, if this grey zone is located inside the boundary layer, i.e. ∆ x or ∆ z is equivalent to the boundary layer thickness δ, a drop of the eddy viscosity in this area can occur. As the mesh is not yet adapted to the LES mode and the eddy viscosity does not decrease fast enough to allow the eddies development, it yields an under-estimation of the Reynolds stress levels.

A.1.3.2 Zonal Detached Eddy Simulation

Inspired from the previous approach, the Zonal Detached Eddy Simulation (ZDES) was first proposed by Deck, and the complete formulation that proposes an efficient solution to prevent delay in the formation of instabilities has been recently published in [START_REF] Deck | Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF].

This approach takes full advantage of its zonal nature, not only to allow the user to specify RANS and LES regions, but also to allow the use of various formulations within the same calculation. ZDES offers thus an attractive flexibility in the treatment of turbulent flows in technical applications and has been applied often with good results over a wide range of Mach numbers and configurations (see [START_REF] Deck | High-Fidelity Simulations of Unsteady Civil Aircraft Aerodynamics: Stakes and Perspectives. Application of Zonal Detached Eddy Simulation[END_REF]).

Three specific hybrid length scale formulations, also called modes, are then optimized to be employed on three typical flow field topologies (see Figure A.14). Mode 1 concerns flows where the separation is triggered by a relatively abrupt variation in the geometry, mode 2 is retained when the location of separation is induced by a pressure gradient on a gently-curved surface, and is the one used in this work, and mode 3 for flows where the separation is strongly influenced by the dynamics of the incoming boundary layer.

A.2 Codes used

Two codes have mainly been used in this study, one for the flow simulation, and one for the postprocessing, following the scheme presented in Figure A. [START_REF] Caruana | Buffet and Buffeting Active Control[END_REF].

The Onera elsA code [START_REF] Cambier | The Onera elsA CFD Software: Input from Research and Feedback from Industry[END_REF] is an object-oriented CFD software. The conservative variables of the flow field, computed and extracted by elsA at each instant, are used as input for the postprocessing tool. The time derivative in particular can be extracted from elsA, but, in the present work, they are computed as a finite difference between the field at a given grid node at two consecutive time steps, and a correction is added in the case of a rigid body rotation or translation of the mesh:

where w and Ω denote the velocity and rotation vector of the grid respectively. More complex motions, with deforming grids for example, can therefore not be treated with this numerical setup so far. The Onera far-field drag post-processing tool dates back to 1993, with the first ffd, for farfield drag, code. A cell-by-cell approach with physical sensors was developed and a correction for spurious vortex diffusion was implemented one year later. Face-based management was adopted in 2003 for the drag extraction to allow the post-processing of flow solutions computed on both structured and unstructured grids. Multiblock structured data are converted into unstructured solutions. This feature is particularly useful for treating numerical solutions from various flow solvers, especially for cell-centered solvers (elsA) and cell-vertex solvers (DLR-TAU [START_REF] Schwamborn | The DLR TAU-Code: Recent Applications in Research and Industry[END_REF]).

The ffd72 code contains about 300 subroutines for 75,000 lines before the present development. The unsteady formulation developed in Chapters I, and II has been implemented in the ffd72 code. The developments represent about 20 new routines, and about 1,000 lines. As far as the practical use of the unsteady method is concerned, the user can specify the time step if the time derivatives are computed manually. The output of ffd72 is then a text file containing the results of the drag breakdown at the given instant. A small Python program has been written in order to run ffd72 for each time step, read the output file, and add the results to a data file containing at the end the evolution in time of all the drag components.

Appendix B

Grid studies

This appendix summarizes the grid studies carried out during this thesis. All test cases could not benefit from such grid studies, however two or three levels of grid density have been considered for five cases.

B.1 Airfoil in a steady transonic inviscid flow

The first convergence study concerns the steady test case used to assess the behavior of the wave drag component (see Chapter II, Section 1.5.1). The flow is inviscid, with a transonic Mach number M ∞ = 0.8. Three levels of grid density are considered here, courtesy of Vassberg and Jameson [START_REF] Vassberg | In Pursuit of Grid Convergence for Two-Dimensional Euler Solutions[END_REF]. The coarse mesh shown in Figure B.17 The convergence of the near-field drag coefficients for the three grid levels is depicted in Figure B.18(a). All three computations are converged to the thousandth of drag count. The fine grid seems to be very close to the converged value for the near-field drag, since there is less than half a drag count between the medium and the fine grid. The far-field estimation of the wave drag, as the downstream extension of the integration surface is set to vary, is shown in Figure B.18(b). We can see again that the fine grid offers a satisfactory level of convergence. What is remarkable is the little dependence on the refinement level of the drag computing using u irr . The computations using u † irr or u * irr are much more dependent. For the coarse and the medium grids, CD † w and CD * w do not converge to a given value but keep increasing with the number of cell layers. The discrepancy between CD † w and CD * w close to the shock wave is also stronger with the coarse mesh. , for which the wing tip vortex dissipation is also dependent on the mesh quality. The near-field and far-field drag coefficients seem however to converge towards a common value about 108 drag counts. An even finer grid would have been better suited. We can however conclude that again CD † i and CD * i are even further from CD i in the vicinity of the trailing edge when the mesh is coarser. CD * i shows also larger discrepancies with CD † i with the coarse grid. The drag breakdown method has been applied on these two refinement levels. The corresponding results can be seen in Figure B.24. We can first notice that all drag components vary little with the grid refinement level. The total far-field drag (in pink) is in particular almost identical, whereas the near-field drag varies more: the spurious drag, in orange, is indeed smaller with the fine grid, which confirms that far-field prediction allows to put aside a part of the numerical artificial drag, even for unsteady cases. A small phase shift can be observed for the viscous (green) and induced drag (dark blue) components, which may be explained by the frontier between viscous and complementary volumes which is less accurate with a coarse grid. The propagation and acoustics component (light blue) sees its amplitude decrease with a finer mesh. The motion drag, in grey, is very slightly affected. The evolution in time of the drag components for both grids is presented in Figure B.26. The grid dependence is here stronger, since the vortex decay for this 3-D flow is more sensitive to grid refinement as we can see with the near-field drag in black. The total far-field drag, in pink, is however remarkably stable or in other words, the spurious drag (orange) is much stronger with the coarse grid. The induced drag component (dark blue) sees its maximal peaks reduced with a coarse grid, which is probably due to the vortex decay. The propagation and acoustics drag (light blue) is accordingly increased. There is actually a transfer from one drag component to the other due to the dissipation of the wingtip vortex. The motion drag (grey) is again very slightly affected.

List of Tables

II.1

Good practice recommendations for the robustness of the method . . . . . . . . . . . Prediction and phenomenological breakdown of drag for unsteady flows Accurate drag prediction is now of a major issue for aircraft designers. Its phenomenological sources need to be identified and quantified for an efficient design process. Far-field methods, which allow such phenomenological drag breakdown, are however restricted to steady flows.

This study consists in developing a far-field drag prediction method aiming at a phenomenological breakdown of drag for unsteady flows. The first step has consisted in generalizing the steady formulation of Van der Vooren to unsteady flows, starting from a new rigorous proof. Axes for the improvement of the robustness and physical background have then been explored. Acoustic contributions have in particular been highlighted and quantified. The resulting five-components formulation has then been applied to simple cases, in order to validate as best as possible the phenomenological breakdown. The behavior of the drag components has proved to be consistent with the physics of the flow. Finally, the method has been applied to complex cases in order to demonstrate its capabilities: a 3-D case and a flow simulated by the ZDES method.

In the future, it would be interesting to further improve the definition of the induced drag component, for example by using velocity-based formulations. As far as the application cases are concerned, the performance evaluation of a Counter-Rotating-Open-Rotor would strongly benefit from such a method. Unsteady optimization of one of the drag component could also be contemplated. Finally, applications in aeroelasticity or flapping flight would be an interesting perspective.

Keywords: DRAG; FAR-FIELD; PERFORMANCE; AERODYNAMICS; CFD; UNSTEADY Prévision et décomposition phénoménologique de la traînée pour des écoulements instationnaires L'estimation précise de la traînée est aujourd'hui un enjeu majeur pour les avionneurs. Il est nécessaire d'identifier et de quantifier ses sources phénoménologiques dans le cadre d'un processus de design efficace. Les méthodes champ lointain, qui permettent une telle décomposition de la traînée, sont cependant limitées aux applications stationnaires.

Cette étude consiste à développer une méthode d'extraction champ lointain destinée à permettre une décomposition phénoménologique de la traînée pour des écoulements instationnaires. La première étape a consisté à généraliser la formulation stationnaire de Van der Vooren aux écoulements instationnaires, en partant d'une nouvelle démonstration rigoureuse de sa méthode. Des axes pour l'amélioration de la robustesse et du contenu physique ont ensuite été explorés. Des contributions acoustiques ont en particulier été mises en évidence et quantifiées. La formulation à cinq composantes ainsi obtenue a ensuite été appliquée à des cas tests simples, dans le but de valider aussi bien que possible la décomposition phénoménologique. Le comportement des composantes de traînée s'est avéré cohérent avec la physique de l'écoulement. Enfin, la méthode a été appliquée à des cas complexes afin de démontrer ses capacités : un cas instationnaire 3D ainsi qu'un écoulement simulé en ZDES.

Dans l'avenir, il serait intéressant de continuer à explorer la définition de la composante de traînée induite, par exemple en utilisant les formulations basées sur le vecteur vitesse. En ce qui concerne les cas d'application, l'évaluation de la performance d'un doublet d'hélices contra-rotatives pourrait fortement bénéficier de l'utilisation d'une méthode comme celle-ci. L'optimisation instationnaire d'une des composantes de traînée pourrait également être envisagée. Enfin, des applications en aéroélasticité ou en vol d'ailes battantes pourraient être d'intéressantes perspectives.

Mots