C Pozzolini
email: cedric.pozzolini@cnes.fr

F Buffe

M Salaun

A LOCALIZATION AND UPDATING STRATEGY OF LARGE FINITE ELEMENT MODELS IN STRUCTURAL DYNAMICS

Keywords: Model updating, Error localization, Visibility, Constitutive Relation Error

The purpose of this paper is to evaluate the application of the error of constitutive law method to the updating of large FE models of space structures using FRF experimental results. First, we briefly recall the theoretical basis of this method in modal and frequency approaches. Then, the notion of visibility is introduced to improve the modelling of localization error and the quality of modal updating, for low frequencies. Finally we propose a global strategy and discuss the results we obtained on satellite JASON2.

Introduction

During the launch of a satellite, the level of excitation is very high. However, the security margins must be reduced in order to remain competitive. So CNES is concerned with the problem of Finite Element Model (FEM) updating in lower frequency to ensure that numerical simulations are predictive enough to maintain a high reliability of spacecraft structures. This problem has been a subject of intense research over the past twenty years, but the introduced methods are still hard to apply on industrial structures. Usually, the validation of a numerical dynamical model is performed by comparing numerical eigendata or Frequency Response Functions (FRF) with natural data measured from sinesweep base excitation vibration tests. In FEM analysis, it is usual to make simplifying assumptions, in the modelling of joints for example, or to neglect some phenomena, such as contact between structures. Therefore, the purpose of modal updating is to modify parameters (such as mass, stiffness and damping of sub-structures or connections between components) in the numerical model in order to obtain better agreement between numerical and experimental data. To select erroneous parameters, a localization criterion is applied. It should be noted that the conditioning of the updating inverse problem is deteriorated when the number of erroneous parameters is too large or when the number of sensors is not large enough. That's why a sensitivity Monte Carlo analysis is used to remove non-significant parameters. Finally, the udapting process is assumed to be correct when the Modal Assurance Criterion (MAC) is accurate enough.

The modelling of error localization uses the existing error of constitutive law. In particular, the applicability of such a method to large industrial test cases such as JASON2 is adressed. In order to improve the quality of modal updating of the low frequency behaviour, a comparison between FRF and modal approaches of the error of constitutive law is made. Error location is extensively studied in the case of the dynamical behaviour of a clamped structure, with realistic limitations in order to find a better linear equivalent model. Design variables are typically thickness of panels and section or number of beams. Particular attention is paid to automatic selection of relevant parameters during updating in order to keep the sensitivity matrix well-conditioned and to obtain an accurate solution for parameters fitting.

Basics on error of constitutive law methods

The classical problem of linear elastodynamics to modelize a structure, which occupies a domain Ω of R 3 , is to find a displacement field u kinematically admissible (K.A.) and a stress σ such as ∀v K.A. ,

Ω f vdΩ = Ω tr σ(u)ε(v) dΩ - Ω ρ ∂ 2 u ∂t 2 vdΩ,
where ρ is the mass density and ε(.) stands for the linearized strain tensor. The material constitutive law is σ(u) = Hε(u). In sine-sweep base forced vibration, boundary conditions are given, on a part Γ of ∂Ω, by

∂u ∂t (., t) = 0 u(., t) = - a ω 2 sin(ωt)
The associated discrete problem reads

[K]u + iω[D]u -ω 2 [M]u = F (1)
where [K], [D] and [M] are respectively the stiffness, damping and mass matrices. For spacecraft structures, the damping is considered as small enough to use the so-called Basile hypothesis :

[D] = α[K] + β[M].
Then, all these matrices depend on geometric and mechanical parameters (sections, thicknesses, Young modulus, Poisson ratio and density). The Constitutive Relation Error (CRE) was initially proposed by Ladavèze et al in the early 1980's as an error estimator for FEM [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF]. An extended version was described by Chouaki et al [START_REF] Chouaki | Updating structural dynamic models with emphasis on the damping properties[END_REF] in order to update mass, stiffness, and damping properties using the FRF data. In [?] is introduced a distance between two admissible displacement fields u and v in order to compare two different dynamical models of Ω. This distance in material constitutive law reads

d(u, v) 2 = Ω tr CL(u, v)H -1 CL(u, v) dΩ where CL(u, v) = σ(u) -Hε(v).
In the industrial framework, we have a discrete numerical model of the structure and the measured data are a set of eigenmodes {ω 2 k , u k } k or a set of displacement fields {u ω } ω∈ [ωmin,ωmax] . Let us remark that these eigenvectors or displacement fields are only given at a few nodes of the mesh, more precisely where the sensors are located. So, the first step of the updating is to "expand" the measures to the full set of mesh nodes. In the case of an eigenmode {ω 2 k , u k } this expansion, say u k , is the displacement field, denoted below by u, which minimizes the following function over the admissible displacements set

J k (u, v, w) = γ 2 (u -v) t [K](u -v) + 1-γ 2 (u -w) t ω 2 k [M](u -w) + (πu -u k ) t [K] π (πu -u k ) under the constraint [K]v = ω 2 k [M]w.
(2) In this equation, πu denotes the restriction of u to the measured DOFs (sensors), while [K] π is the corresponding restriction of the rigidity matrix. Moreover, γ is a scalar parameter which value is generally taken as 0.5.

Remark 2.1. Let φ k be the k th -numerical mode: [K]φ k = ω 2 k [M]φ k . Clearly, if ω k = ω k and u k = φ k , then J k (φ k , φ k , φ k ) = 0 and u k = φ k , as- sociated with v k = w k = φ k , is solution of the previous minimization problem.
For FRF data with damping, the expension of u ω is obtained by minimizing

J ω (u, v, w) = γ 2 (u -v) t [K] + T ω 2 [D] (u -v) + 1-γ 2 (u -w) t ω 2 [M](u -w) + (πu -u ω ) t [G] π (πu -u ω ) such that [K]v + iω[D]v -ω 2 [M]w = F (3) where [G] π = γ 2 [[K] π + T ω 2 [D] π + 1 -γ 2 ω 2 [M] π ,
scalar T being a characteristic time, taken as the time length of measures for example. Then, the updating problem belongs to the large class of inverse problems. To solve it, we assume that modelization errors are located where error of constitutive law is the highest. That's why an iterative algorithm is used, for which each iteration needs two steps: error location and correction. For this, we use a local indicator on each sub-structure of Ω, say Ω i , which reads

ξ 2 i = 1 m k≤m eli∈Ωi ξ 2 λ k ,eli
, in which el i stands for the finite elements belonging to Ω i , the elementary terms being

ξ 2 λ k ,eli = γ 2 (u k -v k ) t [K] eli (u k -v k ) + 1-γ 2 ω 2 k (u k -w k ) t [M] eli (u k -w k )
ξ i can be normalized by the total energy on Ω, the local energy or the volume of Ω i . A drawback of total energy normalization is that sub-structures with relatively small modeling errors but containing a large part of total structural energy may be associated with an indicator ξ i much greater than a really erroneous sub-structure. For example, it can be the case when the number of sensors is relatively small. Then, the residual energy is spread over the structure and correctly modeled sub-structures with high energy may still be localized as erroneous.

The parameter correction process is typically performed by solving the non-linear optimization problem around (u k , v k , w k ) solu-tions of ( 2)

min p∈Pj k γ 2 (u k -v k ) t [K](p)(u k -v k ) + 1-γ 2 (u k -w k ) t ω 2 k [M](p)(u k -w k ) +(πu k -u k ) t [K] π (p)(πu k -u k ) (4)
where the sum is taken on all eigenmodes. p is a vector of active parameters, belonging to the set P j defined at the j th localization step. In a same way, for (u, v, w) solutions of (3), the minimization problem reads

min p∈Pj ω k 1-γ 2 (u -w) t ω 2 k [M](p)(u -w) + γ 2 (u -v) t [K](p) + T ω 2 k [D](p) (u -v) +(πu -u ω ) t [G] π (p)(πu -u ω )
(5) the sum being taken over the set of frequencies ω k which are identified with FRF in [ω min , ω max ]. It is important to note that the extended CRE formulation does not require a a priori pairing between analytical and experimental eigemodes. It is worth noticing that this last method has been tested only for small academic problems and not for industrial structures, such as JASON2.

Sub-structure visibility

This section deals with a generalization of visibility notion, introduced in [START_REF] Michot | Contribution à la validation de modèles en dynamique des structures[END_REF]. Let us just remark that this quantity was introduced in the static case whereas we develop it for dynamics. For a given FEM, let n be the number of DOFs of the mesh, n i the number of DOFs of sub-structure i, n c the number of measured DOFs and n ss the number of sub-structures. Then we introduce space F i , containing load vectors applied to DOFs of sub-structure i: more precisely, F i is a n i -dimensional subspace of R n , spanned by "unit vector fields" which values are 1 for one DOF of sub-structure i and 0 for all the others. Similarly, we define space F c , associated with loads applied to the measured DOFs: here again, F c is spanned by "unit vectors" which values are 1 for one measured DOF and 0 for all the others. So it is clearly n cdimensional. Now, let us introduce space U i of displacements due to any loading on i thsub-structure. Therefore, U i is the image of

F i by [K] -1 .
Our goal is to measure the sensors capability to localize defaults in material constitutive law on any sub-structure. In this respect, we introduce space U ck such that

[Z] k (U ck ) = F c ⊂ R n , with [Z] k := [K] -ω 2 k [M]
. So U ck can be seen as set of displacement fields for which sensors are excited, in the case of the k theigenmode. A measure of the visibility of sub-structure i for mode k is linked with the capability of space U i to be properly represented by space U ck . For practical reasons, we denote by [U ck ] a matrix which columns are vectors of R n generating a basis of space U ck : so [U ck ] is a rectangular matrix. Similarly, [F c ] is matrix of the unit vectors in each measured DOF (sensor). By definition of U ck , we have:

[Z] k [U ck ] = [F c ]. Let us recall that matrix [Z] k is singular. So column q of [U ck ], say [U ck ] q , is such that [Z] k [U ck ] q = [F c ] p (p th -column of [F c ]). It is obtained through the classical formula [U ck ] q := j =k φ ⊤ j [F c ] p ω 2 j -ω 2 k φ j ,
which gives [U ck ] q up to an element of the kernel of [Z] k , which is removed by

[U ck ] q := [U ck ] q -φ ⊤ k [M] [U ck ] q φ k .
Then we set Definition 3.1. The maximal dynamic visibility V ik max associated to i th -sub-structure and k th -mode is given by the ratio of energy induced by a perturbation on sensors DOFs of the i th -sub-structure, by the total residual energy

V ik max = max x∈R n -{0} x ⊤ [A]x x ⊤ [B]x , with [A] = ([U ck ]) ⊤ [K] i [U ck ] and [B] = ([U ck ]) ⊤ [K][U ck ].

Let us remark that the size of matrix [A][B]

-1 is the number of sensors, then can be quickly computed. The visibility value belongs to the interval [0, 1] and depends on sensors locations and on mode energy. A value of V ik max close to 0 ensure that a erroneous parameter is undetectable in the i thsub-structure, and can be neglected in the location error process. Conversely, V ik max close to 1 indicates that this sub-structure contributes significantly to the location process, for the error of constitutive law. Naturally, defining a frontier between "visible" and "unvisible" sub-structures is still an open problem, as it will appear in the following examples.

Remark 3.2. As in general matrix

[A][B] -1 is singular, the so-called minimal dynamic visibility V ik min := min x x ⊤ [A]x
x ⊤ [B]x = 0, which explains this quantity is a priori useless.

Modal Updating of JA-SON2

To illustrate our strategy of localization/correction, we take the example of large FEM of JASON2 (weight 525 kg and 3 meters high), made of plate and shell elements assembled by beam elements, with a total of 50392 nodes and 302352 DOFs, see Fig. 1.

PIM

Fig.1 FEM JASON2

The modal and FRF CRE methods are implemented in AESOP software [START_REF] Cogan | AESOP. Analytical-Experimental Structural Optimization Platform Version 5.0 User Guide[END_REF], which is developed by the LMARC Besançon. But very few numerical expenriences have been tested on the FRF CRE method. AESOP use SQP algorithm of Matlab to solve (4) or [START_REF] Roy | PRIMODAL An Innovative Software tool for Structural Dynamics in Industry User's Manual Version 1[END_REF] , and the NASTRAN DEMAP.

MAC coeff. = (π (φ i ).u k ) 2 π (φ i ) u k
The vibration tests have been performed by ALCATEL in low frequency and around 87 modes were measured from 5Hz to 150Hz. In collaboration with the CNES, a new "multi-dofs" RTMVI (Real Time Modal Vibration Identification) method have been recently developed by TOPMODAL (Toulouse -France) in the software PRIMODAL [START_REF] Roy | PRIMODAL An Innovative Software tool for Structural Dynamics in Industry User's Manual Version 1[END_REF]. Thanks to this method 29 experimental modes have been indentified, and the MAC have been calculated (Fig. 2). We only select in the sum of formula (4) a set of modes with significant modal effective mass. We choose parameter associated to a good visibility zone of the model, to be sure that the sensors permit to see its erroneous. For the mode 1, the first step localization indicate the Pshell 151000 (in NASTRAN Model) lateral panel of PIM cube. But the updating investigation with thickness parameter fails in spite of a 0.64 static visibility. Ploting the error location by elements permits the identification of the most affected CQUAD elements. We noted the modeling error zones are very small and close to the interfaces of all lateral panels of PIM cube. The highest error, and visibility beams localized are properties Pbar 20000 and 18000, corresponding to the links the PROTEUS platform to PIM cube. By a Monte Carlo analysis we tested the sensitivity of area moment of inertia and the section of these beams. We note that mode 1 grows asymptocally as the inertia of the Pbar 20000. Finally an increase nearly ten time the area moment of inertia, three time the section of the Pbar 20000, and 1.8 time the section Pbar 18000, minimizes (4), and improves the correlation for the two firsts paired modes (see Fig. 6 and Fig. 7). The complex geometry of these elements confirmes your choice (see Fig. 4).

PBAR 20 000 PBAR 18 000 Fig. 4 The complex geometry of pbar 20000 and 18000 The second step's localization shows that the error is concentrated around the AMR antenna (see Fig. 3). For the modes 2 to 7 the pbar 3050075 linking the AMR antenna to PIM are indicate, nevertheless they have low static visibility (0.2). Indeed the tickness updating of these beams do not improve the frequencies. To correct the the third mode (global vertical torsion) and the fourteenth mode (X-pumping) we take the thickness of Pshell 3050015 (AMR wire stretcher), it has the highest influence on the CRE, and 0.67 static visibility. A Monte Carlo run shows that the decrease of the thickness minimizes the CRE. We obtained a very good adjustment for the two firsts (Fig. 5) modes changing the Pshell 3050015 thickness from 0.001016m to 0.0008m. However, one's again the error is localized in a small zone around local stiffeners as we can see on Fig. 6. We performed the same updating parameters with CRE FRF method. By comparing the FRF before and after adjustment parameters, we cleary see the improvement around the firts modes of the distance with the experimental function. We illustrate on the following example the difficulty to select a good parameter on the large FEM of JASON2, and We simulate an experimental model with 175 DOFs (sensors). There is no explicit relation between the extended CRE and the output errors that can be evaluated directly on the basis of the paired analytical and experimental eigensolutions. We modify the thikness of Pshell 12000, Pshell 151000 and section of Pbar 20000. Therefore we compare the intial and modified models, using the location process. In order to valid the extended CRE updating process on JASON2, we have tried various strategies to solve the (4) problem. On Fig. 8, Pshell 12000 and 151000 have good static visibility (0.99, both), but only 0.32 for Pbar 20000. In a first time, we have tried to find directly (without sensitivity analysis) the minimum of (4) with the three parameters of pshell 12000, Pshell 151000, and Pbar 20000 : the SQP algorithm did not converge. In a second time, we chosen to update only Pshell 12000 and Pshell 151000. After 27'34 hours running, we only optained a good adjustment for the Pshell 151000 parameter. We applied a Monte Carlo Latin Hypercube sensitivity analysis to the Pshell 12000 parameter. After 22'24 hours, the third updating corrected the three parameters, see Fig. 9. But there exists situation where J k (u, v, w) decreases and distance between φ k and u k grows. To illustrate this phenomenon, we consider the numerical updating problem defined in Example 4.1. However if we choose the thikness of the pshell 13000 to perform the updating, the design objective function (4) decreases while the error between the eigendata increases !

Conclusion

The CRE method is costly in large F.E. model to adjust parameters. The localization-correction cycle have to be repeated until no further improment is possible, but we do not know the global minimum of (4). Visibility function permits to select zone with significant error energy. How to use the visibility dynamic function to estimate the quality of the updating ? How the scalar factors in (4) influence the localization of dominantly erroneous subdomains ? There is not explicit relation between the extended CRE and the eigendata, however the location and number of sensors can damage the condition number of the (4) minimizing problem.

We are now working on this line.

Fig. 3

 3 Fig.3 Static Visibility on 164 DOFs

Fig. 6

 6 Fig.6 Localization : pshell 3050015 stiffener on AMR

Fig. 7

 7 Fig. 7 JASON2 FRF comparison

Example 4 . 1 .

 41 The model updating is usually limited to a correction of modal data, by changing the most sensitive parameters.

Fig. 8

 8 Fig.8 Static Visibility on 175 DOFs Remark 4.2. The error estimator J k (u, v, w) defined in (2) grows when the term (πuu k ) t [K] π (πuu k ) increases.

  Numerical updating test on JASON2 with 175 DOFs.

		Data		Updating 2	Updating 3
	Property	Ini. Val.	Modif. Val.	Result 2	Monte Carlo Val.	Result 3
	Pshell 12000	0.0012	0.004	0.00061	0.0009	0.00121
	Pshell 151000	0.0012	0.003	0.00117	blocked	0.000119
	Pshell 20000	0.0001849	0.0004	blocked	blocked	0.00028
	CEE Ini. Cost : 0.38	CEE Cost 2 : 0.12	CEE Cost 3 : from 0.05 to 0.005
		Fig. 9


Nicolas Roussouly
email: nicolas.roussouly@icam.fr

Michel Salaun
email: michel.salaun@isae.fr

Fabrice Buffe
email: fabrice.buffe@cnes.fr

Anne Carpine
email: anne.carpine@thalesaleniaspace.fr

RELIABILITY APPROACH IN SPACECRAFT STRUCTURE

This paper presents an application of the probabilistic approach with reliability assessment on a spacecraft structure. The adopted strategy uses meta-modeling with first and second order polynomial functions. This method aims at minimizing computational time while giving relevant results. The first part focuses on computational tools employed in the strategy development.

The second part presents a spacecraft application. The purpose is to highlight benefits of the probabilistic approach compared with the current deterministic one. From examples of reliability assessment we show some advantages which could be found in industrial applications.

INTRODUCTION

Up to now, mechanical analysis in spacecraft structure has been deterministic. Numerical values of input variables are specified in order to guarantee structural performance in the worst cases. Structure validation is based on safety margin calculation which contains safety factors. This approach aims at simplifying the problem of uncertainties by makeing sure they are covered but is unable to bring under control neither risk of failure nor oversizing of the structure.

For several years, certain approaches to take into account uncertainties in modeling have arisen. The probabilistic approach is one of them and consists of considering uncertainties as probabilistic ones, i.e with ran-dom variables characterised by their probability density function (pdf). This procedure gives more information like output uncertainties, reliability with respect to a failure mode or sensitivity of input variables. These aspects are very helpful in industrial applications especially in risk analysis and in order to find an optimal design from an economical point of view that deals with technical and financial information.

To formalize the reliability concept we consider X = {X 1 , ..., X n } an input random vector with n random variables with a joint probability density function f X (x) and a performance function G(X) which defines three mechanical domains : (a) domain of success : G(X) > 0, (b) limit state : G(X) = 0 and (c) domain of failure : G(X) < 0. Then probability of failure reads

P f = ˆG(X)<0 f X (x)dx
and reliability, called F , is equal to F = 1 -P f . Therefore, the probability of failure calculation implies knowledge of the joint pdf and evaluation of the performance function. Joint pdf is generally assumed because of the lack of information on marginal laws and variable correlation. In large scale structures, performance function stems from a finite element code and the number of parameters to take into account grow significantly which is increasingly time consuming. Moreover, it is often necessary to determine reliability with respect to several performance functions because of the number of mechanical components and load cases.

Several method have been developed to handle these problems [START_REF] Lemaire | Fiabilité des structures[END_REF], [START_REF] Ditlevsen | Structural reliability methods[END_REF]. Although simulation methods such as Monte Carlo or improved techniques are robust, they are computationally too demanding for assessing low probability of failure, i.e lower than 10 -6 . Other methods, like first and second order reliability methods, commonly called FORM and SORM exist. They use an optimization algorithm to find the most probable failure point and then approximate the limit state (with a first or second order function) in order to assess probability of failure. Such methods are interesting with a regular limit state and a low probability of failure but require an optimization procedure per response. Another idea to achieve propagation of uncertainties is meta-modeling as an approximation of a finite element model. They are widely used in other domains of science and their efficiency in mechanical structure has already been proved. The most commonly used are first and second order response surface, polynomial chaos, kriging or radial basis function but there also exists more complicated kinds of meta-models such as neural networks or support vector machines [START_REF] Barton | Metamodeling : a state of the art review[END_REF]. In industrial cases, advantages can be found in meta-models. They can be treated with a relatively low computational cost by comparison with a finite element one. Moreover, as a result of lack of knowledge about uncertainty of variables, repeating a study with the lowest cost to assess importance of pdf is necessary. Finally, they may be useful in optimization and reliability based optimization.

In this paper we propose to use first or second order polynomial functions which are sufficiently well suited approximations in case of linear elastic behavior under static load with low variation on parameters. In the next section we present some computational tools used. The procedure is then applied to a spacecraft structure in order to obtain reliability results on components which can potentially fail at first. Since the finite element model is not updated after a first failure we are not able to take into account successive failures. Finally, the kind of results and benefits obtained will be discussed.

STRATEGY DEVELPOMENT

Sentivity analysis and variable selection

A finite element model of an industrial structure contains a lot of variables and it is often computationally too demanding to include all of them into the study. Generally, it appears that only a few influence studied responses. Therefore, detecting a couple of the most relevant is valuable and can be done through sensitivity analysis. Some methods exist in computer simulation whose purpose is to rank input variables in terms of their importance. These techniques are known as screening methods and are able to deal with hundreds variables (more information on screening methods can be found in [START_REF] Saltelli | Sensitivity analysis. Series in Probability and Statistics[END_REF]). One-at-a-time (OAT) design proposed by Morris [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF] appeared to be a good compromise between computational cost and the relevance of results [START_REF] Garnier | Etude des méthodes d'analyse de sensibilité, application aux structures spatiales[END_REF]. Morris OAT is a global sensitivity experiment because it covers the entire space over which variables vary. It is used to distinguish (a) negligible effects, (b) linear effects without interactions, and (c) non-linear effects and/or interactions. In practice, results are plotted on a graph which represents non-linear and/or interaction effects versus linear effects. Figure 1 shows an example of a Morris OAT experiment for the displacement at the end of a bending embedded beam in terms of force, young modulus, length, width and height. The number of computer runs is equal to 3(n + 1) where n is the number of variables (it could be more but 3(n + 1) is enough to obtain good information). Morris OAT doesn't really quantify the sensitivity of variables but rather gives a hierarchy. More than in probabilistic approach, this kind of method can help in model understanding, optimization and model updating. When it is used in variable selection, selected variables will be considered as random variables characterized by their distribution while others will be fixed to deterministic values. Therefore, meta-models will only link selected variables to studied responses. The meta-modeling process imply several steps which are presented now.

Specification of experimental design

Experimental design specifies a set of input values for which finite element computation is performed in order to get corresponding outputs. Well known techniques called Design of Experiment (DOE) are generally used in physical experiments. They focus on planning the experiment in order to satisfy some optimality criteria based on physical random error. In our study, the experiment is a computer deterministic simulation and some distinctions such as random error or number of factors and their range of variation can be noticed. In this way, many authors advise the use of "space filling" design in order to best cover the design space. Methods like latin hypercube sampling (LHS) or low discrepancy sequences which imply uniform distribution of points in space are often used in literature [START_REF] Wang | Review of Metamodeling Techniques in Support of Engineering Design Optimization[END_REF][START_REF] Simpson | Metamodels for Computer-based Engineering Design: Survey and recommendations[END_REF][START_REF] Saltelli | Sensitivity analysis. Series in Probability and Statistics[END_REF][START_REF] Sacks | Design and analysis of computer experiments[END_REF]. We will afterward perform LHS technique.

Determination of meta-models

In this paper we use first and second order polynomial functions as an approximation of the finite element model in linear regression. The linear regression equation is given by

Y = Xβ + ε (1)
where

Y
QR(X β) = V ar[X β] -Bias[X β] 2 (2)
If we suppose that for all p, ε p is an independently, normally distributed random variable of zero mean and with a variance noted σ 2 , then QR can be estimated by the well known Mallows Cp criteria. Other criteria exist based on other formulations than QR to check model quality. The most well-known are Akaike Information Criteria (AIC) and Bayesian Information Criteria (BIC) (see more details in [START_REF] Azaïs | Le modèle linéaire par l'exemple[END_REF][START_REF] Cornillon | Régression, Théorie et applications[END_REF]). These criteria can be used to compute the best model from several potential regressors. This is known as stepwise regression. Another interesting measure used in stepwise regression is the R-square coefficient. It represents the fraction of variation of output explained by the model and therefore indicates how well the model reproduces output. R-square coefficient varies between 0 and 1, and is equal to 1 if the model perfectly fits the sample. Actually this coefficient increases in terms of the number of regressors and is equal to 1 when the number of regressors is equal to the number of samples. We will prefer to use the adjusted R-square, noted adjRsquare, which is the R-squared penalized by the number of terms [START_REF] Cornillon | Régression, Théorie et applications[END_REF].

Mallows Cp, AIC and BIC are based on the normal distribution of error ε p in linear regression which is a classical assumption in statistics to approximate physical experiments. Although in our case the experiment is deterministic and terms of error ε p only specify the meta-model lack of fit, these tools can be used as a guideline but is unable to conclude about model validity. The validation step is done without assumptions.

Selection and validation of meta-models

Previous criteria used in stepwise regression are more or less restrictive and result in different models. Selection of the best model must be done with respect to its ability to predict responses from a new input sample. The best way would be to have two samples, one for learning step and another for selection but this process imply a larger computational time. A possible approach is the use of cross-validation techniques. This method consists of the following steps : (a) split the global sample into D samples, (b) leave one sample out (called the validation sample) and build the model with the others, (c) use the validation sample to assess the approximation error, (d) repeat previous steps with all samples. Then, prediction error is estimated as the mean of errors computed on all validation samples. The measure of error used is the classical root mean square error (RMSE). Therefore among all models proposed by stepwise regression, the one with the lowest prediction RMSE is chosen.

To assess model accuracy in the region of interest we also compute the error MAX on validation samples, given by max p∈ [1,P ] ( Y p -Ŷp ) which estimates a more local error than RMSE. The validation of the model is guaranteed if predicted RMSE and predicted MAX are lower than a threshold. In our case, we use as the threshold 5% of maximum variation of the response. This ensures that error due to the meta-model is low with respect to response variation. From a graphical point of view, a good way to assess model accuracy is to plot response of the finite element model versus response of the meta-model with the domain of admissible error (cf. Figure 2). Once the meta-model is built and validated, it is used as a surrogate of the finite element model. It is then employed in probabilistic procedures in order to analyse output uncertainty and to assess reliability information.

Methodology in application

We present here the manner in which the previously described tools are used in application.

The morris OAT method is performed on the finite element model but no variable selection is done. Therefore, we decided to build linear meta-models with all variables from a Latin Hypercube sampling with 3(n + 1) points which are necessary for the validation step (n is the number of variables). For each response, three linear models are established through stepwise regression with Mallows Cp, BIC and adjR-square criteria (Mallows Cp and AIC are equivalent in linear regression [START_REF] Stone | An asymptotic equivalence of choice of model by cross-validation and Akaike's criterion[END_REF]). The best one is chosen from predicted RMSE computed with the cross-validation technique. If predicted MAX error and predicted RMSE are lower than admissible error (5% of maximum variation of response), the linear model is validated. If not, we add quadratic and cross-product regressors of the most significant variables selected with Morris OAT from the value indicating non-linear and/or interaction effects. Since stepwise regression has removed useless linear regressors, added quadratic and cross-product terms do not involve more sampling points. This procedure avoids performing more finite element computation than necessary. It allows a sufficiently accurate model to be obtained in the following application. In other cases, if more regressors are useful, we must add sampling points and thus perform more finite element computation. Finally, if second order polynomial functions are not good enough, more sophisticated metamodels have to be used.

APPLICATION ON A SPACECRAFT

STRUCTURE

The following example involves the TARANIS spacecraft whose platform belongs to the Myriade family developped by CNES. The analysis has been treated under 8 quasi-static load cases which are those of the qualification step.

Finite element model description

The mechanical model is deterministic and has been done with the finite element code MSC NASTRAN. Figure 3 presents the finite element model. It contains about 380 000 degrees of freedom. Loads consist of acceleration of -9.5g in longitudinal direction (X direction) and 5.2g in lateral direction (Y and Z direction). Eight cases are different projection of lateral load on Y and Z axis, they are presented in Table 1. We have to notice that the model has been optimized and resultant values are not those provided by CNES. The methodology described in section 1.5 is applied on the spacecraft. Latin Hypercube sampling and Morris OAT are performed on ±20% of variation on each variable. These two methods each involve 252 finite element computations (see time computation in Table 3).

Uncertainty and sensitivity analysis

Uncertainty of responses can be analysed from a Monte Carlo sampling. It gives statistical measures of response such as its expectation and standard deviation (and higher order moment) but a probability density function can also be determined thanks to a statistical hypothesis test. Another interesting analysis consists of variance decomposition of responses (see more details in [START_REF] Saltelli | Sensitivity analysis. Series in Probability and Statistics[END_REF]). This procedure aims at identifying the fraction of the total variance of a response which is due to any individual input variables. Figure 5 shows the variance decomposition of Von Mises stress in +Z panel (i.e panel whose normal direction is +Z).

Reliability analysis

Reliability assessment implies the definition of a performance function for each response. By comparison with a deterministic approach, failure modes considered are the elastic limit for panels and the sliding limit for interface screws. These two limits involve new random variables (elastic limit, friction coefficient, preload) which are assumed independent (with a normal distribution and a coefficient of variation equal to 10%) for each component and with respect to those already defined in the finite element model. This can be a strong assumption especially between the elastic limit and the young modulus of a material which could be correlated.

As the mechanical model is a system, components depend on each others and only potential first failures can be determined. If one component fails, the finite element model does not enable to know what happen afterward and it is not able to take into account successive failures. The probabilities of failure are assessed through a Monte Carlo analysis with 10 7 samples (see time computation in Table 3).

Probability of failure results

Results are given in Table 2. It only contains the worst load case, i.e the case which involves the largest probability of failure. For this case, components which are able to fail are presented with their probability of failure and their percentage of first failure. 1.9 10 -3 1.87 10 Results show that +Z panel under case 400 is the most critical element. The probability that Von Mises stress exceeds the elastic limit of material is 1.9 10 -3 . Figure 6 shows the probability of failure of +Z panel with the distribution of performance function. For this kind of component and failure mode, such a high value could be problematic. This kind of problem would be missed in a deterministic approach since calculated safety margin is equal to 0.23. Another interesting result is the percentage of first failure which indicates how many times the component fails in a first instance. This information could be helpful in envisaging construction of system failure scenario.

Some examples of the benefits

We present here two examples which emphasize the benefits of a reliability approach compared to a deterministic one.

The first one is an illustration of safety margin calculation versus reliability assessment. We consider here the sliding limit of two interface screws between -Z panel and -X panel under load case 700. We assume that the friction coefficient and pre-load follow a normal distribution with 0.23 and 17260N respectively for the mean and 10% for both coefficient of variation [START_REF] Carpine | Faisabilité d'un dimensionnement basé sur une approche stochastique[END_REF]. In a deterministic approach, the safety margin is calculated with a B-value for friction coefficient and with a value lower than a B-value for pre-load (all other variables are considered at their mean). Moreover a safety coefficient of 1.25 is considered. From this configuration the sliding margin is negative and is equal to -0.022 for one screw and -0.017 for the other. This can create some problems in a qualification step. With a reliability approach random variables are considered, characterized by the previous distributions. Moreover, the safety coefficient is removed. In this case, probabilities that screws slide are 7.2 10 -4 and 6.5 10 -4 respectively. This kind of probability of failure could be considered as acceptable in this context.

The second example illustrates the importance of reliability assessment in order to compare two manufacturing processes of aluminium panel skins. Skin thicknesses are assumed to follow a uniform distribution with a mean of 0.6 mm. With such a low thickness, two manufacturing processes are considered, one expensive with good accuracy which implies a thickness variation of ±5% ; one cheaper with a lower accuracy which implies a thickness variation of ±20%. Probabilities that Von Mises stresses exceed elastic resistances (under case 400) are 1.9 10 -3 and 3.47 10 -4 respectively. Considering manufacturing costs of both processes, a well-founded decision can be taken which deals with the technical risk and financial aspect [START_REF] Buffe | Aide à la prise de décision[END_REF].

CONCLUSION

This paper presents two main issues. On one hand, we attempted to show the relevance of results which arise from a probabistic approach and more precisely in reliability assessment. Although current deterministic approach is firmly fixed in structural analysis, it doesn't enable to improve designs with respect to robustness or reliability linked with economical constraints whihc are always more and more restrictive. A probabilistic approach is one of possible means and is starting to be applied in several industrial domains. Of course, current practices will not be replaced but it is a way to enrich them with more accurate information especially in order to help in decision. Making this kind of ap-

Step Time Morris OAT 9h30 Latin Hypercube Sampling 9h30 Meta-models construction 20 min 10 7 samples on meta-models 1h

Table 3: Time computation of each step in methodology. Desktop computer used is a dual core 2.8 GHz, 1.97 Go RAM proach can bring benefits to the definition of less conservative sizing rules. On the other hand, we presented a resolution strategy based on meta-modeling. It can be employed on quite large structures and is not too computationally time demanding. Table 3 resumes the computational time of each step of methodology. Even if it could appear large relatively to a deterministic approach, it is very fast when compared with a direct Monte Carlo on finite element model which would last 43 years. Finally advantages of meta-models are numerous. Structural behaviour can be represented well even with simple models : linear polynomial functions are often sufficient. Although in our case the structure is linear elastic under quasi-static load with low variations, more sophisticated meta-models could deal with more complicated mechanical behaviour. Meta-models are also re-usable in order to measure the influence of stochastic chosen input parameters. This seems to be an essential point as distribution of random variables are usually assumed because of the lack of databases. Another interest of meta-models can be found in optimization or more precisely reliability based optimization.

Figure 1 :

 1 Figure 1: Graphical result of Morris OAT experiment

Figure 2 :

 2 Figure 2: Meta-model accuracy

Figure 3 : 2 . 2

 322 Figure 3: TARANIS finite element model

Figure 4 :

 4 Figure 4: Bottom panel of TARANIS

Figure 5 :

 5 Figure 5: Variance decomposition on the +Z panel stress response

Figure 6 :

 6 Figure 6: Representation of performance function distribution with reliability result on +Z panel

Table 1 :

 1 Load cases

	Acceleration on axis (g)
	Cases	X	Y	Z
	100	-9.75 5.20	0.00
	200	-9.75 3.68	3.68
	300	-9.75 0.00	5.20
	400	-9.75 -3.68 3.68
	500	-9.75 -5.20 0.00
	600	-9.75 -3.68 -3.68
	700	-9.75 0.00 -5.20
	800	-9.75 3.68 -3.68

Table 2 :

 2 Reliabilty results

		98.94
	-5	0.96
	< 10 -6	0.05
	< 10 -6	0.03
	< 10 -6	0.02