David Lebel

Christian Soize

Christine Fünfschilling

Guillaume Perrin

Statistical inverse identification for nonlinear train dynamics using a surrogate model in a Bayesian framework

teaching and research institutions in France or abroad, or from public or private research centers.   Statistical inverse identification for nonlinear train dynamics using a surrogate model in a Bayesian framework

Introduction

The goal of this work is to perform the inverse identification of the input parameters of an expensive computer code with functional input excitation and output response, under uncertainty, from experimental data. The industrial application motivating the work presented in this paper consists in monitoring high-speed train suspensions using prediction for untested conditions. As long as the quality of the prediction is satisfying, the calibrated value of the parameters may not correspond to the "real" physical value, in the case when the parameters have a physical meaning. The problem studied in this paper differs in this point because the "real" physical value of the parameters is precisely what we are looking for.

Because the problem is affected by various sources of uncertainty (noise on the measurements, model error, lack of knowledge about the model parameters), a probabilistic model of the system response is built based on the simulation code. The Bayesian framework is well suited to combine a probabilistic model with experimental data to obtain information about the model input parameters. A Bayesian calibration approach is adopted. However, except a few particular analytical cases, Bayesian calibration methods generally require the computation of the model response for numerous values of the input parameters. Surrogate modeling is a solution to circumvent the numerical cost induced by the numerous calls to an expensive computer code. It consists in replacing the latter by an efficient algebraic approximation of the output response. Although it is a classical approach for simulations with scalar output, surrogate modeling of functional output remains a complex task. Because we want to keep the whole information provided by the measurements, we ruled out the possibility of condensing the output into a few indicators of interest. The solution we propose in this paper is, instead of relying on a surrogate model of the output response, building a surrogate model of the likelihood function that is at the core of Bayesian calibration. The novelty of this work also lies in the consideration of the uncertainty introduced by the surrogate model.

We start this paper by recalling in Section 2 the main features of MCMC-based Bayesian calibration, followed by two examples of calibration procedures. The approach we propose is then detailed in Section 3, with a particular focus put on including the surrogate model uncertainty in the calibration procedure. The method is then applied on our industrial railway case in Section 4.

Setting the calibration problem in a Bayesian framework

This section recalls the main principles of Bayesian calibration and introduces the formalism that is used in the rest of the paper. Two calibration methods for cases approaching the one treated in this paper are also presented.

Bayesian calibration principle

We consider a system described by the model Y = H(W). Variable W gathers the various input parameters of the model. In general, W can be defined as a vector belonging to an admissible set C W . In the context of Bayesian calibration, the value of these parameters is uncertain, hence vector W is random. The initial information about the parameters is given by the prior probability density function (PDF) p prior W of W. Variable Y is an observable output quantity of the system. The dependence of Y to parameters W is modeled by the function H. Function H is considered random (meaning that H(w 0 ) is a random value even if w 0 is deterministic), which makes Y random as well. Various sources of uncertainty may account for H randomness: model error, approximation errors... A measurement y mes is provided as experimental data that can be affected by measurement noise. It is considered as a realization of random quantity Y.

The goal of the Bayesian calibration procedure is to determine the posterior PDF p post W of parameters W. The latter represents the updated knowledge about the parameters, according to the new information about the system brought by measurement y mes .

Mathematically, p post

W corresponds to the conditional PDF of W knowing Y

p post W (w) = p W | Y (w | y mes ), w ∈ C W . (1)
According to the Bayes formula, the conditional probability can then be decomposed in the following way:

p post W (w) = p Y | W (y mes | w) p prior W (w) p Y (y mes ) (2) ∝ L(w) p prior W (w) (3)
where, for a fixed y mes , the function From Eq. ( 3), Markov Chain Monte Carlo (MCMC) is a classical way to construct independent realization of posterior PDF p post W . This class of algorithms is briefly presented in the following section. It should be noted that calibration methods that do not rely on MCMC have been developed, in order to avoid the numerous difficulties associated with the implementation of MCMC algorithms. Other approaches for constructing surrogate models consist in using polynomial chaos expansions [START_REF] Marzouk | A stochastic collocation approach to bayesian inference in inverse problems[END_REF][START_REF] Arnst | Identification of bayesian posteriors for coefficients of chaos expansions[END_REF][START_REF] Frangos | Surrogate and reduced-order modeling: A comparison of approaches for large-scale statistical inverse problems[END_REF][START_REF] Nagel | Spectral likelihood expansions for Bayesian inference[END_REF].

L : w → p Y | W (y mes | w) is

Estimation of the posterior density with MCMC

The purpose of Markov Chain Monte Carlo (MCMC, see [START_REF] Givens | Computational statistics[END_REF]) is to draw a sample according to a target PDF. It requires the ability to evaluate this PDF or a function proportional to it anywhere in its definition set. MCMC algorithms perform this draw by building a Markov Chain for which the invariant probability distribution is represented by the target PDF. For Bayesian calibration, the target PDF is posterior PDF p post W ; the function injected in the algorithm is w → L(w) p prior W (w). The drawn sample usually needs to be large in order to be distributed as the target PDF. At least one call to the likelihood function is required for each point of the sample.

Once the sample is drawn, the target PDF can be studied by estimating its moments and quantiles for instance. The marginal PDFs can also be plotted using histograms or kernel methods.

Calibration of an expensive computer code with scalar output

Let us first consider the case of a simulation-based model with scalar output, represented for example by the relationship

Y = h(W) + ε (4)
where the simulation is represented by the deterministic function h and ε that is an additive noise representing the measurement noise and the model error, modeled by a Gaussian centered random variable of variance σ 2 ε . For a single measurement y mes , the corresponding likelihood function is

L(w) = p N (y mes ; h(w), σ 2 ε ), (5)
where p N (.; µ, σ 2 ) stands for the normal density of mean µ and variance σ 2 .

Using MCMC requires numerous evaluations of the likelihood function L and consequently numerous runs of the simulation represented by function h. An expensive simulation makes the procedure unaffordable. A classical way of addressing such numerical issue is to rely on surrogate models. In particular, Gaussian process (GP) models is a commonly used class of surrogate models, as it offers closed-form expressions and an estimation of the approximation error. The principle of GP surrogate modeling is to represent a deterministic function by a conditioned Gaussian process. The mean function of the process (also called the Kriging predictor) constitutes the best approximation of the function anywhere in its definition set in a L 2 -sense, while its variance represents the approximation error. More details, formulas and their mathematical justifications are provided in Appendix A.

For the present calibration case, the simulation is replaced by a GP surrogate model:

for any w ∈ C w , the simulation output h(w) is replaced by the random value ĥ(w) + Z(w) where ĥ is the Kriging predictor and Z a centered GP of variance σ 2 Z independent from noise ε. The model then becomes

Y = ĥ(W) + Z(W) + ε (6)
and the associated likelihood function

L(w) = p N (y mes ; ĥ(w), σ 2 Z (w) + σ 2 ε ), (7)
As shown by this last equation, the interest of building a GP surrogate model of the simulation output is that the surrogate model uncertainty can be readily introduced in the likelihood function. The calibration procedure is not modified. The Bayesian formalism is kept, which allows for a correct evaluation of the posterior uncertainty on the parameters W.

Calibration of an expensive computer code with functional output

We now consider a simulation-based model with an output which is no longer scalar but functional. The method presented in the previous section can no longer be applied as such. Indeed, the surrogate modeling of functional quantities is a complex task and remains a subject of current research (see for example [START_REF] Perrin | Adaptive calibration of a computer code with time-series output[END_REF]).

For the calibration of expensive computer codes, the issue of functional outputs is addressed in [START_REF] Ranjan | Inverse Problem for a Time-Series Valued Computer Simulator via Scalarization[END_REF] and [START_REF] Pratola | Fast sequential computer model calibration of large nonstationary spatial-temporal processes[END_REF]. In both papers, the authors perform a scalarization of the problem. They define a distance between the experimental and modeled output, based on a L 2 -norm for [START_REF] Ranjan | Inverse Problem for a Time-Series Valued Computer Simulator via Scalarization[END_REF] or on likelihood ratios for [START_REF] Pratola | Fast sequential computer model calibration of large nonstationary spatial-temporal processes[END_REF]. The distance being scalar, it can be represented by a GP surrogate model. The calibration procedure then consists in an optimization problem. The GP model is used to minimize the defined distance. The goal is to obtain the optimal parameters values that provide the best fit between the experimental data and the model, according to the chosen distance.

Although both methods achieve interesting results, they cannot be considered as Bayesian approaches. They do not take into account the uncertainties of the model to determine an uncertainty on the calibrated parameters. In [START_REF] Pratola | Fast sequential computer model calibration of large nonstationary spatial-temporal processes[END_REF], the authors propose to estimate the distribution of the minimum of the GP model. This analysis allows for evaluating the uncertainty on the parameters stemming from the approximation by a surrogate model, but still ignores the other sources of uncertainty, such as model error or measurement noise.

Calibration with GP surrogate model of the likelihood function

The goal of this paper is the development of a Bayesian calibration method for expensive computer codes with functional output. As shown in [START_REF] Ranjan | Inverse Problem for a Time-Series Valued Computer Simulator via Scalarization[END_REF] and [START_REF] Pratola | Fast sequential computer model calibration of large nonstationary spatial-temporal processes[END_REF], scalarization combined with GP surrogate modeling is an efficient way of addressing high computational costs caused by a simulation-based model. In the Bayesian formalism, a natural scalarization is provided by the likelihood function. Hence we propose to perform the calibration using a GP surrogate model of the likelihood function. The objective is to remain in the Bayesian framework while taking advantage of GP surrogate modeling.

This approach raises various questions that will be addressed in this paper:

• how should a GP model of the likelihood function be built ?

• how to perform MCMC with a random likelihood function ?

• how to take into account the uncertainty associated with the approximation error in the estimation of the parameters posterior probability distribution ?

H.

In the following sections, we propose a method that takes into account the uncertainty of the GP surrogate model in order to estimate more correctly the calibration error. In Section 4.5, we refer to this new approach as the MCT (for "Monte Carlo on the trajectories") method.

Monte Carlo on the GP model trajectories

If the deterministic log-likelihood function L is replaced by the random surrogate model L(.; Θ), the corresponding posterior PDF becomes random as well. Random variable Θ accounts for the randomness of the surrogate model; it is independent of the other "physical" random quantities involved in the calibration problem, such as the system parameters W or the system response Y. The random posterior PDF corresponding to L(.; Θ) can be expressed using the conditioning by Θ as

w → p W | Y,Θ (w | y mes , Θ) . (8)
Yet, the quantity we are looking for remains unchanged: the conditional PDF of parameters W knowing Y. The rule of conditional probabilities states that it is equal to the expectation with respect to Θ of the previous PDF:

p post W (w) = p W | Y (w | y mes ) (9) = E Θ p W | Y,Θ (w | y mes , Θ) (10)
≈ 1 N N i=1 p W | Y,Θ (w | y mes , θ i ). (11)
Eq. ( 11) corresponds to the empirical estimate of the expectation using N realizations

{θ i } 1≤i≤N of random variable Θ. This Monte Carlo approach requires the estimation of N PDF p W | Y,Θ (w | y mes , θ i ) by MCMC using N deterministic trajectories L(.; θ i ) of surrogate model L(.; Θ).
In practice, for each i, the MCMC algorithm provides a sample, subset of

C W , distributed as p W | Y,Θ (w | y mes , θ i ).
If, for all i, these samples have the same number of points, their simple concatenation corresponds to a sample approximately distributed as p post W .

Approximation of the GP model trajectories

The Monte Carlo approach proposed in the previous section implies the ability to draw a trajectory L(.; θ) of the Gaussian surrogate model and to evaluate its value at any point w in admissible set C W . A classic way to manipulate a Gaussian process trajectory is to compute its value on a fine grid spanning the entire admissible set. However, this solution is not adapted to the calibration of more than two or three parameters, because the number of points on the grid increases exponentially with the dimension of the admissible set, equal to the number of parameters to calibrate. Moreover, since the MCMC algorithm randomly draws the points at which the likelihood needs to be evaluated, computing the values of a trajectory in advance is not relevant. Another approach could consist in iteratively conditioning the Gaussian process by the previously drawn points for each step of the MCMC. However, it implies the manipulation of a full covariance matrix whose size would increase at each step.

Instead of exactly computing the trajectory, we propose to approximate it in the following way: first, the value of trajectory L(.; θ) is computed at the points of a con-

ditioning set W c = {w c j } 1≤j≤Nc .
It should contain only a limited number N c of points. The approximation L(.; θ) of L(.; θ) then consists of the expectation of the surrogate model conditioned by the value of the trajectory at the points of W c :

L(.; θ) : w → E Θ L(w; Θ) | L(w c j ; Θ) = L(w c j ; θ), 1 ≤ j ≤ N c . (12)
The deterministic function L(.; θ) can then be injected in the MCMC algorithm to

estimate PDF p W | Y,Θ (. | y mes , θ).
The crucial step of the method lies in the choice of conditioning set W c . The aim of this conditioning approach is to decrease as much as possible the variance of the surrogate model. Indeed, the lower the variance of a stochastic process, the closer any trajectory is to the mean function. We chose to focus on the region of interest of C W , where the relative value of the likelihood function is high. We define this region as the set P where the surrogate model has a probability higher than a tolerance ρ ∈]0, 1[ to be greater than at the point where its mean function is maximum:

P = {w ∈ C W | P (L(w; Θ) > L(w max ; Θ)) ≥ ρ} (13)
with w max = arg max

w∈C W E Θ {L(w; Θ)}.
Conditioning set W c is then built as a discrete set of N c points space-filling in P. First a large sample uniformly distributed in P is drawn using MCMC with the indicator function of P. Then several subsets of size N c are randomly drawn in this sample. The "most space-filling" subset according to the maximin criterion defined in Appendix A.6 is kept as W c . However, relying solely on measurements of the train dynamic response is not sufficient, because of its strong dependence on the track geometric irregularities. The latter constitute the main excitation source of a rolling train and, consequently, have a major influence on its dynamic behavior (see [START_REF] Perrin | Quantification of the influence of the track geometry variability on the train dynamics[END_REF][START_REF] Lestoille | Sensitivity of train stochastic dynamics to long-time evolution of track irregularities[END_REF][START_REF] Lestoille | Stochastic prediction of high-speed train dynamics to long-term evolution of track irregularities[END_REF]). Track geometry is also subject to damage caused by railway traffic (see [START_REF] Bing | Development of railroad track degradation models[END_REF]): it is gradually degraded and regularly maintained. Consequently, the track geometric irregularities evolve through time. Because of the high sensitivity of the system to these irregularities, their evolution has to be taken into account to perform a correct monitoring of the suspensions state. Train dynamics simulation is thus necessary to include the excitation source in the analysis. More precisely, we propose to build a simulation-based model of the train dynamic response and to calibrate its parameters describing the mechanical properties of the suspensions of interest. The experimental data used for this calibration consist of joint measurements of the track geometric irregularities and of the train dynamic response.

Probabilistic model of the high-speed train dynamic response

This section presents how the train response model is built, based on simulation.

This construction is sketched in figure 1. It has been shown in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF] that the track geometric irregularities can be modeled as a nonstationary R 4 -valued random field {X(s), s ∈ [0 , S]} indexed by curvilinear abscissa [0 , S]. Though the irregularities are represented by random field X, the calibration procedure presented here does not make use of the model developed in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF].

X Simulation h(X; W) Y W B
Instead, the available irregularities measurements performed on various track stretches are considered as realizations of random field X.

The measurement method for the track geometric irregularities has shown a very good reproducibility. Consequently, the measurement noise affecting these irregularities is considered negligible compared to the other sources of uncertainty.

Description of the high-speed train model

The train multibody model consists of rigid bodies linked together by mechanical joints (stiffnesses and dampers) with nonlinear behavior. The wheel-rail contact law is also nonlinear. The flexible modes of the different bodies are not taken into account.

More details about the definition of the multibody train model can be found in [START_REF] Kraft | Parameter identification for a TGV model[END_REF].

Such multibody model contains numerous parameters : body masses and inertiae, mechanical properties of the joints and relative positions of the bodies and the joints.

For this application, only q = 7 mechanical parameters are involved in the calibration process. These parameters of interest have been chosen based on the criticity of the corresponding suspension element and on their influence on the train response. Moreover, for each bogie, all the suspensions of the same type depend on one single parameter.

Indeed, we observed that we were unable to distinguish between the contributions of the different elements of the same type with the available sensors. It also allows for significantly decreasing the dimension of the problem.

The q parameters are gathered in the random vector W, belonging to the admis-

sible set C W ⊂ R q . The
It associates the response â = h(x; w) in the frequency domain with the irregularity signal x and vector w of parameters. If the quantities x and w are replaced by their stochastic counterpart X and W, the simulation output h(X; W) becomes stochastic as well.

The multibody modeling and simulation necessarily contain inaccuracies and simplifications compared to the real system. Numerical solving is also a source of errors.

To perform a robust calibration, a train model error has to be introduced, in order to take into account the fact that the model cannot exactly represent the reality. Moreover, the measurements of the train response performed by embedded accelerometers may contain a certain level of measurement noise. These two types of uncertainties (model Frequency ω (Hz) Amplitude (dB)

Y = h(X; W) + B . (16)
4.3. Specificities of the calibration procedure for the applicative railway case

Experimental data

The calibration is performed to determine the suspensions health state at a given date. We are focusing on a single train, for which a model error has been identified. At a given date, the experimental data consist of a set U of independent joint measurements of the track geometric irregularities and of the corresponding train dynamic response, denoted respectively x mes,i and y mes,i , 1 ≤ i ≤ ν on ν different track stretches:

U = {(x mes,i , y mes,i )} 1≤i≤ν . (17)

Calculation of the likelihood function

The calculation of the likelihood function presented in Section 2 is adapted to the present railway case. We need to take into account measurements of the track geometric irregularities. Consequently, the model parameters are not conditioned solely by the output Y but also by the input-output couple (X, Y). It makes a difference to consider the probability of obtaining a certain dynamic response whatever the irregularities, and to obtain the same response knowing the irregularities that have triggered it.

Moreover, since several independent measurements are used to perform the calibration, it is necessary to introduce a set {(X i , Y i )} 1≤i≤ν of ν independent and identically distributed copies of couple (X, Y). Eq. ( 1) giving the expression of the posterior PDF of parameters W then becomes (the PDF arguments are omitted for simplicity):

p post W = p W | X 1 ,Y 1 ...X ν ,Y ν (18) ∝ p Y 1 ...Y ν | X 1 ...X ν ,W p W | X 1 ...X ν (19) ∝ ν i=1 p Y i |X i ,W p prior W . (20)
Eq. ( 20) is obtained by considering that:

• the track geometric irregularities X and the train parameters W are independent, so p W | X = p prior W ;

• for any (i, j) ∈ {1, . . . ν} 2 , X i and Y i are supposed to be independent from X j and Y j if i = j. Indeed, the track has been divided into numerous stretches so that the realizations of the track geometric irregularities and of the train dynamic response can be considered independent on two different stretches.

The likelihood function is then given by

L : w → ν i=1 p Y i |X i ,W (y mes,i |x mes,i , w) . ( 21
)
Using Eq. ( 16) the log-likelihood function can be written

L : w → ν i=1 log p B (y mes,i -h(x mes,i ; w)) . ( 22
)
From the previous expression, one can observe that ν simulation runs are required for each evaluation of log-likelihood L at each points w ∈ C W . The number ν of available measurements can typically be several hundreds in the present case, which explains the high computational cost for each call to the likelihood function.

Parameters of the GP surrogate modeling

For this application, the number of suspension parameters to identify is q = 7 (see Section 4.2.2). The number of parameters gives the dimension of random vector W and the dimension of the definition set of likelihood function L.

The size of training set is 500. The chosen form for the covariance function is Matérn-5 2 [START_REF] Matérn | Spatial variation[END_REF]. Moreover, we chose to build the GP surrogate model considering that the observations of the log-likelihood function were noisy (see Appendix A.5). The variance of this noise was optimized to fit the available data along with the other correlation parameters. This choice was made in order to improve the quality of the surrogate model. Introducing a small noise improved its regularity. Indeed, it offers the surrogate a certain margin of freedom around the observation. On the contrary, forcing the surrogate model to be strictly interpolating by considering that the observations are exact sometimes resulted in unexpected oscillations of the surrogate model.

We tried to refine the Gaussian surrogate around the maximum of the likelihood function following the KGCP policy proposed in [START_REF] Scott | The Correlated Knowledge Gradient for Simulation Optimization of Continuous Parameters using Gaussian Process Regression[END_REF]. However, no significant improvement of the calibration results were observed with the refined surrogate model. It suggests that the initial training set is large enough to correctly know the behavior of the likelihood function. The observations being considered noisy, adding new points to the training set does not significantly increase the accuracy of the maximum location.

Choice of the MCMC algorithm

In our preliminary tests, the posterior PDF appeared to be very peaked, which prevents the classical Metropolis-Hastings algorithm [START_REF] Hastings | Monte Carlo Sampling Methods Using Markov Chains and Their Applications[END_REF] from correctly sampling from it. We also tried to run several chains in parallel, but the final results were too dependent on the initial choice of the chains starting points. Consequently, we chose to use a different algorithm called Transitional MCMC [START_REF] Ching | Transitional Markov Chain Monte Carlo Method for Bayesian Model Updating, Model Class Selection, and Model Averaging[END_REF][START_REF] Betz | Transitional Markov Chain Monte Carlo: Observations and Improvements[END_REF], more adapted to peaked target PDF. Details about the algorithm are given in Appendix C.

Calibration results

This section presents the results obtained with the Bayesian calibration method for different cases:

• numerical experiments, that is to say using simulated data;

• with measurements performed at the reference date, which were used to identify the model error;

• with measurements performed six month after the reference date.

The goals of this section are to:

• validate the proposed method, using the numerical experiments and the measurements at the reference date;

• show the impact of considering the uncertainty introduced by the use of a surrogate model of the likelihood function on the calibration results.

It should be noted that the speed up provided by the use of a surrogate model compared to an estimation of the parameters posterior distribution that would rely on the computation of the exact likelihood is estimated of order 10 6 .

Numerical experiments

In the present case, a numerical experiment consists of simulated train responses that are used as if they were experimental data. They are generated using actual measurements of the track geometric irregularities on several track stretches and known degraded suspension parameters. Moreover, an independent realization of the error B is added to the response signal on each track stretch in order to generate a quantity as close as possible to an actual measurement. The numerical experiments allows for validating the calibration procedure: the procedure is applied on the artificial train response, the calibration results can then be compared to the reference parameters used to generate the response.

At a given date, we suppose that a set of ν 1 track irregularities measurements {x mes,i } 1≤i≤ν1 is available. The validation procedure from a numerical experiment can then be summed up as follows.

1. Choose artificial parameters w 1 ;

2. Run the simulation on the ν 1 track stretches with parameters w 1 ; 3. Generate ν 1 independent realizations {b i } 1≤i≤ν1 of output predictive error B.

Because B is Gaussian, this generation can be performed very easily from normal independent samples using Eq. (B.6) ; 4. Add these realizations to the simulated response to obtain artificial realizations of the train dynamic response

y num,i = h(x mes,i ; w 1 ) + b i , 1 ≤ i ≤ ν 1 ; (23) 5.
Perform the calibration using input data

U 1 = {(x mes,i , y num,i )} 1≤i≤ν1 (24)
to obtain the calibrated (random) parameters W opt 1 ; 6. Compare W opt 1 to w 1 .

The validation procedure has been performed with different values of the artificial parameters w 1 . Figure 6 presents calibration results with w 1 set to the nominal values.

On this graph are displayed the marginal densities of the posterior PDF of W opt 1 . It also compares the results obtained with the KP method that solely uses the Kriging predictor provided by the surrogate model of the likelihood function and with the MCT method that also includes the uncertainty related to the surrogate model.

One can first observe that the distributions are close to the nominal value (corresponding to 0.5 on the graph). Except parameter 2, the difference between the maximum of the marginal PDFs and the nominal value is always lower than 5% of the size of the admissible interval. The dispersion varies from one parameter to another. This can be explained by the initial choice of the admissible intervals, since the results are rescaled according to their size, but more importantly by the sensitivity of the train dynamic response to the different mechanical parameters. The results for parameter 2 are a good illustration of this sensitivity question. Contrary to the other parameter, it corresponds to a suspension that is not located on the same bogie as the sensors, but at the other end of the carbody. Consequently, we expect this parameter to have less influence on the train dynamic response measured by the sensors. This is coherent with the results: the distribution is further away from the nominal value, and its dispersion is greater than for the other parameters. These considerations also highlight the interest of Bayesian calibration: because we measure the uncertainty on the calibrated parameter, we have a way to assess the accuracy of the calibration and the confidence we can put on its results. The fact that the marginal PDFs are rather peaked comes from the large database that is available for the calibration. The comparison of the results obtained using the KP and MCT methods shows that the marginals PDF are more spread with the MCT methods. With the KP method, a source of uncertainty, the error introduced by the approximation of the likelihood function by a surrogate model, is not taken into account. As a consequence, the uncertainty on the calibrated parameters is reduced. Using the KP method thus leads to an overestimation of the calibration accuracy. Nevertheless, the marginals maxima appear to be located at very similar parameter values with the two methods. The KP method seems satisfying to determine only the most probable parameter values. It also allows for estimating the proportion of the calibrated parameters variance that stems from the system uncertainties. The graph layout is the identical to figure 7.

Results using actual measurements of the high-speed train dynamic response

This section presents calibration results obtained with actual measurements. The first set of measurements is the one collected at the reference date, which is used to identify the output predictive error B, denoted as T 0 . The second set gathers measurements performed at a date T 1 , six months after the reference date T 0 .

Figure 8 presents the calibration results obtained for these two dates using the MCT method. As expected, one can observe that the parameter mean values for T 0 are close to the nominal value (0.5 on the normalized scale). Indeed, the output predictive error has been identified considering the parameters value is nominal at the reference date T 0 . This good correspondence is a second way of validating of the method. The results for the date T 1 show a significant evolution compared to the nominal values. Except for the second parameter, whose case has been treated in the previous section, we can observe a high confidence in the calibrated results. To obtain a complete validation of the identification method, experimental tests need to be performed on the isolated suspension elements in order to measured their mechanical characteristics. Such tests have not been conducted yet. to the ones made for figure 6. A difference that can be noted between the numerical experiments and the calibration on the real measurements is the fact that the posterior PDF is more dispersed for the latter. This higher dispersion highlights the limitation of the chosen model for the output predictive error (a Gaussian additive error being independent of the parameter space). For the numerical experiments, by construction, this model exactly matches the error artificially applied on the generated data, which appears to be less true for the measured data. The use of a more complex, non-Gaussian model for the output predictive error could be considered in future works.

Conclusion and perspectives

We have proposed a method for the Bayesian calibration of expensive computer codes. The main particularity consists in the surrogate modeling of the likelihood function to address the computational costs. We then focused on evaluating the consequences of this use of an approximation on the calibration results. The uncertainty associated with the surrogate modeling is taken into account by estimating the parameters posterior PDF from trajectories of the random surrogate model. We thus make full use of the variance of the GP model that evaluates the accuracy of the approximation.

The need to draw GP trajectories also made us develop an approximation method for the latter, relying on a further conditioning of the GP model.

The procedure has been tested on an industrial case, the monitoring of high-speed train suspensions using acceleration measurements. In this application, the complete MCT method appears to have little effect on the mean value of the calibrated parameters compared to the simple use of the KP method. It however significantly increases the variance of the calibrated parameters. This emphasizes the interest of the method if importance is put on correctly estimating the confidence of the calibration.

Function y is supposed to be a particular sample path of an underlying Gaussian second-order real-valued stochastic process {Y (x), x ∈ X }, indexed by X . The surrogate model consists of this GP conditioned by the observations of y at the points of the training set. The underlying process is defined according to a parametric formulation, and denoted

Y | {β, σ, ψ} ∼ GP f (.) t β ; σ 2 R(., . | ψ) (A.1)
where:

• GP(m; Σ) denotes the Gaussian process whose mean function is

x → m(x) = E{Y (x)} from X into R and covariance function is (x, x ) → Σ(x, x ) = E{(Y (x) -m(x))(Y (x ) -m(x ))} from X 2 into R;
• f : X → R p gathers p deterministic regression functions on X ;

• β ∈ R p is the vector of regression coefficients parameterizing the GP mean function;

• R(x, x | ψ) gives the shape of the covariance function. It depends only on x-x and on a vector ψ of parameters (for example, correlation lengths). It is assumed that R(x, x | ψ) = 1. Consequently, σ 2 is the variance of random variable Y (x) for x fixed in X and is thus independent of x.

The mean function of the process is decomposed as a linear combination of regression functions chosen by the user, polynomials for example. The first parameter is thus the vector of regression coefficients β. The second parameter is the GP variance. The Stationarizing consists in averaging on the diagonals of the matrix in order to obtain the where CV j (q) is the coefficient of variation (equal to the standard deviation divided by the mean) of set {L(x j,k ) q-qj } 1≤k≤Ns . If q j > 1, set it to 1.

y n =      y(x 1 ) . . . y(x n )      , Y n =      Y (x 1 ) . . . Y (x n )      . (A.2) Y n | {β, σ,
  Y (x) Y n   | {β, σ, ψ} ∼ N     f (x) t [F ]   β ; σ 2   1 r(x) t r(x) [R]     (A.3) with [F ] =      f (x 1 ) t . . . f (x n ) t      , r(x) =      R(x, x 1 | ψ) . . . R(x, x n | ψ)      (A.4) and [R] ij = R(x i , x j | ψ), 1 ≤ i, j ≤ n. Matrix [R] is assumed to be invertible.

  called the likelihood function. Its computation depends on the relationship Y = H(W). Because of dimensionality problems, the log-likelihood L : w → log(L(w)) is usually used instead.

Figure 1 :

 1 Figure 1: Diagram of the system quantities of interest

Figure 2 :

 2 Figure 2: Track irregularities example. From bottom to top: lateral, gauge, vertical, cross-level irregularities.

Figure 3 :

 3 Figure 3: Vertical (top curve) and lateral (bottom curve) accelerations in a carbody, in the time domain.

Figure 4 :

 4 Figure 4: Vertical (top curve) and lateral (bottom curve) accelerations in a bogie, in the time domain.

Figure 5 :

 5 Figure 5: Lateral acceleration in a carbody, represented in the frequency domain. The solid line represents the mean function, while the filled area represents the 5% and 95% quantiles.

Figure 6 :

 6 Figure 6: Results of the numerical experiment with nominal parameters. The parameter values are alsonormalized so that every parameter varies between 0 and 1; 0.5 then corresponds to the nominal value. The curves represent the normalized posterior PDF marginals for each parameter, obtained using the KP (red or dark gray) and MCT (orange or light gray) methods.

Figure 7 W opt 1 show a good correspondence to the input parameters w 1 .

 711 Figure 7 presents calibration for arbitrary degraded values of the artificial parameters w 1 . Two different locations of the bogie are considered: at the head and the rear of the train. The vertical and lateral accelerations the train undergoes at these two locations are indeed very different. One can also note on the bottom graph that we studied the case of a parameter set at the boundary of the admissible interval. The calibrated parameters W opt 1 show a good correspondence to the input parameters w 1 . The numerical experiments results are conclusive for the validation of both the KP and the MCT methods, with the limits previously noted for the KP method.

Figure 7 : 1 (

 71 Figure 7: Results of the numerical experiments with arbitrary parameters, using the MCT method, for a bogie located at the head (top graph) and rear (bottom graph) of the train. The input parameters values w 1 (black triangles) are compared to the mean of the posterior PDF marginals (orange dots). The orange lines represents the 98% confidence intervals around these calibrated values. The parameter scale normalization is identical to figure 6.

Figure 8 :

 8 Figure 8: Calibration results using the MCT method at time T 0 (blue diamonds) and T 1 (orange circles).

Figure 9

 9 Figure 9 compares the calibration results using the KP and the MCT method for date T 1 . The observations concerning the comparison of the two methods are similar

Figure 9 :

 9 Figure 9: Comparison of calibration results obtained at time T 1 using the KP (red or dark gray curves) and the MCT (orange or light gray curves) methods. The graph layout is identical to figure 6.

  shape of the correlation function R(., .) (exponential or Matérn for instance) is also chosen by the user according to the expected regularity of y. The third parameter is the vector of correlation parameters ψ. They depend on the shape of the correlation function. Parameters (β, σ, ψ) are a priori unknown. They are not set by the user but have to be determined using the information about y provided by the observations on the training set. Appendix A.2. Conditioning by the observations Let y n and Y n be respectively the values of function y and of process Y at points (x i ) 1≤i≤n :

  Figure A.10: 1-D example of GP surrogate model with exact observations: the solid black line represents target function y, the red diamonds the points of the training set, the green dashed line the mean function of the GP model and the green area its standard deviation.

2 . 1 ( 3 .

 213 For 1 ≤ k ≤ N s , compute the weighting coefficient w j,k = L(x j-1,k ) qj -qj-Compute the covariance matrix for the proposal distribution[Σ j ] = β 2 Ns k=1 wj,k (x j-1,k -m j )(x j-1,k -m j ) t (C.6) with m j = Ns k=1 wj,k x j-1,k . 4. Initialize {x c k } 1≤k≤Ns such that x c k = x j-1,k for any k. 5. Build the sample {x j,k } 1≤k≤Ns by iterating as following : for k going from 1 to N s , (a) Select index between 1 and N s , with a probability given by the weight wj, ; (b) Draw a proposal x * from a normal distribution N (x c ; [Σ j ]); (c) Draw r from a uniform distribution on [0, 1]; (d) If r ≤ pj (x * ) pj (x c ) then set x c = x * , otherwise do nothing; (e) Set x j,k = x c . 6. Increase j.

  ψ} is a marginal of process Y | {β, σ, ψ}. As a consequence, for any x ∈ X , the vector concatenating Y (x) | {β, σ, ψ} and Y n | β, σ, ψ is Gaussian:

acceleration measurements for maintenance purposes.This identification problem can be seen as a calibration problem. The main purpose of calibration is to fit a model to experimental data, in order to provide a good Preprint submitted to Journal of Sound and Vibration July 5, 2019

mized along with parameters σ and ψ.

Acknowledgments

This research has been supported by SNCF, the French National Railway Company.

Appendix A. Gaussian process surrogate model

Gaussian process surrogate modeling consists in representing a deterministic scalar target function by a conditioned Gaussian process (GP). Information about the target function is provided by a set of observations of the function value at a few points in the definition set. These points constitute the initial training set. Randomness is introduced as a way to quantify the approximation stemming from the use of a surrogate model.

The objectives of GP surrogate modeling are for example to build an easy-to-compute approximation of an expensive computer code or to model a phenomenon for which all observations are affected by a random noise.

The evaluation of the target function at the points of the training set is generally the computationally expensive step of GP surrogate modeling. Consequently, this set should contain only a limited number of points while maximizing the information provided. Therefore a common choice is to define a space-filling training set to get as much information as possible on the global behavior of the target function on its whole definition set. Appendix A.6 provides an example on how to build a space-filling training set.

After this initial training phase, a second refining step (that will not be detailed in this paper) can be performed in order to explore particular features of the target function. For instance, [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF] proposes the EGO algorithm for the optimization of the target function using an expected improvement criterion. In [START_REF] Picheny | Adaptive Designs of Experiments for Accurate Approximation of a Target Region[END_REF], an adaptive design of experiment is proposed for a target function that needs to be accurately approximated around a certain level.

In this appendix, we present the principle of GP surrogate modeling based on the Bayesian approach presented in [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. Another approach giving equivalent results, based on the minimization of the mean square error, can be found in [START_REF] Sacks | Design and Experiments of Computer Experiments[END_REF].

Appendix A.1. Problem statement

The goal is to build a surrogate model of a given deterministic real-valued function y defined on a subset X of R q , whose value is known only at the n training points (x i ) 1≤i≤n in X . An illustration is provided on figure A.10.

Appendix A.3. Mean function parameter β

In this section, σ and ψ are supposed to be known. Only the regression coefficients β are supposed to be unknown. The fact that no information is a priori available about β is taken into account by following a hierarchical approach. Parameter β is represented by a random vector with non-informative prior: p β ∝ 1. The training 545 29 set is then used to learn about the distribution of β. The principle is to determine the distribution of β | Y n = y n for a fixed value of σ and ψ. In this section, the conditioning on σ, ψ is not systematically repeated for simplicity. Using the Bayes formula, one can write:

. This can be written as

. Variance and correlation parameters 555

Parameters σ and ψ could be estimated using the same hierarchical approach, by putting prior distributions on these variables. However, in this case, no closed form can be determined in general for Y . Instead, parameters σ and ψ are determined according to a criterion assessing how well the Gaussian process is fitting the data provided by the training set. The criterion usually used is the Maximum Likelihood Estimation (MLE).

The principle of MLE criterion is to maximize the density

The value of σ maximizing this density can be determined explicitly:

With this value of σ, the previous maximization is equivalent to minimizing the quan-

1 n in order to determine the optimal value of ψ. In general, no closed form exists for the covariance parameter; this step has to be performed numerically.

Appendix A.6. Space-filling training set

This appendix presents one method for the generation of a space filling training set in a set X consisting of the set product of bounded intervals of R. The chosen method is the optimization of a Latin Hypercube Sample (LHS) according to a maximin distance (see [START_REF] Chen | A review of design and modeling in computer experiments[END_REF]).

To build a N -points LHS, admissible set X must be divided into N cells along each dimension. The points are then scattered in X so that in every dimension, each cell only contains one point of the sample. The interest of a LHS is that if you consider a particular dimension, the sample points are regularly spaced, with no redundancy. 580 However, a LHS is not necessarily space-filling. We choose to measure the spacefilling property with the distance δ that consists of the smaller Euclidian distance between two points of the training set: for a finite discrete subset W of X , distance δ is defined as

The greater δ(W), the most space-filling W is. In practice, numerous LHS candidates are drawn. The one for which the value of criterion δ is the greatest is kept as the best training set in X .

Appendix B. Identification of the output predictive error

In this section, we consider that all the functional quantities are discretized. The µ discretization points (ω i ) 1≤i≤µ span the frequency band Ω. For the processes containing several components, the latter are concatenated to obtain a single column vector.

For instance, the discretized version of process {Y(ω), ω ∈ Ω} is the vector of size nµ

Reference set U 0 contains independent joint measurements of the track geometric irregularities and of the train response, denoted respectively x ref,i and y ref,i , 1 ≤ i ≤ ν 0 , performed on ν 0 track stretches:

For 1 ≤ i ≤ ν 0 , let then b ref,i be the difference between the measured train response and the simulated one (with nominal parameters) on the i th track stretch: Since random vector B is Gaussian, it is completely defined by its mean vector and its covariance matrix.

From the realizations {b ref,i } 1≤i≤ν0 the empirical estimates of the mean function m B and covariance function [ C B ] are computed: to each component or correlation between two components, in the following way: der to consider only the most statistically significant eigenvectors. Indeed, we observed that as the value of the eigenvalues decreases, the eigenvectors display characteristics similar to white noise. We choose m ≤ ν 0 , which has the advantage to make the estimated covariance matrix of the projection coefficients invertible (see Eq. (B.12)).

Coefficient α is determined by minimizing the Leave-One-Out error of projection (see [START_REF] Perrin | A Posteriori error and optimal reduced basis for stochastic processes defined by a finite set of realizations[END_REF]) of the realizations {b ref,i } i on the basis [V ] obtained by the spectral decomposition of [ C B ], which depends on α.

Random vector B can then be represented as

where γ is the centered random vector gathering the projection coefficients of B. For standard PCA, the components of γ are uncorrelated, but it is not the case here because of the stationarization. Thus it becomes necessary to estimate the covariance matrix of γ

where

Finally, we get

Appendix C. Transitional Markov Chain Monte Carlo

The TMCMC algorithm [START_REF] Ching | Transitional Markov Chain Monte Carlo Method for Bayesian Model Updating, Model Class Selection, and Model Averaging[END_REF][START_REF] Betz | Transitional Markov Chain Monte Carlo: Observations and Improvements[END_REF] is a MCMC algorithm designed to sample from the posterior PDF knowing the prior PDF and the likelihood function in a Bayesian framework. It is based on the following equation giving the relationship between these quantities:

The particularity of TMCMC is to start from a sample distributed according the prior PDF, and making it gradually evolve toward a sample distributed according to the posterior PDF. The algorithms work in m steps by sampling successively from the distributions p j defined for 0 ≤ j ≤ m and x ∈ X by

where 0 = q 0 < q 1 < • • • < q m-1 < q m = 1. One can immediately notice that p 0 = p prior and p m = p post .

At the beginning of each step j, 1 ≤ j ≤ m, the sample of size N s is distributed as p j-1 . Each point x j-1,k of the sample is then affected a weighting coefficient w j,k = L(x j-1,k ) qj -qj-1 , 1 ≤ k ≤ N s . A new sample of size N s is drawn from {x j-1,k } k according to these weighting coefficients. This new sample is distributed as p j . In order to avoid the repetition of identical elements in the new sample, MCMC steps are applied to disturb the sample while keeping the same distribution. The Metropolis-Hastings algorithm is used to draw the proposals for the MCMC steps: a Gaussian distribution around the previous point of the Markov Chain. Its covariance matrix is estimated from the sample {x j-1,k } k . A factor β is introduced to control the step size.

The algorithm can be summed up as follows.

Initialization: set j = 1, draw a sample {x 0,k } 1≤k≤Ns of N s distributed as p prior .

Iterations: while q j-1 < 1, 1. Determine the optimal value for q j by solving q j = arg min q∈]qj-1,1]

|CV j-1 (q) -1| (C.3)


Evangéline Capiez-Lernout
email: evangeline.capiez-lernout@u-pem.fr

Christian Soize
email: christian.soize@u-pem.fr

Nonparametric probabilistic approach for uncertainty quantification of geometrically nonlinear mistuned bladed-disks

Nonparametric probabilistic approach for uncertainty quantification of geometrically nonlinear mistuned bladed-disks

Evangéline Capiez-Lernout, Christian Soize

Introduction

Nowadays, a main challenge in engineering structural dynamics concerns the development of computational strategies for constructing stochastic nonlinear reduced-order models that are able to accurately reproduce dynamical structural responses that occur in nonlinear vibrational operating ranges. On one hand, a particular attention has to be paid on the construction of nonlinear deterministic reduced-order models [1,2,3,4,5], that have to be compatible with non intrusive computational strategies regarding existing industrial softwares. On the other hand, uncertainties can be implemented from the mechanical or geometrical parameters of a given deterministic computational model or with nonparametric probabilistic approaches, which are parameterized by hyperparameters directly linked to the considered reduced-order models. We are interested in this latter strategy that allows both model and parameter uncertainties to be considered. In previous researches, uncertainty was introduced from the nonlinear reduced-operators issued from a given mean reduced-order model [6,2,7]. With the proposed method, the uncertainty is introduced by replacing the deterministic reduced-order basis with a stochastic reduced-order basis. This latter one is obtained by using a new nonparametric probabilistic approach [8,9] so that each realization of the random projection basis respects some mathematical properties linked to the available information. Moreover, it is parameterized by a small number of hyperparameters. It is then used as the projection basis for constructing the stochastic nonlinear reduced-order model. With such uncertainty modeling, the numerical effort is focused on the construction of the stochastic nonlinear reduced internal forces, which is explicitly carried out using the stochastic reduced-order basis combined with the finite element method. The paper is organized as follows. In Section 2, the computational methodology is described. In particular, the computation of the mean nonlinear dynamical finite element response is used for constructing an adapted reduced-order basis with the proper orthogonal decomposition method. The stochastic nonlinear computational model is then written and the main steps concerning the computational resolution of the set of stochastic nonlinear differential equations are discussed. Section 3 is devoted to a numerical application consisting in a rotating bladed-disk structure for which the uncertainty propagation of the nonlinear dynamical response is analyzed in detail.

Description of the methodology 2.1. Context of the method

This research concerns the propagation of uncertainties through structural elasto-dynamical systems with finite displacements. In such a context, a complete methodology has already been developed and a dedicated numerical tool that is adapted for large numerical computational models and that is non intrusive with respect to the industrial softwares has been constructed [2,10,7]. More specifically, the mistuning analysis of bladed-disks that undergo large strains/displacements has been considered in [7]. In this case, the model used for taking into account mistuning uncertainties has been issued from the nonparametric probabilistic approach for which a dedicated probability model issued from the maximum entropy principle has been directly implemented on the reduced operators of a mean nonlinear reduced-order model [START_REF] Soize | Stochastic Models of Uncertainties in Computational Mechanics[END_REF]. With the presence of geometric nonlinearities, this strategy requires (i) the knowledge of all the linear, quadratic and cubic stiffness operators of a deterministic nonlinear reduced-order model that can be obtained for instance by judiciously exploring the indirect use of industrial softwares to get all these linear, quadratic and cubic stiffness contributions [3] or in the present case by explicit construction with the finite element method [2]. (ii) the construction of a positive stiffness operator whose entries are all the linear, quadratic and cubic stiffness contributions. (iii) the construction of the corresponding random stiffness operator whose thanks to the nonparametric probabilistic theory defined in [START_REF] Soize | Stochastic Models of Uncertainties in Computational Mechanics[END_REF] combined with a strategy for reducing the size of the random germ [12]. (iv) the extraction of the different random entries in order to compute the stochastic reduced nonlinear forces and the corresponding stochastic tangent operator. The first one is required by the computational algorithms that are used for solving the stochastic nonlinear differential equations. It should be noted that the use of the latter one depends on the nonlinear algorithm used (Newton or arc-length based).

In the present context, a nonparametric probabilistic approach on the model-form recently introduced in [8] is considered. This paper presents a first attempt for evaluating this new method in the context of turbomachinery mistuning. Such strategy requires (i) the knowledge of a deterministic reduced-order basis that is used for constructing a deterministic nonlinear reduced-order model for which the reduced quadratic and cubic stiffness do not need to be constructed.

(ii) the construction of a stochastic reduced-order basis by using this new nonparametric approach of model-form uncertainties. (iii) the explicit construction of the stochastic reduced nonlinear internal forces and its corresponding stochastic tangential operator by using the finite element method. Note that such construction is carried out for each finite element and that the assemblage is performed with a summation over each finite element.

As a consequence, the computational strategy has strongly been modified with respect to the previous researches and has therefore to be adapted. Below, the main steps of this new methodology are summarized, the related novel developments are explained and the advantages and drawbacks are pointed out.

Mean nonlinear finite element model of the structure

The structure under consideration is a bladed-disk structure with a M -order cyclic symmetry that rotates around its symmetry axis with constant angular speed Ω. The dynamical equations are expressed in the rotating frame of an equilibrium configuration. The structure is assumed to be fixed on a part of its boundary (i. e.the disk's inner radius). In the context of the finite element method, the nonlinear finite element computational model that describes the nonlinear dynamical forced response of the considered structure is characterized in the time domain by the following set of nonlinear coupled differential equations such that

[M ] Ü(t) + ([D] + [C G (Ω)]) U(t) + [K(Ω)]U(t) + F NL (U(t)) = F(t) , (1) [B] T U(t) = 0 (2)
in which the R n -vector U(t) is the instantaneous displacement vector. In Eq. ( 2), the real matrices [M ], [D] and [K(Ω)] are the mass, damping, and linear stiffness (n × n) real matrices with symmetry positive definiteness property and the real matrix , [C G (Ω)] is the gyroscopic coupling matrix with antisymmetry property. The R n -vector F(t) is issued from the finite element discretization of the external load and the R n -vector F NL (U(t)) describes the nonlinear finite element internal forces induced by the geometrical nonlinearities. The linear stiffness matrix can be written as

[K(Ω)] = [K e ] + [K g ] + [K c (Ω)] , (3)
in which [K e ] represents the elastic stiffness, where [K g ] describes the geometric stiffness that is calculated from the stress state of the system submitted to centrifugal loads and where [K c (Ω)] is the centrifugal stiffness. In Eq. ( 2), the real (n × n BC ) matrix [B] describes the n BC constraint relations defining the Dirichlet conditions, verifying the relation

[B] T [B] = [I n BC ].
Since we are interested in analyzing the nonlinear vibrations of the structure, the external load is defined in the time domain for t ∈ R and a Fourier transform with respect to the time domain is performed on the nonlinear solution and allows then the nonlinear dynamical response in the frequency domain to be analyzed. In the present research, it is chosen to solve this set of nonlinear differential equations in order to construct the nonlinear dynamical response that will be considered as the reference response. In this case, the nonlinear algorithms require to construct the nonlinear internal finite element forces and the related tangential operator with the finite element model. Note that such computation is particularly time consuming but is carried out once.

Construction of the projection basis

The construction of the mean nonlinear reduced-order model requires a projection basis. In the present work, such projection basis is computed from the nonlinear reference response with the proper-orthogonal decomposition method (POD-method), which is proved to be particularly relevant for nonlinear problems. Let [A] be the (n×n) correlation matrix related to the nonlinear reference dynamical response U(t). It is defined by

[A]= [Y] T [Y] , [Y] ij = U i (t j ) √ δ t , (4)
in which δ t is the constant sampling time step, and where t j denotes the sampling time number j. The projection basis is defined by the eigenvectors α related to the N most contributing eigenvalues λ α , solution of

[A] α = λ α α with T α β = δ αβ , ( 5
)
in which δ αβ is the Kronecker symbol set to 1 if α = β and to 0 otherwise. Note that in practice, for large computational models, the numerical construction of correlation matrix [A] is difficult to achieve. In such a case, [A] is not computed and the eigenvalue problem is replaced by a singular value decomposition of [Y] or [Y] T .

Stochastic nonlinear reduced-order model

In this work, the nonparametric probabilistic approach of model-form uncertainties recently introduced in [8] is used. In such case, the probability model is directly implemented in the projection basis. Let [Φ] be the (n × N ) matrix whose columns are the eigenvectors 1 , • • • , N .

The matrix [Φ] is then replaced by a stochastic matrix [¨( )], whose probability model is described in [8]. It should be noted that the random vectors contained in random matrix [¨] statisfy the boundary conditions and orthonormality conditions

[B] T [¨]= 0 , [¨] T [¨]= [I N ] , (6)
in which [I N ] is the identity matrix with N order. The probability distribution of random matrix [¨] depends on a R m -valued hyper-parameter = {s, β, [σ]}, whose dimension m = 0.5N (N + 1) + 2 has been optimized according to the modeling. More specifically, s is a scalar, which controls the statistical fluctuations around each deterministic vector basis, β is a scalar that controls the correlation between the components of the vector basis, and [σ] is an (N × N ) real upper triangular matrix, which allows for controlling the correlations between the column vector basis. The Monte Carlo numerical procedure for generating such random basis is briefly given below:

[¨] = [Φ] + s[Z] [I N ] + s 2 [Z] T [Z] -1/2 (7) [Z] = [A] - 1 2 [Φ] [Φ] T [A] + [A] T [Φ] ( 8
) [A] = [U] -[B] [B] T [U] ( 9
) [U] = [G(β)][σ] , ( 10
)
in which [G(β)] is a (n × N ) random matrix whose random generator can be found in details in [8].

The stochastic computational model consists then in solving the set of stochastic nonlinear differential equations

[M] Q(t) + ([D] + [C G (Ω)]) Q(t) + [K(Ω)]Q(t) + [¨] T F NL ([¨] Q(t)) = [¨] T F(t) , ( 11
) in which [M] = [¨] T [M ] [¨], [D] = [¨] T [D] [¨], [C G (Ω)] = [¨] T [C G (Ω)] [¨] and [K(Ω)] = [¨] T [K(Ω)] [¨]
are the full (N × N ) random mass, damping, gyroscopic coupling and stiffness matrices, and where Q is the R N -valued random vector of the generalized coordinates from which the R n -valued random vector of the physical solution U is reconstructed by

U = [¨] Q . ( 12
)
Regarding the numerical computation, the set of these stochastic nonlinear differential equations is solved in the time domain by using the Monte Carlo numerical simulation, using an implicit and unconditionally stable integration scheme (Newmark method with the averaging acceleration scheme). For each sampling time, the nonlinearity is solved iteratively by using either the fixed point method or an arc-length-based algorithm, depending on the nonlinearity rate. Concerning the current strategy, the knowledge of the nonlinear reduced-order internal force vector is required by both nonlinear algorithms whereas the knowledge of the related reduced tangent operator (that is time consuming) is only required by the use of the arc-lengthbased algorithm. Practically, the contributions of both nonlinear reduced-order internal force and tangent reduced operator induced by each finite element are explicitly constructed with the finite element method. The assemblage over the finite elements is particularly quick because it only requires to sum each of these contributions. A posterior nonlinear dynamical analysis is performed from Eq.( 12) in the frequency domain by computing the Fourier transform of the nonlinear dynamical response.

The main advantage of this new method is that it allows for considering general nonlinear dynamical systems for which the operators are not necessarily known. It also constitutes a global approach because there is no need to analyze and decide whether the mass, the damping or the stiffness matrices are random or not. In opposite, the real difficulty arises when considering the inverse problem because the hyperparameter that has to be identified has a large dimension that increases with the square of the order of the reduced-order model. This latter point requires some further open researches that allows for reducing the size of the hyperparameter and that are not investigated in the present research.

Numerical application

The structure under consideration is a bladed-disk structure of order M = 24 whose finite element mesh is constructed with hexahedral solid finite elements with 8 nodes and is constituted of 52 290 degrees of freedom. The computational model is issued from [13] for which the boundary conditions have been modified. The structure is in rotation around its revolution axis with a constant velocity Ω = 1 000 rpm. Since the nonlinear dynamic analysis is carried out in the rotating frame of the structure, the rigid body motion due to the rotation of the structure corresponds to a fixed boundary condition located at the inner radius of the structure. The bladed-disk is made up of a homogeneous isotropic linear elastic steel material. A Rayleigh damping model is added for the bladed-disk structure such that [D]= α[M ]+β[K] with α = 39.27, β = 1.53 × 10 -6 such that fundamental natural frequency ν 1 = 383 Hz has a critical damping ξ 1 = 0.01.

Description of the excitation

The cyclic symmetry is first used for solving the linear generalized eigenvalue problem related to the tuned rotating bladed-disk. Figure 1 displays the natural frequencies of the tuned structure ν i with respect to the circumferential wave number h. From this graph, we are interested by the veering corresponding to circumferential wave number h = 4 in the frequency band [2 500 , 3 200] Hz. As a consequence, the spatial repartition of the external load described by the normalized vector F is such that an external point load is applied along all directions at the excitation nodes located at the tip of each blades with a constant phase shift π 3 . In this case, function g(t) is written as

g(t) = 2 π t sin(π ∆ν t) cos(2π s ∆ν t) with ∆ν = 700 Hz , s = 4.0714 , (13)
such that | g(2πν)| = 1 B e (2πν) as represented in Figure 2. Practically, the computation is carried on a truncated time domain T = [t ini , t ini +T ]. The initial load is chosen as t ini = -0.057 s yielding a null initial load. The time duration T is then adjusted so that the system be returned at its equilibrium state within a given numerical tolerance for both linear and nonlinear computations.

Even the fundamental natural frequency does not belong to excitation frequency band B e , it can be indirectly excited through the geometrical nonlinear effects. Time duration is chosen as T = 0.128 s which ensures the system to return to its equilibrium state with a relative tolerance τ = e -2πξ 1 ν 1 T = 4.6% when fundamental natural frequency is excited. The sample frequency ν e and the number of time steps are then chosen as ν e = 8 000 Hz and n t = 1 024 yielding a constant sampling time step δt = 1.25 × 10 -4 s and a constant sampling frequency step δν = 3.907 Hz.

Finally, the external loading can be written as

F(t) = s 0 g(t) F , ( 14
)
in which s 0 is the load intensity. In the present case, the case s 0 = 28.184 N is investigated and corresponds to the case in which geometrical nonlinearities deeply affect the vibrational behavior of the bladed-disk structure.

Computational results

In this section, the nonlinear finite element response described in Eqs. ( 2) and( 2) is computed and taken as a reference response in the frequency band of analysis B a = [350 , 4000] Hz. Such nonlinear reference response is then used in order to construct the projection basis as explained in Section 1. A convergence analysis is then performed in order to set the optimal order of the reduced-order model that can accurately reproduce the reference response. In this case, it is shown that the a good agreement is obtained with N = 42. The computational results are presented for s ∈ {0.01, 0.0125, 0.015, 0.0175}, β = 0.0025 and [σ] an upper triangular matrix with N order for which the non-zeros entries are defined by [σ] ii = 10, [σ] i,i+1 = 5 and [σ] i,i+2 = 2.5.

Note that these results corresponds to a weakly mistuned situation for which amplification effects are still moderate when only considering the linear effects. The results are computed with n s = 70 Monte Carlo numerical realizations. Let u j α (t) be the displacement at the tip of the blade number j along direction e α corresponding to the reference solution. We are then interested by the observation defined by v(2πν)= max j || u j α (2πν)||. Figures 3,4, 5, and 6 display for s = {0.01, 0.0125, 0.015, 0.0175}, the graph of function ν → v(2πν) (black dashed line) and the graph of the confidence region (yellow area) of random function ν → V (2πν) of the observation with a probability level p c = 0.95 related to the nonlinear reduced-order model. In these graphs, the sensitivity of the stochastic nonlinear forced response is analyzed with respect to the uncertainty level. As expected, it can be seen that the structure responds outside B e . More particularly, the resonance located around 2,100 Hz has similar amplitude to the main resonances located in B e and matches well to its deterministic counterpart. The resonance located around 1,000 Hz cannot either be neglected. Its robustness to uncertainty is sensitively weak and a frequency shift is observed with respect to the reference response.

An interesting point concerns the particular unexpected robustness of the main resonances that are located in B e . In the present case, the structure clearly undergoes significant geometrical nonlinear effects since it responds outside the excitation frequency band. One can see that such effects have a tendency to stabilize the dynamical behaviour in B e . On the contrary, the resonances issued from the propagation of the geometrical nonlinearities have a weaker response level accompanied by a very low response level. It can be seen as a "linear" dynamical behaviour that is known to be sensitive to mistuning uncertainties. 0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 10 -10

Conclusions

A computational methodology has been presented for analyzing the propagation of uncertainties in the context of structural dynamics with geometric nonlinearities. A nonparametric probabilistic method for modeling model-form uncertainties, which is implemented from a given projection basis is used for constructing a stochastic computational nonlinear reduced-order model . A numerical application, which consists in a rotating bladed-disk structure, is presented and shows the capability of the stochastic computational model to capture the unexpected resonances induced by the geometrical nonlinearities.

Figure 1 .Figure 2 .

 12 Figure 1. Graph of natural frequencies ν i with respect to its circumferential wave number h. Excitation frequency band B e (light red area)

Figure 3 .Figure 4 .Figure 5 .Figure 6 .

 3456 Figure 3. Graph of functions ν → v(2πν) (black dashed line) computed with the deterministic nonlinear finite element model, graph of the confidence region of random function ν → V (2πν) computed with the nonlinear stochastic reduced-order model (yellow area), localization of the excitation frequency band B e (light red area) for s = 0.01.

Acknowledgement

The authors thanks Pr. Christophe Pierre from the Illinois University concerning the use of the finite element model and also thank Pr. Marc Mignolet from the Arizona State Unversity.


Anthony Picou
email: anthony.picou@u-pem.fr

Evangéline Capiez-Lernout

Christian Soize

Moustapha Mbaye
email: moustapha.mbaye@safrangroup.com

Evang Éline Capiez-Lernout

Mistuning analysis of a

come

INTRODUCTION

In the framework of linear dynamic analysis of bladed-disk, the intentional detuning has been identified as an efficient technological way for reducing the sensitivity of the forced response of bladed-disks to mistuning induced by small variations of the mechanical properties [START_REF] Ewins | The effects of detuning upon the forced vibrations of bladed disks[END_REF][START_REF] Whitehead | The maximum factor by which forced vibration of blades can increase due to mistuning[END_REF]. Various research concerning the modeling of mistuned bladed-disk have been carried out, with the development of reduced-order models (see for instance [START_REF] Castanier | Reduced order modeling technique for mistuned bladed disks[END_REF][START_REF] Yang | A reduced-order model of mistuning using a subset of nominal system modes[END_REF][START_REF] Bladh | Component-mode-based reduced order modeling techniques for mistuned bladed disks-part 1: The-oretical models[END_REF][START_REF] Sinha | Reduced-order model of a bladed rotor with geometric mistuning[END_REF] for linear cases and [START_REF] Capiez-Lernout | Mistuning analysis and uncertainty quantification of an industrial bladed disk with geometrical nonlinearity[END_REF][START_REF] Grolet | Computing multiple periodic solutions of nonlinear vibration problems using the harmonic balance method and groebner bases[END_REF] for nonlinear cases) and of probabilistic approaches [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Mignolet | Nonparametric stochastic modeling of linear systems with prescribed variance of several natural frequencies[END_REF]. The pattern optimization has been proved to be efficient for reducing the response amplitudes in the linear context. Nevertheless, geometrical nonlinear effects can no longer be neglected with the current technological improvements due to the use of more flexible and lighter materials which sometimes leads to large displacements and strains.

The main objective of this paper is to investigate the effects of geometrical nonlinearities on the nonlinear dynamical behavior of a bladed-disk that is both mistuned and intentionally mistuned. The intentional mistuning is introduced by a given pattern that defines the spatial distribution of different types of sectors on the full bladed-disk. The mistuning is caused by all the discrepancies that are encountered in the manufacturing process and that break the cyclic symmetry of the structure which may cause dynamical amplification effects. For the sake of clarity "intentional mistuning" and "mistuning" are respectively denoted by "detuning" and "mistuning". We are then interested in quantifying the nonlinear dynamical behavior of the detuned-mistuned bladed-disk assuming the intensity of the load sufficiently high to rightly consider large displacements and large strains.

The main novelty of this work is that both detuning and mistuning phenomena are taken into account within the computational model. The detuning modeling requires to consider two different generating sectors whose finite element meshes are compatible from one sector to another one and to assemble the full structure according to a given pattern. The mistuning phenomena is taken into account through the nonparametric probabilistic approach of uncertainties [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF] extended to the geometrical nonlinear context as proposed in [START_REF] Capiez-Lernout | Mistuning analysis and uncertainty quantification of an industrial bladed disk with geometrical nonlinearity[END_REF]. An adapted methodology in accordance with those objectives is implemented. We are therefore interested in the physical nonlinear analysis of such detuned and mistuned bladed-disk structure.

The first Section summarizes the main methodological steps of the nonlinear analysis of such detuned and mistuned bladed-disk structures. Details are given on the choice of the projection basis used to build the nonlinear reduced-order model that is required [START_REF] Ewins | The effects of detuning upon the forced vibrations of bladed disks[END_REF] to optimize the dimension of the nonlinear reduced-order model and (2) to implement the mistuning uncertainties according to the nonparametric probabilistic approach keeping a reasonable size of the random germ matrix. The second Section is devoted to a numerical application for which a simple model representing an industrial detuned and mistuned bladed-disk is considered. In particular, a mistuned case and a mistuned-detuned case are considered. The nonlinear dynamical amplification factors resulting from a uniform excitation over a narrow low-frequency band of excitation are then compared for the two mistuned cases, inside and outside of the frequency band of excitation.

METHODOLOGY

This section is devoted to the construction of a stochastic nonlinear reduced-order model. The constitutive equations of the materials are linear elastic. The external loads are sufficiently large so that the structure undergoes geometrical nonlinear effects induced by the large displacements and strains.

The main difficulty is to obtain a nonlinear stochastic computational model that is able to take into account both detuning and mistuning phenomena considering industrial high-dimensional computational model (HDM). This means that the nonlinear stochastic computational model has to be optimally reduced in order to perform nonlinear stochastic computations with a reasonable computational cost. With the current use of the nonparametric probabilistic strategy, a deterministic nonlinear reduced-order computational model has to be constructed for representing the detuned bladed-disk structure. Uncertainties are then added to this mean nonlinear reduced-order model in order to model the mistuning phenomena. The main steps are as follows:

• A first projection modal basis is computed by solving a generalized eigenvalue problem related to the HDM.

• A second basis is obtained using the Proper-Orthogonal Decomposition (POD) method performed on the structural displacements that are computed with the nonlinear reduced-order model resulting from the projection of the HDM on the first modal basis. • Mistuning uncertainties are then added to the mean nonlinear reduced-order model obtained from the double projection method.

Nonlinear finite element computational model of a detuned bladed-disk

The detuning is done according to a given pattern, defining the spatial repartition of several different sector types on the full bladeddisk. The finite element meshes of each sector type are compatible from one sector interface to another one. In the rotating frame, the nonlinear differential equation that describes the computational model of the nonlinear dynamical behavior of a rotating detuned bladed-disk as a function of the rotational speed Ω is written as

[M] ü + ( [D] + [C g (Ω)] ) u + ( [K e ] + [K c (Ω)] + [K g (Ω)] ) u + F NL (u) = F , (1)
in which the R n -vector u(t) is the vector of unknown displacements. In equation ( 1), the

(n × n) real matrices [M], [D], [C g (Ω)], [K e ], [K c (Ω)], [K g (Ω)
] are respectively the mass matrix with positive definiteness properties, the damping matrix with symmetry and positive definiteness properties, the gyroscopic coupling matrix with skew-symmetry properties, the elastic stiffness matrix with symmetry and positive definiteness properties, the centrifugal stiffness matrix with symmetry and negative definiteness properties, the geometrical stiffness matrix with symmetry and positive definiteness properties. It is assumed that rotational velocity is sufficiently large so that

[K e ] + [K c (Ω)] + [K g (Ω)
] is positive definite and there is no static unstability (divergence). In equation ( 1), the R n -vector F results from the discretization of the external forces and the R n -vector F NL (q) describes the nonlinear internal forces issued from the geometrical nonlinearities.

Numerical construction of an adapted projection basis

A possible approach would consist in solving the generalized eigenvalue problem related to the conservative and linear part of the differential equation [START_REF] Ewins | The effects of detuning upon the forced vibrations of bladed disks[END_REF]. Nevertheless, such generalized eigenvalue problem would yield complex eigenvectors because of the skewsymmetry property of the gyroscopic coupling. To avoid this difficulty, the projection basis is built in two steps. First, the previous generalized eigenvalue problem is solved without damping and gyroscopic coupling matrices. Let ϕ 1 , . . . , ϕ m be the eigenmodes related to the m first eigenvalues, 0

< λ 1 λ 2 • • • λ n are stored in the (n × m) real modal matrix [Φ m ].
Such projection modal basis is then used to build the nonlinear reduced-order model related to equation [START_REF] Ewins | The effects of detuning upon the forced vibrations of bladed disks[END_REF]. A convergence analysis is carried out with respect to m. Once convergence is reached, the nonlinear reduced-order model is used for computing the forced time-response denoted q 0 (t). The POD basis is obtained by solving the eigenvalue problem for the correlation matrix related to the nonlinear time response q 0 (t). The eigenvectors related to the N most contributing eigenvalues are retained to form the POD vector basis. Let [W (m,N) ] be the (m × N) real matrix whose columns are these POD vector basis. As explained in [START_REF] Capiez-Lernout | Post-buckling nonlinear static and dynamical analyses of uncertain cylindrical shells and experimental validation[END_REF], this POD vector basis is efficiently generated by computing the N left singular vectors of the real matrix [V ] defined by

[V ] i j = q 0 i (t j ) √ ∆t . (2)
The projection basis of the HDM, which allows for constructing the nonlinear reduced-order model, is represented by the (n × N) real matrix [Ψ] that is written as

[Ψ] = [Φ m ] [W (m,N) ] . (3)
With such reduced-order basis, the nonlinear reduced-order model is written as

u(t) = [Ψ] q(t) , (4)
[M] q(t)

+ ( [D] + [C(Ω)] ) q(t) + [K(Ω)]q(t) + F NL (q(t)) = F(t) , (5)
in which

[M] = [Ψ] T [M][Ψ] (6)
[D] = [Ψ] T [D][Ψ] (7)
[C(Ω)] = [Ψ] T [C g (Ω)][Ψ] (8)
[K] = [Ψ] T ([K e ] + [K c (Ω)] + [K g (Ω)]) [Ψ] (9)
F(t) = [Ψ] T F(t) . ( 10
)
The reduced quadratic and cubic stiffness operators K

i jk , {i, j, k} ∈ {1, 2, . . . , N} 3 and K

i jk , {i, j, k, l} ∈ {1, 2, . . . , N} 4 are explicitly constructed [START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in three-dimensional elasticity[END_REF] and the R N -vector F NL (q) is then calculated using

F NL i (q) = K (2)
i jk q j q k + K

i jk q j q k q , ∀q ∈ R N .

After solving the nonlinear differential equation ( 5), for which q(t) is the R N -vectors of the unknown generalized coordinates, the physical response is retrieved by using equation [START_REF] Yang | A reduced-order model of mistuning using a subset of nominal system modes[END_REF].

Construction of the stochastic reduced-order model

The stochastic nonlinear computational model corresponds to a probabilistic modeling of the mistuning for the detuned rotating bladed-disk. As previously explained, the nonparametric probabilistic approach of uncertainties is used [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF]. Its generalization to the nonlinear geometrical context [START_REF] Mignolet | Stochastic reduced order models for uncertain geometrically nonlinear dynamical systems[END_REF], including a strategy for reducing the size of the germs [START_REF] Capiez-Lernout | An improvement of the uncertainty quantification in computational structural dynamics with nonlinear geometrical effects[END_REF] is used. It involves a positive matrix [K NL ] with order P = N(N + 1) containing entries contain all the cubic, quadratic and linear elastic stiffness contributions. The corresponding random matrix [K NL ] is written as:

[K NL ] = [L K ][G K (δ K )][L K ] T + [∆ K] (12)
in which [L K ] is a (P × Q) real matrix whose columns contain Q, the predominant eigenvectors of matrix [K NL ], where

[∆ K] = [K NL ] - [L K ][L K ]
T and where [G K (δ K )] is the full (Q × Q) random matrix with Q ≪ P, constructed using the Maximum Entropy principle [START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF] and where δ K is the hyperparameter that allows for controlling the level of uncertainties. By construction the mathematical expectation,

E{[K NL ]} of random matrix [K NL ] is the deterministic matrix [K NL ].
All the details concerning the construction of [K NL ] can be found in [START_REF] Mignolet | Stochastic reduced order models for uncertain geometrically nonlinear dynamical systems[END_REF].

NUMERICAL APPLICATION

The numerical application concerns the nonlinear dynamical analysis of a detuned and mistuned rotating bladed-disk. From an industrial point of view, such an application is of industrial interest when exceptional operating ranges such as severe loads or unstable flutter situations yielding low-damping levels are considered.

Description of the finite element model

The bladed-disk comprises 24 blades and is made of an homogeneous and isotropic elastic material (steel) with Young modulus 2 × 10 11 N.m -2 , Poisson's ratio 0.3, and mass density 7650 Kg.m -3 . In the rotating frame, the disk is clamped at the inner radius of the disk sector. The geometrical characteristics of the generating sector are summarized in Table 1. Software ANSYS has been used to construct the computational model of the reference sector, which is constituted of tridimensional solid finite elements: brick elements with 20 nodes, pyramidal elements with 13 nodes, and tetrahedral elements with 10 nodes. The computational model of the tuned bladeddisk (with 24 blades) has been constructed using the computational model of the sector. The characteristics of the finite element model are given in Table 2 1. The mesh of the finite element sector is shown in Figure 2. Eigenfrequencies of the linearized tuned rotating bladed-disk Figure 3 displays the Campbell diagram representing the evolution of the eigenfrequencies according to the rotational velocity Ω for the tuned case (pattern P 1 ). Figure 4 displays the graph of eigenfrequencies ν α = √ λ α as a function of the nodal diameter number h for the linearized tuned rotating bladed-disk (configuration P 1 ) for which the rotational speed is Ω = 465 rad.s -1 .

P 1 = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1] ( 13
) P 2 = [1 1 1 1 1 1 2 2 2 2 2 2 1 1 1 1 1 1 2 2 2 2 2 2] ( 14

Defining the external load

An h = 20 engine order excitation is chosen for which there are two close eigenfrequencies from two different classes of eigenmodes. In the time domain, the load vector F(t) is defined by,

F(t) = s 0 g(t) F . ( 15
)
The time-function g(t) is defined for t ∈ R and is constructed so that its Fourier transform is | g(2πν)| = 1 on the frequency band B e = [1000, 1600] Hz of excitation and zero outside B e . In the numerical process, the signal g(t) is truncated by choosing t ini = -0.065 s such that g(t ini ) = 0 with a time duration T = 0.35 s, ensuring the low frequencies located outside B e to be captured when considering nonlinear dynamical responses. Figures 5 and6 show the graph t → g(t) and its Fourier transform ν → | g(2πν)|. The spatial distribution of the external load is a normalized vector F for which all directions related to the node located at the tip of each blade are simultaneously excited and the load intensity is given by s 0 = 1 N. The dynamical analysis is performed in the time domain. Then a Fourier Transform of the time responses allows for analyzing the nonlinear dynamical responses in the frequency domain.

NUMERICAL PARAMETER CONCERNING THE STOCHASTIC REDUCED-ORDER MODEL

We are interested in analyzing the nonlinear dynamical response in the frequency band of analysis B = [0, 4 000] Hz. For each pattern, the first generalized eigenvalue problem is solved and a convergence analysis of the nonlinear dynamical responses showed that m = 145 eigenmodes were required for the first reduced-order basis. The convergence analysis with respect to the responses computed with the first nonlinear reduced-order model shows that N = 55 vectors in the second reduced-order basis (constructed by the POD method) are sufficient to reproduce the nonlinear dynamical behavior in frequency band B [START_REF] Picou | Effects of geometrical nonlinearities for a rotating intentionally mistuned bladed-disk[END_REF]. In the following, all calculations are carried out with a reduced-order basis characterized by (m, N) = (145, 55). Concerning the size of the random germ, a convergence analysis has shown that a good approximation was obtained for Q = 500 N(N + 1) = 3 080.

Nonlinear dynamical analyses of both tuned and detuned rotating bladed-disk

In this Section, the nonlinear deterministic dynamical response of the tuned and detuned bladed-disk is considered. In the following, the subscripts NL and LIN are used for the nonlinear case and for the linear case (when a quantity is used either for NL or for LIN, the Circumferential wave number subscript is removed). Let u j α (t) be the time displacement of blade j according to local direction e α . For a given pattern, the observation |||u j ||| is defined by

|||u j ||| 2 = max t ||u j (t)|| 2 with ||u j (t)|| 2 = 3 ∑ α=1 |u j α (t)| 2 . ( 16
)
For each pattern, the response related to the most responding blade j 0 in the time domain is considered. We then have 2,NL (t) for pattern (P 1 ) and pattern (P 2 ). Table 3 summarizes the maximum displacements for both linear and nonlinear cases corresponding to pattern P 1 and P 2 Linear Nonlinear It should be underlined that the major effects of the geometrical nonlinearities is to induce an indirect excitation outside the frequency excitation frequency band as shown in [START_REF] Capiez-Lernout | Mistuning analysis and uncertainty quantification of an industrial bladed disk with geometrical nonlinearity[END_REF][START_REF] Picou | Effects of geometrical nonlinearities for a rotating intentionally mistuned bladed-disk[END_REF] in turbomachinery context or in [START_REF] Akkaoui | Nonlinear dynamical analysis of a fluidstructure computational model with sloshing and capillarity[END_REF] in fluide-structure interaction.

j 0 = arg max j |||u j ||| . ( 17

Nonlinear dynamical analysis of the mistuned-detuned rotating bladed-disk

Unlike the previous Section where the observations related to the deterministic case were denoted by lowercase letters, the random observations corresponding to the stochastic case (presence of mistuning) are similarly denoted by capital letters. In this part, we are interested in the nonlinear dynamical stochastic analysis of the mistuned-detuned rotating bladed-disk. Let B NL (2πν) be the random amplification dynamic factor defined by,

B NL (2πν) = || U j 0 NL (2πν)|| ||| u j 0 ,tuned NL ||| , (18)

CONCLUSION

A methodology adapted to the modeling of the mistuned and the detuned rotating bladed-disk has been detailed and applied to a simple model representing an industrial bladed-disk. The results highlight the indirect excitation of the structure through the geometrical nonlinearities, above and below the excitation frequency band. Furthermore, the dynamical analysis of the considered detuned and mistuned rotating bladed-disk has allowed the sensitivity of the structure outside the excitation frequency band to be quantified. The computational implementation of the proposed methodology shows that we are able to propose a nonlinear geometrical analysis of a detuned-mistuned rotating bladed-disk structure. In the present case, the considered detuned pattern is not optimized regarding the dynamical amplification factor. An interesting perspective is to perform an optimization of the detuning with respect to dynamical criteria involving the amplification factor that occur in the excitation frequency band as well as the local amplification factors that occurs outside the excitation frequency band.

ACKNOWLEDGMENT

SafranTech is gratefully acknowledged for the permission to publish this work.

  NOMENCLATUREb NL Deterministic dynamical amplification factor g(t) Load in time domain g(2πν) Fourier transform of g(t) h Nodal diameter number j Blade number j 0 Blade of interest m Number of eigenvectors n Number of degrees-of-freedom s 0 Load intensity t Time t ini Initial time u j α (t) Displacement at the tip of blade number j according to local direction e α u j α (2πν) Fourier transform of u j α (t) B(2πν) Random dynamical amplification factor N Dimension of the nonlinear reduced-order model P Dimension of [K NL ] matrix

FIGURE 1 :

 1 FIGURE 1: Configurations P 1 (left) and P 2 (right)

FIGURE 2 :FIGURE 3 :

 23 FIGURE 2: Mesh of the finite element sector

FIGURE 4 :FIGURE 5 :

 45 FIGURE 4:Graph h → ν α (h) of the eigenfrequencies ν α related to the tuned pattern P 1 of the rotating bladed-disk as a function of the nodal diameter number h.

P 1 1 .TABLE 3 :

 13 48 × 10 -3 m 5.33 × 10 -4 m P 2 1.67 × 10 -3 m 6.18 × 10 -4 m Maximum displacements in the time domain and related to a load intensity s 0 = 1 N. It is clearly seen that this load yields to significant geometrical effects that stiffen the nonlinear dynamical response of the tuned and detuned structures. The observed nonlinear time responses look more irregular, suggesting numerous resonances contributions outside B e . The nonlinear dynamical response is then analyzed in the frequency domain. Let u j 0 α (2πν) be the Fourier Transform of u j 0 α (t).

Figure 11

 11 displays the graph of ν → || u j 0 LIN (2πν)|| (red line) and the graph of ν → || u j 0 NL (2πν)|| (blue line) related to pattern P 1 . Figure 12 displays the graph of ν → || u j 0 LIN (2πν)|| (red line) and the graph of ν → || u j 0 NL (2πν)|| (blue line) corresponding to pattern P 2 . For both patterns, secondary resonances induced by the geometrical nonlinear effects appear below and above excitation frequency band B e .

2 FIGURE 7 : 2 FIGURE 8 :

 2728 FIGURE 7: Linear dynamical analysis in the time domain: graph of t → u j 0 2,LIN (t) corresponding to pattern P 1 .

2 FIGURE 9 : 2 FIGURE 10 :

 29210 FIGURE 9: Nonlinear dynamical analysis in the time domain: graph of t → u j 0 2,NL (t) corresponding to pattern P 1 .

FrequencyFIGURE 11 :

 11 FIGURE 11: Nonlinear dynamical analysis of the deterministic tuned structure (pattern P 1 : graphs of ν → || u j 0 LIN (2πν)|| (red line) and ν → || u j 0 NL (2πν)|| (blue line) and excitation frequency band B e (light blue area).

Frequency

FIGURE 12 :

 12 FIGURE 12: Nonlinear dynamical analysis of the deterministic detuned structure (pattern P 2 : graphs of ν → || u j 0 LIN (2πν)|| (red line) and ν → || u j 0 NL (2πν)|| (blue line) and excitation frequency band B e (light blue area).

Frequency

FIGURE 13 :Frequency

 13 FIGURE 13: Nonlinear mistuning analysis of the detuned structure (pattern P 2 ): graph of ν → b NL (2πν) (thick dashed line), confidence region of B NL (2πν) (yellow area) and excitation frequency band B e (light blue area) for δ = 0.1.

FIGURE 14 :

 14 FIGURE 14: Nonlinear mistuning analysis of the detuned structure (pattern P 2 ): graph of ν → b NL (2πν) (thick dashed line), confidence region of B NL (2πν) (yellow area) and excitation frequency band B e (light blue area) for δ = 0.3.

TABLE 1 :

 1 . A detuned bladed-disk structure is constructed by defining another sector type (denoted 2), which is obtained Geometrical characteristics of the reference sector by decreasing the Young modulus of the blade. The Young modulus of the disk remains equal to 2 × 10 11 N × m -2 . Consequently, two patterns are considered: a cyclic one (tuned) denoted as P 1 and a detuned one denoted as P 2 . The finite element model of tuned configuration P 1 and the detuned configuration P 2 are shown in Figure

	Inner disk Radius	19.8 mm
	Outer disk Radius	100 mm
	Disk width	20 mm
	Blade thickness at root section 4.8 mm
	Blade thickness at tip section	2 mm

  )

	Structure	Elements	Nodes	DOFs
	Blade sector	2,714	6,896	20,688
	Disk sector	836	4,554	13,662
	Full structure	85,200	265,080 787,176

TABLE 2 :

 2 Characteristics of the finite element model

  LIN (t) for both patterns P 1 and P 2 . It can be seen that both tuned and detuned structures yield different responses in the time domain. Figures 9 and 10 respectively display the graphs of t → u

	10 2
	jĝ(2:8)j 10 0
	10 -2
	10 -4	0	1000	2000	3000	4000
			frequency 8 (Hz)
			FIGURE 6: Graph of ν → | g(2πν)|.
	Figures 7 and 8 display the graphs of t → u	j 0 2,j 0
						)