RELIABILITY ASSESSMENT USING A COMBINATION OF POLYNOMIAL CHAOS AND SIMULATIONS : APPLICATION TO NONLINEAR FRACTURE MECHANICS

This paper presents a probabilistic approach based on polynomial chaos expansion, in order to provide accurate explicit approximation of the structural response to be considered in the limit state function. The main difficulties in this approach are related to the calculation of the expansion coefficients which are defined by multi-dimensional integrals. As an alternative to the quadrature methods, Monte-Carlo simulations based on low discrepancy Halton sequence have been used for this issue. The accuracy and the efficiency of the proposed approach have been approved through analytical models. It is shown that the use of low discrepancy sequence provides more rapidly converging estimates. The proposed approach has been applied to assess the integrity of a cracked pipe.


INTRODUCTION
The engineering experience has shown that fracture mechanics is one of the main causes of failure of structures and mechanical components.Form the deterministic point of view, modeling the mechanical behavior of cracked components is not a trivial task, as it involves large singularities in the structural response.The progress of finite element analysis offers solutions for wide varieties of structural behaviors with various complexity levels, allowing for accurate predictions of the integrity of cracked structures and consequently the remaining resistance capacity.The failure mode of cracked structures can be analyzed by applying the deterministic fracture mechanics theory.Although this theory is widely used, the obtained predictions are conservative since the system parameters are considered as deterministic.However, it has been widely recognized that uncertainties can strongly arise in material proprieties, crack geometry and loading parameters.For this reason, stochastic approaches have been developed to take account for various sources of uncertainties in the assessment of structural performance.For instance, the applied reliability approaches are mainly based on second moment reliability methods, such as FORM/SORM [1].In these methods, a closed expression is often required to define the limit state function.Unfortunately, in engineering practice, the limit state function cannot be defined by explicit function of the basic random variables, and only implicit structural response can be available when numerical methods are involved (i.e.finite element model).Although the sampling methods [2] such as Monte-Carlo simulations are very useful, they require high computation effort, especially when the finite element model is complex and the failure probability is low.To cope with these difficulties, the Response Surface Method (RSM) [3] has been developed in the past decades, in order to reduce the number of calls to the mechanical model.The main idea of this approach is to build a second order polynomial approximation of the limit state function based on a limited number of realizations of the implicit function.This polynomial approximation is refined during an iterative process in order to obtain the best closed form.However, the obtained solution is only accurate in the vicinity of the most likely failure point.Consequently, the accuracy is often reduced when estimating the failure probability.In the literature, several RSM approaches have been proposed [4][5][6][7], but all of them have the same backbone, as the main differences lie in the choice of the experiment design points and the convergence criteria.As an alternative to these approaches, the present work proposes a stochastic response surface approach based on the polynomial chaos expansion [8] of the limit state function in the standard space.This approximation is defined as a projection in the basis of multidimensional Hermite polynomials.The unknown coefficients of the polynomial chaos expansion, which are defined by multidimensional integrals, are computed using Monte-Carlo simulations based on low discrepancy random number sequences, which are more efficient than classical pseudo-random numbers [9,10] (i.e. the sampling points are more uniformly distributed in the random space).The specificity of this approach is that the explicit function is constructed by taking into account the whole available statistical information of the input parameters.Consequently, the obtained approximation is close to the target limit state function not only in the vicinity of the most likely failure point, but also in the whole random space.The efficiency and the accuracy of the proposed approach are evaluated through an example dealing with nonlinear fracture mechanics.This paper is divided into three sections.In the first one, the mathematical formulation of the polynomial chaos expansion are presented for implicit response, and completed by the computation of multidimensional integrals using quasi Monte-Carlo simulations.In the second section, the accuracy and the efficiency of the proposed approach are shown through an analytical test model.The third section is dedicated to perform the sensitivity analysis on a cracked pipe under internal pressure, and then to estimate its failure probability.

Construction of the Stochastic Response Surface (SRS)
Let us consider a mechanical model having N uncertain parameters gathered in the vector ൌ ሺ‫ݔ‬ ଵ , … , ‫ݔ‬ ே ሻ.Mathematically speaking, the model response can be represented by the function:

‫ݕ‬ ൌ ݂ሺሻ ሺ1ሻ
where ‫ݕ‬ is the model response, which is considered to be scalar throughout this paper without loss of generality, obtained by explicit (i.e.analytical equations) or implicit (i.e.finite element analysis) representation of the function ݂.Let ܻ ൌ ݂ሺሻ be the stochastic model associated with the model defined by equation (1), where is a N-dimensional random variable modeling the uncertainties in the input parameters , and ܻ a scalar random variable representing the response ‫.ݕ‬In practical problems, the components of the N-dimensional random variable have different probability distributions and can also be correlated.This difficulty can be overcome by applying probabilistic transformations [11,12], which allow us to represent the N-dimensional random variable by N-dimensional standard Gaussian variable (i.e. the components are independent Gaussian variables with zero mean and unit standard deviation).Therefore, the ‫‬ ௧ order polynomial chaos expansion of the random variable ܻ is defined by [8]: where ܶ denotes the probabilistic transformation, ሼܽ ሽ ୀ ିଵ are the unknown deterministic coefficients of the polynomial chaos expansion, and ሼߖ ሺሻሽ ୀ ିଵ are multivariate Hermite polynomials which form an orthonormal basis with respect to Gaussian probability density.The number of terms ܲ in the above summation depends on the stochastic dimension N and the order ‫‬ of the polynomial chaos expansion; it is given by: The expansion defined by equation ( 2) converges in the sense of the ‫ܮ‬ ଶ -norm, that is: where ॱሾ.ሿ is the expectation operator.
As the polynomial chaos expansion is defined in the standard random space, the Ndimensional Hermite polynomial ሼߖ ࢻ ሺܷሻ, ࢻ ‫א‬ Գ ே ሽ can be constructed as the tensor products of N one-dimensional Hermite polynomials ‫ܪ‬ ఈ : where ࢻ is a N-dimensional index, whose components ሼߙ , ݅ ൌ 1, … , ܰሽ represent the orders of the one-dimensional Hermite polynomials ൛‫ܪ‬ ఈ ൟ ୀଵ ே .To build the polynomial chaos expansion of the random variable ܻ, we just need to compute the unknown coefficients ሼܽ ሽ ୀ ିଵ .Classically, two classes of approaches are distinguished for the computation of the expansion coefficients: the intrusive and non-intrusive approaches.In the first one, the unknown coefficients are obtained through the minimization of the residual under the constraint that this later is orthogonal in the polynomial chaos basis.This procedure is often called Galerkin projection scheme, which usually requires the modification of the numerical code.When the mathematical model involves nonlinearities, this procedure can be a challenging task and difficult to implement.To overcome these difficulties, the non-intrusive approaches have been developed.They only require the computation of the deterministic model at some realizations of the input parameters.Among the various types of non-intrusive methods, the projection methods are the most widely applied in the literature.In this paper, Monte-Carlo simulations based on low discrepancy random sequences are proposed for the computation of the expansion coefficients.

Computation of the PCE coefficients by simulations
The projection methods use the orthogonality property of the polynomial chaos basis to define the unknown coefficients using N-dimensional integrals.The projection of the polynomial chaos expansion (Eq.( 2)) on the polynomial basis ሼߖ ࢻ , |ࢻ| ‫‬ሽ, with respect to the inner product ‫,݃ۃ‬ ‫ۄ݄‬ հ ॱሾ݃, ݄ሿ gives: which can be written as: where the expectation ॱሾߖ ࢻ ଶ ሿ can be computed analytically by using the orthogonality of the one-dimensional Hermite polynomials, and ߮ ே denotes the probability density of Ndimensional standard Gaussian variable.To simplify the notations, let us consider the following N-dimensional integral: It is clear that the major difficulty in computing the expansion coefficients lies in the evaluation of the integral ‫ܫ‬ ࢻ ே .To this end, it is proposed herein to use Monte-Carlo simulations for which the integral ‫ܫ‬ ࢻ ே can be estimated by the following sum: where ൛ ൟ ୀଵ ெ are sample points generated according to N-dimensional Gaussian density.Monte-Carlo simulations are robust and converge for any ‫ܮ‬ ଶ -function.The mean square error ߝ ெ of the integral estimate ‫ܫ‬ ࢻ ே is given by: with ॽሾ.ሿ the variance operator.This expression shows the convergence rate of Monte-Carlo simulations, which is only proportional to 1 ‫ܯ√‬ ⁄ .This low convergence rate is mainly due to the use of pseudo-random number generators for which the sampling points ൛ ൟ ୀଵ ெ are not uniformly distributed in the random space.To overcome this problem, other sampling schemes have been proposed in the present work, the low discrepancy works [13,14] have shown significant Schlier [14] has demonstrate simulations based on low discrepancy sequences is convergence rate achieved when using pseudo such low discrepancy sequences are known as these methods have been used in Dai and work, the focus is put on quasi Monte reveals to be efficient in computing The one-dimensional Halton sequence is a part prime number , the components of the one expanding the integers where are integers of the base We can observe that the terms of Halton sequence are distributed in the range Moreover, the terms are separated by schemes have been proposed in the literature, such as, Latin hypercube sampling.low discrepancy sequences appear as an interesting tool significant efficiency in performing N-dimensional integration.demonstrated that the convergence rate associated low discrepancy sequences is which is much faster than achieved when using pseudo-random numbers.ncy sequences are known as quasi Monte-Carlo simulations.have been used in Dai and Wang [10] to perform reliability analysis.Figure 1 shows the distributions of 500 sample points, obtained cube and low discrepancy Halton sample points.constructed by random pairing N one-dimensional Halton [9] has proved prime numbers has a discrepancy , the good uniformity of the sample points in the random space.It is to note that in and consequently the the dimensions, the prime cycles can be very large.dimensional Halton sequence, distributed in the unit hypercube ሾ0,1ሿ ே .It is necessary to transform the later into ‫ܯ‬ realizations ሼ ଵ , … , ெ ሽ of N-dimensional standard Gaussian variable .This can be achieved by the following change of variables: where Φ and ‫ܨ‬ are respectively the cumulative distributions of standard Gaussian variable and uniform variable over ሾ0,1ሿ.The estimation of the unknown coefficients ሼܽ ሽ ୀ ିଵ of the polynomial chaos expansion is obtained by injecting the ‫ܯ‬ realizations ሼ ଵ , … , ெ ሽ in equation (9).Due to the deterministic nature of Halton sequence, it is not possible to obtain the confidence bounds for the estimates.As mentioned above, Monte-Carlo simulations based on low discrepancy sequences are generally more efficient than those based on pseudorandom numbers.However, the efficiency might be considerably reduced for high dimension problems (i.e.ܰ 10).This behavior has been explained by Caflish et al. [16] by introducing the notion of effective dimension, which has also been demonstrated by Kucherenko et al. [15] on the basis of Sobol's sensitivity indices.

Computation of the statistical moments
The statistical moments of the random variable ܻ (i.e.representing the model response) can be analytically obtained by the polynomial chaos expansion, using the orthogonality of Hermite polynomials.The first four statistical moments, namely the mean, the standard deviation, the skewness and the kurtosis respectively, are computed by: where ॱሾ.ሿ and ॽሾ.ሿ denotes resepectively the expectation and the variance.The expectations ॱൣߖ ࢻ ߖ ࢼ ߖ ࢽ ൧ and ॱൣߖ ࢻ ߖ ࢼ ߖ ࢽ ߖ ൧ can be easily derived by using the orthogonality of Hermite polynomials.Note that, most of the terms in the above summand are nil, which greatly simplifies the computations.Having the polynomial chaos expansion of the random variable ܻ, the associated probability density function can be easily constructed by Monte-Carlo simulations.This procedure is accurate and efficient since it requires very low computing time.

VALIDATION EXAMPLE
This analytical example aims at validating the efficiency and the accuracy of the proposed approach, in estimating the statistical characteristics (i.e.statistical moments and probability density) of the model response.Let us consider a three order polynomial function with input Gaussian variables ൌ ሼܺ ሽ ୀଵ ଷ having identical mean ߤ ൌ 1 and standard deviation ߪ: The accuracy of the statistical moments is investigated for Monte-Carlo simulations based on different sample scheme: pseudo-random numbers, Latin Hyper-cube sampling and Halton low discrepancy sequence.The global accuracy of the polynomial chaos expansion is measured by the lack of fit ߝ ைி [17], which is the normalized error in the mean square sense.
In other words, ߝ ைி allows us to evaluate the capacity of the polynomial chaos expansion to reproduce the target model; it is defined as: where ‫ݕ‬ ெ, is the ‫‬ ௧ order polynomial chaos expansion of the response, whose coefficients are obtained by Monte-Carlo simulations based on ‫ܯ‬ sample points.The computation of the quantities in equation ( 19) is a heavy task when the model is not available in explicit form.Therefore, the lack of fit ߝ ைி is estimated using Monte-Carlo simulations, as following: where ߪ ಾ, ଶ is the standard deviation computed using equation (15).
In the cases of moderate uncertainties (ߪ ൌ 0.3) and large uncertainties (ߪ ൌ 0.6), respectively, Figures 2a and 2b show the evolution of the estimate of the lack of fit ߝ̂ ைி in terms of the number of simulations, according to the three sample schemes.It can be noted that, the numerical computation is carried out by adopting 2 nd order polynomial chaos expansion, in order to reproduce the exact model.Consequently, the approximation error is only due to the integration errors when computing the coefficients of the polynomial chaos expansion.These two figures show that Monte sequence converge faster than those ba samples.Moreover, the convergence rate is not affected by other words, the approximation error is practically the same for both moderate and large input uncertainties.Note that Latin Hyper numbers, especially when dealing with large input uncertainties.evolution of the stochastic response surface obtained by polynomial chaos expansion (Eq.( 2)) with respect to each uncertain paramete schemes.show that Monte-Carlo simulations based on Halton low discrepancy than those based on pseudo-random numbers and the convergence rate is not affected by the input uncertainty , the approximation error is practically the same for both moderate and large input atin Hyper-cube samples are better performing than pseudo especially when dealing with large input uncertainties.Figure 3 evolution of the stochastic response surface obtained by polynomial chaos expansion (Eq.( 2)) ect to each uncertain parameter in the range , and for different sample onvergence of the stochastic response surface using different sample schemes As can be seen, the analytical approximation converges rapidly towards the target defined by equation ( 18), when low discrepancy Halton sequence is used, which confirms the 2a and 2b.It is now required to compute the first four statistical moments of the model one hand, the exact solution is obtained by analytical computation the other hand, the statistical moments of the random variable are for different samples Halton low discrepancy and Latin Hyper-cube the input uncertainty level.In , the approximation error is practically the same for both moderate and large input than pseudo-random Figure 3 depicts the evolution of the stochastic response surface obtained by polynomial chaos expansion (Eq.( 2)) and for different sample using different sample schemes rapidly towards the target solution , which confirms the model response, for analytical computation, are estimated using equations ( 14) to (17).As can be solution.The low discrepancy Halton sequence provides more accurate estimates of the statistical moments than the other methods not affected when dealing with large input uncertainties (i.e.The probability density function of the random variable On one hand, the exact solution is obtained by crude Monte equation ( 18), using 10 6 sample function are constructed using 10 figure 4 that the probability density function obtained by Latin Hyper discrepancy Halton sequence are much closer to the exact solution.Consider the axi-symmetrical pressure and axial tension has an internal radius centered circumferential internal crack with length (17).As can be seen in Table 1, all the methods converge he low discrepancy Halton sequence provides more accurate estimates of the than the other methods.Moreover, it can be observed not affected when dealing with large input uncertainties (i.e. ).
: The probability density function of the random variable is estimated by various methods.On one hand, the exact solution is obtained by crude Monte-Carlo simulations applied to samples.On the other hand, the estimates of the probability density function are constructed using 10 6 runs of the polynomial chaos expansion figure 4 that the probability density function obtained by Latin Hyper-cube sampling and low crepancy Halton sequence are much closer to the exact solution.

STOCHASTIC ANALYSIS OF CRACKED PIPE
cracked pipe depicted in figure 5 [18], and subjected to internal and axial tension .The pipe of length , a thickness , and a circumferential internal crack with length .
methods converge to the exact he low discrepancy Halton sequence provides more accurate estimates of the that the accuracy is comparison of the statistical moments  The pipe is made of steel, whose strain relationship [19]: where is the stress, is the strain, respectively, stress, is a dimensionless material parameter stress-strain curve associated with Figure There are several ways to characterize the stress singularity in tip.For instance, when the material crack driving forces.For elastic characterized by the J-integral crack growth is observed when the : Geometry and applied loads of the cracked pipe [18] Due to the boundary conditions at the pipe ends, the internal pressure longitudinal tension in addition to the applied axial tension .The resulting stress the end effects reads: steel, whose behavior follows the well-known Ramberg strain, respectively, is the Young's modulus, imensionless material parameter and is a strain hardening exponent.
with Ramberg-Osgood constitutive law is shown in figure 6.There are several ways to characterize the stress singularity in the neighborhood of the crack tip.For instance, when the material is elastic, the stress intensity factors are often tak elastic-plastic material behavior, the crack driving force is integral, whose formulation is based on energy considerations when the J-integral exceeds the material toughness [18] the internal pressure produces a he resulting stress Ramberg-Osgood stressthe Young's modulus, is a reference is a strain hardening exponent.The law is shown in figure 6.
neighborhood of the crack elastic, the stress intensity factors are often taken as the the crack driving force is , whose formulation is based on energy considerations.The exceeds the material toughness .In general, the exact solution of the J-integral is only known for academic problems.For engineering problems, the computation of the J-integral is only possible through numerical procedures.In order to predict the behavior of the cracked pipe, a finite element model has been developed using the software Cast3m [21].Due to the symmetry of the problem, only one half of the pipe has been modeled.The mesh is composed of 709 quadratic elements and 1553 nodes, as depicted in figure 7.
Figure 7: Mesh of the cracked pipe

Sensitivity analysis
At first, the statistical moments and the probability density function of the mechanical responses, calculated by the J-integral, are evaluated.In this example, the uncertain parameters are the Young's modulus ‫,ܧ‬ the reference stress ߪ ௬ and the two parameters ߙ and ݊ of the Ramberg-Osgood law.These variables are considered as independent, with the parameters given in table 2. As can be seen, accurate estimates of the four moments are obtained when using low discrepancy Halton sequence for computing the expansion coefficients.The estimates of the mean are in good agreement with the reference solution.The higher relative error is observed on estimates of the kurtosis, which is only about 11.76%.As expected, the accuracy of the proposed method decreases with the order of the statistical moment.However, when using pseudo-random numbers to compute the unknown coefficients of the polynomial chaos expansion, the proposed approach failed to estimate the statistical moments.The accuracy is more affected for statistical moments of higher order.The relative error for the skewness is ߝ ൌ 97.4%.In general, the proposed approach offers the best compromise between efficiency and accuracy.

Reliability analysis
Under the applied tensile load ߪ ௧ , the failure occurs when the crack driving force ‫ܬ‬ exceeds the strength ‫ܬ‬ ூ representing the material toughness.The limit state function reads: In addition to the uncertain parameters used in the sensitivity analysis, the material toughness ‫ܬ‬ ூ is modeled as lognormal random variable with mean ߤ ൌ 59 ‫.ܽܲܯ‬ ݉݉ and standard deviation ߪ ൌ 9.5 ‫.ܽܲܯ‬ ݉݉.The failure probability ܲ is estimated by Monte-Carlo simulations applied to the 3 rd order polynomial chaos expansion of the limit state function.
The coefficients of the polynomial chaos expansion are estimated by Monte-Carlo simulations with 10000 samples generated from low discrepancy Halton sequence.The estimates of the failure probability for various levels of the applied tensile load ߪ ௧ are given in table 4. As expected, we can see that the failure probability increases with the applied tensile load ߪ ௧ .An increase of 25% of the applied tensile load ߪ ௧ produces an increase of the failure probability by a factor of 10 5 .The role of the applied tensile load ߪ ௧ is therefore very important for assessing the integrity of the cracked pipe.

CONCLUSION
In this paper, a polynomial chaos expansion-based probabilistic approach is presented to perform sensitivity and reliability analyses of structures.The originality of the approach lies in computing the expansion coefficients by using quasi Monte-Carlo simulation method.This method is based on low discrepancy sequence which ensures a better uniformity of the sample points in the random space than the pseudo-random numbers classically used in Monte-Carlo simulations.On the basis of the validation example, the low discrepancy Halton sequence improves the accuracy and the efficiency of the Monte-Carlo simulations, compared to Latin Hyper-cube samples and pseudo-random numbers.Therefore, the statistical moments of the response can be accurately estimated with low computation cost.The proposed approach is then applied to assess the integrity of cracked pipe.The stress field singularity in the neighborhood of the crack tip and the residual fracture strength can be adequately characterized by the J-integral concept applied to the finite element model with Ramber-Osgood constitutive law.It is shown that the applied tensile load is the most important parameter for the integrity of the cracked pipe.In the future, it would be interesting to use other low discrepancy sequences to enhance the efficiency of the Monte-Carlo simulations in particular when dealing with higher stochastic dimensions.

Figure 1 :
Figure 1: Comparison of the distribution of 500 sample points obtained by numbers generator, Latin Hyper schemes have been proposed in the literature, such as, Latin hypercube sampling.low discrepancy sequences appear as an interesting tool significant efficiency in performing N-dimensional integration.demonstrated that the convergence rate associated low discrepancy sequences is which is much faster than achieved when using pseudo-random numbers.ncy sequences are known as quasi Monte-Carlo simulations.have been used in Dai and Wang[10] to perform reliability analysis.onquasi Monte-Carlo simulation based on Halton sequence which in computing high-dimensional integrals.dimensional Halton sequence is a particular low discrepancy sequence the components of the one-dimensional Helton sequence in the base ; the term of the sequence

Figure 3 :
Figure 3: Convergence of the stochastic response surface

Figure 5 :
Figure 5: Geometry and applied loads of the cracked pipe

Tableau 4 :
Cracked pipe problem -failure probability in terms of the applied tension.
[20]first four statistical moments of the crack driving force are computed from the coefficients of the 3 rd order polynomial chaos expansion, which are estimated by Monte-Carlo simulations based on 10000 samples generated by low discrepancy Halton sequence.Note that the reference solution is obtained by Lagrange polynomial-based stochastic collocation method (L-SRS)[20].The estimates of the statistical moments when using pseudo-random numbers and Halton sequence, are reported in table 3 together with the reference values.
Tableau 2 : Cracked pipe problem -statistical characteristics of the uncertain parameters