Uncertainty Quantification in Small-Timescale Model-Based Fatigue Crack Growth Analysis Using a Stochastic Collocation Method

Due to the uncertainties originating from the underlying physical model, material properties and the measurement data in fatigue crack growth (FCG) processing, the prediction of fatigue crack growth lifetime is still challenging. The objective of this paper was to investigate a methodology for uncertainty quantification in FCG analysis and probabilistic remaining useful life prediction. A small-timescale growth model for the fracture mechanics-based analysis and predicting crack-growth lifetime is studied. A stochastic collocation method is used to alleviate the computational difficulties in the uncertainty quantification in the small-timescale model-based FCG analysis, which is derived from tensor products based on the solution of deterministic FCG problems on sparse grids of collocation point sets in random space. The proposed method is applied to the prediction of fatigue crack growth lifetime of Al7075-T6 alloy plates and verified by fatigue crack-growth experiments. The results show that the proposed method has the advantage of computational efficiency in uncertainty quantification of remaining life prediction of FCG.


Introduction
Fatigue damage is among the most common failure modes in structural safety. It refers to the performance degradation process of structural materials due to cyclic loading. The fatigue damage process accumulates from small crack initiation, propagates to macro cracks and finally leads to structural failure under certain conditions. In practical engineering, the prognosis of fatigue damage plays a critical role in estimating an engineering system's remaining useful life based on the damage evolution model [1]. However, sources of uncertainties widely exist in the process of fatigue propagation analysis, including the variability in material properties, experimental errors, variation in loading conditions and model inaccuracy [2]. Therefore, the prognosis of the fatigue damage should take the form of a probabilistic distribution. This paper focuses on the remaining useful life prognosis taking into account the uncertainties in fatigue crack growth (FCG) analysis.
Numerous stochastic analysis methods have been proposed to deal with uncertainty factors observed in large replicate FCG tests and investigate the uncertainties of FCG analysis [3][4][5]. The Monte Carlo (MC) method, whose implementation is straightforward, has been the most commonly used approach. Liu and Mahadevan proposed a MC simulation-based method for predicting the probabilistic fatigue life [6]. However, the convergence rate of the MC method is slow due to the large number of samples needed, especially for a complex system. Several variants of the MC have been developed to improve the computational efficiency for probabilistic fatigue crack analysis, such as the Markov chain Monte Carlo method (MCMC) and subset simulation [7,8], which have application limitations due to respective additional restrictions. Next, stochastic spectral methods [9] have been been proposed for uncertainty quantification of complex multi-dimensional problems. The stochastic spectral methods are based on the polynomial chaos expansion technique and construct the explicit response surface. Hermite polynomials were usually used as the polynomial chaos basis for modeling the stochastic responses caused by Gaussian variables. In recent years, the generalized polynomial chaos (gPC) method has become the most popular stochastic spectral method, which is an analysis and quantification method for the non-Gaussian random input parameters by the differential equations [10]. Stochastic Galerkin method and stochastic collocation method are most commonly used to solve the coefficients of the polynomial expansion in gPC. The stochastic Galerkin method has many applications for stochastic problems in chemical systems and computational fluid dynamics [11,12]. However, this method is computationally cumbersome and can be inapplicable for some highly complex and nonlinear problems due to the large coupled set of equations. The stochastic collocation method achieves higher resolution for polynomial approximations in random space compared with the stochastic Galerkin method [10]. The functional interpolation by polynomials and integration of functions according to Gaussian quadrature rules are the basic ideas of the stochastic collocation method. In this method, a sparse grid is used for the construction of the set of the interpolation points. The stochastic collocation method has been widely applied in many fields in the last several years, e.g., uncertainty analysis of flow in random porous media [13], threedimensional problems in solid mechanics [14], reliability analysis of structures with stochastic loadings and material properties [15], etc. From the above literature review, this approach shows the ease of implementation in high-dimensional random spaces and maintains the integration accuracy as much as possible. Recently, the stochastic collocation method has also served as an efficient way for fatigue crack propagation and health prognosis [16,17].
A large number of stochastic methods for FCG analysis have been developed to investigate the uncertainty of remaining useful life prediction [18][19][20][21]. In general, the traditional fatigue prognosis methods are performed in a cycle-based manner. This class of methods requires one to transform a realistic loading history by cycle counting before the FCG analysis is performed. Paris proposed the most widely applied fracture mechanics-based model for fatigue crack propagation in the 1960s. The Paris model relates the stress intensity factor range ΔK with the crack growth rate [22]. However, the original Paris model can only be applied to calculate the crack growth rate in a certain range with fixed stress ratios under amplitude loading, which is impractical. Many modifications and extensions of Paris model have been proposed. Forman et al. incorporated the stress ratio effect in the Paris model [23]. In addition, the inclusions of the fracture toughness and the threshold stress intensity range based on the Forman formulation were proposed. Wolf observed the contact between the crack surfaces and introduced the concept of crack closure into crack propagation analysis under cyclic loading for the first time [24]. Wolf proposed to replace ΔK with the stress intensity factor ΔKEFF, which represents the difference between the intensity factor at the crack opening load and maximum stress intensity factor. Many experimental and theoretical studies have been made based on Wolf's research [25,26]. Nevertheless, it has been argued that the concept of crack closure is unable to be applied under some special conditions [27,28]. Another group of approaches including the Willenborg model and Wheeler model investigate a much more complicated problem-the FCG analysis under variable amplitude loading. Both models employ the interaction of the plastic zone to explain the complex FCG behavior.
However, the abovementioned models are all cycle-based, which has many inherent difficulties even for some widely-used models. Due to the temporal scale limitation, a realistic stress history needs to be transformed into a cyclic manner time history, and it is impossible to reduce the time scale in more fundamental investigation. A new FCG formulation, small-timescale FCG model, was proposed to overcome the limitations of previous models [29]. The core issue in the small-timescale model is the calculation of the incremental crack growth at some arbitrary time instant during a loading cycle. Extensive experiments under both constant and variable amplitudes have been implemented to validate the accuracy of the novel small-timescale model for FCG analysis [30][31][32][33][34]. Previous investigations have indicated that the presented model greatly promotes the fatigue prognosis for structural systems. In addition, it has been observed from experimental observations that huge scatter exists in crack growth process [35].
In view of above discussions, the objective of this paper was to investigate the uncertainty quantification in small-timescale model-based FCG analysis using a stochastic collocation method. Firstly, the fundamental basis of the stochastic collocation method was presented. Secondly, the small-timescale model was introduced briefly and the framework of uncertainty quantification based on this model was proposed. Then the experimental data for 7075-T6 aluminum alloy from the literature was adopted to demonstrate the computational accuracy and efficiency of the proposed approach. Finally, Section 5 draws some conclusions.

Stochastic Collocation Method for Uncertainty Quantification
In this section, the uncertainty quantification in the small-timescale model, which is defined by partial differential equations (PDE) with uncertain input parameters, is considered. A class of generalized polynomial chaos (gPC) method known as the stochastic collocation method is a popular choice for these complex systems. The stochastic collocation method transforms the random problems to the corresponding deterministic problems at each collocation point by using the high dimensional polynomial interpolation technique [10]. In the following, this method is briefly introduced after the description of the formulation of stochastic systems.

Stochastic Systems and Interpolation Formulation
The stochastic PDE system is defined in a spatial domain ⊂ D   = 1,2, 3  and a time domain [0,T]: where  is a differential operator,  is the support of Z ,  is a boundary condition operator, 0 u is the initial condition, where =1 u n is the dimension of u for the crack growth system in this paper. For any given x and t, it is feasible to find a numerical approximation ( ) The multivariate gPC expansion is an efficient approximation for the stochastic process and random variable.
be the univariate gPC basis functions with degree up to N . The basis functions satisfy: where γ m is a normalized factor. i 1 =( ,..., ) The Nth degree gPC basis functions for d variates are given by products of gPC polynomials as: Subsequently, the gPC approximation of ( ) u Z is: where î v is the polynomial chaos expansion coefficient. The mean value and variance of ( ) u Z can be approximated by the formulation.
In the collocation method, I is a set of nodes generated from the random space.
For all = 1,..., j M , the stochastic PDE system is transformed at the node ( ) j Z by solving the deterministic problem: , the Lagrange interpolation polynomial can be expressed as: where: are the Lagrange interpolating polynomials. The Lagrange polynomials have a specific property: Several choices exist for the interpolation nodes and the interpolation for univariate problems is well studied. As for the multivariate problem, the Lagrange interpolation can be extended to apply to the entire multidimensional space. The interpolation is then constructed by the full tensor product or the sparse grid collocation introduced in the following sections.

Full Tensor Product Collocation
For the multivariate problem with > 1 d , a straightforward approach is using the tensor product of the Lagrange interpolation operator for each dimension: based on the one-dimensional nodal sets: Then in the entire space ⊂ d Z I  , the interpolation formulation calculated by tensor product is: in a general way is the level of interpolation on each dimension, the level of the interpolation polynomial, and = + q d k . According to the basic property of the tensor product, the detailed full tensor product formulation is: and the corresponding nodal set is: Usually, the same interpolation formulation is used in each dimension and the number of the points is the same in each dimension, i.e., Therefore, the total number of points, = d M m , increases rapidly in high dimensions, which makes the tensor product approach timeconsuming.

Smolyak Sparse Grid Collocation
An alternative to the tensor product, the Smolyak sparse grid collocation method, has wide applications due to its high accuracy and fast convergence. The Smolyak algorithm constructs a linear combination of tensor products to select the set of collocation points. The linear combination uses a relatively small number of points and preserves the interpolation property of = 1 d for > 1 d compared with the full tensor product. The Smolyak algorithm is based on the difference value of the one-dimension interpolation polynomials on each level of interpolation. Let us define: Then the formulation of the Smolyak algorithm is: where is the level of the interpolation polynomial, same as the tensor product. Thus, the number of values of the multi-index i 1 =( ,..., ) The construction of the Smolyak sparse grid can also be expressed as: The corresponding nodal set for the interpolation is a sparse grid: In the Smolyak algorithm, the Clenshaw-Curtis sparse grids are used to construct the onedimensional interpolation formulations. The Clenshaw-Curtis sparse grids are the extreme value of Chebyshev polynomials. For any ≤ ≤ 1 i d , the points are given by:  [36].

Description of Small-Timescale FCG Model
The detailed derivations of the small-timescale model can be found in [29] and a brief introduction of this model is presented here. The basis of this FCG model is the incremental crack growth at any arbitrary small incremental time in a loading cycle. Figure 1 shows the geometric relationship constructed between the instantaneous crack growth kinetics and the crack tip opening displacement (CTOD). Thus, the crack growth rate da is calculated by: where θ is the crack tip opening angle (CTOA), δ is the CTOD, and θ = / 2 C ctg . It is assumed that the crack growth in Equation (20) is infinitesimal. The CTOA changes from °90 to a very small angle beyond a certain length until the final failure. Thus, the CTOA is assumed approximately to be a function of material properties and the applied loading, which can be expressed as: where ΔK is the stress intensity factor range, Δ th K is the intrinsic threshold stress intensity factor, and c K is the fracture toughness value. Δ th K corresponds with the intrinsic material resistance to the FCG. In this study, the FCG problem is a predominantly plane stress problem for surface cracks. Considering the effect of material strain hardening, the CTOD can be expressed as: where π λ σ is a constant of material property, σ is the nominal stress, E is Young's modulus, and σ y is the yield stress. According to Figure 1, when the crack tip grows from O to ' O , the increment of CTOD ignoring the high order terms can be expressed as: Substituting Equation (23) into Equation (20), the instantaneous crack growth rate can be expressed as: The previous discussion can only be applied when the crack grows. Nevertheless, the crack tip plastic zone is critical for FCG and the reverse plastic zone produces a compressive residual stress after unloading. Due to the energy principle and the reverse plastic zone effect, the crack growth may stop during the initial part of the loading cycle as well as during the unloading cycle. Using the hypotheses for the non-uniform crack growth kinetics, the general expression of the instantaneous crack growth rate under arbitrary loading histories is: is the Heaviside function and σ ref is the reference applied stress level.
Crack will only grow beyond a certain stress level σ ref in the loading path. A detailed calculation of σ ref on the basis of the interaction of forward and reversed plastic zone can be found in [29].
Accordingly, the equivalent cycle-based formulation of the small-timescale model can be obtained by solving the integral at both ends of the integral of Equation (24). The length increment Δa during one cycle can be expressed as: where:  (25)). The objective of this study was to adopt the previously discussed small-timescale model for the uncertainty quantification in the FCG analysis. The small-timescale model is fundamentally different from the classical cycle-based FCG and it can easily be utilized on any timescale and over any crack-length range.

Uncertainty Quantification in FCG
The presented stochastic collocation method is addressed accounting for uncertainty quantification in the small-timescale FCG model. The stochastic collocation separates uncertainty quantification from solving differential equations. The framework of the approach includes four steps. First, define the uncertainty sources and the probability space in the study and choose the nodal set Θ d M from the probability space based on the extrema of Chebyshev polynomials and the Smolyak algorithm. Second, solve the deterministic differential equations (Equation (27)) of the small-timescale FCG model at the nodal set. Then, the Lagrange interpolation polynomial is employed to construct the interpolation polynomial for remaining useful life prediction based on the nodal set and the solution set. Finally, estimate the probability characteristics of the prediction results by MC simulation as the post process. The uncertainty quantification framework of the stochastic collocation method is shown in Figure 2.

The Remaining Useful Life Prediction
The remaining useful life prediction is performed at an inspection point when the crack length is estimated as initial crack length 0 a . The component will end in fatigue failure when the crack propagates from 0 a to the critical length c a in the small-timescale FCG model. According to the previous discussion, the crack growth rate under constant amplitude loading is employed here. Let us rewrite the formulation of crack growth rate: where N is the number of loading cycles when the crack length propagates from 0 a to the critical length c a . Then the crack growth process is established by solving the following equation using the Runge-Kutta algorithm: The remaining useful life N at the inspection point can be expressed as: The random input parameters in the fatigue crack propagation can be denoted by The distributions of the random parameters are assumed to be independent and identical. Then the remaining useful life N is denoted by a function of the random parameters Z as ( ) N Z .

Results and Discussion
In order to show the efficiency and accuracy of the proposed method, uncertainty quantification in the crack growth model of a metal plate is implemented with the stochastic collocation method using the experimental data available in previous literature [35]. First, the crack growth curves are established and the uncertainty parameters are defined in the FCG model. Then the stochastic collocation method is applied to predict the remaining useful life accounting for the uncertainty. The prediction results are compared with the experimental data from the literature.

FCG Analysis of A Metal Plate
The experimental data on 7075-T6 aluminum alloy plate specimens with a center through crack ( Figure 3) taken from the above research are collected to investigate the proposed approach of uncertainty quantification for FCG. For easier access to the information in the previous paper, some essential information is presented here: The material and geometry properties of the specimen and the corresponding test data are listed in Table 1 and Table 2, respectively. Figure 4 shows the seven raw data curves under constant amplitude loading. It can be observed that the scatter of raw data curves exists due to the material uncertainty. Based on the test data, the initial crack length is ac = 0.011 m and the critical length is set as ac = 0.0258 m to investigate the remaining useful life prediction.  In this study, it is assumed that the uncertainty comes from the model parameters uncertainty. To generate the simulated crack growth curves, the first step is to quantitate the statistics of the material crack growth parameters. The small-timescale model has four independent underlying parameters to be determined in Equation (27), i.e., the Young's modulus E , the intrinsic threshold stress intensity factor Δ th K , the fracture toughness value c K , and the yield strength σ y . According to the reference [35], the Young's modulus E is assumed to be a constant. The randomness in the crack growth is represented by three random variables Δ th K , c K , and σ y . Then according to the reference [6] and [37], the distribution characteristics of these random variables are determined by fitting crack growth data from the literature. It is found that the lognormal distribution is appropriate for the materials investigated in this paper (as shown in Table 3).
where α = a/W, and the calculation is valid when α ≥ 0.2 . B and W are the thickness and the width of the specimen, respectively. The generated crack growth curves by a directed Monte Carlo simulation (10000 samples) show a general satisfactory agreement compared with the seven raw test data curves ( Figure 5), which indicates the ability of the small-timescale model to describe the fatigue propagation. The negligible inconsistency may be caused by different types of uncertainty including natural variability in local geometry, variant loading conditions, and data uncertainty due to measurement errors during experiment, which are over the range of this study. In this investigation, due to the lack of the exact solution of the problem, the prediction results from the proposed method are compared with the results from Monte Carlo simulation, which is considered as the real distribution of the remaining useful life.

The Stochastic Collocation Method for Remaining Useful Life Prediction
Based on the small-timescale FCG model under constant amplitude loading (Equation (27)), the remaining useful life N of the 7075-T6 aluminum alloy plate can be expressed as Equation (31).
When calculating the remaining useful life in the uncertainty damage prognosis, the material parameters c K , Δ th K and σ y are sampled from the prescribed lognormal distributions. Therefore, the remaining useful life prediction formula for the specimen with uncertainty is: where σ = Δ ( , , ) c t h y Z K K is the set of independent random variables, and Z I is the probability space. The quantization interval of the random variable is defined as μ σ μ σ to ensure the probability of the value outside the interval is less than 0.003. Thus, the probability space is set as: In the stochastic collocation method, the three-dimensional standard probability space − 3 [ 1,1] is required by using the linear transformation of the previous probability space. The core issue for the implementation of the Smolyak sparse grid is the choice of the nodal set depending on the level of interpolation k . To investigate the computational accuracy, the relative cumulative error S err is defined as: and the norm is: where ( ) k N z is the k-level interpolation polynomial. The relative cumulative error decreases rapidly as the level of interpolation increases as shown in Figure 6a. On the other hand, the computational efficiency of the proposed approach is considered by the total number of nodes for the threedimensional interpolation problem, which is also the number of solving the differential equations M (Figure 6b). According to the above discussion, the level of interpolation =5 k is appropriate to balance the computational accuracy and efficiency with . Take for example, the elapsed time of directly solving the differential equations (9.1034s) is more than three times as the elapsed time using 5 ( ) N z (2.7600 s) at 10000 samples. In this study, 5 ( ) N z is selected as the interpolation polynomial of the remaining useful life taking advantage of the efficient sparse grid method.

Results of Remaining Useful Life Prediction
The mean value and the standard deviation of the remaining useful life N are convergent to × × . Due to the enormous number of samples, four types of distribution including uniform distribution, normal distribution, 2-parameter lognormal distribution and 3-parameter lognormal distribution, are used to fit the distribution of N by a Kolmogorov-Smirnov test (K-S test). The K-S test quantifies a distance D between the cumulative distribution function (CDF) of the samples and the CDF of the reference distribution. It is considered that the samples take the specific distribution when α < ( , ) D D n . The null hypothesis is rejected at levelα = 0.05 , so that the critical distance is defined as (10000,0.05) D . Table 4 lists the results of the K-S test indicating that N takes a 3-parameter lognormal distribution. The probability density function can be expressed as: where 0 N , μ and σ are constant parameters. The probability density function (PDF) and the CDF of the samples and the 3-parameter lognormal distribution are shown in Figure 7 and Figure 8, and the MC simulation results from solving the differential equations directly are also included for comparison. Figure 7 and Figure 8 show that the 3-parameter lognormal distribution fit the dataset of the remaining useful life N acceptably, resulting in the following distribution parameters: N0 = 24306, μ = 4.408 and σ = 0.0135 . In this study, due to the lack of the experimental data, the results from MC simulation by solving the differential equations in the small-timescale FCG model is considered as the real distribution of the remaining useful life.
It can be observed that the distribution of remaining useful life prediction using the proposed stochastic collocation method shows the goodness of fit of the MC simulation. Furthermore, the predicted remaining useful life with 95% reliability is = 24610 N , which is among the range of experimental results. In addition, it is noted that the total number of nodes, which is also the number of solving the differential equations is 757 and the number of realizations in direct MC simulation is 10000. The present results indicate the proposed approach has the capability to give an overall efficient and accurate prediction of the FCG taking into account uncertainties of input parameters.

Conclusions
In this paper, a stochastic collocation approach is developed for uncertainty quantification of remaining useful life prediction in small-timescale FCG model. The presented FCG model improves the accuracy of deterministic prediction of − a N curve based on the instantaneous crack growth increment. Then the stochastic collocation method is utilized to evaluate the remaining useful life prediction uncertainty. The stochastic collocation method based on the gPC expansion is able to reduce the repetitious equation solving significantly. The main conclusions of this research are: 1) The FCG curves predicted by the Monte Carlo simulation of the small-timescale model considering the uncertainty parameters agreed well with the experimental data under constant amplitude loading. It can be observed that parameter uncertainties do exist in FCG analysis.
Besides, the small-timescale model shows a continuous relationship between the fatigue crack length and the loading history.
2) The remaining useful life predictions obtained from the proposed approach were compared with the Monte Carlo simulation results for the plate specimen. The comparisons show that the proposed approach can greatly improve the efficiency and accuracy for uncertainty quantification in FCG analysis by utilizing the stochastic collocation method.
3) The significant computation efficiency of this uncertainty quantification approach for FCG can be potentially applied to much more complicated multi-dimensional problems for crack and damage propagations. Moreover, this study focuses on the material uncertainty quantification and other types of uncertainties can be involved in potential applications.