Monsieur Franck

Monsieur Frédéric Druesne

Laurent Dubar

Frédéric Collègues Alexandre

Kossiga

Ruqing, Gorka, Hans Anna Anthony

Xavier, Sylvain, Mariem, Rim, Oussama Ivan

Keywords:

Tout d'abord, je remercie particulièrement le Professeur Thierry Tison, Directeur de thèse qui m'a beaucoup aidé dans le travail et la vie pour me donner les pistes de recherche. Je vous remercie également le Professeur Hakim Naceur, Co-directeur de thèse qui

Context and objectives of thesis

Additive Manufacturing (AM), also commonly known as 3D printing, is a manufacturing process for building almost any physical solid part from its three-dimensional digital model. Nowadays, AM technology has a very fast growing market capability mostly in the aerospace, automotive, healthcare, and consumer products sectors to build 3D complex parts. In this context, the products are built-up by direct deposition of successive controlled thin layers of melted metal on top of each other considering different materials, such as for example titanium, carbon fiber or fiberglass. Thus, this technology facilitates the realization of complex lightweight structures, such as lattice structures which are interesting for energy absorption problematics. These last ones are defined with the aspect ratio, the unit-cell size and shape and can have dimensions up to nanometer.

Although the manufacturing process is more and more well managed with for example the Electron beam melting (EBM) process, some geometric and material manufacturing variability are detected. They are generally due to the built-part orientation, the layer thickness and the melted pool depth. The manufactured structure presents high variability in strut size and roughness and so it is necessary to define an effective strut diameter in comparison to the CAD diameter. In the current work, we focus our discussion on the case of the linear buckling problem under compression loading by considering both experimental and numerical ways. To improve the predictivity of numerical simulations for lattice structures, it is necessary to integrate the observed variability in the solving algorithm by performing specific uncertainty propagation methods.

Indeed, to take into account these uncertainties and tend to reliable and robust predictions, it is necessary to generate multiple numerical simulations by performing sensitivity analyses, designs of experiments, non-deterministic studies or even reliable and robust optimizations. The objective is to simulate the evolution of mechanical responses (uncertainty propagation step) as a function of input parameter variations (uncertainty modeling step) and to detect failures and performance reductions of products (uncertainty management step). Nevertheless, whatever the theory used to model uncertainty, it is important to maintain the computational time as low as possible and define efficient numerical strategies with the help of surrogate or reduced order models.

In this context, the objective of the thesis is to replace the original buckling problem by a suitable approximation or to directly reduce the size of the original problem. Although many works have already been performed in this domain, the efficiency of the approximations is very dependent of the mathematical nature of the studied problem through the definition of input parameters and evolution of output solutions. Thus, the application of the approximation techniques in uncertainty propagation methods for studying the uncertain critical buckling loads of lattice structures request specific research investigations. Accordingly, the aim of this thesis is to integrate material and geometric uncertainty, experimentally quantified, in uncertain buckling analysis of lattice structures, produced by additive manufacturing.

More precisely, this thesis investigates different approximation techniques compatible with both uncertainty propagation methods and buckling analysis. A state of the art allows first to classify the different approximation techniques. It is followed, for the first time, by a comparison with a same general numerical framework to address an assessment of advantages and drawbacks and finally supply the more efficient numerical method. Next, in a second part, lattice structures and additive manufacturing are studied experimentally and numerically in a way to perform correlation considering uncertainty inherent to the process.

Organization of this report

The report is organized in two parts, respectively defined with five and two chapters. The first part is dedicated to the evaluation and the coupling of metamodels and reduced order models by considering a same numerical framework. The aim is first, to analyze the role of key parameters for each method and second, to quantify their effect on the quality of approximation of output variables. The current simulations concern the linear buckling problem and in particular the calculation of the buckling load and associated mode. A numerical strategy coupling a metamodel and a reduced order model is then proposed and so tested.

The first chapter summarizes the main steps relative to the numerical solving of a deterministic linear buckling problem and presents the studied structure and the numerical benchmark.

In this context, a multiparametric analysis is performed considering modifications of material and geometric properties. In the second, third and fourth chapters, different numerical methods (metamodels, reduced order models, projection) based on specific samplings and different techniques (regression, correlation, decomposition, perturbation, derivatives) are first compared.

Second, new numerical methods coupling metamodels and reduced order models are proposed and compared in the fifth chapter.

In the second part of this report, the integration of the proposed approximation in uncertainty propagation method is discussed. The objective is to study the evolution of uncertain buckling loads as a function of uncertainty, identified on material and geometric parameters, for two lattice structures. Uncertainty is mainly due to the additive manufacturing process and is modeled either by the fuzzy sets theory or by the probabilistic theory. A correlation between uncertain numerical data and experimental results is finally proposed.

The sixth chapter is dedicated to the presentation of the two studied lattice structures, defined Introduction by additive manufacturing, and the quantification of the variability by experimentations. Moreover, quasi static tests are performed to identify the experimental buckling loads of each lattice structure. Next, finite element models of lattice structures are proposed and some deterministic buckling simulations are compared to experimental results. The seventh chapter discusses of uncertainty propagation in numerical simulations by considering two uncertain modeling theories, the probabilistic and the fuzzy sets ones. The proposed approximation is successively integrated in Monte Carlo Simulations (MCS), Zadeh's Extension Principle (ZEP) and optimization technique. Uncertain numerical solutions can then be compared to experimental results in term of buckling loads.

Finally, a general conclusion closes this report and provides some perspectives for future works.

2.1

Elastic stability based on the potential energy . . . . . . . . . . . . . . .

2.2

Formulation of the eigenvalue problem and linearization . . . . . . . . . Passing the critical level, the structure undergoes a significant and a prompt deformation which allows it to resist again to the applied loads, but in a totally different geometric configuration than the original undeformed shape.

Instabilities can be classified into two main types: bifurcations and limit points (or snap through). The so-called Euler beams considered straight behave as bifurcations. These are characterized by the appearance of a secondary branch on the load-displacement curve, which is discontinuous with respect to the primary branch as indicated in Fig 1 .2. In the case of Figure 1.2a, where the axial load P is a function of the lateral displacement d of the free end of the column, we note that, if we consider the column perfectly straight (without imperfection), then there will be no transverse displacement before buckling (theoretical blue curve). In this case, we talk about buckling without initial displacements also called EULER buckling.

In the case of Fig 1 .2b, blistering buckling of the structure occurs after a non-negligible displacement. In this case, we talk about buckling with initial displacements also called DUPUIS Buckling. It is therefore clear that the stability of a structure must be studied taking into account geometric nonlinearities. But it is also possible to estimate the critical loads by performing a simplified study. This estimate will be more or less precise depending on the type of instability and the assumptions used in the simplified model. We talk here about eigenvalue analysis, which consists in solving a problem with eigenvalues around a point of equilibrium of the structure at a given load level.

It should be noted that a stability study via eigenvalue calculation is generally limited to the elastic case. We will limit the developments to this case in this thesis.

2 Linear buckling analysis

Elastic stability based on the potential energy

Suppose we can estimate the potential energy Π of a structure, as a function of the displacement field U. Let's develop the functional Π(U) in Taylor series around solution state U 0 , gives

Π(U 0 + δU) = Π(U 0 ) + δΠ(U 0 , δU) + δ 2 Π(U 0 , δU) + .... (1.1)
where δU is a perturbation of the displacement field around the point of equilibrium U 0 . The term δΠ(U 0 , δU) is zero since the structure is under equilibrium, and if we neglect higher terms, we obtain: This equation can be interpreted as follows. If ∆Π is greater than zero, then the potential energy for the displacement field U 0 + δU will be higher than the calculated for the field U 0 .

∆Π = Π(U 0 + δU) -Π(U 0 ) ≈ δ 2 Π(U 0 , δU) (1.
Hence, the structure will be stable, because in order to go from the state U 0 to the state U 0 + δU it is necessary to supply energy (Fig 1.3a).

Contrarily, if the potential energy for the field of displacement U 0 + δU is less than the one calculated for the field U 0 , then ∆Π will be smaller than zero. So the structure will have tendency to restore energy and therefore it will be unstable (Fig 1.3b). In the light of the previous remarks, it is clear that the limit of moving from the stable state to the unstable one corresponds to the case ∆Π = 0 as indicated in figure 1.3c.

Thus, from the second variation of the potential energy, it is possible to characterize the majority of the instabilities by using equation 1.2 (δ 2 Π = 0).

Formulation of the eigenvalue problem and linearization

Let consider the classical form of the internal elastic strain energy of a structure stated as

Π int = 1 2 V ε : σdV (1.3)
where V is the volume of the undeformed structure, σ is the internal symmetric stress tensor and ε is the Green Lagrange symmetric strain tensor, which can be split into a linear part ε l = ε l (U) and a quadratic one ε nl = ε nl (U 2 ), and hence we can write ε = ε l + ε nl .

The second variation of the internal strain energy functional Eq. 1.3 leads to:

δ 2 Π int = V δε l : δσ + δε nl : δσ + δ 2 ε nl : σ dV (1.4)
In Eq. 1.4, δε nl = δε nl (U, δU) is a linear function of the displacements, while δ 2 ε nl (δU 2 ) is constant. The assumption of a linear elastic behavior of the structure can be expressed by introducing the elastic isotropic material tensor H with:

σ = H ε l (1.5)
Finally the expression of the second variation of the internal strain energy functional takes the following form:

δ 2 Π int = V δε l : (H δε l ) + δε l : (H δε nl ) + δε nl : (H δε l ) + δ 2 ε nl : σ dV (1.6)
It has to be noticed that in Eq. 1.6, the contribution of δε nl : (H δε nl ) has been neglected, since we assume that strains induced in the structure during buckling will remain small.

Let consider now a structure under a compressive load F 0 at its equilibrium state given by a displacement vector U 0 and a stress state σ 0 . If we assume that the perturbed unstable solution can be expressed as U = λU 0 , with λ a constant real factor, then we can have for the perturbed state:

δε nl = λ δε nl0 σ = λ σ 0 (1.7)
where δε nl0 = δε nl0 (U 0 , δU) and σ 0 is the stress tensor at the initial stable state before buckling. Substituting now Eq. 1.7 into Eq. 1.6, one can obtain:

δ 2 Π int = V δε l : (H δε l ) +
λ δε l : (H δε nl0 ) + δε nl0 : (H δε l ) + δ 2 ε nl : σ 0 dV (1.8)

For structures under pressure loading p which follows the normals to the deformed structure, it is important to consider the variation of the external energy which can be stated as:

δΠ ext = A p n • δUdA (1.9)
where n = n(U) are the normals to the deformed surface A. The second variation of the external energy functional Eq. 1.9 leads to:

δ 2 Π ext = A p δn • δUdA = A λ p δn 0 • δUdA (1.10)
with n 0 is the normal to the initial state surface before buckling.

By considering the critical equilibrium state of a structure, we can write

δ 2 Π = δ 2 Π int -δ 2 Π ext = 0 ∀ δU = 0 (1.11)
Finally the variational form of equilibrium is stated as:

V δε l : (H δε l ) dV = -λ V δε l : (H δε nl0 ) + δε nl0 : (H δε l ) + δ 2 ε nl : σ 0 dV -A p δn 0 • δUdA ∀ δU = 0 (1.12)

Eq. 1.12 represents a linear eigenvalue problem, for which the first set of solution (λ 1 , U 1 ) will represent respectively, the critical load factor (as F cr = λ 1 F 0 ) and the associated buckling mode shape of the structure.

Finite element discretization

Generally the displacement-based formulation of the Finite Element (FE) method when applied for the discretization of a structure using any kind of finite element model such as 1D rod, 2D plate, or 3D brick, will result in the following general approximation form of the displacement vector U:

U =      U X U Y U Z      = nne i=1      N i U Xi N i U Y i N i U Zi      (1.13)
with U X , U Y , U Z the three cartesian components of the displacement vector U, N i is the shape function at node i, U Xi , U Y i , U Zi are the nodal displacement components at node i and nne number of nodes per element.

Eq. 1.13 can be rewritten in a condensed form as follows:

U = N U n and δU = N δU n (1.14)
where N 3×3nne the FE approximation matrix of size ndim × nne × ndim, with ndim number of dimension (1, 2 or 3) and U n the nodal displacement vector of dimension nne × ndim.

Once the displacement field is approximated, the next step is to express the Green Lagrange strain tensor ε in terms of nodal displacements vector U n , this is stated as

ε = ε l + ε nl = 1 2 ∇U T + ∇U + 1 2 ∇U T ∇U (1.15)
where ∇U is the displacement gradient matrix. By using the engineering notation, one can derive the expression of the linear strain vector as follows

ε l = {U X,X , U Y,Y , U Z,Z , U X,Y + U Y,X , U X,Z + U Z,X , U Y,Z + U Z,Y } T (1.16)
The nonlinear part of the strain vector is

ε nl =                      0.5 U 2 X,X + U 2 Y,X + U 2 Z,X 0.5 U 2 X,Y + U 2 Y,Y + U 2 Z,Y 0.5 U 2 X,Z + U 2 Y,Z + U 2 Z,Z U X,X U X,Y + U Y,X U Y,Y + U Z,X U Z,Y U X,X U X,Z + U Y,X U Y,Z + U Z,X U Z,Z U X,Y U X,Z + U Y,Y U Y,Z + U Z,Y U Z,Z                     
(1.17)

Based on the expressions of linear and nonlinear strains, their variations δε l and δε nl can be written in a condensed matrix form, as follows

δε l = B l δU n δε nl = B nl δU n or δε nl0 = B nl0 δU n (1.18)
with B l the constant strain approximation matrix, B nl = B nl (U) or B nl0 = B nl (U 0 ) is the nonlinear strain approximation matrix which is function of the initial displacement.

The contracted tensor product δ 2 ε nl : σ 0 can be rewritten as a function of the variation of nodal displacement vector δU n as follows:

δ 2 ε nl : σ 0 = δU T n B T g S 0 B g δU n (1.19)
where S 0 = S 0 (σ 0 ) is an initial stress matrix, and B g a geometric constant matrix. The last term concerning the normals first variation δn 0 with the respect to the displacements, which can be approximated using the FE, as:

δn 0 = G δU n (1.20)
Finally by using Eqs.1.18-1.20, then the linear eigenvalue problem Eq. 1.12 reads:

δU T n V B T l H B l dV δU n = -δU T n λ V B T l H B nl0 + B T nl0 H B l + B T g S 0 B g dV - A p G T NdA δU n (1.21)
which can be written in a condensed form:

K l δU n = -λ K u + K σ -K λ δU n (1.22)
where K l , K u , K σ and K λ are respectively the linear stiffness matrix, the initial displacement matrix, the geometric matrix and the follower force matrix. These matrices are defined by:

K l = V B T l H B l dV K u = V B T l H B nl0 + B T nl0 H B l dV K σ = V B T g S 0 B g dV K λ = A p G T NdA
(1.23)

In the case of constant compressive loads (no follower forces such as pressure), K λ is zero.

Also, in the case of bifurcation-type buckling problems, the initial displacement matrix K u may be neglected, Eq. 1.22 simplifies to the so-called EULER buckling problem, stated as:

K l δU n = -λ K σ δU n (1.24)
For the remaining manuscript we will focus on the EULER buckling and therefore Eq. 1.24 will constitute our state problem. Three types of FE models have been developed during our

Computational procedure

The general procedure for buckling simulation of a structure under compressive load is as described in Figure 1.5.

1. Start by meshing the structure.

2. Assign physical properties, apply boundary conditions and the compressive load F 0 .

3. Proceed with the FE assembly to build the global linear stiffness matrix K l and the global external force vector F 0 .

4. Solve the linear static system K l U 0 = F 0 .

5. Proceed with the FE assembly to build the global matrices K σ , K u and K λ .

6. Solve the eigenvalue problem :

K l z = -λ K u + K σ -K λ z
and determine the solution λ cr = λ min and the corresponding eigenmode z min .

7. Post-process the critical buckling load F cr = λ cr F 0 .

8. Post-process the critical buckling mode shape using displacement vector U cr = z.

Figure 1.5 -Computational procedure for linear buckling analysis.

In the following chapters, we will consider the vector s which aggregates both the eigenvector z and the buckling load factor λ cr , such as:

s = z λ cr (1.25)
3 Multiparametric analysis

Numerical model description

In the sequel we will present the academic numerical application, which will serve as a benchmark for validation of different methods. The application consist in a two-bar truss under compression as described in Fig 1 .6. The two-bar structure is pined at both ends whereas an external vertical compressive load F is applied at the top-center of the structure.

The geometric and material parameters of the two-bar structure are summarized in Table 1.1.

In order to evaluate the efficiency of different variability techniques, the two-bar structure is modeled using a 3D mesh based on brick elements as shown in Fig 1 .7. The resulting FE mesh corresponds to 4 × 4 × 40 bricks for each side which results in a total of 1280 bricks.

For simplification reasons, a boundary condition U Z = 0 is imposed on all nodes in order to reduce the solution to a 2D plane case problem. Table 1.1 -Geometric and material characteristics. Figure 1.9 -Two-bars under compression: First buckling mode.

Response surface definition

In order to validate the performances of newly developed alternative mathematical models in the present thesis, variations on four geometric and material parameters have been studied.

Namely the Young's modulus E, the inclination angle α and the two dimensions of the cross section (height h and width b). The limits of variations of the four parameters are summarized in To study the effects of each parameter on the critical buckling load, a full factorial DOE was performed considering seven levels of responses per factor, which led to 2401 samples. The aim here is to generate a significant sets of samples in order to discuss about the associated computational time of numerical methods. Two response surfaces representing the critical load were generated and studied. The first one was defined as a function of Young's modulus E and First, we observe that the dimensions section b and h have larger effects on the buckling load factor compared to the effects of the Young's modulus and the angle. Indeed, the maximal variation on output data for the first response surface is close to 25% whereas, for the second response surface, almost 40% is observed for only 10% of variations imposed on input parameters.

Second, we naturally detect a linear behavior of critical load as a function of the Young's modulus.

For the three other parameters, the behavior is rather quadratic.

As the computational time relative to the evaluation of one sample of parameters is close to 3s, the total time for the multiparametric analysis is close to 2h for a finite element model with only 6075 dof (or 4030 dof active). This application shows clearly the necessity to develop an alternative mathematical model to facilitate the integration of multiparametric analysis in design step and consider larger finite element models. Chapter 2

Metamodels and reduced order models

Introduction

Metamodels and Reduced Order Models (ROM) have largely been integrated these last years in advanced simulations, such as the Design Of Experiments (DOE), the parametric optimization or the uncertainty propagation. The main aim is to reduce the computational cost of simulations, which rely on a multiparametric analysis. Indeed, the original finite element problem must be solved several times for different values of input parameters. Thus, the integration of metamodels and reduced order models allow here to replace the original problem by a suitable approximation or to directly reduce the size of the original problem. To achieve this objective, many methods and algorithms have been developed and applied in multiple scientific domains.

Obviously, it is rather difficult to propose a full overview considering the large number of publications. However, the main works about the metamodels discussed on methods based on regression and/or correlations models, such as polynomial approximations or Response Surface Method (RSM) [START_REF] Myers | Response Surface Methodology: Process and Product Optimization Using Designed Experiments[END_REF], Radial Basis Functions (RBF) [Baxter, 1992] or Kriging [START_REF] Sakata | Structural optimization using kriging approximation[END_REF], Wei et al., 2004, Kaymaz, 2005, Wang et al., 2013]. The idea is to build functions alternative to classical numerical simulations from only a set of specific samples of a DOE. Neural networks [Haykin, 1999] have been introduced in this context too. The main problem, relative to metamodels approach, is to identify the necessary number of samples to build the alternative models with accuracy. Indeed, if only few samples are available, the prediction capability of resulting approximated models would be insufficient. On the contrary, ill-conditioning and expensive cost of additional computations can be observed with large samples.

In parallel, reduced order models, such as Proper Orthogonal Decomposition (POD) [Ryckelynck, 2005, Yvonnet et al., 2007, David et al., 2012], Proper Generalized Decomposition (PGD) [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF], Nouy, 2010, Ladevèze and Chamoin, 2011, David et al., 2012] or more recently Hyperreduction [Ryckelynck, 2005, Rutzmoser and[START_REF] Rutzmoser | [END_REF], have been investigated. The general idea is to approximate a higher dimension system by another one of much lower dimension. Indeed, PGD method has been introduced by Ladevèze [Ladevèze, 1999] and called radial approximation in the Large Time INcrement (LATIN) method framework. Contrary to POD, the PGD is based on a separated representation of the variables and takes directly the nature of the problem into account. In the case of POD, an optimal basis is built considering proper orthogonal modes extracted from a correlation matrix, composed of different snapshots of numerical simulations or experiments.

This chapter is dedicated to the evaluation of the performance of quadratic regression, RBF, Kriging metamodels and ROM-POD considering a same numerical framework. The objective here is to test these methods in term of precision and computational time for a multi-parametric linear buckling problem with material and geometric parameters variations. Finally, considering the advantages and drawbacks of these methods, a new method, coupling ROM-POD and Kriging, will be proposed.

A metamodel is an approximation of a system response build from data available for limited number of selected snapshots or samples. These last years, metamodels have been considered

as valuable tools for investigating multiparametric analysis for a wide range of engineering applications. In this section, the main equations governing the definition of quadratic polynomial model, Radial Basis Function, Kriging and Proper Orthogonal Decomposition are resumed. The objective is to highlight the key parameters associated to each technique, which could affect the quality of the approximation and the associated computational time.

Let's consider n s snapshots to define the metamodels. a buckling problem is solved using the procedure described in Section 2.4 of Chapter 1 and the output eigensolution vector s(p i ) (aggregation of eigenvalue and eigenvector) is calculated. The j th studied component of this vector is noted s j (p i ) with j = 1...(n dof + 1) whereas the vector s j and the matrix S respectively contains all the output results of the j th studied component and the output results for all the components and all the snapshots.

All variables described above are used to build the metamodel in a phase which is generally called the offline step. In the following, the vector noted p represents a set of parameters for a new evaluation of the problem and then correspond to the online step.

Quadratic model

The first approximation classically used in engineering to characterize the behavior of a mechanical solution as a function of input parameters is based on a polynomial model with or without cross terms. Indeed, this model is easy to implement and only requires the calculation of unknown coefficients associated to each order. For the present applications, a quadratic model with a complete basis is considered. This choice is also in adequation with the evolution of buckling loads reported in Section 3.2 of Chapter 1.

The general form of the j th studied component of output vector s j (p) for a parameter vector p is :

s j (p) = a 0 + np k=1 a k pk + np k=1 np l=1 a kl pk pl (2.1)
where a 0 , a k and a kl are the unknown coefficients to be determined by using a least squares

Correlation functions R(θ, d) Linear max{0, 1 -θ|d|} Exponential exp(-θ|d|) Gaussian exp(-θd 2 )
Table 2.1 -Classical correlation functions.

fitting technique. Thus, the minimal number of snapshots must be equal to the number of unknown coefficients.

Radial Basis Functions

Radial Basis Function (RBF), originally proposed by Hardy [Hardy, 1971], has been developed

for scattered multivariate data interpolation [Buhmann, 2004]. This method relies on linear combinations of a radially symmetric function based on specific metric to approximate response functions.

The RBF regression model relative to the approximation of a component of eigensolution vector s j (p) for a parameter vector p is expressed as follows:

s j (p) = ns i=1 w i R(θ, d) (2.2)
where d = ppi is the radial distance from p to pi . w i is the weighted factor of the i th basis function.

The correlation function R(θ, d) can be given by different forms as shown in Table 2.1. The variable θ is a user defined scaling parameter. In the present application, it is fixed to the average distance of studied samples.

The unknown interpolation coefficients w i can be computed by solving the following system of linear equations:

      R 11 R 12 • • • R 1ns R 21 R 22 • • • R 2ns . . . . . . . . . . . . R ns1 R ns2 • • • R nsns             w 1 w 2 . . . w ns       =       s j (p 1 ) s j (p 2 ) . . . s j (p ns )       (2.3)
where R ij = R(θ, pjpi )

Kriging

Kriging is an interpolation method useful for replacing expensive high dimensional model with computationally efficient approximations of nonlinear functions [START_REF] Sacks | Design and analysis of computer experiments[END_REF].

The approximation of the j th component of eigensolution s j (p) can be written as the sum of a regression model and a random contribution:

s j (p) = n f i=1 β i f i (p) + W (p) = f(p) T β + W (p) (2.4)
where f i (p), i = 1...n f are monomial basis functions (constant, linear or quadratic terms) weighted by regression parameters β i . W (p) is a realization of a normally distributed random process defined with zero mean value and a covariance function, such as:

E[W (p 1 )W (p 2 )] = σ 2 R(θ, p 1 , p 2 ) (2.5)
where E[.] is the expectation operator, σ 2 is the process variance and θ > 0 a scaling vector, as previously defined for RBF.

The correlation model R(θ, p 1 , p 2 ) is a monotone function depending on the distance of two known samples p 1 and p 2 and is here defined as a product of stationary one-dimensional correlations (Table 2.1). The correlation function is equal to 1 if the points are identical and 0 for infinitely distant points and so uncorrelated.

R(θ, p 1 , p 2 ) = np j=1 R(θ j , p 1,j -p 2,j ) (2.6)
with p 1,j and p 2,j are respectively the j th parameter of samples p 1 and p 2 .

Considering an unbiased linear predictor, Eq. 2.4 can be rewritten as follows:

s j (p) = f(p) T β + r(p) T γ (2.7) where β = (F T R -1 F) -1 F T R -1 Y is the usual generalized least square estimate of β and γ = R -1 (Y -F β).
To obtain further details about the mathematical developments, the reader can refer to the following reference [START_REF] Lophaven | Dace a matlab kriging toolbox[END_REF].

First, the vector of regression functions f(p) and the vector of correlation r(p) evaluated for a untried set of parameters p are expressed as follows:

f(p) =     f 1 (p) . . . f n f (p)     (2.8) r(p) =     R(θ, p 1 , p) . . . R(θ, p ns , p)     (2.9)
Second, the matrix of regression functions F and the matrix of correlation R evaluated from the different snapshots p i , i = 1...n s are defined by: 2.11) where

F =       f 1 (p 1 ) f 2 (p 1 ) • • • f n f (p 1 ) f 1 (p 2 ) f 2 (p 2 ) • • • f n f (p 2 ) . . . . . . . . . . . . f 1 (p ns ) f 2 (p ns ) • • • f n f (p ns )       (2.10) R =       R 11 R 12 • • • R 1ns R 21 R 22 • • • R 2ns . . . . . . . . . . . . R ns1 R ns2 • • • R nsns       (
R ij = R(θ, p i , p j )
The associated optimal variance σ 2 is obtained by the function:

σ 2 = 1 n s (Y -F β) T R -1 (Y -F β) (2.
12)

The parameters β and γ depend on the scaling parameters vector θ. An optimal choice, noted θ , is determined as the solution of an optimization problem.

θ = arg min θ (det(R) 1 ns σ 2 ) (2.13)
where det(R) is the determinant of R. This equation corresponds to a maximum likelihood estimation.

Proper Orthogonal Decomposition

Proper Orthogonal Decomposition (POD) method, also denoted as the Karhunen-Loeve decomposition, allows to replace the original high-dimensional problem by a low-dimensional approximation by building an optimal linear basis through experimental or numerical data, which describes the full-system behavior. In this section, we focus our discussion on the snapshot form of POD proposed by Sirovich [Sirovich, 1987].

To define the POD model, a snapshot matrix S, as defined for the previous methods, is considered. The approximation of the eigensolution s(p) can be defined as a mean contribution obtained from the data sets S and a fluctuating contribution, described by a linear combination of n m POD modes Ψ k and associated POD coefficients a k , such as:

s j (p) = S + nm k=1 a k Ψ k (2.14)
with n m the number of POD modes.

The POD basis Ψ k of order n m is the orthonormal set which minimizes the sum of squared projection errors over all the orthonormal sets of the subspace composed of the solutions s(p i ).

POD modes Ψ k can be calculated by solving a generalized eigenvalue problem defined as follows:

Q S -S S -S T QΨ k = λ k QΨ k (2.15)
where Q is a symmetric positive definite matrix, such as Q ij = Ω ξ i ξ j dΩ with ξ i is the i th shape function associated to the finite element model.

Considering a Choleski decomposition of matrix Q = L T L, the generalized eigenvalue problem is transformed into a symmetric ordinary eigenvalue problem

YY T Φ k = λ k Φ k (2.16) with Φ k = LΨ k and Y = L S -S .
The size of the problem to be solved is directly linked to the number of studied components of eigensolution vector. To reduce the computational time, the eigenvalue problem can be defined with Y T Y matrix instead of YY T . These two problems have the same eigenvalues and their eigenvectors are related by simple relationships. Moreover, to reduce the size of the studied problem, only POD modes associated to a high eigenvalue participation are retained in the approximation. Thus, POD approximation is efficient if we have

n m << n dof .
As soon as the POD basis is obtained, each POD coefficient a k , k = 1...n m can be defined with the following scalar product for all the snapshots p i , i = 1...n s :

a k (p i ) =< s(p i ), Ψ k > (2.17)
The final step needs to efficiently link the parameter space and the POD coefficients space. A solution, based on a kriging technique, is given in reference [START_REF] Hamdaoui | Pod surrogates for real-time multi-parametric sheet metal forming problems[END_REF].

3 Reduced order model with POD modes

Projection technique

Let us now consider an eigensolution approximation using a projection of the equilibrium equations onto a subspace spanned by the columns of a rectangular projection basis T. Such projection methods have already been applied to approximate data from dynamic systems. The approximate modes of the model, projected on the basis T, are given by z = Tq, where q is the eigenvector of the reduced problem as expressed in Eq. 2.18.

T T K l T + λ cr T T K σ T q = 0 (2.18)
In this section, we propose to build the projection basis T by using the different POD modes

Ψ k . T = [Ψ 1 ... Ψ nm ] (2.19)
The size of this projection matrix is

[n dof × n m ] while the reduced problem is [n m × n m ].

Calculation of POD coefficients with Kriging

With the projection technique, a large amount of computational time is devoted to the assembly of the stiffness matrices K l and K σ (Eq. 2.18). In a way to avoid this step for a new set of parameters (during the online step), we propose here to link the space of parameters and the space of participation factors q with a kriging technique. Indeed, the participation parameters q can be seen as POD coefficients calculated with Eq. 2.17.

The proposed method relies on 4 steps.

Firstly, the POD modes Ψ 1 , ... Ψ nm are determined according to a set of snapshots P, as exposed in Section 2.3, to build the projection matrix T (Eq. 2.19).

Secondly, the reduced linear eigenvalue problem, defined by Eq. 2.18, is solved for each snapshot and the participation parameters q(p i ) are identified. The minimal number of snapshots must be equal to the number of unknowns used to define the second order regression part of kriging.

Thirdly, a new normalisation for modes, defined by Eq. 2.20, is chosen. The objective is to be able to approximate the buckling load factor without the calculation of modified geometric stiffness matrix.

(

Tq(p i )) T K σ (Tq(p i )) = λ cr (p i ) (2.20)
This new normalisation modifies the definition of participation factors for each snapshot. We obtain q :

q (p i ) = q(p i ) λ cr (p i ) z T (p i )z(p i ) (2.21)
In the fourth step, we aggregate all the participation factors q (p i ) corresponding to all the snapshots in a matrix Q (p i ), whose size is [n s × (n + 1)] and we link the two space of parameters by a kriging technique to calculate q (p) for a new sample p.

Finally, an approximation of the mode z(p) is obtained by the following equation:

z(p) = Tq (p) (2.22)
The approximation of buckling load factor λ cr is then calculated using:

λ cr (p) = (Tq (p)) T (Tq (p)) (2.23)
4 Numerical applications: buckling load evaluation of a modified structure

This section is devoted to the evaluation of the methods presented in the previous sections thanks to the test case described in Chapter 1. In this study, the finite element model is supposed to be imperfect with variations on its material and geometric parameters, such as the Young's modulus, the angle value and the two dimensions of the cross section. We globally focus our discussion on the maximal errors of the critical load and associated mode shape, calculated with

M AC (zapp,z ref ) = (z T app z ref ) 2 (z T app z app )(z T ref z ref ) (2.24) err (zapp,z ref ) = z app -z ref z ref (2.25)
where z app and z ref are the approximated and the reference solution sets, calculated with the Sorensen algorithm [C. [START_REF] Sorensen | [END_REF]. As described in Section 3.2, a discretization of 7 values for each input parameter, namely 2401 samples, were used a first for defining the reference solutions. A threshold on the maximal error was fixed at 1% to evaluate the efficiency of the methods. Moreover, to compare all the methods in the same framework, the snapshots or training samples were generated using a DOE based on a Latin Hypercube Sampling (LHS) algorithm.

The calculations were performed ten times to estimate the influence of samples position in the design space on the approximation and the same sets were used for all the comparisons.

In the next sections, six methods were tested. The five first ones, namely Quadratic model, RBF, Kriging, POD and ROM-POD, have already been investigated in several domains by different authors whereas the last one ROM-POD-Kriging, relying on both projection, POD and Kriging, is here proposed.

Case 1: Quadratic model

For this numerical application, a full quadratic regression model has been chosen to approximate the behaviour of eigensolution, presented in Fig. 2.2. As the coefficients of each monomial are identified by solving of linear system considering results of several snapshots, the key parameters of this method is the sample data set, namely the position of different samples in the design space and the number of snapshots. Obviously, the minimal size of snapshot set must be superior or equal to the number of unknown coefficients.

The errors are presented in Fig. For both cases, the error rapidly tends to a limit, namely 0.2% for two parameters and 2% for four parameters, with globally a low dispersion as a function of the set of snapshots. The studied method is efficient for the case with two parameters, for which the behavior is rather linear or quadratic, but shows some limitations if the section parameters are added. In this case, the evolution of parameters is not really quadratic but slightly nonlinear leading to some errors mainly at the bounds of parameters, as shown in Fig. 2.2. Note that these figures highlight the evolution of buckling loads and errors as a function of two parameters while the two others are fixed to the nominal values.

About the buckling modes, the maximal norm error can reach 100% with a MAC close to 0.1.

Once lower and upper bounds of input parameters are used, the level of error increases seriously because the convergence radius of approximation is too low.

On the first hand, the main time consuming step is due to the offline step relative to the solving of n s full buckling problems corresponding to the n s snapshots. Indeed, 100 snapshots were here necessary to tend to the best level of errors. On the other hand, the calculation of unknown coefficients of the regression, relying on a linear system to be solved whose size is [n s × n c ], is not time consuming. The same observation can be performed for the online step, where the approximations were obtained with a matrix vector product. The size of this matrix product is depending on the number of studied output data n o and the number of unknown coefficients n c .

In conclusion, this method is not efficient in terms of precision and computational time as soon as parameters dependency becomes highly nonlinear. Gaussian models led to a good approximation of the buckling load factor with a maximal error of 0.89% and 0.09% respectively. However, the approximation of the associated mode shape with futher precision was more difficult. Indeed, some MAC can be close to 0.5 and the norm error of the order of 200% with poor approximations for specific vector components.

Considering a Gaussian correlation model, we investigated the effect of the number of samples in the offline phase on quality of the online step (Fig. 2.3). We noticed that the maximal errors on critical buckling load are clearly lower than quadratic model even if a small number of snapshots is considered. Indeed, for the first test case with 2 parameters, the maximal error is less than 0.1% once 10 samples were chosen. From 20 samples, the level of errors does not change significantly and is globally constant. However for the second test case with 4 parameters, the maximal error decreases up to 100 samples and increases next with an oscillatory behavior. The same behavior was observed for the error variation, particularly beyond 300 samples. It was observed also, that the introduction of a large sample size does not improve the precision of results and can tend to a less efficient approximation. This observation was confirmed with the eigenvectors, for which the maximal error can reach 200% with a MAC of 0.5. As presented in Fig. 2.4, the quality of the approximation of critical buckling load is homogeneous for all the design spaces with very low level of errors when 50 samples are used. Regarding the computational time, one can address the same conclusion as previously. The main time consuming step is due to the offline step relative to the generation of snapshots. We noticed that the number of snapshots necessary to supply an approximation error up to 1% In conclusion, the use of a low number of snapshots is very interesting to provide a good approximation. Nevertheless, the level of approximation is not homogeneous with the respect of studied solutions, such as buckling load or the elements of mode. The matter does not ensure the use of the present method with a high degree of confidence.

Case 3: Kriging

As exposed in Section 2.3, the accuracy of kriging depends simultaneously on the regression, the correlation models and the definition of input sample data. We propose to successively analyze the effects of these parameters on the approximated solutions by considering, as for RBF, the same 500 random snapshots. In the following, three different orders (zero, first and second) and three correlation functions (linear, exponential and Gaussian) were respectively chosen for the regression and correlation models. The results of simulations are summarized in Table 2.3 in terms of maximal buckling load errors.

We observe a high sensitivity of the approximation quality to the model regression order.

Indeed, whatever the nature of correlation function was chosen, the errors were lower than to 2% for a second order polynomial. For zero and first orders, the level of error was significantly higher.

The maximal error was beyond to 5% for the best case. On the other hand, considering a second order polynomial, the effect of correlation function was less important. The maximal errors vary between 0.76% and 1.86%. The linear correlation model was not sufficiently accurate whereas Gaussian correlation model implied interesting error levels. The same behavior was observed for the mode shape with respect to the correlation model. Although the MAC was closed to 0.95, the norm error was slightly higher with respectively 30.28%, 28.95% and 15.65% for the linear, exponential and Gaussian correlation models respectively. In conclusion, second order polynomial regression and Gaussian correlation models were chosen for the following numerical tests.

The results of maximal errors and associated variations for buckling load factor are presented in Figure 2.5 considering the cases with 2 and 4 parameters. First of all, we observe for the two cases a decrease of both the maximal error levels and associated variations when the number of samples increases. In the first case, the error was very low even if a low number of samples is used. The maximal errors and associated variation are respectively near to 0.6% and 1% for 8 samples. These two data tend to 10 -4 % for 350 samples. However, when four material and geometric parameters were considered, the level of errors and variations are clearly higher and can reach more than 5% for 15 samples. Moreover, the decrease of errors as a function of the number of samples is slower. To obtain a maximal error inferior to 1%, 350 samples were here necessary. The observed errors tend to 0.6% with a very small variation of errors (less than 0.1%). In addition, the norm error for eigenvectors is close to 40% with a MAC close to 0.95. This second investigation highlights that the precision of data is directly linked to the nature and the number of input parameters. Finally, to achieve a good precision with kriging is enough complex although the behavior of the critical load is not highly non linear. The behavior of buckling load is rather well characterized, as it is confirmed in Fig. 2.6 by observing the response surfaces and associated errors, for all the input variable parameters.

As for RBF, the quality of the approximation about buckling loads is homogeneous for all the design space with very low level of errors when 350 samples are used. We have highlighted that the threshold of 1% can be achieved when a large number of snapshots was integrated in the simulation, namely 350 snapshots. Kriging seems to be less efficient than RBF, defined with only correlation functions, to approximate with precision the buckling load. Nevertheless, it is important to note that the buckling modes are well characterized (MAC close to 0.95 with a norm error equal to 40%) on the contrary to RBF or Quadratic regression.

The computational time of this method depends again on the number of snapshots needed in the offline step and so on the solving of n s full buckling problems. Moreover, the calculation of optimal parameters β and γ depends on matricial operations and on the solving of some linear systems of reduced size, equivalent to the regression method for β, and an optimization problem for θ . These steps are not time consuming if a low number of output variables are studied.

Nevertheless, the associated time can rapidly increase as soon as all the degrees of freedom of the Eigenvalue participation [%] 99.8 0.17 0.005 1.4 10 -6 4.9 10 -7 3.3 10 -8

Maximal errors [%] 0.85 0.58 0.55 0.52 0.52 0.52

Table 2.4 -Maximal errors on critical buckling load as a function of the number of POD modes.

eigenvector are investigated. In this application, the online time is globally equivalent to one full buckling calculation.

In conclusion, the decrease of the error is globally low when the number of input parameters increases. Moreover, the global level of approximation is rather better than RBF method, but remains unsatisfactory for the eigenvector. In addition to the time dedicated to offline step, the online step is clearly dependent on the output data to be approximated and can be significant.

Case 4: POD-Kriging

The accuracy of POD method depends on the number of selected modes useful for the approximation and on the definition of input sample data, both for the determination of POD modes and for the definition of kriging, as detailed in Section 2.4.

First of all, we consider 500 random values to define the reference snapshots set and built the POD modes. The evolution of POD eigenvalue participation and the maximal buckling load errors as a function of the number of POD modes are summarized in Table 2.4.

In the present case, we can approximate the buckling load factor with a low number of POD modes. We observe in Table 2.4 a low evolution of the error as a function of the number of modes.

Indeed, as soon as four modes are selected, the level of errors is stable at 0.52%.

Considering now 4 POD modes, we study the evolution of maximal buckling loads as a function of sample size. The results for POD-Kriging are drawn in Fig 2 .7.

The observed errors for POD-Kriging are very closed to those obtained for Kriging only.

Indeed, the lowest maximal errors are respectively for the two numerical tests around 10 -4 % and 0.6%. The response surfaces, presented Fig. 2.8 for 350 samples, confirm the observed trend, namely a good characterization of the global behavior with irreducible errors. The main difference with Kriging only is associated to the level of errors for eigenvectors, which is close to 10% with a MAC of 0.99. The integration in a correlation matrix of different eigenmodes corresponding to the snapshots and the interpolation of POD coefficients allowed improving mainly the quality of the approximation for the eigenmodes.

As previously, the same trends were observed for the computational time. The main time consuming part of the offline step was due to the calculation of reference problems according to the number of snapshots n s . In addition, the time dedicated to the eigenvalue problem, always during the offline step, is not really important as the size of the correlation matrix is dependent of the n s snapshots. Finally, Kriging method was used to approximate only some POD coefficients. During the online step, only some matricial operations were used without increasing the computational time.

In conclusion, this method supplies the best compromise between precision/time in comparison to the previous ones as soon as the eigensolutions are investigated. If the study is limited to the buckling load, no improvement was detected and RBF or Kriging alone could be used.

Case 5: ROM-POD

This section is dedicated to the calculation of the approximation of eigensolutions by considering POD modes in a projection matrix. Considering the results of the previous section, four POD modes are sufficient to obtain a good approximation. The evolution of maximal buckling loads was studied as a function of sample size and the results for ROM-POD method were drawn in Fig 2 .9.

First of all, the maximal error for the buckling load was up to 1% for the two test cases regardless the number of snapshots. When two parameters were considered, the error was constant and close to 10 -4 %. This level of error was previously found for 350 samples whereas the ROM-POD requires only 8 samples. Moreover, the variation of errors was near null. For four parameters, the error was close to 10 -2 %. This level was never obtained in the previous simulations excepted for RBF. In addition, the norm error for eigenvectors was up to 0.5% with a MAC of 1. The behavior of buckling load was captured accurately, as it is confirmed in Fig. 2.10 by observing the response surfaces and associated errors.

Using the POD modes in a projection basis was very efficient in term of precision. The level of errors and the associated variations were undoubtedly low. However, this method implies very large computational time because the modified stiffness matrices must be assembled for each new input parameter vector. The computational time was then globally the same as the reference solution. Thus, this method cannot be used in its standard form.

Case 6: ROM-POD-Kriging

The proposed method ROM-POD-Kriging, exposed in Section 3.2, relies as for the two previous methods on the number of selected modes and on the definition of input sample data. Four POD modes were always considered for the numerical application and the maximal errors and associated variations for the critical buckling load are presented in Figure 2.5.

For the two applications, the level of errors was higher than for ROM-POD but was however interesting. Indeed, we observed for the two cases a decrease of both the maximal error levels and associated variations when the number of samples increases. For two parameters, the error was around 10 -1 % for all the sample size and always less than the threshold of 1%.

When four parameters were considered, the level of errors and variations were closed to kriging one but presenting a better convergence rate. Indeed, only 35 samples were necessary to obtain a maximal error inferior to 1%. Fig. 2.6 shows that the behavior of buckling load is rather well characterized for all the design space. With 350 samples, the error can be close to 10 -1 % without variations of error due to the random generation. In addition, the error norm for eigenvectors is close to 1% with a MAC of 1. This method supplies interesting results for the approximation of the critical buckling load factor and associated buckling mode. Indeed, POD method requires 350 samples whereas ROM-POD-Kriging uses only 35 snapshots for the same level of precision. In addition, the computational time was less important than for the previous methods only based on POD modes. The consuming time offline step was always relative to the calculation of reference solutions for the n s snapshots. As the participation factors were calculated with the same stiffness matrices, the additional time is only relative to the solving of n s reduced buckling problems. This last step is not time consuming. The same assessment can be performed for the kriging interpolation of participation factors for new input data during the online step.

In conclusion, ROM-POD-Kriging method seems to be a best comprise between ROM-POD, which is very precise, and POD or Kriging alone.

Assessment of results

In this section, an assessment of the different methods tested with a same benchmark is given.

First of all, if two parameters, whatever they are, are considered all the method allows obtaining interesting and acceptable results. Indeed, the main differences between the different approximations appears when more parameters are managed. For the case of four parameters, Fig. 2.13 summarizes the maximal error on the critical loads and the needed computational time as a function of the finite element size by considering the optimal key parameters for each method.

The best precise method is obviously the ROM-POD method for which the level of error is close to 10 -4 %. Nevertheless, the associated CPU time is comparable to the reference one, that naturally excludes this numerical method. The Quadratic regression method supplies a poor level of precision in comparison to the other methods. Although the gain factor is near equal to six, this method can not be selected. The other methods, namely RBF, Kriging, POD and ROM-POD-Kriging, globally propose a same level of precision between 0.5% and 1% for the buckling load factor. Concerning the buckling mode, RBF is clearly less precise than the three other methods. About the computational time, Kriging and POD requires the same computational time, the gain factor is close to six as for the Quadratic regression. Finally, RBF and ROM-POD-Kriging allow to reduce the CPU time with a factor close to 72. Considering all the comparisons, the proposed ROM-POD-Kriging method supplies a good compromise between the precision for both buckling load factor and buckling mode and the computational time. This chapter has allowed to compare different metamodels or coupling of metamodels in the case of the approximation of linear buckling solutions and draw some advantages and drawbacks about each metamodel. Finally, the ROM-POD-Kriging method has been proposed.

Introduction

In addition to metamodels and ROM discussed in Chapter 2, a multiparametric problem can be addressed with reanalysis techniques. The main idea of this kind of methods is to perform only a nominal calculation of the full problem and to re-utilize previous results to evaluate the perturbed solutions. Each new sample is then considered as a perturbation of the nominal problem.

These approaches have an intrusive character as they required an adaptation of the equations for each studied problem. Different numerical reanalysis techniques have already been developed, such as Combined Approximation (CA) [START_REF] Kirsch | Procedures for approximate eigenproblem reanalysis of structures[END_REF], Kirsch et al., 2006, Kirsch and Bogomolni, 2007], perturbation and series development [He, 1999, Rong et al., 2003, Massa et al., 2004, Sliva et al., 2010], Padé approximants [START_REF] Damil | An iterative method based upon padé approximants[END_REF], Yang et al., 2001, Duigou et al., 2003, Massa et al., 2008b, Massa et al., 2009a] or different adhoc projections [Balmès, 1996, Bouazzouni et al., 1997, Corus et al., 2006, Boyaval et al., 2009, Nouy, 2010]. The spectrum of development and application is relatively large in static, modal and dynamic. Indeed, Lallemand et al [START_REF] Lallemand | Fuzzy modal finite element analysis of structures with imprecise material properties[END_REF] utilized Neumann series expansion to treat eigenvalue problems including fuzzy parameters. Massa et al [START_REF] Massa | A fuzzy procedure for the static design of imprecise structures[END_REF], Massa et al., 2008a] integrated series development and Padé approximations in fuzzy static and modal problems. Qiu and

Elishakoff [START_REF] Qiu | Antioptimization of structures with large uncertain-but-non-random parameters via interval analysis[END_REF]] propagated uncertainty with combination of the perturbation technique and interval mathematics to determine the region of static response surface. Chen et al [START_REF] Chen | Hybrid uncertain analysis for the prediction of exterior acoustic field with interval and random parameters[END_REF] employed first-order series to study the exterior acoustic field prediction with interval and random variables.

For the solving of multiparametric linear eigenvalue problem, several developments based on homotopy perturbation and projection technique [START_REF] Massa | Structural modal reanalysis methods using homotopy perturbation and projection techniques[END_REF][START_REF] Massa | Multi-level homotopy perturbation and projection techniques for the reanalysis of quadratic eigenvalue problems: The application of stability analysis[END_REF], Massa et al., 2016] have already been successfully applied for vibration problems, where modifications are integrated in mass and stiffness matrices. In this chapter, the performance of these reanalysis techniques, such as modal stability, series development, Padé approximations, homotopy developments or residue iteration, are evaluated for the case of modified linear buckling problems with the numerical framework used in Chapter 2. In the present case, the material and geometric parameters variations are already integrated in linear and geometric stiffness matrices.

After identifying the advantages and drawbacks of these methods, a new method, coupling high order perturbed modes, projection and Kriging, will be exposed.

2 Perturbed linear buckling analysis

Nominal case

Let consider a nominal linear buckling problem described in a finite element context corresponding to a nominal set of parameters noted

p (0) = [p (0) 1 ... p (0)
np ]. This problem is numerically solved considering three main steps. The first one is the solving of a system of linear equations considering the nominal linear stiffness matrix K (0) l and the external load F ext to obtain a nominal displacement vector U (0) :

K (0) l U (0) = F ext (3.1)
The second step is dedicated to the assembling of the nominal geometric stiffness matrix K (0) σ by considering the nominal stress state σ (0) , calculated with the nominal displacement vector

U (0) .
The last step is relative to the calculation of the nominal buckling load factor λ

(0)
cr by solving of a linear eigenvalue problem considering the nominal linear and geometric stiffness matrices.

K (0) l + λ (0) cr K (0) σ z (0) = 0 (3.2)
where λ

(0)
cr is the nominal minimum positive eigenvalue of Eq 3.2, z (0) is a non-zero vector which denotes the nominal eigenvector reflecting the structural buckling shape. The nominal critical buckling load F (0) cr can be obtained by multiplying the initial external load F ext by the nominal buckling load factor λ (0) cr .

F (0) cr = λ (0) cr F ext (3.3)
These two eigensolutions z (0) and λ

(0)
cr are aggregated in a vector s (0) , such as:

s (0) = z (0) λ (0) cr (3.4)

Perturbed case

Let consider now some modified input parameters

p (m) = [p (m) 1 ... p (m)
np ] which define the modified model problem. In this section, the exponent (m) refers to the modified state whereas the exponent (0) identifies the previous nominal one.

The equations associated to the modified buckling problem are summarized as follows:

K (m) l U (m) = F ext (3.5) K (m) l + λ (m) cr K (m) σ z (m) = 0 (3.6) s (m) = z (m) λ (m) cr (3.7)
where

K (m) l , K (m) σ , U (m) , λ (m) cr
and z (m) are respectively the modified linear and geometric stiffness matrices, the modified displacement, the modified buckling load factor and the modified eigenvector.

Perturbation techniques

The perturbation techniques can be independently applied for the different steps of calculation procedure, namely for the modified linear static step and/or for the modified eigenvalue problem.

The different technical solutions are exposed in the next subsections.

Modal stability method

The modal stability method, proposed by Lardeur [START_REF] Arnoult | The modal stability procedure for dynamic and linear finite element analysis with variability[END_REF], has been applied

to different applications domains [START_REF] Druesne | Fast methods based on modal stability procedure to evaluate natural frequency variability for industrial shell-type structures[END_REF], Druesne et al., 2016] and integrated in uncertainty propagation method in the case of modified linear eigenvalue problem. This method relies on the assumption that the mode shapes of the structure are assumed to be weakly sensitive to variations of the input parameters.

The modified eigenvector z (m) can be expressed as the nominal contribution z (0) , defined by the Eq. 3.2, and a perturbed vector, noted ∆z:

z (m) = z (0) + ∆z (3.8)
Considering the expression of the Rayleigh quotient, classically used in modal analysis, for the modified eigenvalue problem and Eq. 3.8, we can derive the expression of the modified buckling load factor λ (m) cr :

λ (m) cr = - (z (0) + ∆z) T K (m) l (z (0) + ∆z) (z (0) + ∆z) T K (m) σ (z (0) + ∆z) (3.9)
Therefore, according to the modal stability assumption, the modified eigenvector z (m) is assumed to be equal to the nominal vector z (0) . Eq. 3.9 is then reduced to:

λ (m) cr = - z (0)T K (m) l z (0) z (0)T K (m) σ z (0) (3.10)
Thus, the perturbations are only present on the two stiffness matrices, which can be calculated for each parametric modifications. To reduce the computational time, the authors propose to replace the numerator and denominator of this fraction by using equivalent energies to avoid the assembling step of modified matrices. In our numerical applications, only Eq. 3.10 will be used to test the efficiency of this kind of approximation.

Homotopy developments

On the contrary to the previous section, the aim is here to take into account the contribution associated to the eigenvector perturbation ∆z in the case of eigenvalue problem, by identifying some high order vectors. The following techniques can be applied both for modified system of linear equations and for modified linear eigenvalue problem.

System of linear equations

The integration of perturbations in the linear stiffness matrix K (m) l

, useful to apply Eq. 3.5

governing a linear static problem, inevitably imply the modification of nodal displacements vector m) , stress state of the studied structure σ (m) and so the geometric stiffness matrix K (m) σ . Thus, for each perturbation, it is necessary to update the geometric stiffness matrix to determine the associated global level of perturbation.

U (
The homotopy perturbation technique expresses a non-linear problem as a set of linear problems, by defining a convergent series solution. The technique consists in introducing an additional unknown parameter ε to highlight the non-linearity. At the end of the formulation, this parameter is set to 1 in order to evaluate the perturbed problem. In the present case, the non-linearity is relative to the introduced perturbations on finite element matrices.

The modified linear stiffness matrix K (m) l can be expressed as nominal contribution K (0) l and a perturbation ∆K l such as:

K (m) l = K (0) l + ε∆K l (3.11)
The modified displacement vector U (m) is developed as a series considering a nominal dis-placement vector U (0) and different perturbed high order vectors U (j) , such as:

U (m) = U (0) + εU (1) + ... + ε (n) U (n) (3.12)
where n is the order of truncation of the series.

After having introduced Eqs 3.11 and 3.12 in Eq. 3.5, we obtain:

(K (0) l + ε∆K)(U (0) + ε (1) U (1) + ... + ε (n) U (n) ) = F ext (3.13)
We can identify the nominal displacement vector and perturbed vectors by considering each order of parameter ε. For the nominal displacement solution U (0) , we detect the equation relative to the nominal static problem, already exposed in Section 2.1 by Eq 3.1.

The perturbed vectors U (j) are then calculated by considering the following recursive equation:

K (0) l U (j) = -∆KU (j-1) (3.14)

Linear eigenvalue problem

For the modified linear eigenvalue problem, the eigensolutions λ (m) and z (m) are expressed, as previously, in series expansion with respect to the homotopy parameter ε, respectively by Eq. 3.15 and 3.16.

λ (m) cr = λ (0) cr + ελ (1) cr + ... + ε (n) λ (n) cr (3.15) z (m) = z (0) + εz (1) + ... + ε (n) z (n) (3.16)
Unlike the modal stability method, the aim is here to take into account some modifications on mode shape between the nominal and perturbed state by integrating higher order contributions

z (j)
Moreover, the matrices K (m) l and K (m) σ (Eq. 3.11 and 3.17) are decomposed into matrices K (0) l and K (0) σ , associated to an initial eigenvalue problem, and perturbed matrices ∆K l and ∆K σ , which include all the modifications identified on the material and geometric parameters.

K (m) σ = K (0) σ + ε∆K σ (3.17)
In complement to the modified eigenvalue problem, a normalization equation, relative to modified data, is added.

z (m)T K (m) σ z (m) = 1 (3.18)
By introducing Eq. 3.11, Eq. 3.15 to 3.17 in Eq. 3.6 and 3.18, it is possible to identify the perturbed eigensolutions by considering each order of parameter ε.

((K (0) l + ε∆K l ) + (λ (0) cr + ελ (1) cr + ... + ε (n) λ (n) cr )(K (0) σ + ε∆K σ )) (z (0) + εz (1) + ... + ε (n) z (n) ) = 0 (3.19) (z (0) + εz (1) + ... + ε (n) z (n) ) T (K (0) σ + ε∆K σ )(z (0) + εz (1) + ... + ε (n) z (n) ) = 1 (3.20)
For 0-order, we obtain the nominal eigenvalue problem, already defined by Eq. 3.2, and associated nominal normalization equation (Eq. 3.21).

z (0)T K (0) σ z (0) = 1 (3.21)
For the j th order, the high order perturbed eigensolutions are calculated with the following linear system of equations, depending on nominal characteristics and different perturbations:

K * (0) s (j) = F * (j-1,...,0) (3.22) with K * (0) = K (0) l + λ (0) cr K (0) σ K (0) σ z (0) z (0)T K (0) σ 0 (3.23) s (j) = z (j) λ (j) cr (3.24) F * (j-1,...,0) =        -∆K l z (j-1) + j-1 k=0 λ (k) cr ∆K σ z (j-k-1) + j-1 k=1 λ (k) cr K (0) σ z (j-k) -0.5 j-1 k=0 z (k)T ∆K σ z (j-k-1) + j-1 k=1 z (k)T K (0) σ z (j-k)        (3.25)
These high order perturbed vectors s (j) can then be integrated in different approximations such as series, Padé approximants or projection techniques, exposed in the next sections.

Approximation with series and Padé approximants

Considering the developments exposed in Section 3.2, the modified static displacement U (m) can be approximated at first by Eq. 3.12 to calculate the modified geometric stiffness matrix

K (m) σ .
Second, considering the perturbations on stiffness matrices, the modified critical buckling load λ (m) cr is approximated by a series described by Eq. 3.15. This solution, already discussed in several papers [Massa et al., 2008b, Massa et al., 2009b, Massa et al., 2011a], can be efficient when the level of variations of studied solutions is very low and the behavior is clearly linear.

For a more general case, more robust techniques must be considered for integrating significant perturbations on material and geometric parameters.

For instance, Padé approximants, have been already used in other numerical methods such as Asymptotic Numerical Method to extend the convergence of the developement, are considered to determine the modified studied solutions as U (m) , z (m) or λ (m) cr . Thus, the series can be replaced by a rational fraction function. In the present case, we consider Padé approximants with a fixed denominator. In the case of displacement vector, we obtain:

U (m) = U (0) + ε D n-2 (ε) D n-1 (ε) U (1) + ... + ε (n-1) D 0 (ε) D n-1 (ε) U (n-1) (3.26)
where (D j (ε)) j=0...n-1 are polynomials of degree j with real coefficients d j , such as D j (ε) = d 0 +εd 1 +...+ε (j) d j . The coefficients d j are identified using an iterative Gram-Schmidt algorithm, recalled in reference [Massa et al., 2008b].

The same procedure can be applied for the critical buckling load λ (m) cr and associated mode z (m) .

Projection techniques

Reduced problems

As in section 3.1, another possibility to calculate modified solutions is to consider a projection of the equilibrium equations onto a subspace spanned by the columns of a rectangular projection bases. In the present case, two projections matrices, noted T U and T z , are respectively considered for the static and linear eigenvalue problems. Thus, the different steps of a buckling problem can be respectively rewritten considering reduced finite element matrices as follows:

T T U K l T U q U = T T U F ext (3.27) U = T U q U (3.28) T T z K l T z + λ cr T T z K σ T z q z = 0 (3.29) z = T z q z (3.30)
where q U and q z are respectively the reduced vector of displacement and mode shape.

The next sections are dedicated to the definition of the reduction matrices T U and T z .

Perturbation projection matrix

Firstly, the matrices T U and T z can be directly composed on nominal data and associated perturbed vectors [Massa et al., 2008a[START_REF] Massa | Structural modal reanalysis methods using homotopy perturbation and projection techniques[END_REF], Do et al., 2016], such as: method. This method will be consider in the next section.

T U = [U (0) U (1) ... U (n) ] (3.31) T z = [z (0) z (1) ... z (n) ] ( 3

Residue iteration

Secondly, Balmes [START_REF] Bobillot | Iterative techniques for eigenvalue solutions of damped structures coupled with fluids[END_REF] proposes a reanalysis technique, based on an error control and an enrichment mechanism, compatible with the static and eigenvalue problems.

The first step is to start this numerical strategy with the projection matrices, T U and T z , respectively initialized to the nominal data U (0) , z (0) and λ

(0)
cr to calculate a first approximation of studied modified solutions, using Eqs. 3.27 and 3.29.

Next, load residuals, R L (U (m) ) and R L (λ (m) cr , z (m) ), are calculated with the first approximated solutions for each problem, such as:

R L (U (m) ) = K (m) l U (m) -F ext (3.33) R L (λ (m) cr , z (m) ) = K (m) l + λ (m) cr K (m) σ z (m) (3.34)
where U (m) and z (m) are calculated with Eqs. 3.28 and 3.30.

To enrich each projection basis, we define displacement residuals using the inverse of nominal linear stiffness matrix K (0) l and the load residuals:

R D (U (m) ) = K (0) l -1 R L (U (m) ) (3.35) R D (λ (m) cr , z (m) ) = K (0) l -1 R L (λ (m) cr , z (m) ) (3.36)
Finally, the projection bases T U and T z are expressed at the j th iteration with displacement residuals as follows:

T (j+1) U = T (j) U R D (U (m) ) (3.37) T (j+1) z = T (j) z R D (λ (m) cr , z (m) ) (j) (3.38)
The number of used residuals, noted n R , is comparable to the order of high order perturbed vectors, described in the previous section, and can be controlled using a relative error, as proposed by the authors in [START_REF] Bobillot | Iterative techniques for eigenvalue solutions of damped structures coupled with fluids[END_REF]. The size of the reduced matrices is so [(n R + 1) ×

(n R + 1)].

Fixed projection basis

Considering Eqs. 3.31 and 3.32, we have seen that it is possible to define a reduced projection space, which allows to calculate the modified buckling load factor as a function of input parameters variations. With this strategy, the computational time is naturally reduced through a classical calculation because only one factorisation of a reference nominal matrices

K (0) l and K * (0)
is necessary (Eqs. 3.14 and 3.22). However, the gain factor remains limited due to the calculation of the necessary high order perturbed vectors for each new modifications and to orthonormalise the projection matrix for each new modified problem. To avoid these limitations, we propose here to build a projection basis compatible with a set of modifications representing the studied design space rather than only one modification [START_REF] Massa | From homotopy perturbation technique to reduced order model for multiparametric modal analysis of large finite element models[END_REF], Do et al., 2017].

Considering the k th modified input parameter p

(m)
k , used to define the modified linear and geometric stiffness matrices, this one can be decomposed as a nominal part and a perturbation multiplied by the homotopy parameter ε:, such as:

p (m) k = p (0) k + ε∆p k (3.39) Each modified parameter p (m) k
generate some perturbations on stiffness matrices, noted

∆K l (p (m) k ) and ∆K σ (p (m)
k ), such as:

K l (p (m) k ) = K (0) l + ε∆K l (p (m) k ) (3.40) K σ (p (m) k ) = K (0) σ + ε∆K σ (p (m) k ) (3.41)
Next, it is so possible to perform a homotopy development for each parameter p As the homotopy development is efficient for a range of values defined around the nominal calculation along the parameter direction p k according to the homotopy parameter ε, only one value of parameter p k is necessary to characterize the high-order vectors per parameter direction. In practice, only the minimal and/or maximal bound of variations of parameter p k will be considered to calculate the high order vectors.

To build the projection matrices, we propose here to concatenate data identified for all the parameter directions p k and suppress the possible collinear vectors during the orthonormalization step. Moreover, the nominal part is always fixed in the definition of the projection basis whatever the modification is introduced.

The projection bases are so defined as follows:

T U (p 1 , ..., p np ) = [U (0) U (1) p 1 ... U (n) p 1 ... U (1) pn p ... U (n) pn p ] (3.42) T z (p 1 , ..., p np ) = [z (0) z (1) p 1 ... z (n) p 1 ... z (1) pn p ... z (n) pn p ] (3.43)
The maximal size of the reduced matrices is [(n p × n + 1) × (n p × n + 1)] because the number of collinear vectors is dependent on the studied problem.

Coupling of homotopy developments and kriging

The method, proposed in Section 3.3, allows to reduce the size of problems to be solved and so the associated computational time, but requires the calculation of modified matrices for each new samples. To improve the previous proposal, we propose here to couple the homotopy developments for the linear eigenvalue problem and kriging techniques, exposed in Section 2.3 for POD modes.

As previously, the proposed method relies on 4 steps.

Firstly, the projection basis T z (p 1 , ..., p np ) (Eq. 3.43) is built with high order perturbed vectors per parameter direction p (m) k . Secondly, the reduced linear eigenvalue problem, defined by Eq. 3.29 is solved for each snapshot and the variable q z (p i ) is identified. The variable q z (p i ) can be seen as a participation factor variable which traduces the participation of each high order perturbed modes to calculate the associated modified solutions. The minimal number of snapshots must be equal to the number of unknowns used to define the second order regression part of kriging. The snapshots are generated by a LHS algorithm as previously.

Thirdly, the normalisation equation Eq. 3.18 associated to the modified eigenvalue problem, is changed such as:

(T z q z (p i )) T K (m) σ (T z q z (p i )) = λ (m) cr (p i ) (3.44)
The objective is to be able to approximate the modified buckling load factor without the knowledge of the modified geometric stiffness matrix. This new normalisation implies the definition of new participation factor q z (p i ) for each snapshot, such as:

q z (p i ) = q z (p i ) λ (m) cr (p i ) z (m)T (p i )z (m) (p i ) (3.45)
In the fourth step, we aggregate all the participation factors q z (p i ) corresponding to all the snapshots in a matrix Q z (p i ), whose size is [n s × (n + 1)] and we link the two space of paramaters by a kriging technique to calculate q z (p) for a new sample p.

Finally, an approximation of the mode z(p) is obtained by the following equation:

z(p) = T z q z (p) (3.46)
The approximation of buckling load factor λ cr is so calculated:

λ (m) cr (p) = (T z q z (p)) T (T z q z (p)) (3.47)
4 Numerical applications: buckling load evaluation of a modified structure

For this numerical section, only the case with four parameters is here considered because we have shown in chapter 2 that the case with two parameters is not discriminant to test the efficiency of the different methods in the same condition. The same numerical framework is considered by keeping the same random snapshots and the other key parameters.

Case 1: Modal Stability

On the contrary to previous exposed methods, modal stability method is not dependent on specific parameters. Only the nominal eigenvector and matrices perturbation are considered. To test the efficiency of this method, specifically dedicated to the eigenvalue problem, the linear static step is calculated with the reference way to avoid the introduction of errors.

The response surfaces and associated errors are presented in Fig. 3.1 for the two sets of parameters, namely for Young's modulus and angle, and width and height geometric parameters.

In Fig. 3.1, we can shown that the error for the approximation of buckling load is enough low for Young's modulus and width of section. On the contrary, the level of errors increases seriously and can achieve 100% for parameters, such as angle or height of the section. In this last case, the perturbations introduced in geometric parameters and so in geometric stiffness matrix do not allow to verify the assumption of stability. Indeed, a MAC of 0.9 with a norm variation of 35% is observed between two reference configurations.

This method, which has already successfully applied to approximate the first mode shapes in modified modal analysis, cannot be considered as an suitable alternative to the reference method for the reanalysis of modified buckling problems with topological parameters.

Case 2: Series and Padé approximants

The approximation of buckling load factor with series or Padé approximants is conditioned with the order of truncation n, as detailed in Section 3.2. Thus, we first study this effect on the evaluation of buckling load by successively applying series or Padé approximants both for the reanalysis of the static solution and the eigenvalue problem. The maximal errors on the buckling loads with series and Padé are respectively exposed in Tables 3.1 and 3.2.

Globally, the quality of the approximation is clearly poor. The least maximal error, equal to 11%, is observed for the seventh order for each approximation. Moreover, when we study the first row of Table 3.1, which is dedicated to the reanalysis of the static step and a calculation of the critical load factor with the reference method, a divergence of the series is highlighted. For a first column, associated to the calculation of the static step with the reference method and the reanalysis of the critical load factor, a slow convergence of the series is detected. In conclusion, a direct operating of the high order solutions in a series is not enough sufficient to answer to our problematic.

To improve the efficiency of the approximation, we introduce the high order data in rational fractions. As it can be seen in Table 3.2, the results are of better quality. First, each approximation is convergent independently. Indeed, we observe, respectively in the first row and column of the Table, some errors of 0.81% for a seventh order for static step and 0.07% for a seventh order for eigenvalue reanalysis. By combining the two approximations, a maximal error of 0.68% is obtained for the buckling load factor.

However, we observe that the level of errors is in general very high when the two approximations are used together. Only a seventh order of truncation can be chosen in the present case. In this case, the MAC value is close to 0.9 with an error of 30%. The results are in line with already published results which discuss of the optimal order of truncation with Padé approximants. The order of truncation is obviously dependent of the studied problem but must generally be high.

This choice is unfortunately incompatible with the researched strategy which is to reduce the computational time too. Indeed, the use of Padé approximants requires the same computational time as the reference method. In conclusion, the use of perturbations in series or in Padé approximants is not efficient either for the precision, or for the computational time. The nature of these mathematical functions does not allow to get a rapid convergence. To achieve our objectives, a discussion of projection techniques, where high order solutions can be introduced in projection matrix, is proposed in the next section.

Case 3: Projection technique with HPP modes

This section is dedicated to the evaluation of the capability of projection techniques with HPP modes. First, the maximal errors for approximation of buckling load factor is exposed in Table 3.3 as a function of the orders of truncation for both static and eigenvalue problems.

Considering the same high order solutions, as in previous section, we globally obtain interesting levels of errors with a best convergence of each approximation. Indeed, for the static step, the errors are close to 3.44% for an order 3 and 0.21% for an order 5. This result is better than with Padé approximants. Moreover, a maximal error of 0.03% is detected only for an order 3 in the case of eigenvalue problem. The use of high order solutions as a projection space clearly improve the quality of the approximation.

After coupling the two projection techniques, the maximal errors of buckling load factors are inferior to 1% as soon as an fifth order is chosen for static step and a third order for eigenvalue problem. In addition, the MAC is equal to 1 for all the configurations with a norm error inferior to 0.5%. This approximation is very efficient to investigate both eigenvalue and eigenvector. It is important to note that the main sensitive step is not the eigenvalue problem which can appear more complicated as static one but rather the approximation of the linear static problem with large precision. Indeed, a poor approximation of displacements introduces an error in the geometric stiffness matrix and the associated perturbation is not well quantified. About the computational time, as the modified stiffness matrices are necessary for each modified set of input parameters, the computational gain factor is not significant. The CPU time is globally equal to the reference one.

Case 4: Residue iteration

Another solution is to use the Residue iteration projection to improve step by step the quality of approximation by successively increasing the size of the projection matrix. Classically, a convergence criterion is used to stop the calculation and obtain with precision the researched solution. In this section, we rather consider the size of projection matrix and so the number of introduced residue as criterion in order to directly compare the efficiency with previously proposed method. A residue is so equivalent to the order of truncation because each method requires the solving of a system of linear equations to calculate the needed data. The maximal errors on buckling load as a function of number of residue is presented in Table 3.4.

First, the behavior of this method is very close to that observed for the projection with HPP modes. Indeed, with 5 residues for the static problem and 3 residues for the eigenvalue one, the maximal errors can reach 0.22%. In the previous case, we have highlighted that an order 5 and an order 3, respectively for the static and the eigenvalue problems, was a good choice. We observe the same convergence too. For the buckling mode, the results are very precise too with a MAC equal to one and a norm errors inferior to 1%.

About the computational time, the observations are the same as for the projection with HPP modes. As the orders of truncation are lower than for Padé approximants, we observe a decrease of the computational time of 15% for each studied samples. The cost associated to the building of the projection matrices and the solving of the reduced problems are less time consuming as the reference method. However, this reduction is not significant because the projection matrix must be updated for each new perturbation.

Case 5: ROM-HPP method

The proposed method ROM-HPP relies on the high order perturbed vectors calculated per direction of parameters for both the static and eigenvalue problems. First, the maximal errors of buckling loads approximation are investigated as a function of the order of truncation for each problem. The results are described in Fig. 3.2. First of all, the convergence of error is clearly different between static and eigenvalue problems.

As already observed in the previous section, the static problem requires an 23 th order to guarantee the maximal error inferior to 1%. The static step is clearly the most important step to maintain the level of error for the global calculation. Indeed, an third order is sufficient for the eigenvalue problem to obtain the targeted threshold. Considering this development in the two projection matrices, it is possible to approximate with success the response surfaces of buckling load as a function of the different parameter, as shown in Fig. 3.3.

With this proposed method, the objective of precision is achieved both for the buckling load and the associated mode. For the last case, the MAC is equal to 0.99 with an error of 2%.

Nevertheless, Figure 3.3 show clearly that the approximation of the buckling loads as a function of angle and height section parameters is a little more difficult in comparison to Young's modulus or width section.

About the computational time, the assessment is obviously less interesting because it is necessary to obtain the modified matrices for each new sample, as for the ROM-POD method, even if the matrix projection is fixed.

Case 6: ROM-HPP-Kriging method

This section is dedicated to ROM-HPP-Kriging method, which is comparable to the proposed method ROM-POD-Kriging, exposed in the chapter 2. This new method depends on the same key parameters, namely the number of modes of the projection matrix, defined as a function of the order of truncation, and on the definition of input sample data.

First, the maximal errors of buckling loads approximation are investigated as a function of the order of truncation for each problem. The results are described in Fig. 3.4.

To tend to the targeted threshold, an third order is sufficient to define the projection matrix by considering perturbed modes per parameter direction. In the present case by considering 4 parameters and a third order, the number of modes is equal to 13. However, this number of We observe that the maximal errors are inferior to 1% as soon as 25 snapshots are considered.

If more snapshots are taken into account, the mean error become stable. However, an increase of the error variation is observed. This approximation is rather efficient to investigate both the buckling loads, as shown in Fig. 3.6, and the associated modes, where the MAC is equal to 0.99 with a norm error inferior to 1%. In conclusion, ROM-HPP-Kriging method is the best comprise to approximate eigensolutions when perturbations techniques are considered. This chapter has allowed to test several methods relying on perturbations and projection techniques for the same numerical benchmark. First, we have shown that the Modal Stability assumption cannot be considered for buckling problems. The perturbations observed on geometric stiffness matrices do not allow to verify this assumption. Second, the techniques based on Series or Padé approximants have a low convergence radius or requires very high order of developments. In this last case, the computational gain is not interesting. Third, the projection techniques, based either on HPP modes or on residues, are very precise for approximating both the critical buckling loads and associated buckling modes. The cost associated to the building of the projection matrices and the solving of the reduced problems are less time consuming as the reference method.

However, this reduction is not significant to be used in multiparametric analyses. Only the ROM-HPP-Kriging method supplies the compromise in terms of precision and computational gain. Indeed, Fig. 3.7 summarizes the maximal error on the critical loads and the needed computational time as a function of the finite element size by considering the optimal key parameters for HPP, ROM-HPP and ROM-HPP-Kriging. The results relative to the kriging method, exposed in chapter 2, are also added on the figure. ROM-HPP-Kriging allows to reduce the CPU time with a factor close to 74 in comparison to the reference. For recall, the direct use of Kriging gives only a decrease of a factor 7. This computational gain is even more interesting if we consider that the precision is inferior to the threshold of 1% whatever the number of degrees of freedom considered for the numerical model.

Introduction

Another way to approximate the modified buckling loads and the associated eigenmodes in a multiparametric linear buckling problems is to calculate the derivatives of eigensolutions and perform Taylor series expansion or projection techniques based on eigenmodes derivatives. In the litterature, several works, mainly for the case of uncertainty propagation, have already integrated these techniques. Concerning probabilistic uncertainty propagation method, Taylor series expansion [START_REF] Benfratello | A perturbation approach for the response of dynamically modified structural systems[END_REF], Chakraborty and Dey, 1998, Nieuwenhof and Coyette, 2002, Falsone and Ferro, 2005] have been considered with the aim to substitute the equations of numerical model and the random functions by the mathematic developments to quantify the first statistical moments (mean, variance) of output solutions [Sudret, 2007].

Massa et al [START_REF] Massa | A fuzzy procedure for the static design of imprecise structures[END_REF], Massa et al., 2008b] integrated Taylor series development in the case of a fuzzy modelisation of the variability for static and modal problems. Qiu and Elishakoff [START_REF] Qiu | Antioptimization of structures with large uncertain-but-non-random parameters via interval analysis[END_REF]] propagated uncertainty with combination of the series expansion and interval mathematics to determine the region of static response surface. Chen et al [START_REF] Chen | Hybrid uncertain analysis for the prediction of exterior acoustic field with interval and random parameters[END_REF] employed first-order Taylor series to study the exterior acoustic field prediction with interval and random variables. Hong et al [START_REF] Hong | Parametric reduced-order models for predicting the vibration response of complex structures with component damage and uncertainties[END_REF] proposed the multiplecomponent parametric reduced-order models for predicting the vibration response of complex structures with parametric variability in multiple components. In this case, the modified component matrix was described by the Taylor series development. All these methods are dependent on the validity of Taylor series expansion and have been mainly applied for linear static and dynamic problems [START_REF] Handa | Application of finite element methods in the statistical analysis of structures[END_REF], Liu et al., 1986, Shinozuka and Yamazaki, 1988] with low output dispersions around the nominal values as shown by Ghanem and Spanos [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF].

To improve the domain of validity of the approximation and the quality of the approximation, Massa et al [START_REF] Massa | Structural modal reanalysis methods using homotopy perturbation and projection techniques[END_REF] proposed to integrate derivatives of eigenvectors in a projection matrix and solve some reduced linear eigenvalue problems in vibration instead of directly using Taylor series expansion. In this chapter, a comparison of different developments, based on derivatives and Taylor series expansion, is firstly proposed using the previous numerical framework.

The aim is to analyze the possibility of extension of previous works, performed in vibration, to the domain of linear buckling.

In addition, as a key step of this strategy is the calculation of derivatives of linear and geometric stiffness matrices, we propose to investigate the approximation of these finite element matrices to obtain a parametric description of the studied problem by using a set of snapshots. Finally, a numerical strategy, comparable to those proposed in the previous chapters, is proposed coupling 2 Decomposition of finite element matrices

General form

In this chapter, we assume that a decomposition of finite element matrices K l (p) and K σ (p) as a function of parameters functions (Equations 4.1 and 4.2) is available considering the following form:

K l (p) = n f k=1 f k (p)K lk (4.1) K σ (p) = n f k=1 f k (p)K σk (4.2)
where f k (p), k = 1...n f are monomial basis functions (constant, linear or quadratic terms) as used in kriging method. K lk and K σk are constant matrices including the coefficients of regression.

Approximation of stiffness terms

To perform multiparametric analyses, it is interesting to have a model description of finite element matrices in terms of input parameters p to describe the evolution from the initial model as a linear combination, as described in Eqs. 4.1 and 4.2. The nonlinear behavior of parameters is thus included in monomial basis functions.

This decomposition is easily performed for specific parameters of linear stiffness matrix, such as Young's modulus or density, possible for a thickness of a plate by decomposing the membrane and bending properties and the coupling effects but not directly available for topological modifications, for example, which modify the coordinates of nodes. For the geometric stiffness matrix, this decomposition is not directly possible whatever the parameter is studied because of the dependency on the strain field.

Instead of using snapshots to directly define the approximation method as in chapters 2 and 3, we propose here to approximate the non null terms of the stiffness matrices using the regression technique, presented in section 2.1 of chapter 2. The order of approximation can be chosen parameter by parameter. For example, a first order is considered for Young's modulus or density whereas a third or fourth order is selected for geometric modifications. To calculate the unknown coefficients of the regression, a set of snapshots is considered as for the metamodels in chapter 2.

For the linear stiffness matrix, the computational time relative to model description can be reduced by considering some family of elements defined with the same geometric properties.

In this case, the parametric analysis can be performed for one elementary stiffness matrix and extended for the other ones. On the contrary, no specific improvement can be performed for the geometric stiffness matrix.

3 Taylor series expansion

Approximation of eigensolutions

In this section, we kept the same exponents (0) and (m) for respectively referencing the nominal and modified state as in chapter 3. Let's consider a Taylor series expansion of eigensolutions s (m) , described Eq. 4.3:

s (m) = s (0) +ε np k 1 =1 ∂s (0) ∂p k 1 (p (m) k 1 -p (0) k 1 )+ ε 2 2 np k 1 =1 np k 2 =1 [ ∂ 2 s (0) ∂p k 1 p k 2 (p (m) k 1 -p (0) k 1 )(p (m) k 2 -p (0) k 2 )+...+R n (p 1 , ..., p np )] (4.3)
with, for example, for the first order:

∂s (0) ∂p k =      ∂z (0) ∂p k ∂λ (0) cr ∂p k      (4.4)
R n (p 1 , ..., p np ) is a residual at the order n. A homotopy parameter ε is introduced for each order of the expansion as in Chapter 3. If ε is equal to zero, Eq. 4.3 defines the nominal eigensolution.

For ε equal to one, we obtain a Taylor series expansion of the variable s (m) . This variable will be used in a next section to identify the high order perturbed eigensolutions, already defined in Chapter 3.

Calculation of derivatives

The eigensolutions derivatives depend on the derivatives of the nominal finite element stiffness matrices, which can be written as follows by considering the form proposed in Eqs. 4.1 and 4.2:

∂K l ∂p k = n f k=1 ∂f k (p (0) ) ∂p k K lk (4.5) ∂K σ ∂p k = n f k=1 ∂f k (p (0) ) ∂p k K σk (4.6)
The calculation of n th derivative of eigensolutions

∂ n s (0)
∂p n k according to p n k is performed with a method comparable to the calculation of high order perturbed vectors (Section 3.2 in chapter 3) and requires the solving of a linear system for each order of derivatives:

K * (0) ∂ n s (0) ∂p n k = ∂ n F (0) ∂p n k (4.7)
For a first order, the vector

∂F (0)
∂p k is given for example by :

∂F (0) ∂p k =         -   ∂K l ∂p k z (0) + ∂K σ ∂p k z (0) λ (0) cr   -0.5   z (0)T ∂K σ ∂p k z (0)           (4.8)
4 Reduced order model with high order derivatives

Reduced order model definition

Considering the specific form of finite element matrices, the reduced eigenvalue problem can be defined as follows:

n f k=1 f k (p) T T z K lk T z + λ cr n f k=1 f k (p) T T z K σk T z q z = 0 (4.9)
with T T z K lk T z and T T z K σk T z are reduced constant matrices calculated only one time during the offline step.

Derivatives in projection matrix

Considering the reduced eigenvalue problem, described in Eq. 4.9, the projection basis can be defined with the nominal eigenvector and the associated high order eigenvector derivatives truncated at the order n, such as:

T z =   z (0) • • • ∂z (0) ∂p k • • • ∂ n z (0)

Calculation of perturbed modes with derivatives

Directly collecting all the eigenvector derivatives in the projection basis T z is not advantageous since the number of parameters is high for a truncation order n. To control the size of the basis and reduce the CPU time, we proposed to express the high-order terms in terms of eigensolution derivatives in order to maintain the reduced basis at the same size as of the ROM-HPP method, described in chapter 3.

After developing Eq. 4.3 for buckling mode and reorganizing with respect to the different orders of ε, we can identify a link between eigenvector derivatives and high order eigenvectors.

For the two first orders, we have:

s (1) = np k 1 =1 ∂s (0) ∂p k 1 (p (m) k 1 -p (0) k 1 ) (4.11) s (2) = 1 2! np k 1 =1 np k 2 =1 [ ∂ 2 s (0) ∂p k 1 p k 2 (p (m) k 1 -p (0) k 1 )(p (m) k 2 -p (0) k 2 )] (4.12)
These data are integrated in the projection matrix T z as explained in Eq. 4.10. As a reminder, the maximal size of the projection matrix is [n dof × (n p × n + 1)] whereas the reduced problem is

[(n p × n + 1) × (n p × n + 1)].

Coupling of high order derivatives with kriging

In this section, the same strategy, already used with success in sections 3.2 of chapter 2 and 3.4 of chapter 3, is considered for the case of eigenderivatives. The space of parameters and the space of participation factors are linked with a kriging technique to calculate q z (p) for a new sample p. In the present case, the projection basis T can be composed either by the derivatives of eigenvectors (Section 4.2), either by the high order perturbed eigenvectors defined with eigenvectors derivatives (Section 4.3).

5 Numerical applications: buckling load evaluation of a modified structure

Finite element matrices approximation

Before evaluating the efficiency of the different approximation methods, this section is dedicated to the validation of the stiffness matrices approximation. In the present application, four parameters are considered as variable, namely the Young's modulus, the angle value and the two dimensions of the structure section.

For the Young's modulus, the order of truncation is naturally fixed to one whereas a third or fourth order is chosen for the other geometric parameters. The same numerical framework is considered concerning the random snapshots. As the determination of regression coefficients relies on the knowledge of reference solutions of a set of snapshot, the maximal errors on buckling loads are drawn in Fig. 4.1, respectively for a third and fourth order of truncation for the geometric parameters. The approximated solutions are calculated with the parametric stiffness matrices and a standard Sorensen algorithm. For a third order, we observe that the maximal errors are always superior to the threshold of 1% whatever the number of snapshots used. On the contrary, the maximal errors are inferior to 1% for a fourth order as soon as the number of snapshots is superior to 60. These results confirm the non linear behavior of the critical load as a function of geometric parameters of the structure.

Considering this decomposition of matrices, the next sections are dedicated to test the efficiency of the different methods based on eigensolutions derivatives.

Case 1: Taylor expansion series

The approximation of buckling load factor with Taylor expansion series is dependent on the order of truncation n, as detailed in Section 3. To study the effect of this parameter, the maximal errors on the buckling loads are exposed in Fig. 4.2 for a first, second and third order. First of all, we observe a decrease of the error for a second and a third order in comparison to the first one. Nevertheless, the errors remain high with a level of error close to 10%. The associated modes are badly calculated too with a norm error superior to 250% and some MAC inferior to 0.2. The convergence radius of this approximation is clearly too low as soon as the bounds of input parameters are considered (Fig. 4.3). The behavior is globally comparable to the series developments, exposed in chapter 3.

Even if the computational time is very low (only some matricial operations), this kind of approximation is not compatible with a multiparametric analysis for buckling. Whatever the order of approximation is chosen, the maximal errors on buckling loads are inferior to 1% for the two methods. The quality of approximation is high for both critical loads and associated modes for the two projections with only a first order, as shown in Figs. 4.5 and 4.6. Moreover, the norm error for modes is inferior to 1% with a MAC equal to 1. About the computational time, the gain factors can be rather interesting. After calculating the approximation of matrices with n s snapshots and the derivatives for a nominal set of parameters during the offline step, the online step relies, for each new input set, on the assembling of reduced matrices (already stored in disk during the offline step) and on the calculation of a reduced buckling problem. Thus, the more the number of degrees of freedom increases, the more the gain factor increases too. The effective gain factor will be presented in Section 6. As previously, a first order is sufficient to obtain a maximal error inferior to 1% as described in Fig. 4.7. However, we observe that ROM-Taylor-HPP-Kriging is not enough robust as soon as the order of truncation increases.

About the sample data size, the use of thirty snapshots (Fig. 4.8) allows to guarantee an acceptable level of error for the buckling loads (Fig. 4.9) in the case of ROM-Taylor-Kriging and ROM-Taylor-HPP-Kriging methods. Nevertheless, the evolution of ROM-Taylor-HPP-Kriging is less stable than ROM-Taylor-Kriging as a function of the sample data size.

For the buckling modes, the assessment is the same as the previous section for ROM-Taylor-Kriging method because the norm error is inferior to 2% with a MAC of 0.99. In conclusion of this section, ROM-Taylor-Kriging method seems to be more efficient than ROM-Taylor-HPP-Kriging method in term of global behavior as a function of key parameters (order of truncation, sample data size). The updating of the projection basis and the calculation of participation factor for each new sample do not guarantee a good behavior.

Assessment of results

In this section, we have highlighted that methods, based on derivatives, can be a possible alternative to the reference method to perform a multiparametric buckling analysis. The main condition for the success is to built a suitable approximation of finite element matrices compatible with the behavior of each input parameter.

Except for a direct approximation with Taylor expansion series, the other proposals supply efficient results as summarized in Fig. 4.10. Indeed, the level of precision verifies the threshold fixed at the beginning of this benchmark both for buckling loads and associated modes. Only the ROM-Taylor-HPP-Kriging method exceeds the threshold for the largest finite element model. Moreover, the computational gain is of the order of 36 for ROM-Taylor, ROM-Taylor-HPP, ROM-Taylor-HPP-Kriging and 39 for ROM-Taylor-Kriging. Indeed, as soon as the parametric description of finite element matrices are available, the calculation of the other steps are not clearly time consuming. The integration of Kriging interpolation do not supply a additional reduction of computational time as in previous chapters. For the following chapter, ROM-Taylor-HPP-Kriging method will not select because we have shown that the behavior is not always stable. The level of precision is globally comparable for the three methods. The maximal errors of the critical loads are always inferior to the threshold of 1% with a MAC close to 1 for the associated eigenvectors. About the computational time, ROM-POD-Kriging and ROM-HPP-Kriging supply the best gain factor in comparison to ROM-Taylor-Kriging. Indeed, the number of snapshots necessary to build the approximation of matrices is superior to the two other ones. The number of snapshots is respectively equal to 35, 25 and 60 for ROM-POD-Kriging, ROM-HPP-Kriging and ROM-Taylor-Kriging methods in Table 5.1. The CPU gain factor will be more important as soon as the number of degrees of freedom increases.

In conclusion, ROM-HPP-Kriging method, summarized in Fig. 5.2, is the best compromise to be used in multiparametric analyses. Part II

Application on lattice structures 1 Introduction

Additive Manufacturing (AM), also commonly known as 3D printing, is a manufacturing process for building almost any physical solid part from its three-dimensional digital model. In AM, 3D complex parts are built-up by direct deposition of successive controlled thin layers of melted metal on top of each other.

In its early days, AM was mainly used for rapid prototyping [START_REF] Wohlers | Wohlers Report 2014: 3D Printing and Additive Manufacturing State of the Industry Annual Worldwide Progress Report[END_REF], but quickly manufacturers discovered the potential of this new manufacturing process. AM was mostly plebiscited (and still) by industrials of strategic sectors such as defense industry, aerospace and aeronautics mainly for lightweight capabilities of the produced parts. However, in the recent years, numerous industrial sectors such as automotive, medical, building,... have started to use AM, because it enables the creation of innovative shapes which were not feasible before, using the standard manufacturing processes. Some of the most used materials for these demanding applications, such as Inconel 718 and Titanium Ti6Al4V, can be printed in 3D.

Nowadays, AM technology has a very fast growing market capability mostly in the aerospace, automotive, healthcare, and consumer products sectors. According to the annual report on Additive Manufacturing Technology by Wohlers [START_REF] Wohlers | Wohlers Report 2016: 3D Printing and Additive Manufacturing State of the Industry : Annual Worldwide Progress Report[END_REF], the AM industry grew by the most interesting, mainly for their lightweight characteristics and energy absorption in specific applications such as impact protection. [Gordon, 2008] studied the mechanical properties of lattice structures under compressive loading. The author used analytical models to predict the properties of the strength based on the critical inelastic buckling stress of an elementary cell of the whole lattice. [START_REF] Campanelli | Manufacturing and characterization of ti6al4v lattice components manufactured by selective laser melting[END_REF] investigated titanium alloy Ti6Al4V lattice structures for the production of lightweight components. Compression tests were carried out in order to evaluate the mechanical strength and the energy absorbed. [START_REF] Choy | Compressive properties of ti-6al-4v lattice structures fabricated by selective laser melting: Design, orientation and density[END_REF] and [START_REF] Sing | Selective laser melting of lattice structures: A statistical approach to manufacturability and mechanical behavior[END_REF] studied Ti6Al4V lattice structures made of different unit cell geometries which were fabricated by AM. Designs were compared in different geometrical orientations and different densities for investigation of deformation behavior and compressive properties (Figure 6.4).

Figure 6.4 -Ti6Al4V made lattice structures produced by AM and their behavior was studied under compression [START_REF] Choy | Compressive properties of ti-6al-4v lattice structures fabricated by selective laser melting: Design, orientation and density[END_REF], [START_REF] Sing | Selective laser melting of lattice structures: A statistical approach to manufacturability and mechanical behavior[END_REF].

[ [START_REF] Sercombe | Failure modes in high strength and stiffness to weight scaffolds produced by selective laser melting[END_REF] studied the failure mechanisms of porous scaffold structures obtained by AM. Lattice structures made of Ti6Al4V were tested under interrupted compression and analyzed by X-Ray Micro Tomography to study the deformation and failure of the scaffolds. Authors concluded that failure was observed in areas that exhibit the greatest tensile stress, while the onset of the commonly observed layered failure occurred afterwards. Paulose et al. [START_REF] Paulose | Selective buckling via states of self-stress in topological metamaterials[END_REF] investigated the design of 3D printed lattice structures to test a new mechanical material that can be manipulated at strategic points for possible use as an insulator for frames as a start. Authors were able to introduce controlled imperfections in some areas of the 3D lattice to let them buckle in a controlled manner allowing the deformed structure to morph into a specific shape due to an increase in temperature or stresses (see Figure 6.5). The authors demonstrate the effectiveness of the topological design through analytical and numerical calculations as well as buckling experiments performed on two and three-dimensional metamaterials built out of stacked kagome lattices.

In 2016 [START_REF] Kaur | Three-dimensionally printed cellular architecture materials: perspectives on fabrication, material advances, and applications[END_REF] made a literature survey about recent developments in AM of smallscale 3D lattice structures with dimensions that span to nanometers. They highlight the advantage of well-defined architecture and small-scale size effect towards a unique class of materials with supreme mechanical strength at ultra-low densities.

More recently [START_REF] Chen | 3d printed hierarchical honeycombs with shape integrity under large compressive deformations[END_REF] investigated the compressive performance of a hierarchical lattice structure created by replacing cell walls in regular honeycombs with triangular lattice configurations. They found that hierarchical lattice honeycombs exhibit a progressive failure mode, along with improved stiffness and energy absorption under uniaxial compression. They reported also that high energy dissipation and shape integrity was noticed at large imposed strains. One can observe clearly from Figure 6.6, that lattice honeycombs with different densities undergo at first a linear buckling with a clearly identified first critical load before collapsing under higher loading pressure (Figure 6.6).

Therefore, it appears that a linear buckling analysis could be an important task as a first attempt, to characterize the first critical buckling load, at the preliminary design stage of such complex 3D lattice structures obtained by AM. This matter will be the subject of the present chapter, where two lattice structures (see Figure 6.7) made by AM and supplied by the G-SCOP laboratory (www.g-scop.grenoble-inp.fr) at Grenoble Alpes University will be studied and the variability of some geometrical and material parameters will be investigated.

Brief description of additive manufacturing process

AM process begins with a CAD model exported to the AM machine in a STL file format.

The AM machine slices the CAD model into a finite number of layers, according to the desired precision. Nowadays, a large variety of AM processes exists and unfortunately still no standardization norm in this fabrication process contrarily to conventional manufacturing processes. The -Material Extrusion (ME) : where beads of thermoplastic filament are molten using a controlled heating power, then extruded and selectively distributed through a nozzle. The technique used here is the Fused Deposition Modeling, and used material are limited to polymers and thermoplastics.

The 3D lattice structures (Figure 6.7) studied in this thesis were obtained using the Electron Beam Melting technique, therefore we will limit description only to this AM process.

Electron beam melting (EBM)

Electron beam melting (EBM) is an advanced AM process in which metal powder bed is fused by a concentrated beam of electrons heated using a controlled energy [Murr, 2015]. Production in a vacuum chamber allows that electrons don't collide with gas molecules and also prevent oxidation which can compromise reactive materials such as titanium [START_REF] Galati | A literature review of powder-based electron beam melting focusing on numerical simulations[END_REF].

EBM is generally used to produce components for several industrial sectors such as aerospace, automotive and medical industry.

The principle of EBM process is presented in its simplest scheme in Figure 6.8. Once a part has been produced, the build wrapper is removed together with other attached objects from the loose powder. Powder clinging to the part and remaining in internal cavities is cleaned and blasted away. Post-processing techniques, like hot isostatic pressing may be employed to release residual stresses and improve mechanical properties. In some cases, machining may be used to achieve required tolerances. Electro-polishing is also often used to improve the final surface quality of the part.

Compatible materials with EBM process

Strategic industrial sectors such as defense and aerospace have been for a very long time attracted by lightweight materials for their higher strength and corrosion resistance such as Titanium [START_REF] Campanelli | Manufacturing and characterization of ti6al4v lattice components manufactured by selective laser melting[END_REF]. It has been reported from several research papers such as [START_REF] Castellanos | A comparative study of manufacturing processes of complex surface parts in titanium ti6al4v[END_REF], [START_REF] Galati | A literature review of powder-based electron beam melting focusing on numerical simulations[END_REF] that lightweight Ti6Al4V parts obtained by AM have high performances in strength and durability, a serious benefit in applications where Figure 6.8 -Principle of Electron beam melting on an Arcam machine [START_REF] Galati | A literature review of powder-based electron beam melting focusing on numerical simulations[END_REF].

weight saving is critical.

In medical industry, biocompatibility of titanium alloy [Sidambe, 2014] permits the production of sophisticated healthcare solutions by means of AM process such as prosthesis, orthopedics and dental implants made of Ti6Al4V.

The mechanical properties of the standard titanium alloy parts made using the EBM process are summarized in (Table 6.1). The Arcam Titanium Ti6Al4V (Grade 5) [Arcam, 2017b] powder has a particle size between 45 and 100 microns. This limit on the minimum particle size ensures safe handling of the powder.

Ti6Al4V ELI (Grade 23) [Arcam, 2017a] is very similar to its parent Ti6Al4V (Grade 5), except that it contains reduced levels of oxygen, nitrogen, carbon and iron. ELI is short for Extra Low Interstitials, and these lower interstitials provide improved ductility and better fracture toughness for the Ti6Al4V ELI material. It is typically used for biomedical implants, cryogenic applications and offshore equipments. The Arcam Titanium Ti6Al4V ELI (Grade 23) powder particle size is the same as for Ti6Al4V. The mechanical properties of the Ti6Al4V ELI parts made using the EBM process are summarized in (Table 6.2)

Other super-alloys can be processed with EBM, such as cobalt chrome [START_REF] Gora | Enhancing surface finish of additively manufactured titanium and cobalt chrome elements using laser based finishing[END_REF] Table 6.1 -Mechanical properties of Ti6Al4V [START_REF] Arcam ; Arcam | Ti6Al4V Titanium Alloy[END_REF].

Table 6.2 -Mechanical properties of Ti6Al4V ELI [Arcam, 2017a].

with its excellent mechanical characteristics in terms of strength and temperature resistance has attracted industrials in the defense, the aerospace where parts operate at very high temperatures. The 316L stainless steel powder is often used to produce parts by means of AM [START_REF] Zhong | Additive manufacturing of 316l stainless steel by electron beam melting for nuclear fusion applications[END_REF]. Obtained parts are known for their excellent corrosion resistance, which makes them attractive for automotive and medical industries.

Inconel 718 is another superalloy which can be processed by EBM [START_REF] Popovich | Functionally graded inconel 718 processed by additive manufacturing: Crystallographic texture, anisotropy of microstructure and mechanical properties[END_REF] because of its excellent mechanical properties, corrosion resistance and top performance at high temperatures. EBM-printed parts fabricated from Inconel 718 are used in racing and aerospace applications. Valves fabricated from this nickel alloy are used in the petrochemical industry where durability and corrosion resistance are important characteristics.

the build orientation and specifying a minimum cross-sectional thickness were ways to reduce the effect of the anisotropy and heterogeneity in material properties. Earlier, [START_REF] Simonelli | Effect of the build orientation on the mechanical properties and fracture modes of slm ti6al4v[END_REF] investigated the tensile properties of Ti6Al4V parts obtained by AM. They conclude that ductility is dependent on the build orientation of the parts. In particular, horizontally orientated AM parts typically showed a higher mechanical strength as compared to vertically orientated ones.

Similar results were reported by [START_REF] Carroll | Anisotropic tensile behavior of ti6al4v components fabricated with directed energy deposition additive manufacturing[END_REF] on the influence of the anisotropic tensile behavior of Ti6Al4V parts obtained by AM. Generally in EBM process, the final surface quality has to be processed because in addition to the staircase irregularities, the EBM process also faces other issues. When the powder is melted, heat exchanges cause some of the bed powder to be partially sintered and sticked to the part surface. The surface sintering effect not only degrades the geometrical tolerancing accuracy of the part, but also creates surface defects which constitutes crack initiation sites. Therefore, it appears that the surface quality improvement of EBM produced parts is of critical importance.

Within the EBM process, numerous post-process surface treatment techniques exist, to improve the accuracy of the final parts. The most used are : machining, mechanical-polishing, abrasivepolishing, electro-polishing, chemical milling and electroplating. Figure 6.11 shows the EBM manufacturing steps of cylindrical lattice structure which will be analyzed in our work.

Advantages and drawbacks of EBM

The EBM process has several advantages and also some drawbacks, as with every manufacturing process.

Advantages

Among the advantages of EBM, we can cite : -EBM process is more competitive economically than other processes.

-EBM process offers more flexibility for the design of complex parts.

-EBM produces parts are of better quality than those obtained by casting.

-EBM parts request less supports during forming than other processes.

-EBM process is 95 % energy-efficient, 5∼10 times better than Selective Laser Melting (SLM)

process.

-EBM parts have less residual stresses than in other processes.

Drawbacks

EBM process has some drawbacks, among :

-EBM parts have surfaces with low quality, it generally requires post-processing.

-EBM process needs a significant amount of calibration.

-EBM requires the use of pure powder metals.

-EBM machines require important preventative maintenance.

-EBM parts quality depends on the position of the part during the forming.

-EBM parts generally have a lower strength and a lower durability compared to SLM parts.

3 Description of the studied 3D printed lattice structures

Lattice structures for energy absorption

As seen in the previous sections, lattice structures are one of the best choices when it comes to energy absorption during impact. Several research groups have investigated the influence of different patterns on cellular structures performances. [START_REF] Labeas | Investigation on the static response and failure process of metallic open lattice cellular structures[END_REF] studied the response of three different cellular core types, suitable for manufacturing crashworthiness sandwich cellular structures. They showed that the structural response is strongly influenced by the aspect ratio (radius/length), as well as the unit-cell size and shape. [START_REF] Yang | Mechanical properties of 3d re-entrant honeycomb auxetic structures realized via additive manufacturing[END_REF] investigated experimentally and analytically the behavior of a 3D re-entrant honeycomb auxetic cellular structure shown in Figure 6.12. They found that re-entrant auxetic structures have relatively little sensitivity to the number of unit cells.

(a) (b) Figure 6.12 -3D re-entrant honeycomb auxetic cellular structure a) Re-entrant auxetic samples by EBM. b) Auxetic unit cell [START_REF] Yang | Mechanical properties of 3d re-entrant honeycomb auxetic structures realized via additive manufacturing[END_REF].

[ [START_REF] Warmuth | Fabrication and characterisation of a fully auxetic 3d lattice structure via selective electron beam melting[END_REF] investigated the three-dimensional chiral cellular structure of Figure 6.13 which exhibits auxetic behavior in all three directions caused by negative Poissons ratio.

Samples were fabricated from Ti6Al4V powder via selective electron beam melting. The authors studied the influence of the strut thickness on the mechanical properties and the deformation behavior of cellular structures. They showed that the deformation changes from stretching to bending dominated by an appropriate choice of the amplitude of the strut.

[ [START_REF] Ozdemir | Energy absorption in lattice structures in dynamics: Nonlinear fe simulations[END_REF] undertook an experimental study of Ti6Al4V re-entrant and diamond lattice structures across a range of loading rates in Figure 6.14. The authors found that, whilst re-entrant cube specimens made up of multiple layers of unit cells were load rate sensitive, the mechanical properties of lattice structure cell layers were insensitive to load rate. This is an important characteristic since the later type of lattice structure could be used as an energy absorber in impact problems.

(a) (b) Figure 6.14 -Representative unit cells of (a) diamond and (b) re-entrant cube lattice structures [START_REF] Ozdemir | Energy absorption in lattice structures in dynamics: Nonlinear fe simulations[END_REF].

With AM process one can theoretically produce an infinite variety of lattice structures using different patterns within a design space. It is now possible to find through the literature, a wide variety of lattice shapes (Figure 6.15) having more or less complex patterns [START_REF] Ahmadi | Additively manufactured open-cell porous biomaterials made from six different space-filling unit cells: The mechanical and morphological properties[END_REF].

They have all in common a lightweight and high energy absorption in impacts compared to classical parts.

In the present thesis work, two lattice structures made of Ti6Al4V ELI and supplied by the G-SCOP laboratory at Grenoble Alpes University will be studied and the variability of some geometrical and material parameters will be investigated. The two lattice structures are designed as new solutions for the system of shock absorber for aircraft seats as shown in Figure 6.16.

Octagonal lattice structure

The first lattice structure is 3D truss made of titanium alloy Ti6Al4V ELI. The lattice structure has an octagonal shape and is made of four-levels and four symmetrically linked branches in the horizontal plane. Levels are linked together by branches of equal height L 2 . The cross section size of the branches is circular of the same diameter as depicted in Figure 6.17.

(a) (b) (c) (d) (e) (f)
All dimensions of the lattice structure are given in the 2D layout of the CAD geometry given in Figure 6.18. The mechanical properties of Ti6Al4V ELI Titanium Alloy, may be found in [Arcam, 2017a] and a wide range of values may be found in the literature. The main geometric and material characteristics of the first studied octagonal lattice structure, are summarized in Table 6.3.

Re-entrant cube lattice structure

The second lattice structure is 3D truss made also of titanium alloy Ti6Al4V ELI. The lattice structure has a cubic shape and is generated by duplicating three re-entrant cube unit cells along three axes, as shown in Figure 6.19.

Parameter Value Unit

Length L 1 12. The re-entrant cube unit cell has dimensions a×a×a and all branches are of the same circular cross section with a constant diameter d as depicted in Figure 6.20. The geometric and material characteristics of the second lattice structure, are summarized in Table 6 4 Uncertainty quantification and experimental buckling loads

Geometrical variability associated to EBM process

The geometrical variability of 3D lattices structures obtained by means of EBM process is well established. This is generally due to several parameters among, the built-part orientation, the layer thickness and the melted pool depth. Some authors reported that 3D printed lattice parts produced by EBM may present high variability in strut size and roughness [START_REF] Yang | Modeling of uniaxial compression in a 3d periodic re-entrant lattice structure[END_REF],

and they introduce the so-called effective strut diameter d which falls entirely within the strut, as shown in Figure 6.21. From the microscopic observations, [START_REF] Antonysamy | Effect of build geometry on the β-grain structure and texture in additive manufacture of ti6al4v by selective electron beam melting[END_REF] analyzed the effect of build geometry on the β-grain structure in Ti6Al4V parts, produced by EBM. They showed that nucleation of the β-grains always occurred heterogeneously from boundary layers. This is resulting from partially melted powder in the surrounding bed. The authors found that in bulk sections, coarse columnar grains develop a strong texture parallel to the build direction.

In contrast, the contour produced a complex skin grain structure which dominated thin sections less than 2mm thick. (which is the reference for our two lattice structures). From Figure 6.22, we can understand that powders of 80µm may partially increase the geometrical diameter of our parts, without necessarily increasing the mechanical properties. In addition, knowing that Arcam AB had announced that the roughness of parts produced by EBM is about 25 -30µm. By thus, accumulating this value with the error induced by fused powders is of 30% on diameter, an error of about 90µm can be obtained on the working radius of the strut.

[ [START_REF] Suard | Towards stiffness prediction of cellular structures made by electron beam melting (ebm)[END_REF], Suard, 2015] studied struts geometry based on X-ray micro-tomography.

They start from 3D images and develop a criterion of "mechanically efficient volume" for stiffness prediction. They showed that the effective volume ratio, defined as the ratio of volume of the inscribed cylinder over the volume of the strut estimated from 3D image analysis (6.23), leads to a slight underestimation of the stiffness. Thus, Suard proposes a correction of the CAD diameter d CAD as a function of the manufacturing orientation α considering a linear regression model.

The effective diameter d EF F is described as follows: with k 0 = 8.65 10 -1 , k 1 = 1.02 10 -3 , k 2 = -1.16 10 -2 , k 3 = -3.37 10 -3 the regression coefficients of the polynomial. The author precises that this identification leads to a maximal error between the measured values and the mathematical model of 5%. This mathematical model and associated identified errors will be used in the next chapter relative to uncertainty management.

d EF F = (k 0 + k 1 α)d CAD + (k 2 + k 3 α) (6.1) (a) (b)

Material variability associated to EBM process

Quasi-static traction tests are presented for 3 or 4 sets of raw specimens made of Ti6Al4V by ARCAM EBM machine with three orientations as shown in Figure 6.24. Variation of material 123

Experimental setup variability

In the present investigation, a SINTECH 20/D machine was used for the quasi-static measurements. Besides the standard experimental measurement errors which usually exist in all devices, it has to be noticed that the upper clamping claw of the machine was free to rotate before it stabilizes and starts pressing the specimen. Figure 6.25 shows the possible rotation of the upper clamping claw of the SINTECH machine. Therefore it has to be mentioned that the clamping claw may take an initial rotation angle which may have a maximal amplitude of β max = ±2 • . This value has to be taken into account as a variability parameter of the experimental setup itself. Moreover, one has to add to this variability parameter the standard error measurement on the load and displacements results.

Experimental buckling tests

Quasi-static experiments have been conducted to determine the buckling loads of the two studied lattice structures (octagonal and re-entrant cube lattice structures). Compression tests have been achieved using a SINTECH 20/D machine with a maximum loading capacity of 100 kN, available at LAMIH as shown in Figure 6.26. The specimen is placed between two lower and upper claws and the load was applied very slowly with a speed of 2mm/min in such a way that the structure will deform slowly at very low strain rate. The observation of the buckling process of the octagonal lattice structure during compression is shown in Figure 6.27. As it can be observed experimentally, the deformed lattice structure behaves as expected showing the first buckling mode shape consistent with the theory. The observation of the buckling process of the re-entrant cube lattice structure during com- The linear buckling formulation allows to obtain two matrices: a classical linear stiffness matrix and a geometric stiffness one. The former can be written as Figure 6.31 -A 3D beam element with two nodes and 6 dof per node.

K l = L 0 B T l H B l dx (6.2)
where L 0 is the initial length of the beam, H is the 6 × 6 generalized linear material law which depends for isotropic materials on two parameters E and ν. B l is the 6 × 12 constant finite element strain matrix.

Taking into account the membrane, bending and sheering forces, the geometric stiffness matrix reads:

K σ = L 0 B T g σ B g dx (6.3)
where σ is a generalized stress matrix, which is function of the generalized forces N x , T y , T z and moments M x , M y and M z . B g is a constant matrix defined in terms of the first derivatives of the shape functions. Only one quadrature Gauss point is necessary for full integration of both linear and geometric stiffness matrices.

In the following, the linear buckling problem will be solved as explained previously in chapter 1 using these two matrices. It has to be noticed that the choice of a beam element has been justified by the high accuracy which can be obtained using much less elements compared to the solid hexahedral elements.

FE modeling of the octagonal lattice structure

The octagonal lattice structure, shown in Figure 6.32, is defined with 2 strut types corresponding to two orientations 0 • and 90 • . Thus, two elastic modulus E 0 • and E 90 • and two diameters of section d 0 • and d 90 • are so considered for the parameterization. Hence, an inclination angle β is added as shown in Figure 6.32.

Figure 6.32 -Choice of variability parameters for the octagonal lattice structure.

The octagonal lattice structure is submitted to a total compression load F = 37kN which is distributed on 9 nodes at the top of the structure as indicated in Figure 6.33. At the bottom all nodes are considered as clamped. 6.6. These different data will be used in the next section for the correlation with experimental results.

The linear buckling of the structure has to be validated using CAD nominal. To this purpose, ABAQUS © software has been used as a reference solution where the B31 two-node beam model was adopted. For our modeling using MATLAB programming, the linear buckling analysis has been done by means of eigenvalue solver IRA/Sorensen which is valid for very large number of degrees of freedom. The first buckling mode shape is depicted in Figure 6.34.

Figure 6.34 -First buckling mode shape of the octagonal lattice structure using ABAQUS © .

For the computation of the critical buckling load, different finite element mesh refinements have been used and our solution obtained using MATLAB was compared to the one obtained using ABAQUS © software for exactly the same mesh. Figure 6.35 shows the comparison between the two solutions in function of the number of dof. Both behaviors of our solution and the one obtained using ABAQUS © , shown that the critical loads reduces and tends to be unchanged after a certain refinement. Considering the input values of Suard's case, the critical load changes from 24162N to 22253N for ABAQUS © , and it changes from 28908N to 22188N in our model when refining mesh from 819 dof to 57267 dof.

As it can be observed, the critical buckling load becomes stable after 5427 dof. The related accuracy of critical load between ABAQUS © and our model is 0.3% up to 5427 dof. Hence, we can consider that our developed 3D beam element is validated. As previously, the different description of parameters values are resumed in Table 6.7.

Parameter CAD description Mean description Suard description

Young's Modulus E 0 ABAQUS © and MATLAB buckling analyses were performed in the same conditions that in the previous case. After resolution, the first buckling mode shape was obtained and is depicted in Figure 6.38.

Figure 6.39 shows the comparison between the two solutions in function of the number of dof.

Both behaviors of our solution and the one obtained using ABAQUS © , shown that the critical load reduces and tends to be unchanged after a certain refinement. Considering the input values of Suard's case, the critical load changes from 22486N to 20252N for ABAQUS © , and it changes from 24825N to 19947N in our model when refining mesh from 5196 dof to 75180 dof.

As it can be observed, the critical buckling load becomes stable after 36300 dof. The related accuracy of critical load between ABAQUS © and our model is 1.2% up to 36300 dof.

Deterministic numerical and experimental comparisons

The aim of this section is to compare numerical deterministic results of critical buckling loads to those experimentally measured. For recall, three descriptions are here considered corresponding to different values of diameters as a function of the orientation. Tables 6.8 and 6.9 summarize buckling loads respectively obtained for the octagonal and re-entrant cube lattice structures.

Conclusions

The main objective of this chapter was to propose an assessment of the variability associated to the additive manufacturing of the lattice structures and to perform a first comparison between numerical and experimental buckling loads.

After the brief description of the EBM method, two studied lattice structures with octagonal and re-entrant cube shapes, are presented. Different experiments have been performed to highlight the variations on the buckling loads and the elastic modulus as a function of the manufacturing orientations. Next, we have proposed to investigate the Suardś works proposing a correction of CAD diameters to consider the effective active volume of the built structures. Finally, two finite element models have been developed and validated thanks to an industrial software to study the buckling loads of the lattice structures.

The comparison of numerical and experimental results have highlighted that the correction, proposed by Suard, is a good starting point to simulate the behavior of our lattice structures.

However, to numerically capture the variable experimental behaviors, the integration of uncertainty in the numerical model seems to be interesting. This way will be developed in the next chapter by considering both the different description to model the diameter variability and the uncertainty associated to elastic modulus and angle inclination.

Chapter 7

Additive manufacturing of lattice structures: Uncertainty propagation In spite of the developments of efficient numerical simulations, the correlation between a deterministic simulation and experimentations are not obvious in many cases. Indeed, during the manufacturing of mechanical structures, it is not uncommon to observe some uncertainties resulting in product variability either on material properties, on geometric characteristics and on boundary conditions. These observed variability necessarily affect the response of structures and the efficiency of the studied system. Numerically, to take into account these uncertainties and tend to reliable and robust predictions, a current industrial trend involves making multiple numerical simulations by performing sensitivity analyses, designs of experiments [Fisher, 1935, Taguchi et al., 2005, Kleijnen, 2010] non-deterministic studies [Moore, 1966, Shinozuka, 1972, Zadeh, 1978, Ben-Haim and Elishakoff, 1990] or even reliable and robust optimizations [Goldberg, 1989, Kennedy andEberhart, 1995]. Thus, the idea is to simulate the evolution of mechanical responses as a function of input parameter variations and to detect failures and performance reductions of products.

In mechanical engineering, several theories, such as probabilistic theory [START_REF] Metropolis | The monte carlo method[END_REF], interval theory [Moore, 1966], and fuzzy set theory [Zadeh, 1965],

have already been used. Whatever the used theory, the introduction of uncertainty in numerical simulations relies on several steps, namely uncertainty modeling step, uncertainty propagation step and finally uncertainty management step.

Two major axes have generally been chosen considering intrusive or non-intrusive propagation methods. With intrusive approaches, the uncertain parameters are replaced in the governing equations by appropriate expansions or approximations. To implement this category of approach in a classical simulation code, it is necessary to rewrite a large part of the software. On the contrary, the main objective of non-intrusive methods is to obtain the perturbed output quantities without making any modifications to the deterministic software. This second approach treats the deterministic code as a black-box and approximates the unknown coefficients with formulas based on deterministic code evaluations. The strategy for the selection of samples in the design space and their number depend on the chosen non-intrusive technique. Several intrusive or non-intrusive alternative propagation methods have been developed in order to overcome inconveniences of reference intensive samplings methods. In probabilistic approach, Monte Carlo Simulations (MCS) [START_REF] Metropolis | The monte carlo method[END_REF], represent the numerical reference method to propagate uncertainty.

This method is based on random sampling of input parameters and on deterministic calcula-problems to avoid the significant difficulties of implementation for non-linear problems. The second way, non-intrusive, relies on either regression methods or spectral projection methods. The first category calculates the stochastic modes by minimizing the gap between random function and its approximation on the polynomial chaos basis. On the contrary, the Non-Intrusive Spectral Projection (NISP) methods [START_REF] Nechak | Prediction of Random Self Frictioninduced Vibrations in Uncertain Dry Friction Systems Using a Multi-Element Generalized Polynomial Chaos Approach[END_REF], Nechak et al., 2013a] project the solution of problem on the basis of chaos. The initial problem is transformed into a calculation of integrals [START_REF] Crestaux | Polynominal chaos expansion for sensitivity analysis[END_REF] by the collocation approach for example [START_REF] Babuska | A stochastic collocation method for elliptic partial differential equations with random input data[END_REF], Loeven et al., 2007]. To avoid some limitations encountered in the case of large numbers of random variables or large range of variations, several alternatives have been proposed such as sparse generalized polynomial chaos [START_REF] Xiu | High-order collocation methods for differential equations with random inputs[END_REF], Ganapathysubramanian and Zabaras, 2007, Blatman and Sudret, 2008], multi-element generalized polynomial chaos [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF], Wiener-Haar chaos [START_REF] Maître | Uncertainty propagation using wiener-haar expansions[END_REF], Nechak et al., 2013b] or Wiener-Fourier chaos [START_REF] Millman | Airfoil pitch-and-plunge bifurcation behavior with fourier chaos expansions[END_REF]. As the parametric probabilistic approach are not sufficient to take modeling errors and modeling uncertainties into account, a second category, called the non-parametric methods [Soize, 2000] has been developed

and are based on the use of a reduced-order model and the random matrix theory. The random matrix theory [Mehta, 1991] is used to construct the prior probability distribution of the random matrices modeling the uncertain operators of the mean computational model. On the same way, the non-parametric method consists in constructing directly the stochastic modeling of the operators of the mean computational model. In practice, the generalized matrices issued from the nominal finite element model of the structure are replaced by random matrices. The uncertainties are introduced directly on global matrices of the model by using a dispersion parameter and its prior probability distribution. The dispersion parameter is estimated according to the Maximum Entropy Principle [Jaynes, 1957] for which the constraints are defined by the available information [Soize, 2005a[START_REF] Soize ; Soize | Random matrix theory for modeling uncertainties in computational mechanics[END_REF]. Many works have been published in order to validate the non-parametric probabilistic approach with experimental results [Chebli andSoize, 2004, Chen et al., 2006] and to extend the applicability of the theory to different areas in structural mechanics [START_REF] Capiez-Lernout | Nonparametric modeling of random uncertainties for dynamic response of mistuned bladed disks[END_REF]. Recently, an improvement of the non-parametric approach, called the generalized probabilistic approach of uncertainties, has been proposed [Soize, 2010] and allows the prior stochastic model of each type of uncertainties (uncertainty on the parameters of the model and modeling errors) to be separately constructed.

Other developments are dedicated to quantify the role of each kind of uncertainty (parametric and non-parametric) and to integrate them simultaneously in uncertain applications.

The numerical developments proposed in interval and fuzzy set theories are very close and can be explained together since a fuzzy problem can be decomposed into a set of interval prob-lems [START_REF] Massa | A fuzzy procedure for the static design of imprecise structures[END_REF], Massa et al., 2008a[START_REF] Massa | Multi-level homotopy perturbation and projection techniques for the reanalysis of quadratic eigenvalue problems: The application of stability analysis[END_REF]. The reference method used to solve problems described by the fuzzy formalism is based on Zadeh's Extension Principle (ZEP) [Zadeh, 1975a[START_REF] Zadeh ; Zadeh | The concept of a linguistic variable and its application to approximate reasoning-2[END_REF][START_REF] Zadeh ; Zadeh | The concept of a linguistic variable and its application to approximate reasoning-3[END_REF], which extends general operations for real numbers to the corresponding operations for fuzzy numbers. The Extension Principle states that the degree of membership of one combination is equal to the smallest degree of membership of the independent parameters of this combination. In the case of multiple occurrences of a solution, the final membership degree is equal to the maximum membership degree of the different solutions. In practice, fuzzy numbers can be discretized either according to the support (membership value equal to 0) or, more conveniently, according to the membership function. Formally, implementation of the method leads to a combinatorial process which can be compared to MCS in the sense that it presents the same drawback in terms of computational time. In a similar way as with the probabilistic theory, many researches have been performed to propose alternatives considering successively interval arithmetic [Moore, 1966], perturbation or series expansion [Qiu andElishakoff, 1998, Lallemand et al., 1999] and specific samplings or designs of experiments [Hanss, 2002, Hanss, 2003, Donders et al., 2005, Hanss and Turrin, 2010, Haag et al., 2010]. However, these methods have rapidly shown some limitations such as the overestimation of the output set with interval arithmetic or a bad approximation of non-linear problems with perturbation techniques. The best results have been obtained using either sensibility analyses or optimization algorithms. The first works have been initiated by McWilliam [McWilliam, 2001] and relied on monotonic evolutions of response surfaces. Bounds of solutions are determined with bounds of the intervals associated to each parameter. More recently, an interesting solution has been proposed by transforming interval problems as a min-max optimization problem to limit computational time of the uncertainty propagation step and thus to be compatible with an industrial context [Massa et al., 2009a]. Indeed, this approach relies on iterative search for global optima (local gradient algorithm) in a search space, the size of which is gradually increasing when the degree of membership is decreasing [START_REF] Massa | A fuzzy procedure for the static design of imprecise structures[END_REF], Massa et al., 2008b, Degrauwe et al., 2009]. For the crisp values, the output quantities and their first sensitivities are determined for each fuzzy parameter. The signs of the first-order sensitivities indicate the functional dependence of the response function and define the combinations of discrete fuzzy parameter values for the following level of degree of membership, which then can supply the minimum and maximum variations. For each α-cut level, the first derivatives of the output quantities are evaluated for the combinations of discrete fuzzy parameter values determined at the previous level. The comparison of the signs of the derivatives with those obtained at the previous level indicates the search path. In this context, the uncertainty propagation have already been studied in different linear problems such as static [START_REF] Massa | A fuzzy procedure for the static design of imprecise structures[END_REF],

modal [Massa et al., 2008b], frequency responses [Ruffin-Mourier, 2008, Massa et al., 2009b] and transient [START_REF] Rao | Transient response of structures with uncertain structural parameters[END_REF]. Moreover, to control the computational time, the proposed methods have been coupled with different approximations or reanalyses to limit the number of calls to finite element calculations.

The last step is the uncertainty management. After having calculated the non-deterministic solutions, it is possible to integrate them either on sensitivity analyses, to identify the most sensitive parameters or in optimization process, to study and improve the reliability and the robustness of the design. The aim of sensitivity analysis [START_REF] Saltelli | Structural reliability methods[END_REF], Cao et al., 2013, Zhang and Pandey, 2014] is to quantify the influence of each uncertain input parameter of a model. In this context, a global sensitivity is considered rather than a local deterministic sensitivity based on the computation of gradients of the response with respect to input parameters around a specific value. First, the aim is to take into account both the input interactions and the information for all possible ranges of variation of the uncertain parameters with the aim to estimate their influence [START_REF] Gauger | A new uncertainty analysis for the transformation method[END_REF]. Secondly, the idea is to avoid the drawback linked to the choice of a reference point where local gradients are generally calculated. In non-probabilistic context, a state-of-the-art of available methods is proposed in [START_REF] Saltelli | Structural reliability methods[END_REF]. Regressionbased methods or variance-based methods are mainly considered to measure the effects of inputs on outputs. In the case of variance methods, generally called ANalysis Of VAriance (ANOVA), the variances of the output are decomposed as a sum of contributions of each input variable.

Next, the reliability analysis [START_REF] Ditlevsen | Sensitivity Analysis[END_REF], Melchers, 1999, Lemaire, 2005] consists in computing the probability of failure of the system which can be described by a failure criterion. This mathematic criterion is directly linked to uncertain input parameters and mechanical output responses and is represented by a limit state function which separates safe and failure domain. The probability of failure is defined as the integral of the failure domain of the joint probability function of random vector and cannot be computed analytically. It is necessary to use numerical methods such as MCS or approximation methods such as First Order Reliability Method (FORM) or Second Order Reliability Method (SORM) [START_REF] Kiureghian | Second-order reliability approximations[END_REF] to evaluate the joint probability function. With the two last methods, the idea is to approximate the probability of failure and to limit the computational time compared to MCS. In current industrial applications, the objective is to optimize the system and to guaranty that 95%, even 99%, of the samplings stay in the safe domain. The recent developments are focussed on the integration of efficient approximation to define a more suitable approximated probabilistic problem. Finally, the robustness analysis traduces the ability of the optimum to be insensitive with respect to the input uncertainty. The aim is to guarantee in operation, low modifications of the proposed design behaviour. Robustness measures can be directly incorporates in the objective functions of non-probabilistic approaches of the non-deterministic optimization problem. For example, Kelesoglu [Kelesoglu, 2007] has proposed a method for fuzzy multi-objective optimization of space trusses, and Ramík [Ramík, 2007] has developed a class of fuzzy optimization problems with an objective function depending on fuzzy parameters. Massa et al [Massa et al., 2011a] have proposed a robust fuzzy optimization of a lower wall of micro satellite. The aim was to minimize the mass of the component while maximizing the robustness of the design parameters. In probabilistic theory, the robustness problem can be seen as a minimization problem of the variance of output data. Shah et al [START_REF] Shah | Multifidelity robust aerodynamic design optimization under mixed uncertainty[END_REF] have proposed a robust and computationally efficient optimization algorithm for air foil design under mixed uncertainty. The algorithm exploits stochastic expansions derived from the Non-Intrusive Polynomial Chaos techniques. Li et al [START_REF] Li | Multiobjective robust optimization for crashworthiness design of foam filled thin-walled structures with random and interval uncertainties[END_REF] proposed a multi-objective robust optimization procedure mixing both random and interval variables in crashworthiness design of a foam filled column. These advanced developments require very expensive computational times, which must be reduced, or at least maintained, using mathematical approximation as an alternative to the evaluation of large systems of equations (due to the finite element discretization for example).

Probabilistic approach

Many methods have been proposed in the literature for uncertainty propagation when uncertain parameters are represented by random variables. The following sections summarized the two methods used in this chapter for the numerical tests, namely Monte Carlo Simulation and polynomial chaos expansion.

In the probabilistic theory, uncertainty is described either by random variables or by random fields if spatial variation of the model properties is of concern. Construction of the associated probability distributions requires a statistical study of measured data from specific tests. This particular phase of uncertainty quantification is primordial to introduce realistic variations. When few experimental data are available from some observations of the system, these experimental data can be used to identify an optimal prior stochastic model of uncertainties using the maximum likelihood method. In mechanical engineering, Gaussian, log-normal and uniform laws are preferred in most applications.

Monte Carlo Simulation

Monte Carlo Simulation (MCS) [START_REF] Metropolis | The monte carlo method[END_REF] is a general mathematical technique which allows to propagate uncertainty on input parameters through a deterministic model.

MCS is a non-intrusive method and can be used in a wide variety of fields such as in quantitative finance, computational biology or physical science. The principle of the MCS is to simulate a large

A p,d = {α ∈ N p : |α| < d}cardA p,d = r = (p + d)! p!d! (7.9)
To limit the number of polynomial coefficients, one can consider the hyperbolic truncation scheme proposed by Blatman [Blatman, 2009].

This scheme is formally given by A p,d,q = {α ∈ A p,d : α q < d} (7.10)

where:

α q = p i=1 α q i 1 q (7.11)
For q=1, the hyperbolic truncation scheme correspond exactly to the standard one. For values inferior to 1, the number of interaction terms decrease as q decrease. An example of the behaviour of the hyperbolic truncation is given in Fig 7 .1 for a problem with 2 uncertain parameters and a maximum order of polynomial basis d=3.

Figure 7.1 -Hyperbolic truncation scheme for 2 uncertain parameters p1, p2 expand at degree 3.

Once the polynomial basis is created, the last step of the method concerns the evaluation of the coefficients y α of the PCE. Several methods have been proposed in the literature. They are either intrusive or non-intrusive. The non-intrusive ordinary least-square minimization method [START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF] has been chosen and implemented. The non-intrusive aspect and the fact that an arbitrary number of snapshots can be used to calculate the coefficients are the main advantages of the method.

Given an experimental design of N input random vectors X ed = {x (1) , • • • , x (N ) } T and their corresponding model responses Y ed = {y (1) , • • • , y (N ) } T , the vector y of coefficients y α is calculated with:

y = A T A A T Y ed (7.12)
where the experimental matrix A which contains the values of all basis polynomials in the experimental design points is given by:

A ij = Ψ x (i) with i = 1, • • • , N ; j = 1, • • • , r (7.13)
In practice, it is preferable that the size of the experimental design be superior to the rank r of the polynomial basis (note that rank 1 corresponds to degree 0). Of course, snapshots must be representative of the random input vector and are often evaluated using a Latin Hypercube Sampling.

With the coefficients in hand it is then straightforward to calculate the two first moments of the PCE and so of the outputs. The mean value of an output is given by:

µ Y = y 1 (7.14)
where y 1 is the coefficient of the constant basis term Ψ 1 = 1. The variance of an output is given by all over coefficients of the non-constant basis terms using:

σ 2 Y = α∈A,α =1 y 2 α (7.15)
Finally, one can obtain a posteriori error estimation of the PCE using the leave-one-out error [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF] which uses cross-validation. The technique consists in building N different PCE meta-models where one snapshot of the experimental design has been suppressed

and then verify the prediction of this snapshot with the real value. In practice, when the results of the least-square minimization are available, there is no need to explicitly calculate the PCE meta-models, which can be time consuming. In this case, leave-one-out error can be calculated using Eq. 7.16.

ε LOO = N i=1 M(x (i) )-M P CE (x (i) ) 1-h i 2 N i=1 (M (x (i) ) -µ Y ) 2 (7.16)
where h i is the i th component of the vector calculated from the experimental matrix defined by Eq 7.13. h = diag A A T A -1 A T (7.17)

Interval and fuzzy sets approaches

In interval theory [Moore, 1966], the uncertain data are represented by intervals of variation.

The inferior and superior bounds represent the minimal and maximal observed data and each value has the same confidence degree. Fuzzy set theory, introduced by Zadeh [Zadeh, 1978], can be seen as an extension of the classical set theory and a generalization of the interval theory. Zadeh's Extension Principle is outlined below:

Given a function ϕ that maps from X = X 1 × X 2 × . . . × X N to universe Y , such that y = ϕ(x 1 , x 2 , . . . , x n ) where y ∈ Y and x i ∈ X i , ∀i and considering fuzzy subsets A 1 , A 2 , . . . , A n signs, the response function cannot be considered monotonic, giving rise to an extremum between these two α-cut levels. The combination nearest the extremum is chosen and the search is stopped for this variation. Finally, the fuzzy output data are then built α-cut by α-cut.

To reduce the computational time associated with uncertainty propagation, we propose now to substitute the classical solver for stability problem (which can be solved several times in the Fuzzy numbers related diameters of struts are described using common triangular fuzzy numbers where typical values are defined as follows:

• The maximum of confidence corresponding to the crisp value (α = 1) is given to the diameter calculated from Suard's response surface model in Eq. 6.1 of Chapter 6.

• The lower bound of the support (α = 0) corresponds to 95% of the crisp value. In fact, a maximum gap of 5% was observed by Suard between the response surface and the measured values.

• The upper bound of the support is fixed to the CAD value expectation which cannot be attained, justifying a degree of confidence of 0.

Fuzzy numbers related to Young modulus of struts are described using common fuzzy intervals where typical values are defined as follows:

• Bounds of the interval at maximum of confidence are fixed with the extremum of measured Young moduli in Table 6.5 of Chapter 6

• A variability of 10% is generally adopted to quantify uncertainty on elastic modulus. This percentage is then shared between the lower and upper values of the previous interval which gives the bounds of the support.

A triangular fuzzy number is adopted for loading angle β with a crisp value at 90 • . A variation of 2 • around this value was estimated and allows to fix the bounds of the supports.

Tables 7.2 and 7.3 summarized the description of uncertain parameters for respectively octagonal lattice structure and re-entrant cube lattice structure. In both case, two orientations of struts are considered leading to two uncertain diameters and two uncertain Young moduli. Note that for reentrant cube lattice structure, a Young modulus measured for a 0 • manufacturing orientation has been applied for the 17.5 • struts orientation because material characterization was not realized for this orientation. A fifth parameter, namely the loading angle finish the parameterization of the uncertain buckling problem.

Definition of the Reduced Order Models

The ROM-HPP-Kriging method is used to approximate the critical loads by using the strategy summarized in Chapter 5. First, the projection matrix is defined considering the nominal data and the high order perturbed modes calculated per parameter direction. In the present case, five calculations are necessary to capture the evolution of the five variable parameters. For each homotopy development, a third order of truncation is selected. Thus, the number of modes in the projection matrix is equal to 16 to analyze the first buckling mode. About the kriging parameters, we consider a second order regression model and a Gaussian correlation function. The sample data size is equal to 50 for the first lattice structure and 50 for the second one. To highlight the efficiency of the ROM-HPP-Kriging method, we propose to compare it to kriging method as in the first part of this report. The results for the octogonal lattice are presented in Fig. 7.6 as a function of the number of degrees of freedom to test the robustness of the method too. First, kriging method does not allow to achieve the, certainly demanding, threshold fixed to 1% whatever the size of data samples used. Indeed, with only kriging, we cannot obtain an error less than 2.6% for 50 samples and less than 1.5% for 400 samples. Nevertheless, depending on the context, it is clear that these last values can be considered as acceptable. On the contrary, the maximal error for ROM-HPP-Kriging method is inferior to 1% with a DOE of 50 samples for all finite element models. Moreover, the computational gain for the proposed method is close to 70 in comparison to the reference calculation. Notice that the computational gain associated to kriging method is only to 7 as resumed in Table 7.4 for the case of 57267 dofs. The same efficiency is obtained for the second lattice structure. This particularity of the ROM-HPP-Kriging method will be highlighted again in sections 2.4 and 2.5 with the comparison of uncertain buckling loads.

Probabilistic simulations

The same procedure is applied for both lattice structures. In every cases, mean value and coefficient of variation of the buckling load as well as data relative to the numerical cost of the problem are stored. First, a MCS is performed using the full model to get the probabilistic reference results. To avoid any outlier value of the buckling load during the simulation, the current buckling mode is compared to the nominal one using a MAC criterion threshold of 70%.

Negative eigenvalues that may arise are also suppressed of the cumulative MCS process. For the studied problems, none of these cases occurred. Second, a MCS is performed using the ROM-HPP-Kriging described in section 2.3 and third, some PCE calculations was done in a way to find the parameterization giving the best results compared to the MCS reference ones. Relative errors on mean and CV of buckling load gives the first criterion to evaluate the quality of the different reduce models. They are completed with a criterion which gives the maximum error between buckling values given by the reference MCS and buckling approximations given by the reduce models. This last criterion, called verification, error is formalized by expression Eq. 7.22.

Err Y = Y i M CSref -M(X i M CSref ) Y i M CSref (7.22)
where M is either the ROM-HPP-Kriging or PCE model.

For a given number of uncertain parameter, the control parameters of the PCE model are the maximum order of polynomials basis (d) and the truncation scheme (q). These control parameters have a direct impact on the number of coefficients of the polynomial chaos and thus on the number of necessary snapshots of the experimental design. The snapshots themselves have also an impact on the accuracy of the PCE model because they are evaluated for random evaluations of the uncertain parameters using a LHS. As previously mentioned, the accuracy of the PCE can be a posteriori controlled using the leave-one-out error (ε LOO , Eq. 7.16) which behaviour is correlated with the verification error. However, note that in an effective use of the method, the verification error cannot be evaluated since it requires MCS reference. Thus, evaluation of the precision of the PCE will be performed using only the leave-one-out error in the sequel.

The nature relatively smooth of the studied problems (in terms of stochastic complexity) does not need high order polynomials. In a way to evaluate the minimal experimental design needed, different PCE configurations have been tested using either orders 2 or 3 with standard truncation and hyperbolic truncation with q=0.75 for order 3. With 5 uncertain parameters, the minimum size of the experimental designs needed to avoid underdetermined linear system in eq. 7.10 are N min = 21 for order 2, N min = 56 for order 3 (q=1) and N min = 26 for order 3 (q=0.75).

The Figs. 7.7 and 7.8 show the behaviour of the leave-one-out error thanks to the number of samples used to build the PCE for respectively the octagonal and the re-entrant cube lattice structures. Simulations was performed using a third order expansion and a standard truncation scheme starting from the minimum number of samples and ending by doubling this number. Note that a good evaluation of the leave-one-out error in the minimum case needs in fact N min + 1 samples (57 samples in these cases). For each sample set, 20 LHS have been performed to evaluate the influence of the random sampling. The verification error is also plotted for each case. On these figures, a blue dot corresponds to one PCE evaluation while the red cross gives the mean of the 20 evaluations per experimental design size. A comparison between the two figures show that it is not obvious to a priori set a threshold for the leave-one-out error since the minimal value that can be reached is structure dependent. For the selection of the optimal size of the experimental design, it is preferable to reason in terms of slope of the leave-one-out error. A second point that can be mentioned concerns the behaviour of the two errors which is well correlated. It can also be seen on both cases that the quality of the PCE are very sensitive to the random sampling when the number of samples is set to the minimum. From the observations of the mean values of the leave-one-out errors, one can consider that a number of 75 samples are necessary for the octagonal lattice structure while 70 samples can be selected for the re-entrant cube lattice structure in a way to get a good and stabilized quality of the PCE. In fact, for these selections good convergences of the leave-one-out curves are reached.

Following this last results, a similar study is performed considering a reduced order with a standard truncation scheme (order 2, q=1) and the same order with a hyperbolic truncation scheme (order 3, q=0.75). In both cases, sizes of experimental designs were fixed to a proportional value corresponding to the previous selection. These values correspond to an increase of the minimum experimental design size of 34% for the first structure and 25% for the second one. Mean results about the PCE studies are reported in Table 7.5 for the octagonal lattice structure and Table 7.6 for the re-entrant cube lattice structure. In both cases, one can notice that the two reduce basis cases does not allow to reach the same Loo error order than the previous simulation.

The verification errors follow the same rule. They are inferior to 1%, as expected, only in the cases corresponding to order 3 with a standard truncation scheme which confirms the choice of the PCE model parameterizations. a comparison between the MCS reference results performed with the full models, MCS results performed with the ROM-HPP-Kriging models and the chosen PCE models (order 3, q=1). As a measure of numerical cost, numbers of full model evaluations or their equivalent values are also reported in Table 7.7 for each method. Gains finally complete measures of efficiency.

Relative errors on the mean value and coefficient of variation of the first buckling factor calculated with both alternative methods are very small and of the same order. The same remark can be formulated for the verification error which reveal also a very good precision for the reduce

Fuzzy simulations

To perform the calculation with the Extension Principle and the optimization strategy, the fuzzy numbers are discretized with 4 α-cuts. The Figs. 7.11 and 7.12 present the fuzzy input parameters and the fuzzy buckling load used for the octagonal lattice structure whereas Figs. 7.13 and 7.14 are relative to the re-entrant cube structure. be selected. This last one requires a number of calculations clearly inferior to the EP method.

Indeed, a EP method relies in this case on 21952 calculations whereas a fifty of calculations are necessary for the OPT method.

For the two lattice structures, the experimental results are included in the bounds of fuzzy numbers. The intersection between the experimental values and the α-cut level give respectively for the two lattice structures 0.9 and 1. As these values are enough close to one, the adding of a slight observed variability on variable parameters is sufficient to enclose the experimental data. This section has highlighted that uncertain propagation techniques allows to obtain more confident results and more generally improve the predictivity of numerical simulations.

This chapter has allowed to investigate the uncertain linear buckling of lattice structures which was one of the objectives defined in the introduction of this thesis. Indeed, material and geometric uncertainty, experimentally identified in previous chapter and associated to the octagonal and reentrant cube lattice structures, have been modelled using fuzzy sets and probabilistic theories. An examination of all uncertain results clearly shows that the choice of uncertainty modelling has a strong influence on the analysis that can be made concerning the design of the structures.

Although all uncertain results enclosed the experimental ones, engineers take generally a decision thanks to a specific value such as the mean for the probabilistic theory or a defuzzified value for the fuzzy sets theory. For example, one can consider that mean values of the probabilistic result are more representative of the experimental ones for the octagonal lattice than for the re-entrant cube lattice (see figures 7.9 and 7.10). On the contrary fuzzy results are better representative of the experimental ones for the re-entrant cube lattice since the interval with the maximum degree of confidence perfectly enclose all experimental values (see figures 7.12 and 7.14). We can even add that the median value of this interval perfectly match the experimental value. Moreover, focussing on fuzzy buckling load of octagonal lattice (figure 7.12), the intersection between experimental results and the fuzzy number give a membership value of 0.9 which corresponds to a high level of confidence which relativizes the previous analysis. From these considerations it appears that the fuzzy uncertainty representation is more appropriate in this early design phase. Although this situation is not a general case this verdict is not completely surprising. Of course, in a more advance project phase, an ad hoc quantity of information to describe probabilistic uncertainty should certainly give better results than simple uniform distributions on relatively large intervals.

Anyway in all cases, it appears that uncertain modelling gives more insight that makes the designer more confident on the quality of its models. This fact would certainly be reinforced if experimental buckling tests have been performed on a more representative family of manufactured specimens.

General conclusions and perspectives

General conclusions

This thesis has allowed to investigate several numerical methods compatible with the linear buckling analysis of lattice structures described with uncertain material and non-topologic as well as topologic geometric properties. The observed variability is associated to the additive manufacturing, largely used nowadays to design highly complex part shapes. The uncertain The first chapter has summarized the main steps relative to the numerical solving of a deterministic linear buckling problem and presented the numerical application used in the first part of this thesis. The uncertain parameters are associated to material and geometric properties, such as Youngs modulus, angle or section of the structure. Finally, the response surfaces highlight the effect of section modifications on the buckling load values and the nonlinear evolution due to the topological modification of node coordinates.

The second chapter has allowed to compare different metamodels in the case of the approximation of linear buckling solutions and extract some advantages and drawbacks about each metamodel. We have highlighted a good precision between 0.5% and 1% of errors for RBF, Kriging, POD techniques to approximate the buckling loads. However, RBF method is not adapted to the approximation of the buckling mode whereas the quadratic regression is not efficient for the eigensolutions. Finally, the ROM-POD-Kriging has been proposed to couple the advantages of POD in a projection technique and Kriging. For this proposal, the level of precision (error inferior to 1% for the buckling load and MAC superior to 0.99 for the buckling mode) and the computational gain (close to 72 for the proposed multi-parametric analysis) is very interesting for large finite element models.

Several methods based on perturbations and projection techniques have been investigated

in the third chapter considering the same previous benchmark. First, we have shown that the modal stability assumption is not verified for the buckling problem. Second, Series or Pade approximants have a low convergence radius or require very high order of developments that implies important computational times. Third, the projection techniques, based either on HPP modes or on residues, are very precise for approximating both the critical buckling loads and associated buckling modes. However, the associated computational time is close to the reference one and so the computational gain is not significant for multiparametric analyses. Finally, a new method, called ROM-HPP-Kriging, has been developed and supplies an excellent compromise in terms of precision and computational gain. Indeed, the maximal error on buckling load is inferior to 1% and a MAC of 0.99 with a computational gain close to 74.

Several approximation methods, based on derivatives, are studied in the fourth chapter. These last ones rely on the availability of a suitable parameterization of finite element matrices as a function of input parameters. Except for a direct approximation with Taylor expansion series, the methods based on projection and eigensolutions derivatives supply interesting results in term of precision (error inferior to 1% for the buckling load and MAC superior to 0.99 for the buckling mode) for the studied benchmark. Considering the same strategy as in the two previous chapters, the ROM-Taylor-Kriging has been tested giving precision level similar to the projection techniques with derivatives. At present, the computational gain is close to 39 and is less interesting than the previous methods. This reduced performance is due to the computational time needed for the approximation of finite element matrices.

The fifth chapter proposes an assessment of numerical results proposed in the three previous chapters. ROM-POD-Kriging, ROM-HPP-Kriging and ROM-Taylor-Kriging performances are compared by considering different size of finite element models. Finally, we have highlighted that the ROM-HPP-Kriging gives the best compromise between precision and computational gain and is chosen for the second part of this thesis. This method will be introduced in uncertainty propagation methods. A flowchart, resuming the main steps of the ROM-HPP-Kriging method, is then proposed.

The sixth chapter is firstly dedicated to the presentation of the additive manufacturing and two studied lattice structures. Secondly, the variability associated to the EBM process is quantified thanks to a set of experiments. The variability is identified on Youngś modulus, diameters and angle. Moreover, quasi static tests are performed to identify the experimental buckling loads of each lattice structure. Next, we have proposed a finite element model of lattice structures and some deterministic buckling simulations have been compared to experimental results. The proposed analysis has confirmed the results developed in the Suardś thesis about the correction of the diameters and shown the necessity to integrate uncertainty in numerical simulations.

The seventh chapter discusses of the uncertainty propagation in numerical simulations. Two ways are here investigated considering two uncertainty modeling theories, namely fuzzy sets and probability. The ROM-HPP-Kriging method is integrated in Monte Carlo Simulations for the probabilistic formalism, and in Extension Principle and optimization technique for the fuzzy one.

The aim is to calculate the membership or the probability density functions of the buckling load and to compare them to the experimental results. Finally, the non-deterministic correlation highlights first the interest to consider uncertainty to improve the predictivity of simulations.

Second, we have shown a great flexibility of use for the proposed ROM.

Perspectives

This thesis has allowed to investigate several numerical methods on a same benchmark to facilitate the comparisons in terms of precision and computational time. All these methods depend on several key parameters that can certainly be optimized for future applications. For example, different correlations functions, such as Matérn one, could be tested for RBF or Kriging methods.

The scaling parameter for RBF could perhaps be updated to reduce the error on buckling mode.

For ROM-POD-Kriging method, as a decrease of the level of precision in comparison to ROM-POD is observed, it will be interesting to study the parametrization of Kriging in a way to improve the interpolation of the participation factors. Next, for ROM-POD-Kriging method, it will be interesting to study the convergence radius of this proposal as a function of the kind of input parameters and the sampling offline step. Indeed, the generation of samples and the associated size is another key parameter for this kind of study. Finally, the definition of samples with a Design Of Experiments or an iterative strategy based on observed errors can improve the quality of output results. Considering an improvement of these key parameters, another coupling could be investigated for example ROM-HPP-RBF.

We have highlighted the efficiency of the proposed reduced order model for studying uncertain buckling problems with fuzzy or probabilistic descriptions. Now, it is possible to investigate different designs of lattice structures with an uncertain description which generalizes the Suardś model. In this context, an experimental Design Of Experiments (DOE) can be proposed for studying more precisely the different uncertainty linked to the additive manufacturing and by considering different conditions of manufacturing. A focus can be performed about the quantification and the modelisation of the uncertainty.

As the solving of the linear buckling problem relies on a linear eigenvalue problem, it is first possible to apply the ROM-HPP-Kriging method to vibration problem and study the frequency and associated eigenmodes of mechanical structures. Second, we can extend the proposed method to the quadratic eigenvalue problem and define a reduced order model to investigate, for example, the stability problem of rubbing systems. Another way is to extend the current works to the case of nonlinear buckling problem by taking into account the advantages and drawback of the studied methods. The starting point is to analyze the behavior of Asymptotic Numerical Method, successfully used for buckling problems, and study the possibility to introduce a reduced model.

Finally, about the uncertainty propagation in a probabilistic context, it will be interesting to futher study the non-intrusive Polynomial Chaos Expansion method again and the possibility to integrate the proposed reduced order model or another technique, instead of using a classical regression metamodel, to improve the precision of results.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Linear buckling analysis . . . . . . . . . . . . . . . . . . . . . . . . . .

2 . 3

 23 Figure 1.1 -Buckling damage of : a) an engine connecting rod b) a vinification vat due to wine fermentation.

Figure 1 . 2 -

 12 Figure 1.2 -Buckling types : a) bifurcation of a column b) limit points or snap through.

Figure 1 . 3 -

 13 Figure 1.3 -Definition of Equilibrium states.

Figure 1 . 4 -

 14 Figure 1.4 -Finite element models developed during the thesis: a) rod b) beam c) brick.

Figure 1 . 6 -

 16 Figure 1.6 -Two-bars under compression: problem definition.

Figure 1

 1 Figure 1.7 -Two-bars under compression: 3D finite element mesh.

Figure 1

 1 Figure 1.8 -Two-bars under compression: static compressive deflection.

  Figure 1.10 -Two bars under compression: a) first mode λ cr = 11.09, b) second mode λ cr = 24.64, c) third mode λ cr = 39.59, d) fourth mode λ cr = 42.61.

  Figure 1.11 -Critical load response surface : (a) function of the Young's modulus and the angle, (b) function of the cross section (width and height).

Figure 2 .

 2 Figure 2.1 -Evolution of maximal errors on buckling load for quadratic model as a function of sample data size (a) two input parameters E and α, (b) four input parameters E, α, b and h.

Figure 2 . 2 -

 22 Figure 2.2 -Response surface and associated errors on buckling load for quadratic model with 4 parameters and 350 samples (a,b) parameters E and α, (c,d) parameters b and h.

Figure 2 . 3 -

 23 Figure 2.3 -Evolution of maximal errors on buckling load for RBF as a function of sample data size (a) two input parameters E and α, (b) four input parameters E, α, b and h.

Figure 2 . 4 -

 24 Figure 2.4 -Response surface and associated errors of buckling load for RBF model with 4 parameters and 50 samples (a,b) parameters E and α, (c,d) parameters b and h.

Figure 2 . 5 -

 25 Figure 2.5 -Evolution of maximal errors on buckling load for Kriging as a function of sample data size (a) two input parameters E and α, (b) four input parameters E, α, b and h.

Figure 2 . 6 -

 26 Figure 2.6 -Response surface and associated errors on buckling load for Kriging with 4 parameters and 350 samples (a,b) parameters E and α, (c,d) parameters b and h.

Figure 2

 2 Figure 2.7 -Evolution of maximal errors on buckling load for POD-Kriging as a function of sample data size (a) two input parameters E and α, (b) four input parameters E, α, b and h.

Figure 2 . 8 -

 28 Figure 2.8 -Response surface and associated errors on buckling load for POD-Kriging with 4 parameters and 350 samples (a,b) parameters E and α, (c,d) parameters b and h.

Figure 2

 2 Figure 2.9 -Evolution of maximal errors on buckling load for ROM-POD as a function of sample data size (a) two input parameters E and α, (b) four input parameters E, α, b and h.

Figure 2 .

 2 Figure 2.10 -Response surface and associated errors on buckling load for ROM-POD with 4 parameters and 15 samples (a,b) parameters E and α, (c,d) parameters b and h.

Figure 2 .

 2 Figure 2.11 -Evolution of maximal errors on buckling load for ROM-POD-Kriging as a function of sample data size (a) two input parameters E and α, (b) four input parameters E, α, b and h.

Figure 2 .

 2 Figure 2.12 -Response surface and associated errors on buckling load for ROM-POD-Kriging with 4 parameters and 35 samples (a,b) parameters E and α, (c,d) parameters b and h.

Figure 2 . 2 Case 2 : 4 :

 2224 Figure 2.13 -(a) Evolution of maximal errors on buckling load and (b) of the computational time as a function of number of degrees of freedom for 4 variables parameters.

  .32) The projection bases, whose size is [n dof × (n + 1)] with n << n dof , are orthonormalized with the iterative Gram Schmidt algorithm. This projection technique allows to decrease the size of the initial problem to the number of retained projection vectors , thus reducing the CPU time required to solve the modified problems, namely [(n + 1) × (n + 1)]. Combination of the homotopy development in perturbation technique or Homotopy Perturbation and Projection, is named HPP

  the associated high-order perturbed vectors.

Figure 3 . 1 -

 31 Figure 3.1 -Response surface and associated errors on buckling load for Modal Stability (a,b) parameters E and α, (c,d) parameters b and h.

Figure 3 . 2 -

 32 Figure 3.2 -Evolution of the maximal errors as a function of the order of approximation.

Figure 3 . 3 -

 33 Figure 3.3 -Response surface and associated errors on buckling load for ROM-HPP method with respectively a 23 th and a 3 rd order for static and eigenvalue problems (a,b) parameters E and α, (c,d) parameters b and h.

Figure 3 . 4 -

 34 Figure 3.4 -Evolution of the maximal errors as a function of the order of approximation.

  Maximal errors on buckling loads [%]

Figure 3 . 5 -

 35 Figure 3.5 -Evolution of the maximal errors as a function of the sample data size.

Figure 3 . 6 -

 36 Figure 3.6 -Response surface and associated errors on buckling load for ROM-HPP-Kriging method with 25 samples (a,b) parameters E and α, (c,d) parameters b and h.

Figure 3

 3 Figure 3.7 -(a) Evolution of maximal errors on buckling load (b) and of the computational time as a function of number of degrees of freedom.

  k = 1 • • • n p . The size of this projection matrix is [n dof × (n d + 1)] with n d the number of derivatives whereas the reduced problem is [(n d + 1) × (n d + 1)].

Figure 4 . 1 -

 41 Figure 4.1 -Evolution of the maximal errors on buckling load as a function of the sample data size and the order of truncation.

  Maximal errors on buckling loads[%]

Figure 4 . 2 -

 42 Figure 4.2 -Evolution of the maximal errors for Taylor expansion series as a function of the order of truncation.

Figure 4 . 3 -

 43 Figure 4.3 -Response surface and associated errors on buckling load for Taylor expansion series (a,b) parameters E and α, (c,d) parameters b and h.

Figure 4 . 4 -

 44 Figure 4.4 -Evolution of the maximal errors for ROM-Taylor as a function of the order of approximation.

Figure 4 . 5 -

 45 Figure 4.5 -Response surface and associated errors on buckling load for ROM-Taylor (a,b) parameters E and α, (c,d) parameters b and h.

Figure 4 . 6 -

 46 Figure 4.6 -Response surface and associated errors on buckling load for ROM-Taylor-HPP (a) parameters E and α, (b) parameters b and h.

Figure 4 .

 4 Figure 4.7 -Evolution of the maximal errors on buckling load for ROM-Taylor-Kriging and ROM-Taylor-HPP-Kriging as a function of the order of approximation.

Figure 4 . 8 -

 48 Figure 4.8 -Evolution of the maximal errors on buckling load for ROM-Taylor-Kriging as a function of the sample data size.

Figure 4 .

 4 Figure 4.9 -Response surface and associated errors on buckling load for ROM-Taylor-Kriging (a,b) parameters E and α, (c,d) parameters b and h.

Figure 4 .

 4 Figure 4.10 -(a) Evolution of maximal errors on buckling load (b) and of the computational time as a function of number of degrees of freedom.

  The three previous chapters have allowed to propose a new strategy to build a reduced order model relying on the interpolation of participation factors with Kriging. These participation factors are respectively associated to POD modes in Chapter 2, HPP modes in Chapter 3 and eigenvectors derivatives in Chapter 4. Three methods, named ROM-POD-Kriging, ROM-HPP-Kriging and ROM-Taylor-Kriging, have been proposed. The obtained results in the case of four parameters are summarized Fig.5.1 in term of precision and computational time.

Figure 5

 5 Figure 5.1 -(a) Evolution of maximal errors on buckling load (b) and of the computational time as a function of number of degrees of freedom.

Figure 5 . 2 -

 52 Figure 5.2 -Main steps of ROM-HPP-Kriging method.

  25.9% to $5.16 billion in 2015. The annual growth rate for the previous three years (2012-2014) was 33.8% and by 2021, the market is expected to grow to more than $26 billion. According to Frost & Sullivan Institute in their report of 2015 [www.frost.com], the aerospace, automotive, and medical industries are expected to account for 51% of the AM market by 2025. Figure 6.1 shows the future of Additive Manufacturing with a schematic of revenue generation in manufacturing sectors for the period of 2015-2025.

Figure 6 . 1 -

 61 Figure 6.1 -Additive Manufacturing market potential according to the Frost & Sullivan Institute [www.frost.com].

Figure 6 . 5 -

 65 Figure 6.5 -Controlled buckling of a 3D printed lattice structure : a) Top view of the sample under compression where controlled buckling zone is highlighted with red dots. b) Side view along the compression axis at 20% of compression [Paulose et al., 2015].

Figure 6 . 6 -Figure 6 . 7 -

 6667 Figure 6.6 -Stress-strain response of lattice honeycombs under compression with different relative densities, respectively. [Chen et al., 2018]. The lattice undergoes buckling with a clearly identified first critical load before collapsing.

Figure 6 .

 6 Figure 6.10 -Part build orientation in AM.

Figure 6 .

 6 Figure 6.11 -Manufacturing steps of a cylindrical lattice structure obtained by EBM process.

Figure 6 .

 6 Figure 6.13 -A cubic chiral lattice structure (a) with a high amplitude, (b) with a low amplitude [Warmuth et al., 2017].

Figure 6 .

 6 Figure 6.15 -Lattice structures based on different types of unit cells: (a) Cubic; (b) Diamond; (c) Truncated cube; (d) Truncated cuboctahedron; (e) Rhombic dodecahedron; (f) Rhombicuboctahedron [Ahmadi et al., 2015].

Figure 6 .

 6 Figure 6.16 -Aircraft seat absorber.

Figure 6 .

 6 Figure 6.17 -Description of the studied octagonal lattice structure : (a) EBM real part (b) CAD geometry.

Figure 6 .

 6 Figure 6.18 -2D layout of the CAD geometry of the studied octagonal lattice structure.

Figure 6 .

 6 Figure 6.19 -Description of the re-entrant cube lattice structure : (a) EBM real part (b) CAD geometry.

Figure 6 .

 6 Figure 6.20 -2D layout of the CAD geometry of the re-entrant cube lattice structure.

Figure 6 .

 6 Figure 6.21 -Strut size of the parts made by EBM [Yang et al., 2013].

Figure 6 .

 6 Figure 6.22 describes the concept of the geometrical diameter and the effectively working diameter for specimens manufactured by EBM. The respect of the desired geometry described by the CAD depends significantly on the size of the powders used during the building process. The dimensions of the powders are given by a Gaussian distribution centered on the value of 80µm

Figure 6 .

 6 Figure 6.22 -Illustration of EBM built-part contour roughness : (a) During the EBM process (b) After cooling and solidification.

Figure 6 .

 6 Figure 6.23 -Inscribed and circumscribed surface : (a) Tomographic reconstruction of a 1mm strut (b): Inscribed cylinder (red) inside the real strut (green)[START_REF] Suard | Towards stiffness prediction of cellular structures made by electron beam melting (ebm)[END_REF].

Figure 6 .

 6 Figure 6.25 -Uncertainty of the experimental setup.

Figure 6 .

 6 Figure 6.26 -SINTECH 20/D machine used for buckling experiments.

Figure 6 .

 6 Figure 6.27 -Experimental buckling of the octagonal lattice structure.

Figure 6 .

 6 Figure 6.29 -Mechanical response of the octagonal lattice structure.

Figure 6 .

 6 Figure 6.30 -Mechanical response of the re-entrant cube lattice structure.

Figure 6 .

 6 Figure 6.33 -Description of loading and boundary conditions of the octagonal lattice structure.

Figure 6 .

 6 Figure 6.35 -Comparison of the critical buckling load vs. total number of dof.

Figure 6 .

 6 Figure 6.37 -Description of loading and boundary conditions of the re-entrant cube lattice structure.

Figure 6 .Figure 6 . 39 -

 6639 Figure 6.38 -First buckling mode shape of the re-entrant cube lattice structure using ABAQUS © .

Fuzzy

  set introduces a degree of membership, represented by the membership function. This membership function describes the grade of membership for each element in the domain defined by the fuzzy set. The concept allows thus membership values between zero and one. Different forms of membership functions can be chosen, depending on the kind of imperfection considered. Triangular membership functions are the simplest and common choice to represent imprecisely known quantities. Indeed, one can fix a membership degree of 1 to the classical used value while membership values of 0 are affected to an estimate of the bounds of this value. However, with additional information, other membership functions can be chosen (π-shape, trapezoidal ....) to model the perception and the expert judgment of the user. Triangular or trapezoidal membership functions represented, Fig. 7.2, are a common choice when no specific information is available about imperfections.

Figure 7 . 2 -

 72 Figure 7.2 -Triangular and trapezoidal membership functions.

Figure 7 . 4 -

 74 Figure 7.4 -Propagation of uncertainty for fuzzy data.

Figure 7 . 6 -

 76 Figure 7.6 -(a) Evolution of maximal errors on buckling load (b) and of the computational time as a function of number of degrees of freedom for octagonal lattice structure.

Figure 7

 7 Figure 7.7 -Evolution of (a) leave-one-out and (b) verification error for the octagonal lattice structure.

Figure 7

 7 Figure 7.8 -Evolution of (a) leave-one-out and (b) verification error for the re-entrant cube lattice structure.

Figure 7

 7 Figure 7.9 -Histogram of buckling load for octagonal lattice structure.

Figure 7 Figure 7

 77 Figure 7.11 -Fuzzy input parameters of octagonal lattice structure.

Figure 7

 7 Figure 7.14 -Fuzzy buckling loads of re-entrant cube structure.

  Next, considering several uncertainty propagation methods and ROM-HPP-Kriging proposal, we have respectively studied the membership functions and the probability densities of the buckling loads. The results with proposed ROM have been confronted to those obtained with a classical Monte Carlo Simulation or Polynomial Chaos Expansion. Finally, the non-deterministic results have been successfully compared to experimental buckling loads. The different analysis and results have highlighted the good behavior of the PCE and the great flexibility of use for the proposed ROM in a multiparametric analysis context.

  numerical behaviors are finally compared to those obtained experimentally during quasi-static tests.The proposed work has been organized in two parts. The first one has allowed to test different metamodels or ROM methods, proposed in the literature, considering a same numerical framework and develop a new numerical method coupling homotopy perturbation technique, projection and kriging. This last proposal is then integrated in fuzzy and probabilistic uncertain propagation methods in the second part to calculate the membership function and the probability density function of the critical buckling loads. Two different lattice structures, for which experimental data are available, are considered in numerical applications. These structures are designed for energy absorption application and consider either octagonal or re-entrant cube unit cells.

  Table 1.2.

	Parameters	Lower bound Upper bound
	Young's Modulus E [GPa]	189	231
	Angle α [ • ]	10.5	19.5
	Width of section b [m]	0.09	0.11
	Height of section h [m]	0.09	0.11
	Table 1.2 -Inferior and superior bound of parameters.

  The i th sample is defined from a vector of n p design variables, noted p i = [p 1 ... p np ]. The matrix, noted P, aggregates all the vectors p

i such as the size is equal to [n s × n p ]. If the design variables are normalized, the vector and matrix of parameters are respectively noted pi = [p 1 ... pnp ] and P. Moreover, for each sample,

  the norm error err (zapp,z ref ) and Modal Assurance Criterion (MAC) M AC (zapp,z ref ) . The associated formula are recalled in Eqs. 2.24 and 2.25.

Table 3 .

 3 1 -Maximal errors on buckling load as a function of order of truncation for series.

	Maximal errors		Static step
	on buckling load [%]	Reference Order 1 Order 3 Order 5 Order 7
		Reference	-	57	32	18	0.81
		Order 1	1644	1415	1085	1653	1445
	Eigenvalue problem	Order 3	232	215	53	456	167
		Order 5	4	2904	641	29	3.50
		Order 7	0.07	1954	18	21	0.68

Table 3 .

 3

2 -Maximal errors on buckling load as a function of order of truncation for Padé.

Table 3 .

 3 3 -Maximal errors on buckling load as a function of order of truncation for perturbation projection.

	Maximal errors		Static step
	on buckling load [%]	Reference Order 1 Order 3 Order 5 Order 7
		Reference	-	9.4	3.44	0.21	0.08
		Order 1	10.70	37.60	14.60	10.90	10.70
	Eigenvalue problem	Order 3	0.03	9.20	3.48	0.24	0.05
		Order 5	0.004	9.23	3.45	0.22	0.08
		Order 7	0.004	9.20	3.46	0.22	0.08

Table 3 .

 3 4 -Maximal errors on buckling load as a function of number of residue for Residue iteration projection.

	Maximal errors		Static step
	on buckling load [%]	Reference Order 1 Order 3 Order 5 Order 7
		Reference	-	9.20	3.20	0.19	0.08
		Order 1	9.50	10.20	13.50	9	8
	Eigenvalue problem	Order 3	0.03	8.10	3.50	0.22	0.05
		Order 5	0.004	7.90	3.20	0.19	0.04
		Order 7	0.003	7.80	3.10	0.18	0.04

Table 5 .

 5 1 -Decomposition of CPU times for Kriging and ROMs with 48598 dof.

Table 6 . 3 -

 63 Geometric and material characteristics of the first studied octagonal lattice structure.

  .4.

	Parameter	Value Unit
	Length a	7.0 mm
	Length b	1.56 mm
	Angle β	90 deg
	Diameter d CAD	1.12 mm
	Young's modulus E	113.8-120 GPa
	Poisson ratio ν	0.342
	Density	4.43 g/cm 3

Table 6 . 4 -

 64 Geometric and material characteristics of the re-entrant cube lattice structure.

Table 6 . 5 -

 65 Material Elastic modulus: Influence of manufacturing direction.

	Type	Orientation	Elastic Modulus [GPa]
		0 •	109.2 ±6.3
	Machined material	45 •	121.5 ±9.2
		90 •	107.1 ±11.6
	Average		112.6 ±9.03

Table 6 . 6 -

 66 Significant set of values for the 5 selected parameters.

Table 6 .

 6 7 -Significant set of values for the 5 selected parameters.

Table 6 .

 6 8 -Comparison of critical buckling loads for octagonal lattice structure: experiments vs. numerical.

	Experiments		Test 1 [N]	Test 2 [N]
			26971	26751
	Numerical	CAD [N] Mean for ROM [N] [Suard, 2015] [N]
		44878	29867	22498

Table 6 .

 6 9 -Comparison of critical buckling loads for re-entrant cube lattice structure: experiments vs.For the two studied lattice structures, the buckling loads calculated with CAD description are clearly higher than those observed during experiments as explained by Suard. The correction of the CAD diameter is mandatory for our lattice applications. For the octagonal lattice structure, the mean and Suard's, descriptions, namely 29867N and 22498N , bound those obtained experimentally, corresponding to 26834N in mean. On the contrary for the re-entrant cube lattice structure where the diameter is closer, Suard's description, namely 19591N , is enough close to the experimental ones, either 19635N in mean.

	numerical.

  Towards reliable and robust designs . . . . . . . . . . . . . . . . . . . . 1.2 Probabilistic approach . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	1 Uncertainty propagation in mechanical engineering
	1.1 Towards reliable and robust designs
	Contents
	1	Uncertainty propagation in mechanical engineering . . . . . . . . . .
		1.1
		1.3	Interval and fuzzy sets approaches . . . . . . . . . . . . . . . . . . . . .

  Table7.5 -Mean results of the PCE studies for octagonal lattice structure.As a final evaluation for the probabilistic methods, Table7.7 gives, for the two structures,

	Order Truncation [q] N o of PCE	Size of	ε LOO	Err Y
			coefficients experimental
			N min	design
	3	1	56	75	5.12 10 -5 0.76
	3	0.75	26	35	7.57 10 -4 2.50
	2	1	21	29	9.10 10 -4 2.76

  Table 7.6 -Mean results of the PCE studies for re-entrant cube lattice structure.

	Order Truncation [q] N o of PCE	Size of	ε LOO	Err Y
			coefficients experimental
			N min	design
	3	1	56	70	2.60 10 -6 0.24
	3	0.75	26	33	5.95 10 -4 2.65
	2	1	21	27	9.10 10 -4 3.50

Remerciements

As reported by recent publications [START_REF] Javaid | Additive manufacturing applications in medical cases: A literature based review[END_REF], [START_REF] Brandicourt | Reconstruction of cranial vault defect with polyetheretherketone implants[END_REF], [START_REF] Berretta | Additive manufacture of peek cranial implants: Manufacturing considerations versus accuracy and mechanical performance[END_REF] and [START_REF] Bose | Surface modification of biomaterials and biomedical devices using additive manufacturing[END_REF] with its capability of producing complex shapes and lightweight parts, metal AM is becoming highly suited to the production of medical implants, orthopaedics and dental devices in titanium, providing innovative healthcare solutions (Figure 6.2). [ [START_REF] Bin | Customized ca-p/phbv nanocomposite scaffolds for bone tissue engineering: design, fabrication, surface modification and sustained release of growth factor[END_REF], [START_REF] Berretta | Additive manufacture of peek cranial implants: Manufacturing considerations versus accuracy and mechanical performance[END_REF].

More recently, researchers have found ways to print composite materials, such as carbon fiber and fiberglass, which improve the durability and the strength of 3D-printed products. For instance, [START_REF] Compton | 3d-printing of lightweight cellular composites[END_REF] Among lightweight structures obtained by AM, cellular or 3D lattice structures are probably

Influence of the part orientation in EBM process

In AM, the part building is done by stacking layers of material on top of each other producing staircase effect which leads to poor surface finish of the part. Therefore part build orientation appears to be a crucial parameter which affects the part quality. [START_REF] Das | Optimum part build orientation in additive manufacturing for minimizing part errors and support structures[END_REF] proposed a procedure to find optimal part build orientation which minimize the volume of support structures while meeting the specified geometric tolerancing criteria of the part.

[ [START_REF] Leutenecker-Twelsiek | Considering part orientation in design for additive manufacturing[END_REF] proposed a procedure for the early determination of the part orientation as depicted in Figure 6.9. Their method is based on dividing the part design into several elements which are analyzed separately. During this analysis the orientation's effects on the quality features of the elements are evaluated and rated. The final orientation of the part is determined based on this rating. Figure 6.9 -Determination of part orientation in AM to meet geometric tolerancing criteria [START_REF] Leutenecker-Twelsiek | Considering part orientation in design for additive manufacturing[END_REF].

Besides to the geometric tolerancing requirements for parts obtained by AM, the part build directions (Figure 6.10) during the forming process can have a direct influence of the mechanical characteristics of the part itself. For instance [START_REF] Kok | Anisotropy and heterogeneity of microstructure and mechanical properties in metal additive manufacturing: A critical review[END_REF] reported that among the influencing factors for the anisotropy of parts obtained by AM are : lack-of-fusion defects; phase transformation; layer banding and microstructural coarsening. They concluded that adjusting parameters obtained from the experiments for Ti6Al4V material are given in Table 6.5. We observe a non linear behavior as a function of the orientation during the manufacturing.

The highest values are obtained for an orientation of 45 • whereas the highest variation about elastic modulus is detected for an orientation of 90 • . The order of observed variation is classical to this kind of experimental test bench.

pression is shown in Figure 6.28. Here, we observed experimentally that the deformed lattice structure stays stable and collapse vertically when reaching the buckling load. This phenomenon was expected because the re-entrant cube lattice structure belongs to the auxetic family of lattice structures which are known to be stable because of their negative Poisson ratio. The mechanical response in terms of load-displacement behavior is shown in Figure 6.29.

The two blue and red curves are the results of two quasi-static tests obtained from two different specimens. The critical buckling load values are defined by the first limit point on the loaddisplacement curve. These values correspond to 26.917kN and 26.751kN respectively.

On the same manner, the mechanical response of the re-entrant cube lattice structure is shown in Figure 6.30. The critical buckling load values are 19.557kN , 20.066kN and 19.281kN corresponding to the first limit points on the curves.

5 Deterministic simulations and correlation

3D Beam element for buckling analysis

In order to perform the numerical simulation of the buckling analysis of the two studied lattice structures, the 3-dimensional Timoshenko beam finite element developed during this thesis.

Consider a 3D beam finite element given in Figure 6.31, we assume a two-node straight line beam of length l having six degrees of freedom (dof) per node i as u i , v i , w i , θ xi , θ yi , θ zi and hence 12 dof in total.

tions associated with each combination of generated parameters [Schueller, 1997, Schueller, 2001].

Simple to implement, this method presents guarantee in terms of convergence but can be proved prohibitive when the number of uncertain parameters increase or the deterministic method itself is time consuming. Many improvements have been proposed to reduce this drawback, such as the parallelization [START_REF] Johnson | Parallel implementation of monte carlo simulation -comparative studies from stochastic structural dynamics[END_REF], Papadrakakis and Kotsopulos, 1999, Johnson et al., 2003], the sampling technique [Cambier et al., 2002, Helton andDavis, 2003] or the used of meta-model [START_REF] Gayton | CQ2RS: a new statistical approach to the response surface method for reliability analysis[END_REF], Lew et al., 2006] which will be discussed in the next subsection. If a classification is considered as a function of uncertainty nature, two main categories can be built namely parametric and non-parametric methods. The parametric methods [Sudret, 2007] are the most popular and allow to take into account uncertainties on the parameters of the model. Amongst numerous proposals, the first methods were based on Taylor series expansion or perturbation techniques [START_REF] Benfratello | A perturbation approach for the response of dynamically modified structural systems[END_REF], Chakraborty and Dey, 1998, Nieuwenhof and Coyette, 2002, Falsone and Ferro, 2005] with the aim to substitute the equations of numerical model and the random functions by the mathematical developments to quantify the first statistical moments (mean, variance) of output solutions [Sudret, 2007]. These methods are dependent on the validity of Taylor series expansion and can only be applied for linear static and dynamic problems [START_REF] Handa | Application of finite element methods in the statistical analysis of structures[END_REF], Liu et al., 1986, Shinozuka and Yamazaki, 1988] with low output dispersions around the nominal values as shown by Ghanem and Spanos [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF]]. These last years, authors were particularly interested in methods relying on polynomial chaos expansion developed by Wiener [Wiener, 1938a]. The first works in mechanical engineering are due to Ghanem and Spanos [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF] who have combined Wiener-Hermite developments with finite element method to model and propagate uncertainty in mechanical structures. A more general extension, called generalized polynomial chaos, was proposed by Xiu et al. [Xiu and Karniadakis, 2002a]. The method is based on the correspondence between the probability density functions of certain random variables and the weight functions of orthogonal polynomials of the Askey scheme. The polynomial chaos method allows the stochastic contribution of a random variables and the deterministic contribution to be dissociated.

The random part is then decomposed on a basis of orthogonal polynomials (for example Hermite polynomials for Gaussian variable, Legendre for uniform variable). Intrusive and nonintrusive methods have been successively applied to calculate stochastic modes useful to uncertainty propagation. The intrusive method is based on the Galerkin projection to expand the uncertain model on the polynomial chaos basis as a linear system of deterministic coupled equations [Ghanem andSpanos, 1991b, Babuska et al., 2004]. Nevertheless, this way is limited to a low number of random parameters and mainly applied to the case of linear number of samples of input random parameters then compute for each sample the response quantity under consideration through the deterministic model. The last step of the method consists in performing a statistical treatment of the output samples.

Suppose the response of a system is modelled by a mathematical function or an algorithm M : R p → R and x is the vector of input parameters. Uncertainty on the input parameters is modelled by a random vector X ∈ R p described by the joint probability function f X (x). The response Y is also a random variable and is given by:

Assume that a sample set of N input vectors has been generated, unbiased estimates of the mean value µ Y and standard deviation σ Y of a response quantity Y are given by:

Additionally, the coefficient of variation of a response quantity Y which gives a measure of its relative dispersion is given by:

In practice, the quantities of interest in Eq 7.2 and Eq 7.3 are updated continuously during the process after each run. The accuracy of the MCS can be estimated by the coefficient of variation of the empirical mean given by equation Eq 7.5. Generally, a value of 1% or 5% is chosen as a threshold for the convergence loop.

The precision of the method is mainly dependent on the number of simulations which can be very time consuming. Nevertheless, in many applications, MCS serves as a reference for the evaluation of random response quantities given by alternative methods.

Polynomial chaos expansion

Polynomial Chaos Expansion (PCE), firstly introduced by Norbert Wiener [Wiener, 1938b] is a powerful metamodeling technique to evaluate uncertainty of a dynamical system thanks to a probabilistic description of the system parameters. A first application of PCE in structural mechanics is due to Ghanem and Spanos [Ghanem and Spanos, 1991a] who proposed an efficient construction of the PCE for representing second-order stochastic processes and random fields.

They used PCE for solving boundary value problems with a stochastic finite elements method.

PCE was initially studied for Gaussian random variables. In a way to extend the method to any random distributions, Xiu and Kaniakakis [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF] have proposed a generalization of the original Wiener chaos expansion. Given a random vector with independent components X ∈ R p described by the joint probability density function f X (x) and a model M : R p → R, from the Wiener theory and the generalized Cameron-Martin theorem, any second order random process can be expanded in a polynomial function series, such as:

where Ψ α (X) are multivariate polynomials orthonormal with respect to f X , α is a multiindex that identifies the components of the multivariate polynomials and real values y α are the coefficients of the polynomial chaos. In practice, the sum in Eq. 7.6 is truncated which leads to the truncated PCE given by:

The polynomial basis Ψ α (X) is classically built from a set of univariate orthonormal polynomials of degree k, φ

k (x i ) associated to each input parameters (i) and is given by:

.1 -List of common random variable and their corresponding univariate polynomials.

The correspondence between families of univariate polynomials and common random variables to which they are orthonormal is given in Table 7.1. Considering a maximum order of truncation d, the standard truncation scheme corresponds to all polynomials in the p input variables of total degree less than or equal to d: defined for reference sets X 1 , X 2 , . . . , X n , the Extension Principle defines a fuzzy subset B of Y using data from fuzzy subsets A 1 , A 2 , . . . , A n of X. A fuzzy characterization of the membership function in Y is written as follows:

Zadeh's Extension Principle has already been implemented in its discrete form to propagate uncertainties in mechanical engineering applications and is called the Transformation Method [START_REF] Moens | Non-probabilistic finite element analysis for parametric uncertainty treatment in applied mechanics: Recent advances[END_REF]. In practice, this approach initially requires a discretization of input membership functions following the level of confidence in Fig 7 .3 rather than the support. This strategy defines a set of intervals x inf i ; x sup i α for which a level of confidence is also associated and thus transforms the fuzzy problem in Eq 7.19 into several interval problems in Eq 7.20.

[ỹ 1 , . . . , ỹm ] = f(x 1 , . . In this chapter, ZEP is considered to calculate the membership function of buckling loads.

Next, a calculation of deterministic solutions corresponding to all combinations of fuzzy input parameter values is performed to determine all possible evolutions for output data. Finally, an evaluation of the degree of confidence for all solutions is considered to build the membership functions of output data.

Since the implementation of the Extension Principle is time consuming to propagate the uncertainty with precision in a finite element context, several methods have already been proposed in the literature and are known as the fuzzy finite element method. These methods are based either on specific combinatorial samplings [START_REF] Moens | Non-probabilistic finite element analysis for parametric uncertainty treatment in applied mechanics: Recent advances[END_REF], on interval arithmetic [Massa et al., 2008a], on fuzzy logic [Massa et al., 2015a], on series development [START_REF] Massa | Fuzzy eigensolutions of mechanical structures[END_REF] or on optimization problems [Massa et al., 2011a]. All these methods have a common goal, which is to precisely determine all modifications of behaviours of studied solutions in a time frame compatible with a mechanical design step.

To reduce the conservatism of the output solutions, we propose to build the output intervals by considering an optimization technique rather than interval arithmetic which can seriously overestimate the output set because of the same parameter [Moens andHanss, 2011, Massa et al., 2008a].

The interval problem, defined for each α-cut level, is transformed as a min-max optimization problem in Eq 7.21 for each output data y i , defined as follows:

where f is the function representative of the deterministic problem to solve.

The choice of optimization technique is mainly guided by the functional dependence of the solutions, the kind of studied uncertain parameters and the mathematical problem to solve (linear system for static analysis, eigenvalue problem for modal analysis, time or frequency dependence for dynamic analysis). For certain parameters, simple rules can be defined because the functional dependence of the solutions is monotonic. For example, eigenvalues tend to increase as Young's modulus increases and decrease as density increases. However, for other parameters, the functional dependence is not always monotonic, particularly with large variations. Moreover, functional dependence can be different for each component of the solution vector.

In the present case, the functional dependence of the solutions are not always monotonic and we propose to consider a local gradient algorithm adapted to the management of the different α-cut levels in Fig 7 .4.

optimization scheme) with a reanalysis technique. The proposed strategy is described in the next section.

2 Buckling analysis with uncertainty for lattice structures

Organization of calculations

The two lattice structures, described in the previous chapter, are here studied by considering the five uncertain parameters. For each lattice structure, a ROM is built by using ROM-HPP-Kriging method and integrated in fuzzy sets and probabilistic uncertainty propagation methods.

All the investigated non-deterministic simulations are summarized in Fig 7 .5.

Figure 7.5 -Organization of non-deterministic calculations.

To evaluate the efficiency of proposed methods based on ROM, reference calculations are respectively performed with MCS and Extension Principle implementation. Moreover, a polynomial chaos is calculated for the probabilistic way whereas an uncertainty propagation by optimization is proposed for the fuzzy sets theory.

Uncertain parameters description

Description of uncertain parameters is done using the same rules for both lattice structures.

The main idea is to maximize the use of information in hand for the description of fuzzy numbers. In a way to compare the results and considering the weak quantity of experimental data, probabilistic uncertainty is built considering uniform laws where bounds are those of the fuzzy number supports. Considering that samples are different, comparison of the histograms for the same structure confirms again the quality of the evaluation with the reduce models. Concerning the experimental values, in both cases they are enclosed in the variation domains given by uncertainty propagation simulations. However, one must recognize that mean values are much closer to the experimental ones for the octagonal lattice structure than for the re-entrant cube lattice structure. Nevertheless, this section has highlighted that uncertain propagation techniques allows to obtain more confident results and more generally improve the predictivity of numerical simulations.