domain: phys spi

text:
Clément Caillaud
email: clement.caillaud@cea.fr

Mathieu Lugrin
email: mathieu.lugrin@onera.fr

Sébastien Esquieu
email: sebastien.esquieu@cea.fr

Cédric Content
email: cedric.content@onera.fr

Global stability analysis of a hypersonic cone-cylinder-flare geometry

Characterizing the boundary layer transition to turbulence around realistic hypersonic vehicles is a challenging task due to the numerous parameters that affect the process. To address this challenge, the conecylinder-flare (CCF) geometry has been designed to provide a flow topology that captures various transition mechanisms observed on reentry objects [6], such as absolute and convective instabilities, which are dependent on the free stream conditions. In this study, a global linear stability analysis is performed on the CCF model at M ∞ = 6.0 to investigate and map the dominant instabilities at wind tunnel flow conditions. We examine optimal responses and forcings computed using resolvent analysis, as well as global modes originating from the recirculation bubble at the cylinder flare junction. The effects of bluntness are assessed through analyses of both blunt and sharp configurations. Our results shed light on the linear flow mechanisms that promote the transition to turbulence around such hypersonic objects.

INTRODUCTION

The current CCF geometry consists in a 5 • sharp or blunt cone, followed by a cylindrical section which terminates to a flare with an angle of θ f = 10 • (CCF10). The geometry is depicted in Fig. 1 and can be divided into two primary regions: the cone section, which provides a canonical convectively unstable boundary layer that supports first and second mode waves, and the cylinder-flare section, which in-cludes a separation bubble generated by the adverse pressure gradient caused by the shock-wave at the cylinder-flare junction. This flow topology is similar to that described in previous works by [START_REF] Mathieu Lugrin | Transition scenario in hypersonic axisymmetrical compression ramp flow[END_REF][START_REF] Anubhav Dwivedi | Oblique transition in hypersonic double-wedge flow[END_REF]. In this second region, the separation bubble has been found to be globally unstable for flare angles θ f ≥ 8 • and Reynolds numbers above Re = 11.5 × 10 6 [START_REF] Li | Nonlinear Evolution of Instabilities in a Laminar Separation Bubble at Hypersonic Mach Number[END_REF]. At this critical Reynolds number, convectively unstable modes such as first-mode and second mode waves can also exist along the recirculation bubble and/or the reattachment point. For cases where either the convective, the global instabilities or a combination of both, lead to transition, the non-linear terms will induce a strong modification of the bubble. Therefore making the baseflow stability analysis irrelevant and leading to a mean flow analysis [START_REF] Marxen | The effect of small-amplitude convective disturbances on the size and bursting of a laminar separation bubble[END_REF][START_REF] Mathieu Lugrin | Transition scenario in hypersonic axisymmetrical compression ramp flow[END_REF]. However, all the cases treated in this article stay fully laminar in the whole domain and thus base flow stability analysis is justified.

For these flare angles, the previous experimental campaigns in quiet and noisy tunnels, associated with numerical investigations using PSE, revealed that both first and second-mode waves are amplified along the sharp cone region. On the cylinder-flare region, second mode instabilities were not found to be growing until the bubble reattachment whereas first-mode waves were seen to be continuously amplified along the bubble [START_REF] Benitez | Instability Measurements on an Axisymmetric Separation Bubble at Mach 6[END_REF][START_REF] Paredes | Boundary-Layer Instabilities over a Cone-Cylinder-Flare Model at Mach 6[END_REF]. Additionally, infrared thermography revealed the presence of steady streaks originating at the reattachment and amplifying along the flare at an azimuthal wavenumber of m = 36. The origin of these streaks remain poorly understood but is suggested to be the byproduct of non-linear interaction in the bubble [START_REF] Li | Nonlinear Evolution of Instabilities in a Laminar Separation Bubble at Hypersonic Mach Number[END_REF]. Considering these previous results, we employ a global linear stability framework to explore the transition dynamics of the CCF10 object. Our aims are twofold: first, characterise the base flow and identify the instabilities present in the separation bubble, and second, analyse the forced response of the linear dynamics through resolvent analysis. To achieve these goals, we organise this paper as follows: in Section 2 describes the numerical stability tools and the base flow, along with a validation. In Section 3 are discussed the global modes found in the bubble, and the Section 4 investigates the forced response of the linear dynamics. Finally, Section 5 draw conclusions and discuss the future directions of this research.

LINEAR STABILITY

The analysis is supported by the stability toolbox BROADCAST [START_REF] Poulain | BROAD-CAST: A high-order compressible CFD toolbox for stability and sensitivity using Algorithmic Differentiation[END_REF]. It uses high order finite-volume schemes and algorithmic differentiation to compute fixed points of the compressible Navier-Stokes equations and extract the associated direct and adjoint global linear operators, along with its derivatives up to arbitrary order if required. Using this toolbox, the methodology is introduced in the next paragraphs. For further details about the framework presented in this section, the reader is referred to the comprehensive review of [START_REF] Sipp | Dynamics and Control of Global Instabilities in Open-Flows: A Linearized Approach[END_REF].

Considering a domain in cylindrical coordinates (x, r, θ) discretised in a structured fashion with N p points. The flow dynamics around the CCF geome-try are governed by the following discrete non-linear dynamical system,

∂q q q ∂t = R(q q q) + f e f e f e , (1)
where q q q is the conservative state vector of N v variables, R R R is the compressible Navier-Stokes operator and f e f e f e is a harmonic exogenous forcing of small amplitude. Supposing the existence of a fixed-point q q q 0 such that R R R(q q q 0 ) = 0, the non-linear system can be linearised around this steady state in order to retrieve the linear dynamics of small disturbances q q q ′ , with q q q = q q q 0 + ϵq q q ′ , ϵ ≪ 1. The linear disturbances equations reads, ∂q q q ′ ∂t = Lq q q ′ + f e f e f e , with L = ∂R R R(q q q) ∂q q q q q q0 ,

with L, the Jacobian of the non linear operator around q q q 0 . To analyse the disturbances dynamics, we make the hypothesis of a zero-angle of attack, and consider the problem to be 2D axisymmetric.

Hence the disturbances and forcing vectors are expressed in terms of their frequency ω and azimuthal Fourier modes m, resulting in the harmonic ansatzs, q q q ′ (x, r, θ, t) = q q q(x, r)e i(mθ+ωt) ,

f f f e (x, r, θ, t) = f f f (x, r)e i(mθ+ωt) . (3)
Using these vectors, the linear dynamics of Eq. 2 will be studied in two ways. First, by considering the autonomous system, i.e f e f e f e = 0. In Fourier space, the global stability of the baseflow at a given frequency and wavenumber m (with this formalism, the Jacobian operator becomes dependent on m) can be studied by solving the eigenvalue problem, L(m)q q q = iωq q q, ω ∈ C,

with the growth-rate of an eigenfunction given by ℜ(ω) and its frequency by ℑ(ω).

Next, the stability of the forced system, i.e f e f e f e ̸ = 0, can be studied. Due to the non-normality of the jacobian operator L, the exogenous forcing can trigger the non-modal amplification of responses even if the system is stable [START_REF] Peter | Nonmodal Stability Theory[END_REF]. These non-normal mechanisms are given by the resolvent operator defined as R R R = (iωI -L) -1 . For noise-amplifier flows, this operator yields an input-output relation at a given frequency and wavenumber between forcings and responses over the baseflow, such that,

q q q = R R R(ω, m) f f f , ω ∈ R. (6)
To find the most amplified convective instabilities at a pair (ω, m), an optimal decomposition of the resolvent matrix can be computed and provides an orthogonal basis of optimal forcings and responses ranked by energy. Starting from Eq. 6, such functions can be found by solving an optimisation problem. Using the discrete norms ||q q q|| E = q q q * W E q q q and || f

f f || F = f f f * P * W f P f f f
, where W E is the Chu energy weight matrix [START_REF] Chu | On the energy transfer to small disturbances in fluid flow (Part I)[END_REF] and W F is taken as the identity along with a restriction matrix P used to impose the forcing on specific regions or variables, a Rayleigh quotient can be obtained,

µ 2 0 = sup f f f ||q q q|| E || f f f || F = sup f f f ||R f f f || E || f f f || F . (7)
The optimal forcing that satisfies equation 7 may then be found by solving the eigenvalue problem :

P * R R R * W E R R RP f f f i = µ 2 i W f f f f i . (8)
Where (•) * is the hermitian transpose and the eigenvalues µ 2 0 > ... > µ 2 i > µ 2 i+1 > ... of Eq. 8 are the optimal gains ranked by energy and their associated eigenfunctions f f f opt i . Solving Eq. 8 for various values of (ω, m) allows to map the system resonance peaks and characterise the receptivity of the baseflow introduced in the next section.

BASEFLOW COMPUTATION

The laminar flow around the CCF geometry comprises a detached shock, expansion waves, separation and reattachment shocks around the bubble.

Case

Re These features poses a challenge to the computation of fixed-points for different Reynolds numbers. This section outlines the flow conditions and details the procedure used to compute fixed points. Finally the baseflows are described and validated.

R n T ∞ ρ ∞ R-R n m -1 mm K kg/m 3 R01
The flow around CCF10 is studied at a Mach number M ∞ = 6.0 for a unit Reynolds number of Re = 11.5 × 10 6 . The domain extends from the wall to the free-stream and from the stagnation point to the end of the flare. An isothermal wall at T w = 300K is considered. These conditions relate to experiments performed in quiet and noisy wind tunnels [START_REF] Benitez | Instability Measurements on an Axisymmetric Separation Bubble at Mach 6[END_REF] and are summarised in Table 1 The computation of a fixed point uses the following steps : first, a preliminary flow is computed on a coarse structured mesh which is not shock-aligned. Second, once this coarse solution reaches a sufficient convergence level, the shock boundary is extracted and a fine shock-aligned grid is generated. The coarse solution is linearly interpolated on the fine grid and serves as the initial guess to a pseudo transient continuation technique [START_REF] Crivellini | An implicit matrixfree Discontinuous Galerkin solver for viscous and turbulent aerodynamic simulations[END_REF]. The solution vector update q q q n+1 = q q q n + δq q q n is given by the computation of the step δq as, The residuals R(q q q n ) are computed using a 7-th order FE-MUSCL scheme with a high order dissipation and Jameson-like shock capturing controlled by a shock sensor, resulting in low dissipation levels [START_REF] Sciacovelli | Assessment of a high-order shock-capturing central-difference scheme for hypersonic turbulent flow simulations[END_REF]. The jacobian is computed with the same numerical methods through algorithmic differentiation. Between 5 and 20 pseudo-Newton iterations are required to converge the fixed-point to a normalised residual value of ||R(q q q 0 )|| 2 ≤ 10 -14 .

I ∆t + Lq q q n δq n = -R(q q q n ). (9)
Baseflows q q q 0 are shown in Fig. 1 for cases R01 and R5. Observing the flow features for R n = 0.1mm & R n = 5mm reveals two main differences. First, starting from the nosetip, the R5 case shows a thicker entropy layer up to the cone-cylinder junction induced by the detached shock. Whereas for the R01 case, the entropy layer remains thin and is actually following closely the boundary layer height. The entropy layer effects are known to affect the growth of instabilities on hypersonic blunt bodies [START_REF] Stetson | Nosetip bluntness effects on cone frustum boundary layer transition in hypersonic flow[END_REF] and will be investigated in the coming sections. The second noticeable difference between the two cases is the length of the separation bubble, increasing the bluntness leads to a bigger separation length induced by the reduced edge Mach number M e and the thicker boundary layer.

The computed fixed point for R n = 0.1mm is validated against the previous study of [START_REF] Paredes | Boundary-Layer Instabilities over a Cone-Cylinder-Flare Model at Mach 6[END_REF] for the same conditions but different numerical methods. Fig. 2 shows a perfect agreement between the two solutions hence validating both computations of the hypersonic laminar baseflow.

For each fixed-points, the jacobian is computed and is used to probe the linear stability of this two flows.

GLOBAL STABILITY

This section aims at solving the eigenvalue problem introduced in Eq. 5 for successive discrete azimuthal wavenumbers m. For case R01 some modes of the spectrum are compared and validated against the previous findings of [START_REF] Paredes | Boundary-Layer Instabilities over a Cone-Cylinder-Flare Model at Mach 6[END_REF]. On the other hand the discussion on the global stability case R5 constitutes a new set of results.

Each spectrum is obtained from the jacobian using an iterative LU+Arnoldi algorithm as depicted in [START_REF] Poulain | BROAD-CAST: A high-order compressible CFD toolbox for stability and sensitivity using Algorithmic Differentiation[END_REF], the computations uses the PETSc/SLEPc routines along with the MUMPS library. This iterative solver is set to converge the 5 eigenvalues closest to a target τ = 1000.0 + 0.0i using a shift-invert strategy. As the unstable modes and spectrum computation are dependent on the baseflow discretisation, a spectrum convergence is carried out beforehand. The Fig. 3a displays this convergence in the wavenumbergrowth-rate frame. A reference computation from [START_REF] Paredes | Boundary-Layer Instabilities over a Cone-Cylinder-Flare Model at Mach 6[END_REF] is used to plot the modes branches. Note that since the domain is axisymmetric, only discrete values of m can physically be allowed to exist. For the spectrum convergence, three grids with increasing number of streamwise and vertical points are considered. All of them yield consistent results with the reference, therefore confirming that even a coarse grid resolution allows to properly converge these global structures. Although, as it will be required later for the forced analysis, a grid of N x × N y = 4030 × 550 will be used for the subsequent computations.

Aside from the validation analysis, comparing the Figs. 3a & 3b unveils significant differences between the global dynamics of the bubble between R n = 0.1mm and R n = 5mm. The spectrum of the blunt configuration reveals three branches of unstable modes. All the eigen values are coloured by their frequency to distinguish unsteady dynamics. Looking at the spectrum of Fig. 3b reveals two steady branches for m ∈ [START_REF] Chu | On the energy transfer to small disturbances in fluid flow (Part I)[END_REF][START_REF] Sidharth | Onset of three-dimensionality in supersonic flow over a slender double wedge[END_REF] and m ∈ [START_REF] Poulain | BROAD-CAST: A high-order compressible CFD toolbox for stability and sensitivity using Algorithmic Differentiation[END_REF]26]. This latter branch dominates the growth-rates of the spectrum. A third unsteady branch is visible for m ∈ [START_REF] Hartman | Nonlinear transition mechanism on a blunt cone at Mach 6: Oblique breakdown[END_REF][START_REF] Paredes | Boundary-Layer Instabilities over a Cone-Cylinder-Flare Model at Mach 6[END_REF], these modes have a frequency between f = 60Hz and f = 90Hz, suggesting quasi-steady dynamics in comparison to the flow characteristic time scales.

All the unstable modes arising in Fig. 3 appear to be structures modulating the bubble shape with wavenumbers m > 0. Therefore, the growth of these structures leads to a three-dimensionalisation of the baseflow. Similar mechanisms were found by [START_REF] Sidharth | Onset of three-dimensionality in supersonic flow over a slender double wedge[END_REF][START_REF] Cao | Transition to turbulence in hypersonic flow over a compression ramp due to intrinsic instability[END_REF][START_REF] Mathieu Lugrin | Multi-scale study of the transitional shock-wave boundary layer interaction in hypersonic flow[END_REF]. The eigenfunctions of the dominating modes for each unstable branch of Fig. 3b are plotted in Fig. 4. The mode at (f, m) = (0, 5) is representative of bubble breathing mode. Its structure is similar for R n = 0.1mm and R n = 5mm, suggesting that it comes from the same underlying physical mechanism in both cases. The eigenfunction for (f, m) = (81, 13) seem to be mainly amplified near the reattachment point, displaying a strong transverse motion at this position. Finally, the instability for (f, m) = (0, 19) also seem to be a reattachment mode and show a transverse motion very similar to the mode (f, m) = (60, 13). Further analysis have to be made with the eigenfunctions in order understand the physical origins of these similarities.

The global dynamics obtained from the linear stability analysis of Eq. 5 revealed the presence of multiple steady and unsteady structures susceptible to alter the laminar base flow and even lead to turbulence. In the non-linear regime, these global mode will saturate and may interact with the convective instabilities originating from the external disturbances, hence, the next section aims at characterising these convective instabilities.

OPTIMAL FORCING

The response of the baseflow to external disturbances is analysed using Eq. 8. Maps of the optimal gain in (ω, m) are computed for each case and allow to see the most amplified mechanism at a given frequency and wavenumber along the object for both nose radii.

In order to quantify the effects of the separation region, these maps are computed for both the full geometry and a subdomain restricted to the cone. Indeed, as the baseflow can be globally unstable on the cylinder-flare portion, an input-output analysis on the steady-homogenous fixed point q q q 0 computed through Eq. 9 could be ill-posed as the system bifurcated. Hence, solving Eq. 8 for the cone region, which is globally stable and independent of the bifurcated state due to the supersonic regime, allows for a consistent optimal response computation and gives a linear estimation of the wave content that might enter the bubble in a non-linear scenario.

Figure 5 regroups the gain maps for cases R01 and R5. First, a baseline analysis is made for the R01 case before showing the variations induced by the blunt nose of case R5. Figure 5a shows the optimal gains map for both the full geometry and the region restricted to the cone only. In summary, three main resonance regions appear : the dominating resonance peak is found near (f, m) = (75kHz, 20) and seems to be related to the family of first mode instabilities. Next, a resonance is visible for (f, m) = (225kHz, 10) and is experimentally verified to be related to the family of second-mode instabilities. Finally, a third resonance peak is found for steady waves at (f, m) = (0kHz, 70) identified as streamwise streaks.

Looking at the cone-only gain map highlights different dynamics, as the dominating peak is now the second mode at (f, m) = (250kHz, 10), which is consistent with previous literature results on canonical sharp cone boundary layers at M ∞ = 6.0. The first mode broadband peak remains present for (f, m) = ([50, 100]kHz, 20) but shows an amplifica- tion at least one order of magnitude smaller that the second mode. Finally, the main difference with a resolvent computation on the full geometry resides in the absence of a clear streaks amplification peak at the previously observed (f, m) = (0kHz, 70) region. This does not necessarily imply the total absence of streaks as these can still have a marginal amplification from non-normality on the cone (µ 0 ≈ 10 2 ). Additionally, the maximum value of m chosen for the gain map might not be high enough for capturing the higher wavenumber streaks that might appear close to the nosetip where the boundary layer height is smaller.

Comparing the previous results with the gain maps of the R5 case in Fig. 5b shows that the bluntness induces substantial differences in the baseflow response to an external forcing. Looking at the full object gain map, the peaks of the first and second mode instabilities observed for case R01 have vanished. The R5 baseflow displays a weak non-normal amplification for frequencies above f = 50kHz. The only observed peak is related to the steady streaks having a maximum gain : µ 0 = 1.4 × 10 3 at (f, m) = (0, 50). The gain map restricted to the blunt cone of case R5 shows the same absence of strongly amplified peaks for the frequencies and wavenumbers corresponding to the first and second modes found in Fig. 5a.

As such, these latter frequency-wavenumber ranges can no-longer be related to first and second mode waves as no clear resonance peak emerges. Instead, the gain map region for f ∈ [25, 300]kHz and m ∈ [0, 150] is most likely related to weak nonnormal amplification mechanisms. This hypothesis has to be further investigated in following studies by looking at the rank of the resolvent operator in this region. Nonetheless, the behaviour of the waves at (f, m) = (75kHz, 20) and (f, m) = (225kHz, 10) will also be studied for case R5 in order to have a comparison basis between the sharp and blunt dynamics.

These differences in the optimal gains maps be- tween the sharp and blunt nosetips can be related to the transition reversal phenomenon encountered for hypersonic boundary layers [START_REF] Stetson | Nosetip bluntness effects on cone frustum boundary layer transition in hypersonic flow[END_REF]. For the blunt nose, the entropy layer added to the decrease in edge Mach number induced by the strong shock might suppress the support of the unstable waves found for R01. It should be precised that the gain maps for case R5 contains some unresolved gains value (white regions), these arise from difficulties in solving Eq. 8 in presence of the strong detached shock. The discontinuity also induces a gain peak at (f, m) = (0, 0). This gain peak can not be directly interpreted as a convective instability and seem related to be related to an infinitesimal shock motion induced by the forcing of the momentum equations which reproduce the effect of infinitesimal variations in the Mach number.

Further insights can be gained by looking at the spatial growth along the object of the structures corresponding to the three peaks discussed above. In supplement to the global information provided by the gain maps, it helps to understand how the op- timal forcings and responses evolve around the different regions. Considering the mesh to be quasiwall normal, up to a small approximation error, the integrated energy of the modes is computed in the wall-normal grid direction η as,

E y = η q q q * W E q q qdη. (10)
Where, q q q and W E represents the resolvent response mode and its associated energy weighting matrix as stated in Sec. 2. A similar procedure holds for the optimal forcings as well.

Figure 6 displays the integrated energy of the optimal responses in thick lines and the associated forcing norm in dotted lines. For case R01 in Fig. 6a, it can be immediately noticed that the overall energy of the first and second modes remains two or three order of magnitude higher than the steady streaks at m = 70 along the whole geometry. These former two modes displays similar energy levels on the cone up to x ≈ 0.28m. From that point, the (100, 20) wave shows a steady amplification, whereas the (225, 10) wave exhibits an increased growth that leads to its Figure 8: Comparing optimal responses q q q and optimal forcings f f f at (f, m) = (100kHz, 20) for cases R01 (left) and (R5). Normalised contours of (ρv) ′ ; ( ) entropy layer ; ( ) boundary layer edge (h t /h t∞ = 0.995) ; ( ) detached shock.

dominance at the end of the cone. This increased growth of the second mode wave is explained by the dependency of this family of instabilities to the boundary layer height. For case R5, the same modes are plotted for unsteady waves as a comparison basis. These waves undergo a weak amplification consistent with the low gain observed in Fig. 5b. Further discussions on the differences of these non-normal mechanisms for R01 and R5 will be provided in a subsequent section.

A second point of interest is the expansion fan of the cone-cylinder junction after x = 0.4, it displays a stabilising effect for the three considered waves in case R01 and R5, especially for the streaks. Oppositely, the expansion fan is a favourable region for receptivity as it can be seen by the sharp and localised increase in the forcings norm. Next, the instabilities are convected through the separated flow. For this region, the curves of Fig. 6a are supplemented by the corresponding optimal response eigenfunctions of case R01 shown in Fig. 7. Focusing on Fig. 6a, along the cylinder, the first-mode and the streaks instabilities remain amplified, especially the streaks which shows a sudden growth rate increase at the onset of separation. On the contrary, the second mode waves get slightly damped, suggesting a stabilising effect of the separated region on these convective instabilities. Near the reattachment, all modes display an increase in energy, most probably linked with the strong re-compression, and keep amplifying along the flare. Especially, the first mode mode wave has an energetic content at the end of the flare similar to the second mode wave. This states that both mode have to be equally considered when studying a transition scenario on this geometry at similar conditions. The response modes shown in Fig. 7 indicate that the spatial support of the streaks is mainly located within the boundary layer and after the separation region. On the other hand both first and second mode waves show some fluctuations along and under the shear layer of the bubble around x = 0.52m.

Focusing on the forcings location along the full object, some comments can be made on the receptivity support of these instabilities. For case R01, the first and second mode forcings shown in Fig. 6a indicate that the region close to the nosetip is the most sensitive to the receptivity process. After x = 0.1m, the forcings energy steadily decreases, offering lesser projection support for external disturbances. A similar behaviour is observed for the streaks forcing vector, although, its energy remains quite constant up to the end of the cone, suggesting that the disturbances leading to the streaks on the flare are actually growing all along the geometry, a similar observation is made for the streaks of case R5.

In relation to the forcing spatial support, additional information in the difference of optimal gains between the sharp and blunt case can be gained by looking at the resolvent optimal eigenfunction around the nosetip. Figure 8 provides this comparison for the wave at (f, m) = (100kHz, 20). This frequency and eigenvalue corresponds to the peak of the first mode instability for case R01 and to a weak and broadband non-normal amplification mechanism for case R5. For both geometries, the entropy layer, the boundary layer and the detached shock-wave of the respective baseflows are plotted in order to illustrate the differences in spatial support of the non-normal mechanisms. Consistently with previous literature results, for case R01, the optimal forcings are located upstream on the cone, close to the nosetip and within the boundary layer. The associated response is amplifying downstream around the boundary layer edge. For the wave of R5 at the same frequency and wavenumber, different spatial supports are observed. The forcing is no longer restrained close to the nosetip and spread along the cone. Moreover, its most energetic region is mainly located within the entropy layer and does not attain the boundary layer. The associated response is located downstream on the cone but is also mainly amplifying along the entropy layer edge and does not seem to fluctuate within the boundary layer.

Such differences support the hypothesis of a different mechanism at the origin of the wave amplification on the blunt cone. Similar entropy layer modes were found by previous DNS studies for moderately blunt cones [START_REF] Hartman | Nonlinear transition mechanism on a blunt cone at Mach 6: Oblique breakdown[END_REF][START_REF] Paredes | Mechanism for frustum transition over blunt cones at hypersonic speeds[END_REF]. 9: Comparison of the R01 case second-mode pressure response at the wall with wind tunnel measurements from unsteady pressure sensors for the same conditions [START_REF] Paredes | Boundary-Layer Instabilities over a Cone-Cylinder-Flare Model at Mach 6[END_REF].

Finally, an assessment of the resolvent ability to predict hypersonic instabilities growth around the CCF10 object is made in Fig. 9 by comparing the optimal second-mode (f, m) = (225kHz, 10) pressure at the wall, to experimental results [START_REF] Paredes | Boundary-Layer Instabilities over a Cone-Cylinder-Flare Model at Mach 6[END_REF]. The computation is shown along with the spectra peak value at f ≈ 225kHz acquired during a wind tunnel test at the same free-stream conditions. These spectra were obtained with wall mounted PCB and Kulite pressure sensors placed at different stations at the end of the cone and along the flare. Since the computation is linear, an initial amplitude A 0 has to be set for the resolvent mode to have a physically meaningful amplitude, this value is calibrated to match the first PCB sensor energy at x = 0.387m, leading to A 0 = 1.2 × 10 -20 . The computation and the sensors wall pressure shows a close agreement in the reattachment region, the resolvent pressure amplitude function is matching the experimentally observed trends in pressure fluctuations. Such a result confirms the ability of the current linear framework to predict accurate hypersonic instabilities evolution in presence of complex laminar flow features.

CONCLUSION

The present study provides a comprehensive review of an ongoing numerical investigation aiming at uncovering the linear instabilities of the hypersonic laminar flow around the CCF10 geometry. This research works towards providing a deeper understanding of the transition to turbulence under experimental and flight conditions by mapping the various linear instabilities at play.

To achieve this objective, the numerical framework was first validated against previous results, followed by a mapping and comparison of the linear mechanisms on both sharp and blunt configurations. The validation procedure demonstrated the accuracy of the BROADCAST toolbox in providing fixed points as a starting point for the global stability investigation. Additionally, the obtained global modes for the sharp case were found to compare well with the literature's spectrum for the same conditions. Further analysis revealed differences in the laminar bubble dynamics for the sharp and blunt global spectra, with an increased number of unstable branches and the presence of unsteady bubble modes for the R5 case. The origin of these different dynamics is yet to be explained and may be investigated through the structural sensitivity of the baseflow.

The study also highlighted the differences in amplifying mechanisms between the sharp and blunt configurations. The sharp configuration displayed clear amplification peaks for the first and second Mack modes, whereas the blunt configuration only showed an amplification peak for steady streaks.

Investigating further these entropy layer modes allowed to observe their associated optimal forcings structures. Such forcing waves were scarcely discussed in the literature and provide a good starting point for explaining the origin of the observed entropy layer instabilities observed in the aforementioned DNS studies.

Finally, a direct comparison with experimental data obtained from the Purdue BAM6QT quiet wind tunnel confirmed the framework ability to accurately estimate convective instabilities evolution along the object.

Overall, this preliminary study validated the stability framework and provided insights into the linear dynamics around the Cone-Cylinder-Flare geometry. Future work will focus on explaining the physical mechanisms at the origin of the discrepancies observed between the sharp and blunt cone dynamics. This set of results will then help to better understand future non-linear simulations and experimental campaigns.

Figure 1 :

 1 Figure 1: Baseflow for the R n = 0.1mm (top) & R n = 5mm (bottom) CCF10 geometries. Greyscale : normalised numerical shadowgraphy with : entropy layer ; : recirculation bubble with separation and reattachment points (•) ; : line along which the profiles of Fig. 2 are extracted.

Figure 2 :

 2 Figure 2: Validation of the base-flow by comparison with a reference solution at the cylinder-flare junction [14].

Figure 3 :

 3 Figure 3: Eigen values growth-rates as a function of m for both nose radius. The left spectrum shows a grid convergence study. The right spectrum eigenvalues are coloured by their frequencies (Hz)

Figure 4 :

 4 Figure 4: Eigenfunctions of the leading global modes for R n = 5mm. Normalised contours of disturbances transverse momentum ; ( ) boundary layer ; ( ) separated region.

  (a) Rn = 0.1mm (b) Rn = 5mm

Figure 5 :

 5 Figure 5: Maps of optimal gain for the R01 case for the full geometry (left) and the cone only (right). Dominating peaks : (•). Light grey dots show the points actually computed for building the maps.

  m) = (0kHz,70) (f, m) = (100kHz,20) (f, m) = (225kHz,[START_REF] Mathieu Lugrin | Multi-scale study of the transitional shock-wave boundary layer interaction in hypersonic flow[END_REF]

Figure 6 :

 6 Figure6: Integrated energy of optimal responses q q q ( ) and forcings f f f ( ) for the different peaks observed in Figs. 5a & 5b. • : separation and reattachment points. In the legend from left to right : streaks, first mode, second mode, frequencies are in kHz. The object shape is depicted at the bottom for clarity.

Figure 7 :

 7 Figure 7: Optimal responses corresponding to the leading resolvent modes for R n = 0.1mm. Normalised contours of disturbances momentum, saturated at 10% of the normalised value ; ( ) boundary layer ; ( ) separated region.

Figure

  Figure9: Comparison of the R01 case second-mode pressure response at the wall with wind tunnel measurements from unsteady pressure sensors for the same conditions[START_REF] Paredes | Boundary-Layer Instabilities over a Cone-Cylinder-Flare Model at Mach 6[END_REF].


domain: spi phys

text:
Camille Fournis
email: camille.fournis@onera.fr

Ilias Petropoulos

Fulvio Sartor

Unsteady Far-Field Drag Analyses of Transonic Buffet over the NASA Common Research Model

This paper presents an analysis of the physics of transonic buffet by means of phenomenological far-field decomposition of drag. The analysis relies on the unsteady farfield drag decomposition developed by Toubin et al. ("Improved Unsteady Far-Field Drag Breakdown Method and Application to Complex Cases," AIAA Journal, Vol. 54, No. 6, 2016) so far applied to pitching airfoils and wings, vortex-shedding and two-dimensional transonic buffet. It provides a decomposition of drag into viscous, wave, lift-induced and acoustic components. First of all, the formulation of Toubin et al. is presented and a physical interpretation of the terms of the decomposition is given. Then, the main buffet characteristics are shown on the wing-body configuration of the NASA Common Research Model 1 (CRM) investigated in the frame of the 7 th AIAA Drag Prediction Workshop (DPW-7). Later on, the unsteady drag exerted on the aircraft is decomposed and the evolution of the various drag components over time is investigated.

INTRODUCTION

Given the inevitable economical and environmental challenges that have to be addressed, the aviation sector is ever more keen on reducing its operating costs and carbon footprint. One big aspect of this task is to minimize the fuel consumption of aircraft.

That is why accurate drag prediction has always been of paramount importance when designing an aircraft. Indeed, fuel consumption is a direct consequence of the amount of drag exerted on the aircraft, hence drag reduction implies lower fuel consumption. But before reducing drag, it is necessary to accurately measure and quantify it. To do so, multiple scientists have been de-veloping drag prediction methods throughout the twentieth century. Some of those methods rely on the farfield approach, which consists in considering the momentum balance in a control volume of fluid surrounding the aircraft. It differs from the classical near-field (n f ) approach, which mechanically decomposes drag into pressure and friction components on the aircraft skin:

D n f = D p + D f = S a (p -p ∞ ) n x dS - S a τ x • n dS (1)
Indeed, the far-field ( f f ) approach traditionally decomposes drag into phenomenological contributions: liftinduced drag D i , viscous drag D v and wave drag D w .

D f f = D i + D v + D w (2)
This breakdown is crucial for aircraft designers because it provides an insight on the physical sources of drag and helps to identify a potential for the improvement of the aircraft geometry. The lift-induced drag was first identified at the beginning of the twentieth century by Prandtl in his famous lifting-line theory [START_REF] Prandtl | Theory of lifting surfaces[END_REF] and its definition was amended in the seventies by Maskell [START_REF] Maskell | Progress towards a method for the measurement of the components of the drag of a wing of finite span[END_REF]. Methods for viscous and wave drag computation were then defined using thermodynamic approaches introduced by Oswatitsch [START_REF] Oswatitsch | Der luftwiderstand als integral des entropiestromes[END_REF] in 1945 and later by van der Vooren and Sloof [START_REF] Van Der Vooren | CFD-based drag prediction: state-of-the-art, theory, prospects[END_REF] in the nineties. In further developments, Destarac and van der Vooren [START_REF] Destarac | Far-field/near-field drag balance and applications of drag extraction in CFD[END_REF][START_REF] Destarac | Drag/thrust analysis of jet-propelled transonic transport aircraft; definition of physical drag components[END_REF] defined the lift-induced drag indirectly, by merely subtracting the irreversible drag (viscous plus wave drag) from the total far-field drag. In spite of the insights provided by those developments, these drag prediction and decomposition methods were only applicable to steady flows.

In the recent years, several formulations aimed at assessing and decomposing drag in unsteady flows have been developed [START_REF] Gariépy | Generalization of the far-field drag decomposition method to unsteady flows[END_REF][START_REF] Ostieri | Linear and nonlinear decomposition of aerodynamic force acting on an oscillating plate[END_REF][START_REF] Ostieri | Aerodynamic force and Lamb vector field in compressible unsteady flows[END_REF][START_REF] Toubin | Development and application of a new unsteady far-field drag decomposition method[END_REF][START_REF] Toubin | Improved unsteady far-field drag breakdown method and application to complex cases[END_REF]. Gariépy et al. [START_REF] Gariépy | Generalization of the far-field drag decomposition method to unsteady flows[END_REF] were the first to propose a generalization of the Destarac-van der Vooren approach [START_REF] Destarac | Far-field/near-field drag balance and applications of drag extraction in CFD[END_REF][START_REF] Destarac | Drag/thrust analysis of jet-propelled transonic transport aircraft; definition of physical drag components[END_REF] in unsteady flows. However, as shown by Toubin and Bailly [START_REF] Toubin | Development and application of a new unsteady far-field drag decomposition method[END_REF], their decomposition is very sensitive to the size of the control volumes used for the integration of the various drag components: an extension of the integration volumes in the wake creates an unphysical phase shift on the time evolution of the farfield drag components. Indeed, it is necessary that the phase of all far-field drag components be invariant to the size of the integration domain, in order to guarantee an objective physical unsteady drag breakdown. Moreover, it ensures that the farfield drag integrated in the flowfield surrounding the aircraft (hence comprising all the unsteady phenomena occuring therein) is consistent at all times with the near-field drag integrated on the skin of the aircraft. This particular aspect has been one of the major successes of Toubin's doctoral thesis [START_REF] Toubin | Prediction and phenomenological breakdown of drag for unsteady flows[END_REF]. Toubin et al. [START_REF] Toubin | Prediction and phenomenological breakdown of drag for unsteady flows[END_REF][START_REF] Toubin | Development and application of a new unsteady far-field drag decomposition method[END_REF][START_REF] Toubin | Improved unsteady far-field drag breakdown method and application to complex cases[END_REF] developed an unsteady drag breakdown based on the Destarac-van der Vooren formulation [START_REF] Destarac | Drag/thrust analysis of jet-propelled transonic transport aircraft; definition of physical drag components[END_REF]. They validated this method on a number of academic cases and demonstrated its robustness with respect to the integration volume extension in the wake.

The interest in unsteady phenomenological drag breakdown has so far been limited to relatively academic cases. Recent drag-related studies have however begun to emerge on the more applied NASA CRM aircraft in transonic buffet conditions. Steady investigations of the CRM have been the main focus of recent AIAA Drag Prediction Workshops (DPWs) [START_REF] Coder | Contributions to the sixth drag prediction workshop using structured, overset grid methods[END_REF][START_REF] Tinoco | Summary data from the sixth AIAA CFD drag prediction workshop: CRM cases[END_REF]. As an extension of these studies, the latest edition (DPW-7) featured a task dedicated to "Beyond RANS" analyses, aiming at investigating flight conditions at greater Reynolds numbers and higher angles of attack. ONERA took part in this DPW-7 task on the basis of unsteady Reynolds-averaged Navier-Stokes (URANS) computations. In particular, when considering a wing in cruise conditions and increasing the angle of attack, periodic shockwave motions known as transonic buffet occur, and the flow is unsteady [START_REF] Lee | Self-sustained shock oscillations on airfoils at transonic speeds[END_REF]. The unsteadiness predominantly consists of a periodic selfsustained motion of the shock as well as unsteady separation of the boundary layer, consequences in the flowfield that will be analyzed in the present work.

The aim of this study is to investigate whether an unsteady drag decomposition can establish links between the evolution of drag components and unsteady phenomena occurring within the flow field. This work has been carried out by using the far-field decomposition of Toubin et al. [START_REF] Toubin | Improved unsteady far-field drag breakdown method and application to complex cases[END_REF] to analyze URANS computations of the DPW-7 CRM configuration in transonic buffet regime. The study is structured as follows. The formulation employed is first presented in Section 1, followed by Section 2 which gives an overview of ONERA computational studies of buffet prediction in the frame of the DPW-7. Finally, Section 3 presents investigations on the unsteady drag breakdown of the CRM case.

PRESENTATION OF THE

UNSTEADY FAR-FIELD DRAG DECOMPOSITION OF TOUBIN ET AL.

Original formulation

To overcome the phase shift issue raised by the formulation of Gariépy et al. [START_REF] Gariépy | Generalization of the far-field drag decomposition method to unsteady flows[END_REF], Toubin and Bailly [START_REF] Toubin | Development and application of a new unsteady far-field drag decomposition method[END_REF] adopted a different approach, based on a drag balance in streamtubes, in order to account for the time delay and propagation of flow perturbations created at the source of drag by shockwave motion, vortex-shedding or a pitching airfoil or wing for instance. By doing so, they could eliminate the phase shift entailed by an extension of the control volumes in the wake. Their formulation is presented below:

D f f = D w + D v + D ui + D m (3)
with:

D w = - S w ρ (U irr -U ∞ ) (q • n) dS - V w ∂ ρ (U -U ∞ ) ∂t dV - V wd ∂ ρ (U -U irr ) ∂t + 1 U irr ∂ p ∂t dV (4)
D v = S v (-ρ (U irr -U ∞ ) q + τ x ) • n dS - V v ∂ ρ (U -U ∞ ) ∂t dV (5)
D ui = S e (-ρ (U -U irr ) q -(p -p ∞ ) e x ) • n dS - V c ∂ ρ (U -U irr ) ∂t + 1 U irr ∂ p ∂t dV (6)
D m = - S a ρ (U -U ∞ ) (q • n) dS (7)
where D ui is the unsteady lift-induced drag and D m is the drag generated by solid body motion, henceforth referred to as motion drag. The volumes and surfaces are illustrated in Fig. 1 and S a is the skin surface of the aircraft.

Here q is the velocity vector, τ x = τ • e x and U irr is the irreversible axial velocity defect [START_REF] Destarac | Far-field/near-field drag balance and applications of drag extraction in CFD[END_REF][START_REF] Destarac | Drag/thrust analysis of jet-propelled transonic transport aircraft; definition of physical drag components[END_REF] defined by:

U irr = U ∞ 1 + 2∆H U 2 ∞ - 2 (γ -1) M 2 ∞ e ∆s/C p -1 (8)
and H = C p T + q 2 /2 is the specific stagnation enthalpy, ∆H = H -H ∞ is the specific stagnation enthalpy variation, γ is the ratio of specific heats, C p is the specific heat at constant pressure, and ∆s = ss ∞ is the specific entropy variation. Later on, Toubin et al. [START_REF] Toubin | Improved unsteady far-field drag breakdown method and application to complex cases[END_REF] further decomposed D ui into an acoustic propagation drag D pa and the lift-induced drag D i as follows:

D ui = D i + D pa (9)
Figure 1: Control volumes used in the integration of the various drag components (from [START_REF] Toubin | Improved unsteady far-field drag breakdown method and application to complex cases[END_REF])

with:

D i = S e (-ρ (U -U irr ) q -(p -p ∞ ) e x ) • n dS + V c ρ (U -U irr ) 1 a + 1 U irr ∂U ∂t dV ( 10
)
D pa = - V c ρ (U -U irr ) a ∂ R + ∂t dV ( 11
)
where a is the speed of sound and

R + = U + 2a/ (γ -1)
is the Riemann invariant on characteristic curves of the (x,t) plane defined by dx/dt = U + a. In summary, the unsteady far-field drag can be decomposed into five contributions:

D f f = D w + D v + D i + D pa + D m (12)
First of all, it is important to note that each drag component includes an unsteady contribution in the formulation of Toubin et al. [START_REF] Toubin | Improved unsteady far-field drag breakdown method and application to complex cases[END_REF], whereas in the formulation of Gariépy et al. [START_REF] Gariépy | Generalization of the far-field drag decomposition method to unsteady flows[END_REF], unsteady effects are grouped into a single unsteady drag component. Secondly, in steady flows, all the time derivatives vanish and D w , D v and D i recover their steady expressions given in [START_REF] Destarac | Far-field/near-field drag balance and applications of drag extraction in CFD[END_REF][START_REF] Destarac | Drag/thrust analysis of jet-propelled transonic transport aircraft; definition of physical drag components[END_REF], while D pa is eliminated. Moreover, D m also disappears in steady regime given the impermeable wall condition q • n = 0 which holds for both viscous and inviscid flows.

Single-vector formulation

In the frame of this work, a different but mathematically equivalent formulation has been derived in order to establish an unsteady counterpart of the single-vector version of the steady formulation of Destarac [START_REF] Destarac | Far-field/near-field drag balance and applications of drag extraction in CFD[END_REF]. The derivations rely on the following property in unsteady regime:

∇ • f = ∂ ρ (U -U ∞ ) ∂t (13)
with:

f = -ρ (U -U ∞ ) q -(p -p ∞ ) e x + τ x ( 14
)
By merely decomposing f as:

f = f * vw + f * i ( 15
)
with:

f * vw = -ρ (U irr -U ∞ ) q + τ x (16) f * i = -ρ (U -U irr ) q -(p -p ∞ ) e x ( 17
)
and after some algebraic manipulations, it is possible to re-express D w , D v , D ui and D i using only the f * i vector:

D w = - V wd ∂ ρ (U -U irr ) ∂t + 1 U irr ∂ p ∂t dV - S w f * i • n dS ( 18
)
D v = - S v f * i • n dS -D m + D n f ( 19
)
D ui = - V c ∂ ρ (U -U irr ) ∂t + 1 U irr ∂ p ∂t dV + S e f * i • n dS ( 20
)
D i = V c ρ (U -U irr ) 1 a + 1 U irr ∂U ∂t dV + S e f * i • n dS (21)
The following far-field drag results have all been obtained by post-processing the CFD solutions with the latter single-vector formulation.

ANALYSIS OF THE MAIN BUFFET CHARACTERISTICS ON THE NASA CRM

This section is devoted to the presentation and analysis of the main features of the buffet phenomenon predicted by URANS computations. They were performed using the ONERA-SAFRAN elsA solver [START_REF] Cambier | The ONERA elsA CFD software: input from research and feedback from industry[END_REF], in which the turbulence was modeled using the Spalart-Allmaras model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF], and the space discretization was done using the Jameson-Schmidt-Turkel numerical scheme for the convective flux [START_REF] Jameson | Numerical solution of the Euler equations by finite volume methods using Runge Kutta time stepping schemes[END_REF]. The time algorithm chosen is Gear, providing a global time stepping. Those time-accurate simulations correspond to the work carried out by ONERA [START_REF] Hue | Steady and unsteady computations of the CRM aircraft at different Reynolds numbers[END_REF] as part of the "Beyond RANS" optional test-case of the DPW-7. The participants were invited to use alternative methods to analyze a flow at M ∞ = 0.85, with high angle of attack (α = 4 • and 4.25 • ) and possibly buffet phenomenon for a Reynolds number of 20 million over the NASA CRM (the geometry is shown in Fig. 2). ONERA opted for URANS simulations since this approach can provide an accurate description of the unsteadiness of the shockwave while keeping the cost of the simulation reasonable. Contrary to high-fidelity methods such as DNS or ZDES, URANS is able to provide useful information about the unsteady dynamics of the flow in buffet conditions in a short amount of time and with limited resources. However, the turbulence remains completely modeled as for the steady-state RANS simulations and only the large-scale unsteadiness can be reproduced. The choice of the time step is the first question that has been addressed. A convergence study was carried out, starting from the converged RANS solution, using different time steps. At the beginning, the Tiny grid level of the DPW-7 (5M points) was used in order to have a quick estimation of the buffet phenomenon, localize the position of the unsteadiness on the wing, its intensity and, most importantly, the buffet period. While the y + = 0.4 discretization is sufficiently fine and the buffet phenomenon could be reproduced with shockwave motions observed on the wing, the spanwise spatial discretization in the buf-fet region is too coarse to validate the result when compared to the Coarse grid (17M points) and the Medium grid (41M points), as depicted in Figs. 3 and4. According to the data obtained with the Tiny grid, the buffet phenomenon has a frequency of 5.90 Hz (see Fig. 5). For this case, a time step of ∆t = 10 -4 s, corresponding to 1700 steps per buffet period, was sufficient to reproduce the phenomenon. Thus, it was concluded that the Tiny grid can reproduce the main features of buffet, despite poor spanwise discretization, but the unsteady behavior of the shockwave might not be correctly captured. Hence the Medium grid has been considered in the following studies. For this grid level, the time step is ∆t = 5 × 10 -5 s, corresponding to roughly 1600 steps per period, with each step converged during 8 inner iterations.

The three images in Fig. 6

INVESTIGATION AND EVOLUTION OF FAR-FIELD DRAG COMPONENTS OVER A PERIOD

The present section is dedicated to a deeper phenomenological analysis of the unsteadiness entailed by the buffet and its impact on the level of viscous, wave and liftinduced drag over time. First of all, it is necessary to investigate the effect of the size of the integration domain on the various drag components. To do so, the viscous, shock and induced volumes are extended in the wake of the aircraft wing, between x W = 50 m and x W = 60 m (see Fig. 7). In Fig. 8 and Tab. 1 (top rows), each value is obtained sition allowing to reduce as much as possible the effect of numerical dissipation. On the contrary, the averaged unsteady wave drag coefficient C D w remains constant at a rather significant level (around 90 counts), while the time-averaged acoustic propagation drag C D pa has a practically negligible contribution to overall drag. It is also interesting to note that the steady drag decomposition of the time-averaged solution is very close to the time-averaged unsteady drag decomposition obtained at x W = 50 m: in the present case, it is because the unsteadiness entailed by the buffet phenomenon is very mild. Nevertheless, Toubin [START_REF] Toubin | Prediction and phenomenological breakdown of drag for unsteady flows[END_REF] clearly demonstrated that values of steady drag decompositions applied to time-averaged flow solutions can be strikingly imprecise with respect to averaged values of unsteady drag decompositions.

At this stage, it should be pointed out that, in Fig. 7, there are some cells in the wake which are located in between the viscous and shock volumes, but do not belong to either of them. Yet, the wake of the shockwave is to be gradually entrained into the wake of the boundary layer, therefore those untagged cells are not expected. In fact, the extensions of the viscous and shock volumes are defined here using a physical sensor based on the level of viscous stresses and the vorticity modulus, but the latter is dissipated by the rapid coarsening of the grid. Hence, the integration volume issue might be solved by reducing the vorticity threshold used in this criterion. This would lead to the same total far-field drag, but to an eventually slightly different drag decomposition. Besides, some cells below the viscous volume are detected as part of the extended shock volume. There is no shockwave on the pressure side of the wing, therefore this feature is not expected either. Here, increasing the entropy threshold would likely remove those cells from this volume, although their contribution to the overall wave drag value is expected to be small. An important following step is to verify that the evolution of the far-field drag versus time is synchronized with that of the near-field drag at each time step. In other terms, it consists in checking that the far-field and nearfield drag forces are in phase. For this purpose, the evolution of far-field drag versus time has been plotted for all x W positions and compared to the evolution of near-field drag (see Fig. 9). It is clear that far-field and near-field drag are in phase for all x W positions, which confirms the consistency of the far-field drag computed by the formulation of Toubin et al. [START_REF] Toubin | Improved unsteady far-field drag breakdown method and application to complex cases[END_REF]. Regarding viscous and wave drag, the sensitivity of their phase with respect to x W is plotted in Figs. 10 and11: for the viscous drag a slight phase shift is observed (perhaps linked to the presence of cells not belonging to either the viscous or shock volume), while for the wave drag all curves are in phase. However, the C D v curves corresponding to the two most downstream positions x W = 58 m and x W = 60 m show almost no phase shift. In Fig. 12, the phase shift of the unsteady induced drag D ui seems slightly more pronounced than for the viscous drag, although it remains reasonably small. On the contrary, the phase shift becomes even more visible for the lift-induced drag D i (see Fig. 13). Actually, D i is computed from D ui after having distinguished the acoustic effects materialized by R + and involved in D pa . This phase shift in D i is indeed directly caused by the great sensitivity of the phase of D pa with respect to x W (see Fig. 14): although its variation is of a very low amplitude and therefore has a minimal contribution to total drag, it still has a direct impact on the phase of its counterpart D i . Similarly to the present case, Toubin [START_REF] Toubin | Prediction and phenomenological breakdown of drag for unsteady flows[END_REF] also observed such a phase shift on D i for the two-dimensional buffet case over the OAT15A airfoil but noticed that the phase shift tends to zero when integrating even further in the wake. The author explained that this is caused by a drastic decrease of the contribution of D pa , an effect which is clearly depicted in Fig. 14 for the present case. A gradual decrease in the magnitude of the dephasing D pa contribution would thus lead to a reciprocal decrease in the phase shift of D i (visible between the x w = 58 m and x w = 60 m curves in Fig. 13). Despite this observation, the choice of a sufficiently downstream x W position to force the invariance of the phase of D i is not satisfactory (at least in the present case), due to the influence of numerical dissipation in the wake in liftinduced drag prediction (cf. above disussion). Hence, the x W = 50 m position was chosen for the following analyses in order to ensure, as much as possible, that the computed viscous drag represents the outcome of viscous phenomena occurring on the wing, rather than that of numerical dissipation in the wing wake region.

x W (in m) C D i C D v C D w C D pa C D sp
As shown in Fig. 15, the maximum in viscous drag corresponds to a minimum in wave drag, which is consistent with Toubin's results in the case of two-dimensional buffet [START_REF] Toubin | Prediction and phenomenological breakdown of drag for unsteady flows[END_REF]. In a two-dimensional buffet case, the wave drag is minimum when the shockwave strength is the weakest, i.e. when the stall cell created by the boundary layer separation downstream of the shockwave foot pushes the shockwave itself at its most upstream buffeting position. This corresponds to the instant when the stall cell reaches its maximum size, hence producing maximum viscous pressure drag D vp (the part of the viscous drag that is not caused by skin friction). The interpretation is however more complex in the current three-dimensional buffet case, because the shockwave oscillates both in the chordwise and the spanwise directions, behaving similarly to a wave propagating towards the wingtip. This is better illustrated in Figs. 16 and17. At t = t 2 , the shockwave at the y = 20 m section (see Fig. 2) reaches its downstream limit while the shockwave at the y = 22.5 m section reaches its upstream limit. The maximum wave drag at t = t 1 could be explained by the fact that the shockwave at the y = 22.5 m section is at its most downstream position and overcompensates the contribution of the y = 20 m section, where the shockwave is at its most upstream position.

Figs. 18, 19 and 20 respectively display approximate spanwise distributions of the wave and viscous pressure drag coefficients C D w and C D vp computed by the ONERA FFDπ software at t = t 1 , t ≈ (t 1 + t 2 ) /2 and t = t 2 , comprising a fine discretization in the spanwise direction. The symbols correspond to the time-average of the spanwise distributions and is shown in order to highlight the deviation of the unsteady distribution. A contour of negative streamwise skin friction coefficient C f x on the wing suction side is also added on each figure in order to visualize the temporal evolution of the size of the separation region caused by the shockwave, allowing to corre-late the unsteady flow field variation on the wing to phenomenological drag sources. On average, the presence of the separation region entails a drastic increase in the viscous pressure drag between y = 10 m and y = 22.5 m. This also coincides with the region of the wing where the wave drag generation is highest. Indeed, the shockwave is strongest around y = 12 m, where the flow starts separating downstream. An investigation of the same figures illustrates the convection of the stall cell towards the tip of the wing by the displacement of a region of negative streamwise skin friction. This region matches the displacement of the peak in the spanwise distribution of C D vp . This peak in C D vp is very closely preceded by a dip in the spanwise distribution of C D w . The imperfect synchronization between the peak in dC D vp /dy and the dip in dC D w /dy is likely related to the very slight phase shift between the maximum reached by C D vp and the minimum reached by C D w illustrated in Fig. 15.

CONCLUSION

In the present study, the unsteady far-field drag decomposition method developed by Toubin et al. [START_REF] Toubin | Improved unsteady far-field drag breakdown method and application to complex cases[END_REF] has been adapted to a single-vector formulation and applied to the three-dimensional transonic buffet phenomenon over the NASA Common Research Model. To the authors' knowledge, this constitutes the first attempt at decomposing the drag in three-dimensional buffet conditions using a farfield technique.

The flight conditions investigated in this study enabled to effectively observe the onset of a very mild buffet phenomenon where the unsteadiness was confined to a region of limited size on the wing suction side. The Strouhal number of the observed buffet phenomenon appeared to be in line with previous results from the literature. In spite of this unsteadiness being very mild, the unsteady far-field drag decomposition method could capture the oscillations in the wave, viscous and lift-induced drag exerted on the aircraft. It was observed that the wave drag was maximum when the viscous drag was minimum and vice versa, as it was previously shown by Toubin et al. [START_REF] Toubin | Improved unsteady far-field drag breakdown method and application to complex cases[END_REF] in two space dimensions. Sensitivity analyses of the phase shift with respect to the wake extension of the integration volumes were also proposed, suggesting that the phase of the wave drag and the viscous drag is almost insensitive to this parameter, whereas the phase of the lift-induced drag showed a more pronounced sensitivity to the investigated integration volume wake boundaries, as these limits remained rather close to the wing.

In the end, some local analyses were carried out with the aim to link the unsteady variation of drag components to flow field variations. It was observed that the convection of the stall cell and the oscillation of the shockwave on the wing suction side are well depicted on the intensity of unsteady spanwise generation of wave and viscous pressure drag. The stall cell is characterized by a bump in the viscous pressure drag distribution and a dip in the wave drag distribution, whereas the shock oscillation is merely translated into oscillations on both distributions.

In the future, further investigations will be conducted with the aim to perform grid convergence studies. Then, potential improvements to physical sensors can be made (essentially the threshold robustness). Finally, additionnal simulations (URANS and/or ZDES) and drag analyses on a stronger buffet case can be carried out.

Figure 2 :

 2 Figure 2: Geometry of the NASA CRM.

Figure 3 :

 3 Figure 3: Pressure contour on the wing (Tiny grid).

  present a comparison of the pressure fluctuations obtained at three angles of attack. Here, the color-map has a logarithmic scale, allowing for a broader comparison of pressure fluctuations. On the left-hand-side, very low values of pressure fluctuations indicate that the flow is steady at C L = 0.58 (α = 2.755 • ). Very small fluctuations can be observed on the shockwave foot, but without any link to a buffet phenomenon. In the center, the pressure fluctuations at α = 4 • show that the flow is in buffet conditions. However, contrary to what is observed on the right-hand-side, where α = 4.25 • , some unsteadiness can be observed in the inner part of the wing, where the so-called "λ shockwave foot" is visible because of the interaction between the wing and the fuselage. The frequencies of the buffet phenomenon observed when α = 4 • and α = 4.25 • are 11.9 Hz and 12.2 Hz, respectively. These values correspond to a Strouhal number St = 0.464, in fair agreement with what was found by Paladini et al. [13] on a similar configuration, where the Strouhal number based on the mean aerodynamic chord for the FLIRET geometry (a half wing-body configuration) was St = 0.48.

  The results have been obtained by post-processing the CFD solutions at M ∞ = 0.85, α = 4.25 • , Re = 20 million of the Medium grid with the ONERA FFD72 and FFDπ far-field drag analysis codes.

Figure 4 :Figure 5 :

 45 Figure 4: Pressure fluctuations (linear colormap) for different grids.

Figure 6 :

 6 Figure 6: Pressure fluctuations (logarithmic colormap) for different wing loadings at the Medium grid level.

Figure 7 :

 7 Figure 7: Extension of the viscous and shock volumes for x W = 50 m (in red) and x W = 60 m (in orange) respectively.

Figure 8 :

 8 Figure 8: Relative evolution of each drag component averaged in time by the extension of the integration volumes in the wake.

Figure 9 :

 9 Figure 9: Temporal evolution of the far-field drag versus the near-field drag for various x W positions.

Figure 10 :

 10 Figure 10: Temporal evolution of the viscous drag normalized by the time-averaged value for various x W positions.

Figure 11 :

 11 Figure 11: Temporal evolution of the wave drag normalized by the time-averaged value for various x W positions.

Figure 12 :

 12 Figure 12: Temporal evolution of the unsteady induced drag normalized by the time-averaged value for various x W positions.

Figure 13 :

 13 Figure 13: Temporal evolution of the lift-induced drag for various x W positions.

Figure 14 :

 14 Figure 14: Temporal evolution of the acoustic propagation drag for various x W positions.

Figure 15 :

 15 Figure 15: Temporal evolution of C D v and C D w (normalized by their average in time) at x W = 50 m.

Figure 16 :

 16 Figure 16: Pressure coefficient at spanwise station y = 20 m for t = t 1 and t = t 2 .

Figure 17 :

 17 Figure 17: Pressure coefficient at spanwise station y = 22.5 m for t = t 1 and t = t 2 .

Figure 18 :

 18 Figure 18: Spanwise distributions of C D vp and C D w for t = t 1 .

Figure 19 :

 19 Figure 19: Spanwise distributions of C D vp and C D w for t ≈ (t 1 + t 2 ) /2.

Figure 20 :

 20 Figure 20: Spanwise distributions of C D vp and C D w for t = t 2 .

Table 1 :

 1

		Time-averaged drag -Unsteady formulation
	50	200.57 216.78 89.28 0.00 -0.23
	52	199.89 218.01 89.27 0.00 -0.77
	54	198.86 219.65 89.26 0.01 -1.38
	56	197.54 221.28 89.26 0.02 -1.69
	58	196.30 222.59 89.26 0.02 -1.76
	60	195.35 223.66 89.27 0.01 -1.90
		Time-averaged solution -Steady formulation
	-	201.40 215.96 89.02	-	0.03

Time-averaged value of each drag component (in counts, i.e. 10 -4 ) for x W ∈ [50 m, 60 m] and drag breakdown of the time-averaged solution.

CRM website: http://commonresearchmodel.larc.nasa.gov/

ACKNOWLEDGEMENTS

The studies presented have used the elsA software, whose development is partially funded by its two coowners: ONERA and SAFRAN.

The authors are grateful to Didier Bailly, ONERA, for the fruitful discussions on the unsteady far-field drag decomposition and to David Hue, ONERA, for his suggestions and for providing the steady flow solutions on the NASA CRM, from which the unsteady computations were restarted.