Uncertainties in structural dynamics : overview and comparative analysis of methods

When studying mechanical systems, engineers usually consider that mathematical models and structural parameters are deterministic. However, experimental results show that these elements are uncertain in most cases, due to natural variability or lack of knowledge. Therefore, engineers are becoming more and more interested in uncertainty quantification. In order to improve the predictability and robustness of numerical models, a variety of methods and techniques have been developed. In this work we propose to review the main probabilistic approaches used to model and propagate uncertainties in structural mechanics. Then we present the Lack-Of-Knowledge theory that was recently developed to take into account all sources of uncertainties. Finally, a comparative analysis of different parametric probabilistic methods and the Lack-Of-Knowledge theory in terms of accuracy and computation time provides useful information on modeling and propagating uncertainties in structural dynamics.


Introduction
Industrial structures are mainly assemblies of many parts with complex geometries and non-linear characteristics.In most cases, parameters of the mathematical-mechanical model linked to geometry, boundary conditions and material properties can neither be identified nor modeled accurately.For that, there has been a growing interest in modeling uncertainties in mechanical structures [1][2][3].
In computational mechanics, random uncertainties are due to data uncertainties and model uncertainties.Data uncertainties concern the parameters of the mathematical-mechanical model and can be taken into account by parametric approaches.Table 1 gives a quick overview of some of the main probabilistic and nonprobabilistic parametric methods.
Statistical methods, such as Monte Carlo methods [4], seem to be the simplest way to evaluate the response of a a Corresponding author: sami.daouk@lmt.ens-cachan.frmodel with uncertain parameters but their convergence is time-consuming.Sampling methods and variance reduction techniques help overcome this difficulty.Stochastic Finite Element methods result from the combination of finite elements with probabilities.The main variants are the perturbation method [5], and the spectral stochastic finite element method [6].Model reduction techniques aim to approximate the response of a complex model by the response of a surrogate model built through a projection on a reduced basis, such as the Proper Orthogonal Decomposition [7] and the Proper Generalized Decomposition [8].Nevertheless, uncertainties are sometimes epistemic due to imprecise or incomplete information, and that is why non-probabilistic approaches were developed [9].Among the first developments was the interval theory, usually applied when the focus lies in identifying the bounds of the interval in which the output of interest varies, such as to problems in mechanics [10,11] and robotics [12].The fuzzy set theory is another nonprobabilistic approach that uses the fuzzy logic which is based on human reasoning rather than rigid calculations [13,14].More recently, the info-gap decision theory was developed in order to solve decision problems under severe uncertainty [15].Other non-probabilistic techniques, less developed in structural mechanics, are used to model uncertainties, such as the generalized information theory [16] and the evidence theory, also called Dempster-Shafer theory [17].
Due to the introduction of simplifications in the modeling procedure, model uncertainties should be considered but cannot be taken into account by the parameters of the mathematical-mechanical model.A description of the nonparametric probabilistic approach for dynamical systems is presented in [18,19].Model uncertainties are considered in a global way when considering the matrices of the dynamical model as random matrices, built using the maximum entropy approach.Recently, efforts have been made to develop a modeling technique taking into account different sources and types of uncertainty.The Lack-Of-Knowledge (LOK) theory [20,21] was proposed in the intention to qualify and quantify the difference between numerical models and real structures.
The aim of this paper is to provide engineers with practical information on the main stochastic methods used to model and propagate uncertainties in structural mechanics.Firstly, an overview of stochastic parametric and nonparametric approaches is presented.Then, a comparative analysis of different parametric probabilistic methods and the LOK theory in terms of accuracy and computation time will provide useful information on modeling and propagating small uncertainties in structural dynamics.

Reference problem and notations
For many physical problems studied by engineers, the conceptual model can be written in terms of stochastic partial differential equations.Uncertainties are then modeled in a suitable probability space and the response of the model is a random variable that satisfies a set of equations formally written where A is a differential operator, b is associated with the source terms and ζ = ζ(θ) are the uncertain parameters modeled as random variables, θ denoting randomness.For instance, when analyzing free vibrations of a structure, Equation (1) becomes the eigenvalue problem of linear dynamics.

Monte Carlo simulation
Monte Carlo simulation (MCS) is the basic approach that characterizes the response variability of a structure through the use of repeated deterministic experiments [4].After defining a variation domain of the inputs (i.e.uncertain parameters), n pseudo-random samples of variance σ 2 are generated from a probabilistic distribution over the domain.For each sample, a deterministic structure is defined and computed in the framework of the FE method.Finally, a statistical study of the output of interest provides first and second moments, probability density function, and failure probability.
Easy to implement and robust, the convergence rate of MCS is in the order of σ/ √ n.Parallel computation and processing can accelerate the convergence [22].Variance reduction and sampling techniques were also proposed in order to reduce the computational effort needed to obtain accurate results.

Variance reduction techniques
The optimization of the location of the samples is a first way to improve MCS.Variance reduction techniques [23], such as importance sampling, are based on concentrating the distribution of the generated samples in zones of most importance, namely the parts which mainly contribute to the statistical estimate of the output of interest.

Sampling techniques
While MCS uses a pseudorandom sequence, Quasi-Monte Carlo methods consist in choosing the samples using low-discrepancy sequences (e.g.Halton and Sobol sequences), also called quasi-random sequences [24,25].That yields a convergence rate in the order of 1/n.Latin Hypercube Sampling allows a reduction of the number of samples required [26].For each of the p variables of the problem, its range is divided into N intervals, in each one sample is generated.The N samples of each variable are then randomly matched together until obtaining N ptuples.The number N of p-tuples is significantly smaller than n with the same accuracy level.

Reliability methods
The aim of reliability methods [27] is to evaluate the probability of failure of a system containing uncertainties.Let S = S(ζ) denote a vector of load effects (e.g.displacements, stresses, etc.) that is related to the uncertain parameters and whose values define the failure of the system.To assess the reliability of a structure, a limit state function g is defined, such as g(S) ≤ 0 defines the failure state.If f S (s) denotes the joint probability density function (PDF) of S, then the probability of failure is given by: Computing this PDF is not possible with MCS due to the large number of samples required for an accurate evaluation (10 m+2 samples for a precision of 10 −m , with m ≥ 4).In addition, f S (s) is usually unknown.To overcome these difficulties, a probabilistic transformation Y = Y(ζ) is applied (e.g. the Nataf or Rosenblatt transformations [27]).The limit state function is mapped onto the standard normal space by using Y and the failure probability can be rewritten as:

ϕ(y)dy
where ϕ(y) denotes the standard normal PDF.Due to the form of ϕ(y), the most likely failure point in the standard normal space is the nearest one from the origin, and is the solution of the following optimization problem: After computing y * , reliability methods consist in approximating the failure domain by a simpler one whose probability can be estimated analytically.The First Order Reliability Method (FORM) replaces the integration domain by the hyperplane that passes through the design point and which is orthogonal to vector y * .A better approximation can be provided by the Second Order Reliability Method (SORM) when replacing the limit state surface by a quadratic surface whose probabilistic content is known analytically.Lately, a more efficient and accurate approach, the Mean-Value First Order Saddlepoint Approximation (MVFOSA) method, was proposed [28].

Stochastic finite element methods
The combination of the probability theory with the FE method resulted in the development of stochastic finite element methods [2,29,30].The FE method is used in the spatial discretization.When the uncertain parameters are random but spatially constant, they are modeled as random variables.To account for the spatial variability of uncertain quantities (e.g.material property), a characterization in terms of a random field is usually employed.Through a process of discretization, it is possible to represent the random field by a vector of random variables.The discretization methods can be divided into three groups, i.e. point discretization methods, average discretization methods and series expansion methods [31].Stochastic FE methods are based on a series representation of the solution.

Perturbation method
The perturbation method is one of the most popular techniques used for analyzing random systems in stochastic engineering [5,32].After the discretization step, the solution of ( 1) is approximated by a Taylor series expansion around the mean value μ ζ of the random variables {ζ k (θ)} 1 k n as follows: where

By writing the operator A(•, ζ(θ)
) and the right-hand side b(ζ) similarly and by injecting them in Equation ( 1), one obtains that the coefficients of the expansion of u are solution of a sequence of deterministic problems.
Due to limitations in numerical differentiation, usually expansions of degree two or smaller are used, and only the first and second order statistics are computed.Therefore, the perturbation method is accurate in the case of small uncertainties.The mean and covariance can be simply expressed in terms of the expansion coefficients and the statistics of the random variables: where C ζiζj denotes the covariance of variables ζ i and ζ j .

Neumann decomposition
The idea of the Neumann decomposition method is to expand in a Neumann series the inverse of the random operator A(•, ζ) [33][34][35].It starts by writing the operator as: where A 0 is the mean of A, A −1 0 its inverse, and I denotes the identity operator.The inverse of A can be then expanded in a kind of geometric series, known as the Neumann series: The series converges only if the absolute values of all eigenvalues of A −1 0 A * are less than 1, which can always be satisfied after introducing a scalar factor as proposed in [34].The solution of (1) can then be represented by the series as: where the terms u i (ζ) are solutions of the following problems: The Neumann decomposition method has the advantage that only the deterministic part of the random operator has to be inverted and only once, and that partial derivatives with respect to the random variables are not required.However, computing the expansion terms of u requires to solve deterministic problems with random right-hand sides, which can be expensive.In practice, the series is truncated and the first moments of the approximated solution are evaluated.

Spectral stochastic finite element method
The Spectral Stochastic Finite Element method aims at discretizing the random dimension in an efficient way using two procedures.Firstly, the input random field is discretized using the truncated Karhunen-Loève expansion.Then, the solution is represented by its coordinates in an appropriate basis of the space of random variables, "the polynomial chaos" [6].
The Karhunen-Loève (KL) expansion is a representation of a stochastic field (or process) as an infinite linear combination of orthogonal functions, that are the eigenvalues of a Fredholm integral equation of the second kind, considering the autocovariance function as kernel.Considering a stochastic field h(x, θ) modeling an uncertainty in the system (e.g.material Young's modulus), this expansion is defined as follows: where x is a vector containing space variables, μ(x) is the mean of the field, λ i and Φ i (x) are the eigenvalues and eigenfunctions of the autocovariance function Cov[x 1 , x 2 ], and {ξ i (θ)} forms a set of uncorrelated random variables.If the random field h(x, θ) is Gaussian, then {ξ i (θ)} are independent standard normal random variables.Otherwise the joint distribution of {ξ i (θ)} is almost impossible to obtain.Hence, the KL expansion is mainly applicable to the discretization of Gaussian fields.
The solution u(θ) is then searched under the form of an expansion over the polynomial chaos: where {Ψ j (θ)} ∞ j=0 is a complete set of orthogonal random polynomials defined by means of {ξ i (θ)} ∞ i=1 forming an orthonormal basis, satisfying: The correspondence of the type of polynomial chaos {Ψ j (ξ i )} and the probability distribution of variables {ξ i } is discussed in [36] (e.g.Hermite polynomials are associated with the Gaussian distribution).
In order to evaluate the solution, the coordinates {u j } ∞ j=0 , that are deterministic vectors having as much components as the number of degrees of freedom, must be computed.Then a truncation of both expansions is performed: M + 1 terms in the Karhunen-Loève expansion and P terms in the polynomial chaos expansion are taken (P + 1 = (M+p)!M!p! , p is chaos polynomials order).Upon substituting the truncated expansion of the solution into the governing Equation (1) with ξ 0 (θ) ≡ 1, one finds: A Galerkin projection of Equation ( 5) onto each polynomial {Ψ j (θ)} is then performed in order to ensure that the error is orthogonal to the functional space spanned by the finite-dimensional basis {Ψ j (θ)} P −1 j=0 .By using the polynomial chaos expansion of u(θ), the randomness is effectively transferred into the basis polynomials.Thus, using orthogonality of the polynomial basis leads to a set of P coupled equations for the coefficients u j that are deterministic.The mean and covariance matrix can be then obtained, along with the probability density function of the components of u j : While mathematically elegant, this approach suffers from the curse of dimensionality.In addition, higher order approximations are necessary to more accurately capture higher order response statistics.To address such difficulties, a non intrusive method that builds a sparse PC expansion was introduced [37].Another work proposes a model reduction technique for chaos representations that mitigates the computational cost without reducing the accuracy of the results [38].

Model reduction techniques
Modeling any physical phenomenon on large domains results in inadequate computing resources, both in terms of speed and storage capacity.The methods described previously focus on the equations to solve rather than the model itself.Model reduction techniques seek to reduce the studied model through approximations and projections without modifying the mathematical/physical model.
The basic idea is to build an approximation of the solution, which is a multivariable function (e.g.space and structural uncertain parameters), as a finite sum of products of functions where the {ψ p (x)} p=1... n and {λ pk (ζ k )} p=1... n form low-dimensional reduced bases of deterministic vectors and stochastic functions of the random parameters {ζ k } k=1... n k respectively.Reduced bases can be defined in many ways leading to an optimal decomposition of the solution for a given order n of decomposition.A statistical study provides first moments and probability density distribution of the solution.

A posteriori techniques
A posteriori model reduction techniques build the reduced basis using some knowledge (at least partial) on the solution (called snapshot) of the reference problem.Then the coefficients in the linear combination (6) are computed by solving the equations of the initial problem projected in the reduced basis.The offline determination of the reduced basis usually has a high computational cost, while the cost of the online solving of the reduced model is significantly smaller.
A first popular technique is the proper orthogonal decomposition (POD) [7,39] where random values of the uncertain parameters are considered in order to evaluate the snapshots.The basis vectors are chosen among the snapshots technique such as the separated representation minimizes the distance to the solution with respect to a given metric.Another way to build the reduced basis is given by the reduced basis (RB) method [40], especially developed for parametrized problems.In this case, the uncertain parameters are chosen in an iterative optimal way to compute the snapshots.
The simplest approximation of the properties of dynamic systems is the projection on a truncated modal basis, also referred to as the subspace reduction (SR) method [41].An accurate approximation of the response is assumed to be found in a subspace spanned by the columns of a rectangular matrix T (N rows and N R N columns, where N is the number of degrees of freedom of the initial problem).The approximate eigenmodes of a model are given by Φ i = TΦ iR , where Φ iR are the eigenvectors of the projected (reduced) eigenvalue problem In practice, the columns of matrix T can be the truncated eigenmodes of several models, each given for different values of the uncertain parameters.This approach can be very useful when analyzing uncertainty propagation to some eigenfrequencies and eigenmodes of a model.

A priori techniques
A priori model reduction techniques [42] do not require any knowledge of the solution and operate in an iterative way, where a set of simple problems, that can be seen as pseudo eigenvalue problems, need to be solved.They are based on the proper generalized decomposition (PGD) [8,43] where the deterministic vectors and stochastic functions in ( 6) are initially unknown and computed at the same time.This computation phase is usually very expensive, but done only once.The solution of the problem for a variation of any parameter can be evaluated in a significantly fast and cheap online computation.The PGD has been successfully applied to stochastic problems [44,45] and to solve multidimensional models [46].

Nonparametric methods
In computational mechanics, random uncertainties in model predictions are due to data uncertainties and to model uncertainties.Data uncertainties concern the parameters of the mathematical-mechanical model and can be taken into account by parametric probabilistic approaches, such as the ones previously presented.However, for complex mechanical systems, the constructed model cannot be considered representative due to the introduction of approximations and because some details are not accurately known, such as in joining areas.Clearly, model uncertainties should be introduced but cannot be taken into account by the parameters of the mathematicalmechanical model under consideration.A description of the nonparametric probabilistic approach for dynamical systems can be found in [18,19].Model uncertainties are considered in a global way when considering the matrices of the dynamical model as random matrices, built using the maximum entropy approach.Developments for taking into account data and model uncertainties in linear dynamical systems are presented in [47].Recently, a generalized probabilistic approach that models both modelparameter uncertainties and modeling errors in structural dynamics, for linear and nonlinear problems, has been introduced and validated [48].

The lack-of-knowledge theory
Parametric probabilistic approaches are usually used to model intrinsic uncertainties resulting from natural variabilities.In case of epistemic uncertainties resulting from a lack of knowledge, it seems interesting to look for the solution in an interval rather than conducting a global stochastic search.The concept of lack-ofknowledge [20,21] is based on the idea of globalizing all sources of uncertainty for each substructure through scalar parameters that belong to an interval whose boundaries are random variables.

Basic lack-of-knowledge
The starting point of the LOK theory is to consider a theoretical deterministic model, to which a lack-ofknowledge model is added.Let Ω be a family of real and quasi-identical structures, each being modeled as an assembly of several substructures E. To each substructure E two positive scalar random variables m − E (θ) and m + E (θ) are associated, called basic LOKs and defined in [20,49] as follows: where K E and K E are the stiffness matrices of E, for the real structure and the theoretical deterministic model respectively.The basic LOKs are stochastic variables characterized by probability laws defined using the deterministic interval [−m − E ; m + E ]. Equation ( 7) may seem sometimes difficult to use because of the usual size of the stiffness matrix of industrial structures.Therefore the basic LOKs are expressed in practice using scalar quantities related to the stiffness, namely the strain energies: where e E (U, θ) = 1 2 U T K E (θ)U is the strain energy of a real substructure taken from the family of structures Ω, and e E (U) = 1 2 U T K E U is the strain energy of the theoretical deterministic model of the substructure.

Effective LOK of an output
From the basic LOKs, the first step is to establish a general procedure that leads to the evaluation of the dispersion of any variable of interest α.The second step is to write the difference: between the value α LOK of the variable of interest given by the LOK model, and α the one from the theoretical deterministic model.Then Δα can be expressed as a function of the stiffness or the strain energy.Using in Equations ( 7) and ( 8) leads to the propagation of the intervals ([m − E (θ); m + E (θ)]) E∈Ω throughout the stochastic model.In the LOK model, one associates to each generated sample of the basic LOKs two bounds Δα − LOK et Δα + LOK satisfying: where As long as the probability laws of the basic LOKs are known, the dispersion of these bounds can be determined.
After generating the basic LOKs, the probability density functions of Δα − LOK (θ) and Δα + LOK (θ) are found.When carrying out a comparative analysis, representative quantities of the dispersion of Δα are extracted from α LOK , namely Δα − LOK (P ) and Δα + LOK (P ), called effective LOKs on α LOK .These quantities are the bounds defining the smallest interval containing P % of the values of Δα.
As an application of this procedure, let the quantities of interest be the angular eigenfrequencies in the case of free vibrations where K = E∈Ω K E is the random global stiffness matrix, M = E∈Ω M E is the random global mass matrix, Φ i is the ith eigenvector and ω i the corresponding eigenvalue.When the basic LOKs are small enough to approximate, one writes the difference on the square of ω i as: where the modes are normalized with respect to the mass matrix.Through Equation (10), the definition (8) enables the propagation of the LOK intervals.Therefore the lower and upper bounds are found as: where In the case of small uncertainties, the effective LOKs on the eigenmodes Φ i can be similarly evaluated.Details are presented in [20] along with an extension of the LOK theory to consider the presence of a large basic LOK.

Applications and results
Among the methods cited previously, four were implemented and used to evaluate eigenvalues of problem ( 9) on two different assemblies.Each chosen method propagates uncertainty differently: MCS is a statistical approach, the perturbation method is based on the stochastic finite elements, the SR method is a model reduction technique and the LOK theory, new to engineers, considers different sources and types of uncertainty.These methods are lightly intrusive and have been applied to industrial problems.For the study cases that will be presented, the Young's modulus of the material of each substructure was assumed to be a random variable such as: where η is a uniform random variable defined in the interval [-1; +1].In the case of the LOK theory, the basic LOKs m E (θ) were considered as uniform random variables in [-0.08; +0.08].
For the statistical studies, 20 000 samples were generated and some eigenfrequencies of the global structure were evaluated for each sample.As a common representation, the bounds for the probability P = 0.99 were extracted from the PDFs given by each method.That means that the bounds define the smallest interval containing 99% of the values of the eigenfrequencies.

First assembly
In order to validate the implementation of the considered methods, a planar truss formed of four pin-jointed bars was considered, as shown in Figure 1.The connections between the structure and the base are assumed to be perfectly rigid.All beams are made of aluminum (E = 72 GPa, ρ = 2700 kg.m −3 ) with a circular crosssection of 10 −4 m 2 .
The deterministic structure has four tensilecompressive eigenmodes, and the associated eigenfrequencies were evaluated.The 99%-intervals presented in Table 2 show the influence of stiffness uncertainties on the eigenfrequencies of the assembly.It can be noticed that the intervals are nearly identical and that the ones given by the LOK theory are included in the ones resulting from Monte Carlo simulations.An important question can then be asked whether the same accuracy can be obtained when considering a larger assembly, representative of a real industrial structure.The next study aims to provide an answer in this matter.

Second assembly
The assembly shown in Figure 2 was inspired from the geometry of the booster pump studied in the framework of the international benchmark SICODYN [50,51].The main goal of this benchmark was to measure the effective variability on structural dynamics computations and then quantify the confidence in numerical models used either in design purpose or in expertise purpose and finally to ensure robust predictions.The structure is an assembly of a cone and a tube with a wedge between these two elements.The whole set is made of steel (E = 210 GPa, ρ = 7800 kg.m −3 ) and the largest side of the cone is clamped.The mesh is composed of 3240 finite elements (∼20 000 DOFs), which seems representative of industrial components.
In this case, the quantities of interest were the first three eigenfrequencies of the assembly.Figure 3 shows the shapes of the corresponding eigenmodes for the reference structure and Table 3 gives the 99%-intervals resulting from the statistical study of uncertainty propagation.
The Monte Carlo method was considered as reference for the comparative analysis of the methods in terms of quality of the results.For that purpose, Table 3 shows the highest algebraic percent error for both interval bounds and the central processing unit (CPU) computation time.
The results show that all implemented methods converge to similar results.The LOK theory seems to be less conservative than MC, leading to reduced intervals.In   addition, the LOK theory and the perturbation method reveal to be the best in terms of computation time with a substantial advantage compared to other methods.The more numerous the degrees of freedom of the studied model are, the greater this advantage will be.

Conclusion
This paper aimed to provide engineers with practical information on the main probabilistic methods used to model and propagate uncertainties in structural mechanics.A brief overview of probabilistic approaches and the LOK theory was presented.Four lightly intrusive methods, modeling and propagating parametric uncertainties differently, were used to evaluate eigenfrequencies of two different assemblies.In the case of small uncertainties on stiffness parameters, the comparative analysis presented accurate results with a significant difference in terms of computation time.MCS and model reduction techniques are time consuming.Stochastic FE methods provide results more quickly but are relatively intrusive.As for the LOK theory, it globalizes all sources of uncertainties, related to data and model, and seems to be less conservative.This modeling technique bounds the quantity of interest in a stochastic interval.The main advantage is the low implementation and computation time which can be handy when modeling real industrial assemblies.
In the framework of the SICODYN Project [50,51], initiated in 2012 and carried out till 2016, the comparative analysis will be extended to the case of a booster pump studied within its industrial environment.One of the goals is to evaluate the ability of some methods, such as the LOK theory and the generalized probabilistic approach of uncertainties [52], to quantify, not only data uncertainties, but also model uncertainties, especially in the case of high values of uncertainties.

Table 3 .
99%-intervals for the first three eigenfrequencies of the 3D assembly, with percent error and computation time (20 000 samples).