Advanced computational technique based on kriging and Polynomial Chaos Expansion for structural stability of mechanical systems with uncertainties

In this paper, a numerical strategy based on the combination of the kriging approach and the Polynomial Chaos Expansion (PCE) is proposed for the prediction of buckling loads due to random geometric imperfections and fluctuations in material properties of a mechanical system. The original computational approach is applied on a beam simply supported at both ends by rigid supports and by one punctual spring whose location and stiffness vary. The beam is subjected to a deterministic axial compression load. The PCE-kriging meta-modelling approach is employed to efficiently perform a parametric analysis with random geometrical and material properties. The approach proved to be computationally efficient in terms of number of model evaluations and in terms of computational time to predict accurately the buckling loads of a beam system. It is demonstrated that the buckling loads are substantially impacted not only by both the location and the stiffness of the spring, but also by the random parameters.


Introduction
The problem of structural stability of mechanical systems [1], such as buckling of beam structures under axial compressive loading has attracted the attention of many researchers. Generally speaking, the first step in engineering for the prediction of buckling is to consider idealized structures with mathematically exact geometries. In reality, structures have imperfections due to the manufacturing process and cannot be considered as homogeneous or geometrically exact. So, the effect of material and geometric imperfections needs to be considered for a robust prediction of buckling and for the estimation of the main differences with the idealized buckling problem. Also, many researches carried out in the field of structural stability theory and practice for buckling of beams were focused E. Denimal Univ. Gustave Eiffel, Inria, COSYS/SII, I4S, Campus de Beaulieu, 35042 Rennes, France e-mail: enora.denimal@inria.fr to get deep insights in the role played by different parameters, which is not possible with the different approaches presented previously.
The objective of the present work is then to predict the buckling behaviour of a mechanical system when the two kinds of uncertainties are present, namely epistemic and aleatory. The considered system is a beam where all geometrical and material properties are considered as uncertain and can vary. This uncertainty is modelled through a PDF for each parameter. The buckling properties of the beam system are studied when a spring is added to the latter. The equation of motion of the system is obtained with the Rayleigh-Ritz method. The spring location and stiffness can vary, modifying strongly the buckling properties. Indeed, the buckling loads as well as the buckling mode shapes are impacted by this additional stiffness. Maximum buckling loads are observed at specific location of the spring and a saturation effect is observable when the spring stiffness increases. By considering the uncertainties related to the beam model, the stochastic properties of the buckling loads are analysed by using the PCE-kriging meta-model. It is demonstrated that with a low number of model evaluations and a drastic reduction of the computational time, the meta-model is able to predict accurately the different buckling loads and their stochastic moments. Finally, a sensitivity analysis based on the Sobol indices directly assessed from the PCE-kriging formulation is achieved to investigate in depth the influence of the different beam model parameters.

Hybrid surrogate model presentation
In the present work, random variables are employed to describe the aleatory uncertainty as it is one of the most traditional and used way to model it. The epistemic uncertainty is modelled through parametric variables that can take the values in a given interval. They will be denoted as parametric and random uncertainties or variables in the following for clarity considerations. Each of them is treated in a different manner to take into account their characteristics.
The general problem consists then in approximating the function f that depends on x and on ξ , where x is the vector of parametric variables of dimension n and ξ the vector of random variables of dimension m.

Kriging
For epistemic uncertainties, surface response approaches are well adapted due to their non-intrusive characteristics. The most well-known one being the polynomial regression, despite its limited efficiency especially when nonlinear behaviours or when functions with discontinuities are approximated. Another approach is kriging, which is the best unbiased predictor [16]. It has proven its efficiency in many applications, even for non-linearities and/or discontinuities [15], and will be used in the following.
Let h denote the function to be approximated from a limited number of evaluations that depends on x ∈ R n , then a kriging approximation of h is given by [16,29] where (g k ) is a basis of q orthogonal regression functions, often chosen as polynomials of low order (up to order 2), β the vector of the regression coefficients to be determined and Z a zero-mean Gaussian process of variance σ 2 . The covariance matrix is given by where R ∈ [0, 1] is the spatial correlation function and E the expectation operator. The correlation function is usually unknown but constructed from a family of kernel functions parameterized by the hyper-parameter θ , and depends on the distance between x and x . If n > 1, then the correlation function is constructed by taking the tensor product of a univariate kernel function family k: where θ j is the value of θ in the j-th dimension, and d j = x j − x j the distance between x and x in the j-th dimension. For isotropic kriging, θ j = θ is a constant, and in anisotropic kriging θ j is different in each direction.
In Table 1, examples of kernel functions in dimension 1 are given [30].
To construct a kriging meta-model, N evaluations of the model are required, i.e. a set of inputs x (1) , . . . , x (N ) and outputs y = y (1) , . . . , y (N ) , called learning set or experimental design (ED). From this set, the regression matrix G of coefficients G i j = g j (x (i) ) can be constructed as well as the matrix R of coefficients Then, three parameters must be determined, namely the regression vector β, the process variance σ 2 and the hyperparameter θ. The first step consists in the determination of θ by maximizing the likelihood function L [16,29]: where C = σ 2 R and the estimation of β is β = (G T RG) −1 G T R −1 y and the estimation of the process variance is The kriging predictor at a new point x 0 is then given by [16,29] h( where g 0 is the vector of g j (x 0 ) and r 0 is the vector of R(θ , x (i) , x 0 ).

Polynomial Chaos Expansion
For random variables, stochastic modelling is particularly adapted as it is based on the stochastic properties of the random variable. According to the Polynomial Chaos theory, the second-order random-parameter function Y (ξ ), which depends on ξ the vector of m independent random variables (ξ i ) i∈ [1,m] , can be approximated by a convergent expansion [19,20,31,32]: where Ψ J are the multivariate orthogonal polynomials and (a J ) J ∈[0,P−1] are the weighting coefficients to be determined. Each random variable (ξ i ) is described by its PDF f i and the joined PDF of ξ is f = m i=1 f i . The multivariate polynomials Ψ J are obtained by the tensorization of orthogonal mono-variate polynomials (P (i) j ) of order j and orthogonal to the PDF f i . It writes [19,20,31,32] The inner product between the two polynomials P i and P j that depends on a random variable X of PDF f X is given by If c o is the chaos order, i.e. the maximum order of polynomials kept in the expansion, then P = m+c o c o . The sets of polynomials (P (i) j ) are chosen based on the Askey scheme that gives a correspondence between classical PDF and the associated orthogonal polynomial family [20,31]. This is illustrated in Table 2 for a few families. The number of terms in the PCE can increase drastically when high chaos order is considered or when the number of random parameters increases. To limit the number of terms in the expansion and avoid convergence issues, some truncation strategies are possible [33], here the hyperbolic truncation norm is adopted to select the polynomials to keep [34]. It is based on the idea that models and physics are often driven by low-order interaction effects. The polynomials that are kept in Eq. (6) verify where γ ∈ (0, 1] is the chosen hyperbolic norm coefficient. The coefficients (a J ) J ∈[0,P−1] can be computed based on intrusive and non-intrusive approaches. The non-intrusive approach based on the regression method has proved to be efficient [32,34] and is used here. The (a J ) J ∈[0,P−1] are computed by minimizing in the least-square sense the distance between the random function Y (ξ ) and and its PCE evaluations at M points.

Hybrid approach coupling kriging and PCE
In this section, the hybrid method that couples kriging and PCE to deal with both random and parametric parameters is presented [28]. As previously explained, the objective is to create a surrogate function x is a vector of parametric variables of dimension n and ξ the vector of random parameters of dimension m. First, the random part is decomposed using the PCE approximation: where the chaos coefficients (a J ) J ∈[0,P−1] now depend on the parametric variables x. Each PCE coefficient is then approximated based on a kriging meta-model: Finally, it comes

Exploitation of the coefficients
The PCE decomposition gives a direct access to many statistical information without any additional simulations, that are directly dependent on the parametric variables here, which is the main interest of the approach. Indeed, if a unique kriging meta-model were build, then it would have been necessary to perform an MCS simulation. Thus, the average of f at the points x is given by and the variance is given by The Sobol indices can also be directly estimated from the PCE decomposition due to its uniqueness. If V ( f (x)) is the variance of the output at the point x, then the first-order Sobol index S i related to the variable ξ i is [35] where The latter is directly deduced from the PCE coefficients [32]: where α i is the set of the multivariate indices for which only polynomials related to the variable ξ i are present.

Equations of motion
The mechanical system under study is shown in Fig. 2. It is composed of a beam of length L beam and of rectangular cross section with inner sides b i and h i and external sides b o and h o . The system is simply supported at both ends by rigid supports and at position x = L spring by a spring of stiffness k spring . The beam is loaded by a uniform axial force P which acts on the beam neutral axis. The beam deflection is described by a deflection function v(x, t) (i.e. amount of vertical displacement at position x on the beam from the left end). In the following, the classical Rayleigh-Ritz method is proposed to formulate the mathematical problem and obtain the equations of motion of the studied system. As a reminder, the Rayleigh-Ritz approach that corresponds to an approximate energy-based method has been widely used to solve problems of solid mechanics and more particularly to investigate the dynamics (such as the natural frequencies and mode shapes) as well as the buckling behaviour of vibrating beam structures [36]. Using the Rayleigh-Ritz method, the displacement function v(x, t) of the beam is approximated as a linear combination of n trial functions where λ i (t) are arbitrary coefficients. Each trial function φ i should satisfy at least the essential boundary conditions (i.e. geometrical boundary conditions). In the present case, these two boundary conditions are given by zero displacement at the right-hand and left-hand supports: v(0) = 0 and v (L beam ) = 0. It is to be noted that the use of polynomials or trigonometric series is classically used in the modelling of many physical problems such as buckling and vibration. Finally, the original variational problem of extremizing the functional integral F [v (x, t)] is replaced by the problem of finding the coefficients λ i that extremizes F (λ i , . . . , λ n ). So minimizing the energy functional F can be done by taking partial derivatives with respect to λ i such as ∂ F ∂λ i = 0. The convergence of the Rayleigh-Ritz process implies that the approximation given in Eq. (17) will tend towards the exact deflection function v(x, t) that extremizes the integral F [v (x, t)] when n tends to infinity. For the beam system under study, the maximum kinetic energy T is given by where ρ and S define the beam density and the cross-sectional area (i.e. S = b e h e − b i h i ), respectively. The maximum strain energy V is composed of the contributions of the beam and the punctual spring k spring such as where E and I define the Young's modulus of the beam and the second moment of area of the beam cross section , respectively. Then the amount of work W P performed by the external concentrated load P from the initial position to the bent situation is equal to where v P corresponds to the horizontal displacement generated by the load P. By using a Taylor expansion and neglecting the higher order terms, the horizontal displacement v P can be approximated by [1] v P = 1 2 By equating the total potential energy Π = −T + V + W P and substituting v(x, t) by the approximate admissible displacement provided in Eq. (17), the equation of motion with respect to the generalized coordinates λ i is given by with for i, j = 1, . . . , n.
Considering the essential boundary conditions to be verified, a linear combination of sinusoidal functions are applied for approximating the displacement function v(x, t): Substituting Eq. (27) in Eqs. (24), (25) and (26), the mass matrix M and the stiffness matrices K and K P are given by Assuming a solution of the form Λ(t) = ue iωt (where ω corresponds to the angular frequency of the system and u the associated eigenvector in the generalized coordinates λ i ), system (22) becomes Eigenvalues are determined by solving the associated characteristic equation (i.e. det −ω 2 M + K − K P = 0). It is recalled that each angular frequency ω i resulting from the discretization of the displacement variational principle by the Rayleigh-Ritz method is an upper bound to the corresponding exact angular frequency It is often very difficult to determine the exact critical buckling load in structures using the classical Euler's theory. Therefore, the estimation of buckling loads is often approximated using energy conservation and referred to as an energy method in structural analysis. The principle of stationary potential energy can be used to find the critical loads by considering the total potential energy given by Π = V + W P (i.e. at the critical load, the neutral equilibrium is possible, therefore, Π P = 0). It corresponds to the simplification of the eigenproblem previously discussed (see Eq. (31)) to its restriction to the sudden change in shape (deformation) of the beam under load. For a linear buckling analysis, the eigenvalue problem to be solved to get the jth buckling load multiplier α j and the associated buckling mode shapes Ψ j is defined by K − α j K P Ψ j = 0. The jth buckling load is then provided by P b j = α j P. It is to be noted that buckling load multiplier is an indicator of the factor of safety against buckling: if α j > 1 the initial applied load is less that the estimated critical buckling load P b j and buckling are not predicted. The natural frequencies of the mechanical system with the axial applied load P are provided by Eq. (31). If α j < 1 the initial applied load exceeds the estimated critical buckling load P b j and buckling will occur. For α j = 1 the applied load is exactly the critical buckling load P b j and buckling is expected. The critical loads obtained will be an approximate value of the exact critical loads due to the fact that the exact deflected shape is not available. The approximated deflected shape assumed is provided by Eq. (17). Moreover, it is worth mentioning that in the present case, the presence of a spring with variable position and intensity makes impossible an analytical expression of the buckling loads and therefore requires numerical processing. Therefore, taking into account uncertainties (i.e. random geometric imperfections and fluctuations in material properties) for the prediction of buckling loads requires an efficient numerical methodology.

Preliminary deterministic study
First of all, evolutions of the first three natural frequencies of the beam system subjected to the load P are displayed in Fig. 3a for different values of both the stiffness k spring and position of the spring L spring . The values of the physical parameters used for the beam system are indicated in Table 3. It appears very clearly that each frequency has a potentially very different evolution from the others depending on the two parameters of the spring k spring , L spring . Density ρ (kg/m 3 ) 7800 One can more particularly observe in certain cases a crossing of two frequencies when the value of the applied load P increases (see crossing between the first and second frequencies (second case, in blue) and crossing between the second and third frequencies (third case, in red)). This well-known effect in the field of vibrations is due to the position of the spring L spring in regard to the nodes and anti-nodes of mode shapes: if the spring is close to a node its influence will be less important than if it is located close to an anti-node (location where the amplitude is maximum). Of course, for a given position L spring the increase of the stiffness k spring leads to an increase of the frequencies (except if the spring is located on a node of the associated mode shape). It is well known that gradually increasing the load P decreases the natural frequencies of the beam system. When the load reaches a critical level, the beam system may suddenly change its shape. This phenomenon called buckling corresponds to the fact that a natural frequency of the system has reached a zero value. For example, the first three buckling loads P b i (for i = 1, . . . , 3) are identified by a circle •, a square and a triangle Δ in Fig. 3a. The associated buckling mode shapes are also given in Fig. 3b. Once again these results clearly demonstrate the impact of the two spring parameters k spring , L spring on the value of the buckling loads and also on the fundamental and second buckling modes. In conclusion, the eigenfrequencies, the associated eigenmodes as well as the buckling loads of the mechanical system are dependent on the combination P, k spring , L spring .
In order to have a more global overview of the influence of the two parameters k spring , L spring on the first two buckling loads P b 1 and P b 2 , simulations are conducted for 100 × 100 combinations of these two deterministic inputs  (17)). A convergence study has also been conducted in order to verify the convergence of the results on the estimation of the first two buckling loads as a function of the number of functions selected in the Rayleigh-Ritz method. Figure 5 give the two values Δ i,n max and Δ i,n mean defined by where P b i,ref defines the ith reference buckling load (chosen for n = 50) and P b i,n corresponds to the ith approximate buckling load based on n trial functions (with n ≤ 30). These two indicators Δ i,n max and Δ i,n mean are estimated by considering the 100 × 100 calculations in regard to the two deterministic inputs k spring , L spring . As shown in Fig. 5, using 10 functions φ i (x) allows to have a mean error Δ i,n mean of the global estimate lower than 0.02% and a maximum error Δ i,n max lower than 0.1% for both the first and second buckling loads P b 1 and P b 2 . Coming back briefly to Fig. 4, the buckling loads evolutions are symmetric about L spring = 0.5 m. The first buckling load is maximal when the spring is located at the centre of the beam, whereas for the second buckling load the maximum is reached when the spring is at the third or two-third of the length of the beam. The buckling loads tend to increase when the spring stiffness increases as well until a threshold is reached. The latter is reached at about 2.5×10 8 N/m for the first buckling load, and at about 6.4 × 10 8 N/m for the second buckling load.
Finally, it is important to remember that many physical parameters of the beam system can have a significant influence on the buckling loads. Some parameters may vary throughout the life of the mechanical system or some uncertainties have to be taken into account during the manufacturing process. All of these facts can affect the estimation of buckling loads and their impacts must therefore be quantified. However, conducting such a study requires a considerable number of simulations, resulting in excessive calculation times. As a result, it is necessary to deploy new strategies in order to considerably reduce the number of calculations and to respond to this current problem. The main objective of the rest of study is to promote the use of PCE-kriging surrogate model and to  illustrate its efficiency to accurately predict the first two buckling loads P b 1 and P b 2 of the beam system with multiple uncertainties.

Application and numerical examples
In this section, the PCE-kriging approach is employed to predict the two first buckling loads P b 1 and P b 2 when the uncertainties present in the model are considered. In the present study, all the model parameters are considered as uncertain, and are described by a normal distribution. It represents a total of six uncertain parameters: the Young modulus E, the beam length L beam , the inner base b i , the inner height h i , the outer base b o and the outer height h o . Their characteristics (mean and standard deviation) are given in Table 4. The choice was done to ensure that 95% of the samples represent a 5% variation around the mean value for each parameter. It is worth noticing that the density is not considered as random, as it has no influence on P b 1 and P b 2 . To be noted that the effect of boundary conditions is not investigated in the proposed study. The system is thus subjected to two types of varying parameters: a parametric variation of the spring location L spring and stiffness k spring , and the random aspect of all the parameters that characterize the beam. To treat with these two natures of parameter, the PCE-kriging approach is employed.

Construction of the PCE-kriging meta-models
As a first step, the experimental design is constructed. As the system is symmetric about L spring , the location of L spring is taken between 0 and 0.5L beam . A set of 50 points generated by an LHS of dimension 6 is taken for the PCE construction, and a set of 51 points constructed by a LHS of dimension 2 that satisfies a maximin criterion is created for the kriging. Finally, the buckling loads P b 1 and P b 2 are computed for the 50 × 51 = 2550 configurations of the system.
For each value of L spring and k spring , a PCE is constructed for each buckling load. Based on a convergence study realized for different spring properties and not presented here for the sake of concision, the chaos order is taken equal to 5 and a 0.5-hyperbolic norm is chosen, so the PCE basis is composed of 40 terms. The PCE coefficients are computed based on the regression approach for each buckling load and for each value of L spring and k spring . At this stage, a first validation step has been undertaken. On a 11 × 11 grid of the spring properties (stiffness and location), at each grid point a PCE for each buckling load was created. In addition, for each point of the grid, 1000 random samples were also generated. The PCE predictions were compared to the reference values. In each case, the PCEs have shown good prediction quality, not shown here for the sake of concision.
The second step consists in the construction of the kriging meta-models for each PCE coefficient, i.e. for each of the 40 coefficients of the PCE basis for each load. For the first buckling load P b 1 , a spline correlation with a zero-order polynomial regression is taken, whereas a matérn 3/2 correlation with a zero-order regression kriging is taken for the second buckling load P b 2 . The final PCE-kriging meta-models are validated by comparing the relative error on the mean and the variance of each buckling load on a 26 × 51 grid between the meta-model and reference values. The reference mean and variance are determined from a LHS of 1000 points, which has reached convergence. The evolution of the relative error for the mean and the variance of the first buckling load P b 1 (resp. P b 2 ) are given in Fig. 6 (resp. Fig. 7) on the right column. The red points correspond to the 51 points of the kriging-experimental design, i.e. points where PCE are constructed. As a visual indicator of the quality of the meta-models, the average value and variance of each buckling loads of the PCE-kriging meta-models and of the reference are also displayed on the left column. The reference results obtained from Monte Carlo Simulations are given in purple and the PCE-kriging predictions are given in green.
Considering the first buckling load, the relative error is always inferior to 5% for the average prediction, as illustrated in Fig. 6b, d, demonstrating the high quality of the prediction. Looking at Fig. 6a one can see the very good agreement between the predictions and the reference results. In average, the buckling load P b 1 increases when L spring increases as well. When k spring increases, P b 1 increases as well but reaches a threshold. Maximal values of the buckling load P b 1 are then reached when the spring is fixed at the centre of the beam and when k spring is larger than 2.5 × 10 8 N/m. For the variance, except a peak of 50%, the relative error is always inferior to 10% which shows here again the high quality of the PCE-kriging meta-model. The error peak is reached when L spring /L beam = 0.5 and for low values of k spring . Comparing to Fig. 6c, it corresponds to a zone with a large and brutal variation of the variance. Higher quality of the meta-model could be reached by adding more points in this area. Despite this peak, one can see in Fig. 6c that the variance prediction remains substantially satisfactory. The variance tends to increase when both the spring location and stiffness increase. When the spring is near the base or when the stiffness is low, the variance is low and the buckling load P b 1 is not significantly influenced by the uncertainties. However, when the spring is getting closer to the centre of the beam or when the spring stiffness increases, the variance increases and the bucking load is significantly impacted by the uncertainty present in the model.
Considering the second buckling load and Fig. 7, the relative error on the average is always inferior to 3%, as illustrated in Fig. 7b, d, which demonstrates again the accurate prediction of the PCE-kriging meta-model. Looking at the comparison between the predictions of the average buckling load P b 2 and the reference in Fig. 7a, one can see the good agreement between the two surfaces. The maximum buckling load is reached when the spring is located at the one-third of the beam length (i.e. L spring /L beam = 0.33 or L spring /L beam = 0.66 due to the symmetry). As for the first buckling load, the second one tends to remain constant when the spring is near the base or when the stiffness is low. Similarly, a threshold of the maximum buckling load is reached when k spring increases. The latter is reached later compared to the first buckling load (2.5 × 10 8 N/m versus 6.4 × 10 8 N/m). Similar analysis can be done for the variance, the relative error is globally low (inferior to 10%) except at two peaks that correspond to brutal variation of the variance. The variance keeps low value when the spring is near the base or has a low stiffness, and reaches high values when the buckling load is maximal.
It is worth noticing here that the standard deviation corresponds to the square root of the variance and that in the case of a Gaussian distribution, 95% of the data are in the interval i with i = 1, 2 and σ i defined the associated standard deviation. Thus here, it can represent a variation of the buckling load up to about 8.95 × 10 6 N and 1.90 × 10 7 N for P b 1 and P b 2 around their average value, respectively, which corresponds to the the maximum of (2σ i ) on the two-dimensional deterministic space inputs (i.e. the couple values L spring /L beam , k spring ).
It is worth saying a few words here about the computational time related to the approach. First, considering the number of evaluation of the model, with the PCE-kriging meta-model, 2550 evaluations have been done to construct the surrogate model. To obtain the surface response of the reference presented in Figs. 6a, c and 7a, c, 26 × 51 × 1000 = 1326000 model evaluations have been performed. Thus, using this meta-model has divided the computational cost by 520. Considering the computational time, realizing 1000 simulations with the full model takes 30.2528 s, and so the evaluation of the 26 × 51 grid takes 4.0115 × 10 4 s=11.14 h. Once the PCE-kriging meta-model is constructed, the construction of the surface response for the same grid takes 0.7770 s as no LHS is required with the formulation. The computational time is then divided here by 5.1628 × 10 4 . The computational interest of using such approach is clearly demonstrated here, even on a small model whose evaluation time is small. This is mainly due to the fact that all the stochastic information is retained in the PCE coefficients parameterized with the kriging, avoiding the realization of costly MCS or LHS to get the stochastic moments.

Sensitivity analysis
To go further in the analysis, the Sobol indices are directly extracted from the PCE coefficients, they are given in Fig. 8 for both buckling loads, and the name of each variable is given at the top of each sub- figure. Only the first-order Sobol indices are displayed here as the total order Sobol indices appear to be equal. This demonstrates the low coupling between the different parameters, and justify the use of a hyperbolic norm as a truncation scheme for the PCE here. As a complement, in Table 5 the minimum and maximum values observed on the grid for each first-order Sobol index are given for each buckling load. The Sobol indices vary a lot from one parameter to another and depends on the considered buckling load. As a reminder, when the Sobol index of a parameter i is close to 1, then this parameter has a large influence on the output variance, whereas if the Sobol index is near 0, then this parameter i has a low influence. Sobol indices are a convenient tool to evaluate and compare the influence of different parameters on a quantity of interest. Considering the first buckling load, the outer height h o has the highest influence as its Sobol index remains always constant near 0.7. A jump to 0.85 around L spring /L beam = 0.5 and k spring is however observed, but it corresponds to the area where the error of the meta-model is the highest and so this higher value has to be taken with attention. The second parameter that has the highest influence is the beam length with a Sobol index equal to 0.2, which is much lower than the outer height of the beam. The other parameters have low Sobol indices values, showing their low influence on the first buckling load. Thus, if one wants to reduce the variance on the buckling load, it has to control the uncertainty on the beam length and its outer height in priority.
Considering the second buckling load, similar analysis can be done. The parameter with the highest influence is the outer height h o with a Sobol index of about 0.7. The latter increases up to 0.8 when the buckling load reaches its threshold value. The length of the beam is the second parameter that influences the second buckling load, with a lower influence as the Sobol index is only equal to about 0.2. The latter drops to about 0.02, having an evolution complementary to the Sobol index of h o .
Sobol indices are usually used to compare and assess the influence of each parameter on the output variance. By ranking the influence of the random parameters on the buckling loads here, the following conclusion is drawn: -high influence: it corresponds to the parameters with the highest Sobol indices, here h o as the Sobol indices are always superior to 0.62 for both buckling loads. -medium influence: L beam can also be considered as moderately influential as the Sobol indices can reach up to 0.2 in some cases. However, L beam reaches low values in some cases (see minimum equal to 0.0092 and 0.0011) so its influence depends on the spring location and stiffness. -low influence: all the other parameters, namely E, b o , b i and h i , have a low influence as their Sobol indices are always inferior to 0.1 and much lower than the other Sobol indices. Their effect can even be neglected based on this analysis.

Conclusion
In this work, the efficiency of an advanced computational technique based on kriging and Polynomial Chaos Expansion for structural stability of mechanical systems with uncertainties is illustrated. This methodology is more specifically tested in order to predict the buckling behaviour of a beam subjected to random geometrical and material properties and to a deterministic axial compression. In addition, the system includes a punctual spring whose location and stiffness vary, which drastically changes the values of buckling loads. The use of a surrogate model allows reducing the computational burden and avoid a computationally expensive MCS/scanning method. It is demonstrated that the method predicts accurately the buckling loads and reduces drastically the number of model evaluations and the computational time. Finally, it is also illustrated that a sensitivity analysis, without additional computational cost, can be performed to rank the influence of different random parameters. material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.