C T Nitschke

P Cinnella

D Lucor

J.-C Chassaing
email: jean-camille.chassaing@upmc.fr

Model-form and predictive uncertainty quantification in linear aeroelasticity

Keywords:

Bayesian techniques are employed to quantify model-form uncertainty in aeroelasticity.

• Variabilities in the flutter speed of a typical airfoil section are investigated.

• Inference is performed using a stochastic formulation of the analytical lift function.

• An efficient adjusted model is obtained by considering a possible bias in the random error term.

Introduction

Considerable attention has been paid during the last decade to the development of probabilistic methods for uncertainty quantification in aeroelasticity [START_REF] Beran | Uncertainty quantification in aeroelasticity[END_REF][START_REF] Georgia | Uncertainty quantification of aeroelastic stability[END_REF]. Accounting for inherent variabilities in the aeromechanical parameters is of utmost importance for achieving robust aeroelastic designs in regards with the determination of the probability of observing undesirable fluid-structure phenomena like flutter and limit-cycle-oscillations [START_REF] Lee | Nonlinear aeroelastic analysis of airfoils: bifurcation and chaos[END_REF]. Forward uncertainty quantification of aleatory uncertainty traditionally consists in characterizing the probability density functions (pdf) of input random parameters which are propagated through an aeroelastic model to predict the stochastic response of the quantity of interest. Various sources of uncertainty can lead to important variations in the critical

Nomenclature m

Total moving mass of the system k h

Translational Mean value of the error term flutter speed [START_REF] Pettit | Uncertainty quantification in aeroelasticity: Recent results and research challenges[END_REF][START_REF] Verhoosel | Uncertainty and reliability analysis of fluid-structure stability boundaries[END_REF]. Boosted by the advances in the development of non-intrusive statistical approaches, the scope of forward propagation of aleatory uncertainty in aeroelasticity has been widened substantially [START_REF] Lindsley | Nonlinear plate aeroelastic response with uncertain stiffness and boundary conditions[END_REF][START_REF] Marques | Transonic aeroelastic stability predictions under the influence of structural variability[END_REF]. However, a complete description of the variability in the predictions from an aeroelastic model must account for the quantification of model assumptions and predictive uncertainties [START_REF] Clyde | Model uncertainty[END_REF][START_REF] Park | A Bayesian approach for quantification of model uncertainty[END_REF][START_REF] Rojas | Application of a multimodel approach to account for conceptual model and scenario uncertainties in groundwater modelling[END_REF]. In many aeroelastic problems, such epistemic uncertainties, which results from a lack of knowledge, may induce greater variability in the stochastic response than real physical randomness [START_REF] Thomas | An investigation of the sensitivity of f-16 fighter flutter onset and limit cycle oscillations to uncertainties[END_REF]Zhang et al.). Model-form uncertainty is related to the fact that several plausible models may be employed to evaluate the flutter boundary. These uncertainties can be attributed to the selection of a class of aeroelastic models, like for instance the choice of low-or high-fidelity structural and aerodynamic operators. It is also well known that for a given model class, the selection of most adequate or accurate models for computing the aeroelastic response can be a very delicate point since the sensitivity of these models may strongly vary depending on the physical scenario of interest. Moreover, achieving reliable aeroelastic designs requires the quantification of the predictive uncertainty of a given aeroelastic model. This type of uncertainty may be viewed as the error between the flutter velocity prediction given by an aeroelastic model and the true physical value of the flutter speed. Bayesian model averaging (BMA) [START_REF] Leamer | Specification Searches[END_REF][START_REF] Hoeting | Bayesian model averaging: A tutorial[END_REF] and adjustment factors methods [START_REF] Mosleh | The assessment of probability distributions from expert opinions with an application to seismic fragility curves[END_REF][START_REF] Zio | Two methods for the structured assessment ofmodel uncertainty by experts in performance assessments of radioactive waste repositories[END_REF] are two popular mathematical tools which can be employed to account for model-form and predictive uncertainty [START_REF] Beck | Updating models and their uncertainties. I: Bayesian statistical framework[END_REF][START_REF] Marzouk | A stochastic collocation approach to Bayesian inference in inverse problems[END_REF]Park et al., 2010a;[START_REF] Zhang | A comprehensive bayesian approach for model updating and quantification of modeling errors[END_REF]. Both approaches rely on the framework of Bayesian inference but the latter uses expert opinion when no experimental data are available for the construction of an adjusted model. Basically, the application of the Bayesian inference method for parameter calibration and model updating involves the following steps. First, a set of parametrized stochastic models is built where prior distributions of model parameters are specified according to physical considerations. Second, a statistical model, which relates the model to the observations, is constructed by defining the likelihood function. This represents the probability of observing the data given model coefficients. Next, the joint posterior pdf of model parameters is determined using Bayes' theorem. This calibration step requires the use of Markov Chain Monte Carlo (MCMC) sampling algorithms for high-dimensional parameter vectors. Finally, an adjusted model is constructed by weighting the predictive distribution of each of the previously calibrated individual models by using their corresponding posterior model probability.

Although Bayesian updating is being increasingly used by the Computational Fluid Dynamics (CFD) community to estimate the parameters of turbulence models in both RANS [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF][START_REF] Oliver | Bayesian uncertainty quantification applied to RANS turbulence models[END_REF]Edeling et al., 2014b, a;[START_REF] Guillas | Bayesian calibration of the constants of the kε turbulence model for a CFD model of street canyon flow[END_REF][START_REF] Papadimitriou | Bayesian uncertainty quantification of turbulence models based on high-order adjoint[END_REF][START_REF] Prudhomme | Adaptive surrogate modeling for response surface approximations with application to bayesian inference[END_REF] and LES modeling [START_REF] Tran | A sparse grid method for Bayesian uncertainty quantification with application to large eddy simulation turbulence models[END_REF], few works deal with the study of epistemic uncertainty in aeroelastic flutter boundaries. [START_REF] Dwight | Reducing uncertainty in aeroelastic boundaries using experimental data[END_REF] applied a Bayesian parameter estimation of structural uncertainties to the Goland wing problem by considering a highfidelity CFD/Finite Element aeroelastic solver. This study was performed using a probabilistic collocation approach as surrogate model for the propagation step. Substantial reduction in both the uncertainty in the flutter speed and in the parameters variability were observed using few observation data. In [START_REF] Riley | Quantification of model-form and predictive uncertainty for multi-physics simulation[END_REF] constructed an adjusted statistical model for the prediction of the flutter velocity of the AGARD 445.6 weakened wing by using two adjustment factor approaches and BMA. Model-form uncertainty quantification relied on the use of three commercially available aeroelastic solvers. The highest computed highest posterior model probability corresponded to the model considered as the more realistic model among candidates according to expert judgement. A complete uncertainty quantification process of this probabilistic aeroelastic problem was presented in [START_REF] Riley | Quantification of Model-form, Predictive, and Parametric Uncertainties in Simulation-based Design[END_REF], where parametric uncertainty are propagated by coupling the Fast Fourier Transforms technique with a weighted-Stack Response Surface method. A parameter identification study about the characterization of the free decay response of a two-degrees-of-freedom typical section model, was performed in [START_REF] Khalil | Application of Bayesian inference to the flutter margin method: New developments[END_REF] by applying Bayesian inference to flutter speed estimation. A Bayesian parameter estimation was performed by [START_REF] Khalil | Probabilistic parameter estimation of a fluttering aeroelastic system in the transitional Reynolds number regime[END_REF] and [START_REF] Sandhu | Bayesian model selection for nonlinear aeroelastic systems using wind-tunnel data[END_REF] to investigate self-sustained oscillations in pure pitching motion of a NACA0012 airfoil in the transitional Reynolds number regime by using a simplified analytical aerodynamic model. To this end, a Markov Chain Monte Carlo algorithm was combined with the extended Kalman filter technique to build the joint posterior probability density function in the presence of noisy and dense experimental data [START_REF] Khalil | Probabilistic parameter estimation of a fluttering aeroelastic system in the transitional Reynolds number regime[END_REF]. This methodology was extended to BMA in [START_REF] Sandhu | Bayesian model selection for nonlinear aeroelastic systems using wind-tunnel data[END_REF] by computing an estimate of the evidence using the [START_REF] Chib | Marginal likelihood from the metropolis hastings output[END_REF] method. Moreover, a parallel adaptive MCMC sampling algorithm was proposed as a computationally more efficient alternative to the conventional MCMC method for the evaluation of the joint pdf. The Bayesian Model Selection approach was applied in Sandhu et al. (2016b) to the calibration of fully unsteady nonlinear aerodynamic model using windtunnel test data at various numbers. Finally, an application of the Bayesian analysis to the Zimmerman and Weissenburger flutter margin method [START_REF] Zimmerman | Prediction of flutter onset speed based on flight testing at subcritical speeds[END_REF] was proposed in Sandhu et al. (2016a).

In the present study, a Bayesian statistical framework is employed to quantify and reduce model-form uncertainty in the canonical flutter problem of an incompressible two-degree-of-freedom typical airfoil section. The model-form uncertainty is here a generic designation that encompasses both the model selection choice and the inherent parametric uncertainty intrinsic to the chosen models. In our case, it emerges from the use of several available approximations of the Theodorsen's lift function that is employed to compute the unsteady aerodynamic forces acting on the moving structure. The objective of this work is to construct an adjusted stochastic model for the flutter speed by probabilistically combining the predictive distribution of individual models whose uncertain coefficients are calibrated using experimental values of the flutter speed obtained from different scenarios. To this end, a special attention is paid to the computation of the hyperparameters that describe the predictive error term. Note that the present work, which can be seen as a proof of concept, is restricted to the use of small size experimental data sets and to the calibration of low-order aeroelastic models. Nevertheless the proposed strategy may potentially be extended to more complex configurations in future work.

The outline of this paper is as follows. Section 2 describes the key ingredients involved in the derivation of the studied loworder aeroelastic system, including a brief description of several approximations of the Theodorsen's circulation function and surveys available experimental data. The stochastic parametrization of the aerodynamic model is introduced in Section 3. Section 4 is devoted to the Bayesian calibration of the aerodynamical models, including the construction of an adjusted stochastic model using Bayesian model averaging. Section 5 illustrates the benefit of incorporating systematic errors of the deterministic aeroelastic model in the discrepancy model. Concluding remarks, including a discussion about possible extensions of the present methodology to more complex aeroelastic configurations are presented in Section 6. Finally, the effect of increasing the number of data points in the calibration step is investigated in Appendix using synthetic data.

Deterministic aeroelastic model

We consider the canonical flutter problem of an oscillating typical airfoil section in an incompressible flow. Although the Theodorsen's lift model [START_REF] Theodorsen | General Theory of Aerodynamic Instability and the Mechanism of Flutter[END_REF]Bisplinghoff et al., 1955;[START_REF] Fung | An Introduction to the Theory of Aeroelasticity[END_REF][START_REF] Leishman | Principles of Helicopter Aerodynamics[END_REF] was developed three quarters of a century ago, this frequency-domain based aerodynamical operator is still very popular since it provides a simple but accurate description of the unsteady aerodynamic loads [START_REF] Olsman | Acoustic forcing to simulate the plunging motion of an airfoil[END_REF][START_REF] Amandolese | Low speed flutter and limit cycle oscillations of a two-degree-of-freedom flat plate in a wind tunnel[END_REF][START_REF] Sterenborg | Experimental benchmark of a free plunging wing with imposed flap oscillations[END_REF]. Moreover, state-space representations for this lift model can be derived to develop feedback control laws and optimal controllers by using rational function approximations to the Theodorsen's transfer function [START_REF] Jones | Operational Treatment of the Nonuniform-lift Theory in Airplane Dynamics[END_REF][START_REF] Vepa | Finite state modeling of aeroelastic systems[END_REF][START_REF] Coller | Structural non-linearities and the failure of the classic flutter instability[END_REF][START_REF] Leishman | State-space representation of unsteady airfoil behavior[END_REF][START_REF] Brunton | Empirical state-space representations for Theodorsens lift model[END_REF]. The purpose of the present work originates in the model-form uncertainty associated to the choice of these function approximations and the construction of an adjusted model for robust stochastic predictions of the flutter speed.

Governing equations

We consider a typical airfoil section model subject to a two-degree-of-freedom motion with plunge displacement of the elastic axis h and angular pitch displacement α, again about the elastic axis. The additional mechanical parameters, which are represented in Fig. 1, are the static unbalance x α , and the non-dimensional distance between the mid-chord and the elastic axis a h .

The coupled bending-torsion equations of motion can be written in non-dimensional form as follows [START_REF] Theodorsen | General Theory of Aerodynamic Instability and the Mechanism of Flutter[END_REF] -

ḧ b -x α α -ω 2 h h b + L mb = 0 (1) -r 2 α α -x α ḧ b -r 2 α ω 2 α α + M α mb 2 = 0 (2)
where r α is the radius of gyration about the elastic axis, ω h and ω α are the uncoupled plunge and pitch frequencies of the undamped motion, b represents the half chord, m is the mass, and the upper points denote double derivation with respect to time. The aerodynamic operator is described through the lift L and the pitching moment M α about the elastic axis. Here the unsteady aerodynamic forces are modeled assuming a small amplitude harmonic motion of the 2DOF airfoil subject to an incompressible flow [START_REF] Fung | An Introduction to the Theory of Aeroelasticity[END_REF])

L = -ρb 2 ( V π α + π ḧ -πba h α) -2π ρVbC(k) { V α + ḣ + b ( 1 2 -a h
) α}

(3)

M α = -ρb 2 [ π ( 1 2 -a h ) Vb α + π b 2 ( 1 8 + a 2 h ) α -a h π b ḧ] + 2ρVb 2 π ( 1 2 + a h ) C (k) { V α + ḣ + b ( 1 2 -a h ) α} (4)
where V is the airspeed, k = ωb/V denotes the reduced frequency and ω is the flutter angular frequency. This aerodynamical model includes added-mass and quasi steady effects and the influence of wake vorticity by means of Theodorsen's transfer function C (k) [START_REF] Theodorsen | General Theory of Aerodynamic Instability and the Mechanism of Flutter[END_REF])

C (k) = J 1 (k) -ıY 1 (k) (J 1 (k) + Y 0 (k)) -ı(Y 1 (k) -J 0 (k)) (5)
where J n and Y n are the Bessel functions of first and second kind, respectively. Note that this expression of C (k) is valid for harmonic motion only.

Flutter solution method

Assuming time-harmonic oscillations of the structure, the governing equations can be reformulated in the frequency domain and written in matrix form as where ω denotes the frequency of the motion, [M] is the mass matrix, [K] represents the stiffness matrix, [A e ] is the aerodynamic matrix, which is a function of the reduced frequency k and {r} is the vector of generalized coordinates. The iterative frequency matching V -g method (Bisplinghoff et al., 1955) is used in this work to evaluate the critical flutter velocity V f . This velocity-damping approach consists in incorporating an unknown damping coefficient g into Eq. ( 6) leading to the following eigenvalue problem

( [K] -ω 2 [M] ) {r} = ω 2 [A e ] {r} (6)
([M] + [A e ] -λ [K]) {r} = {0} (7)
where λ = (1 + ıg)/ω 2 . The frequency ω, the artificial damping g and the nondimensional velocity V * can be computed from each complex eigenvalue as

ω = 1 √ Re(λ j ) , g = Im(λ j ) Re(λ j ) , V * = V bω α . (8)
In practice, the determination of the nondimensional linear critical flutter velocity V * f involves the following steps: (i) start the iterative procedure by selecting a value for the reduced frequency k, (ii) solve the eigenvalue problem (Eq. ( 7)), (iii) for each eigenvalue, compute ω, g, V * and update the V -g and Vω plots, (iv) decrease the reduced frequency k and repeat steps (ii)-(iv) until a single mode changes from negative to positive modal damping. The critical flutter speed V * f is found to be the lowest value of V * for which g > 0, as illustrated in Fig. 2.

A bisection method is used in the velocity-damping solver to further refine the root of the artificial damping. This requires to specify a tolerance criterion on the value of the critical flutter velocity. A typical value of 10 -5 is used to minimize spurious numerical noise.

Approximations to Theodorsen's lift function

Alongside the analytical lift function by [START_REF] Theodorsen | General Theory of Aerodynamic Instability and the Mechanism of Flutter[END_REF], multiple approximated formulae for C (k) have been proposed in the literature. Earlier approximations proposed by [START_REF] Jones | Operational Treatment of the Nonuniform-lift Theory in Airplane Dynamics[END_REF] and [START_REF] Jones | Aerodynamic Forces on Wings in Nonuniform Motion[END_REF] aimed at quantifying the aerodynamic behavior of wings due to instationary loads generated as a response to wind gusts. Later approximations, mostly used in control theory, were directly derived from Theodorsen's function either by Padé approximation [START_REF] Vepa | Finite state modeling of aeroelastic systems[END_REF] or by a balanced truncation development of state-space realizations of the system [START_REF] Brunton | Empirical state-space representations for Theodorsens lift model[END_REF]. [START_REF] Riley | Quantification of Model-form, Predictive, and Parametric Uncertainties in Simulation-based Design[END_REF] provides some examples of common approaches in use. The approximations considered in this paper are summarized in Table 1. These functions can be divided in two families according to the number of states they are derived from and, as a consequence, to the number of terms they involve. Models of a given family share the same mathematical form, but involve different coefficients. In this work we consider either two-state or four-state models.

The real and imaginary parts of the approximations of C (k) presented in Table 1 are plotted in Fig. 3 as a function of the reduced frequency, for each one of the families considered in this study.

Although the different models are analytic approximations of the same function, they exhibit differences in phase shift driven by the imaginary part of the function. For instance, the phase shifts of models [START_REF] Jones | Operational Treatment of the Nonuniform-lift Theory in Airplane Dynamics[END_REF][START_REF] Jones | Aerodynamic Forces on Wings in Nonuniform Motion[END_REF] differ significantly from that of Theodorsen's original function. In contrast, the return to the stationary system response is relatively similar.

Note that the rounded version of R. T. Jones' approximation given by [START_REF] Riley | Quantification of Model-form, Predictive, and Parametric Uncertainties in Simulation-based Design[END_REF], deviates significantly from the original function and the other approximations, both in phase shift as well as decline to stationary response, despite its coefficients exhibiting only slight differences with respect to the original model.

Obviously, it is expected that the different forms of the function affect the V -g analysis. In particular, differences in the phase shift between the models lead to different estimates of the artificial damping and, consequently, to different critical flutter velocities.

C (k) ≈ 1.0 -0.165k k-0.0455ı -0.335k k-0.3ı
Jones ( 1945)

C (k) ≈ 1.0 -0.165k k-0.041ı -0.335k k-0.32ı
Riley ( 2011)

C (k) ≈ (1.0+10.61ık)(1.0+1.774ık) (1.0+13.51ık)(1.0+2.745ık)
Jones rounded [START_REF] Riley | Quantification of Model-form, Predictive, and Parametric Uncertainties in Simulation-based Design[END_REF])

C (k) ≈ 0.015+0.3ık-0.5k 2 0.015+0.35ık-k * * 2
Four states Brunton and Rowley ( 2013) 32510 -4 k 4 -1.158ık 3 -0.3052k 2 +0.02028ık+2.32510 -4 Vepa ( 1977) 551ı10 -3 k+9.55710 -6 2k 4 -1.064ık 3 -0.1134k 2 +2.617ı10 -3 k+9.55710 - 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.1 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.1 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.1 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Fig. 3. Plots of typical approximations of Theodorsen's function taken from [START_REF] Riley | Quantification of Model-form, Predictive, and Parametric Uncertainties in Simulation-based Design[END_REF] and [START_REF] Brunton | Empirical state-space representations for Theodorsens lift model[END_REF] as described in Table 1.

C (k) ≈ 0.5k 4 -0.703ık 3 -0.2393k 2 +0.01894ık+2.
C (k) ≈ k 4 -0.761ık 3 -0.1021k 2 +2.

Table 2

Experimental data for V * f , corresponding to four values of ω h /ω α . The other aeroelastic parameters are taken fixed, namely r α = 0.5, x α = 0.2, a h = -0.4, µ = 400 [START_REF] Theodorsen | General Theory of Aerodynamic Instability and the Mechanism of Flutter[END_REF].

Scenario A B C D ω h /ω α 0.33 0.5 0.83 1. V * f 10.67
9.19 6.41 7.30

Experimental data

In this section we describe the experimental data obtained by [START_REF] Theodorsen | General Theory of Aerodynamic Instability and the Mechanism of Flutter[END_REF] for four aeromechanical configurations, also called scenarios in the following, characterized by different values of the plunge and pitch frequency ratio. The scenarios and their corresponding measured critical flutter speed indices are illustrated in Table 2. [START_REF] Theodorsen | General Theory of Aerodynamic Instability and the Mechanism of Flutter[END_REF] for the four scenarios of Table 2. The other symbols refer to the model results. The blue dashed curve is a least squares polynomial fit of order seven to the experimental data.

Deterministic analysis

The aeroelastic responses of the deterministic models corresponding to different approximations of C (k) (see Table 1) are computed for each scenario in Table 2. The results obtained in terms of critical flutter velocity index V * f , are presented in Fig. 4.

The different deterministic models provide relatively similar results. They all tend to underestimate the flutter speed, except for scenario C , where the model results are closely distributed around the experimental value. The highest spread is observed for scenario B, for which the 1945 model from [START_REF] Jones | Aerodynamic Forces on Wings in Nonuniform Motion[END_REF] behaves as an outlier and fortuitously provides a solution significantly closer to the experimental values than the other models, as shown in Fig. 4.

In general, four-state models, which approximate Theodorsen's function more accurately, provide results for V * f that are also closer to those of the exact model, as shown in Table 3.

In summary, changing the lift function model introduces a dispersion of the prediction results, both because of the mathematical form of the function in use (e.g. two-state or four-state) but also because of the coefficients associated to a given class of functions. As the number of examples shows, different choices are possible for the form of the function and for the values of their parameters. For this reason, in the following we introduce a statistical methodology for quantifying and reducing model uncertainties in aeroelastic problems.

Stochastic framework

In this Section, the six deterministic lift function models under investigation are casted into two families of stochastic models. A sensitivity analysis is performed on their uncertain parameters and the impact of both parametric and model-form uncertainties on predictions of the critical flutter velocity is quantified.

Stochastic models

In order to derive the stochastic models, the two families of function approximations in Table 1 are brought to the following form where N is the number of states of the models. N = 2 and N = 4 will be considered in the present work. The resulting stochastic models will be referred to as C 2 (k) and C 4 (k) in the following. Corresponding values for the coefficients α j and β j for the deterministic models are given in Table 4.

C (k) ≈ 1 - N ∑ j=1 α j k k -β j ı (9)
To highlight the physical role of the parameters, the inverse Laplace transform of the two-states approximation is considered, which relates function C(k) and its approximations to the time domain.

L -1 {C} = 1 - N ∑ j=1 α j exp(-β j τ ). (10)
In the time space, we can identify coefficients α j as influence factors, whereas the β j represent time constants.

Given the variety of values that the α j and β j coefficients can take according to the specific model in use, in the following we model them as random variables, described by a probability density function (pdf ) and introducing a parametric uncertainty on the solution outcome, which also becomes random.

To construct the a priori probability density functions, the following requirements are made to ensure the physical admissibility of the stochastic models:

(1) When the wing is not in oscillatory motion, the response is required to fall back to the stationary solution. This means that for a reduced frequency approaching zero, the real part of the function value has to tend to one and the imaginary

part to zero lim k→0 ℜe(C(k)) → 1, lim k→0 ℑm(C(k)) → 0 (2) It is also required that lim k→∞ ℜe(C(k)) → 1 2 .
(3) The imaginary part of the function must always be smaller than or equal to zero ℑm(C(k)) < 0, meaning that there is a positive phase shift between the flow movement and the reaction of the wing.

(4) The real part of the function is required to decrease monotonically between the previously defined limits ℜe(C(k)) ∈ [ 1 2 , 1].

To meet requirements 1, 3 and 4, it is sufficient to restrict the parameters to positive values. To respond to requirement 2, the influence factors have to sum up to 0.5. For this reason, we replace the α j by the normalized coefficients αj such that αj =

α j ∑ N j=1 α j . (11)
In the following, we treat the whole set of coefficients αj and β j as parametric uncertainties of the function C(k) and we calibrate them using the available experimental values for the critical velocity by means of a Bayesian inference approach.

As a first step, we assign to coefficients αj and β j an a priori probability density function (pdf), characterizing their random behavior. These pdf are then updated in the calibration step.

The normalized influence factors αj are assumed to be independent and uniformly distributed in the interval [0, 1]. The time constants β j are also modeled as independent and uniformly distributed random variables, with a lower bound equal to 0 and an upper bound equal to 0.9, which is about 10% higher than the highest time constant found in the present models [START_REF] Brunton | Empirical state-space representations for Theodorsens lift model[END_REF]. In summary, we consider the following distributions for αj and β j αj ∼ U [0, 1] , β j ∼ U [0, 0.9] , j = 1, . . . , N.

(12)

Let θ = [ α; β] ∈ R 2×N be the vector of all uncertain coefficients associated to a given Nstate stochastic model. With the aforementioned assumption on coefficients αj and β j , plus the additional hypothesis of statistical independence of the parameters, the joint pdf for θ i is just the product of the preceding uniform pdfs. Note that the choice of considering independent model coefficients for the definition of the prior pdf in Eq. ( 12) does not exclude the possibility of the presence of correlation in the posterior pdf, since the calibration process inherently accounts for dependencies between the parameters (see Section 4.1). We propagate in Fig. 5 the described uncertainty through the aeroelastic model by means of Monte Carlo sampling with 10 7 samples. The mean and the 50% maximum credibility interval obtained for the two stochastic models of the approximation of C (k). For comparison purpose, the exact deterministic value of C (k) is also represented. We observe that the 2-and 4-states stochastic model are driven by extreme outliers. As a consequence, their mean values lay outside the 50% confidence intervals for large parts of the considered frequency range. However, the maximum credibility intervals and the plots of the realizations show that the functions are complying to the physical requirements. The spread is relatively wide. Results suggest that calibration might be needed, because already Theodorsen's original function is outside the 50% credibility interval I 0.5 .

The stochastic representations of C (k) are used to estimate the output pdf of the predicted critical flutter velocity index. The results of direct propagation of the prior parameter distribution are shown in Fig. 6. The models perform well for scenarios B, C and D, the value σ being at about 0.5. For scenarios A, C and D, the experimental data can be found within an interval of one standard deviation σ . However, the spread of V * f values for scenario A is very high, suggesting extremely poor accuracy of the models for these conditions, even if the result found by experiment is within 1 • σ of the mean response.

Sensitivity analysis

A variance-based sensitivity analysis of the critical flutter velocity is carried out by means of the analysis of total Sobol coefficients [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF].

We consider the two stochastic models of C (k) introduced in Section 3.1. Hereafter, we call M 1 the two-state model and M 2 the four-state model. We note y the model output, obtained by using the aeroelastic solver supplemented by a lift function model M i :

y = y(x, θ i , M i )
where x is a vector of explanatory variables, e.g. variables that are not considered for calibration, and θ i is the vector of uncertain parameters associated to model M i :

θ i = {( αj ) 1,...,N , (β j ) 1,...,N }.
In the present work, the model output y for V * f is directly obtained by the V -g solver. The explanatory variables are represented by the aeromechanical parameters of the aeroelastic configuration under consideration, i.e.

x = {x α , r α , ω h , ω α , a h , µ}.

(13) Several studies have shown the importance of taking into account uncertainties associated to the definition of such parameters in aeroelastic systems. However, since in this work we wish to focus on uncertainties associated to different formulations of the lift function, the explanatory variables x are considered as deterministic in the following.

The Sobol coefficients can be seen as partial variances of V f associated with the uncertain model coefficients θ i . For instance, the total Sobol index of the kth component of θ i reads

S θ i k = Var(y) -Var (E(y|θ i (m))) Var(y) (14)
where m are all subsets of indices that do not contain k. The total Sobol indices are approximated by means of Monte Carlo sampling of the prior pdf, with 10 7 samples. It was found that model parameters of the same nature (influence factors or time factors) have the same behavior in terms of Sobol indices. For this reason, in the following we only typical results for one of the α parameters and one of the β parameters for each class of stochastic models, as well as one typical second-order index, generically denoted as S αβ , which accounts for coupled influences of coefficients of both kinds.

By inspecting Fig. 7, it can be noted that the time factors are the more influential parameters. Thus, it is expected that this family of parameters can be more easily updated by the subsequent calibration step. We also note that second order influences are significant. The sensitivity to a combination of time factors is lower than that to an influence/time factor pair, indicating that it is the coupling of the two that changes the flutter result. However, since the ratios of the influence factors α decide on the contribution of the different time scales, the total sensitivity to combinations of the influence factors is particularly high.

Bayesian approach

To reduce the uncertainty associated with the choice of model parameters, experimental data are assimilated to the stochastic lift function models introduced in the previous Section by means of Bayesian inference. To this end, the concept of likelihood is introduced to model the dispersion of model predictions around the observed data. Subsequently, the marginal likelihood is used to evaluate and compare different model classes.

Calibration method

In this Section, we introduce a statistical methodology for calibrating model coefficients of a given model M i using the available experimental observations. As models M i , we consider again the two stochastic models of C(k) introduced in Section 3.1, where M 1 denotes the two-state model C 2 (k) and M 2 the four-state model C 4 (k). Using these models, we aim to make predictions of V * f -the quantity of interest -which will be denoted q in the following. Contrary to deterministic approaches, it is assumed that the model output which will be noted y in general does not match the quantity of interest exactly. This discrepancy can be attributed to either model inadequacy or measurement error. The quantity of interest is thus modeled as the output of the deterministic model plus a random error term ε

q = y(x, θ i , M i ) + ε(µ i , σ i , M i ) (15)
where x are the so-called explicative parameters which define the problem, and θ i are the random model parameters subject to epistemic uncertainties as developed in Section 3.1. In many cases, the error term ε is chosen to be Gaussian [START_REF] Carlin | Bayesian Methods for Data Analysis[END_REF] because the final result is assumed to have a fixed mean offset and a known variance, in which case a Gaussian is the distribution with the highest information entropy. Moreover, in many cases, the dominant error in the discrepancy term is the measurement error which is usually modeled using a normal law. The mean µ i and the standard deviation σ i are the hyperparameters that describe the error term. They can be different from model to model.

Experimental data points appear in this framework as realizations of the stochastic quantity of interest. Under the hypothesis that the mean of the error term is zero, and assuming furthermore that the realizations of the error term for different observed experimental data points d j are independent, we can write the likelihood that the model outputs correspond to the observed data points as Beck and Katafygiotis (1998)

f N (D|θ i , σ i , M i ) = Π n d j=1 1 √ 2π σ 2 i exp ( - (d j -y(x j , θ i , M i )) 2 2σ 2 i ) (16)
with D the set of individual data points d j , n d the number of elements of D and index j is used to indicate the individual scenarios associated to the data. Note that there is a specific vector of explanatory parameters x j associated, also being indexed to match the scenario. In contrast, x without index j is used to define a new scenario for the prediction step. The individual models are then updated using the likelihood function in Bayes' rule (see e.g. Beck and Katafygiotis, 1998)

p(θ i |D, M i ) ∝ f N (D|θ i , M i )p(θ i |M i ) (17)
where p(θ i |M i ) represents the joint prior probability density of the uncertain parameters associated to model M i , f N is the likelihood function, i.e. the probability of observing the data D given a model M i and a set of parameters θ i , and p(θ i |D, M i ) is the joint posterior distribution of the model parameters, i.e. the output of the calibration procedure. The data D corresponds to observations of the quantity of interest q, i.e. to a set of experimentally observed values for V * f . To this end, we will take the data from the scenarios presented in Section 2.4. The next step in the application of the Bayesian framework consists in describing the statistics in the error term ε in Eq. ( 15). In this work, two statistical modeling of the error terms are analyzed in Sections 4.2 and 4.3.

Calibration using the evidence method

Up to now, the descriptive parameters of the error term, which we call hyperparameters in the following because they are not part of the physical problem, are not given yet. A common approach is to set the error mean to zero and use the standard deviation of the experimental noise as the error standard deviation [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF]. However, such information is not available in the present case. Hereafter, the calibration step of the Bayesian framework is performed using the evidence method. This approach is depicted in Section 4.2.1 and corresponding calibration results are presented in 4.2.2.

Choice of the hyperparameter

A possibility to determine hyperparameter σ is to optimize the marginal evidence, also called marginal or integrated likelihood, with respect to the parameter that is looked for. The marginal likelihood of a model is given by Cheung et al. ( 2011)

P(D|M i ) = ∫ f N (D|θ i , M i )p(θ i |M i )dθ i . (18)
Optimizing this term with respect to the hyperparameters is equivalent to picking the value which gives the highest posterior model probability.

For better legibility of the results, we use the Bayesian Information Criterion (BIC), which is often used in the context of Bayesian Model Selection and is defined as Lebarbier and Mary-Huard ( 2004) Usually, a Laplace approximation of the marginal likelihood is used to compute the BIC. We avoid the Laplace approximation because of the small number of data points in D. Instead, P(D|M i ) is approximated by Monte Carlo integration, propagating the prior model parameter distribution through the likelihood function as

BIC (M i ) = -2 ln(P(D|M i )). (19)
P(D|M i ) ≈ 1 n s ns ∑ ℓ=1 f N (D|θ i (ℓ), M i ) (20)
where θ i (ℓ)are the ℓth samples of p(θ i |M i ). The direct Monte Carlo approach is affordable here because the aeroelastic solver is mostly analytic, thus calls to it are computationally cheap. Throughout this paper, 10 6 samples are used for the evaluation of the integrated likelihoods. Note that the same values for the integrated likelihood will also be used for the application of the model averaging step in Eq. ( 26).

In the following, all four data points from Table 2 are employed to compute the BIC. Results on the Bayesian information criteria are summarized in Table 5 and plotted in Fig. 8.

In the present case, the values of σ i minimizing this coefficient are at σ 2 ≈ 0.6 for M 2 and at σ 1 ≈ 0.75 for M 1 . For too high values of σ i , the BIC coefficients increase. This is due to the fact that the allowed discrepancy margin is too wide to optimally exploit the information provided by the calibration data. Note that the BIC coefficient also rises quickly for too low values of σ i , indicating that the model is not able to capture the data within the given variation and does not benefit from the data in an optimal way.

Calibration results

In this Section, we apply our statistical framework to an aeroelastic 2DOF pitch and plunge system defined by r α = 0.5, x α = 0.2, a h = -0.4 and µ = 400. The model-form uncertainties are given by the different number of states in M 1 and M 2 . The parametric uncertainty is described by the random coefficients θ = [ α; β] introduced in Section 3.1. The experimental data set

D = {d A , d B , d C , d D } is given in Table 2.
Posterior distributions of the model parameters are computed using the calibration procedure described in Section 4.1, using all of the scenarios D = {d A , d B , d C , d D }. Initially, the hyperparameter σ i in the likelihood model is taken fixed and the calibrations are repeated for different choices of this hyperparameter. Samples from the joint posterior model parameter distributions are drawn from Eq. ( 17) using the Markov Chain Monte Carlo algorithm with Metropolis-Hasting sampling. A total of 250 000 samples, of which 50 000 are used for burn-in and are thus ignored. Fig. 9 shows the posterior of model coefficients α1 and β 1 where the uniform prior distributions of model parameters are also represented for comparison. Note that the other αj and β j parameters exhibit similar posteriors as α1 and β 1 respectively.

In general, the parameters are weakly informed by the calibration data. We observe that the probability density of the influence factor α1 for model M 1 decreases below values of about 0.25. This can be interpreted in a way that the calibrated models attribute weight to all time scales available. This remark holds for model M 2 (Fig. 9(c)).

At the same time, the time factors β 1 show a sharp peak below values of about 0.1, (Fig. 9(b),(d)), indicating that small frequency scales are important for the calibration. This confirms a tendency visible in the deterministic models, which all have at least one very low β value. We also see that the time factor probability density function there is a drop right of the first sharp peak for σ = 0.4, which is as pronounced for the higher σ values. This indicates that there might be separate time scale regimes, however, the evidence is not strong enough to confirm that.

It can also be noticed that in general, parameters are better informed when using lower values of σ i . The calibrated model is forced to remain closer to the data, which is only possible for some narrow ranges of the parameters. Moreover, the fact that the distributions stretch to the limits of the prior uniform distribution means that the final distributions of the parameters are not independent of the chosen prior. This is most probably due to the small number of data used for the calibration. Fig. 10 shows the joint probability density functions for influence/time factor pairs ( α1 , β 1 ). The two-state model M 1 exhibits a region of high probability density for small time factors of under 0.1 at low values of the influence factor α1 . A second region of high probability density is observed for values of β 1 between 0.3 and 0.5 and α1 > 0.6.

For the four-state model M 2 , however, the whole parameter space is more evenly used. The higher probability density field for low time factor values stretches now to influence factor values down to 0.1. The observed pairings coincide with the values used in the deterministic model. It is also noticed that having more time/influence factor groups leads to a higher exploitation of the parameter space.

The results of forward propagation of the posterior samples are shown in Figs. 11 and 12 along with the forward propagation results of the prior parameter distributions for comparison. For both models, the standard deviation is set to σ = 0.6, which represents the optimum in terms of BIC for the four-state model M 2 . Strictly speaking, this value is too low for the two-state model M 1 , but we choose the same value for both models to comparison easier. For both models, the output in critical velocity gets higher after calibration for all the data points, which means that parts of their distributions shift to the right. However, in all the cases, a peak observed in the a priori propagation left of the observed values persist, indicating that the models are not informed ideally. Note that for point d C , there is a probability density mass above the observed value for V * f , although this point was already calculated with relatively good accuracy before. The correction that the model has to compensate an underestimation of V * f might be reflected in that result. Table 6 summarizes the first two statistical moments of the prior and posterior sample distributions of V * f . The mean value of the critical flutter speed becomes closer to the reference value in almost all cases. Moreover, the calibration process succeeds in reducing the value of the standard deviation σ V * f for both models.

Calibration using hyperparameter inference

Up to now, we used the BIC to decide on a deterministic value for σ i . However, this can be difficult because its computation is either expensive or obtained using approximations with relatively strong assumptions [START_REF] Lebarbier | Le critère BIC : fondements théoriques et interprétation[END_REF]).

An alternative consists in considering σ i as an uncertain hyperparameter that is calibrated together with the uncertain model parameters during the inference process (Eq. ( 17)). In the following, we perform the calibration with σ i taken as random parameter. The posterior distribution Eq. ( 17) is thus changed to

p(θ i , σ i |D, M i ) ∝ f N (D|θ i , σ i , M i )p(θ i , σ i |M i ). (21)
To include the least information possible, we take a uniform distribution for σ . We still assume that the models are capable of capturing the point. We thus choose the lower limit of the distribution just above zero. We then assume that the models deliver at least predictions with the most likely offset, which is the σ found using the BIC, or maximally just above. We thus set the upper limit of the distribution to the σ found by the BIC plus a good 10%, which gives the distribution σ ∼ U [0.01, 0.7].

The corresponding posterior distributions of σ i are shown in Fig. 13. It is observed that the probability density function of σ i drops for σ i → 0, indicating that the models are not capable of directly reproducing the results without the addition of a discrepancy term. The higher probability density for σ 2 below of σ = 0.3 means that it is possible to obtain more accurate results with the four-state model M 2 .

We also note that the pdf s of σ i do not differ much from a uniform distribution, meaning that the hyperparameters are not informed well enough to establish an upper error bound within the prior. However, considering σ i as a hyperparameter might be a good compromise in the absence of known values. We clearly observe that adding σ i in the inference process attributes higher probability density to the searched value. The behavior is similar for both models.

Table 7 shows that the hyperparameter inference over σ i gives more reliable values of µ V * f compared to the use of σ i based on the analysis of the BIC coefficients (Table 6).

Finally, we examine the influence of the number of data points used in the calibration step. All computations are performed using hyperparameter inference. The posterior parameter distributions obtained using different numbers of data points in the calibration are presented in Fig. 15. Adding data leads to sharper posteriors, especially for parameters β i . So better informed posteriors can be obtained as soon as additional data become available.

Table 7

Posterior sample distribution moments using hyperparameter inference.

M 1 M 2 A B C D A B C

Robust predictions

The posterior distribution of the random parameters derived in Section 4.1 can be used to predict an updated estimate of the quantity of interest. According to [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF], the marginal posterior predictive distribution for the quantity of interest reads

p(q|x, D, M i ) = ∫ p(q|x, θ i , σ i , M i )p(θ i , σ i |D, M i )dθ i dσ i (22)
where p(q|x, θ i , σ i , M i ) is the predictive distribution for a given set of parameters and hyperparameters. It is constructed using the definition of the quantity of interest including the error term (15), so the error is added to the model output. The predictive distribution for one point in the parameter space reads consequently Practically, the approximation of the marginalization in Eq. ( 22) is then performed using Monte Carlo integration as follows

p(q|x, θ i , σ i , M i ) = 1 √ 2π σ 2 i exp ( - (q -y(x, θ i , M i )) 2 2σ 2 i ) . (23)
p(q|D, x, M i ) ≈ 1 n s ns ∑ ℓ=1 p(q|x, θ i (ℓ), σ i (ℓ), M i ) (24)
where θ i (ℓ) and σ i (ℓ) are the ℓth sample of p(θ i , σ i |D, M i ). Note that hyperparameter inference is used during the calibration process and consequently the distribution of σ is also marginalized over.

In the following, we will take interest in the prediction of a new point for which no data is available. For comparison purposes, only three of the points of data set D = {d A , d C , d D } will be used for calibration, while point d B , defined by ω h ωα = 0.5 and V * f = 9.19 is used for the validation of the prediction step.

Fig. 16 shows the predictive distributions for the two-and four state models using different numbers of data points for calibration.

The reference is indicated by the black vertical bar. Despite the fact that differences between the two models are small it can be seen that adding data points improves the accuracy of the predictions in the sense that the probability to recover the experimental value of V * f at point d B increases. The local hump observed for the pdf of V * f for M 1 when D = {d A , d c } can be interpreted as the tentative of the calibration step to account for the extreme V * f values in the experimental set. This hump thus disappears when adding case d d in the experimental set, resulting in a smoother predictive distribution of V * f . However, once the first two data points have been added, the resulting information gain is small in the present case. Adding more data points, and specifically close to the prediction point could help to further improve the results. The possible benefits of considering larger data sets are demonstrated in Appendix, where synthetic data points are added to the calibration process.

Bayesian model averaging

In order to account for the uncertainty in the selection of the model itself, Bayesian Model Averaging (BMA) is used to construct a stochastic model from a set of aerodynamic models involved in the computation of the critical flutter speed. Here, the discrepancy between the model output and the realizations of the critical flutter velocity V * f is modeled as a Gaussian distribution as in Eq. ( 15). The output distribution of the BMA is an average of the posterior predictive distributions of the m models weighted by the posterior model probability of each model [START_REF] Hoeting | Bayesian model averaging: A tutorial[END_REF])

p(q|x, D) = m ∑ i=1 p(q|x, D, M i )P(M i |D) (25)
where p(q|x, D) is the total predictive distribution for the variable of interest q averaged over the different models given an experimental data set D, p(q|D, x, M i ) corresponds to the robust or posterior predictive distribution of model M i and P(M i |D) represents the posterior model probability. Here, the vector of explanatory deterministic parameters x is explicitly mentioned for the sake of clarity, but we drop it in the following to simplify the notation.

Again, Bayesian inference is employed to evaluate the model probability P(M i |D) used in the BMA (Eq. ( 25)) [START_REF] Hoeting | Bayesian model averaging: A tutorial[END_REF])

P(M i |D) = P(D|M i )P(M i ) ∑ m j=1 P(D|M j )P(M j ) (26)
where P(M i ) represents the prior model probability assumed to follow a discrete uniform distribution in the present work. P(D|M i ) is again the marginal likelihood as given in Eq. ( 18). The equation is adapted to hyperparameter inference in the following. For a given model M i , the corresponding marginal likelihood P(D|M i ) is computed as

P(D|M i ) = ∫ f N (D|θ i , σ i , M i )p(θ i |M i )p(σ i |M i )dθ i dσ i ( 27
)
where p(θ i |M i ) represents the prior density which is defined based on expert opinion as described in Section 3.1 and p(σ i |M i ) is the prior probability density distribution of the hyperparameter σ i .

In the following, the averaging step of Eq. ( 25) will be applied. To achieve this, the posterior model probabilities are computed, which are listed in Table 8. We notice that the more complex model M 2 is assigned a higher model probability than M 1 . Similar trends are obtained for both a fixed σ i and the hyperparameter model. Now we compare the total predictive distribution p(q|x, D) for point d B obtained using BMA performed on the stochastic models M 1 and M 2 with the total predictive distribution p(q|x, D) determ. based on BMA of the deterministic models listed in Table 1. To this end, the model coefficients are set to their nominal values as given in Table 4. Individual calibrations are applied by hyperparameter inference on σ i over data set D = {d A , d C , d D }. Fig. 17 clearly illustrates the benefit of calibrating model coefficients in addition to the model averaging process. Indeed, considering deterministic aerodynamic models M i gives a lower probability density value for the searched value of V * f than those obtained using BMA based on stochastic models.

Moreover, differences among the individual predictions of stochastic models M 1 and M 2 and their BMA counterpart remain small since results given by the individuals model after calibration are quite similar.

Attempts are conducted in the next Section to reduce confidence intervals of the prediction.

Considering a bias

As already observed in Fig. 4, Theodorsen's original model systematically underestimates the flutter speed index. This might introduce a systematical underestimation for the other deterministic models as well, as at least some of them were created to mimic the original function's behavior. In this context, supposing a zero-mean error for the stochastic model might not be an optimal choice. In the following, an error representation with non-zero mean is adopted to take into account a possible systematical error.

Consequently, the error term used in the Bayesian inference is defined as ϵ(µ i , σ i , M i ) where µ i denotes the mean of the random error term for model M i .

The likelihood function is then

f N (D|θ i , σ i , M i ) = Π n d j=1 1 √ 2π σ 2 i exp ( - (d j -y(x j , θ i , M i ) -µ i ) 2 2σ 2 i ) (28)
and the predictive distribution reads

p(q|x, θ i , σ i , M i ) = 1 √ 2π σ 2 i exp ( - (q -y(x, θ i , M i ) -µ i ) 2 2σ 2 i ) . (29)
Both σ i and µ i are treated as stochastic hyperparameters during the inference process. In order to select reasonable ranges for the prior distributions of µ i and σ i , we first carry out a preliminary optimization of these coefficients by performing an analysis of the BIC. As shown in Fig. 18, minimal BIC values are obtained for values of about µ i = 0.7 and an error standard deviation of about 0.35. According to these results, we adopt hereafter the following distributions for the hyperparameters: σ i ∼ U [0.01, 0.5] and µ i ∼ U [0.55, 0.85]. The joint marginal posterior distributions of (µ i , σ i ) are shown in Fig. 19.

The probability density maximum for µ i is situated around µ 2 = 0.7 for M 2 , whereas M 1 tends to prefer higher values of µ 1 .

The maximum of the probability density of the standard deviation of the error is obtained for values of σ i ranging from 0.3 to 0.4. This represents a significant reduction in error standard deviation in comparison to the calibrations without bias as given in Fig. 13. The difference with respect to the BIC analysis (Fig. 18) can be explained by the fact that the posterior distributions obtained for µ i and σ i become correlated due to the inference, whereas the BIC is a global criterion calculated using uncorrelated prior distributions.

As shown in Figs. 20 and 21, the prediction based on calibration using non-zero error succeeds much better in capturing data points d A , d B and d D than its unbiased counterpart. The uncertainty is also reduced significantly. This can be explained by the fact that, since the mean of the posterior predictive distribution cannot shift for µ = 0 during the calibration step, the adjusted model captures the data by increasing the variance in p(y). On the other hand, considering a bias in the definition of the error terms introduces an additional stochastic degree of freedom, leading to a reduction in the variance in p(y). In addition, predictions for the d B data point are improved for both model because the probability to observe the reference value is higher.

The plot of the total posterior predictive distribution p(q|D) in Fig. 22 as obtained by BMA shows that the support of possible realizations of V * f is noticeably reduced for calibrations including the bias effects compared to BMA based on zero mean error p(q|D) µ=0 . In addition, the probability density for the experimental value predicted is much higher. Moreover, the maximum of the probability density function is much closer to the value as well. We can conclude, for the present case, that the calibrated model with ϵ(µ i , σ i , M i ) shows better robust predictions than BMA with a zero mean error.

Concluding remarks

In this work, Bayesian calibration is performed on a set of low-order aeroelastic models in order achieve a better level of confidence in prediction of the critical flutter boundary of a pitching and plunging airfoil due to model choice uncertainties related to the approximation of the Theodorsen's lift function. Instead of considering these model-form uncertainties from a deterministic point of view, stochastic models are constructed from a set of mathematical formulations using physically admissible uncertain coefficients. Afterwards, Bayesian model averaging is performed on these stochastic models and the standard deviation characterizing the random error term is considered either as a deterministic or as stochastic hyperparameter.

Results show that the predictive distribution of a new flutter speed scenario experimentally assessed (for which validation experimental data are available), based on the calibration of both stochastic models, is more robust than the one obtained from the application of conventional BMA with deterministic models. It was also observed that performing the calibration over multiple scenarios avoids overfitting problems that would be obtained by the calibration of individual points. In addition, the present methodology provides each individual model with a probability associated to its contribution on the global model-form uncertainty. Such information is of crucial importance since it can be directly used to discard or to confirm the need to consider a particular model for constructing a combined stochastic model for reliable prediction of the flutter speed of the aeroelastic system of interest. In the present case, the four-state model has the higher model probability compared to its two-states counterpart. This seems to confirm that this more sophisticated model with higher number of parameters, captures more physics to properly adapt to the experimental results.

Although we acknowledge that, both the use of a limited set of data per experimental scenario and the study of loworder aeroelastic models, limits the generality of the present results, the proposed Bayesian calibration framework could be employed to predict reliable flutter margins from subcritical flight test data. To end, future work should focus on the uncertainty quantification and Bayesian inference of the aeromechanical damping of high-fidelity aeroelastic models in the presence of combined parametric and epistemic uncertainties. In such context, more efficient strategies than conventional MCMC algorithm should be employed to sample the posterior pdf due to the expensive computational models. Potential candidates are, among others, the Hamiltonian Monte Carlo [START_REF] Cheung | New Bayesian updating methodology for model validation and robust predictions based on data from hierarchical subsystem tests[END_REF] method, the Importance sampling algorithm [START_REF] Li | An adaptive importance sampling algorithm for Bayesian inversion with multimodal distributions[END_REF] or parallel adaptive MCMC (Sandhu et al., 2016a,b). Additionally, computational costs can be greatly alleviated by using statistical response surface models to explore the state space of the model parameters instead of the costly computer code. Surrogate models, like for instance Kriging based response surface methods, the stochastic spectral projection or the stochastic collocation approach, which are employed for uncertainty quantification in Computational Fluid Dynamics [START_REF] Jeroen Witteveen | Subcell resolution in simplex stochastic collocation for spatial discontinuities[END_REF][START_REF] De Baar | Efficient uncertainty quantification with gradient-enhanced kriging: Applications in FSI[END_REF][START_REF] Chassaing | Stochastic investigation of flows about airfoils at transonic speeds[END_REF], and could be coupled with Bayesian calibration of high-fidelity aeroelastic systems [START_REF] Dwight | Reducing uncertainty in aeroelastic boundaries using experimental data[END_REF]. Eventually, the reduction of the random space dimension can be carried out by performing a variance-based sensitivity analysis of the uncertain model coefficient before performing the calibration step. Recently, such approaches were successfully used in Bayesian applications [START_REF] Sraj | Coordinate transformation and Polynomial Chaos for the Bayesian inference of a Gaussian process with parametrized prior covariance function[END_REF].

Appendix. Performing a BMA study with additional interpolated data

The aim of this Appendix is to investigate the effect of increasing the number of data points in the calibration step of the present BMA framework. In the absence of available additional experimental data, we consider synthetic data which are obtained from the interpolation of the experimental data set of V * f for different values of ω h /ω α , as depicted in Fig. 4. From a practical point of view, n d = 20 equally-spaced observation data are extracted from the range ω h /ω α = [1/3, 1].

A.1. Calibration

The calibration step is based on the use of the evidence theory as described in Section 4.2 and an analysis of the Bayesian Information Criterion (BIC) is performed similar to those presented in Table 5 and in Fig. 8. The optimal value of the hyperparameter σ which minimizes the BIC is σ = 0.25. Then the posterior distributions of the model coefficients are computed using Importance Sampling with n IS = 10 6 samples. factor p( α1 |D, M 1 ) is observed for n d = 20, meaning that this model has a tendency to choose its weight parameters equal. On the contrary, this distribution remains mainly flat for n d = 4, meaning that no global information on the parameters can be observed. Moreover, the posterior distribution of the time factor p( β1 |D, M 1 ) obtained using a calibration over 20 observation data show a bi-modal behavior which was not observed with a small data set.

As far as calibration of the coefficients of the four-state model M 2 is concerned, the bi-modal shape of p( β1 |D, M 2 ) obtained using n d = 20 is also visible although less pronounced than for p( β1 |D, M 1 ). Finally, we remark that the posterior distribution p( α1 |D, M 2 ) is nearly insensitive to the use of a larger data set for the calibration.

A.1.1. Robust predictions

Next, BMA is applied to compute the total predictive distribution p(q|D) for point d B as done in Section 4.5 for D = {d A , d C , d D }. The four-state model has the higher model probability with about 70%, which means that the model benefits from the higher number of parameters, being more capable to properly adapt to the experimental results. It is interesting to remark that performing calibration with few data is sufficient to give consistent values of P(q|x, D, M i ).

Fig. 1 .

 1 Fig. 1. (a) Two-degree-of-freedom pitch-and-plunge airfoil aeroelastic model.

Fig. 2 .

 2 Fig.2. Plots of modal damping g versus non-dimensional velocity V * and reduced frequency k using the aerodynamics of[START_REF] Theodorsen | General Theory of Aerodynamic Instability and the Mechanism of Flutter[END_REF] and r α = 0.5, x α = 0.2, a h = -0.4, µ = 400, ω h /ω α = 0.83. The critical flutter conditions are given by k f = 0.13 and V * f = 6.4 as indicated by the black vertical bar.

Fig. 4 .

 4 Fig. 4. Linear flutter boundary as a function of the ratio of uncoupled natural frequencies in pitch and plunge. Crosses show the experimental values[START_REF] Theodorsen | General Theory of Aerodynamic Instability and the Mechanism of Flutter[END_REF] for the four scenarios of Table2. The other symbols refer to the model results. The blue dashed curve is a least squares polynomial fit of order seven to the experimental data.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Mean and 0.5 confidence intervals of approximation families of Theodorsen's function, with 10 realizations each, compared to the analytical function.

Fig. 7 .

 7 Fig. 7. Total Sobol indices for scenario B (ω h /ω α = 0.5).

Fig. 9 .

 9 Fig. 9. Posterior of model coefficients for different values of the hyperparameter σ i of model M i .

Fig. 10 .

 10 Fig. 10. Heatmap of a influence/time factor pair for M 1 and M 2 .

Fig. 11 .

 11 Fig. 11. Prior and posterior sample distributions for V * f due to parametric uncertainties in two-state model M 1 given D = {d A , d B , d C , d D } and computed for all points in D using σ 1 = 0.6.

Fig. 12 .

 12 Fig. 12. Prior and posterior sample distributions for V * f due to parametric uncertainties in four-state model M 2 given D = {d A , d B , d C , d D } and computed for all points in D using σ 2 = 0.6.

  Fig. 14 presents a comparison of the posterior distributions of V * f obtained using σ i = 0.6 with those resulting from the hyperparameter inference. The calibration is performed using D = {d A , d B , d C , d D } and the pdfs are plotted for data point d B .

Fig. 13 .

 13 Fig. 13. Kernel density estimations of the posterior distributions of σ i for models M 1 and M 2 . The prior distribution is σ i = U [0.01, 0.7].

Fig. 14 .

 14 Fig. 14. Posterior sample distribution of V * f at d B given D = {d A , d B , d C , d D } for σ i = 0.6 and σ i ∼ U [0.01, 0.7].

Fig. 15 .

 15 Fig. 15. Influence of size of the calibration data set D on the posterior of model coefficients computed using hyperparameter inference with σ i ∼ U [0.01, 0.7].

Fig. 16 .

 16 Fig. 16. Influence of the number of experimental data on the predictive distributions of V * f at point d B .

Fig. 17 .Fig. 18 .

 1718 Fig. 17. Combined model from BMA, predicting d B using data points D = {d A , d C , d D } and σ i ∼ U [0.01, 0.7].

Fig. 19 .Fig. 20 .

 1920 Fig. 19. Heatmap of error mean and error standard deviation joint marginal probability density distribution for M 1 and M 2 , with data D = {d A , d B , d C , d D } for calibration.

Fig. 21 .

 21 Fig. 21. Effect of considering a bias in the random error term on the posterior predictive distributions of V * f due to model-form uncertainties in M 2 given D = {d A , d B , d C , d D }. Hyperparameters are defined by σ 2 ∼ U [0.01, 0.5] and µ 2 ∼ U [0.55, 0.85].

Fig. 22 .

 22 Fig. 22. Combined statistical model for robust predictions the of the distribution of the critical flutter speed computed using BMA with σ i ∼ U [0.01, 0.5] and µ i ∼ U [0.55, 0.85].

  Fig. A.1 compares the posterior distribution of model coefficients α1 and β 1 obtained using n d = 20 data points with those obtained using a calibration over D = {d A , d B , d C , d D } and σ i = 0.6 (Fig. 9). The examination of the results of the two-state model M 1 shows that the parameters of this model better inform in the presence of additional data, since the range of the support of most probable values is more narrow compared to the calibration over D = {d A , d B , d C , d D } only. A preference for mid-interval values in the posterior distribution of the influence

Fig. A. 1 .

 1 Fig. A.1. Calibration step: Posterior of model coefficients using n d 4 and n d = 20 where n d is the number of scenarios which are considered for the calibration.

Fig. A. 2 .

 2 Fig. A.2. Prediction step: results of combined statistical models for robust prediction at scenario B constructed using calibration using 3 and 20 data points.

  Fig. A.2 shows that accounting for a larger data set significantly increases the probability of obtaining realization of the flutter speed index V * f close to the reference value of the reference data point d b (V * f = 9.19). At the same time, the support of possible realizations of V * f is noticeably reduced compared to the BMA performed with D = {d A , d C , d D }.

Table 1

 1 Approximations of Theodorsen's function used to construct the stochastic lift functions.

	Nb. of states	Reference	Function definition
	Two states	Jones (1938)

Table 3

 3 Results for critical flutter velocity V *

	ω h ωα	Exp.	Theodorsen	Jones (1938)	Jones (1945)	Riley	Jones rounded	Vepa	Brunton
	0.33	10.67	9.967	9.952	10.018	9.805	9.943	9.937	9.965
	0.5	9.19	8.029	8.107	8.797	7.843	8.171	7.973	8.024
	0.83	6.41	6.312	6.402	6.499	6.317	6.261	6.293	6.321
	1.0	7.30	6.960	7.255	7.044	6.973	6.805	6.949	6.966

f with the different lift function models presented in Table

1

.

Table 4

 4 Coefficients of the approximations of C (k) according to Eq. (9).

		α 1	α 2	α 3	α 4	β 1	β 2	β 3	β 4
	Jones (1938)	0.165	0.335			0.0455	0.3
	Jones (1945)	0.165	0.335			0.041	0.32
	Riley and Grandhi (2011)	0.2346	0.2664			0.074	0.3643
	Jones rounded	0.1	0.4			0.05	0.3
	Brunton and Rowley (2013)	0.0396	0.1555	0.2438	0.0612	0.0144	0.0786	0.2522	0.8128
	Vepa (1977)	0.0128	0.0333	0.2279	0.2259	0.0045	0.0257	0.1042	0.3976

Table 5

 5 Computation of the Bayesian information criteria (BIC) for values of the standard deviation σ of the error term ranging from 0.3 to 1.0. Influence of the standard deviation in ϵ(σ i , M i ) on the calibration process for M 1 and M 2 , obtained using {d A , d B , d C , d D }.

	σ	0.3	0.4	0.5		0.6		0.7		0.8		0.85	0.9	0.95	1.0
	BIC(M 1 )	12.15	11.57	11.19		10.78		10.55		10.51		10.54	10.60	10.69	10.80
	BIC(M 2 )	9.60	9.14	9.07		9.06		9.15		9.35		9.49	9.63	9.79	9.97
			12
			11.5
			11
			10.5
			10
			9.5
			9 0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
		Fig. 8.

Table 6

 6 Prior and posterior sample distribution moments.

		Scenario	M 1				M 2
			A	B	C	D	A	B	C	D
	Prior	µ	10.071 2.806	8.067 0.569	6.231 0.455	6.524 0.605	9.755 1.238	8.218 0.450	6.358 0.404	6.643 0.522
	Posterior	µ σ	10.051 0.469	8.535 0.495	6.507 0.352	6.807 0.345	10.113 0.370	8.561 0.400	6.591 0.287	6.890 0.312

Table 8

 8 Posterior model probabilities P(q|x, D, M i ) for calibration over , d C , d D } obtained from BMA with different modeling of the random error term ϵ(σ i , M i ).

	Model error	P(M 1 |D)	P(M 2 |D)
	σ i = 0.6	0.3526	0.6474
	σ i ∼ U [0.01, 0.7]	0.311	0.689

D

= {d A

Table A .1

 A Individual posterior model probabilities obtained for calibrations with different size of data set. Table A.1 summarizes the posterior model probabilities P(M|D) for calibration over D = {d A , d C , d D } (Table 8) and the extended data set with n d = 20.

	Observation data	Model error	P(M 1 |D)	P(M 2 |D)
	n d = 20	σ i = 0.25	0.3077	0.6923
	n d = 3	σ i = 0.6	0.3526	0.6474
	n d = 3	σ i ∼ U [0.01, 0.7]	0.311	0.689
	Finally,