domain: spi chim info math phys

text:
Frédéric Plante

Julien Dandois

Samir Beneddine

Éric Laurendeau

Denis Sipp

Link between subsonic stall and transonic buffet on swept and unswept wings: from global stability analysis to nonlinear dynamics

This paper examines the three-dimensional cellular patterns appearing on wings in subsonic stall and transonic buffet conditions. Unsteady Reynolds-Averaged Navier-Stokes simulations are carried out for three-dimensional infinite swept configurations closed by periodic boundary conditions in the spanwise direction. In both flow conditions the occurrence of stall/buffet cells is observed, as well as their convection at a speed proportional to the sweep angle. In transonic buffet conditions, this phenomenon is superimposed to the well-documented two-dimensional buffet instability. These results indicate that the discrepancies between two-dimensional and three-dimensional buffet are caused by the occurrence of buffet cells and that this phenomenon is similar to the one observed at low speed. These phenomena are then studied using global linear stability analysis with the assumption of a periodic flow in the spanwise direction. From these analyses a mode coherent with the two-dimensional buffet is obtained, as well as a mode coherent with two-dimensional vortex shedding in stall conditions. In addition, in both flow conditions an unstable mode reminiscent of stall/buffet cells is observed.

Introduction

Transonic buffet and subsonic stall are two boundaries of the flight envelope of civil aircraft. Transonic buffet is an oscillation of the shock wave position over a wing caused by an interaction between a separated boundary layer and the shock wave. This oscillation induces variation of the aerodynamic forces which can be detrimental to aircraft handling. A thorough review of the state of knowledge about transonic buffet is presented by [START_REF] Giannelis | A review of recent developments in the understanding of transonic shock buffet[END_REF]. On the other hand, subsonic stall occurs when flow separation causes a reduction of the lift forces as the wing angle of attack is increased. These phenomena involve complex physics, part of which is the occurrence of three-dimensional flow features which are named buffet cells and stall cells respectively. The stall cells are identified as a spanwise variation of the flow separation, while the buffet cells are a spanwise variation of the shock wave position, and thus of the shock induced flow separation.

Stall cells were observed experimentally [START_REF] Moss | Two-Dimensional Low-Speed Tunnel Tests on the NACA 0012 section Including Measurements Made During Pitching Oscillations at the Stall[END_REF]; [START_REF] Gregory | Progress Report on Observations of Three-Dimensional Flow Patterns obtained during Stall Development on and on the Problem of Measuring Two-Dimensional Characteristics[END_REF]) since the 1970's. Most of these studies were carried out for wing which spans the entire wind tunnel width, a setup which is often considered as "two-dimensional". However, the flow becomes strongly three-dimensional when stall cells are present. Early studies suggested that the phenomenon might be caused by the side walls. However, the experimental results of [START_REF] Winkelmann | Flowfield Model for a Rectangular Planform Wing beyond Stall[END_REF] contradict this conclusion since the stall cells are observed over the entire span of wings of aspect ratios up to twelve and with free tips. They reported that the number of cells is proportional to the aspect ratio. [START_REF] Schewe | Reynolds-number effects in flow around more-or-less bluff bodies[END_REF] also reports the effect of the aspect ratio on the number of stall cells. Experiments by [START_REF] Broeren | Spanwise Variation in the Unsteady Stalling Flowfields of Two-Dimensional Airfoil Models[END_REF] of aerofoil with trailing edge, leading edge and thin aerofoil type of stall suggest that trailing edge separation is a necessary condition for the occurrence of stall cells. [START_REF] Yon | Study of the Unsteady Flow Features on a Stalled Wing[END_REF] obtained stall cells over a narrow range of angles of attack and discussed the unsteady nature of the stall cells. More recently, Dell'Orso & Amitay (2018) carried out a parametric study of the effect of the angle of attack (α) and the Reynolds number Re = ρ ∞ U ∞ c µ ∞ . They identified eight types of flow topologies from a full span separation (2-D) to the presence of one or two stall cells. From these results, the stall cells are observed above a given critical Reynolds number and angle of attack. They also observed conditions for which the topology oscillates in time between two types of flow topology. This might explain the fact that steady and unsteady stall cells are reported in the literature. Hence, the stall cells seem to be very sensitive to the flow conditions and are observed in the presence of a trailing edge separation. The stall cells can be steady or unsteady depending on the experiments. However, the measurement methods might also have a time-averaging effect and impact the visualisation of an unsteady behaviour.

Stall cells were computed by [START_REF] Bertagnolio | New Insight Into the Flow Around a Wind Turbine Airfoil Section[END_REF] for a wind turbine aerofoil with Reynolds-Averaged Navier-Stokes (RANS) and Unsteady Reynolds-Averaged Navier-Stokes (URANS) simulations. [START_REF] Manni | Numerical study of airfoil stall cells using a very wide computational domain[END_REF] studied the stall cells over a wing with an aspect ratio equal to ten at a Reynolds number of one million with Unsteady RANS (URANS) and Delayed Detached Eddy Simulation. Their URANS results showed stall cells with spanwise spacing of 1.4 to 1.8 chord length over a narrow range of angles of attack. [START_REF] Liu | Numerical analysis on the oscillation of stall cells over a NACA 0012 aerofoil[END_REF] studied the unsteady behaviour of stall cells over a NACA0012 at a Reynolds number of 1.35 × 10 5 and 1 × 10 6 . Previous work by the authors (Plante et al. (2019b)), based on URANS simulations of a NACA4412 aerofoil at Reynolds 350 000, showed the convection of the stall cells in the spanwise direction when the wing is swept. To the authors' knowledge this is the first analysis of the stall cells phenomenon on swept wings.

Multiple models have been proposed for the origin of the phenomenon. Analyses based on the lifting line theory by [START_REF] Spalart | Prediction of Lift Cells for Stalling Wings by Lifting-Line Theory[END_REF] and [START_REF] Gross | Criterion for Spanwise Spacing of Stall Cells[END_REF] suggest that a negative slope of the lift versus the angle of attack is necessary to observe the stall cells. Analyses with lifting surface models (Vortex Lattice Method) coupled with RANS data by [START_REF] Gallay | Nonlinear Generalized Lifting-Line Coupling Algorithms for Pre/Poststall Flows[END_REF] and [START_REF] Paul | Iteration schemes for rapid post-stall aerodynamic prediction of wings using a decambering approach[END_REF] exhibits similar lift distributions. [START_REF] Kitsios | Biglobal stability analysis in curvilinear coordinates of massively separated lifting bodies[END_REF] studied an ellipse and a NACA0015 aerofoil in stall conditions at a Reynolds number of 200 with global linear stability analyses. They found a non-oscillatory unstable mode with a spanwise wavenumber β = 2π/L z = 1.0. [START_REF] Rodríguez | On the birth of stall cells on airfoils[END_REF] linked this non-oscillating unstable global mode to the stall cells phenomenon. However, the existence of this mode has then been questioned by the later work of [START_REF] He | Linear instability of low Reynolds number massively separated flow around three NACA airfoils[END_REF], who were unable to recover the unstable global modes, and this with two different numerical approaches. [START_REF] Zhang | Biglobal linear stability analysis on low-re flow past an airfoil at high angle of attack[END_REF] carried out similar stability analyses on a NACA0012 aerofoil at Reynolds numbers in a range from 400 to 1000. Besides the classical oscillatory modes, they identified an unstable non-oscillatory threedimensional (β = 2 and 4) mode on an aerofoil in laminar regime but for slightly higher Reynolds numbers (between 800 and 1000). They associated this mode to a centrifugal instability of the recirculation bubble, present when the backflow velocity is larger than 20% of the inflow velocity. In these articles, flows with a low Reynolds number of a few hundreds to a thousand were studied. As such, the literature lacks a proper relation between stall cells and a global instability for high Reynolds turbulent flows. The present paper has the purpose of providing this link using RANS models. This is done by carrying out three-dimensional URANS simulations and global stability analyses.

Similar flow instabilities are reported for separation bubbles at low Reynolds numbers on various geometries. Separation bubbles on flat plates were analysed by [START_REF] Theofilis | On the origins of unsteadiness and three-dimensionality in a laminar separation bubble[END_REF] and [START_REF] Rodríguez | Structural changes of laminar separation bubbles induced by global linear instability[END_REF]. In the latter study a link to the structure of stall cells was proposed in relation with the above mentioned questionable computations in [START_REF] Rodríguez | On the birth of stall cells on airfoils[END_REF]. [START_REF] Barkley | Three-dimensional instability in flow over a backward-facing step[END_REF] investigated a backward facing step and found a threedimensional instability. They studied flows at Reynolds numbers within a range between 450 to 1050 based on the height of the step, and associated the three-dimensional unstable mode to a centrifugal instability thanks to the criterion of [START_REF] Bayly | Instability mechanism in shear-flow transition[END_REF]. [START_REF] Gallaire | Three-dimensional transverse instabilities in detached boundary layers[END_REF] studied the flow over a bump at a Reynolds number of 400. They found a non-oscillatory mode causing the flow to become three-dimensional. A similar instability was found by [START_REF] Marquet | Direct and adjoint global modes of a recirculation bubble: lift-up and convective non-normalities[END_REF] for an S-shaped duct and by [START_REF] Picella | Successive bifurcations in a fully three-dimensional open cavity flow[END_REF] in open-cavity flows. Finally, a three-dimensional non-oscillatory unstable mode was found for a shockwave boundary layer interaction by [START_REF] Robinet | Bifurcations in shock-wave/laminar-boundary-layer interaction: global instability approach[END_REF]. These studies show that three-dimensional instabilities such as those found on aerofoils at low Reynolds numbers are also found on other configurations displaying recirculation bubbles.

Concerning transonic buffet, the two-dimensional phenomenon was investigated for a NACA0012 aerofoil by [START_REF] Mcdevitt | Static and dynamic pressure measurements on a NACA 0012 airfoil in the Ames High Reynolds Number Facility[END_REF] for a wide range of Mach number and angle of attacks at a Reynolds number of 10 million. Strouhal numbers (St = f c U ∞ ) in a range from 0.06 to 0.07 were identified, and the buffet onset boundary in the Mach number -angle of attack plane was identified. Similar buffet frequencies were obtained by [START_REF] Benoit | Buffeting prediction for transport aircraft applications based on unsteady pressure measurements[END_REF] for a RA16SC1 aerofoil and [START_REF] Jacquin | Experimental Study of Shock Oscillation over a Transonic Supercritical Profile[END_REF] for the OAT15A aerofoil. The latter study is now widely used to validate numerical models. Recently, [START_REF] Brion | Experimental analysis of the shock dynamics on a transonic laminar airfoil[END_REF] studied the OALT25, an aerofoil designed to promote laminar flow, with free transition and tripped boundary layer. They obtained Strouhal numbers of the order of 0.07 in the turbulent boundary layer case. From these experiments, turbulent transonic buffet is identified as a large amplitude oscillation of the shock position with a well-identified Strouhal number around 0.07.

Many numerical studies of two-dimensional transonic buffet have been carried out. Le Balleur & Girodroux-Lavigne (1989) and [START_REF] Edwards | Transonic shock oscillations calculated with a new interactive boundary layer coupling method[END_REF] used viscous-inviscid coupling strategies. However, the bulk of the numerical simulations reported in the literature is carried out with Unsteady Reynolds Averaged Navier-Stokes models (URANS) [START_REF] Goncalves | Turbulence model and numerical scheme assessment for buffet computations[END_REF]; [START_REF] Thiery | Numerical prediction of shock induced oscillations over a 2D airfoil: Influence of turbulence modelling and test section walls[END_REF]; [START_REF] Iovnovich | Reynolds-Averaged Navier-Stokes Study of the Shock-Buffet Instability Mechanism[END_REF]; [START_REF] Giannelis | Influence of Mach number and angle of attack on the two-dimensional transonic buffet phenomenon[END_REF]; [START_REF] Grossi | Prediction of Transonic Buffet by Delayed Detached-Eddy Simulation[END_REF]; Plante & Laurendeau (2019)). These simulations identified frequencies in the range of those obtained experimentally. However, the simulations are found to be very sensitive to the numerical scheme and turbulence model. Large Eddy Simulations (LES) were carried out by [START_REF] Garnier | Large-Eddy Simulation of Transonic Buffet over a Supercritical Airfoil[END_REF] and [START_REF] Fukushima | Wall-Modeled Large-Eddy Simulation of Transonic Airfoil Buffet at High Reynolds Number[END_REF] with excellent agreement to the experiment of [START_REF] Jacquin | Experimental Study of Shock Oscillation over a Transonic Supercritical Profile[END_REF]. However, such simulations are computationally very expensive. Hence, hybrid RANS-LES simulations were attempted [START_REF] Deck | Numerical Simulation of Transonic Buffet over a Supercritical Airfoil[END_REF]; [START_REF] Huang | Simulation of shock wave buffet and its suppression on an OAT15A supercritical airfoil by IDDES[END_REF]; [START_REF] Grossi | Prediction of Transonic Buffet by Delayed Detached-Eddy Simulation[END_REF]; [START_REF] Ishida | Transonic Buffet Simulation over Supercritical Airfoil by Unsteady-FaSTAR Code[END_REF]; [START_REF] Deck | Towards an enhanced protection of attached boundary layers in hybrid RANS/LES methods[END_REF]) and [START_REF] Ribeiro | Buffet Simulations with a Lattice-Boltzmann based Transonic Solver[END_REF] used the Lattice Boltzmann Method. These simulations predicted transonic buffet, but discrepancies to the experimental pressure distributions were found. Hence, URANS simulations reproduce the main feature of transonic buffet. However, the turbulence modelling has a strong effect on the quantitative results.

Three-dimensional transonic buffet is more complex. Most of the studies were carried out on half-wing body aircraft configuration. A broadband frequency content with a dominant frequency higher than the one of two-dimensional buffet was found experimentally [START_REF] Roos | The buffeting pressure field of a high-aspect-ratio swept wing[END_REF], [START_REF] Dandois | Experimental study of transonic buffet phenomenon on a 3D swept wing[END_REF], Paladini et al. (2019b), [START_REF] Masini | Analysis of a civil aircraft wing transonic shock buffet experiment[END_REF], [START_REF] Koike | Unsteady pressure measurement of transonic buffet on nasa common research model[END_REF]). On such configurations, the convection of flow perturbations towards the wing tip was observed and buffet cells can be seen in the pressure-sensitive paint measurements of [START_REF] Sugioka | Experimental analysis of transonic buffet on a 3D swept wing using fast-response pressure-sensitive paint[END_REF]. Numerical simulations of such configurations were carried out with URANS [START_REF] Sartor | Mach number effects on buffeting flow on a half wing-body configuration[END_REF]) and hybrid RANS-LES simulations [START_REF] Brunet | Zonal-Detached Eddy Simulation of Transonic Buffet on a Civil Aircraft Type Configuration[END_REF], [START_REF] Sartor | Delayed Detached-Eddy Simulation of Shock Buffet on Half Wing-Body Configuration[END_REF], [START_REF] Ishida | Transonic Buffet Simulation over NASA-CRM by Unsteady-FaSTAR Code[END_REF]). [START_REF] Ohmichi | Modal Decomposition Analysis of Three-Dimensional Transonic Buffet Phenomenon on a Swept Wing[END_REF] used Dynamic Mode Decomposition and Proper Orthogonal Decomposition to identify two dominant modes. One high-frequency mode (St ≈ 0.2 -0.6) associated to buffet cells and one at St ≈ 0.06, coherent with two-dimensional buffet. [START_REF] Iovnovich | Numerical Study of Shock Buffet on Three-Dimensional Wings[END_REF] studied simplified swept wings closed by a symmetry plane and an extrapolation boundary condition to emulate an infinite swept wing. They observed buffet cells and an effect of the sweep angle on the buffet amplitude and frequency. Previous work by the authors (Plante et al. (2019b)) investigated infinite swept wings closed by periodic boundary conditions. Such simulations produced much more regular buffet cells and showed that buffet cells occur without any three-dimensional disturbance from the physical setup, thus making this phenomenon a candidate for an explanation based on global stability analysis.

Transonic buffet was described as a feedback loop between the shock wave motion and acoustic waves generated near the trailing edge by [START_REF] Lee | Oscillatory shock motion caused by transonic shock boundary-layer interaction[END_REF]. The description of the transonic buffet as a globally unstable mode was introduced by [START_REF] Crouch | Origin of transonic buffet on aerofoils[END_REF] and further investigated by [START_REF] Iorio | Direct and adjoint global stability analysis of turbulent transonic flows over a NACA0012 profile[END_REF] and [START_REF] Sartor | Stability, Receptivity, and Sensitivity Analyses of Buffeting Transonic Flow over a Profile[END_REF]. Such analyses produce an accurate prediction of the buffet frequency and onset angle of attack angle. Three-dimensional global stability analyses were attempted by [START_REF] Iorio | Direct and adjoint global stability analysis of turbulent transonic flows over a NACA0012 profile[END_REF], but an unstable mode specific to three-dimensional buffet was not found. Recently, global stability analysis was applied to infinite swept wings using a spanwise periodicity assumption [START_REF] Crouch | Global Instability Analysis of Unswept and Swept-Wing Transonic Buffet Onset[END_REF] 2020)). Stability analysis has also been applied to the NASA Common Research Model by [START_REF] Timme | Global instability of wing shock buffet[END_REF][START_REF] Timme | Global Shock Buffet Instability on NASA Common Research Model[END_REF][START_REF] Timme | Global instability of wing shock-buffet onset[END_REF]. These analyses exhibit unstable modes coherent with the occurrence of buffet cells, thus emphasising the phenomenological differences between twodimensional and three-dimensional buffet.

This paper aims at showing that the buffet cell phenomenon occurring in transonic flow has the same origin as the stall cell one occurring during stall at low speed. In particular, we will show that a similar unstable global mode is at the origin of both phenomena and highlight discrepancies between linear dynamics predictions provided by global stability analyses and saturated nonlinear dynamics given by URANS simulations. In addition, the unstable global mode corresponding to the buffet/stall cell phenomenon is followed from low speed to transonic conditions to establish a close link between these two phenomena. First, the numerical methods are presented. Then, numerical solutions for low-speed stall and transonic buffet conditions are analysed. Linear global stability analyses around steady base-flow are presented and the transition from the linear to the nonlinear dynamic is investigated. Finally, the link between the two flow regimes is discussed.

Physical and numerical set-up

This section first describes the physical setup studied in this paper. Then the governing equations and the numerical methods used to solve them are presented.

Configurations

In this study, infinite swept wings are considered. Fig. 1 shows the numerical setup. Usually, aerodynamic analyses are carried out in a referential (x, y, z) where the span of the wing is defined as the length L z along the z direction, and the sweep angle δ is the angle between the z-axis and the leading edge of the wing. With these definitions, the velocity in the z direction is null and the free-stream velocity is defined by its norm V ∞,3D and the angle of attack α 3D (see fig. 1). However, it is convenient for the present analysis to define another referential (x , y , z ) with x normal to the leading edge, y = y and z aligned with the leading edge. In this referential, the sweep angle becomes a sideslip angle applied to the free-stream velocity V ∞,3D .

Since the sweep angle will be changed in this article, we define the 2-D flow quantities as the conditions in the plane normal to the leading edge i.e. in the plane (x , y ). Hence we get the relations

x z y L z = L z c o s δ Periodicity δ V ∞,3D δ x z y = y S p a n -Lz α 3D α 2D c 3D = 1/cosδ c = c 2D = 1 Figure 1: Infinite swept wing configuration M 3D = V ∞,3D γP ∞ /ρ ∞ = M 2D cos(δ )
(2.1)

α 3D = arctan [tan(α 2D ) cos(δ )] (2.2) Re 2D = ρ ∞ V ∞,2D c 2D µ ∞ (2.3) Re 3D = ρ ∞ V ∞,3D c 3D µ ∞ = Re 2D cos(δ ) 2 (2.4)
For this study, the 2-D flow conditions are kept constant when the sweep angle changes to maintain the similarity. This means that the Mach number M 2D , the Reynolds number Re 2D and the angle of attack α 2D are kept constant. For convenience the 2-D subscripts are dropped (in the following M = M 2D , Re = Re 2D , α = α 2D , etc.) and the reference quantities are defined as

V ∞ = V ∞,2D = 1.0, c = c 2D = 1.0. For example, the Strouhal number is defined as St = f c 2D V ∞,2D
= f . This is done to help the comparison of non-dimensional quantities based on these reference quantities when the sweep angle varies. The aspect ratio AR is defined as AR = L z c 2D = L z and it should be noted that the baseline aerofoil is always retrieved in the plane (x , y ). The (x , y , z ) referential is more adapted for the present study, and therefore, it will be the only referential used throughout this paper. As such, the prime symbol is dropped for the sake of conciseness.

The present paper focuses on two regimes: the subsonic stall regime and the transonic buffet regime. The subsonic study is based on the NACA4412 aerofoil at a Reynolds number Re = 350 000 and a Mach number M = 0.2 such that compressibility effects are negligible. This case is selected because it exhibits a trailing edge type of stall. The transonic case focuses on the ONERA OALT25 aerofoil at M = 0.7352 and Re = 3 × 10 6 , such that the experimental 2-D results of [START_REF] Brion | Experimental analysis of the shock dynamics on a transonic laminar airfoil[END_REF] can be used for comparison. This aerofoil was designed to promote laminarity and has been studied experimentally both with free transition and a boundary layer tripped at 7% on both sides. With the tripped boundary layer, the buffet instability is similar to that of other turbulent aerofoils. In the present study, the laminar region is neglected and simulations are carried out in fully turbulent conditions. Hence, the experimental data with a tripped boundary layer are used to validate the numerical set up for two-dimensional buffet computations. The final part of the paper is dedicated to a third aerofoil, the NACA0012, at a Reynolds number Re = 1 × 10 7 and various Mach numbers and angles of attack. This aerofoil is used to link the results from the two previous cases, since it allows to track phenomena from the low-Mach stall regime to the transonic buffet regime. This particular test case is selected since the experiments of [START_REF] Mcdevitt | Static and dynamic pressure measurements on a NACA 0012 airfoil in the Ames High Reynolds Number Facility[END_REF] showed the transonic buffet in the high Mach range and this aerofoil is suitable to have a trailing edge type stall behaviour in the low Mach range. This makes it suitable to observe both the stall cells and the buffet cells. The aerofoils considered in this paper have a trailing edge thickness of 0.0025, 0.005 and 0.0025 chord units for the NACA4412, OALT25 and NACA0012 aerofoils respectively.

Governing equations

Non-linear model

In this study, the Reynolds-Averaged Navier-Stokes (RANS) equations are used to model the fluid flow. The one-equation Spalart-Allmaras turbulence model with the Edwards-Chandra modification is used to close the RANS equations system. This system reads in integral form

V V V ∂W W W ∂t + R R R( ( (W W W ) ) ) = 0 (2.5)
where W W W = [ρ , ρu , ρv , ρw , ρe , ρ ν] t is the conservative variables vector, V V V is a diagonal matrix whose non-zero coefficients are the mesh cells volume and R R R is the summation of the convective and viscous fluxes, and the source terms. In this article, z-invariant solutions (all the derivatives of the flow variables with respect to z are null) with w = 0 will be referred to as a 2-D solution, while z-invariant solutions with w = 0 will be referred to as 2.5-D solutions [START_REF] Ghasemi | A Two-Dimensional/Infinite Swept Wing Navier-Stokes Solver[END_REF][START_REF] Bourgault-Côté | Extension of a Two-Dimensional Navier-Stokes Solver for Infinite Swept Flow[END_REF]).

Linear model

Steady solutions W W W 0 of equation (2.5), also known as base-flows, are computed to carry out global linear stability analyses [START_REF] Sipp | Dynamics and Control of Global Instabilities in Open-Flows: A Linearized Approach[END_REF]; [START_REF] Theofilis | Global Linear Instability[END_REF]; [START_REF] Taira | Modal Analysis of Fluid Flows: An Overview[END_REF]. By definition, they satisfy R R R( ( (W W W 0 0 0 ) ) ) = 0, and under a first-order approximation, infinitesimal perturbations W W W about a base-flow W W W 0 are governed by the linearisation of equation (2.5), which reads

∂W W W ∂t = A A AW W W = -V V V -1 ∂ R R R ∂W W W W W W 0 W W W (2.6) where ∂ R R R ∂W W W W W W 0
is the Jacobian of the discretised RANS equations evaluated about the base-flow, including the turbulence model. Solutions in the form of normal modes are sought.

W W W = e λt Ŵ W W , (2.7)
such that eq. 2.6 reduces to the eigenproblem 1 reports the number of points on the aerofoil surface (n i ), in the direction normal to the aerofoil surface (n j ), on the pressures side (n p.s. ), on the suction side (n s.s. ), on the trailing edge (n t.e. ), the wall spacing (∆ y w ), the distance from the aerofoil surface to the farfield boundary condition (∆ y f . f . ) and the spacing across the shockwave (∆ x s. f . ).

A A A Ŵ W W = λ Ŵ W W . ( 2

Solution of the non-linear model

In this paper, the ONERA-Airbus-Safran elsA [START_REF] Cambier | The Onera elsA CFD software: Input from research and feedback from industry[END_REF]) software is used to solve the (U)RANS equations. For this study, a cell-centred structured finite volume method is used. A second order cell-centred scheme with scalar numerical dissipation is used for the discretisation of the convective fluxes. Convergence towards steady-state solutions is accelerated using local time steps and geometrical multigrid with a LU-SSOR pseudo-time integration scheme. In some case, Selective Frequency Damping (SFD) (Åkervik et al. (2006); [START_REF] Jordi | Encapsulated formulation of the selective frequency damping method[END_REF][START_REF] Jordi | An adaptive selective frequency damping method[END_REF]; [START_REF] Richez | Selective frequency damping method for steady RANS solutions of turbulent separated flows around an airfoil at stall[END_REF]) or a Newton solver [START_REF] Busquet | Global stability analysis of turbulent flows around an airfoil near stall[END_REF], [START_REF] Wales | Numerical continuation of high Reynolds number external flows[END_REF]) is used to force the convergence to a steady state. Time-accurate solutions are obtained by using a second order dual time stepping approach.

Solution of the linear model

This study involves three-dimensional geometries for which A A A has a high dimension and a large bandwidth. Therefore, solving the eigenproblem is often too costly. But in our cases, the base-flows are always homogeneous in the spanwise direction, such that the method proposed by [START_REF] Schmid | Stability analysis for n -periodic arrays of fluid systems[END_REF] and used by Paladini et al. (2019a) and Plante et al. (2019a) for nperiodic arrays of fluid systems may be used to significantly reduce the cost of the eigenvalues computation. This method exploits the block-circulant structure of matrix A A A, which takes the form

A A A =        A A A 0 A A A 1 A A A 2 . . . A A A n-1 A A A n-1 A A A 0 A A A 1 . . . A A A n-2 A A A n-2 A A A n-1 A A A 0 . . . A A A n-3 . . . . . . . . . . . . . . . A A A 1 A A A 2 A A A 3 . . . A A A 0        (2.9)
when the degrees of freedom are properly indexed. As shown by [START_REF] Schmid | Stability analysis for n -periodic arrays of fluid systems[END_REF], A A A can be transformed into a block diagonal matrix

 A A =        A A 0 0 • • • 0 0  A A 1 . . . . . . . . . . . . . . . 0 0 • • • 0  A A n -1       (2.10) with  A A j = n-1 ∑ k=0 ρ k j A A A k (2.11)
and

ρ j = e i j 2π n (2.12)
Eigenvalues of  A A are also eigenvalues of A A A and each block matrix  A A j j j can be treated separately, resulting in n separated eigenproblems of the size of a two-dimensional problem

 A A j j j v v v j = λ j v v v j .
(2.13)

The eigenvectors of A A A are then retrieved as

Ŵ W W = v v v j , ρ j v v v j , ρ 2 j v v v j , . . . , ρ n-1 j v v v j t (2.14)
This method implies to solve the eigenproblem for all the  A A j matrices to get the stability result for all the n possible wavenumbers. However, the only difference between each  A A j matrix is the coefficients ρ j such that changing the value of j or n has the effect of changing the wavenumber. Hence, solving for the matrice  A A 2 is equivalent to solving  A A 1 with n divided by two. Therefore, in this study, the problem is only solved for the matrix Â

A A 1 with various values of n. This allows to get the results over a continuous range of wavenumbers.

This analysis yields results similar to the biGlobal stability analysis carried out by [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF] and [START_REF] He | Linear instability of low Reynolds number massively separated flow around three NACA airfoils[END_REF]. However, this approach presents several advantages. In particular, it does not require tedious mathematical developments to include the periodicity assumption in the Jacobian matrix. This makes this approach suitable for the use of a black box CFD solver, in this case elsA.

The linear operator A A A is obtained with the fully discrete approach proposed by [START_REF] Mettot | Computation of eigenvalue sensitivity to base flow modifications in a discrete framework: Application to open-loop control[END_REF]. This Jacobian is extracted with a second-order finite difference scheme as

A A Au u u ≈ 1 2ε (R R R(W W W 0 0 0 + εu u u) -R R R(W W W 0 0 0 -εu u u)) (2.15)
where ε is a small parameter and u u u are vectors chosen to efficiently compute every non-zero coefficient of the Jacobian matrix (see [START_REF] Beneddine | Characterization of unsteady flow behavior by linear stability analysis[END_REF] for details). For the second-order finite volume scheme used in this study, the numerical stencils have a width of five grid cells. Hence,

A A A 3 = A A A 4 = . . . = A A A n-3 = 0,
and all the non-zero block matrices A A A j (eq. 2.9) can be retrieved from a three-dimensional grid with only five grid cells in the spanwise direction (see Paladini et al. (2019a); Plante et al. (2019a)). Hence, the base-flows are computed as two-dimensional solutions and the grid is extruded to get the three-dimensional mesh for the extraction of the 3-D Jacobian A A A. The spanwise spacing ∆ z must be imposed. Refinement study of this parameter will be presented, but no extra computational cost is associated with change in ∆ z, which can therefore be chosen arbitrarily small. This Jacobian matrix is then reduced as  A A 1 1 1 by assuming a given wavenumber in eq. (2.11). The eigenproblems are then solved with the Arnoldi iteration method [START_REF] Sorensen | Implicit Application of Polynomial Filters in a k -Step Arnoldi Method[END_REF]) coupled with a shift-and-invert strategy [START_REF] Christodoulou | Finding leading modes of a viscous free surface flow: An asymmetric generalized eigenproblem[END_REF]) to extract inner eigenvalues.

Unsteady Reynolds-Averaged Simulations

This section presents the results of URANS simulations for the low-speed and transonic cases. These simulations are carried out to observe stall and buffet cells as well as the superposition between the 2-D transonic buffet mode and the buffet cells.

Subsonic Stall

In this section the 2-D and 3-D simulations of the NACA4412 aerofoil at a Reynolds number of 350 000 are compared. Then the simulations of swept wings are presented.

3.1.1. Unswept δ = 0 • case
Let us first focus on the flow over an unswept wing in the post-stall regime. Fig. 2 shows the lift polar in the unswept case (δ = 0 • ) for two kinds of simulations: two-dimensional solutions and three-dimensional ones with an aspect ratio AR = 6 and 113 spanwise grid points. The maximum lift coefficient of 1.57 is obtained for an angle of attack around α = 14 • in both cases. For angles of attack higher than 16 • the flow solver failed to converge. This is caused by the occurrence of an instability linked to a vortex shedding. This is further discussed in section 4.1.4. In order to obtain steady-state two-dimensional solutions in this regime, the SFD method presented by [START_REF] Richez | Selective frequency damping method for steady RANS solutions of turbulent separated flows around an airfoil at stall[END_REF] and the Newton algorithm have been used. The SFD is used with a local time-stepping algorithm and the cut-off frequency and damping coefficient were based on a parametric study to obtain solutions converged within machine accuracy. A value of the filter width ∆ = 1/74 and the proportional controller coefficient χ = 50 allowed the solution to converge and the CFL is set to 50 to compute the local timestep in the LU-SSOR implicit scheme. To ensure consistency, it has been checked that both methods provide the same solution for several cases. Successive solutions are obtained by initialising each computation from the previous angle of attack. The lower branch of the lift polar is obtained by increasing the angle of attack up to values for which a converged solution could no longer be obtained. Then, by decreasing the angle of attack starting from these unconverged flow solution, the solution method is able to converge on the lower branch of the lift curve. Every solution has been converged within machine accuracy. Two distinct solutions are found for some angles of attack. Such hysteresis phenomena of the discretised RANS equations were obtained by [START_REF] Wales | Numerical continuation of high Reynolds number external flows[END_REF] and [START_REF] Busquet | Global stability analysis of turbulent flows around an airfoil near stall[END_REF] for the NACA0012 and the OA209 aerofoils, respectively. Grid convergence has been ensured by carrying out a two-dimensional simulation at an angle of attack of 15 • with a grid refined by a factor 2 in both directions (1024 by 256 grid cells), resulting in non-significant differences in the pressure distribution.

For angles of attack over 14 • the three-dimensional lift coefficient is lower than the twodimensional one. These solutions are obtained using a global time step and the SFD technique. The lower lift coefficient is associated with the presence of stall cells, as seen in fig. 3 for an angle of attack of 15 • . Four steady stall cells are observed and the skin friction lines exhibit the "owlface" shape reported in the literature. The left part of this figure shows the result with a spanwise mesh refinement of a factor 2. No significant change in the result is noticeable, ensuring the grid convergence of the results. In this simulation, the wavenumber of the stall cell is 4.2. However, the periodic boundary conditions constrain the possible wavenumbers since the number of cells must be an integer. Hence, the only possible wavenumbers that may appear in the simulation are necessarily multiples of β 0 = 2π/6. There might therefore exist a discrepancy with respect to a truly infinite configuration, with AR = ∞. A more accurate estimation of the most amplified wavenumber may be obtained by considering a larger aspect ratio. For this reason, a simulation with AR = 12 has also been run and 8 stall cells were obtained. Therefore, the most amplified wavenumber is narrowed down to a value between 14π/12 ≈ 3.7 and 18π/12 ≈ 4.7. This range will be compared to the result of linear stability analyses in section 4.1.

Since the 2-D solutions are particular cases where every derivative in the spanwise direction are null, they are also solutions of the three-dimensional discretised RANS equations. Hence these results demonstrate the existence of multiple solutions, the two-dimensional and the three-dimensional ones, for a given flow conditions. Such results are in agreement with those from [START_REF] Kamenetskiy | Numerical Evidence of Multiple Solutions for the Reynolds-Averaged Navier-Stokes Equations[END_REF].

3.1.2. Swept cases δ > 0 • , α = 15 •
A fixed value α = 15 • is now considered and a sweep angle is added. In these conditions, the flow is unsteady and time-accurate simulations with a non-dimensional time step of 0.008 are carried out. Fig. 4 shows the instantaneous pressure coefficient and the skin friction lines over wings swept at 10 • , 20 • and 30 • . Four stall cells are observed and the wavenumber is the same as for the unswept wing (β = 4.2). The topology of the stall cells also changes with the sweep angle, and the skin friction lines are now aligned in the cross-flow direction. This is associated with a spanwise convection of the stall cells. This induces a frequency proportional to the wavelength of the cells and their convection speed (St = V C/λ z ). Table 2 summarises all the wavelength and convection speed obtained for sweep angles of 0 • , 10 • , 20 • and 30 • . Fig. 5 shows the definition of these quantities. For the simulations at a constant aspect ratio, the wavenumber is always 4.2. These results suggest that the sweep angle has no effect on the wavelength of the stall cells. However, as mentioned previously, this result is constrained by the wing aspect ratio. Finally, the

Transonic Buffet

This section analyses the transonic buffet on wings based on the OALT25 aerofoil at a Reynolds number Re = 3 × 10 6 , a Mach number of 0.7352 and an angle of attack of 4 • . The non-dimensional time step is set to 0.0105. Simulations with the 2.5-D assumptions are first analysed. Then URANS results for 3-D infinite swept wings are presented.

3.2.1. Spanwise invariant solutions (2.5-D) Fig. 6 shows the time-averaged wall pressure coefficient of 2-D and 2.5-D simulations. The pressure coefficient on the pressure side and the trailing edge agree with the experiments of [START_REF] Brion | Experimental analysis of the shock dynamics on a transonic laminar airfoil[END_REF] and the pressure of the supersonic plateau is also the same. This indicates that the Mach number and angle of attack are the same between the experiments and the simulations. However, the shock wave of the simulations is located further downstream compared to the experimental field. It is seen that the pressure distribution is the same for every sweep angle, except in the shock wave region. This shows that the flowfields in the (x, y) plane remain close due to the fact that the upstream normal inflow conditions are kept constant when changing the sweep angle. The time evolution of the lift coefficient yields a Strouhal number of 0.075. The shock wave motion amplitude can be evaluated from the slope of the time-averaged pressure coefficient in the vicinity of the shock wave position. The largest buffet amplitude is observed for the unswept case and the buffet amplitude decreases as the sweep angle is increased. This relation between the sweep angle and the buffet amplitude is investigated with global stability analyses in section 4.2.3.

δ ( • ) AR N cells λ z β V C St 0 6.

superposition between 2-D buffet mode and 3-D buffet cells

Fig. 7 shows the instantaneous pressure coefficient and skin friction lines on swept wings of aspect ratio AR = 6 meshed with 168 spanwise grid cells. For the unswept wing and the wing swept at 30 • , a variation of the shock wave position is observed, forming the so-called buffet cells. For the unswept wing (δ = 0 • ), spanwise structures are observed. However, the amplitude of the buffet cells is irregular and the number of buffet cells varies in time. As such, the wavenumber varies in time and the buffet cells appear and disappear at random spanwise location since there is no convection. This results in a flow that can seem chaotic but a dominant Strouhal number of 0.06 is observed. This is related to the frequency of a chordwise variation of the shock wave position and a variation of the buffet cells amplitude. The same frequency is observed for swept wings (δ > 0 • ). However, in these cases the number of buffet cells is constant and they are convected in the spanwise direction. Two cells are observed for a sweep angle of δ = 30 • . This result is compared to the stability analysis results in section 4.2. Fig. 7 also shows the instantaneous solution with 112 and 168 spanwise grid cells. The comparison of the two fields shows that the grid spacing is sufficient for this case since no significant differences are noticeable in the skin friction lines and the pressure coefficient.

Fig. 8 shows the power spectral density of the sectional lift coefficient with 112, 168 and 224 spanwise grid cells. The sectional lift coefficient is extracted on a given chordwise wing section, the choice of this section has no effect since the buffet cells are convected in the spanwise direction. The sectional lift coefficient is sampled at every time-step (0.0105 nondimensional time-unit) over a time frame of 630 nondimensional time-units. A periodogram with boxcar window and a single block is used to estimate the power spectral density, which yields a Strouhal number resolution of 0.00158. This method is used since the signals for the grid refinement study are too short to use a sufficient number of averaging blocks. This yields noisier spectra and estimation error on the amplitude. Nevertheless, one can observe that the peak frequencies are the same. A dashed line is added to indicate the Strouhal number around 0.06 (the frequency observed for every sweep angle) and a dashed-dotted line indicates the Strouhal number of 0.135. The latter is associated with the convection of buffet cells, as it will be shown in this section. Another peak at a frequency of 0.135 could correspond to the nonlinear interaction between the two dominant modes (St = 0.135 + 0.055 = 0.19).

Fig. 9 shows the power spectral density of the sectional and global lift coefficient with sweep angles of 10 • , 20 • and 30 • . These spectra are evaluated using the Welch method with a Hann window and a 50% overlap between the blocks. The signal length is of 1 260 nondimensional time-units and nine blocks are used. The Welch method is used since the time series is longer than in the mesh refinement study. The sampling is done at every time step (0.0105 nondimensional time-units) for the sectional lift and every ten time step for the global lift. This yields a Strouhal number resolution of 0.0039. One can notice the frequency between St = 0.055 and 0.06 indicated by a dashed line for every sweep angle in the global and sectional lift signals. Some harmonics of this frequency are also observed. This frequency is associated with a chordwise oscillation of the shock wave position which modulates the buffet cells amplitude (2-D oscillation). However, a second dominant frequency is observed in the sectional lift coefficient. The latter increases with the sweep angle and is due to the buffet cells convection in the spanwise direction. This second frequency is not observed in the spectrum of the global lift coefficient since the number of buffet cells is constant. For this reason, the variations of the lift coefficient associated with the convection are cancelled out when the lift coefficient is averaged along the span. Like in fig. 8 further discussed in the next paragraph. Some harmonics of the buffet cells convection frequency are also observed. For instance the convection Strouhal number of 0.045 for a sweep angle of 10 • has a harmonic at 0.09, which is also the convection frequency for a sweep angle of 20 • .

To visualize the unsteady behaviour of the flow, the pressure has been extracted along a line parallel to the leading edge. Fig. 10 is a (z,t) diagram showing this extraction of the pressure coefficient near x /c = 0.43 for the sweep angle of 30 • . Most of the time, the pressure is equal to the one on the supersonic plateau. However two buffet cells periodically cross the extraction line because the amplitude of the cells changes in time. This occurs with a period T = 18.0. This results in a Strouhal number equal to 0.055 (around the Strouhal number of 0.06 shown in fig. 9). A Strouhal number of this order is observed for every sweep angle, as previously shown, thus suggesting that the variation of the buffet cells amplitude relates to a two-dimensional buffet mode. In addition, this figure allows to follow the spanwise location of the buffet cells in time. Hence, the convection speed of the buffet cells can be extracted as the slope of an isocontour of pressure on this pressure map (here V C = 0.4) and the associated frequency can be computed with the wavelength of the buffet cells (St = V C/λ z ). These quantities are schematised in fig. 5 and reported for multiple sweep angles in Table 3. Fig. 11 shows the relation between the sweep angle and the convection speed. The cross-flow velocity varies in tan δ , meaning that the convection speed is proportional to the far field velocity projected along the leading edge. This proportionality is true for low-speed flow conditions as well. A fit coefficient of 0.7 is found for the buffet cells and 0.8 for the stall cells. More details are provided in Plante et al. (2019b). This relation between the sweep angle and the convection speed of the buffet cells is further analysed using stability analyses in section 4 for the stall and buffet conditions.

Partial conclusion

Three-dimensional flow features in the form of stall cells and buffet cells are observed with URANS modelling, in the low-speed stall regime and in the high-speed buffet regime respectively. These phenomena can be obtained with a numerical setup in which no disturbances are imposed by the initial conditions and the boundary conditions. This suggests that there is an instability inherent to these flow conditions. Moreover, the spanwise length of the cells remains approximately unchanged when the aspect ratio is increased. This suggests that there is a preferential wavenumber for these flow structures. However, finding a very accurate estimation of this wavenumber with URANS simulations would be computationally expensive since the aspect ratio constrains the number of cells which can be represented, and AR needs to be increased to reduce the uncertainty with respect to the most amplified wavelength. As it will be shown in the next section, stability analyses are a more appropriate tool to address this point. Finally, both types of cells are convected with a velocity proportional to the free stream cross-flow velocity. This causes an unsteady behaviour which is superimposed to a two-dimensional unsteadiness in the case of the transonic buffet simulations. These flow behaviours are analysed with global linear stability analyses in the next section.

Global Stability Analyses

The second part of this study uses global stability analyses to identify linearly unstable modes which could explain the stall cells and buffet cells observed in section 3. Results for the lowspeed and transonic flow conditions are presented and discussed. Parameters characterising the size and intensity of the recirculation zones in the base-flows are provided as supplementary material to this article.

Subsonic Stall

Steady 2.5-D and 2-D solutions of the NACA4412 aerofoil at a Reynolds number of 350 000, a Mach number of 0.2, and multiple angles of attack and sweep angles are converged within machine accuracy. These solutions are used as base-flows for the stability analyses. For consistency with the previous section, most of the analyses are carried out for an angle of attack of 15 • . The following focuses on the case δ = 0 • and α = 15 • . Fig. 12 presents the eigenvalues spectra obtained for two values of β . For both values, there is an unstable eigenvalue (positive real part) located along the real axis. This is a non-oscillatory unstable mode. A mode similar to this one has been observed by [START_REF] Zhang | Biglobal linear stability analysis on low-re flow past an airfoil at high angle of attack[END_REF] ensure that modes with non-zero imaginary part are not missed. Fig. 13 shows the growth rate of the steady unstable mode for a range of wavenumbers β . Multiple values of ∆ z used in the extrusion of the grid are plotted to ensure that the grid is sufficiently refined. These results can be considered converged for every β with ∆ z = 0.002. One can observe that the grid convergence is poorer in the high β range associated with small structures. This indicates that a finer grid is required to compute small flow structures, as expected. The grid spacing used in the URANS simulations (∆ z = 0.054) is well converged up to β = 6, covering the case with the maximum growth rate and the wavenumbers obtained in the URANS simulations. This is consistent with the grid convergence shown in fig. 3. The maximum growth rate is obtained for a spanwise structure with β ≈ 4.4. This corresponds to a wavelength λ z ≈ 1.43, which is consistent with the value observed in the URANS simulations (β = 4.2). Indeed, in the simulations, the periodic boundary conditions constrained the possible wavenumbers to multiples of β 0 = 2π/AR. Therefore, the wavelength of the URANS simulation is as close as possible to that predicted by the stability analysis, showing that this unstable mode is responsible for the formation of stall cells. Fig. 14 shows the base-flow and the real part of the eigenvector ρe reconstructed in three dimensions with eq. 2.14, and extracted on the aerofoil surface and the plane z = 0. The instability is located in the shear layer and upstream of the separation line. The shape of this mode is consistent with the structure of the stall cells. 4.1.2. Case δ = 0 • and varying α.

These stability analyses are carried out for several angles of attack for which a solution converged to machine accuracy can be obtained (fig. 2). The three-dimensional unstable mode can be identified on the upper and lower branch of the lift polar hysteresis. Fig. 15 shows the growth rate of this mode for the most amplified wavenumber of the two branches. The mode is unstable from an angle of attack of 14.15 • (the maximum lift) up to 20.5 • (the last converged base-flow). On the lower branch, the growth rates are higher than for the upper one, but the mode is only unstable for angles of attack lower than 16.9 • . Above this angle, the mode becomes stable again, even though the flow is massively separated. Fig. 16 shows the shape of the modes on the upper and lower branches. One can observe that the unstable mode is located around the region where the flow separates. For α = 18 • the flow separates far from the leading edge. However, for the other angles of attack the separation is close to the leading edge. Hence, it is possible that the mode is no longer unstable past α = 16.9 • on the lower branch because the flow separation becomes fixed on the leading edge. This tends to validate the idea that stall cells are only observed with a trailing edge type of stall. Works based on the lifting line theory by [START_REF] Spalart | Prediction of Lift Cells for Stalling Wings by Lifting-Line Theory[END_REF] and [START_REF] Gross | Criterion for Spanwise Spacing of Stall Cells[END_REF] suggested that a negative slope of the C Lα curve is a necessary condition for the occurrence of stall cells. Fig. 15 shows the growth rate obtained in these flow conditions for different values of α, and for the upper branch of the lift hysteresis, the most unstable eigenvalue crosses the imaginary axis at an angle of attack slightly below 14.15 • . This coincides with the maximum lift angle of attack in fig. 2, which is therefore consistent with the negative C Lα slope criterion. However, the stability analysis of the lower branch of the lift hysteresis shows that the stall cell modes are stable for α > 16.9 • , while the slope of the C Lα curve is negative between 17.0 and 17.4 • (see fig. 2), which contradicts this criterion. Another conclusion of the analysis of [START_REF] Gross | Criterion for Spanwise Spacing of Stall Cells[END_REF] is that the wavelength of the cells should be λ z = -0.5 ∂C L/∂ α with a discrete model or λ z =π /4 ∂C L/∂ α with a continuous model. Since the linear instability occurs around the twodimensional solution, ∂C L/∂ α is taken as the slope of the two-dimensional lift polar around an angle of attack of 15.0 • (-2.52rad -1 ). With this slope, the wavelength should be 1.26 or 1.98 (β = 3. the mode becomes stable. The value of β for which the growth rate is maximum also slightly decreases as the sweep angle increases. The stability results are consistent with the URANS simulations from section 3. For example, the wavenumber observed in fig. 3 in the unswept case is 4.2 with 4 stall cells. By accounting for the numerical constraint on the wavenumber (the domain size imposes the possible values of β ), the actual wavenumber on an infinite domain is expected to be between 3.7 and 4.7. The modes with the maximum growth rate have β between 4.0 and 4.4, hence within the expected range.

Fig. 17b also shows the frequency obtained with URANS simulations (red points) and the frequency extrapolated from the trend of the convection speed (see fig. 11) and the wavenumber (red lines). The frequencies observed in the URANS computations are well aligned with the frequency obtained with the global linear stability analysis. Hence, the linear stability analyses predict the behaviour of the stall cells. 4.1.4. Two-dimensional wake instability mode.

These analyses are now carried for β = 0, which corresponds to two-dimensional analyses. Such analyses are performed for several angles of attack. Fig. 18 shows the part of the spectra where unstable modes are found. The mode becomes unstable for an angle of attack between 16 • and 17 • for the upper branch (UB) of the lift hysteresis, and is unstable for the entire lower branch (LB). This is consistent with the need to use the SFD technique or a Newton solver to converge the base-flow for angles of attack above 16 • , where it was observed that the local time stepping technique failed to converge. Fig. 19 shows the eigenmodes for the upper and lower branch of the lift hysteresis at an angle of attack of 18 • . The shape of these modes corresponds to a meandering of the wake which would lead to shedding of vortices in the nonlinear regime. The wavelength of this mode is larger for the lower branch, where the size of the separation is larger. As the angle of attack increases, the size of the separation region becomes larger and the Strouhal number based on the chord length (St c ) decreases. However, the Strouhal number based on the frontal height of the separation zone (St l = St c l = St c (x rx s ) sin α) stays constant around 0.2. This Strouhal number definition based on the separation height is inspired from the one of Yon & Katz (1998) and [START_REF] Zaman | A natural low-frequency oscillation of the flow over an airfoil near stalling conditions[END_REF], but adapted to the case where the separation point x s is not located at the leading edge but somewhere between the leading and the trailing edge. Here the flow re-attaches at the trailing edge (x r = 1.0). This shows that the Strouhal number based on the separation height is the correct one and its value is close to the expected one for vortex shedding [START_REF] Chabert | Experimental closed-loop control of separated-flow over a plain flap using extremum seeking[END_REF]). Hence, for α 16.5 • , the three-dimensional non-oscillating mode coexists with a two-dimensional wake instability mode, such that an unsteady behaviour should be observed in the URANS simulations. On the other hand, the case α = 15 • was analysed in section 3.1, and consistently with the stability results (only the three-dimensional steady mode is unstable for this value of α), the unsteady computation converges toward a steady state. Hence, once again, the results of section 3 are consistent with the stability analyses. The same stability analyses are now carried out for sweep angles of 0 • , 10 • , 20 • and 30 • and an angle of attack of 18 • . Fig. 20 shows the value of the growth rate and the Strouhal number of the wake instability mode for values of β around zero. This unstable mode exists for large wavenumbers β . In these figures, modes with a negative value of β are shown. These modes can also be seen as modes with positive β and a negative frequency. One can observe that, as expected, the spectrum is symmetric with respect to the imaginary axis when the sweep angle is null. In this case, the base-flow indeed exhibits a zero spanwise velocity and perturbations may equivalently propagate (in terms of phase, the phase-speed being v φ = ω/β ) towards the inboard (negative z, v φ < 0) and outboard (positive z, v φ > 0) directions. If a positive sweep is considered, then a positive base-flow velocity in the z direction is added and the phase of the perturbation is generally preferentially advected in this direction (in the case of a constant W (x, y) = W base-flow velocity, one can show that the frequency ω (resp. phase-speed v φ ) of a mode just shifts by an amount of βW (resp. W ), due to the Doppler-shift). Therefore, a positive sweep generally tends to increase the frequency of 3-D perturbations and their phase speed. For exactly two-dimensional oscillating perturbations (as obtained for the wake instability mode), the phase velocity is infinite (positive or negative), even if sweep is added. It is also seen that, if the perturbations become slightly 3-D (β = 0), then the outboard propagating 3-D perturbations are slightly more unstable than the inboard propagating ones. Note that the discussion would have been different if the group velocity (v g = dω/dβ ) was considered, since the group velocity determines the speed of the envelope and fronts of a wave-packet and therefore better accounts for actual propagation of the perturbation: for this a convective/absolute stability analysis should be conducted following [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF]. But this is outside the scope of this paper.

Transonic Buffet

The stability of base-flows obtained in the transonic regime for the OALT25 aerofoil at a Mach number of 0.7352 and a Reynolds number of 3 million is now studied. To begin with, the unswept case is analysed for several angles of attack. Then the effect of the sweep angle is analysed for a constant angle of attack of 4 • . Finally, the two-dimensional case is analysed.

Effect of the angle of attack on the unswept case

Fig. 21 shows the spectra obtained for β = 2π and β = 0 for an angle of attack of 4 • , and a sweep angle of 0 • . The cases where β tends to 0 are equivalent to the two-dimensional cases. As one can see, an unstable global mode is found at a Strouhal number of 0.075 when β = 0 (unsteady two-dimensional) and a Strouhal number of zero when β = 2π (non-oscillating threedimensional). Hence, two unstable eigenmodes are found for these flow conditions. This result is analogous to that at low speed where both a stall cells and a wake instability mode were observed. It should be noted that the spectra are symmetric with respect to the real axis, due to the symmetry of the problem in the non-swept case (a wave may propagate in the positive and negative z directions). Thus, an unsteady mode with a negative frequency is also observed. Following this analysis, fig. 22a shows the effect of the wavenumber β on the growth rate of the non-oscillating unstable mode for several angles of attack. All the analyses are done with a spanwise grid spacing of 0.002. An analysis similar to the one shown in fig. 13 has been carried out. It shows that refining the grid by a factor 100 does not change significantly the results for β up to 70 and for α = 4.0 and 8.0. It should also be noted that the spanwise grid spacing used in the URANS simulations (∆ z = 0.035) is properly converged for the range of unstable wavenumbers at an angle of attack of 4 • . As for the subsonic stall case, a bump in the growth rate is observed for a medium wavelength around β = 6.0 for the low angle of attack. However, a second bump with a maximum growth rate at higher β is obtained as the angle of attack increases. For the particular angle of attack of 4.25 • , the two bumps can clearly be distinguished. These results are similar to those obtained by [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF] and Paladini et al. (2019a) on the OAT15A aerofoil, where an intermediate-wavelength mode and a shortwavelength one were reported. [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF] suspected that the short-wavelength might be too small to be relevant to a RANS framework. Fig. 22b shows the spectra obtained for several β in between the intermediate-wavelength and short-wavelength bumps at an angle of attack of 4.25 • . One can notice that there is only one unstable eigenvalue for each β . We do not observe a 2-mode switching, where the originally unstable mode would become stable while a previously stable mode would become unstable. This shows that these two bumps are in fact the same mode and not two separate modes which could have been both unstable for some angles of attack. Although this does not discard arguments on the physical relevance of the low-wavelength mode, the physical mechanism responsible for both these modes is the same. Fig. 23 shows the C Lα curve of the steady base-flow of the OALT25. As it is the case for the NACA4412, the angles of attack for which the linear stability shows an unstable three-dimensional mode corresponds to angles of attack where the C Lα slope is negative, in agreement with [START_REF] Gross | Criterion for Spanwise Spacing of Stall Cells[END_REF]. However, the slope of this curve around α = 4 • and the most unstable wavelength does not match the predictions of the model of [START_REF] Gross | Criterion for Spanwise Spacing of Stall Cells[END_REF]. If one takes the slope between 5 • and 9 • , wavenumbers of 6.5 and 4.13 are predicted with the discrete and continuous models respectively. The unstable mode has a wavenumber β = 6.0, which is in the predicted range. However, as it was the case for subsonic mode, changes in the C Lα slope are not associated with changes in the wavelength predicted by stability analyses. 4.2.2. Sweep angle effect at α = 4 • . Fig. 24a shows the effect of the sweep angle on the growth rate of the three-dimensional mode. As for the subsonic stall, the growth rate of the three-dimensional mode decreases as the sweep angle increases. From this trend, the mode becomes stable for a sweep angle slightly higher than 40 • and is stable at 50 • . Also, the wavenumber for which the mode is the most unstable slightly decreases when the sweep angle increases. However, the wavenumber obtained in the swept URANS (β = 2.09) simulations is not unstable in the stability analysis. It will be shown in section 5 that this discrepancy is due to nonlinear effects. Fig. 24b shows the Strouhal number of the three-dimensional mode for sweep angles δ of 0 • , 10 • , 20 • , 30 • and 40 • for a range of β . As expected, the frequency increases with β (smaller spanwise structures) and with the sweep angles (higher convection speed). The URANS results (red points) and the frequencies extrapolated from the convection speed of the URANS simulations (red lines) are plotted on this graph as well. The frequencies obtained by [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF] on the OAT15A aerofoil are also plotted. For the latter, the values have been corrected for the fact that the reference speed is taken as the speed normal to the leading edge instead of the free stream velocity. These curves match the stability results, even though the work of [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF] was performed on a different aerofoil. Hence, the linear stability is able to predict the convection speed of the cells, but because of the discrepancy on the wavenumber with respect to the URANS simulations, the frequency of the simulations is incorrectly predicted by the stability analyses. This is due to nonlinear effects that shift the wavelength of the cells (see section 5).

Fig. 25 shows the base-flow with a sweep angle of 30 • , and the real part of the ρe twodimensional and three-dimensional unstable modes. The two unstable modes are mostly located in the vicinity of the shock wave and in the shear layer. The two-dimensional mode displays the shape of the unstable mode reported in two-dimensional buffet studies of [START_REF] Crouch | Origin of transonic buffet on aerofoils[END_REF] and [START_REF] Sartor | Stability, Receptivity, and Sensitivity Analyses of Buffeting Transonic Flow over a Profile[END_REF]. The three-dimensional mode has spanwise fluctuations in the shock foot and trailing edge region. Those are regions of flow separation. The shape of this mode is in line with the observations of [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF] and Paladini et al. (2019a).

2-D buffet mode

Fig. 26 shows the portion of the spectrum where the unsteady two-dimensional mode is located. One interesting point is that the unstable unsteady mode (two-dimensional buffet) is only observed over a range of angles of attack between 3.75 • and 6.0 • . The three-dimensional buffet mode becomes unstable at an angle of attack around 3.75 • as well. However, it does not become stable again at higher angles of attack in the range studied in this paper. Hence, this mode is observed for transonic stall conditions, where the two-dimensional buffet is stable. Although this flow condition is an extreme one, which is not expected to be reached for a civil aircraft, this result suggests that the 3-D buffet can exist without the 2-D buffet phenomenon occurring. This is similar to the fact that the stall cells can be obtained without the wake instability mode.

Finally, it is possible to analyse the effect of β on the two-dimensional buffet mode. Fig. 27 shows this effect on the growth rate and the frequency for several sweep angles. These are modes with very large wavelengths of at least 10 chord lengths. It should be noted that the spectrum is symmetric with respect to the real axis for δ = 0 • . But it is not the case for swept wings. Hence, the effect of β on the growth rate and frequency is not the same for the positive and negative frequency. The two-dimensional buffet mode is more stable as the sweep angle increases. This is consistent with the fact that the buffet amplitude is lower when the sweep angle increases as shown in fig. 6. Also, β has an effect on the growth rate. Such results are similar to those of [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF] for the OAT15A aerofoil. However, they did not report the peak in the growth rate between β = ±0.5 and β = ±0.6. As verification, the case of [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF] has been executed with the numerical setup of the present paper with similar results (see appendix A). Hence, these peaks are not caused by any problem related to the stability method or its implementation, and their existence remains unexplained for the OALT25 aerofoil. Also, the inboard propagating modes (β < 0) are more unstable than the inboard propagating modes when the wing is swept. This behaviour is opposite to the tendency of the wake instability mode (fig. 20).

Partial conclusion

The global linear stability analyses predict two unstable modes in the subsonic stall conditions. One which can be related to the stall cells (3-D) and the other to the shedding of vortices (2-D). The wavenumber, convection speed and frequency of the 3-D mode are consistent with the URANS results presented in section 3.1. The 2-D mode has the expected frequency with a Strouhal number based on the separation height of 0.2. However, multiple frequency and size of the vortices are identified for some flow conditions since multiple steady-state base-flows are found.

Two unstable modes are also identified in the transonic buffet regime. The first one is associated to the 2-D buffet instability at a frequency around 0.07. The second one leads to spanwise flow structures which can be associated to buffet cells. There is a discrepancy between the wavenumber predicted by the stability analyses and the URANS simulations. However, the frequency of this mode is that of 3-D structures convected at the same speed as those in the URANS simulations.

For the case of the NACA4412, two solution branches have been analysed. Such multiple solutions have been observed in the literature before [START_REF] Wales | Numerical continuation of high Reynolds number external flows[END_REF]; [START_REF] Busquet | Global stability analysis of turbulent flows around an airfoil near stall[END_REF]). Multiple base-flows were not found for the OALT25 aerofoil, however, they cannot be ruled out. Nevertheless such multiple solutions do not seem to be related to the occurrence of stall cells since an unstable stall cell mode is found at angles of attack where only one base-flow is found.

From linear to non-linear regime

Section 3 analysed the non-linear saturated regime. This was done by starting the solution from a non-converged RANS solution. Section 4 studied the behaviour of the solution in the linear regime around a fully converged RANS solution. This section investigates the transient between the two regimes. To do so, the steady base-flows are extruded to get a three-dimensional initial condition, and the URANS solver is run from this initial solution. Fig. 28 shows how the solution diverges from the steady two-dimensional base-flow towards the solution shown in section 3. As one can see, the flow diverges with four regular stall cells and the flow topology converges to the one shown previously. This corresponds to a wavenumber β = 4.2 which is the wavenumber with the maximum growth rate allowable on this grid (see fig. 17a), since the number of cells is constrained by the size of the domain. This result is consistent with the stability analysis. Similar analyses can be done for the transonic buffet case. Fig. 29 shows three time samples of the URANS simulation with the extruded base-flow as the initial state and a sweep angle of 20 • . One can observe that three-dimensionality occurs with the appearance of 6 buffet cells (β = 6.28), in agreement with the largest growth rate allowable for this computational domain (see fig. 24a). However, the wavelength progressively increases toward the value observed in section 3. This shows that the stability analysis predicts the initial bifurcation from the base-flow to the presence of buffet cells. However, some non-linear effects are responsible for the shift of wavelength observed such that the saturated regime is different from the value predicted by the linear model.

Link between stall and buffet cells

Previous solutions were obtained for two different aerofoils, at two Reynolds numbers. To clearly establish a link between buffet cells and stall cells, the stability analyses are now carried out for a NACA0012 at a Reynolds number of 10 million and the unstable eigenvalues are tracked from the low-speed regime to the transonic one. The NACA0012 aerofoil at a high Reynolds number is selected since a trailing edge stall is expected at low Mach number and transonic buffet at higher speed [START_REF] Mcdevitt | Static and dynamic pressure measurements on a NACA 0012 airfoil in the Ames High Reynolds Number Facility[END_REF]. As the Mach number increases, the angle of attack is decreased in order to recover the stall and buffet onset conditions. Fig. 30 shows the base-flow obtained for Mach numbers from 0.2 to 0.72. The flow topology at a Mach number of 0.2 is the one of a subsonic stall. Similar topology, but with a stronger flow acceleration and a weak shock wave near the leading edge are obtained for Mach numbers of 0.3 and 0.4. However, for Mach numbers higher than 0.5 the shock wave is stronger and flow separation occurs at the shock foot. Fig. 31a shows the growth rate of the three-dimensional unstable mode. For the low-speed cases, the maximum growth rate is found for β around 4. This result is similar to the previous observations on the NACA4412 aerofoil, even though the Reynolds number is by far larger. For Mach numbers between 0.4 and 0.5, the flow topology changes significantly and so does the value of the most amplified β . In this case β is greater than 170. These are three-dimensional modes with a very small wavelength. As the Mach number increases, the value of β for which the three-dimensional mode is found decreases to reach a value of 100 for a Mach number of 0.6 and 50 for a Mach number of 0.72. This is consistent with the results of [START_REF] Crouch | Global Instability Analysis of Unswept and Swept-Wing Transonic Buffet Onset[END_REF], who reported only a short-wavelength mode for the NACA0012 aerofoil and not the intermediatewavelength like on the OAT15A and OALT25. One major difference between the NACA0012 and the OALT25 aerofoils is the fact that there is only a flow separation in the shock foot area for the NACA0012, while the flow over the OALT25 also separates at the trailing edge. Figs. 31(b) shows the growth rate of the three-dimensional mode normalized by the height of the separation h s . By using this normalization, the maximum growth rates are located around β h s ≈ 0.25 for every Mach number, which shows that this is the correct scaling length for this three-dimensional mode. The amplification factors are normalized by the same length scale. Two other scalings were attempted, namely a scaling by the separation length l s and by the projection of the separation length in the direction normal to the inflow l (see supplementary material to this article). By using the separation length, a proper scaling is obtained for the subsonic and transonic cases separately. But this scaling does not hold between these two flow regimes. Fig. 32 displays the eigenvectors at these flow conditions, indicating that the unstable modes are all located in the vicinity of the flow separation. For a Mach number larger than 0.5, a shock wave occurs in the base-flow. In these cases, the eigenmode reaches high values at the separation point and further downstream, but it does not propagate in the supersonic region upstream of the shock wave, contrary to subsonic cases.

In the case of the OALT25 aerofoil, the three-dimensional and two-dimensional modes become unstable around the same angle of attack. This may suggest that the two unstable modes are somewhat linked. For this reason, the angle of attack at which an unstable 2D buffet mode appears has been searched for the NACA0012 at M = 0.72. Fig. 33 shows the region of the spectra around which the unstable mode is expected. An unstable mode is present at an angle of attack of 4.5 • , but this mode is stable for lower angles of attack. This indicates that the three-dimensional buffet mode can be unstable (see unstable branch at α = 3.5 • in fig. 31) without an underlying twodimensional mode (fig. 33).

As it is the case for the NACA4412 and the OALT25 aerofoils, 3-D modes with a non-zero frequency are obtained when a sweep angle is added. Fig. 34 shows the frequency of the 3-D modes with sweep angles of δ = 0 • , 10 • , 20 • and 30 • , and the scaling found previously. The angles of attack are sligthly increased in comparison to previous results on the NACA0012 because the sweep angle has a stabilizing effect. As one can see, the frequency of these unstable modes is still proportional to the spanwise velocity. Moreover, the proportionality coefficients found for the NACA4412 and the OALT25 still hold for the NACA0012 in subsonic and transonic regimes respectively. Hence it is possible to follow a single and same unstable mode that represents stall cells in the subsonic regime and buffet cells in the transonic flow regime. This branch of modes exhibits the same features as those found separately for the NACA4412 and the OALT25 aerofoils. This also shows that the stall cell mode is present both at an intermediate (350 000) and high (10 7 ) Reynolds number. Moreover, a scaling by the height of the separation bubble is proposed to link the size of the stall cells to that of the buffet cells.

Conclusion

Three-dimensional flows over infinite swept configurations were obtained using URANS simulations for low-speed stall and transonic buffet conditions. The observed stall/buffet cells are steady for unswept wings, and convected in the spanwise direction when the wing is swept. This induces a frequency which can be related to the distance between two cells and their convection speed. In the case of the subsonic stall, this unsteadiness dominates the flow. However, it is superimposed to the two-dimensional buffet phenomenon at high-speed.

Global linear stability analyses of 2.5-D base-flows were carried out for subsonic stall and transonic buffet. An unstable three-dimensional mode was observed in both flow conditions. The frequency of the latter is null for unswept wings and becomes unsteady as the sweep angle is increased. The frequencies and wavelength of stall cells obtained with stability analysis are in line with the observation of the URANS simulations. However, discrepancies between the wavelength of buffet cells and the global mode with the largest growth rate are observed. Nevertheless, the frequency of the unstable modes is consistent with that of structures convected at the convection speed observed in the URANS simulations. A mode consistent with a two-dimensional vortexshedding in subsonic stall conditions was found, as well as a two-dimensional buffet mode in transonic buffet conditions. The two-dimensional vortex shedding is not observed in the URANS simulations because the angle of attack is lower than the one for the onset of this unstable mode. A similar condition where the buffet cell mode is present but not the two-dimensional buffet mode might allow to obtain steady buffet cells on an unswept wing.

Nonlinear flow simulations starting from the base-flows effectively shows that the flows diverge following the shapes identified in the stability analyses. In the case of the subsonic stall, the number of cells predicted in the stability analyses stay constant in the nonlinear regime. However, the solution of transonic buffet over swept wings diverges towards buffet cells with a lower wavenumber after the appearance of spanwise flow structures of the length identified with the stability analyses, possibly because of non-linear effects.

Finally, results of the stability analyses allow to follow the angle of onset of a threedimensional mode from the low-speed to the transonic flow regime. This study was done on a NACA0012 aerofoil at a Reynolds number of 10 million. Although the wavelength of this mode varies greatly with the Mach number, using the separation height as the reference length allows to centre the non-dimensional wavenumber around a value β h s ≈ 0.25. This result reinforces the conclusion that the two phenomena are the same and are related to the occurrence of flow separation. Also, the three-dimensional mode is identified at an angle of attack lower than the two-dimensional buffet onset angle. This shows that the buffet cells associated with three-dimensional buffet can be observed without the presence of the two-dimensional buffet.

Hence, this paper provides evidences that the buffet cells and the stall cells are in fact the same flow instability. Even though it is convenient to describe the three-dimensional mode in buffet condition as a three-dimensional buffet phenomenon in the context of applied aerodynamics, one can show that this mode is oscillatory only when the wings are swept. Moreover, the physical mechanisms involved in the two-dimensional buffet and the buffet cells are different. As such, one can question the accuracy of naming the so-called buffet cells as three-dimensional buffet and not just as stall cells.

  ,Crouch et al. (2019); Paladini et al. (2019a); Plante et al. (2019a)) and fully three-dimensional analyses (Paladini (2018); He & Timme (

Figure 2 :

 2 Figure 2: Lift curve for the 2-D and 3-D steady solutions (NACA4412, Re = 350 000, M = 0.2, δ = 0 • ).

Figure 3 :

 3 Figure 3: Surface pressure coefficient and skin friction lines of 3-D steady solution with 224 (left) and 112 (right) spanwise grid cells (NACA4412, M = 0.2, Re = 350 000, α = 15 • , δ = 0 • , L z = 6).

  (a) Sweep 10 • . (b) Sweep 20 • . (c) Sweep 30 • .

Figure 4 :

 4 Figure 4: Instantaneous surface pressure coefficient and skin friction lines over a swept wing (NACA4412, M = 0.2, Re = 350 000, α = 15 • , L z = 6).

Figure 5 :

 5 Figure 5: Stall and buffet cells convection model.

Figure 6 :

 6 Figure 6: Time-averaged wall pressure coefficient for spanwise invariant solution with several sweep angle (OALT25, M = 0.7352, Re = 3 × 10 6 , α = 4 • ).

  (a) δ = 0 • . (b) δ = 30 • (right 112 spanwise grid cells, left 168 spanwise grid cells).

Figure 7 :

 7 Figure 7: Instantaneous surface pressure coefficient and skin friction lines for URANS simulations of infinite swept wings with two sweep angles (OALT25, M = 0.7352, Re = 3 × 10 6 , α = 4 • , L z = 6).

Figure 8 :Figure 9 :

 89 Figure 8: Sectional lift coefficient power spectral density (OALT25, M = 0.7352, Re = 3 × 10 6 , α = 4 • , δ = 30 • ).

Figure 10 :

 10 Figure 10: Extraction of the pressure on a line parallel to the leading edge near x /c = 0.43 (OALT25, M = 0.7352, Re = 3 × 10 6 , α = 4 • , δ = 30 • ).

  δ ( • ) AR N cells λ z β Buffet cells convection frequency for several sweep angles (OALT25, M = 0.7352, Re = 3 × 10 6 , α = 4 • ).

FFigure 11 :

 11 Figure 11: Convection speed of stall (NACA4412, M = 0.2, Re = 350 000, α = 15 • ) and buffet (OALT25, M = 0.7352, Re = 3 × 10 6 , α = 4 • ) cells.

  4.1.1. Case δ = 0 • and α = 15 •

Figure 12 :

 12 Figure 12: Eigenspectrum (NACA4412, M = 0.2, Re = 350 000, α = 15.0 • , δ = 0.0 • ).

Figure 13 :

 13 Figure 13: Effect of spanwise grid density and the wavenumber on the growth rate of the unstable mode (NACA4412, M = 0.2, Re = 350 000, α = 15.0 • , δ = 0.0 • ).

  (a) ρe /ρe ∞ base-flow. (b) Real part of the ρe /ρe ∞ eigenmode (β = 4.4) -St = 0.0.

Figure 14 :

 14 Figure 14: Flow visualisation of the (a): base-flow and (b): 3-D unstable mode (NACA4412, M = 0.2, Re = 350 000, α = 15.0 • , δ = 0.0 • ).

  Figure 15: Effect of the angle of attack on the growth rate of the unstable mode (NACA4412, M = 0.2, Re = 350 000, δ = 0.0 • ).

  (a) Real part of the ρe /ρe ∞ eigenmode (α = 18 • , upper branch, β = 5) -St = 0.0. (b) Real part of the ρe /ρe ∞ eigenmode (α = 20.5 • , upper branch, β = 7.5) -St = 0.0. (c) Real part of the ρe /ρe ∞ eigenmode (α = 16.78 • , lower branch, β = 8) -St = 0.0. (d) Real part of the ρe /ρe ∞ eigenmode (α = 16.9 • , lower branch, β = 5.5) -St = 0.0.

Figure 16 :

 16 Figure 16: Flow visualisation of the 3-D unstable mode for several angle of attack on the two branches of the lift hysteresis (NACA4412, M = 0.2, Re = 350 000, δ = 0.0 • ).

Figure 17 :

 17 Figure 17: Effect of the sweep angle on the growth rate and frequency of the 3-D unstable mode (NACA4412, M = 0.2, Re = 350 000, α = 15.0 • , black: stability analysis, red: URANS simulations).

Figure 18 :

 18 Figure 18: Effect of the angle of attack on the spectra of the 2-D unstable mode (NACA4412, M = 0.2, Re = 350 000, δ = 0.0 • , β = 0).

Figure 19 :

 19 Figure 19: Flow visualisation of the (a,b): base-flow and (c,d): 2-D unstable mode on the two branches of the lift hysteresis (NACA4412, M = 0.2, Re = 350 000, α = 18.0 • , δ = 0.0 • , β = 0.0).

Figure 20 :

 20 Figure 20: Effect of the wavenumber on the 2-D wake instability mode (NACA4412, M = 0.2, Re = 350 000, α = 18.0 • ).

Figure 21 :

 21 Figure 21: Eigenspectrum in 2-D and 3-D (OALT25, M = 0.7352, Re = 3 × 10 6 , α = 4.0 • , δ = 0.0 • ).

Figure 22 :Figure 23 :

 2223 Figure 22: Effect of the wavenumber β on the 3-D mode (OALT25, M = 0.7352, Re = 3 × 10 6 , δ = 0.0 • )

Figure 24 :

 24 Figure 24: Effect of the sweep angle on the 3-D mode (OALT25, M = 0.7352, Re = 3 × 10 6 , α = 4.0 • , black: stability analysis, red: URANS simulations).

  (a) ρe /ρe ∞ base-flow. (b) Real part of the 2-D ρe /ρe ∞ eigenmode (β = 0.0) -St = 0.075. (c) Real part of the 3-D ρe /ρe ∞ eigenmode (β = 5.0) -St = 0.356.

Figure 25 :

 25 Figure 25: Flow visualisation of the (a): base-flow, (b,c) 2-D and 3-D unstable modes (OALT25, M = 0.7352, Re = 3 × 10 6 , α = 4.0 • , δ = 30 • ).

Figure 27 :

 27 Figure 27: Effect of the wavenumber on the 2-D buffet mode (OALT25, M = 0.7352, Re = 3 × 10 6 , α = 4.0 • ).

  (a) t = 68 (b) t = 85. (c) t = 102.

Figure 28 :

 28 Figure 28: Snapshots of the surface pressure coefficient and skin friction lines of the URANS simulations (NACA4412, M = 0.2, Re = 350 000, α = 15.0 • , δ = 0 • , L z = 6).

  Figure 29: Snapshots of the surface pressure coefficient and skin friction lines of the URANS simulations (OALT25, M = 0.7352, Re = 3 × 10 6 , α = 4.0 • , δ = 20 • , L z = 6).

  Figure 30: Base-flow Mach number field at several Mach number and angle of attack (NACA0012, Re = 10 7 , δ = 0 • )

Figure 31 :

 31 Figure 31: Effect of the Mach number and angle of attack on the growth rate of the 3-D mode (NACA0012, Re = 10 7 , δ = 0 • ).

  (a) M = 0.2 and α = 18.2 • -St = 0.0. (b) M = 0.3 and α = 15.0 • -St = 0.0. (c) M = 0.4 and α = 12.0 • -St = 0.0. (d) M = 0.5 and α = 7.9 • -St = 0.0. (e) M = 0.6 and α = 6.15 • -St = 0.0. (f) M = 0.72 and α = 3.5 • -St = 0.0.

Figure 32 :

 32 Figure 32: Flow visualisation of the 3-D ρe /ρe ∞ unstable mode (NACA0012, Re = 10 7 , δ = 0 • ).

  Subsonic (NACA0012, M = 0.2, Re = 10 7 , α = 18.5 • ). Transonic (NACA0012, M = 0.72, Re = 10 7 , α = 4.5 • ).

Figure 34 :

 34 Figure 34: Effect of the sweep angle on the 3-D mode frequency.

Figure 35 :

 35 Figure 35: Effect of the wavenumber of the 2-D buffet mode (OAT15A, Re = 3 × 10 6 , α = 3.5 • , δ = 0 • ).

Table 1 :

 1 Characteristics of the grids.

	.8)

Table 2 :

 2 Stall cells convection frequency for several wing spans and sweep angles (NACA4412, M = 0.2, Re = 350 000, α = 15 • ).

	00	4	1.50 4.20 0.00 0.000
	0 12.00	8	1.50 4.20 0.00 0.000
	10 6.00	4	1.50 4.20 0.14 0.033
	20 6.00	4	1.50 4.20 0.29 0.069
	30 6.00	4	1.50 4.2 0.46 0.110

† Email address for correspondence: frederic.plante@polymtl.ca

Acknowledgment

Parts of this paper were presented at the 54th 3AF International Conference AERO2019 Conference (FP63-AERO2019-plante). This work is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). Part of this work made use of the GENCI CINES facilities (Grant DARI No.A0042A10423).

Declaration of interests

The authors report no conflict of interest.

Appendix A. OAT15A analysis

As verification for the implementation of the method proposed by [START_REF] Schmid | Stability analysis for n -periodic arrays of fluid systems[END_REF], we have carried out a global linear stability analysis of transonic buffet of an OAT15A aerofoil in the same regime as in [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF]. A grid with the same property and the same numerical setup as the ones used for the OALT25 aerofoil are used. Fig. 35 shows the growth rate and the frequency of the 2-D mode with respect to β . The results of [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF] are also presented. One can observe that the frequencies and growth rates are close to the values found by [START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF]. The trend with respect to the wavenumber is also the same. This agreement verifies the method used in this paper.