Statistical inference of equivalent initial flaw size with complicated structural geometry and multi-axial variable amplitude loading

International Journal of Fatigue 32 (2010) 1689–1700

Contents lists available at ScienceDirect

International Journal of Fatigue journal homepage: www.elsevier.com/locate/ijfatigue

Statistical inference of equivalent initial flaw size with complicated structural geometry and multi-axial variable amplitude loading Shankar Sankararaman, You Ling, Sankaran Mahadevan * Department of Civil and Environmental Engineering, Vanderbilt University, Nashville, Tennessee 37235, USA

a r t i c l e

i n f o

Article history: Received 17 August 2009 Received in revised form 1 February 2010 Accepted 25 March 2010 Available online 3 April 2010 Keywords: Fatigue life prediction Initial flaw size Crack growth Multi-axial loading Bayesian inference

a b s t r a c t This paper presents several efficient statistical inference techniques to calibrate the equivalent initial flaw size (EIFS) of fatigue cracks for mechanical components with complicated geometry and multi-axial, variable amplitude loading. Finite element analysis is used to address the complicated geometry and calculate the stress intensity factors. Multi-modal stress intensity factors due to multi-axial loading are combined to calculate an equivalent stress intensity factor using a characteristic plane approach. During cycle-by-cycle integration of the crack growth law, a Gaussian process surrogate model is used to replace the expensive finite element analysis, resulting in rapid computation. Experimental data (crack size after a particular number of loading cycles) and statistical methods are used to calibrate the EIFS. The methods of least squares and maximum likelihood method are extended to evaluate the entire probability distribution of EIFS. Bayesian techniques are also implemented for this purpose. A fast numerical integration technique is developed as an efficient alternative to the expensive Markov Chain Monte Carlo sampling approach in the Bayesian analysis. An application problem of cracking in a cylindrical structure is used to illustrate the proposed methods. Ó 2010 Elsevier Ltd. All rights reserved.

1. Motivation Mechanical components of engineering systems are often subjected to cyclic loads, leading to fatigue, crack initiation and crack growth. System models are coupled with crack growth models to study crack growth propagation. This, in turn, can be used in life prediction and reliability evaluation of the system. A major problem in fracture mechanics-based life prediction is that the initial crack size in the structural component is unknown. This issue is further complicated by the fact that small crack growth propagation is anomalous in nature [1]. This problem was overcome by the introduction of an equivalent initial flaw size (EIFS) nearly thirty years ago [2]. The concept of EIFS was introduced to by-pass small crack growth analysis and to use EIFS as the initial crack size in long crack growth models such as Paris’ law. However, EIFS does not represent any physical quantity and hence is impossible to measure through experiments. Initially, empirical crack lengths in the range of 0.25–1 mm for metals were assumed [3–5]. Since then, several rigorous procedures have been developed to estimate the EIFS. Several researchers have estimated EIFS using back extrapolation of experimental data [6,7]. The underlying idea of this approach is to use a crack growth model and experimental data * Corresponding author. Tel.: +1 615 322 3040; fax: +1 615 322 3365. E-mail address: [email protected] (S. Mahadevan). 0142-1123/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijfatigue.2010.03.012

(number of load cycles required to reach a final crack size) to estimate the initial flaw size. Sometimes, this requires inversion of crack growth models, which may be expensive. A disadvantage of the back-extrapolation method is that the results are dependent on the stress level [2]. Hence, the EIFS estimate with this method becomes problem specific and not suitable for extrapolation to other structural configurations and loading, thus requiring a new experimental data set for each loading case. The concept of EIFS was modified by Barter [8] using an equivalent pre-crack size (EPS). This method involves little extrapolation and the solution is independent of stress levels and load amplitude. This concept of EPS has been pursued by other researchers as well. White et al. [9] proposed a probabilistic fracture approach to derive the equivalent pre-crack size based on the back-extrapolation method. Molent et al. [10] derived the EPS for Al 7050 using the back projection of the experimental crack growth curve. These methods have been illustrated for simple structures under constant amplitude uniaxial loading. Liu and Mahadevan [2] developed a new technique to calculate EIFS based on the Kitagawa–Takahashi [11] diagram and the ElHaddad [12] model. This method is independent of the load levels and uses the concepts of fatigue limit and threshold stress intensity factor. This line of research work has so far considered only smooth plate specimens under uniaxial [13] and multi-axial [14] loading conditions. The methodology requires the calculation of a geometry factor (Y) to estimate the equivalent initial flaw size. This

1690

S. Sankararaman et al. / International Journal of Fatigue 32 (2010) 1689–1700

geometry factor is a function of the geometry of the structure as well as the crack configuration. The calculation of Y is not straightforward for structures with complicated crack configurations. The estimation of EIFS has been viewed as a model calibration problem by Makeev et al. [15] and Cross et al. [16] and has been solved using the methods of maximum likelihood and Bayesian inference respectively. In this approach, EIFS is treated as a model parameter that is calibrated through observation of crack sizes at different numbers of cycles. These methods have been previously illustrated for simple structural elements subjected to constant amplitude loading. This paper extends the work of Makeev et al. [15] and Cross et al. [16] to structures with complicated geometry and multi-axial variable amplitude loading. It develops an efficient and accurate approach (within the general framework of [15,16]) to combine several types of analysis – finite element analysis, crack growth analysis, surrogate model construction, Bayesian inference, and probability integration – in implementing EIFS calibration methodology to realistic structures with complicated geometry and multi-axial variable amplitude loading. The calculation of stress intensity factor is computationally difficult in the presence of complicated geometry. Look-up tables are available for the calculation of stress intensity factor for a few wellknown standard structures and certain known crack configurations. Examples include cylinders, plate with a central hole, etc. with standard loading conditions such as uniform loading, uniaxial loading, etc. For structures with complicated geometry, and unseen crack configurations, and under multi-axial variable amplitude loading, finite element analysis may be necessary to calculate the three stress intensity factors (mode-I, mode-II, and mode-III) to account for mixed mode crack growth. However, it is computationally expensive to incorporate a finite element model into the above methods, since they would require repeated use of finite element analysis, in every loading cycle; therefore, an accurate surrogate modeling approach is used in this paper. Experimental observations in the form of crack size of specimens after a particular number of load cycles are used to infer the EIFS. A suitable long crack growth model (modified Paris’ law) is chosen and is used in conjunction with Wheeler’s retardation model to account for variable amplitude loading, for the sake of illustration. Multi-modal stress intensity factors due to multiaxial loading are combined to calculate an equivalent stress intensity factor using a characteristic plane approach. Starting with the EIFS, the growth of the crack in each load cycle is calculated. In each cycle, the stress intensity factor needs to be computed as a function of the current crack size. Thus, the number of load cycles to reach a particular crack size can be calculated using the crack growth model. The equivalent initial flaw size can be computed by comparing model predictions and corresponding experimental observations. Note that this procedure is significantly different from the methods used for EIFS calculation using inspection data, which are not based on statistics and use back extrapolation of the crack growth curve to calculate the EIFS. It is essential to observe that fatigue cracking is a stochastic process, i.e. there are several kinds of uncertainties involved in crack growth. For example, the loading conditions are random, experimental observations are noisy and the crack growth model is subjected to uncertainties and model errors. Hence methods for estimation of EIFS should account for these sources of uncertainty and calculate a probability distribution for EIFS instead of just point estimates. In this paper, three different statistical techniques – least squares approach, maximum likelihood, Bayesian inference – are investigated to estimate the probability distribution of EIFS. Conventionally, the least squares framework and the method of maximum likelihood have been used to obtain point estimates [17] of the inferred parameter (EIFS, in this case). However, this paper extends these methods to obtain the entire probability distri-

bution of EIFS. If the distribution parameters of EIFS are uncertain, then they can also be inferred using the proposed methods. The proposed approach also solves the difficulties of expensive computations in two different steps of the EIFS calibration process. First, the calculation of stress intensity factor used in crack growth analysis involves the remeshing of the finite element model and use of the expensive finite element analysis in each loading cycle. Second, statistical inference using Bayesian updating is commonly done using Markov Chain Monte Carlo (MCMC) sampling which requires numerous model evaluations (of the order of 106 [16]). This paper proposes methods to reduce the computational effort in both types of calculations. The first issue is addressed through the use of an accurate Gaussian process surrogate model to replace the finite element analysis. A few finite element analysis runs are used to train the surrogate model, which can then predict the stress intensity factor to be used in crack growth analysis in each cycle. The second issue, i.e. the integration involved in Bayesian updating is solved numerically through the use of an advanced quadrature technique instead of MCMC sampling. Thus, the proposed methods drastically reduce the computational costs in the estimation of EIFS and facilitate application to practical engineering problems. Section 2 explains the various steps in the proposed methodology for the inference of EIFS. First, it introduces the concept of EIFS in more detail and explains its use in crack growth analysis. The proposed method can be applied with any long crack growth law and this paper uses a modified Paris’ law (with Wheeler’s retardation model) for illustration purposes only. Section 3 introduces the concept of a Gaussian process model as a surrogate for finite element analysis, in order to calculate the stress intensity factor for use in crack growth analysis in each cycle. Section 4 develops several different inference techniques to quantify the statistics of EIFS. Section 5 develops an efficient integration method to replace the MCMC simulation. Section 6 illustrates the proposed techniques using a numerical example, and Section 7 provides the concluding remarks. 2. Proposed framework for the inference of EIFS This section describes the various steps involved in the proposed methodology for the estimation of statistics of EIFS. The basic concept of EIFS is introduced and the use of EIFS in fatigue crack growth analysis is discussed. Although the following discussion uses Paris’ law, any crack growth model can be used to implement the proposed methodology. This crack growth model is used to simulate the growth of the crack and hence predict the number of cycles to reach a given crack size. At each cycle, the stress intensity factor is to be calculated using a Gaussian process surrogate model, which is in turn trained using a few runs of finite element analysis. A crack propagation framework is established using which it is possible to estimate the number of loading cycles used to reach a particular final crack size. The model predictions are compared with the experimental data, and the statistics of EIFS are calculated using a model calibration approach. The various steps involved are shown schematically in Fig. 1. The following sub-sections explain in detail each of the steps involved in the inference of EIFS. 2.1. Basic concept of EIFS A rigorous approach to fatigue life prediction would be to perform crack growth analysis starting from the actual initial flaw, accounting for imperfections, voids and non-metallic inclusions. If the initial crack size is large, then long crack growth models such as Paris’ law can be used directly. However, this is not the case in

1691

S. Sankararaman et al. / International Journal of Fatigue 32 (2010) 1689–1700

Multi-axial, Variable Amplitude Loading

Surrogate Model Trained Using FE Model (Section 3)

EIFS (Section 2.1)

Measurement Data

Crack Propagation Analysis (Section 2.3)

Crack Size after N cycles

Predict Number of Cycles to Reach Given Final Crack Size (Section 2.3)

Infer Statistics of EIFS (Section 4)

Fig. 1. Proposed model calibration framework for the inference of EIFS.

1/g(a)

Equal area 1/g l (a) 1/g s(a)

IFS

EIFS

ac

a

Fig. 3. Actual IFS vs. equivalent IFS. Fig. 2. Small crack growth vs. long crack growth.

most materials. Hence the long crack growth model cannot be used directly. A schematic plot of the long crack and short crack growth curves is given in Fig. 2. Consider the crack growth model given in Eq. (1). N represents the number of cycles, a represents the crack size and DK represents the stress intensity factor. f may be any function used to describe the relationship between da/dN and DK, such as the Paris law, NASGRO equation, etc.

da ¼ f ðDKÞ dN

ð1Þ

The stress intensity factor can be expressed as a function of the crack size and loading (L).

da ¼ gða; LÞ dN

ð2Þ

As mentioned earlier, the rigorous approach to calculate the fatigue life would be to use a crack growth law (gs), which accounts for both the short crack region and the long crack region. This implementation is shown in Eq. (3), where ‘IFS’ stands for the actual initial flaw size, and A stands for the crack size after N cycles.

N¼

Z

dN ¼

Z

A IFS

1 da g s ða; LÞ

N¼

dN ¼

Z

A h

1 da g l ða; LÞ

An illustration of Eq. (3) and (4) is shown in Fig. 3.

2.2. Crack growth analysis and EIFS Inference First, a crack growth analysis procedure for complicated structures subjected to multi-axial random amplitude loading is discussed. Second, EIFS is treated to be a model parameter in this methodology and is inferred using experimental observations. The method for crack growth analysis can use any long crack growth model and a modified Paris law is used in this paper for illustration purposes only. Further, Wheeler’s retardation model [18,19] for fatigue crack growth is also included to account for the effects of variable amplitude loading. The expression in Eq. (1) is replaced by Eq. (5).

  da DK th m ¼ ur CðDKÞn 1  dN DK

ð4Þ

ð5Þ

In Eq. (5), ur refers to the retardation parameter [19], which can be defined as:

ð3Þ

An alternative approach would be to use an equivalent initial flaw size (EIFS, denoted by h) and only a long crack growth model (gl) to calculate the fatigue life, as shown in Eq. (4).

Z

The EIFS can be estimated by equating areas under the curves 1/ gs(a) and 1/gl(a) shown in Fig. 3. The choice of the crack growth law will affect the value of EIFS. As mentioned earlier, EIFS is not a physical quantity; it is a model calibration parameter estimated from experimental observations.

r

u ¼

8h < :

ik r p;i aOL þr p;OL ai

1

if ai þ r p;i < aOL þ r p;OL

ð6Þ

if ai þ r p;i P aOL þ r p;OL

In Eq. (6), aOL is the crack length at which the overload is applied, ai is the current crack length, rp,OL is the size of the plastic zone produced by the overload at aOL, rp,i is the size of the plastic zone produced at the current crack length ai, and k is the curve fitting parameter for the original Wheeler model termed the shaping

1692

S. Sankararaman et al. / International Journal of Fatigue 32 (2010) 1689–1700

exponent [18]. Sheu et al. [19] and Song. et al. [20] observed that crack growth retardation actually takes place within an effective plastic zone. Hence the size of the plastic zone can be calculated in terms of the applied stress intensity factor (K) and yield strength (r) as:

 2 K rp ¼ a

ð7Þ

r

In Eq. (7), a is known as the effective plastic zone size constant which is calculated experimentally [18]. The expressions in Eq. (6) and (7) can be combined with Eq. (5) and used in crack growth analysis. Then the expression in Eq. (5) can be integrated to calculate the number of cycles (N) required to reach a final crack size (A), as:

N¼

Z

dN ¼

Z

A h

u

r CðDKÞn

1  m da 1  DDKKth

ð8Þ

The stress intensity factor DK in Eq. (6) depends on the loading conditions, geometry and crack size; it can be expressed as a closed form function of the crack size for specimens with simple geometry subjected to constant amplitude loading. However, in many mechanical components, the loading is multi-axial (for example, simultaneous tension, torsion and bending), and the stress intensity factors corresponding to multiple cracking modes need to be taken into account. This can be accomplished using an equivalent stress intensity factor. For example, if KI, KII, KIII represent the mode-I, mode-II and mode-III stress intensity factors respectively, then the equivalent stress intensity factor Keqv has been calculated using a characteristic plane approach proposed by Liu and Mahadevan [21]. The use of the characteristic plane approach for crack growth prediction under multi-axial, variable amplitude loading has been validated earlier with several data sets [22].During each cycle of loading, if the crack grows, then, the stress intensity factor needs to be reevaluated at the new crack size for the loading in the next cycle. Hence, it becomes necessary to integrate the expression in Eq. (8) through a cycle-by-cycle procedure. Each cycle involves the computation of DK using a finite element analysis represented by W.

DK eqv ¼ Wða; LÞ

ð9Þ

and the principle of maximum of likelihood usually yield point estimates, they are extended here to calculate the entire distribution of the EIFS (h). These approaches are explained in detail in Section 4.

3. Gaussian process modeling As mentioned in the previous section, the expression in Eq. (9) used to calculate the stress intensity factor may be evaluated through finite element analyses for structures with complicated geometry. Repeated evaluation of the stress intensity factor at different crack sizes is impractical with an expensive, time-consuming finite element model. Hence, the finite element model needs to be replaced by a surrogate model trained using input (crack size, loading conditions) and output (stress intensity factor) variables. Different types of surrogate modeling techniques (conventional polynomial response surface, polynomial chaos expansion [23], support vector regression [24], relevance vector regression [25], and Gaussian Process interpolation [26]) were investigated, and the Gaussian process modeling technique was found to be most effective. Gaussian process (GP) modeling is a powerful technique based on spatial statistics for interpolating data and is increasingly being used to build surrogates to expensive computer simulations for the purposes of optimization and uncertainty propagation [26– 28]. The procedure used to construct a GP model is summarized in Fig. 4. The basic idea of the GP model is that the response values Y (Keqv in this case) evaluated at different values of the input variables X (crack size and loading conditions in this case), are modeled as a Gaussian random field, with a defined mean and covariance function. There are three important reasons why a GP model has been used in this research work. i. The GP model is capable of capturing highly nonlinear relationships that exist between input and output variables without the need for an explicit functional form. Hence, a closed form expression (as in polynomial type regression methods) need not be assumed. ii. For a non-parametric interpolation technique, this method requires fewer sample points (usually 30 or less) as against methods such as kernel estimation and non-parametric multiplicative regression. iii. A GP model provides a direct estimate of the variance in the output prediction.

Repeated mesh generation and execution of the finite element analysis in Eq. (9) makes the aforementioned cycle-by-cycle integration extremely expensive, perhaps impossible in some cases. Hence, it is necessary to substitute the finite element model by an inexpensive surrogate model. Different kinds of surrogate models have been explored and the Gaussian process modeling technique has been employed in this paper. A few runs of the finite element analysis are used to train this surrogate model and then, this model is used to predict the stress intensity factor for other crack sizes and loading magnitudes during the cycle-by-cycle crack growth analysis. The construction of a Gaussian process model is explained in Section 3. The integral in Eq. (8) can be viewed as a black box G as shown in Eq. (10), where N represents the number of cycles to reach a final crack size A from a given equivalent initial flaw size h.

Suppose that there are m training points, x1, x2, x3 . . . xm of a ddimensional input variable vector (the input variables being the crack size and loading conditions in this case), yielding the resultant observed random vector Y(x1), Y(x2), Y(x3) . . . Y(xm) (the stress intensity factor, in this case). R is the m  m matrix of correlations among input variables at the training points. Under the assumption that the parameters governing both the trend function (fT(xi), at each training point) and the covariance (k) are known, the expected value and the variance of the Gaussian process at an untested location x is calculated as in Eq. (11) and (12) respectively.

N ¼ Gðh; AÞ

Y  ¼ E½Yðx ÞjY ¼ f T ðx Þb þ r T ðx ÞR1 ðY  FbÞ

ð11Þ

r2Y  ¼ Var½Yðx ÞjY ¼ kð1  rT R1 rÞ

ð12Þ

ð10Þ

Experimental data is collected in the form of ‘Crack size Ai after Ni cycles’, and used to infer the EIFS, i.e. h. Thus, this problem can be viewed as a model calibration problem in which the model parameter h needs to be evaluated. In this paper, well known statistical techniques such as the method of maximum likelihood, principle of least squares, and Bayesian inference are investigated for this purpose. While methods like the least squares technique

In Eq. (11) and (12), F is a matrix containing the rows of trend functions fT(xi) (i = 1 to m), r is the vector of correlations between x and each of the training points, b represents the coefficients of the regression trend. See McFarland [28] for details of implementation of this method.

1693

S. Sankararaman et al. / International Journal of Fatigue 32 (2010) 1689–1700

More Training Points

Finite Element Analysis (generate training points) No

Stress Intensity Factors KI, KII, KIII

Input Crack Size Tension Load Torsion Load Bending Load

Characteristic Plane

Accurate?

(Training Points) Keqv Gaussian Process (GP) (Predict Keqv)

Yes

Use GP in Crack Growth Model (Predict Keqv) Fig. 4. Flowchart for Gaussian process surrogate model construction.

A recently developed adaptive technique can be used to select training points for the GP model, in order to construct the response surface to a desired level of accuracy [28]. In order to accurately capture the functional relationship between input and output quantities, it was stated earlier that the Gaussian Process model must be ‘‘trained” using a set of observed inputs and outputs. In order to determine which points within the range of possible values (design space) would be most advantageous to use as training points, a greedy point selection technique has been proposed by McFarland [28]. Since the GP model is capable of estimating the variance in model output, the greedy point algorithm identifies the next training point at the input variable value which corresponds to the largest variance. By repeating this algorithm, training points are adaptively identified until the estimated variance is below a desired threshold. This adaptive method eliminates the subjectivity associated with choosing training points and ensures that the identified training points minimize the variance in the surrogate model output. See McFarland [28] for details of this procedure. Once the surrogate model is constructed using this approach, it can be used instead of the finite element model to predict the stress intensity factor, as shown earlier in Fig. 4. 4. Statistical Inference techniques for EIFS This section explores three techniques – a least squares based approach, a likelihood based approach and a Bayesian approach – for statistical inference of the equivalent initial flaw size. As mentioned earlier in Section 1, Makeev [15] and Cross [16] have used concepts of likelihood and Bayesian inference respectively to infer the probability distribution of equivalent initial flaw size. However, this earlier work addresses the inference of EIFS for structures with simple geometry and uniaxial, constant amplitude loading conditions. Two limitations of the earlier studies are addressed in this paper. i. A closed form expression for stress intensity factor is not available for complicated structures with multi-axial loading conditions. This problem is overcome by the use of finite element analysis, which is used to train an inexpensive surrogate model.

ii. It is not practical to use sampling-based techniques such as Monte Carlo simulation and Markov Chain Monte Carlo Simulation, which require repeated evaluation (of the order of 106) of the model G in Eq. (10) (each evaluation in turn requiring a cycle-by-cycle integration of the crack growth law and hence, a stress intensity factor calculation in each cycle). In lieu of sampling, an alternative fast numerical integration technique is proposed in this paper. The problem of estimation of EIFS is reduced to a model calibration problem as explained in Section 2.3. To recapitulate, this method evaluates the number of cycles needed to reach a final crack size (A) from a given equivalent initial flaw size (h) using the model G. A modified Paris’ law (Eq. (5)) along with Wheeler’s retardation model is used in crack growth analysis, and the stress intensity factor in each crack growth cycle is calculated as a function of crack size and loading through a Gaussian process model, as explained in Section 2. As explained earlier in Section 1, crack growth is a stochastic process and several sources of uncertainty – physical variability, data uncertainty and modeling errors – contribute to the uncertainty in crack growth. The loading conditions, physical parameters (such as yield stress, effective plastic zone size constant, etc.) are variable in nature. Data uncertainty refers to the sparseness of data and measurement inaccuracies. Modeling errors include uncertainty in the parameters of the modified Paris model, Wheeler model, model-form errors, etc. Hence, the model prediction is different from experimental observations and the aforementioned sources of uncertainty contribute to this difference, treated as a random variable (e) in this paper.

N ¼ Gðh; AÞ þ e

ð13Þ

Makeev et al. [15] and Cross et al. [16] defined a lognormal noise in the target data through the probability of observing a target value for a given input. The term e in Eq. (13) is the same as the aforementioned lognormal noise, the difference being that e is now assumed to be a normal random variable. The normal distribution assumption regarding the noise term in Eq. (13) is equivalent to lognormal distribution assumption of Makeev et al. [15] and Cross et al.[16], since the noise term is additive in Eq. (13), whereas it is multiplicative in [15,16]. Several classical statistics-based regression techniques [17] assume normally distributed noise, thereby

1694

S. Sankararaman et al. / International Journal of Fatigue 32 (2010) 1689–1700

allowing the use of both classical statistics-based methods and Bayesian inference for EIFS calibration. The various sources of uncertainty leading to this noise were mentioned in the paragraph preceding Eq. (13). It is assumed that experimental data is available in the form of (N, A) where N represents the number of cycles to reach a final crack size (A). The distribution of EIFS needs to be inferred using these data. The various methods used to calculate EIFS in this research work are listed below, under two types of estimation. i. Inference of EIFS. ii. The probability distribution of EIFS (h) can be calculated can be calculated using one of three different methods. a. Extended likelihood method. b. Bayesian method. c. Extended least squares method. iii. Inference of distribution parameters of EIFS. If the distribution parameters of EIFS are uncertain, then the probability distributions of these distribution parameters can be calculated using one of two different methods. a. Extended likelihood method. b. Bayesian method. These methods are represented through a schematic in Fig. 5. The following sub-sections discuss each of these inference methods in detail. 4.1. Inference on EIFS This subsection focuses on obtaining the distribution of equivalent initial flaw size (h) from data collected through experiments. Three different approaches are presented here to tackle this problem. The first approach is the likelihood method explained in Section 4.1.1. This approach is ideal to use when the analyst is completely ignorant about the distribution of h. The second approach is based on the principle of Bayesian updating. In some cases, the analyst might possess some information about the distribution of h. This may be used to construct a prior probability distribution and then updated using Bayes’ theorem, in presence of experimental evidence. This approach is presented in Section 4.1.2. The third approach is based on the well known least squares technique. In the literature, the least squares technique has been commonly employed to give only point estimates for the variable

for inference (h, in this case). However, this method is extended here to determine the entire distribution of h by calculating confidence intervals at different significance levels, as shown in Section 4.1.3. 4.1.1. Extended likelihood method The likelihood of a random variable is defined as the probability of observing an experimental result. Suppose that an experimental measurement suggests that it takes Ni number of cycles to reach a crack size Ai. Given h (EIFS) and Ai (final crack size), Eq. (11) can be used to evaluate the probability density of N = Ni. Suppose that there were m such independent experimental measurements, then the likelihood of h, i.e. L(h) can be calculated as shown in Eq. (14).

LðhÞ /

m Y

12



e

Maximize Likelihood L(θ )

rexp

ð14Þ

i¼1

This likelihood evaluated in Eq.14 is proportional to the probability density of h. The most likely value of h can be calculated by maximizing the expression in Eq. (14). This approach is known as the method of maximum likelihood. The likelihood approach can be extended to calculate the probability distribution of h by normalizing the expression in Eq. (14).

LðhÞ fh ðhÞ ¼ R LðhÞ

ð15Þ

The integral in Eq. (15) is a one-dimensional integral and does not need the use of sophisticated integration techniques such as Monte Carlo Integration or MCMC sampling. Therefore, this paper uses an advanced integration technique known as Adaptive Recursive Simpson’s Quadrature, which is described in Section 5. 4.1.2. Bayesian Method In this approach, a prior distribution is assumed for the uncertain quantity (h) and after collecting data through experiments, the distribution of h can be updated to estimate the posterior distribution h. If the observed data is represented by D, Bayesian updating of h can be expressed as:

f ðhjDÞ ¼ R

f ðDjhÞf ðhÞ f ðDjhÞf ðhÞdh

ð16Þ

The term f(D|h) represents the probability of observing data D, given h. (This is identical to L(h) in Eq. (14).) The term f(h) is the prior probability density of h. The posterior probability density f(h|D) is proportional to the product of prior times the likelihood.

EIFS (θ ) ~ LN (λ, ζ) Inference on λ, ζ

Inference on EIFS (θ ) Minimize Sum of Squares S(θ )

2

N i Gðh;Ai Þ

Black Box N = G (θ , A)+ε

Experimental Data Ni, Ai (i = 1 to m)

Bayesian Inference L(θ )*Prior Fig. 5. Schematic of inference methods.

Likelihood L(λ, ζ)

Bayesian Inference L(λ , ζ)*Prior

S. Sankararaman et al. / International Journal of Fatigue 32 (2010) 1689–1700

The denominator is a normalizing constant; this is also a onedimensional integral and can be easily evaluated using the integration technique described in Section 5. [Note: The extended likelihood method in Section 4.1.1 can be visualized as a limiting case of the Bayesian method. Consider a prior f(h) defined to be a uniform probability distribution on the interval [h/2, h/2]. Applying the limit h ? 1 to the expression in Eq. (16), f ðhjDÞ ¼ limh!1 R f ðDjhÞf ðhÞ ¼ limh!1 R LðhÞh ¼ RLðhÞ , f ðDjhÞf ðhÞdh

LðhÞhdh

LðhÞ

which is identical to the expression in Eq. (15) for the extended likelihood method.]

ments, the distributions of k and f can be inferred statistically. This can be done by two different methods, as explained below. 4.2.1. Extended likelihood method Consider m experimental measurements, (Ni, Ai; i = 1–m) Given k and f (parameters of the lognormal distribution of EIFS) and Ai (final crack size), Eq. (13) can be used to evaluate the probability density of N = Ni. The likelihood of (k, f), denoted by L(k, f) can be expressed as:

Lðk; fÞ / 4.1.3. Extended Least Squares Technique This well known technique has been commonly used to estimate unknown parameters in regression analysis. The model in Eq. (13) can be treated as a regression model and the least squares estimate of h (EIFS) can be evaluated. If there are m measurements, then a ‘total sum of squares error’, S(h) can be associated with any value of h as: m X SðhÞ ¼ ðNi  Gðh; ai ÞÞ2

ð17Þ

i¼1

The least squares estimate h is obtained by minimizing S(h) with respect to h. However this method can be extended to obtain the entire distribution of h, as explained below. A well known first step [29] is to estimate the confidence bounds on h at a particular significance level a. For a given significance level a, a target error value Starget can be defined as:

Starget

 ¼ Sðh Þ 1 þ 

p Fa n  p p;mp

 ð18Þ

In Eq. (18), F refers to the F-statistic evaluated at significance level a [29]. p refers to the number of parameters being inferred (p = 1 here), and m is the number of measurements. The upper and lower bound values of h satisfying Eq. (18), are estimated through an optimization technique. This is illustrated in Fig. 6. By choosing different values of a, the entire cumulative probability distribution can be obtained. For example, a = 0.1 yields 90% confidence bounds that correspond to cumulative distribution function (CDF) values of 0.05 and 0.95. 4.2. Inference on distribution parameters of EIFS If the distribution parameters of EIFS are uncertain, then these parameters can be calculated using the likelihood and Bayesian methods. Suppose expert opinion suggests that EIFS follows a lognormal distribution, but there is uncertainty about the values of the distribution parameters, k and f. Hence the distribution of EIFS can be denoted as fh(h|k, n). After collecting data through measure-

1695

m Z Y

e

12



Ni Gðh;Ai Þ

rexp

2 fh ðhjk; fÞdh

ð19Þ

i¼1

Using principles of conditional probability, the individual likelihoods of k and f can be expressed as:

LðkÞ / LðfÞ /

Z Z

Lðk; fÞdf

ð20Þ

Lðk; fÞdk

ð21Þ

Similar to that in Section 4.1.1, the most likely values of k and f can be estimated by maximizing the likelihoods in Eqs. (20) and (21). The probability distributions of k and f can also be obtained by normalizing the expressions as shown in Eq. (20) and (21).

R Lðk; fÞdf LðkÞ ¼ R R LðkÞdk Lðk; fÞdf dk R Lðk; fÞdk LðfÞ ¼ R R ff ðfÞ ¼ R LðfÞdf Lðk; fÞdk df

fk ðkÞ ¼ R

ð22Þ ð23Þ

The expressions in Eqs. (22) and (23) consist of one-dimensional integrals in the numerators and nested one-dimensional integrals in the denominator. Each of these can be evaluated through the adaptive quadrature technique described in Section 5. 4.2.2. Bayesian inference This method is similar to that in Section 4.1.2, the difference being that parameters being inferred here are k and f. It was already explained that the posterior distribution (obtained after Bayesian updating) is proportional to the product of prior times the likelihood. The joint likelihood of the two parameters k and f was expressed as L(k, f) in Eq. (19). Thus, the joint posterior distribution of k and f can be expressed as:

f ðk; fjDÞ / Lðk; fÞf ðk; fÞ

ð24Þ

Normalizing the expression in Eq. (24) would lead to a two-dimensional integral. There is an alternate route to evaluate the posterior distribution of k and f individually. First, the individual likelihoods are calculated. These are different from the expressions in Eqs. (20) and (21) to account for the prior information.

LðkÞ / LðfÞ /

Z Z

Lðk; fÞf ðfÞdf

ð25Þ

Lðk; fÞf ðfÞdk

ð26Þ

The posterior distributions of k and f can then be estimated as:

R Lðk; fÞf ðfÞdf f ðkÞ LðkÞf ðkÞ  R R ¼ LðkÞf ðkÞdk Lðk; fÞf ðfÞdf f ðkÞdk R Lðk; fÞf ðkÞdk f ðfÞ LðfÞf ðfÞ ¼ R R ffjD ðfÞ ¼ R LðfÞf ðfÞdf Lðk; fÞf ðkÞdk f ðfÞdf

fkjD ðkÞ ¼ R

Fig. 6. Confidence bounds at a-level.

ð27Þ ð28Þ

The expressions in Eqs. (27) and (28) consist of one-dimensional integrals in the numerators and nested one-dimensional integrals in the denominators. Thus, the two-dimensional integral required

1696

S. Sankararaman et al. / International Journal of Fatigue 32 (2010) 1689–1700

to normalize Eq. (24) has been simplified into one-dimensional integrals which are evaluated through an adaptive quadrature technique described in Section 5. 5. Method of integration This section develops an advanced integration technique to replace the expensive Markov Chain Monte Carlo (MCMC) sampling commonly used to evaluate the integrals in Eq. (15) through Eq. (28). While one-dimensional integrals can be evaluated directly with the proposed technique, multi-dimensional integrals are first transformed to nested one-dimensional integrals and then evaluated through the adaptive quadrature technique. Note that these methods can be applied for all problems involving Bayesian inference where the likelihood (L) can be expressed as a function of (implicit or explicit) the parameters being inferred and is not limited to the model calibration problem (Eq. (13)) in this paper. Bayesian updating often encounters evaluation of integrals of the form:

K¼

Z

6. Numerical example

Lðx; yÞf ðxÞf ðyÞdxdy

ð29Þ

The integral in Eq. (29) is a two-dimensional integral. However, by the inherent definition of likelihood, this can be converted into two nested one-dimensional integrals. Recalling the definition of likelihood, the likelihood of x and y, L(x, y) is proportional to the probability of observing data conditioned on x and y. Hence, the likelihood of x, L(x) can be calculated based on the principle of total probability as

LðxÞ / ProbðDjxÞ /

Z

ProbðDjx; yÞProbðyÞdy

ð30Þ

According to the definition of likelihood, Prob(D|x, y) is proportional to L(x, y) and hence Eq. (30) can be rewritten as:

LðxÞ /

Z

Lðx; yÞf ðyÞdy

ð31Þ

Substituting the expression in Eq. (31) into Eq. (29),

k¼

Z

LðxÞf ðxÞdx ¼

Z Z

 Lðx; yÞf ðyÞdy f ðxÞdx

ð32Þ

Thus, the two-dimensional integral in Eq. (29) has been converted into two nested one-dimensional integrals in Eq. (32). These one-dimensional integrals are evaluated using an advanced numerical algorithm known as Adaptive Recursive Simpson’s Quadrature [30]. Consider the one-dimensional integral and its approximation using Simpson’s rule as:

Z a

recursive quadrature algorithm requires about 50 evaluations of each nested integral. Hence, if there are only two parameters of interest, this requires about 2500 evaluations, which is extremely efficient in comparison with the MCMC technique. It is acknowledged that this method of integration (splitting the integrals and using adaptive recursive Simpson’s quadrature) is efficient only when the number of variables is small (less than 5). If the number of variables is greater than 5, then the MCMC sampling technique is more efficient than the proposed technique. In this paper, the objective is only to infer either a single variable (EIFS) or its distribution parameters. Most probability distributions have four parameters or less (only two in many cases, especially in EIFS calculations). Therefore, the numerical integration technique can be used in lieu of MCMC sampling to considerably reduce the computational expense. This method of integration is implemented in Section 6 to calculate the probability distribution of EIFS as well as the probability distributions of the distribution parameters (k, f) of the lognormal distribution of EIFS.

b

f ðxÞdx 

    ba aþb f ðaÞ þ 4f þ f ðbÞ ¼ Sða; bÞ 6 2

ð33Þ

The adaptive recursive quadrature algorithm calls for subdividing the interval of integration (a, b) into two sub-intervals ((a, c) and (c, b), a < c < b) and then, Simpson’s rule is applied to each sub-interval. The error in the estimate of the integral is calculated by comparing the integral values before and after splitting. The criterion for determining when to stop dividing a particular interval depends on the tolerance level e. The tolerance level for stopping [30] may be chosen as:

jSða; cÞ þ Sðc; bÞ  Sða; bÞj < 15e

ð34Þ

It is observed that while the MCMC sampling method requires about a million evaluations of the model G in Eq. (13) [15] to evaluate the integral in Eq. (29), the implementation of the adaptive

This section illustrates the three statistical methods proposed for the inference of EIFS through a numerical example. The following sub-sections present the details of the problem followed by inference results. 6.1. Description of the problem An embedded crack in a hollow cylindrical component is used to illustrate the proposed methodology. The component is subjected to multi-axial loading, namely tension, torsion, and bending. The details of the geometry and material properties are given in the next subsection. A cylindrical structure has been chosen for illustration purposes only. The stress intensity factor cannot be simply expressed as a function of crack size and loading conditions using a closed form function. A finite element model is constructed using the commercially available software ANSYS (version 11.0) and is used to generate training points (to predict Keqv) for the surrogate model. The methods in this paper can be applied to any practical structural component that requires extensive finite element analysis (or any similar expensive computation) to calculate the stress intensity factor. 6.2. Finite element model The finite element model of the cylindrical structural component with a circular embedded crack is shown in Fig. 7. In the finite element model, the crack faces are modeled as surface to surface contact elements in order to prevent the penetration of the crack’s upper and lower surfaces. The augmented Lagrangian [31] method is used for contact simulation. Additionally, friction effect is included in the material properties of the contact element, in which a Coulomb friction model is used. The Coulomb friction model defines an equivalent shear stress which is proportional to the contact pressure and the friction coefficient. Friction coefficients between two crack faces are difficult to measure and are generally assumed to vary between 0 and 0.5 [32]. The friction coefficient, l, used within this study is assumed to be 0.1. Since the primary quantity of interest is the stress intensity factor at the crack tip, the volume along the crack front is subdivided into many smaller blocks, which allows for better mesh control and enables SIF evaluation at various locations along the crack front.

1697

S. Sankararaman et al. / International Journal of Fatigue 32 (2010) 1689–1700 Table 3 Load spectra. Load type

Cycle magnitude

Load level (L)

Bending

Largest Medium Small

3 2 1

Torsion

Largest Medium Small

3 2 1

Tension

Largest Medium Small

3 2 1

sumed to be equally likely to occur within each block. It is essential to note that this particular loading history has been considered for illustration purposes only and any loading history can be analyzed using the proposed approach. The following section describes the calculation of stress intensity factor as a function of loading and crack size. 6.4. Calculation of stress intensity factor

Fig. 7. Embedded crack in a cylindrical structure.

The mesh around the crack location (at the crack front and surrounding areas) is refined in order to obtain a more accurate solution and avoid convergence problems. The stress intensity factor is calculated all along the crack front and used in crack growth analysis. The material properties and the geometrical properties of the cylindrical component are shown in Tables 1 and 2 respectively. 6.3. Load spectrum A variable-amplitude multi-axial loading scenario is considered in this paper. The cylinder shown in Fig. 7 is subjected to bending, torsion and tensile loads. Three different load combinations have been considered and numbered as L = 1, 2 and 3. As seen from Table 3, L = 1 represents the combination of smallest magnitudes in bending, torsion and tension while L = 3 represents the combination of largest magnitudes in bending, tension and torsion. A block loading history is shown in Fig. 8, where the number of load cycles in each block is treated as a random variable and any of the three load cases are as-

Table 1 Material properties. Material properties: 4340 Steel E N

ry rult

29000 ksi 0.32 163 ksi 200 ksi

Table 2 Geometrical properties.

A surrogate model is constructed to replace the finite element analysis explained in Section 6.2. Twenty six different combinations of crack size and loading levels were considered and the effective stress intensity factor was calculated through finite element analyses. First, a few standard response surfaces were tested and the Gaussian process model was selected after careful analysis and rigorous scrutiny. First, a conventional polynomial response surface (second order, without interaction terms) was used to fit the data. 10 data points were selected for fitting and 16 more data points were used to quantify the error in the surrogate model. The error sum of squares (sum taken over the remaining 16 points) was observed to be 58.1 units (measured in terms of the predicted stress intensity factor). Advanced response surface methods such as polynomial chaos expansion [23], support vector regression [24], relevance vector regression [25], and Gaussian Process interpolation [26] yielded error sum of squares values of 57.7 units, 55.1 units, 54.9 units, and 50.3 units respectively. Hence, the Gaussian process model (least sum of squares error) was selected for replacing the finite element model, and was trained again using all the available 26 points. The error analysis of Gaussian process interpolation is graphically shown in Fig. 9. 6.5. Simulation of experimental data For the sake of illustration, an experimental data set is simulated using an assumed EIFS distribution and crack growth according to modified Paris law including Wheeler’s retardation effects (Eq. (5)). The ‘true’ EIFS follows a lognormal distribution with k = 0.642 and f = 0.067 (corresponding to mean = 1.905 mm and standard deviation = 0.127 mm). Random samples of ‘true’ EIFS are generated, and using Eq. (5) experimental data points describing the number of life cycles (N) to reach a final crack size (A) are simulated. A noise term (10% Gaussian white noise) is added to each N so as to account for experimental errors. In order to examine the sensitivity of the calibrated EIFS to the data size, three different sets of data are simulated – 20 ‘‘measurements”, 10 ‘‘measurements” and 5 ‘‘measurements”.

Geometry Inside radius, ri Outside radius, ro Height, h Crack location

0.345 in. 0.568 in. 6 in. Mid-height

6.6. Statistical Inference The various methods of EIFS inference developed in Section 4 are now employed to estimate the statistics of EIFS. First, the dis-

1698

S. Sankararaman et al. / International Journal of Fatigue 32 (2010) 1689–1700

Fig. 8. Sample loading history.

25 Training Points GP Interpolation Validation Points

Keq (in Pa. m

1/ 2

)

20 15 10 5 0 -5 2.5

2

3

1.5

2.5 1

2 0.5

Crack Size (in mm)

1.5 0

1

Load Level

Fig. 9. Gaussian process surrogate model for DK = W(a, L).

tribution of EIFS (h) is calculated directly. Second, EIFS is assumed to follow a lognormal distribution (h  LN(k, f)) as done by Makeev et al. [15] and Cross et al. [16] and the probability distributions of k and f are calculated using the methods of maximum likelihood and Bayesian inference. The most likely values of these estimates are used to calculate the mean and standard deviation of h. The results (mean and standard deviation of h calculated using five different methods) are summarized in Table 4.

It is seen from Table 4 that all the methods give good estimates of the mean of equivalent initial flaw size. The maximum error is about 0.5% when inferred from 20 points, and this increases to 1.8% when inferred from 5 data points. This suggests that even with a small number of data points (in this case 5), a good estimate of the mean of EIFS can be obtained. However, the estimate of variance varies from method to method. It is observed that the standard deviation obtained using Method

Table 4 Summary of EIFS inference results. No.

Method

1 2 3 4 5

Likelihood of EIFS Bayesian inference of EIFS Least squares technique Likelihood of k and f Bayesian inference of k and f

Distribution of h (5 data points)

Distribution of h (10 data points)

Distribution of h (20 data points)

Mean (mm)

Standard deviation (mm)

Mean (mm)

Standard deviation (mm)

Mean (mm)

Standard deviation (mm)

1.901 1.900 1.897 1.936 1.934

0.069 0.071 0.209 0.201 0.198

1.895 1.899 1.911 1.889 1.891

0.045 0.048 0.105 0.179 0.175

1.907 1.906 1.905 1.910 1.913

0.025 0.021 0.069 0.128 0.128

Note: True distribution of EIFS assumed to be lognormal (mean = 1.905 mm and standard deviation = 0.127 mm).

S. Sankararaman et al. / International Journal of Fatigue 32 (2010) 1689–1700

1699

Fig. 10. Bayesian Inference.

1 through Method 3 decreases with increase in number of data points. It is difficult to calculate the ‘true’ standard deviation of EIFS, as the estimate is found keep on decreasing with more number of points. Method 4 and Method 5, i.e. inference on the distribution parameters of EIFS yield more accurate estimates of standard deviation. In particular, when using 5 data points, the standard deviation is higher than when using 20 data points. Hence, with more data points, the standard deviation can be expected to reach the ‘true’ standard deviation asymptotically. This is shown in Fig. 10, where the inference using 20 points is found to represent the ‘true’ distribution with reasonable accuracy.

7. Summary This study investigated different methods of statistical inference to estimate the distribution of equivalent initial flaw size (EIFS). The methods developed can be applied to structural systems with complicated geometry, subjected to multi-axial variable amplitude loading. The calculation of effective stress intensity factor is facilitated through the use of a Gaussian process model, which replaces repeated finite element analyses. The various inference techniques are based on the concepts of likelihood, Bayesian updating and least squares. The probability distribution of EIFS is obtained in each of these methods. The principles of likelihood and Bayesian inference can also be used to estimate the probability distributions of the distribution parameters of EIFS, assuming that EIFS follows a particular distribution type, with unknown parameters. The methods are illustrated using a numerical example. Experimental data is simulated using an assumed EIFS distribution and the various methods developed are employed to estimate this distribution. The inferred probability distributions are compared with the ‘true’ distribution and it is seen that satisfactory numerical results are obtained. Future work needs to explicitly consider multiple sources of uncertainty in the inference of EIFS. It was mentioned earlier that fatigue crack growth is a stochastic process with several uncertainties – physical variability, data uncertainty and modeling errors. In this paper, the contributions of all these different sources of uncertainty were combined into a single random variable (e) in Eq. (13). The study of individual contributions from different sources uncertainty would help in deciding trade-offs between model-refinement and experimental data collection, in order to reduce the overall uncertainty in probabilistic fatigue life prediction.

Acknowledgements The research reported in this paper was supported in part by the NASA ARMD/AvSP IVHM project under NRA Award No. NNX09AY54A (Project Manager: Dr. K. Goebel, NASA Ames Research Center) through subcontract to Clarkson University (No. 375-32531, Principal Investigator: Dr. Y. Liu). The support is gratefully acknowledged. References [1] Lankford J, Hudak SJ. Relevance of the small crack problem to lifetime prediction in gas turbines. Int J Fatigue 1987;9(2):87–93. [2] Liu Y, Mahadevan S. Probabilistic fatigue life prediction using an equivalent initial flaw size distribution. Int J Fatigue 2009;31(3):476–87. [3] Joint service specification guide aircraft structures, JSSG-2006. United States of America: Department of Defense; 1998. [4] Gallagher JP, Berens AP, Engle Jr RM. USAF damage tolerant design handbook: guidelines for the analysis and design of damage tolerant aircraft structures. Final report; 1984. [5] Merati A, Eastaugh G. Determination of fatigue related discontinuity state of 7000 series of aerospace aluminum alloys. Eng Failure Anal 2007;14(4): 673–85. [6] Heller RA, Yang JN. Statistical distribution of time to crack initiation and initial crack size using service data, final report NASA NSG.1099 (Supplement 1); 1977. [7] Moreira PMGP, de Matos PFP, de Castro PMST. Fatigue striation spacing and equivalent initial flaw size in Al 2024-T3 riveted specimens. Theor Appl Fract Mech 2005;43(1):89–99. [8] Barter SA, Sharp PK, Holden G, Clark G. Initiation and early growth of fatigue cracks in an aerospace aluminum alloy. Fatigue Fract Eng Mater 2002;25(2): 111–25. [9] White P, Molent L, Barter S. Interpreting fatigue test results using a probabilistic fracture approach. Int J Fatigue 2005;27(7):752–67. [10] Molent L, Sun Q, Green A. Characterization of equivalent initial flaw sizes in 7050 aluminium alloy. J Fatigue Fract Eng Mater Struct 2006;29:916–37. [11] Kitagawa H, Takahashi S. Applicability of fracture mechanics to vary small cracks or cracks in early stage. In: Proceedings of the 2nd international conference on mechanical behavior of materials. USA (OH): ASM International; 1976. [12] El Haddad MH, Topper TH, Smith KN. Prediction of nonpropagating cracks. Eng Fract Mech 1979;11:573–84. [13] Xiang Y, Lu Z, Liu Y. Crack growth-based fatigue life prediction using an equivalent initial flaw model. Part I: uniaxial loading. International Journal of Fatigue 2010;32(2):341–9. [14] Lu Z, Xiang Y, Liu Y. Crack growth-based fatigue life prediction using an equivalent initial flaw model. Part II: multiaxial loading. Intl J Fatigue 2010;32(2):376–81. [15] Makeev A, Nikishkov Y, Armanios E. A concept for quantifying equivalent initial flaw size distribution in fracture mechanics based life prediction models. Int J Fatigue 2006;29(1). [16] Cross R, Makeev A, Armanios E. Simultaneous uncertainty quantification of fracture mechanics based life prediction model parameters. Int J Fatigue 2007;29(8). [17] Haldar A, Mahadevan S. Probability, reliability and statistical methods in engineering design. New York: Wiley; 2000.

1700

S. Sankararaman et al. / International Journal of Fatigue 32 (2010) 1689–1700

[18] Yuen BKC, Taheri F. Proposed modifications to the Wheeler retardation model for multiple overloading fatigue life prediction. Int J Fatigue, 0142-1123 2006;28(12):1803–19. doi:10.1016/j.ijfatigue.2005.12.007. [19] Sheu BC, Song PS, Hwang S. Shaping exponent in wheeler model under a single overload. Eng Fract Mech 1995;51(1):135–43. [20] Song PS, Sheu BC, Chang L. A modified wheeler model to improve predictions of crack growth following a single overload. JSME Int J Ser A 2001;44(1): 117–22. [21] Liu Y, Mahadevan S. Multiaxial high-cycle fatigue criterion and life prediction for metals. Int J Fatigue 2005;27(7):790–800. [22] Liu Y, Mahadevan S. A unified multiaxial fatigue model for isotropic and anisotropic materials. Int J Fatigue 2007;29:347–59. [23] Ghanem R, Spanos J. Polynomial chaos in stochastic finite elements. Appl Mech 1990;57:197. doi:10.1115/1.2888303. [24] Boser B, Guyon I, Vapnik VN. A training algorithm for optimal margin classifiers. In: Proceedings of the Fifth annual workshop on computational learning theory; 1992. p. 144–52.

[25] Tipping M. Sparse Bayesian learning and the relevance vector machine. J Mach Learn Res 2001;1:211–44. [26] Rasmussen, C. Evaluation of Gaussian processes and other methods for nonlinear regression. PhD thesis, University of Toronto; 1996. [27] Santner TJ, Williams BJ, Noltz WI. The design and analysis of computer experiments. New York: Springer-Verlag; 2003. [28] McFarland J. Uncertainty analysis for computer simulations through validation and calibration. Ph D. Dissertation, Vanderbilt University, 2008. [29] Seber GAF, Wild CJ. Nonlinear regression. Wiley Series in Probability and Mathematical Statistics; 1989. [30] McKeeman MW. Algorithm 145: adaptive numerical integration by Simpson’s rule. Commun ACM 1962;5(12). [31] ANSYS. ANSYS theory reference, release 11.0. ANSYS Inc., 2007. [32] Liu Y, Liu L, Mahadevan S. Analysis of subsurface crack propagation under rolling contact loading in railroad wheels using FEM. Eng Fract Mech 2007;74:2659–74.