I. Introduction

I commercial transonic aviation, the interaction between shock-waves and boundary-layers has significant consequences on aerothermodynamic performances and aircraft emissions, and represented one of the most important topics within the aeronautical scientific community for the past seventy years [START_REF] Dolling | Fifty years of shock-wave/boundary-layer interaction research: what next?[END_REF]. A specific type of shock-wave and turbulent boundary-layer interaction is the so-called transonic buffet, whose associated unsteady loads can originate wing vibrations, or buffeting. Degradation of performances, potential failure to fatigue and pilots and passengers' discomfort are some of the main repercussions of the buffeting conditions. To circumvent these complications, standard certification procedures impose the use of safety design margin that effectively limit the flight envelope. For both industrial and academic reasons, understanding, predicting and modeling transonic buffet are therefore of fundamental importance.

Turbulent transonic buffet consists of two distinct instabilities: a two-dimensional (2D) instability, related to the streamwise oscillations of the shock-wave forming on the wing suction side, and a three-dimensional (3D) one, characterized by the cross-flow propagation of the commonly named buffet cells.

Early investigations on turbulent transonic buffet mostly focused the attention on the 2D streamwise shock-oscillations.

The experiments of McDevitt and Okuno [START_REF] Mcdevitt | Static and dynamic pressure measurements on a NACA 0012 airfoil in the Ames High Reynolds Number Facility[END_REF] and Benoit and Legrain [START_REF] Benoit | Buffeting prediction for transport aircraft applications based on unsteady pressure measurements[END_REF], later confirmed by Jacquin et al. [START_REF] Jacquin | Experimental study of shock oscillation over a transonic supercritical profile[END_REF], identified characteristic Strouhal numbers ( , based on the free-stream velocity ∞ and chord length ) between 0.06 and 0.07.

Confirming the robustness of this 2D turbulent instability, the same frequencies were recovered when the boundary-layer over the laminar OALT25 airfoil was tripped to induce transition to turbulence [START_REF] Brion | Experimental analysis of the shock dynamics on a transonic laminar airfoil[END_REF]. From a numerical point of view, the literature is prevalently based on 2D Reynolds-Averaged Navier-Stokes equations (RANS) studies [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] Crouch | Origin of transonic buffet on aerofoils[END_REF][START_REF] Iovnovich | Reynolds-Averaged Navier-Stokes study of the shock-buffet instability mechanism[END_REF][START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF][START_REF] Giannelis | Influence of Mach number and angle of attack on the two-dimensional transonic buffet phenomenon[END_REF][START_REF] Paladini | Various approaches to determine active regions in an unstable global mode: application to transonic buffet[END_REF][START_REF] Plante | Simulation of transonic buffet using a time-spectral method[END_REF]. Although a considerable sensitivity to numerical schemes and turbulence models exists [START_REF] Goncalves | Turbulence model and numerical scheme assessment for buffet computations[END_REF], the 2D RANS calculations can recover the experimental frequencies associated with the chordwise shock-oscillations. The extensive reviews on turbulent transonic buffet [START_REF] Lee | Self-sustained shock oscillations on airfoils at transonic speeds[END_REF][START_REF] Giannelis | A review of recent developments in the understanding of transonic shock buffet[END_REF] reveal two widely accepted but conflicting interpretations of this 2D mechanism. On the one hand, these oscillations are believed to be the consequence of a feedback loop involving downstream traveling pressure waves emanated from the shock foot and upstream traveling pressure waves generated at the trailing edge that provide energy back to the shock [START_REF] Jacquin | Experimental study of shock oscillation over a transonic supercritical profile[END_REF][START_REF] Lee | Oscillatory shock motion caused by transonic shock boundary-layer interaction[END_REF]. On the other hand, the streamwise shock oscillations are interpreted as a global instability localized on the shock-foot that produces pressure perturbations in the direction normal to the profile and directly transfers energy along the shock [START_REF] Crouch | Origin of transonic buffet on aerofoils[END_REF][START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF]. Despite the large body of both experimental and numerical studies, a definite explanation is still lacking.

Due to the complexity of the phenomenon and thanks to modern experimental measurement techniques, more attention could only recently be given to describing the 3D aspects of buffet [START_REF] Dandois | Experimental study of transonic buffet phenomenon on a 3D swept wing[END_REF][START_REF] Koike | Unsteady pressure measurement of transonic buffet on NASA common research model[END_REF][START_REF] Sugioka | Experimental analysis of transonic buffet on a 3D swept wing using fast-response pressure-sensitive paint[END_REF][START_REF] Paladini | Analysis and comparison of transonic buffet phenomenon over several three-dimensional wings[END_REF][START_REF] Masini | Analysis of a civil aircraft wing transonic shock buffet experiment[END_REF]. Wing-body aircraft configurations showed that besides the characteristic frequencies of 2D shock-oscillations, a broadband energy content exists in the = 0.2 -0.6 range. This energy content has been associated with the convection of buffet cells towards the wingtip. With the increasing access to large computational resources, also the numerical interest could shift towards the investigation of fully 3D geometries. Made exception of a few large-eddy simulations (LES) [START_REF] Garnier | Large-Eddy Simulation of transonic buffet over a supercritical airfoil[END_REF][START_REF] Fukushima | Wall-modeled Large-Eddy Simulation of transonic airfoil buffet at high Reynolds number[END_REF] and direct numerical simulations [START_REF] Zauner | Direct Numerical Simulations of transonic flow around an airfoil at moderate Reynolds numbers[END_REF][START_REF] Zauner | Modal analysis of a laminar-flow airfoil under buffet conditions at Re = 500,000[END_REF][START_REF] Zauner | Wide domain simulations of flow over an unswept laminar wing section undergoing transonic buffet[END_REF], RANS or hybrid RANS/LES remain nonetheless the most popular approaches due to their relative inexpensiveness. Calculations on both 2.5D (extruded wings with spanwise periodicity) [START_REF] Deck | Numerical simulation of transonic buffet over a supercritical airfoil[END_REF][START_REF] Grossi | Prediction of transonic buffet by Delayed Detached-Eddy Simulation[END_REF][START_REF] Iovnovich | Numerical study of shock buffet on three-dimensional wings[END_REF][START_REF] Crouch | Global instability in the onset of transonic-wing buffet[END_REF][START_REF] Paladini | Transonic buffet instability: From two-dimensional airfoils to three-dimensional swept wings[END_REF][START_REF] Plante | Similarities between cellular patterns occurring in transonic buffet and subsonic stall[END_REF][START_REF] Deck | Towards an enhanced protection of attached boundary layers in hybrid RANS/LES methods[END_REF][START_REF] Plante | Link between subsonic stall and transonic buffet on swept and unswept wings: from global stability analysis to nonlinear dynamics[END_REF] and full-aircraft [START_REF] Sartor | Mach number effects on buffeting flow on a half wing-body configuration[END_REF][START_REF] Ishida | Transonic buffet simulation over supercritical airfoil by unsteady-FaSTAR code[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][START_REF] Timme | Global instability of wing shock-buffet onset[END_REF][START_REF] Masini | Scale-resolving simulations of a civil aircraft wing transonic shock-buffet experiment[END_REF][START_REF] Ehrle | Simulation of transonic bufet with an automated zonal DES approach[END_REF] configurations are proliferating and, although a strong sensitivity to the numerical grid and turbulence modeling prevents satisfactory quantitative comparisons, new light is being shed on 3D buffet.

While the mechanisms related to both instabilities need further scrutiny, it is worth to bring attention to the associated dynamics. By considering a 2D flow and using a RANS approach to model turbulence, the present contribution aims at showing that, despite the complexity of the physical mechanisms coming into play, the dynamics of transonic buffet within such a framework is relatively simple. For this purpose, techniques from system identification can be used to obtain low-dimensional models that describe the full-state dynamics in a computationally efficient way.

Reduced-order modeling (ROM) is nowadays applied to a wide spectrum of scientific and industrial applications as well as to economic processes. The techniques can be very diverse and range from physical reductions [START_REF] Noack | Reduced-Order Modelling for Flow Control[END_REF] to data-driven algorithms [START_REF] Brunton | Closed-loop turbulence control: progress and challenges[END_REF]. The limitations of the most commonly used approaches, like those based on proper orthogonal decomposition(POD)-Galerkin methods [START_REF] Noack | Reduced-Order Modelling for Flow Control[END_REF][START_REF] Carlberg | Preserving Lagrangian structure in nonlinear model reduction with application to structural dynamics[END_REF] or dynamic mode decomposition (DMD) [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF], pushed the scientific community to develop new models applicable to a wider range of flow conditions. In the present work, a system identification method is trained on a dataset consisting of the lift coefficient time-series, ( ), obtained at different angles of attack (AoA). The identification method makes use of the sparse identification of nonlinear dynamics (SINDy) [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Loiseau | Sparse reduced-order modelling: sensor-based dynamics to full-state estimation[END_REF] to identify a parametrized minimal-order model that best describes the system dynamics over the range of operating conditions considered. These models are then supplemented with an equally parametrized DMD representation of the instantaneous state vector to enable on-the-fly estimation of the whole flow field around the airfoil directly from its lift measurements.

The paper is organized as follows: details of the 2D RANS simulations performed to calculate the lift coefficient timehistories along with a description of the system identification and DMD methods are given in Sec. II. A comprehensive description of the training dataset and the results obtained by system identification are presented in Sec. III. The identified system is tested over a case not present in the training data set in Sec. IV, and the physical implications of the findings are also discussed. Conclusions and perspectives on the future work are summarized in Sec. V.

II. Numerical Methodology

This section presents the numerical details of the unsteady RANS simulations used to obtain the training data set consisting of the time series of the lift coefficient at various angles of attack. In a second time, the system identification framework, consisting of the joint use of a preprocessing of the lift coefficient time-series using the Broomhead-King embedding technique and the SINDy method for the subsequent model identification, is presented. Finally, a brief recall on DMD used for the full-state estimation is given.

A. 2D Steady and Unsteady RANS Simulations

For a compressible perfect gas, the 2D unsteady Reynolds-Averaged Navier-Stokes equations can be written as

q = (q), (1)
where is the differential nonlinear Navier-Stokes operator in Cartesian coordinates and q = [ , u, , ] is the conservative quantities vector with, , u, and being the fluid density, velocity vector, total energy and turbulent viscosity, respectively. All numerical simulations in this paper are run with an in-house code solving the compressible Navier-Stokes equations on multi-block fully-parallelized structured grids in a finite-volume fashion. While the mean convective fluxes are calculated using a Jameson-Schmidt-Turkel third-order scheme with artificial dissipation, a second-order scheme with Harten's correction is used for the turbulent ones. All viscous terms are differentiated with a second-order centered scheme. The temporal discretization is based on a dual time-stepping method [START_REF] Jameson | Time-dependent calculations using multigrid with applications to unsteady flows past airfoils and wings[END_REF] using a second-order extrapolation for the derivatives with respect to the physical time. The boundary conditions used for both the steady or unsteady nonlinear calculations are: no-slip velocity, adiabatic temperature, and pressure extrapolation on the profile walls; the AoA is applied on the imposed uniform velocity at the inflow of the numerical domain; characteristic boundary conditions are set at the domain lateral boundaries and outflow to minimize wave reflections.

The Spalart-Allmaras turbulence model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] with Edwards modification [START_REF] Edwards | Comparison of Eddy viscosity-transport turbulence models for three-dimensional, shockseparated flowfields[END_REF] and quadratic constitutive relation (QCR)

2000 version [START_REF] Spalart | Strategies for turbulence modelling ands imulation[END_REF] is used to close the RANS equations. Each simulation is started with the solution of the steady RANS equations and evolved in time until it reaches a fully-established nonlinear regime. During this whole evolution, the lift coefficient is monitored and will serve as an input for the system identification task.

B. Feature extraction

This work aims at identifying low-order models of the transonic buffet using solely the time-evolution of the lift coefficient and instantaneous flow fields obtained from URANS simulations. To do so, the measurement first needs to be augmented into a higher-dimensional latent space where the dynamics is easy to model. Leveraging Taken's embedding theorem [START_REF] Takens | Detecting strange attractors in turbulence[END_REF], such an embedding can be obtained using the algorithm proposed by Broomhead and King [START_REF] Broomhead | Extracting qualitative dynamics from experimental data[END_REF]. For each angle of attack considered herein, a × ( -) Hankel matrix is first constructed based on the recorded evolution of the lift coefficient, i.e.

H =                 0 1 2 • • • - 1 2 3 • • • -+1 . . . . . . . . . . . . . . . +1 +2 • • •                 (2)
where = ( Δ ) is the value of the lift coefficient recorded at the th time step and is the length of the recorded signal. The number of time delays considered to construct this Hankel matrix defines the embedding window. It is usually taken to be of the order of the dominant time-scale in the signal.

As a second step, this Hankel matrix is decomposed via its economy-sized singular value decomposition, i.e.

H = UΣV

with U ∈ R × a unitary matrix, Σ ∈ R × a diagonal matrix whose entries are the singular values of H and V ∈ R ( -)× .

From a physical point of view, the th column of U can be understood as a temporal convolution kernel extracting specific features of the original signal. The singular value then characterizes how important this feature is to reconstruct the original signal while the th column of V corresponds to the (scaled) signal resulting from this convolution. As for principal component analysis, the optimal dimension of the desired latent space can be inferred from the distribution of the singular values . Because of the close connection of this time-delay embedding with principal component analysis, the signals in V are linearly uncorrelated and sometimes called principal component trajectories. For more details about this decomposition, interested readers are referred to the recent work of Kamb et al. [START_REF] Kamb | Time-delay observables for Koopman: Theory and applications[END_REF].

C. SINDy -System Identification

Having lifted the single measurement signal into a higher-dimensional latent space, a model of the dynamics in this latent space can be obtained using system identification techniques. To do so, the sparse identification of nonlinear dynamics (SINDy) framework originally proposed by Brunton et al. [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF] is considered. For the sake of simplicity, the presentation of the method is restricted to the case of a 2D state vector a( ) = [ 1 ( ), 2 ( )] (as is the case in the rest of this paper) albeit the method readily generalizes to higher-dimensional systems. SINDy relies on the observation that, for many systems modeled as da d

= f (a) (4)
the unknown function f (a) : R → R is sparse in the space of possible right-hand side functions.

To identify f (a), time-series data is first collected and formed into the data matrix

A = a( 1 ) a( 2 ) • • • a( ) . (5)
A similar matrix A of time derivatives is formed. As a second step, a possibly over-complete library (or dictionary) Θ of candidate nonlinear functions is constructed, e.g.

Θ(A) = 1 A 2 (A) • • • (A) . (6)
Here (A) denotes a matrix with columns vectors given by all possible time-series of th degree monomials in the state a, e.g.

2 ( 1 , 2 ) = 2 1 1 2 2 2 . ( 7
)
Note that any basis function may be included in the library Θ albeit polynomials have proven to work well for fluid dynamics [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Loiseau | Sparse reduced-order modelling: sensor-based dynamics to full-state estimation[END_REF][START_REF] Loiseau | Constrained sparse Galerkin regression[END_REF][START_REF] Guan | Sparse nonlinear models of chaotic electroconvection[END_REF][START_REF] Loiseau | Data-driven modeling of the chaotic thermal convection in an annular thermosyphon[END_REF]. The unknown dynamical system may now be represented in terms of the data matrices as

A = Θ(A)Ξ. (8)
Each column Ξ is a vector of coefficients determining the active terms in the equation governing the dynamics of ( ).

A parsimonious model will provide an accurate model fit with as few non-zero terms as possible in Ξ. Identifying this sparsity pattern may be formulated as the following optimization problem minimize

Ξ card(Ξ ) subject to A -Θ(A)Ξ 2 2 ≤ (9)
where card(Ξ ) = Ξ 0 is the cardinality (or ℓ 0 pseudo-norm) of Ξ , i.e. its number of non-zero entries, and the constraint A -Θ(A)Ξ 2 2 ≤ quantifies the fidelity of the model with respect to the data. The above minimization problem is however a combinatorially complex optimization problem. Various convex relaxations to this problem have been proposed over the years, most of them replacing the ℓ 0 with the ℓ 1 norm such as in LASSO regression. In this work, the simple hard-thresholded least-squares algorithm proposed in [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF] was however found to be sufficient. Note finally that, since its introduction, numerous extensions of SINDy have been proposed, see for instance [START_REF] Loiseau | Constrained sparse Galerkin regression[END_REF][START_REF] Mangan | Inferring biological networks by sparse identification of nonlinear dynamics[END_REF][START_REF] Kaiser | Sparse identification of nonlinear dynamics for model predictive control in the low-data limit[END_REF][START_REF] Rudy | Data-driven discovery of partial differential equations[END_REF][START_REF] Mangan | Model selection for dynamical systems via sparse regression and information criteria[END_REF][START_REF] Brunton | Chaos as an intermittently forced linear system[END_REF][START_REF] Callaham | Nonlinear stochastic modeling with Langevin regression[END_REF][START_REF] Kaheman | SINDy-PI: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics[END_REF][START_REF] Kaheman | Automatic differentiation to simultaneously identify nonlinear dynamics and extract noise probability distributions from data[END_REF][START_REF] Kaptanoglu | Promoting global stability in data-driven models of quadratic nonlinear dynamics[END_REF][START_REF] Cortiella | Sparse identification of nonlinear dynamical systems via reweighted ℓ 1 -regularized least squares[END_REF] and references therein.

D. Dynamic Mode Decomposition

Fluid flows are characterized by high-dimensional nonlinear dynamics giving rise to rich structures. Despite this apparent complexity, the dynamics often evolve on a low-dimensional attractor defined by a few dominant coherent structures [START_REF] Holmes | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF]. Over the years, being able to leverage these coherent structures has proven instrumental for real-time control or estimation of these complex spatio-temporal systems. Since the introduction of POD in the fluid dynamics community by Lumley [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF] in the late sixties, various other such techniques have been proposed to extract these coherent structures. In the present work, DMD is selected and here briefly discussed.

Given a sequence of evenly sampled high-dimensional snapshots { ( , )} =1, related by +1 = ( ) where = ( , ) and : R → R (with 1), DMD aims at finding a low-rank operator ∈ R × that best approximates the unknown function in the least-squares sense. While the original idea has been proposed by Rowley et al. [START_REF] Rowley | Spectral analysis of nonlinear flows[END_REF] and Schmid [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF], the well-posed optimization problem leading to the optimal identification of this low-rank operator has only been recently put in the limelight by Héas and Herzet [START_REF] Héas | Low-Rank Dynamic Mode Decomposition: Optimal Solution in Polynomial-Time[END_REF] and the proof of its optimality has been simplified in Loiseau [START_REF] Loiseau | Data-driven modeling of the chaotic thermal convection in an annular thermosyphon[END_REF]. Introducing the data matrix

= 1 2 • • • -1
and its time-shifted counterpart where is the × identity matrix. As shown in [START_REF] Héas | Low-Rank Dynamic Mode Decomposition: Optimal Solution in Polynomial-Time[END_REF] and [START_REF] Loiseau | Data-driven modeling of the chaotic thermal convection in an annular thermosyphon[END_REF], the matrix solution to this problem is formed by leading eigenvectors of the following eigenvalue problem

= -1
while is solution to

= -1
where = is covariance of and is the cross-covariance matrix. Once and have been computed, the low-rank DMD operator can easily be factorized as = where , ∈ C × are the left and right DMD eigenmodes and ∈ C × is the DMD eigenvalues matrix (i.e.

= ). These eigenmodes and associated eigenvalues may provide valuable insights into the dynamics of the spatio-temporal coherent structures identified. Alternatively, the low-rank DMD operator may also be factorized using its singular value decomposition = .

It must be noted that other variants of DMD have been proposed in the literature [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF][START_REF] Rowley | Spectral analysis of nonlinear flows[END_REF][START_REF] Tu | On dynamic mode decomposition: Theory and applications[END_REF][START_REF] Jovanović | Sparsity-promoting dynamic mode decomposition[END_REF][START_REF] Williams | A Data-Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition[END_REF][START_REF] Kutz | Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems[END_REF][START_REF] Dawson | Characterizing and correcting for the effect of sensor noise in the dynamic mode decomposition[END_REF][START_REF] Noack | Recursive dynamic mode decomposition of transient and post-transient wake flows[END_REF][START_REF] Le Clainche | Higher Order Dynamic Mode Decomposition[END_REF][START_REF] Hirsh | Centering Data Improves the Dynamic Mode Decomposition[END_REF], each with its pros and cons. To the best of the authors' knowledge, however, none of these techniques can easily be recast in this seemingly natural framework of reduced-rank regressions.

III. Results

Steady and unsteady RANS solutions forming the training set for the system identification are described in detail in Sec. III.A along with some comparisons with the literature. Results on the system identification are instead discussed in Sec. III.B.

A. Training Dataset: Unsteady RANS Lift Coefficient Time-Histories

The training dataset consists of the lift coefficient time-histories at different AoA in the transonic buffet regime. For each AoA, the lift coefficient time evolution covers the whole transient and post-transients dynamics, from the steady solution to the fully-developed buffet shock-oscillations. Thus, the results concerning both the steady and unsteady calculations are presented in this section and compared with some studies present in the literature.

Numerical Setup

RANS Solutions

The pseudo-unsteady temporal integration with a Courant-Friedrichs-Lewy (CFL) number equal to 10 is used to

Unsteady RANS Solutions

For each AoA, an unsteady RANS calculation is started from its corresponding steady solution with a physical time step fixed to Δ * = 6.85 × 10 -7 . For the AoA at which transonic buffet occurs, the time evolution of the lift coefficient is recorded from the steady solution up to non-linear saturation, when the shock-oscillations due to buffet are fully developed. In agreement with Sartor et al. [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF], buffet onset occurs around AoA = 3.5 o . Some differences exist in terms of buffet offset (AoA > 5.00 o in the present study against AoA > 6.25 o in Sartor et al. [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF]) but some scattering in the literature [START_REF] Giannelis | Influence of Mach number and angle of attack on the two-dimensional transonic buffet phenomenon[END_REF] suggests that cases at high AoA may present separations too large to be correctly simulated in a RANS framework. The lift coefficient time evolution for AoA = [3.50 o , 4.75 o ] and 0.25 o step is reported in figure 4, where it is possible to appreciate the exponential instability as the flow transitions from a steady solution to a limit cycle. In agreement with Sartor et al. [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF], for increasing incidence, the buffet frequency does not vary significantly (between 74 and 77 ), the mean (time-averaged) value decreases, and the oscillation amplitude increases. This ensemble of lift coefficient time-histories is selected as a database for the training of the system identification method presented in the following sections. An important feature to discuss concerning the evolution of the lift coefficient in the case of transonic buffet is the "asymmetric" nature of the oscillations around its corresponding steady value. The baseflow-subtracted lift coefficient evolution normalized with the maximum of the oscillations amplitude, 1 , is plotted against time and its temporal most upstream shock position (blue), and most downstream shock position (red). The phase plot shows that the limit cycle has a non-circular shape because the shock slows down when reaching its most upstream position, causing the high-lift phase to be at a higher absolute value with respect to the low-lift counterpart. The nearly steady (left plot), most upstream position (middle plot), and most downstream position (right plot) states are reported in figure 6, where the sonic (black solid line) and zero-streamwise velocity (white solid line) iso-lines are superimposed to the contours of the streamwise velocity.

DMD Analysis

For each AoA, 512 snapshots are collected at saturation. About 16 snapshots per shock-oscillation period were sampled, for a total of 32 periods. The snapshots are 2D flow fields consisting of density, velocities, energy, and turbulent viscosity for each domain cell.

While some differences exist between each AoA, the DMD results are qualitatively similar and only those for AoA = 4.75 o are shown here for illustrative purposes. Figure 7 (left plot) reports the imaginary and real part of the DMD eigenvalues. The four eigenvalues, and their complex conjugates, resting on the unit circle (dashed gray line) are indicated by black full symbols and the streamwise velocity contours of the corresponding eigenfunctions are reported on the right part of the figure (the correspondence between eigenvalues and eigenfunctions is denoted by the letters indicated). While the top left subplot in figure 7 corresponds to DMD mode at the buffet frequency, the remaining plots are simply its higher harmonics. The eigenfunction is localized on the shock and shear-layer regions, as already seen in the DMD analysis by Poplingher et al. [START_REF] Poplingher | Modal analysis of transonic shock buffet on 2D airfoil[END_REF]. It is interesting to emphasize that a similar eigenmode structure is given by linear global stability analysis [START_REF] Crouch | Origin of transonic buffet on aerofoils[END_REF][START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF], suggesting that the mechanisms that govern buffet are fundamentally linear. While linear global stability analysis is performed around a RANS (steady base flow solution) and the eigen-perturbations are concentrated on the shock, the DMD analysis is done around the mean flow and the eigenfunction appears to be smeared around the shock.

B. System Identification

Given the time series of the lift coefficient at various AoA, the aim is now to identify a relatively simple dynamical system capable of describing accurately the observations made in the previous section, namely the exponential instability of the base flow and the subsequent transition to a limit cycle, as well as the asymmetry in the oscillation (i.e. high-lift vs low-lift).

Latent space

As discussed in Sec. II.B, time-series of the lift coefficient obtained from URANS simulations are augmented lifted into a higher-dimensional latent space using the method proposed by Broomhead and King [START_REF] Broomhead | Extracting qualitative dynamics from experimental data[END_REF]. For that purpose, a window of Δ ≈ 0.013 is considered, i.e. approximately one period of oscillation. Figure 8 depicts the corresponding distribution of singular values, with all cases superimposed. Evidently, these Hankel matrices can be well approximated by a rank 2 model, consistent with the fact that, for the range of parameters considered, the flow is characterized by an oscillatory behavior. Figure 9 shows the convolution kernels learned by the method while figure 10 presents the trajectory of the system in the corresponding 2D latent space for the different angles of attack. Looking at the latter figure, it can be seen that the oscillatory behavior of the transonic buffet is nicely captured by this embedding.

Additionally, comparing these embeddings with the naïve one in figure 5, it can be seen that azimuthal invariance is restored in this latent space. In the next section, we will leverage this azimuthal invariance to obtain a simple dynamical model expressed in terms of the oscillation's amplitude ( ) and its phase ( ). are consistent with recent results available in the literature [START_REF] Crouch | Origin of transonic buffet on aerofoils[END_REF][START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF].

Estimating the lift coefficient from the latent space

The dynamical systems identified in the previous section can be augmented with a measurement equation relating the instantaneous state of the system in the latent space to the corresponding lift measurement. Based on physical arguments, this measurement equation can be expressed as The constant term 0 models the contribution to the lift of the linearly unstable base flow while that of the oscillatory dynamics is captured by 2 cos( ). The term 3 2 cos(2 ) captures the contributions resulting from nonlinear interactions giving rise to the buffet's second harmonic. Finally, the term proportional to is included as being the simplest term that could break the symmetry of the reconstructed signal. Including this term is of crucial importance to capture the asymmetry between the high-lift and low-lift phases clearly visible in figure 5. These coefficients are estimated using a simple linear least-squares approach for all of the angles of attack considered herein. The identified numerical values of these coefficients are reported in table 2. Given the true time-series of ( ) and ( ) for an angle of attack of 3.75 o , figure 13 (top plot) compares the true measurements of the lift coefficient against those reconstructed from the measurement equation Eq. [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF]. Similar results are reported in figure 13 (bottom plot) for an angle of attack of 4.75 o for which capturing the asymmetry of the lift signal can be of utmost importance. In both cases, excellent agreement between the true time series and the estimated one is obtained.

( , ) = 0 + 1 + 2 cos( ) + 3 2 cos(2 ). (10)
It should be noted finally that the evolution of these parameters with respect to the angle of attack is monotonic. As for the parameters of the identified dynamical systems, the dependence of the measurement equation's parameters with respect to the angle of attack can easily be captured by a simple interpolation scheme such as cubic splines.

Estimating the full state vector from the latent space

The last step in the data-driven analysis is to reconstruct the high-dimensional state vector from the 2D latent space presented in III.B.1. Various approaches have been proposed in the literature to perform this task, from a simple locally linear embedding [START_REF] Loiseau | Sparse reduced-order modelling: sensor-based dynamics to full-state estimation[END_REF] to advanced decoder networks [START_REF] Nair | Leveraging reduced-order models for state estimation using deep learning[END_REF]. In this work, a simple map will be identified to relate the low-dimensional latent state vector to the representation of the flow in the DMD basis.

The projection of the k th snapshot onto the span of the DMD modes is given by

( Δ ) = † ( Δ )
where ( Δ ) is the state vector from the URANS simulation, † is the Moore-Penrose pseudoinverse of the DMD eigenbasis and ∈ C 4 is the complex-valued amplitude of the four complex-conjugate pairs of DMD modes retained.

As before, this complex-valued nature is leveraged to rewrite each entry of in polar coordinates, i.e.

= exp( ) ∀ = 1, • • • , 4.
The objective is now to find the map from the polar representation of the latent space state vector to that of the DMD representation . Since the dynamics of the flow in the saturated stage are strictly periodic, one can easily show that the simplest map is given by = = ± + ( )

0
where and are the coordinates of the latent space state vector in polar coordinates, and and ( ) 0 are the map's parameters for the th DMD mode. These parameters are once again identified using a simple linear least-squares procedure and are reported in table 3. The predicted evolution of the DMD modes is compared against the ground truth in figure 14. Only the configuration at an angle of attack of 4 o is presented for the sake of conciseness albeit similar results have been obtained for all other cases. Excellent agreement is obtained for the first two pairs of complex conjugate DMD modes while the agreement slightly deteriorates for the higher-order ones (possibly due to errors resulting from the sampling of the dynamics used to compute these DMD modes). The orbits in these phase portraits are reminiscent of Lissajous curves further highlighting that the higher-order DMD modes are simply harmonics of the fundamental one.

It needs to be noted finally that this DMD-based estimation of the full velocity field provides excellent agreement with the ground truth data only once the flow has reached its limit cycle. During the initial stages of transition, the spatial support of the perturbation would be better captured by the instability modes [START_REF] Crouch | Origin of transonic buffet on aerofoils[END_REF][START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF] which can differ quite substantially from the DMD modes. Using arguments from dynamical systems theory, one can however argue that these 3 4 2.12 0.65 0.34 0.23 instability modes are continuously deformed into the DMD modes as their amplitude increases. As such, one could use the technique proposed in [START_REF] Loiseau | 9 From the POD-Galerkin method to sparse manifold models[END_REF] wherein a parameterized reduced basis is constructed to continuously transform the instability modes to the POD/DMD modes based on interpolation on the Grassman matrix manifold [START_REF] Amsallem | Interpolation on manifolds of CFD-based fluid and finite element-based structural reduced-order models for on-line aeroelastic predictions[END_REF][START_REF] Amsallem | Interpolation Method for Adapting Reduced-Order Models and Application to Aeroelasticity[END_REF]. A similar procedure can be used to interpolate the DMD modes as a function of the angle of attack.

IV. Discussions

A. Identified system predictive capabilities and computational cost

To test the predictive capabilities of the identified system, an angle of attack of 4.15 o , not included in the training set, has been selected to compare the time-histories of true and predicted lift coefficients. As previously done in section III.A, a RANS solution is first obtained and used as initial condition for a URANS calculation to generate the ground truth data.

Starting from the steady solution, the flow undergoes a linear transient and eventually saturates nonlinearly towards a limit cycle characterized by periodic self-sustained shock oscillations.

The parameters of both the dynamical system and measurement equation identified in the previous section are interpolated at the corresponding angle of attack using cubic splines. The resulting models are then used to generate a prediction of the lift coefficient time-history with no additional information about the flow. The Stuart-Landau oscillator model is first integrated forward in time starting from an initial condition close to its fixed point to generate the time series of the latent space vector. This 2D time series is then processed by the measurement equation to generate the estimated evolution of the true lift coefficient. A comparison of the true (solid black line) and predicted (gray dashed While this phase difference is unavoidable due to the arbitrary choice of the initial condition, the shock oscillation amplitudes, the linear transient growth rate (bottom left figure), and the signal frequency energy contents (bottom right figure) remarkably match. Although the identified system is in general able to reproduce the asymmetry of the shock oscillations in the nonlinearly saturated stage (see figure 13), in this case, the asymmetry is minimal. Note moreover that, since only up to quadratic terms were included in the equation [START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF] to account for nonlinear interactions, the predicted power spectral density (PSD) is not able to reproduce third or higher harmonics. For the prediction of the lift coefficient time history, the contribution of these higher harmonics is however negligible.

Although not presented for the sake of conciseness, the DMD modes can be equally interpolated at this particular angle of attack (along with an interpolation of the mean flow) and the simple map identified in section III.B.4 can be used to obtain an estimation of the whole velocity field around the airfoil.

One final thought concerns the computational effort required to obtained both the true and predicted lift coefficient time histories. While the cost of the CFD simulation is of the order of a few hundred core hours, the identified model generates the predicted lift on a single process in a matter of seconds with remarkable accuracy.

B. Physical implications

The physical mechanisms governing the turbulent transonic buffet are at the center of many academic and industrial discussions. Being high-speed buffet a limiting cause of the flight envelope, accurate predictions are of utmost importance during design, certification, and operational phases. Even though a better understanding of the physical mechanisms underlying transonic buffet remains critical, this work highlights the duality of the physical interpretation of the governing mechanisms versus the dynamical behavior of interest for engineering applications. Despite the complexity of the problem, where turbulent boundary layers, separation, and unsteady shock waves coexist, the dynamics of the system (as modeled in a URANS framework) is very simple. The identified system is in fact no different than what was found for the canonical 2D cylinder flow at low Reynolds numbers [START_REF] Loiseau | Sparse reduced-order modelling: sensor-based dynamics to full-state estimation[END_REF]. To the authors' knowledge, the model proposed herein is currently the simplest one able to accurately predict 2D transonic buffet.

One might argue that many aspects of high Reynolds number turbulence are neglected in the current 2D RANS framework. While this is certainly true, and the consideration of 3D effects and corresponding appearance of buffet cells [START_REF] Dandois | Experimental study of transonic buffet phenomenon on a 3D swept wing[END_REF][START_REF] Koike | Unsteady pressure measurement of transonic buffet on NASA common research model[END_REF][START_REF] Sugioka | Experimental analysis of transonic buffet on a 3D swept wing using fast-response pressure-sensitive paint[END_REF][START_REF] Paladini | Analysis and comparison of transonic buffet phenomenon over several three-dimensional wings[END_REF][START_REF] Masini | Analysis of a civil aircraft wing transonic shock buffet experiment[END_REF] further complexifies the dynamics of the flow, experimental and higher fidelity numerical studies have shown that only two instabilities are dominant. At buffeting conditions, these instabilities are also robust over a relatively large range of Mach and Reynolds numbers. A high-frequency instability (3D spanwise traveling buffet cells) superimposed onto a low-frequency one (2D streamwise shock oscillations) are the main characters at play. In particular, looking at the time-series generated by 3D LES computations (i.e., see figure 2a in [START_REF] Zauner | Modal analysis of a laminar-flow airfoil under buffet conditions at Re = 500,000[END_REF]), it can be speculated that the dynamics of the lift coefficient could be modeled using the same simple Stuart-Landau oscillator augmented with a second component to account for high-frequency buffet cells and a third stochastic one for the small scale turbulent fluctuations. Such a modeling strategy is currently undertaken by Callaham et al. [START_REF] Callaham | An empirical mean-field model of symmetry-breaking in a turbulent wake[END_REF] to model the temporal evolution of the center of pressure in an experimental turbulent wake using an extension of SINDy for stochastic systems [START_REF] Callaham | Nonlinear stochastic modeling with Langevin regression[END_REF]. If these speculations are confirmed by further investigations, the use of these simple models could represent a very powerful tool for the assessment of high-speed buffet in the flight envelope. For a given geometry, only a relatively modest number of numerical simulations or experiments would be needed to adjust the coefficients of the identified system and have an accurate prediction in the complete buffet regime.

V. Conclusions

The sparse identification of nonlinear dynamics (SINDy) method is used to obtain a low-dimensional model to describe the nonlinear dynamics of a 2D supercritical profile at transonic buffet conditions within a 2D RANS framework.

A simple parametrized Stuart-Landau oscillator model is identified and used to reconstruct lift time histories and full states at different angles of attack within the buffet condition range. The identified model is tested on one case not included in the training dataset and the predicted lift time-history from the RANS solution to nonlinear saturation (attainment of the periodic shock oscillations in the limit cycle) agrees remarkably with the corresponding CFD simulation.

The findings presented here offer a few points for reflection. While the scientific efforts on the elucidation of the buffet governing mechanisms should not diminish, the modeling of the associated flow response might prove to be relatively less complex. The dynamics identified in this work (and the shared similarities with 3D wings and full aircraft configurations) suggest that an accurate representation of the flow response might be achieved with simple models.

Additionally, the accuracy of the model and its rapidity in reproducing the buffet dynamics make the identified system an attractive and powerful tool for real-time simulations in the context of design, certification and operational phases in the life cycle of commercial aircrafts.

While further investigations should be carried out by considering 3D configurations and different levels of CFD high-fidelity datasets, the present work poses the basis for accurate real-time buffet dynamics prediction for both industrial and academic purposes.

=

  shock-wave dynamics on a two-dimensional aerofoil at transonic buffet conditions is characterized by time-periodic oscillations. Given the time series of the lift coefficient at different angles of attack for the OAT15A supercritical profile, the sparse identification of nonlinear dynamics (SINDy) technique is used to extract a parametrized, interpretable and minimal-order description of this dynamics. For all of the operating conditions considered, SINDy infers that the dynamics in the lift coefficient time series can be modeled by a simple parametrized Stuart-Landau oscillator, reducing the computation time from hundreds of core hours to seconds. The identified models are then supplemented with equally parametrized measurement equations and low-rank DMD representation of the instantaneous state vector to reconstruct the true lift signal and enable real-time estimation of the whole flow field. Simplicity, accuracy and interpretability make the identified model a very attractive tool towards the construction of real-time systems to be used during the design, certification and operational phases of the aircraft life cycle. Nomenclature a = [ 1 , 2 ] = features space vector; features space components AoA = angle of attack = airfoil chord = parameters of the identified lift measurement equation q = [ , , , , ] = conservative variables vector = time u = [ , ] = velocity vector; streamwise and wall-normal velocity components = complex DMD mode amplitudes ,

  =where the rank of the desired approximation. Because of the rank constraint, this optimization problem is highly non-convex. It admits nonetheless a closed-form solution. Factorizing the unknown low-rank operator as = with and ∈ R × two rank-matrices, this optimization problem can be recast as minimize ,

  Eleven simulations have been carried out by spanning the AoA from 3.00 o to 5.50 o by an increment of 0.25 o . The flow conditions are set to correspond to those in Sartor et al. [10], where an OAT15A airfoil with a sharp trailing edge and chord of = 0.23 was considered. The same geometry used in Sartor et al. [10] was kindly provided by Dr. Sartor, ONERA. The stagnation pressure and temperature are 101325 and 300 , respectively. The Mach number is = 0.73 and the Reynolds number based on the chord length is = 3.2 × 10 6 .

Fig. 1 Fig. 2

 12 Fig. 1 Numerical grid used for the steady/unsteady RANS calculations. Whole domain showing the C-type grid topology and far-field boundaries (left) and zoom around the airfoil (right). The airfoil section is colored in red in both plots.

  filter the unsteadiness. Steady solutions are converged until the residuals of the state variables in the ℓ 2 -norm are lower than 10 -8 . Distribution of the pressure coefficient for the AoA = 3.50, 4.50 and 5.50 o are compared against those of Deck [27] and [10] in figure 2 (left plot) and show good agreement on both airfoil sides and in terms of shock positions. Contours of the dimensional streamwise velocity for the steady solution at AoA = 4.50 o are plotted in figure 2 (right plot) along with the sonic (black solid line) and zero-streamwise velocity (white solid line) iso-lines, showing the supersonic flow region, shock position, and separation, respectively. Figure 3 shows the pressure, , and skin-friction, , coefficient distributions for all the AoA examined. Some transparency has been added to the monochromatic lines representing the distributions. Hence, regions where the pressure and skin friction distributions become darker indicate a weak dependency on the incidence angle. The distribution on the upper side of the airfoil is flat upstream of the shock position and very weak dependency on the AoA can be seen on the lower airfoil surface. The shows that the region between normal shock and trailing edge is fully separated for all cases examined. The separation point moves upstream for increasing AoA and, for high incidence (above AoA = 4.25 o ) the distribution is nearly independent of the AoA for / > 0.5.

Fig. 3

 3 Fig. 3 RANS pressure (left plot) and skin-friction (right plot) coefficient distributions for all cases between AoA = 3.00 : 0.25 : 5.50 o . The arrows indicate the evolution of the AoA.

Fig. 4

 4 Fig. 4 Lift coefficient time-histories at different AoA constituting the training dataset for the system identification method.

Fig. 5

 5 Fig. 5 Baseflow-subtracted lift coefficient normalized with the maximum of the oscillation amplitude against time (left) and its time derivative (right) for the AoA = 4.75 o . The colored dots correspond to different flow states: nearly-steady state (green), transient (magenta), most upstream shock position (blue) and most downstream shock position (red).

Fig. 6

 6 Fig. 6 Sonic (black solid line) and zero-streamwise velocity (white solid line) iso-lines superimposed to the dimensional streamwise-velocity contours for the nearly steady (left plot), most upstream position (middle plot) and most downstream position (right plot) states at AoA = 4.75 o . The different colored dots correspond to the same states indicated in figure 5.

Fig. 7

 7 Fig. 7 DMD eigenvalues and corresponding streamwise eigen-velocity contours (indicated by the letters) for the AoA = 4.75 o case.

Fig. 8 Fig. 9

 89 Fig. 8 Singular value distribution of the Hankel matrix. For all cases, these distributions indicate that a 2D latent space is sufficient.

Fig. 10 Fig. 11 Fig. 12

 101112 Fig. 10 Trajectories in the latent space. Only a subset of the different AoA considered herein is shown.

Fig. 13 Table 2

 132 Fig. 13 Comparison of the ground truth and modeled evolution of the lift coefficient time history for AoA = 3.75 o (top plot) and 4.75 o (bottom plot).

3

 3 Fig. 14 Projection of the dynamics onto the span of the first four DMD eigenmodes. The black lines represent the projection of the true snapshot onto these modes while the light gray one depicts the evolution directly inferred from the latent state vector.

Fig. 15

 15 Fig. 15 Comparisons between the true (black solid line) and predicted (gray dashed line) lift coefficients for AoA=4.15 o .

Table 1 Parameters of the models identified using SINDy for the various angles of attack considered.

 1

	AoA	3.5 o	3.75 o	4 o	4.25 o	4.5 o	4.75 o
		4.34	18.3	25.93	22.24	17.71	11.55
		377.23 264.86 205.68 129.86 86.05	52.35
		465.78 465.62 466.31 466.60 467.30 469.34

Dynamical modeling using SINDy

As stated, the key features of the transonic buffet's dynamics (when modeled using URANS) can be captured by a 2D latent space. Moreover, in this latent space, figure 10 shows that the corresponding phase portraits are invariant in the azimuthal direction. In order to identify a simple model, we can thus work in the polar coordinates

where the phase ( ), initially in [-, [ is additionally unwrapped. Figure 11 shows the corresponding time-series (in black) that will be used for the system identification.

The oscillation's amplitude being independent of its phase, one can identify a low-order model independently for both of them. To do so, two libraries are constructed for each angle of attack, namely

and

SINDy models are then fitted for each angle of attack. In all cases, the structure of the identified model is

The coefficients , , and are reported in table 1 for each AoA. Figure 11 provides a comparison of the ground truth evolution (black solid lines) of the oscillations' amplitude and phase with that predicted by the identified models (gray dashed lines). Figure 12 provides a comparison of the growth rate and oscillation frequency predicted by the models against those obtained in the linear transient of the URANS equations. For all angles of attack considered herein, excellent agreement is obtained. This comparison can serve to infer a physical meaning to the identified parameters.

While and correspond to the growth rate and oscillation frequency, the parameter is the cause for the nonlinear saturation and attainment of the limit cycle.

Despite the complex physical mechanisms at play, the identified models thus show that, from a dynamical point of view, the transonic buffet can be reduced to a simple Stuart-Landau oscillator. It should be noted that this finding is consistent with what concluded by Bouhadji and Braza [START_REF] Bouhadji | Organised modes and shock-vortex interaction in unsteady viscous transonic flows around an aerofoil: Part II: Reynolds number effect[END_REF] for the Reynolds number effect on zero-incidence transonic aerofoils. Additionally, we have > 0 and > 0 for all angles of attack, highlighting that the oscillatory dynamics of the transonic buffet arise due to a supercritical Hopf bifurcation of the underlying steady base flow. These observations