# Efficient parametric study of a stochastic airfoil system based on hybrid surrogate modeling with advanced automatic kriging construction

Enora Denimal, Jean-Jacques Sinou

► To cite this version:

Enora Denimal, Jean-Jacques Sinou. Efficient parametric study of a stochastic airfoil system based on hybrid surrogate modeling with advanced automatic kriging construction. European Journal of Mechanics - A/Solids, 2023, 99, pp.104926. 10.1016/j.euromechsol.2023.104926. hal-04093089

HAL Id: hal-04093089

https://inria.hal.science/hal-04093089

Submitted on 9 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Denimal, E., and J-J. Sinou. "Efficient parametric study of a stochastic airfoil system based on hybrid surrogate modelling with advanced automatic kriging construction." *European Journal of Mechanics-A/Solids* (2023): 104926. https://doi.org/10.1016/j.euromechsol.2023.104926

# Efficient parametric study of a stochastic airfoil system based on hybrid surrogate modelling with advanced automatic kriging construction

Enora Denimal<fnref n="1" />,\* and Jean-Jacques Sinou<fnref n="2" />,<fnref n="3" />

## Abstract

Flutter is one of the most important aeroelastic instability phenomena that arises from the interaction between the structural dynamics of the mechanical airfoil system and the surrounding airflow. This instability phenomenon can lead not only to a reduction in aircraft performance but also to catastrophic structural failure.

Therefore, one of the major challenges is to perform parametric and sensitivity studies on the stability behaviour of a wing system subject to many random uncertainties in order to achieve a thorough understanding and reliable estimation of the role played by each parameter in the flutter phenomenon. To carry out such a study, an advanced surrogate modelling technique based on kriging and polynomial chaos expansion (PCE) is proposed for the prediction of flutter instability. In addition, a methodology based on hybrid surrogate modelling with advanced automatic kriging construction is discussed to promote an efficient parametric study of the airfoil system with uncertainties subjected to flutter. The Sobol indices highlight that the role played by each random parameter depends strongly on the flow speed and airfoil geometry with complex behaviours, giving valuable insights into the physics and the complexity of flutter.

**Keywords** stability analysis, airfoil flutter, hybrid uncertainties, kriging, polynomial chaos expansion, sensitivity analysis

## 1 Introduction

Flutter is a well-known instability problem in aero-elasticity [1, 2]. It is a dangerous phenomenon which can lead to lower performances or even to dramatic structural failures. This phenomenon comes from an interaction between the structure and the surrounding air flow that affects the effective structural mass, damping and stiffness matrices which become asymmetric. For a critical speed of the airflow, i.e. the flutter speed, the airfoil starts to absorb energy from the surrounding airflow leading to a dynamic instability called flutter. To suppress vibration and instability, active and passive solutions have been extensively investigated [3–7]. However, designing an effective structure robust to environmental and manufacturing tolerances that limits or avoids the appearance of flutter in a first stage is of major importance. Yet, the complexity of the flutter phenomenon and the high number of parameters involved in the modelling make this task complex. Indeed, when designing such structure, one needs to identify the best set of some design parameters to ensure the stability of the system under the high number of model uncertainties present, that might come from environmental considerations, material variability or manufacturing tolerances for example. The latter are often considered as *aleatory* parameters and are modelled with random variables, whereas the design parameters can be seen as *epistemic* uncertainties as increasing the number of samplings reduces the uncertainty [8]. Due to the presence of aleatory uncertainties, the response of the airfoil is stochastic, but this stochastic response varies when the design parameters vary as well. In this context, it is necessary to deal with different natures of uncertainties, and so advanced and complex numerical tools are required. Moreover, the presence of a large number of uncertain parameters make

<footnote n="1">¹Univ. Gustave Eiffel, Inria, COSYS/SII, I4S, Campus de Beaulieu, 35042 Rennes, France</footnote>

<footnote n="2">²Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, École Centrale de Lyon, 36 avenue Guy de Collongue 69134 Écully Cedex, France</footnote>

<footnote n="3">³Institut Universitaire de France, 75005 Paris, France</footnote>

<footnote>\* Corresponding authors: enora.denimal@inria.fr</footnote>
any parametric study or sensitivity study challenging due to the high numerical time involved, even for simple
academic models. But getting deep insights in the role and influence of each parameter is essential to gain a
better understanding of the phenomenon and the role played by each parameter for different flow speeds. Thus,
the main objective of the paper is to propose a deep analysis of the stability of an airfoil model by taking into
consideration both types of uncertainties. As it may represent a high numerical cost, advanced numerical tools
based on hybrid surrogate modelling are developed to ensure an affordable computational cost and extensive
parametric study of the stochastic system and parametric sensitivity analysis are performed.

Most studies are focused only on the propagation of one type of uncertainties. Different methods are available
to propagate them. If the uncertain parameters can be modelled as random variables, the most classical and
robust method remains the Monte Carlo Sampling (MCS) strategy which approximates the stochastic response
by realising a large number of runs of the model. Even though it is robust, the convergence rate is low and the
computational cost is too expensive to consider real applications. Advanced sampling strategies also exist which
reduces the number of samples, as the Latin Hypercube Sampling (LHS) [9], however this approach remains
expensive. MCS has been used for example in [10] to solve the random flutter equations of an airfoil model under
turbulent flow conditions. Or in [11] to study the influence of the boundary conditions on a two-dimensional
panel flutter system in the absence or presence of random pressure. More advanced and efficient techniques
as Polynomial Chaos Expansion (PCE) [12] for example, are possible. PCE has been extensively considered,
with the generalized PCE (gPCE) [13] or arbitrary PCE (aPCE) [14] to deal with different probability density
functions (PDF). In [15], the PCE approach is used to represent the periodic response of an airfoil when the
structural stiffness are random. In [14], the aPCE is employed to study the stability of an airfoil model when
the natural frequency or the structural nonlinearity are uncertain. In [16], the bifurcating behaviour of an
airfoil model is predicted with PCE when three uncertain parameters are considered. In [17], an intrusive
PCE approach is developed to predict the flutter instability of an airfoil system when the structural stiffness
is random. Other approaches are based on the use of Response Surface Methods (RSM). The latter is cheap
to evaluate and can be exploited to perform MCS studies for uncertainty quantification (UQ) for example.
Methods as polynomial regression, kriging [18], support vector regression [19], etc. are in this category. These
methods are also often used for parametric studies and design space exploration, and are particularly relevant for
epistemic uncertainty. It consists in the approximation of the model with an equivalent analytical mathematical
expression. A comparison of PCE, Singular Value Decomposition (SVD), Fuzzy and MCS for flutter prediction
is proposed in [20]. Advanced multi-fidelity models to couple low-fidelity and high-fidelity models to predict
 flutter have also been studied [21, 22]. In cases where both nature of uncertainties are present, i.e. aleatory and
epistemic, adopting a unique modelling strategy might not be the best numerical strategy. Recent works [23]
have proposed hybrid surrogate modelling to deal with both epistemic and aleatory uncertainties based on the
association of PCE and kriging to model both uncertainties. By considering that the PCE coefficients depend
on the epistemic uncertainties, they can then be approximated with a kriging surrogate model. This hybrid
surrogate model is particularly efficient when a parametric study of a stochastic model is required, but also
for a parametric sensitivity analysis. Previous works have demonstrated the robustness of the approach for
systems with large number of uncertain parameters and/or large and complex mechanical systems and can also
be improved by integrating the physical properties in the surrogate modelling formulation [24–26]. However, it
is known that the choice of the kriging properties (regression and correlation functions) have a large influence on
its performances [27]. These properties are often chosen by hand based on previous knowledge, or by comparing
kriging predictions to reference points. In the hybrid surrogate model, as each PCE coefficient is modelled with
a kriging model, such strategies are not possible due to the large number of surrogate models. It justifies the
necessity to find and identify efficient criterion to construct automatically these kriging surrogate models with
the best properties, which is one objective of this study. More precisely, different criterion are considered as it
is possible to consider criterion based on the full hybrid surrogate model or construction criterion based on the
error of each kriging surrogate model (i.e. on each PCE coefficients).

Finally, the objective of the paper is to perform a parametric study and a sensitivity study of the stability
behaviour of an airfoil system subjected to numerous random uncertainties, i.e. when two different types of
uncertainties are present. By doing this design space exploration and this parametric sensitivity analysis, the
goal is to gain deep insights in the role played by each parameter in the flutter phenomenon. An advanced
hybrid surrogate model that couples kriging and PCE is employed to deal with both natures of uncertainties.
A secondary objective of the study, is to propose criterion that can be used for an automatic construction of
the different kriging surrogate models to improve the performances of the hybrid surrogate model.

The paper is structured as follows. First, the model of the airfoil and the stability analysis for flutter detection
is presented. Second, the uncertain parameters in the model are presented. Third, the hybrid surrogate model
and the different automatic construction criterion considered are presented. Finally, the results are presented.
In a first time, the different criterion are compared in terms of convergence rate and efficiency. In a second
time, once the final hybrid surrogate model is validated, it is exploited to analyse deeply the stability of the airfoil in an uncertain context and a sensitivity analysis is perform. Complex behaviours and dependencies are shown, justifying the use of such advanced tools.

## 2 Airfoil model and basic stability analysis

The proposed airfoil model corresponds to a two-degrees-of-freedom system which are the plunging deflection $h$ (positive in the downward direction) and the pitching angle $\alpha$ about the elastic axis (positive nose up), as illustrated in Figure 1. This system is subjected to the quasi-steady aerodynamic lift $L_{qs}$ and moment $T_{qs}$ that can be expressed for a two-dimensional and incompressible flow such as [5, 7]:

$$L_{qs} = \pi \rho_a b s_p \left( b \left( \ddot{h} + V \dot{\alpha} - b \bar{a} \dot{\alpha} \right) + 2V \left( V \alpha + \dot{h} + b \left( \frac{1}{2} - \bar{a} \right) \dot{\alpha} \right) \right) \quad (1)$$

$$T_{qs} = \pi \rho_a b^2 s_p \left( b \bar{a} \ddot{h} - V b \left( \frac{1}{2} - \bar{a} \right) \dot{\alpha} - b^2 \left( \frac{1}{8} + \bar{a}^2 \right) \ddot{\alpha} + 2V \left( \frac{1}{2} + \bar{a} \right) \left( V \alpha + \dot{h} + b \left( \frac{1}{2} - \bar{a} \right) \dot{\alpha} \right) \right) \quad (2)$$

The definitions of all the parameters are given in Table 1. Then the equation of the airfoil model can be defined by:

$$\mathbf{M}_s \ddot{\mathbf{x}} + \mathbf{C}_s \dot{\mathbf{x}} + \mathbf{K}_s \mathbf{x} = \mathbf{F} \quad (3)$$

where $\mathbf{x} = [h \ \alpha]^T$ and the vector $\mathbf{F}$ includes the contribution of the quasi-steady aerodynamic lift and moment such as $\mathbf{F} = [-L_{qs} \ T_{qs}]^T$. The mass, damping and stiffness matrices of the structural mechanical system ($\mathbf{M}_s, \mathbf{C}_s, \mathbf{K}_s$) are defined by

$$\mathbf{M}_s = \begin{bmatrix} m & mx_\alpha \\ mx_\alpha & I_\alpha \end{bmatrix} \quad ; \quad \mathbf{C}_s = \begin{bmatrix} c_h & 0 \\ 0 & c_\alpha \end{bmatrix} \quad ; \quad \mathbf{K}_s = \begin{bmatrix} k_h & 0 \\ 0 & k_\alpha \end{bmatrix} \quad (4)$$

Finally, the system can be rewritten in the following form

$$(\mathbf{M}_s + \mathbf{M}_{aero}) \ddot{\mathbf{x}} + (\mathbf{C}_s + \mathbf{C}_{aero}) \dot{\mathbf{x}} + (\mathbf{K}_s + \mathbf{K}_{aero}) \mathbf{x} = \mathbf{0} \quad (5)$$

where the matrices $\mathbf{M}_{aero}$, $\mathbf{C}_{aero}$ and $\mathbf{K}_{aero}$ correspond to the additional contributions due to quasi-steady aerodynamic flow conditions previously given in the vector $\mathbf{F} = [-L_{qs} \ T_{qs}]^T$. $\mathbf{M}_{aero}$, $\mathbf{C}_{aero}$ and $\mathbf{K}_{aero}$ are given by:

$$\mathbf{M}_{aero} = \begin{bmatrix} \pi \rho_a b^2 s_p & -\pi \rho_a b^3 s_p \bar{a} \\ -\pi \rho_a b^3 s_p \bar{a} & \pi \rho_a b^4 s_p \left(\frac{1}{8} + \bar{a}^2\right) \end{bmatrix} \quad (6)$$

$$\mathbf{C}_{aero} = \begin{bmatrix} 2\pi\rho_a b s_p V & \pi\rho_a b^2 s_p V \left(1 + 2\left(\frac{1}{2} - \bar{a}\right)\right) \\ -2\pi\rho_a b^2 s_p V \left(\frac{1}{2} + \bar{a}\right) & -2\pi\rho_a b^3 s_p V \bar{a} \left(\frac{1}{2} - \bar{a}\right) \end{bmatrix} \quad (7)$$

$$\mathbf{K}_{aero} = \begin{bmatrix} 0 & 2\pi\rho_a b s_p V^2 \\ 0 & -2\pi\rho_a b^2 s_p V^2 \left(\frac{1}{2} + \bar{a}\right) \end{bmatrix} \quad (8)$$

The usual approach to determine the system stability is the Complex Eigenmodes Analysis (CEA). The stability is determined by considering the complex eigenvalues $\lambda_k$ and the associated eigenvectors $\Phi_k$ defined as the solution of:

$$(\lambda_k^2 (\mathbf{M}_s + \mathbf{M}_{aero}) + \lambda_k (\mathbf{C}_s + \mathbf{C}_{aero}) + (\mathbf{K}_s + \mathbf{K}_{aero})) \Phi_k = \mathbf{0} \quad (9)$$

If all eigenvalues have negative real parts, the system is stable. If at least one eigenvalue has a positive real part, the system is unstable. In this second case, if the imaginary part of the associated positive eigenvalue is not equal to zero, the system encounters flutter and the imaginary part corresponds to the angular frequency of the associated unstable mode. To be noted that in the case where the associated imaginary part is equal to zero, the system encounters divergence.

In order to illustrate the problem that will be addressed in the rest of the paper, a first brief study is presented on the notion of stability for the airfoil system. Figure 2(a) displays the evolutions of the real parts and angular frequencies as a function of the flow speed $V$. Figure 2(b) illustrates the appearance of instability in the complex plane. Initially, at $V = 0$ the system is stable with the two modes at angular frequencies equal to $\omega_{I,0} = 10.94$ rad/s $\omega_{II,0} = 37.70$ rad/s, respectively. The be noted that the associated real parts are equal to $\Re(\lambda_{I,0})=-0.318$ and $\Re(\lambda_{II,0})=-7.066$, respectively. Increasing the flow speed $V$ induces an increase of the real part $\Re(\lambda_{I,V})$ and a decrease of the real part $\Re(\lambda_{II,V})$. When the flow speed $V$ reaches 23.46 m/s the system's stability switches. As illustrated in Figure 2(b), one complex eigenvalue crosses the complex plane imaginary axis (i.e. $\Re(\lambda) = 0$ and $\Im(\lambda) \neq 0$). The flutter angular frequency $\omega_f$ at the Hopf bifurcation point (marked by
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic mechanical system of a two-dimensional airfoil [7]</figcaption></figure>

<table><caption>Table 1: Geometrical and physical parameters of the airfoil system</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Value</th></tr></thead><tbody><tr><td>Span of the airfoil</td><td>s<sub>p</sub></td><td>1 m</td></tr><tr><td>Airfoil mass</td><td>m</td><td>2.049 kg</td></tr><tr><td>Dimensionless distance of the elastic axis from midchord</td><td>a&#773;</td><td>-0.6847</td></tr><tr><td>Semichord of the airfoil section</td><td>b</td><td>0.135 m</td></tr><tr><td>Torsional stiffness</td><td>k<sub>&alpha;</sub></td><td>6.833 N m/rad</td></tr><tr><td>Plunging stiffness</td><td>k<sub>h</sub></td><td>2844.4 N/m</td></tr><tr><td>Damping coefficient in pitching</td><td>c<sub>&alpha;</sub></td><td>0.036 N s/rad</td></tr><tr><td>Damping coefficient in plunging</td><td>c<sub>h</sub></td><td>27.43 N s/m</td></tr><tr><td>Air density</td><td>&rho;<sub>a</sub></td><td>1.225 kg/m<sup>3</sup></td></tr><tr><td>Mass moment of inertia of the airfoil about the elastic axis</td><td>I<sub>&alpha;</sub></td><td>mx<sub>&alpha;</sub><sup>2</sup> + 0.0517 kg m<sup>2</sup></td></tr><tr><td>Distance of the elastic axis from center of mass</td><td>x<sub>&alpha;</sub></td><td>(0.0873 - (1 + a&#773;) b) m</td></tr></tbody></table>

a red cross in Figure 2(b)) is equal to 24.32 rad/s. After the Hopf bifurcation, the airfoil system is unstable and the real parts of the two modes are opposite: one mode is stable (i.e. the associated real part is negative and decreases versus the increase of the flow speed) whereas the other one is unstable (i.e. the associated real part is positive and increases versus the increase of the flow speed).

## 3 Description of the parametric and uncertain parameters

In the airfoil model, different parameters are assumed to be uncertain and could vary due to environmental considerations, manufacturing tolerances, wear etc. In this work, it is considered that five parameters, namely $k_{\alpha}$, $c_{\alpha}$, $k_{h}$, $c_{h}$, $s_{p}$, are affected by such variations. A classical modelling choice is to use random variables to model them. They are characterized by a probability law, and the airfoil response becomes stochastic. It is assumed in this work that the random parameters follow a uniform distribution and can vary up to ±5%. It is summarized in Table 2.

However, it is still of prior importance to study the stability behaviour of such stochastic structure when considering a parametric variation of different parameters that drive the system. Indeed, it could be interesting in a design process to select a robust design, or for an extensive parametric study to understand deeply the system behaviour. In this work, it is chosen to study the stability of the stochastic airfoil when the distance of the elastic axis to mid-chord $a$, the flow speed $V$ and the semi-chord length $b$ vary. So for each value of the vector $\mathbf{p} = [a, V, b]$, one is looking for the stochastic response of the airfoil system. The variation of these three parameters are given in Table 3.

The variation of these parameters has an impact on the other parameters of the airfoil, such as the mass or the inertia. The airfoil is assumed to have the same thickness and camber ratio, denoted $\tau$ and $\eta$ respectively (thickness camber and over chord length), during the design and operations. It is known that such system have an area $A$ and a bending inertia $I$ equal to:

$$A = K_A b^2 \tau \qquad (10)$$

$$I = K_I (\tau^2 + \eta^2) b^4 \qquad (11)$$
<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of the stability analysis of the airfoil system – (a) real parts and angular frequencies versus the flow speed V – (b) eigenvalues evolution in the complex plan</figcaption></figure>

<table><caption>Table 2: Random parameters description</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Average value</th><th>Law</th><th>% variation</th></tr></thead><tbody><tr><td>Torsional stiffness</td><td>k<sub>α</sub></td><td>6.833 N m/rad</td><td>Uniform</td><td>± 5 %</td></tr><tr><td>Damping coefficient in pitching</td><td>c<sub>α</sub></td><td>0.036 N s/rad</td><td>Uniform</td><td>± 5 %</td></tr><tr><td>Plunging stiffness</td><td>k<sub>h</sub></td><td>2844.4 N/m</td><td>Uniform</td><td>± 5 %</td></tr><tr><td>Damping coefficient in plunging</td><td>c<sub>h</sub></td><td>27.43 N s/m</td><td>Uniform</td><td>± 5 %</td></tr><tr><td>Span of the airfoil</td><td>s<sub>p</sub></td><td>1 m</td><td>Uniform</td><td>± 5 %</td></tr></tbody></table>

with $K_A$ and $K_I$ constants. The mass of the airfoil is $m = As_p\rho = K_A b^2 \tau s_p \rho = \beta b^2 s_p$ with $\beta = K_a \tau \rho$ a constant, and the bending moment of inertia is $I = \gamma b^4$ with $\gamma = K_I(\tau^2 + \eta^2)$ a constant. From these, the distance of the elastic axis from center of mass is $x_\alpha = 0.0873 - (1 + \bar{a})b$ and the mass moment of inertia of the airfoil about the elastic axis is $I_\alpha = mx_\alpha^2 + I$.

## 4 Hybrid surrogate model and automatic construction criterion

The parameters to predict are the real parts and the imaginary parts of the airfoil eigenvalues, in this context $\lambda_1 = \kappa_1 + i\omega_1$ and $\lambda_2 = \kappa_2 + i\omega_2$, with $\kappa_\bullet$ the real part and $\omega_\bullet$ the imaginary part. For the sake of clarity and concision, the parameter to predict will be denoted $\lambda$ in the following. It can represent $\kappa_1, \kappa_2, \omega_1, \omega_2$. As each variable depends on the vector of design parameter $\mathbf{p} = [\bar{a}, V, B]$ and the vector of random parameters $\xi = [k_{\alpha\rho}, c_\alpha, k_h, c_h, s_p]$, it writes $\lambda(\mathbf{p}, \xi)$.

<table><caption>Table 3: Interval variation of the parametric variables</caption><thead><tr><th>Parameter</th><th>Notation</th><th>[min,max]</th></tr></thead><tbody><tr><td>Dimensionless distance of the elastic axis from midchord</td><td>a&#773;</td><td>[-0.5,-0.8]</td></tr><tr><td>Semichord of the airfoil section</td><td>b</td><td>[0.13, 0.16] m</td></tr><tr><td>Flow speed</td><td>V</td><td>[0, 40] m/s</td></tr></tbody></table>
## 4.1 Polynomial Chaos Expansion

The uncertainty related to the random part is modelled and propagated with PCE. According to the PCE theory, any random parameter Y can be approximated by a convergent, in the $L^2$ sense, expansion [12, 13, 28, 29]:

$$Y(\xi) = \sum_{k=0}^{P-1} \mu_k \Psi_{\pi_k}(\xi) \quad (12)$$

with $(\mu_k)$ the weighting coefficients and $(\Psi_{\pi_k})$ the multivariate polynomials. The expansion is truncated and $P$ terms are retained. Only polynomials of order inferior to $o_{pce}$ are retained when no additional truncation scheme is retained. In this case, the number of term is $P = (o_{pce}^{r+1})$, with $r$ the number of random parameters.

The multivariate polynomial basis $(\Psi_{\pi_k})$ is constructed by tensorisation of monovariate polynomial basis $(\Phi_{\pi}^{(j)})$, with $\pi$ the polynomial order and $j \in \{k_{alpha}, c_{\alpha}, k_h, c_h, s_p\}$. Thus, each random parameter is characterised by a distribution and by a family of polynomials orthogonal w.r.t. this distribution. The correspondence between the distribution and the polynomial family is given by the Askey scheme [13, 28]. As uniform laws are considered here, Legendre polynomials are used. If $\pi_k$ denotes the multivariate index $\pi_k = [\pi_{k_{alpha}}, \pi_{c_{\alpha}}, \pi_{k_h}, \pi_{c_h}, \pi_{s_p}]$, then [12, 13, 28, 29]

$$\Psi_{\pi} = \Phi_{\pi_{k_{\alpha}}}^{(k_{\alpha})} \times \dots \times \Phi_{\pi_{s_p}}^{(s_p)} \quad (13)$$

A hyperbolic truncation scheme of parameter $q \in (0, 1]$ is adopted to reduce the size of the PCE and to select a subset of polynomials with low order interactions [30]. Only multi-index $\pi$ that verify

$$||\pi|| \left( \sum_{j \in [k_{alpha}, c_{\alpha}, k_h, c_h, s_p]} \pi_j^q \right)^{1/q} \le o_{pce} \quad (14)$$

are kept. The lower is $q$, the smaller is the PCE basis. It translates the idea that the main effects are often driven by low order interactions between the different input parameters.

The coefficients $(\mu_k)$ can be computed by intrusive or non-intrusive methods [29]. A non-intrusive method based on regression is employed here as it is the best compromise between numerical cost, accuracy of results and implementation difficulty. The coefficients $(\mu_k)$ are the solution of a least-square minimisation problem between the random function $Y(\xi)$ and its PCE evaluation at $N$ points. The input points of experimental design are generated based on a Latin Hypercube Sampling (LHS) [9].

## 4.2 Kriging

The uncertainty related to the parametric parameters $\mathbf{p} = [\bar{a}, V, b]$ will be modelled and propagated with the kriging method. It approximates a function $\mu(\mathbf{p})$ by a response surface based on a given set of evaluations of the method [18, 31–33]. The kriging approach approximates the function $\mu$ as:

$$\mu(\mathbf{p}) = g(\mathbf{p})^T \beta + Z(\mathbf{p}) \quad (15)$$

where $g$ is a set of $n_r$ regressive functions, often taken as polynomial of low order (up to order 2), the vector $\beta$ is the weighting coefficients and are a solution of a least square problem. $Z$ is a zero mean Gaussian process of variance $\sigma^2$, whose covariance is

$$\mathbb{E}[Z(\mathbf{p}), Z(\mathbf{p}')] = \sigma^2 \mathcal{R}(\theta, \mathbf{p}, \mathbf{p}'). \quad (16)$$

$\mathcal{R}$ is the spatial correlation function of scaling parameter $\theta$, $\mathbf{p}$ and $\mathbf{p}'$ are two points of the input space. Anisotropic kriging is considered here, and so $\theta$ has the same size as $\mathbf{p}$, here 3. The correlation function $\mathcal{R}$ is built based on a product of mono-variate correlation families [18, 31–33], in dimension 3 it writes:

$$\mathcal{R}(\theta, \mathbf{p}, \mathbf{p}') = \prod_{j=1}^{3} \mathcal{R}_j(\theta_j, d_j) \quad (17)$$

with $\mathcal{R}_j$ a 1D-correlation family taken from Table 4, $\theta_j$ the value of $\theta$ in the j-th dimension, and $d_j = p_j - p'_j$ the distance between $\mathbf{p}$ and $\mathbf{p}'$ in the j-th dimension.

The construction of a kriging surrogate model is based on $Q$ evaluations of the function $\mu$. Thus, from $Q$ inputs $(\mathbf{p}^{(k)})_{k \in [1, Q]}$, $Q$ outputs $(\mu^{(k)} = \mu(\mathbf{p}^{(k)})_{k \in [1, Q]}$ are obtained. From this set, the regression matrix $\mathbf{G}$ of coefficients $G_{ij} = g_j(\mathbf{p}^{(i)})$ and the regression matrix $\mathbf{R}$ of coefficients $R_{ij} = \mathcal{R}(\theta, \mathbf{p}^{(i)}, \mathbf{p}^{(j)})$ are computed. The
and the variance by

$$ \sigma_{\lambda(\mathbf{p})}^2 = \sum_{k=1}^{P-1} \mu_k(\mathbf{p})^2 \|\Psi_{\pi_k}\|^2 = \sum_{k=1}^{P-1} (\mathbf{g}^{(k)}\beta^{(k)} + Z^{(k)}(\mathbf{p}))^2 \|\Psi_{\pi_k}\|^2 \quad (22) $$

The Sobol indices [34] also depend on the parametric vector and can also be directly deduced from the formulation as they depend on the PCE coefficients [29]. Let $V(\lambda(\mathbf{p}))$ denotes the variance of $\lambda$ at $\mathbf{p}$. The first Sobol index $S_i$ related to the random variable $i \in \{k_{\alpha\alpha}, c_{\alpha}, k_h, c_h, s_p\}$ is defined as [34]:

$$ S_i(\mathbf{p}) = \frac{V_i(\lambda(\mathbf{p}))}{V(\lambda(\mathbf{p}))} \quad (23) $$

where $V_i(\lambda(\mathbf{p})) = V(E(\lambda(\mathbf{p})|X_i)$ is obtained from the PCE coefficients [29]:

$$ V_i(\lambda(\mathbf{p})) = \sum_{k \in \mathbf{v}_i} \mu_k^2 \|\Psi_{\pi_k}\|^2 \quad (24) $$

with $\mathbf{v}_i$ the set of multivariate indices composed only of polynomials related to the variable $i$.

## 4.5 Automatic construction criterion of the kriging surrogate models

As the hybrid surrogate model requires one kriging surrogate model per PCE coefficient, it requires many kriging surrogate models. Usually, the choice the correlation and regression function to be used is based on a trial-and-error process or based on previous knowledge [23, 25, 27]. However, this approach is hardly applicable in the context of the hybrid surrogate model as many kriging surrogate models must be constructed. So different criterion for an automatic selection of the best kriging properties have been compared. Two options are possible in terms of construction of the hybrid surrogate model, either all PCE coefficients are approximated with the same kriging properties, or each coefficient has its own kriging properties. With the first option, error criterion based on the error of the hybrid surrogate model are preferred, whereas with the second option it is necessary to use a criterion for each PCE coefficient (and so for each kriging surrogate model).

Two methods are adopted for the computation of the error criterion. The first one is based on the comparison of the predicted values from a surrogate model and reference values. These reference values are obtained from a set of $N_v$ validations point (and so it requires additional simulations). The second option is based on the well known Leave-One-Out (LOO) criterion [30]. This approach does not require additional simulations and exploits the training set composed of $Q$ points. From this set, one point $i$ is removed and a new surrogate model is constructed with the new training set. The prediction of the new surrogate model is compared to the reference value. This process is repeated for each point of the training set.

The different cases retained for the automatic construction of the kriging surrogate models are:

* case 1: it is a baseline case which correspond to the case where all the kriging have the same properties, and a classical choice is taken, namely a Matérn 5/2 correlation function and a constant regression part. This option is extremely convenient in terms of implementation and in terms of complexity. This case will be to see how much improvement one can get by using advanced criterion.

* case 2: all the kriging surrogate models have the same properties, and the correlation and regression functions are selected to minimise the maximum of the absolute relative error of the hybrid surrogate model on a set of reference validation points. This error can be written $e_{r,glob} = \max_{k \in [1, N_v]} \left| \frac{\lambda^{(k)} - \tilde{\lambda}^{(k)}}{\lambda^{(k)}} \right|$, with $\lambda^{(k)}$ the exact value of the eigenvalue at the k-th point of the validation set and $\tilde{\lambda}^{(k)}$ the hybrid surrogate model prediction at the k-th point of the validation set.

* case 3: all the kriging surrogate models have the same properties, and the correlation and regression functions are selected to minimise the average of the absolute relative error of the hybrid surrogate model on a set of reference validation points. This error can be written $\tilde{e}_{r,glob} = \frac{1}{N_v} \sum_{k=1}^{N_v} \left| \frac{\lambda^{(k)} - \tilde{\lambda}^{(k)}}{\lambda^{(k)}} \right|$, with $\lambda^{(k)}$ the exact value of the eigenvalue at the k-th point of the validation set and $\tilde{\lambda}^{(k)}$ the hybrid surrogate model prediction at the k-th point of the validation set.

* case 4: each kriging surrogate model has its own properties. The correlation and regression function of a kriging surrogate model for a given PCE coefficient are chosen to minimise the maximum of the absolute relative LOO error on the PCE coefficient. This error can be written: $e_{LOO,coeff} = \max_{k \in [1,Q]} \left| \frac{\mu_j^{(k)} - \tilde{\mu}_j^{(\sim k)}}{\mu_j^{(k)}} \right|$, with $\mu_j^{(k)}$ the reference value of the j-th PCE coefficient at the k-th point of the training set and $\tilde{\mu}_j^{(\sim k)}$ the prediction of the kriging surrogate model constructed by removing the k-th point.
* case 5: each kriging surrogate model has its own properties. The correlation and regression function of a kriging surrogate model for a given PCE coefficient are chosen to minimise the average of the absolute relative LOO error on the PCE coefficient. This error can be written: $\bar{e}_{LOO,coeff} = \frac{1}{Q} \sum_{k=1}^{Q} \left| \frac{\mu_j^{(k)} - \tilde{\mu}_j^{(k)}}{\mu_j^{(k)}} \right|$, with $\mu_j^{(k)}$ the reference value of the j-th PCE coefficient at the k-th point of the training set and $\tilde{\mu}_j^{(k)}$ the prediction of the kriging surrogate model constructed by removing the k-th point.

* case 6: all the kriging surrogate models have the same properties, and the correlation and regression functions are selected to minimise the maximum of the absolute relative LOO error of the hybrid surrogate model. This error can be written: $e_{LOO,glob} = \max_{k \in [1,Q]} \left| \frac{\lambda^{(k)} - \tilde{\lambda}^{(\sim k)}}{\lambda^{(k)}} \right|$, with $\lambda^{(k)}$ the exact value of the eigenvalue at the k-th point of the training set and $\tilde{\lambda}^{(\sim k)}$ the prediction of the hybrid surrogate model constructed by removing the k-th point of the training set.

* case 7: all the kriging surrogate models have the same properties, and the correlation and regression functions are selected to minimise the average of the absolute relative LOO error of the hybrid surrogate model. This error can be written: $\bar{e}_{LOO,glob} = \frac{1}{Q} \sum_{k=1}^{Q} \left| \frac{\lambda^{(k)} - \tilde{\lambda}^{(\sim k)}}{\lambda^{(k)}} \right|$, with $\lambda^{(k)}$ the exact value of the eigenvalue at the k-th point of the training set and $\tilde{\lambda}^{(\sim k)}$ the prediction of the hybrid surrogate model constructed by removing the k-th point of the training set.

# 5 Results

In this section, the hybrid surrogate model is employed to surrogate the stability behaviour of the airfoil, i.e. to predict its complex eigenvalues. The first part of the section is dedicated to the construction of the hybrid surrogate model. More precisely, the PCE properties are first chosen based on a convergence study. Then, the different automatic construction criterion for the construction of the kriging surrogate models of the PCE coefficients are compared. Finally, the final surrogate model is validated. The second part of the section is dedicated to the exploitation of the surrogate model to make an extensive study of the stability of the airfoil considering the uncertainties. The stochastic eigenvalues are predicted over the parametric space in a first time. In a second time, a sensitivity analysis using the Sobol indices is performed.

## 5.1 Hybrid surrogate model construction with automatic criterion

### 5.1.1 PCE construction

The goal of this section is to set up the PCE properties, i.e. to chose the DoE size and the PCE size with the chaos order $o_{pce}$ and the truncation norm $q$. A convergence study over the DoE size is done for different PCE sizes. The different DoE sizes are [50, 75, 100, 125, 150, 175, 200, 250, 300, 350, 400, 500], and for each size 10 simulations are done to get an estimation of the robustness of PCE properties to the DoE. The different chaos order considered are [1, 2, 3, 5] and two truncation norms are considered, namely 1 and 0.5. For each case, the PCE predictions are compared to 500 reference values from a validation set, and the relative error is computed. As it is required to have more training points than PCE coefficients, for each PCE basis, the error is computed only when there are enough points in the training set. In Figure 3, the evolution of the average relative errors w.r.t. the DoE size are given for each eigenvalue and for the different PCE sizes. The errorbar limits represent the average values $\pm$ the standard deviation.

From Figure 3, one can see that the error decreases when the DoE sizes increase until it reaches a plateau. Also, it can be noticed that the error dispersion tends to decrease when the DoE size increases (see the errorbars that are less spread). The final error is strongly related to the PCE basis size. Indeed, cases with the highest error have the smaller bases (7 terms for ($o_{pce}$ = 1, $q$ = 0.5)) than for basis with more terms (84 terms for ($o_{pce}$ = 3, $q$ = 1)) or even 462 terms for ($o_{pce}$ = 5, $q$ = 1). However, the largest the PCE basis, the more points in the DoE are needed to compute the PCE coefficients. However, in all cases, the level of error remains low as well as the dispersion of the error. Indeed, the average relative error is always inferior to $10e-2$ for the real parts and inferior to $10e-3$ for the angular frequencies.

Finally, a chaos order of $m = 2$ with a truncation norm of $q = 1$ are taken and a DoE of 100 points is chosen. This choice is based on the fact that the error level is satisfactory for the eigenvalue prediction. Indeed, it gives an error of about 2e-3 for $Re(\lambda_1)$, 2.5e-4 for $Re(\lambda_2)$, 2.5e-4 for $\omega_1$ and around 2e-4 for $\omega_2$. The case ($o_{pce}$ = 3, $q$ = 1) gives lower error but represents a PCE basis of 84 terms, whereas the case ($o_{pce}$ = 2, $q$ = 1) represents a basis of 28 terms. Considering that one kriging surrogate model is constructed per PCE coefficients, a compromise between the PCE basis size and the PCE accuracy is done.
<figure><img src="image_4.png" /><figcaption>Figure 3: Evolution of the error for different PCE size versus the size of the experimental design for the real (a,b) and imaginary (c,d) parts of the first (a,c) and second (b,d) eigenvalues</figcaption></figure>

5.1.2 Comparison of the automatic construction criterion

The different criterion for the automatic construction of the kriging surrogate models is presented here. The
change in the kriging properties (correlation and regression) has a strong impact on its accuracy and so these
properties must be well tuned. This can be done by hand, but in a context where numerous kriging surrogate
models must be created, i.e. one per PCE coefficient, this strategy cannot be used anymore and automatic
construction criterion must be investigated. Seven different strategies are considered and compared.

Thus, for each eigenvalue real part and angular frequency, hybrid surrogate models are constructed with different DoE sizes for the kriging part, namely [10, 25, 50, 100, 150, 200, 500] points. For each DoE size, 20 LHS are generated. The final error of each hybrid surrogate models is computed by computing the relative error of the hybrid surrogate model prediction at reference values from a validation set composed of 50 × 125 = 6250 points. From these 6250 values, one get the average relative error $\bar{e}_r$ and the standard deviation of this error $\sigma_{e_r}$. For each hybrid surrogate model and for each DoE size, the average and the standard deviation over the 20 LHS are computed. They are displayed in Figure 4 for the evolution of the average relative error and in Figure 5 for the standard deviation of the relative error. For each figure, on the left column the average of the quantity of interest is given (i.e. $E[\bar{e}_r]$ or $E[\sigma_r]$) and in the right column, the standard deviation of the quantity is given (i.e. $Std[\bar{e}_r]$ or $Std[\sigma_r]$). The different automatic strategies for the kriging construction can be compared in terms of performances, convergence speed and robustness.

As a first general comment, the evolutions of the errors are similar and the average and the standard deviation
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 4: Average relative error - mean (left) and std (right) - Evolution of the average error of the average relative error on 20 LHS</figcaption></figure>
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Figure 5: Std relative error - mean (left) and std (right) - Evolution of the average error of the average relative error on 20 LHS</figcaption></figure>
of the relative errors decrease in mean and variance when the DoE size increases. Looking at the evolution of
the average relative error $\bar{e}_r$ in Figure 4 and of its standard deviation $\sigma_r$ in Figure 5, some differences between
the different cases is visible. In all cases the indicators that are global and that are based on the maximum
error computed (orange and light blue lines) always underperform compared to the others. Indeed, the level
of error is always higher and the quality of the hybrid surrogate models varies between two different DoE. For
example, for $E[\bar{e}_r^{Re}(\lambda_1)]$, for a DoE of 25 points, the mean of the average relative error is two times higher than
for the other cases. Similarly, the standard deviation of the average relative errors $\bar{e}_r$ are also always higher for
these two cases, and sometimes they even reach a plateau (see orange line in Figure 4(d)). It demonstrates that
the quality of the hybrid surrogate models vary a lot from one DoE to another and show a lack of robustness.
On the other sides, the indicators that are based on the average value tend to be better with lower error levels
(see the yellow, green and red curves). They also have performances that are more constant than the indicators
based on a max. More particularly, the criterion that is based on the average global relative error (yellow
curves) present the best performances in terms of average error and levels of deviations. However, the criterion
based on a LOO process on each PCE coefficients and on the global LOO present similar performances. As
a first conclusion, indicators that are based on the average error in the kriging construction rather than on
the maximum absolute error have better performances in terms of error levels and robustness to the DoE.
Compared to a classic construction (black blue line) with constant and predetermined kriging characteristics,
one can see that the error is slightly reduced with these advanced criterion based on an average of the error.
But the improvement in terms of dispersion is important, especially for the real parts. So using constant and
predetermined properties is very efficient in terms of numerical time (no need to optimise the kriging properties)
and implementation, but the quality of the hybrid surrogate model vary a lot.

A last comment, would be that the criteria based on the average of the global relative error gives the best performances. However, it requires an additional set of validation points, which is not always available. In this context, criterion based on LOO are preferable, especially as they demonstrate similar performances. Using the LOO process for each PCE coefficient, and so having different kriging properties for each PCE coefficient, gives better performances than using the same kriging properties and a global LOO process. However, the process of the selection of the best kriging properties for each PCE coefficient takes much more time and should be considered. So using the same kriging properties for all PCE coefficients tends to be the best compromise in terms of prediction quality and computational time.

5.1.3 Validation of the surrogate model

The final hybrid surrogate model is constructed with 100 points in the Kriging DoE and is constructed using
the average global relative error. It is finally validated by comparing its prediction to reference values taken
from a set of 5000 validation points. The comparison for each hybrid surrogate model is displayed in Figure 6,
where the validation points are in red and the reference points in black. The good accuracy of the different
surrogate models is clearly visible. The view in the complex plan is given in Figure 7. It is clear that the
stability behaviour of the airfoil is accurately predicted with the different surrogate models, so they can be
employed for large parametric, stochastic and sensitivity studies.

## 5.2 Eigenvalue prediction

The hybrid surrogate model can now be used and exploited to analyse the stability behaviour of the airfoil.
The evolution of the average of the real parts and of the angular frequencies of the two modes are given in
Figures 8 and 9 respectively. The evolutions of the standard deviations are given in the Figures 10 and 11
for the real parts and angular frequencies respectively. The first mode is the mode that can become unstable
(positive real part) and the stability limit of the mode is represented with a black line. It corresponds to the
case where the average real part is null. The Table 5 summarizes the minimum and maximum values observed
on each surface for each flow speed.

Considering the first mode, the impact of the wind on the stability limit is clear. Under a speed of 20 m/s,
the mode is always stable. At 25 m/s and above, the mode becomes more and more unstable. The stability
limit is mostly driven by the distance of the elastic axis to the midchord $\bar{a}$ and the semichord length $b$ has a
low impact (see the black line that is almost horizontal). When the flow speed increases, the stability limit is
reached for higher values of the distance of the elastic axis to the midchord: at 25 m/s, the system is unstable
when $\bar{a}$ is inferior to about 0.65, but at 40 m/s the system is unstable when $\bar{a}$ is inferior to about 0.55. When
the semichord length $b$ increases, the stability limit is slightly moved and the system is a bit more stable, i.e.
lower values of $\bar{a}$ must be reached to have an unstable airfoil. Also, at low flow speed, the real part is almost
constant, but when the flow speed increases, the real part varies and reaches its highest values when $\bar{a}$ and $b$ are
minimal. Considering the corresponding standard deviation (see Figure 11(a)), it is almost always null when
<figure><img src="image_21.png" /><figcaption>Figure 6: Comparison of the hybrid surrogate model predictions (●) and reference values (o) of the first (a,c) and second (c,d) eigenvalues at the validation points for the real part (a,b) and the angular frequency (c,d)</figcaption></figure>

the flow speed is low and so the uncertainties of the model don't impact much the real part of the eigenvalue. However, when the flow speed increases, the standard deviation also increases. It tends to remain low in the stable area (see the blue areas above the stability limit), but when the system is unstable the standard deviation tends to increase showing that the uncertainties impact the eigenvalue. Larger dispersions are also observed around the stability limit.

Looking at the angular frequency of the first mode (see Figure 9(a)), it increases when the flow speed increases. At very low speeds (10 m/s and under), the angular frequency increases with *b* and does not really depend on *ā*. At 15 m/s and above, the angular frequency is maximal for the lowest values of *ā* and *b* and decreases when *ā* or when *b* increases. Looking at the corresponding standard deviation (see Figure 11(a)), it is always low and constant (around 0.1) for wind velocities of 10 m/s or under. In terms of standard deviation (see Figure 11), at low velocities, i.e. at 10 m/s or under, the standard deviation is constant and low (between 0.11 and 0.18 according to Table 5). So the impact of the uncertainties on the angular frequency is low and constant. It has a higher impact on the angular frequency than on the real part (where the standard deviation was equal to 0). When the flow speed increases, the standard deviation increases as well up to 0.30. At 15 m/s and 20 m/s, the location of the highest standard deviations depend mostly on *ā* and the maximum is observed for lower *ā*. At 20 m/s, a peak of the standard deviation at 0.57 around ($ā = -0.7, b = 0.135$) is observed. It means that a large dispersion of the angular frequency is observed and that the model uncertainties have a strong impact. This increase in the standard deviation is coupled to an increase of the average angular frequency. When the flow speed increases again, the maximum of the standard deviation tends to move towards larger *b* and *ā*, and the larger values are always observed in the stable area of the mode. At flow speeds that are larger than 35 m/s, the standard deviation decreases (a maximum of 0.39 at 35 m/s and of 0.3 at 40 m/s) and tends to be more
<figure><img src="image_22.png" /><figcaption>Figure 7: Comparison of the hybrid surrogate model predictions (●) and reference values (●) of the eigenvalues at the validation points in the complex plan</figcaption></figure>

constant (vary between 0.12 and 0.30).

Considering the second mode, it is always stable (see Figure 8(b) and Table 5) as the real part of the eigenvalue is always negative. The behaviour of the second mode is somewhat complementary to the behaviour of the first mode: when the average real part of the first eigenvalue increases, the average real part of the second value decreases. At low flow speed, the average real part does not vary a lot (between -8 and -4 at 0 m/s). When the flow speed increases, the variation of the average real part increases (it varies between -11 and -33 at 40 m/s) and becomes more and more negative (the maximum average real part at 40 m/s is equal to -11.57, whereas it was equal to -4.32 at 0 m/s). The minimum average real part is always reached for the minimum values of *b* and *ā*. The standard deviation of the real part of the second mode is higher than the one of the first mode. Indeed, the minimum standard deviation is equal to 0.1 and reaches 0.22 at 30 m/s, whereas the maximum value of the minimum standard deviation for the real part of the first eigenvalue was 0.08. It demonstrates that the model uncertainties have a higher influence on the second eigenvalue real part than on the first one. The maximum standard deviation tends to increase with the flow speed (up 0.77 at 40 m/s). The variations of the standard deviation depend a lot of the flow speed. Indeed, at 15 m/s, the minimum standard deviation is observed for *b* = 0.16, whereas at 30 m/s, the minimum standard deviation is observed at *ā* = -0.53. A peak of the standard deviation is observed at 40 m/s for the lower values of *b* and maximum values of *ā*. These different observations demonstrate that the increase in the flow speed tends to increase the dispersion of the real part of the second eigenvalue. However, the impact on this dispersion depends also on *b* and *ā* and the location of the largest dispersion depends on the flow speed, the semichord length and the distance of the elastic axis to the midchord. These evolutions are complex, which illustrates the necessity to conduct such uncertainty propagation and quantification studies.

Looking at the evolution of the second angular frequency (see Figure 9(b), Figure 11(b) and Table 5), it tends to decrease when the flow speed increases. Indeed, at 0 m/s, $\omega_2$ varies between 30.62 rad/s and 40.67 rad/s, whereas it varies between 13.81 rad/s and 24.46 rad/s at 40 m/s. The evolution of the average angular frequency $\omega_2$ is very smooth. At lower flow speeds, the maximum is reached for lower values of *b* and *ā*, but when the flow speed increases, the maximum is reached for the highest values of *b* and the lowest values of *ā*. In terms of standard deviation, the second angular frequency is the most sensitive to the model uncertainties as the standard deviations are the highest (always superior to 0.46 and up to 1.78). At low flow speeds, the standard deviation is almost constant (between 0.62 and 0.77 at 0 m/s). But when the flow speed increases, the standard deviation varies more. The minimum decreases and the maximum standard deviation increases. For example, at 40 m/s, it varies between 0.46 and 1.78. The location of the maximum varies with the flow speed and shows complex evolutions. Indeed, at 20 m/s, an increase of the standard deviation is observed for $b \in [0.13, 0.14]$ and $ā \in [-0.77, -0.65]$. At 35 m/s, this maximum is less spread and is reached for higher values of $ā$ and at 30 m/s, this maximum is reached for $ā = -0.55$ and doesn't depend on *b* any more. Finally, at 35 m/s and 40 m/s, higher values of the standard deviation are reached and they are located for low values of *b* and high values of $ā$. This demonstrates that the impact of the model uncertainties of the angular speeds depends a lot on the flow speed, but also on the semichord length *b* and the distance of the elastic axis to the midchord $ā$.

Finally, the probability of instability depending on the flow speed *V*, the semichord length *b* and the distance of the elastic axis to the midchord $ā$ is given in Figure 12. For different values of the flow speed and over $(b, ā)$,
<table><caption>Table 5: Minimum and maximum values of the average and standard deviation over the considered values of *b* and *ā* for different flow speeds and for the different modes</caption><thead><tr><th rowspan="2">Flow speed (m/s)</th><th colspan="4">Re(λ₁)</th><th colspan="4">ω₁</th></tr><tr><th>min mean</th><th>max mean</th><th>min std</th><th>max std</th><th>min mean</th><th>max mean</th><th>min std</th><th>max std</th></tr></thead><tbody><tr><td>0</td><td>-0.40</td><td>-0.04</td><td>4.5E-3</td><td>0.08</td><td>7.76</td><td>12.29</td><td>0.11</td><td>0.18</td></tr><tr><td>5</td><td>-1.21</td><td>-0.55</td><td>6.1E-3</td><td>0.08</td><td>7.78</td><td>12.31</td><td>0.10</td><td>0.18</td></tr><tr><td>10</td><td>-2.14</td><td>-0.67</td><td>6.2E-3</td><td>0.15</td><td>7.88</td><td>14.88</td><td>0.11</td><td>0.18</td></tr><tr><td>15</td><td>-2.83</td><td>-1.10</td><td>0.02</td><td>0.37</td><td>7.98</td><td>18.75</td><td>0.12</td><td>0.30</td></tr><tr><td>20</td><td>-2.61</td><td>-0.48</td><td>0.03</td><td>0.63</td><td>8.12</td><td>24.94</td><td>0.12</td><td>0.57</td></tr><tr><td>25</td><td>-2.73</td><td>3.53</td><td>0.04</td><td>0.59</td><td>7.92</td><td>29.23</td><td>0.16</td><td>0.42</td></tr><tr><td>30</td><td>-4.01</td><td>6.31</td><td>0.06</td><td>0.54</td><td>7.43</td><td>31.44</td><td>0.15</td><td>0.65</td></tr><tr><td>35</td><td>-4.79</td><td>9.50</td><td>0.08</td><td>0.60</td><td>7.89</td><td>33.18</td><td>0.14</td><td>0.39</td></tr><tr><td>40</td><td>-4.03</td><td>12.19</td><td>0.06</td><td>0.86</td><td>8.40</td><td>34.70</td><td>0.12</td><td>0.30</td></tr></tbody><thead><tr><th rowspan="2"></th><th colspan="4">Re(λ₂)</th><th colspan="4">ω₂</th></tr><tr><th>min mean</th><th>max mean</th><th>min std</th><th>max std</th><th>min mean</th><th>max mean</th><th>min std</th><th>max std</th></tr></thead><tbody><tr><td>0</td><td>-8.25</td><td>-4.62</td><td>0.18</td><td>0.31</td><td>30.62</td><td>40.67</td><td>0.62</td><td>0.77</td></tr><tr><td>5</td><td>-8.88</td><td>-5.48</td><td>0.19</td><td>0.32</td><td>30.45</td><td>39.81</td><td>0.62</td><td>0.80</td></tr><tr><td>10</td><td>-9.53</td><td>-6.20</td><td>0.20</td><td>0.33</td><td>30.02</td><td>38.43</td><td>0.64</td><td>0.84</td></tr><tr><td>15</td><td>-10.23</td><td>-6.72</td><td>0.20</td><td>0.40</td><td>27.86</td><td>35.31</td><td>0.66</td><td>0.99</td></tr><tr><td>20</td><td>-13.68</td><td>-8.20</td><td>0.18</td><td>0.44</td><td>24.64</td><td>33.06</td><td>0.67</td><td>1.27</td></tr><tr><td>25</td><td>-20.05</td><td>-8.57</td><td>0.18</td><td>0.57</td><td>23.08</td><td>29.98</td><td>0.54</td><td>1.13</td></tr><tr><td>30</td><td>-24.55</td><td>-8.09</td><td>0.22</td><td>0.51</td><td>22.57</td><td>27.71</td><td>0.51</td><td>1.16</td></tr><tr><td>35</td><td>-29.20</td><td>-9.19</td><td>0.21</td><td>0.41</td><td>19.80</td><td>26.12</td><td>0.46</td><td>1.35</td></tr><tr><td>40</td><td>-32.94</td><td>-11.57</td><td>0.10</td><td>0.77</td><td>13.81</td><td>24.46</td><td>0.46</td><td>1.78</td></tr></tbody></table>

the percent of unstable cases considering the variation of the random parameters is displayed. One can clearly see the impact of the different parameters on the airfoil stability. Thus, if the flow speed is inferior to 20 m/s, then the airfoil is always stable. Similarly, if *ā* is equal to -0.5 then the system is always stable. In terms of impact on the bifurcation point, an increase in *b* tends to stabilise the system, i.e. the Hopf bifurcation point is pushed towards higher flow speeds. Similarly, an increase in *ā* tends to stabilise the system. One can also see that for *V* = 25 m/s, the impact of the uncertainties on the bifurcation point is important as the transition zone between 0 % of instability case and 100 % of instability cases is larger than for higher flow speeds. Moreover, the higher the flow speed, the smaller this transition zone and so the smaller the impact of the random parameters on the bifurcation point.

## 5.3 Sensitivity analysis

In the previous section, the stochastic output was analysed in terms of average and variance over the parametric space. In this section, the impact of each random parameter is assessed in details thanks to a sensitivity based on a Sobol analysis. As a reminder, a Sobol index is an indicator on the influence of a parameter on a quantity of interest. Thus, if the Sobol index $S_i$ related to a parameter *i* is close to 1, then this parameter has a strong influence on the output. On the opposite, if the Sobol index is close to 0, then this parameter *i* has a low influence. The Sobol indices are a common and convenient way to rank and compare the influence of different parameters on a quantity of interest. The first order Sobol indices are given in Figure 13 and in Figure 14 for the real part of the eigenvalues and the angular frequencies, respectively. From a look at the figures, one can see that the Sobol indices vary a lot with the flow speed, are very different from one parameter to another and strongly vary with *b* and *ā*. These evolutions are complex and demonstrate a complex evolution of the role played by each parameter.

Considering the real part of the first eigenvalue, some parameters have always a negligible impact, namely $k_α$ and $c_h$ as their Sobol indices are always close to 0. $c_α$ has a strong influence for some values of *b* and *ā* at flow speeds inferior to 10 m/s. $c_h$ has a medium impact on the first eigenvalue regardless of the flow speed. At some values of *b* and *ā*, $c_h$ has a strong influence (Sobol index equal to about 1). For 30 m/s and above, it corresponds to a sort of arc of circle. The latter is always located just above the stability limit of the mode. Finally, $s_p$ is the parameter with the highest impact on the real part of the first mode, and so on the instability behaviour of the airfoil. There is always one area where $s_p$ has almost no influence (Sobol index close to 0), which is complementary to the influence of $c_h$. For flow speeds inferior to 20 m/s, the real part is driven by $s_p$,
<figure><img src="image_23.png" /><figcaption>Figure 8: Evolution of the average real part of the first (a) and second (b) modes versus the semichord *b* and the distance of the elastic axis to the midchord *ā* for different values of the flow speed - (–) stability limit</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure 9: Evolution of the average angular frequency of the first (a) and second (b) modes versus the semichord $b$ and the distance of the elastic axis to the midchord $\bar{a}$ for different values of the flow speed - (–) stability limit</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 10: Evolution of the standard deviation of the real part of the first (a) and second (b) modes versus the semichord *b* and the distance of the elastic axis to the midchord *ā* for different values of the flow speed - (–) stability limit</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 11: Evolution of the standard deviation of the angular frequency of the first (a) and second (b) modes versus the semichord *b* and the distance of the elastic axis to the midchord *ā* for different values of the flow speed - (–) stability limit</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 12: % of unstable cases over the considered values of b and $\bar{a}$ for different flow speeds</figcaption></figure>

$k_h$ and $c_\alpha$. But for flow speeds at 30 m/s or above, there is almost only $s_p$ that has an influence.

For the real part of the second eigenvalue (see Figure 13(b)), $k_\alpha$ and $c_\alpha$ have no influence as their Sobol indices are always null for all flow speeds, all semichord length and all distance of the elastic axis to the semichord. $k_h$ has some influence on the real part only when the flow speed is superior to 20 m/s. At low flow speeds (0 m/s and 10 m/s), $c_h$ and $s_p$ are the two parameters that have an influence on the real part of the eigenvalue. Their influence is similar (Sobol indices around 0.5) and is constant over $b$ and $\bar{a}$. For larger flow speed, $k_h$ also becomes influential. The influence of these three parameters depend on $b$ and $\bar{a}$ and the flow speed. For $s_p$, there is a line at constant $\bar{a}$ for which the Sobol index $S_{s_p}$ is equal to 0. The location of this line changes with the flow speed and is reached for higher $\bar{a}$ (see the line at $\bar{a} = -0.7$ at 20 m/s, at $\bar{a} = -0.55$ at 30 m/s, etc.). When $S_{s_p}$ reaches a zero value (and so $s_p$ has no influence), $c_h$ becomes the parameter with the highest influence (see the lines with higher Sobol index). Thus in this context, to control the variability of the second eigenvalue, the variability of $c_h$ and $s_p$ must be controlled at low flow speeds for higher flow speeds $k_h$ must also be considered. However, depending on the value of $b$ and $\bar{a}$, $s_p$ might have no influence.

Considering the first angular frequency (see Figure 14(a)), the influential parameters are different. Indeed, $c_\alpha$ and $c_h$ have almost no impact on $\omega_1$ (Sobol indices always equal to 0). For the other parameters, their influence depend on the flow speed. Thus, $k_\alpha$ has a high influence for low flow speeds. At 0 m/s, it is even the only parameter that impacts the first angular frequency. But when the flow speed increases, its influence decreases and becomes almost null at high flow speeds. At medium flow speed, $\omega_1$ is mostly driven by the airfoil span $s_p$, except at $\bar{a} \simeq -0.5$ where $k_\alpha$ drives the angular frequency. Finally, when the flow speed increases again, $k_h$ becomes the parameter that drives $\omega_1$. So for the first angular frequency, the parameter that has the highest impact depends on the flow speed. It justifies the interest of performing sensitivity analysis to identify and assess this important shifts in terms of main parameters with high influence.

For the second angular frequency (see Figure 14(b)), the impact of $k_\alpha$, $c_\alpha$ and $s_p$ is almost null as their Sobol indices are equal to zero. On the other side, the Sobol indices of $k_h$ and $s_p$ are almost constant over $V$, $b$ and $\bar{a}$ and equal to about 0.5 so their influence are similar.

As a conclusion, $c_\alpha$ has no influence on the complex eigenvalues and could be removed from future studies and considered as deterministic. $k_\alpha$ impacts only the first angular for $V \le 10$ m/s, but otherwise it doesn't impact the other parameters. Similarly, $c_h$ influences only the real part of the second eigenvalue. $k_h$ impacts slightly the real parts, impacts $\omega_1$ only for $V \ge 30$ m/s and always has an influence of $\omega_2$. Finally, $s_p$ has a strong influence on the real parts of the eigenvalues and $\omega_2$, but it impacts $\omega_1$ only for $20 \text{ m/s} \le V \le 30 \text{ m/s}$.
a

<figure><img src="image_28.png" /></figure>

ā

<figure><img src="image_29.png" /><figcaption>Figure 13: First order Sobol index of the real part of the first (a) and second (b) eigenvalues versus the semichord *b* and the distance of the elastic axis to the midchord ā for different values of the flow speed</figcaption></figure>
<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /><figcaption>Figure 14: First order Sobol index of the angular frequency of the first (a) and second (b) modes versus the semichord *b* and the distance of the elastic axis to the midchord *ā* for different values of the flow speed</figcaption></figure>
# 6 Conclusion

The objective of the present work was to study the stability behaviour of a stochastic airfoil model when considering a parametric variation of several parameters. The stability of the 2-dof airfoil is assessed based on the analysis of the real part of the complex eigenvalues with a CEA. To limit the computational time, a hybrid surrogate model that couples polynomial chaos expansion and kriging is employed. It allows to consider both kind of parameters and perform efficiently a parametric study of the stochastic system. In the formulation, each PCE coefficient is approximated with a kriging surrogate models. Different criterion are compared to construct efficiently the kriging surrogate models. It is shown that error criterion that minimises the average relative error, instead of the maximum, have better performances in terms of accuracy, convergence and robustness towards the DoE. Once the hybrid surrogate models for each real part and each imaginary part of the eigenvalues are constructed, they are exploited to perform a deep parametric study of the stochastic system and investigate the stability behaviour of the airfoil in the presence of uncertainties. Finally, an extensive sensitivity analysis based on Sobol indices is done. Deep insights in the role played by each random parameter are gained and it is shown that the parameters that drive the airfoil eigenvalues change a lot over the parametric domain.

# Acknowledgments

J.-J. Sinou acknowledges the support of the Institut Universitaire de France. E. Denimal acknowledges the support of Rennes Metropole.

# References

[1] J. R. Wright and J. E. Cooper. *Introduction to aircraft aeroelasticity and loads*, volume 20. John Wiley & Sons, 2008.

[2] D. Hodges and G. Pierce. *Introduction to structural dynamics and aeroelasticity*, volume 15. Cambridge University Press, 2011.

[3] V. Mukhopadhyay. Historical Perspective on Analysis and Control of Aeroelastic Responses. *Journal of Guidance, Control, and Dynamics*, 26(5):673–684, September 2003.

[4] J. Mottershead and Y. Ram. Inverse eigenvalue problems in vibration absorption: passive modification and active control. *Mechanical systems and signal processing*, 20(1):5–44, 2006.

[5] Y.H. Zhao. Stability of a two-dimensional airfoil with time-delayed feedback control. *Journal of Fluids and Structures*, 25:1–25, 2009.

[6] K. V. Singh. Active aeroelastic control with time delay for targeted flutter modes. *Aerospace science and technology*, 43:281–288, 2015.

[7] Y.H. Zhao. Flutter instability and active aeroelastic control with time delay for a two-dimensional airfoil. *European Journal of Mechanics - A/Solids*, 92:104465, 2022.

[8] C. J Roy and W. Oberkampf. A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing. *Computer methods in applied mechanics and engineering*, 200(25-28):2131–2144, 2011.

[9] M. McKay, R. Beckman, and W. Conover. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. *Technometrics*, 42(1):55–61, 2000.

[10] D. Poirel and S.J. Price. Random binary (coalescence) flutter of a two-dimensional linear airfoil. *Journal of Fluids and Structures*, 18(1):23–42, 2003.

[11] RA Ibrahim, DM Beloiu, and CL Pettit. Influence of joint relaxation on deterministic and stochastic panel flutter. *AIAA journal*, 43(7):1444–1454, 2005.

[12] N. Wiener. The homogeneous chaos. *American Journal of Mathematics*, 60(4):897–936, 1938.

[13] D. Xiu and G. Karniadakis. The Wiener-Askey polynomial chaos for stochastic differential equations. Technical Report 2, 2002.

[14] J. Witteveen, S. Sarkar, and H. Bijl. Modeling physical uncertainties in dynamic stall induced fluid-structure interaction of turbine blades using arbitrary polynomial chaos. *Computers & structures*, 85(11-14):866–878, 2007.

[15] C. Pettit and P. Beran. Polynomial chaos expansion applied to airfoil limit cycle oscillations. In *45th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics & Materials Conference*, page 1691, 2004.

[16] S. Sarkar, JAS Witteveen, A Loeven, and H Bijl. Effect of uncertainty on the bifurcation behavior of pitching airfoil stall flutter. *Journal of Fluids and Structures*, 25(2):304–320, 2009.

[17] M Ghommem, MR Hajj, and AH Nayfeh. Uncertainty analysis near bifurcation of an aeroelastic system. *Journal of Sound and Vibration*, 329(16):3335–3347, 2010.

[18] D. Jones, M. Schonlau, and W. Welch. Efficient Global Optimization of Expensive Black-Box Functions. Technical report, 1998.

[19] C. Cortes and V. Vapnik. Support-vector networks. *Machine learning*, 20(3):273–297, 1995.

[20] I. Tartaruga, J. Cooper, G. Georgiou, and H. Khodaparast. Flutter uncertainty quantification for the s4t model. In *55th AIAA Aerospace Sciences Meeting*, page 1653, 2017.

[21] M. Rumpfkeil and P. Beran. Multi-fidelity surrogate models for flutter database generation. *Computers & Fluids*, 197:104372, 2020.

[22] A. Thelen, L. Leifsson, and P. Beran. Multifidelity flutter prediction using regression cokriging with adaptive sampling. *Journal of Fluids and Structures*, 97:103081, 2020.

[23] E. Denimal, L. Nechak, J-J. Sinou, and S. Nacivet. A novel hybrid surrogate model and its application on a mechanical system subjected to friction-induced vibration. *Journal of Sound and Vibration*, 434:456–474, 2018.

[24] E. Denimal and J-J. Sinou. Advanced kriging-based surrogate modelling and sensitivity analysis for rotordynamics with uncertainties. *European Journal of Mechanics-A/Solids*, 90:104331, 2021.
[25] E. Denimal and J-J. Sinou. Advanced computational technique based on kriging and polynomial chaos expansion for structural stability of mechanical systems with uncertainties. *Journal of Engineering Mathematics*, 130(1):1–19, 2021.

[26] J-J. Sinou and E. Denimal. Reliable crack detection in a rotor system with uncertainties via advanced simulation models based on kriging and polynomial chaos expansion. *European Journal of Mechanics-A/Solids*, 92:104451, 2022.

[27] E. Denimal, L. Nechak, J-J. Sinou, and S Nacivet. Kriging surrogate models for predicting the complex eigenvalues of mechanical systems subjected to friction-induced vibration. *Shock and vibration*, 2016, 2016.

[28] D. Xiu and G. Karniadakis. Modeling uncertainty in flow simulations via generalized polynomial chaos. *Journal of Computational Physics*, 187(1):137–167, may 2003.

[29] B. Sudret. Global sensitivity analysis using polynomial chaos expansions. *Reliability engineering & system safety*, 93(7):964–979, 2008.

[30] G. Blatman and B. Sudret. Adaptive sparse polynomial chaos expansion based on least angle regression. *Journal of computational Physics*, 230(6):2345–2367, 2011.

[31] S. Lophaven, H. Nielsen, and J. Søndergaard. DACE: a Matlab kriging toolbox, volume 2. Citeseer, 2002.

[32] C. Rasmussen and C. Williams. Gaussian processes for machine learning the mit press. *Cambridge, MA*, 2006.

[33] J. Kleijnen. Kriging metamodeling in simulation: A review. *European journal of operational research*, 192(3):707–716, 2009.

[34] I. M. Sobol and S. S. Kucherenko. Global sensitivity indices for nonlinear mathematical models. Review. *Wilmott*, 2005(1):56–61, jan 2005.